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
21 import java.io.Serializable;
22
23 /**
24 * This class defines a sphere.
25 */
26 @SuppressWarnings("WeakerAccess")
27 public class Sphere implements Serializable {
28
29 /**
30 * Constant defining minimum allowed radius.
31 */
32 public static final double MIN_RADIUS = 0.0;
33
34 /**
35 * Constant defining default threshold value used when none is provided.
36 */
37 public static final double DEFAULT_THRESHOLD = 1e-9;
38
39 /**
40 * Constant defining minimum allowed threshold.
41 */
42 public static final double MIN_THRESHOLD = 0.0;
43
44 /**
45 * Constant defining machine precision.
46 */
47 public static final double EPS = 1e-12;
48
49 /**
50 * Center of sphere.
51 */
52 private Point3D center;
53
54 /**
55 * Radius of sphere.
56 */
57 private double radius;
58
59 /**
60 * Empty constructor.
61 * Creates sphere located at space origin (0,0) with radius 1.0.
62 */
63 public Sphere() {
64 center = Point3D.create();
65 radius = 1.0;
66 }
67
68 /**
69 * Constructor.
70 * Sets center and radius of sphere.
71 *
72 * @param center Center of sphere.
73 * @param radius Radius of sphere.
74 * @throws IllegalArgumentException Raised if provided radius is negative.
75 */
76 public Sphere(final Point3D center, final double radius) {
77 setCenterAndRadius(center, radius);
78 }
79
80 /**
81 * Constructor.
82 * Computes a sphere by using 4 points that must belong to its locus.
83 *
84 * @param point1 point 1.
85 * @param point2 point 2.
86 * @param point3 point 3.
87 * @param point4 point 4.
88 * @throws CoplanarPointsException if provided set of points are coincident
89 * or coplanar (form a single plane). In such cases a singularity occurs
90 * since a sphere having an infinite radius would be required to contain all
91 * four points in its locus.
92 */
93 public Sphere(final Point3D point1, final Point3D point2, final Point3D point3, final Point3D point4)
94 throws CoplanarPointsException {
95 setParametersFromPoints(point1, point2, point3, point4);
96 }
97
98 /**
99 * Constructor.
100 * Computes a sphere from a valid quadric corresponding to a sphere
101 *
102 * @param quadric a quadric to create a sphere from.
103 * @throws IllegalArgumentException if provided quadric is not a sphere.
104 */
105 public Sphere(final Quadric quadric) {
106 setFromQuadric(quadric);
107 }
108
109 /**
110 * Returns center of sphere.
111 *
112 * @return Center of sphere.
113 */
114 public Point3D getCenter() {
115 return center;
116 }
117
118 /**
119 * Sets center of sphere.
120 *
121 * @param center Center of sphere.
122 * @throws NullPointerException Raised if provided center is null.
123 */
124 public void setCenter(final Point3D center) {
125 if (center == null) {
126 throw new NullPointerException();
127 }
128
129 this.center = center;
130 }
131
132 /**
133 * Returns radius of sphere.
134 *
135 * @return Radius of sphere.
136 */
137 public double getRadius() {
138 return radius;
139 }
140
141 /**
142 * Sets radius of sphere.
143 *
144 * @param radius Radius of sphere.
145 * @throws IllegalArgumentException Raised if provided radius is negative.
146 */
147 public void setRadius(final double radius) {
148 if (radius < MIN_RADIUS) {
149 throw new IllegalArgumentException();
150 }
151
152 this.radius = radius;
153 }
154
155 /**
156 * Sets center and radius of this sphere.
157 *
158 * @param center Center to be set.
159 * @param radius Radius to be set.
160 * @throws IllegalArgumentException Raised if provided radius is negative.
161 * @throws NullPointerException Raised if provided center is null.
162 */
163 public final void setCenterAndRadius(final Point3D center, final double radius) {
164 setRadius(radius);
165 setCenter(center);
166 }
167
168 /**
169 * Sets parameters of a sphere by using four points that must belong to its
170 * locus.
171 *
172 * @param point1 point 1.
173 * @param point2 point 2.
174 * @param point3 point 3.
175 * @param point4 point 4.
176 * @throws CoplanarPointsException if provided set of points are coincident
177 * or coplanar (form a single plane). In such cases a singularity occurs
178 * since a sphere having an infinite radius would be required to contain all
179 * four points in its locus.
180 */
181 public final void setParametersFromPoints(
182 final Point3D point1, final Point3D point2, final Point3D point3, final Point3D point4)
183 throws CoplanarPointsException {
184
185 // normalize points to increase accuracy
186 point1.normalize();
187 point2.normalize();
188 point3.normalize();
189 point4.normalize();
190
191 try {
192 final var m = new Matrix(4, 4);
193 final var b = new double[4];
194
195 // 1st point
196 var x = point1.getHomX();
197 var y = point1.getHomY();
198 var z = point1.getHomZ();
199 var w = point1.getHomW();
200 m.setElementAt(0, 0, 2.0 * x * w);
201 m.setElementAt(0, 1, 2.0 * y * w);
202 m.setElementAt(0, 2, 2.0 * z * w);
203 m.setElementAt(0, 3, w * w);
204 b[0] = -x * x - y * y - z * z;
205
206 // 2nd point
207 x = point2.getHomX();
208 y = point2.getHomY();
209 z = point2.getHomZ();
210 w = point2.getHomW();
211 m.setElementAt(1, 0, 2.0 * x * w);
212 m.setElementAt(1, 1, 2.0 * y * w);
213 m.setElementAt(1, 2, 2.0 * z * w);
214 m.setElementAt(1, 3, w * w);
215 b[1] = -x * x - y * y - z * z;
216
217 // 3rd point
218 x = point3.getHomX();
219 y = point3.getHomY();
220 z = point3.getHomZ();
221 w = point3.getHomW();
222 m.setElementAt(2, 0, 2.0 * x * w);
223 m.setElementAt(2, 1, 2.0 * y * w);
224 m.setElementAt(2, 2, 2.0 * z * w);
225 m.setElementAt(2, 3, w * w);
226 b[2] = -x * x - y * y - z * z;
227
228 // 4th point
229 x = point4.getHomX();
230 y = point4.getHomY();
231 z = point4.getHomZ();
232 w = point4.getHomW();
233 m.setElementAt(3, 0, 2.0 * x * w);
234 m.setElementAt(3, 1, 2.0 * y * w);
235 m.setElementAt(3, 2, 2.0 * z * w);
236 m.setElementAt(3, 3, w * w);
237 b[3] = -x * x - y * y - z * z;
238
239 // normalize each row to increase accuracy
240 final var row = new double[4];
241 double rowNorm;
242
243 for (var j = 0; j < 4; j++) {
244 m.getSubmatrixAsArray(j, 0, j, 3, row);
245 rowNorm = com.irurueta.algebra.Utils.normF(row);
246 for (var i = 0; i < 4; i++) {
247 m.setElementAt(j, i, m.getElementAt(j, i) / rowNorm);
248 }
249 b[j] /= rowNorm;
250 }
251
252 final var params = com.irurueta.algebra.Utils.solve(m, b);
253
254 // g = -cx
255 final var g = params[0];
256 // h = -cy
257 final var h = params[1];
258 // i = -cz
259 final var i = params[2];
260 // j = cx^2 + cy^2 + cz^2 - R^2
261 final var j = params[3];
262
263 // compute center
264 final var inhomCx = -g;
265 final var inhomCy = -h;
266 final var inhomCz = -i;
267 final var c = new InhomogeneousPoint3D(inhomCx, inhomCy, inhomCz);
268
269 // compute radius
270 final var r = Math.sqrt(inhomCx * inhomCx + inhomCy * inhomCy + inhomCz * inhomCz - j);
271
272 setCenterAndRadius(c, r);
273 } catch (final AlgebraException e) {
274 throw new CoplanarPointsException(e);
275 }
276 }
277
278 /**
279 * Returns volume of a sphere having provided radius.
280 *
281 * @param radius Radius of a sphere.
282 * @return Volume of a sphere having provided radius.
283 * @throws IllegalArgumentException Raised if provided radius is negative.
284 */
285 public static double volume(final double radius) {
286 if (radius < MIN_RADIUS) {
287 throw new IllegalArgumentException();
288 }
289 return 4.0 / 3.0 * Math.PI * radius * radius * radius;
290 }
291
292 /**
293 * Returns volume of this sphere.
294 *
295 * @return Volume of this sphere.
296 */
297 public double getVolume() {
298 return volume(radius);
299 }
300
301 /**
302 * Returns surface of a sphere having provided radius.
303 *
304 * @param radius Radius of a sphere.
305 * @return Surface of a sphere having provided radius.
306 * @throws IllegalArgumentException Raised if provided radius is negative.
307 */
308 public static double surface(final double radius) {
309 if (radius < MIN_RADIUS) {
310 throw new IllegalArgumentException();
311 }
312 return 4.0 * Math.PI * radius * radius;
313 }
314
315 /**
316 * Returns surface of this sphere.
317 *
318 * @return Surface of this sphere.
319 */
320 public double getSurface() {
321 return surface(radius);
322 }
323
324 /**
325 * Determines if provided point is inside this sphere or not up to a certain
326 * threshold.
327 * If provided threshold is positive, the sphere behaves as if it was a
328 * larger sphere increased by threshold amount, if provided threshold is
329 * negative, the sphere behaves as if it was a smaller sphere decreased by
330 * threshold amount in radius.
331 *
332 * @param point Point to be checked.
333 * @param threshold Threshold to determine if point is inside or not
334 * @return True if point is considered to be inside this sphere, false
335 * otherwise.
336 */
337 public boolean isInside(final Point3D point, final double threshold) {
338 return point.distanceTo(center) - threshold <= radius;
339 }
340
341 /**
342 * Determines if provided point is inside this sphere or not.
343 *
344 * @param point Point to be checked.
345 * @return True if point is considered to be inside this sphere, false
346 * otherwise.
347 */
348 public boolean isInside(final Point3D point) {
349 return isInside(point, 0.0);
350 }
351
352 /**
353 * Returns distance from provided point to the closest point located in the
354 * sphere boundary.
355 * Returned distance will be negative when point is inside of sphere, and
356 * positive otherwise.
357 *
358 * @param point Point to be checked.
359 * @return Distance from point to sphere boundary.
360 */
361 public double getSignedDistance(final Point3D point) {
362 return signedDistance(this, point);
363 }
364
365 /**
366 * Returns distance from provided point to the closest point located in
367 * provided sphere boundary.
368 * Returned distance will be negative when point is inside of sphere, and
369 * positive otherwise.
370 *
371 * @param sphere A sphere.
372 * @param point Point to be checked.
373 * @return Distance from point to provided sphere boundary.
374 */
375 public static double signedDistance(final Sphere sphere, final Point3D point) {
376 return point.distanceTo(sphere.getCenter()) - sphere.getRadius();
377 }
378
379 /**
380 * Returns distance from provided point to the closest point located in the
381 * sphere boundary.
382 *
383 * @param point Point to be checked.
384 * @return Distance from point to sphere boundary.
385 */
386 public double getDistance(final Point3D point) {
387 return Math.abs(getSignedDistance(point));
388 }
389
390 /**
391 * Returns distance from provided point to the closest point located in
392 * provided sphere boundary.
393 *
394 * @param sphere A sphere.
395 * @param point Point to be checked.
396 * @return Distance from point to provided sphere boundary.
397 */
398 public static double distance(final Sphere sphere, final Point3D point) {
399 return Math.abs(signedDistance(sphere, point));
400 }
401
402 /**
403 * Returns closest point to provided point that is located in this sphere
404 * boundary.
405 *
406 * @param point A point to be checked.
407 * @return Closest point laying in sphere boundary.
408 * @throws UndefinedPointException Raised if provided point is at sphere
409 * center or very close to it.
410 */
411 public Point3D getClosestPoint(final Point3D point) throws UndefinedPointException {
412 final var result = Point3D.create();
413 closestPoint(point, result);
414 return result;
415 }
416
417 /**
418 * Computes closest point to provided point that is located in this sphere
419 * boundary and stores the result in provided result instance.
420 *
421 * @param point A point to be checked.
422 * @param result Instance where result will be stored.
423 * @throws UndefinedPointException Raised if provided point is at sphere
424 * center or very close to ti.
425 */
426 public void closestPoint(final Point3D point, final Point3D result) throws UndefinedPointException {
427
428 var directionX = point.getInhomX() - center.getInhomX();
429 var directionY = point.getInhomY() - center.getInhomY();
430 var directionZ = point.getInhomZ() - center.getInhomZ();
431 // normalize direction and multiply by radius to set result as locus of
432 // circle
433 final var norm = Math.sqrt(directionX * directionX + directionY * directionY + directionZ * directionZ);
434
435 // check if point is at center or very close to center, in that case the
436 // closest point cannot be found (would be all points of a circle)
437 if (norm < EPS) {
438 throw new UndefinedPointException();
439 }
440
441 directionX *= radius / norm;
442 directionY *= radius / norm;
443 directionZ *= radius / norm;
444
445 result.setInhomogeneousCoordinates(center.getInhomX() + directionX,
446 center.getInhomY() + directionY, center.getInhomZ() + directionZ);
447 }
448
449 /**
450 * Determines whether provided point lies at sphere boundary or not up to
451 * a certain threshold.
452 *
453 * @param point Point to be checked.
454 * @param threshold A small threshold to determine whether point lies at
455 * sphere boundary.
456 * @return True if point lies at sphere boundary, false otherwise.
457 * @throws IllegalArgumentException Raised if provided threshold is negative.
458 */
459 public boolean isLocus(final Point3D point, final double threshold) {
460 if (threshold < MIN_THRESHOLD) {
461 throw new IllegalArgumentException();
462 }
463
464 return Math.abs(point.distanceTo(center) - radius) <= threshold;
465 }
466
467 /**
468 * Determines whether provided point lies at sphere boundary or not.
469 *
470 * @param point Point to be checked.
471 * @return True if point lies at sphere boundary, false otherwise.
472 */
473 public boolean isLocus(final Point3D point) {
474 return isLocus(point, DEFAULT_THRESHOLD);
475 }
476
477 /**
478 * Returns a plane tangent to this sphere at provided point. Provided point
479 * must be locus of this sphere, otherwise a NotLocusException will be
480 * thrown.
481 *
482 * @param point a locus point of this sphere.
483 * @return a 3D plane tangent to this sphere at provided point.
484 * @throws NotLocusException if provided point is not locus of this sphere
485 * up to DEFAULT_THRESHOLD.
486 */
487 public Plane getTangentPlaneAt(final Point3D point) throws NotLocusException {
488 return getTangentPlaneAt(point, DEFAULT_THRESHOLD);
489 }
490
491 /**
492 * Returns a plane tangent to this sphere at provided point. Provided point
493 * must be locus of this sphere, otherwise a NotLocusException will be
494 * thrown.
495 *
496 * @param point a locus point of this sphere.
497 * @param threshold threshold to determine if provided point is locus.
498 * @return a 3D plane tangent to this circle at provided point.
499 * @throws NotLocusException if provided point is not locus of this sphere
500 * up to provided threshold.
501 * @throws IllegalArgumentException if provided threshold is negative.
502 */
503 public Plane getTangentPlaneAt(final Point3D point, final double threshold) throws NotLocusException {
504 final var plane = new Plane();
505 tangentPlaneAt(point, plane, threshold);
506 return plane;
507 }
508
509 /**
510 * Computes a plane tangent to this sphere at provided point. Provided point
511 * must be locus of this sphere, otherwise a NotLocusException will be
512 * thrown.
513 *
514 * @param point a locus point of this sphere.
515 * @param plane instance of a 3D plane where result will be stored.
516 * @param threshold threshold to determine if provided point is locus.
517 * @throws NotLocusException if provided point is not locus of this sphere.
518 * @throws IllegalArgumentException if provided threshold is negative.
519 */
520 public void tangentPlaneAt(final Point3D point, final Plane plane, final double threshold)
521 throws NotLocusException {
522 if (!isLocus(point, threshold)) {
523 throw new NotLocusException();
524 }
525
526 point.normalize();
527 center.normalize();
528
529 // Q = [a d f g]
530 // [d b e h]
531 // [f e c i]
532 // [g h i j]
533
534 // Q = [1 0 0 -cx]
535 // [0 1 0 -cy]
536 // [0 0 1 -cz]
537 // [-cx -cy -cz cx^2 + cy^2 + cz^2 - r^2]
538
539 // a = b = c = 1.0
540 // d = e = f = 0.0
541 // g = -cx, h = -cy, i = -cz
542 // j = cx^2 + cy^2 + cz^2 - r^2
543
544 // Hence plane is P = Q * p, where Q is the sphere quadric and p is
545 // a point in the locus of the sphere
546
547 final var homX = point.getHomX();
548 final var homY = point.getHomY();
549 final var homZ = point.getHomZ();
550 final var homW = point.getHomW();
551 final var cx = center.getInhomX();
552 final var cy = center.getInhomY();
553 final var cz = center.getInhomZ();
554 final var quadricG = -cx;
555 final var quadricH = -cy;
556 final var quadricI = -cz;
557 final var quadricJ = cx * cx + cy * cy + cz * cz - radius * radius;
558
559 final var planeA = homX + quadricG * homW;
560 final var planeB = homY + quadricH * homW;
561 final var planeC = homZ + quadricI * homW;
562 final var planeD = quadricG * homX + quadricH * homY + quadricI * homZ + quadricJ * homW;
563 plane.setParameters(planeA, planeB, planeC, planeD);
564 }
565
566 /**
567 * Converts this sphere into a quadric.
568 * Quadrics are a more general representation of spheres.
569 *
570 * @return A quadric representing this circle.
571 */
572 public Quadric toQuadric() {
573 center.normalize();
574 // use inhomogeneous center coordinates
575 final var cx = center.getInhomX();
576 final var cy = center.getInhomY();
577 final var cz = center.getInhomZ();
578
579 final var a = 1.0;
580 final var b = 1.0;
581 final var c = 1.0;
582 final var d = 0.0;
583 final var e = 0.0;
584 final var f = 0.0;
585 final var g = -cx;
586 final var h = -cy;
587 final var i = -cz;
588 final var j = cx * cx + cy * cy + cz * cz - radius * radius;
589
590 return new Quadric(a, b, c, d, e, f, g, h, i, j);
591 }
592
593 /**
594 * Sets parameters of this sphere from a valid quadric corresponding to a
595 * sphere.
596 *
597 * @param quadric quadric to set parameters from.
598 * @throws IllegalArgumentException if provided quadric is not a sphere.
599 */
600 public final void setFromQuadric(final Quadric quadric) {
601 final var isSphere = quadric.getA() == quadric.getB() && quadric.getB() == quadric.getC()
602 && quadric.getA() != 0.0 && quadric.getD() == 0.0 && quadric.getE() == 0.0 && quadric.getF() == 0.0;
603
604 if (!isSphere) {
605 throw new IllegalArgumentException();
606 }
607
608 quadric.normalize();
609
610 final var a = quadric.getA();
611 // normalize parameters so that
612 // a = b = c = 1.0
613 // d = e = f = 0.0
614 // g = -cx, h = -cy, i = -cz
615 // j = cx^2 + cy^2 + cz^2 - r^2
616 final var normG = quadric.getG() / a;
617 final var normH = quadric.getH() / a;
618 final var normI = quadric.getI() / a;
619 final var normJ = quadric.getJ() / a;
620
621 final var cx = -normG;
622 final var cy = -normH;
623 final var cz = -normI;
624 final var r = Math.sqrt(cx * cx + cy * cy + cz * cz - normJ);
625
626 center = new InhomogeneousPoint3D(cx, cy, cz);
627 this.radius = r;
628 }
629 }