Skip to content

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

Merged
asmorkalov merged 3 commits intoopencv:4.xfrom
s-trinh:feat/getClosestEllipsePoints_2
Jun 3, 2025
Merged

Add getClosestEllipsePoints() function to get the closest point on an ellipse#26299
asmorkalov merged 3 commits intoopencv:4.xfrom
s-trinh:feat/getClosestEllipsePoints_2

Conversation

@s-trinh
Copy link
Copy Markdown
Contributor

@s-trinh s-trinh commented Oct 12, 2024

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

image

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 s-trinh mentioned this pull request Oct 12, 2024
Closed
6 tasks
@asmorkalov asmorkalov added this to the 4.11.0 milestone Oct 14, 2024
@asmorkalov asmorkalov self-requested a review October 14, 2024 06:51
@s-trinh 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
Copy link
Copy Markdown
Contributor Author

s-trinh commented Oct 15, 2024

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 s-trinh force-pushed the feat/getClosestEllipsePoints_2 branch 4 times, most recently from c25415c to 35fb4ea Compare October 19, 2024 20:44
@asmorkalov asmorkalov modified the milestones: 4.11.0, 4.12.0 Dec 19, 2024
@s-trinh s-trinh mentioned this pull request Jan 3, 2025
Open
6 tasks
@s-trinh
Copy link
Copy Markdown
Contributor Author

s-trinh commented Feb 4, 2025 •
edited
Loading

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

@asmorkalov asmorkalov force-pushed the feat/getClosestEllipsePoints_2 branch from f817d9f to d4114dc Compare March 7, 2025 06:29
vpisarev
vpisarev reviewed May 30, 2025
View reviewed changes
@vpisarev
Copy link
Copy Markdown
Contributor

vpisarev commented May 30, 2025

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

sturkmen72
sturkmen72 reviewed May 30, 2025
View reviewed changes
@vpisarev vpisarev requested review from sturkmen72 and vpisarev June 3, 2025 13:53
vpisarev
vpisarev approved these changes Jun 3, 2025
View reviewed changes
Copy link
Copy Markdown
Contributor

@vpisarev vpisarev left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good! Thank you for the contribution!

@asmorkalov asmorkalov merged commit e258f25 into opencv:4.x Jun 3, 2025
26 of 28 checks passed
@asmorkalov asmorkalov mentioned this pull request Jun 11, 2025
Merged
@asmorkalov asmorkalov mentioned this pull request Jul 21, 2025
Merged
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

@vpisarev vpisarev vpisarev approved these changes

@asmorkalov asmorkalov Awaiting requested review from asmorkalov

@sturkmen72 sturkmen72 Awaiting requested review from sturkmen72

Assignees

@vpisarev vpisarev

Labels

category: imgproc feature test

Projects

None yet

Milestone

4.12.0

Development

Successfully merging this pull request may close these issues.

4 participants

@s-trinh @vpisarev @sturkmen72 @asmorkalov