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.geometry.ellipse.AlgebraicEllipse;
013import imagingbook.common.math.Matrix;
014import org.apache.commons.math3.linear.RealMatrix;
015
016import static imagingbook.common.math.Arithmetic.sqr;
017
018/**
019 * Interface for algebraic ellipse fits.
020 * 
021 * @author WB
022 */
023public interface EllipseFitAlgebraic {
024        
025        public enum FitType {
026                FitzgibbonNaive,
027                FitzgibbonOriginal,
028                FitzgibbonStable,
029                Taubin1,
030                Taubin2
031        }
032        
033        public static EllipseFitAlgebraic getFit(FitType type, Pnt2d[] points, Pnt2d xref) {
034                switch (type) {
035                case FitzgibbonNaive:           return new EllipseFitFitzgibbonNaive(points);
036                case FitzgibbonOriginal:        return new EllipseFitFitzgibbonOriginal(points);
037                case FitzgibbonStable:          return new EllipseFitFitzgibbonStable(points, xref);
038                case Taubin1:                           return new EllipseFitTaubin1(points, xref);
039                case Taubin2:                           return new EllipseFitTaubin2(points, xref);
040                }
041                throw new RuntimeException("unknown algebraic fit type: " + type);
042        }
043
044        /**
045         * Returns a vector of algebraic ellipse parameters: {@code (a,b,c,d,e,f)}, representing the ellipse by
046         * {@code a x^2 + b x y + c y^2 + d x + e y + f = 0}.
047         *
048         * @return the ellipse parameters
049         */
050        public double[] getParameters();
051
052        /**
053         * Returns a algebraic ellipse constructed from the estimated ellipse parameters.
054         *
055         * @return an {@link AlgebraicEllipse} instance
056         */
057        public default AlgebraicEllipse getEllipse() {
058                double[] p = getParameters();
059                return (p == null) ? null : new AlgebraicEllipse(p);
060        }
061        
062        public default boolean isEllipse() {
063                double[] p = getParameters();   // p = (a,b,c,d,e,f)
064                if (p == null) {
065                        return false;
066                }
067                else {
068                        double a = p[0];
069                        double b = p[1];
070                        double c = p[2];
071                        return (sqr(b) - 4*a*c) < 0;
072                }
073        }
074
075        /**
076         * <p>
077         * Returns a 6x6 matrix (U) to convert ellipse parameters qr = (ar,...,fr) calculated for data centered at (xr, yr)
078         * to ellipse parameters q = (a,...,f) for the original non-centered data: q = U * qr. See [1, Sec. 11.2.1] (e.g.,
079         * Alg. 11.6) for additional information.
080         * </p>
081         * <p>
082         * [1] W. Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic Introduction</em>, 3rd ed,
083         * Springer (2022).
084         * </p>
085         *
086         * @param xr data reference point (x)
087         * @param yr data reference point (y)
088         * @return the offset correction matrix (U)
089         */
090        public default RealMatrix getDataOffsetCorrectionMatrix(double xr, double yr) {
091                return Matrix.makeRealMatrix(6, 6,
092                                1,       0,     0,        0,   0,  0 ,
093                                0,       1,     0,        0,   0,  0 ,
094                                0,       0,     1,        0,   0,  0 ,
095                           -2*xr,   -yr,    0,        1,   0,  0 ,
096                                0,      -xr,   -2*yr,     0,   1,  0 ,
097                                sqr(xr), xr*yr, sqr(yr), -xr, -yr, 1 );
098        }
099        
100        // ----------------------------------------------------------------------
101        
102        public static void main(String[] args) {
103                Pnt2d p1 = Pnt2d.from(31, 90);
104                Pnt2d p2 = Pnt2d.from(79, 25);
105                Pnt2d p3 = Pnt2d.from(159, 31);
106                Pnt2d p4 = Pnt2d.from(174, 55);
107                Pnt2d p5 = Pnt2d.from(173, 84);
108                Pnt2d p6 = Pnt2d.from(135, 114);
109                
110                Pnt2d[] pts = {p1, p2, p3, p4, p5, p6};
111                Pnt2d xr = Pnt2d.from(100,100);
112                
113                for (FitType type : FitType.values()) {
114                        EllipseFitAlgebraic fit = EllipseFitAlgebraic.getFit(type, pts, xr);
115                        System.out.println(type);
116                        if (fit != null) {
117                                System.out.println(fit.getEllipse());
118                        }
119                }
120                
121        }
122
123
124}