1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.geometry.euclidean.threed.shape;
18
19 import java.text.MessageFormat;
20 import java.util.List;
21 import java.util.stream.Collectors;
22 import java.util.stream.Stream;
23
24 import org.apache.commons.geometry.core.partitioning.bsp.RegionCutRule;
25 import org.apache.commons.geometry.euclidean.AbstractNSphere;
26 import org.apache.commons.geometry.euclidean.threed.Plane;
27 import org.apache.commons.geometry.euclidean.threed.Planes;
28 import org.apache.commons.geometry.euclidean.threed.RegionBSPTree3D;
29 import org.apache.commons.geometry.euclidean.threed.RegionBSPTree3D.RegionNode3D;
30 import org.apache.commons.geometry.euclidean.threed.Vector3D;
31 import org.apache.commons.geometry.euclidean.threed.line.Line3D;
32 import org.apache.commons.geometry.euclidean.threed.line.LineConvexSubset3D;
33 import org.apache.commons.geometry.euclidean.threed.line.LinecastPoint3D;
34 import org.apache.commons.geometry.euclidean.threed.line.Linecastable3D;
35 import org.apache.commons.geometry.euclidean.threed.mesh.SimpleTriangleMesh;
36 import org.apache.commons.geometry.euclidean.threed.mesh.TriangleMesh;
37 import org.apache.commons.numbers.core.Precision;
38
39 /** Class representing a 3 dimensional sphere in Euclidean space.
40 */
41 public final class Sphere extends AbstractNSphere<Vector3D> implements Linecastable3D {
42
43 /** Message used when requesting a sphere approximation with an invalid subdivision number. */
44 private static final String INVALID_SUBDIVISION_MESSAGE =
45 "Number of sphere approximation subdivisions must be greater than or equal to zero; was {0}";
46
47 /** Constant equal to {@code 4 * pi}. */
48 private static final double FOUR_PI = 4.0 * Math.PI;
49
50 /** Constant equal to {@code (4/3) * pi}. */
51 private static final double FOUR_THIRDS_PI = FOUR_PI / 3.0;
52
53 /** Construct a new sphere from its component parts.
54 * @param center the center of the sphere
55 * @param radius the sphere radius
56 * @param precision precision context used to compare floating point numbers
57 * @throws IllegalArgumentException if center is not finite or radius is not finite or is
58 * less than or equal to zero as evaluated by the given precision context
59 */
60 private Sphere(final Vector3D center, final double radius, final Precision.DoubleEquivalence precision) {
61 super(center, radius, precision);
62 }
63
64 /** {@inheritDoc} */
65 @Override
66 public double getSize() {
67 final double r = getRadius();
68 return FOUR_THIRDS_PI * r * r * r;
69 }
70
71 /** {@inheritDoc} */
72 @Override
73 public double getBoundarySize() {
74 final double r = getRadius();
75 return FOUR_PI * r * r;
76 }
77
78 /** {@inheritDoc} */
79 @Override
80 public Vector3D project(final Vector3D pt) {
81 return project(pt, Vector3D.Unit.PLUS_X);
82 }
83
84 /** Build an approximation of this sphere using a {@link RegionBSPTree3D}. The approximation is constructed by
85 * taking an octahedron (8-sided polyhedron with triangular faces) inscribed in the sphere and subdividing each
86 * triangular face {@code subdivisions} number of times, each time projecting the newly created vertices onto the
87 * sphere surface. Each triangle subdivision produces 4 triangles, meaning that the total number of triangles
88 * inserted into tree is equal to \(8 \times 4^s\), where \(s\) is the number of subdivisions. For
89 * example, calling this method with {@code subdivisions} equal to {@code 3} will produce a tree having
90 * \(8 \times 4^3 = 512\) triangular facets inserted. See the table below for other examples. The returned BSP
91 * tree also contains structural cuts to reduce the overall height of the tree.
92 *
93 * <table>
94 * <caption>Subdivisions to Triangle Counts</caption>
95 * <thead>
96 * <tr>
97 * <th>Subdivisions</th>
98 * <th>Triangles</th>
99 * </tr>
100 * </thead>
101 * <tbody>
102 * <tr><td>0</td><td>8</td></tr>
103 * <tr><td>1</td><td>32</td></tr>
104 * <tr><td>2</td><td>128</td></tr>
105 * <tr><td>3</td><td>512</td></tr>
106 * <tr><td>4</td><td>2048</td></tr>
107 * <tr><td>5</td><td>8192</td></tr>
108 * </tbody>
109 * </table>
110 *
111 * <p>Care must be taken when using this method with large subdivision numbers so that floating point errors
112 * do not interfere with the creation of the planes and triangles in the tree. For example, if the number of
113 * subdivisions is too high, the subdivided triangle points may become equivalent according to the sphere's
114 * {@link #getPrecision() precision context} and plane creation may fail. Or plane creation may succeed but
115 * insertion of the plane into the tree may fail for similar reasons. In general, it is best to use the lowest
116 * subdivision number practical for the intended purpose.</p>
117 * @param subdivisions the number of triangle subdivisions to use when creating the tree; the total number of
118 * triangular facets inserted into the returned tree is equal to \(8 \times 4^s\), where \(s\) is the number
119 * of subdivisions
120 * @return a BSP tree containing an approximation of the sphere
121 * @throws IllegalArgumentException if {@code subdivisions} is less than zero
122 * @throws IllegalStateException if tree creation fails for the given subdivision count
123 * @see #toTriangleMesh(int)
124 */
125 public RegionBSPTree3D toTree(final int subdivisions) {
126 if (subdivisions < 0) {
127 throw new IllegalArgumentException(MessageFormat.format(INVALID_SUBDIVISION_MESSAGE, subdivisions));
128 }
129 return new SphereTreeApproximationBuilder(this, subdivisions).build();
130 }
131
132 /** Build an approximation of this sphere using a {@link TriangleMesh}. The approximation is constructed by
133 * taking an octahedron (8-sided polyhedron with triangular faces) inscribed in the sphere and subdividing each
134 * triangular face {@code subdivisions} number of times, each time projecting the newly created vertices onto the
135 * sphere surface. Each triangle subdivision produces 4 triangles, meaning that the total number of triangles
136 * in the returned mesh is equal to \(8 \times 4^s\), where \(s\) is the number of subdivisions. For
137 * example, calling this method with {@code subdivisions} equal to {@code 3} will produce a mesh having
138 * \(8 \times 4^3 = 512\) triangular facets inserted. See the table below for other examples.
139 *
140 * <table>
141 * <caption>Subdivisions to Triangle Counts</caption>
142 * <thead>
143 * <tr>
144 * <th>Subdivisions</th>
145 * <th>Triangles</th>
146 * </tr>
147 * </thead>
148 * <tbody>
149 * <tr><td>0</td><td>8</td></tr>
150 * <tr><td>1</td><td>32</td></tr>
151 * <tr><td>2</td><td>128</td></tr>
152 * <tr><td>3</td><td>512</td></tr>
153 * <tr><td>4</td><td>2048</td></tr>
154 * <tr><td>5</td><td>8192</td></tr>
155 * </tbody>
156 * </table>
157 *
158 * <p><strong>BSP Tree Conversion</strong></p>
159 * <p>Inserting the boundaries of a sphere mesh approximation directly into a BSP tree will invariably result
160 * in poor performance: since the region is convex the constructed BSP tree degenerates into a simple linked
161 * list of nodes. If a BSP tree is needed, users should prefer the {@link #toTree(int)} method, which creates
162 * balanced tree approximations directly, or the {@link RegionBSPTree3D.PartitionedRegionBuilder3D} class,
163 * which can be used to insert the mesh faces into a pre-partitioned tree.
164 * </p>
165 * @param subdivisions the number of triangle subdivisions to use when creating the mesh; the total number of
166 * triangular faces in the returned mesh is equal to \(8 \times 4^s\), where \(s\) is the number
167 * of subdivisions
168 * @return a triangle mesh approximation of the sphere
169 * @throws IllegalArgumentException if {@code subdivisions} is less than zero
170 * @see #toTree(int)
171 */
172 public TriangleMesh toTriangleMesh(final int subdivisions) {
173 if (subdivisions < 0) {
174 throw new IllegalArgumentException(MessageFormat.format(INVALID_SUBDIVISION_MESSAGE, subdivisions));
175 }
176 return new SphereMeshApproximationBuilder(this, subdivisions).build();
177 }
178
179 /** Get the intersections of the given line with this sphere. The returned list will
180 * contain either 0, 1, or 2 points.
181 * <ul>
182 * <li><strong>2 points</strong> - The line is a secant line and intersects the sphere at two
183 * distinct points. The points are ordered such that the first point in the list is the first point
184 * encountered when traveling in the direction of the line. (In other words, the points are ordered
185 * by increasing abscissa value.)
186 * </li>
187 * <li><strong>1 point</strong> - The line is a tangent line and only intersects the sphere at a
188 * single point (as evaluated by the sphere's precision context).
189 * </li>
190 * <li><strong>0 points</strong> - The line does not intersect the sphere.</li>
191 * </ul>
192 * @param line line to intersect with the sphere
193 * @return a list of intersection points between the given line and this sphere
194 */
195 public List<Vector3D> intersections(final Line3D line) {
196 return intersections(line, Line3D::abscissa, Line3D::distance);
197 }
198
199 /** Get the first intersection point between the given line and this sphere, or null
200 * if no such point exists. The "first" intersection point is the first such point
201 * encountered when traveling in the direction of the line from infinity.
202 * @param line line to intersect with the sphere
203 * @return the first intersection point between the given line and this instance or
204 * null if no such point exists
205 */
206 public Vector3D firstIntersection(final Line3D line) {
207 return firstIntersection(line, Line3D::abscissa, Line3D::distance);
208 }
209
210 /** {@inheritDoc} */
211 @Override
212 public List<LinecastPoint3D> linecast(final LineConvexSubset3D subset) {
213 return getLinecastStream(subset)
214 .collect(Collectors.toList());
215 }
216
217 /** {@inheritDoc} */
218 @Override
219 public LinecastPoint3D linecastFirst(final LineConvexSubset3D subset) {
220 return getLinecastStream(subset)
221 .findFirst()
222 .orElse(null);
223 }
224
225 /** Get a stream containing the linecast intersection points of the given
226 * line subset with this instance.
227 * @param subset line subset to intersect against this instance
228 * @return a stream containing linecast intersection points
229 */
230 private Stream<LinecastPoint3D> getLinecastStream(final LineConvexSubset3D subset) {
231 return intersections(subset.getLine()).stream()
232 .filter(subset::contains)
233 .map(pt -> new LinecastPoint3D(pt, getCenter().directionTo(pt), subset.getLine()));
234 }
235
236 /** Construct a sphere from a center point and radius.
237 * @param center the center of the sphere
238 * @param radius the sphere radius
239 * @param precision precision context used to compare floating point numbers
240 * @return a sphere constructed from the given center point and radius
241 * @throws IllegalArgumentException if center is not finite or radius is not finite or is
242 * less than or equal to zero as evaluated by the given precision context
243 */
244 public static Sphere from(final Vector3D center, final double radius, final Precision.DoubleEquivalence precision) {
245 return new Sphere(center, radius, precision);
246 }
247
248 /** Internal class used to construct hyperplane-bounded approximations of spheres as BSP trees. The class
249 * begins with an octahedron inscribed in the sphere and then subdivides each triangular face a specified
250 * number of times.
251 */
252 private static final class SphereTreeApproximationBuilder {
253
254 /** Threshold used to determine when to stop inserting structural cuts and begin adding facets. */
255 private static final int PARTITION_THRESHOLD = 2;
256
257 /** The sphere that an approximation is being created for. */
258 private final Sphere sphere;
259
260 /** The number of triangular subdivisions to use. */
261 private final int subdivisions;
262
263 /** Construct a new builder for creating a BSP tree approximation of the given sphere.
264 * @param sphere the sphere to create an approximation of
265 * @param subdivisions the number of triangle subdivisions to use in tree creation
266 */
267 SphereTreeApproximationBuilder(final Sphere sphere, final int subdivisions) {
268 this.sphere = sphere;
269 this.subdivisions = subdivisions;
270 }
271
272 /** Build the sphere approximation BSP tree.
273 * @return the sphere approximation BSP tree
274 * @throws IllegalStateException if tree creation fails for the configured subdivision count
275 */
276 RegionBSPTree3D build() {
277 final RegionBSPTree3D tree = RegionBSPTree3D.empty();
278
279 final Vector3D center = sphere.getCenter();
280 final double radius = sphere.getRadius();
281 final Precision.DoubleEquivalence precision = sphere.getPrecision();
282
283 // insert the primary split planes
284 final Plane plusXPlane = Planes.fromPointAndNormal(center, Vector3D.Unit.PLUS_X, precision);
285 final Plane plusYPlane = Planes.fromPointAndNormal(center, Vector3D.Unit.PLUS_Y, precision);
286 final Plane plusZPlane = Planes.fromPointAndNormal(center, Vector3D.Unit.PLUS_Z, precision);
287
288 tree.insert(plusXPlane.span(), RegionCutRule.INHERIT);
289 tree.insert(plusYPlane.span(), RegionCutRule.INHERIT);
290 tree.insert(plusZPlane.span(), RegionCutRule.INHERIT);
291
292 // create the vertices for the octahedron
293 final double cx = center.getX();
294 final double cy = center.getY();
295 final double cz = center.getZ();
296
297 final Vector3D maxX = Vector3D.of(cx + radius, cy, cz);
298 final Vector3D minX = Vector3D.of(cx - radius, cy, cz);
299
300 final Vector3D maxY = Vector3D.of(cx, cy + radius, cz);
301 final Vector3D minY = Vector3D.of(cx, cy - radius, cz);
302
303 final Vector3D maxZ = Vector3D.of(cx, cy, cz + radius);
304 final Vector3D minZ = Vector3D.of(cx, cy, cz - radius);
305
306 // partition and subdivide the face triangles and insert them into the tree
307 final RegionNode3D root = tree.getRoot();
308
309 try {
310 partitionAndInsert(root.getMinus().getMinus().getMinus(), minX, minZ, minY, 0);
311 partitionAndInsert(root.getMinus().getMinus().getPlus(), minX, minY, maxZ, 0);
312
313 partitionAndInsert(root.getMinus().getPlus().getMinus(), minX, maxY, minZ, 0);
314 partitionAndInsert(root.getMinus().getPlus().getPlus(), minX, maxZ, maxY, 0);
315
316 partitionAndInsert(root.getPlus().getMinus().getMinus(), maxX, minY, minZ, 0);
317 partitionAndInsert(root.getPlus().getMinus().getPlus(), maxX, maxZ, minY, 0);
318
319 partitionAndInsert(root.getPlus().getPlus().getMinus(), maxX, minZ, maxY, 0);
320 partitionAndInsert(root.getPlus().getPlus().getPlus(), maxX, maxY, maxZ, 0);
321 } catch (final IllegalStateException | IllegalArgumentException exc) {
322 // standardize any tree construction failure as an IllegalStateException
323 throw new IllegalStateException("Failed to construct sphere approximation with subdivision count " +
324 subdivisions + ": " + exc.getMessage(), exc);
325 }
326
327 return tree;
328 }
329
330 /** Recursively insert structural BSP tree cuts into the given node and then insert subdivided triangles
331 * when a target subdivision level is reached. The structural BSP tree cuts are used to help reduce the
332 * overall depth of the BSP tree.
333 * @param node the node to insert into
334 * @param p1 first triangle point
335 * @param p2 second triangle point
336 * @param p3 third triangle point
337 * @param level current subdivision level
338 */
339 private void partitionAndInsert(final RegionNode3D node,
340 final Vector3D p1, final Vector3D p2, final Vector3D p3, final int level) {
341
342 if (subdivisions - level > PARTITION_THRESHOLD) {
343 final int nextLevel = level + 1;
344
345 final Vector3D center = sphere.getCenter();
346
347 final Vector3D m1 = sphere.project(p1.lerp(p2, 0.5));
348 final Vector3D m2 = sphere.project(p2.lerp(p3, 0.5));
349 final Vector3D m3 = sphere.project(p3.lerp(p1, 0.5));
350
351 RegionNode3D curNode = node;
352
353 checkedCut(curNode, createPlane(m3, m2, center), RegionCutRule.INHERIT);
354 partitionAndInsert(curNode.getPlus(), m3, m2, p3, nextLevel);
355
356 curNode = curNode.getMinus();
357 checkedCut(curNode, createPlane(m2, m1, center), RegionCutRule.INHERIT);
358 partitionAndInsert(curNode.getPlus(), m1, p2, m2, nextLevel);
359
360 curNode = curNode.getMinus();
361 checkedCut(curNode, createPlane(m1, m3, center), RegionCutRule.INHERIT);
362 partitionAndInsert(curNode.getPlus(), p1, m1, m3, nextLevel);
363
364 partitionAndInsert(curNode.getMinus(), m1, m2, m3, nextLevel);
365 } else {
366 insertSubdividedTriangles(node, p1, p2, p3, level);
367 }
368 }
369
370 /** Recursively insert subdivided triangles into the given node. Each triangle is inserted into the minus
371 * side of the previous triangle.
372 * @param node the node to insert into
373 * @param p1 first triangle point
374 * @param p2 second triangle point
375 * @param p3 third triangle point
376 * @param level the current subdivision level
377 * @return the node representing the inside of the region after insertion of all triangles
378 */
379 private RegionNode3D insertSubdividedTriangles(final RegionNode3D node,
380 final Vector3D p1, final Vector3D p2, final Vector3D p3,
381 final int level) {
382
383 if (level >= subdivisions) {
384 // base case
385 checkedCut(node, createPlane(p1, p2, p3), RegionCutRule.MINUS_INSIDE);
386 return node.getMinus();
387 } else {
388 final int nextLevel = level + 1;
389
390 final Vector3D m1 = sphere.project(p1.lerp(p2, 0.5));
391 final Vector3D m2 = sphere.project(p2.lerp(p3, 0.5));
392 final Vector3D m3 = sphere.project(p3.lerp(p1, 0.5));
393
394 RegionNode3D curNode = node;
395 curNode = insertSubdividedTriangles(curNode, p1, m1, m3, nextLevel);
396 curNode = insertSubdividedTriangles(curNode, m1, p2, m2, nextLevel);
397 curNode = insertSubdividedTriangles(curNode, m3, m2, p3, nextLevel);
398 curNode = insertSubdividedTriangles(curNode, m1, m2, m3, nextLevel);
399
400 return curNode;
401 }
402 }
403
404 /** Create a plane from the given points, using the precision context of the sphere.
405 * @param p1 first point
406 * @param p2 second point
407 * @param p3 third point
408 * @return a plane defined by the given points
409 */
410 private Plane createPlane(final Vector3D p1, final Vector3D p2, final Vector3D p3) {
411 return Planes.fromPoints(p1, p2, p3, sphere.getPrecision());
412 }
413
414 /** Insert the cut into the given node, throwing an exception if no portion of the cutter intersects
415 * the node.
416 * @param node node to cut
417 * @param cutter plane to use to cut the node
418 * @param cutRule cut rule to apply
419 * @throws IllegalStateException if no portion of the cutter plane intersects the node
420 */
421 private void checkedCut(final RegionNode3D node, final Plane cutter, final RegionCutRule cutRule) {
422 if (!node.insertCut(cutter, cutRule)) {
423 throw new IllegalStateException("Failed to cut BSP tree node with plane: " + cutter);
424 }
425 }
426 }
427
428 /** Internal class used to construct geodesic mesh sphere approximations. The class begins with an octahedron
429 * inscribed in the sphere and then subdivides each triangular face a specified number of times.
430 */
431 private static final class SphereMeshApproximationBuilder {
432
433 /** The sphere that an approximation is being created for. */
434 private final Sphere sphere;
435
436 /** The number of triangular subdivisions to use. */
437 private final int subdivisions;
438
439 /** Mesh builder object. */
440 private final SimpleTriangleMesh.Builder builder;
441
442 /** Construct a new builder for creating a mesh approximation of the given sphere.
443 * @param sphere the sphere to create an approximation of
444 * @param subdivisions the number of triangle subdivisions to use in mesh creation
445 */
446 SphereMeshApproximationBuilder(final Sphere sphere, final int subdivisions) {
447 this.sphere = sphere;
448 this.subdivisions = subdivisions;
449 this.builder = SimpleTriangleMesh.builder(sphere.getPrecision());
450 }
451
452 /** Build the mesh approximation of the configured sphere.
453 * @return the mesh approximation of the configured sphere
454 */
455 public SimpleTriangleMesh build() {
456 final Vector3D center = sphere.getCenter();
457 final double radius = sphere.getRadius();
458
459 // create the vertices for the octahedron
460 final double cx = center.getX();
461 final double cy = center.getY();
462 final double cz = center.getZ();
463
464 final Vector3D maxX = Vector3D.of(cx + radius, cy, cz);
465 final Vector3D minX = Vector3D.of(cx - radius, cy, cz);
466
467 final Vector3D maxY = Vector3D.of(cx, cy + radius, cz);
468 final Vector3D minY = Vector3D.of(cx, cy - radius, cz);
469
470 final Vector3D maxZ = Vector3D.of(cx, cy, cz + radius);
471 final Vector3D minZ = Vector3D.of(cx, cy, cz - radius);
472
473 addSubdivided(minX, minZ, minY, 0);
474 addSubdivided(minX, minY, maxZ, 0);
475
476 addSubdivided(minX, maxY, minZ, 0);
477 addSubdivided(minX, maxZ, maxY, 0);
478
479 addSubdivided(maxX, minY, minZ, 0);
480 addSubdivided(maxX, maxZ, minY, 0);
481
482 addSubdivided(maxX, minZ, maxY, 0);
483 addSubdivided(maxX, maxY, maxZ, 0);
484
485 return builder.build();
486 }
487
488 /** Recursively subdivide and add triangular faces between the given outer boundary points.
489 * @param p1 first point
490 * @param p2 second point
491 * @param p3 third point
492 * @param level recursion level; counts up
493 */
494 private void addSubdivided(final Vector3D p1, final Vector3D p2, final Vector3D p3, final int level) {
495 if (level >= subdivisions) {
496 // base case
497 builder.addFaceUsingVertices(p1, p2, p3);
498 } else {
499 // subdivide
500 final int nextLevel = level + 1;
501
502 final Vector3D m1 = sphere.project(p1.lerp(p2, 0.5));
503 final Vector3D m2 = sphere.project(p2.lerp(p3, 0.5));
504 final Vector3D m3 = sphere.project(p3.lerp(p1, 0.5));
505
506 addSubdivided(p1, m1, m3, nextLevel);
507 addSubdivided(m1, p2, m2, nextLevel);
508 addSubdivided(m3, m2, p3, nextLevel);
509 addSubdivided(m1, m2, m3, nextLevel);
510 }
511 }
512 }
513 }