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 org.apache.commons.math3.linear.MatrixUtils; 014import org.apache.commons.math3.linear.RealMatrix; 015import org.apache.commons.math3.linear.SingularMatrixException; 016import org.apache.commons.math3.linear.SingularValueDecomposition; 017 018import static imagingbook.common.math.Arithmetic.sqr; 019 020/** 021 * This is an improved implementation of the Kåsa [1] circle fitting algorithm described in [2, Sec. 5.2] (Eq. 5.12). It 022 * is based on the Moore-Penrose pseudo-inverse which is applied to the full data matrix (i.e, no 3x3 scatter matrix is 023 * mounted). See also [3, Sec. 11.1.2] and {@link CircleFitKasaA} for the original version. 024 * <p> 025 * This algorithm is assumed to be numerically more stable than solutions based on solving a 3x3 system. The 026 * pseudo-inverse is obtained by singular-value decomposition (SVD). However, the significant bias on points sampled 027 * from a small circle segment remains. Fits to exactly 3 (non-collinear) points are handled properly. No data centering 028 * (which should improve numerical stability) is used. 029 * </p> 030 * <p> 031 * [1] I. Kåsa. "A circle fitting procedure and its error analysis", <em>IEEE Transactions on Instrumentation and 032 * Measurement</em> <strong>25</strong>(1), 8–14 (1976). 033 * <br> 034 * [2] N. Chernov. "Circular and Linear Regression: Fitting Circles and Lines by Least Squares". Monographs on 035 * Statistics and Applied Probability. Taylor & Francis (2011). 036 * <br> 037 * [3] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 3rd ed, Springer 038 * (2022). 039 * </p> 040 * 041 * @author WB 042 */ 043public class CircleFitKasaC implements CircleFitAlgebraic { 044 045 private final double[] q; // q = (B,C,D) circle parameters, A=1 046 047 /** 048 * Constructor. The centroid of the sample points is used as the reference point. 049 * 050 * @param points sample points 051 */ 052 public CircleFitKasaC(Pnt2d[] points) { 053 this(points, null); 054 } 055 056 /** 057 * Constructor. The centroid of the sample points is used as the reference point for data centering if {@code null} 058 * is passed for {@code xref}. 059 * 060 * @param points sample points 061 * @param xref reference point or {@code null} 062 */ 063 public CircleFitKasaC(Pnt2d[] points, Pnt2d xref) { 064 this.q = fit(points, xref); 065 } 066 067 @Override 068 public double[] getParameters() { 069 return q; 070 } 071 072 private double[] fit(Pnt2d[] pts, Pnt2d xref) { 073 final int n = pts.length; 074 if (n < 3) { 075 throw new IllegalArgumentException("at least 3 points are required"); 076 } 077 078 if (xref == null) { 079 xref = PntUtils.centroid(pts); 080 } 081 final double xr = xref.getX(); 082 final double yr = xref.getY(); 083 084 085 final double[] z = new double[n]; 086 final double[][] Xa = new double[n][]; 087 for (int i = 0; i < n; i++) { 088 final double x = pts[i].getX() - xr; 089 final double y = pts[i].getY() - yr; 090 Xa[i] = new double[] {x, y, 1}; 091 z[i] = -(sqr(x) + sqr(y)); 092 } 093 094 RealMatrix X = MatrixUtils.createRealMatrix(Xa); 095 096 RealMatrix Xi = null; 097 try { 098 SingularValueDecomposition svd = new SingularValueDecomposition(X); 099 Xi = svd.getSolver().getInverse(); // get (3,N) pseudoinverse of X 100 } catch (SingularMatrixException e) { } 101 102 if (Xi == null) { 103 return null; 104 } 105 else { 106 double[] qq = Xi.operate(z); // solution vector qq = X^-1 * z = (B, C, D) 107 // re-adjust for data centering 108 RealMatrix M = CircleFitAlgebraic.getDecenteringMatrix(xr, yr); 109 return M.operate(new double[] {1, qq[0], qq[1], qq[2]}); // q = (A,B,C,D) 110 } 111 } 112 113}