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.ellipse.algebraic;
010
011import imagingbook.common.geometry.basic.Pnt2d;
012import imagingbook.common.math.Matrix;
013import imagingbook.common.math.PrintPrecision;
014import org.apache.commons.math3.linear.Array2DRowRealMatrix;
015import org.apache.commons.math3.linear.ArrayRealVector;
016import org.apache.commons.math3.linear.DecompositionSolver;
017import org.apache.commons.math3.linear.LUDecomposition;
018import org.apache.commons.math3.linear.RealMatrix;
019import org.apache.commons.math3.linear.RealVector;
020import org.apache.commons.math3.linear.SingularMatrixException;
021
022import java.util.Random;
023
024import static imagingbook.common.math.Arithmetic.sqr;
025
026/**
027 * Performs an exact ellipse fit to 5 given points. If the fit is unsuccessful, {@link #getParameters()} returns
028 * {@code null}. The underlying algorithm is described in Section F.3.3 of [1].
029 * <p>
030 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic Introduction</em>, 3rd ed, Springer
031 * (2022).
032 * </p>
033 *
034 * @author WB
035 * @version 2022/11/17
036 */
037public class EllipseFit5Points implements EllipseFitAlgebraic {
038        
039        private final double[] p;       // p = (A,B,C,D,E,F) ellipse parameters
040        
041        public EllipseFit5Points(Pnt2d[] points) {
042                this.p = fit(points);
043        }
044
045        @Override
046        public double[] getParameters() {
047                return this.p;
048        }
049        
050        private double[] fit(Pnt2d[] points) {
051                if (points.length < 5) {
052                        throw new IllegalArgumentException("5 points are needed for " + this.getClass().getSimpleName());
053                }
054                final int n = points.length;
055
056                // create design matrix X:
057                RealMatrix X = new Array2DRowRealMatrix(n, 5);
058                RealVector b = new ArrayRealVector(n);
059                for (int i = 0; i < n; i++) {
060                        double x = points[i].getX();
061                        double y = points[i].getY();
062                        X.setEntry(i, 0, sqr(x) - sqr(y));
063                        X.setEntry(i, 1, x * y);
064                        X.setEntry(i, 2, x);
065                        X.setEntry(i, 3, y);
066                        X.setEntry(i, 4, 1);
067                        b.setEntry(i, - sqr(y));
068                }
069                
070//              System.out.println("X = \n" + Matrix.toString(X));
071//              System.out.println("b = " + Matrix.toString(b));
072                
073                LUDecomposition decomp = new LUDecomposition(X); 
074                // see https://commons.apache.org/proper/commons-math/userguide/linear.html
075                 
076                DecompositionSolver solver = decomp.getSolver();
077                RealVector x = null;
078                try {
079                        x = solver.solve(b);
080                } catch(SingularMatrixException e) { }
081
082                if (x == null) {
083                        return null;
084                }
085                
086                double A = x.getEntry(0);
087                double B = x.getEntry(1);
088                double C = 1 - A;
089                double D = x.getEntry(2);
090                double E = x.getEntry(3);
091                double F = x.getEntry(4);
092                
093                boolean isEllipse = (4 * A * C - sqr(B)) > 0;
094                
095                return (isEllipse) ? new double[] {A,B,C,D,E,F} : null;
096        }
097        
098        // --------------------------------------------------------------------------------------------
099        
100        public static void main(String[] args) {
101                PrintPrecision.set(9);
102                Pnt2d p0 = Pnt2d.from(40, 53);
103                Pnt2d p1 = Pnt2d.from(107, 20);
104                Pnt2d p2 = Pnt2d.from(170, 26);
105                Pnt2d p3 = Pnt2d.from(186, 55);
106                Pnt2d p4 = Pnt2d.from(135, 103);
107                
108                Pnt2d[] points = {p0, p1, p2, p3, p4};
109                EllipseFit5Points fit = new EllipseFit5Points(points);
110                System.out.println("fit parameters = " + Matrix.toString(fit.getParameters()));
111                System.out.println("fit ellipse = " + fit.getEllipse());
112                System.out.println("fit ellipse = " +  Matrix.toString(fit.getEllipse().getParameters()));
113                
114                // create random 5-point sets and try to fit ellipses, counting null results:
115                Random rg = new Random(17);
116                int N = 1000;
117                int nullCnt = 0;
118                for (int k = 0; k < N; k++) {
119                        for (int i = 0; i < points.length; i++) {
120                                double x = rg.nextInt(200);
121                                double y = rg.nextInt(200);
122                                points[i] = Pnt2d.from(x, y);
123                        }
124                        fit = new EllipseFit5Points(points);
125                        if (fit.getEllipse() == null) {
126                                nullCnt++;
127                        }
128//                      System.out.println("fit ellipse = " + fit.getEllipse());
129                        
130                }
131                
132                System.out.println(nullCnt + " null results out of " + N);
133                
134                // fit ellipse = {0.317325319, 0.332954818, 0.875173557, -89.442594143, -150.574265066, 7886.192568730}
135        }
136}