diff --git a/modules/imgproc/src/shapedescr.cpp b/modules/imgproc/src/shapedescr.cpp index 6107596560..efc11b09a5 100644 --- a/modules/imgproc/src/shapedescr.cpp +++ b/modules/imgproc/src/shapedescr.cpp @@ -771,202 +771,6 @@ cvContourArea( const void *array, CvSlice slice, int oriented ) } -/* for now this function works bad with singular cases - You can see in the code, that when some troubles with - matrices or some variables occur - - box filled with zero values is returned. - However in general function works fine. -*/ -static void -icvFitEllipse_F( CvSeq* points, CvBox2D* box ) -{ - cv::Ptr D; - double S[36], C[36], T[36]; - - int i, j; - double eigenvalues[6], eigenvectors[36]; - double a, b, c, d, e, f; - double x0, y0, idet, scale, offx = 0, offy = 0; - - int n = points->total; - CvSeqReader reader; - int is_float = CV_SEQ_ELTYPE(points) == CV_32FC2; - - CvMat matS = cvMat(6,6,CV_64F,S), matC = cvMat(6,6,CV_64F,C), matT = cvMat(6,6,CV_64F,T); - CvMat _EIGVECS = cvMat(6,6,CV_64F,eigenvectors), _EIGVALS = cvMat(6,1,CV_64F,eigenvalues); - - /* create matrix D of input points */ - D = cvCreateMat( n, 6, CV_64F ); - - cvStartReadSeq( points, &reader ); - - /* shift all points to zero */ - for( i = 0; i < n; i++ ) - { - if( !is_float ) - { - offx += ((CvPoint*)reader.ptr)->x; - offy += ((CvPoint*)reader.ptr)->y; - } - else - { - offx += ((CvPoint2D32f*)reader.ptr)->x; - offy += ((CvPoint2D32f*)reader.ptr)->y; - } - CV_NEXT_SEQ_ELEM( points->elem_size, reader ); - } - - offx /= n; - offy /= n; - - // fill matrix rows as (x*x, x*y, y*y, x, y, 1 ) - for( i = 0; i < n; i++ ) - { - double x, y; - double* Dptr = D->data.db + i*6; - - if( !is_float ) - { - x = ((CvPoint*)reader.ptr)->x - offx; - y = ((CvPoint*)reader.ptr)->y - offy; - } - else - { - x = ((CvPoint2D32f*)reader.ptr)->x - offx; - y = ((CvPoint2D32f*)reader.ptr)->y - offy; - } - CV_NEXT_SEQ_ELEM( points->elem_size, reader ); - - Dptr[0] = x * x; - Dptr[1] = x * y; - Dptr[2] = y * y; - Dptr[3] = x; - Dptr[4] = y; - Dptr[5] = 1.; - } - - // S = D^t*D - cvMulTransposed( D, &matS, 1 ); - cvSVD( &matS, &_EIGVALS, &_EIGVECS, 0, CV_SVD_MODIFY_A + CV_SVD_U_T ); - - for( i = 0; i < 6; i++ ) - { - double a = eigenvalues[i]; - a = a < DBL_EPSILON ? 0 : 1./sqrt(sqrt(a)); - for( j = 0; j < 6; j++ ) - eigenvectors[i*6 + j] *= a; - } - - // C = Q^-1 = transp(INVEIGV) * INVEIGV - cvMulTransposed( &_EIGVECS, &matC, 1 ); - - cvZero( &matS ); - S[2] = 2.; - S[7] = -1.; - S[12] = 2.; - - // S = Q^-1*S*Q^-1 - cvMatMul( &matC, &matS, &matT ); - cvMatMul( &matT, &matC, &matS ); - - // and find its eigenvalues and vectors too - //cvSVD( &matS, &_EIGVALS, &_EIGVECS, 0, CV_SVD_MODIFY_A + CV_SVD_U_T ); - cvEigenVV( &matS, &_EIGVECS, &_EIGVALS, 0 ); - - for( i = 0; i < 3; i++ ) - if( eigenvalues[i] > 0 ) - break; - - if( i >= 3 /*eigenvalues[0] < DBL_EPSILON*/ ) - { - box->center.x = box->center.y = - box->size.width = box->size.height = - box->angle = 0.f; - return; - } - - // now find truthful eigenvector - _EIGVECS = cvMat( 6, 1, CV_64F, eigenvectors + 6*i ); - matT = cvMat( 6, 1, CV_64F, T ); - // Q^-1*eigenvecs[0] - cvMatMul( &matC, &_EIGVECS, &matT ); - - // extract vector components - a = T[0]; b = T[1]; c = T[2]; d = T[3]; e = T[4]; f = T[5]; - - ///////////////// extract ellipse axes from above values //////////////// - - /* - 1) find center of ellipse - it satisfy equation - | a b/2 | * | x0 | + | d/2 | = |0 | - | b/2 c | | y0 | | e/2 | |0 | - - */ - idet = a * c - b * b * 0.25; - idet = idet > DBL_EPSILON ? 1./idet : 0; - - // we must normalize (a b c d e f ) to fit (4ac-b^2=1) - scale = sqrt( 0.25 * idet ); - - if( scale < DBL_EPSILON ) - { - box->center.x = (float)offx; - box->center.y = (float)offy; - box->size.width = box->size.height = box->angle = 0.f; - return; - } - - a *= scale; - b *= scale; - c *= scale; - d *= scale; - e *= scale; - f *= scale; - - x0 = (-d * c + e * b * 0.5) * 2.; - y0 = (-a * e + d * b * 0.5) * 2.; - - // recover center - box->center.x = (float)(x0 + offx); - box->center.y = (float)(y0 + offy); - - // offset ellipse to (x0,y0) - // new f == F(x0,y0) - f += a * x0 * x0 + b * x0 * y0 + c * y0 * y0 + d * x0 + e * y0; - - if( fabs(f) < DBL_EPSILON ) - { - box->size.width = box->size.height = box->angle = 0.f; - return; - } - - scale = -1. / f; - // normalize to f = 1 - a *= scale; - b *= scale; - c *= scale; - - // extract axis of ellipse - // one more eigenvalue operation - S[0] = a; - S[1] = S[2] = b * 0.5; - S[3] = c; - - matS = cvMat( 2, 2, CV_64F, S ); - _EIGVECS = cvMat( 2, 2, CV_64F, eigenvectors ); - _EIGVALS = cvMat( 1, 2, CV_64F, eigenvalues ); - cvSVD( &matS, &_EIGVALS, &_EIGVECS, 0, CV_SVD_MODIFY_A + CV_SVD_U_T ); - - // exteract axis length from eigenvectors - box->size.width = (float)(2./sqrt(eigenvalues[0])); - box->size.height = (float)(2./sqrt(eigenvalues[1])); - - // calc angle - box->angle = (float)(180 - atan2(eigenvectors[2], eigenvectors[3])*180/CV_PI); -} - - CV_IMPL CvBox2D cvFitEllipse2( const CvArr* array ) { @@ -993,12 +797,11 @@ cvFitEllipse2( const CvArr* array ) n = ptseq->total; if( n < 5 ) CV_Error( CV_StsBadSize, "Number of points should be >= 5" ); -#if 1 - icvFitEllipse_F( ptseq, &box ); -#else + /* * New fitellipse algorithm, contributed by Dr. Daniel Weiss */ + CvPoint2D32f c = {0,0}; double gfp[5], rp[5], t; CvMat A, b, x; const double min_eps = 1e-6; @@ -1015,6 +818,23 @@ cvFitEllipse2( const CvArr* array ) cvStartReadSeq( ptseq, &reader ); is_float = CV_SEQ_ELTYPE(ptseq) == CV_32FC2; + + for( i = 0; i < n; i++ ) + { + CvPoint2D32f p; + if( is_float ) + p = *(CvPoint2D32f*)(reader.ptr); + else + { + p.x = (float)((int*)reader.ptr)[0]; + p.y = (float)((int*)reader.ptr)[1]; + } + CV_NEXT_SEQ_ELEM( sizeof(p), reader ); + c.x += p.x; + c.y += p.y; + } + c.x /= n; + c.y /= n; for( i = 0; i < n; i++ ) { @@ -1027,6 +847,8 @@ cvFitEllipse2( const CvArr* array ) p.y = (float)((int*)reader.ptr)[1]; } CV_NEXT_SEQ_ELEM( sizeof(p), reader ); + p.x -= c.x; + p.y -= c.y; bd[i] = 10000.0; // 1.0? Ad[i*5] = -(double)p.x * p.x; // A - C signs inverted as proposed by APP @@ -1065,6 +887,8 @@ cvFitEllipse2( const CvArr* array ) p.y = (float)((int*)reader.ptr)[1]; } CV_NEXT_SEQ_ELEM( sizeof(p), reader ); + p.x -= c.x; + p.y -= c.y; bd[i] = 1.0; Ad[i * 3] = (p.x - rp[0]) * (p.x - rp[0]); Ad[i * 3 + 1] = (p.y - rp[1]) * (p.y - rp[1]); @@ -1086,8 +910,8 @@ cvFitEllipse2( const CvArr* array ) if( rp[3] > min_eps ) rp[3] = sqrt(2.0 / rp[3]); - box.center.x = (float)rp[0]; - box.center.y = (float)rp[1]; + box.center.x = (float)rp[0] + c.x; + box.center.y = (float)rp[1] + c.y; box.size.width = (float)(rp[2]*2); box.size.height = (float)(rp[3]*2); if( box.size.width > box.size.height ) @@ -1100,7 +924,6 @@ cvFitEllipse2( const CvArr* array ) box.angle += 360; if( box.angle > 360 ) box.angle -= 360; -#endif return box; } diff --git a/modules/imgproc/test/test_convhull.cpp b/modules/imgproc/test/test_convhull.cpp index f8998c865e..84fa0cf5fe 100644 --- a/modules/imgproc/test/test_convhull.cpp +++ b/modules/imgproc/test/test_convhull.cpp @@ -1134,6 +1134,8 @@ int CV_FitEllipseTest::validate_test_results( int test_case_idx ) _exit_: #if 0 + if( code < 0 ) + { cvNamedWindow( "test", 0 ); IplImage* img = cvCreateImage( cvSize(cvRound(low_high_range*4), cvRound(low_high_range*4)), 8, 3 ); @@ -1154,6 +1156,7 @@ _exit_: cvShowImage( "test", img ); cvReleaseImage( &img ); cvWaitKey(0); + } #endif if( code < 0 ) @@ -1164,6 +1167,36 @@ _exit_: } +class CV_FitEllipseSmallTest : public cvtest::BaseTest +{ +public: + CV_FitEllipseSmallTest() {} + ~CV_FitEllipseSmallTest() {} +protected: + void run(int) + { + Size sz(50, 50); + vector > c; + c.push_back(vector()); + int scale = 1; + Point ofs = Point(0,0);//sz.width/2, sz.height/2) - Point(4,4)*scale; + c[0].push_back(Point(2, 0)*scale+ofs); + c[0].push_back(Point(0, 2)*scale+ofs); + c[0].push_back(Point(0, 6)*scale+ofs); + c[0].push_back(Point(2, 8)*scale+ofs); + c[0].push_back(Point(6, 8)*scale+ofs); + c[0].push_back(Point(8, 6)*scale+ofs); + c[0].push_back(Point(8, 2)*scale+ofs); + c[0].push_back(Point(6, 0)*scale+ofs); + + RotatedRect e = fitEllipse(c[0]); + CV_Assert( fabs(e.center.x - 4) <= 1. && + fabs(e.center.y - 4) <= 1. && + fabs(e.size.width - 9) <= 1. && + fabs(e.size.height - 9) <= 1. ); + } +}; + /****************************************************************************************\ * FitLine Test * \****************************************************************************************/ @@ -1664,6 +1697,7 @@ TEST(Imgproc_FitEllipse, accuracy) { CV_FitEllipseTest test; test.safe_run(); } TEST(Imgproc_FitLine, accuracy) { CV_FitLineTest test; test.safe_run(); } TEST(Imgproc_ContourMoments, accuracy) { CV_ContourMomentsTest test; test.safe_run(); } TEST(Imgproc_ContourPerimeterSlice, accuracy) { CV_PerimeterAreaSliceTest test; test.safe_run(); } +TEST(Imgproc_FitEllipse, small) { CV_FitEllipseSmallTest test; test.safe_run(); } /* End of file. */