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 imagingbook.common.math.Arithmetic; 014import imagingbook.common.math.Matrix; 015import org.apache.commons.math3.fitting.leastsquares.LeastSquaresFactory; 016import org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer; 017import org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.Optimum; 018import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem; 019import org.apache.commons.math3.fitting.leastsquares.LevenbergMarquardtOptimizer; 020import org.apache.commons.math3.fitting.leastsquares.MultivariateJacobianFunction; 021import org.apache.commons.math3.linear.Array2DRowRealMatrix; 022import org.apache.commons.math3.linear.ArrayRealVector; 023import org.apache.commons.math3.linear.MatrixUtils; 024import org.apache.commons.math3.linear.RealMatrix; 025import org.apache.commons.math3.linear.RealVector; 026import org.apache.commons.math3.optim.SimpleVectorValueChecker; 027import org.apache.commons.math3.util.Pair; 028 029import java.util.LinkedList; 030import java.util.List; 031 032import static imagingbook.common.math.Arithmetic.sqr; 033import static java.lang.Math.sqrt; 034import static org.apache.commons.math3.fitting.leastsquares.LeastSquaresFactory.evaluationChecker; 035import static org.apache.commons.math3.linear.MatrixUtils.createRealVector; 036 037/** 038 * <p> 039 * "Coordinate-based" geometric circle fitter using a nonlinear least-squares (Levenberg-Marquart) optimizer. See [1, 040 * Sec. 11.1.3] for a detailed description (Alg. 11.4). 041 * </p> 042 * <p> 043 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 3rd ed, Springer 044 * (2022). 045 * </p> 046 * 047 * @author WB 048 * @version 2022/11/17 049 */ 050public class CircleFitGeometricCoord implements CircleFitGeometric { 051 052 public static final int DefaultMaxIterations = 1000; 053 public static final double DefaultTolerance = 1e-6; 054 055 private final Pnt2d[] pts; 056 private final double[] V; 057 private final double[][] J; 058 059 private final Optimum solution; 060 private final List<double[]> history = new LinkedList<>(); 061 062 public CircleFitGeometricCoord(Pnt2d[] pts, GeometricCircle initCircle) { 063 this(pts, initCircle, DefaultMaxIterations, DefaultMaxIterations, DefaultTolerance); 064 } 065 066 public CircleFitGeometricCoord(Pnt2d[] pts, GeometricCircle initCircle, int maxEvaluations, int maxIterations, double tolerance) { 067 this.pts = pts; 068 this.V = new double[2 * pts.length]; 069 this.J = new double[2 * pts.length][3]; 070 071 LeastSquaresProblem problem = 072 LeastSquaresFactory.create( 073 new AnalyticModel(), // model(V, J), 074 makeTargetVector(pts), 075 createRealVector(initCircle.getParameters()), 076 evaluationChecker(new SimpleVectorValueChecker(tolerance, tolerance)), 077 maxEvaluations, maxIterations); 078 079 LeastSquaresOptimizer optimizer = new LevenbergMarquardtOptimizer(); // new GaussNewtonOptimizer(); 080 this.solution = optimizer.optimize(problem); 081 } 082 083 // -------------------------------------------------------------------------- 084 085 @Override 086 public double[] getParameters() { 087 return solution.getPoint().toArray(); 088 } 089 090 @Override 091 public int getIterations() { 092 return solution.getIterations(); 093 } 094 095 @Override 096 public List<double[]> getHistory() { 097 return history; 098 } 099 100 101 private RealVector makeTargetVector(final Pnt2d[] X) { 102 final int n = X.length; 103 final double[] target = new double[2*n]; 104 for (int i = 0; i < n; i++) { 105 target[2*i] = X[i].getX(); // x_i 106 target[2*i+1] = X[i].getY(); // y_i 107 } 108 return MatrixUtils.createRealVector(target); 109 } 110 111 // -------------------------------------------------------------------------- 112 113 /** 114 * Defines function {@link #value(RealVector)} which returns the vector of model values and the associated Jacobian 115 * matrix for a given parameter point. 116 */ 117 class AnalyticModel implements MultivariateJacobianFunction { 118 int k = 0; 119 120 @Override 121 public Pair<RealVector, RealMatrix> value(RealVector point) { 122 k++; 123 final double[] p = point.toArray(); 124 if (RecordHistory) { 125 history.add(p.clone()); 126 } 127 final double xc = p[0], yc = p[1], r = p[2]; 128 for (int i = 0; i < pts.length; i++) { 129 final double dx = pts[i].getX() - xc; 130 final double dy = pts[i].getY() - yc; 131 final double ri2 = sqr(dx) + sqr(dy); // r_i^2 132 final double ri = sqrt(ri2); // r_i 133 final double ri3 = ri2 * ri; // r_i^3 134 135 if (Arithmetic.isZero(ri)) { 136 throw new ArithmeticException("zero radius ri encountered"); 137 } 138 139 V[2*i] = xc + dx * r / ri; 140 V[2*i+1] = yc + dy * r / ri; 141 142 // Jacobian: 143 J[2*i][0] = 1 + (r / ri) * (sqr(dx) / ri2 - 1); // 1 + (r * sqr(dx) / ri3) - (r / ri); 144 J[2*i][1] = r * dx * dy / ri3; 145 J[2*i][2] = dx / ri; 146 147 J[2*i+1][0] = r * dx * dy / ri3; 148 J[2*i+1][1] = 1 + (r / ri) * (sqr(dy) / ri2 - 1); // 1 + (r * sqr(dy) / ri3) - (r / ri); 149 J[2*i+1][2] = dy / ri; 150 } 151 if (k == 1) { 152 System.out.println("V = " + Matrix.toString(V)); 153 System.out.println("J = \n" + Matrix.toString(J)); 154 } 155 RealVector VV = new ArrayRealVector(V, false); 156 RealMatrix JJ = new Array2DRowRealMatrix(J, false); 157 return new Pair<>(VV, JJ); 158 } 159 } 160 161 // ------------------------------------------------------------------- 162 163// public static void main(String[] args) { 164// PrintPrecision.set(3); 165//// CircleFitGeometric.RecordHistory = true; 166// 167// Pnt2d[] points = { 168// Pnt2d.from(15,9), 169// Pnt2d.from(68,33), 170// Pnt2d.from(35,69), 171// Pnt2d.from(17,51), 172// Pnt2d.from(90,54) 173// }; 174// 175// GeometricCircle estimA = new CircleFitPratt(points).getGeometricCircle(); 176// System.out.println("estimA: " + estimA.toString()); 177// System.out.println("estimA error = " + estimA.getMeanSquareError(points)); 178// 179// GeometricCircle init = new GeometricCircle(45, 40, 30); // Example (a) 180//// GeometricCircle init = new GeometricCircle(75, 75, 12); // Example (b) 181// //GeometricCircle init = estimA; 182//// GeometricCircle init = estimA.disturb(0, 0, 0); 183// System.out.println(" init: " + init.toString()); 184// System.out.println(" init error = " + init.getMeanSquareError(points)); 185// 186// CircleFitGeometricCoord geomfitter = new CircleFitGeometricCoord(points, init); 187// GeometricCircle circleG = geomfitter.getCircle(); 188//// 189//// //Circle2D refined = Doube.levenMarqFull(points, init); 190//// 191// System.out.println("circleG: " + circleG.toString()); 192// System.out.println("iterations: " + geomfitter.getIterations()); 193// System.out.println("final error = " + circleG.getMeanSquareError(points)); 194// for (double[] p : geomfitter.getHistory()) { 195// System.out.println(Matrix.toString(p)); 196// } 197// } 198 199}