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 org.apache.commons.math3.linear.MatrixUtils; 015import org.apache.commons.math3.linear.RealMatrix; 016import org.apache.commons.math3.linear.SingularValueDecomposition; 017 018import static imagingbook.common.math.Arithmetic.sqr; 019 020/** 021 * This is an implementation of the algebraic circle fitting algorithm by Taubin [1], following the description in [2] 022 * (Sec. 5.9-5.10). 023 * <p> 024 * The algorithm uses singular-value decomposition (SVD). Fits to exactly 3 (non-collinear) points are handled properly. 025 * Data centering is used to improve numerical stability (alternatively, a reference point can be specified). 026 * </p> 027 * <p> 028 * [1] G. Taubin. "Estimation of planar curves, surfaces, and nonplanar space curves defined by implicit equations with 029 * applications to edge and range image segmentation". <em>IEEE Transactions on Pattern Analysis and Machine 030 * Intelligence</em> <strong>13</strong>(11), 1115–1138 (1991). 031 * <br> 032 * [2] N. Chernov. "Circular and Linear Regression: Fitting Circles and Lines by Least Squares". Monographs on 033 * Statistics and Applied Probability. Taylor & Francis (2011). 034 * </p> 035 * 036 * @author WB 037 */ 038public class CircleFitTaubin implements CircleFitAlgebraic { 039 040 private final double[] q; // p = (A,B,C,D) circle parameters 041 042 /** 043 * Constructor. The centroid of the sample points is used as the reference point. 044 * 045 * @param points sample points 046 */ 047 public CircleFitTaubin(Pnt2d[] points) { 048 this(points, null); 049 } 050 051 /** 052 * Constructor. The centroid of the sample points is used as the reference point for data centering if {@code null} 053 * is passed for {@code xref}. 054 * 055 * @param points sample points 056 * @param xref reference point or {@code null} 057 */ 058 public CircleFitTaubin(Pnt2d[] points, Pnt2d xref) { 059 this.q = fit(points, xref); 060 } 061 062 private double[] fit(Pnt2d[] pts, Pnt2d xref) { 063 int n = pts.length; 064 if (n < 3) 065 throw new IllegalArgumentException("at least 3 points are required"); 066 067 if (xref == null) { 068 xref = PntUtils.centroid(pts); 069 } 070 double xr = xref.getX(); 071 double yr = xref.getY(); 072 073 double[][] Xa = new double[n][3]; // Xa[i] = (z, x, y) 074 double zs = 0; 075 for (int i = 0; i < n; i++) { 076 double x = pts[i].getX() - xr; // = x_i 077 double y = pts[i].getY() - yr; // = y_i 078 double z = sqr(x) + sqr(y); // = z_i 079 zs = zs + z; 080 Xa[i][0] = z; 081 Xa[i][1] = x; 082 Xa[i][2] = y; 083 } 084 085 double zm = zs / n; 086 double zmr = Math.sqrt(zm); 087 088 // normalize z to zero mean, unit variance: 089 for (int i = 0; i < n; i++) { 090 Xa[i][0] = (Xa[i][0] - zm) / (2 * zmr); 091 } 092 093 RealMatrix X = MatrixUtils.createRealMatrix(Xa); 094 SingularValueDecomposition svd = new SingularValueDecomposition(X); 095 RealMatrix V = svd.getV(); 096 097 // get the column vector of V associated with the smallest singular value: 098 double[] svals = svd.getSingularValues(); 099 int minIdx = Matrix.idxMin(svals); 100 double[] a = V.getColumn(minIdx); 101 102 double[] qq = new double[4]; 103 qq[0] = a[0] / (2 * zmr); 104 qq[1] = a[1]; 105 qq[2] = a[2]; 106 qq[3] = -a[0] * zmr / 2; 107 108 // re-adjust for data centering: 109 RealMatrix M = CircleFitAlgebraic.getDecenteringMatrix(xr, yr); 110 return M.operate(qq); // q = (A,B,C,D) 111 } 112 113 @Override 114 public double[] getParameters() { 115 return this.q; 116 } 117 118}