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.math.Arithmetic; 013import imagingbook.common.math.Matrix; 014 015import java.util.Locale; 016 017import static imagingbook.common.math.Arithmetic.sqr; 018import static imagingbook.common.math.Matrix.multiply; 019import static imagingbook.common.math.Matrix.normL2; 020import static java.lang.Math.cos; 021import static java.lang.Math.sin; 022import static java.lang.Math.sqrt; 023 024/** 025 * <p> 026 * Represents an algebraic ellipse with the implicit equation A x^2 + B x y + C y^2 + D x + E y + F = 0. Parameters A, 027 * ..., F are normalized such that B^2 - 4 A C = -1. Instances are immutable. See Secs. 11.2.1 and F.3.1 for details. 028 * </p> 029 * <p> 030 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 3rd ed, Springer 031 * (2022). 032 * </p> 033 * 034 * @author WB 035 * @version 2022/11/17 036 */ 037public class AlgebraicEllipse { 038 039 /** Ellipse parameters. */ 040 private final double A, B, C, D, E, F; 041 042 /** 043 * Constructor. Creates a {@link AlgebraicEllipse} instance by normalizing the supplied parameters [A,...,F] such 044 * that d = B^2 - 4 A C = -1. Throws an exception if d is non-negative. 045 * 046 * @param A ellipse parameter A 047 * @param B ellipse parameter B 048 * @param C ellipse parameter C 049 * @param D ellipse parameter D 050 * @param E ellipse parameter E 051 * @param F ellipse parameter F 052 */ 053 public AlgebraicEllipse(double A, double B, double C, double D, double E, double F) { 054 double d = sqr(B) - 4 * A * C; 055 if (d > -Arithmetic.EPSILON_DOUBLE) { // discriminant (d) must be negative! 056 throw new IllegalArgumentException("illegal ellipse parameters, non-negative discriminant " + d); 057 } 058 // d < 0 059 // normalize to B^2 - 4 A C = -1 060 double s = 1 / sqrt(-d); 061 if (A >= 0) { 062 this.A = A * s; 063 this.B = B * s; 064 this.C = C * s; 065 this.D = D * s; 066 this.E = E * s; 067 this.F = F * s; 068 } 069 else { 070 this.A = -A * s; 071 this.B = -B * s; 072 this.C = -C * s; 073 this.D = -D * s; 074 this.E = -E * s; 075 this.F = -F * s; 076 } 077 } 078 079 /** 080 * Constructor. Creates a {@link AlgebraicEllipse} instance from the specified parameter vector [A,...,F]. 081 * 082 * @param p algebraic ellipse parameters 083 */ 084 public AlgebraicEllipse(double[] p) { 085 this(p[0], p[1], p[2], p[3], p[4], p[5]); 086 } 087 088 /** 089 * Constructor. Creates a {@link AlgebraicEllipse} instance from a {@link GeometricEllipse}. 090 * 091 * @param ge a {@link GeometricEllipse} 092 */ 093 public AlgebraicEllipse(GeometricEllipse ge) { 094 this(getAlgebraicEllipseParameters(ge)); 095 } 096 097 private static double[] getAlgebraicEllipseParameters(GeometricEllipse ge) { 098 double ra = ge.ra; 099 double rb = ge.rb; 100 double xc = ge.xc; 101 double yc = ge.yc; 102 double theta = ge.theta; 103 104 double cosT = cos(theta); 105 double sinT = sin(theta); 106 107 double A = sqr(ra * sinT) + sqr(rb * cosT); 108 double B = 2 * (sqr(rb) - sqr(ra)) * sinT * cosT; 109 double C = sqr(ra * cosT) + sqr(rb * sinT); 110 double D = -2 * A * xc - B * yc; 111 double E = -2 * C * yc - B * xc; 112 double F = A * sqr(xc) + B * xc * yc + C * sqr(yc) - sqr(ra * rb); 113 return new double[] {A, B, C, D, E, F}; 114 } 115 116 /** 117 * Return a vector of parameters for this ellipse. The length of the vector and the meaning of the parameters 118 * depends on the concrete ellipse type. 119 * 120 * @return a vector of parameters [A, B, C, D, E, F] 121 */ 122 public double[] getParameters() { 123 return new double[] {A, B, C, D, E, F}; 124 } 125 126 public double getAlgebraicDistance(Pnt2d p) { 127 double x = p.getX(); 128 double y = p.getY(); 129 return A * sqr(x) + B * x * y + C * sqr(y) + D * x + E * y + F; 130 } 131 132 public double getSampsonDistance(Pnt2d p) { 133 double x = p.getX(); 134 double y = p.getY(); 135 double ad = A * sqr(x) + B * x * y + C * sqr(y) + D * x + E * y + F; // algebraic distance 136 double nGrad = Math.hypot(2*A*x + B*y + D, 2*C*y + B*x + E); 137 return ad / nGrad; 138 } 139 140 public double getGoncharovaDistance(Pnt2d p) { 141 double x = p.getX(); 142 double y = p.getY(); 143 double G = A * sqr(x) + B * x * y + C * sqr(y) + D * x + E * y + F; // algebraic distance 144 145 double[] gV = {2*A*x + B*y + D, 2*C*y + B*x + E}; 146 double ngV = normL2(gV); 147 double sd = G / ngV; 148 149 double[][] H = {{ 2*A, B }, { B, 2*C }}; 150 double gd = sd * (1 + G * Matrix.dotProduct(multiply(gV, H), gV) / (2 * Math.pow(ngV, 4))); 151 return gd; 152 } 153 154 // -------------------------------------------------------------------------- 155 156 @Override 157 public boolean equals(Object other) { 158 if (this == other) { 159 return true; 160 } 161 if (!(other instanceof AlgebraicEllipse)) { 162 return false; 163 } 164 return this.equals((AlgebraicEllipse) other, Arithmetic.EPSILON_DOUBLE); 165 } 166 167 public boolean equals(AlgebraicEllipse other, double tolerance) { 168 return 169 Arithmetic.equals(A, other.A, tolerance) && 170 Arithmetic.equals(B, other.B, tolerance) && 171 Arithmetic.equals(C, other.C, tolerance) && 172 Arithmetic.equals(D, other.D, tolerance) && 173 Arithmetic.equals(E, other.E, tolerance) && 174 Arithmetic.equals(F, other.F, tolerance) ; 175 } 176 177 // -------------------------------------------------------------------------- 178 179 public AlgebraicEllipse duplicate() { 180 return new AlgebraicEllipse(this.getParameters()); 181 } 182 183 @Override 184 public String toString() { 185 return String.format(Locale.US, "%s %s", 186 this.getClass().getSimpleName(), 187 Matrix.toString(getParameters())); 188 } 189 190 // -------------------------------------------------------------------------- 191 192// public static void main(String[] args) { 193// PrintPrecision.set(9); 194// GeometricEllipse eg = new GeometricEllipse(150, 50, 150, 150, 0.5); 195// System.out.println("eg = " + eg); 196// AlgebraicEllipse ea = new AlgebraicEllipse(eg); 197// System.out.println("ea = " + ea); 198// 199// Pnt2d p = Pnt2d.from(300, 300); 200// System.out.println("p = " + p); 201// Pnt2d pc = eg.getClosestPoint(p); 202// System.out.println("pc = " + pc); 203// 204// System.out.println("Sampson dist(p) = " + ea.getSampsonDistance(p)); 205// System.out.println("Gonch dist(p) = " + ea.getGoncharovaDistance(p)); 206// 207// } 208 209}