From 4ebaf2dbb4ad1c1420a1c69858441540654e7ab0 Mon Sep 17 00:00:00 2001 From: gideok Kim Date: Tue, 10 Mar 2026 23:18:02 +0900 Subject: [PATCH] imgproc: fix fitEllipseDirect/AMS determinant threshold for near-circular data MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit fitEllipseDirect: replace det(M) threshold with eigenvector quality check - Add Ts ≈ 0 guard to avoid division by zero in Schur complement - Move eigenNonSymmetric inside perturbation loop - Validate eigenvector with 4ac-b² > 1e-6*||v||² to filter garbage from complex eigenvalues fitEllipseAMS: scale threshold by 1/n^5 to match det(M) magnitude Add Imgproc_FitEllipseDirect_NearCircular regression test --- modules/imgproc/src/shapedescr.cpp | 32 +++++++++++++------ modules/imgproc/test/test_fitellipse_ams.cpp | 20 ++++++++++++ .../imgproc/test/test_fitellipse_direct.cpp | 22 +++++++++++++ 3 files changed, 64 insertions(+), 10 deletions(-) diff --git a/modules/imgproc/src/shapedescr.cpp b/modules/imgproc/src/shapedescr.cpp index 6ddd855e32..4952a9722e 100644 --- a/modules/imgproc/src/shapedescr.cpp +++ b/modules/imgproc/src/shapedescr.cpp @@ -579,7 +579,9 @@ cv::RotatedRect cv::fitEllipseAMS( InputArray _points ) M(4,3)=DM(3,4); M(4,4)=DM(4,4); - if (fabs(cv::determinant(M)) > 1.0e-10) { + double npow = (double)n * (double)n; + npow = npow * npow * (double)n; // n^5 + if (fabs(cv::determinant(M)) > 1.0e-10 / npow) { break; } @@ -703,6 +705,8 @@ cv::RotatedRect cv::fitEllipseDirect( InputArray _points ) Matx pVec; double x0, y0, a, b, theta, Ts; + Mat eVal, eVec; + double cond[3]; double s = 0; for( i = 0; i < n; i++ ) @@ -763,6 +767,11 @@ cv::RotatedRect cv::fitEllipseDirect( InputArray _points ) Ts=(-(DM(3,5)*DM(4,4)*DM(5,3)) + DM(3,4)*DM(4,5)*DM(5,3) + DM(3,5)*DM(4,3)*DM(5,4) - \ DM(3,3)*DM(4,5)*DM(5,4) - DM(3,4)*DM(4,3)*DM(5,5) + DM(3,3)*DM(4,4)*DM(5,5)); + if (fabs(Ts) < DBL_EPSILON) { + eps = (float)(s/(n*2)*1e-2); + continue; + } + M(0,0) = (DM(2,0) + (DM(2,3)*TM(0,0) + DM(2,4)*TM(1,0) + DM(2,5)*TM(2,0))/Ts)/2.; M(0,1) = (DM(2,1) + (DM(2,3)*TM(0,1) + DM(2,4)*TM(1,1) + DM(2,5)*TM(2,1))/Ts)/2.; M(0,2) = (DM(2,2) + (DM(2,3)*TM(0,2) + DM(2,4)*TM(1,2) + DM(2,5)*TM(2,2))/Ts)/2.; @@ -773,18 +782,9 @@ cv::RotatedRect cv::fitEllipseDirect( InputArray _points ) M(2,1) = (DM(0,1) + (DM(0,3)*TM(0,1) + DM(0,4)*TM(1,1) + DM(0,5)*TM(2,1))/Ts)/2.; M(2,2) = (DM(0,2) + (DM(0,3)*TM(0,2) + DM(0,4)*TM(1,2) + DM(0,5)*TM(2,2))/Ts)/2.; - double det = cv::determinant(M); - if (fabs(det) > 1.0e-10) - break; - eps = (float)(s/(n*2)*1e-2); - } - - if( iter < 2 ) { - Mat eVal, eVec; eigenNonSymmetric(M, eVal, eVec); // Select the eigen vector {a,b,c} which satisfies 4ac-b^2 > 0 - double cond[3]; cond[0]=(4.0 * eVec.at(0,0) * eVec.at(0,2) - eVec.at(0,1) * eVec.at(0,1)); cond[1]=(4.0 * eVec.at(1,0) * eVec.at(1,2) - eVec.at(1,1) * eVec.at(1,1)); cond[2]=(4.0 * eVec.at(2,0) * eVec.at(2,2) - eVec.at(2,1) * eVec.at(2,1)); @@ -793,6 +793,18 @@ cv::RotatedRect cv::fitEllipseDirect( InputArray _points ) } else { i = (cond[0](i,0), v1 = eVec.at(i,1), v2 = eVec.at(i,2); + double vnorm2 = v0*v0 + v1*v1 + v2*v2; + if (cond[i] > 1e-6 * vnorm2) + break; + } + eps = (float)(s/(n*2)*1e-2); + } + + if( iter < 2 ) { double norm = std::sqrt(eVec.at(i,0)*eVec.at(i,0) + eVec.at(i,1)*eVec.at(i,1) + eVec.at(i,2)*eVec.at(i,2)); if (((eVec.at(i,0)<0.0 ? -1 : 1) * (eVec.at(i,1)<0.0 ? -1 : 1) * (eVec.at(i,2)<0.0 ? -1 : 1)) <= 0.0) { norm=-1.0*norm; diff --git a/modules/imgproc/test/test_fitellipse_ams.cpp b/modules/imgproc/test/test_fitellipse_ams.cpp index f2c9d1f793..b88660e1f3 100644 --- a/modules/imgproc/test/test_fitellipse_ams.cpp +++ b/modules/imgproc/test/test_fitellipse_ams.cpp @@ -337,6 +337,26 @@ TEST(Imgproc_FitEllipseAMS_Issue_7, accuracy) { EXPECT_TRUE(checkEllipse(ellipseAMSTest, ellipseAMSTrue, tol)); } +TEST(Imgproc_FitEllipseAMS_NearCircular, accuracy) +{ + std::vector points; + double cx = 27.0, cy = 27.0, a = 17.0, b = 16.5; + for (int i = 0; i < 360; i++) { + double theta = 2.0 * CV_PI * i / 360.0; + points.push_back(cv::Point2f( + (float)(cx + a * cos(theta)), + (float)(cy + b * sin(theta)))); + } + + cv::RotatedRect ams = cv::fitEllipseAMS(points); + + // AMS should produce a valid result close to ground truth + EXPECT_NEAR(ams.center.x, 27.0, 0.5); + EXPECT_NEAR(ams.center.y, 27.0, 0.5); + EXPECT_NEAR(std::max(ams.size.width, ams.size.height), 34.0, 1.0); + EXPECT_NEAR(std::min(ams.size.width, ams.size.height), 33.0, 1.0); +} + TEST(Imgproc_FitEllipseAMS_HorizontalLine, accuracy) { vector pts({{-300, 100}, {-200, 100}, {-100, 100}, {0, 100}, {100, 100}, {200, 100}, {300, 100}}); const RotatedRect el = fitEllipseAMS(pts); diff --git a/modules/imgproc/test/test_fitellipse_direct.cpp b/modules/imgproc/test/test_fitellipse_direct.cpp index e41c52764d..81d05a9280 100644 --- a/modules/imgproc/test/test_fitellipse_direct.cpp +++ b/modules/imgproc/test/test_fitellipse_direct.cpp @@ -337,6 +337,28 @@ TEST(Imgproc_FitEllipseDirect_Issue_7, accuracy) { EXPECT_TRUE(checkEllipse(ellipseDirectTest, ellipseDirectTrue, tol)); } +TEST(Imgproc_FitEllipseDirect_NearCircular, accuracy) +{ + // 360 points on a near-circular ellipse (a=17, b=16.5) + // This data previously triggered unnecessary fallback to fitEllipseNoDirect + std::vector points; + double cx = 27.0, cy = 27.0, a = 17.0, b = 16.5; + for (int i = 0; i < 360; i++) { + double theta = 2.0 * CV_PI * i / 360.0; + points.push_back(cv::Point2f( + (float)(cx + a * cos(theta)), + (float)(cy + b * sin(theta)))); + } + + cv::RotatedRect direct = cv::fitEllipseDirect(points); + + // Direct should produce a valid result close to ground truth + EXPECT_NEAR(direct.center.x, 27.0, 0.1); + EXPECT_NEAR(direct.center.y, 27.0, 0.1); + EXPECT_NEAR(std::max(direct.size.width, direct.size.height), 34.0, 0.5); + EXPECT_NEAR(std::min(direct.size.width, direct.size.height), 33.0, 0.5); +} + TEST(Imgproc_FitEllipseDirect_HorizontalLine, accuracy) { vector pts({{-300, 100}, {-200, 100}, {-100, 100}, {0, 100}, {100, 100}, {200, 100}, {300, 100}}); const RotatedRect el = fitEllipseDirect(pts);