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 //                        Intel License Agreement
11 //                For Open Source Computer Vision Library
12 //
13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
14 // Third party copyrights are property of their respective owners.
15 //
16 // Redistribution and use in source and binary forms, with or without modification,
17 // are permitted provided that the following conditions are met:
18 //
19 //   * Redistribution's of source code must retain the above copyright notice,
20 //     this list of conditions and the following disclaimer.
21 //
22 //   * Redistribution's in binary form must reproduce the above copyright notice,
23 //     this list of conditions and the following disclaimer in the documentation
24 //     and/or other materials provided with the distribution.
25 //
26 //   * The name of Intel Corporation may not be used to endorse or promote products
27 //     derived from this software without specific prior written permission.
28 //
29 // This software is provided by the copyright holders and contributors "as is" and
30 // any express or implied warranties, including, but not limited to, the implied
31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
32 // In no event shall the Intel Corporation or contributors be liable for any direct,
33 // indirect, incidental, special, exemplary, or consequential damages
34 // (including, but not limited to, procurement of substitute goods or services;
35 // loss of use, data, or profits; or business interruption) however caused
36 // and on any theory of liability, whether in contract, strict liability,
37 // or tort (including negligence or otherwise) arising in any way out of
38 // the use of this software, even if advised of the possibility of such damage.
39 //
40 //M*/
41 
42 #include "test_precomp.hpp"
43 #include "opencv2/core/core_c.h"
44 #include "opencv2/calib3d/calib3d_c.h"
45 
46 namespace cvtest {
47 
cvTsRodrigues(const CvMat * src,CvMat * dst,CvMat * jacobian)48 static int cvTsRodrigues( const CvMat* src, CvMat* dst, CvMat* jacobian )
49 {
50     int depth;
51     int i;
52     float Jf[27];
53     double J[27];
54     CvMat _Jf, matJ = cvMat( 3, 9, CV_64F, J );
55 
56     depth = CV_MAT_DEPTH(src->type);
57 
58     if( jacobian )
59     {
60         assert( (jacobian->rows == 9 && jacobian->cols == 3) ||
61                 (jacobian->rows == 3 && jacobian->cols == 9) );
62     }
63 
64     if( src->cols == 1 || src->rows == 1 )
65     {
66         double r[3], theta;
67         CvMat _r = cvMat( src->rows, src->cols, CV_MAKETYPE(CV_64F,CV_MAT_CN(src->type)), r);
68 
69         assert( dst->rows == 3 && dst->cols == 3 );
70 
71         cvConvert( src, &_r );
72 
73         theta = sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
74         if( theta < DBL_EPSILON )
75         {
76             cvSetIdentity( dst );
77 
78             if( jacobian )
79             {
80                 memset( J, 0, sizeof(J) );
81                 J[5] = J[15] = J[19] = 1;
82                 J[7] = J[11] = J[21] = -1;
83             }
84         }
85         else
86         {
87             // omega = r/theta (~[w1, w2, w3])
88             double itheta = 1./theta;
89             double w1 = r[0]*itheta, w2 = r[1]*itheta, w3 = r[2]*itheta;
90             double alpha = cos(theta);
91             double beta = sin(theta);
92             double gamma = 1 - alpha;
93             double omegav[] =
94             {
95                 0, -w3, w2,
96                 w3, 0, -w1,
97                 -w2, w1, 0
98             };
99             double A[] =
100             {
101                 w1*w1, w1*w2, w1*w3,
102                 w2*w1, w2*w2, w2*w3,
103                 w3*w1, w3*w2, w3*w3
104             };
105             double R[9];
106             CvMat _omegav = cvMat(3, 3, CV_64F, omegav);
107             CvMat matA = cvMat(3, 3, CV_64F, A);
108             CvMat matR = cvMat(3, 3, CV_64F, R);
109 
110             cvSetIdentity( &matR, cvRealScalar(alpha) );
111             cvScaleAdd( &_omegav, cvRealScalar(beta), &matR, &matR );
112             cvScaleAdd( &matA, cvRealScalar(gamma), &matR, &matR );
113             cvConvert( &matR, dst );
114 
115             if( jacobian )
116             {
117                 // m3 = [r, theta]
118                 double dm3din[] =
119                 {
120                     1, 0, 0,
121                     0, 1, 0,
122                     0, 0, 1,
123                     w1, w2, w3
124                 };
125                 // m2 = [omega, theta]
126                 double dm2dm3[] =
127                 {
128                     itheta, 0, 0, -w1*itheta,
129                     0, itheta, 0, -w2*itheta,
130                     0, 0, itheta, -w3*itheta,
131                     0, 0, 0, 1
132                 };
133                 double t0[9*4];
134                 double dm1dm2[21*4];
135                 double dRdm1[9*21];
136                 CvMat _dm3din = cvMat( 4, 3, CV_64FC1, dm3din );
137                 CvMat _dm2dm3 = cvMat( 4, 4, CV_64FC1, dm2dm3 );
138                 CvMat _dm1dm2 = cvMat( 21, 4, CV_64FC1, dm1dm2 );
139                 CvMat _dRdm1 = cvMat( 9, 21, CV_64FC1, dRdm1 );
140                 CvMat _dRdm1_part;
141                 CvMat _t0 = cvMat( 9, 4, CV_64FC1, t0 );
142                 CvMat _t1 = cvMat( 9, 4, CV_64FC1, dRdm1 );
143 
144                 // m1 = [alpha, beta, gamma, omegav; A]
145                 memset( dm1dm2, 0, sizeof(dm1dm2) );
146                 dm1dm2[3] = -beta;
147                 dm1dm2[7] = alpha;
148                 dm1dm2[11] = beta;
149 
150                 // dm1dm2(4:12,1:3) = [0 0 0 0 0 1 0 -1 0;
151                 //                     0 0 -1 0 0 0 1 0 0;
152                 //                     0 1 0 -1 0 0 0 0 0]'
153                 //                     -------------------
154                 //                     0 0 0  0 0 0 0 0 0
155                 dm1dm2[12 + 6] = dm1dm2[12 + 20] = dm1dm2[12 + 25] = 1;
156                 dm1dm2[12 + 9] = dm1dm2[12 + 14] = dm1dm2[12 + 28] = -1;
157 
158                 double dm1dw[] =
159                 {
160                     2*w1, w2, w3, w2, 0, 0, w3, 0, 0,
161                     0, w1, 0, w1, 2*w2, w3, 0, w3, 0,
162                     0, 0, w1, 0, 0, w2, w1, w2, 2*w3
163                 };
164 
165                 CvMat _dm1dw = cvMat( 3, 9, CV_64FC1, dm1dw );
166                 CvMat _dm1dm2_part;
167 
168                 cvGetSubRect( &_dm1dm2, &_dm1dm2_part, cvRect(0,12,3,9) );
169                 cvTranspose( &_dm1dw, &_dm1dm2_part );
170 
171                 memset( dRdm1, 0, sizeof(dRdm1) );
172                 dRdm1[0*21] = dRdm1[4*21] = dRdm1[8*21] = 1;
173 
174                 cvGetCol( &_dRdm1, &_dRdm1_part, 1 );
175                 cvTranspose( &_omegav, &_omegav );
176                 cvReshape( &_omegav, &_omegav, 1, 1 );
177                 cvTranspose( &_omegav, &_dRdm1_part );
178 
179                 cvGetCol( &_dRdm1, &_dRdm1_part, 2 );
180                 cvReshape( &matA, &matA, 1, 1 );
181                 cvTranspose( &matA, &_dRdm1_part );
182 
183                 cvGetSubRect( &_dRdm1, &_dRdm1_part, cvRect(3,0,9,9) );
184                 cvSetIdentity( &_dRdm1_part, cvScalarAll(beta) );
185 
186                 cvGetSubRect( &_dRdm1, &_dRdm1_part, cvRect(12,0,9,9) );
187                 cvSetIdentity( &_dRdm1_part, cvScalarAll(gamma) );
188 
189                 matJ = cvMat( 9, 3, CV_64FC1, J );
190 
191                 cvMatMul( &_dRdm1, &_dm1dm2, &_t0 );
192                 cvMatMul( &_t0, &_dm2dm3, &_t1 );
193                 cvMatMul( &_t1, &_dm3din, &matJ );
194 
195                 _t0 = cvMat( 3, 9, CV_64FC1, t0 );
196                 cvTranspose( &matJ, &_t0 );
197 
198                 for( i = 0; i < 3; i++ )
199                 {
200                     _t1 = cvMat( 3, 3, CV_64FC1, t0 + i*9 );
201                     cvTranspose( &_t1, &_t1 );
202                 }
203 
204                 cvTranspose( &_t0, &matJ );
205             }
206         }
207     }
208     else if( src->cols == 3 && src->rows == 3 )
209     {
210         double R[9], A[9], I[9], r[3], W[3], U[9], V[9];
211         double tr, alpha, beta, theta;
212         CvMat matR = cvMat( 3, 3, CV_64F, R );
213         CvMat matA = cvMat( 3, 3, CV_64F, A );
214         CvMat matI = cvMat( 3, 3, CV_64F, I );
215         CvMat _r = cvMat( dst->rows, dst->cols, CV_MAKETYPE(CV_64F, CV_MAT_CN(dst->type)), r );
216         CvMat matW = cvMat( 1, 3, CV_64F, W );
217         CvMat matU = cvMat( 3, 3, CV_64F, U );
218         CvMat matV = cvMat( 3, 3, CV_64F, V );
219 
220         cvConvert( src, &matR );
221         cvSVD( &matR, &matW, &matU, &matV, CV_SVD_MODIFY_A + CV_SVD_U_T + CV_SVD_V_T );
222         cvGEMM( &matU, &matV, 1, 0, 0, &matR, CV_GEMM_A_T );
223 
224         cvMulTransposed( &matR, &matA, 0 );
225         cvSetIdentity( &matI );
226 
227         if( cvNorm( &matA, &matI, CV_C ) > 1e-3 ||
228             fabs( cvDet(&matR) - 1 ) > 1e-3 )
229             return 0;
230 
231         tr = (cvTrace(&matR).val[0] - 1.)*0.5;
232         tr = tr > 1. ? 1. : tr < -1. ? -1. : tr;
233         theta = acos(tr);
234         alpha = cos(theta);
235         beta = sin(theta);
236 
237         if( beta >= 1e-5 )
238         {
239             double dtheta_dtr = -1./sqrt(1 - tr*tr);
240             double vth = 1/(2*beta);
241 
242             // om1 = [R(3,2) - R(2,3), R(1,3) - R(3,1), R(2,1) - R(1,2)]'
243             double om1[] = { R[7] - R[5], R[2] - R[6], R[3] - R[1] };
244             // om = om1*vth
245             // r = om*theta
246             double d3 = vth*theta;
247 
248             r[0] = om1[0]*d3; r[1] = om1[1]*d3; r[2] = om1[2]*d3;
249             cvConvert( &_r, dst );
250 
251             if( jacobian )
252             {
253                 // var1 = [vth;theta]
254                 // var = [om1;var1] = [om1;vth;theta]
255                 double dvth_dtheta = -vth*alpha/beta;
256                 double d1 = 0.5*dvth_dtheta*dtheta_dtr;
257                 double d2 = 0.5*dtheta_dtr;
258                 // dvar1/dR = dvar1/dtheta*dtheta/dR = [dvth/dtheta; 1] * dtheta/dtr * dtr/dR
259                 double dvardR[5*9] =
260                 {
261                     0, 0, 0, 0, 0, 1, 0, -1, 0,
262                     0, 0, -1, 0, 0, 0, 1, 0, 0,
263                     0, 1, 0, -1, 0, 0, 0, 0, 0,
264                     d1, 0, 0, 0, d1, 0, 0, 0, d1,
265                     d2, 0, 0, 0, d2, 0, 0, 0, d2
266                 };
267                 // var2 = [om;theta]
268                 double dvar2dvar[] =
269                 {
270                     vth, 0, 0, om1[0], 0,
271                     0, vth, 0, om1[1], 0,
272                     0, 0, vth, om1[2], 0,
273                     0, 0, 0, 0, 1
274                 };
275                 double domegadvar2[] =
276                 {
277                     theta, 0, 0, om1[0]*vth,
278                     0, theta, 0, om1[1]*vth,
279                     0, 0, theta, om1[2]*vth
280                 };
281 
282                 CvMat _dvardR = cvMat( 5, 9, CV_64FC1, dvardR );
283                 CvMat _dvar2dvar = cvMat( 4, 5, CV_64FC1, dvar2dvar );
284                 CvMat _domegadvar2 = cvMat( 3, 4, CV_64FC1, domegadvar2 );
285                 double t0[3*5];
286                 CvMat _t0 = cvMat( 3, 5, CV_64FC1, t0 );
287 
288                 cvMatMul( &_domegadvar2, &_dvar2dvar, &_t0 );
289                 cvMatMul( &_t0, &_dvardR, &matJ );
290             }
291         }
292         else if( tr > 0 )
293         {
294             cvZero( dst );
295             if( jacobian )
296             {
297                 memset( J, 0, sizeof(J) );
298                 J[5] = J[15] = J[19] = 0.5;
299                 J[7] = J[11] = J[21] = -0.5;
300             }
301         }
302         else
303         {
304             r[0] = theta*sqrt((R[0] + 1)*0.5);
305             r[1] = theta*sqrt((R[4] + 1)*0.5)*(R[1] >= 0 ? 1 : -1);
306             r[2] = theta*sqrt((R[8] + 1)*0.5)*(R[2] >= 0 ? 1 : -1);
307             cvConvert( &_r, dst );
308 
309             if( jacobian )
310                 memset( J, 0, sizeof(J) );
311         }
312 
313         if( jacobian )
314         {
315             for( i = 0; i < 3; i++ )
316             {
317                 CvMat t = cvMat( 3, 3, CV_64F, J + i*9 );
318                 cvTranspose( &t, &t );
319             }
320         }
321     }
322     else
323     {
324         assert(0);
325         return 0;
326     }
327 
328     if( jacobian )
329     {
330         if( depth == CV_32F )
331         {
332             if( jacobian->rows == matJ.rows )
333                 cvConvert( &matJ, jacobian );
334             else
335             {
336                 _Jf = cvMat( matJ.rows, matJ.cols, CV_32FC1, Jf );
337                 cvConvert( &matJ, &_Jf );
338                 cvTranspose( &_Jf, jacobian );
339             }
340         }
341         else if( jacobian->rows == matJ.rows )
342             cvCopy( &matJ, jacobian );
343         else
344             cvTranspose( &matJ, jacobian );
345     }
346 
347     return 1;
348 }
349 
350 
Rodrigues(const Mat & src,Mat & dst,Mat * jac)351 /*extern*/ void Rodrigues(const Mat& src, Mat& dst, Mat* jac)
352 {
353     if(src.rows == 1 || src.cols == 1)
354         dst.create(3, 3, src.depth());
355     else
356         dst.create(3, 1, src.depth());
357     CvMat _src = cvMat(src), _dst = cvMat(dst), _jac;
358     if( jac )
359         _jac = cvMat(*jac);
360     cvTsRodrigues(&_src, &_dst, jac ? &_jac : 0);
361 }
362 
363 } // namespace
364 namespace opencv_test {
365 
test_convertHomogeneous(const Mat & _src,Mat & _dst)366 static void test_convertHomogeneous( const Mat& _src, Mat& _dst )
367 {
368     Mat src = _src, dst = _dst;
369     int i, count, sdims, ddims;
370     int sstep1, sstep2, dstep1, dstep2;
371 
372     if( src.depth() != CV_64F )
373         _src.convertTo(src, CV_64F);
374 
375     if( dst.depth() != CV_64F )
376         dst.create(dst.size(), CV_MAKETYPE(CV_64F, _dst.channels()));
377 
378     if( src.rows > src.cols )
379     {
380         count = src.rows;
381         sdims = src.channels()*src.cols;
382         sstep1 = (int)(src.step/sizeof(double));
383         sstep2 = 1;
384     }
385     else
386     {
387         count = src.cols;
388         sdims = src.channels()*src.rows;
389         if( src.rows == 1 )
390         {
391             sstep1 = sdims;
392             sstep2 = 1;
393         }
394         else
395         {
396             sstep1 = 1;
397             sstep2 = (int)(src.step/sizeof(double));
398         }
399     }
400 
401     if( dst.rows > dst.cols )
402     {
403         CV_Assert( count == dst.rows );
404         ddims = dst.channels()*dst.cols;
405         dstep1 = (int)(dst.step/sizeof(double));
406         dstep2 = 1;
407     }
408     else
409     {
410         assert( count == dst.cols );
411         ddims = dst.channels()*dst.rows;
412         if( dst.rows == 1 )
413         {
414             dstep1 = ddims;
415             dstep2 = 1;
416         }
417         else
418         {
419             dstep1 = 1;
420             dstep2 = (int)(dst.step/sizeof(double));
421         }
422     }
423 
424     double* s = src.ptr<double>();
425     double* d = dst.ptr<double>();
426 
427     if( sdims <= ddims )
428     {
429         int wstep = dstep2*(ddims - 1);
430 
431         for( i = 0; i < count; i++, s += sstep1, d += dstep1 )
432         {
433             double x = s[0];
434             double y = s[sstep2];
435 
436             d[wstep] = 1;
437             d[0] = x;
438             d[dstep2] = y;
439 
440             if( sdims >= 3 )
441             {
442                 d[dstep2*2] = s[sstep2*2];
443                 if( sdims == 4 )
444                     d[dstep2*3] = s[sstep2*3];
445             }
446         }
447     }
448     else
449     {
450         int wstep = sstep2*(sdims - 1);
451 
452         for( i = 0; i < count; i++, s += sstep1, d += dstep1 )
453         {
454             double w = s[wstep];
455             double x = s[0];
456             double y = s[sstep2];
457 
458             w = w ? 1./w : 1;
459 
460             d[0] = x*w;
461             d[dstep2] = y*w;
462 
463             if( ddims == 3 )
464                 d[dstep2*2] = s[sstep2*2]*w;
465         }
466     }
467 
468     if( dst.data != _dst.data )
469         dst.convertTo(_dst, _dst.depth());
470 }
471 
472 namespace {
473 
474 void
test_projectPoints(const Mat & _3d,const Mat & Rt,const Mat & A,Mat & _2d,RNG * rng,double sigma)475 test_projectPoints( const Mat& _3d, const Mat& Rt, const Mat& A, Mat& _2d, RNG* rng, double sigma )
476 {
477     CV_Assert( _3d.isContinuous() );
478 
479     double p[12];
480     Mat P( 3, 4, CV_64F, p );
481     gemm(A, Rt, 1, Mat(), 0, P);
482 
483     int i, count = _3d.cols;
484 
485     Mat noise;
486     if( rng )
487     {
488         if( sigma == 0 )
489             rng = 0;
490         else
491         {
492             noise.create( 1, _3d.cols, CV_64FC2 );
493             rng->fill(noise, RNG::NORMAL, Scalar::all(0), Scalar::all(sigma) );
494         }
495     }
496 
497     Mat temp( 1, count, CV_64FC3 );
498 
499     for( i = 0; i < count; i++ )
500     {
501         const double* M = _3d.ptr<double>() + i*3;
502         double* m = temp.ptr<double>() + i*3;
503         double X = M[0], Y = M[1], Z = M[2];
504         double u = p[0]*X + p[1]*Y + p[2]*Z + p[3];
505         double v = p[4]*X + p[5]*Y + p[6]*Z + p[7];
506         double s = p[8]*X + p[9]*Y + p[10]*Z + p[11];
507 
508         if( !noise.empty() )
509         {
510             u += noise.at<Point2d>(i).x*s;
511             v += noise.at<Point2d>(i).y*s;
512         }
513 
514         m[0] = u;
515         m[1] = v;
516         m[2] = s;
517     }
518 
519     test_convertHomogeneous( temp, _2d );
520 }
521 
522 
523 /********************************** Rodrigues transform ********************************/
524 
525 class CV_RodriguesTest : public cvtest::ArrayTest
526 {
527 public:
528     CV_RodriguesTest();
529 
530 protected:
531     int read_params( const cv::FileStorage& fs );
532     void fill_array( int test_case_idx, int i, int j, Mat& arr );
533     int prepare_test_case( int test_case_idx );
534     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
535     double get_success_error_level( int test_case_idx, int i, int j );
536     void run_func();
537     void prepare_to_validation( int );
538 
539     bool calc_jacobians;
540     bool test_cpp;
541 };
542 
543 
CV_RodriguesTest()544 CV_RodriguesTest::CV_RodriguesTest()
545 {
546     test_array[INPUT].push_back(NULL);  // rotation vector
547     test_array[OUTPUT].push_back(NULL); // rotation matrix
548     test_array[OUTPUT].push_back(NULL); // jacobian (J)
549     test_array[OUTPUT].push_back(NULL); // rotation vector (backward transform result)
550     test_array[OUTPUT].push_back(NULL); // inverse transform jacobian (J1)
551     test_array[OUTPUT].push_back(NULL); // J*J1 (or J1*J) == I(3x3)
552     test_array[REF_OUTPUT].push_back(NULL);
553     test_array[REF_OUTPUT].push_back(NULL);
554     test_array[REF_OUTPUT].push_back(NULL);
555     test_array[REF_OUTPUT].push_back(NULL);
556     test_array[REF_OUTPUT].push_back(NULL);
557 
558     element_wise_relative_error = false;
559     calc_jacobians = false;
560 
561     test_cpp = false;
562 }
563 
564 
read_params(const cv::FileStorage & fs)565 int CV_RodriguesTest::read_params( const cv::FileStorage& fs )
566 {
567     int code = cvtest::ArrayTest::read_params( fs );
568     return code;
569 }
570 
571 
get_test_array_types_and_sizes(int,vector<vector<Size>> & sizes,vector<vector<int>> & types)572 void CV_RodriguesTest::get_test_array_types_and_sizes(
573     int /*test_case_idx*/, vector<vector<Size> >& sizes, vector<vector<int> >& types )
574 {
575     RNG& rng = ts->get_rng();
576     int depth = cvtest::randInt(rng) % 2 == 0 ? CV_32F : CV_64F;
577     int i, code;
578 
579     code = cvtest::randInt(rng) % 3;
580     types[INPUT][0] = CV_MAKETYPE(depth, 1);
581 
582     if( code == 0 )
583     {
584         sizes[INPUT][0] = cvSize(1,1);
585         types[INPUT][0] = CV_MAKETYPE(depth, 3);
586     }
587     else if( code == 1 )
588         sizes[INPUT][0] = cvSize(3,1);
589     else
590         sizes[INPUT][0] = cvSize(1,3);
591 
592     sizes[OUTPUT][0] = cvSize(3, 3);
593     types[OUTPUT][0] = CV_MAKETYPE(depth, 1);
594 
595     types[OUTPUT][1] = CV_MAKETYPE(depth, 1);
596 
597     if( cvtest::randInt(rng) % 2 )
598         sizes[OUTPUT][1] = cvSize(3,9);
599     else
600         sizes[OUTPUT][1] = cvSize(9,3);
601 
602     types[OUTPUT][2] = types[INPUT][0];
603     sizes[OUTPUT][2] = sizes[INPUT][0];
604 
605     types[OUTPUT][3] = types[OUTPUT][1];
606     sizes[OUTPUT][3] = cvSize(sizes[OUTPUT][1].height, sizes[OUTPUT][1].width);
607 
608     types[OUTPUT][4] = types[OUTPUT][1];
609     sizes[OUTPUT][4] = cvSize(3,3);
610 
611     calc_jacobians = cvtest::randInt(rng) % 3 != 0;
612     if( !calc_jacobians )
613         sizes[OUTPUT][1] = sizes[OUTPUT][3] = sizes[OUTPUT][4] = cvSize(0,0);
614 
615     for( i = 0; i < 5; i++ )
616     {
617         types[REF_OUTPUT][i] = types[OUTPUT][i];
618         sizes[REF_OUTPUT][i] = sizes[OUTPUT][i];
619     }
620     test_cpp = (cvtest::randInt(rng) & 256) == 0;
621 }
622 
623 
get_success_error_level(int,int,int j)624 double CV_RodriguesTest::get_success_error_level( int /*test_case_idx*/, int /*i*/, int j )
625 {
626     return j == 4 ? 1e-2 : 1e-2;
627 }
628 
629 
fill_array(int test_case_idx,int i,int j,Mat & arr)630 void CV_RodriguesTest::fill_array( int test_case_idx, int i, int j, Mat& arr )
631 {
632     if( i == INPUT && j == 0 )
633     {
634         double r[3], theta0, theta1, f;
635         Mat _r( arr.rows, arr.cols, CV_MAKETYPE(CV_64F,arr.channels()), r );
636         RNG& rng = ts->get_rng();
637 
638         r[0] = cvtest::randReal(rng)*CV_PI*2;
639         r[1] = cvtest::randReal(rng)*CV_PI*2;
640         r[2] = cvtest::randReal(rng)*CV_PI*2;
641 
642         theta0 = sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
643         theta1 = fmod(theta0, CV_PI*2);
644 
645         if( theta1 > CV_PI )
646             theta1 = -(CV_PI*2 - theta1);
647 
648         f = theta1/(theta0 ? theta0 : 1);
649         r[0] *= f;
650         r[1] *= f;
651         r[2] *= f;
652 
653         cvtest::convert( _r, arr, arr.type() );
654     }
655     else
656         cvtest::ArrayTest::fill_array( test_case_idx, i, j, arr );
657 }
658 
659 
prepare_test_case(int test_case_idx)660 int CV_RodriguesTest::prepare_test_case( int test_case_idx )
661 {
662     int code = cvtest::ArrayTest::prepare_test_case( test_case_idx );
663     return code;
664 }
665 
666 
run_func()667 void CV_RodriguesTest::run_func()
668 {
669     cv::Mat v = test_mat[INPUT][0], M = test_mat[OUTPUT][0], v2 = test_mat[OUTPUT][2];
670     cv::Mat M0 = M, v2_0 = v2;
671     if( !calc_jacobians )
672     {
673         cv::Rodrigues(v, M);
674         cv::Rodrigues(M, v2);
675     }
676     else
677     {
678         cv::Mat J1 = test_mat[OUTPUT][1], J2 = test_mat[OUTPUT][3];
679         cv::Mat J1_0 = J1, J2_0 = J2;
680         cv::Rodrigues(v, M, J1);
681         cv::Rodrigues(M, v2, J2);
682         if( J1.data != J1_0.data )
683         {
684             if( J1.size() != J1_0.size() )
685                 J1 = J1.t();
686             J1.convertTo(J1_0, J1_0.type());
687         }
688         if( J2.data != J2_0.data )
689         {
690             if( J2.size() != J2_0.size() )
691                 J2 = J2.t();
692             J2.convertTo(J2_0, J2_0.type());
693         }
694     }
695     if( M.data != M0.data )
696         M.reshape(M0.channels(), M0.rows).convertTo(M0, M0.type());
697     if( v2.data != v2_0.data )
698         v2.reshape(v2_0.channels(), v2_0.rows).convertTo(v2_0, v2_0.type());
699 }
700 
701 
prepare_to_validation(int)702 void CV_RodriguesTest::prepare_to_validation( int /*test_case_idx*/ )
703 {
704     const Mat& vec = test_mat[INPUT][0];
705     Mat& m = test_mat[REF_OUTPUT][0];
706     Mat& vec2 = test_mat[REF_OUTPUT][2];
707     Mat* v2m_jac = 0, *m2v_jac = 0;
708     double theta0, theta1;
709 
710     if( calc_jacobians )
711     {
712         v2m_jac = &test_mat[REF_OUTPUT][1];
713         m2v_jac = &test_mat[REF_OUTPUT][3];
714     }
715 
716 
717     cvtest::Rodrigues( vec, m, v2m_jac );
718     cvtest::Rodrigues( m, vec2, m2v_jac );
719     cvtest::copy( vec, vec2 );
720 
721     theta0 = cvtest::norm( vec2, CV_L2 );
722     theta1 = fmod( theta0, CV_PI*2 );
723 
724     if( theta1 > CV_PI )
725         theta1 = -(CV_PI*2 - theta1);
726     vec2 *= theta1/(theta0 ? theta0 : 1);
727 
728     if( calc_jacobians )
729     {
730         //cvInvert( v2m_jac, m2v_jac, CV_SVD );
731         double nrm = cvtest::norm(test_mat[REF_OUTPUT][3], CV_C);
732         if( FLT_EPSILON < nrm && nrm < 1000 )
733         {
734             gemm( test_mat[OUTPUT][1], test_mat[OUTPUT][3],
735                   1, Mat(), 0, test_mat[OUTPUT][4],
736                   v2m_jac->rows == 3 ? 0 : CV_GEMM_A_T + CV_GEMM_B_T );
737         }
738         else
739         {
740             setIdentity(test_mat[OUTPUT][4], Scalar::all(1.));
741             cvtest::copy( test_mat[REF_OUTPUT][2], test_mat[OUTPUT][2] );
742         }
743         setIdentity(test_mat[REF_OUTPUT][4], Scalar::all(1.));
744     }
745 }
746 
747 
748 /********************************** fundamental matrix *********************************/
749 
750 class CV_FundamentalMatTest : public cvtest::ArrayTest
751 {
752 public:
753     CV_FundamentalMatTest();
754 
755 protected:
756     int read_params( const cv::FileStorage& fs );
757     void fill_array( int test_case_idx, int i, int j, Mat& arr );
758     int prepare_test_case( int test_case_idx );
759     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
760     double get_success_error_level( int test_case_idx, int i, int j );
761     void run_func();
762     void prepare_to_validation( int );
763 
764     int method;
765     int img_size;
766     int cube_size;
767     int dims;
768     int f_result;
769     double min_f, max_f;
770     double sigma;
771     bool test_cpp;
772 };
773 
774 
CV_FundamentalMatTest()775 CV_FundamentalMatTest::CV_FundamentalMatTest()
776 {
777     // input arrays:
778     //   0, 1 - arrays of 2d points that are passed to %func%.
779     //          Can have different data type, layout, be stored in homogeneous coordinates or not.
780     //   2 - array of 3d points that are projected to both view planes
781     //   3 - [R|t] matrix for the second view plane (for the first one it is [I|0]
782     //   4, 5 - intrinsic matrices
783     test_array[INPUT].push_back(NULL);
784     test_array[INPUT].push_back(NULL);
785     test_array[INPUT].push_back(NULL);
786     test_array[INPUT].push_back(NULL);
787     test_array[INPUT].push_back(NULL);
788     test_array[INPUT].push_back(NULL);
789     test_array[TEMP].push_back(NULL);
790     test_array[TEMP].push_back(NULL);
791     test_array[OUTPUT].push_back(NULL);
792     test_array[OUTPUT].push_back(NULL);
793     test_array[REF_OUTPUT].push_back(NULL);
794     test_array[REF_OUTPUT].push_back(NULL);
795 
796     element_wise_relative_error = false;
797 
798     method = 0;
799     img_size = 10;
800     cube_size = 10;
801     dims = 0;
802     min_f = 1;
803     max_f = 3;
804     sigma = 0;//0.1;
805     f_result = 0;
806 
807     test_cpp = false;
808 }
809 
810 
read_params(const cv::FileStorage & fs)811 int CV_FundamentalMatTest::read_params( const cv::FileStorage& fs )
812 {
813     int code = cvtest::ArrayTest::read_params( fs );
814     return code;
815 }
816 
817 
get_test_array_types_and_sizes(int,vector<vector<Size>> & sizes,vector<vector<int>> & types)818 void CV_FundamentalMatTest::get_test_array_types_and_sizes( int /*test_case_idx*/,
819                                                 vector<vector<Size> >& sizes, vector<vector<int> >& types )
820 {
821     RNG& rng = ts->get_rng();
822     int pt_depth = cvtest::randInt(rng) % 2 == 0 ? CV_32F : CV_64F;
823     double pt_count_exp = cvtest::randReal(rng)*6 + 1;
824     int pt_count = cvRound(exp(pt_count_exp));
825 
826     dims = cvtest::randInt(rng) % 2 + 2;
827     method = 1 << (cvtest::randInt(rng) % 4);
828 
829     if( method == CV_FM_7POINT )
830         pt_count = 7;
831     else
832     {
833         pt_count = MAX( pt_count, 8 + (method == CV_FM_8POINT) );
834         if( pt_count >= 8 && cvtest::randInt(rng) % 2 )
835             method |= CV_FM_8POINT;
836     }
837 
838     types[INPUT][0] = CV_MAKETYPE(pt_depth, 1);
839 
840     sizes[INPUT][0] = cvSize(dims, pt_count);
841     if( cvtest::randInt(rng) % 2 )
842     {
843         types[INPUT][0] = CV_MAKETYPE(pt_depth, dims);
844         if( cvtest::randInt(rng) % 2 )
845             sizes[INPUT][0] = cvSize(pt_count, 1);
846         else
847             sizes[INPUT][0] = cvSize(1, pt_count);
848     }
849 
850     sizes[INPUT][1] = sizes[INPUT][0];
851     types[INPUT][1] = types[INPUT][0];
852 
853     sizes[INPUT][2] = cvSize(pt_count, 1 );
854     types[INPUT][2] = CV_64FC3;
855 
856     sizes[INPUT][3] = cvSize(4,3);
857     types[INPUT][3] = CV_64FC1;
858 
859     sizes[INPUT][4] = sizes[INPUT][5] = cvSize(3,3);
860     types[INPUT][4] = types[INPUT][5] = CV_MAKETYPE(CV_64F, 1);
861 
862     sizes[TEMP][0] = cvSize(3,3);
863     types[TEMP][0] = CV_64FC1;
864     sizes[TEMP][1] = cvSize(pt_count,1);
865     types[TEMP][1] = CV_8UC1;
866 
867     sizes[OUTPUT][0] = sizes[REF_OUTPUT][0] = cvSize(3,1);
868     types[OUTPUT][0] = types[REF_OUTPUT][0] = CV_64FC1;
869     sizes[OUTPUT][1] = sizes[REF_OUTPUT][1] = cvSize(pt_count,1);
870     types[OUTPUT][1] = types[REF_OUTPUT][1] = CV_8UC1;
871 
872     test_cpp = (cvtest::randInt(rng) & 256) == 0;
873 }
874 
875 
get_success_error_level(int,int,int)876 double CV_FundamentalMatTest::get_success_error_level( int /*test_case_idx*/, int /*i*/, int /*j*/ )
877 {
878     return 1e-2;
879 }
880 
881 
fill_array(int test_case_idx,int i,int j,Mat & arr)882 void CV_FundamentalMatTest::fill_array( int test_case_idx, int i, int j, Mat& arr )
883 {
884     double t[12]={0};
885     RNG& rng = ts->get_rng();
886 
887     if( i != INPUT )
888     {
889         cvtest::ArrayTest::fill_array( test_case_idx, i, j, arr );
890         return;
891     }
892 
893     switch( j )
894     {
895     case 0:
896     case 1:
897         return; // fill them later in prepare_test_case
898     case 2:
899         {
900         double* p = arr.ptr<double>();
901         for( i = 0; i < arr.cols*3; i += 3 )
902         {
903             p[i] = cvtest::randReal(rng)*cube_size;
904             p[i+1] = cvtest::randReal(rng)*cube_size;
905             p[i+2] = cvtest::randReal(rng)*cube_size + cube_size;
906         }
907         }
908         break;
909     case 3:
910         {
911         double r[3];
912         Mat rot_vec( 3, 1, CV_64F, r );
913         Mat rot_mat( 3, 3, CV_64F, t, 4*sizeof(t[0]) );
914         r[0] = cvtest::randReal(rng)*CV_PI*2;
915         r[1] = cvtest::randReal(rng)*CV_PI*2;
916         r[2] = cvtest::randReal(rng)*CV_PI*2;
917 
918         cvtest::Rodrigues( rot_vec, rot_mat );
919         t[3] = cvtest::randReal(rng)*cube_size;
920         t[7] = cvtest::randReal(rng)*cube_size;
921         t[11] = cvtest::randReal(rng)*cube_size;
922         Mat( 3, 4, CV_64F, t ).convertTo(arr, arr.type());
923         }
924         break;
925     case 4:
926     case 5:
927         t[0] = t[4] = cvtest::randReal(rng)*(max_f - min_f) + min_f;
928         t[2] = (img_size*0.5 + cvtest::randReal(rng)*4. - 2.)*t[0];
929         t[5] = (img_size*0.5 + cvtest::randReal(rng)*4. - 2.)*t[4];
930         t[8] = 1.;
931         Mat( 3, 3, CV_64F, t ).convertTo( arr, arr.type() );
932         break;
933     }
934 }
935 
936 
prepare_test_case(int test_case_idx)937 int CV_FundamentalMatTest::prepare_test_case( int test_case_idx )
938 {
939     int code = cvtest::ArrayTest::prepare_test_case( test_case_idx );
940     if( code > 0 )
941     {
942         const Mat& _3d = test_mat[INPUT][2];
943         RNG& rng = ts->get_rng();
944         double Idata[] = { 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0 };
945         Mat I( 3, 4, CV_64F, Idata );
946         int k;
947 
948         for( k = 0; k < 2; k++ )
949         {
950             const Mat& Rt = k == 0 ? I : test_mat[INPUT][3];
951             const Mat& A = test_mat[INPUT][k == 0 ? 4 : 5];
952             Mat& _2d = test_mat[INPUT][k];
953 
954             test_projectPoints( _3d, Rt, A, _2d, &rng, sigma );
955         }
956     }
957 
958     return code;
959 }
960 
run_func()961 void CV_FundamentalMatTest::run_func()
962 {
963     // cvFindFundamentalMat calls cv::findFundamentalMat
964     cv::Mat _input0 = test_mat[INPUT][0], _input1 = test_mat[INPUT][1];
965     cv::Mat& F = test_mat[TEMP][0], &mask = test_mat[TEMP][1];
966     F = cv::findFundamentalMat( _input0, _input1, method, MAX(sigma*3, 0.01), 0, mask );
967     f_result = !F.empty();
968 }
969 
prepare_to_validation(int test_case_idx)970 void CV_FundamentalMatTest::prepare_to_validation( int test_case_idx )
971 {
972     const Mat& Rt = test_mat[INPUT][3];
973     const Mat& A1 = test_mat[INPUT][4];
974     const Mat& A2 = test_mat[INPUT][5];
975     double f0[9], f[9];
976     Mat F0(3, 3, CV_64FC1, f0), F(3, 3, CV_64F, f);
977 
978     Mat invA1, invA2, R=Rt.colRange(0, 3), T;
979 
980     cv::invert(A1, invA1, CV_SVD);
981     cv::invert(A2, invA2, CV_SVD);
982 
983     double tx = Rt.at<double>(0, 3);
984     double ty = Rt.at<double>(1, 3);
985     double tz = Rt.at<double>(2, 3);
986 
987     double _t_x[] = { 0, -tz, ty, tz, 0, -tx, -ty, tx, 0 };
988 
989     // F = (A2^-T)*[t]_x*R*(A1^-1)
990     cv::gemm( invA2, Mat( 3, 3, CV_64F, _t_x ), 1, Mat(), 0, T, CV_GEMM_A_T );
991     cv::gemm( R, invA1, 1, Mat(), 0, invA2 );
992     cv::gemm( T, invA2, 1, Mat(), 0, F0 );
993     F0 *= 1./f0[8];
994 
995     uchar* status = test_mat[TEMP][1].ptr();
996     double err_level = get_success_error_level( test_case_idx, OUTPUT, 1 );
997     uchar* mtfm1 = test_mat[REF_OUTPUT][1].ptr();
998     uchar* mtfm2 = test_mat[OUTPUT][1].ptr();
999     double* f_prop1 = test_mat[REF_OUTPUT][0].ptr<double>();
1000     double* f_prop2 = test_mat[OUTPUT][0].ptr<double>();
1001 
1002     int i, pt_count = test_mat[INPUT][2].cols;
1003     Mat p1( 1, pt_count, CV_64FC2 );
1004     Mat p2( 1, pt_count, CV_64FC2 );
1005 
1006     test_convertHomogeneous( test_mat[INPUT][0], p1 );
1007     test_convertHomogeneous( test_mat[INPUT][1], p2 );
1008 
1009     Mat Fsrc = test_mat[TEMP][0];
1010     if( Fsrc.rows > 3 )
1011         Fsrc = Fsrc.rowRange(0, 3);
1012 
1013     cvtest::convert(Fsrc, F, F.type());
1014 
1015     if( method <= CV_FM_8POINT )
1016         memset( status, 1, pt_count );
1017 
1018     for( i = 0; i < pt_count; i++ )
1019     {
1020         double x1 = p1.at<Point2d>(i).x;
1021         double y1 = p1.at<Point2d>(i).y;
1022         double x2 = p2.at<Point2d>(i).x;
1023         double y2 = p2.at<Point2d>(i).y;
1024         double n1 = 1./sqrt(x1*x1 + y1*y1 + 1);
1025         double n2 = 1./sqrt(x2*x2 + y2*y2 + 1);
1026         double t0 = fabs(f0[0]*x2*x1 + f0[1]*x2*y1 + f0[2]*x2 +
1027                    f0[3]*y2*x1 + f0[4]*y2*y1 + f0[5]*y2 +
1028                    f0[6]*x1 + f0[7]*y1 + f0[8])*n1*n2;
1029         double t = fabs(f[0]*x2*x1 + f[1]*x2*y1 + f[2]*x2 +
1030                    f[3]*y2*x1 + f[4]*y2*y1 + f[5]*y2 +
1031                    f[6]*x1 + f[7]*y1 + f[8])*n1*n2;
1032         mtfm1[i] = 1;
1033         mtfm2[i] = !status[i] || t0 > err_level || t < err_level;
1034     }
1035 
1036     f_prop1[0] = 1;
1037     f_prop1[1] = 1;
1038     f_prop1[2] = 0;
1039 
1040     f_prop2[0] = f_result != 0;
1041     f_prop2[1] = f[8];
1042     f_prop2[2] = cv::determinant( F );
1043 }
1044 /******************************* find essential matrix ***********************************/
1045 class CV_EssentialMatTest : public cvtest::ArrayTest
1046 {
1047 public:
1048     CV_EssentialMatTest();
1049 
1050 protected:
1051     int read_params( const cv::FileStorage& fs );
1052     void fill_array( int test_case_idx, int i, int j, Mat& arr );
1053     int prepare_test_case( int test_case_idx );
1054     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
1055     double get_success_error_level( int test_case_idx, int i, int j );
1056     void run_func();
1057     void prepare_to_validation( int );
1058 
1059 #if 0
1060     double sampson_error(const double* f, double x1, double y1, double x2, double y2);
1061 #endif
1062 
1063     int method;
1064     int img_size;
1065     int cube_size;
1066     int dims;
1067     double min_f, max_f;
1068     double sigma;
1069 };
1070 
1071 
CV_EssentialMatTest()1072 CV_EssentialMatTest::CV_EssentialMatTest()
1073 {
1074     // input arrays:
1075     //   0, 1 - arrays of 2d points that are passed to %func%.
1076     //          Can have different data type, layout, be stored in homogeneous coordinates or not.
1077     //   2 - array of 3d points that are projected to both view planes
1078     //   3 - [R|t] matrix for the second view plane (for the first one it is [I|0]
1079     //   4 - intrinsic matrix for both camera
1080     test_array[INPUT].push_back(NULL);
1081     test_array[INPUT].push_back(NULL);
1082     test_array[INPUT].push_back(NULL);
1083     test_array[INPUT].push_back(NULL);
1084     test_array[INPUT].push_back(NULL);
1085     test_array[TEMP].push_back(NULL);
1086     test_array[TEMP].push_back(NULL);
1087     test_array[TEMP].push_back(NULL);
1088     test_array[TEMP].push_back(NULL);
1089     test_array[TEMP].push_back(NULL);
1090     test_array[OUTPUT].push_back(NULL); // Essential Matrix singularity
1091     test_array[OUTPUT].push_back(NULL); // Inliers mask
1092     test_array[OUTPUT].push_back(NULL); // Translation error
1093     test_array[OUTPUT].push_back(NULL); // Positive depth count
1094     test_array[REF_OUTPUT].push_back(NULL);
1095     test_array[REF_OUTPUT].push_back(NULL);
1096     test_array[REF_OUTPUT].push_back(NULL);
1097     test_array[REF_OUTPUT].push_back(NULL);
1098 
1099     element_wise_relative_error = false;
1100 
1101     method = 0;
1102     img_size = 10;
1103     cube_size = 10;
1104     dims = 0;
1105     min_f = 1;
1106     max_f = 3;
1107     sigma = 0;
1108 }
1109 
1110 
read_params(const cv::FileStorage & fs)1111 int CV_EssentialMatTest::read_params( const cv::FileStorage& fs )
1112 {
1113     int code = cvtest::ArrayTest::read_params( fs );
1114     return code;
1115 }
1116 
1117 
get_test_array_types_and_sizes(int,vector<vector<Size>> & sizes,vector<vector<int>> & types)1118 void CV_EssentialMatTest::get_test_array_types_and_sizes( int /*test_case_idx*/,
1119                                                 vector<vector<Size> >& sizes, vector<vector<int> >& types )
1120 {
1121     RNG& rng = ts->get_rng();
1122     int pt_depth = cvtest::randInt(rng) % 2 == 0 ? CV_32F : CV_64F;
1123     double pt_count_exp = cvtest::randReal(rng)*6 + 1;
1124     int pt_count = MAX(5, cvRound(exp(pt_count_exp)));
1125 
1126     dims = cvtest::randInt(rng) % 2 + 2;
1127     dims = 2;
1128     method = CV_LMEDS << (cvtest::randInt(rng) % 2);
1129 
1130     types[INPUT][0] = CV_MAKETYPE(pt_depth, 1);
1131 
1132     sizes[INPUT][0] = cvSize(dims, pt_count);
1133     if( cvtest::randInt(rng) % 2 )
1134     {
1135         types[INPUT][0] = CV_MAKETYPE(pt_depth, dims);
1136         if( cvtest::randInt(rng) % 2 )
1137             sizes[INPUT][0] = cvSize(pt_count, 1);
1138         else
1139             sizes[INPUT][0] = cvSize(1, pt_count);
1140     }
1141 
1142     sizes[INPUT][1] = sizes[INPUT][0];
1143     types[INPUT][1] = types[INPUT][0];
1144 
1145     sizes[INPUT][2] = cvSize(pt_count, 1 );
1146     types[INPUT][2] = CV_64FC3;
1147 
1148     sizes[INPUT][3] = cvSize(4,3);
1149     types[INPUT][3] = CV_64FC1;
1150 
1151     sizes[INPUT][4] = cvSize(3,3);
1152     types[INPUT][4] = CV_MAKETYPE(CV_64F, 1);
1153 
1154     sizes[TEMP][0] = cvSize(3,3);
1155     types[TEMP][0] = CV_64FC1;
1156     sizes[TEMP][1] = cvSize(pt_count,1);
1157     types[TEMP][1] = CV_8UC1;
1158     sizes[TEMP][2] = cvSize(3,3);
1159     types[TEMP][2] = CV_64FC1;
1160     sizes[TEMP][3] = cvSize(3, 1);
1161     types[TEMP][3] = CV_64FC1;
1162     sizes[TEMP][4] = cvSize(pt_count,1);
1163     types[TEMP][4] = CV_8UC1;
1164 
1165     sizes[OUTPUT][0] = sizes[REF_OUTPUT][0] = cvSize(3,1);
1166     types[OUTPUT][0] = types[REF_OUTPUT][0] = CV_64FC1;
1167     sizes[OUTPUT][1] = sizes[REF_OUTPUT][1] = cvSize(pt_count,1);
1168     types[OUTPUT][1] = types[REF_OUTPUT][1] = CV_8UC1;
1169     sizes[OUTPUT][2] = sizes[REF_OUTPUT][2] = cvSize(1,1);
1170     types[OUTPUT][2] = types[REF_OUTPUT][2] = CV_64FC1;
1171     sizes[OUTPUT][3] = sizes[REF_OUTPUT][3] = cvSize(1,1);
1172     types[OUTPUT][3] = types[REF_OUTPUT][3] = CV_8UC1;
1173 
1174 }
1175 
1176 
get_success_error_level(int,int,int)1177 double CV_EssentialMatTest::get_success_error_level( int /*test_case_idx*/, int /*i*/, int /*j*/ )
1178 {
1179     return 1e-2;
1180 }
1181 
1182 
fill_array(int test_case_idx,int i,int j,Mat & arr)1183 void CV_EssentialMatTest::fill_array( int test_case_idx, int i, int j, Mat& arr )
1184 {
1185     double t[12]={0};
1186     RNG& rng = ts->get_rng();
1187 
1188     if( i != INPUT )
1189     {
1190         cvtest::ArrayTest::fill_array( test_case_idx, i, j, arr );
1191         return;
1192     }
1193 
1194     switch( j )
1195     {
1196     case 0:
1197     case 1:
1198         return; // fill them later in prepare_test_case
1199     case 2:
1200         {
1201         double* p = arr.ptr<double>();
1202         for( i = 0; i < arr.cols*3; i += 3 )
1203         {
1204             p[i] = cvtest::randReal(rng)*cube_size;
1205             p[i+1] = cvtest::randReal(rng)*cube_size;
1206             p[i+2] = cvtest::randReal(rng)*cube_size + cube_size;
1207         }
1208         }
1209         break;
1210     case 3:
1211         {
1212         double r[3];
1213         Mat rot_vec( 3, 1, CV_64F, r );
1214         Mat rot_mat( 3, 3, CV_64F, t, 4*sizeof(t[0]) );
1215         r[0] = cvtest::randReal(rng)*CV_PI*2;
1216         r[1] = cvtest::randReal(rng)*CV_PI*2;
1217         r[2] = cvtest::randReal(rng)*CV_PI*2;
1218 
1219         cvtest::Rodrigues( rot_vec, rot_mat );
1220         t[3] = cvtest::randReal(rng)*cube_size;
1221         t[7] = cvtest::randReal(rng)*cube_size;
1222         t[11] = cvtest::randReal(rng)*cube_size;
1223         Mat( 3, 4, CV_64F, t ).convertTo(arr, arr.type());
1224         }
1225         break;
1226     case 4:
1227         t[0] = t[4] = cvtest::randReal(rng)*(max_f - min_f) + min_f;
1228         t[2] = (img_size*0.5 + cvtest::randReal(rng)*4. - 2.)*t[0];
1229         t[5] = (img_size*0.5 + cvtest::randReal(rng)*4. - 2.)*t[4];
1230         t[8] = 1.;
1231         Mat( 3, 3, CV_64F, t ).convertTo( arr, arr.type() );
1232         break;
1233     }
1234 }
1235 
1236 
prepare_test_case(int test_case_idx)1237 int CV_EssentialMatTest::prepare_test_case( int test_case_idx )
1238 {
1239     int code = cvtest::ArrayTest::prepare_test_case( test_case_idx );
1240     if( code > 0 )
1241     {
1242         const Mat& _3d = test_mat[INPUT][2];
1243         RNG& rng = ts->get_rng();
1244         double Idata[] = { 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0 };
1245         Mat I( 3, 4, CV_64F, Idata );
1246         int k;
1247 
1248         for( k = 0; k < 2; k++ )
1249         {
1250             const Mat& Rt = k == 0 ? I : test_mat[INPUT][3];
1251             const Mat& A = test_mat[INPUT][4];
1252             Mat& _2d = test_mat[INPUT][k];
1253 
1254             test_projectPoints( _3d, Rt, A, _2d, &rng, sigma );
1255         }
1256     }
1257 
1258     return code;
1259 }
1260 
1261 
run_func()1262 void CV_EssentialMatTest::run_func()
1263 {
1264     Mat _input0(test_mat[INPUT][0]), _input1(test_mat[INPUT][1]);
1265     Mat K(test_mat[INPUT][4]);
1266     double focal(K.at<double>(0, 0));
1267     cv::Point2d pp(K.at<double>(0, 2), K.at<double>(1, 2));
1268 
1269     RNG& rng = ts->get_rng();
1270     Mat E, mask1(test_mat[TEMP][1]);
1271     E = cv::findEssentialMat( _input0, _input1, focal, pp, method, 0.99, MAX(sigma*3, 0.0001), mask1 );
1272     if (E.rows > 3)
1273     {
1274         int count = E.rows / 3;
1275         int row = (cvtest::randInt(rng) % count) * 3;
1276         E = E.rowRange(row, row + 3) * 1.0;
1277     }
1278 
1279     E.copyTo(test_mat[TEMP][0]);
1280 
1281     Mat R, t, mask2;
1282     recoverPose( E, _input0, _input1, R, t, focal, pp, mask2 );
1283     R.copyTo(test_mat[TEMP][2]);
1284     t.copyTo(test_mat[TEMP][3]);
1285     mask2.copyTo(test_mat[TEMP][4]);
1286 }
1287 
1288 #if 0
1289 double CV_EssentialMatTest::sampson_error(const double * f, double x1, double y1, double x2, double y2)
1290 {
1291     double Fx1[3] = {
1292         f[0] * x1 + f[1] * y1 + f[2],
1293         f[3] * x1 + f[4] * y1 + f[5],
1294         f[6] * x1 + f[7] * y1 + f[8]
1295     };
1296     double Ftx2[3] = {
1297         f[0] * x2 + f[3] * y2 + f[6],
1298         f[1] * x2 + f[4] * y2 + f[7],
1299         f[2] * x2 + f[5] * y2 + f[8]
1300     };
1301     double x2tFx1 = Fx1[0] * x2 + Fx1[1] * y2 + Fx1[2];
1302 
1303     double error = x2tFx1 * x2tFx1 / (Fx1[0] * Fx1[0] + Fx1[1] * Fx1[1] + Ftx2[0] * Ftx2[0] + Ftx2[1] * Ftx2[1]);
1304     error = sqrt(error);
1305     return error;
1306 }
1307 #endif
1308 
prepare_to_validation(int test_case_idx)1309 void CV_EssentialMatTest::prepare_to_validation( int test_case_idx )
1310 {
1311     const Mat& Rt0 = test_mat[INPUT][3];
1312     const Mat& A = test_mat[INPUT][4];
1313     double f0[9], f[9], e[9];
1314     Mat F0(3, 3, CV_64FC1, f0), F(3, 3, CV_64F, f);
1315     Mat E(3, 3, CV_64F, e);
1316 
1317     Mat invA, R=Rt0.colRange(0, 3), T1, T2;
1318 
1319     cv::invert(A, invA, CV_SVD);
1320 
1321     double tx = Rt0.at<double>(0, 3);
1322     double ty = Rt0.at<double>(1, 3);
1323     double tz = Rt0.at<double>(2, 3);
1324 
1325     double _t_x[] = { 0, -tz, ty, tz, 0, -tx, -ty, tx, 0 };
1326 
1327     // F = (A2^-T)*[t]_x*R*(A1^-1)
1328     cv::gemm( invA, Mat( 3, 3, CV_64F, _t_x ), 1, Mat(), 0, T1, CV_GEMM_A_T );
1329     cv::gemm( R, invA, 1, Mat(), 0, T2 );
1330     cv::gemm( T1, T2, 1, Mat(), 0, F0 );
1331     F0 *= 1./f0[8];
1332 
1333     uchar* status = test_mat[TEMP][1].ptr();
1334     double err_level = get_success_error_level( test_case_idx, OUTPUT, 1 );
1335     uchar* mtfm1 = test_mat[REF_OUTPUT][1].ptr();
1336     uchar* mtfm2 = test_mat[OUTPUT][1].ptr();
1337     double* e_prop1 = test_mat[REF_OUTPUT][0].ptr<double>();
1338     double* e_prop2 = test_mat[OUTPUT][0].ptr<double>();
1339     Mat E_prop2 = Mat(3, 1, CV_64F, e_prop2);
1340 
1341     int i, pt_count = test_mat[INPUT][2].cols;
1342     Mat p1( 1, pt_count, CV_64FC2 );
1343     Mat p2( 1, pt_count, CV_64FC2 );
1344 
1345     test_convertHomogeneous( test_mat[INPUT][0], p1 );
1346     test_convertHomogeneous( test_mat[INPUT][1], p2 );
1347 
1348     cvtest::convert(test_mat[TEMP][0], E, E.type());
1349     cv::gemm( invA, E, 1, Mat(), 0, T1, CV_GEMM_A_T );
1350     cv::gemm( T1, invA, 1, Mat(), 0, F );
1351 
1352     for( i = 0; i < pt_count; i++ )
1353     {
1354         double x1 = p1.at<Point2d>(i).x;
1355         double y1 = p1.at<Point2d>(i).y;
1356         double x2 = p2.at<Point2d>(i).x;
1357         double y2 = p2.at<Point2d>(i).y;
1358 //        double t0 = sampson_error(f0, x1, y1, x2, y2);
1359 //        double t = sampson_error(f, x1, y1, x2, y2);
1360         double n1 = 1./sqrt(x1*x1 + y1*y1 + 1);
1361         double n2 = 1./sqrt(x2*x2 + y2*y2 + 1);
1362         double t0 = fabs(f0[0]*x2*x1 + f0[1]*x2*y1 + f0[2]*x2 +
1363                    f0[3]*y2*x1 + f0[4]*y2*y1 + f0[5]*y2 +
1364                    f0[6]*x1 + f0[7]*y1 + f0[8])*n1*n2;
1365         double t = fabs(f[0]*x2*x1 + f[1]*x2*y1 + f[2]*x2 +
1366                    f[3]*y2*x1 + f[4]*y2*y1 + f[5]*y2 +
1367                    f[6]*x1 + f[7]*y1 + f[8])*n1*n2;
1368         mtfm1[i] = 1;
1369         mtfm2[i] = !status[i] || t0 > err_level || t < err_level;
1370     }
1371 
1372     e_prop1[0] = sqrt(0.5);
1373     e_prop1[1] = sqrt(0.5);
1374     e_prop1[2] = 0;
1375 
1376     e_prop2[0] = 0;
1377     e_prop2[1] = 0;
1378     e_prop2[2] = 0;
1379     SVD::compute(E, E_prop2);
1380 
1381 
1382 
1383     double* pose_prop1 = test_mat[REF_OUTPUT][2].ptr<double>();
1384     double* pose_prop2 = test_mat[OUTPUT][2].ptr<double>();
1385     double terr1 = cvtest::norm(Rt0.col(3) / cvtest::norm(Rt0.col(3), NORM_L2) + test_mat[TEMP][3], NORM_L2);
1386     double terr2 = cvtest::norm(Rt0.col(3) / cvtest::norm(Rt0.col(3), NORM_L2) - test_mat[TEMP][3], NORM_L2);
1387     Mat rvec(3, 1, CV_32F);
1388     cvtest::Rodrigues(Rt0.colRange(0, 3), rvec);
1389     pose_prop1[0] = 0;
1390     // No check for CV_LMeDS on translation. Since it
1391     // involves with some degraded problem, when data is exact inliers.
1392     pose_prop2[0] = method == CV_LMEDS || pt_count == 5 ? 0 : MIN(terr1, terr2);
1393 
1394 
1395 //    int inliers_count = countNonZero(test_mat[TEMP][1]);
1396 //    int good_count = countNonZero(test_mat[TEMP][4]);
1397     test_mat[OUTPUT][3] = true; //good_count >= inliers_count / 2;
1398     test_mat[REF_OUTPUT][3] = true;
1399 
1400 
1401 }
1402 
1403 
1404 /********************************** convert homogeneous *********************************/
1405 
1406 class CV_ConvertHomogeneousTest : public cvtest::ArrayTest
1407 {
1408 public:
1409     CV_ConvertHomogeneousTest();
1410 
1411 protected:
1412     int read_params( const cv::FileStorage& fs );
1413     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
1414     void fill_array( int test_case_idx, int i, int j, Mat& arr );
1415     double get_success_error_level( int test_case_idx, int i, int j );
1416     void run_func();
1417     void prepare_to_validation( int );
1418 
1419     int dims1, dims2;
1420     int pt_count;
1421 };
1422 
1423 
CV_ConvertHomogeneousTest()1424 CV_ConvertHomogeneousTest::CV_ConvertHomogeneousTest()
1425 {
1426     test_array[INPUT].push_back(NULL);
1427     test_array[OUTPUT].push_back(NULL);
1428     test_array[REF_OUTPUT].push_back(NULL);
1429     element_wise_relative_error = false;
1430 
1431     pt_count = dims1 = dims2 = 0;
1432 }
1433 
1434 
read_params(const cv::FileStorage & fs)1435 int CV_ConvertHomogeneousTest::read_params( const cv::FileStorage& fs )
1436 {
1437     int code = cvtest::ArrayTest::read_params( fs );
1438     return code;
1439 }
1440 
1441 
get_test_array_types_and_sizes(int,vector<vector<Size>> & sizes,vector<vector<int>> & types)1442 void CV_ConvertHomogeneousTest::get_test_array_types_and_sizes( int /*test_case_idx*/,
1443                                                 vector<vector<Size> >& sizes, vector<vector<int> >& types )
1444 {
1445     RNG& rng = ts->get_rng();
1446     int pt_depth1 = cvtest::randInt(rng) % 2 == 0 ? CV_32F : CV_64F;
1447     int pt_depth2 = pt_depth1;//cvtest::randInt(rng) % 2 == 0 ? CV_32F : CV_64F;
1448     double pt_count_exp = cvtest::randReal(rng)*6 + 1;
1449     int t;
1450 
1451     pt_count = cvRound(exp(pt_count_exp));
1452     pt_count = MAX( pt_count, 5 );
1453 
1454     dims1 = 2 + (cvtest::randInt(rng) % 2);
1455     dims2 = dims1 + 1;
1456 
1457     if( cvtest::randInt(rng) % 2 )
1458         CV_SWAP( dims1, dims2, t );
1459 
1460     types[INPUT][0] = CV_MAKETYPE(pt_depth1, 1);
1461 
1462     sizes[INPUT][0] = cvSize(dims1, pt_count);
1463     if( cvtest::randInt(rng) % 2 )
1464     {
1465         types[INPUT][0] = CV_MAKETYPE(pt_depth1, dims1);
1466         if( cvtest::randInt(rng) % 2 )
1467             sizes[INPUT][0] = cvSize(pt_count, 1);
1468         else
1469             sizes[INPUT][0] = cvSize(1, pt_count);
1470     }
1471 
1472     types[OUTPUT][0] = CV_MAKETYPE(pt_depth2, dims2);
1473     sizes[OUTPUT][0] = cvSize(1, pt_count);
1474 
1475     types[REF_OUTPUT][0] = types[OUTPUT][0];
1476     sizes[REF_OUTPUT][0] = sizes[OUTPUT][0];
1477 }
1478 
1479 
get_success_error_level(int,int,int)1480 double CV_ConvertHomogeneousTest::get_success_error_level( int /*test_case_idx*/, int /*i*/, int /*j*/ )
1481 {
1482     return 1e-5;
1483 }
1484 
1485 
fill_array(int,int,int,Mat & arr)1486 void CV_ConvertHomogeneousTest::fill_array( int /*test_case_idx*/, int /*i*/, int /*j*/, Mat& arr )
1487 {
1488     Mat temp( 1, pt_count, CV_MAKETYPE(CV_64FC1,dims1) );
1489     RNG& rng = ts->get_rng();
1490     CvScalar low = cvScalarAll(0), high = cvScalarAll(10);
1491 
1492     if( dims1 > dims2 )
1493         low.val[dims1-1] = 1.;
1494 
1495     cvtest::randUni( rng, temp, low, high );
1496     test_convertHomogeneous( temp, arr );
1497 }
1498 
1499 
run_func()1500 void CV_ConvertHomogeneousTest::run_func()
1501 {
1502     cv::Mat _input = test_mat[INPUT][0], &_output = test_mat[OUTPUT][0];
1503     if( dims1 > dims2 )
1504         cv::convertPointsFromHomogeneous(_input, _output);
1505     else
1506         cv::convertPointsToHomogeneous(_input, _output);
1507 }
1508 
1509 
prepare_to_validation(int)1510 void CV_ConvertHomogeneousTest::prepare_to_validation( int /*test_case_idx*/ )
1511 {
1512     test_convertHomogeneous( test_mat[INPUT][0], test_mat[REF_OUTPUT][0] );
1513 }
1514 
1515 
1516 /************************** compute corresponding epipolar lines ************************/
1517 
1518 class CV_ComputeEpilinesTest : public cvtest::ArrayTest
1519 {
1520 public:
1521     CV_ComputeEpilinesTest();
1522 
1523 protected:
1524     int read_params( const cv::FileStorage& fs );
1525     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
1526     void fill_array( int test_case_idx, int i, int j, Mat& arr );
1527     double get_success_error_level( int test_case_idx, int i, int j );
1528     void run_func();
1529     void prepare_to_validation( int );
1530 
1531     int which_image;
1532     int dims;
1533     int pt_count;
1534 };
1535 
1536 
CV_ComputeEpilinesTest()1537 CV_ComputeEpilinesTest::CV_ComputeEpilinesTest()
1538 {
1539     test_array[INPUT].push_back(NULL);
1540     test_array[INPUT].push_back(NULL);
1541     test_array[OUTPUT].push_back(NULL);
1542     test_array[REF_OUTPUT].push_back(NULL);
1543     element_wise_relative_error = false;
1544 
1545     pt_count = dims = which_image = 0;
1546 }
1547 
1548 
read_params(const cv::FileStorage & fs)1549 int CV_ComputeEpilinesTest::read_params( const cv::FileStorage& fs )
1550 {
1551     int code = cvtest::ArrayTest::read_params( fs );
1552     return code;
1553 }
1554 
1555 
get_test_array_types_and_sizes(int,vector<vector<Size>> & sizes,vector<vector<int>> & types)1556 void CV_ComputeEpilinesTest::get_test_array_types_and_sizes( int /*test_case_idx*/,
1557                                                 vector<vector<Size> >& sizes, vector<vector<int> >& types )
1558 {
1559     RNG& rng = ts->get_rng();
1560     int fm_depth = cvtest::randInt(rng) % 2 == 0 ? CV_32F : CV_64F;
1561     int pt_depth = cvtest::randInt(rng) % 2 == 0 ? CV_32F : CV_64F;
1562     int ln_depth = pt_depth;
1563     double pt_count_exp = cvtest::randReal(rng)*6;
1564 
1565     which_image = 1 + (cvtest::randInt(rng) % 2);
1566 
1567     pt_count = cvRound(exp(pt_count_exp));
1568     pt_count = MAX( pt_count, 1 );
1569     bool few_points = pt_count < 5;
1570 
1571     dims = 2 + (cvtest::randInt(rng) % 2);
1572 
1573     types[INPUT][0] = CV_MAKETYPE(pt_depth, 1);
1574 
1575     sizes[INPUT][0] = cvSize(dims, pt_count);
1576     if( cvtest::randInt(rng) % 2 || few_points )
1577     {
1578         types[INPUT][0] = CV_MAKETYPE(pt_depth, dims);
1579         if( cvtest::randInt(rng) % 2 )
1580             sizes[INPUT][0] = cvSize(pt_count, 1);
1581         else
1582             sizes[INPUT][0] = cvSize(1, pt_count);
1583     }
1584 
1585     types[INPUT][1] = CV_MAKETYPE(fm_depth, 1);
1586     sizes[INPUT][1] = cvSize(3, 3);
1587 
1588     types[OUTPUT][0] = CV_MAKETYPE(ln_depth, 3);
1589     sizes[OUTPUT][0] = cvSize(1, pt_count);
1590 
1591     types[REF_OUTPUT][0] = types[OUTPUT][0];
1592     sizes[REF_OUTPUT][0] = sizes[OUTPUT][0];
1593 }
1594 
1595 
get_success_error_level(int,int,int)1596 double CV_ComputeEpilinesTest::get_success_error_level( int /*test_case_idx*/, int /*i*/, int /*j*/ )
1597 {
1598     return 1e-5;
1599 }
1600 
1601 
fill_array(int test_case_idx,int i,int j,Mat & arr)1602 void CV_ComputeEpilinesTest::fill_array( int test_case_idx, int i, int j, Mat& arr )
1603 {
1604     RNG& rng = ts->get_rng();
1605 
1606     if( i == INPUT && j == 0 )
1607     {
1608         Mat temp( 1, pt_count, CV_MAKETYPE(CV_64FC1,dims) );
1609         cvtest::randUni( rng, temp, cvScalar(0,0,1), cvScalarAll(10) );
1610         test_convertHomogeneous( temp, arr );
1611     }
1612     else if( i == INPUT && j == 1 )
1613         cvtest::randUni( rng, arr, cvScalarAll(0), cvScalarAll(10) );
1614     else
1615         cvtest::ArrayTest::fill_array( test_case_idx, i, j, arr );
1616 }
1617 
1618 
run_func()1619 void CV_ComputeEpilinesTest::run_func()
1620 {
1621     cv::Mat _points = test_mat[INPUT][0], _F = test_mat[INPUT][1], &_lines = test_mat[OUTPUT][0];
1622     cv::computeCorrespondEpilines( _points, which_image, _F, _lines );
1623 }
1624 
1625 
prepare_to_validation(int)1626 void CV_ComputeEpilinesTest::prepare_to_validation( int /*test_case_idx*/ )
1627 {
1628     Mat pt( 1, pt_count, CV_MAKETYPE(CV_64F, 3) );
1629     Mat lines( 1, pt_count, CV_MAKETYPE(CV_64F, 3) );
1630     double f[9];
1631     Mat F( 3, 3, CV_64F, f );
1632 
1633     test_convertHomogeneous( test_mat[INPUT][0], pt );
1634     test_mat[INPUT][1].convertTo(F, CV_64F);
1635     if( which_image == 2 )
1636         cv::transpose( F, F );
1637 
1638     for( int i = 0; i < pt_count; i++ )
1639     {
1640         double* p = pt.ptr<double>() + i*3;
1641         double* l = lines.ptr<double>() + i*3;
1642         double t0 = f[0]*p[0] + f[1]*p[1] + f[2]*p[2];
1643         double t1 = f[3]*p[0] + f[4]*p[1] + f[5]*p[2];
1644         double t2 = f[6]*p[0] + f[7]*p[1] + f[8]*p[2];
1645         double d = sqrt(t0*t0 + t1*t1);
1646         d = d ? 1./d : 1.;
1647         l[0] = t0*d; l[1] = t1*d; l[2] = t2*d;
1648     }
1649 
1650     test_convertHomogeneous( lines, test_mat[REF_OUTPUT][0] );
1651 }
1652 
TEST(Calib3d_Rodrigues,accuracy)1653 TEST(Calib3d_Rodrigues, accuracy) { CV_RodriguesTest test; test.safe_run(); }
TEST(Calib3d_FindFundamentalMat,accuracy)1654 TEST(Calib3d_FindFundamentalMat, accuracy) { CV_FundamentalMatTest test; test.safe_run(); }
TEST(Calib3d_ConvertHomogeneoous,accuracy)1655 TEST(Calib3d_ConvertHomogeneoous, accuracy) { CV_ConvertHomogeneousTest test; test.safe_run(); }
TEST(Calib3d_ComputeEpilines,accuracy)1656 TEST(Calib3d_ComputeEpilines, accuracy) { CV_ComputeEpilinesTest test; test.safe_run(); }
TEST(Calib3d_FindEssentialMat,accuracy)1657 TEST(Calib3d_FindEssentialMat, accuracy) { CV_EssentialMatTest test; test.safe_run(); }
1658 
TEST(Calib3d_FindFundamentalMat,correctMatches)1659 TEST(Calib3d_FindFundamentalMat, correctMatches)
1660 {
1661     double fdata[] = {0, 0, 0, 0, 0, -1, 0, 1, 0};
1662     double p1data[] = {200, 0, 1};
1663     double p2data[] = {170, 0, 1};
1664 
1665     Mat F(3, 3, CV_64F, fdata);
1666     Mat p1(1, 1, CV_64FC2, p1data);
1667     Mat p2(1, 1, CV_64FC2, p2data);
1668     Mat np1, np2;
1669 
1670     correctMatches(F, p1, p2, np1, np2);
1671 
1672     cout << np1 << endl;
1673     cout << np2 << endl;
1674 }
1675 
1676 }} // namespace
1677 /* End of file. */
1678