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.algebraic; 010 011import imagingbook.common.geometry.basic.Pnt2d; 012import imagingbook.common.geometry.basic.PntUtils; 013import imagingbook.common.math.Matrix; 014import imagingbook.common.math.eigen.EigenDecompositionJama; 015import imagingbook.common.util.PrimitiveSortMap; 016import org.apache.commons.math3.linear.MatrixUtils; 017import org.apache.commons.math3.linear.RealMatrix; 018import org.apache.commons.math3.linear.RealVector; 019import org.apache.commons.math3.linear.SingularValueDecomposition; 020 021import static imagingbook.common.math.Arithmetic.sqr; 022 023 024/** 025 * This is an implementation of the algebraic circle fitting algorithm by Pratt [1], as described in [2] (Sec. 5.5-5.6). 026 * The algorithm uses singular-value decomposition (SVD) and eigen-decomposition. See [3, Alg. 11.2] for additional 027 * details. 028 * <p> 029 * Fits to exactly 3 (non-collinear) points are handled properly. Data centering is used to improve numerical stability 030 * (alternatively, a reference point can be specified). 031 * </p> 032 * <p> 033 * [1] V. Pratt. "Direct least-squares fitting of algebraic surfaces". <em>ACM SIGGRAPH Computer Graphics</em> 034 * <strong>21</strong>(4), 145–152 (July 1987). 035 * <br> 036 * [2] N. Chernov. "Circular and Linear Regression: Fitting Circles and Lines by Least Squares". Monographs on 037 * Statistics and Applied Probability. Taylor & Francis (2011). 038 * <br> 039 * [3] 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 */ 045public class CircleFitPratt implements CircleFitAlgebraic { 046 047 private static final RealMatrix Ci = // inverse of constraint matrix C 048 MatrixUtils.createRealMatrix(new double[][] { 049 { 0, 0, 0, -0.5 }, 050 { 0, 1, 0, 0 }, 051 { 0, 0, 1, 0 }, 052 { -0.5, 0, 0, 0 }}); 053 054 private final double[] q; // q = (A,B,C,D) circle parameters 055 056 @Override 057 public double[] getParameters() { 058 return this.q; 059 } 060 061 /** 062 * Constructor. The centroid of the sample points is used as the reference point. 063 * 064 * @param points sample points 065 */ 066 public CircleFitPratt(Pnt2d[] points) { 067 this(points, null); 068 } 069 070 /** 071 * Constructor. The centroid of the sample points is used as the reference point for data centering if {@code null} 072 * is passed for {@code xref}. 073 * 074 * @param points sample points 075 * @param xref reference point or {@code null} 076 */ 077 public CircleFitPratt(Pnt2d[] points, Pnt2d xref) { 078 this.q = fit(points, xref); 079 } 080 081 private double[] fit(Pnt2d[] pts, Pnt2d xref) { 082 final int n = pts.length; 083 if (n < 3) { 084 throw new IllegalArgumentException("at least 3 points are required"); 085 } 086 087 if (xref == null) { 088 xref = PntUtils.centroid(pts); 089 } 090 final double xr = xref.getX(); 091 final double yr = xref.getY(); 092 093 double[][] Xa = new double[Math.max(n, 4)][4]; // Xa must have at least 4 rows! 094 for (int i = 0; i < pts.length; i++) { 095 double x = pts[i].getX() - xr; // = x_i 096 double y = pts[i].getY() - yr; // = y_i 097 Xa[i][0] = sqr(x) + sqr(y); // = z_i 098 Xa[i][1] = x; 099 Xa[i][2] = y; 100 Xa[i][3] = 1; 101 } 102 // if nPoints = 3 (special case) the last row of the 103 // 4x4 matrix contains all zeros (make X singular)! 104 105 RealMatrix X = MatrixUtils.createRealMatrix(Xa); 106 107 SingularValueDecomposition svd = new SingularValueDecomposition(X); 108 RealMatrix S = svd.getS(); 109 RealMatrix V = svd.getV(); 110 double[] svals = svd.getSingularValues(); // note: singular values are all positive (>= 0) 111 112 int k = Matrix.idxMin(svals); 113 double smin = svals[k]; 114 double smax = Matrix.max(svals); 115 116 RealVector qq = null; // = \dot{q} solution vector (algebraic circle parameters) 117 118 double icond = smin / smax; 119 if (icond < 1e-12) { // smin/smax = inverse condition number of X, 1e-12 120 // singular case (X is rank deficient) 121 qq = V.getColumnVector(k); 122 } else { 123 // regular (non-singular) case 124 125 // Version1 (seems to create smaller roundoff errors, better matrix symmetry): 126 RealMatrix Y = V.multiply(S).multiply(V.transpose()); 127 RealMatrix Z = Y.multiply(Ci).multiply(Y); // = Y * Ci * Y 128 129 // Version2: 130// RealMatrix Y = V.multiply(S); 131// RealMatrix Z = Y.transpose().multiply(Ci).multiply(Y); // = Y^T * Ci * Y 132 133 EigenDecompositionJama ed = new EigenDecompositionJama(Z); 134 double[] evals = ed.getRealEigenvalues(); 135 int l = new PrimitiveSortMap(evals).getIndex(1); // index of the 2nd-smallest eigenvalue 136 RealVector el = ed.getEigenvector(l); 137 138 // Version1 --------------------------------------------------- 139// qq = Matrix.solve(S.multiply(svd.getVT()), el); // solve S * V^T * p = el 140 141 // Version2 --------------------------------------------------- 142 qq = V.operate(MatrixUtils.inverse(S).operate(el)); // simpler since S is diagonal (i.e., easy to invert) 143 } 144 145 RealMatrix M = CircleFitAlgebraic.getDecenteringMatrix(xr, yr); 146 return M.operate(qq).toArray(); // q = (A,B,C,D) 147 } 148 149}