Add getClosestEllipsePoints() function to get the closest point on an ellipse #26299

s-trinh · 2024-10-12T10:55:53Z

Sorry, for some reason the GitHub CI was not run in the previous PR #26237
I have opened this one to have at least results from the GitHub CI + pullrequest.opencv.org

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:

Demo code:

code

#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;
}

Pull Request Readiness Checklist

See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request

I agree to contribute to the project under Apache 2 License.
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
The PR is proposed to the proper branch
There is a reference to the original bug report and related work
There is accuracy test, performance test and test data in opencv_extra repository, if applicable
Patch to opencv_extra has the same branch name.
The feature is well documented and sample code can be built with the project CMake

s-trinh · 2024-10-15T21:54:52Z

I have added additional tests to check the code correctness when handling "minor/major axes swapping" cases.

Python code to validate ellipse points sampling:

code

#!/usr/bin/env python3
import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
def getEllipsePointsCV(ellipse):
    center, axes, angle = ellipse
    ellipse_points = cv.ellipse2Poly((int(center[0]), int(center[1])), (int(axes[0] // 2), int(axes[1] // 2)), int(angle), 0, 360, 1)
    return ellipse_points
# https://github.com/opencv/opencv/pull/26299
def main():
    (xc, yc) = (1442.97900390625, 662.1879272460938)
    (a, b) = (579.5570678710938/2, 730.834228515625/2)
    theta_deg = 20.190902709960938
    theta = np.deg2rad(theta_deg)
    ellipse_params = ((float(xc), float(yc)),(float(a*2), float(b*2)), float(theta_deg))
    ellipse_points = getEllipsePointsCV(ellipse_params)
    xs = []
    ys = []
    ts = np.linspace(0, 360, 60)
    for t in ts:
        # https://fr.wikipedia.org/wiki/Ellipse_(math%C3%A9matiques)#%C3%89quation_param%C3%A9trique
        ax = a * np.cos(theta)
        ay = a * np.sin(theta)
        bx = -b * np.sin(theta)
        by = b * np.cos(theta)
        cos_t = np.cos(np.deg2rad(t))
        sin_t = np.sin(np.deg2rad(t))
        x = xc + ax*cos_t + bx*sin_t
        y = yc + ay*cos_t + by*sin_t
        xs.append(x)
        ys.append(y)
    xs_ = np.expand_dims(np.array(xs), axis=1)
    ys_ = np.expand_dims(np.array(ys), axis=1)
    input_pts = np.hstack((xs_, ys_)).astype(np.float32)
    closest_pts = cv.getClosestEllipsePoints(ellipse_params, input_pts).reshape((input_pts.shape[0], 2))
    error_pts = closest_pts - input_pts
    print(f"input_pts={input_pts}")
    print(f"closest_pts={closest_pts}")
    print(f"error_pts={error_pts}")
    mse = (error_pts**2).mean(axis=1).mean()
    print(f"mse={mse}")
    _, axs = plt.subplots(figsize=(8, 8))
    axs.scatter(xs, ys, c='blue', label='Points')
    axs.scatter(closest_pts[:,0], closest_pts[:,1], c='red', label='Closest points')
    ellipse_polygon = plt.Polygon(ellipse_points, fill=None, edgecolor='green', \
                                      label='Ellipse from points')
    axs.add_patch(ellipse_polygon)
    axs.set_aspect('equal', adjustable='box')
    axs.legend()
    plt.show()
if __name__ == '__main__':
    main()

s-trinh · 2025-02-04T20:56:01Z

Additional information if it can help someone.

In skimage.measure.EllipseModel, there is the residuals() function to compute the ellipse fitting residuals.
From a very quick glance to the source code, the idea seems to be:

use the parametric form to describe the ellipse
from an initial guess t value
optimize the parametric t (see optimize.leastsq(fun, t0[i], args=(xi, yi))) to compute the value yielding the lowest distance

…int belonging to a considered ellipse equation.

modules/imgproc/src/shapedescr.cpp

vpisarev · 2025-05-30T05:37:30Z

@s-trinh, sorry for late reply. The functionality is useful. Except for the very slow and unnecessary Mat allocation, PR looks good.

modules/imgproc/test/test_fitellipse.cpp

… typo.

vpisarev

Looks good! Thank you for the contribution!

s-trinh mentioned this pull request Oct 12, 2024

Add getClosestEllipsePoints() function to get the closest point on an ellipse #26237

Closed

6 tasks

asmorkalov added test feature category: imgproc labels Oct 14, 2024

asmorkalov added this to the 4.11.0 milestone Oct 14, 2024

asmorkalov self-requested a review October 14, 2024 06:51

s-trinh force-pushed the feat/getClosestEllipsePoints_2 branch 2 times, most recently from 291c1cf to a5f19d2 Compare October 15, 2024 21:49

s-trinh force-pushed the feat/getClosestEllipsePoints_2 branch 4 times, most recently from c25415c to 35fb4ea Compare October 19, 2024 20:44

asmorkalov modified the milestones: 4.11.0, 4.12.0 Dec 19, 2024

s-trinh mentioned this pull request Jan 3, 2025

Fix and add tests for ellipse fitting #26694

Open

6 tasks

s-trinh added 2 commits March 7, 2025 09:29

Add function getClosestEllipsePoints() in order to get the closest po…

9e6083d

…int belonging to a considered ellipse equation.

Prefer to use https links.

d4114dc

asmorkalov force-pushed the feat/getClosestEllipsePoints_2 branch from f817d9f to d4114dc Compare March 7, 2025 06:29

vpisarev reviewed May 30, 2025

View reviewed changes

modules/imgproc/src/shapedescr.cpp Outdated Show resolved Hide resolved

sturkmen72 reviewed May 30, 2025

View reviewed changes

modules/imgproc/test/test_fitellipse.cpp Outdated Show resolved Hide resolved

Adress review comments: use Matx instead of Mat in the loop code, fix…

a26fa7b

… typo.

vpisarev requested review from vpisarev and sturkmen72 June 3, 2025 13:53

vpisarev approved these changes Jun 3, 2025

View reviewed changes

asmorkalov assigned vpisarev Jun 3, 2025

asmorkalov merged commit e258f25 into opencv:4.x Jun 3, 2025
26 of 28 checks passed

asmorkalov mentioned this pull request Jun 11, 2025

5.x merge 4.x #27435

Merged

asmorkalov mentioned this pull request Jul 21, 2025

5.x merge 4.x #27566

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Uh oh!

Add getClosestEllipsePoints() function to get the closest point on an ellipse #26299

Add getClosestEllipsePoints() function to get the closest point on an ellipse #26299

Uh oh!

s-trinh commented Oct 12, 2024

Uh oh!

s-trinh commented Oct 15, 2024

Uh oh!

s-trinh commented Feb 4, 2025 •

edited

Loading

Uh oh!

Uh oh!

vpisarev commented May 30, 2025

Uh oh!

Uh oh!

vpisarev left a comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Add getClosestEllipsePoints() function to get the closest point on an ellipse #26299

Add getClosestEllipsePoints() function to get the closest point on an ellipse #26299

Uh oh!

Conversation

s-trinh commented Oct 12, 2024

Pull Request Readiness Checklist

Uh oh!

s-trinh commented Oct 15, 2024

Uh oh!

s-trinh commented Feb 4, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

vpisarev commented May 30, 2025

Uh oh!

Uh oh!

vpisarev left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

s-trinh commented Feb 4, 2025 •

edited

Loading