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.eigen.EigenDecompositionJama;
014import org.apache.commons.math3.linear.MatrixUtils;
015import org.apache.commons.math3.linear.RealMatrix;
016
017import static imagingbook.common.math.Arithmetic.sqr;
018
019/**
020 * <p>
021 * Algebraic ellipse fit using Fitzgibbon's original method [1], based on simple matrix inversion. See [2, Sec. 11.2.1]
022 * for a detailed description. Note: This implementation does not use data centering nor accepts a specific reference
023 * point. With exactly 5 input points (generally sufficient for ellipse fitting) the scatter matrix X is singular and in
024 * this case matrix S has no inverse. Thus at least 6 distinct input points are required (i.e., no duplicate points are
025 * allowed).
026 * </p>
027 * <p>
028 * [1] A. W. Fitzgibbon, M. Pilu, and R. B. Fisher. Direct least- squares fitting of ellipses. IEEE Transactions on
029 * Pattern Analysis and Machine Intelligence 21(5), 476-480 (1999).<br> [2] W. Burger, M.J. Burge, <em>Digital Image
030 * Processing &ndash; An Algorithmic Introduction</em>, 3rd ed, Springer (2022).
031 * </p>
032 *
033 * @author WB
034 * @version 2022/11/17
035 */
036public class EllipseFitFitzgibbonNaive implements EllipseFitAlgebraic {
037        
038        // constraint matrix
039        private static final RealMatrix C = Matrix.makeRealMatrix(6, 6,
040                                0, 0, 2, 0, 0, 0,
041                                0,-1, 0, 0, 0, 0,
042                                2, 0, 0, 0, 0, 0,
043                                0, 0, 0, 0, 0, 0,
044                                0, 0, 0, 0, 0, 0,
045                                0, 0, 0, 0, 0, 0 );
046        
047        private final double[] q;       // = (A,B,C,D,E,F) ellipse parameters
048        
049        public EllipseFitFitzgibbonNaive(Pnt2d[] points) {
050                if (points.length < 6) {
051                        throw new IllegalArgumentException("fitter requires at least 6 sample points instead of " + points.length);
052                }
053                this.q = fit(points);
054        }
055
056        @Override
057        public double[] getParameters() {
058                return this.q;
059        }
060        
061        private double[] fit(Pnt2d[] points) {
062                final int n = points.length;
063
064                // create design matrix X:
065                RealMatrix X = MatrixUtils.createRealMatrix(n, 6);
066                for (int i = 0; i < n; i++) {
067                        double x = points[i].getX();
068                        double y = points[i].getY();
069                        X.setEntry(i, 0, sqr(x));
070                        X.setEntry(i, 1, x * y);
071                        X.setEntry(i, 2, sqr(y));
072                        X.setEntry(i, 3, x);
073                        X.setEntry(i, 4, y);
074                        X.setEntry(i, 5, 1);
075                }
076                
077                // scatter matrix S:
078                RealMatrix S = X.transpose().multiply(X);       
079                RealMatrix Si = MatrixUtils.inverse(S, 1e-15);
080                
081//              EigenDecomposition ed = new EigenDecomposition(Si.multiply(C));
082                EigenDecompositionJama ed = new EigenDecompositionJama(Si.multiply(C));
083
084                double[] evals = ed.getRealEigenvalues();
085                int k = Matrix.idxMax(evals);                           // index of the largest eigenvalue
086                return ed.getEigenvector(k).toArray();
087        }
088        
089}