}
-/* 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<CvMat> 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 )
{
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;
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++ )
{
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
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]);
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 )
box.angle += 360;
if( box.angle > 360 )
box.angle -= 360;
-#endif
return box;
}
_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 );
cvShowImage( "test", img );
cvReleaseImage( &img );
cvWaitKey(0);
+ }
#endif
if( code < 0 )
}
+class CV_FitEllipseSmallTest : public cvtest::BaseTest
+{
+public:
+ CV_FitEllipseSmallTest() {}
+ ~CV_FitEllipseSmallTest() {}
+protected:
+ void run(int)
+ {
+ Size sz(50, 50);
+ vector<vector<Point> > c;
+ c.push_back(vector<Point>());
+ 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 *
\****************************************************************************************/
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. */
projects / profile / ivi / opencv.git / commitdiff