This method builds this intervals set as an ordered list of + * {@link Interval Interval} elements. If the intervals set has no + * lower limit, the first interval will have its low bound equal to + * {@code Double.NEGATIVE_INFINITY}. If the intervals set has + * no upper limit, the last interval will have its upper bound equal + * to {@code Double.POSITIVE_INFINITY}. An empty tree will + * build an empty list while a tree representing the whole real line + * will build a one element list with both bounds beeing + * infinite.
+ * @return a new ordered list containing {@link Interval Interval} + * elements + */ + public ListAn hyperplane in 1D is a simple point, its orientation being a + * boolean.
+ *Instances of this class are guaranteed to be immutable.
+ * @version $Revision$ $Date$ + */ +public class OrientedPoint implements Hyperplane { + + /** Dummy region returned by the {@link #wholeHyperplane} method. */ + private static final Region DUMMY_REGION = new DummyRegion(); + + /** Point location. */ + private Point1D location; + + /** Orientation. */ + private boolean direct; + + /** Simple constructor. + * @param location location of the hyperplane + * @param direct if true, the plus side of the hyperplane is towards + * abscissae greater than {@code location} + */ + public OrientedPoint(final Point1D location, final boolean direct) { + this.location = location; + this.direct = direct; + } + + /** Copy the instance. + *Since instances are immutable, this method directly returns + * the instance.
+ * @return the instance itself + */ + public Hyperplane copySelf() { + return this; + } + + /** Get the offset (oriented distance) of a point to the hyperplane. + * @param point point to check + * @return offset of the point + */ + public double getOffset(final Point point) { + final double delta = ((Point1D) point).getAbscissa() - location.getAbscissa(); + return direct ? delta : -delta; + } + + /** Transform a space point into a sub-space point. + *Since this class represent zero dimension spaces which does + * not have lower dimension sub-spaces, this method cannot be + * supported here. It always throws a {@code RuntimeException} + * when called.
+ * @param point n-dimension point of the space + * @return (n-1)-dimension point of the sub-space corresponding to + * the specified space point + * @see #toSpace + */ + public Point toSubSpace(final Point point) { + throw new MathUnsupportedOperationException(LocalizedFormats.NOT_SUPPORTED_IN_DIMENSION_N, 1); + } + + /** Transform a sub-space point into a space point. + *Since this class represent zero dimension spaces which does + * not have lower dimension sub-spaces, this method cannot be + * supported here. It always throws a {@code RuntimeException} + * when called.
+ * @param point (n-1)-dimension point of the sub-space + * @return n-dimension point of the space corresponding to the + * specified sub-space point + * @see #toSubSpace + */ + public Point toSpace(final Point point) { + throw new MathUnsupportedOperationException(LocalizedFormats.NOT_SUPPORTED_IN_DIMENSION_N, 1); + } + + /** Build the sub-space shared by the instance and another hyperplane. + *Since this class represent zero dimension spaces which does + * not have lower dimension sub-spaces, this method cannot be + * supported here. It always throws a {@code RuntimeException} + * when called.
+ * @param other other sub-space (must have the same dimension as the + * instance) + * @return a sub-space at the intersection of the instance and the + * other sub-space (it has a dimension one unit less than the + * instance) + */ + public SubSpace intersection(final Hyperplane other) { + throw new MathUnsupportedOperationException(LocalizedFormats.NOT_SUPPORTED_IN_DIMENSION_N, 1); + } + + /** Build a region covering the whole hyperplane. + *Since this class represent zero dimension spaces which does + * not have lower dimension sub-spaces, this method returns a dummy + * implementation of a {@link Region Region} (always the same + * instance). This implementation is only used to allow the {@link + * SubHyperplane SubHyperplane} class implementation to work + * properly, it should not be used otherwise.
+ * @return a dummy region + */ + public Region wholeHyperplane() { + return DUMMY_REGION; + } + + /** Build a region covering the whole space. + * @return a region containing the instance (really an {@link + * IntervalsSet IntervalsSet} instance) + */ + public Region wholeSpace() { + return new IntervalsSet(); + } + + /** Check if the instance has the same orientation as another hyperplane. + *This method is expected to be called on parallel hyperplanes + * (i.e. when the {@link #side side} method would return {@link + * org.apache.commons.math.geometry.partitioning.Hyperplane.Side#HYPER} + * for some sub-hyperplane having the specified hyperplane + * as its underlying hyperplane). The method should not + * re-check for parallelism, only for orientation, typically by + * testing something like the sign of the dot-products of + * normals.
+ * @param other other hyperplane to check against the instance + * @return true if the instance and the other hyperplane have + * the same orientation + */ + public boolean sameOrientationAs(final Hyperplane other) { + return !(direct ^ ((OrientedPoint) other).direct); + } + + /** Compute the relative position of a sub-hyperplane with respect + * to the instance. + * @param sub sub-hyperplane to check + * @return one of {@link org.apache.commons.math.geometry.partitioning.Hyperplane.Side#PLUS PLUS}, + * {@link org.apache.commons.math.geometry.partitioning.Hyperplane.Side#MINUS MINUS} + * or {@link org.apache.commons.math.geometry.partitioning.Hyperplane.Side#HYPER HYPER} + * (in dimension 1, this method never returns {@link + * org.apache.commons.math.geometry.partitioning.Hyperplane.Side#BOTH BOTH}) + * + */ + public Side side(final SubHyperplane sub) { + final double global = getOffset(((OrientedPoint) sub.getHyperplane()).location); + return (global < -1.0e-10) ? Side.MINUS : ((global > 1.0e-10) ? Side.PLUS : Side.HYPER); + } + + /** Split a sub-hyperplane in two parts by the instance. + * @param sub sub-hyperplane to split + * @return an object containing both the part of the sub-hyperplane + * on the plus side of the instance and the part of the + * sub-hyperplane on the minus side of the instance + */ + public SplitSubHyperplane split(final SubHyperplane sub) { + final double global = getOffset(((OrientedPoint) sub.getHyperplane()).location); + return (global < -1.0e-10) ? new SplitSubHyperplane(null, sub) : new SplitSubHyperplane(sub, null); + } + + /** Get the hyperplane location on the real line. + * @return the hyperplane location + */ + public Point1D getLocation() { + return location; + } + + /** Check if the hyperplane orientation is direct. + * @return true if the plus side of the hyperplane is towards + * abscissae greater than hyperplane location + */ + public boolean isDirect() { + return direct; + } + + /** Revert the instance. + */ + public void revertSelf() { + direct = !direct; + } + + /** Dummy region representing the whole set of reals. */ + private static class DummyRegion extends Region { + + /** Simple constructor. + */ + public DummyRegion() { + super(); + } + + /** {@inheritDoc} */ + public Region buildNew(final BSPTree tree) { + return this; + } + + /** {@inheritDoc} */ + protected void computeGeometricalProperties() { + setSize(0); + setBarycenter(Point1D.ZERO); + } + } + +} diff --git a/src/main/java/org/apache/commons/math/geometry/euclidean/oneD/Point1D.java b/src/main/java/org/apache/commons/math/geometry/euclidean/oneD/Point1D.java new file mode 100644 index 000000000..cb607ac6f --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/oneD/Point1D.java @@ -0,0 +1,53 @@ +/* + * 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.math.geometry.euclidean.oneD; + +import org.apache.commons.math.geometry.partitioning.Point; + +/** This class represents a 1D point. + *Instances of this class are guaranteed to be immutable.
+ * @version $Revision$ $Date$ + */ +public class Point1D implements Point { + + /** Point at 0.0 abscissa. */ + public static final Point1D ZERO = new Point1D(0.0); + + /** Point at 1.0 abscissa. */ + public static final Point1D ONE = new Point1D(1.0); + + /** Point at undefined (NaN) abscissa. */ + public static final Point1D UNDEFINED = new Point1D(Double.NaN); + + /** Abscissa of the point. */ + private double x; + + /** Simple constructor. + * @param x abscissa of the point + */ + public Point1D(final double x) { + this.x = x; + } + + /** Get the abscissa of the point. + * @return abscissa of the point + */ + public double getAbscissa() { + return x; + } + +} diff --git a/src/main/java/org/apache/commons/math/geometry/euclidean/oneD/package.html b/src/main/java/org/apache/commons/math/geometry/euclidean/oneD/package.html new file mode 100644 index 000000000..2ac354fe7 --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/oneD/package.html @@ -0,0 +1,24 @@ + + + + ++This package provides basic 1D geometry components. +
+ + diff --git a/src/main/java/org/apache/commons/math/geometry/CardanEulerSingularityException.java b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/CardanEulerSingularityException.java similarity index 96% rename from src/main/java/org/apache/commons/math/geometry/CardanEulerSingularityException.java rename to src/main/java/org/apache/commons/math/geometry/euclidean/threeD/CardanEulerSingularityException.java index 1125687ac..c6498fb65 100644 --- a/src/main/java/org/apache/commons/math/geometry/CardanEulerSingularityException.java +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/CardanEulerSingularityException.java @@ -15,7 +15,7 @@ * limitations under the License. */ -package org.apache.commons.math.geometry; +package org.apache.commons.math.geometry.euclidean.threeD; import org.apache.commons.math.MathException; import org.apache.commons.math.exception.util.LocalizedFormats; diff --git a/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/Line.java b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/Line.java new file mode 100644 index 000000000..0604e5286 --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/Line.java @@ -0,0 +1,153 @@ +/* + * 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.math.geometry.euclidean.threeD; + +import org.apache.commons.math.geometry.euclidean.oneD.Point1D; +import org.apache.commons.math.geometry.partitioning.Point; +import org.apache.commons.math.geometry.partitioning.SubSpace; +import org.apache.commons.math.util.FastMath; + +/** The class represent lines in a three dimensional space. + + *Each oriented line is intrinsically associated with an abscissa + * wich is a coordinate on the line. The point at abscissa 0 is the + * orthogonal projection of the origin on the line, another equivalent + * way to express this is to say that it is the point of the line + * which is closest to the origin. Abscissa increases in the line + * direction.
+ + * @version $Revision$ $Date$ + */ +public class Line implements SubSpace { + + /** Line direction. */ + private Vector3D direction; + + /** Line point closest to the origin. */ + private Point3D zero; + + /** Build a line from a point and a direction. + * @param p point belonging to the line (this can be any point) + * @param direction direction of the line + * @exception IllegalArgumentException if the direction norm is too small + */ + public Line(final Vector3D p, final Vector3D direction) { + reset(p, direction); + } + + /** Reset the instance as if built from a point and a normal. + * @param p point belonging to the line (this can be any point) + * @param dir direction of the line + * @exception IllegalArgumentException if the direction norm is too small + */ + public void reset(final Vector3D p, final Vector3D dir) { + final double norm = dir.getNorm(); + if (norm == 0.0) { + throw new IllegalArgumentException("null norm"); + } + this.direction = new Vector3D(1.0 / norm, dir); + zero = new Point3D(1.0, p, -Vector3D.dotProduct(p, this.direction), this.direction); + } + + /** Revert the line direction. + */ + public void revertSelf() { + direction = direction.negate(); + } + + /** Get the normalized direction vector. + * @return normalized direction vector + */ + public Vector3D getDirection() { + return direction; + } + + /** Get the abscissa of a point with respect to the line. + *The abscissa is 0 if the projection of the point and the + * projection of the frame origin on the line are the same + * point.
+ * @param point point to check (must be a {@link Vector3D Vector3D} + * instance) + * @return abscissa of the point (really a + * {org.apache.commons.math.geometry.euclidean.oneD.Point1D Point1D} instance) + */ + public Point toSubSpace(final Point point) { + final double x = Vector3D.dotProduct(((Vector3D) point).subtract(zero), direction); + return new Point1D(x); + } + + /** Get one point from the line. + * @param point desired abscissa for the point (must be a + * {org.apache.commons.math.geometry.euclidean.oneD.Point1D Point1D} instance) + * @return one point belonging to the line, at specified abscissa + * (really a {@link Vector3D Vector3D} instance) + */ + public Point toSpace(final Point point) { + return new Point3D(1.0, zero, ((Point1D) point).getAbscissa(), direction); + } + + /** Check if the instance is similar to another line. + *Lines are considered similar if they contain the same + * points. This does not mean they are equal since they can have + * opposite directions.
+ * @param line line to which instance should be compared + * @return true if the lines are similar + */ + public boolean isSimilarTo(final Line line) { + final double angle = Vector3D.angle(direction, line.direction); + return ((angle < 1.0e-10) || (angle > (FastMath.PI - 1.0e-10))) && contains(line.zero); + } + + /** Check if the instance contains a point. + * @param p point to check + * @return true if p belongs to the line + */ + public boolean contains(final Vector3D p) { + return distance(p) < 1.0e-10; + } + + /** Compute the distance between the instance and a point. + * @param p to check + * @return distance between the instance and the point + */ + public double distance(final Vector3D p) { + final Vector3D d = p.subtract(zero); + final Vector3D n = new Vector3D(1.0, d, -Vector3D.dotProduct(d, direction), direction); + return n.getNorm(); + } + + /** Compute the shortest distance between the instance and another line. + * @param line line to check agains the instance + * @return shortest distance between the instance and the line + */ + public double distance(final Line line) { + + final Vector3D normal = Vector3D.crossProduct(direction, line.direction); + if (normal.getNorm() < 1.0e-10) { + // lines are parallel + return distance(line.zero); + } + + // separating middle plane + final Plane middle = new Plane(new Vector3D(0.5, zero, 0.5, line.zero), normal); + + // the lines are at the same distance on either side of the plane + return 2 * FastMath.abs(middle.getOffset(zero)); + + } + +} diff --git a/src/main/java/org/apache/commons/math/geometry/NotARotationMatrixException.java b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/NotARotationMatrixException.java similarity index 96% rename from src/main/java/org/apache/commons/math/geometry/NotARotationMatrixException.java rename to src/main/java/org/apache/commons/math/geometry/euclidean/threeD/NotARotationMatrixException.java index 910e40157..2b5737b9f 100644 --- a/src/main/java/org/apache/commons/math/geometry/NotARotationMatrixException.java +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/NotARotationMatrixException.java @@ -15,7 +15,7 @@ * limitations under the License. */ -package org.apache.commons.math.geometry; +package org.apache.commons.math.geometry.euclidean.threeD; import org.apache.commons.math.MathException; import org.apache.commons.math.exception.util.Localizable; diff --git a/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/OutlineExtractor.java b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/OutlineExtractor.java new file mode 100644 index 000000000..842945297 --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/OutlineExtractor.java @@ -0,0 +1,250 @@ +/* + * 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.math.geometry.euclidean.threeD; + +import org.apache.commons.math.geometry.euclidean.twoD.Point2D; +import org.apache.commons.math.geometry.euclidean.twoD.PolygonsSet; +import org.apache.commons.math.geometry.partitioning.BSPTree; +import org.apache.commons.math.geometry.partitioning.BSPTreeVisitor; +import org.apache.commons.math.geometry.partitioning.Region; +import org.apache.commons.math.geometry.partitioning.SubHyperplane; +import org.apache.commons.math.util.FastMath; + +import java.util.ArrayList; + +/** Extractor for {@link PolygonsSet polyhedrons sets} outlines. + *This class extracts the 2D outlines from {{@link PolygonsSet + * polyhedrons sets} in a specified projection plane.
+ * @version $Revision$ $Date$ + */ +public class OutlineExtractor { + + /** Abscissa axis of the projection plane. */ + private Vector3D u; + + /** Ordinate axis of the projection plane. */ + private Vector3D v; + + /** Normal of the projection plane (viewing direction). */ + private Vector3D w; + + /** Build an extractor for a specific projection plane. + * @param u abscissa axis of the projection point + * @param v ordinate axis of the projection point + */ + public OutlineExtractor(final Vector3D u, final Vector3D v) { + this.u = u; + this.v = v; + w = Vector3D.crossProduct(u, v); + } + + /** Extract the outline of a polyhedrons set. + * @param polyhedronsSet polyhedrons set whose outline must be extracted + * @return an outline, as an array of loops. + */ + public Point2D[][] getOutline(final PolyhedronsSet polyhedronsSet) { + + // project all boundary facets into one polygons set + final BoundaryProjector projector = new BoundaryProjector(); + polyhedronsSet.getTree(true).visit(projector); + final PolygonsSet projected = projector.getProjected(); + + // Remove the spurious intermediate vertices from the outline + final Point2D[][] outline = projected.getVertices(); + for (int i = 0; i < outline.length; ++i) { + final Point2D[] rawLoop = outline[i]; + int end = rawLoop.length; + int j = 0; + while (j < end) { + if (pointIsBetween(rawLoop, end, j)) { + // the point should be removed + for (int k = j; k < (end - 1); ++k) { + rawLoop[k] = rawLoop[k + 1]; + } + --end; + } else { + // the point remains in the loop + ++j; + } + } + if (end != rawLoop.length) { + // resize the array + outline[i] = new Point2D[end]; + System.arraycopy(rawLoop, 0, outline[i], 0, end); + } + } + + return outline; + + } + + /** Check if a point is geometrically between its neighbour in an array. + *The neighbours are computed considering the array is a loop + * (i.e. point at index (n-1) is before point at index 0)
+ * @param loop points array + * @param n number of points to consider in the array + * @param i index of the point to check (must be between 0 and n-1) + * @return true if the point is exactly between its neighbours + */ + private boolean pointIsBetween(final Point2D[] loop, final int n, final int i) { + final Point2D previous = loop[(i + n - 1) % n]; + final Point2D current = loop[i]; + final Point2D next = loop[(i + 1) % n]; + final double dx1 = current.x - previous.x; + final double dy1 = current.y - previous.y; + final double dx2 = next.x - current.x; + final double dy2 = next.y - current.y; + final double cross = dx1 * dy2 - dx2 * dy1; + final double dot = dx1 * dx2 + dy1 * dy2; + final double d1d2 = FastMath.sqrt((dx1 * dx1 + dy1 * dy1) * (dx2 * dx2 + dy2 * dy2)); + return (FastMath.abs(cross) <= (1.0e-6 * d1d2)) && (dot >= 0.0); + } + + /** Visitor projecting the boundary facets on a plane. */ + private class BoundaryProjector implements BSPTreeVisitor { + + /** Projection of the polyhedrons set on the plane. */ + private PolygonsSet projected; + + /** Simple constructor. + */ + public BoundaryProjector() { + projected = new PolygonsSet(new BSPTree(Boolean.FALSE)); + } + + /** {@inheritDoc} */ + public Order visitOrder(final BSPTree node) { + return Order.MINUS_SUB_PLUS; + } + + /** {@inheritDoc} */ + public void visitInternalNode(final BSPTree node) { + final Region.BoundaryAttribute attribute = + (Region.BoundaryAttribute) node.getAttribute(); + if (attribute.getPlusOutside() != null) { + addContribution(attribute.getPlusOutside(), false); + } + if (attribute.getPlusInside() != null) { + addContribution(attribute.getPlusInside(), true); + } + } + + /** {@inheritDoc} */ + public void visitLeafNode(final BSPTree node) { + } + + /** Add he contribution of a boundary facet. + * @param facet boundary facet + * @param reversed if true, the facet has the inside on its plus side + */ + private void addContribution(final SubHyperplane facet, final boolean reversed) { + + // extract the vertices of the facet + final Plane plane = (Plane) facet.getHyperplane(); + Point2D[][] vertices = + ((PolygonsSet) facet.getRemainingRegion()).getVertices(); + + final double scal = Vector3D.dotProduct(plane.getNormal(), w); + if (FastMath.abs(scal) > 1.0e-3) { + + if ((scal < 0) ^ reversed) { + // the facet is seen from the inside, + // we need to invert its boundary orientation + final Point2D[][] newVertices = new Point2D[vertices.length][]; + for (int i = 0; i < vertices.length; ++i) { + final Point2D[] loop = vertices[i]; + final Point2D[] newLoop = new Point2D[loop.length]; + if (loop[0] == null) { + newLoop[0] = null; + for (int j = 1; j < loop.length; ++j) { + newLoop[j] = loop[loop.length - j]; + } + } else { + for (int j = 0; j < loop.length; ++j) { + newLoop[j] = loop[loop.length - (j + 1)]; + } + } + newVertices[i] = newLoop; + } + + // use the reverted vertices + vertices = newVertices; + + } + + // compute the projection of the facet in the outline plane + final ArrayListThe plane is oriented in the direction of + * {@code (p2-p1) ^ (p3-p1)}
+ * @param p1 first point belonging to the plane + * @param p2 second point belonging to the plane + * @param p3 third point belonging to the plane + * @exception IllegalArgumentException if the points do not constitute a plane + */ + public Plane(final Vector3D p1, final Vector3D p2, final Vector3D p3) { + this(p1, Vector3D.crossProduct(p2.subtract(p1), p3.subtract(p1))); + } + + /** Copy constructor. + *The instance created is completely independant of the original + * one. A deep copy is used, none of the underlying object are + * shared.
+ * @param plane plane to copy + */ + public Plane(final Plane plane) { + originOffset = plane.originOffset; + origin = plane.origin; + u = plane.u; + v = plane.v; + w = plane.w; + } + + /** Copy the instance. + *The instance created is completely independant of the original + * one. A deep copy is used, none of the underlying objects are + * shared (except for immutable objects).
+ * @return a new hyperplane, copy of the instance + */ + public Hyperplane copySelf() { + return new Plane(this); + } + + /** Reset the instance as if built from a point and a normal. + * @param p point belonging to the plane + * @param normal normal direction to the plane + */ + public void reset(final Vector3D p, final Vector3D normal) { + setNormal(normal); + originOffset = -Vector3D.dotProduct(p, w); + setFrame(); + } + + /** Reset the instance from another one. + *The updated instance is completely independant of the original + * one. A deep reset is used none of the underlying object is + * shared.
+ * @param original plane to reset from + */ + public void reset(final Plane original) { + originOffset = original.originOffset; + origin = original.origin; + u = original.u; + v = original.v; + w = original.w; + } + + /** Set the normal vactor. + * @param normal normal direction to the plane (will be copied) + * @exception IllegalArgumentException if the normal norm is too small + */ + private void setNormal(final Vector3D normal) { + final double norm = normal.getNorm(); + if (norm < 1.0e-10) { + throw new IllegalArgumentException("null norm"); + } + w = new Vector3D(1.0 / norm, normal); + } + + /** Reset the plane frame. + */ + private void setFrame() { + origin = new Point3D(-originOffset, w); + u = w.orthogonal(); + v = Vector3D.crossProduct(w, u); + } + + /** Get the origin point of the plane frame. + *The point returned is the orthogonal projection of the + * 3D-space origin in the plane.
+ * @return the origin point of the plane frame (point closest to the + * 3D-space origin) + */ + public Point3D getOrigin() { + return origin; + } + + /** Get the normalized normal vector. + *The frame defined by ({@link #getU getU}, {@link #getV getV}, + * {@link #getNormal getNormal}) is a rigth-handed orthonormalized + * frame).
+ * @return normalized normal vector + * @see #getU + * @see #getV + */ + public Vector3D getNormal() { + return w; + } + + /** Get the plane first canonical vector. + *The frame defined by ({@link #getU getU}, {@link #getV getV}, + * {@link #getNormal getNormal}) is a rigth-handed orthonormalized + * frame).
+ * @return normalized first canonical vector + * @see #getV + * @see #getNormal + */ + public Vector3D getU() { + return u; + } + + /** Get the plane second canonical vector. + *The frame defined by ({@link #getU getU}, {@link #getV getV}, + * {@link #getNormal getNormal}) is a rigth-handed orthonormalized + * frame).
+ * @return normalized second canonical vector + * @see #getU + * @see #getNormal + */ + public Vector3D getV() { + return v; + } + + /** Revert the plane. + *Replace the instance by a similar plane with opposite orientation.
+ *The new plane frame is chosen in such a way that a 3D point that had + * {@code (x, y)} in-plane coordinates and {@code z} offset with + * respect to the plane and is unaffected by the change will have + * {@code (y, x)} in-plane coordinates and {@code -z} offset with + * respect to the new plane. This means that the {@code u} and {@code v} + * vectors returned by the {@link #getU} and {@link #getV} methods are exchanged, + * and the {@code w} vector returned by the {@link #getNormal} method is + * reversed.
+ */ + public void revertSelf() { + final Vector3D tmp = u; + u = v; + v = tmp; + w = w.negate(); + originOffset = -originOffset; + } + + /** Transform a 3D space point into an in-plane point. + * @param point point of the space (must be a {@link Vector3D + * Vector3D} instance) + * @return in-plane point (really a {@link + * org.apache.commons.math.geometry.euclidean.twoD.Point2D Point2D} instance) + * @see #toSpace + */ + public Point toSubSpace(final Point point) { + final Vector3D p3D = (Vector3D) point; + return new Point2D(Vector3D.dotProduct(p3D, u), + Vector3D.dotProduct(p3D, v)); + } + + /** Transform an in-plane point into a 3D space point. + * @param point in-plane point (must be a {@link + * org.apache.commons.math.geometry.euclidean.twoD.Point2D Point2D} instance) + * @return 3D space point (really a {@link Vector3D Vector3D} instance) + * @see #toSubSpace + */ + public Point toSpace(final Point point) { + final Point2D p2D = (Point2D) point; + return new Point3D(p2D.x, u, p2D.y, v, -originOffset, w); + } + + /** Get one point from the 3D-space. + * @param inPlane desired in-plane coordinates for the point in the + * plane + * @param offset desired offset for the point + * @return one point in the 3D-space, with given coordinates and offset + * relative to the plane + */ + public Vector3D getPointAt(final Point2D inPlane, final double offset) { + return new Vector3D(inPlane.x, u, inPlane.y, v, offset - originOffset, w); + } + + /** Check if the instance is similar to another plane. + *Planes are considered similar if they contain the same + * points. This does not mean they are equal since they can have + * opposite normals.
+ * @param plane plane to which the instance is compared + * @return true if the planes are similar + */ + public boolean isSimilarTo(final Plane plane) { + final double angle = Vector3D.angle(w, plane.w); + return ((angle < 1.0e-10) && (FastMath.abs(originOffset - plane.originOffset) < 1.0e-10)) || + ((angle > (FastMath.PI - 1.0e-10)) && (FastMath.abs(originOffset + plane.originOffset) < 1.0e-10)); + } + + /** Rotate the plane around the specified point. + *The instance is not modified, a new instance is created.
+ * @param center rotation center + * @param rotation vectorial rotation operator + * @return a new plane + */ + public Plane rotate(final Vector3D center, final Rotation rotation) { + + final Vector3D delta = origin.subtract(center); + final Plane plane = new Plane(center.add(rotation.applyTo(delta)), + rotation.applyTo(w)); + + // make sure the frame is transformed as desired + plane.u = rotation.applyTo(u); + plane.v = rotation.applyTo(v); + + return plane; + + } + + /** Translate the plane by the specified amount. + *The instance is not modified, a new instance is created.
+ * @param translation translation to apply + * @return a new plane + */ + public Plane translate(final Vector3D translation) { + + final Plane plane = new Plane(origin.add(translation), w); + + // make sure the frame is transformed as desired + plane.u = u; + plane.v = v; + + return plane; + + } + + /** Get the intersection of a line with the instance. + * @param line line intersecting the instance + * @return intersection point between between the line and the + * instance (null if the line is parallel to the instance) + */ + public Point3D intersection(final Line line) { + final Vector3D direction = line.getDirection(); + final double dot = Vector3D.dotProduct(w, direction); + if (FastMath.abs(dot) < 1.0e-10) { + return null; + } + final Vector3D point = (Vector3D) line.toSpace(Point1D.ZERO); + final double k = -(originOffset + Vector3D.dotProduct(w, point)) / dot; + return new Point3D(1.0, point, k, direction); + } + + /** Build the line shared by the instance and another plane. + * @param other other plane + * @return line at the intersection of the instance and the + * other plane (really a {@link Line Line} instance) + */ + public SubSpace intersection(final Hyperplane other) { + final Plane otherP = (Plane) other; + final Vector3D direction = Vector3D.crossProduct(w, otherP.w); + if (direction.getNorm() < 1.0e-10) { + return null; + } + return new Line(intersection(this, otherP, new Plane(direction)), + direction); + } + + /** Get the intersection point of three planes. + * @param plane1 first plane1 + * @param plane2 second plane2 + * @param plane3 third plane2 + * @return intersection point of three planes, null if some planes are parallel + */ + public static Vector3D intersection(final Plane plane1, final Plane plane2, final Plane plane3) { + + // coefficients of the three planes linear equations + final double a1 = plane1.w.getX(); + final double b1 = plane1.w.getY(); + final double c1 = plane1.w.getZ(); + final double d1 = plane1.originOffset; + + final double a2 = plane2.w.getX(); + final double b2 = plane2.w.getY(); + final double c2 = plane2.w.getZ(); + final double d2 = plane2.originOffset; + + final double a3 = plane3.w.getX(); + final double b3 = plane3.w.getY(); + final double c3 = plane3.w.getZ(); + final double d3 = plane3.originOffset; + + // direct Cramer resolution of the linear system + // (this is still feasible for a 3x3 system) + final double a23 = b2 * c3 - b3 * c2; + final double b23 = c2 * a3 - c3 * a2; + final double c23 = a2 * b3 - a3 * b2; + final double determinant = a1 * a23 + b1 * b23 + c1 * c23; + if (FastMath.abs(determinant) < 1.0e-10) { + return null; + } + + final double r = 1.0 / determinant; + return new Vector3D( + (-a23 * d1 - (c1 * b3 - c3 * b1) * d2 - (c2 * b1 - c1 * b2) * d3) * r, + (-b23 * d1 - (c3 * a1 - c1 * a3) * d2 - (c1 * a2 - c2 * a1) * d3) * r, + (-c23 * d1 - (b1 * a3 - b3 * a1) * d2 - (b2 * a1 - b1 * a2) * d3) * r); + + } + + /** Build a region covering the whole hyperplane. + * @return a region covering the whole hyperplane + */ + public Region wholeHyperplane() { + return new PolygonsSet(); + } + + /** Build a region covering the whole space. + * @return a region containing the instance (really a {@link + * PolyhedronsSet PolyhedronsSet} instance) + */ + public Region wholeSpace() { + return new PolyhedronsSet(); + } + + /** Check if the instance contains a point. + * @param p point to check + * @return true if p belongs to the plane + */ + public boolean contains(final Point3D p) { + return FastMath.abs(getOffset(p)) < 1.0e-10; + } + + /** Get the offset (oriented distance) of a parallel plane. + *This method should be called only for parallel planes otherwise + * the result is not meaningful.
+ *The offset is 0 if both planes are the same, it is + * positive if the plane is on the plus side of the instance and + * negative if it is on the minus side, according to its natural + * orientation.
+ * @param plane plane to check + * @return offset of the plane + */ + public double getOffset(final Plane plane) { + return originOffset + (sameOrientationAs(plane) ? -plane.originOffset : plane.originOffset); + } + + /** Get the offset (oriented distance) of a point. + *The offset is 0 if the point is on the underlying hyperplane, + * it is positive if the point is on one particular side of the + * hyperplane, and it is negative if the point is on the other side, + * according to the hyperplane natural orientation.
+ * @param point point to check + * @return offset of the point + */ + public double getOffset(final Point point) { + return Vector3D.dotProduct((Vector3D) point, w) + originOffset; + } + + /** Check if the instance has the same orientation as another hyperplane. + * @param other other hyperplane to check against the instance + * @return true if the instance and the other hyperplane have + * the same orientation + */ + public boolean sameOrientationAs(final Hyperplane other) { + return Vector3D.dotProduct(((Plane) other).w, w) > 0.0; + } + + /** Compute the relative position of a sub-hyperplane with respect + * to the instance. + * @param sub sub-hyperplane to check + * @return one of {@link org.apache.commons.math.geometry.partitioning.Hyperplane.Side#PLUS PLUS}, + * {@link org.apache.commons.math.geometry.partitioning.Hyperplane.Side#MINUS MINUS}, + * {@link org.apache.commons.math.geometry.partitioning.Hyperplane.Side#BOTH BOTH}, + * {@link org.apache.commons.math.geometry.partitioning.Hyperplane.Side#HYPER HYPER} + */ + public Side side(final SubHyperplane sub) { + + final Plane otherPlane = (Plane) sub.getHyperplane(); + final Line inter = (Line) intersection(otherPlane); + + if (inter == null) { + // the hyperplanes are parallel, + // any point can be used to check their relative position + final double global = getOffset(otherPlane); + return (global < -1.0e-10) ? Side.MINUS : ((global > 1.0e-10) ? Side.PLUS : Side.HYPER); + } + + // create a 2D line in the otherPlane canonical 2D frame such that: + // - the line is the crossing line of the two planes in 3D + // - the line splits the otherPlane in two half planes with an + // orientation consistent with the orientation of the instance + // (i.e. the 3D half space on the plus side (resp. minus side) + // of the instance contains the 2D half plane on the plus side + // (resp. minus side) of the 2D line + Point2D p = (Point2D) otherPlane.toSubSpace(inter.toSpace(Point1D.ZERO)); + Point2D q = (Point2D) otherPlane.toSubSpace(inter.toSpace(Point1D.ONE)); + if (Vector3D.dotProduct(Vector3D.crossProduct(inter.getDirection(), + otherPlane.getNormal()), + w) < 0) { + final Point2D tmp = p; + p = q; + q = tmp; + } + final Hyperplane line2D = new org.apache.commons.math.geometry.euclidean.twoD.Line(p, q); + + // check the side on the 2D plane + return sub.getRemainingRegion().side(line2D); + + } + + /** Split a sub-hyperplane in two parts by the instance. + * @param sub sub-hyperplane to split + * @return an object containing both the part of the sub-hyperplane + * on the plus side of the instance and the part of the + * sub-hyperplane on the minus side of the instance + */ + public SplitSubHyperplane split(final SubHyperplane sub) { + + final Plane otherPlane = (Plane) sub.getHyperplane(); + final Line inter = (Line) intersection(otherPlane); + + if (inter == null) { + // the hyperplanes are parallel + final double global = getOffset(otherPlane); + return (global < -1.0e-10) ? new SplitSubHyperplane(null, sub) : new SplitSubHyperplane(sub, null); + } + + // the hyperplanes do intersect + Point2D p = (Point2D) otherPlane.toSubSpace(inter.toSpace(Point1D.ZERO)); + Point2D q = (Point2D) otherPlane.toSubSpace(inter.toSpace(Point1D.ONE)); + if (Vector3D.dotProduct(Vector3D.crossProduct(inter.getDirection(), + otherPlane.getNormal()), + w) < 0) { + final Point2D tmp = p; + p = q; + q = tmp; + } + final SubHyperplane l2DMinus = + new SubHyperplane(new org.apache.commons.math.geometry.euclidean.twoD.Line(p, q)); + final SubHyperplane l2DPlus = + new SubHyperplane(new org.apache.commons.math.geometry.euclidean.twoD.Line(q, p)); + + final BSPTree splitTree = + sub.getRemainingRegion().getTree(false).split(l2DMinus); + final BSPTree plusTree = Region.isEmpty(splitTree.getPlus()) ? + new BSPTree(Boolean.FALSE) : + new BSPTree(l2DPlus, new BSPTree(Boolean.FALSE), + splitTree.getPlus(), null); + + final BSPTree minusTree = Region.isEmpty(splitTree.getMinus()) ? + new BSPTree(Boolean.FALSE) : + new BSPTree(l2DMinus, new BSPTree(Boolean.FALSE), + splitTree.getMinus(), null); + + return new SplitSubHyperplane(new SubHyperplane(otherPlane.copySelf(), + new PolygonsSet(plusTree)), + new SubHyperplane(otherPlane.copySelf(), + new PolygonsSet(minusTree))); + + } + +} diff --git a/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/Point3D.java b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/Point3D.java new file mode 100644 index 000000000..127ae8d9d --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/Point3D.java @@ -0,0 +1,112 @@ +/* + * 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.math.geometry.euclidean.threeD; + +import org.apache.commons.math.geometry.partitioning.Point; + +/** This class represents a 3D point. + *Instances of this class are guaranteed to be immutable.
+ * @version $Revision$ $Date$ + */ +public class Point3D extends Vector3D implements Point { + + /** Point at undefined (NaN) coordinates. */ + public static final Point3D UNDEFINED = new Point3D(Double.NaN, Double.NaN, Double.NaN); + + /** Serializable UID. */ + private static final long serialVersionUID = 9128130934224884451L; + + /** Simple constructor. + * Build a vector from its coordinates + * @param x abscissa + * @param y ordinate + * @param z height + * @see #getX() + * @see #getY() + * @see #getZ() + */ + public Point3D(final double x, final double y, final double z) { + super(x, y, z); + } + + /** Simple constructor. + * Build a vector from its azimuthal coordinates + * @param alpha azimuth (α) around Z + * (0 is +X, π/2 is +Y, π is -X and 3π/2 is -Y) + * @param delta elevation (δ) above (XY) plane, from -π/2 to +π/2 + * @see #getAlpha() + * @see #getDelta() + */ + public Point3D(final double alpha, final double delta) { + super(alpha, delta); + } + + /** Multiplicative constructor + * Build a vector from another one and a scale factor. + * The vector built will be a * u + * @param a scale factor + * @param u base (unscaled) vector + */ + public Point3D(final double a, final Vector3D u) { + super(a, u); + } + + /** Linear constructor + * Build a vector from two other ones and corresponding scale factors. + * The vector built will be a1 * u1 + a2 * u2 + * @param a1 first scale factor + * @param u1 first base (unscaled) vector + * @param a2 second scale factor + * @param u2 second base (unscaled) vector + */ + public Point3D(final double a1, final Vector3D u1, final double a2, final Vector3D u2) { + super(a1, u1, a2, u2); + } + + /** Linear constructor + * Build a vector from three other ones and corresponding scale factors. + * The vector built will be a1 * u1 + a2 * u2 + a3 * u3 + * @param a1 first scale factor + * @param u1 first base (unscaled) vector + * @param a2 second scale factor + * @param u2 second base (unscaled) vector + * @param a3 third scale factor + * @param u3 third base (unscaled) vector + */ + public Point3D(final double a1, final Vector3D u1, final double a2, final Vector3D u2, + final double a3, final Vector3D u3) { + super(a1, u1, a2, u2, a3, u3); + } + + /** Linear constructor + * Build a vector from four other ones and corresponding scale factors. + * The vector built will be a1 * u1 + a2 * u2 + a3 * u3 + a4 * u4 + * @param a1 first scale factor + * @param u1 first base (unscaled) vector + * @param a2 second scale factor + * @param u2 second base (unscaled) vector + * @param a3 third scale factor + * @param u3 third base (unscaled) vector + * @param a4 fourth scale factor + * @param u4 fourth base (unscaled) vector + */ + public Point3D(final double a1, final Vector3D u1, final double a2, final Vector3D u2, + final double a3, final Vector3D u3, final double a4, final Vector3D u4) { + super(a1, u1, a2, u2, a3, u3, a4, u4); + } + +} diff --git a/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/PolyhedronsSet.java b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/PolyhedronsSet.java new file mode 100644 index 000000000..a3629397a --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/PolyhedronsSet.java @@ -0,0 +1,413 @@ +/* + * 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.math.geometry.euclidean.threeD; + +import java.awt.geom.AffineTransform; +import java.util.Arrays; +import java.util.Collection; + +import org.apache.commons.math.geometry.euclidean.twoD.Point2D; +import org.apache.commons.math.geometry.partitioning.BSPTree; +import org.apache.commons.math.geometry.partitioning.BSPTreeVisitor; +import org.apache.commons.math.geometry.partitioning.Hyperplane; +import org.apache.commons.math.geometry.partitioning.Point; +import org.apache.commons.math.geometry.partitioning.Region; +import org.apache.commons.math.geometry.partitioning.SubHyperplane; +import org.apache.commons.math.geometry.partitioning.Transform; +import org.apache.commons.math.util.FastMath; + +/** This class represents a 3D region: a set of polyhedrons. + * @version $Revision$ $Date$ + */ +public class PolyhedronsSet extends Region { + + /** Build a polyhedrons set representing the whole real line. + */ + public PolyhedronsSet() { + super(); + } + + /** Build a polyhedrons set from a BSP tree. + *The leaf nodes of the BSP tree must have a + * {@code Boolean} attribute representing the inside status of + * the corresponding cell (true for inside cells, false for outside + * cells). In order to avoid building too many small objects, it is + * recommended to use the predefined constants + * {@code Boolean.TRUE} and {@code Boolean.FALSE}
+ * @param tree inside/outside BSP tree representing the region + */ + public PolyhedronsSet(final BSPTree tree) { + super(tree); + } + + /** Build a polyhedrons set from a Boundary REPresentation (B-rep). + *The boundary is provided as a collection of {@link + * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the + * interior part of the region on its minus side and the exterior on + * its plus side.
+ *The boundary elements can be in any order, and can form + * several non-connected sets (like for example polyhedrons with holes + * or a set of disjoints polyhedrons considered as a whole). In + * fact, the elements do not even need to be connected together + * (their topological connections are not used here). However, if the + * boundary does not really separate an inside open from an outside + * open (open having here its topological meaning), then subsequent + * calls to the {@link Region#checkPoint(Point) checkPoint} method will + * not be meaningful anymore.
+ *If the boundary is empty, the region will represent the whole + * space.
+ * @param boundary collection of boundary elements, as a + * collection of {@link SubHyperplane SubHyperplane} objects + */ + public PolyhedronsSet(final CollectionThe instance is not modified, a new instance is created.
+ * @param center rotation center + * @param rotation vectorial rotation operator + * @return a new instance representing the rotated region + */ + public PolyhedronsSet rotate(final Vector3D center, final Rotation rotation) { + return (PolyhedronsSet) applyTransform(new RotationTransform(center, rotation)); + } + + /** 3D rotation as a Transform. */ + private static class RotationTransform implements Transform { + + /** Center point of the rotation. */ + private Vector3D center; + + /** Vectorial rotation. */ + private Rotation rotation; + + /** Cached original hyperplane. */ + private Hyperplane cachedOriginal; + + /** Cached 2D transform valid inside the cached original hyperplane. */ + private Transform cachedTransform; + + /** Build a rotation transform. + * @param center center point of the rotation + * @param rotation vectorial rotation + */ + public RotationTransform(final Vector3D center, final Rotation rotation) { + this.center = center; + this.rotation = rotation; + } + + /** {@inheritDoc} */ + public Point apply(final Point point) { + final Vector3D delta = ((Vector3D) point).subtract(center); + return new Point3D(1.0, center, 1.0, rotation.applyTo(delta)); + } + + /** {@inheritDoc} */ + public Hyperplane apply(final Hyperplane hyperplane) { + return ((Plane) hyperplane).rotate(center, rotation); + } + + /** {@inheritDoc} */ + public SubHyperplane apply(final SubHyperplane sub, + final Hyperplane original, final Hyperplane transformed) { + if (original != cachedOriginal) { + // we have changed hyperplane, reset the in-hyperplane transform + + final Plane oPlane = (Plane) original; + final Plane tPlane = (Plane) transformed; + final Vector3D p00 = oPlane.getOrigin(); + final Vector3D p10 = (Vector3D) oPlane.toSpace(new Point2D(1.0, 0.0)); + final Vector3D p01 = (Vector3D) oPlane.toSpace(new Point2D(0.0, 1.0)); + final Point2D tP00 = (Point2D) tPlane.toSubSpace(apply((Point) p00)); + final Point2D tP10 = (Point2D) tPlane.toSubSpace(apply((Point) p10)); + final Point2D tP01 = (Point2D) tPlane.toSubSpace(apply((Point) p01)); + final AffineTransform at = + new AffineTransform(tP10.getX() - tP00.getX(), tP10.getY() - tP00.getY(), + tP01.getX() - tP00.getX(), tP01.getY() - tP00.getY(), + tP00.getX(), tP00.getY()); + + cachedOriginal = original; + cachedTransform = org.apache.commons.math.geometry.euclidean.twoD.Line.getTransform(at); + + } + return sub.applyTransform(cachedTransform); + } + + } + + /** Translate the region by the specified amount. + *The instance is not modified, a new instance is created.
+ * @param translation translation to apply + * @return a new instance representing the translated region + */ + public PolyhedronsSet translate(final Vector3D translation) { + return (PolyhedronsSet) applyTransform(new TranslationTransform(translation)); + } + + /** 3D translation as a transform. */ + private static class TranslationTransform implements Transform { + + /** Translation vector. */ + private Vector3D translation; + + /** Cached original hyperplane. */ + private Hyperplane cachedOriginal; + + /** Cached 2D transform valid inside the cached original hyperplane. */ + private Transform cachedTransform; + + /** Build a translation transform. + * @param translation translation vector + */ + public TranslationTransform(final Vector3D translation) { + this.translation = translation; + } + + /** {@inheritDoc} */ + public Point apply(final Point point) { + return new Point3D(1.0, (Vector3D) point, 1.0, translation); + } + + /** {@inheritDoc} */ + public Hyperplane apply(final Hyperplane hyperplane) { + return ((Plane) hyperplane).translate(translation); + } + + /** {@inheritDoc} */ + public SubHyperplane apply(final SubHyperplane sub, + final Hyperplane original, final Hyperplane transformed) { + if (original != cachedOriginal) { + // we have changed hyperplane, reset the in-hyperplane transform + + final Plane oPlane = (Plane) original; + final Plane tPlane = (Plane) transformed; + final Point2D shift = (Point2D) tPlane.toSubSpace(apply((Point) oPlane.getOrigin())); + final AffineTransform at = + AffineTransform.getTranslateInstance(shift.getX(), shift.getY()); + + cachedOriginal = original; + cachedTransform = + org.apache.commons.math.geometry.euclidean.twoD.Line.getTransform(at); + + } + + return sub.applyTransform(cachedTransform); + + } + + } + +} diff --git a/src/main/java/org/apache/commons/math/geometry/Rotation.java b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/Rotation.java similarity index 99% rename from src/main/java/org/apache/commons/math/geometry/Rotation.java rename to src/main/java/org/apache/commons/math/geometry/euclidean/threeD/Rotation.java index babc4a2a6..7f4f443c9 100644 --- a/src/main/java/org/apache/commons/math/geometry/Rotation.java +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/Rotation.java @@ -15,7 +15,7 @@ * limitations under the License. */ -package org.apache.commons.math.geometry; +package org.apache.commons.math.geometry.euclidean.threeD; import java.io.Serializable; diff --git a/src/main/java/org/apache/commons/math/geometry/RotationOrder.java b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/RotationOrder.java similarity index 98% rename from src/main/java/org/apache/commons/math/geometry/RotationOrder.java rename to src/main/java/org/apache/commons/math/geometry/euclidean/threeD/RotationOrder.java index f6aae1972..844ddefd5 100644 --- a/src/main/java/org/apache/commons/math/geometry/RotationOrder.java +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/RotationOrder.java @@ -15,7 +15,7 @@ * limitations under the License. */ -package org.apache.commons.math.geometry; +package org.apache.commons.math.geometry.euclidean.threeD; /** * This class is a utility representing a rotation order specification diff --git a/src/main/java/org/apache/commons/math/geometry/Vector3D.java b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/Vector3D.java similarity index 99% rename from src/main/java/org/apache/commons/math/geometry/Vector3D.java rename to src/main/java/org/apache/commons/math/geometry/euclidean/threeD/Vector3D.java index 2d915e570..5b50c2265 100644 --- a/src/main/java/org/apache/commons/math/geometry/Vector3D.java +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/Vector3D.java @@ -15,7 +15,7 @@ * limitations under the License. */ -package org.apache.commons.math.geometry; +package org.apache.commons.math.geometry.euclidean.threeD; import java.io.Serializable; diff --git a/src/main/java/org/apache/commons/math/geometry/Vector3DFormat.java b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/Vector3DFormat.java similarity index 99% rename from src/main/java/org/apache/commons/math/geometry/Vector3DFormat.java rename to src/main/java/org/apache/commons/math/geometry/euclidean/threeD/Vector3DFormat.java index ff63a0e94..2aff0925e 100644 --- a/src/main/java/org/apache/commons/math/geometry/Vector3DFormat.java +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/Vector3DFormat.java @@ -15,7 +15,7 @@ * limitations under the License. */ -package org.apache.commons.math.geometry; +package org.apache.commons.math.geometry.euclidean.threeD; import java.text.FieldPosition; import java.text.NumberFormat; diff --git a/src/main/java/org/apache/commons/math/geometry/package.html b/src/main/java/org/apache/commons/math/geometry/euclidean/threeD/package.html similarity index 100% rename from src/main/java/org/apache/commons/math/geometry/package.html rename to src/main/java/org/apache/commons/math/geometry/euclidean/threeD/package.html diff --git a/src/main/java/org/apache/commons/math/geometry/euclidean/twoD/Line.java b/src/main/java/org/apache/commons/math/geometry/euclidean/twoD/Line.java new file mode 100644 index 000000000..5e8754e10 --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/twoD/Line.java @@ -0,0 +1,512 @@ +/* + * 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.math.geometry.euclidean.twoD; + +import java.awt.geom.AffineTransform; + +import org.apache.commons.math.exception.MathIllegalArgumentException; +import org.apache.commons.math.exception.util.LocalizedFormats; +import org.apache.commons.math.geometry.euclidean.oneD.IntervalsSet; +import org.apache.commons.math.geometry.euclidean.oneD.OrientedPoint; +import org.apache.commons.math.geometry.euclidean.oneD.Point1D; +import org.apache.commons.math.geometry.partitioning.BSPTree; +import org.apache.commons.math.geometry.partitioning.Hyperplane; +import org.apache.commons.math.geometry.partitioning.Point; +import org.apache.commons.math.geometry.partitioning.Region; +import org.apache.commons.math.geometry.partitioning.SubHyperplane; +import org.apache.commons.math.geometry.partitioning.SubSpace; +import org.apache.commons.math.geometry.partitioning.Transform; +import org.apache.commons.math.util.FastMath; +import org.apache.commons.math.util.MathUtils; + +/** This class represents an oriented line in the 2D plane. + + *An oriented line can be defined either by prolongating a line + * segment between two points past these points, or by one point and + * an angular direction (in trigonometric orientation).
+ + *Since it is oriented the two half planes at its two sides are + * unambiguously identified as a left half plane and a right half + * plane. This can be used to identify the interior and the exterior + * in a simple way by local properties only when part of a line is + * used to define part of a polygon boundary.
+ + *A line can also be used to completely define a reference frame + * in the plane. It is sufficient to select one specific point in the + * line (the orthogonal projection of the original reference frame on + * the line) and to use the unit vector in the line direction and the + * orthogonal vector oriented from left half plane to right half + * plane. We define two coordinates by the process, the + * abscissa along the line, and the offset across + * the line. All points of the plane are uniquely identified by these + * two coordinates. The line is the set of points at zero offset, the + * left half plane is the set of points with negative offsets and the + * right half plane is the set of points with positive offsets.
+ + * @version $Revision$ $Date$ + */ +public class Line implements Hyperplane { + + /** Angle with respect to the abscissa axis. */ + private double angle; + + /** Cosine of the line angle. */ + private double cos; + + /** Sine of the line angle. */ + private double sin; + + /** Offset of the frame origin. */ + private double originOffset; + + /** Build a line from two points. + *The line is oriented from p1 to p2
+ * @param p1 first point + * @param p2 second point + */ + public Line(final Point2D p1, final Point2D p2) { + reset(p1, p2); + } + + /** Build a line from a point and an angle. + * @param p point belonging to the line + * @param angle angle of the line with respect to abscissa axis + */ + public Line(final Point2D p, final double angle) { + reset(p, angle); + } + + /** Build a line from its internal characteristics. + * @param angle angle of the line with respect to abscissa axis + * @param cos cosine of the angle + * @param sin sine of the angle + * @param originOffset offset of the origin + */ + private Line(final double angle, final double cos, final double sin, final double originOffset) { + this.angle = angle; + this.cos = cos; + this.sin = sin; + this.originOffset = originOffset; + } + + /** Copy constructor. + *The created instance is completely independant from the + * original instance, it is a deep copy.
+ * @param line line to copy + */ + public Line(final Line line) { + angle = MathUtils.normalizeAngle(line.angle, FastMath.PI); + cos = FastMath.cos(angle); + sin = FastMath.sin(angle); + originOffset = line.originOffset; + } + + /** {@inheritDoc} */ + public Hyperplane copySelf() { + return new Line(this); + } + + /** Reset the instance as if built from two points. + *The line is oriented from p1 to p2
+ * @param p1 first point + * @param p2 second point + */ + public void reset(final Point2D p1, final Point2D p2) { + final double dx = p2.x - p1.x; + final double dy = p2.y - p1.y; + final double d = FastMath.hypot(dx, dy); + if (d == 0.0) { + angle = 0.0; + cos = 1.0; + sin = 0.0; + originOffset = p1.y; + } else { + angle = FastMath.PI + FastMath.atan2(-dy, -dx); + cos = FastMath.cos(angle); + sin = FastMath.sin(angle); + originOffset = (p2.x * p1.y - p1.x * p2.y) / d; + } + } + + /** Reset the instance as if built from a line and an angle. + * @param p point belonging to the line + * @param alpha angle of the line with respect to abscissa axis + */ + public void reset(final Point2D p, final double alpha) { + this.angle = MathUtils.normalizeAngle(alpha, FastMath.PI); + cos = FastMath.cos(this.angle); + sin = FastMath.sin(this.angle); + originOffset = cos * p.y - sin * p.x; + } + + /** Revert the instance. + */ + public void revertSelf() { + if (angle < FastMath.PI) { + angle += FastMath.PI; + } else { + angle -= FastMath.PI; + } + cos = -cos; + sin = -sin; + originOffset = -originOffset; + } + + /** Get the reverse of the instance. + *Get a line with reversed orientation with respect to the + * instance. A new object is built, the instance is untouched.
+ * @return a new line, with orientation opposite to the instance orientation + */ + public Line getReverse() { + return new Line((angle < FastMath.PI) ? (angle + FastMath.PI) : (angle - FastMath.PI), + -cos, -sin, -originOffset); + } + + /** Transform a 2D space point into a line point. + * @param point 2D point (must be a {@link Point2D Point2D} + * instance) + * @return line point corresponding to the 2D point (really a {@link + * org.apache.commons.math.geometry.euclidean.oneD.Point1D Point1D} instance) + * @see #toSpace + */ + public Point toSubSpace(final Point point) { + final Point2D p2D = (Point2D) point; + return new Point1D(cos * p2D.x + sin * p2D.y); + } + + /** Get one point from the line. + * @param point desired abscissa for the point (must be a {@link + * org.apache.commons.math.geometry.euclidean.oneD.Point1D Point1D} instance) + * @return line point at specified abscissa (really a {@link Point2D + * Point2D} instance) + */ + public Point toSpace(final Point point) { + final double abscissa = ((Point1D) point).getAbscissa(); + return new Point2D(abscissa * cos - originOffset * sin, + abscissa * sin + originOffset * cos); + } + + /** Get the intersection point of the instance and another line. + * @param other other line + * @return intersection point of the instance and the other line + * (really a {@link Point2D Point2D} instance) + */ + public SubSpace intersection(final Hyperplane other) { + final Line otherL = (Line) other; + final double d = sin * otherL.cos - otherL.sin * cos; + if (FastMath.abs(d) < 1.0e-10) { + return null; + } + return new Point2D((cos * otherL.originOffset - otherL.cos * originOffset) / d, + (sin * otherL.originOffset - otherL.sin * originOffset) / d); + } + + /** Build a region covering the whole hyperplane. + * @return a region covering the whole hyperplane + */ + public Region wholeHyperplane() { + return new IntervalsSet(); + } + + /** Build a region covering the whole space. + * @return a region containing the instance (really a {@link + * PolygonsSet PolygonsSet} instance) + */ + public Region wholeSpace() { + return new PolygonsSet(); + } + + /** Get the offset (oriented distance) of a parallel line. + *This method should be called only for parallel lines otherwise + * the result is not meaningful.
+ *The offset is 0 if both lines are the same, it is + * positive if the line is on the right side of the instance and + * negative if it is on the left side, according to its natural + * orientation.
+ * @param line line to check + * @return offset of the line + */ + public double getOffset(final Line line) { + return originOffset + + ((cos * line.cos + sin * line.sin > 0) ? -line.originOffset : line.originOffset); + } + + /** Get the offset (oriented distance) of a point to the line. + *The offset is 0 if the point belongs to the line, it is + * positive if the point is on the right side of the line and + * negative if it is on the left side, according to its natural + * orientation.
+ * @param point point to check (must be a {@link Point2D Point2D} instance) + * @return offset of the point + */ + public double getOffset(final Point point) { + final Point2D p2D = (Point2D) point; + return sin * p2D.x - cos * p2D.y + originOffset; + } + + /** Check if the instance has the same orientation as another hyperplane. + *This method is expected to be called on parallel hyperplanes + * (i.e. when the {@link #side side} method would return {@link + * org.apache.commons.math.geometry.partitioning.Hyperplane.Side#HYPER HYPER} + * for some sub-hyperplane having the specified hyperplane + * as its underlying hyperplane). The method should not + * re-check for parallelism, only for orientation, typically by + * testing something like the sign of the dot-products of + * normals.
+ * @param other other hyperplane to check against the instance + * @return true if the instance and the other hyperplane have + * the same orientation + */ + public boolean sameOrientationAs(final Hyperplane other) { + final Line otherL = (Line) other; + return (sin * otherL.sin + cos * otherL.cos) >= 0.0; + } + + /** Get one point from the plane. + * @param abscissa desired abscissa for the point + * @param offset desired offset for the point + * @return one point in the plane, with given abscissa and offset + * relative to the line + */ + public Point2D getPointAt(final Point1D abscissa, final double offset) { + final double x = abscissa.getAbscissa(); + final double dOffset = offset - originOffset; + return new Point2D(x * cos + dOffset * sin, x * sin - dOffset * cos); + } + + /** Check if the line contains a point. + * @param p point to check + * @return true if p belongs to the line + */ + public boolean contains(final Point2D p) { + return FastMath.abs(getOffset(p)) < 1.0e-10; + } + + /** Check the instance is parallel to another line. + * @param line other line to check + * @return true if the instance is parallel to the other line + * (they can have either the same or opposite orientations) + */ + public boolean isParallelTo(final Line line) { + return FastMath.abs(sin * line.cos - cos * line.sin) < 1.0e-10; + } + + /** Translate the line to force it passing by a point. + * @param p point by which the line should pass + */ + public void translateToPoint(final Point2D p) { + originOffset = cos * p.y - sin * p.x; + } + + /** Get the angle of the line. + * @return the angle of the line with respect to the abscissa axis + */ + public double getAngle() { + return MathUtils.normalizeAngle(angle, FastMath.PI); + } + + /** Set the angle of the line. + * @param angle new angle of the line with respect to the abscissa axis + */ + public void setAngle(final double angle) { + this.angle = MathUtils.normalizeAngle(angle, FastMath.PI); + cos = FastMath.cos(this.angle); + sin = FastMath.sin(this.angle); + } + + /** Get the offset of the origin. + * @return the offset of the origin + */ + public double getOriginOffset() { + return originOffset; + } + + /** Set the offset of the origin. + * @param offset offset of the origin + */ + public void setOriginOffset(final double offset) { + originOffset = offset; + } + + /** Compute the relative position of a sub-hyperplane with respect + * to the instance. + * @param sub sub-hyperplane to check + * @return one of {@link org.apache.commons.math.geometry.partitioning.Hyperplane.Side#PLUS PLUS}, + * {@link org.apache.commons.math.geometry.partitioning.Hyperplane.Side#MINUS MINUS}, + * {@link org.apache.commons.math.geometry.partitioning.Hyperplane.Side#BOTH BOTH}, + * {@link org.apache.commons.math.geometry.partitioning.Hyperplane.Side#HYPER HYPER} + */ + public Side side(final SubHyperplane sub) { + + final Hyperplane otherHyp = sub.getHyperplane(); + final Point2D crossing = (Point2D) intersection(otherHyp); + + if (crossing == null) { + // the lines are parallel, + final double global = getOffset((Line) otherHyp); + return (global < -1.0e-10) ? Side.MINUS : ((global > 1.0e-10) ? Side.PLUS : Side.HYPER); + } + + // the lines do intersect + final boolean direct = FastMath.sin(((Line) otherHyp).angle - angle) < 0; + final Point1D x = (Point1D) otherHyp.toSubSpace(crossing); + return sub.getRemainingRegion().side(new OrientedPoint(x, direct)); + + } + + /** Split a sub-hyperplane in two parts by the instance. + * @param sub sub-hyperplane to split + * @return an object containing both the part of the sub-hyperplane + * on the plus side of the instance and the part of the + * sub-hyperplane on the minus side of the instance + */ + public SplitSubHyperplane split(final SubHyperplane sub) { + + final Line otherLine = (Line) sub.getHyperplane(); + final Point2D crossing = (Point2D) intersection(otherLine); + + if (crossing == null) { + // the lines are parallel + final double global = getOffset(otherLine); + return (global < -1.0e-10) ? + new SplitSubHyperplane(null, sub) : + new SplitSubHyperplane(sub, null); + } + + // the lines do intersect + final boolean direct = FastMath.sin(otherLine.angle - angle) < 0; + final Point1D x = (Point1D) otherLine.toSubSpace(crossing); + final SubHyperplane subPlus = new SubHyperplane(new OrientedPoint(x, !direct)); + final SubHyperplane subMinus = new SubHyperplane(new OrientedPoint(x, direct)); + + final BSPTree splitTree = + sub.getRemainingRegion().getTree(false).split(subMinus); + final BSPTree plusTree = Region.isEmpty(splitTree.getPlus()) ? + new BSPTree(Boolean.FALSE) : + new BSPTree(subPlus, new BSPTree(Boolean.FALSE), + splitTree.getPlus(), null); + final BSPTree minusTree = Region.isEmpty(splitTree.getMinus()) ? + new BSPTree(Boolean.FALSE) : + new BSPTree(subMinus, new BSPTree(Boolean.FALSE), + splitTree.getMinus(), null); + + return new SplitSubHyperplane(new SubHyperplane(otherLine.copySelf(), + new IntervalsSet(plusTree)), + new SubHyperplane(otherLine.copySelf(), + new IntervalsSet(minusTree))); + + } + + /** Get a {@link org.apache.commons.math.geometry.partitioning.Transform + * Transform} embedding an affine transform. + * @param transform affine transform to embed (must be inversible + * otherwise the {@link + * org.apache.commons.math.geometry.partitioning.Transform#apply(Hyperplane) + * apply(Hyperplane)} method would work only for some lines, and + * fail for other ones) + * @return a new transform that can be applied to either {@link + * Point2D Point2D}, {@link Line Line} or {@link + * org.apache.commons.math.geometry.partitioning.SubHyperplane + * SubHyperplane} instances + * @exception MathIllegalArgumentException if the transform is non invertible + */ + public static Transform getTransform(final AffineTransform transform) throws MathIllegalArgumentException { + return new LineTransform(transform); + } + + /** Class embedding an affine transform. + *This class is used in order to apply an affine transform to a + * line. Using a specific object allow to perform some computations + * on the transform only once even if the same transform is to be + * applied to a large number of lines (for example to a large + * polygon)./
+ */ + private static class LineTransform implements Transform { + + // CHECKSTYLE: stop JavadocVariable check + private double cXX; + private double cXY; + private double cX1; + private double cYX; + private double cYY; + private double cY1; + + private double c1Y; + private double c1X; + private double c11; + // CHECKSTYLE: resume JavadocVariable check + + /** Build an affine line transform from a n {@code AffineTransform}. + * @param transform transform to use (must be invertible otherwise + * the {@link LineTransform#apply(Hyperplane)} method would work + * only for some lines, and fail for other ones) + * @exception MathIllegalArgumentException if the transform is non invertible + */ + public LineTransform(final AffineTransform transform) throws MathIllegalArgumentException { + + final double[] m = new double[6]; + transform.getMatrix(m); + cXX = m[0]; + cXY = m[2]; + cX1 = m[4]; + cYX = m[1]; + cYY = m[3]; + cY1 = m[5]; + + c1Y = cXY * cY1 - cYY * cX1; + c1X = cXX * cY1 - cYX * cX1; + c11 = cXX * cYY - cYX * cXY; + + if (FastMath.abs(c11) < 1.0e-20) { + throw new MathIllegalArgumentException(LocalizedFormats.NON_INVERTIBLE_TRANSFORM); + } + + } + + /** {@inheritDoc} */ + public Point apply(final Point point) { + final Point2D p2D = (Point2D) point; + final double x = p2D.getX(); + final double y = p2D.getY(); + return new Point2D(cXX * x + cXY * y + cX1, + cYX * x + cYY * y + cY1); + } + + /** {@inheritDoc} */ + public Hyperplane apply(final Hyperplane hyperplane) { + final Line line = (Line) hyperplane; + final double rOffset = c1X * line.cos + c1Y * line.sin + c11 * line.originOffset; + final double rCos = cXX * line.cos + cXY * line.sin; + final double rSin = cYX * line.cos + cYY * line.sin; + final double inv = 1.0 / FastMath.sqrt(rSin * rSin + rCos * rCos); + return new Line(FastMath.PI + FastMath.atan2(-rSin, -rCos), + inv * rCos, inv * rSin, + inv * rOffset); + } + + /** {@inheritDoc} */ + public SubHyperplane apply(final SubHyperplane sub, + final Hyperplane original, final Hyperplane transformed) { + final OrientedPoint op = (OrientedPoint) sub.getHyperplane(); + final Point1D newLoc = + (Point1D) transformed.toSubSpace(apply(original.toSpace(op.getLocation()))); + return new SubHyperplane(new OrientedPoint(newLoc, op.isDirect())); + } + + } + +} diff --git a/src/main/java/org/apache/commons/math/geometry/euclidean/twoD/NestedLoops.java b/src/main/java/org/apache/commons/math/geometry/euclidean/twoD/NestedLoops.java new file mode 100644 index 000000000..e9be044f2 --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/twoD/NestedLoops.java @@ -0,0 +1,190 @@ +/* + * 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.math.geometry.euclidean.twoD; + +import java.util.ArrayList; +import java.util.Arrays; +import java.util.Iterator; + +import org.apache.commons.math.exception.MathIllegalArgumentException; +import org.apache.commons.math.exception.util.LocalizedFormats; +import org.apache.commons.math.geometry.euclidean.oneD.OrientedPoint; +import org.apache.commons.math.geometry.euclidean.oneD.Point1D; +import org.apache.commons.math.geometry.partitioning.Hyperplane; +import org.apache.commons.math.geometry.partitioning.Region; +import org.apache.commons.math.geometry.partitioning.SubHyperplane; + +/** This class represent a tree of nested 2D boundary loops. + + *
This class is used during Piece instances construction. + * Beams are built using the outline edges as + * representative of facets, the orientation of these facets are + * meaningful. However, we want to allow the user to specify its + * outline loops without having to take care of this orientation. This + * class is devoted to correct mis-oriented loops.
+ + *
Orientation is computed assuming the piece is finite, i.e. the + * outermost loops have their exterior side facing points at infinity, + * and hence are oriented counter-clockwise. The orientation of + * internal loops is computed as the reverse of the orientation of + * their immediate surrounding loop.
+ + * @version $Revision$ $Date$ + */ +class NestedLoops { + + /** Boundary loop. */ + private Point2D[] loop; + + /** Surrounded loops. */ + private ArrayListBuild an empty tree of nested loops. This instance will become + * the root node of a complete tree, it is not associated with any + * loop by itself, the outermost loops are in the root tree child + * nodes.
+ */ + public NestedLoops() { + surrounded = new ArrayListBuild a tree node with neither parent nor children
+ * @param loop boundary loop (will be reversed in place if needed) + * @exception MathIllegalArgumentException if an outline has an open boundary loop + */ + private NestedLoops(final Point2D[] loop) throws MathIllegalArgumentException { + + if (loop[0] == null) { + throw new MathIllegalArgumentException(LocalizedFormats.OUTLINE_BOUNDARY_LOOP_OPEN); + } + + this.loop = loop; + surrounded = new ArrayListThis is this method that really inverts the loops that where + * provided through the {@link #add(Point2D[]) add} method if + * they are mis-oriented
+ */ + public void correctOrientation() { + for (NestedLoops child : surrounded) { + child.setClockWise(true); + } + } + + /** Set the loop orientation. + * @param clockwise if true, the loop should be set to clockwise + * orientation + */ + private void setClockWise(final boolean clockwise) { + + if (originalIsClockwise ^ clockwise) { + // we need to inverse the original loop + int min = -1; + int max = loop.length; + while (++min < --max) { + final Point2D tmp = loop[min]; + loop[min] = loop[max]; + loop[max] = tmp; + } + } + + // go deeper in the tree + for (final NestedLoops child : surrounded) { + child.setClockWise(!clockwise); + } + + } + +} diff --git a/src/main/java/org/apache/commons/math/geometry/euclidean/twoD/Point2D.java b/src/main/java/org/apache/commons/math/geometry/euclidean/twoD/Point2D.java new file mode 100644 index 000000000..eb8fa2da9 --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/twoD/Point2D.java @@ -0,0 +1,72 @@ +/* + * 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.math.geometry.euclidean.twoD; + +import org.apache.commons.math.geometry.partitioning.Point; +import org.apache.commons.math.geometry.partitioning.SubSpace; + +/** This class represents a 2D point. + *Instances of this class are guaranteed to be immutable.
+ * @version $Revision$ $Date$ + */ +public class Point2D extends java.awt.geom.Point2D.Double implements Point, SubSpace { + + /** Point at undefined (NaN) coordinates. */ + public static final Point2D UNDEFINED = new Point2D(java.lang.Double.NaN, java.lang.Double.NaN); + + /** Serializable UID. */ + private static final long serialVersionUID = 8883702098988517151L; + + /** Build a point with default coordinates. + */ + public Point2D() { + } + + /** Build a point from its coordinates. + * @param x abscissa + * @param y ordinate + */ + public Point2D(final double x, final double y) { + super(x, y); + } + + /** Build a point from a java awt point. + * @param point java awt point + */ + public Point2D(final java.awt.geom.Point2D.Double point) { + super(point.x, point.y); + } + + /** Transform a 2D space point into a sub-space point. + * @param point 2D point of the space + * @return always return null + * @see #toSpace + */ + public Point toSubSpace(final Point point) { + return null; + } + + /** Transform a sub-space point into a space point. + * @param point ignored parameter + * @return always return the instance + * @see #toSubSpace + */ + public Point toSpace(final Point point) { + return this; + } + +} diff --git a/src/main/java/org/apache/commons/math/geometry/euclidean/twoD/PolygonsSet.java b/src/main/java/org/apache/commons/math/geometry/euclidean/twoD/PolygonsSet.java new file mode 100644 index 000000000..728a357bb --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/twoD/PolygonsSet.java @@ -0,0 +1,343 @@ +/* + * 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.math.geometry.euclidean.twoD; + +import java.util.ArrayList; +import java.util.Arrays; +import java.util.Collection; +import java.util.List; + +import org.apache.commons.math.geometry.euclidean.oneD.Point1D; +import org.apache.commons.math.geometry.partitioning.BSPTree; +import org.apache.commons.math.geometry.partitioning.Hyperplane; +import org.apache.commons.math.geometry.partitioning.Region; +import org.apache.commons.math.geometry.partitioning.SubHyperplane; +import org.apache.commons.math.geometry.partitioning.utilities.AVLTree; +import org.apache.commons.math.util.FastMath; + +/** This class represents a 2D region: a set of polygons. + * @version $Revision$ $Date$ + */ +public class PolygonsSet extends Region { + + /** Vertices organized as boundary loops. */ + private Point2D[][] vertices; + + /** Build a polygons set representing the whole real line. + */ + public PolygonsSet() { + super(); + } + + /** Build a polygons set from a BSP tree. + *The leaf nodes of the BSP tree must have a + * {@code Boolean} attribute representing the inside status of + * the corresponding cell (true for inside cells, false for outside + * cells). In order to avoid building too many small objects, it is + * recommended to use the predefined constants + * {@code Boolean.TRUE} and {@code Boolean.FALSE}
+ * @param tree inside/outside BSP tree representing the region + */ + public PolygonsSet(final BSPTree tree) { + super(tree); + } + + /** Build a polygons set from a Boundary REPresentation (B-rep). + *The boundary is provided as a collection of {@link + * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the + * interior part of the region on its minus side and the exterior on + * its plus side.
+ *The boundary elements can be in any order, and can form + * several non-connected sets (like for example polygons with holes + * or a set of disjoints polyhedrons considered as a whole). In + * fact, the elements do not even need to be connected together + * (their topological connections are not used here). However, if the + * boundary does not really separate an inside open from an outside + * open (open having here its topological meaning), then subsequent + * calls to the {@link + * Region#checkPoint(org.apache.commons.math.geometry.partitioning.Point) + * checkPoint} method will not be meaningful anymore.
+ *If the boundary is empty, the region will represent the whole + * space.
+ * @param boundary collection of boundary elements, as a + * collection of {@link SubHyperplane SubHyperplane} objects + */ + public PolygonsSet(final CollectionThe polygon boundary can be represented as an array of loops, + * each loop being itself an array of vertices.
+ *In order to identify open loops which start and end by + * infinite edges, the open loops arrays start with a null point. In + * this case, the first non null point and the last point of the + * array do not represent real vertices, they are dummy points + * intended only to get the direction of the first and last edge. An + * open loop consisting of a single infinite line will therefore be + * represented by a three elements array with one null point + * followed by two dummy points. The open loops are always the first + * ones in the loops array.
+ *If the polygon has no boundary at all, a zero length loop + * array will be returned.
+ *All line segments in the various loops have the inside of the + * region on their left side and the outside on their right side + * when moving in the underlying line direction. This means that + * closed loops surrounding finite areas obey the direct + * trigonometric orientation.
+ * @return vertices of the polygon, organized as oriented boundary + * loops with the open loops first (the returned value is guaranteed + * to be non-null) + */ + public Point2D[][] getVertices() { + if (vertices == null) { + if (getTree(false).getCut() == null) { + vertices = new Point2D[0][]; + } else { + + // sort the segmfinal ents according to their start point + final SegmentsBuilder visitor = new SegmentsBuilder(); + getTree(true).visit(visitor); + final AVLTree+ * The object built is not a real segment, only the sorting key is used to + * allow searching in the neighborhood of a point. This is an horrible hack ... + *
+ * @param start start point of the segment + * @param dx abscissa offset from the start point + * @param dy ordinate offset from the start point + */ + public Segment(final Point2D start, final double dx, final double dy) { + this.start = null; + this.end = null; + this.line = null; + sortingKey = new OrderedTuple(start.x + dx, start.y + dy); + } + + /** Get the start point of the segment. + * @return start point of the segment + */ + public Point2D getStart() { + return start; + } + + /** Get the end point of the segment. + * @return end point of the segment + */ + public Point2D getEnd() { + return end; + } + + /** Get the line containing the segment. + * @return line containing the segment + */ + public Line getLine() { + return line; + } + + /** {@inheritDoc} */ + public int compareTo(final Segment o) { + return sortingKey.compareTo(o.sortingKey); + } + + /** {@inheritDoc} */ + @Override + public boolean equals(final Object other) { + if (this == other) { + return true; + } else if (other instanceof Segment) { + return compareTo((Segment) other) == 0; + } else { + return false; + } + } + + /** {@inheritDoc} */ + @Override + public int hashCode() { + return start.hashCode() ^ end.hashCode() ^ line.hashCode() ^ sortingKey.hashCode(); + } + +} diff --git a/src/main/java/org/apache/commons/math/geometry/euclidean/twoD/SegmentBuilder.java b/src/main/java/org/apache/commons/math/geometry/euclidean/twoD/SegmentBuilder.java new file mode 100644 index 000000000..bf77c2448 --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/euclidean/twoD/SegmentBuilder.java @@ -0,0 +1,90 @@ +/* + * 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.math.geometry.euclidean.twoD; + +import java.util.List; + +import org.apache.commons.math.geometry.euclidean.oneD.Interval; +import org.apache.commons.math.geometry.euclidean.oneD.IntervalsSet; +import org.apache.commons.math.geometry.euclidean.oneD.Point1D; +import org.apache.commons.math.geometry.partitioning.BSPTree; +import org.apache.commons.math.geometry.partitioning.BSPTreeVisitor; +import org.apache.commons.math.geometry.partitioning.Region.BoundaryAttribute; +import org.apache.commons.math.geometry.partitioning.SubHyperplane; +import org.apache.commons.math.geometry.partitioning.utilities.AVLTree; + +/** Visitor building segments. + * @version $Revision$ $Date$ + */ +class SegmentsBuilder implements BSPTreeVisitor { + + /** Sorted segments. */ + private AVLTree+This package provides basic 2D geometry components. +
+ + diff --git a/src/main/java/org/apache/commons/math/geometry/partitioning/BSPTree.java b/src/main/java/org/apache/commons/math/geometry/partitioning/BSPTree.java new file mode 100644 index 000000000..18c105309 --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/partitioning/BSPTree.java @@ -0,0 +1,631 @@ +/* + * 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.math.geometry.partitioning; + +import org.apache.commons.math.geometry.partitioning.Hyperplane.Side; +import org.apache.commons.math.util.FastMath; + +/** This class represent a Binary Space Partition tree. + + *BSP trees are an efficient way to represent space partitions and + * to associate attributes with each cell. Each node in a BSP tree + * represents a convex region which is partitioned in two convex + * sub-regions at each side of a cut hyperplane. The root tree + * contains the complete space.
+ + *The main use of such partitions is to use a boolean attribute to + * define an inside/outside property, hence representing arbitrary + * polytopes (line segments in 1D, polygons in 2D and polyhedrons in + * 3D) and to operate on them.
+ + *Another example would be to represent Voronoi tesselations, the + * attribute of each cell holding the defining point of the cell.
+ + *The application-defined attributes are shared among copied + * instances and propagated to split parts. These attributes are not + * used by the BSP-tree algorithms themselves, so the application can + * use them for any purpose. Since the tree visiting method holds + * internal and leaf nodes differently, it is possible to use + * different classes for internal nodes attributes and leaf nodes + * attributes. This should be used with care, though, because if the + * tree is modified in any way after attributes have been set, some + * internal nodes may become leaf nodes and some leaf nodes may become + * internal nodes.
+ + *One of the main sources for the development of this package was + * Bruce Naylor, John Amanatides and William Thibault paper Merging + * BSP Trees Yields Polyhedral Set Operations Proc. Siggraph '90, + * Computer Graphics 24(4), August 1990, pp 115-124, published by the + * Association for Computing Machinery (ACM).
+ + * @version $Revision$ $Date$ + */ +public class BSPTree { + + /** Cut sub-hyperplane. */ + private SubHyperplane cut; + + /** Tree at the plus side of the cut hyperplane. */ + private BSPTree plus; + + /** Tree at the minus side of the cut hyperplane. */ + private BSPTree minus; + + /** Parent tree. */ + private BSPTree parent; + + /** Application-defined attribute. */ + private Object attribute; + + /** Build a tree having only one root cell representing the whole space. + */ + public BSPTree() { + cut = null; + plus = null; + minus = null; + parent = null; + attribute = null; + } + + /** Build a tree having only one root cell representing the whole space. + * @param attribute attribute of the tree (may be null) + */ + public BSPTree(final Object attribute) { + cut = null; + plus = null; + minus = null; + parent = null; + this.attribute = attribute; + } + + /** Build a BSPTree from its underlying elements. + *This method does not perform any verification on + * consistency of its arguments, it should therefore be used only + * when then caller knows what it is doing.
+ *This method is mainly useful kto build trees + * bottom-up. Building trees top-down is realized with the help of + * method {@link #insertCut insertCut}.
+ * @param cut cut sub-hyperplane for the tree + * @param plus plus side sub-tree + * @param minus minus side sub-tree + * @param attribute attribute associated with the node (may be null) + * @see #insertCut + */ + public BSPTree(final SubHyperplane cut, final BSPTree plus, final BSPTree minus, + final Object attribute) { + this.cut = cut; + this.plus = plus; + this.minus = minus; + this.parent = null; + this.attribute = attribute; + plus.parent = this; + minus.parent = this; + } + + /** Insert a cut sub-hyperplane in a node. + *The sub-tree starting at this node will be completely + * overwritten. The new cut sub-hyperplane will be built from the + * intersection of the provided hyperplane with the cell. If the + * hyperplane does intersect the cell, the cell will have two + * children cells with {@code null} attributes on each side of + * the inserted cut sub-hyperplane. If the hyperplane does not + * intersect the cell then no cut hyperplane will be + * inserted and the cell will be changed to a leaf cell. The + * attribute of the node is never changed.
+ *This method is mainly useful when called on leaf nodes + * (i.e. nodes for which {@link #getCut getCut} returns + * {@code null}), in this case it provides a way to build a + * tree top-down (whereas the {@link #BSPTree(SubHyperplane, + * BSPTree, BSPTree, Object) 4 arguments constructor} is devoted to + * build trees bottom-up).
+ * @param hyperplane hyperplane to insert, it will be chopped in + * order to fit in the cell defined by the parent nodes of the + * instance + * @return true if a cut sub-hyperplane has been inserted (i.e. if + * the cell now has two leaf child nodes) + * @see #BSPTree(SubHyperplane, BSPTree, BSPTree, Object) + */ + public boolean insertCut(final Hyperplane hyperplane) { + + if (cut != null) { + plus.parent = null; + minus.parent = null; + } + + final SubHyperplane chopped = fitToCell(new SubHyperplane(hyperplane)); + if (chopped.getRemainingRegion().isEmpty()) { + cut = null; + plus = null; + minus = null; + return false; + } + + cut = chopped; + plus = new BSPTree(); + plus.parent = this; + minus = new BSPTree(); + minus.parent = this; + return true; + + } + + /** Copy the instance. + *The instance created is completely independant of the original + * one. A deep copy is used, none of the underlying objects are + * shared (except for the nodes attributes and immutable + * objects).
+ * @return a new tree, copy of the instance + */ + public BSPTree copySelf() { + + if (cut == null) { + return new BSPTree(attribute); + } + + return new BSPTree(cut.copySelf(), plus.copySelf(), minus.copySelf(), + attribute); + + } + + /** Get the cut sub-hyperplane. + * @return cut sub-hyperplane, null if this is a leaf tree + */ + public SubHyperplane getCut() { + return cut; + } + + /** Get the tree on the plus side of the cut hyperplane. + * @return tree on the plus side of the cut hyperplane, null if this + * is a leaf tree + */ + public BSPTree getPlus() { + return plus; + } + + /** Get the tree on the minus side of the cut hyperplane. + * @return tree on the minus side of the cut hyperplane, null if this + * is a leaf tree + */ + public BSPTree getMinus() { + return minus; + } + + /** Get the parent node. + * @return parent node, null if the node has no parents + */ + public BSPTree getParent() { + return parent; + } + + /** Associate an attribute with the instance. + * @param attribute attribute to associate with the node + * @see #getAttribute + */ + public void setAttribute(final Object attribute) { + this.attribute = attribute; + } + + /** Get the attribute associated with the instance. + * @return attribute associated with the node or null if no + * attribute has been explicitly set using the {@link #setAttribute + * setAttribute} method + * @see #setAttribute + */ + public Object getAttribute() { + return attribute; + } + + /** Visit the BSP tree nodes. + * @param visitor object visiting the tree nodes + */ + public void visit(final BSPTreeVisitor visitor) { + if (cut == null) { + visitor.visitLeafNode(this); + } else { + switch (visitor.visitOrder(this)) { + case PLUS_MINUS_SUB: + plus.visit(visitor); + minus.visit(visitor); + visitor.visitInternalNode(this); + break; + case PLUS_SUB_MINUS: + plus.visit(visitor); + visitor.visitInternalNode(this); + minus.visit(visitor); + break; + case MINUS_PLUS_SUB: + minus.visit(visitor); + plus.visit(visitor); + visitor.visitInternalNode(this); + break; + case MINUS_SUB_PLUS: + minus.visit(visitor); + visitor.visitInternalNode(this); + plus.visit(visitor); + break; + case SUB_PLUS_MINUS: + visitor.visitInternalNode(this); + plus.visit(visitor); + minus.visit(visitor); + break; + case SUB_MINUS_PLUS: + visitor.visitInternalNode(this); + minus.visit(visitor); + plus.visit(visitor); + break; + default: + throw new RuntimeException("internal error"); + } + + } + } + + /** Fit a sub-hyperplane inside the cell defined by the instance. + *Fitting is done by chopping off the parts of the + * sub-hyperplane that lie outside of the cell using the + * cut-hyperplanes of the parent nodes of the instance.
+ * @param sub sub-hyperplane to fit + * @return a new sub-hyperplane, gueranteed to have no part outside + * of the instance cell + */ + private SubHyperplane fitToCell(final SubHyperplane sub) { + SubHyperplane s = sub; + for (BSPTree tree = this; tree.parent != null; tree = tree.parent) { + if (tree == tree.parent.plus) { + s = tree.parent.cut.getHyperplane().split(s).getPlus(); + } else { + s = tree.parent.cut.getHyperplane().split(s).getMinus(); + } + } + return s; + } + + /** Get the cell to which a point belongs. + *If the returned cell is a leaf node the points belongs to the + * interior of the node, if the cell is an internal node the points + * belongs to the node cut sub-hyperplane.
+ * @param point point to check + * @return the tree cell to which the point belongs (can be + */ + public BSPTree getCell(final Point point) { + + if (cut == null) { + return this; + } + + // position of the point with respect to the cut hyperplane + final double offset = cut.getHyperplane().getOffset(point); + + if (FastMath.abs(offset) < 1.0e-10) { + return this; + } else if (offset <= 0) { + // point is on the minus side of the cut hyperplane + return minus.getCell(point); + } else { + // point is on the plus side of the cut hyperplane + return plus.getCell(point); + } + + } + + /** Perform condensation on a tree. + *The condensation operation is not recursive, it must be called + * explicitely from leaves to root.
+ */ + private void condense() { + if ((cut != null) && (plus.cut == null) && (minus.cut == null) && + (((plus.attribute == null) && (minus.attribute == null)) || + ((plus.attribute != null) && plus.attribute.equals(minus.attribute)))) { + attribute = (plus.attribute == null) ? minus.attribute : plus.attribute; + cut = null; + plus = null; + minus = null; + } + } + + /** Merge a BSP tree with the instance. + *All trees are modified (parts of them are reused in the new + * tree), it is the responsibility of the caller to ensure a copy + * has been done before if any of the former tree should be + * preserved, no such copy is done here!
+ *The algorithm used here is directly derived from the one + * described in the Naylor, Amanatides and Thibault paper (section + * III, Binary Partitioning of a BSP Tree).
+ * @param tree other tree to merge with the instance (will be + * unusable after the operation, as well as the + * instance itself) + * @param leafMerger object implementing the final merging phase + * (this is where the semantic of the operation occurs, generally + * depending on the attribute of the leaf node) + * @return a new tree, result ofinstance <op>
+ * tree, this value can be ignored if parentTree is not null
+ * since all connections have already been established
+ */
+ public BSPTree merge(final BSPTree tree, final LeafMerger leafMerger) {
+ return merge(tree, leafMerger, null, false);
+ }
+
+ /** Merge a BSP tree with the instance.
+ * @param tree other tree to merge with the instance (will be
+ * unusable after the operation, as well as the
+ * instance itself)
+ * @param leafMerger object implementing the final merging phase
+ * (this is where the semantic of the operation occurs, generally
+ * depending on the attribute of the leaf node)
+ * @param parentTree parent tree to connect to (may be null)
+ * @param isPlusChild if true and if parentTree is not null, the
+ * resulting tree should be the plus child of its parent, ignored if
+ * parentTree is null
+ * @return a new tree, result of instance <op>
+ * tree, this value can be ignored if parentTree is not null
+ * since all connections have already been established
+ */
+ private BSPTree merge(final BSPTree tree, final LeafMerger leafMerger,
+ final BSPTree parentTree, final boolean isPlusChild) {
+ if (cut == null) {
+ // cell/tree operation
+ return leafMerger.merge(this, tree, parentTree, isPlusChild, true);
+ } else if (tree.cut == null) {
+ // tree/cell operation
+ return leafMerger.merge(tree, this, parentTree, isPlusChild, false);
+ } else {
+ // tree/tree operation
+ final BSPTree merged = tree.split(cut);
+ if (parentTree != null) {
+ merged.parent = parentTree;
+ if (isPlusChild) {
+ parentTree.plus = merged;
+ } else {
+ parentTree.minus = merged;
+ }
+ }
+
+ // merging phase
+ plus.merge(merged.plus, leafMerger, merged, true);
+ minus.merge(merged.minus, leafMerger, merged, false);
+ merged.condense();
+ if (merged.cut != null) {
+ merged.cut =
+ merged.fitToCell(new SubHyperplane(merged.cut.getHyperplane()));
+ }
+
+ return merged;
+
+ }
+ }
+
+ /** This interface gather the merging operations between a BSP tree
+ * leaf and another BSP tree.
+ * As explained in Bruce Naylor, John Amanatides and William + * Thibault paper Merging + * BSP Trees Yields Polyhedral Set Operations, + * the operations on {@link BSPTree BSP trees} can be expressed as a + * generic recursive merging operation where only the final part, + * when one of the operand is a leaf, is specific to the real + * operation semantics. For example, a tree representing a region + * using a boolean attribute to identify inside cells and outside + * cells would use four different objects to implement the final + * merging phase of the four set operations union, intersection, + * difference and symmetric difference (exclusive or).
+ * @version $Revision$ $Date$ + */ + public static interface LeafMerger { + + /** Merge a leaf node and a tree node. + *This method is called at the end of a recursive merging + * resulting from a {@code tree1.merge(tree2, leafMerger)} + * call, when one of the sub-trees involved is a leaf (i.e. when + * its cut-hyperplane is null). This is the only place where the + * precise semantics of the operation are required. For all upper + * level nodes in the tree, the merging operation is only a + * generic partitioning algorithm.
+ *Since the final operation may be non-commutative, it is + * important to know if the leaf node comes from the instance tree + * ({@code tree1}) or the argument tree + * ({@code tree2}). The third argument of the method is + * devoted to this. It can be ignored for commutative + * operations.
+ *The {@link BSPTree#insertInTree BSPTree.insertInTree} method + * may be useful to implement this method.
+ * @param leaf leaf node (its cut hyperplane is guaranteed to be + * null) + * @param tree tree node (its cut hyperplane may be null or not) + * @param parentTree parent tree to connect to (may be null) + * @param isPlusChild if true and if parentTree is not null, the + * resulting tree should be the plus child of its parent, ignored if + * parentTree is null + * @param leafFromInstance if true, the leaf node comes from the + * instance tree ({@code tree1}) and the tree node comes from + * the argument tree ({@code tree2}) + * @return the BSP tree resulting from the merging (may be one of + * the arguments) + */ + BSPTree merge(BSPTree leaf, BSPTree tree, + BSPTree parentTree, boolean isPlusChild, + boolean leafFromInstance); + + } + + /** Split a BSP tree by an external sub-hyperplane. + *Split a tree in two halves, on each side of the + * sub-hyperplane. The instance is not modified.
+ *The tree returned is not upward-consistent: despite all of its + * sub-trees cut sub-hyperplanes (including its own cut + * sub-hyperplane) are bounded to the current cell, it is not + * attached to any parent tree yet. This tree is intended to be + * later inserted into an higher level tree.
+ *The algorithm used here is the one given in Naylor, Amanatides + * and Thibault paper (section III, Binary Partitioning of a BSP + * Tree).
+ * @param sub partitioning sub-hyperplane, must be already clipped + * to the convex region represented by the instance, will be used as + * the cut sub-hyperplane of the returned tree + * @return a tree having the specified sub-hyperplane as its cut + * sub-hyperplane, the two parts of the split instance as its two + * sub-trees and a null parent + */ + public BSPTree split(final SubHyperplane sub) { + + if (cut == null) { + return new BSPTree(sub, copySelf(), new BSPTree(attribute), null); + } + + final Hyperplane cHyperplane = cut.getHyperplane(); + final Hyperplane sHyperplane = sub.getHyperplane(); + switch (cHyperplane.side(sub)) { + case PLUS : + { // the partitioning sub-hyperplane is entirely in the plus sub-tree + final BSPTree split = plus.split(sub); + if (sHyperplane.side(cut) == Side.PLUS) { + split.plus = new BSPTree(cut.copySelf(), + split.plus, minus.copySelf(), attribute); + split.plus.condense(); + split.plus.parent = split; + } else { + split.minus = new BSPTree(cut.copySelf(), + split.minus, minus.copySelf(), attribute); + split.minus.condense(); + split.minus.parent = split; + } + return split; + } + case MINUS : + { // the partitioning sub-hyperplane is entirely in the minus sub-tree + final BSPTree split = minus.split(sub); + if (sHyperplane.side(cut) == Side.PLUS) { + split.plus = new BSPTree(cut.copySelf(), + plus.copySelf(), split.plus, attribute); + split.plus.condense(); + split.plus.parent = split; + } else { + split.minus = new BSPTree(cut.copySelf(), + plus.copySelf(), split.minus, attribute); + split.minus.condense(); + split.minus.parent = split; + } + return split; + } + case BOTH : + { + final Hyperplane.SplitSubHyperplane cutParts = sHyperplane.split(cut); + final Hyperplane.SplitSubHyperplane subParts = cHyperplane.split(sub); + final BSPTree split = new BSPTree(sub, + plus.split(subParts.getPlus()), + minus.split(subParts.getMinus()), + null); + split.plus.cut = cutParts.getPlus(); + split.minus.cut = cutParts.getMinus(); + final BSPTree tmp = split.plus.minus; + split.plus.minus = split.minus.plus; + split.plus.minus.parent = split.plus; + split.minus.plus = tmp; + split.minus.plus.parent = split.minus; + split.plus.condense(); + split.minus.condense(); + return split; + } + default : + return cHyperplane.sameOrientationAs(sHyperplane) ? + new BSPTree(sub, plus.copySelf(), minus.copySelf(), attribute) : + new BSPTree(sub, minus.copySelf(), plus.copySelf(), attribute); + } + + } + + /** Insert the instance into another tree. + *The instance itself is modified so its former parent should + * not be used anymore.
+ * @param parentTree parent tree to connect to (may be null) + * @param isPlusChild if true and if parentTree is not null, the + * resulting tree should be the plus child of its parent, ignored if + * parentTree is null + * @see LeafMerger + */ + public void insertInTree(final BSPTree parentTree, final boolean isPlusChild) { + + // set up parent/child links + parent = parentTree; + if (parentTree != null) { + if (isPlusChild) { + parentTree.plus = this; + } else { + parentTree.minus = this; + } + } + + // make sure the inserted tree lies in the cell defined by its parent nodes + if (cut != null) { + + // explore the parent nodes from here towards tree root + for (BSPTree tree = this; tree.parent != null; tree = tree.parent) { + + // this is an hyperplane of some parent node + final Hyperplane hyperplane = tree.parent.cut.getHyperplane(); + + // chop off the parts of the inserted tree that extend + // on the wrong side of this parent hyperplane + if (tree == tree.parent.plus) { + cut = hyperplane.split(cut).getPlus(); + plus.chopOffMinus(hyperplane); + minus.chopOffMinus(hyperplane); + } else { + cut = hyperplane.split(cut).getMinus(); + plus.chopOffPlus(hyperplane); + minus.chopOffPlus(hyperplane); + } + + } + + // since we may have drop some parts of the inserted tree, + // perform a condensation pass to keep the tree structure simple + condense(); + + } + + } + + /** Chop off parts of the tree. + *The instance is modified in place, all the parts that are on + * the minus side of the chopping hyperplane are disgarded, only the + * parts on the plus side remain.
+ * @param hyperplane chopping hyperplane + */ + private void chopOffMinus(final Hyperplane hyperplane) { + if (cut != null) { + cut = hyperplane.split(cut).getPlus(); + plus.chopOffMinus(hyperplane); + minus.chopOffMinus(hyperplane); + } + } + + /** Chop off parts of the tree. + *The instance is modified in place, all the parts that are on + * the plus side of the chopping hyperplane are disgarded, only the + * parts on the minus side remain.
+ * @param hyperplane chopping hyperplane + */ + private void chopOffPlus(final Hyperplane hyperplane) { + if (cut != null) { + cut = hyperplane.split(cut).getMinus(); + plus.chopOffPlus(hyperplane); + minus.chopOffPlus(hyperplane); + } + } + +} diff --git a/src/main/java/org/apache/commons/math/geometry/partitioning/BSPTreeVisitor.java b/src/main/java/org/apache/commons/math/geometry/partitioning/BSPTreeVisitor.java new file mode 100644 index 000000000..6cfd38cbd --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/partitioning/BSPTreeVisitor.java @@ -0,0 +1,110 @@ +/* + * 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.math.geometry.partitioning; + +/** This interface is used to visit {@link BSPTree BSP tree} nodes. + + *Navigation through {@link BSPTree BSP trees} can be done using + * two different point of views:
+ *-
+ *
- + * the first one is in a node-oriented way using the {@link + * BSPTree#getPlus}, {@link BSPTree#getMinus} and {@link + * BSPTree#getParent} methods. Terminal nodes without associated + * {@link SubHyperplane sub-hyperplanes} can be visited this way, + * there is no constraint in the visit order, and it is possible + * to visit either all nodes or only a subset of the nodes + * + *
- + * the second one is in a sub-hyperplane-oriented way using + * classes implementing this interface which obeys the visitor + * design pattern. The visit order is provided by the visitor as + * each node is first encountered. Each node is visited exactly + * once. + * + *
Before attempting to visit an internal node, this method is + * called to determine the desired ordering of the visit. It is + * guaranteed that this method will be called before {@link + * #visitInternalNode visitInternalNode} for a given node, it will be + * called exactly once for each internal node.
+ * @param node BSP node guaranteed to have a non null cut sub-hyperplane + * @return desired visit order, must be one of + * {@link Order#PLUS_MINUS_SUB}, {@link Order#PLUS_SUB_MINUS}, + * {@link Order#MINUS_PLUS_SUB}, {@link Order#MINUS_SUB_PLUS}, + * {@link Order#SUB_PLUS_MINUS}, {@link Order#SUB_MINUS_PLUS} + */ + Order visitOrder(BSPTree node); + + /** Visit a BSP tree node node having a non-null sub-hyperplane. + *It is guaranteed that this method will be called after {@link + * #visitOrder visitOrder} has been called for a given node, + * it wil be called exactly once for each internal node.
+ * @param node BSP node guaranteed to have a non null cut sub-hyperplane + * @see #visitLeafNode + */ + void visitInternalNode(BSPTree node); + + /** Visit a leaf BSP tree node node having a null sub-hyperplane. + * @param node leaf BSP node having a null sub-hyperplane + * @see #visitInternalNode + */ + void visitLeafNode(BSPTree node); + +} diff --git a/src/main/java/org/apache/commons/math/geometry/partitioning/Characterization.java b/src/main/java/org/apache/commons/math/geometry/partitioning/Characterization.java new file mode 100644 index 000000000..a96b5ba5b --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/partitioning/Characterization.java @@ -0,0 +1,90 @@ +/* + * 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.math.geometry.partitioning; + +/** Characterization of a sub-hyperplane. + * @version $Revision$ $Date$ + */ +class Characterization { + + /** Parts of the sub-hyperplane that have inside cells on the tested side. */ + private SubHyperplane in; + + /** Parts of the sub-hyperplane that have outside cells on the tested side. */ + private SubHyperplane out; + + /** Create an empty characterization of a sub-hyperplane. + */ + public Characterization() { + in = null; + out = null; + } + + /** Check if the sub-hyperplane that have inside cells on the tested side. + * @return true if the sub-hyperplane that have inside cells on the tested side + */ + public boolean hasIn() { + return (in != null) && (!in.getRemainingRegion().isEmpty()); + } + + /** Get the parts of the sub-hyperplane that have inside cells on the tested side. + * @return parts of the sub-hyperplane that have inside cells on the tested side + */ + public SubHyperplane getIn() { + return in; + } + + /** Check if the sub-hyperplane that have outside cells on the tested side. + * @return true if the sub-hyperplane that have outside cells on the tested side + */ + public boolean hasOut() { + return (out != null) && (!out.getRemainingRegion().isEmpty()); + } + + /** Get the parts of the sub-hyperplane that have outside cells on the tested side. + * @return parts of the sub-hyperplane that have outside cells on the tested side + */ + public SubHyperplane getOut() { + return out; + } + + /** Add a part of the sub-hyperplane known to have inside or outside cell on the tested side. + * @param sub part of the sub-hyperplane to add + * @param inside if true, the part added as an inside cell on the tested side, otherwise + * it has an outside cell on the tested side + */ + public void add(final SubHyperplane sub, final boolean inside) { + if (inside) { + if (in == null) { + in = sub; + } else { + in = new SubHyperplane(in.getHyperplane(), + Region.union(in.getRemainingRegion(), + sub.getRemainingRegion())); + } + } else { + if (out == null) { + out = sub; + } else { + out = new SubHyperplane(out.getHyperplane(), + Region.union(out.getRemainingRegion(), + sub.getRemainingRegion())); + } + } + } + +} diff --git a/src/main/java/org/apache/commons/math/geometry/partitioning/Hyperplane.java b/src/main/java/org/apache/commons/math/geometry/partitioning/Hyperplane.java new file mode 100644 index 000000000..dd8fa775a --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/partitioning/Hyperplane.java @@ -0,0 +1,159 @@ +/* + * 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.math.geometry.partitioning; + +/** This interface represents an hyperplane of a space. + + *The most prominent place where hyperplane appears in space + * partitioning is as cutters. Each partitioning node in a {@link + * BSPTree BSP tree} has a cut {@link SubHyperplane sub-hyperplane} + * which is either an hyperplane or a part of an hyperplane. In an + * n-dimensions euclidean space, an hyperplane is an (n-1)-dimensions + * hyperplane (for example a traditional plane in the 3D euclidean + * space). They can be more exotic objects in specific fields, for + * example a circle on the surface of the unit sphere.
+ + * @version $Revision$ $Date$ + */ +public interface Hyperplane extends SubSpace { + + /** Enumerate for specifying sides of the hyperplane. */ + enum Side { + + /** Code for the plus side of the hyperplane. */ + PLUS, + + /** Code for the minus side of the hyperplane. */ + MINUS, + + /** Code for elements crossing the hyperplane from plus to minus side. */ + BOTH, + + /** Code for the hyperplane itself. */ + HYPER; + + } + + /** Copy the instance. + *The instance created is completely independant of the original + * one. A deep copy is used, none of the underlying objects are + * shared (except for immutable objects).
+ * @return a new hyperplane, copy of the instance + */ + Hyperplane copySelf(); + + /** Get the offset (oriented distance) of a point. + *The offset is 0 if the point is on the underlying hyperplane, + * it is positive if the point is on one particular side of the + * hyperplane, and it is negative if the point is on the other side, + * according to the hyperplane natural orientation.
+ * @param point point to check + * @return offset of the point + */ + double getOffset(Point point); + + /** Check if the instance has the same orientation as another hyperplane. + *This method is expected to be called on parallel hyperplanes + * (i.e. when the {@link #side side} method would return {@link + * Side#HYPER} for some sub-hyperplane having the specified hyperplane + * as its underlying hyperplane). The method should not + * re-check for parallelism, only for orientation, typically by + * testing something like the sign of the dot-products of + * normals.
+ * @param other other hyperplane to check against the instance + * @return true if the instance and the other hyperplane have + * the same orientation + */ + boolean sameOrientationAs(Hyperplane other); + + /** Build the sub-space shared by the instance and another hyperplane. + * @param other other hyperplane + * @return a sub-space at the intersection of the instance and the + * other sub-space (it has a dimension one unit less than the + * instance) + */ + SubSpace intersection(Hyperplane other); + + /** Build a region covering the whole hyperplane. + *The region build is restricted to the sub-space defined by the + * hyperplane. This means that the regions points are consistent + * with the argument of the {@link SubSpace#toSpace toSpace} method + * and with the return value of the {@link SubSpace#toSubSpace + * toSubSpace} method.
+ * @return a region covering the whole hyperplane + */ + Region wholeHyperplane(); + + /** Build a region covering the whole space. + * @return a region containing the instance + */ + Region wholeSpace(); + + /** Compute the relative position of a sub-hyperplane with respect + * to the instance. + * @param sub sub-hyperplane to check + * @return one of {@link Side#PLUS}, {@link Side#MINUS}, {@link Side#BOTH}, + * {@link Side#HYPER} + */ + Side side(SubHyperplane sub); + + /** Split a sub-hyperplane in two parts by the instance. + * @param sub sub-hyperplane to split + * @return an object containing both the part of the sub-hyperplane + * on the plus side of the instance and the part of the + * sub-hyperplane on the minus side of the instance + */ + SplitSubHyperplane split(SubHyperplane sub); + + /** Class holding the results of the {@link Hyperplane#split Hyperplane.split} + * method. */ + class SplitSubHyperplane { + + /** Part of the sub-hyperplane on the plus side of the splitting hyperplane. */ + private final SubHyperplane plus; + + /** Part of the sub-hyperplane on the minus side of the splitting hyperplane. */ + private final SubHyperplane minus; + + /** Build a SplitSubHyperplane from its parts. + * @param plus part of the sub-hyperplane on the plus side of the + * splitting hyperplane + * @param minus part of the sub-hyperplane on the minus side of the + * splitting hyperplane + */ + public SplitSubHyperplane(final SubHyperplane plus, final SubHyperplane minus) { + this.plus = plus; + this.minus = minus; + } + + /** Get the part of the sub-hyperplane on the plus side of the splitting hyperplane. + * @return part of the sub-hyperplane on the plus side of the splitting hyperplane + */ + public SubHyperplane getPlus() { + return plus; + } + + /** Get the part of the sub-hyperplane on the minus side of the splitting hyperplane. + * @return part of the sub-hyperplane on the minus side of the splitting hyperplane + */ + public SubHyperplane getMinus() { + return minus; + } + + } + +} diff --git a/src/main/java/org/apache/commons/math/geometry/partitioning/Point.java b/src/main/java/org/apache/commons/math/geometry/partitioning/Point.java new file mode 100644 index 000000000..ad9d3daea --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/partitioning/Point.java @@ -0,0 +1,29 @@ +/* + * 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.math.geometry.partitioning; + +/** This interface represents a generic point to be used in a space partition. + *
Points are completely virtual entities with no specification at + * all, so this class is essentially a marker interface with no + * methods. This allows to perform partition in traditional euclidean + * n-dimensions spaces, but also in more exotic universes like for + * example the surface of the unit sphere.
+ * @version $Revision$ $Date$ + */ +public interface Point { + // nothing here, this is only a marker interface +} diff --git a/src/main/java/org/apache/commons/math/geometry/partitioning/Region.java b/src/main/java/org/apache/commons/math/geometry/partitioning/Region.java new file mode 100644 index 000000000..b28c42361 --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/partitioning/Region.java @@ -0,0 +1,1069 @@ +/* + * 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.math.geometry.partitioning; + +import java.util.Collection; +import java.util.TreeSet; +import java.util.Comparator; +import java.util.Iterator; +import java.util.ArrayList; + +/** This class represent a region of a space as a partition. + + *Region are subsets of a space, they can be infinite (whole + * space, half space, infinite stripe ...) or finite (polygons in 2D, + * polyhedrons in 3D ...). Their main characteristic is to separate + * points that are considered to be inside the region from + * points considered to be outside of it. In between, there + * may be points on the boundary of the region.
+ + *This implementation is limited to regions for which the boundary + * is composed of several {@link SubHyperplane sub-hyperplanes}, + * including regions with no boundary at all: the whole space and the + * empty region. They are not necessarily finite and not necessarily + * path-connected. They can contain holes.
+ + *Regions can be combined using the traditional sets operations : + * union, intersection, difference and symetric difference (exclusive + * or) for the binary operations, complement for the unary + * operation.
+ + * @version $Revision$ $Date$ + */ +public abstract class Region { + + /** Enumerate for the location of a point with respect to the region. */ + public static enum Location { + /** Code for points inside the partition. */ + INSIDE, + + /** Code for points outside of the partition. */ + OUTSIDE, + + /** Code for points on the partition boundary. */ + BOUNDARY; + } + + /** Inside/Outside BSP tree. */ + private BSPTree tree; + + /** Size of the instance. */ + private double size; + + /** Barycenter. */ + private Point barycenter; + + /** Build a region representing the whole space. + */ + protected Region() { + tree = new BSPTree(Boolean.TRUE); + } + + /** Build a region from an inside/outside BSP tree. + *The leaf nodes of the BSP tree must have a + * {@code Boolean} attribute representing the inside status of + * the corresponding cell (true for inside cells, false for outside + * cells). In order to avoid building too many small objects, it is + * recommended to use the predefined constants + * {@code Boolean.TRUE} and {@code Boolean.FALSE}. The + * tree also must have either null internal nodes or + * internal nodes representing the boundary as specified in the + * {@link #getTree getTree} method).
+ * @param tree inside/outside BSP tree representing the region + */ + protected Region(final BSPTree tree) { + this.tree = tree; + } + + /** Build a Region from a Boundary REPresentation (B-rep). + *The boundary is provided as a collection of {@link + * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the + * interior part of the region on its minus side and the exterior on + * its plus side.
+ *The boundary elements can be in any order, and can form + * several non-connected sets (like for example polygons with holes + * or a set of disjoints polyhedrons considered as a whole). In + * fact, the elements do not even need to be connected together + * (their topological connections are not used here). However, if the + * boundary does not really separate an inside open from an outside + * open (open having here its topological meaning), then subsequent + * calls to the {@link #checkPoint(Point) checkPoint} method will not be + * meaningful anymore.
+ *If the boundary is empty, the region will represent the whole + * space.
+ * @param boundary collection of boundary elements, as a + * collection of {@link SubHyperplane SubHyperplane} objects + */ + protected Region(final CollectionThis method allow to create new instances without knowing + * exactly the type of the region. It is an application of the + * prototype design pattern.
+ *The leaf nodes of the BSP tree must have a + * {@code Boolean} attribute representing the inside status of + * the corresponding cell (true for inside cells, false for outside + * cells). In order to avoid building too many small objects, it is + * recommended to use the predefined constants + * {@code Boolean.TRUE} and {@code Boolean.FALSE}. The + * tree also must have either null internal nodes or + * internal nodes representing the boundary as specified in the + * {@link #getTree getTree} method).
+ * @param newTree inside/outside BSP tree representing the new region + * @return the built region + */ + public abstract Region buildNew(BSPTree newTree); + + /** Recursively build a tree by inserting cut sub-hyperplanes. + * @param node current tree node (it is a leaf node at the beginning + * of the call) + * @param boundary collection of edges belonging to the cell defined + * by the node + */ + private void insertCuts(final BSPTree node, final CollectionThe instance created is completely independant of the original + * one. A deep copy is used, none of the underlying objects are + * shared (except for the underlying tree {@code Boolean} + * attributes and immutable objects).
+ * @return a new region, copy of the instance + */ + public Region copySelf() { + return buildNew(tree.copySelf()); + } + + /** Check if the instance is empty. + * @return true if the instance is empty + */ + public boolean isEmpty() { + return isEmpty(tree); + } + + /** Check if the sub-tree starting at a given node is empty. + * @param node root node of the sub-tree (must have {@link + * Region Region} tree semantics, i.e. the leaf nodes must have + * {@code Boolean} attributes representing an inside/outside + * property) + * @return true if the sub-tree starting at the given node is empty + */ + public static boolean isEmpty(final BSPTree node) { + + // we use a recursive function rather than the BSPTreeVisitor + // interface because we can stop visiting the tree as soon as we + // have found an inside cell + + if (node.getCut() == null) { + // if we find an inside node, the region is not empty + return !isInside(node); + } + + // check both sides of the sub-tree + return isEmpty(node.getMinus()) && isEmpty(node.getPlus()); + + } + + /** Check a leaf node inside attribute. + * @param node leaf node to check + * @return true if the leaf node is an inside node + */ + private static boolean isInside(final BSPTree node) { + return (Boolean) node.getAttribute(); + } + + /** Check if the instance entirely contains another region. + * @param region region to check against the instance + * @return true if the instance contains the specified tree + */ + public boolean contains(final Region region) { + return difference(region, this).isEmpty(); + } + + /** Check a point with respect to the region. + * @param point point to check + * @return a code representing the point status: either {@link + * Location#INSIDE}, {@link Location#OUTSIDE} or {@link Location#BOUNDARY} + */ + public Location checkPoint(final Point point) { + return checkPoint(tree, point); + } + + /** Check a point with respect to the region starting at a given node. + * @param node root node of the region + * @param point point to check + * @return a code representing the point status: either {@link + * Location#INSIDE}, {@link Location#OUTSIDE} or {@link Location#BOUNDARY} + */ + protected Location checkPoint(final BSPTree node, final Point point) { + final BSPTree cell = node.getCell(point); + if (cell.getCut() == null) { + // the point is in the interior of a cell, just check the attribute + return isInside(cell) ? Location.INSIDE : Location.OUTSIDE; + } + + // the point is on a cut-sub-hyperplane, is it on a boundary ? + final Location minusCode = checkPoint(cell.getMinus(), point); + final Location plusCode = checkPoint(cell.getPlus(), point); + return (minusCode == plusCode) ? minusCode : Location.BOUNDARY; + + } + + /** Get the complement of the region (exchanged interior/exterior). + *The instance is not modified, a new region is built.
+ * @return a new region, complement of the instance + */ + public Region getComplement() { + return buildNew(recurseComplement(tree)); + } + + /** Recursively build the complement of a BSP tree. + * @param node current node of the original tree + * @return new tree, complement of the node + */ + private static BSPTree recurseComplement(final BSPTree node) { + if (node.getCut() == null) { + return new BSPTree(isInside(node) ? Boolean.FALSE : Boolean.TRUE); + } + + BoundaryAttribute attribute = (BoundaryAttribute) node.getAttribute(); + if (attribute != null) { + final SubHyperplane plusOutside = + (attribute.plusInside == null) ? null : attribute.plusInside.copySelf(); + final SubHyperplane plusInside = + (attribute.plusOutside == null) ? null : attribute.plusOutside.copySelf(); + attribute = new BoundaryAttribute(plusOutside, plusInside); + } + + return new BSPTree(node.getCut().copySelf(), + recurseComplement(node.getPlus()), + recurseComplement(node.getMinus()), + attribute); + + } + + /** Get the underlying BSP tree. + + *Regions are represented by an underlying inside/outside BSP + * tree whose leaf attributes are {@code Boolean} instances + * representing inside leaf cells if the attribute value is + * {@code true} and outside leaf cells if the attribute is + * {@code false}. These leaf attributes are always present and + * guaranteed to be non null.
+ + *In addition to the leaf attributes, the internal nodes which + * correspond to cells split by cut sub-hyperplanes may contain + * {@link BoundaryAttribute BoundaryAttribute} objects representing + * the parts of the corresponding cut sub-hyperplane that belong to + * the boundary. When the boundary attributes have been computed, + * all internal nodes are guaranteed to have non-null + * attributes, however some {@link BoundaryAttribute + * BoundaryAttribute} instances may have their {@link + * BoundaryAttribute#plusInside plusInside} and {@link + * BoundaryAttribute#plusOutside plusOutside} fields both null if + * the corresponding cut sub-hyperplane does not have any parts + * belonging to the boundary.
+ + *Since computing the boundary is not always required and can be + * time-consuming for large trees, these internal nodes attributes + * are computed using lazy evaluation only when required by setting + * the {@code includeBoundaryAttributes} argument to + * {@code true}. Once computed, these attributes remain in the + * tree, which implies that in this case, further calls to the + * method for the same region will always include these attributes + * regardless of the value of the + * {@code includeBoundaryAttributes} argument.
+ + * @param includeBoundaryAttributes if true, the boundary attributes + * at internal nodes are guaranteed to be included (they may be + * included even if the argument is false, if they have already been + * computed due to a previous call) + * @return underlying BSP tree + * @see BoundaryAttribute + */ + public BSPTree getTree(final boolean includeBoundaryAttributes) { + if (includeBoundaryAttributes && (tree.getCut() != null) && (tree.getAttribute() == null)) { + // we need to compute the boundary attributes + recurseBuildBoundary(tree); + } + return tree; + } + + /** Class holding boundary attributes. + *This class is used for the attributes associated with the + * nodes of region boundary shell trees returned by the {@link + * Region#getTree Region.getTree}. It contains the + * parts of the node cut sub-hyperplane that belong to the + * boundary.
+ *This class is a simple placeholder, it does not provide any + * processing methods.
+ * @see Region#getTree + */ + public static class BoundaryAttribute { + + /** Part of the node cut sub-hyperplane that belongs to the + * boundary and has the outside of the region on the plus side of + * its underlying hyperplane (may be null). + */ + private final SubHyperplane plusOutside; + + /** Part of the node cut sub-hyperplane that belongs to the + * boundary and has the inside of the region on the plus side of + * its underlying hyperplane (may be null). + */ + private final SubHyperplane plusInside; + + /** Simple constructor. + * @param plusOutside part of the node cut sub-hyperplane that + * belongs to the boundary and has the outside of the region on + * the plus side of its underlying hyperplane (may be null) + * @param plusInside part of the node cut sub-hyperplane that + * belongs to the boundary and has the inside of the region on the + * plus side of its underlying hyperplane (may be null) + */ + public BoundaryAttribute(final SubHyperplane plusOutside, + final SubHyperplane plusInside) { + this.plusOutside = plusOutside; + this.plusInside = plusInside; + } + + /** Get the part of the node cut sub-hyperplane that belongs to the + * boundary and has the outside of the region on the plus side of + * its underlying hyperplane. + * @return part of the node cut sub-hyperplane that belongs to the + * boundary and has the outside of the region on the plus side of + * its underlying hyperplane + */ + public SubHyperplane getPlusOutside() { + return plusOutside; + } + + /** Get the part of the node cut sub-hyperplane that belongs to the + * boundary and has the inside of the region on the plus side of + * its underlying hyperplane. + * @return part of the node cut sub-hyperplane that belongs to the + * boundary and has the inside of the region on the plus side of + * its underlying hyperplane + */ + public SubHyperplane getPlusInside() { + return plusInside; + } + + + } + + /** Recursively build the boundary shell tree. + * @param node current node in the inout tree + */ + private void recurseBuildBoundary(final BSPTree node) { + if (node.getCut() != null) { + + SubHyperplane plusOutside = null; + SubHyperplane plusInside = null; + + // characterize the cut sub-hyperplane, + // first with respect to the plus sub-tree + final Characterization plusChar = new Characterization(); + characterize(node.getPlus(), node.getCut().copySelf(), plusChar); + + if (plusChar.hasOut()) { + // plusChar.out corresponds to a subset of the cut + // sub-hyperplane known to have outside cells on its plus + // side, we want to check if parts of this subset do have + // inside cells on their minus side + final Characterization minusChar = new Characterization(); + characterize(node.getMinus(), plusChar.getOut(), minusChar); + if (minusChar.hasIn()) { + plusOutside = minusChar.getIn(); + } + } + + if (plusChar.hasIn()) { + // plusChar.in corresponds to a subset of the cut + // sub-hyperplane known to have inside cells on its plus + // side, we want to check if parts of this subset do have + // outside cells on their minus side + final Characterization minusChar = new Characterization(); + characterize(node.getMinus(), plusChar.getIn(), minusChar); + if (minusChar.hasOut()) { + plusInside = minusChar.getOut(); + } + } + + node.setAttribute(new BoundaryAttribute(plusOutside, plusInside)); + recurseBuildBoundary(node.getPlus()); + recurseBuildBoundary(node.getMinus()); + + } + } + + /** Filter the parts of an hyperplane belonging to the boundary. + *The filtering consist in splitting the specified + * sub-hyperplane into several parts lying in inside and outside + * cells of the tree. The principle is to call this method twice for + * each cut sub-hyperplane in the tree, once one the plus node and + * once on the minus node. The parts that have the same flag + * (inside/inside or outside/outside) do not belong to the boundary + * while parts that have different flags (inside/outside or + * outside/inside) do belong to the boundary.
+ * @param node current BSP tree node + * @param sub sub-hyperplane to characterize + * @param characterization placeholder where to put the characterized parts + */ + private static void characterize(final BSPTree node, final SubHyperplane sub, + final Characterization characterization) { + if (node.getCut() == null) { + // we have reached a leaf node + final boolean inside = (Boolean) node.getAttribute(); + characterization.add(sub, inside); + } else { + final Hyperplane hyperplane = node.getCut().getHyperplane(); + switch (hyperplane.side(sub)) { + case PLUS: + characterize(node.getPlus(), sub, characterization); + break; + case MINUS: + characterize(node.getMinus(), sub, characterization); + break; + case BOTH: + final Hyperplane.SplitSubHyperplane split = hyperplane.split(sub); + characterize(node.getPlus(), split.getPlus(), characterization); + characterize(node.getMinus(), split.getMinus(), characterization); + break; + default: + // this should not happen + throw new RuntimeException("internal error"); + } + } + } + + /** Get the size of the boundary. + * @return the size of the boundary (this is 0 in 1D, a length in + * 2D, an area in 3D ...) + */ + public double getBoundarySize() { + final BoundarySizeVisitor visitor = new BoundarySizeVisitor(); + getTree(true).visit(visitor); + return visitor.getSize(); + } + + /** Visitor computing the boundary size. */ + private static class BoundarySizeVisitor implements BSPTreeVisitor { + + /** Size of the boundary. */ + private double boundarySize; + + /** Simple constructor. + */ + public BoundarySizeVisitor() { + boundarySize = 0; + } + + /** {@inheritDoc}*/ + public Order visitOrder(final BSPTree node) { + return Order.MINUS_SUB_PLUS; + } + + /** {@inheritDoc}*/ + public void visitInternalNode(final BSPTree node) { + final BoundaryAttribute attribute = (BoundaryAttribute) node.getAttribute(); + if (attribute.plusOutside != null) { + boundarySize += attribute.plusOutside.getRemainingRegion().getSize(); + } + if (attribute.plusInside != null) { + boundarySize += attribute.plusInside.getRemainingRegion().getSize(); + } + } + + /** {@inheritDoc}*/ + public void visitLeafNode(final BSPTree node) { + } + + /** Get the size of the boundary. + * @return size of the boundary + */ + public double getSize() { + return boundarySize; + } + + } + + /** Get the size of the instance. + * @return the size of the instance (this is a length in 1D, an area + * in 2D, a volume in 3D ...) + */ + public double getSize() { + if (barycenter == null) { + computeGeometricalProperties(); + } + return size; + } + + /** Set the size of the instance. + * @param size size of the instance + */ + protected void setSize(final double size) { + this.size = size; + } + + /** Get the barycenter of the instance. + * @return an object representing the barycenter + */ + public Point getBarycenter() { + if (barycenter == null) { + computeGeometricalProperties(); + } + return barycenter; + } + + /** Set the barycenter of the instance. + * @param barycenter barycenter of the instance + */ + protected void setBarycenter(final Point barycenter) { + this.barycenter = barycenter; + } + + /** Compute some geometrical properties. + *The properties to compute are the barycenter and the size.
+ */ + protected abstract void computeGeometricalProperties(); + + /** Transform a region. + *Applying a transform to a region consist in applying the + * transform to all the hyperplanes of the underlying BSP tree and + * of the boundary (and also to the sub-hyperplanes embedded in + * these hyperplanes) and to the barycenter. The instance is not + * modified, a new instance is built.
+ * @param transform transform to apply + * @return a new region, resulting from the application of the + * transform to the instance + */ + public Region applyTransform(final Transform transform) { + + // transform the BSP tree + final Region tRegion = buildNew(recurseTransform(tree, transform)); + + // transform the barycenter + if (barycenter != null) { + tRegion.size = size; + tRegion.barycenter = transform.apply(barycenter); + } + + return tRegion; + + } + + /** Recursively transform an inside/outside BSP-tree. + * @param node current BSP tree node + * @param transform transform to apply + * @return a new tree + */ + private BSPTree recurseTransform(final BSPTree node, final Transform transform) { + + if (node.getCut() == null) { + return new BSPTree(node.getAttribute()); + } + + final SubHyperplane sub = node.getCut(); + final SubHyperplane tSub = sub.applyTransform(transform); + BoundaryAttribute attribute = (BoundaryAttribute) node.getAttribute(); + if (attribute != null) { + final SubHyperplane tPO = + (attribute.getPlusOutside() == null) ? null : attribute.getPlusOutside().applyTransform(transform); + final SubHyperplane tPI = + (attribute.getPlusInside() == null) ? null : attribute.getPlusInside().applyTransform(transform); + attribute = new BoundaryAttribute(tPO, tPI); + } + + return new BSPTree(tSub, + recurseTransform(node.getPlus(), transform), + recurseTransform(node.getMinus(), transform), + attribute); + + } + + /** Compute the relative position of the instance with respect to an + * hyperplane. + * @param hyperplane reference hyperplane + * @return one of {@link Hyperplane.Side#PLUS Hyperplane.Side.PLUS}, {@link + * Hyperplane.Side#MINUS Hyperplane.Side.MINUS}, {@link Hyperplane.Side#BOTH + * Hyperplane.Side.BOTH} or {@link Hyperplane.Side#HYPER Hyperplane.Side.HYPER} + * (the latter result can occur only if the tree contains only one + * cut hyperplane) + */ + public Hyperplane.Side side(final Hyperplane hyperplane) { + final Sides sides = new Sides(); + recurseSides(tree, new SubHyperplane(hyperplane), sides); + return sides.plusFound() ? + (sides.minusFound() ? Hyperplane.Side.BOTH : Hyperplane.Side.PLUS) : + (sides.minusFound() ? Hyperplane.Side.MINUS : Hyperplane.Side.HYPER); + } + + /** Search recursively for inside leaf nodes on each side of the given hyperplane. + + *The algorithm used here is directly derived from the one + * described in section III (Binary Partitioning of a BSP + * Tree) of the Bruce Naylor, John Amanatides and William + * Thibault paper Merging + * BSP Trees Yields Polyhedral Set Operations Proc. Siggraph + * '90, Computer Graphics 24(4), August 1990, pp 115-124, published + * by the Association for Computing Machinery (ACM)..
+ + * @param node current BSP tree node + * @param sub sub-hyperplane + * @param sides object holding the sides found + */ + private void recurseSides(final BSPTree node, final SubHyperplane sub, final Sides sides) { + + if (node.getCut() == null) { + if (isInside(node)) { + // this is an inside cell expanding across the hyperplane + sides.rememberPlusFound(); + sides.rememberMinusFound(); + } + return; + } + + final Hyperplane hyperplane = node.getCut().getHyperplane(); + switch (hyperplane.side(sub)) { + case PLUS : + // the sub-hyperplane is entirely in the plus sub-tree + if (sub.getHyperplane().side(node.getCut()) == Hyperplane.Side.PLUS) { + if (!isEmpty(node.getMinus())) { + sides.rememberPlusFound(); + } + } else { + if (!isEmpty(node.getMinus())) { + sides.rememberMinusFound(); + } + } + if (!(sides.plusFound() && sides.minusFound())) { + recurseSides(node.getPlus(), sub, sides); + } + break; + case MINUS : + // the sub-hyperplane is entirely in the minus sub-tree + if (sub.getHyperplane().side(node.getCut()) == Hyperplane.Side.PLUS) { + if (!isEmpty(node.getPlus())) { + sides.rememberPlusFound(); + } + } else { + if (!isEmpty(node.getPlus())) { + sides.rememberMinusFound(); + } + } + if (!(sides.plusFound() && sides.minusFound())) { + recurseSides(node.getMinus(), sub, sides); + } + break; + case BOTH : + // the sub-hyperplane extends in both sub-trees + final Hyperplane.SplitSubHyperplane split = hyperplane.split(sub); + + // explore first the plus sub-tree + recurseSides(node.getPlus(), split.getPlus(), sides); + + // if needed, explore the minus sub-tree + if (!(sides.plusFound() && sides.minusFound())) { + recurseSides(node.getMinus(), split.getMinus(), sides); + } + break; + default : + // the sub-hyperplane and the cut sub-hyperplane share the same hyperplane + if (node.getCut().getHyperplane().sameOrientationAs(sub.getHyperplane())) { + if ((node.getPlus().getCut() != null) || isInside(node.getPlus())) { + sides.rememberPlusFound(); + } + if ((node.getMinus().getCut() != null) || isInside(node.getMinus())) { + sides.rememberMinusFound(); + } + } else { + if ((node.getPlus().getCut() != null) || isInside(node.getPlus())) { + sides.rememberMinusFound(); + } + if ((node.getMinus().getCut() != null) || isInside(node.getMinus())) { + sides.rememberPlusFound(); + } + } + } + + } + + /** Utility class holding the already found sides. */ + private static final class Sides { + + /** Indicator of inside leaf nodes found on the plus side. */ + private boolean plusFound; + + /** Indicator of inside leaf nodes found on the plus side. */ + private boolean minusFound; + + /** Simple constructor. + */ + public Sides() { + plusFound = false; + minusFound = false; + } + + /** Remember the fact that inside leaf nodes have been found on the plus side. + */ + public void rememberPlusFound() { + plusFound = true; + } + + /** Check if inside leaf nodes have been found on the plus side. + * @return true if inside leaf nodes have been found on the plus side + */ + public boolean plusFound() { + return plusFound; + } + + /** Remember the fact that inside leaf nodes have been found on the minus side. + */ + public void rememberMinusFound() { + minusFound = true; + } + + /** Check if inside leaf nodes have been found on the minus side. + * @return true if inside leaf nodes have been found on the minus side + */ + public boolean minusFound() { + return minusFound; + } + + } + + /** Get the parts of a sub-hyperplane that are contained in the region. + *The parts of the sub-hyperplane that belong to the boundary are + * not included in the resulting parts.
+ * @param sub sub-hyperplane traversing the region + * @return filtered sub-hyperplane + */ + public SubHyperplane intersection(final SubHyperplane sub) { + return recurseIntersection(tree, sub); + } + + /** Recursively compute the parts of a sub-hyperplane that are + * contained in the region. + * @param node current BSP tree node + * @param sub sub-hyperplane traversing the region + * @return filtered sub-hyperplane + */ + private SubHyperplane recurseIntersection(final BSPTree node, final SubHyperplane sub) { + + if (node.getCut() == null) { + return isInside(node) ? sub.copySelf() : null; + } + + final Hyperplane hyperplane = node.getCut().getHyperplane(); + switch (hyperplane.side(sub)) { + case PLUS : + return recurseIntersection(node.getPlus(), sub); + case MINUS : + return recurseIntersection(node.getMinus(), sub); + case BOTH : + final Hyperplane.SplitSubHyperplane split = hyperplane.split(sub); + final SubHyperplane plus = recurseIntersection(node.getPlus(), split.getPlus()); + final SubHyperplane minus = recurseIntersection(node.getMinus(), split.getMinus()); + if (plus == null) { + return minus; + } else if (minus == null) { + return plus; + } else { + return new SubHyperplane(plus.getHyperplane(), + Region.union(plus.getRemainingRegion(), + minus.getRemainingRegion())); + } + default : + return recurseIntersection(node.getPlus(), + recurseIntersection(node.getMinus(), sub)); + } + + } + + /** Compute the union of two regions. + * @param region1 first region (will be unusable after the operation as + * parts of it will be reused in the new region) + * @param region2 second region (will be unusable after the operation as + * parts of it will be reused in the new region) + * @return a new region, result of {@code region1 union region2} + */ + public static Region union(final Region region1, final Region region2) { + final BSPTree tree = region1.tree.merge(region2.tree, new UnionMerger()); + tree.visit(new InternalNodesCleaner()); + return region1.buildNew(tree); + } + + /** Compute the intersection of two regions. + * @param region1 first region (will be unusable after the operation as + * parts of it will be reused in the new region) + * @param region2 second region (will be unusable after the operation as + * parts of it will be reused in the new region) + * @return a new region, result of {@code region1 intersection region2} + */ + public static Region intersection(final Region region1, final Region region2) { + final BSPTree tree = region1.tree.merge(region2.tree, new IntersectionMerger()); + tree.visit(new InternalNodesCleaner()); + return region1.buildNew(tree); + } + + /** Compute the symmetric difference (exclusive or) of two regions. + * @param region1 first region (will be unusable after the operation as + * parts of it will be reused in the new region) + * @param region2 second region (will be unusable after the operation as + * parts of it will be reused in the new region) + * @return a new region, result of {@code region1 xor region2} + */ + public static Region xor(final Region region1, final Region region2) { + final BSPTree tree = region1.tree.merge(region2.tree, new XORMerger()); + tree.visit(new InternalNodesCleaner()); + return region1.buildNew(tree); + } + + /** Compute the difference of two regions. + * @param region1 first region (will be unusable after the operation as + * parts of it will be reused in the new region) + * @param region2 second region (will be unusable after the operation as + * parts of it will be reused in the new region) + * @return a new region, result of {@code region1 minus region2} + */ + public static Region difference(final Region region1, final Region region2) { + final BSPTree tree = region1.tree.merge(region2.tree, new DifferenceMerger()); + tree.visit(new InternalNodesCleaner()); + return region1.buildNew(tree); + } + + /** Leaf node / tree merger for union operation. */ + private static final class UnionMerger implements BSPTree.LeafMerger { + /** {@inheritDoc} */ + public BSPTree merge(final BSPTree leaf, final BSPTree tree, + final BSPTree parentTree, final boolean isPlusChild, + final boolean leafFromInstance) { + if (isInside(leaf)) { + // the leaf node represents an inside cell + leaf.insertInTree(parentTree, isPlusChild); + return leaf; + } + // the leaf node represents an outside cell + tree.insertInTree(parentTree, isPlusChild); + return tree; + } + }; + + /** Leaf node / tree merger for intersection operation. */ + private static final class IntersectionMerger implements BSPTree.LeafMerger { + /** {@inheritDoc} */ + public BSPTree merge(final BSPTree leaf, final BSPTree tree, + final BSPTree parentTree, final boolean isPlusChild, + final boolean leafFromInstance) { + if (isInside(leaf)) { + // the leaf node represents an inside cell + tree.insertInTree(parentTree, isPlusChild); + return tree; + } + // the leaf node represents an outside cell + leaf.insertInTree(parentTree, isPlusChild); + return leaf; + } + }; + + /** Leaf node / tree merger for xor operation. */ + private static final class XORMerger implements BSPTree.LeafMerger { + /** {@inheritDoc} */ + public BSPTree merge(final BSPTree leaf, final BSPTree tree, + final BSPTree parentTree, final boolean isPlusChild, + final boolean leafFromInstance) { + BSPTree t = tree; + if (isInside(leaf)) { + // the leaf node represents an inside cell + t = recurseComplement(t); + } + t.insertInTree(parentTree, isPlusChild); + return t; + } + }; + + /** Leaf node / tree merger for difference operation. + *The algorithm used here is directly derived from the one + * described in section III (Binary Partitioning of a BSP + * Tree) of the Naylor, Amanatides and Thibault paper. An error + * was detected and corrected in the figure 5.1 of the article for + * merging leaf nodes with complete trees. Contrary to what is said + * in the figure, the {@code ELSE} part of if is not the same + * as the first part with {@code T1} and {@codeT2} + * swapped. {@code T1} and {@codeT2} must be swapped + * everywhere except in the {@code RETURN} part of the + * {@code DIFFERENCE} operation: if {@codeT2} is an + * in-cell, we must return {@code Complement_Bspt(T2)}, not + * {@code Complement_Bspt(T1)}, and if {@codeT2} is an + * out-cell, we must return {@code T1}, not {@codeT2}
+ */ + private static final class DifferenceMerger implements BSPTree.LeafMerger { + /** {@inheritDoc} */ + public BSPTree merge(final BSPTree leaf, final BSPTree tree, + final BSPTree parentTree, final boolean isPlusChild, + final boolean leafFromInstance) { + if (isInside(leaf)) { + // the leaf node represents an inside cell + final BSPTree argTree = recurseComplement(leafFromInstance ? tree : leaf); + argTree.insertInTree(parentTree, isPlusChild); + return argTree; + } + // the leaf node represents an outside cell + final BSPTree instanceTree = leafFromInstance ? leaf : tree; + instanceTree.insertInTree(parentTree, isPlusChild); + return instanceTree; + } + }; + + /** Visitor removing internal nodes attributes. */ + private static final class InternalNodesCleaner implements BSPTreeVisitor { + + /** {@inheritDoc} */ + public Order visitOrder(final BSPTree node) { + return Order.PLUS_SUB_MINUS; + } + + /** {@inheritDoc} */ + public void visitInternalNode(final BSPTree node) { + node.setAttribute(null); + } + + /** {@inheritDoc} */ + public void visitLeafNode(final BSPTree node) { + } + + } + +} diff --git a/src/main/java/org/apache/commons/math/geometry/partitioning/SubHyperplane.java b/src/main/java/org/apache/commons/math/geometry/partitioning/SubHyperplane.java new file mode 100644 index 000000000..759eeaa94 --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/partitioning/SubHyperplane.java @@ -0,0 +1,138 @@ +/* + * 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.math.geometry.partitioning; + +/** This interface represents the remaining parts of an hyperplane after + * other parts have been chopped off. + + *sub-hyperplanes are obtained when parts of an {@link + * Hyperplane hyperplane} are chopped off by other hyperplanes that + * intersect it. The remaining part is a convex region. Such objects + * appear in {@link BSPTree BSP trees} as the intersection of a cut + * hyperplane with the convex region which it splits, the chopping + * hyperplanes are the cut hyperplanes closer to the tree root.
+ + * @version $Revision$ $Date$ + */ +public class SubHyperplane { + + /** Underlying hyperplane. */ + private final Hyperplane hyperplane; + + /** Remaining region of the hyperplane. */ + private final Region remainingRegion; + + /** Build a chopped hyperplane that is not chopped at all. + * @param hyperplane underlying hyperplane + */ + public SubHyperplane(final Hyperplane hyperplane) { + this.hyperplane = hyperplane; + remainingRegion = hyperplane.wholeHyperplane(); + } + + /** Build a sub-hyperplane from an hyperplane and a region. + * @param hyperplane underlying hyperplane + * @param remainingRegion remaining region of the hyperplane + */ + public SubHyperplane(final Hyperplane hyperplane, final Region remainingRegion) { + this.hyperplane = hyperplane; + this.remainingRegion = remainingRegion; + } + + /** Copy the instance. + *The instance created is completely independant of the original + * one. A deep copy is used, none of the underlying objects are + * shared (except for the nodes attributes and immutable + * objects).
+ * @return a new sub-hyperplane, copy of the instance + */ + public SubHyperplane copySelf() { + return new SubHyperplane(hyperplane.copySelf(), remainingRegion.copySelf()); + } + + /** Get the underlying hyperplane. + * @return underlying hyperplane + */ + public Hyperplane getHyperplane() { + return hyperplane; + } + + /** Get the remaining region of the hyperplane. + *The returned region is expressed in the canonical hyperplane + * frame and has the hyperplane dimension. For example a chopped + * hyperplane in the 3D euclidean is a 2D plane and the + * corresponding region is a convex 2D polygon.
+ * @return remaining region of the hyperplane + */ + public Region getRemainingRegion() { + return remainingRegion; + } + + /** Apply a transform to the instance. + *The instance must be a (D-1)-dimension sub-hyperplane with + * respect to the transform not a (D-2)-dimension + * sub-hyperplane the transform knows how to transform by + * itself. The transform will consist in transforming first the + * hyperplane and then the all region using the various methods + * provided by the transform.
+ * @param transform D-dimension transform to apply + * @return the transformed instance + */ + public SubHyperplane applyTransform(final Transform transform) { + final Hyperplane tHyperplane = transform.apply(hyperplane); + final BSPTree tTree = + recurseTransform(remainingRegion.getTree(false), tHyperplane, transform); + return new SubHyperplane(tHyperplane, remainingRegion.buildNew(tTree)); + } + + /** Recursively transform a BSP-tree from a sub-hyperplane. + * @param node current BSP tree node + * @param transformed image of the instance hyperplane by the transform + * @param transform transform to apply + * @return a new tree + */ + private BSPTree recurseTransform(final BSPTree node, final Hyperplane transformed, + final Transform transform) { + if (node.getCut() == null) { + return new BSPTree(node.getAttribute()); + } + + Region.BoundaryAttribute attribute = + (Region.BoundaryAttribute) node.getAttribute(); + if (attribute != null) { + final SubHyperplane tPO = (attribute.getPlusOutside() == null) ? + null : + transform.apply(attribute.getPlusOutside(), + hyperplane, transformed); + final SubHyperplane tPI = (attribute.getPlusInside() == null) ? + null : + transform.apply(attribute.getPlusInside(), + hyperplane, transformed); + attribute = new Region.BoundaryAttribute(tPO, tPI); + } + + return new BSPTree(transform.apply(node.getCut(), + hyperplane, transformed), + recurseTransform(node.getPlus(), transformed, + transform), + recurseTransform(node.getMinus(), transformed, + transform), + attribute); + + } + +} diff --git a/src/main/java/org/apache/commons/math/geometry/partitioning/SubSpace.java b/src/main/java/org/apache/commons/math/geometry/partitioning/SubSpace.java new file mode 100644 index 000000000..148a24329 --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/partitioning/SubSpace.java @@ -0,0 +1,50 @@ +/* + * 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.math.geometry.partitioning; + + +/** This interface represents a sub-space of a space. + + *Sub-spaces are the lower dimensions subsets of a n-dimensions + * space. The (n-1)-dimension sub-spaces are specific sub-spaces known + * as {@link Hyperplane hyperplanes}.
+ + *In the 3D euclidean space, hyperplanes are 2D planes, and the 1D + * sub-spaces are lines.
+ + * @see Hyperplane + * @version $Revision$ $Date$ + */ +public interface SubSpace { + + /** Transform a space point into a sub-space point. + * @param point n-dimension point of the space + * @return (n-1)-dimension point of the sub-space corresponding to + * the specified space point + * @see #toSpace + */ + Point toSubSpace(Point point); + + /** Transform a sub-space point into a space point. + * @param point (n-1)-dimension point of the sub-space + * @return n-dimension point of the space corresponding to the + * specified sub-space point + * @see #toSubSpace + */ + Point toSpace(Point point); + +} diff --git a/src/main/java/org/apache/commons/math/geometry/partitioning/Transform.java b/src/main/java/org/apache/commons/math/geometry/partitioning/Transform.java new file mode 100644 index 000000000..a55dab567 --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/partitioning/Transform.java @@ -0,0 +1,74 @@ +/* + * 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.math.geometry.partitioning; + + +/** This interface represents an inversible affine transform in a space. + *Inversible affine transform include for example scalings, + * translations, rotations.
+ + *Transforms are dimension-specific. The consistency rules between + * the three {@code apply} methods are the following ones for a + * transformed defined for dimension D:
+ *-
+ *
- + * the transform can be applied to a point in the + * D-dimension space using its {@link #apply(Point)} + * method + * + *
- + * the transform can be applied to a (D-1)-dimension + * hyperplane in the D-dimension space using its + * {@link #apply(Hyperplane)} method + * + *
- + * the transform can be applied to a (D-2)-dimension + * sub-hyperplane in a (D-1)-dimension hyperplane using + * its {@link #apply(SubHyperplane, Hyperplane, Hyperplane)} + * method + * + *
+{@link org.apache.commons.math.geometry.partitioning.BSPTree BSP trees} +are an efficient way to represent parts of space and in particular +polytopes (line segments in 1D, polygons in 2D and polyhedrons in 3D) +and to operate on them. The main principle is to recursively subdivide +the space using simple hyperplanes (points in 1D, lines in 2D, planes +in 3D). +
+ ++We start with a tree composed of a single node without any cut +hyperplane: it represents the complete space, which is a convex +part. If we add a cut hyperplane to this node, this represents a +partition with the hyperplane at the node level and two half spaces at +each side of the cut hyperplane. These half-spaces are represented by +two child nodes without any cut hyperplanes associated, the plus child +which represents the half space on the plus side of the cut hyperplane +and the minus child on the other side. Continuing the subdivisions, we +end up with a tree having internal nodes that are associated with a +cut hyperplane and leaf nodes without any hyperplane which correspond +to convex parts. +
+ ++When BSP trees are used to represent polytopes, the convex parts are +known to be completely inside or outside the polytope as long as there +is no facet in the part (which is obviously the case if the cut +hyperplanes have been chosen as the underlying hyperplanes of the +facets (this is called an autopartition) and if the subdivision +process has been continued until all facets have been processed. It is +important to note that the polytope is not defined by a +single part, but by several convex ones. This is the property that +allows BSP-trees to represent non-convex polytopes despites all parts +are convex. The {@link +org.apache.commons.math.geometry.partitioning.Region Region} class is +devoted to this representation, it is build on top of the {@link +org.apache.commons.math.geometry.partitioning.BSPTree BSPTree} class using +boolean objects as the leaf nodes attributes to represent the +inside/outside property of each leaf part, and also adds various +methods dealing with boundaries (i.e. the separation between the +inside and the outside parts). +
+ ++Rather than simply associating the internal nodes with an hyperplane, +we consider sub-hyperplanes which correspond to the part of +the hyperplane that is inside the convex part defined by all the +parent nodes (this implies that the sub-hyperplane at root node is in +fact a complete hyperplane, because there is no parent to bound +it). Since the parts are convex, the sub-hyperplanes are convex, in +3D the convex parts are convex polyhedrons, and the sub-hyperplanes +are convex polygons that cut these polyhedrons in two +sub-polyhedrons. Using this definition, a BSP tree completely +partitions the space. Each point either belongs to one of the +sub-hyperplanes in an internal node or belongs to one of the leaf +convex parts. +
+ +
+In order to determine where a point is, it is sufficient to check its
+position with respect to the root cut hyperplane, to select the
+corresponding child tree and to repeat the procedure recursively,
+until either the point appears to be exactly on one of the hyperplanes
+in the middle of the tree or to be in one of the leaf parts. For
+this operation, it is sufficient to consider the complete hyperplanes,
+there is no need to check the points with the boundary of the
+sub-hyperplanes, because this check has in fact already been realized
+by the recursive descent in the tree. This is very easy to do and very
+efficient, especially if the tree is well balanced (the cost is
+O(log(n)) where n is the number of facets)
+or if the first tree levels close to the root discriminate large parts
+of the total space.
+
+One of the main sources for the development of this package was Bruce +Naylor, John Amanatides and William Thibault paper Merging +BSP Trees Yields Polyhedral Set Operations Proc. Siggraph '90, +Computer Graphics 24(4), August 1990, pp 115-124, published by the +Association for Computing Machinery (ACM). The same paper can also be +found here. +
+ + + diff --git a/src/main/java/org/apache/commons/math/geometry/partitioning/utilities/AVLTree.java b/src/main/java/org/apache/commons/math/geometry/partitioning/utilities/AVLTree.java new file mode 100644 index 000000000..e1e5dc4f7 --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/partitioning/utilities/AVLTree.java @@ -0,0 +1,631 @@ +/* + * 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.math.geometry.partitioning.utilities; + +/** This class implements AVL trees. + + *The purpose of this class is to sort elements while allowing + * duplicate elements (i.e. such that {@code a.equals(b)} is + * true). The {@code SortedSet} interface does not allow this, so + * a specific class is needed. Null elements are not allowed.
+ + *Since the {@code equals} method is not sufficient to + * differentiate elements, the {@link #delete delete} method is + * implemented using the equality operator.
+ + *In order to clearly mark the methods provided here do not have + * the same semantics as the ones specified in the + * {@code SortedSet} interface, different names are used + * ({@code add} has been replaced by {@link #insert insert} and + * {@code remove} has been replaced by {@link #delete + * delete}).
+ + *This class is based on the C implementation Georg Kraml has put + * in the public domain. Unfortunately, his page seems not + * to exist any more.
+ + * @paramThe element is deleted only if there is a node {@code n} + * containing exactly the element instance specified, i.e. for which + * {@code n.getElement() == element}. This is purposely + * different from the specification of the + * {@code java.util.Set} {@code remove} method (in fact, + * this is the reason why a specific class has been developed).
+ * @param element element to delete (silently ignored if null) + * @return true if the element was deleted from the tree + */ + public boolean delete(final T element) { + if (element != null) { + for (Node node = getNotSmaller(element); node != null; node = node.getNext()) { + // loop over all elements neither smaller nor larger + // than the specified one + if (node.element == element) { + node.delete(); + return true; + } else if (node.element.compareTo(element) > 0) { + // all the remaining elements are known to be larger, + // the element is not in the tree + return false; + } + } + } + return false; + } + + /** Check if the tree is empty. + * @return true if the tree is empty + */ + public boolean isEmpty() { + return top == null; + } + + + /** Get the number of elements of the tree. + * @return number of elements contained in the tree + */ + public int size() { + return (top == null) ? 0 : top.size(); + } + + /** Get the node whose element is the smallest one in the tree. + * @return the tree node containing the smallest element in the tree + * or null if the tree is empty + * @see #getLargest + * @see #getNotSmaller + * @see #getNotLarger + * @see Node#getPrevious + * @see Node#getNext + */ + public Node getSmallest() { + return (top == null) ? null : top.getSmallest(); + } + + /** Get the node whose element is the largest one in the tree. + * @return the tree node containing the largest element in the tree + * or null if the tree is empty + * @see #getSmallest + * @see #getNotSmaller + * @see #getNotLarger + * @see Node#getPrevious + * @see Node#getNext + */ + public Node getLargest() { + return (top == null) ? null : top.getLargest(); + } + + /** Get the node whose element is not smaller than the reference object. + * @param reference reference object (may not be in the tree) + * @return the tree node containing the smallest element not smaller + * than the reference object or null if either the tree is empty or + * all its elements are smaller than the reference object + * @see #getSmallest + * @see #getLargest + * @see #getNotLarger + * @see Node#getPrevious + * @see Node#getNext + */ + public Node getNotSmaller(final T reference) { + Node candidate = null; + for (Node node = top; node != null;) { + if (node.element.compareTo(reference) < 0) { + if (node.right == null) { + return candidate; + } + node = node.right; + } else { + candidate = node; + if (node.left == null) { + return candidate; + } + node = node.left; + } + } + return null; + } + + /** Get the node whose element is not larger than the reference object. + * @param reference reference object (may not be in the tree) + * @return the tree node containing the largest element not larger + * than the reference object (in which case the node is guaranteed + * not to be empty) or null if either the tree is empty or all its + * elements are larger than the reference object + * @see #getSmallest + * @see #getLargest + * @see #getNotSmaller + * @see Node#getPrevious + * @see Node#getNext + */ + public Node getNotLarger(final T reference) { + Node candidate = null; + for (Node node = top; node != null;) { + if (node.element.compareTo(reference) > 0) { + if (node.left == null) { + return candidate; + } + node = node.left; + } else { + candidate = node; + if (node.right == null) { + return candidate; + } + node = node.right; + } + } + return null; + } + + /** Enum for tree skew factor. */ + private static enum Skew { + /** Code for left high trees. */ + LEFT_HIGH, + + /** Code for right high trees. */ + RIGHT_HIGH, + + /** Code for Skew.BALANCED trees. */ + BALANCED; + } + + /** This class implements AVL trees nodes. + *AVL tree nodes implement all the logical structure of the + * tree. Nodes are created by the {@link AVLTree AVLTree} class.
+ *The nodes are not independant from each other but must obey + * specific balancing constraints and the tree structure is + * rearranged as elements are inserted or deleted from the tree. The + * creation, modification and tree-related navigation methods have + * therefore restricted access. Only the order-related navigation, + * reading and delete methods are public.
+ * @see AVLTree + */ + public class Node { + + /** Element contained in the current node. */ + private T element; + + /** Left sub-tree. */ + private Node left; + + /** Right sub-tree. */ + private Node right; + + /** Parent tree. */ + private Node parent; + + /** Skew factor. */ + private Skew skew; + + /** Build a node for a specified element. + * @param element element + * @param parent parent node + */ + Node(final T element, final Node parent) { + this.element = element; + left = null; + right = null; + this.parent = parent; + skew = Skew.BALANCED; + } + + /** Get the contained element. + * @return element contained in the node + */ + public T getElement() { + return element; + } + + /** Get the number of elements of the tree rooted at this node. + * @return number of elements contained in the tree rooted at this node + */ + int size() { + return 1 + ((left == null) ? 0 : left.size()) + ((right == null) ? 0 : right.size()); + } + + /** Get the node whose element is the smallest one in the tree + * rooted at this node. + * @return the tree node containing the smallest element in the + * tree rooted at this node or null if the tree is empty + * @see #getLargest + */ + Node getSmallest() { + Node node = this; + while (node.left != null) { + node = node.left; + } + return node; + } + + /** Get the node whose element is the largest one in the tree + * rooted at this node. + * @return the tree node containing the largest element in the + * tree rooted at this node or null if the tree is empty + * @see #getSmallest + */ + Node getLargest() { + Node node = this; + while (node.right != null) { + node = node.right; + } + return node; + } + + /** Get the node containing the next smaller or equal element. + * @return node containing the next smaller or equal element or + * null if there is no smaller or equal element in the tree + * @see #getNext + */ + public Node getPrevious() { + + if (left != null) { + final Node node = left.getLargest(); + if (node != null) { + return node; + } + } + + for (Node node = this; node.parent != null; node = node.parent) { + if (node != node.parent.left) { + return node.parent; + } + } + + return null; + + } + + /** Get the node containing the next larger or equal element. + * @return node containing the next larger or equal element (in + * which case the node is guaranteed not to be empty) or null if + * there is no larger or equal element in the tree + * @see #getPrevious + */ + public Node getNext() { + + if (right != null) { + final Node node = right.getSmallest(); + if (node != null) { + return node; + } + } + + for (Node node = this; node.parent != null; node = node.parent) { + if (node != node.parent.right) { + return node.parent; + } + } + + return null; + + } + + /** Insert an element in a sub-tree. + * @param newElement element to insert + * @return true if the parent tree should be re-Skew.BALANCED + */ + boolean insert(final T newElement) { + if (newElement.compareTo(this.element) < 0) { + // the inserted element is smaller than the node + if (left == null) { + left = new Node(newElement, this); + return rebalanceLeftGrown(); + } + return left.insert(newElement) ? rebalanceLeftGrown() : false; + } + + // the inserted element is equal to or greater than the node + if (right == null) { + right = new Node(newElement, this); + return rebalanceRightGrown(); + } + return right.insert(newElement) ? rebalanceRightGrown() : false; + + } + + /** Delete the node from the tree. + */ + public void delete() { + if ((parent == null) && (left == null) && (right == null)) { + // this was the last node, the tree is now empty + element = null; + top = null; + } else { + + Node node; + Node child; + boolean leftShrunk; + if ((left == null) && (right == null)) { + node = this; + element = null; + leftShrunk = node == node.parent.left; + child = null; + } else { + node = (left != null) ? left.getLargest() : right.getSmallest(); + element = node.element; + leftShrunk = node == node.parent.left; + child = (node.left != null) ? node.left : node.right; + } + + node = node.parent; + if (leftShrunk) { + node.left = child; + } else { + node.right = child; + } + if (child != null) { + child.parent = node; + } + + while (leftShrunk ? node.rebalanceLeftShrunk() : node.rebalanceRightShrunk()) { + if (node.parent == null) { + return; + } + leftShrunk = node == node.parent.left; + node = node.parent; + } + + } + } + + /** Re-balance the instance as left sub-tree has grown. + * @return true if the parent tree should be reSkew.BALANCED too + */ + private boolean rebalanceLeftGrown() { + switch (skew) { + case LEFT_HIGH: + if (left.skew == Skew.LEFT_HIGH) { + rotateCW(); + skew = Skew.BALANCED; + right.skew = Skew.BALANCED; + } else { + final Skew s = left.right.skew; + left.rotateCCW(); + rotateCW(); + switch(s) { + case LEFT_HIGH: + left.skew = Skew.BALANCED; + right.skew = Skew.RIGHT_HIGH; + break; + case RIGHT_HIGH: + left.skew = Skew.LEFT_HIGH; + right.skew = Skew.BALANCED; + break; + default: + left.skew = Skew.BALANCED; + right.skew = Skew.BALANCED; + } + skew = Skew.BALANCED; + } + return false; + case RIGHT_HIGH: + skew = Skew.BALANCED; + return false; + default: + skew = Skew.LEFT_HIGH; + return true; + } + } + + /** Re-balance the instance as right sub-tree has grown. + * @return true if the parent tree should be reSkew.BALANCED too + */ + private boolean rebalanceRightGrown() { + switch (skew) { + case LEFT_HIGH: + skew = Skew.BALANCED; + return false; + case RIGHT_HIGH: + if (right.skew == Skew.RIGHT_HIGH) { + rotateCCW(); + skew = Skew.BALANCED; + left.skew = Skew.BALANCED; + } else { + final Skew s = right.left.skew; + right.rotateCW(); + rotateCCW(); + switch (s) { + case LEFT_HIGH: + left.skew = Skew.BALANCED; + right.skew = Skew.RIGHT_HIGH; + break; + case RIGHT_HIGH: + left.skew = Skew.LEFT_HIGH; + right.skew = Skew.BALANCED; + break; + default: + left.skew = Skew.BALANCED; + right.skew = Skew.BALANCED; + } + skew = Skew.BALANCED; + } + return false; + default: + skew = Skew.RIGHT_HIGH; + return true; + } + } + + /** Re-balance the instance as left sub-tree has shrunk. + * @return true if the parent tree should be reSkew.BALANCED too + */ + private boolean rebalanceLeftShrunk() { + switch (skew) { + case LEFT_HIGH: + skew = Skew.BALANCED; + return true; + case RIGHT_HIGH: + if (right.skew == Skew.RIGHT_HIGH) { + rotateCCW(); + skew = Skew.BALANCED; + left.skew = Skew.BALANCED; + return true; + } else if (right.skew == Skew.BALANCED) { + rotateCCW(); + skew = Skew.LEFT_HIGH; + left.skew = Skew.RIGHT_HIGH; + return false; + } else { + final Skew s = right.left.skew; + right.rotateCW(); + rotateCCW(); + switch (s) { + case LEFT_HIGH: + left.skew = Skew.BALANCED; + right.skew = Skew.RIGHT_HIGH; + break; + case RIGHT_HIGH: + left.skew = Skew.LEFT_HIGH; + right.skew = Skew.BALANCED; + break; + default: + left.skew = Skew.BALANCED; + right.skew = Skew.BALANCED; + } + skew = Skew.BALANCED; + return true; + } + default: + skew = Skew.RIGHT_HIGH; + return false; + } + } + + /** Re-balance the instance as right sub-tree has shrunk. + * @return true if the parent tree should be reSkew.BALANCED too + */ + private boolean rebalanceRightShrunk() { + switch (skew) { + case RIGHT_HIGH: + skew = Skew.BALANCED; + return true; + case LEFT_HIGH: + if (left.skew == Skew.LEFT_HIGH) { + rotateCW(); + skew = Skew.BALANCED; + right.skew = Skew.BALANCED; + return true; + } else if (left.skew == Skew.BALANCED) { + rotateCW(); + skew = Skew.RIGHT_HIGH; + right.skew = Skew.LEFT_HIGH; + return false; + } else { + final Skew s = left.right.skew; + left.rotateCCW(); + rotateCW(); + switch (s) { + case LEFT_HIGH: + left.skew = Skew.BALANCED; + right.skew = Skew.RIGHT_HIGH; + break; + case RIGHT_HIGH: + left.skew = Skew.LEFT_HIGH; + right.skew = Skew.BALANCED; + break; + default: + left.skew = Skew.BALANCED; + right.skew = Skew.BALANCED; + } + skew = Skew.BALANCED; + return true; + } + default: + skew = Skew.LEFT_HIGH; + return false; + } + } + + /** Perform a clockwise rotation rooted at the instance. + *The skew factor are not updated by this method, they + * must be updated by the caller
+ */ + private void rotateCW() { + + final T tmpElt = element; + element = left.element; + left.element = tmpElt; + + final Node tmpNode = left; + left = tmpNode.left; + tmpNode.left = tmpNode.right; + tmpNode.right = right; + right = tmpNode; + + if (left != null) { + left.parent = this; + } + if (right.right != null) { + right.right.parent = right; + } + + } + + /** Perform a counter-clockwise rotation rooted at the instance. + *The skew factor are not updated by this method, they + * must be updated by the caller
+ */ + private void rotateCCW() { + + final T tmpElt = element; + element = right.element; + right.element = tmpElt; + + final Node tmpNode = right; + right = tmpNode.right; + tmpNode.right = tmpNode.left; + tmpNode.left = left; + left = tmpNode; + + if (right != null) { + right.parent = this; + } + if (left.left != null) { + left.left.parent = left; + } + + } + + } + +} diff --git a/src/main/java/org/apache/commons/math/geometry/partitioning/utilities/OrderedTuple.java b/src/main/java/org/apache/commons/math/geometry/partitioning/utilities/OrderedTuple.java new file mode 100644 index 000000000..71bfec062 --- /dev/null +++ b/src/main/java/org/apache/commons/math/geometry/partitioning/utilities/OrderedTuple.java @@ -0,0 +1,417 @@ +/* + * 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.math.geometry.partitioning.utilities; + +import java.util.Arrays; + +import org.apache.commons.math.util.FastMath; + +/** This class implements an ordering operation for T-uples. + + *Ordering is done by encoding all components of the T-uple into a + * single scalar value and using this value as the sorting + * key. Encoding is performed using the method invented by Georg + * Cantor in 1877 when he proved it was possible to establish a + * bijection between a line and a plane. The binary representations of + * the components of the T-uple are mixed together to form a single + * scalar. This means that the 2k bit of component 0 is + * followed by the 2k bit of component 1, then by the + * 2k bit of component 2 up to the 2k bit of + * component {@code t}, which is followed by the 2k-1 + * bit of component 0, followed by the 2k-1 bit of + * component 1 ... The binary representations are extended as needed + * to handle numbers with different scales and a suitable + * 2p offset is added to the components in order to avoid + * negative numbers (this offset is adjusted as needed during the + * comparison operations).
+ + *The more interesting property of the encoding method for our + * purpose is that it allows to select all the points that are in a + * given range. This is depicted in dimension 2 by the following + * picure:
+ + *
+
+ * This picture shows a set of 100000 random 2-D pairs having their + * first component between -50 and +150 and their second component + * between -350 and +50. We wanted to extract all pairs having their + * first component between +30 and +70 and their second component + * between -120 and -30. We built the lower left point at coordinates + * (30, -120) and the upper right point at coordinates (70, -30). All + * points smaller than the lower left point are drawn in red and all + * points larger than the upper right point are drawn in blue. The + * green points are between the two limits. This picture shows that + * all the desired points are selected, along with spurious points. In + * this case, we get 15790 points, 4420 of which really belonging to + * the desired rectangle. It is possible to extract very small + * subsets. As an example extracting from the same 100000 points set + * the points having their first component between +30 and +31 and + * their second component between -91 and -90, we get a subset of 11 + * points, 2 of which really belonging to the desired rectangle.
+ + *the previous selection technique can be applied in all + * dimensions, still using two points to define the interval. The + * first point will have all its components set to their lower bounds + * while the second point will have all its components set to their + * upper bounds.
+ + *T-uples with negative infinite or positive infinite components + * are sorted logically.
+ + *Since the specification of the {@code Comparator} interface + * allows only {@code ClassCastException} errors, some arbitrary + * choices have been made to handle specific cases. The rationale for + * these choices is to keep regular and consistent T-uples + * together.
+ *-
+ *
- instances with different dimensions are sorted according to + * their dimension regardless of their components values + *
- instances with {@code Double.NaN} components are sorted + * after all other ones (even after instances with positive infinite + * components + *
- instances with both positive and negative infinite components + * are considered as if they had {@code Double.NaN} + * components + *
The ordering method is detailed in the general description of + * the class. Its main property is to be consistent with distance: + * geometrically close T-uples stay close to each other when stored + * in a sorted collection using this comparison method.
+ + *T-uples with negative infinite, positive infinite are sorted + * logically.
+ + *Some arbitrary choices have been made to handle specific + * cases. The rationale for these choices is to keep + * normal and consistent T-uples together.
+ *-
+ *
- instances with different dimensions are sorted according to + * their dimension regardless of their components values + *
- instances with {@code Double.NaN} components are sorted + * after all other ones (evan after instances with positive infinite + * components + *
- instances with both positive and negative infinite components + * are considered as if they had {@code Double.NaN} + * components + *
+This package provides multidimensional ordering features for partitioning. +
+ + diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index 7cc1ffd7a..bd2043c8b 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -52,6 +52,16 @@ TheThe geometry package provides classes useful for many physical simulations - in the real 3D space, namely vectors and rotations. + in Euclidean spaces, like vectors and rotations in 3D, as well as a general + implentation of Binary Space Partitioning Trees (BSP trees).
+ + Interval and + IntervalsSet represent one dimensional regions. All classical set operations are available + for intervals sets: union, intersection, symmetric difference (exclusive or), difference, complement, + as well as region predicates (point inside/outside/on boundary, emptiness, other region contained). + It is also possible to compute geometrical properties like size, barycenter or boundary size. + Intervals sets can be built by constructive geometry (union, intersection ...) or from a boundary + representation. +
++ + PolygonsSet represent two dimensional regions. All classical set operations are available + for polygons sets: union, intersection, symmetric difference (exclusive or), difference, complement, + as well as region predicates (point inside/outside/on boundary, emptiness, other region contained). + It is also possible to compute geometrical properties like size, barycenter or boundary size and + to extract the vertices. Polygons sets can be built by constructive geometry (union, intersection ...) + or from a boundary representation. +
++ + PolyhedronsSet represent three dimensional regions. All classical set operations are available + for polyhedrons sets: union, intersection, symmetric difference (exclusive or), difference, complement, + as well as region predicates (point inside/outside/on boundary, emptiness, other region contained). + It is also possible to compute geometrical properties like size, barycenter or boundary size and + to extract the vertices. Polyhedrons sets can be built by constructive geometry (union, intersection ...) + or from a boundary representation. +
- + Vector3D provides a simple vector type. One important feature is that instances of this class are guaranteed to be immutable, this greatly simplifies modelling dynamical systems @@ -57,14 +86,12 @@ is of course also implemented.
- + Vector3DFormat is a specialized format for formatting output or parsing input with text representation of 3D vectors.
-- + Rotation represents 3D rotations. Rotation instances are also immutable objects, as Vector3D instances.
@@ -136,6 +163,41 @@applyInverseTo(Rotation).
+ + BSP trees are an efficient way to represent space partitions and + to associate attributes with each cell. Each node in a BSP tree + represents a convex region which is partitioned in two convex + sub-regions at each side of a cut hyperplane. The root tree + contains the complete space. +
+ ++ The main use of such partitions is to use a boolean attribute to + define an inside/outside property, hence representing arbitrary + polytopes (line segments in 1D, polygons in 2D and polyhedrons in + 3D) and to operate on them. +
+ ++ Another example would be to represent Voronoi tesselations, the + attribute of each cell holding the defining point of the cell. +
+ ++ The application-defined attributes are shared among copied + instances and propagated to split parts. These attributes are not + used by the BSP-tree algorithms themselves, so the application can + use them for any purpose. Since the tree visiting method holds + internal and leaf nodes differently, it is possible to use + different classes for internal nodes attributes and leaf nodes + attributes. This should be used with care, though, because if the + tree is modified in any way after attributes have been set, some + internal nodes may become leaf nodes and some leaf nodes may become + internal nodes. +
+