Sphere.java
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.geometry.euclidean.threed.shape;
import java.text.MessageFormat;
import java.util.List;
import java.util.stream.Collectors;
import java.util.stream.Stream;
import org.apache.commons.geometry.core.partitioning.bsp.RegionCutRule;
import org.apache.commons.geometry.euclidean.AbstractNSphere;
import org.apache.commons.geometry.euclidean.threed.Plane;
import org.apache.commons.geometry.euclidean.threed.Planes;
import org.apache.commons.geometry.euclidean.threed.RegionBSPTree3D;
import org.apache.commons.geometry.euclidean.threed.RegionBSPTree3D.RegionNode3D;
import org.apache.commons.geometry.euclidean.threed.Vector3D;
import org.apache.commons.geometry.euclidean.threed.line.Line3D;
import org.apache.commons.geometry.euclidean.threed.line.LineConvexSubset3D;
import org.apache.commons.geometry.euclidean.threed.line.LinecastPoint3D;
import org.apache.commons.geometry.euclidean.threed.line.Linecastable3D;
import org.apache.commons.geometry.euclidean.threed.mesh.SimpleTriangleMesh;
import org.apache.commons.geometry.euclidean.threed.mesh.TriangleMesh;
import org.apache.commons.numbers.core.Precision;
/** Class representing a 3 dimensional sphere in Euclidean space.
*/
public final class Sphere extends AbstractNSphere<Vector3D> implements Linecastable3D {
/** Message used when requesting a sphere approximation with an invalid subdivision number. */
private static final String INVALID_SUBDIVISION_MESSAGE =
"Number of sphere approximation subdivisions must be greater than or equal to zero; was {0}";
/** Constant equal to {@code 4 * pi}. */
private static final double FOUR_PI = 4.0 * Math.PI;
/** Constant equal to {@code (4/3) * pi}. */
private static final double FOUR_THIRDS_PI = FOUR_PI / 3.0;
/** Construct a new sphere from its component parts.
* @param center the center of the sphere
* @param radius the sphere radius
* @param precision precision context used to compare floating point numbers
* @throws IllegalArgumentException if center is not finite or radius is not finite or is
* less than or equal to zero as evaluated by the given precision context
*/
private Sphere(final Vector3D center, final double radius, final Precision.DoubleEquivalence precision) {
super(center, radius, precision);
}
/** {@inheritDoc} */
@Override
public double getSize() {
final double r = getRadius();
return FOUR_THIRDS_PI * r * r * r;
}
/** {@inheritDoc} */
@Override
public double getBoundarySize() {
final double r = getRadius();
return FOUR_PI * r * r;
}
/** {@inheritDoc} */
@Override
public Vector3D project(final Vector3D pt) {
return project(pt, Vector3D.Unit.PLUS_X);
}
/** Build an approximation of this sphere using a {@link RegionBSPTree3D}. The approximation is constructed by
* taking an octahedron (8-sided polyhedron with triangular faces) inscribed in the sphere and subdividing each
* triangular face {@code subdivisions} number of times, each time projecting the newly created vertices onto the
* sphere surface. Each triangle subdivision produces 4 triangles, meaning that the total number of triangles
* inserted into tree is equal to \(8 \times 4^s\), where \(s\) is the number of subdivisions. For
* example, calling this method with {@code subdivisions} equal to {@code 3} will produce a tree having
* \(8 \times 4^3 = 512\) triangular facets inserted. See the table below for other examples. The returned BSP
* tree also contains structural cuts to reduce the overall height of the tree.
*
* <table>
* <caption>Subdivisions to Triangle Counts</caption>
* <thead>
* <tr>
* <th>Subdivisions</th>
* <th>Triangles</th>
* </tr>
* </thead>
* <tbody>
* <tr><td>0</td><td>8</td></tr>
* <tr><td>1</td><td>32</td></tr>
* <tr><td>2</td><td>128</td></tr>
* <tr><td>3</td><td>512</td></tr>
* <tr><td>4</td><td>2048</td></tr>
* <tr><td>5</td><td>8192</td></tr>
* </tbody>
* </table>
*
* <p>Care must be taken when using this method with large subdivision numbers so that floating point errors
* do not interfere with the creation of the planes and triangles in the tree. For example, if the number of
* subdivisions is too high, the subdivided triangle points may become equivalent according to the sphere's
* {@link #getPrecision() precision context} and plane creation may fail. Or plane creation may succeed but
* insertion of the plane into the tree may fail for similar reasons. In general, it is best to use the lowest
* subdivision number practical for the intended purpose.</p>
* @param subdivisions the number of triangle subdivisions to use when creating the tree; the total number of
* triangular facets inserted into the returned tree is equal to \(8 \times 4^s\), where \(s\) is the number
* of subdivisions
* @return a BSP tree containing an approximation of the sphere
* @throws IllegalArgumentException if {@code subdivisions} is less than zero
* @throws IllegalStateException if tree creation fails for the given subdivision count
* @see #toTriangleMesh(int)
*/
public RegionBSPTree3D toTree(final int subdivisions) {
if (subdivisions < 0) {
throw new IllegalArgumentException(MessageFormat.format(INVALID_SUBDIVISION_MESSAGE, subdivisions));
}
return new SphereTreeApproximationBuilder(this, subdivisions).build();
}
/** Build an approximation of this sphere using a {@link TriangleMesh}. The approximation is constructed by
* taking an octahedron (8-sided polyhedron with triangular faces) inscribed in the sphere and subdividing each
* triangular face {@code subdivisions} number of times, each time projecting the newly created vertices onto the
* sphere surface. Each triangle subdivision produces 4 triangles, meaning that the total number of triangles
* in the returned mesh is equal to \(8 \times 4^s\), where \(s\) is the number of subdivisions. For
* example, calling this method with {@code subdivisions} equal to {@code 3} will produce a mesh having
* \(8 \times 4^3 = 512\) triangular facets inserted. See the table below for other examples.
*
* <table>
* <caption>Subdivisions to Triangle Counts</caption>
* <thead>
* <tr>
* <th>Subdivisions</th>
* <th>Triangles</th>
* </tr>
* </thead>
* <tbody>
* <tr><td>0</td><td>8</td></tr>
* <tr><td>1</td><td>32</td></tr>
* <tr><td>2</td><td>128</td></tr>
* <tr><td>3</td><td>512</td></tr>
* <tr><td>4</td><td>2048</td></tr>
* <tr><td>5</td><td>8192</td></tr>
* </tbody>
* </table>
*
* <p><strong>BSP Tree Conversion</strong></p>
* <p>Inserting the boundaries of a sphere mesh approximation directly into a BSP tree will invariably result
* in poor performance: since the region is convex the constructed BSP tree degenerates into a simple linked
* list of nodes. If a BSP tree is needed, users should prefer the {@link #toTree(int)} method, which creates
* balanced tree approximations directly, or the {@link RegionBSPTree3D.PartitionedRegionBuilder3D} class,
* which can be used to insert the mesh faces into a pre-partitioned tree.
* </p>
* @param subdivisions the number of triangle subdivisions to use when creating the mesh; the total number of
* triangular faces in the returned mesh is equal to \(8 \times 4^s\), where \(s\) is the number
* of subdivisions
* @return a triangle mesh approximation of the sphere
* @throws IllegalArgumentException if {@code subdivisions} is less than zero
* @see #toTree(int)
*/
public TriangleMesh toTriangleMesh(final int subdivisions) {
if (subdivisions < 0) {
throw new IllegalArgumentException(MessageFormat.format(INVALID_SUBDIVISION_MESSAGE, subdivisions));
}
return new SphereMeshApproximationBuilder(this, subdivisions).build();
}
/** Get the intersections of the given line with this sphere. The returned list will
* contain either 0, 1, or 2 points.
* <ul>
* <li><strong>2 points</strong> - The line is a secant line and intersects the sphere at two
* distinct points. The points are ordered such that the first point in the list is the first point
* encountered when traveling in the direction of the line. (In other words, the points are ordered
* by increasing abscissa value.)
* </li>
* <li><strong>1 point</strong> - The line is a tangent line and only intersects the sphere at a
* single point (as evaluated by the sphere's precision context).
* </li>
* <li><strong>0 points</strong> - The line does not intersect the sphere.</li>
* </ul>
* @param line line to intersect with the sphere
* @return a list of intersection points between the given line and this sphere
*/
public List<Vector3D> intersections(final Line3D line) {
return intersections(line, Line3D::abscissa, Line3D::distance);
}
/** Get the first intersection point between the given line and this sphere, or null
* if no such point exists. The "first" intersection point is the first such point
* encountered when traveling in the direction of the line from infinity.
* @param line line to intersect with the sphere
* @return the first intersection point between the given line and this instance or
* null if no such point exists
*/
public Vector3D firstIntersection(final Line3D line) {
return firstIntersection(line, Line3D::abscissa, Line3D::distance);
}
/** {@inheritDoc} */
@Override
public List<LinecastPoint3D> linecast(final LineConvexSubset3D subset) {
return getLinecastStream(subset)
.collect(Collectors.toList());
}
/** {@inheritDoc} */
@Override
public LinecastPoint3D linecastFirst(final LineConvexSubset3D subset) {
return getLinecastStream(subset)
.findFirst()
.orElse(null);
}
/** Get a stream containing the linecast intersection points of the given
* line subset with this instance.
* @param subset line subset to intersect against this instance
* @return a stream containing linecast intersection points
*/
private Stream<LinecastPoint3D> getLinecastStream(final LineConvexSubset3D subset) {
return intersections(subset.getLine()).stream()
.filter(subset::contains)
.map(pt -> new LinecastPoint3D(pt, getCenter().directionTo(pt), subset.getLine()));
}
/** Construct a sphere from a center point and radius.
* @param center the center of the sphere
* @param radius the sphere radius
* @param precision precision context used to compare floating point numbers
* @return a sphere constructed from the given center point and radius
* @throws IllegalArgumentException if center is not finite or radius is not finite or is
* less than or equal to zero as evaluated by the given precision context
*/
public static Sphere from(final Vector3D center, final double radius, final Precision.DoubleEquivalence precision) {
return new Sphere(center, radius, precision);
}
/** Internal class used to construct hyperplane-bounded approximations of spheres as BSP trees. The class
* begins with an octahedron inscribed in the sphere and then subdivides each triangular face a specified
* number of times.
*/
private static final class SphereTreeApproximationBuilder {
/** Threshold used to determine when to stop inserting structural cuts and begin adding facets. */
private static final int PARTITION_THRESHOLD = 2;
/** The sphere that an approximation is being created for. */
private final Sphere sphere;
/** The number of triangular subdivisions to use. */
private final int subdivisions;
/** Construct a new builder for creating a BSP tree approximation of the given sphere.
* @param sphere the sphere to create an approximation of
* @param subdivisions the number of triangle subdivisions to use in tree creation
*/
SphereTreeApproximationBuilder(final Sphere sphere, final int subdivisions) {
this.sphere = sphere;
this.subdivisions = subdivisions;
}
/** Build the sphere approximation BSP tree.
* @return the sphere approximation BSP tree
* @throws IllegalStateException if tree creation fails for the configured subdivision count
*/
RegionBSPTree3D build() {
final RegionBSPTree3D tree = RegionBSPTree3D.empty();
final Vector3D center = sphere.getCenter();
final double radius = sphere.getRadius();
final Precision.DoubleEquivalence precision = sphere.getPrecision();
// insert the primary split planes
final Plane plusXPlane = Planes.fromPointAndNormal(center, Vector3D.Unit.PLUS_X, precision);
final Plane plusYPlane = Planes.fromPointAndNormal(center, Vector3D.Unit.PLUS_Y, precision);
final Plane plusZPlane = Planes.fromPointAndNormal(center, Vector3D.Unit.PLUS_Z, precision);
tree.insert(plusXPlane.span(), RegionCutRule.INHERIT);
tree.insert(plusYPlane.span(), RegionCutRule.INHERIT);
tree.insert(plusZPlane.span(), RegionCutRule.INHERIT);
// create the vertices for the octahedron
final double cx = center.getX();
final double cy = center.getY();
final double cz = center.getZ();
final Vector3D maxX = Vector3D.of(cx + radius, cy, cz);
final Vector3D minX = Vector3D.of(cx - radius, cy, cz);
final Vector3D maxY = Vector3D.of(cx, cy + radius, cz);
final Vector3D minY = Vector3D.of(cx, cy - radius, cz);
final Vector3D maxZ = Vector3D.of(cx, cy, cz + radius);
final Vector3D minZ = Vector3D.of(cx, cy, cz - radius);
// partition and subdivide the face triangles and insert them into the tree
final RegionNode3D root = tree.getRoot();
try {
partitionAndInsert(root.getMinus().getMinus().getMinus(), minX, minZ, minY, 0);
partitionAndInsert(root.getMinus().getMinus().getPlus(), minX, minY, maxZ, 0);
partitionAndInsert(root.getMinus().getPlus().getMinus(), minX, maxY, minZ, 0);
partitionAndInsert(root.getMinus().getPlus().getPlus(), minX, maxZ, maxY, 0);
partitionAndInsert(root.getPlus().getMinus().getMinus(), maxX, minY, minZ, 0);
partitionAndInsert(root.getPlus().getMinus().getPlus(), maxX, maxZ, minY, 0);
partitionAndInsert(root.getPlus().getPlus().getMinus(), maxX, minZ, maxY, 0);
partitionAndInsert(root.getPlus().getPlus().getPlus(), maxX, maxY, maxZ, 0);
} catch (final IllegalStateException | IllegalArgumentException exc) {
// standardize any tree construction failure as an IllegalStateException
throw new IllegalStateException("Failed to construct sphere approximation with subdivision count " +
subdivisions + ": " + exc.getMessage(), exc);
}
return tree;
}
/** Recursively insert structural BSP tree cuts into the given node and then insert subdivided triangles
* when a target subdivision level is reached. The structural BSP tree cuts are used to help reduce the
* overall depth of the BSP tree.
* @param node the node to insert into
* @param p1 first triangle point
* @param p2 second triangle point
* @param p3 third triangle point
* @param level current subdivision level
*/
private void partitionAndInsert(final RegionNode3D node,
final Vector3D p1, final Vector3D p2, final Vector3D p3, final int level) {
if (subdivisions - level > PARTITION_THRESHOLD) {
final int nextLevel = level + 1;
final Vector3D center = sphere.getCenter();
final Vector3D m1 = sphere.project(p1.lerp(p2, 0.5));
final Vector3D m2 = sphere.project(p2.lerp(p3, 0.5));
final Vector3D m3 = sphere.project(p3.lerp(p1, 0.5));
RegionNode3D curNode = node;
checkedCut(curNode, createPlane(m3, m2, center), RegionCutRule.INHERIT);
partitionAndInsert(curNode.getPlus(), m3, m2, p3, nextLevel);
curNode = curNode.getMinus();
checkedCut(curNode, createPlane(m2, m1, center), RegionCutRule.INHERIT);
partitionAndInsert(curNode.getPlus(), m1, p2, m2, nextLevel);
curNode = curNode.getMinus();
checkedCut(curNode, createPlane(m1, m3, center), RegionCutRule.INHERIT);
partitionAndInsert(curNode.getPlus(), p1, m1, m3, nextLevel);
partitionAndInsert(curNode.getMinus(), m1, m2, m3, nextLevel);
} else {
insertSubdividedTriangles(node, p1, p2, p3, level);
}
}
/** Recursively insert subdivided triangles into the given node. Each triangle is inserted into the minus
* side of the previous triangle.
* @param node the node to insert into
* @param p1 first triangle point
* @param p2 second triangle point
* @param p3 third triangle point
* @param level the current subdivision level
* @return the node representing the inside of the region after insertion of all triangles
*/
private RegionNode3D insertSubdividedTriangles(final RegionNode3D node,
final Vector3D p1, final Vector3D p2, final Vector3D p3,
final int level) {
if (level >= subdivisions) {
// base case
checkedCut(node, createPlane(p1, p2, p3), RegionCutRule.MINUS_INSIDE);
return node.getMinus();
} else {
final int nextLevel = level + 1;
final Vector3D m1 = sphere.project(p1.lerp(p2, 0.5));
final Vector3D m2 = sphere.project(p2.lerp(p3, 0.5));
final Vector3D m3 = sphere.project(p3.lerp(p1, 0.5));
RegionNode3D curNode = node;
curNode = insertSubdividedTriangles(curNode, p1, m1, m3, nextLevel);
curNode = insertSubdividedTriangles(curNode, m1, p2, m2, nextLevel);
curNode = insertSubdividedTriangles(curNode, m3, m2, p3, nextLevel);
curNode = insertSubdividedTriangles(curNode, m1, m2, m3, nextLevel);
return curNode;
}
}
/** Create a plane from the given points, using the precision context of the sphere.
* @param p1 first point
* @param p2 second point
* @param p3 third point
* @return a plane defined by the given points
*/
private Plane createPlane(final Vector3D p1, final Vector3D p2, final Vector3D p3) {
return Planes.fromPoints(p1, p2, p3, sphere.getPrecision());
}
/** Insert the cut into the given node, throwing an exception if no portion of the cutter intersects
* the node.
* @param node node to cut
* @param cutter plane to use to cut the node
* @param cutRule cut rule to apply
* @throws IllegalStateException if no portion of the cutter plane intersects the node
*/
private void checkedCut(final RegionNode3D node, final Plane cutter, final RegionCutRule cutRule) {
if (!node.insertCut(cutter, cutRule)) {
throw new IllegalStateException("Failed to cut BSP tree node with plane: " + cutter);
}
}
}
/** Internal class used to construct geodesic mesh sphere approximations. The class begins with an octahedron
* inscribed in the sphere and then subdivides each triangular face a specified number of times.
*/
private static final class SphereMeshApproximationBuilder {
/** The sphere that an approximation is being created for. */
private final Sphere sphere;
/** The number of triangular subdivisions to use. */
private final int subdivisions;
/** Mesh builder object. */
private final SimpleTriangleMesh.Builder builder;
/** Construct a new builder for creating a mesh approximation of the given sphere.
* @param sphere the sphere to create an approximation of
* @param subdivisions the number of triangle subdivisions to use in mesh creation
*/
SphereMeshApproximationBuilder(final Sphere sphere, final int subdivisions) {
this.sphere = sphere;
this.subdivisions = subdivisions;
this.builder = SimpleTriangleMesh.builder(sphere.getPrecision());
}
/** Build the mesh approximation of the configured sphere.
* @return the mesh approximation of the configured sphere
*/
public SimpleTriangleMesh build() {
final Vector3D center = sphere.getCenter();
final double radius = sphere.getRadius();
// create the vertices for the octahedron
final double cx = center.getX();
final double cy = center.getY();
final double cz = center.getZ();
final Vector3D maxX = Vector3D.of(cx + radius, cy, cz);
final Vector3D minX = Vector3D.of(cx - radius, cy, cz);
final Vector3D maxY = Vector3D.of(cx, cy + radius, cz);
final Vector3D minY = Vector3D.of(cx, cy - radius, cz);
final Vector3D maxZ = Vector3D.of(cx, cy, cz + radius);
final Vector3D minZ = Vector3D.of(cx, cy, cz - radius);
addSubdivided(minX, minZ, minY, 0);
addSubdivided(minX, minY, maxZ, 0);
addSubdivided(minX, maxY, minZ, 0);
addSubdivided(minX, maxZ, maxY, 0);
addSubdivided(maxX, minY, minZ, 0);
addSubdivided(maxX, maxZ, minY, 0);
addSubdivided(maxX, minZ, maxY, 0);
addSubdivided(maxX, maxY, maxZ, 0);
return builder.build();
}
/** Recursively subdivide and add triangular faces between the given outer boundary points.
* @param p1 first point
* @param p2 second point
* @param p3 third point
* @param level recursion level; counts up
*/
private void addSubdivided(final Vector3D p1, final Vector3D p2, final Vector3D p3, final int level) {
if (level >= subdivisions) {
// base case
builder.addFaceUsingVertices(p1, p2, p3);
} else {
// subdivide
final int nextLevel = level + 1;
final Vector3D m1 = sphere.project(p1.lerp(p2, 0.5));
final Vector3D m2 = sphere.project(p2.lerp(p3, 0.5));
final Vector3D m3 = sphere.project(p3.lerp(p1, 0.5));
addSubdivided(p1, m1, m3, nextLevel);
addSubdivided(m1, p2, m2, nextLevel);
addSubdivided(m3, m2, p3, nextLevel);
addSubdivided(m1, m2, m3, nextLevel);
}
}
}
}