Merge pull request #26299 from s-trinh:feat/getClosestEllipsePoints_2

s-trinh · web-flow · commit e258f2595e84 · 2025-06-03T17:12:16.000+03:00
Add getClosestEllipsePoints() function to get the closest point on an ellipse #26299 Following #26078, I was thinking that a function to get for a considered 2d point the corresponding closest point (or maybe directly the distance?) on an ellipse could be useful. This would allow computing the fitting error with `fitEllipse()` for instance. Code is based from: - https://stackoverflow.com/questions/22959698/distance-from-given-point-to-given-ellipse/46007540#46007540 - https://blog.chatfield.io/simple-method-for-distance-to-ellipse/ - https://github.com/0xfaded/ellipse_demo --- Demo code: <details> <summary>code</summary> ```cpp #include <iostream> #include <opencv2/opencv.hpp> namespace { void scaleApplyColormap(const cv::Mat &img_float, cv::Mat &img) { cv::Mat img_scale = cv::Mat::zeros(img_float.size(), CV_8UC3); double min_val = 0, max_val = 0; cv::minMaxLoc(img_float, &min_val, &max_val); std::cout << "min_val=" << min_val << " ; max_val=" << max_val << std::endl; if (max_val - min_val > 1e-2) { float a = 255 / (max_val - min_val); float b = -a * min_val; cv::convertScaleAbs(img_float, img_scale, a, b); cv::applyColorMap(img_scale, img, cv::COLORMAP_TURBO); } else { std::cerr << "max_val - min_val <= 1e-2" << std::endl; } } cv::Mat drawEllipseDistanceMap(const cv::RotatedRect &ellipse_params) { float bb_rect_w = ellipse_params.center.x + ellipse_params.size.width; float bb_rect_h = ellipse_params.center.y + ellipse_params.size.height; std::vector<cv::Point2f> points_list; points_list.resize(1); cv::Mat pointsf; cv::Mat closest_pts; cv::Mat dist_map = cv::Mat::zeros(bb_rect_h*1.5, bb_rect_w*1.5, CV_32F); for (int i = 0; i < dist_map.rows; i++) { for (int j = 0; j < dist_map.cols; j++) { points_list[0].x = j; points_list[0].y = i; cv::Mat(points_list).convertTo(pointsf, CV_32F); cv::getClosestEllipsePoints(ellipse_params, pointsf, closest_pts); dist_map.at<float>(i, j) = std::hypot(closest_pts.at<cv::Point2f>(0).x-j, closest_pts.at<cv::Point2f>(0).y-i); } } cv::Mat dist_map_8u; scaleApplyColormap(dist_map, dist_map_8u); return dist_map_8u; } } int main() { std::vector<cv::Point2f> points_list; // [1434, 308], [1434, 309], [1433, 310], [1427, 310], [1427, 312], [1426, 313], [1422, 313], [1422, 314], points_list.push_back(cv::Point2f(1434, 308)); points_list.push_back(cv::Point2f(1434, 309)); points_list.push_back(cv::Point2f(1433, 310)); points_list.push_back(cv::Point2f(1427, 310)); points_list.push_back(cv::Point2f(1427, 312)); points_list.push_back(cv::Point2f(1426, 313)); points_list.push_back(cv::Point2f(1422, 313)); points_list.push_back(cv::Point2f(1422, 314)); // [1421, 315], [1415, 315], [1415, 316], [1414, 317], [1408, 317], [1408, 319], [1407, 320], [1403, 320], points_list.push_back(cv::Point2f(1421, 315)); points_list.push_back(cv::Point2f(1415, 315)); points_list.push_back(cv::Point2f(1415, 316)); points_list.push_back(cv::Point2f(1414, 317)); points_list.push_back(cv::Point2f(1408, 317)); points_list.push_back(cv::Point2f(1408, 319)); points_list.push_back(cv::Point2f(1407, 320)); points_list.push_back(cv::Point2f(1403, 320)); // [1403, 321], [1402, 322], [1396, 322], [1396, 323], [1395, 324], [1389, 324], [1389, 326], [1388, 327], points_list.push_back(cv::Point2f(1403, 321)); points_list.push_back(cv::Point2f(1402, 322)); points_list.push_back(cv::Point2f(1396, 322)); points_list.push_back(cv::Point2f(1396, 323)); points_list.push_back(cv::Point2f(1395, 324)); points_list.push_back(cv::Point2f(1389, 324)); points_list.push_back(cv::Point2f(1389, 326)); points_list.push_back(cv::Point2f(1388, 327)); // [1382, 327], [1382, 328], [1381, 329], [1376, 329], [1376, 330], [1375, 331], [1369, 331], [1369, 333], points_list.push_back(cv::Point2f(1382, 327)); points_list.push_back(cv::Point2f(1382, 328)); points_list.push_back(cv::Point2f(1381, 329)); points_list.push_back(cv::Point2f(1376, 329)); points_list.push_back(cv::Point2f(1376, 330)); points_list.push_back(cv::Point2f(1375, 331)); points_list.push_back(cv::Point2f(1369, 331)); points_list.push_back(cv::Point2f(1369, 333)); // [1368, 334], [1362, 334], [1362, 335], [1361, 336], [1359, 336], [1359, 1016], [1365, 1016], [1366, 1017], points_list.push_back(cv::Point2f(1368, 334)); points_list.push_back(cv::Point2f(1362, 334)); points_list.push_back(cv::Point2f(1362, 335)); points_list.push_back(cv::Point2f(1361, 336)); points_list.push_back(cv::Point2f(1359, 336)); points_list.push_back(cv::Point2f(1359, 1016)); points_list.push_back(cv::Point2f(1365, 1016)); points_list.push_back(cv::Point2f(1366, 1017)); // [1366, 1019], [1430, 1019], [1430, 1017], [1431, 1016], [1440, 1016], [1440, 308] points_list.push_back(cv::Point2f(1366, 1019)); points_list.push_back(cv::Point2f(1430, 1019)); points_list.push_back(cv::Point2f(1430, 1017)); points_list.push_back(cv::Point2f(1431, 1016)); points_list.push_back(cv::Point2f(1440, 1016)); points_list.push_back(cv::Point2f(1440, 308)); cv::Mat pointsf; cv::Mat(points_list).convertTo(pointsf, CV_32F); cv::RotatedRect ellipse_params = cv::fitEllipseAMS(pointsf); std::cout << "ellipse_params, center=" << ellipse_params.center << " ; size=" << ellipse_params.size << " ; angle=" << ellipse_params.angle << std::endl; cv::TickMeter tm; tm.start(); cv::Mat dist_map_8u = drawEllipseDistanceMap(ellipse_params); tm.stop(); std::cout << "Elapsed time: " << tm.getAvgTimeSec() << " sec" << std::endl; cv::Point center(ellipse_params.center.x, ellipse_params.center.y); cv::Point axis(ellipse_params.size.width/2, ellipse_params.size.height/2); std::vector<cv::Point> ellipse_pts_list; cv::ellipse2Poly(center, axis, ellipse_params.angle, 0, 360, 1, ellipse_pts_list); cv::polylines(dist_map_8u, ellipse_pts_list, false, cv::Scalar(0, 0, 0), 3); // Points to be fitted cv::Mat closest_pts; cv::getClosestEllipsePoints(ellipse_params, pointsf, closest_pts); for (int i = 0; i < closest_pts.rows; i++) { cv::Point pt; pt.x = closest_pts.at<cv::Point2f>(i).x; pt.y = closest_pts.at<cv::Point2f>(i).y; cv::circle(dist_map_8u, pt, 8, cv::Scalar(0, 0, 255), 2); } cv::imwrite("dist_map_8u.png", dist_map_8u); return EXIT_SUCCESS; } ``` </details> ![image](https://github.com/user-attachments/assets/3345cc86-ba83-44f9-ac78-74058a33a7dc) --- ### Pull Request Readiness Checklist See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request - [x] I agree to contribute to the project under Apache 2 License. - [x] To the best of my knowledge, the proposed patch is not based on a code under GPL or another license that is incompatible with OpenCV - [x] The PR is proposed to the proper branch - [x] There is a reference to the original bug report and related work - [x] There is accuracy test, performance test and test data in opencv_extra repository, if applicable Patch to opencv_extra has the same branch name. - [x] The feature is well documented and sample code can be built with the project CMake
diff --git a/doc/opencv.bib b/doc/opencv.bib
@@ -301,6 +301,13 @@ @article{ChambolleEtAl
   publisher = {Walter de Gruyter},
   url = {https://hal.science/hal-00437581v1}
 }
+@misc{Chatfield2017,
+  author = {Chatfield, Carl},
+  title = {A Simple Method for Distance to Ellipse},
+  year = {2017},
+  publisher = {GitHub},
+  howpublished = {\url{https://blog.chatfield.io/simple-method-for-distance-to-ellipse/}},
+}
 @article{Chaumette06,
   author = {Chaumette, Fran{\c c}ois and Hutchinson, S.},
   title = {{Visual servo control, Part I: Basic approaches}},
diff --git a/modules/imgproc/include/opencv2/imgproc.hpp b/modules/imgproc/include/opencv2/imgproc.hpp
@@ -118,7 +118,7 @@ This module offers a comprehensive suite of image processing functions, enabling
     coordinates needs to be retrieved. In the simplest case, the coordinates can be just rounded to the
     nearest integer coordinates and the corresponding pixel can be used. This is called a
     nearest-neighbor interpolation. However, a better result can be achieved by using more
-    sophisticated [interpolation methods](https://en.wikipedia.org/wiki/Multivariate_interpolation) ,
+    sophisticated [interpolation methods](https://en.wikipedia.org/wiki/Multivariate_interpolation) ,
     where a polynomial function is fit into some neighborhood of the computed pixel \f$(f_x(x,y),
     f_y(x,y))\f$, and then the value of the polynomial at \f$(f_x(x,y), f_y(x,y))\f$ is taken as the
     interpolated pixel value. In OpenCV, you can choose between several interpolation methods. See
@@ -1467,7 +1467,7 @@ CV_EXPORTS_W void getDerivKernels( OutputArray kx, OutputArray ky,
 /** @brief Returns Gabor filter coefficients.
 
 For more details about gabor filter equations and parameters, see: [Gabor
-Filter](https://en.wikipedia.org/wiki/Gabor_filter).
+Filter](https://en.wikipedia.org/wiki/Gabor_filter).
 
 @param ksize Size of the filter returned.
 @param sigma Standard deviation of the gaussian envelope.
@@ -1549,7 +1549,7 @@ CV_EXPORTS_W void GaussianBlur( InputArray src, OutputArray dst, Size ksize,
 /** @brief Applies the bilateral filter to an image.
 
 The function applies bilateral filtering to the input image, as described in
-https://www.dai.ed.ac.uk/CVonline/LOCAL_COPIES/MANDUCHI1/Bilateral_Filtering.html
+https://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/MANDUCHI1/Bilateral_Filtering.html
 bilateralFilter can reduce unwanted noise very well while keeping edges fairly sharp. However, it is
 very slow compared to most filters.
 
@@ -1659,7 +1659,7 @@ stackBlur can generate similar results as Gaussian blur, and the time consumptio
 It creates a kind of moving stack of colors whilst scanning through the image. Thereby it just has to add one new block of color to the right side
 of the stack and remove the leftmost color. The remaining colors on the topmost layer of the stack are either added on or reduced by one,
 depending on if they are on the right or on the left side of the stack. The only supported borderType is BORDER_REPLICATE.
-Original paper was proposed by Mario Klingemann, which can be found https://underdestruction.com/2004/02/25/stackblur-2004.
+Original paper was proposed by Mario Klingemann, which can be found https://underdestruction.com/2004/02/25/stackblur-2004.
 
 @param src input image. The number of channels can be arbitrary, but the depth should be one of
 CV_8U, CV_16U, CV_16S or CV_32F.
@@ -1874,7 +1874,7 @@ Check @ref tutorial_canny_detector "the corresponding tutorial" for more details
 The function finds edges in the input image and marks them in the output map edges using the
 Canny algorithm. The smallest value between threshold1 and threshold2 is used for edge linking. The
 largest value is used to find initial segments of strong edges. See
-<https://en.wikipedia.org/wiki/Canny_edge_detector>
+<https://en.wikipedia.org/wiki/Canny_edge_detector>
 
 @param image 8-bit input image.
 @param edges output edge map; single channels 8-bit image, which has the same size as image .
@@ -2143,7 +2143,7 @@ An example using the Hough line detector
 /** @brief Finds lines in a binary image using the standard Hough transform.
 
 The function implements the standard or standard multi-scale Hough transform algorithm for line
-detection. See <https://homepages.inf.ed.ac.uk/rbf/HIPR2/hough.htm> for a good explanation of Hough
+detection. See <https://homepages.inf.ed.ac.uk/rbf/HIPR2/hough.htm> for a good explanation of Hough
 transform.
 
 @param image 8-bit, single-channel binary source image. The image may be modified by the function.
@@ -2986,13 +2986,13 @@ CV_EXPORTS_W void accumulateWeighted( InputArray src, InputOutputArray dst,
 
 The operation takes advantage of the Fourier shift theorem for detecting the translational shift in
 the frequency domain. It can be used for fast image registration as well as motion estimation. For
-more information please see <https://en.wikipedia.org/wiki/Phase_correlation>
+more information please see <https://en.wikipedia.org/wiki/Phase_correlation>
 
 Calculates the cross-power spectrum of two supplied source arrays. The arrays are padded if needed
 with getOptimalDFTSize.
 
 The function performs the following equations:
-- First it applies a Hanning window (see <https://en.wikipedia.org/wiki/Hann_function>) to each
+- First it applies a Hanning window (see <https://en.wikipedia.org/wiki/Hann_function>) to each
 image to remove possible edge effects. This window is cached until the array size changes to speed
 up processing time.
 - Next it computes the forward DFTs of each source array:
@@ -3022,7 +3022,7 @@ CV_EXPORTS_W Point2d phaseCorrelate(InputArray src1, InputArray src2,
 
 /** @brief This function computes a Hanning window coefficients in two dimensions.
 
-See (https://en.wikipedia.org/wiki/Hann_function) and (https://en.wikipedia.org/wiki/Window_function)
+See (https://en.wikipedia.org/wiki/Hann_function) and (https://en.wikipedia.org/wiki/Window_function)
 for more information.
 
 An example is shown below:
@@ -3507,7 +3507,7 @@ An example using the GrabCut algorithm
 
 /** @brief Runs the GrabCut algorithm.
 
-The function implements the [GrabCut image segmentation algorithm](https://en.wikipedia.org/wiki/GrabCut).
+The function implements the [GrabCut image segmentation algorithm](https://en.wikipedia.org/wiki/GrabCut).
 
 @param img Input 8-bit 3-channel image.
 @param mask Input/output 8-bit single-channel mask. The mask is initialized by the function when
@@ -3834,7 +3834,7 @@ CV_EXPORTS_W Moments moments( InputArray array, bool binaryImage = false );
 /** @brief Calculates seven Hu invariants.
 
 The function calculates seven Hu invariants (introduced in @cite Hu62; see also
-<https://en.wikipedia.org/wiki/Image_moment>) defined as:
+<https://en.wikipedia.org/wiki/Image_moment>) defined as:
 
 \f[\begin{array}{l} hu[0]= \eta _{20}+ \eta _{02} \\ hu[1]=( \eta _{20}- \eta _{02})^{2}+4 \eta _{11}^{2} \\ hu[2]=( \eta _{30}-3 \eta _{12})^{2}+ (3 \eta _{21}- \eta _{03})^{2} \\ hu[3]=( \eta _{30}+ \eta _{12})^{2}+ ( \eta _{21}+ \eta _{03})^{2} \\ hu[4]=( \eta _{30}-3 \eta _{12})( \eta _{30}+ \eta _{12})[( \eta _{30}+ \eta _{12})^{2}-3( \eta _{21}+ \eta _{03})^{2}]+(3 \eta _{21}- \eta _{03})( \eta _{21}+ \eta _{03})[3( \eta _{30}+ \eta _{12})^{2}-( \eta _{21}+ \eta _{03})^{2}] \\ hu[5]=( \eta _{20}- \eta _{02})[( \eta _{30}+ \eta _{12})^{2}- ( \eta _{21}+ \eta _{03})^{2}]+4 \eta _{11}( \eta _{30}+ \eta _{12})( \eta _{21}+ \eta _{03}) \\ hu[6]=(3 \eta _{21}- \eta _{03})( \eta _{21}+ \eta _{03})[3( \eta _{30}+ \eta _{12})^{2}-( \eta _{21}+ \eta _{03})^{2}]-( \eta _{30}-3 \eta _{12})( \eta _{21}+ \eta _{03})[3( \eta _{30}+ \eta _{12})^{2}-( \eta _{21}+ \eta _{03})^{2}] \\ \end{array}\f]
 
@@ -4072,7 +4072,7 @@ CV_EXPORTS_W void findContoursLinkRuns(InputArray image, OutputArrayOfArrays con
 
 The function cv::approxPolyDP approximates a curve or a polygon with another curve/polygon with less
 vertices so that the distance between them is less or equal to the specified precision. It uses the
-Douglas-Peucker algorithm <https://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm>
+Douglas-Peucker algorithm <https://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm>
 
 @param curve Input vector of a 2D point stored in std::vector or Mat
 @param approxCurve Result of the approximation. The type should match the type of the input curve.
@@ -4326,6 +4326,9 @@ ellipse/rotatedRect data contains negative indices, due to the data points being
 border of the containing Mat element.
 
 @param points Input 2D point set, stored in std::vector\<\> or Mat
+
+@note Input point types are @ref Point2i or @ref Point2f and at least 5 points are required.
+@note @ref getClosestEllipsePoints function can be used to compute the ellipse fitting error.
  */
 CV_EXPORTS_W RotatedRect fitEllipse( InputArray points );
 
@@ -4363,6 +4366,9 @@ CV_EXPORTS_W RotatedRect fitEllipse( InputArray points );
  \f}
 
  @param points Input 2D point set, stored in std::vector\<\> or Mat
+
+ @note Input point types are @ref Point2i or @ref Point2f and at least 5 points are required.
+ @note @ref getClosestEllipsePoints function can be used to compute the ellipse fitting error.
  */
 CV_EXPORTS_W RotatedRect fitEllipseAMS( InputArray points );
 
@@ -4408,9 +4414,26 @@ CV_EXPORTS_W RotatedRect fitEllipseAMS( InputArray points );
  The scaling factor guarantees that  \f$A^T C A =1\f$.
 
  @param points Input 2D point set, stored in std::vector\<\> or Mat
+
+ @note Input point types are @ref Point2i or @ref Point2f and at least 5 points are required.
+ @note @ref getClosestEllipsePoints function can be used to compute the ellipse fitting error.
  */
 CV_EXPORTS_W RotatedRect fitEllipseDirect( InputArray points );
 
+/** @brief Compute for each 2d point the nearest 2d point located on a given ellipse.
+
+ The function computes the nearest 2d location on a given ellipse for a vector of 2d points and is based on @cite Chatfield2017 code.
+ This function can be used to compute for instance the ellipse fitting error.
+
+ @param ellipse_params Ellipse parameters
+ @param points Input 2d points
+ @param closest_pts For each 2d point, their corresponding closest 2d point located on a given ellipse
+
+ @note Input point types are @ref Point2i or @ref Point2f
+ @see fitEllipse, fitEllipseAMS, fitEllipseDirect
+ */
+CV_EXPORTS_W void getClosestEllipsePoints( const RotatedRect& ellipse_params, InputArray points, OutputArray closest_pts );
+
 /** @brief Fits a line to a 2D or 3D point set.
 
 The function fitLine fits a line to a 2D or 3D point set by minimizing \f$\sum_i \rho(r_i)\f$ where
@@ -4429,7 +4452,7 @@ of the following:
 - DIST_HUBER
 \f[\rho (r) =  \fork{r^2/2}{if \(r < C\)}{C \cdot (r-C/2)}{otherwise} \quad \text{where} \quad C=1.345\f]
 
-The algorithm is based on the M-estimator ( <https://en.wikipedia.org/wiki/M-estimator> ) technique
+The algorithm is based on the M-estimator ( <https://en.wikipedia.org/wiki/M-estimator> ) technique
 that iteratively fits the line using the weighted least-squares algorithm. After each iteration the
 weights \f$w_i\f$ are adjusted to be inversely proportional to \f$\rho(r_i)\f$ .
 
diff --git a/modules/imgproc/src/shapedescr.cpp b/modules/imgproc/src/shapedescr.cpp
@@ -877,6 +877,121 @@ cv::RotatedRect cv::fitEllipseDirect( InputArray _points )
     return box;
 }
 
+namespace cv
+{
+// @misc{Chatfield2017,
+//   author = {Chatfield, Carl},
+//   title = {A Simple Method for Distance to Ellipse},
+//   year = {2017},
+//   publisher = {GitHub},
+//   howpublished = {\url{https://blog.chatfield.io/simple-method-for-distance-to-ellipse/}},
+// }
+// https://github.com/0xfaded/ellipse_demo/blob/master/ellipse_trig_free.py
+static void solveFast(float semi_major, float semi_minor, const cv::Point2f& pt, cv::Point2f& closest_pt)
+{
+    float px = std::abs(pt.x);
+    float py = std::abs(pt.y);
+
+    float tx = 0.707f;
+    float ty = 0.707f;
+
+    float a = semi_major;
+    float b = semi_minor;
+
+    for (int iter = 0; iter < 3; iter++)
+    {
+        float x = a * tx;
+        float y = b * ty;
+
+        float ex = (a*a - b*b) * tx*tx*tx / a;
+        float ey = (b*b - a*a) * ty*ty*ty / b;
+
+        float rx = x - ex;
+        float ry = y - ey;
+
+        float qx = px - ex;
+        float qy = py - ey;
+
+        float r = std::hypotf(rx, ry);
+        float q = std::hypotf(qx, qy);
+
+        tx = std::min(1.0f, std::max(0.0f, (qx * r / q + ex) / a));
+        ty = std::min(1.0f, std::max(0.0f, (qy * r / q + ey) / b));
+        float t = std::hypotf(tx, ty);
+        tx /= t;
+        ty /= t;
+    }
+
+    closest_pt.x = std::copysign(a * tx, pt.x);
+    closest_pt.y = std::copysign(b * ty, pt.y);
+}
+} // namespace cv
+
+void cv::getClosestEllipsePoints( const RotatedRect& ellipse_params, InputArray _points, OutputArray closest_pts )
+{
+    CV_INSTRUMENT_REGION();
+
+    Mat points = _points.getMat();
+    int n = points.checkVector(2);
+    int depth = points.depth();
+    CV_Assert(depth == CV_32F || depth == CV_32S);
+    CV_Assert(n > 0);
+
+    bool is_float = (depth == CV_32F);
+    const Point* ptsi = points.ptr<Point>();
+    const Point2f* ptsf = points.ptr<Point2f>();
+
+    float semi_major = ellipse_params.size.width / 2.0f;
+    float semi_minor = ellipse_params.size.height / 2.0f;
+    float angle_deg = ellipse_params.angle;
+    if (semi_major < semi_minor)
+    {
+        std::swap(semi_major, semi_minor);
+        angle_deg += 90;
+    }
+
+    Matx23f align_T_ori_f32;
+    float theta_rad = static_cast<float>(angle_deg * M_PI / 180);
+    float co = std::cos(theta_rad);
+    float si = std::sin(theta_rad);
+    float shift_x = ellipse_params.center.x;
+    float shift_y = ellipse_params.center.y;
+
+    align_T_ori_f32(0,0) = co;
+    align_T_ori_f32(0,1) = si;
+    align_T_ori_f32(0,2) = -co*shift_x - si*shift_y;
+    align_T_ori_f32(1,0) = -si;
+    align_T_ori_f32(1,1) = co;
+    align_T_ori_f32(1,2) = si*shift_x - co*shift_y;
+
+    Matx23f ori_T_align_f32;
+    ori_T_align_f32(0,0) = co;
+    ori_T_align_f32(0,1) = -si;
+    ori_T_align_f32(0,2) = shift_x;
+    ori_T_align_f32(1,0) = si;
+    ori_T_align_f32(1,1) = co;
+    ori_T_align_f32(1,2) = shift_y;
+
+    std::vector<Point2f> closest_pts_list;
+    closest_pts_list.reserve(n);
+    for (int i = 0; i < n; i++)
+    {
+        Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
+        Matx31f pmat(p.x, p.y, 1);
+
+        Matx21f X_align = align_T_ori_f32 * pmat;
+        Point2f closest_pt;
+        solveFast(semi_major, semi_minor, Point2f(X_align(0,0), X_align(1,0)), closest_pt);
+
+        pmat(0,0) = closest_pt.x;
+        pmat(1,0) = closest_pt.y;
+        Matx21f closest_pt_ori = ori_T_align_f32 * pmat;
+        closest_pts_list.push_back(Point2f(closest_pt_ori(0,0), closest_pt_ori(1,0)));
+    }
+
+    cv::Mat(closest_pts_list).convertTo(closest_pts, CV_32F);
+}
+
 ////////////////////////////////////////////// C API ///////////////////////////////////////////
 
 CV_IMPL int
diff --git a/modules/imgproc/test/test_fitellipse.cpp b/modules/imgproc/test/test_fitellipse.cpp
