1 /*
2 * Copyright (C) 2017 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.geometry;
17
18 import com.irurueta.algebra.AlgebraException;
19 import com.irurueta.algebra.Matrix;
20 import com.irurueta.algebra.SingularValueDecomposer;
21
22 import java.io.Serializable;
23
24 /**
25 * This class defines an ellipse.
26 * This class uses formulas as defined at:
27 * <a href="https://en.wikipedia.org/wiki/Ellipse">https://en.wikipedia.org/wiki/Ellipse</a>
28 */
29 @SuppressWarnings("DuplicatedCode")
30 public class Ellipse implements Serializable {
31
32 /**
33 * Constant defining default threshold value used when none is provided.
34 */
35 public static final double DEFAULT_THRESHOLD = 1e-9;
36
37 /**
38 * Constant defining minimum allowed threshold.
39 */
40 public static final double MIN_THRESHOLD = 0.0;
41
42 /**
43 * Center of ellipse.
44 */
45 private Point2D center;
46
47 /**
48 * Semi-major axis length (a).
49 */
50 private double semiMajorAxis;
51
52 /**
53 * Semi-minor axis length (b).
54 */
55 private double semiMinorAxis;
56
57 /**
58 * Rotation angle.
59 */
60 private double rotationAngle;
61
62 /**
63 * Empty constructor.
64 * Creates an ellipse equal to a circle located at space origin (0,0) with
65 * radius 1.0.
66 */
67 public Ellipse() {
68 center = Point2D.create();
69 semiMajorAxis = semiMinorAxis = 1.0;
70 rotationAngle = 0.0;
71 }
72
73 /**
74 * Constructor.
75 *
76 * @param center center of ellipse.
77 * @param semiMajorAxis semi-major axis length.
78 * @param semiMinorAxis semi-minor axis length.
79 * @param rotationAngle rotation angle expressed in radians.
80 */
81 public Ellipse(final Point2D center, final double semiMajorAxis, final double semiMinorAxis,
82 final double rotationAngle) {
83 setCenterAxesAndRotation(center, semiMajorAxis, semiMinorAxis, rotationAngle);
84 }
85
86 /**
87 * Constructor.
88 *
89 * @param center center of ellipse.
90 * @param semiMajorAxis semi-major axis length.
91 * @param semiMinorAxis semi-minor axis length.
92 * @param rotation 2D rotation.
93 */
94 public Ellipse(final Point2D center, final double semiMajorAxis, final double semiMinorAxis,
95 final Rotation2D rotation) {
96 setCenterAxesAndRotation(center, semiMajorAxis, semiMinorAxis, rotation);
97 }
98
99 /**
100 * Constructor from 2 points, ellipse center and rotation.
101 *
102 * @param point1 1st point.
103 * @param point2 2nd point.
104 * @param center center of ellipse.
105 * @param theta rotation angle expressed in radians.
106 * @throws ColinearPointsException if points are in a co-linear or degenerate
107 * configuration.
108 */
109 public Ellipse(final Point2D point1, final Point2D point2, final Point2D center, final double theta)
110 throws ColinearPointsException {
111 setParametersFromPointsCenterAndRotation(point1, point2, center, theta);
112 }
113
114 /**
115 * Constructor from 5 points.
116 *
117 * @param point1 1st point.
118 * @param point2 2nd point.
119 * @param point3 3rd point.
120 * @param point4 4th point.
121 * @param point5 5th point.
122 * @throws ColinearPointsException if points are in a co-linear or degenerate
123 * configuration.
124 */
125 public Ellipse(final Point2D point1, final Point2D point2, final Point2D point3, final Point2D point4,
126 final Point2D point5) throws ColinearPointsException {
127 setParametersFromPoints(point1, point2, point3, point4, point5);
128 }
129
130 /**
131 * Constructor from 5 points.
132 *
133 * @param point1 1st point.
134 * @param point2 2nd point.
135 * @param point3 3rd point.
136 * @param point4 4th point.
137 * @param point5 5th point.
138 * @param threshold threshold to determine whether points form an ellipse.
139 * This is usually a very small value
140 * @throws ColinearPointsException if points are in a co-linear or degenerate
141 * configuration.
142 */
143 public Ellipse(final Point2D point1, final Point2D point2, final Point2D point3, final Point2D point4,
144 final Point2D point5, double threshold) throws ColinearPointsException {
145 setParametersFromPoints(point1, point2, point3, point4, point5, threshold);
146 }
147
148 /**
149 * Constructor setting parameters of canonical equation of an ellipse, which
150 * is:
151 * a*x^2 + b*x*y + c*y^2 + d*x + e*y + f = 0
152 *
153 * @param a a parameter.
154 * @param b b parameter.
155 * @param c c parameter.
156 * @param d d parameter.
157 * @param e e parameter.
158 * @param f f parameter.
159 * @throws IllegalArgumentException if parameters do not follow
160 * b^2 - 4*a*c < 0.0
161 */
162 public Ellipse(final double a, final double b, final double c, final double d, final double e, final double f) {
163 setParameters(a, b, c, d, e, f);
164 }
165
166 /**
167 * Constructor setting parameters of canonical equation of an ellipse, which
168 * is:
169 * a*x^2 + b*x*y + c*y^2 + d*x + e*y + f = 0
170 *
171 * @param a a parameter.
172 * @param b b parameter.
173 * @param c c parameter.
174 * @param d d parameter.
175 * @param e e parameter.
176 * @param f f parameter.
177 * @param threshold threshold to determine whether parameters are valid due
178 * to machine precision.
179 * @throws IllegalArgumentException if parameters do not follow
180 * b^2 - 4*A*c < threshold
181 */
182 public Ellipse(final double a, final double b, final double c, final double d, final double e, final double f,
183 final double threshold) {
184 setParameters(a, b, c, d, e, f, threshold);
185 }
186
187 /**
188 * Constructor.
189 *
190 * @param conic conic to build ellipse from.
191 * @throws IllegalArgumentException if provided conic is not an ellipse.
192 */
193 public Ellipse(final Conic conic) {
194 setFromConic(conic);
195 }
196
197 /**
198 * Constructor.
199 *
200 * @param circle a circle to set parameters from.
201 */
202 public Ellipse(final Circle circle) {
203 setFromCircle(circle);
204 }
205
206
207 /**
208 * Returns center of ellipse.
209 *
210 * @return center of ellipse.
211 */
212 public Point2D getCenter() {
213 return center;
214 }
215
216 /**
217 * Sets center of ellipse.
218 *
219 * @param center center of ellipse.
220 * @throws NullPointerException raised if provided center is null.
221 */
222 public void setCenter(final Point2D center) {
223 if (center == null) {
224 throw new NullPointerException();
225 }
226 this.center = center;
227 }
228
229 /**
230 * Gets semi-major axis length.
231 *
232 * @return semi-major axis length.
233 */
234 public double getSemiMajorAxis() {
235 return semiMajorAxis;
236 }
237
238 /**
239 * Sets semi-major axis length.
240 *
241 * @param semiMajorAxis semi-major axis length.
242 */
243 public void setSemiMajorAxis(final double semiMajorAxis) {
244 this.semiMajorAxis = semiMajorAxis;
245 }
246
247 /**
248 * Gets semi-minor axis length.
249 *
250 * @return semi-minor axis length.
251 */
252 public double getSemiMinorAxis() {
253 return semiMinorAxis;
254 }
255
256 /**
257 * Sets semi-minor axis length.
258 *
259 * @param semiMinorAxis semi-minor axis length.
260 */
261 public void setSemiMinorAxis(final double semiMinorAxis) {
262 this.semiMinorAxis = semiMinorAxis;
263 }
264
265 /**
266 * Gets rotation angle expressed in radians.
267 *
268 * @return rotation angle expressed in radians.
269 */
270 public double getRotationAngle() {
271 return rotationAngle;
272 }
273
274 /**
275 * Sets rotation angle expressed in radians.
276 *
277 * @param rotationAngle rotation angle expressed in radians.
278 */
279 public void setRotationAngle(final double rotationAngle) {
280 this.rotationAngle = rotationAngle;
281 }
282
283 /**
284 * Gets 2D rotation.
285 *
286 * @return 2D rotation.
287 */
288 public Rotation2D getRotation() {
289 return new Rotation2D(rotationAngle);
290 }
291
292 /**
293 * Sets 2D rotation.
294 *
295 * @param rotation 2D rotation to be set.
296 */
297 public void setRotation(final Rotation2D rotation) {
298 rotationAngle = rotation.getTheta();
299 }
300
301 /**
302 * Gets parameter A of canonical ellipse equation:
303 * A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
304 *
305 * @return parameter A.
306 */
307 public double getA() {
308 final var sint = Math.sin(rotationAngle);
309 final var cost = Math.cos(rotationAngle);
310
311 final var a = semiMajorAxis;
312 final var b = semiMinorAxis;
313
314 final var sint2 = sint * sint;
315 final var cost2 = cost * cost;
316
317 final var a2 = a * a;
318 final var b2 = b * b;
319
320 return a2 * sint2 + b2 * cost2;
321 }
322
323 /**
324 * Sets parameter A of canonical ellipse equation.
325 * A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
326 *
327 * @param a parameter A to be set.
328 */
329 public void setA(final double a) {
330 setParameters(a, getB(), getC(), getD(), getE(), getF());
331 }
332
333 /**
334 * Gets parameter B of canonical ellipse equation:
335 * A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
336 *
337 * @return parameter B.
338 */
339 public double getB() {
340 final var sint = Math.sin(rotationAngle);
341 final var cost = Math.cos(rotationAngle);
342 final var a = semiMajorAxis;
343 final var b = semiMinorAxis;
344 final var a2 = a * a;
345 final var b2 = b * b;
346
347 return 2.0 * (b2 - a2) * sint * cost;
348 }
349
350 /**
351 * Sets parameter B of canonical ellipse equation:
352 * A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
353 *
354 * @param b parameter B to be set.
355 */
356 public void setB(final double b) {
357 setParameters(getA(), b, getC(), getD(), getE(), getF());
358 }
359
360 /**
361 * Gets parameter C of canonical ellipse equation:
362 * A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
363 *
364 * @return parameter C.
365 */
366 public double getC() {
367 final var sint = Math.sin(rotationAngle);
368 final var cost = Math.cos(rotationAngle);
369
370 final var a = semiMajorAxis;
371 final var b = semiMinorAxis;
372
373 final var sint2 = sint * sint;
374 final var cost2 = cost * cost;
375
376 final var a2 = a * a;
377 final var b2 = b * b;
378
379 return a2 * cost2 + b2 * sint2;
380 }
381
382 /**
383 * Sets parameter C of canonical ellipse equation:
384 * A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
385 *
386 * @param c parameter C to be set.
387 */
388 public void setC(final double c) {
389 setParameters(getA(), getB(), c, getD(), getE(), getF());
390 }
391
392 /**
393 * Gets parameter D of canonical ellipse equation:
394 * A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
395 *
396 * @return parameter D.
397 */
398 public double getD() {
399 final var xc = center.getInhomX();
400 final var yc = center.getInhomY();
401
402 return -2.0 * getA() * xc - getB() * yc;
403 }
404
405 /**
406 * Sets parameter D of canonical ellipse equation:
407 * A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
408 *
409 * @param d parameter D to be set.
410 */
411 public void setD(final double d) {
412 setParameters(getA(), getB(), getC(), d, getE(), getF());
413 }
414
415 /**
416 * Gets parameter E of canonical ellipse equation:
417 * A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
418 *
419 * @return parameter E.
420 */
421 public double getE() {
422 final var xc = center.getInhomX();
423 final var yc = center.getInhomY();
424
425 return -getB() * xc - 2.0 * getC() * yc;
426 }
427
428 /**
429 * Sets parameter E of canonical ellipse equation:
430 * A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
431 *
432 * @param e parameter E to be set.
433 */
434 public void setE(final double e) {
435 setParameters(getA(), getB(), getC(), getD(), e, getF());
436 }
437
438 /**
439 * Gets parameter F of canonical ellipse equation:
440 * A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
441 *
442 * @return parameter F.
443 */
444 public double getF() {
445 final var xc = center.getInhomX();
446 final var yc = center.getInhomY();
447
448 final var a = semiMajorAxis;
449 final var b = semiMinorAxis;
450
451 final var xc2 = xc * xc;
452 final var yc2 = yc * yc;
453
454 final var a2 = a * a;
455 final var b2 = b * b;
456
457 return getA() * xc2 + getB() * xc * yc + getC() * yc2 - a2 * b2;
458 }
459
460 /**
461 * Sets parameter F of canonical ellipse equation.
462 * A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
463 *
464 * @param f parameter F to be set.
465 */
466 public void setF(final double f) {
467 setParameters(getA(), getB(), getC(), getD(), getE(), f);
468 }
469
470 /**
471 * Sets ellipse parameters.
472 *
473 * @param center center of ellipse.
474 * @param semiMajorAxis semi-major axis length.
475 * @param semiMinorAxis semi-minor axis length.
476 * @param rotationAngle rotation angle expressed in radians.
477 */
478 public final void setCenterAxesAndRotation(
479 final Point2D center, final double semiMajorAxis, final double semiMinorAxis, final double rotationAngle) {
480 this.center = center;
481 this.semiMajorAxis = Math.max(semiMajorAxis, semiMinorAxis);
482 this.semiMinorAxis = Math.min(semiMinorAxis, semiMajorAxis);
483 this.rotationAngle = rotationAngle;
484 }
485
486 /**
487 * Sets ellipse parameters.
488 *
489 * @param center center of ellipse.
490 * @param semiMajorAxis semi-major axis length.
491 * @param semiMinorAxis semi-minor axis length.
492 * @param rotation 2D rotation.
493 */
494 public final void setCenterAxesAndRotation(
495 final Point2D center, final double semiMajorAxis, final double semiMinorAxis, final Rotation2D rotation) {
496 setCenterAxesAndRotation(center, semiMajorAxis, semiMinorAxis, rotation.getTheta());
497 }
498
499 /**
500 * Sets parameters from 2 points, ellipse center and rotation.
501 *
502 * @param point1 1st point.
503 * @param point2 2nd point.
504 * @param center center of ellipse.
505 * @param theta rotation angle expressed in radians.
506 * @throws ColinearPointsException if points are in a co-linear or degenerate
507 * configuration.
508 */
509 public final void setParametersFromPointsCenterAndRotation(
510 final Point2D point1, final Point2D point2, final Point2D center, final double theta)
511 throws ColinearPointsException {
512 // unknowns: semi-major axis (a) and semi-minor axis (b)
513
514 // equation of an ellipse follows:
515 // A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
516
517 // where
518 // sint = sin(theta)
519 // cost = cos(theta)
520
521 // A = a^2*sint^2 + b^2*cost^2
522 // B = 2*(b^2 - a^2)*sint*cost
523 // C = a^2*cost^2 + b^2*sint^2
524 // D = -2*A*xc - B*yc = -2*(a^2*sint^2 + b^2*cost^2)*xc - 2*(b^2 - a^2)*sint*cost*yc
525 // E = -B*xc - 2*C*yc = -2*(b^2 - a^2)*sint*cost*xc - 2*(a^2*cost^2 + b^2*sint^2)*yc
526 // F = A*xc^2 + B*xc*yc + C*yc^2 - a^2*b^2 = (a^2*sint^2 + b^2*cost^2)*xc^2 + 2*(b^2 - a^2)*sint*cost*xc*yc + (a^2*cost^2 + b^2*sint^2)*yc^2 - a^2*b^2
527
528 // (a^2*sint^2 + b^2*cost^2)*x^2 +
529 // 2*(b^2 - a^2)*sint*cost*x*y +
530 // (a^2*cost^2 + b^2*sint^2)*y^2 +
531 // (-2*(a^2*sint^2 + b^2*cost^2)*xc - 2*(b^2 - a^2)*sint*cost*yc)*x +
532 // (-2*(b^2 - a^2)*sint*cost*xc - 2*(a^2*cost^2 + b^2*sint^2)*yc)*y +
533 // (a^2*sint^2 + b^2*cost^2)*xc^2 + 2*(b^2 - a^2)*sint*cost*xc*yc + (a^2*cost^2 + b^2*sint^2)*yc^2 - a^2*b^2 = 0
534
535 // a^2*sint^2*x^2 + b^2*cost^2*x^2 +
536 // 2*b^2*sint*cost*x*y - 2*a^2*sint*cost*x*y +
537 // a^2*cost^2*y^2 + b^2*sint^2*y^2 +
538 // -2*a^2*sint^2*xc*x -2*b^2*cost^2*xc*x - 2*b^2*sint*cost*yc*x + 2*a^2*sint*cost*yc*x +
539 // -2*b^2*sint*cost*xc*y + 2*a^2*sint*cost*xc*y -2*a^2*cost^2*yc*y -2*b^2*sint^2*yc*y +
540 // a^2*sint^2*xc^2 + b^2*cost^2*xc^2 + 2*b^2*sint*cost*xc*yc -2*a^2*sint*cost*xc*yc + a^2*cost^2*yc^2 + b^2*sint^2*yc^2 - a^2*b^2 = 0
541
542 // Divide by b^2
543
544 // (a^2/b^2)*sint^2*x^2 + cost^2*x^2 +
545 // 2*sint*cost*x*y - (a^2/b^2)*2*sint*cost*x*y +
546 // (a^2/b^2)*cost^2*y^2 + sint^2*y^2 +
547 // -(a^2/b^2)*2*sint^2*xc*x - 2*cost^2*xc*x - 2*sint*cost*yc*x + (a^2/b^2)*2*sint*cost*yc*x +
548 // -2*sint*cost*xc*y + (a^2/b^2)*2*sint*cost*xc*y - (a^2/b^2)*2*cost^2*yc*y - 2*sint^2*yc*y +
549 // (a^2/b^2)*sint^2*xc^2 + cost^2*xc^2 + 2*sint*cost*xc*yc - (a^2/b^2)*2*sint*cost*xc*yc + (a^2/b^2)*cost^2*yc^2 + sint^2*yc^2 - a^2 = 0
550
551 // (a^2/b^2)*(sint^2*x^2 - 2*sint*cost*x*y + cost^2*y^2 - 2*sint^2*xc*x + 2*sint*cost*yc*x + 2*sint*cost*xc*y - 2*cost^2*yc*y + sint^2*xc^2 - 2*sint*cost*xc*yc + cost^2*yc^2)
552 // - a^2
553 // = - cost^2*x^2 - 2*sint*cost*x*y - sint^2*y^2 + 2*cost^2*xc*x + 2*sint*cost*yc*x + 2*sint*cost*xc*y + 2*sint^2*yc*y - cost^2*xc^2 - 2*sint*cost*xc*yc - sint^2*yc^2
554
555 // Unknowns are:
556 // a^2/b^2 and a^2
557
558 try {
559 final var sint = Math.sin(theta);
560 final var cost = Math.cos(theta);
561
562 final var sint2 = sint * sint;
563 final var cost2 = cost * cost;
564 final var sintcost = sint * cost;
565
566 final var xc = center.getInhomX();
567 final var yc = center.getInhomY();
568 final var xc2 = xc * xc;
569 final var yc2 = yc * yc;
570
571 final var m = new Matrix(2, 2);
572 final var b = new double[2];
573
574 final var x = new double[]{
575 point1.getInhomX(),
576 point2.getInhomX(),
577 };
578
579 final var y = new double[]{
580 point1.getInhomY(),
581 point2.getInhomY(),
582 };
583
584 double x2;
585 double y2;
586 double rowNorm;
587 for (var i = 0; i < 2; i++) {
588 x2 = x[i] * x[i];
589 y2 = y[i] * y[i];
590
591 final var tmp = 2.0 * sintcost * x[i] * y[i];
592
593 m.setElementAt(i, 0, sint2 * x2 - tmp + cost2 * y2 - 2.0 * sint2 * xc * x[i]
594 + 2.0 * sintcost * yc * x[i] + 2.0 * sintcost * xc * y[i] - 2.0 * cost2 * yc * y[i]
595 + sint2 * xc2 - 2.0 * sintcost * xc * yc + cost2 * yc2);
596 m.setElementAt(i, 1, -1.0);
597
598 b[i] = -cost2 * x2 - tmp - sint2 * y2 + 2.0 * cost2 * xc * x[i] + 2.0 * sintcost * yc * x[i]
599 + 2.0 * sintcost * xc * y[i] + 2.0 * sint2 * yc * y[i] - cost2 * xc2 - 2.0 * sintcost * xc * yc
600 - sint2 * yc2;
601
602 // normalize row to increase accuracy
603 rowNorm = 0.0;
604 for (var j = 0; j < 2; j++) {
605 rowNorm += Math.pow(m.getElementAt(i, j), 2.0);
606 }
607 rowNorm = Math.sqrt(rowNorm);
608
609 for (var j = 0; j < 2; j++) {
610 m.setElementAt(i, j, m.getElementAt(i, j) / rowNorm);
611 }
612
613 b[i] /= rowNorm;
614 }
615
616 final var params = com.irurueta.algebra.Utils.solve(m, b);
617
618 // params[1] = a^2
619 final var sMajorAxis = Math.sqrt(Math.abs(params[1]));
620
621 // params[0] = a^2/b^2 --> b^2 = a^2/params[0] = params[1]/params[0]
622 final var sMinorAxis = Math.sqrt(Math.abs(params[1] / params[0]));
623
624 setCenterAxesAndRotation(center, sMajorAxis, sMinorAxis, theta);
625 } catch (final AlgebraException e) {
626 throw new ColinearPointsException(e);
627 }
628 }
629
630 /**
631 * Sets parameters from 5 points.
632 *
633 * @param point1 1st point.
634 * @param point2 2nd point.
635 * @param point3 3rd point.
636 * @param point4 4th point.
637 * @param point5 5th point.
638 * @throws ColinearPointsException if points are in a co-linear or degenerate
639 * configuration.
640 */
641 public final void setParametersFromPoints(
642 final Point2D point1, final Point2D point2, final Point2D point3, final Point2D point4,
643 final Point2D point5) throws ColinearPointsException {
644 setParametersFromPoints(point1, point2, point3, point4, point5, 0.0);
645 }
646
647 /**
648 * Sets parameters from 5 points.
649 *
650 * @param point1 1st point.
651 * @param point2 2nd point.
652 * @param point3 3rd point.
653 * @param point4 4th point.
654 * @param point5 5th point.
655 * @param threshold threshold to determine whether points form an ellipse.
656 * This is usually a very small value
657 * @throws ColinearPointsException if points are in a co-linear or degenerate
658 * configuration.
659 */
660 public final void setParametersFromPoints(
661 final Point2D point1, final Point2D point2, final Point2D point3, final Point2D point4,
662 final Point2D point5, final double threshold) throws ColinearPointsException {
663 // normalize points to increase accuracy
664 point1.normalize();
665 point2.normalize();
666 point3.normalize();
667 point4.normalize();
668 point5.normalize();
669
670 try {
671 // each point belonging to a conic follows equation:
672 // p' * C * p = 0 ==>
673 // x^2 + y^2 + w^2 + 2*x*y + 2*x*w + 2*y*w = 0
674 final var m = new Matrix(5, 6);
675 var x = point1.getHomX();
676 var y = point1.getHomY();
677 var w = point1.getHomW();
678 m.setElementAt(0, 0, x * x);
679 m.setElementAt(0, 1, 2.0 * x * y);
680 m.setElementAt(0, 2, y * y);
681 m.setElementAt(0, 3, 2.0 * x * w);
682 m.setElementAt(0, 4, 2.0 * y * w);
683 m.setElementAt(0, 5, w * w);
684 x = point2.getHomX();
685 y = point2.getHomY();
686 w = point2.getHomW();
687 m.setElementAt(1, 0, x * x);
688 m.setElementAt(1, 1, 2.0 * x * y);
689 m.setElementAt(1, 2, y * y);
690 m.setElementAt(1, 3, 2.0 * x * w);
691 m.setElementAt(1, 4, 2.0 * y * w);
692 m.setElementAt(1, 5, w * w);
693 x = point3.getHomX();
694 y = point3.getHomY();
695 w = point3.getHomW();
696 m.setElementAt(2, 0, x * x);
697 m.setElementAt(2, 1, 2.0 * x * y);
698 m.setElementAt(2, 2, y * y);
699 m.setElementAt(2, 3, 2.0 * x * w);
700 m.setElementAt(2, 4, 2.0 * y * w);
701 m.setElementAt(2, 5, w * w);
702 x = point4.getHomX();
703 y = point4.getHomY();
704 w = point4.getHomW();
705 m.setElementAt(3, 0, x * x);
706 m.setElementAt(3, 1, 2.0 * x * y);
707 m.setElementAt(3, 2, y * y);
708 m.setElementAt(3, 3, 2.0 * x * w);
709 m.setElementAt(3, 4, 2.0 * y * w);
710 m.setElementAt(3, 5, w * w);
711 x = point5.getHomX();
712 y = point5.getHomY();
713 w = point5.getHomW();
714 m.setElementAt(4, 0, x * x);
715 m.setElementAt(4, 1, 2.0 * x * y);
716 m.setElementAt(4, 2, y * y);
717 m.setElementAt(4, 3, 2.0 * x * w);
718 m.setElementAt(4, 4, 2.0 * y * w);
719 m.setElementAt(4, 5, w * w);
720
721 // normalize each row to increase accuracy
722 final var row = new double[6];
723 double rowNorm;
724
725 for (var j = 0; j < 5; j++) {
726 m.getSubmatrixAsArray(j, 0, j, 5, row);
727 rowNorm = com.irurueta.algebra.Utils.normF(row);
728 for (var i = 0; i < 6; i++) {
729 m.setElementAt(j, i, m.getElementAt(j, i) / rowNorm);
730 }
731 }
732
733 final var decomposer = new SingularValueDecomposer(m);
734 decomposer.decompose();
735
736 if (decomposer.getRank() < 5) {
737 throw new ColinearPointsException();
738 }
739
740 // the right null-space of m contains the parameters a, b, c, d, e ,f
741 // of the conic
742 final var v = decomposer.getV();
743
744 final var aPrime = v.getElementAt(0, 5);
745 final var bPrime = v.getElementAt(1, 5);
746 final var cPrime = v.getElementAt(2, 5);
747 final var dPrime = v.getElementAt(3, 5);
748 final var ePrime = v.getElementAt(4, 5);
749 final var fPrime = v.getElementAt(5, 5);
750
751 // an ellipse follows the generic conic equation
752 // A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
753
754 // Or in matrix form
755 // [A B/2 D/2 ] [A' B' D']
756 // [B/2 C E/2 ] = [B' C' E']
757 // [D/2 E/2 F ] [D' E' F']
758
759 final var b = 2.0 * bPrime;
760 final var d = 2.0 * dPrime;
761 final var e = 2.0 * ePrime;
762
763 setParameters(aPrime, b, cPrime, d, e, fPrime, threshold);
764 } catch (final AlgebraException | IllegalArgumentException ex) {
765 throw new ColinearPointsException(ex);
766 }
767 }
768
769 /**
770 * Sets parameters of canonical equation of an ellipse, which is:
771 * a*x^2 + b*x*y + c*y^2 + d*x + e*y + f = 0
772 *
773 * @param a a parameter.
774 * @param b b parameter.
775 * @param c c parameter.
776 * @param d d parameter.
777 * @param e e parameter.
778 * @param f f parameter.
779 * @throws IllegalArgumentException if parameters do not follow
780 * b^2 - 4*a*c < 0.0.
781 */
782 public final void setParameters(
783 final double a, final double b, final double c, final double d, final double e, final double f) {
784 setParameters(a, b, c, d, e, f, 0.0);
785 }
786
787 /**
788 * Sets parameters of canonical equation of an ellipse, which is:
789 * a*x^2 + b*x*y + c*y^2 + d*x + e*y + f = 0
790 *
791 * @param a a parameter.
792 * @param b b parameter.
793 * @param c c parameter.
794 * @param d d parameter.
795 * @param e e parameter.
796 * @param f f parameter.
797 * @param threshold threshold to determine whether parameters are valid due
798 * to machine precision.
799 * @throws IllegalArgumentException if parameters do not follow
800 * b^2 - 4*A*c < threshold.
801 */
802 public final void setParameters(
803 final double a, final double b, final double c, final double d, final double e, final double f,
804 final double threshold) {
805 final var discriminant = b * b - 4.0 * a * c;
806 if (discriminant >= threshold) {
807 // not an ellipse
808 throw new IllegalArgumentException();
809 }
810
811 final var tmp1 = a * e * e + c * d * d - b * d * e + discriminant * f;
812 final var tmp2 = Math.sqrt((a - c) * (a - c) + b * b);
813 final var sMajorAxis = -Math.sqrt(2.0 * tmp1 * (a + c + tmp2)) / discriminant;
814
815 final var sMinorAxis = -Math.sqrt(2.0 * tmp1 * (a + c - tmp2)) / discriminant;
816
817 final var xc = (2.0 * c * d - b * e) / discriminant;
818 final var yc = (2.0 * a * e - b * d) / discriminant;
819 final var centerPoint = new InhomogeneousPoint2D(xc, yc);
820
821 double alpha;
822 if (Math.abs(b) <= threshold && a < c) {
823 alpha = 0.0;
824 } else {
825 alpha = Math.atan2(c - a - tmp2, b);
826 }
827
828 setCenterAxesAndRotation(centerPoint, sMajorAxis, sMinorAxis, alpha);
829 }
830
831 /**
832 * Gets focus distance of ellipse.
833 * Focus determines the distance respect to the center of the ellipse
834 * where the two focus points are located.
835 * The sum of the distances from any point P = P(x,y) on the ellipse to
836 * those two foci is constant and equal to the major axis length.
837 *
838 * @return focus distance.
839 */
840 public double getFocus() {
841 return Math.sqrt(Math.pow(semiMajorAxis, 2.0) - Math.pow(semiMinorAxis, 2.0));
842 }
843
844 /**
845 * Gets semi major axis x,y coordinates.
846 *
847 * @param coords array where x, y coordinates of semi major axis will be
848 * stored.
849 */
850 public void getSemiMajorAxisCoordinates(final double[] coords) {
851 coords[0] = semiMajorAxis * Math.cos(rotationAngle);
852 coords[1] = semiMajorAxis * Math.sin(rotationAngle);
853 }
854
855 /**
856 * Gets semi major axis x,y coordinates.
857 *
858 * @return array containing x, y coordinates of semi major axis.
859 */
860 public double[] getSemiMajorAxisCoordinates() {
861 final var coords = new double[Point2D.POINT2D_INHOMOGENEOUS_COORDINATES_LENGTH];
862 getSemiMajorAxisCoordinates(coords);
863 return coords;
864 }
865
866 /**
867 * Sets semi major axis coordinates.
868 * This method updates rotation angle and semi major axis length.
869 *
870 * @param coords coordinates of semi major axis.
871 */
872 public void setSemiMajorAxisCoordinates(final double[] coords) {
873 rotationAngle = Math.atan2(coords[1], coords[0]);
874 semiMajorAxis = com.irurueta.algebra.Utils.normF(coords);
875 }
876
877 /**
878 * Gets semi minor axis x,y coordinates.
879 *
880 * @param coords array where x, y coordinates of semi minor axis will be
881 * stored.
882 */
883 public void getSemiMinorAxisCoordinates(final double[] coords) {
884 coords[0] = -semiMinorAxis * Math.sin(rotationAngle);
885 coords[1] = semiMinorAxis * Math.cos(rotationAngle);
886 }
887
888 /**
889 * Gets semi minor axis x,y coordinates.
890 *
891 * @return array containing x,y coordinates of semi minor axis.
892 */
893 public double[] getSemiMinorAxisCoordinates() {
894 final var coords = new double[Point2D.POINT2D_INHOMOGENEOUS_COORDINATES_LENGTH];
895 getSemiMinorAxisCoordinates(coords);
896 return coords;
897 }
898
899 /**
900 * Sets semi minor axis coordinates.
901 * This method updates rotation angle and semi minor axis length.
902 *
903 * @param coords coordinates of semi minor axis.
904 */
905 public void setSemiMinorAxisCoordinates(final double[] coords) {
906 rotationAngle = Math.atan2(-coords[0], coords[1]);
907 semiMinorAxis = com.irurueta.algebra.Utils.normF(coords);
908 }
909
910 /**
911 * Gets 1st focus point.
912 * The sum of the distances from any point P = P(x,y) on the ellipse to
913 * this focus is constant and equal to the major axis length.
914 *
915 * @param focusPoint1 1st focus point.
916 */
917 public void getFocusPoint1(final Point2D focusPoint1) {
918 final var focus = getFocus();
919
920 focusPoint1.setInhomogeneousCoordinates(center.getInhomX() - focus * Math.cos(rotationAngle),
921 center.getInhomY() - focus * Math.sin(rotationAngle));
922 }
923
924 /**
925 * Gets 1st focus point.
926 * The sum of the distances from any point P = P(x,y) on the ellipse to
927 * this focus is constant and equal to the major axis length.
928 *
929 * @return 1st focus point.
930 */
931 public Point2D getFocusPoint1() {
932 final var result = Point2D.create();
933 getFocusPoint1(result);
934 return result;
935 }
936
937 /**
938 * Gets 2nd focus point.
939 * The sum of the distances from any point P = P(x,y) on the ellipse to
940 * this focus is constant and equal to the major axis length.
941 *
942 * @param focusPoint2 2nd focus point.
943 */
944 public void getFocusPoint2(final Point2D focusPoint2) {
945 final var focus = getFocus();
946
947 focusPoint2.setInhomogeneousCoordinates(center.getInhomX() + focus * Math.cos(rotationAngle),
948 center.getInhomY() + focus * Math.sin(rotationAngle));
949 }
950
951 /**
952 * Gets 2nd focus point.
953 * The sum of the distances from any point P = P(x,y) on the ellipse to
954 * this focus is constant and equal to the major axis length.
955 *
956 * @return 2nd focus point.
957 */
958 public Point2D getFocusPoint2() {
959 final var result = Point2D.create();
960 getFocusPoint2(result);
961 return result;
962 }
963
964 /**
965 * Sets focus points.
966 *
967 * @param focusPoint1 1st focus point.
968 * @param focusPoint2 2nd focus point.
969 * @param keepSemiMinorAxis true indicates that semi-minor axis is kept,
970 * false indicates that semi-major axis is kept instead.
971 */
972 public void setFocusPoints(final Point2D focusPoint1, final Point2D focusPoint2, final boolean keepSemiMinorAxis) {
973 final var c = new InhomogeneousPoint2D((focusPoint1.getInhomX() + focusPoint2.getInhomX()) / 2.0,
974 (focusPoint1.getInhomY() + focusPoint2.getInhomY()) / 2.0);
975 final var alpha = Math.atan2(focusPoint2.getInhomY() - focusPoint1.getInhomY(),
976 focusPoint2.getInhomX() - focusPoint1.getInhomX());
977
978 final var f = focusPoint1.distanceTo(c);
979 final var f2 = f * f;
980 final double a2;
981 final double b2;
982 if (keepSemiMinorAxis) {
983 b2 = semiMinorAxis * semiMinorAxis;
984 a2 = f2 + b2;
985 } else {
986 a2 = semiMajorAxis * semiMajorAxis;
987 b2 = a2 - f2;
988 }
989
990 final var sMajorAxis = Math.sqrt(a2);
991 final var sMinorAxis = Math.sqrt(b2);
992
993 setCenterAxesAndRotation(c, sMajorAxis, sMinorAxis, alpha);
994 }
995
996 /**
997 * Gets eccentricity of ellipsis.
998 *
999 * @return eccentricity of ellipsis.
1000 */
1001 public double getEccentricity() {
1002 return getFocus() / semiMajorAxis;
1003 }
1004
1005 /**
1006 * Returns area of this ellipse.
1007 *
1008 * @return area of this ellipse.
1009 */
1010 public double getArea() {
1011 return Math.PI * semiMajorAxis * semiMinorAxis;
1012 }
1013
1014 /**
1015 * Returns perimeter of this ellipse.
1016 *
1017 * @return Perimeter of this ellipse.
1018 */
1019 public double getPerimeter() {
1020 final var a = semiMajorAxis;
1021 final var b = semiMinorAxis;
1022
1023 return Math.PI * (3.0 * (a + b) - Math.sqrt((3.0 * a + b) * (a + 3.0 * b)));
1024 }
1025
1026 /**
1027 * Gets curvature of ellipse at provided point.
1028 *
1029 * @param point point to be checked.
1030 * @return curvature of ellipse.
1031 */
1032 public double getCurvature(final Point2D point) {
1033 final var a = semiMajorAxis;
1034 final var b = semiMinorAxis;
1035
1036 final var a2 = a * a;
1037 final var b2 = b * b;
1038
1039 final var a4 = a2 * a2;
1040 final var b4 = b2 * b2;
1041
1042 final var x = point.getInhomX() - center.getInhomX();
1043 final var y = point.getInhomY() - center.getInhomY();
1044
1045 final var x2 = x * x;
1046 final var y2 = y * y;
1047
1048 return Math.pow(x2 / a4 + y2 / b4, -3.0 / 2.0) / (a2 * b2);
1049 }
1050
1051 /**
1052 * Converts this circle into a conic.
1053 * Conics are a more general representation of circles.
1054 *
1055 * @return A conic representing this circle
1056 */
1057 public Conic toConic() {
1058 center.normalize();
1059 // use inhomogeneous center coordinates
1060 final var xc = center.getInhomX();
1061 final var yc = center.getInhomY();
1062
1063 final var sint = Math.sin(rotationAngle);
1064 final var cost = Math.cos(rotationAngle);
1065
1066 final var a = semiMajorAxis;
1067 final var b = semiMinorAxis;
1068
1069 final var xc2 = xc * xc;
1070 final var yc2 = yc * yc;
1071
1072 final var sint2 = sint * sint;
1073 final var cost2 = cost * cost;
1074
1075 final var a2 = a * a;
1076 final var b2 = b * b;
1077
1078 final var aParam = a2 * sint2 + b2 * cost2;
1079 final var bParam = 2.0 * (b2 - a2) * sint * cost;
1080 final var cParam = a2 * cost2 + b2 * sint2;
1081 final var dParam = -2.0 * aParam * xc - bParam * yc;
1082 final var eParam = -bParam * xc - 2.0 * cParam * yc;
1083 final var fParam = aParam * xc2 + bParam * xc * yc + cParam * yc2 - a2 * b2;
1084
1085 final var bConic = bParam / 2.0;
1086 final var dConic = dParam / 2.0;
1087 final var eConic = eParam / 2.0;
1088
1089 return new Conic(aParam, bConic, cParam, dConic, eConic, fParam);
1090 }
1091
1092 /**
1093 * Set parameters of this circle from a valid conic corresponding to a
1094 * circle.
1095 *
1096 * @param conic conic to set parameters from.
1097 * @throws IllegalArgumentException if provided conic is not an ellipse.
1098 */
1099 public final void setFromConic(final Conic conic) {
1100 if (conic.getConicType() != ConicType.ELLIPSE_CONIC_TYPE
1101 && conic.getConicType() != ConicType.CIRCLE_CONIC_TYPE) {
1102 throw new IllegalArgumentException();
1103 }
1104
1105 conic.normalize();
1106
1107 final var aConic = conic.getA();
1108 final var bConic = conic.getB();
1109 final var cConic = conic.getC();
1110 final var dConic = conic.getD();
1111 final var eConic = conic.getE();
1112 final var fConic = conic.getF();
1113
1114 // an ellipse follows the generic conic equation
1115 // A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
1116
1117 // Or in matrix form
1118 // [A B/2 D/2 ] [A' B' D']
1119 // [B/2 C E/2 ] = [B' C' E']
1120 // [D/2 E/2 F ] [D' E' F']
1121
1122 final var b = 2.0 * bConic;
1123 final var d = 2.0 * dConic;
1124 final var e = 2.0 * eConic;
1125
1126 setParameters(aConic, b, cConic, d, e, fConic);
1127 }
1128
1129 /**
1130 * Sets parameters of this ellipse from a circle.
1131 *
1132 * @param circle a circle to set parameters from.
1133 */
1134 public final void setFromCircle(final Circle circle) {
1135 center = circle.getCenter();
1136 semiMajorAxis = semiMinorAxis = circle.getRadius();
1137 rotationAngle = 0.0;
1138 }
1139
1140 /**
1141 * Determines if provided point is inside this ellipse or not up to a
1142 * certain threshold.
1143 * If provided threshold is positive, the ellipse behaves as if it was a
1144 * larger ellipse increased by threshold amount, if provided threshold is
1145 * negative, the ellipse behaves as if it was a smaller ellipse decreased by
1146 * threshold amount in radius.
1147 *
1148 * @param point Point to be checked.
1149 * @param threshold Threshold to determine if point is inside or not.
1150 * @return True if point is considered to be inside this circle, false
1151 * otherwise.
1152 */
1153 public boolean isInside(final Point2D point, final double threshold) {
1154 center.normalize();
1155 // use inhomogeneous center coordinates
1156 final var xc = center.getInhomX();
1157 final var yc = center.getInhomY();
1158
1159 final var sint = Math.sin(rotationAngle);
1160 final var cost = Math.cos(rotationAngle);
1161
1162 final var a = semiMajorAxis;
1163 final var b = semiMinorAxis;
1164
1165 final var xc2 = xc * xc;
1166 final var yc2 = yc * yc;
1167
1168 final var sint2 = sint * sint;
1169 final var cost2 = cost * cost;
1170
1171 final var a2 = a * a;
1172 final var b2 = b * b;
1173
1174 final var aParam = a2 * sint2 + b2 * cost2;
1175 final var bParam = 2.0 * (b2 - a2) * sint * cost;
1176 final var cParam = a2 * cost2 + b2 * sint2;
1177 final var dParam = -2.0 * aParam * xc - bParam * yc;
1178 final var eParam = -bParam * xc - 2.0 * cParam * yc;
1179 final var fParam = aParam * xc2 + bParam * xc * yc + cParam * yc2 - a2 * b2;
1180
1181 final var x = point.getInhomX();
1182 final var y = point.getInhomY();
1183
1184 return aParam * x * x + bParam * x * y + cParam * y * y + dParam * x + eParam * y + fParam <= threshold;
1185 }
1186
1187 /**
1188 * Determines if provided point is inside this circle or not.
1189 *
1190 * @param point Point to be checked.
1191 * @return True if point is considered to be inside this circle, false
1192 * otherwise.
1193 */
1194 public boolean isInside(final Point2D point) {
1195 return isInside(point, 0.0);
1196 }
1197
1198 /**
1199 * Determines whether provided point lies at ellipse boundary or not up to
1200 * a certain threshold.
1201 *
1202 * @param point Point to be checked.
1203 * @param threshold A small threshold to determine whether point lies at
1204 * ellipse boundary.
1205 * @return True if point lies at ellipse boundary, false otherwise.
1206 * @throws IllegalArgumentException Raised if provided threshold is
1207 * negative.
1208 */
1209 public boolean isLocus(final Point2D point, final double threshold) {
1210 if (threshold < MIN_THRESHOLD) {
1211 throw new IllegalArgumentException();
1212 }
1213
1214 center.normalize();
1215 // use inhomogeneous center coordinates
1216 final var xc = center.getInhomX();
1217 final var yc = center.getInhomY();
1218
1219 final var sint = Math.sin(rotationAngle);
1220 final var cost = Math.cos(rotationAngle);
1221
1222 final var a = semiMajorAxis;
1223 final var b = semiMinorAxis;
1224
1225 final var xc2 = xc * xc;
1226 final var yc2 = yc * yc;
1227
1228 final var sint2 = sint * sint;
1229 final var cost2 = cost * cost;
1230
1231 final var a2 = a * a;
1232 final var b2 = b * b;
1233
1234 final var aParam = a2 * sint2 + b2 * cost2;
1235 final var bParam = 2.0 * (b2 - a2) * sint * cost;
1236 final var cParam = a2 * cost2 + b2 * sint2;
1237 final var dParam = -2.0 * aParam * xc - bParam * yc;
1238 final var eParam = -bParam * xc - 2.0 * cParam * yc;
1239 final var fParam = aParam * xc2 + bParam * xc * yc + cParam * yc2 - a2 * b2;
1240
1241 final var x = point.getInhomX();
1242 final var y = point.getInhomY();
1243
1244 return Math.abs(aParam * x * x + bParam * x * y + cParam * y * y + dParam * x + eParam * y + fParam) <= threshold;
1245 }
1246
1247 /**
1248 * Determines whether provided point lies at ellipse boundary or not.
1249 *
1250 * @param point Point to be checked.
1251 * @return True if point lies at ellipse boundary, false otherwise.
1252 */
1253 public boolean isLocus(final Point2D point) {
1254 return isLocus(point, DEFAULT_THRESHOLD);
1255 }
1256
1257 /**
1258 * Returns a line tangent to this ellipse at provided point. Provided point
1259 * must be locus of this ellipse, otherwise a NotLocusException will be
1260 * thrown.
1261 *
1262 * @param point a locus point of this ellipse.
1263 * @return a 2D line tangent to this ellipse at provided point.
1264 * @throws NotLocusException if provided point is not locus of this ellipse
1265 * up to DEFAULT_THRESHOLD.
1266 */
1267 public Line2D getTangentLineAt(final Point2D point) throws NotLocusException {
1268 return getTangentLineAt(point, DEFAULT_THRESHOLD);
1269 }
1270
1271 /**
1272 * Returns a line tangent to this ellipse at provided point. Provided point
1273 * must be locus of this ellipse, otherwise a NotLocusException will be
1274 * thrown.
1275 *
1276 * @param point a locus point of this circle.
1277 * @param threshold threshold to determine if provided point is locus.
1278 * @return a 2D line tangent to this circle at provided point.
1279 * @throws NotLocusException if provided point is not locus of this circle
1280 * up to provided threshold.
1281 * @throws IllegalArgumentException if provided threshold is negative.
1282 */
1283 public Line2D getTangentLineAt(final Point2D point, final double threshold) throws NotLocusException {
1284 final var line = new Line2D();
1285 tangentLineAt(point, line, threshold);
1286 return line;
1287 }
1288
1289 /**
1290 * Computes a line tangent to this circle at provided point. Provided point
1291 * must be locus of this circle, otherwise a NotLocusException will be
1292 * thrown.
1293 *
1294 * @param point a locus point of this circle.
1295 * @param line instance of a 2D line where result will be stored.
1296 * @param threshold threshold to determine if provided point is locus.
1297 * @throws NotLocusException if provided point is not locus of this circle
1298 * up to provided threshold.
1299 * @throws IllegalArgumentException if provided threshold is negative.
1300 */
1301 public void tangentLineAt(final Point2D point, final Line2D line, final double threshold) throws NotLocusException {
1302 if (!isLocus(point, threshold)) {
1303 throw new NotLocusException();
1304 }
1305
1306 final var c = toConic();
1307 c.tangentLineAt(point, line, threshold);
1308 }
1309 }