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.RealMatrix;
015
016import static imagingbook.common.math.Arithmetic.sqr;
017
018/**
019 * This is an implementation of the modified Kåsa [1] circle fitting algorithm described in [2, Sec. 5.1]. A description
020 * of the concrete algorithm can be found in [3, Alg. 11.1]. See {@link CircleFitKasaA} for the original version.
021 * <p>
022 * Compared to the original Kåsa algorithm, this variant also solves a 3x3 linear system but uses a slightly different
023 * setup of the scatter matrix (using only powers of 2 instead of 3). A numerical solver is used for this purpose. The
024 * algorithm is fast but shares the same numerical instabilities and bias when sample points are taken from a small
025 * circle segment. It fails when matrix M becomes singular. Fits to exactly 3 (non-collinear) points are handled
026 * properly. Data centering is used to improve numerical stability (alternatively, a reference point can be specified).
027 * </p>
028 * <p>
029 * [1] I. Kåsa. "A circle fitting procedure and its error analysis",
030 * <em>IEEE Transactions on Instrumentation and Measurement</em> <strong>25</strong>(1),
031 * 8–14 (1976).
032 * <br>
033 * [2] N. Chernov. "Circular and Linear Regression: Fitting Circles and Lines by Least Squares". Monographs on
034 * Statistics and Applied Probability. Taylor &amp; Francis (2011).
035 * <br>
036 * [3] W. Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic Introduction</em>, 3rd ed, Springer
037 * (2022).
038 * </p>
039 *
040 * @author WB
041 * @see CircleFitKasaA
042 * @see CircleFitKasaC
043 */
044public class CircleFitKasaB implements CircleFitAlgebraic {
045
046        private final double[] q;       // q = (B,C,D) circle parameters, A=1
047
048        /**
049         * Constructor. The centroid of the sample points is used as the reference point.
050         *
051         * @param points sample points
052         */
053        public CircleFitKasaB(Pnt2d[] points) {
054                this(points, null);
055        }
056
057        /**
058         * Constructor. The centroid of the sample points is used as the reference point for data centering if {@code null}
059         * is passed for {@code xref}.
060         *
061         * @param points sample points
062         * @param xref reference point or {@code null}
063         */
064        public CircleFitKasaB(Pnt2d[] points, Pnt2d xref) {
065                this.q = fit(points, xref);
066        }
067        
068        
069        @Override
070        public double[] getParameters() {
071                return q;
072        }
073        
074        private double[] fit(Pnt2d[] pts, Pnt2d xref) {
075                final int n = pts.length;
076                if (n < 3) {
077                        throw new IllegalArgumentException("at least 3 points are required");
078                }
079                
080                if (xref == null) {
081                        xref = PntUtils.centroid(pts);
082                }
083                final double xr = xref.getX();
084                final double yr = xref.getY();
085
086                // calculate elements of scatter matrix
087                double sx = 0, sy = 0, sz = 0;
088                double sxy = 0, sxx = 0, syy = 0, sxz = 0, syz = 0;
089                for (int i = 0; i < n; i++) {
090                        final double x = pts[i].getX() - xr;
091                        final double y = pts[i].getY() - yr;
092                        final double x2 = sqr(x);
093                        final double y2 = sqr(y);
094                        final double z = x2 + y2;
095                        sx  += x;
096                        sy  += y;
097                        sz  += z;
098                        sxx += x2;
099                        syy += y2;
100                        sxy += x * y;   
101                        sxz += x * z;
102                        syz += y * z;
103                }
104                
105                double[][] X = {                                // scatter matrix X
106                                {sxx, sxy, sx},
107                                {sxy, syy, sy},
108                                { sx,  sy,  n}};
109            
110                double[] b = {-sxz, -syz, -sz};  // RHS vector
111                double[] qq = Matrix.solve(X, b); // solve M * qq = b (exact), for parameter vector qq = (B, C, D)
112                if (qq == null) {
113                        return null;    // M is singular, no solution
114                }
115                else {
116                        // re-adjust for data centering
117                        RealMatrix M = CircleFitAlgebraic.getDecenteringMatrix(xr, yr);
118                        return M.operate(new double[] {1, qq[0], qq[1], qq[2]});        // q = (A,B,C,D)
119                }
120        }
121
122}