Tìm vòng tròn/hình elip bên trong của đám mây điểm 2D

Question

Tôi có một loạt điểm 2D. Bạn có thể nhìn thấy chúng trên hình bên trái. Họ tạo thành một số loại nhẫn với một vài tai thỏ. Mục tiêu của tôi là tìm chu trình/hình elip bên trong lớn, bạn có thể thấy ở bên phải.

Những loại thuật toán hữu ích cho một tình huống như vậy.

Tôi đã thử một biến thể của thuật toán RANSAC (lấy 5 điểm ngẫu nhiên, tạo thành một hình elip, xác định điểm và lặp lại). Tôi thiết kế chức năng chấm điểm theo một cách, mà điểm bên trong hình elip nhận được rất nhiều điểm tiêu cực, và điểm bên ngoài, nhưng rất gần nhận được rất nhiều điểm tích cực. Nhưng kết quả không hề hứa hẹn. Các thuật toán tìm thấy chiếc nhẫn, nhưng tôi nhận được một số hình elip ngẫu nhiên của chiếc nhẫn, không phải là hình elip bên trong lớn tôi muốn.

Có chiến lược hay nào không?

Answer 1

Dưới đây là một số mã Java cho vòng tròn phù hợp nhất.

https://www.spaceroots.org/documents/circle/CircleFitter.java

// Copyright (c) 2005-2007, Luc Maisonobe 
// All rights reserved. 
// 
// Redistribution and use in source and binary forms, with 
// or without modification, are permitted provided that 
// the following conditions are met: 
// 
// Redistributions of source code must retain the 
// above copyright notice, this list of conditions and 
// the following disclaimer. 
// Redistributions in binary form must reproduce the 
// above copyright notice, this list of conditions and 
// the following disclaimer in the documentation 
// and/or other materials provided with the 
// distribution. 
// Neither the names of spaceroots.org, spaceroots.com 
// nor the names of their contributors may be used to 
// endorse or promote products derived from this 
// software without specific prior written permission. 
// 
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND 
// CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED 
// WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 
// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL 
// THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY 
// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF 
// USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 
// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER 
// IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE 
// USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
// POSSIBILITY OF SUCH DAMAGE. 

package org.spaceroots; 

import java.io.Reader; 
import java.io.BufferedReader; 
import java.io.InputStreamReader; 
import java.io.FileReader; 
import java.io.IOException; 
import java.util.ArrayList; 
import java.util.Iterator; 
import java.util.Locale; 
import java.text.DecimalFormat; 
import java.text.DecimalFormatSymbols; 
import java.awt.geom.Point2D; 

/** Class fitting a circle to a set of points. 
* <p>This class implements the fitting algorithms described in the 
* paper <a 
* href="http://www.spaceroots.org/documents/circle/circle-fitting.pdf"> 
* Finding the circle that best fits a set of points</a></p> 
* @author Luc Maisonobe 
*/ 
public class CircleFitter { 

    /** Test program entry point. 
    * @param args command line arguments 
    */ 
    public static void main(String[] args) { 
    try { 

     BufferedReader br = null; 
     switch (args.length) { 
     case 0 : 
     br = new BufferedReader((new InputStreamReader(System.in))); 
     break; 
     case 1: 
     br = new BufferedReader(new FileReader(args[0])); 
     break; 
     default : 
     System.err.println("usage: java CircleFitter [file]"); 
     System.exit(1); 
     } 

     // read the points, ignoring blank lines and comment lines 
     ArrayList list = new ArrayList(); 
     int l = 0; 
     for (String line = br.readLine(); line != null; line = br.readLine()) { 
     ++l; 
     line = line.trim(); 
     if ((line.length() > 0) && (! line.startsWith("#"))) { 
      // this is a data line, we expect two numerical fields 
      String[] fields = line.split("\\s+"); 
      if (fields.length != 2) { 
      throw new LocalException("syntax error at line " + l + ": " + line 
            + "(expected two fields, found" 
            + fields.length + ")"); 
      } 

      // parse the fields and add the point to the list 
      list.add(new Point2D.Double(Double.parseDouble(fields[0]), 
             Double.parseDouble(fields[1]))); 

     } 
     } 
     Point2D.Double[] points = 
     (Point2D.Double[]) list.toArray(new Point2D.Double[list.size()]); 

     DecimalFormat format = 
     new DecimalFormat("000.00000000", 
          new DecimalFormatSymbols(Locale.US)); 

     // fit a circle to the test points 
     CircleFitter fitter = new CircleFitter(); 
     fitter.initialize(points); 
     System.out.println("initial circle: " 
         + format.format(fitter.getCenter().x) 
         + " "  + format.format(fitter.getCenter().y) 
         + " "  + format.format(fitter.getRadius())); 

     // minimize the residuals 
     int iter = fitter.minimize(100, 0.1, 1.0e-12); 
     System.out.println("converged after " + iter + " iterations"); 
     System.out.println("final circle: " 
         + format.format(fitter.getCenter().x) 
         + " "  + format.format(fitter.getCenter().y) 
         + " "  + format.format(fitter.getRadius())); 

    } catch (IOException ioe) { 
     System.err.println(ioe.getMessage()); 
     System.exit(1); 
    } catch (LocalException le) { 
     System.err.println(le.getMessage()); 
     System.exit(1); 
    } 
    } 

    /** Build a new instance with a default current circle. 
    */ 
    public CircleFitter() { 
    center = new Point2D.Double(0.0, 0.0); 
    rHat = 1.0; 
    points = null; 
    } 


    /** Initialize an approximate circle based on all triplets. 
    * @param points circular ring sample points 
    * @exception LocalException if all points are aligned 
    */ 
    public void initialize(Point2D.Double[] points) 
    throws LocalException { 

    // store the points array 
    this.points = points; 

    // analyze all possible points triplets 
    center.x = 0.0; 
    center.y = 0.0; 
    int n = 0; 
    for (int i = 0; i < (points.length - 2); ++i) { 
     Point2D.Double p1 = (Point2D.Double) points[i]; 
     for (int j = i + 1; j < (points.length - 1); ++j) { 
     Point2D.Double p2 = (Point2D.Double) points[j]; 
     for (int k = j + 1; k < points.length; ++k) { 
      Point2D.Double p3 = (Point2D.Double) points[k]; 

      // compute the triangle circumcenter 
      Point2D.Double cc = circumcenter(p1, p2, p3); 
      if (cc != null) { 
      // the points are not aligned, we have a circumcenter 
      ++n; 
      center.x += cc.x; 
      center.y += cc.y; 
      } 
     } 
     } 
    } 

    if (n == 0) { 
     throw new LocalException("all points are aligned"); 
    } 

    // initialize using the circumcenters average 
    center.x /= n; 
    center.y /= n; 
    updateRadius(); 

    } 

    /** Update the circle radius. 
    */ 
    private void updateRadius() { 
    rHat = 0; 
    for (int i = 0; i < points.length; ++i) { 
     double dx = points[i].x - center.x; 
     double dy = points[i].y - center.y; 
     rHat += Math.sqrt(dx * dx + dy * dy); 
    } 
    rHat /= points.length; 
    } 

    /** Compute the circumcenter of three points. 
    * @param pI first point 
    * @param pJ second point 
    * @param pK third point 
    * @return circumcenter of pI, pJ and pK or null if the points are aligned 
    */ 
    private Point2D.Double circumcenter(Point2D.Double pI, 
             Point2D.Double pJ, 
             Point2D.Double pK) { 

    // some temporary variables 
    Point2D.Double dIJ = new Point2D.Double(pJ.x - pI.x, pJ.y - pI.y); 
    Point2D.Double dJK = new Point2D.Double(pK.x - pJ.x, pK.y - pJ.y); 
    Point2D.Double dKI = new Point2D.Double(pI.x - pK.x, pI.y - pK.y); 
    double sqI = pI.x * pI.x + pI.y * pI.y; 
    double sqJ = pJ.x * pJ.x + pJ.y * pJ.y; 
    double sqK = pK.x * pK.x + pK.y * pK.y; 

    // determinant of the linear system: 0 for aligned points 
    double det = dJK.x * dIJ.y - dIJ.x * dJK.y; 
    if (Math.abs(det) < 1.0e-10) { 
     // points are almost aligned, we cannot compute the circumcenter 
     return null; 
    } 

    // beware, there is a minus sign on Y coordinate! 
    return new Point2D.Double(
      (sqI * dJK.y + sqJ * dKI.y + sqK * dIJ.y)/(2 * det), 
      -(sqI * dJK.x + sqJ * dKI.x + sqK * dIJ.x)/(2 * det)); 

    } 

    /** Minimize the distance residuals between the points and the circle. 
    * <p>We use a non-linear conjugate gradient method with the Polak and 
    * Ribiere coefficient for the computation of the search direction. The 
    * inner minimization along the search direction is performed using a 
    * few Newton steps. It is worthless to spend too much time on this inner 
    * minimization, so the convergence threshold can be rather large.</p> 
    * @param maxIter maximal iterations number on the inner loop (cumulated 
    * across outer loop iterations) 
    * @param innerThreshold inner loop threshold, as a relative difference on 
    * the cost function value between the two last iterations 
    * @param outerThreshold outer loop threshold, as a relative difference on 
    * the cost function value between the two last iterations 
    * @return number of inner loop iterations performed (cumulated 
    * across outer loop iterations) 
    * @exception LocalException if we come accross a singularity or if 
    * we exceed the maximal number of iterations 
    */ 
    public int minimize(int iterMax, 
         double innerThreshold, double outerThreshold) 
    throws LocalException { 

    computeCost(); 
    if ((J < 1.0e-10) || (Math.sqrt(dJdx * dJdx + dJdy * dJdy) < 1.0e-10)) { 
     // we consider we are already at a local minimum 
     return 0; 
    } 

    double previousJ = J; 
    double previousU = 0.0, previousV = 0.0; 
    double previousDJdx = 0.0, previousDJdy = 0.0; 
    for (int iterations = 0; iterations < iterMax;) { 

     // search direction 
     double u = -dJdx; 
     double v = -dJdy; 
     if (iterations != 0) { 
     // Polak-Ribiere coefficient 
     double beta = 
      (dJdx * (dJdx - previousDJdx) + dJdy * (dJdy - previousDJdy)) 
     /(previousDJdx * previousDJdx + previousDJdy * previousDJdy); 
     u += beta * previousU; 
     v += beta * previousV; 
     } 
     previousDJdx = dJdx; 
     previousDJdy = dJdy; 
     previousU = u; 
     previousV = v; 

     // rough minimization along the search direction 
     double innerJ; 
     do { 
     innerJ = J; 
     double lambda = newtonStep(u, v); 
     center.x += lambda * u; 
     center.y += lambda * v; 
     updateRadius(); 
     computeCost(); 
     } while ((++iterations < iterMax) 
       && ((Math.abs(J - innerJ)/J) > innerThreshold)); 

     // global convergence test 
     if ((Math.abs(J - previousJ)/J) < outerThreshold) { 
     return iterations; 
     } 
     previousJ = J; 

    } 

    throw new LocalException("unable to converge after " 
          + iterMax + " iterations"); 

    } 

    /** Compute the cost function and its gradient. 
    * <p>The results are stored as instance attributes.</p> 
    */ 
    private void computeCost() throws LocalException { 
    J = 0; 
    dJdx = 0; 
    dJdy = 0; 
    for (int i = 0; i < points.length; ++i) { 
     double dx = points[i].x - center.x; 
     double dy = points[i].y - center.y; 
     double di = Math.sqrt(dx * dx + dy * dy); 
     if (di < 1.0e-10) { 
     throw new LocalException("cost singularity:" 
           + " point at the circle center"); 
     } 
     double dr = di - rHat; 
     double ratio = dr/di; 
     J += dr * (di + rHat); 
     dJdx += dx * ratio; 
     dJdy += dy * ratio; 
    } 
    dJdx *= 2.0; 
    dJdy *= 2.0; 
    } 

    /** Compute the length of the Newton step in the search direction. 
    * @param u abscissa of the search direction 
    * @param v ordinate of the search direction 
    * @return value of the step along the search direction 
    */ 
    private double newtonStep(double u, double v) { 

    // compute the first and second derivatives of the cost 
    // along the specified search direction 
    double sum1 = 0, sum2 = 0, sumFac = 0, sumFac2R = 0; 
    for (int i = 0; i < points.length; ++i) { 
     double dx  = center.x - points[i].x; 
     double dy  = center.y - points[i].y; 
     double di  = Math.sqrt(dx * dx + dy * dy); 
     double coeff1 = (dx * u + dy * v)/di; 
     double coeff2 = di - rHat; 
     sum1   += coeff1 * coeff2; 
     sum2   += coeff2/di; 
     sumFac  += coeff1; 
     sumFac2R  += coeff1 * coeff1/di; 
    } 

    // step length attempting to nullify the first derivative 
    return -sum1/((u * u + v * v) * sum2 
        - sumFac * sumFac/points.length 
        + rHat * sumFac2R); 

    } 

    /** Get the circle center. 
    * @return circle center 
    */ 
    public Point2D.Double getCenter() { 
    return center; 
    } 

    /** Get the circle radius. 
    * @return circle radius 
    */ 
    public double getRadius() { 
    return rHat; 
    } 

    /** Local exception class for algorithm errors. */ 
    public static class LocalException extends Exception { 
    /** Build a new instance with the supplied message. 
    * @param message error message 
    */ 
    public LocalException(String message) { 
     super(message); 
    } 
    } 

    /** Current circle center. */ 
    private Point2D.Double center; 

    /** Current circle radius. */ 
    private double rHat; 

    /** Circular ring sample points. */ 
    private Point2D.Double[] points; 

    /** Current cost function value. */ 
    private double J; 

    /** Current cost function gradient. */ 
    private double dJdx; 
    private double dJdy; 

} 

EDIT Một khi bạn có một vòng tròn phù hợp nhất bạn có thể thu nhỏ nó cho đến khi một số tỷ lệ nội đến các điểm ngoài được đáp ứng. Hoặc bạn có thể loại bỏ tất cả các điểm bên ngoài vòng tròn và chạy thuật toán vòng tròn phù hợp nhất một lần nữa, lặp lại quy trình này cho đến khi bạn có câu trả lời thỏa đáng.