1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
5 // By downloading, copying, installing or using the software you agree to this license.
6 // If you do not agree to this license, do not download, install,
7 // copy or use the software.
8 //
9 //
10 // License Agreement
11 // For Open Source Computer Vision Library
12 //
13 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
14 // Copyright (C) 2009, Willow Garage Inc., all rights reserved.
15 // Third party copyrights are property of their respective owners.
16 //
17 // Redistribution and use in source and binary forms, with or without modification,
18 // are permitted provided that the following conditions are met:
19 //
20 // * Redistribution's of source code must retain the above copyright notice,
21 // this list of conditions and the following disclaimer.
22 //
23 // * Redistribution's in binary form must reproduce the above copyright notice,
24 // this list of conditions and the following disclaimer in the documentation
25 // and/or other materials provided with the distribution.
26 //
27 // * The name of the copyright holders may not be used to endorse or promote products
28 // derived from this software without specific prior written permission.
29 //
30 // This software is provided by the copyright holders and contributors "as is" and
31 // any express or implied warranties, including, but not limited to, the implied
32 // warranties of merchantability and fitness for a particular purpose are disclaimed.
33 // In no event shall the Intel Corporation or contributors be liable for any direct,
34 // indirect, incidental, special, exemplary, or consequential damages
35 // (including, but not limited to, procurement of substitute goods or services;
36 // loss of use, data, or profits; or business interruption) however caused
37 // and on any theory of liability, whether in contract, strict liability,
38 // or tort (including negligence or otherwise) arising in any way out of
39 // the use of this software, even if advised of the possibility of such damage.
40 //
41 //M*/
42
43 #include "precomp.hpp"
44 #include "opencv2/ximgproc/edge_filter.hpp"
45
46 /* Disable "from double to float" and "from size_t to int" warnings.
47 * Fixing these would make the code look ugly by introducing explicit cast all around.
48 * Here these warning are pointless anyway.
49 */
50 #ifdef _MSC_VER
51 #pragma warning( disable : 4305 4244 4267 4838 )
52 #endif
53 #ifdef __clang__
54 #pragma clang diagnostic ignored "-Wshorten-64-to-32"
55 #endif
56
57 namespace cv
58 {
59 namespace optflow
60 {
61 namespace
62 {
63
64 #ifndef M_SQRT2
65 const float M_SQRT2 = 1.41421356237309504880;
66 #endif
67
mathSign(T val)68 template <typename T> inline int mathSign( T val ) { return ( T( 0 ) < val ) - ( val < T( 0 ) ); }
69
70 /* Stable symmetric Householder reflection that gives c and s such that
71 * [ c s ][a] = [d],
72 * [ s -c ][b] [0]
73 *
74 * Output:
75 * c -- cosine(theta), where theta is the implicit angle of rotation
76 * (counter-clockwise) in a plane-rotation
77 * s -- sine(theta)
78 * r -- two-norm of [a; b]
79 */
symOrtho(double a,double b,double & c,double & s,double & r)80 inline void symOrtho( double a, double b, double &c, double &s, double &r )
81 {
82 if ( b == 0 )
83 {
84 c = mathSign( a );
85 s = 0;
86 r = std::abs( a );
87 }
88 else if ( a == 0 )
89 {
90 c = 0;
91 s = mathSign( b );
92 r = std::abs( b );
93 }
94 else if ( std::abs( b ) > std::abs( a ) )
95 {
96 const double tau = a / b;
97 s = mathSign( b ) / std::sqrt( 1 + tau * tau );
98 c = s * tau;
99 r = b / s;
100 }
101 else
102 {
103 const double tau = b / a;
104 c = mathSign( a ) / std::sqrt( 1 + tau * tau );
105 s = c * tau;
106 r = a / c;
107 }
108 }
109
110 /* Iterative LSQR algorithm for solving least squares problems.
111 *
112 * [1] Paige, C. C. and M. A. Saunders,
113 * LSQR: An Algorithm for Sparse Linear Equations And Sparse Least Squares
114 * ACM Trans. Math. Soft., Vol.8, 1982, pp. 43-71.
115 *
116 * Solves the following problem:
117 * argmin_x ||Ax - b|| + damp||x||
118 *
119 * Output:
120 * x -- approximate solution
121 */
solveLSQR(const Mat & A,const Mat & b,OutputArray xOut,const double damp=0.0,const unsigned iter_lim=10)122 void solveLSQR( const Mat &A, const Mat &b, OutputArray xOut, const double damp = 0.0, const unsigned iter_lim = 10 )
123 {
124 const int n = A.size().width;
125 CV_Assert( A.size().height == b.size().height );
126 CV_Assert( A.type() == CV_32F );
127 CV_Assert( b.type() == CV_32F );
128 xOut.create( n, 1, CV_32F );
129
130 Mat v( n, 1, CV_32F, 0.0f );
131 Mat u = b;
132 Mat x = xOut.getMat();
133 x = Mat::zeros( x.size(), x.type() );
134 double alfa = 0;
135 double beta = cv::norm( u, NORM_L2 );
136 Mat w( n, 1, CV_32F, 0.0f );
137 const Mat AT = A.t();
138
139 if ( beta > 0 )
140 {
141 u *= 1 / beta;
142 v = AT * u;
143 alfa = cv::norm( v, NORM_L2 );
144 }
145
146 if ( alfa > 0 )
147 {
148 v *= 1 / alfa;
149 w = v.clone();
150 }
151
152 double rhobar = alfa;
153 double phibar = beta;
154 if ( alfa * beta == 0 )
155 return;
156
157 for ( unsigned itn = 0; itn < iter_lim; ++itn )
158 {
159 u *= -alfa;
160 u += A * v;
161 beta = cv::norm( u, NORM_L2 );
162
163 if ( beta > 0 )
164 {
165 u *= 1 / beta;
166 v *= -beta;
167 v += AT * u;
168 alfa = cv::norm( v, NORM_L2 );
169 if ( alfa > 0 )
170 v *= 1 / alfa;
171 }
172
173 double rhobar1 = sqrt( rhobar * rhobar + damp * damp );
174 double cs1 = rhobar / rhobar1;
175 phibar = cs1 * phibar;
176
177 double cs, sn, rho;
178 symOrtho( rhobar1, beta, cs, sn, rho );
179
180 double theta = sn * alfa;
181 rhobar = -cs * alfa;
182 double phi = cs * phibar;
183 phibar = sn * phibar;
184
185 double t1 = phi / rho;
186 double t2 = -theta / rho;
187
188 x += t1 * w;
189 w *= t2;
190 w += v;
191 }
192 }
193
_cpu_fillDCTSampledPoints(float * row,const Point2f & p,const Size & basisSize,const Size & size)194 inline void _cpu_fillDCTSampledPoints( float *row, const Point2f &p, const Size &basisSize, const Size &size )
195 {
196 for ( int n1 = 0; n1 < basisSize.width; ++n1 )
197 for ( int n2 = 0; n2 < basisSize.height; ++n2 )
198 row[n1 * basisSize.height + n2] =
199 cosf( ( n1 * CV_PI / size.width ) * ( p.x + 0.5 ) ) * cosf( ( n2 * CV_PI / size.height ) * ( p.y + 0.5 ) );
200 }
201
202 ocl::ProgramSource _ocl_fillDCTSampledPointsSource(
203 "__kernel void fillDCTSampledPoints(__global const uchar* features, int fstep, int foff, __global "
204 "uchar* A, int Astep, int Aoff, int fs, int bsw, int bsh, int sw, int sh) {"
205 "const int i = get_global_id(0);"
206 "const int n1 = get_global_id(1);"
207 "const int n2 = get_global_id(2);"
208 "if (i >= fs || n1 >= bsw || n2 >= bsh) return;"
209 "__global const float2* f = (__global const float2*)(features + (fstep * i + foff));"
210 "__global float* a = (__global float*)(A + (Astep * i + Aoff + (n1 * bsh + n2) * sizeof(float)));"
211 "const float2 p = f[0];"
212 "const float pi = 3.14159265358979323846;"
213 "a[0] = cos((n1 * pi / sw) * (p.x + 0.5)) * cos((n2 * pi / sh) * (p.y + 0.5));"
214 "}" );
215
applyCLAHE(UMat & img,float claheClip)216 void applyCLAHE( UMat &img, float claheClip )
217 {
218 Ptr<CLAHE> clahe = createCLAHE();
219 clahe->setClipLimit( claheClip );
220 clahe->apply( img, img );
221 }
222
reduceToFlow(const Mat & w1,const Mat & w2,Mat & flow,const Size & basisSize)223 void reduceToFlow( const Mat &w1, const Mat &w2, Mat &flow, const Size &basisSize )
224 {
225 const Size size = flow.size();
226 Mat flowX( size, CV_32F, 0.0f );
227 Mat flowY( size, CV_32F, 0.0f );
228
229 const float mult = sqrt( static_cast<float>(size.area()) ) * 0.5;
230
231 for ( int i = 0; i < basisSize.width; ++i )
232 for ( int j = 0; j < basisSize.height; ++j )
233 {
234 flowX.at<float>( j, i ) = w1.at<float>( i * basisSize.height + j ) * mult;
235 flowY.at<float>( j, i ) = w2.at<float>( i * basisSize.height + j ) * mult;
236 }
237 for ( int i = 0; i < basisSize.height; ++i )
238 {
239 flowX.at<float>( i, 0 ) *= M_SQRT2;
240 flowY.at<float>( i, 0 ) *= M_SQRT2;
241 }
242 for ( int i = 0; i < basisSize.width; ++i )
243 {
244 flowX.at<float>( 0, i ) *= M_SQRT2;
245 flowY.at<float>( 0, i ) *= M_SQRT2;
246 }
247
248 dct( flowX, flowX, DCT_INVERSE );
249 dct( flowY, flowY, DCT_INVERSE );
250 for ( int i = 0; i < size.height; ++i )
251 for ( int j = 0; j < size.width; ++j )
252 flow.at<Point2f>( i, j ) = Point2f( flowX.at<float>( i, j ), flowY.at<float>( i, j ) );
253 }
254 }
255
findSparseFeatures(UMat & from,UMat & to,std::vector<Point2f> & features,std::vector<Point2f> & predictedFeatures) const256 void OpticalFlowPCAFlow::findSparseFeatures( UMat &from, UMat &to, std::vector<Point2f> &features,
257 std::vector<Point2f> &predictedFeatures ) const
258 {
259 Size size = from.size();
260 const unsigned maxFeatures = size.area() * sparseRate;
261 goodFeaturesToTrack( from, features, maxFeatures * retainedCornersFraction, 0.005, 3 );
262
263 // Add points along the grid if not enough features
264 if ( maxFeatures > features.size() )
265 {
266 const unsigned missingPoints = maxFeatures - features.size();
267 const unsigned blockSize = sqrt( (float)size.area() / missingPoints );
268 for ( int x = blockSize / 2; x < size.width; x += blockSize )
269 for ( int y = blockSize / 2; y < size.height; y += blockSize )
270 features.push_back( Point2f( x, y ) );
271 }
272 std::vector<uchar> predictedStatus;
273 std::vector<float> predictedError;
274 calcOpticalFlowPyrLK( from, to, features, predictedFeatures, predictedStatus, predictedError );
275
276 size_t j = 0;
277 for ( size_t i = 0; i < features.size(); ++i )
278 {
279 if ( predictedStatus[i] )
280 {
281 features[j] = features[i];
282 predictedFeatures[j] = predictedFeatures[i];
283 ++j;
284 }
285 }
286 features.resize( j );
287 predictedFeatures.resize( j );
288 }
289
removeOcclusions(UMat & from,UMat & to,std::vector<Point2f> & features,std::vector<Point2f> & predictedFeatures) const290 void OpticalFlowPCAFlow::removeOcclusions( UMat &from, UMat &to, std::vector<Point2f> &features,
291 std::vector<Point2f> &predictedFeatures ) const
292 {
293 std::vector<uchar> predictedStatus;
294 std::vector<float> predictedError;
295 std::vector<Point2f> backwardFeatures;
296 calcOpticalFlowPyrLK( to, from, predictedFeatures, backwardFeatures, predictedStatus, predictedError );
297
298 size_t j = 0;
299 const float threshold = occlusionsThreshold * sqrt( static_cast<float>(from.size().area()) );
300 for ( size_t i = 0; i < predictedFeatures.size(); ++i )
301 {
302 if ( predictedStatus[i] )
303 {
304 Point2f flowDiff = features[i] - backwardFeatures[i];
305 if ( flowDiff.dot( flowDiff ) <= threshold )
306 {
307 features[j] = features[i];
308 predictedFeatures[j] = predictedFeatures[i];
309 ++j;
310 }
311 }
312 }
313 features.resize( j );
314 predictedFeatures.resize( j );
315 }
316
getSystem(OutputArray AOut,OutputArray b1Out,OutputArray b2Out,const std::vector<Point2f> & features,const std::vector<Point2f> & predictedFeatures,const Size size)317 void OpticalFlowPCAFlow::getSystem( OutputArray AOut, OutputArray b1Out, OutputArray b2Out,
318 const std::vector<Point2f> &features, const std::vector<Point2f> &predictedFeatures,
319 const Size size )
320 {
321 AOut.create( features.size(), basisSize.area(), CV_32F );
322 b1Out.create( features.size(), 1, CV_32F );
323 b2Out.create( features.size(), 1, CV_32F );
324 if ( useOpenCL )
325 {
326 UMat A = AOut.getUMat();
327 Mat b1 = b1Out.getMat();
328 Mat b2 = b2Out.getMat();
329
330 ocl::Kernel kernel( "fillDCTSampledPoints", _ocl_fillDCTSampledPointsSource );
331 CV_Assert(basisSize.width > 0 && basisSize.height > 0);
332 size_t globSize[] = {features.size(), (size_t)basisSize.width, (size_t)basisSize.height};
333 kernel
334 .args( cv::ocl::KernelArg::ReadOnlyNoSize( Mat( features ).getUMat( ACCESS_READ ) ),
335 cv::ocl::KernelArg::WriteOnlyNoSize( A ), (int)features.size(), (int)basisSize.width,
336 (int)basisSize.height, (int)size.width, (int)size.height )
337 .run( 3, globSize, 0, true );
338
339 for ( size_t i = 0; i < features.size(); ++i )
340 {
341 const Point2f flow = predictedFeatures[i] - features[i];
342 b1.at<float>( i ) = flow.x;
343 b2.at<float>( i ) = flow.y;
344 }
345 }
346 else
347 {
348 Mat A = AOut.getMat();
349 Mat b1 = b1Out.getMat();
350 Mat b2 = b2Out.getMat();
351
352 for ( size_t i = 0; i < features.size(); ++i )
353 {
354 _cpu_fillDCTSampledPoints( A.ptr<float>( i ), features[i], basisSize, size );
355 const Point2f flow = predictedFeatures[i] - features[i];
356 b1.at<float>( i ) = flow.x;
357 b2.at<float>( i ) = flow.y;
358 }
359 }
360 }
361
getSystem(OutputArray A1Out,OutputArray A2Out,OutputArray b1Out,OutputArray b2Out,const std::vector<Point2f> & features,const std::vector<Point2f> & predictedFeatures,const Size size)362 void OpticalFlowPCAFlow::getSystem( OutputArray A1Out, OutputArray A2Out, OutputArray b1Out, OutputArray b2Out,
363 const std::vector<Point2f> &features, const std::vector<Point2f> &predictedFeatures,
364 const Size size )
365 {
366 CV_Assert( prior->getBasisSize() == basisSize.area() );
367
368 A1Out.create( features.size() + prior->getPadding(), basisSize.area(), CV_32F );
369 A2Out.create( features.size() + prior->getPadding(), basisSize.area(), CV_32F );
370 b1Out.create( features.size() + prior->getPadding(), 1, CV_32F );
371 b2Out.create( features.size() + prior->getPadding(), 1, CV_32F );
372
373 if ( useOpenCL )
374 {
375 UMat A = A1Out.getUMat();
376 Mat b1 = b1Out.getMat();
377 Mat b2 = b2Out.getMat();
378
379 ocl::Kernel kernel( "fillDCTSampledPoints", _ocl_fillDCTSampledPointsSource );
380 CV_Assert(basisSize.width > 0 && basisSize.height > 0);
381 size_t globSize[] = {features.size(), (size_t)basisSize.width, (size_t)basisSize.height};
382 kernel
383 .args( cv::ocl::KernelArg::ReadOnlyNoSize( Mat( features ).getUMat( ACCESS_READ ) ),
384 cv::ocl::KernelArg::WriteOnlyNoSize( A ), (int)features.size(), (int)basisSize.width,
385 (int)basisSize.height, (int)size.width, (int)size.height )
386 .run( 3, globSize, 0, true );
387
388 for ( size_t i = 0; i < features.size(); ++i )
389 {
390 const Point2f flow = predictedFeatures[i] - features[i];
391 b1.at<float>( i ) = flow.x;
392 b2.at<float>( i ) = flow.y;
393 }
394 }
395 else
396 {
397 Mat A1 = A1Out.getMat();
398 Mat b1 = b1Out.getMat();
399 Mat b2 = b2Out.getMat();
400
401 for ( size_t i = 0; i < features.size(); ++i )
402 {
403 _cpu_fillDCTSampledPoints( A1.ptr<float>( i ), features[i], basisSize, size );
404 const Point2f flow = predictedFeatures[i] - features[i];
405 b1.at<float>( i ) = flow.x;
406 b2.at<float>( i ) = flow.y;
407 }
408 }
409
410 Mat A1 = A1Out.getMat();
411 Mat A2 = A2Out.getMat();
412 Mat b1 = b1Out.getMat();
413 Mat b2 = b2Out.getMat();
414
415 memcpy( A2.ptr<float>(), A1.ptr<float>(), features.size() * basisSize.area() * sizeof( float ) );
416 prior->fillConstraints( A1.ptr<float>( features.size(), 0 ), A2.ptr<float>( features.size(), 0 ),
417 b1.ptr<float>( features.size(), 0 ), b2.ptr<float>( features.size(), 0 ) );
418 }
419
calc(InputArray I0,InputArray I1,InputOutputArray flowOut)420 void OpticalFlowPCAFlow::calc( InputArray I0, InputArray I1, InputOutputArray flowOut )
421 {
422 const Size size = I0.size();
423 CV_Assert( size == I1.size() );
424
425 UMat from, to;
426 if ( I0.channels() == 3 )
427 {
428 cvtColor( I0, from, COLOR_BGR2GRAY );
429 from.convertTo( from, CV_8U );
430 }
431 else
432 {
433 I0.getMat().convertTo( from, CV_8U );
434 }
435 if ( I1.channels() == 3 )
436 {
437 cvtColor( I1, to, COLOR_BGR2GRAY );
438 to.convertTo( to, CV_8U );
439 }
440 else
441 {
442 I1.getMat().convertTo( to, CV_8U );
443 }
444
445 CV_Assert( from.channels() == 1 );
446 CV_Assert( to.channels() == 1 );
447
448 const Mat fromOrig = from.getMat( ACCESS_READ ).clone();
449 useOpenCL = flowOut.isUMat() && ocl::useOpenCL();
450
451 applyCLAHE( from, claheClip );
452 applyCLAHE( to, claheClip );
453
454 std::vector<Point2f> features, predictedFeatures;
455 findSparseFeatures( from, to, features, predictedFeatures );
456 removeOcclusions( from, to, features, predictedFeatures );
457
458 flowOut.create( size, CV_32FC2 );
459 Mat flow = flowOut.getMat();
460
461 Mat w1, w2;
462 if ( prior.get() )
463 {
464 Mat A1, A2, b1, b2;
465 getSystem( A1, A2, b1, b2, features, predictedFeatures, size );
466 solveLSQR( A1, b1, w1, dampingFactor * size.area() );
467 solveLSQR( A2, b2, w2, dampingFactor * size.area() );
468 }
469 else
470 {
471 Mat A, b1, b2;
472 getSystem( A, b1, b2, features, predictedFeatures, size );
473 solveLSQR( A, b1, w1, dampingFactor * size.area() );
474 solveLSQR( A, b2, w2, dampingFactor * size.area() );
475 }
476 Mat flowSmall( ( size / 8 ) * 2, CV_32FC2 );
477 reduceToFlow( w1, w2, flowSmall, basisSize );
478 resize( flowSmall, flow, size, 0, 0, INTER_LINEAR );
479 ximgproc::fastGlobalSmootherFilter( fromOrig, flow, flow, 500, 2 );
480 }
481
OpticalFlowPCAFlow(Ptr<const PCAPrior> _prior,const Size _basisSize,float _sparseRate,float _retainedCornersFraction,float _occlusionsThreshold,float _dampingFactor,float _claheClip)482 OpticalFlowPCAFlow::OpticalFlowPCAFlow( Ptr<const PCAPrior> _prior, const Size _basisSize, float _sparseRate,
483 float _retainedCornersFraction, float _occlusionsThreshold,
484 float _dampingFactor, float _claheClip )
485 : prior( _prior ), basisSize( _basisSize ), sparseRate( _sparseRate ),
486 retainedCornersFraction( _retainedCornersFraction ), occlusionsThreshold( _occlusionsThreshold ),
487 dampingFactor( _dampingFactor ), claheClip( _claheClip ), useOpenCL( false )
488 {
489 CV_Assert( sparseRate > 0 && sparseRate <= 0.1 );
490 CV_Assert( retainedCornersFraction >= 0 && retainedCornersFraction <= 1.0 );
491 CV_Assert( occlusionsThreshold > 0 );
492 }
493
collectGarbage()494 void OpticalFlowPCAFlow::collectGarbage() {}
495
createOptFlow_PCAFlow()496 Ptr<DenseOpticalFlow> createOptFlow_PCAFlow() { return makePtr<OpticalFlowPCAFlow>(); }
497
PCAPrior(const char * pathToPrior)498 PCAPrior::PCAPrior( const char *pathToPrior )
499 {
500 FILE *f = fopen( pathToPrior, "rb" );
501 CV_Assert( f );
502
503 unsigned n = 0, m = 0;
504 CV_Assert( fread( &n, sizeof( n ), 1, f ) == 1 );
505 CV_Assert( fread( &m, sizeof( m ), 1, f ) == 1 );
506
507 L1.create( n, m, CV_32F );
508 L2.create( n, m, CV_32F );
509 c1.create( n, 1, CV_32F );
510 c2.create( n, 1, CV_32F );
511
512 CV_Assert( fread( L1.ptr<float>(), n * m * sizeof( float ), 1, f ) == 1 );
513 CV_Assert( fread( L2.ptr<float>(), n * m * sizeof( float ), 1, f ) == 1 );
514 CV_Assert( fread( c1.ptr<float>(), n * sizeof( float ), 1, f ) == 1 );
515 CV_Assert( fread( c2.ptr<float>(), n * sizeof( float ), 1, f ) == 1 );
516
517 fclose( f );
518 }
519
fillConstraints(float * A1,float * A2,float * b1,float * b2) const520 void PCAPrior::fillConstraints( float *A1, float *A2, float *b1, float *b2 ) const
521 {
522 memcpy( A1, L1.ptr<float>(), L1.size().area() * sizeof( float ) );
523 memcpy( A2, L2.ptr<float>(), L2.size().area() * sizeof( float ) );
524 memcpy( b1, c1.ptr<float>(), c1.size().area() * sizeof( float ) );
525 memcpy( b2, c2.ptr<float>(), c2.size().area() * sizeof( float ) );
526 }
527 }
528 }
529