Skip to content

O
opencv_contrib
Commit 65d01574 authored by Vladislav Samsonov's avatar Vladislav Samsonov
Browse files
Options

Some improvements

parent 25b2958e
Hide whitespace changes
Inline Side-by-side
Showing with 260 additions and 153 deletions
modules/optflow/include/opencv2/optflow/pcaflow.hpp
View file @ 65d01574
......@@ -47,6 +47,123 @@ namespace cv
{
namespace optflow
{
/*
class PCAFlowBasis
{
public:
Size size;
PCAFlowBasis( Size basisSize = Size( 0, 0 ) ) : size( basisSize ) {}
virtual ~PCAFlowBasis(){};
virtual int getNumberOfComponents() const = 0;
virtual void getBasisAtPoint( const Point2f &p, const Size &maxSize, float *outX, float *outY ) const = 0;
virtual Point2f reduceAtPoint( const Point2f &p, const Size &maxSize, const float *w1, const float *w2 ) const = 0;
};*/
/*
* Orthogonal basis from Discrete Cosine Transform.
* Can be used without any learning or assumptions about flow structure for general purpose.
* Gives low quality estimation.
*/
/*class PCAFlowGeneralBasis : public PCAFlowBasis
{
public:
PCAFlowGeneralBasis( Size basisSize = Size( 18, 14 ) ) : PCAFlowBasis( basisSize ) {}
int getNumberOfComponents() const { return size.area(); }
void getBasisAtPoint( const Point2f &p, const Size &maxSize, float *outX, float *outY ) const
{
for ( int n1 = 0; n1 < size.width; ++n1 )
for ( int n2 = 0; n2 < size.height; ++n2 )
outX[n1 * size.height + n2] =
cosf( ( n1 * M_PI / maxSize.width ) * ( p.x + 0.5 ) ) * cosf( ( n2 * M_PI / maxSize.height ) * ( p.y + 0.5 )
);
memcpy( outY, outX, getNumberOfComponents() * sizeof( *outY ) );
}
Point2f reduceAtPoint( const Point2f &p, const Size &maxSize, const float *w1, const float *w2 ) const
{
Point2f res( 0, 0 );
for ( int n1 = 0; n1 < size.width; ++n1 )
for ( int n2 = 0; n2 < size.height; ++n2 )
{
const float c =
cosf( ( n1 * M_PI / maxSize.width ) * ( p.x + 0.5 ) ) * cosf( ( n2 * M_PI / maxSize.height ) * ( p.y + 0.5 )
);
res.x += c * w1[n1 * size.height + n2];
res.y += c * w2[n1 * size.height + n2];
}
return res;
}
};*/
/*
class PCAFlowLearnedBasis : public PCAFlowBasis
{
private:
float *basisData;
unsigned numberOfComponents;
public:
PCAFlowLearnedBasis( const char *filename )
{
basisData = 0;
FILE *f = fopen( filename, "r" );
CV_Assert( f );
numberOfComponents = 0;
CV_Assert( fread( &numberOfComponents, sizeof( numberOfComponents ), 1, f ) == 1 );
CV_Assert( fread( &size.height, sizeof( size.height ), 1, f ) == 1 );
CV_Assert( fread( &size.width, sizeof( size.width ), 1, f ) == 1 );
CV_Assert( ( numberOfComponents > 0 ) && ( numberOfComponents % 2 == 0 ) );
basisData = new float[size.width * size.height * numberOfComponents];
CV_Assert( fread( basisData, size.width * size.height * sizeof( *basisData ), numberOfComponents, f ) ==
numberOfComponents );
fclose( f );
numberOfComponents /= 2;
}
~PCAFlowLearnedBasis()
{
if ( basisData )
delete[] basisData;
}
int getNumberOfComponents() const { return numberOfComponents; }
void getBasisAtPoint( const Point2f &p, const Size &maxSize, float *outX, float *outY ) const
{
const size_t chunk = size.width * size.height;
size_t offset = size_t( p.y * float(size.height) / maxSize.height ) * size.width + size_t( p.x * float(size.width) /
maxSize.width );
for ( unsigned i = 0; i < numberOfComponents; ++i )
outX[i] = basisData[i * chunk + offset];
offset += numberOfComponents * chunk;
for ( unsigned i = 0; i < numberOfComponents; ++i )
outY[i] = basisData[i * chunk + offset];
}
Point2f reduceAtPoint( const Point2f &p, const Size &maxSize, const float *w1, const float *w2 ) const
{
Point2f res( 0, 0 );
const size_t chunk = size.width * size.height;
const size_t offset = size_t( p.y * float(size.height) / maxSize.height ) * size.width + size_t( p.x *
float(size.width) / maxSize.width );
for ( unsigned i = 0; i < numberOfComponents; ++i )
{
const float c = basisData[i * chunk + offset];
res.x += c * w1[i];
res.y += c * w2[i];
}
return res;
}
};*/
class OpticalFlowPCAFlow : public DenseOpticalFlow
{
......@@ -57,8 +174,8 @@ protected:
const float occlusionsThreshold;
public:
OpticalFlowPCAFlow( Size _basisSize = Size( 18, 14 ), float _sparseRate = 0.02, float _retainedCornersFraction = 1.0,
float _occlusionsThreshold = 0.00002 );
OpticalFlowPCAFlow( const Size _basisSize = Size( 18, 14 ), float _sparseRate = 0.02,
float _retainedCornersFraction = 1.0, float _occlusionsThreshold = 0.00002 );
void calc( InputArray I0, InputArray I1, InputOutputArray flow );
void collectGarbage();
......
modules/optflow/src/pcaflow.cpp
View file @ 65d01574
......@@ -50,7 +50,7 @@ namespace cv
namespace optflow
{
OpticalFlowPCAFlow::OpticalFlowPCAFlow( Size _basisSize, float _sparseRate, float _retainedCornersFraction,
OpticalFlowPCAFlow::OpticalFlowPCAFlow( const Size _basisSize, float _sparseRate, float _retainedCornersFraction,
float _occlusionsThreshold )
: basisSize( _basisSize ), sparseRate( _sparseRate ), retainedCornersFraction( _retainedCornersFraction ),
occlusionsThreshold( _occlusionsThreshold )
......@@ -69,95 +69,6 @@ inline float eDistSq( const Point2f &p1, const Point2f &p2 )
inline float eNormSq( const Point2f &v ) { return v.x * v.x + v.y * v.y; }
void OpticalFlowPCAFlow::findSparseFeatures( Mat &from, Mat &to, std::vector<Point2f> &features,
std::vector<Point2f> &predictedFeatures ) const
{
Size size = from.size();
const unsigned maxFeatures = size.area() * sparseRate;
goodFeaturesToTrack( from, features, maxFeatures * retainedCornersFraction, 0.005, 3 );
// Add points along the grid if not enough features
if ( maxFeatures > features.size() )
{
const unsigned missingPoints = maxFeatures - features.size();
const unsigned blockSize = sqrt( (float)size.area() / missingPoints );
for ( int x = blockSize / 2; x < size.width; x += blockSize )
for ( int y = blockSize / 2; y < size.height; y += blockSize )
features.push_back( Point2f( x, y ) );
}
std::vector<uchar> predictedStatus;
std::vector<float> predictedError;
calcOpticalFlowPyrLK( from, to, features, predictedFeatures, predictedStatus, predictedError );
size_t j = 0;
for ( size_t i = 0; i < features.size(); ++i )
{
if ( predictedStatus[i] )
{
features[j] = features[i];
predictedFeatures[j] = predictedFeatures[i];
++j;
}
}
features.resize( j );
predictedFeatures.resize( j );
}
void OpticalFlowPCAFlow::removeOcclusions( Mat &from, Mat &to, std::vector<Point2f> &features,
std::vector<Point2f> &predictedFeatures ) const
{
std::vector<uchar> predictedStatus;
std::vector<float> predictedError;
std::vector<Point2f> backwardFeatures;
calcOpticalFlowPyrLK( to, from, predictedFeatures, backwardFeatures, predictedStatus, predictedError );
size_t j = 0;
const float threshold = occlusionsThreshold * from.size().area();
for ( size_t i = 0; i < predictedFeatures.size(); ++i )
{
if ( predictedStatus[i] )
{
Point2f flowDiff = features[i] - backwardFeatures[i];
if ( eNormSq( flowDiff ) < threshold )
{
features[j] = features[i];
predictedFeatures[j] = predictedFeatures[i];
++j;
}
}
}
features.resize( j );
predictedFeatures.resize( j );
}
void OpticalFlowPCAFlow::getSystem( OutputArray AOut, OutputArray b1Out, OutputArray b2Out,
const std::vector<Point2f> &features, const std::vector<Point2f> &predictedFeatures,
const Size size )
{
AOut.create( features.size(), basisSize.area(), CV_32F );
b1Out.create( features.size(), 1, CV_32F );
b2Out.create( features.size(), 1, CV_32F );
Mat A = AOut.getMat();
Mat b1 = b1Out.getMat();
Mat b2 = b2Out.getMat();
const Point2f scale =
Point2f( (float)basisSize.width / (float)size.width, (float)basisSize.height / (float)size.height );
for ( size_t i = 0; i < features.size(); ++i )
{
const Point2f p = Point2f( features[i].x * scale.x, features[i].y * scale.y );
for ( int n1 = 0; n1 < basisSize.width; ++n1 )
for ( int n2 = 0; n2 < basisSize.height; ++n2 )
{
const float c = cos( ( n1 * M_PI / basisSize.width ) * ( p.x + 0.5 ) ) *
cos( ( n2 * M_PI / basisSize.height ) * ( p.y + 0.5 ) );
A.at<float>( i, n1 * basisSize.height + n2 ) = c;
}
const Point2f flow = predictedFeatures[i] - features[i];
b1.at<float>( i ) = flow.x;
b2.at<float>( i ) = flow.y;
}
}
template <typename T> static inline int mathSign( T val ) { return ( T( 0 ) < val ) - ( val < T( 0 ) ); }
static inline void symOrtho( double a, double b, double &c, double &s, double &r )
......@@ -200,15 +111,6 @@ static void solveLSQR( const Mat &A, const Mat &b, OutputArray xOut, const doubl
CV_Assert( b.type() == CV_32F );
xOut.create( n, 1, CV_32F );
double anorm = 0;
const double dampsq = damp * damp;
double ddnorm = 0;
double res2 = 0;
double xxnorm = 0;
double z = 0;
double cs2 = -1;
double sn2 = 0;
Mat v( n, 1, CV_32F, 0.0f );
Mat u = b;
Mat x = xOut.getMat();
......@@ -216,11 +118,12 @@ static void solveLSQR( const Mat &A, const Mat &b, OutputArray xOut, const doubl
double alfa = 0;
double beta = cv::norm( u, NORM_L2 );
Mat w( n, 1, CV_32F, 0.0f );
const Mat AT = A.t();
if ( beta > 0 )
{
u *= 1 / beta;
v = A.t() * u;
v = AT * u;
alfa = cv::norm( v, NORM_L2 );
}
......@@ -232,10 +135,7 @@ static void solveLSQR( const Mat &A, const Mat &b, OutputArray xOut, const doubl
double rhobar = alfa;
double phibar = beta;
double rnorm = beta;
double r1norm = rnorm;
double arnorm = alfa * beta;
if ( arnorm == 0 )
if ( alfa * beta == 0 )
return;
for ( unsigned itn = 0; itn < iter_lim; ++itn )
......@@ -246,8 +146,7 @@ static void solveLSQR( const Mat &A, const Mat &b, OutputArray xOut, const doubl
if ( beta > 0 )
{
u *= 1 / beta;
anorm = sqrt( anorm * anorm + alfa * alfa + beta * beta + damp * damp );
v = A.t() * u - beta * v;
v = AT * u - beta * v;
alfa = cv::norm( v, NORM_L2 );
if ( alfa > 0 )
v = ( 1 / alfa ) * v;
......@@ -255,8 +154,6 @@ static void solveLSQR( const Mat &A, const Mat &b, OutputArray xOut, const doubl
double rhobar1 = sqrt( rhobar * rhobar + damp * damp );
double cs1 = rhobar / rhobar1;
double sn1 = damp / rhobar1;
double psi = sn1 * phibar;
phibar = cs1 * phibar;
double cs, sn, rho;
......@@ -266,37 +163,140 @@ static void solveLSQR( const Mat &A, const Mat &b, OutputArray xOut, const doubl
rhobar = -cs * alfa;
double phi = cs * phibar;
phibar = sn * phibar;
double tau = sn * phi;
double t1 = phi / rho;
double t2 = -theta / rho;
Mat dk = ( 1 / rho ) * w;
x = x + t1 * w;
w = v + t2 * w;
ddnorm += cv::norm( dk, NORM_L2SQR );
double delta = sn2 * rho;
double gambar = -cs2 * rho;
double rhs = phi - delta * z;
double gamma = sqrt( gambar * gambar + theta * theta );
cs2 = gambar / gamma;
sn2 = theta / gamma;
z = rhs / gamma;
xxnorm = xxnorm + z * z;
double res1 = phibar * phibar;
res2 = res2 + psi * psi;
rnorm = sqrt( res1 + res2 );
arnorm = alfa * std::abs( tau );
double r1sq = rnorm * rnorm - dampsq * xxnorm;
r1norm = sqrt( std::abs( r1sq ) );
if ( r1sq < 0 )
r1norm = -r1norm;
x += t1 * w;
w *= t2;
w += v;
}
}
void OpticalFlowPCAFlow::findSparseFeatures( Mat &from, Mat &to, std::vector<Point2f> &features,
std::vector<Point2f> &predictedFeatures ) const
{
Size size = from.size();
const unsigned maxFeatures = size.area() * sparseRate;
goodFeaturesToTrack( from, features, maxFeatures * retainedCornersFraction, 0.005, 3 );
// Add points along the grid if not enough features
if ( maxFeatures > features.size() )
{
const unsigned missingPoints = maxFeatures - features.size();
const unsigned blockSize = sqrt( (float)size.area() / missingPoints );
for ( int x = blockSize / 2; x < size.width; x += blockSize )
for ( int y = blockSize / 2; y < size.height; y += blockSize )
features.push_back( Point2f( x, y ) );
}
std::vector<uchar> predictedStatus;
std::vector<float> predictedError;
calcOpticalFlowPyrLK( from, to, features, predictedFeatures, predictedStatus, predictedError );
size_t j = 0;
for ( size_t i = 0; i < features.size(); ++i )
{
if ( predictedStatus[i] )
{
features[j] = features[i];
predictedFeatures[j] = predictedFeatures[i];
++j;
}
}
features.resize( j );
predictedFeatures.resize( j );
}
void OpticalFlowPCAFlow::removeOcclusions( Mat &from, Mat &to, std::vector<Point2f> &features,
std::vector<Point2f> &predictedFeatures ) const
{
std::vector<uchar> predictedStatus;
std::vector<float> predictedError;
std::vector<Point2f> backwardFeatures;
calcOpticalFlowPyrLK( to, from, predictedFeatures, backwardFeatures, predictedStatus, predictedError );
size_t j = 0;
const float threshold = occlusionsThreshold * from.size().area();
for ( size_t i = 0; i < predictedFeatures.size(); ++i )
{
if ( predictedStatus[i] )
{
Point2f flowDiff = features[i] - backwardFeatures[i];
if ( eNormSq( flowDiff ) < threshold )
{
features[j] = features[i];
predictedFeatures[j] = predictedFeatures[i];
++j;
}
}
}
features.resize( j );
predictedFeatures.resize( j );
}
void OpticalFlowPCAFlow::getSystem( OutputArray AOut, OutputArray b1Out, OutputArray b2Out,
const std::vector<Point2f> &features, const std::vector<Point2f> &predictedFeatures,
const Size size )
{
AOut.create( features.size(), basisSize.area(), CV_32F );
b1Out.create( features.size(), 1, CV_32F );
b2Out.create( features.size(), 1, CV_32F );
Mat A = AOut.getMat();
Mat b1 = b1Out.getMat();
Mat b2 = b2Out.getMat();
for ( size_t i = 0; i < features.size(); ++i )
{
const Point2f &p = features[i];
float *row = A.ptr<float>( i );
for ( int n1 = 0; n1 < basisSize.width; ++n1 )
for ( int n2 = 0; n2 < basisSize.height; ++n2 )
row[n1 * basisSize.height + n2] =
cosf( ( n1 * M_PI / size.width ) * ( p.x + 0.5 ) ) * cosf( ( n2 * M_PI / size.height ) * ( p.y + 0.5 ) );
const Point2f flow = predictedFeatures[i] - features[i];
b1.at<float>( i ) = flow.x;
b2.at<float>( i ) = flow.y;
}
}
static void applyCLAHE( Mat &img )
{
Ptr<CLAHE> clahe = createCLAHE();
clahe->setClipLimit( 8 );
clahe->apply( img, img );
}
static void reduceToFlow( const Mat &w1, const Mat &w2, Mat &flow, const Size &basisSize )
{
const Size size = flow.size();
Mat flowX( size, CV_32F, 0.0f );
Mat flowY( size, CV_32F, 0.0f );
const float mult = sqrt( size.area() ) * 0.5;
for ( int i = 0; i < basisSize.width; ++i )
for ( int j = 0; j < basisSize.height; ++j )
{
flowX.at<float>( j, i ) = w1.at<float>( i * basisSize.height + j ) * mult;
flowY.at<float>( j, i ) = w2.at<float>( i * basisSize.height + j ) * mult;
}
for ( int i = 0; i < basisSize.height; ++i )
{
flowX.at<float>( i, 0 ) *= M_SQRT2;
flowY.at<float>( i, 0 ) *= M_SQRT2;
}
for ( int i = 0; i < basisSize.width; ++i )
{
flowX.at<float>( 0, i ) *= M_SQRT2;
flowY.at<float>( 0, i ) *= M_SQRT2;
}
dct( flowX, flowX, DCT_INVERSE );
dct( flowY, flowY, DCT_INVERSE );
for ( int i = 0; i < size.height; ++i )
for ( int j = 0; j < size.width; ++j )
flow.at<Point2f>( i, j ) = Point2f( flowX.at<float>( i, j ), flowY.at<float>( i, j ) );
}
void OpticalFlowPCAFlow::calc( InputArray I0, InputArray I1, InputOutputArray flowOut )
{
const Size size = I0.size();
......@@ -325,6 +325,9 @@ void OpticalFlowPCAFlow::calc( InputArray I0, InputArray I1, InputOutputArray fl
CV_Assert( from.channels() == 1 );
CV_Assert( to.channels() == 1 );
// applyCLAHE(from);
// applyCLAHE(to);
std::vector<Point2f> features, predictedFeatures;
findSparseFeatures( from, to, features, predictedFeatures );
removeOcclusions( from, to, features, predictedFeatures );
......@@ -340,26 +343,13 @@ void OpticalFlowPCAFlow::calc( InputArray I0, InputArray I1, InputOutputArray fl
Mat A, b1, b2, w1, w2;
getSystem( A, b1, b2, features, predictedFeatures, size );
// solve( A, b1, w1, DECOMP_CHOLESKY | DECOMP_NORMAL );
// solve( A, b2, w2, DECOMP_CHOLESKY | DECOMP_NORMAL );
// solve( A1, b1, w1, DECOMP_CHOLESKY | DECOMP_NORMAL );
// solve( A2, b2, w2, DECOMP_CHOLESKY | DECOMP_NORMAL );
solveLSQR( A, b1, w1, 2 );
solveLSQR( A, b2, w2, 2 );
Mat flowSmall( basisSize, CV_32FC2 );
for ( int y = 0; y < basisSize.height; ++y )
for ( int x = 0; x < basisSize.width; ++x )
{
float sumX = 0, sumY = 0;
for ( int n1 = 0; n1 < basisSize.width; ++n1 )
for ( int n2 = 0; n2 < basisSize.height; ++n2 )
{
const float c = cos( ( n1 * M_PI / basisSize.width ) * ( x + 0.5 ) ) *
cos( ( n2 * M_PI / basisSize.height ) * ( y + 0.5 ) );
sumX += c * w1.at<float>( n1 * basisSize.height + n2 );
sumY += c * w2.at<float>( n1 * basisSize.height + n2 );
}
flowSmall.at<Point2f>( y, x ) = Point2f( sumX, sumY );
}
resize( flowSmall, flow, size, 0, 0, INTER_CUBIC );
Mat flowSmall( basisSize * 16, CV_32FC2 );
reduceToFlow( w1, w2, flowSmall, basisSize );
resize( flowSmall, flow, size, 0, 0, INTER_LINEAR );
}
void OpticalFlowPCAFlow::collectGarbage() {}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment