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
011
012import imagingbook.common.geometry.basic.Pnt2d;
013import imagingbook.common.geometry.basic.PntUtils;
014import imagingbook.common.math.Matrix;
015import imagingbook.common.math.eigen.EigenDecompositionJama;
016import imagingbook.common.util.PrimitiveSortMap;
017import org.apache.commons.math3.linear.MatrixUtils;
018import org.apache.commons.math3.linear.RealMatrix;
019import org.apache.commons.math3.linear.RealVector;
020import org.apache.commons.math3.linear.SingularValueDecomposition;
021
022import static imagingbook.common.math.Arithmetic.sqr;
023
024/**
025 * <p>
026 * This is an implementation of the "hyperaccurate" algebraic circle fit proposed by Al-Sharadqah and Chernov in [1],
027 * also described in [2] (Sec. 7.5, Eq. 744). This method combines the Pratt and Taubin fits to eliminate the essential
028 * bias. This version is "optimized for stability, not for speed" (see https://people.cas.uab.edu/~mosya/cl/HyperSVD.m).
029 * Singular-value decomposition is used for testing singularity of the data matrix and finding a solution in the
030 * singular case. Fits to exactly 3 (non-collinear) points are handled properly. Data centering is used to improve
031 * numerical stability (alternatively, a reference point can be specified).
032 * </p>
033 * <p>
034 * [1] A. Al-Sharadqah and N. Chernov. "Error analysis for circle fitting algorithms".
035 * <em>Electronic Journal of Statistics</em>, 3:886–911 (2009).
036 * <br>
037 * [2] N. Chernov. "Circular and Linear Regression: Fitting Circles and Lines by Least Squares". Monographs on
038 * Statistics and Applied Probability. Taylor &amp; Francis (2011).
039 * </p>
040 *
041 * @author WB
042 */
043public class CircleFitHyperStable implements CircleFitAlgebraic {
044        
045        private final double[] q;       // q = (A,B,C,D) circle parameters
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 CircleFitHyperStable(Pnt2d[] points) {
053                this.q = fit(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 CircleFitHyperStable(Pnt2d[] points, Pnt2d xref) {
064                this.q = fit(points, xref);
065        }
066        
067        @Override
068        public double[] getParameters() {
069                return this.q;
070        }
071        
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                double[][] Xa = new double[n][];
087                double xs = 0;
088                double ys = 0;
089                double zs = 0;
090                for (int i = 0; i < n; i++) {
091                        double x = pts[i].getX() - xr;
092                        double y = pts[i].getY() - yr;
093                        double z = sqr(x) + sqr(y);
094                        Xa[i] = new double[] {z, x, y, 1};
095                        xs = xs + x;
096                        ys = ys + y;
097                        zs = zs + z;
098                }
099
100                RealMatrix X = MatrixUtils.createRealMatrix(Xa);
101                SingularValueDecomposition svd = new SingularValueDecomposition(X);
102                RealMatrix S = svd.getS();      
103                RealMatrix V = svd.getV();      
104                double[] svals = svd.getSingularValues();       // note: singular values are all positive
105                
106                int k = Matrix.idxMin(svals);
107                double smin = svals[k];
108                double smax = svals[Matrix.idxMax(svals)];
109                double icond = smin / smax;                     // inverse condition number of X        
110                RealVector qq = null;                           // solution vector (circle parameters)
111
112                if (icond < 1e-12) {                            // singular case                
113                        qq = V.getColumnVector(k);              // take the vector for smallest singular value as the solution
114                } 
115                else {                                                          // regular (non-singular) case
116                        double xm = xs / n;
117                        double ym = ys / n;                     
118                        double zm = zs / n;
119                        // data-dependent constraint matrix:
120//                      RealMatrix C = MatrixUtils.createRealMatrix(new double[][]
121//                              {{ 8 * zm, 4 * xm, 4 * ym, 2 },
122//                               { 4 * xm, 1,      0,      0 },
123//                               { 4 * ym, 0,      1,      0 },
124//                               { 2,      0,      0,      0 }});
125//                      RealMatrix Ci = MatrixUtils.inverse(C);
126                        
127                        // define the inverse constraint matrix directly:
128                        RealMatrix Ci = MatrixUtils.createRealMatrix(new double[][]
129                                        {{ 0,   0, 0, 0.5 }, 
130                                         { 0,   1, 0, -2 * xm }, 
131                                         { 0,   0, 1, -2 * ym },
132                                         { 0.5, -2 * xm, -2 * ym, 4 * (sqr(xm) + sqr(ym)) - 2 * zm }});
133
134                        RealMatrix W = V.multiply(S).multiply(V.transpose());
135                        RealMatrix Z = W.multiply(Ci).multiply(W);
136
137//                      EigenDecomposition ed = new EigenDecomposition(Z);
138                        EigenDecompositionJama ed = new EigenDecompositionJama(Z);
139                        double[] evals = ed.getRealEigenvalues();
140                        
141                        int l = new PrimitiveSortMap(evals).getIndex(1);        // index of the 2nd-smallest eigenvalue
142                        RealVector el = ed.getEigenvector(l);
143                        
144//                      qq = Matrix.solve(W, el);       // may return null!
145                        qq = MatrixUtils.inverse(W).operate(el);                                                                                // alt. 1 (qq = Y^-1 * el)
146//                      qq = V.multiply(MatrixUtils.inverse(S)).multiply(V.transpose()).operate(el);    // alt. 2 (S is diagonal!)
147                }
148                
149                // re-adjust circle parameters for data centering:
150                RealMatrix M = CircleFitAlgebraic.getDecenteringMatrix(xr, yr);
151                return M.operate(qq).toArray(); // q = (A, B, C, D)             
152        }
153        
154}