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;
012import imagingbook.common.math.Arithmetic;
013
014import static imagingbook.common.math.Arithmetic.isZero;
015import static imagingbook.common.math.Arithmetic.sqr;
016
017/**
018 * <p>
019 * Calculates the closest point on the ellipse for a given 2D point inside or outside the ellipse, using orthogonal
020 * 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
021 * details. In contrast to the Newton-based algorithm (see {@link OrthogonalEllipseProjectorNewton}, this version uses
022 * the bisection method for root finding and returns valid results for points close to the x- and y-axis but requires
023 * significantly more iterations to converge.
024 * </p>
025 * <p>
026 * [1] D. Eberly: "Distance from a point to an ellipse, an ellipsoid, or a hyperellipsoid", Technical Report, Geometric
027 * Tools, www.geometrictools.com, Redmont, WA (June 2013). <br> [2] W. Burger, M.J. Burge, <em>Digital Image Processing
028 * &ndash; An Algorithmic Introduction</em>, 3rd ed, Springer (2022).
029 * </p>
030 *
031 * @author WB
032 * @version 2022/04/09
033 * @see OrthogonalEllipseProjectorNewton
034 * @see ConfocalConicEllipseProjector
035 */
036public class OrthogonalEllipseProjector extends EllipseProjector {
037        
038        private static final int MaxIterations = 150;
039        private static final double Epsilon = 1e-6;
040        
041        private final double ra, rb, rab;
042        private int lastIterationCount = 0;     // number of root-finding iterations performed in last projection
043        
044        public OrthogonalEllipseProjector(GeometricEllipse ellipse) {
045                super(ellipse);
046                this.ra = ellipse.ra;
047                this.rb = ellipse.rb;
048                this.rab = sqr(ra / rb);        // ratio of squared axes lengths
049        }
050        
051        @Override
052        protected double[] projectCanonical(double[] u1) {
053                // coordinates of p (mapped to first quadrant of canonical coordinates)
054                final double u = u1[0]; // u,v are both positive
055                final double v = u1[1];
056                
057                double[] ub = null;     // the ellipse contact point (in canonical coordinates)
058                lastIterationCount = 0;
059                
060                if (v > 0) {
061                        if (u > 0) {
062                                double uu = u / ra;
063                                double vv = v / rb;
064                                double ginit = sqr(uu) + sqr(vv) - 1;
065                                if (!isZero(ginit)) {
066                                        double s = getRoot(uu, vv, ginit);
067                                        ub = new double[] {rab * u / (s + rab), v / (s + 1)};
068                                }
069                                else {
070                                        ub = new double[] {u, v};
071                                }
072                        }
073                        else {  // u = 0
074                                ub = new double[] {0, rb};
075                        }
076                }       
077                else {  // v = 0
078                        double numer0 = ra * u;
079                        double denom0 = sqr(ra) - sqr(rb);
080                        if (numer0 < denom0) {
081                                double xde0 = numer0 / denom0;
082                                ub = new double[] {ra * xde0, rb * Math.sqrt(1 - sqr(xde0))};
083                        }
084                        else {
085                                ub = new double[] {ra, 0};
086                        }
087                }
088                return ub;
089        }
090
091        // Find the root of function
092        // G(s) = [(rab * uu) / (s + rab)]^2 + [vv / (s + 1)]^2 - 1
093        // using the bisection method.
094        private double getRoot(double uu, double vv, double g0) {
095                double s0 = vv - 1;
096                double s1 = (g0 < 0) ? 0 : Math.hypot(rab * uu, vv) - 1;
097                double s = 0;
098                
099                int i;
100                for (i = 0; i < MaxIterations; i++) {
101                        s = (s0 + s1) / 2;
102                        if (Arithmetic.equals(s, s0) || Arithmetic.equals(s, s1)) {
103                                break;
104                        }
105                        double g = sqr((rab * uu)/(s + rab)) + sqr(vv/(s + 1)) - 1; // = G(s)
106                        if (g > Epsilon) {                      // G(s) is positive
107                                s0 = s;
108                        }
109                        else if (g < -Epsilon) {        // G(s) is negative
110                                s1 = s;
111                        }
112                        else {                                          // G(s) ~ 0 (done)
113                                break;
114                        }
115                }
116                if (i >= MaxIterations) {
117                        throw new RuntimeException(this.getClass().getSimpleName() + 
118                                        ": max. iteration count exceeded");
119                }
120                lastIterationCount = i;
121                return s;
122        }
123        
124        // for statistics only
125        public int getLastIterationCount() {
126                return this.lastIterationCount;
127        }
128
129}