001/******************************************************************************* 002 * This software is provided as a supplement to the authors' textbooks on digital 003 * image processing published by Springer-Verlag in various languages and editions. 004 * Permission to use and distribute this software is granted under the BSD 2-Clause 005 * "Simplified" License (see http://opensource.org/licenses/BSD-2-Clause). 006 * Copyright (c) 2006-2023 Wilhelm Burger, Mark J. Burge. All rights reserved. 007 * Visit https://imagingbook.com for additional details. 008 ******************************************************************************/ 009package imagingbook.common.geometry.fitting.circle.geometric; 010 011import imagingbook.common.geometry.basic.Pnt2d; 012import imagingbook.common.geometry.circle.GeometricCircle; 013import org.apache.commons.math3.fitting.leastsquares.LeastSquaresFactory; 014import org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer; 015import org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.Optimum; 016import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem; 017import org.apache.commons.math3.fitting.leastsquares.LevenbergMarquardtOptimizer; 018import org.apache.commons.math3.fitting.leastsquares.MultivariateJacobianFunction; 019import org.apache.commons.math3.linear.Array2DRowRealMatrix; 020import org.apache.commons.math3.linear.ArrayRealVector; 021import org.apache.commons.math3.linear.RealMatrix; 022import org.apache.commons.math3.linear.RealVector; 023import org.apache.commons.math3.optim.SimpleVectorValueChecker; 024import org.apache.commons.math3.util.Pair; 025 026import java.util.LinkedList; 027import java.util.List; 028 029import static imagingbook.common.math.Arithmetic.isZero; 030import static java.lang.Math.hypot; 031import static org.apache.commons.math3.fitting.leastsquares.LeastSquaresFactory.evaluationChecker; 032 033/** 034 * <p> 035 * "Distance-based" geometric circle fitter using a nonlinear least-squares (Levenberg-Marquart) optimizer. See [1, Sec. 036 * 11.1.3] for a detailed description (Alg. 11.3). 037 * </p> 038 * <p> 039 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 3rd ed, Springer 040 * (2022). 041 * </p> 042 * 043 * @author WB 044 * @version 2022/11/17 045 */ 046public class CircleFitGeometricDist implements CircleFitGeometric { 047 048 public static final int DefaultMaxIterations = 1000; 049 public static final double DefaultTolerance = 1e-6; 050 051 private final Pnt2d[] pts; 052 private final double[] V; 053 private final double[][] J; 054 private final Optimum solution; 055 private final List<double[]> history = new LinkedList<>(); 056 057 public CircleFitGeometricDist(Pnt2d[] pts, GeometricCircle initCircle) { 058 this(pts, initCircle, DefaultMaxIterations, DefaultMaxIterations, DefaultTolerance); 059 } 060 061 public CircleFitGeometricDist(Pnt2d[] pts, GeometricCircle initCircle, int maxEvaluations, int maxIterations, double tolerance) { 062 this.pts = pts; 063 this.V = new double[pts.length]; 064 this.J = new double[pts.length][3]; 065 066 LeastSquaresProblem problem = 067 LeastSquaresFactory.create( 068 new AnalyticModel(), // model(V, J), 069 new ArrayRealVector(pts.length), // zero vector 070 new ArrayRealVector(initCircle.getParameters()), 071 evaluationChecker(new SimpleVectorValueChecker(tolerance, tolerance)), 072 maxEvaluations, maxIterations); 073 074 LeastSquaresOptimizer optimizer = new LevenbergMarquardtOptimizer(); // new GaussNewtonOptimizer(); 075 this.solution = optimizer.optimize(problem); 076 } 077 078 079 // -------------------------------------------------------------------------- 080 081 @Override 082 public double[] getParameters() { 083 return solution.getPoint().toArray(); 084 } 085 086 @Override 087 public int getIterations() { 088 return solution.getIterations(); 089 } 090 091 @Override 092 public List<double[]> getHistory() { 093 return history; 094 } 095 096 // -------------------------------------------------------------------------- 097 098 /** 099 * Defines function {@link #value(RealVector)} which returns the vector of model values and the associated Jacobian 100 * matrix for a given parameter point. 101 */ 102 class AnalyticModel implements MultivariateJacobianFunction { 103 104 @Override 105 public Pair<RealVector, RealMatrix> value(RealVector point) { 106 final double[] p = point.toArray(); 107 if (RecordHistory) { 108 history.add(p.clone()); 109 } 110 111 final double xc = p[0], yc = p[1], r = p[2]; 112 for (int i = 0; i < pts.length; i++) { 113 double dx = xc - pts[i].getX(); 114 double dy = yc - pts[i].getY(); 115 double ri = hypot(dx, dy); 116 117 // values: 118 V[i] = ri - r; 119 120 // Jacobian: 121 if (isZero(ri)) { // (xi,yi) is at the circle center 122 J[i][0] = 0; // (xc - xi) = 0 too 123 J[i][1] = 0; // (yc - yi) = 0 too 124 } else { 125 J[i][0] = dx / ri; 126 J[i][1] = dy / ri; 127 } 128 J[i][2] = -1; 129 } 130 RealVector VV = new ArrayRealVector(V, false); 131 RealMatrix JJ = new Array2DRowRealMatrix(J, false); 132 return new Pair<>(VV, JJ); 133 } 134 135 } 136 137 138 139 // ------------------------------------------------------------------- 140 141// public static void main(String[] args) { 142// PrintPrecision.set(9); 143// //Pnt2d[] points = CircleMaker.makeTestCircle(XC, YC, R, 100, Angle0, Angle1, SigmaNoise); 144// //Pnt2d[] points = CircleMaker.makeTestCircle(XC, YC, R, 4, Angle0, Angle1, 0.1); 145// //Pnt2d[] points = CircleMaker.makeTestGander(30); 146// //Pnt2d[] points = CircleSampler.make3Points(30); 147// 148// Pnt2d[] points = { 149// Pnt2d.from(15,9), 150// Pnt2d.from(68,33), 151// Pnt2d.from(35,69), 152// Pnt2d.from(17,51), 153// Pnt2d.from(90,54) 154// }; 155// 156//// Circle2D real = new Circle2D(XC, YC, R); 157//// IJ.log(" real: " + real.toString()); 158// 159//// Pnt2d[] points = makeTestPoints(); 160// 161// GeometricCircle estimA = new CircleFitPratt(points).getGeometricCircle(); 162// IJ.log("estimA: " + estimA.toString()); 163// IJ.log("estimA error = " + estimA.getMeanSquareError(points)); 164// 165// //Circle2D init = new Circle2D(XC, YC, R); 166// //Circle2D init = estimA; 167// GeometricCircle init = disturbCircle(estimA, -1.5, 0.5, 10); 168// IJ.log(" init: " + init.toString()); 169// IJ.log(" init error = " + init.getMeanSquareError(points)); 170// 171// CircleFitGeometricDist geomfitter = new CircleFitGeometricDist(points, init); 172// GeometricCircle refined = geomfitter.getCircle(); 173// 174// //Circle2D refined = Doube.levenMarqFull(points, init); 175// 176// IJ.log("refin: " + refined.toString()); 177// IJ.log("iterations: " + geomfitter.getIterations()); 178// 179// IJ.log("final error = " + refined.getMeanSquareError(points)); 180// } 181// 182// private static GeometricCircle disturbCircle (GeometricCircle circle, double dxc, double dyc, double dr) { 183// return new GeometricCircle(circle.getXc() + dxc, circle.getYc() + dyc, circle.getR() + dr); 184// } 185}