📄 cvoptflowlk.cpp
字号:
GradY = ConvY - MemY[memYline][0];
GradX = ConvX - MemX[1][L2];
MemY[memYline][0] = ConvY;
MemX[1][L2] = 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;
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];
GradX = ConvX - MemX[(j - 1) & 1][L2];
MemY[memYline][j] = ConvY;
MemX[(j - 1) & 1][L2] = 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;
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];
GradX = ConvX - MemX[(imageWidth - 2) & 1][L2];
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;
address++;
/* End of derivatives for line */
/****************************************************************************************/
/* ---------Calculating horizontal convolution of processed line----------------------- */
/****************************************************************************************/
address -= BufferWidth;
/* process first HorRadius pixels */
for( j = 0; j < HorRadius; j++ )
{
int jj;
WII[address].xx = 0;
WII[address].xy = 0;
WII[address].yy = 0;
WII[address].xt = 0;
WII[address].yt = 0;
for( jj = -j; jj <= HorRadius; jj++ )
{
float Ker = KerX[jj];
WII[address].xx += II[address + jj].xx * Ker;
WII[address].xy += II[address + jj].xy * Ker;
WII[address].yy += II[address + jj].yy * Ker;
WII[address].xt += II[address + jj].xt * Ker;
WII[address].yt += II[address + jj].yt * Ker;
}
address++;
}
/* process inner part of line */
for( j = HorRadius; j < imageWidth - HorRadius; j++ )
{
int jj;
float Ker0 = KerX[0];
WII[address].xx = 0;
WII[address].xy = 0;
WII[address].yy = 0;
WII[address].xt = 0;
WII[address].yt = 0;
for( jj = 1; jj <= HorRadius; jj++ )
{
float Ker = KerX[jj];
WII[address].xx += (II[address - jj].xx + II[address + jj].xx) * Ker;
WII[address].xy += (II[address - jj].xy + II[address + jj].xy) * Ker;
WII[address].yy += (II[address - jj].yy + II[address + jj].yy) * Ker;
WII[address].xt += (II[address - jj].xt + II[address + jj].xt) * Ker;
WII[address].yt += (II[address - jj].yt + II[address + jj].yt) * Ker;
}
WII[address].xx += II[address].xx * Ker0;
WII[address].xy += II[address].xy * Ker0;
WII[address].yy += II[address].yy * Ker0;
WII[address].xt += II[address].xt * Ker0;
WII[address].yt += II[address].yt * Ker0;
address++;
}
/* process right side */
for( j = imageWidth - HorRadius; j < imageWidth; j++ )
{
int jj;
WII[address].xx = 0;
WII[address].xy = 0;
WII[address].yy = 0;
WII[address].xt = 0;
WII[address].yt = 0;
for( jj = -HorRadius; jj < imageWidth - j; jj++ )
{
float Ker = KerX[jj];
WII[address].xx += II[address + jj].xx * Ker;
WII[address].xy += II[address + jj].xy * Ker;
WII[address].yy += II[address + jj].yy * Ker;
WII[address].xt += II[address + jj].xt * Ker;
WII[address].yt += II[address + jj].yt * Ker;
}
address++;
}
}
/****************************************************************************************/
/* Calculating velocity line */
/****************************************************************************************/
if( PixelLine >= 0 )
{
int USpace;
int BSpace;
int address;
if( PixelLine < VerRadius )
USpace = PixelLine;
else
USpace = VerRadius;
if( PixelLine >= imageHeight - VerRadius )
BSpace = imageHeight - PixelLine - 1;
else
BSpace = VerRadius;
address = ((PixelLine - USpace) % BufferHeight) * BufferWidth;
for( j = 0; j < imageWidth; j++ )
{
int addr = address;
A1B2 = 0;
A2 = 0;
B1 = 0;
C1 = 0;
C2 = 0;
for( i = -USpace; i <= BSpace; i++ )
{
A2 += WII[addr + j].xx * KerY[i];
A1B2 += WII[addr + j].xy * KerY[i];
B1 += WII[addr + j].yy * KerY[i];
C2 += WII[addr + j].xt * KerY[i];
C1 += WII[addr + j].yt * KerY[i];
addr += BufferWidth;
addr -= ((addr >= BufferSize) ? 0xffffffff : 0) & BufferSize;
}
/****************************************************************************************\
* Solve Linear System *
\****************************************************************************************/
{
float delta = (A1B2 * A1B2 - A2 * B1);
if( delta )
{
/* system is not singular - solving by Kramer method */
float deltaX;
float deltaY;
float Idelta = 8 / delta;
deltaX = -(C1 * A1B2 - C2 * B1);
deltaY = -(A1B2 * C2 - A2 * C1);
velocityX[j] = deltaX * Idelta;
velocityY[j] = deltaY * Idelta;
}
else
{
/* singular system - find optical flow in gradient direction */
float Norm = (A1B2 + A2) * (A1B2 + A2) + (B1 + A1B2) * (B1 + A1B2);
if( Norm )
{
float IGradNorm = 8 / Norm;
float temp = -(C1 + C2) * IGradNorm;
velocityX[j] = (A1B2 + A2) * temp;
velocityY[j] = (B1 + A1B2) * temp;
}
else
{
velocityX[j] = 0;
velocityY[j] = 0;
}
}
}
/****************************************************************************************\
* End of Solving Linear System *
\****************************************************************************************/
} /*for */
velocityX += velStep;
velocityY += velStep;
} /*for */
PixelLine++;
ConvLine++;
}
/* Free memory */
for( k = 0; k < 2; k++ )
{
cvFree( &MemX[k] );
cvFree( &MemY[k] );
}
cvFree( &II );
cvFree( &WII );
return CV_OK;
} /*icvCalcOpticalFlowLK_8u32fR*/
/*F///////////////////////////////////////////////////////////////////////////////////////
// Name: cvCalcOpticalFlowLK
// Purpose: Optical flow implementation
// Context:
// Parameters:
// srcA, srcB - source image
// velx, vely - destination image
// Returns:
//
// Notes:
//F*/
CV_IMPL void
cvCalcOpticalFlowLK( const void* srcarrA, const void* srcarrB, CvSize winSize,
void* velarrx, void* velarry )
{
CV_FUNCNAME( "cvCalcOpticalFlowLK" );
__BEGIN__;
CvMat stubA, *srcA = (CvMat*)srcarrA;
CvMat stubB, *srcB = (CvMat*)srcarrB;
CvMat stubx, *velx = (CvMat*)velarrx;
CvMat stuby, *vely = (CvMat*)velarry;
CV_CALL( srcA = cvGetMat( srcA, &stubA ));
CV_CALL( srcB = cvGetMat( srcB, &stubB ));
CV_CALL( velx = cvGetMat( velx, &stubx ));
CV_CALL( vely = cvGetMat( vely, &stuby ));
if( !CV_ARE_TYPES_EQ( srcA, srcB ))
CV_ERROR( CV_StsUnmatchedFormats, "Source images have different formats" );
if( !CV_ARE_TYPES_EQ( velx, vely ))
CV_ERROR( CV_StsUnmatchedFormats, "Destination images have different formats" );
if( !CV_ARE_SIZES_EQ( srcA, srcB ) ||
!CV_ARE_SIZES_EQ( velx, vely ) ||
!CV_ARE_SIZES_EQ( srcA, velx ))
CV_ERROR( CV_StsUnmatchedSizes, "" );
if( CV_MAT_TYPE( srcA->type ) != CV_8UC1 ||
CV_MAT_TYPE( velx->type ) != CV_32FC1 )
CV_ERROR( CV_StsUnsupportedFormat, "Source images must have 8uC1 type and "
"destination images must have 32fC1 type" );
if( srcA->step != srcB->step || velx->step != vely->step )
CV_ERROR( CV_BadStep, "source and destination images have different step" );
IPPI_CALL( icvCalcOpticalFlowLK_8u32fR( (uchar*)srcA->data.ptr, (uchar*)srcB->data.ptr,
srcA->step, cvGetMatSize( srcA ), winSize,
velx->data.fl, vely->data.fl, velx->step ));
__END__;
}
/* End of file. */
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -