我有一堆2D点，你可以在左边的图片上看到它们。它们形成了一些带有一些兔耳朵的环状物。我的目标是找到右侧所示的大内部循环/椭圆。
我尝试了一种变化的RANSAC算法（取5个随机点，形成椭圆，确定一个分数并重复）。我设计了评分函数，使椭圆内的点得到很多负分，而椭圆外但非常接近的点得到很多正分。但结果并不令人满意。该算法可以找到环形，但我得到的是一些随机的环形椭圆，而不是我想要的大内部椭圆。
有没有更好的策略？
寻找二维点云的内部圆/椭圆
8
我尝试了一种变化的RANSAC算法（取5个随机点，形成椭圆，确定一个分数并重复）。我设计了评分函数，使椭圆内的点得到很多负分，而椭圆外但非常接近的点得到很多正分。但结果并不令人满意。该算法可以找到环形，但我得到的是一些随机的环形椭圆，而不是我想要的大内部椭圆。
有没有更好的策略？
- Jakube
6
1个回答
2
这是一些Java代码，用于实现“最佳拟合圆”。
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;
}
编辑：一旦你得到了一个最佳拟合的圆，你可以将其收缩，直到满足一定比例的内部和外部点。或者，你可以移除圆外的所有点并再次运行最佳拟合圆算法，重复这个过程直到得到满意的答案。
- Louis Ricci
2
1这不会在点的内部返回一个圆，但更可能有一半的点在圆内，另一半的点在圆外。 - SirGuy
@GuyGreer - 一旦你有了最佳拟合圆，你可以将其缩小，直到任何当前未定义的点百分比在其外部。这个问题还没有充分阐述，无法给出确切答案。 - Louis Ricci
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接
R(c-Ni/N)的东西。 - user1196549