Line2 = -imgStep; address = 0;
LastLine = imgStep * (imageHeight - 1); while( ConvLine < imageHeight ) {
/*Here we calculate derivatives for line of image */ int memYline = (ConvLine + 1) & 1;
Line2 += imgStep;
Line1 = Line2 - ((Line2 == 0) ? 0 : imgStep);
Line3 = Line2 + ((Line2 == LastLine) ? 0 : imgStep);
/* Process first pixel */
ConvX = CONV( imgA[Line1 + 1], imgA[Line2 + 1], imgA[Line3 + 1] );
ConvY = CONV( imgA[Line3], imgA[Line3], imgA[Line3 + 1] );
GradY = (ConvY - MemY[memYline][0]) * 0.125f; GradX = (ConvX - MemX[1][ConvLine]) * 0.125f;
MemY[memYline][0] = ConvY; MemX[1][ConvLine] = ConvX;
GradT = (float) (imgB[Line2] - imgA[Line2]);
II[address].xx = GradX * GradX; II[address].xy = GradX * GradY; II[address].yy = GradY * GradY; II[address].xt = GradX * GradT; II[address].yt = GradY * GradT;
II[address].alpha = 1 / (Ilambda + II[address].xx + II[address].yy); address++;
/* Process middle of line */
for( j = 1; j < imageWidth - 1; j++ ) {
ConvX = CONV( imgA[Line1 + j + 1], imgA[Line2 + j + 1], imgA[Line3 + j + 1] );
ConvY = CONV( imgA[Line3 + j - 1], imgA[Line3 + j], imgA[Line3 + j + 1] );
GradY = (ConvY - MemY[memYline][j]) * 0.125f;
GradX = (ConvX - MemX[(j - 1) & 1][ConvLine]) * 0.125f;
MemY[memYline][j] = ConvY;
MemX[(j - 1) & 1][ConvLine] = ConvX;
GradT = (float) (imgB[Line2 + j] - imgA[Line2 + j]);
II[address].xx = GradX * GradX; II[address].xy = GradX * GradY; II[address].yy = GradY * GradY; II[address].xt = GradX * GradT; II[address].yt = GradY * GradT;
II[address].alpha = 1 / (Ilambda + II[address].xx + II[address].yy);
address++; }
/* Process last pixel of line */
ConvX = CONV( imgA[Line1 + imageWidth - 1], imgA[Line2 + imageWidth - 1],
imgA[Line3 + imageWidth - 1] );
ConvY = CONV( imgA[Line3 + imageWidth - 2], imgA[Line3 + imageWidth - 1],
imgA[Line3 + imageWidth - 1] );
GradY = (ConvY - MemY[memYline][imageWidth - 1]) * 0.125f; GradX = (ConvX - MemX[(imageWidth - 2) & 1][ConvLine]) * 0.125f;
MemY[memYline][imageWidth - 1] = ConvY;
GradT = (float) (imgB[Line2 + imageWidth - 1] - imgA[Line2 + imageWidth - 1]);
II[address].xx = GradX * GradX; II[address].xy = GradX * GradY; II[address].yy = GradY * GradY; II[address].xt = GradX * GradT; II[address].yt = GradY * GradT;
II[address].alpha = 1 / (Ilambda + II[address].xx + II[address].yy); address++;
ConvLine++; }
/****************************************************************************************\\ * Prepare initial approximation *
\\****************************************************************************************/ if( !usePrevious ) {
float *vx = velocityX; float *vy = velocityY;
for( i = 0; i < imageHeight; i++ ) {
memset( vx, 0, imageWidth * sizeof( float )); memset( vy, 0, imageWidth * sizeof( float ));
vx += velStep; vy += velStep; } }
/****************************************************************************************\\ * Perform iterations *
\\****************************************************************************************/ iter = 0; Stop = 0;
LastLine = velStep * (imageHeight - 1); while( !Stop ) {
float Eps = 0; address = 0;
iter++;
/****************************************************************************************\\ * begin scan velocity and update it *
\\****************************************************************************************/ Line2 = -velStep;
for( i = 0; i < imageHeight; i++ ) {
/* Here average velocity */
float averageX; float averageY; float tmp;
Line2 += velStep;
Line1 = Line2 - ((Line2 == 0) ? 0 : velStep);
Line3 = Line2 + ((Line2 == LastLine) ? 0 : velStep); /* Process first pixel */
averageX = (velocityX[Line2] +
velocityX[Line2 + 1] + velocityX[Line1] + velocityX[Line3]) / 4;
averageY = (velocityY[Line2] +
velocityY[Line2 + 1] + velocityY[Line1] + velocityY[Line3]) / 4;
VelBufX[i & 1][0] = averageX - (II[address].xx * averageX +
II[address].xy * averageY + II[address].xt) * II[address].alpha;
VelBufY[i & 1][0] = averageY - (II[address].xy * averageX +
II[address].yy * averageY + II[address].yt) * II[address].alpha;
/* update Epsilon */
if( criteria.type & CV_TERMCRIT_EPS ) {
tmp = (float)fabs(velocityX[Line2] - VelBufX[i & 1][0]); Eps = MAX( tmp, Eps );
tmp = (float)fabs(velocityY[Line2] - VelBufY[i & 1][0]); Eps = MAX( tmp, Eps ); }
address++;
/* Process middle of line */
for( j = 1; j < imageWidth - 1; j++ ) {
averageX = (velocityX[Line2 + j - 1] + velocityX[Line2 + j + 1] +
velocityX[Line1 + j] + velocityX[Line3 + j]) / 4;
averageY = (velocityY[Line2 + j - 1] + velocityY[Line2 + j + 1] +
velocityY[Line1 + j] + velocityY[Line3 + j]) / 4;
VelBufX[i & 1][j] = averageX - (II[address].xx * averageX +
II[address].xy * averageY + II[address].xt) * II[address].alpha;
VelBufY[i & 1][j] = averageY - (II[address].xy * averageX +
II[address].yy * averageY + II[address].yt) * II[address].alpha;
/* update Epsilon */
if( criteria.type & CV_TERMCRIT_EPS ) {
tmp = (float)fabs(velocityX[Line2 + j] - VelBufX[i & 1][j]);
Eps = MAX( tmp, Eps );
tmp = (float)fabs(velocityY[Line2 + j] - VelBufY[i & 1][j]);
Eps = MAX( tmp, Eps ); }
address++; }
/* Process last pixel of line */
averageX = (velocityX[Line2 + imageWidth - 2] + velocityX[Line2 + imageWidth - 1] + velocityX[Line1 + imageWidth - 1] + velocityX[Line3 + imageWidth - 1]) / 4;
averageY = (velocityY[Line2 + imageWidth - 2] + velocityY[Line2 + imageWidth - 1] + velocityY[Line1 + imageWidth - 1] + velocityY[Line3 + imageWidth - 1]) / 4;
VelBufX[i & 1][imageWidth - 1] = averageX - (II[address].xx * averageX +
II[address].xy * averageY + II[address].xt) * II[address].alpha;