1 /*
2 * Copyright (C) 2012 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 import com.irurueta.algebra.WrongSizeException;
22
23 import java.io.Serializable;
24
25 /**
26 * This class contains the implementation of a quadric.
27 */
28 @SuppressWarnings("DuplicatedCode")
29 public class Quadric extends BaseQuadric implements Serializable {
30
31 /**
32 * Constructor.
33 */
34 public Quadric() {
35 super();
36 }
37
38 /**
39 * Constructor of this class. This constructor accepts every parameter
40 * describing a quadric (parameters a, b, c, d, e, f, g, h, i, j).
41 *
42 * @param a Parameter A of the quadric.
43 * @param b Parameter B of the quadric.
44 * @param c Parameter C of the quadric.
45 * @param d Parameter D of the quadric.
46 * @param e Parameter E of the quadric.
47 * @param f Parameter F of the quadric.
48 * @param g Parameter G of the quadric.
49 * @param h Parameter H of the quadric.
50 * @param i Parameter I of the quadric.
51 * @param j Parameter J of the quadric.
52 */
53 public Quadric(final double a, final double b, final double c, final double d, final double e, final double f,
54 final double g, final double h, final double i, final double j) {
55 super(a, b, c, d, e, f, g, h, i, j);
56 }
57
58 /**
59 * Constructor of this class. This constructor accepts a Matrix describing
60 * a quadric.
61 *
62 * @param m Matrix describing a quadric 4x4 Matrix describing the quadric.
63 * @throws IllegalArgumentException Raised when the size of the matrix is
64 * not 4x4.
65 * @throws NonSymmetricMatrixException Raised when the quadric matrix is not
66 * symmetric.
67 */
68 public Quadric(final Matrix m) throws NonSymmetricMatrixException {
69 super(m);
70 }
71
72 /**
73 * Creates quadric where provided points are contained (are locus).
74 *
75 * @param point1 1st point.
76 * @param point2 2nd point.
77 * @param point3 3rd point.
78 * @param point4 4th point.
79 * @param point5 5th point.
80 * @param point6 6th point.
81 * @param point7 7th point.
82 * @param point8 8th point.
83 * @param point9 9th point.
84 * @throws CoincidentPointsException Raised if points are coincident or
85 * produce a degenerated configuration.
86 */
87 public Quadric(final Point3D point1, final Point3D point2, final Point3D point3, final Point3D point4,
88 final Point3D point5, final Point3D point6, final Point3D point7, final Point3D point8,
89 final Point3D point9) throws CoincidentPointsException {
90 setParametersFromPoints(point1, point2, point3, point4, point5, point6, point7, point8, point9);
91 }
92
93 /**
94 * Checks if the given point is locus (lies within) this quadric.
95 *
96 * @param point Point to be checked.
97 * @param threshold Threshold of distance to determine whether the point
98 * is locus of the quadric or not. Threshold might be needed because of
99 * machine precision. If not provided DEFAULT_LOCUS_THRESHOLD will be used
100 * instead.
101 * @return True if the point lies within this quadric, false otherwise.
102 * @throws IllegalArgumentException Raised if threshold is negative.
103 */
104 public boolean isLocus(final Point3D point, final double threshold) {
105 if (threshold < MIN_THRESHOLD) {
106 throw new IllegalArgumentException();
107 }
108
109 try {
110 normalize();
111 final var q = asMatrix();
112 final var homPoint = new Matrix(Point3D.POINT3D_HOMOGENEOUS_COORDINATES_LENGTH, 1);
113 point.normalize();
114 homPoint.setElementAt(0, 0, point.getHomX());
115 homPoint.setElementAt(1, 0, point.getHomY());
116 homPoint.setElementAt(2, 0, point.getHomZ());
117 homPoint.setElementAt(3, 0, point.getHomW());
118 final var locusMatrix = homPoint.transposeAndReturnNew();
119 locusMatrix.multiply(q);
120 locusMatrix.multiply(homPoint);
121
122 return Math.abs(locusMatrix.getElementAt(0, 0)) < threshold;
123 } catch (final WrongSizeException e) {
124 return false;
125 }
126 }
127
128 /**
129 * Checks if the given point is locus (lies within) this quadric.
130 *
131 * @param point Point to be checked.
132 * @return True if the point lies within this conic, false otherwise
133 * @see #isLocus(Point3D, double)
134 */
135 public boolean isLocus(final Point3D point) {
136 return isLocus(point, DEFAULT_LOCUS_THRESHOLD);
137 }
138
139 /**
140 * Computes the angle between two 3D points using this quadric as a geometry
141 * base.
142 *
143 * @param pointA First point.
144 * @param pointB Second point.
145 * @return Angle between provided points given in radians.
146 */
147 public double angleBetweenPoints(final Point3D pointA, final Point3D pointB) {
148 try {
149 // retrieve quadric as matrix
150 final var q = asMatrix();
151 final var transHomPointA = new Matrix(1, Point3D.POINT3D_HOMOGENEOUS_COORDINATES_LENGTH);
152 pointA.normalize();
153 transHomPointA.setElementAt(0, 0, pointA.getHomX());
154 transHomPointA.setElementAt(0, 1, pointA.getHomY());
155 transHomPointA.setElementAt(0, 2, pointA.getHomZ());
156 transHomPointA.setElementAt(0, 3, pointA.getHomW());
157
158 final var tmp = transHomPointA.multiplyAndReturnNew(q);
159 tmp.multiply(transHomPointA.transposeAndReturnNew()); //This is
160 // homPointA' * Q * homPointA
161
162 final var normA = tmp.getElementAt(0, 0);
163
164 final var homPointB = new Matrix(Point3D.POINT3D_HOMOGENEOUS_COORDINATES_LENGTH, 1);
165 pointB.normalize();
166 homPointB.setElementAt(0, 0, pointB.getHomX());
167 homPointB.setElementAt(1, 0, pointB.getHomY());
168 homPointB.setElementAt(2, 0, pointB.getHomZ());
169 homPointB.setElementAt(3, 0, pointB.getHomW());
170
171 homPointB.transpose(tmp);
172 tmp.multiply(q);
173 tmp.multiply(homPointB);
174
175 final var normB = tmp.getElementAt(0, 0);
176
177 transHomPointA.multiply(q);
178 transHomPointA.multiply(homPointB);
179 // This is homPointA' * Q * homPointB
180
181 final var angleNumerator = transHomPointA.getElementAt(0, 0);
182
183 final var cosTheta = angleNumerator / Math.sqrt(normA * normB);
184 return Math.acos(cosTheta);
185 } catch (final WrongSizeException ignore) {
186 // This will never happen
187 return 0.0;
188 }
189 }
190
191 /**
192 * Checks if two points are perpendicular in the geometry base generated by
193 * this quadric.
194 *
195 * @param pointA First point.
196 * @param pointB Second point.
197 * @param threshold Threshold to determine whether the points are
198 * perpendicular or not. If the dot product between provided points and this
199 * quadric is greater than provided threshold, then points won't be assumed
200 * to be perpendicular. Threshold is provided because of machine precision
201 * limits, if not provided DEFAULT_PERPENDICULAR_THRESHOLD will be used
202 * instead.
203 * @return True if points are perpendicular, false otherwise.
204 * @throws IllegalArgumentException Raised if threshold is negative.
205 */
206 public boolean arePerpendicularPoints(final Point3D pointA, final Point3D pointB, final double threshold) {
207 try {
208 // retrieve quadric as matrix
209 final var transHomPointA = new Matrix(1, Point3D.POINT3D_HOMOGENEOUS_COORDINATES_LENGTH);
210 pointA.normalize();
211 transHomPointA.setElementAt(0, 0, pointA.getHomX());
212 transHomPointA.setElementAt(0, 1, pointA.getHomY());
213 transHomPointA.setElementAt(0, 2, pointA.getHomZ());
214 transHomPointA.setElementAt(0, 3, pointA.getHomW());
215
216 final var homPointB = new Matrix(Point3D.POINT3D_HOMOGENEOUS_COORDINATES_LENGTH, 1);
217 pointB.normalize();
218 homPointB.setElementAt(0, 0, pointB.getHomX());
219 homPointB.setElementAt(1, 0, pointB.getHomY());
220 homPointB.setElementAt(2, 0, pointB.getHomZ());
221 homPointB.setElementAt(3, 0, pointB.getHomW());
222
223 normalize();
224 final var q = asMatrix();
225 transHomPointA.multiply(q);
226 transHomPointA.multiply(homPointB);
227 // This is homPointA' * Q * homPointB
228
229 final var perpend = transHomPointA.getElementAt(0, 0);
230
231 return Math.abs(perpend) < threshold;
232 } catch (final WrongSizeException ignore) {
233 // This will never happen
234 return false;
235 }
236 }
237
238 /**
239 * Checks if two points are perpendicular in the geometry base generated by
240 * this quadric.
241 *
242 * @param pointA First point.
243 * @param pointB Second point.
244 * @return True if points are perpendicular, false otherwise.
245 * @see #arePerpendicularPoints(Point3D, Point3D, double)
246 */
247 public boolean arePerpendicularPoints(final Point3D pointA, final Point3D pointB) {
248 return arePerpendicularPoints(pointA, pointB, DEFAULT_PERPENDICULAR_THRESHOLD);
249 }
250
251 /**
252 * Sets the values of the dual quadric corresponding to this quadric
253 * instance into provided dualQuadric instance.
254 * The dual quadric is equal to the inverse of the quadric matrix.
255 *
256 * @param dualQuadric Dual quadric instance where the values of the dual
257 * quadric of this quadric instance will be stored.
258 * @throws DualQuadricNotAvailableException Raised if the dual quadric does
259 * not exist because this quadric instance is degenerate (its inverse
260 * cannot be computed).
261 */
262 public void dualQuadric(final DualQuadric dualQuadric) throws DualQuadricNotAvailableException {
263 final var quadricMatrix = asMatrix();
264 try {
265 final var invMatrix = com.irurueta.algebra.Utils.inverse(quadricMatrix);
266
267 // ensure that resulting matrix after inversion is symmetric
268 // by computing the mean of off-diagonal elements
269 final var a = invMatrix.getElementAt(0, 0);
270 final var b = invMatrix.getElementAt(1, 1);
271 final var c = invMatrix.getElementAt(2, 2);
272 final var d = 0.5 * (invMatrix.getElementAt(0, 1) + invMatrix.getElementAt(1, 0));
273 final var e = 0.5 * (invMatrix.getElementAt(2, 1) + invMatrix.getElementAt(1, 2));
274 final var f = 0.5 * (invMatrix.getElementAt(2, 0) + invMatrix.getElementAt(0, 2));
275 final var g = 0.5 * (invMatrix.getElementAt(3, 0) + invMatrix.getElementAt(0, 3));
276 final var h = 0.5 * (invMatrix.getElementAt(3, 1) + invMatrix.getElementAt(1, 3));
277 final var i = 0.5 * (invMatrix.getElementAt(3, 2) + invMatrix.getElementAt(2, 3));
278 final var j = invMatrix.getElementAt(3, 3);
279 dualQuadric.setParameters(a, b, c, d, e, f, g, h, i, j);
280 } catch (final AlgebraException e) {
281 throw new DualQuadricNotAvailableException(e);
282 }
283 }
284
285 /**
286 * Computes the dual quadric of this quadric.
287 * The dual quadric is equal to the inverse of the quadric matrix
288 *
289 * @return A new DualQuadric corresponding to the dual quadric of this
290 * instance.
291 * @throws DualQuadricNotAvailableException Raised if the dual quadric does
292 * not exist because this quadric instance is degenerate (its inverse cannot
293 * be computed).
294 */
295 public DualQuadric getDualQuadric() throws DualQuadricNotAvailableException {
296 final var dualQuadric = new DualQuadric();
297 dualQuadric(dualQuadric);
298 return dualQuadric;
299 }
300
301 /**
302 * Sets parameters of this quadric so that provided points lie within it
303 * (are locus).
304 *
305 * @param point1 1st point.
306 * @param point2 2nd point.
307 * @param point3 3rd point.
308 * @param point4 4th point.
309 * @param point5 5th point.
310 * @param point6 6th point.
311 * @param point7 7th point.
312 * @param point8 8th point.
313 * @param point9 9th point.
314 * @throws CoincidentPointsException Raised if points are coincident or
315 * produce a degenerated configuration.
316 */
317 public final void setParametersFromPoints(
318 final Point3D point1, final Point3D point2, final Point3D point3, final Point3D point4,
319 final Point3D point5, final Point3D point6, final Point3D point7, final Point3D point8,
320 final Point3D point9) throws CoincidentPointsException {
321
322 // normalize points to increase accuracy
323 point1.normalize();
324 point2.normalize();
325 point3.normalize();
326 point4.normalize();
327 point5.normalize();
328 point6.normalize();
329 point7.normalize();
330 point8.normalize();
331 point9.normalize();
332
333 try {
334 // each point belonging to a quadric follows equation:
335 // p' * Q * p = 0 ==>
336 // x^2 + y^2 + z^2 + 2*x*y + 2*x*z + 2*y*z + 2*x*w + 2*y*w +
337 // 2*z*w + w^2 = 0
338
339 final var m = new Matrix(9, 10);
340 var x = point1.getHomX();
341 var y = point1.getHomY();
342 var z = point1.getHomZ();
343 var w = point1.getHomW();
344 m.setElementAt(0, 0, x * x);
345 m.setElementAt(0, 1, y * y);
346 m.setElementAt(0, 2, z * z);
347 m.setElementAt(0, 3, 2.0 * x * y);
348 m.setElementAt(0, 4, 2.0 * x * z);
349 m.setElementAt(0, 5, 2.0 * y * z);
350 m.setElementAt(0, 6, 2.0 * x * w);
351 m.setElementAt(0, 7, 2.0 * y * w);
352 m.setElementAt(0, 8, 2.0 * z * w);
353 m.setElementAt(0, 9, w * w);
354 x = point2.getHomX();
355 y = point2.getHomY();
356 z = point2.getHomZ();
357 w = point2.getHomW();
358 m.setElementAt(1, 0, x * x);
359 m.setElementAt(1, 1, y * y);
360 m.setElementAt(1, 2, z * z);
361 m.setElementAt(1, 3, 2.0 * x * y);
362 m.setElementAt(1, 4, 2.0 * x * z);
363 m.setElementAt(1, 5, 2.0 * y * z);
364 m.setElementAt(1, 6, 2.0 * x * w);
365 m.setElementAt(1, 7, 2.0 * y * w);
366 m.setElementAt(1, 8, 2.0 * z * w);
367 m.setElementAt(1, 9, w * w);
368 x = point3.getHomX();
369 y = point3.getHomY();
370 z = point3.getHomZ();
371 w = point3.getHomW();
372 m.setElementAt(2, 0, x * x);
373 m.setElementAt(2, 1, y * y);
374 m.setElementAt(2, 2, z * z);
375 m.setElementAt(2, 3, 2.0 * x * y);
376 m.setElementAt(2, 4, 2.0 * x * z);
377 m.setElementAt(2, 5, 2.0 * y * z);
378 m.setElementAt(2, 6, 2.0 * x * w);
379 m.setElementAt(2, 7, 2.0 * y * w);
380 m.setElementAt(2, 8, 2.0 * z * w);
381 m.setElementAt(2, 9, w * w);
382 x = point4.getHomX();
383 y = point4.getHomY();
384 z = point4.getHomZ();
385 w = point4.getHomW();
386 m.setElementAt(3, 0, x * x);
387 m.setElementAt(3, 1, y * y);
388 m.setElementAt(3, 2, z * z);
389 m.setElementAt(3, 3, 2.0 * x * y);
390 m.setElementAt(3, 4, 2.0 * x * z);
391 m.setElementAt(3, 5, 2.0 * y * z);
392 m.setElementAt(3, 6, 2.0 * x * w);
393 m.setElementAt(3, 7, 2.0 * y * w);
394 m.setElementAt(3, 8, 2.0 * z * w);
395 m.setElementAt(3, 9, w * w);
396 x = point5.getHomX();
397 y = point5.getHomY();
398 z = point5.getHomZ();
399 w = point5.getHomW();
400 m.setElementAt(4, 0, x * x);
401 m.setElementAt(4, 1, y * y);
402 m.setElementAt(4, 2, z * z);
403 m.setElementAt(4, 3, 2.0 * x * y);
404 m.setElementAt(4, 4, 2.0 * x * z);
405 m.setElementAt(4, 5, 2.0 * y * z);
406 m.setElementAt(4, 6, 2.0 * x * w);
407 m.setElementAt(4, 7, 2.0 * y * w);
408 m.setElementAt(4, 8, 2.0 * z * w);
409 m.setElementAt(4, 9, w * w);
410 x = point6.getHomX();
411 y = point6.getHomY();
412 z = point6.getHomZ();
413 w = point6.getHomW();
414 m.setElementAt(5, 0, x * x);
415 m.setElementAt(5, 1, y * y);
416 m.setElementAt(5, 2, z * z);
417 m.setElementAt(5, 3, 2.0 * x * y);
418 m.setElementAt(5, 4, 2.0 * x * z);
419 m.setElementAt(5, 5, 2.0 * y * z);
420 m.setElementAt(5, 6, 2.0 * x * w);
421 m.setElementAt(5, 7, 2.0 * y * w);
422 m.setElementAt(5, 8, 2.0 * z * w);
423 m.setElementAt(5, 9, w * w);
424 x = point7.getHomX();
425 y = point7.getHomY();
426 z = point7.getHomZ();
427 w = point7.getHomW();
428 m.setElementAt(6, 0, x * x);
429 m.setElementAt(6, 1, y * y);
430 m.setElementAt(6, 2, z * z);
431 m.setElementAt(6, 3, 2.0 * x * y);
432 m.setElementAt(6, 4, 2.0 * x * z);
433 m.setElementAt(6, 5, 2.0 * y * z);
434 m.setElementAt(6, 6, 2.0 * x * w);
435 m.setElementAt(6, 7, 2.0 * y * w);
436 m.setElementAt(6, 8, 2.0 * z * w);
437 m.setElementAt(6, 9, w * w);
438 x = point8.getHomX();
439 y = point8.getHomY();
440 z = point8.getHomZ();
441 w = point8.getHomW();
442 m.setElementAt(7, 0, x * x);
443 m.setElementAt(7, 1, y * y);
444 m.setElementAt(7, 2, z * z);
445 m.setElementAt(7, 3, 2.0 * x * y);
446 m.setElementAt(7, 4, 2.0 * x * z);
447 m.setElementAt(7, 5, 2.0 * y * z);
448 m.setElementAt(7, 6, 2.0 * x * w);
449 m.setElementAt(7, 7, 2.0 * y * w);
450 m.setElementAt(7, 8, 2.0 * z * w);
451 m.setElementAt(7, 9, w * w);
452 x = point9.getHomX();
453 y = point9.getHomY();
454 z = point9.getHomZ();
455 w = point9.getHomW();
456 m.setElementAt(8, 0, x * x);
457 m.setElementAt(8, 1, y * y);
458 m.setElementAt(8, 2, z * z);
459 m.setElementAt(8, 3, 2.0 * x * y);
460 m.setElementAt(8, 4, 2.0 * x * z);
461 m.setElementAt(8, 5, 2.0 * y * z);
462 m.setElementAt(8, 6, 2.0 * x * w);
463 m.setElementAt(8, 7, 2.0 * y * w);
464 m.setElementAt(8, 8, 2.0 * z * w);
465 m.setElementAt(8, 9, w * w);
466
467 // normalize each row to increase accuracy
468 final var row = new double[10];
469 double rowNorm;
470 for (var j = 0; j < 9; j++) {
471 m.getSubmatrixAsArray(j, 0, j, 9, row);
472 rowNorm = com.irurueta.algebra.Utils.normF(row);
473 for (var i = 0; i < 10; i++) {
474 m.setElementAt(j, i, m.getElementAt(j, i) / rowNorm);
475 }
476 }
477
478 final var decomposer = new SingularValueDecomposer(m);
479 decomposer.decompose();
480
481 if (decomposer.getRank() < 9) {
482 throw new CoincidentPointsException();
483 }
484
485 // the right null-space of m contains the parameters a, b, c, d, e ,f
486 // of the conic
487 final var v = decomposer.getV();
488
489 final var a = v.getElementAt(0, 9);
490 final var b = v.getElementAt(1, 9);
491 final var c = v.getElementAt(2, 9);
492 final var d = v.getElementAt(3, 9);
493
494 final var f = v.getElementAt(4, 9);
495 final var e = v.getElementAt(5, 9);
496
497 final var g = v.getElementAt(6, 9);
498 final var h = v.getElementAt(7, 9);
499 final var i = v.getElementAt(8, 9);
500 final var j = v.getElementAt(9, 9);
501
502 setParameters(a, b, c, d, e, f, g, h, i, j);
503 } catch (final AlgebraException ex) {
504 throw new CoincidentPointsException(ex);
505 }
506 }
507
508 /**
509 * Returns a plane tangent to this quadric at provided point, as long as
510 * the provided point is locus of this quadric.
511 *
512 * @param point point where plane must be tangent to quadric.
513 * @return a plane tangent to this quadric.
514 * @throws NotLocusException if provided point is not locus of this quadric.
515 */
516 public Plane getTangentPlaneAt(final Point3D point) throws NotLocusException {
517 final var plane = new Plane();
518 tangentPlaneAt(point, plane, DEFAULT_LOCUS_THRESHOLD);
519 return plane;
520
521 }
522
523 /**
524 * Computes a plane tangent to this quadric at provided point, as long as
525 * the provided point is locus of this quadric up to provided threshold.
526 *
527 * @param point point where plane must be tangent to quadric.
528 * @param plane plane where computed result will be stored.
529 * @param threshold threshold to determine if provided point is locus or not
530 * of this quadric. Usually this is a small value close to zero.
531 * @throws NotLocusException if provided point is not locus of this quadric.
532 * @throws IllegalArgumentException if provided threshold is negative.
533 */
534 public void tangentPlaneAt(final Point3D point, final Plane plane, final double threshold)
535 throws NotLocusException {
536
537 if (!isLocus(point, threshold)) {
538 throw new NotLocusException();
539 }
540
541 point.normalize();
542 normalize();
543
544 final var q = asMatrix();
545
546 try {
547 final var p = new Matrix(Point3D.POINT3D_HOMOGENEOUS_COORDINATES_LENGTH, 1);
548 p.setElementAt(0, 0, point.getHomX());
549 p.setElementAt(1, 0, point.getHomY());
550 p.setElementAt(2, 0, point.getHomZ());
551 p.setElementAt(3, 0, point.getHomW());
552
553 q.multiply(p);
554 } catch (final WrongSizeException ignore) {
555 // never happens
556 }
557
558 plane.setParameters(q.getElementAt(0, 0), q.getElementAt(1, 0),
559 q.getElementAt(2, 0), q.getElementAt(3, 0));
560 }
561
562 /**
563 * Intersects this quadric with provided plane.
564 * Notice that result of intersection is expressed on original quadric
565 * coordinates.
566 * If resulting conic needs to be expressed in terms of plane coordinates,
567 * then the plane and the conic must be rotated so that the plane becomes
568 * an xy-plane.
569 *
570 * @param plane plane to intersect this quadric with.
571 * @param result instance where resulting intersection will be stored.
572 */
573 public void intersectWith(final Plane plane, final Conic result) {
574 // A plane follows expression: A*x + B*y + C*z + D*w = 0
575
576 // A quadric has the following matrix form:
577 // Q = [A D F G]
578 // [D B E H]
579 // [F E C I]
580 // [G H I J]
581 // [x y z w][A D F G][x] = [x y z w][A*x + D*y + F*z + G*w] =
582 // [D B E H][y] [D*x + B*y + E*z + H*w]
583 // [F E C I][z] [F*x + E*y + C*z + I*w]
584 // [G H I J][w] [G*x + H*y + I*z + J*w]
585 // = A*x^2 + D*x*y + F*x*z + G*x*w + D*x*y + B*y^2 + E*y*z + H*y*w +
586 // F*x*z + E*y*z + C*z^2 + I*z*w + G*x*w + H*y*w + I*z*w + J*w^2 =
587 // = A*x^2 + B*y^2 + C*z^2 + 2*D*x*y + 2*E*y*z + 2*F*x*z + 2*G*x*w + 2*H*y*w + 2*I*z*w + J*w^2
588 // which follows expression:
589 // A*x^2 + B*y^2 + C*z^2 + 2*D*x*y + 2*E*y*z + 2*F*x*z + 2*G*x*w + 2*H*y*w + 2*I*z*w + J*w^2 = 0
590
591 // A conic has the following matrix form:
592 // C = [A B D]
593 // [B C E]
594 // [D E F]
595 // [x y w][A B D][x] = [x y w][A*x + B*y + D*w] =
596 // [B C E][y] [B*x + C*y + E*w]
597 // [D E F][w] [D*x + E*y + F*w]
598 // = A*x^2 + B*x*y + D*x*w + B*x*y + C*y^2 + E*y*w + D*x*w + E*y*w + F*w^2 =
599 // = A*x^2 + 2*B*x*y + C*y^2 + 2*D*x*w + 2*E*y*w + F*w^2
600 // which follows expression:
601 // A*x^2 + 2*B*x*y + C*y^2 + 2*D*x*w + 2*E*y*w + F*w^2 = 0
602
603 // Quadric parameters
604 final var aQ = getA();
605 final var bQ = getB();
606 final var cQ = getC();
607 final var dQ = getD();
608 final var eQ = getE();
609 final var fQ = getF();
610 final var gQ = getG();
611 final var hQ = getH();
612 final var iQ = getI();
613 final var jQ = getJ();
614
615 // Plane parameters
616 final var aP = plane.getA();
617 final var bP = plane.getB();
618 final var cP = plane.getC();
619 final var dP = plane.getD();
620
621 // we solve the following system of equations:
622 // aQ*x^2 + bQ*y^2 + cQ*z^2 + 2*dQ*x*y + 2*eQ*y*z + 2*fQ*x*z + 2*gQ*x*w + 2*hQ*y*w + 2*iQ*z*w + jQ*w^2 = 0
623 // aP*x + bP*y + cP*z + dP*w = 0
624
625 // Isolating z in plane equation:
626 // z = (- aP*x - bP*y - dP*w)/cP
627
628 // and substituting in quadric equation:
629 // aQ*x^2 + bQ*y^2 + cQ*(- aP*x - bP*y - dP*w)^2/cP^2 + 2*dQ*x*y +
630 // 2*eQ*y*(- aP*x - bP*y - dP*w)/cP + 2*fQ*x*(- aP*x - bP*y - dP*w)/cP +
631 // 2*gQ*x*w + 2*hQ*y*w + 2*iQ*(- aP*x - bP*y - dP*w)/cP*w + jQ*w^2 = 0
632
633 // aQ*x^2 + bQ*y^2 + cQ*(- aP*x - bP*y - dP*w)^2/cP^2 + 2*dQ*x*y +
634 // -2*eQ*aP/cP*x*y -2*eQ*bP/cP*y^2 -2*eQ*dP/cP*y*w +
635 // -2*fQ*aP/cP*x^2 -2*fQ*bP/cP*x*y -2*fQ*dP/cP*x*w +
636 // 2*gQ*x*w + 2*hQ*y*w +
637 // -2*iQ*aP/cP*x*w -2*iQ*bP/cP*y*w -2*iQ*dP/cP*w^2 +
638 // jQ*w^2 = 0
639
640 // aQ*x^2 + bQ*y^2 +
641 // aP^2*cQ/cP^2*x^2 + 2*aP*bP*cQ/cP^2*x*y + 2*aP*cQ*dP/cP^2*x*w +
642 // bP^2*cQ/cP^2*y^2 + 2*bP*cQ*dP/cP^2*y*w + cQ*dP^2/cP^2*w^2 +
643 // 2*dQ*x*y +
644 // -2*eQ*aP/cP*x*y -2*eQ*bP/cP*y^2 -2*eQ*dP/cP*y*w +
645 // -2*fQ*aP/cP*x^2 -2*fQ*bP/cP*x*y -2*fQ*dP/cP*x*w +
646 // 2*gQ*x*w + 2*hQ*y*w +
647 // -2*iQ*aP/cP*x*w -2*iQ*bP/cP*y*w -2*iQ*dP/cP*w^2 +
648 // jQ*w^2 = 0
649
650 // (aQ + aP^2*cQ/cP^2 -2*fQ*aP/cP)*x^2 +
651 // (2*aP*bP*cQ/cP^2 + 2*dQ -2*eQ*aP/cP -2*fQ*bP/cP)*x*y +
652 // (bQ + bP^2*cQ/cP^2 -2*eQ*bP/cP)*y^2 +
653 // (2*aP*cQ*dP/cP^2 -2*fQ*dP/cP + 2*gQ -2*iQ*aP/cP)*x*w +
654 // (2*bP*cQ*dP/cP^2 -2*eQ*dP/cP + 2*hQ -2*iQ*bP/cP)*y*w +
655 // (cQ*dP^2/cP^2 -2*iQ*dP/cP + jQ)*w^2 = 0
656
657
658 // (aQ - 2*aP*fQ/cP + cQ*aP^2/cP^2)*x^2 +
659 // 2*(dQ - bP*fQ/cP - eQ*aP/cP + aP*bP*cQ/cP^2)*x*y +
660 // 2*(gQ - dP*fQ/cP - aP*iQ/cP + cQ*aP*dP/cP^2)*x*w +
661 // (bQ - 2*bP*eQ/cP + cQ*bP^2/cP^2)*y^2 +
662 // 2*(hQ - bP*iQ/cP - eQ*dP/cP + cQ*bP*dP/cP^2)*y*w +
663 // (jQ - 2*dP*iQ/cP + cQ*dP^2/cP^2)*w^2 = 0
664
665
666 // Comparing with conic equation:
667 // aC*x^2 + 2*bC*x*y + 2*dC*x*w + cC*y^2 + 2*eC*y*w + fC*w^2 = 0
668
669 // then conic parameters become:
670 // aC = aQ - 2*aP*fQ/cP + cQ*aP^2/cP^2
671 // bC = dQ - bP*fQ/cP - eQ*aP/cP + aP*bP*cQ/cP^2
672 // cC = bQ - 2*bP*eQ/cP + cQ*bP^2/cP^2
673 // dC = gQ - dP*fQ/cP - aP*iQ/cP + cQ*aP*dP/cP^2
674 // eC = hQ - bP*iQ/cP - eQ*dP/cP + cQ*bP*dP/cP^2
675 // fC = jQ - 2*dP*iQ/cP + cQ*dP^2/cP^2
676
677 final var aP2 = aP * aP;
678 final var bP2 = bP * bP;
679 final var cP2 = cP * cP;
680 final var dP2 = dP * dP;
681
682 final var aC = aQ - 2.0 * aP * fQ / cP + cQ * aP2 / cP2;
683 final var bC = dQ - bP * fQ / cP - eQ * aP / cP + aP * bP * cQ / cP2;
684 final var cC = bQ - 2.0 * bP * eQ / cP + cQ * bP2 / cP2;
685 final var dC = gQ - dP * fQ / cP - aP * iQ / cP + cQ * aP * dP / cP2;
686 final var eC = hQ - bP * iQ / cP - eQ * dP / cP + cQ * bP * dP / cP2;
687 final var fC = jQ - 2.0 * dP * iQ / cP + cQ * dP2 / cP2;
688
689 result.setParameters(aC, bC, cC, dC, eC, fC);
690 }
691
692 /**
693 * Intersects this quadric with provided plane.
694 *
695 * @param plane plane to intersect this quadric with.
696 * @return conic resulting from the intersection.
697 */
698 public Conic intersectWith(final Plane plane) {
699 final var result = new Conic();
700 intersectWith(plane, result);
701 return result;
702 }
703
704 //TODO: shorted distance of point to quadric
705 //TODO: closest point to quadric
706 }