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()

…int belonging to a considered ellipse equation.

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 3 times, most recently from ebd4219 to c25415c Compare October 16, 2024 19:40

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

35fb4ea

…int belonging to a considered ellipse equation.

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

Prefer to use https links.

f817d9f

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

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

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

s-trinh commented Oct 12, 2024

s-trinh commented Oct 15, 2024

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

Are you sure you want to change the base?

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

Conversation

s-trinh commented Oct 12, 2024

Pull Request Readiness Checklist

s-trinh commented Oct 15, 2024