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; 010 011import imagingbook.common.geometry.basic.Pnt2d; 012import imagingbook.common.geometry.basic.Primitive2d; 013import imagingbook.common.geometry.ellipse.project.OrthogonalEllipseProjector; 014import imagingbook.common.geometry.shape.ShapeProducer; 015import imagingbook.common.math.Arithmetic; 016 017import java.awt.Shape; 018import java.awt.geom.AffineTransform; 019import java.awt.geom.Ellipse2D; 020import java.awt.geom.Path2D; 021import java.util.Locale; 022 023import static imagingbook.common.math.Arithmetic.sqr; 024import static java.lang.Math.PI; 025import static java.lang.Math.cos; 026import static java.lang.Math.sin; 027import static java.lang.Math.sqrt; 028 029/** 030 * <p> 031 * Represents an ellipse with major axis length ra, minor axis length rb, center point (xc, yc), and orientation theta. 032 * Instances are immutable. See Secs. 11.2.2 and F.3.1 for details. 033 * </p> 034 * <p> 035 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 3rd ed, Springer 036 * (2022). 037 * </p> 038 * 039 * @author WB 040 * @version 2022/11/17 041 */ 042public class GeometricEllipse implements ShapeProducer, Primitive2d { 043 044 /** 045 * Ellipse parameters. 046 */ 047 public final double ra, rb, xc, yc, theta; 048 private final OrthogonalEllipseProjector projector; 049 050 /** 051 * Constructor. Axis lengths may be exchanged to ensure ra ≥ rb. 052 * 053 * @param ra major axis length 054 * @param rb minor axis length 055 * @param xc center point (x) 056 * @param yc center point (y) 057 * @param theta orientation 058 */ 059 public GeometricEllipse(double ra, double rb, double xc, double yc, double theta) { 060 this.xc = xc; 061 this.yc = yc; 062 if (ra >= rb) { // make sure ra is always the longer axis! 063 this.ra = ra; 064 this.rb = rb; 065 this.theta = Arithmetic.mod(theta, PI); // theta = 0,...,pi always 066 } 067 else { 068 this.ra = rb; // swap ra/rb 069 this.rb = ra; 070 this.theta = Arithmetic.mod(theta + PI/2, PI); // theta = 0,...,pi always 071 } 072 this.projector = new OrthogonalEllipseProjector(this); 073 074 if (this.ra < this.rb) { 075 throw new RuntimeException(this.toString()); 076 } 077 } 078 079 /** 080 * Constructor, short version for axis-aligned ellipses. Axis lengths may be exchanged to ensure ra ≥ rb. 081 * 082 * @param ra major axis length 083 * @param rb minor axis length 084 * @param xc center point (x) 085 * @param yc center point (y) 086 */ 087 public GeometricEllipse(double ra, double rb, double xc, double yc) { 088 this(ra, rb, xc, yc, 0.0); 089 } 090 091 /** 092 * Constructor. 093 * 094 * @param p parameter vector [ra, rb, xc, yc, theta]. 095 * @see #getParameters() 096 */ 097 public GeometricEllipse(double[] p) { 098 this(p[0], p[1], p[2], p[3], p[4]); 099 } 100 101 /** 102 * Constructor. Creates a new {@link GeometricEllipse} from a {@link AlgebraicEllipse} instance. 103 * 104 * @param ae a {@link AlgebraicEllipse} instance 105 */ 106 public GeometricEllipse(AlgebraicEllipse ae) { 107 this(getGeometricEllipseParameters(ae)); 108 } 109 110 /** 111 * Calculates and returns the geometric ellipse parameters from a given algebraic ellipse (see Eqns. 19-23 at 112 * <a href="http://mathworld.wolfram.com/Ellipse.html">http://mathworld.wolfram.com/Ellipse.html</a>). 113 * 114 * @param ae a {@linkplain AlgebraicEllipse} instance with parameters (A,...,F) 115 * @return the geometric ellipse parameters (ra, rb, xc, yc, theta) 116 * @see AlgebraicEllipse 117 */ 118 public static double[] getGeometricEllipseParameters(AlgebraicEllipse ae) { 119 // see Eq. 19-23 at http://mathworld.wolfram.com/Ellipse.html 120 double[] params = ae.getParameters(); 121 final double A = params[0]; // ae.A; 122 final double B = params[1]; // ae.B; 123 final double C = params[2]; // ae.C; 124 final double D = params[3]; // ae.D; 125 final double E = params[4]; // ae.E; 126 final double F = params[5]; // ae.F; 127 final double d = sqr(B) - 4*A*C; 128 129 if (d >= 0) { 130 throw new IllegalArgumentException("B^2 - 4AC must be negative for an ellipse"); 131 } 132 133 final double q = sqrt(sqr(A-C) + sqr(B)); 134 final double s = 2*(A*sqr(E) + C*sqr(D) + F*sqr(B) - B*D*E - 4*A*C*F); 135 136 double xc = (2 * C * D - B * E) / d; 137 double yc = (2 * A * E - B * D) / d; 138 double ra = sqrt(s / (d * (-A - C + q))); 139 double rb = sqrt(s / (d * (-A - C - q))); 140 double theta = 0.5 * Math.atan2(-B, C - A); // theta = -pi/2,...,+pi/2 141 142 return (ra >= rb) ? // make sure ra >= rb (ra is the major axis) 143 new double[] {ra, rb, xc, yc, theta} : 144 new double[] {rb, ra, xc, yc, theta + PI/2} ; 145 } 146 147 // --------------------------------------- 148 149 /** 150 * Returns a vector of parameters for this ellipse. The length of the vector and the meaning of the parameters 151 * depends on the concrete ellipse type. 152 * 153 * @return a vector of parameters [ra, rb, xc, yc, theta] 154 */ 155 public double[] getParameters() { 156 return new double[] {ra, rb, xc, yc, theta}; 157 } 158 159 public Pnt2d getCenter() { 160 return Pnt2d.from(xc, yc); 161 } 162 163 public double getArea() { 164 return this.ra * this.rb * Math.PI; 165 } 166 167 public double getAlgebraicDistance(Pnt2d p) { 168 return new AlgebraicEllipse(this).getAlgebraicDistance(p); 169 } 170 171 /** 172 * Returns the ellipse point closest to the specified point. To perform this calculation for multiple points on the 173 * same ellipse use {@link OrthogonalEllipseProjector}. 174 * 175 * @param pnt some 2D point 176 * @return the closest ("contact") point on the ellipse 177 * @see OrthogonalEllipseProjector 178 */ 179 public Pnt2d getClosestPoint(Pnt2d pnt) { 180 return projector.project(pnt); 181 } 182 183 public double[] getClosestPoint(double[] pnt) { 184 return projector.project(pnt); 185 } 186 187 // --------------------------------------------------------------------------------- 188 189 /** 190 * Returns the mean squared error between this ellipse and a set of 2D points. 191 * 192 * @param points a set of sample points (usually the points used for fitting) 193 * @return the mean squared error 194 */ 195 public double getMeanSquareError(Pnt2d[] points) { 196 double sum2 = 0; 197 for (Pnt2d p : points) { 198 sum2 = sum2 + p.distanceSq(getClosestPoint(p)); 199 } 200 return sum2 / points.length; 201 } 202 203 /** 204 * Returns the closest (geometric) distance of the specified point to this ellipse. 205 * 206 * @param p some 2D point 207 * @return the distance to the closest point on the ellipse 208 */ 209 @Override 210 public double getDistance(Pnt2d p) { 211 return p.distance(getClosestPoint(p)); 212 } 213 214 // ------------------------------------------------------------------------------------------ 215 216 @Override 217 public Shape getShape(double scale) { 218 return this.getOuterShape(); 219 } 220 221 @Override 222 public Shape[] getShapes(double scale) { 223 return new Shape[] { 224 getOuterShape(), 225 getCenterShape(scale), 226 getAxisShape() 227 }; 228 } 229 230 private Shape getOuterShape() { 231 Ellipse2D oval = new Ellipse2D.Double(-ra, -rb, 2 * ra, 2 * rb); 232 AffineTransform trans = new AffineTransform(); 233 trans.translate(xc, yc); 234 trans.rotate(theta); 235 return trans.createTransformedShape(oval); 236 } 237 238 private Shape getCenterShape(double radius) { 239 double dxa = 2 * radius * cos(theta); // major axis is drawn longer 240 double dya = 2 * radius * sin(theta); 241 double dxb = 1 * radius * cos(theta + PI/2); 242 double dyb = 1 * radius * sin(theta + PI/2); 243 Path2D path = new Path2D.Double(); 244 path.moveTo(xc - dxa, yc - dya); 245 path.lineTo(xc + dxa, yc + dya); 246 path.moveTo(xc - dxb, yc - dyb); 247 path.lineTo(xc + dxb, yc + dyb); 248 return path; 249 } 250 251 private Shape getAxisShape() { 252 double dxa = ra * cos(theta); 253 double dya = ra * sin(theta); 254 double dxb = rb * cos(theta + PI/2); 255 double dyb = rb * sin(theta + PI/2); 256 Path2D path = new Path2D.Double(); 257 path.moveTo(xc - dxa, yc - dya); 258 path.lineTo(xc + dxa, yc + dya); 259 path.moveTo(xc - dxb, yc - dyb); 260 path.lineTo(xc + dxb, yc + dyb); 261 return path; 262 } 263 264 @Deprecated // used for book figures only 265 public Pnt2d getLeftAxisPoint() { 266 return Pnt2d.from(xc - ra * cos(theta), yc - ra * sin(theta)); 267 } 268 269 // ------------------------------- 270 271 public GeometricEllipse duplicate() { 272 return new GeometricEllipse(this.getParameters()); 273 } 274 275 @Override 276 public boolean equals(Object other) { 277 if (this == other) { 278 return true; 279 } 280 if (!(other instanceof GeometricEllipse)) { 281 return false; 282 } 283 return this.equals((GeometricEllipse) other, Arithmetic.EPSILON_DOUBLE); 284 } 285 286 public boolean equals(GeometricEllipse other, double tolerance) { 287 return 288 Arithmetic.equals(xc, other.xc, tolerance) && 289 Arithmetic.equals(yc, other.yc, tolerance) && 290 Arithmetic.equals(ra, other.ra, tolerance) && 291 Arithmetic.equals(rb, other.rb, tolerance) && 292 Arithmetic.equals(Math.cos(theta), Math.cos(other.theta), tolerance) && 293 Arithmetic.equals(Math.sin(theta), Math.sin(other.theta), tolerance) ; 294 } 295 296 @Override 297 public String toString() { 298 return String.format(Locale.US, "%s [ra=%.8f, rb=%.8f, xc=%.8f, yc=%.8f, theta=%.8f]", 299 this.getClass().getSimpleName(), ra, rb, xc, yc, theta); 300 } 301 302}