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.ellipse.project;
010
011import imagingbook.common.geometry.ellipse.GeometricEllipse;
012
013import static imagingbook.common.math.Arithmetic.sqr;
014import static java.lang.Math.abs;
015import static java.lang.Math.max;
016import static java.lang.Math.pow;
017
018/**
019 * <p>
020 * Calculates the closest point on the ellipse for a given 2D point inside or outside the ellipse, using orthogonal
021 * projection of points onto the ellipse. This is a robust algorithm based on [1]. See Sec.11.2.2 (Alg. 11.9) of [2] for
022 * details. This version uses the Newton-method for root finding, which is quick but may fail to return a valid result
023 * if the target point is close to the x- or y-axis. See {@link OrthogonalEllipseProjector} for a robust solution or
024 * {@link ConfocalConicEllipseProjector} for an approximate but non-iterative (i.e., fast) alternative.
025 * </p>
026 * <p>
027 * [1] D. Eberly: "Distance from a point to an ellipse, an ellipsoid, or a hyperellipsoid", Technical Report, Geometric
028 * Tools, www.geometrictools.com, Redmont, WA (June 2013). <br> [2] W. Burger, M.J. Burge, <em>Digital Image Processing
029 * &ndash; An Algorithmic Introduction</em>, 3rd ed, Springer (2022).
030 * </p>
031 *
032 * @author WB
033 * @version 2022/04/09
034 * @see OrthogonalEllipseProjector
035 * @see ConfocalConicEllipseProjector
036 */
037public class OrthogonalEllipseProjectorNewton extends EllipseProjector {
038        
039        private final static double NewtonMinStep = 1e-6;
040        private final static int MaxIterations = 100;
041        
042        private final double ra, rb, ra2, rb2;
043        private int lastIterationCount = 0;     // number of root-finding iterations performed in last projection
044        
045        public OrthogonalEllipseProjectorNewton(GeometricEllipse ellipse) {
046                super(ellipse);
047                this.ra = ellipse.ra;
048                this.rb = ellipse.rb;
049                this.ra2 = sqr(ra);
050                this.rb2 = sqr(rb);
051        }
052        
053        @Override
054        protected double[] projectCanonical(double[] u1) {
055                // coordinates of p (mapped to first quadrant)
056                final double u = u1[0]; 
057                final double v = u1[1]; 
058                
059                double[] ub = null;     // the unknown ellipse point
060
061                if (u + v < 1e-6) {     // (u,v) is very close to the ellipse center; u,v >= 0
062                        ub = new double[] {0, rb};
063                }
064                else {                                          
065                        double t = max(ra * u - ra2, rb * v - rb2);
066                        double gprev = Double.POSITIVE_INFINITY;
067                        double deltaT, deltaG;
068                        int k = 0;
069                        do {
070                                k = k + 1;
071                                double g  = sqr((ra * u) / (t + ra2)) + sqr((rb * v) / (t + rb2)) - 1;
072                                double dg = 2 * (sqr(ra * u) / pow(t + ra2, 3) + sqr(rb * v) / pow(t + rb2, 3));
073                                deltaT = g / dg;
074                                t = t + deltaT;                         // Newton iteration
075                                
076                                // in rare cases g(t) is very flat and checking deltaT is not enough for convergence!
077                                deltaG = g - gprev;                     // change of g value
078                                gprev = g;      
079                                
080                        }  while(abs(deltaT) > NewtonMinStep && abs(deltaG) > NewtonMinStep && k < MaxIterations);
081                        
082                        lastIterationCount = k;         // remember iteration count
083                        
084                        if (k >= MaxIterations) {
085                                throw new RuntimeException("max. mumber of iterations exceeded");
086                        }
087                        
088                        ub = new double[] {ra2 * u / (t + ra2), rb2 * v / (t + rb2)};
089                }
090                return ub;
091        }
092
093        // for statistics only
094        
095        public int getLastIterationCount() {
096                return this.lastIterationCount;
097        }
098
099}