diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/geometry/partitioning')
24 files changed, 4877 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java new file mode 100644 index 0000000..d901ab4 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java @@ -0,0 +1,540 @@ +/* + * 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.math3.geometry.partitioning; + +import java.util.ArrayList; +import java.util.Collection; +import java.util.Comparator; +import java.util.HashMap; +import java.util.Iterator; +import java.util.Map; +import java.util.TreeSet; + +import org.apache.commons.math3.geometry.Point; +import org.apache.commons.math3.geometry.Space; +import org.apache.commons.math3.geometry.Vector; + +/** Abstract class for all regions, independently of geometry type or dimension. + + * @param <S> Type of the space. + * @param <T> Type of the sub-space. + + * @since 3.0 + */ +public abstract class AbstractRegion<S extends Space, T extends Space> implements Region<S> { + + /** Inside/Outside BSP tree. */ + private BSPTree<S> tree; + + /** Tolerance below which points are considered to belong to hyperplanes. */ + private final double tolerance; + + /** Size of the instance. */ + private double size; + + /** Barycenter. */ + private Point<S> barycenter; + + /** Build a region representing the whole space. + * @param tolerance tolerance below which points are considered identical. + */ + protected AbstractRegion(final double tolerance) { + this.tree = new BSPTree<S>(Boolean.TRUE); + this.tolerance = tolerance; + } + + /** Build a region from an inside/outside BSP tree. + * <p>The leaf nodes of the BSP tree <em>must</em> 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 <em>must</em> have either null internal nodes or + * internal nodes representing the boundary as specified in the + * {@link #getTree getTree} method).</p> + * @param tree inside/outside BSP tree representing the region + * @param tolerance tolerance below which points are considered identical. + */ + protected AbstractRegion(final BSPTree<S> tree, final double tolerance) { + this.tree = tree; + this.tolerance = tolerance; + } + + /** Build a Region from a Boundary REPresentation (B-rep). + * <p>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.</p> + * <p>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.</p> + * <p>If the boundary is empty, the region will represent the whole + * space.</p> + * @param boundary collection of boundary elements, as a + * collection of {@link SubHyperplane SubHyperplane} objects + * @param tolerance tolerance below which points are considered identical. + */ + protected AbstractRegion(final Collection<SubHyperplane<S>> boundary, final double tolerance) { + + this.tolerance = tolerance; + + if (boundary.size() == 0) { + + // the tree represents the whole space + tree = new BSPTree<S>(Boolean.TRUE); + + } else { + + // sort the boundary elements in decreasing size order + // (we don't want equal size elements to be removed, so + // we use a trick to fool the TreeSet) + final TreeSet<SubHyperplane<S>> ordered = new TreeSet<SubHyperplane<S>>(new Comparator<SubHyperplane<S>>() { + /** {@inheritDoc} */ + public int compare(final SubHyperplane<S> o1, final SubHyperplane<S> o2) { + final double size1 = o1.getSize(); + final double size2 = o2.getSize(); + return (size2 < size1) ? -1 : ((o1 == o2) ? 0 : +1); + } + }); + ordered.addAll(boundary); + + // build the tree top-down + tree = new BSPTree<S>(); + insertCuts(tree, ordered); + + // set up the inside/outside flags + tree.visit(new BSPTreeVisitor<S>() { + + /** {@inheritDoc} */ + public Order visitOrder(final BSPTree<S> node) { + return Order.PLUS_SUB_MINUS; + } + + /** {@inheritDoc} */ + public void visitInternalNode(final BSPTree<S> node) { + } + + /** {@inheritDoc} */ + public void visitLeafNode(final BSPTree<S> node) { + if (node.getParent() == null || node == node.getParent().getMinus()) { + node.setAttribute(Boolean.TRUE); + } else { + node.setAttribute(Boolean.FALSE); + } + } + }); + + } + + } + + /** Build a convex region from an array of bounding hyperplanes. + * @param hyperplanes array of bounding hyperplanes (if null, an + * empty region will be built) + * @param tolerance tolerance below which points are considered identical. + */ + public AbstractRegion(final Hyperplane<S>[] hyperplanes, final double tolerance) { + this.tolerance = tolerance; + if ((hyperplanes == null) || (hyperplanes.length == 0)) { + tree = new BSPTree<S>(Boolean.FALSE); + } else { + + // use the first hyperplane to build the right class + tree = hyperplanes[0].wholeSpace().getTree(false); + + // chop off parts of the space + BSPTree<S> node = tree; + node.setAttribute(Boolean.TRUE); + for (final Hyperplane<S> hyperplane : hyperplanes) { + if (node.insertCut(hyperplane)) { + node.setAttribute(null); + node.getPlus().setAttribute(Boolean.FALSE); + node = node.getMinus(); + node.setAttribute(Boolean.TRUE); + } + } + + } + + } + + /** {@inheritDoc} */ + public abstract AbstractRegion<S, T> buildNew(BSPTree<S> newTree); + + /** Get the tolerance below which points are considered to belong to hyperplanes. + * @return tolerance below which points are considered to belong to hyperplanes + */ + public double getTolerance() { + return tolerance; + } + + /** 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<S> node, final Collection<SubHyperplane<S>> boundary) { + + final Iterator<SubHyperplane<S>> iterator = boundary.iterator(); + + // build the current level + Hyperplane<S> inserted = null; + while ((inserted == null) && iterator.hasNext()) { + inserted = iterator.next().getHyperplane(); + if (!node.insertCut(inserted.copySelf())) { + inserted = null; + } + } + + if (!iterator.hasNext()) { + return; + } + + // distribute the remaining edges in the two sub-trees + final ArrayList<SubHyperplane<S>> plusList = new ArrayList<SubHyperplane<S>>(); + final ArrayList<SubHyperplane<S>> minusList = new ArrayList<SubHyperplane<S>>(); + while (iterator.hasNext()) { + final SubHyperplane<S> other = iterator.next(); + final SubHyperplane.SplitSubHyperplane<S> split = other.split(inserted); + switch (split.getSide()) { + case PLUS: + plusList.add(other); + break; + case MINUS: + minusList.add(other); + break; + case BOTH: + plusList.add(split.getPlus()); + minusList.add(split.getMinus()); + break; + default: + // ignore the sub-hyperplanes belonging to the cut hyperplane + } + } + + // recurse through lower levels + insertCuts(node.getPlus(), plusList); + insertCuts(node.getMinus(), minusList); + + } + + /** {@inheritDoc} */ + public AbstractRegion<S, T> copySelf() { + return buildNew(tree.copySelf()); + } + + /** {@inheritDoc} */ + public boolean isEmpty() { + return isEmpty(tree); + } + + /** {@inheritDoc} */ + public boolean isEmpty(final BSPTree<S> 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 !((Boolean) node.getAttribute()); + } + + // check both sides of the sub-tree + return isEmpty(node.getMinus()) && isEmpty(node.getPlus()); + + } + + /** {@inheritDoc} */ + public boolean isFull() { + return isFull(tree); + } + + /** {@inheritDoc} */ + public boolean isFull(final BSPTree<S> 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 outside cell + + if (node.getCut() == null) { + // if we find an outside node, the region does not cover full space + return (Boolean) node.getAttribute(); + } + + // check both sides of the sub-tree + return isFull(node.getMinus()) && isFull(node.getPlus()); + + } + + /** {@inheritDoc} */ + public boolean contains(final Region<S> region) { + return new RegionFactory<S>().difference(region, this).isEmpty(); + } + + /** {@inheritDoc} + * @since 3.3 + */ + public BoundaryProjection<S> projectToBoundary(final Point<S> point) { + final BoundaryProjector<S, T> projector = new BoundaryProjector<S, T>(point); + getTree(true).visit(projector); + return projector.getProjection(); + } + + /** Check a point with respect to the region. + * @param point point to check + * @return a code representing the point status: either {@link + * Region.Location#INSIDE}, {@link Region.Location#OUTSIDE} or + * {@link Region.Location#BOUNDARY} + */ + public Location checkPoint(final Vector<S> point) { + return checkPoint((Point<S>) point); + } + + /** {@inheritDoc} */ + public Location checkPoint(final Point<S> 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 + * Region.Location#INSIDE INSIDE}, {@link Region.Location#OUTSIDE + * OUTSIDE} or {@link Region.Location#BOUNDARY BOUNDARY} + */ + protected Location checkPoint(final BSPTree<S> node, final Vector<S> point) { + return checkPoint(node, (Point<S>) 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 + * Region.Location#INSIDE INSIDE}, {@link Region.Location#OUTSIDE + * OUTSIDE} or {@link Region.Location#BOUNDARY BOUNDARY} + */ + protected Location checkPoint(final BSPTree<S> node, final Point<S> point) { + final BSPTree<S> cell = node.getCell(point, tolerance); + if (cell.getCut() == null) { + // the point is in the interior of a cell, just check the attribute + return ((Boolean) cell.getAttribute()) ? 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; + + } + + /** {@inheritDoc} */ + public BSPTree<S> getTree(final boolean includeBoundaryAttributes) { + if (includeBoundaryAttributes && (tree.getCut() != null) && (tree.getAttribute() == null)) { + // compute the boundary attributes + tree.visit(new BoundaryBuilder<S>()); + } + return tree; + } + + /** {@inheritDoc} */ + public double getBoundarySize() { + final BoundarySizeVisitor<S> visitor = new BoundarySizeVisitor<S>(); + getTree(true).visit(visitor); + return visitor.getSize(); + } + + /** {@inheritDoc} */ + 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; + } + + /** {@inheritDoc} */ + public Point<S> getBarycenter() { + if (barycenter == null) { + computeGeometricalProperties(); + } + return barycenter; + } + + /** Set the barycenter of the instance. + * @param barycenter barycenter of the instance + */ + protected void setBarycenter(final Vector<S> barycenter) { + setBarycenter((Point<S>) barycenter); + } + + /** Set the barycenter of the instance. + * @param barycenter barycenter of the instance + */ + protected void setBarycenter(final Point<S> barycenter) { + this.barycenter = barycenter; + } + + /** Compute some geometrical properties. + * <p>The properties to compute are the barycenter and the size.</p> + */ + protected abstract void computeGeometricalProperties(); + + /** {@inheritDoc} */ + @Deprecated + public Side side(final Hyperplane<S> hyperplane) { + final InsideFinder<S> finder = new InsideFinder<S>(this); + finder.recurseSides(tree, hyperplane.wholeHyperplane()); + return finder.plusFound() ? + (finder.minusFound() ? Side.BOTH : Side.PLUS) : + (finder.minusFound() ? Side.MINUS : Side.HYPER); + } + + /** {@inheritDoc} */ + public SubHyperplane<S> intersection(final SubHyperplane<S> 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<S> recurseIntersection(final BSPTree<S> node, final SubHyperplane<S> sub) { + + if (node.getCut() == null) { + return (Boolean) node.getAttribute() ? sub.copySelf() : null; + } + + final Hyperplane<S> hyperplane = node.getCut().getHyperplane(); + final SubHyperplane.SplitSubHyperplane<S> split = sub.split(hyperplane); + if (split.getPlus() != null) { + if (split.getMinus() != null) { + // both sides + final SubHyperplane<S> plus = recurseIntersection(node.getPlus(), split.getPlus()); + final SubHyperplane<S> minus = recurseIntersection(node.getMinus(), split.getMinus()); + if (plus == null) { + return minus; + } else if (minus == null) { + return plus; + } else { + return plus.reunite(minus); + } + } else { + // only on plus side + return recurseIntersection(node.getPlus(), sub); + } + } else if (split.getMinus() != null) { + // only on minus side + return recurseIntersection(node.getMinus(), sub); + } else { + // on hyperplane + return recurseIntersection(node.getPlus(), + recurseIntersection(node.getMinus(), sub)); + } + + } + + /** Transform a region. + * <p>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.</p> + * @param transform transform to apply + * @return a new region, resulting from the application of the + * transform to the instance + */ + public AbstractRegion<S, T> applyTransform(final Transform<S, T> transform) { + + // transform the tree, except for boundary attribute splitters + final Map<BSPTree<S>, BSPTree<S>> map = new HashMap<BSPTree<S>, BSPTree<S>>(); + final BSPTree<S> transformedTree = recurseTransform(getTree(false), transform, map); + + // set up the boundary attributes splitters + for (final Map.Entry<BSPTree<S>, BSPTree<S>> entry : map.entrySet()) { + if (entry.getKey().getCut() != null) { + @SuppressWarnings("unchecked") + BoundaryAttribute<S> original = (BoundaryAttribute<S>) entry.getKey().getAttribute(); + if (original != null) { + @SuppressWarnings("unchecked") + BoundaryAttribute<S> transformed = (BoundaryAttribute<S>) entry.getValue().getAttribute(); + for (final BSPTree<S> splitter : original.getSplitters()) { + transformed.getSplitters().add(map.get(splitter)); + } + } + } + } + + return buildNew(transformedTree); + + } + + /** Recursively transform an inside/outside BSP-tree. + * @param node current BSP tree node + * @param transform transform to apply + * @param map transformed nodes map + * @return a new tree + */ + @SuppressWarnings("unchecked") + private BSPTree<S> recurseTransform(final BSPTree<S> node, final Transform<S, T> transform, + final Map<BSPTree<S>, BSPTree<S>> map) { + + final BSPTree<S> transformedNode; + if (node.getCut() == null) { + transformedNode = new BSPTree<S>(node.getAttribute()); + } else { + + final SubHyperplane<S> sub = node.getCut(); + final SubHyperplane<S> tSub = ((AbstractSubHyperplane<S, T>) sub).applyTransform(transform); + BoundaryAttribute<S> attribute = (BoundaryAttribute<S>) node.getAttribute(); + if (attribute != null) { + final SubHyperplane<S> tPO = (attribute.getPlusOutside() == null) ? + null : ((AbstractSubHyperplane<S, T>) attribute.getPlusOutside()).applyTransform(transform); + final SubHyperplane<S> tPI = (attribute.getPlusInside() == null) ? + null : ((AbstractSubHyperplane<S, T>) attribute.getPlusInside()).applyTransform(transform); + // we start with an empty list of splitters, it will be filled in out of recursion + attribute = new BoundaryAttribute<S>(tPO, tPI, new NodesSet<S>()); + } + + transformedNode = new BSPTree<S>(tSub, + recurseTransform(node.getPlus(), transform, map), + recurseTransform(node.getMinus(), transform, map), + attribute); + } + + map.put(node, transformedNode); + return transformedNode; + + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractSubHyperplane.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractSubHyperplane.java new file mode 100644 index 0000000..396b795 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractSubHyperplane.java @@ -0,0 +1,191 @@ +/* + * 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.math3.geometry.partitioning; + +import java.util.HashMap; +import java.util.Map; + +import org.apache.commons.math3.geometry.Space; + +/** This class implements the dimension-independent parts of {@link SubHyperplane}. + + * <p>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.</p> + + * @param <S> Type of the embedding space. + * @param <T> Type of the embedded sub-space. + + * @since 3.0 + */ +public abstract class AbstractSubHyperplane<S extends Space, T extends Space> + implements SubHyperplane<S> { + + /** Underlying hyperplane. */ + private final Hyperplane<S> hyperplane; + + /** Remaining region of the hyperplane. */ + private final Region<T> remainingRegion; + + /** Build a sub-hyperplane from an hyperplane and a region. + * @param hyperplane underlying hyperplane + * @param remainingRegion remaining region of the hyperplane + */ + protected AbstractSubHyperplane(final Hyperplane<S> hyperplane, + final Region<T> remainingRegion) { + this.hyperplane = hyperplane; + this.remainingRegion = remainingRegion; + } + + /** Build a sub-hyperplane from an hyperplane and a region. + * @param hyper underlying hyperplane + * @param remaining remaining region of the hyperplane + * @return a new sub-hyperplane + */ + protected abstract AbstractSubHyperplane<S, T> buildNew(final Hyperplane<S> hyper, + final Region<T> remaining); + + /** {@inheritDoc} */ + public AbstractSubHyperplane<S, T> copySelf() { + return buildNew(hyperplane.copySelf(), remainingRegion); + } + + /** Get the underlying hyperplane. + * @return underlying hyperplane + */ + public Hyperplane<S> getHyperplane() { + return hyperplane; + } + + /** Get the remaining region of the hyperplane. + * <p>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.</p> + * @return remaining region of the hyperplane + */ + public Region<T> getRemainingRegion() { + return remainingRegion; + } + + /** {@inheritDoc} */ + public double getSize() { + return remainingRegion.getSize(); + } + + /** {@inheritDoc} */ + public AbstractSubHyperplane<S, T> reunite(final SubHyperplane<S> other) { + @SuppressWarnings("unchecked") + AbstractSubHyperplane<S, T> o = (AbstractSubHyperplane<S, T>) other; + return buildNew(hyperplane, + new RegionFactory<T>().union(remainingRegion, o.remainingRegion)); + } + + /** Apply a transform to the instance. + * <p>The instance must be a (D-1)-dimension sub-hyperplane with + * respect to the transform <em>not</em> 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.</p> + * @param transform D-dimension transform to apply + * @return the transformed instance + */ + public AbstractSubHyperplane<S, T> applyTransform(final Transform<S, T> transform) { + final Hyperplane<S> tHyperplane = transform.apply(hyperplane); + + // transform the tree, except for boundary attribute splitters + final Map<BSPTree<T>, BSPTree<T>> map = new HashMap<BSPTree<T>, BSPTree<T>>(); + final BSPTree<T> tTree = + recurseTransform(remainingRegion.getTree(false), tHyperplane, transform, map); + + // set up the boundary attributes splitters + for (final Map.Entry<BSPTree<T>, BSPTree<T>> entry : map.entrySet()) { + if (entry.getKey().getCut() != null) { + @SuppressWarnings("unchecked") + BoundaryAttribute<T> original = (BoundaryAttribute<T>) entry.getKey().getAttribute(); + if (original != null) { + @SuppressWarnings("unchecked") + BoundaryAttribute<T> transformed = (BoundaryAttribute<T>) entry.getValue().getAttribute(); + for (final BSPTree<T> splitter : original.getSplitters()) { + transformed.getSplitters().add(map.get(splitter)); + } + } + } + } + + return buildNew(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 + * @param map transformed nodes map + * @return a new tree + */ + private BSPTree<T> recurseTransform(final BSPTree<T> node, + final Hyperplane<S> transformed, + final Transform<S, T> transform, + final Map<BSPTree<T>, BSPTree<T>> map) { + + final BSPTree<T> transformedNode; + if (node.getCut() == null) { + transformedNode = new BSPTree<T>(node.getAttribute()); + } else { + + @SuppressWarnings("unchecked") + BoundaryAttribute<T> attribute = (BoundaryAttribute<T>) node.getAttribute(); + if (attribute != null) { + final SubHyperplane<T> tPO = (attribute.getPlusOutside() == null) ? + null : transform.apply(attribute.getPlusOutside(), hyperplane, transformed); + final SubHyperplane<T> tPI = (attribute.getPlusInside() == null) ? + null : transform.apply(attribute.getPlusInside(), hyperplane, transformed); + // we start with an empty list of splitters, it will be filled in out of recursion + attribute = new BoundaryAttribute<T>(tPO, tPI, new NodesSet<T>()); + } + + transformedNode = new BSPTree<T>(transform.apply(node.getCut(), hyperplane, transformed), + recurseTransform(node.getPlus(), transformed, transform, map), + recurseTransform(node.getMinus(), transformed, transform, map), + attribute); + } + + map.put(node, transformedNode); + return transformedNode; + + } + + /** {@inheritDoc} */ + @Deprecated + public Side side(Hyperplane<S> hyper) { + return split(hyper).getSide(); + } + + /** {@inheritDoc} */ + public abstract SplitSubHyperplane<S> split(Hyperplane<S> hyper); + + /** {@inheritDoc} */ + public boolean isEmpty() { + return remainingRegion.isEmpty(); + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java new file mode 100644 index 0000000..1f1a6ea --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java @@ -0,0 +1,821 @@ +/* + * 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.math3.geometry.partitioning; + +import java.util.ArrayList; +import java.util.List; + +import org.apache.commons.math3.exception.MathIllegalStateException; +import org.apache.commons.math3.exception.MathInternalError; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.geometry.Point; +import org.apache.commons.math3.geometry.Space; +import org.apache.commons.math3.geometry.Vector; +import org.apache.commons.math3.util.FastMath; + +/** This class represent a Binary Space Partition tree. + + * <p>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.</p> + + * <p>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.</p> + + * <p>Another example would be to represent Voronoi tesselations, the + * attribute of each cell holding the defining point of the cell.</p> + + * <p>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.</p> + + * <p>One of the main sources for the development of this package was + * Bruce Naylor, John Amanatides and William Thibault paper <a + * href="http://www.cs.yorku.ca/~amana/research/bsptSetOp.pdf">Merging + * BSP Trees Yields Polyhedral Set Operations</a> Proc. Siggraph '90, + * Computer Graphics 24(4), August 1990, pp 115-124, published by the + * Association for Computing Machinery (ACM).</p> + + * @param <S> Type of the space. + + * @since 3.0 + */ +public class BSPTree<S extends Space> { + + /** Cut sub-hyperplane. */ + private SubHyperplane<S> cut; + + /** Tree at the plus side of the cut hyperplane. */ + private BSPTree<S> plus; + + /** Tree at the minus side of the cut hyperplane. */ + private BSPTree<S> minus; + + /** Parent tree. */ + private BSPTree<S> 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. + * <p>This method does <em>not</em> perform any verification on + * consistency of its arguments, it should therefore be used only + * when then caller knows what it is doing.</p> + * <p>This method is mainly useful to build trees + * bottom-up. Building trees top-down is realized with the help of + * method {@link #insertCut insertCut}.</p> + * @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<S> cut, final BSPTree<S> plus, final BSPTree<S> 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. + * <p>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 <em>no</em> cut hyperplane will be + * inserted and the cell will be changed to a leaf cell. The + * attribute of the node is never changed.</p> + * <p>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).</p> + * @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<S> hyperplane) { + + if (cut != null) { + plus.parent = null; + minus.parent = null; + } + + final SubHyperplane<S> chopped = fitToCell(hyperplane.wholeHyperplane()); + if (chopped == null || chopped.isEmpty()) { + cut = null; + plus = null; + minus = null; + return false; + } + + cut = chopped; + plus = new BSPTree<S>(); + plus.parent = this; + minus = new BSPTree<S>(); + minus.parent = this; + return true; + + } + + /** Copy the instance. + * <p>The instance created is completely independent of the original + * one. A deep copy is used, none of the underlying objects are + * shared (except for the nodes attributes and immutable + * objects).</p> + * @return a new tree, copy of the instance + */ + public BSPTree<S> copySelf() { + + if (cut == null) { + return new BSPTree<S>(attribute); + } + + return new BSPTree<S>(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<S> 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<S> 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<S> getMinus() { + return minus; + } + + /** Get the parent node. + * @return parent node, null if the node has no parents + */ + public BSPTree<S> 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<S> 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 MathInternalError(); + } + + } + } + + /** Fit a sub-hyperplane inside the cell defined by the instance. + * <p>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.</p> + * @param sub sub-hyperplane to fit + * @return a new sub-hyperplane, guaranteed to have no part outside + * of the instance cell + */ + private SubHyperplane<S> fitToCell(final SubHyperplane<S> sub) { + SubHyperplane<S> s = sub; + for (BSPTree<S> tree = this; tree.parent != null && s != null; tree = tree.parent) { + if (tree == tree.parent.plus) { + s = s.split(tree.parent.cut.getHyperplane()).getPlus(); + } else { + s = s.split(tree.parent.cut.getHyperplane()).getMinus(); + } + } + return s; + } + + /** Get the cell to which a point belongs. + * <p>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.</p> + * @param point point to check + * @return the tree cell to which the point belongs + * @deprecated as of 3.3, replaced with {@link #getCell(Point, double)} + */ + @Deprecated + public BSPTree<S> getCell(final Vector<S> point) { + return getCell((Point<S>) point, 1.0e-10); + } + + /** Get the cell to which a point belongs. + * <p>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.</p> + * @param point point to check + * @param tolerance tolerance below which points close to a cut hyperplane + * are considered to belong to the hyperplane itself + * @return the tree cell to which the point belongs + */ + public BSPTree<S> getCell(final Point<S> point, final double tolerance) { + + 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) < tolerance) { + return this; + } else if (offset <= 0) { + // point is on the minus side of the cut hyperplane + return minus.getCell(point, tolerance); + } else { + // point is on the plus side of the cut hyperplane + return plus.getCell(point, tolerance); + } + + } + + /** Get the cells whose cut sub-hyperplanes are close to the point. + * @param point point to check + * @param maxOffset offset below which a cut sub-hyperplane is considered + * close to the point (in absolute value) + * @return close cells (may be empty if all cut sub-hyperplanes are farther + * than maxOffset from the point) + */ + public List<BSPTree<S>> getCloseCuts(final Point<S> point, final double maxOffset) { + final List<BSPTree<S>> close = new ArrayList<BSPTree<S>>(); + recurseCloseCuts(point, maxOffset, close); + return close; + } + + /** Get the cells whose cut sub-hyperplanes are close to the point. + * @param point point to check + * @param maxOffset offset below which a cut sub-hyperplane is considered + * close to the point (in absolute value) + * @param close list to fill + */ + private void recurseCloseCuts(final Point<S> point, final double maxOffset, + final List<BSPTree<S>> close) { + if (cut != null) { + + // position of the point with respect to the cut hyperplane + final double offset = cut.getHyperplane().getOffset(point); + + if (offset < -maxOffset) { + // point is on the minus side of the cut hyperplane + minus.recurseCloseCuts(point, maxOffset, close); + } else if (offset > maxOffset) { + // point is on the plus side of the cut hyperplane + plus.recurseCloseCuts(point, maxOffset, close); + } else { + // point is close to the cut hyperplane + close.add(this); + minus.recurseCloseCuts(point, maxOffset, close); + plus.recurseCloseCuts(point, maxOffset, close); + } + + } + } + + /** Perform condensation on a tree. + * <p>The condensation operation is not recursive, it must be called + * explicitly from leaves to root.</p> + */ + 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. + * <p>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, <em>no</em> such copy is done here!</p> + * <p>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).</p> + * @param tree other tree to merge with the instance (will be + * <em>unusable</em> 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 of <code>instance <op> + * tree</code>, this value can be ignored if parentTree is not null + * since all connections have already been established + */ + public BSPTree<S> merge(final BSPTree<S> tree, final LeafMerger<S> 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 + * <em>unusable</em> 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 <code>instance <op> + * tree</code>, this value can be ignored if parentTree is not null + * since all connections have already been established + */ + private BSPTree<S> merge(final BSPTree<S> tree, final LeafMerger<S> leafMerger, + final BSPTree<S> 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<S> 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(merged.cut.getHyperplane().wholeHyperplane()); + } + + return merged; + + } + } + + /** This interface gather the merging operations between a BSP tree + * leaf and another BSP tree. + * <p>As explained in Bruce Naylor, John Amanatides and William + * Thibault paper <a + * href="http://www.cs.yorku.ca/~amana/research/bsptSetOp.pdf">Merging + * BSP Trees Yields Polyhedral Set Operations</a>, + * 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).</p> + * @param <S> Type of the space. + */ + public interface LeafMerger<S extends Space> { + + /** Merge a leaf node and a tree node. + * <p>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.</p> + * <p>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.</p> + * <p>The {@link BSPTree#insertInTree BSPTree.insertInTree} method + * may be useful to implement this method.</p> + * @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<S> merge(BSPTree<S> leaf, BSPTree<S> tree, BSPTree<S> parentTree, + boolean isPlusChild, boolean leafFromInstance); + + } + + /** This interface handles the corner cases when an internal node cut sub-hyperplane vanishes. + * <p> + * Such cases happens for example when a cut sub-hyperplane is inserted into + * another tree (during a merge operation), and is split in several parts, + * some of which becomes smaller than the tolerance. The corresponding node + * as then no cut sub-hyperplane anymore, but does have children. This interface + * specifies how to handle this situation. + * setting + * </p> + * @param <S> Type of the space. + * @since 3.4 + */ + public interface VanishingCutHandler<S extends Space> { + + /** Fix a node with both vanished cut and children. + * @param node node to fix + * @return fixed node + */ + BSPTree<S> fixNode(BSPTree<S> node); + + } + + /** Split a BSP tree by an external sub-hyperplane. + * <p>Split a tree in two halves, on each side of the + * sub-hyperplane. The instance is not modified.</p> + * <p>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 <em>not</em> + * attached to any parent tree yet. This tree is intended to be + * later inserted into an higher level tree.</p> + * <p>The algorithm used here is the one given in Naylor, Amanatides + * and Thibault paper (section III, Binary Partitioning of a BSP + * Tree).</p> + * @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<S> split(final SubHyperplane<S> sub) { + + if (cut == null) { + return new BSPTree<S>(sub, copySelf(), new BSPTree<S>(attribute), null); + } + + final Hyperplane<S> cHyperplane = cut.getHyperplane(); + final Hyperplane<S> sHyperplane = sub.getHyperplane(); + final SubHyperplane.SplitSubHyperplane<S> subParts = sub.split(cHyperplane); + switch (subParts.getSide()) { + case PLUS : + { // the partitioning sub-hyperplane is entirely in the plus sub-tree + final BSPTree<S> split = plus.split(sub); + if (cut.split(sHyperplane).getSide() == Side.PLUS) { + split.plus = + new BSPTree<S>(cut.copySelf(), split.plus, minus.copySelf(), attribute); + split.plus.condense(); + split.plus.parent = split; + } else { + split.minus = + new BSPTree<S>(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<S> split = minus.split(sub); + if (cut.split(sHyperplane).getSide() == Side.PLUS) { + split.plus = + new BSPTree<S>(cut.copySelf(), plus.copySelf(), split.plus, attribute); + split.plus.condense(); + split.plus.parent = split; + } else { + split.minus = + new BSPTree<S>(cut.copySelf(), plus.copySelf(), split.minus, attribute); + split.minus.condense(); + split.minus.parent = split; + } + return split; + } + case BOTH : + { + final SubHyperplane.SplitSubHyperplane<S> cutParts = cut.split(sHyperplane); + final BSPTree<S> split = + new BSPTree<S>(sub, plus.split(subParts.getPlus()), minus.split(subParts.getMinus()), + null); + split.plus.cut = cutParts.getPlus(); + split.minus.cut = cutParts.getMinus(); + final BSPTree<S> 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<S>(sub, plus.copySelf(), minus.copySelf(), attribute) : + new BSPTree<S>(sub, minus.copySelf(), plus.copySelf(), attribute); + } + + } + + /** Insert the instance into another tree. + * <p>The instance itself is modified so its former parent should + * not be used anymore.</p> + * @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 + * @deprecated as of 3.4, replaced with {@link #insertInTree(BSPTree, boolean, VanishingCutHandler)} + */ + @Deprecated + public void insertInTree(final BSPTree<S> parentTree, final boolean isPlusChild) { + insertInTree(parentTree, isPlusChild, new VanishingCutHandler<S>() { + /** {@inheritDoc} */ + public BSPTree<S> fixNode(BSPTree<S> node) { + // the cut should not be null + throw new MathIllegalStateException(LocalizedFormats.NULL_NOT_ALLOWED); + } + }); + } + + /** Insert the instance into another tree. + * <p>The instance itself is modified so its former parent should + * not be used anymore.</p> + * @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 vanishingHandler handler to use for handling very rare corner + * cases of vanishing cut sub-hyperplanes in internal nodes during merging + * @see LeafMerger + * @since 3.4 + */ + public void insertInTree(final BSPTree<S> parentTree, final boolean isPlusChild, + final VanishingCutHandler<S> vanishingHandler) { + + // 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<S> tree = this; tree.parent != null; tree = tree.parent) { + + // this is an hyperplane of some parent node + final Hyperplane<S> 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 = cut.split(hyperplane).getPlus(); + plus.chopOffMinus(hyperplane, vanishingHandler); + minus.chopOffMinus(hyperplane, vanishingHandler); + } else { + cut = cut.split(hyperplane).getMinus(); + plus.chopOffPlus(hyperplane, vanishingHandler); + minus.chopOffPlus(hyperplane, vanishingHandler); + } + + if (cut == null) { + // the cut sub-hyperplane has vanished + final BSPTree<S> fixed = vanishingHandler.fixNode(this); + cut = fixed.cut; + plus = fixed.plus; + minus = fixed.minus; + attribute = fixed.attribute; + if (cut == null) { + break; + } + } + + } + + // since we may have drop some parts of the inserted tree, + // perform a condensation pass to keep the tree structure simple + condense(); + + } + + } + + /** Prune a tree around a cell. + * <p> + * This method can be used to extract a convex cell from a tree. + * The original cell may either be a leaf node or an internal node. + * If it is an internal node, it's subtree will be ignored (i.e. the + * extracted cell will be a leaf node in all cases). The original + * tree to which the original cell belongs is not touched at all, + * a new independent tree will be built. + * </p> + * @param cellAttribute attribute to set for the leaf node + * corresponding to the initial instance cell + * @param otherLeafsAttributes attribute to set for the other leaf + * nodes + * @param internalAttributes attribute to set for the internal nodes + * @return a new tree (the original tree is left untouched) containing + * a single branch with the cell as a leaf node, and other leaf nodes + * as the remnants of the pruned branches + * @since 3.3 + */ + public BSPTree<S> pruneAroundConvexCell(final Object cellAttribute, + final Object otherLeafsAttributes, + final Object internalAttributes) { + + // build the current cell leaf + BSPTree<S> tree = new BSPTree<S>(cellAttribute); + + // build the pruned tree bottom-up + for (BSPTree<S> current = this; current.parent != null; current = current.parent) { + final SubHyperplane<S> parentCut = current.parent.cut.copySelf(); + final BSPTree<S> sibling = new BSPTree<S>(otherLeafsAttributes); + if (current == current.parent.plus) { + tree = new BSPTree<S>(parentCut, tree, sibling, internalAttributes); + } else { + tree = new BSPTree<S>(parentCut, sibling, tree, internalAttributes); + } + } + + return tree; + + } + + /** Chop off parts of the tree. + * <p>The instance is modified in place, all the parts that are on + * the minus side of the chopping hyperplane are discarded, only the + * parts on the plus side remain.</p> + * @param hyperplane chopping hyperplane + * @param vanishingHandler handler to use for handling very rare corner + * cases of vanishing cut sub-hyperplanes in internal nodes during merging + */ + private void chopOffMinus(final Hyperplane<S> hyperplane, final VanishingCutHandler<S> vanishingHandler) { + if (cut != null) { + + cut = cut.split(hyperplane).getPlus(); + plus.chopOffMinus(hyperplane, vanishingHandler); + minus.chopOffMinus(hyperplane, vanishingHandler); + + if (cut == null) { + // the cut sub-hyperplane has vanished + final BSPTree<S> fixed = vanishingHandler.fixNode(this); + cut = fixed.cut; + plus = fixed.plus; + minus = fixed.minus; + attribute = fixed.attribute; + } + + } + } + + /** Chop off parts of the tree. + * <p>The instance is modified in place, all the parts that are on + * the plus side of the chopping hyperplane are discarded, only the + * parts on the minus side remain.</p> + * @param hyperplane chopping hyperplane + * @param vanishingHandler handler to use for handling very rare corner + * cases of vanishing cut sub-hyperplanes in internal nodes during merging + */ + private void chopOffPlus(final Hyperplane<S> hyperplane, final VanishingCutHandler<S> vanishingHandler) { + if (cut != null) { + + cut = cut.split(hyperplane).getMinus(); + plus.chopOffPlus(hyperplane, vanishingHandler); + minus.chopOffPlus(hyperplane, vanishingHandler); + + if (cut == null) { + // the cut sub-hyperplane has vanished + final BSPTree<S> fixed = vanishingHandler.fixNode(this); + cut = fixed.cut; + plus = fixed.plus; + minus = fixed.minus; + attribute = fixed.attribute; + } + + } + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTreeVisitor.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTreeVisitor.java new file mode 100644 index 0000000..3d09939 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTreeVisitor.java @@ -0,0 +1,114 @@ +/* + * 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.math3.geometry.partitioning; + +import org.apache.commons.math3.geometry.Space; + +/** This interface is used to visit {@link BSPTree BSP tree} nodes. + + * <p>Navigation through {@link BSPTree BSP trees} can be done using + * two different point of views:</p> + * <ul> + * <li> + * 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 + * </li> + * <li> + * 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. + * </li> + * </ul> + + * @param <S> Type of the space. + + * @see BSPTree + * @see SubHyperplane + + * @since 3.0 + */ +public interface BSPTreeVisitor<S extends Space> { + + /** Enumerate for visit order with respect to plus sub-tree, minus sub-tree and cut sub-hyperplane. */ + enum Order { + /** Indicator for visit order plus sub-tree, then minus sub-tree, + * and last cut sub-hyperplane. + */ + PLUS_MINUS_SUB, + + /** Indicator for visit order plus sub-tree, then cut sub-hyperplane, + * and last minus sub-tree. + */ + PLUS_SUB_MINUS, + + /** Indicator for visit order minus sub-tree, then plus sub-tree, + * and last cut sub-hyperplane. + */ + MINUS_PLUS_SUB, + + /** Indicator for visit order minus sub-tree, then cut sub-hyperplane, + * and last plus sub-tree. + */ + MINUS_SUB_PLUS, + + /** Indicator for visit order cut sub-hyperplane, then plus sub-tree, + * and last minus sub-tree. + */ + SUB_PLUS_MINUS, + + /** Indicator for visit order cut sub-hyperplane, then minus sub-tree, + * and last plus sub-tree. + */ + SUB_MINUS_PLUS; + } + + /** Determine the visit order for this node. + * <p>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.</p> + * @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<S> node); + + /** Visit a BSP tree node node having a non-null sub-hyperplane. + * <p>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.</p> + * @param node BSP node guaranteed to have a non null cut sub-hyperplane + * @see #visitLeafNode + */ + void visitInternalNode(BSPTree<S> 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<S> node); + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryAttribute.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryAttribute.java new file mode 100644 index 0000000..dad884c --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryAttribute.java @@ -0,0 +1,116 @@ +/* + * 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.math3.geometry.partitioning; + +import org.apache.commons.math3.geometry.Space; + +/** Class holding boundary attributes. + * <p>This class is used for the attributes associated with the + * nodes of region boundary shell trees returned by the {@link + * Region#getTree(boolean) Region.getTree(includeBoundaryAttributes)} + * when the boolean {@code includeBoundaryAttributes} parameter is + * set to {@code true}. It contains the parts of the node cut + * sub-hyperplane that belong to the boundary.</p> + * <p>This class is a simple placeholder, it does not provide any + * processing methods.</p> + * @param <S> Type of the space. + * @see Region#getTree + * @since 3.0 + */ +public class BoundaryAttribute<S extends Space> { + + /** 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<S> 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<S> plusInside; + + /** Sub-hyperplanes that were used to split the boundary part. */ + private final NodesSet<S> splitters; + + /** 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) + * @deprecated as of 3.4, the constructor has been replaced by a new one + * which is not public anymore, as it is intended to be used only by + * {@link BoundaryBuilder} + */ + @Deprecated + public BoundaryAttribute(final SubHyperplane<S> plusOutside, + final SubHyperplane<S> plusInside) { + this(plusOutside, plusInside, null); + } + + /** 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) + * @param splitters sub-hyperplanes that were used to + * split the boundary part (may be null) + * @since 3.4 + */ + BoundaryAttribute(final SubHyperplane<S> plusOutside, + final SubHyperplane<S> plusInside, + final NodesSet<S> splitters) { + this.plusOutside = plusOutside; + this.plusInside = plusInside; + this.splitters = splitters; + } + + /** 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<S> 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<S> getPlusInside() { + return plusInside; + } + + /** Get the sub-hyperplanes that were used to split the boundary part. + * @return sub-hyperplanes that were used to split the boundary part + */ + public NodesSet<S> getSplitters() { + return splitters; + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryBuilder.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryBuilder.java new file mode 100644 index 0000000..cea4de3 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryBuilder.java @@ -0,0 +1,95 @@ +/* + * 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.math3.geometry.partitioning; + +import org.apache.commons.math3.geometry.Space; + +/** Visitor building boundary shell tree. + * <p> + * The boundary shell is represented as {@link BoundaryAttribute boundary attributes} + * at each internal node. + * </p> + * @param <S> Type of the space. + * @since 3.4 + */ +class BoundaryBuilder<S extends Space> implements BSPTreeVisitor<S> { + + /** {@inheritDoc} */ + public Order visitOrder(BSPTree<S> node) { + return Order.PLUS_MINUS_SUB; + } + + /** {@inheritDoc} */ + public void visitInternalNode(BSPTree<S> node) { + + SubHyperplane<S> plusOutside = null; + SubHyperplane<S> plusInside = null; + NodesSet<S> splitters = null; + + // characterize the cut sub-hyperplane, + // first with respect to the plus sub-tree + final Characterization<S> plusChar = new Characterization<S>(node.getPlus(), node.getCut().copySelf()); + + if (plusChar.touchOutside()) { + // plusChar.outsideTouching() 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<S> minusChar = new Characterization<S>(node.getMinus(), plusChar.outsideTouching()); + if (minusChar.touchInside()) { + // this part belongs to the boundary, + // it has the outside on its plus side and the inside on its minus side + plusOutside = minusChar.insideTouching(); + splitters = new NodesSet<S>(); + splitters.addAll(minusChar.getInsideSplitters()); + splitters.addAll(plusChar.getOutsideSplitters()); + } + } + + if (plusChar.touchInside()) { + // plusChar.insideTouching() 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<S> minusChar = new Characterization<S>(node.getMinus(), plusChar.insideTouching()); + if (minusChar.touchOutside()) { + // this part belongs to the boundary, + // it has the inside on its plus side and the outside on its minus side + plusInside = minusChar.outsideTouching(); + if (splitters == null) { + splitters = new NodesSet<S>(); + } + splitters.addAll(minusChar.getOutsideSplitters()); + splitters.addAll(plusChar.getInsideSplitters()); + } + } + + if (splitters != null) { + // the parent nodes are natural splitters for boundary sub-hyperplanes + for (BSPTree<S> up = node.getParent(); up != null; up = up.getParent()) { + splitters.add(up); + } + } + + // set the boundary attribute at non-leaf nodes + node.setAttribute(new BoundaryAttribute<S>(plusOutside, plusInside, splitters)); + + } + + /** {@inheritDoc} */ + public void visitLeafNode(BSPTree<S> node) { + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjection.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjection.java new file mode 100644 index 0000000..03526e4 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjection.java @@ -0,0 +1,83 @@ +/* + * 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.math3.geometry.partitioning; + +import org.apache.commons.math3.geometry.Point; +import org.apache.commons.math3.geometry.Space; + +/** Class holding the result of point projection on region boundary. + * <p>This class is a simple placeholder, it does not provide any + * processing methods.</p> + * <p>Instances of this class are guaranteed to be immutable</p> + * @param <S> Type of the space. + * @since 3.3 + * @see AbstractRegion#projectToBoundary(Point) + */ +public class BoundaryProjection<S extends Space> { + + /** Original point. */ + private final Point<S> original; + + /** Projected point. */ + private final Point<S> projected; + + /** Offset of the point with respect to the boundary it is projected on. */ + private final double offset; + + /** Constructor from raw elements. + * @param original original point + * @param projected projected point + * @param offset offset of the point with respect to the boundary it is projected on + */ + public BoundaryProjection(final Point<S> original, final Point<S> projected, final double offset) { + this.original = original; + this.projected = projected; + this.offset = offset; + } + + /** Get the original point. + * @return original point + */ + public Point<S> getOriginal() { + return original; + } + + /** Projected point. + * @return projected point, or null if there are no boundary + */ + public Point<S> getProjected() { + return projected; + } + + /** Offset of the point with respect to the boundary it is projected on. + * <p> + * The offset with respect to the boundary is negative if the {@link + * #getOriginal() original point} is inside the region, and positive otherwise. + * </p> + * <p> + * If there are no boundary, the value is set to either {@code + * Double.POSITIVE_INFINITY} if the region is empty (i.e. all points are + * outside of the region) or {@code Double.NEGATIVE_INFINITY} if the region + * covers the whole space (i.e. all points are inside of the region). + * </p> + * @return offset of the point with respect to the boundary it is projected on + */ + public double getOffset() { + return offset; + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjector.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjector.java new file mode 100644 index 0000000..486bbf1 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjector.java @@ -0,0 +1,200 @@ +/* + * 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.math3.geometry.partitioning; + +import java.util.ArrayList; +import java.util.List; + +import org.apache.commons.math3.geometry.Point; +import org.apache.commons.math3.geometry.Space; +import org.apache.commons.math3.geometry.partitioning.Region.Location; +import org.apache.commons.math3.util.FastMath; + +/** Local tree visitor to compute projection on boundary. + * @param <S> Type of the space. + * @param <T> Type of the sub-space. + * @since 3.3 + */ +class BoundaryProjector<S extends Space, T extends Space> implements BSPTreeVisitor<S> { + + /** Original point. */ + private final Point<S> original; + + /** Current best projected point. */ + private Point<S> projected; + + /** Leaf node closest to the test point. */ + private BSPTree<S> leaf; + + /** Current offset. */ + private double offset; + + /** Simple constructor. + * @param original original point + */ + BoundaryProjector(final Point<S> original) { + this.original = original; + this.projected = null; + this.leaf = null; + this.offset = Double.POSITIVE_INFINITY; + } + + /** {@inheritDoc} */ + public Order visitOrder(final BSPTree<S> node) { + // we want to visit the tree so that the first encountered + // leaf is the one closest to the test point + if (node.getCut().getHyperplane().getOffset(original) <= 0) { + return Order.MINUS_SUB_PLUS; + } else { + return Order.PLUS_SUB_MINUS; + } + } + + /** {@inheritDoc} */ + public void visitInternalNode(final BSPTree<S> node) { + + // project the point on the cut sub-hyperplane + final Hyperplane<S> hyperplane = node.getCut().getHyperplane(); + final double signedOffset = hyperplane.getOffset(original); + if (FastMath.abs(signedOffset) < offset) { + + // project point + final Point<S> regular = hyperplane.project(original); + + // get boundary parts + final List<Region<T>> boundaryParts = boundaryRegions(node); + + // check if regular projection really belongs to the boundary + boolean regularFound = false; + for (final Region<T> part : boundaryParts) { + if (!regularFound && belongsToPart(regular, hyperplane, part)) { + // the projected point lies in the boundary + projected = regular; + offset = FastMath.abs(signedOffset); + regularFound = true; + } + } + + if (!regularFound) { + // the regular projected point is not on boundary, + // so we have to check further if a singular point + // (i.e. a vertex in 2D case) is a possible projection + for (final Region<T> part : boundaryParts) { + final Point<S> spI = singularProjection(regular, hyperplane, part); + if (spI != null) { + final double distance = original.distance(spI); + if (distance < offset) { + projected = spI; + offset = distance; + } + } + } + + } + + } + + } + + /** {@inheritDoc} */ + public void visitLeafNode(final BSPTree<S> node) { + if (leaf == null) { + // this is the first leaf we visit, + // it is the closest one to the original point + leaf = node; + } + } + + /** Get the projection. + * @return projection + */ + public BoundaryProjection<S> getProjection() { + + // fix offset sign + offset = FastMath.copySign(offset, (Boolean) leaf.getAttribute() ? -1 : +1); + + return new BoundaryProjection<S>(original, projected, offset); + + } + + /** Extract the regions of the boundary on an internal node. + * @param node internal node + * @return regions in the node sub-hyperplane + */ + private List<Region<T>> boundaryRegions(final BSPTree<S> node) { + + final List<Region<T>> regions = new ArrayList<Region<T>>(2); + + @SuppressWarnings("unchecked") + final BoundaryAttribute<S> ba = (BoundaryAttribute<S>) node.getAttribute(); + addRegion(ba.getPlusInside(), regions); + addRegion(ba.getPlusOutside(), regions); + + return regions; + + } + + /** Add a boundary region to a list. + * @param sub sub-hyperplane defining the region + * @param list to fill up + */ + private void addRegion(final SubHyperplane<S> sub, final List<Region<T>> list) { + if (sub != null) { + @SuppressWarnings("unchecked") + final Region<T> region = ((AbstractSubHyperplane<S, T>) sub).getRemainingRegion(); + if (region != null) { + list.add(region); + } + } + } + + /** Check if a projected point lies on a boundary part. + * @param point projected point to check + * @param hyperplane hyperplane into which the point was projected + * @param part boundary part + * @return true if point lies on the boundary part + */ + private boolean belongsToPart(final Point<S> point, final Hyperplane<S> hyperplane, + final Region<T> part) { + + // there is a non-null sub-space, we can dive into smaller dimensions + @SuppressWarnings("unchecked") + final Embedding<S, T> embedding = (Embedding<S, T>) hyperplane; + return part.checkPoint(embedding.toSubSpace(point)) != Location.OUTSIDE; + + } + + /** Get the projection to the closest boundary singular point. + * @param point projected point to check + * @param hyperplane hyperplane into which the point was projected + * @param part boundary part + * @return projection to a singular point of boundary part (may be null) + */ + private Point<S> singularProjection(final Point<S> point, final Hyperplane<S> hyperplane, + final Region<T> part) { + + // there is a non-null sub-space, we can dive into smaller dimensions + @SuppressWarnings("unchecked") + final Embedding<S, T> embedding = (Embedding<S, T>) hyperplane; + final BoundaryProjection<T> bp = part.projectToBoundary(embedding.toSubSpace(point)); + + // back to initial dimension + return (bp.getProjected() == null) ? null : embedding.toSpace(bp.getProjected()); + + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundarySizeVisitor.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundarySizeVisitor.java new file mode 100644 index 0000000..054838b --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundarySizeVisitor.java @@ -0,0 +1,65 @@ +/* + * 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.math3.geometry.partitioning; + +import org.apache.commons.math3.geometry.Space; + +/** Visitor computing the boundary size. + * @param <S> Type of the space. + * @since 3.0 + */ +class BoundarySizeVisitor<S extends Space> implements BSPTreeVisitor<S> { + + /** Size of the boundary. */ + private double boundarySize; + + /** Simple constructor. + */ + BoundarySizeVisitor() { + boundarySize = 0; + } + + /** {@inheritDoc}*/ + public Order visitOrder(final BSPTree<S> node) { + return Order.MINUS_SUB_PLUS; + } + + /** {@inheritDoc}*/ + public void visitInternalNode(final BSPTree<S> node) { + @SuppressWarnings("unchecked") + final BoundaryAttribute<S> attribute = + (BoundaryAttribute<S>) node.getAttribute(); + if (attribute.getPlusOutside() != null) { + boundarySize += attribute.getPlusOutside().getSize(); + } + if (attribute.getPlusInside() != null) { + boundarySize += attribute.getPlusInside().getSize(); + } + } + + /** {@inheritDoc}*/ + public void visitLeafNode(final BSPTree<S> node) { + } + + /** Get the size of the boundary. + * @return size of the boundary + */ + public double getSize() { + return boundarySize; + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/Characterization.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/Characterization.java new file mode 100644 index 0000000..f8ec2f9 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/Characterization.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.math3.geometry.partitioning; + +import java.util.ArrayList; +import java.util.List; + +import org.apache.commons.math3.exception.MathInternalError; +import org.apache.commons.math3.geometry.Space; + +/** Cut sub-hyperplanes characterization with respect to inside/outside cells. + * @see BoundaryBuilder + * @param <S> Type of the space. + * @since 3.4 + */ +class Characterization<S extends Space> { + + /** Part of the cut sub-hyperplane that touch outside cells. */ + private SubHyperplane<S> outsideTouching; + + /** Part of the cut sub-hyperplane that touch inside cells. */ + private SubHyperplane<S> insideTouching; + + /** Nodes that were used to split the outside touching part. */ + private final NodesSet<S> outsideSplitters; + + /** Nodes that were used to split the outside touching part. */ + private final NodesSet<S> insideSplitters; + + /** Simple constructor. + * <p>Characterization consists in splitting the specified + * sub-hyperplane into several parts lying in inside and outside + * cells of the tree. The principle is to compute characterization + * twice for each cut sub-hyperplane in the tree, once on 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.</p> + * @param node current BSP tree node + * @param sub sub-hyperplane to characterize + */ + Characterization(final BSPTree<S> node, final SubHyperplane<S> sub) { + outsideTouching = null; + insideTouching = null; + outsideSplitters = new NodesSet<S>(); + insideSplitters = new NodesSet<S>(); + characterize(node, sub, new ArrayList<BSPTree<S>>()); + } + + /** Filter the parts of an hyperplane belonging to the boundary. + * <p>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 on 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.</p> + * @param node current BSP tree node + * @param sub sub-hyperplane to characterize + * @param splitters nodes that did split the current one + */ + private void characterize(final BSPTree<S> node, final SubHyperplane<S> sub, + final List<BSPTree<S>> splitters) { + if (node.getCut() == null) { + // we have reached a leaf node + final boolean inside = (Boolean) node.getAttribute(); + if (inside) { + addInsideTouching(sub, splitters); + } else { + addOutsideTouching(sub, splitters); + } + } else { + final Hyperplane<S> hyperplane = node.getCut().getHyperplane(); + final SubHyperplane.SplitSubHyperplane<S> split = sub.split(hyperplane); + switch (split.getSide()) { + case PLUS: + characterize(node.getPlus(), sub, splitters); + break; + case MINUS: + characterize(node.getMinus(), sub, splitters); + break; + case BOTH: + splitters.add(node); + characterize(node.getPlus(), split.getPlus(), splitters); + characterize(node.getMinus(), split.getMinus(), splitters); + splitters.remove(splitters.size() - 1); + break; + default: + // this should not happen + throw new MathInternalError(); + } + } + } + + /** Add a part of the cut sub-hyperplane known to touch an outside cell. + * @param sub part of the cut sub-hyperplane known to touch an outside cell + * @param splitters sub-hyperplanes that did split the current one + */ + private void addOutsideTouching(final SubHyperplane<S> sub, + final List<BSPTree<S>> splitters) { + if (outsideTouching == null) { + outsideTouching = sub; + } else { + outsideTouching = outsideTouching.reunite(sub); + } + outsideSplitters.addAll(splitters); + } + + /** Add a part of the cut sub-hyperplane known to touch an inside cell. + * @param sub part of the cut sub-hyperplane known to touch an inside cell + * @param splitters sub-hyperplanes that did split the current one + */ + private void addInsideTouching(final SubHyperplane<S> sub, + final List<BSPTree<S>> splitters) { + if (insideTouching == null) { + insideTouching = sub; + } else { + insideTouching = insideTouching.reunite(sub); + } + insideSplitters.addAll(splitters); + } + + /** Check if the cut sub-hyperplane touches outside cells. + * @return true if the cut sub-hyperplane touches outside cells + */ + public boolean touchOutside() { + return outsideTouching != null && !outsideTouching.isEmpty(); + } + + /** Get all the parts of the cut sub-hyperplane known to touch outside cells. + * @return parts of the cut sub-hyperplane known to touch outside cells + * (may be null or empty) + */ + public SubHyperplane<S> outsideTouching() { + return outsideTouching; + } + + /** Get the nodes that were used to split the outside touching part. + * <p> + * Splitting nodes are internal nodes (i.e. they have a non-null + * cut sub-hyperplane). + * </p> + * @return nodes that were used to split the outside touching part + */ + public NodesSet<S> getOutsideSplitters() { + return outsideSplitters; + } + + /** Check if the cut sub-hyperplane touches inside cells. + * @return true if the cut sub-hyperplane touches inside cells + */ + public boolean touchInside() { + return insideTouching != null && !insideTouching.isEmpty(); + } + + /** Get all the parts of the cut sub-hyperplane known to touch inside cells. + * @return parts of the cut sub-hyperplane known to touch inside cells + * (may be null or empty) + */ + public SubHyperplane<S> insideTouching() { + return insideTouching; + } + + /** Get the nodes that were used to split the inside touching part. + * <p> + * Splitting nodes are internal nodes (i.e. they have a non-null + * cut sub-hyperplane). + * </p> + * @return nodes that were used to split the inside touching part + */ + public NodesSet<S> getInsideSplitters() { + return insideSplitters; + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/Embedding.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/Embedding.java new file mode 100644 index 0000000..74e2c00 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/Embedding.java @@ -0,0 +1,68 @@ +/* + * 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.math3.geometry.partitioning; + +import org.apache.commons.math3.geometry.Point; +import org.apache.commons.math3.geometry.Space; + +/** This interface defines mappers between a space and one of its sub-spaces. + + * <p>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}. This interface can be used regardless + * of the dimensions differences. As an example, {@link + * org.apache.commons.math3.geometry.euclidean.threed.Line Line} in 3D + * implements Embedding<{@link + * org.apache.commons.math3.geometry.euclidean.threed.Vector3D Vector3D}, {link + * org.apache.commons.math3.geometry.euclidean.oned.Vector1D Vector1D>, i.e. it + * maps directly dimensions 3 and 1.</p> + + * <p>In the 3D euclidean space, hyperplanes are 2D planes, and the 1D + * sub-spaces are lines.</p> + + * <p> + * Note that this interface is <em>not</em> intended to be implemented + * by Apache Commons Math users, it is only intended to be implemented + * within the library itself. New methods may be added even for minor + * versions, which breaks compatibility for external implementations. + * </p> + + * @param <S> Type of the embedding space. + * @param <T> Type of the embedded sub-space. + + * @see Hyperplane + * @since 3.0 + */ +public interface Embedding<S extends Space, T extends Space> { + + /** 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<T> toSubSpace(Point<S> 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<S> toSpace(Point<T> point); + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/Hyperplane.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/Hyperplane.java new file mode 100644 index 0000000..f90c3bc --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/Hyperplane.java @@ -0,0 +1,98 @@ +/* + * 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.math3.geometry.partitioning; + +import org.apache.commons.math3.geometry.Point; +import org.apache.commons.math3.geometry.Space; + +/** This interface represents an hyperplane of a space. + + * <p>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.</p> + + * <p> + * Note that this interface is <em>not</em> intended to be implemented + * by Apache Commons Math users, it is only intended to be implemented + * within the library itself. New methods may be added even for minor + * versions, which breaks compatibility for external implementations. + * </p> + + * @param <S> Type of the space. + + * @since 3.0 + */ +public interface Hyperplane<S extends Space> { + + /** Copy the instance. + * <p>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).</p> + * @return a new hyperplane, copy of the instance + */ + Hyperplane<S> copySelf(); + + /** Get the offset (oriented distance) of a point. + * <p>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.</p> + * @param point point to check + * @return offset of the point + */ + double getOffset(Point<S> point); + + /** Project a point to the hyperplane. + * @param point point to project + * @return projected point + * @since 3.3 + */ + Point<S> project(Point<S> point); + + /** Get the tolerance below which points are considered to belong to the hyperplane. + * @return tolerance below which points are considered to belong to the hyperplane + * @since 3.3 + */ + double getTolerance(); + + /** Check if the instance has the same orientation as another hyperplane. + * <p>This method is expected to be called on parallel hyperplanes. The + * method should <em>not</em> re-check for parallelism, only for + * orientation, typically by testing something like the sign of the + * dot-products of normals.</p> + * @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<S> other); + + /** Build a sub-hyperplane covering the whole hyperplane. + * @return a sub-hyperplane covering the whole hyperplane + */ + SubHyperplane<S> wholeHyperplane(); + + /** Build a region covering the whole space. + * @return a region containing the instance + */ + Region<S> wholeSpace(); + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/InsideFinder.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/InsideFinder.java new file mode 100644 index 0000000..b1db90a --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/InsideFinder.java @@ -0,0 +1,150 @@ +/* + * 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.math3.geometry.partitioning; + +import org.apache.commons.math3.geometry.Space; + +/** Utility class checking if inside nodes can be found + * on the plus and minus sides of an hyperplane. + * @param <S> Type of the space. + * @since 3.4 + */ +class InsideFinder<S extends Space> { + + /** Region on which to operate. */ + private final Region<S> region; + + /** 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. + * @param region region on which to operate + */ + InsideFinder(final Region<S> region) { + this.region = region; + plusFound = false; + minusFound = false; + } + + /** Search recursively for inside leaf nodes on each side of the given hyperplane. + + * <p>The algorithm used here is directly derived from the one + * described in section III (<i>Binary Partitioning of a BSP + * Tree</i>) of the Bruce Naylor, John Amanatides and William + * Thibault paper <a + * href="http://www.cs.yorku.ca/~amana/research/bsptSetOp.pdf">Merging + * BSP Trees Yields Polyhedral Set Operations</a> Proc. Siggraph + * '90, Computer Graphics 24(4), August 1990, pp 115-124, published + * by the Association for Computing Machinery (ACM)..</p> + + * @param node current BSP tree node + * @param sub sub-hyperplane + */ + public void recurseSides(final BSPTree<S> node, final SubHyperplane<S> sub) { + + if (node.getCut() == null) { + if ((Boolean) node.getAttribute()) { + // this is an inside cell expanding across the hyperplane + plusFound = true; + minusFound = true; + } + return; + } + + final Hyperplane<S> hyperplane = node.getCut().getHyperplane(); + final SubHyperplane.SplitSubHyperplane<S> split = sub.split(hyperplane); + switch (split.getSide()) { + case PLUS : + // the sub-hyperplane is entirely in the plus sub-tree + if (node.getCut().split(sub.getHyperplane()).getSide() == Side.PLUS) { + if (!region.isEmpty(node.getMinus())) { + plusFound = true; + } + } else { + if (!region.isEmpty(node.getMinus())) { + minusFound = true; + } + } + if (!(plusFound && minusFound)) { + recurseSides(node.getPlus(), sub); + } + break; + case MINUS : + // the sub-hyperplane is entirely in the minus sub-tree + if (node.getCut().split(sub.getHyperplane()).getSide() == Side.PLUS) { + if (!region.isEmpty(node.getPlus())) { + plusFound = true; + } + } else { + if (!region.isEmpty(node.getPlus())) { + minusFound = true; + } + } + if (!(plusFound && minusFound)) { + recurseSides(node.getMinus(), sub); + } + break; + case BOTH : + // the sub-hyperplane extends in both sub-trees + + // explore first the plus sub-tree + recurseSides(node.getPlus(), split.getPlus()); + + // if needed, explore the minus sub-tree + if (!(plusFound && minusFound)) { + recurseSides(node.getMinus(), split.getMinus()); + } + 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) || ((Boolean) node.getPlus().getAttribute())) { + plusFound = true; + } + if ((node.getMinus().getCut() != null) || ((Boolean) node.getMinus().getAttribute())) { + minusFound = true; + } + } else { + if ((node.getPlus().getCut() != null) || ((Boolean) node.getPlus().getAttribute())) { + minusFound = true; + } + if ((node.getMinus().getCut() != null) || ((Boolean) node.getMinus().getAttribute())) { + 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; + } + + /** 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; + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/NodesSet.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/NodesSet.java new file mode 100644 index 0000000..688279a --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/NodesSet.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.math3.geometry.partitioning; + +import java.util.ArrayList; +import java.util.Iterator; +import java.util.List; + +import org.apache.commons.math3.geometry.Space; + +/** Set of {@link BSPTree BSP tree} nodes. + * @see BoundaryAttribute + * @param <S> Type of the space. + * @since 3.4 + */ +public class NodesSet<S extends Space> implements Iterable<BSPTree<S>> { + + /** List of sub-hyperplanes. */ + private List<BSPTree<S>> list; + + /** Simple constructor. + */ + public NodesSet() { + list = new ArrayList<BSPTree<S>>(); + } + + /** Add a node if not already known. + * @param node node to add + */ + public void add(final BSPTree<S> node) { + + for (final BSPTree<S> existing : list) { + if (node == existing) { + // the node is already known, don't add it + return; + } + } + + // the node was not known, add it + list.add(node); + + } + + /** Add nodes if they are not already known. + * @param iterator nodes iterator + */ + public void addAll(final Iterable<BSPTree<S>> iterator) { + for (final BSPTree<S> node : iterator) { + add(node); + } + } + + /** {@inheritDoc} */ + public Iterator<BSPTree<S>> iterator() { + return list.iterator(); + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/Region.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/Region.java new file mode 100644 index 0000000..9ff3946 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/Region.java @@ -0,0 +1,221 @@ +/* + * 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.math3.geometry.partitioning; + +import org.apache.commons.math3.geometry.Space; +import org.apache.commons.math3.geometry.Point; + +/** This interface represents a region of a space as a partition. + + * <p>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 <em>inside</em> the region from + * points considered to be <em>outside</em> of it. In between, there + * may be points on the <em>boundary</em> of the region.</p> + + * <p>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.</p> + + * <p>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.</p> + + * <p> + * Note that this interface is <em>not</em> intended to be implemented + * by Apache Commons Math users, it is only intended to be implemented + * within the library itself. New methods may be added even for minor + * versions, which breaks compatibility for external implementations. + * </p> + + * @param <S> Type of the space. + + * @since 3.0 + */ +public interface Region<S extends Space> { + + /** Enumerate for the location of a point with respect to the region. */ + 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; + } + + /** Build a region using the instance as a prototype. + * <p>This method allow to create new instances without knowing + * exactly the type of the region. It is an application of the + * prototype design pattern.</p> + * <p>The leaf nodes of the BSP tree <em>must</em> 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 <em>must</em> have either null internal nodes or + * internal nodes representing the boundary as specified in the + * {@link #getTree getTree} method).</p> + * @param newTree inside/outside BSP tree representing the new region + * @return the built region + */ + Region<S> buildNew(BSPTree<S> newTree); + + /** Copy the instance. + * <p>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 underlying tree {@code Boolean} + * attributes and immutable objects).</p> + * @return a new region, copy of the instance + */ + Region<S> copySelf(); + + /** Check if the instance is empty. + * @return true if the instance is empty + */ + boolean isEmpty(); + + /** Check if the sub-tree starting at a given node is empty. + * @param node root node of the sub-tree (<em>must</em> 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 + */ + boolean isEmpty(final BSPTree<S> node); + + /** Check if the instance covers the full space. + * @return true if the instance covers the full space + */ + boolean isFull(); + + /** Check if the sub-tree starting at a given node covers the full space. + * @param node root node of the sub-tree (<em>must</em> 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 covers the full space + */ + boolean isFull(final BSPTree<S> node); + + /** 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 + */ + boolean contains(final Region<S> region); + + /** 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} + */ + Location checkPoint(final Point<S> point); + + /** Project a point on the boundary of the region. + * @param point point to check + * @return projection of the point on the boundary + * @since 3.3 + */ + BoundaryProjection<S> projectToBoundary(final Point<S> point); + + /** Get the underlying BSP tree. + + * <p>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.</p> + + * <p>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#getPlusInside() getPlusInside} and {@link + * BoundaryAttribute#getPlusOutside() getPlusOutside} methods both + * returning null if the corresponding cut sub-hyperplane does not + * have any parts belonging to the boundary.</p> + + * <p>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.</p> + + * @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 + */ + BSPTree<S> getTree(final boolean includeBoundaryAttributes); + + /** 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 ...) + */ + double getBoundarySize(); + + /** 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 ...) + */ + double getSize(); + + /** Get the barycenter of the instance. + * @return an object representing the barycenter + */ + Point<S> getBarycenter(); + + /** Compute the relative position of the instance with respect to an + * hyperplane. + * @param hyperplane reference hyperplane + * @return one of {@link Side#PLUS Side.PLUS}, {@link Side#MINUS + * Side.MINUS}, {@link Side#BOTH Side.BOTH} or {@link Side#HYPER + * Side.HYPER} (the latter result can occur only if the tree + * contains only one cut hyperplane) + * @deprecated as of 3.6, this method which was only intended for + * internal use is not used anymore + */ + @Deprecated + Side side(final Hyperplane<S> hyperplane); + + /** Get the parts of a sub-hyperplane that are contained in the region. + * <p>The parts of the sub-hyperplane that belong to the boundary are + * <em>not</em> included in the resulting parts.</p> + * @param sub sub-hyperplane traversing the region + * @return filtered sub-hyperplane + */ + SubHyperplane<S> intersection(final SubHyperplane<S> sub); + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java new file mode 100644 index 0000000..688ffde --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java @@ -0,0 +1,378 @@ +/* + * 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.math3.geometry.partitioning; + +import java.util.HashMap; +import java.util.Map; + +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.geometry.Point; +import org.apache.commons.math3.geometry.Space; +import org.apache.commons.math3.geometry.partitioning.BSPTree.VanishingCutHandler; +import org.apache.commons.math3.geometry.partitioning.Region.Location; +import org.apache.commons.math3.geometry.partitioning.SubHyperplane.SplitSubHyperplane; + +/** This class is a factory for {@link Region}. + + * @param <S> Type of the space. + + * @since 3.0 + */ +public class RegionFactory<S extends Space> { + + /** Visitor removing internal nodes attributes. */ + private final NodesCleaner nodeCleaner; + + /** Simple constructor. + */ + public RegionFactory() { + nodeCleaner = new NodesCleaner(); + } + + /** Build a convex region from a collection of bounding hyperplanes. + * @param hyperplanes collection of bounding hyperplanes + * @return a new convex region, or null if the collection is empty + */ + public Region<S> buildConvex(final Hyperplane<S> ... hyperplanes) { + if ((hyperplanes == null) || (hyperplanes.length == 0)) { + return null; + } + + // use the first hyperplane to build the right class + final Region<S> region = hyperplanes[0].wholeSpace(); + + // chop off parts of the space + BSPTree<S> node = region.getTree(false); + node.setAttribute(Boolean.TRUE); + for (final Hyperplane<S> hyperplane : hyperplanes) { + if (node.insertCut(hyperplane)) { + node.setAttribute(null); + node.getPlus().setAttribute(Boolean.FALSE); + node = node.getMinus(); + node.setAttribute(Boolean.TRUE); + } else { + // the hyperplane could not be inserted in the current leaf node + // either it is completely outside (which means the input hyperplanes + // are wrong), or it is parallel to a previous hyperplane + SubHyperplane<S> s = hyperplane.wholeHyperplane(); + for (BSPTree<S> tree = node; tree.getParent() != null && s != null; tree = tree.getParent()) { + final Hyperplane<S> other = tree.getParent().getCut().getHyperplane(); + final SplitSubHyperplane<S> split = s.split(other); + switch (split.getSide()) { + case HYPER : + // the hyperplane is parallel to a previous hyperplane + if (!hyperplane.sameOrientationAs(other)) { + // this hyperplane is opposite to the other one, + // the region is thinner than the tolerance, we consider it empty + return getComplement(hyperplanes[0].wholeSpace()); + } + // the hyperplane is an extension of an already known hyperplane, we just ignore it + break; + case PLUS : + // the hyperplane is outside of the current convex zone, + // the input hyperplanes are inconsistent + throw new MathIllegalArgumentException(LocalizedFormats.NOT_CONVEX_HYPERPLANES); + default : + s = split.getMinus(); + } + } + } + } + + return region; + + } + + /** 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 Region<S> union(final Region<S> region1, final Region<S> region2) { + final BSPTree<S> tree = + region1.getTree(false).merge(region2.getTree(false), new UnionMerger()); + tree.visit(nodeCleaner); + 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 Region<S> intersection(final Region<S> region1, final Region<S> region2) { + final BSPTree<S> tree = + region1.getTree(false).merge(region2.getTree(false), new IntersectionMerger()); + tree.visit(nodeCleaner); + 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 Region<S> xor(final Region<S> region1, final Region<S> region2) { + final BSPTree<S> tree = + region1.getTree(false).merge(region2.getTree(false), new XorMerger()); + tree.visit(nodeCleaner); + 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 Region<S> difference(final Region<S> region1, final Region<S> region2) { + final BSPTree<S> tree = + region1.getTree(false).merge(region2.getTree(false), new DifferenceMerger(region1, region2)); + tree.visit(nodeCleaner); + return region1.buildNew(tree); + } + + /** Get the complement of the region (exchanged interior/exterior). + * @param region region to complement, it will not modified, a new + * region independent region will be built + * @return a new region, complement of the specified one + */ + /** Get the complement of the region (exchanged interior/exterior). + * @param region region to complement, it will not modified, a new + * region independent region will be built + * @return a new region, complement of the specified one + */ + public Region<S> getComplement(final Region<S> region) { + return region.buildNew(recurseComplement(region.getTree(false))); + } + + /** Recursively build the complement of a BSP tree. + * @param node current node of the original tree + * @return new tree, complement of the node + */ + private BSPTree<S> recurseComplement(final BSPTree<S> node) { + + // transform the tree, except for boundary attribute splitters + final Map<BSPTree<S>, BSPTree<S>> map = new HashMap<BSPTree<S>, BSPTree<S>>(); + final BSPTree<S> transformedTree = recurseComplement(node, map); + + // set up the boundary attributes splitters + for (final Map.Entry<BSPTree<S>, BSPTree<S>> entry : map.entrySet()) { + if (entry.getKey().getCut() != null) { + @SuppressWarnings("unchecked") + BoundaryAttribute<S> original = (BoundaryAttribute<S>) entry.getKey().getAttribute(); + if (original != null) { + @SuppressWarnings("unchecked") + BoundaryAttribute<S> transformed = (BoundaryAttribute<S>) entry.getValue().getAttribute(); + for (final BSPTree<S> splitter : original.getSplitters()) { + transformed.getSplitters().add(map.get(splitter)); + } + } + } + } + + return transformedTree; + + } + + /** Recursively build the complement of a BSP tree. + * @param node current node of the original tree + * @param map transformed nodes map + * @return new tree, complement of the node + */ + private BSPTree<S> recurseComplement(final BSPTree<S> node, + final Map<BSPTree<S>, BSPTree<S>> map) { + + final BSPTree<S> transformedNode; + if (node.getCut() == null) { + transformedNode = new BSPTree<S>(((Boolean) node.getAttribute()) ? Boolean.FALSE : Boolean.TRUE); + } else { + + @SuppressWarnings("unchecked") + BoundaryAttribute<S> attribute = (BoundaryAttribute<S>) node.getAttribute(); + if (attribute != null) { + final SubHyperplane<S> plusOutside = + (attribute.getPlusInside() == null) ? null : attribute.getPlusInside().copySelf(); + final SubHyperplane<S> plusInside = + (attribute.getPlusOutside() == null) ? null : attribute.getPlusOutside().copySelf(); + // we start with an empty list of splitters, it will be filled in out of recursion + attribute = new BoundaryAttribute<S>(plusOutside, plusInside, new NodesSet<S>()); + } + + transformedNode = new BSPTree<S>(node.getCut().copySelf(), + recurseComplement(node.getPlus(), map), + recurseComplement(node.getMinus(), map), + attribute); + } + + map.put(node, transformedNode); + return transformedNode; + + } + + /** BSP tree leaf merger computing union of two regions. */ + private class UnionMerger implements BSPTree.LeafMerger<S> { + /** {@inheritDoc} */ + public BSPTree<S> merge(final BSPTree<S> leaf, final BSPTree<S> tree, + final BSPTree<S> parentTree, + final boolean isPlusChild, final boolean leafFromInstance) { + if ((Boolean) leaf.getAttribute()) { + // the leaf node represents an inside cell + leaf.insertInTree(parentTree, isPlusChild, new VanishingToLeaf(true)); + return leaf; + } + // the leaf node represents an outside cell + tree.insertInTree(parentTree, isPlusChild, new VanishingToLeaf(false)); + return tree; + } + } + + /** BSP tree leaf merger computing intersection of two regions. */ + private class IntersectionMerger implements BSPTree.LeafMerger<S> { + /** {@inheritDoc} */ + public BSPTree<S> merge(final BSPTree<S> leaf, final BSPTree<S> tree, + final BSPTree<S> parentTree, + final boolean isPlusChild, final boolean leafFromInstance) { + if ((Boolean) leaf.getAttribute()) { + // the leaf node represents an inside cell + tree.insertInTree(parentTree, isPlusChild, new VanishingToLeaf(true)); + return tree; + } + // the leaf node represents an outside cell + leaf.insertInTree(parentTree, isPlusChild, new VanishingToLeaf(false)); + return leaf; + } + } + + /** BSP tree leaf merger computing symmetric difference (exclusive or) of two regions. */ + private class XorMerger implements BSPTree.LeafMerger<S> { + /** {@inheritDoc} */ + public BSPTree<S> merge(final BSPTree<S> leaf, final BSPTree<S> tree, + final BSPTree<S> parentTree, final boolean isPlusChild, + final boolean leafFromInstance) { + BSPTree<S> t = tree; + if ((Boolean) leaf.getAttribute()) { + // the leaf node represents an inside cell + t = recurseComplement(t); + } + t.insertInTree(parentTree, isPlusChild, new VanishingToLeaf(true)); + return t; + } + } + + /** BSP tree leaf merger computing difference of two regions. */ + private class DifferenceMerger implements BSPTree.LeafMerger<S>, VanishingCutHandler<S> { + + /** Region to subtract from. */ + private final Region<S> region1; + + /** Region to subtract. */ + private final Region<S> region2; + + /** Simple constructor. + * @param region1 region to subtract from + * @param region2 region to subtract + */ + DifferenceMerger(final Region<S> region1, final Region<S> region2) { + this.region1 = region1.copySelf(); + this.region2 = region2.copySelf(); + } + + /** {@inheritDoc} */ + public BSPTree<S> merge(final BSPTree<S> leaf, final BSPTree<S> tree, + final BSPTree<S> parentTree, final boolean isPlusChild, + final boolean leafFromInstance) { + if ((Boolean) leaf.getAttribute()) { + // the leaf node represents an inside cell + final BSPTree<S> argTree = + recurseComplement(leafFromInstance ? tree : leaf); + argTree.insertInTree(parentTree, isPlusChild, this); + return argTree; + } + // the leaf node represents an outside cell + final BSPTree<S> instanceTree = + leafFromInstance ? leaf : tree; + instanceTree.insertInTree(parentTree, isPlusChild, this); + return instanceTree; + } + + /** {@inheritDoc} */ + public BSPTree<S> fixNode(final BSPTree<S> node) { + // get a representative point in the degenerate cell + final BSPTree<S> cell = node.pruneAroundConvexCell(Boolean.TRUE, Boolean.FALSE, null); + final Region<S> r = region1.buildNew(cell); + final Point<S> p = r.getBarycenter(); + return new BSPTree<S>(region1.checkPoint(p) == Location.INSIDE && + region2.checkPoint(p) == Location.OUTSIDE); + } + + } + + /** Visitor removing internal nodes attributes. */ + private class NodesCleaner implements BSPTreeVisitor<S> { + + /** {@inheritDoc} */ + public Order visitOrder(final BSPTree<S> node) { + return Order.PLUS_SUB_MINUS; + } + + /** {@inheritDoc} */ + public void visitInternalNode(final BSPTree<S> node) { + node.setAttribute(null); + } + + /** {@inheritDoc} */ + public void visitLeafNode(final BSPTree<S> node) { + } + + } + + /** Handler replacing nodes with vanishing cuts with leaf nodes. */ + private class VanishingToLeaf implements VanishingCutHandler<S> { + + /** Inside/outside indocator to use for ambiguous nodes. */ + private final boolean inside; + + /** Simple constructor. + * @param inside inside/outside indicator to use for ambiguous nodes + */ + VanishingToLeaf(final boolean inside) { + this.inside = inside; + } + + /** {@inheritDoc} */ + public BSPTree<S> fixNode(final BSPTree<S> node) { + if (node.getPlus().getAttribute().equals(node.getMinus().getAttribute())) { + // no ambiguity + return new BSPTree<S>(node.getPlus().getAttribute()); + } else { + // ambiguous node + return new BSPTree<S>(inside); + } + } + + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/Side.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/Side.java new file mode 100644 index 0000000..c9a1357 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/Side.java @@ -0,0 +1,37 @@ +/* + * 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.math3.geometry.partitioning; + +/** Enumerate representing the location of an element with respect to an + * {@link Hyperplane hyperplane} of a space. + * @since 3.0 + */ +public 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; + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/SubHyperplane.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/SubHyperplane.java new file mode 100644 index 0000000..2069f6f --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/SubHyperplane.java @@ -0,0 +1,155 @@ +/* + * 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.math3.geometry.partitioning; + +import org.apache.commons.math3.geometry.Space; + +/** This interface represents the remaining parts of an hyperplane after + * other parts have been chopped off. + + * <p>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.</p> + + * <p> + * Note that this interface is <em>not</em> intended to be implemented + * by Apache Commons Math users, it is only intended to be implemented + * within the library itself. New methods may be added even for minor + * versions, which breaks compatibility for external implementations. + * </p> + + * @param <S> Type of the embedding space. + + * @since 3.0 + */ +public interface SubHyperplane<S extends Space> { + + /** Copy the instance. + * <p>The instance created is completely independent of the original + * one. A deep copy is used, none of the underlying objects are + * shared (except for the nodes attributes and immutable + * objects).</p> + * @return a new sub-hyperplane, copy of the instance + */ + SubHyperplane<S> copySelf(); + + /** Get the underlying hyperplane. + * @return underlying hyperplane + */ + Hyperplane<S> getHyperplane(); + + /** Check if the instance is empty. + * @return true if the instance is empty + */ + boolean isEmpty(); + + /** 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 ...) + */ + double getSize(); + + /** Compute the relative position of the instance with respect + * to an hyperplane. + * @param hyperplane hyperplane to check instance against + * @return one of {@link Side#PLUS}, {@link Side#MINUS}, {@link Side#BOTH}, + * {@link Side#HYPER} + * @deprecated as of 3.6, replaced with {@link #split(Hyperplane)}.{@link SplitSubHyperplane#getSide()} + */ + @Deprecated + Side side(Hyperplane<S> hyperplane); + + /** Split the instance in two parts by an hyperplane. + * @param hyperplane splitting hyperplane + * @return an object containing both the part of the instance + * on the plus side of the hyperplane and the part of the + * instance on the minus side of the hyperplane + */ + SplitSubHyperplane<S> split(Hyperplane<S> hyperplane); + + /** Compute the union of the instance and another sub-hyperplane. + * @param other other sub-hyperplane to union (<em>must</em> be in the + * same hyperplane as the instance) + * @return a new sub-hyperplane, union of the instance and other + */ + SubHyperplane<S> reunite(SubHyperplane<S> other); + + /** Class holding the results of the {@link #split split} method. + * @param <U> Type of the embedding space. + */ + class SplitSubHyperplane<U extends Space> { + + /** Part of the sub-hyperplane on the plus side of the splitting hyperplane. */ + private final SubHyperplane<U> plus; + + /** Part of the sub-hyperplane on the minus side of the splitting hyperplane. */ + private final SubHyperplane<U> 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<U> plus, + final SubHyperplane<U> 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<U> 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<U> getMinus() { + return minus; + } + + /** Get the side of the split sub-hyperplane with respect to its splitter. + * @return {@link Side#PLUS} if only {@link #getPlus()} is neither null nor empty, + * {@link Side#MINUS} if only {@link #getMinus()} is neither null nor empty, + * {@link Side#BOTH} if both {@link #getPlus()} and {@link #getMinus()} + * are neither null nor empty or {@link Side#HYPER} if both {@link #getPlus()} and + * {@link #getMinus()} are either null or empty + * @since 3.6 + */ + public Side getSide() { + if (plus != null && !plus.isEmpty()) { + if (minus != null && !minus.isEmpty()) { + return Side.BOTH; + } else { + return Side.PLUS; + } + } else if (minus != null && !minus.isEmpty()) { + return Side.MINUS; + } else { + return Side.HYPER; + } + } + + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/Transform.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/Transform.java new file mode 100644 index 0000000..ba0c1dd --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/Transform.java @@ -0,0 +1,80 @@ +/* + * 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.math3.geometry.partitioning; + +import org.apache.commons.math3.geometry.Point; +import org.apache.commons.math3.geometry.Space; + + +/** This interface represents an inversible affine transform in a space. + * <p>Inversible affine transform include for example scalings, + * translations, rotations.</p> + + * <p>Transforms are dimension-specific. The consistency rules between + * the three {@code apply} methods are the following ones for a + * transformed defined for dimension D:</p> + * <ul> + * <li> + * the transform can be applied to a point in the + * D-dimension space using its {@link #apply(Point)} + * method + * </li> + * <li> + * the transform can be applied to a (D-1)-dimension + * hyperplane in the D-dimension space using its + * {@link #apply(Hyperplane)} method + * </li> + * <li> + * 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 + * </li> + * </ul> + + * @param <S> Type of the embedding space. + * @param <T> Type of the embedded sub-space. + + * @since 3.0 + */ +public interface Transform<S extends Space, T extends Space> { + + /** Transform a point of a space. + * @param point point to transform + * @return a new object representing the transformed point + */ + Point<S> apply(Point<S> point); + + /** Transform an hyperplane of a space. + * @param hyperplane hyperplane to transform + * @return a new object representing the transformed hyperplane + */ + Hyperplane<S> apply(Hyperplane<S> hyperplane); + + /** Transform a sub-hyperplane embedded in an hyperplane. + * @param sub sub-hyperplane to transform + * @param original hyperplane in which the sub-hyperplane is + * defined (this is the original hyperplane, the transform has + * <em>not</em> been applied to it) + * @param transformed hyperplane in which the sub-hyperplane is + * defined (this is the transformed hyperplane, the transform + * <em>has</em> been applied to it) + * @return a new object representing the transformed sub-hyperplane + */ + SubHyperplane<T> apply(SubHyperplane<T> sub, Hyperplane<S> original, Hyperplane<S> transformed); + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/package-info.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/package-info.java new file mode 100644 index 0000000..6e63c73 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/package-info.java @@ -0,0 +1,114 @@ +/* + * 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. + */ +/** + * + * This package provides classes to implement Binary Space Partition trees. + * + * <p> + * {@link org.apache.commons.math3.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). + * </p> + * + * <p> + * 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. + * </p> + * + * <p> + * 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 <em>not</em> 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.math3.geometry.partitioning.Region Region} class is + * devoted to this representation, it is build on top of the {@link + * org.apache.commons.math3.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). + * </p> + * + * <p> + * Rather than simply associating the internal nodes with an hyperplane, + * we consider <em>sub-hyperplanes</em> 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. + * </p> + * + * <p> + * 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 + * <code>O(log(n))</code> where <code>n</code> is the number of facets) + * or if the first tree levels close to the root discriminate large parts + * of the total space. + * </p> + * + * <p> + * One of the main sources for the development of this package was Bruce + * Naylor, John Amanatides and William Thibault paper <a + * href="http://www.cs.yorku.ca/~amana/research/bsptSetOp.pdf">Merging + * BSP Trees Yields Polyhedral Set Operations</a> 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 <a + * href="http://www.cs.utexas.edu/users/fussell/courses/cs384g/bsp_treemerge.pdf">here</a>. + * </p> + * + * <p> + * Note that the interfaces defined in this package are <em>not</em> intended to + * be implemented by Apache Commons Math users, they are only intended to be + * implemented within the library itself. New methods may be added even for + * minor versions, which breaks compatibility for external implementations. + * </p> + * + */ +package org.apache.commons.math3.geometry.partitioning; diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/AVLTree.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/AVLTree.java new file mode 100644 index 0000000..00c9d3e --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/AVLTree.java @@ -0,0 +1,634 @@ +/* + * 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.math3.geometry.partitioning.utilities; + +/** This class implements AVL trees. + * + * <p>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.</p> + * + * <p>Since the {@code equals} method is not sufficient to + * differentiate elements, the {@link #delete delete} method is + * implemented using the equality operator.</p> + * + * <p>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}).</p> + * + * <p>This class is based on the C implementation Georg Kraml has put + * in the public domain. Unfortunately, his <a + * href="www.purists.org/georg/avltree/index.html">page</a> seems not + * to exist any more.</p> + * + * @param <T> the type of the elements + * + * @since 3.0 + * @deprecated as of 3.4, this class is not used anymore and considered + * to be out of scope of Apache Commons Math + */ +@Deprecated +public class AVLTree<T extends Comparable<T>> { + + /** Top level node. */ + private Node top; + + /** Build an empty tree. + */ + public AVLTree() { + top = null; + } + + /** Insert an element in the tree. + * @param element element to insert (silently ignored if null) + */ + public void insert(final T element) { + if (element != null) { + if (top == null) { + top = new Node(element, null); + } else { + top.insert(element); + } + } + } + + /** Delete an element from the tree. + * <p>The 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 + * <em>different</em> 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).</p> + * @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 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. + * <p>AVL tree nodes implement all the logical structure of the + * tree. Nodes are created by the {@link AVLTree AVLTree} class.</p> + * <p>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.</p> + * @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. + * <p>The skew factor are not updated by this method, they + * <em>must</em> be updated by the caller</p> + */ + 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. + * <p>The skew factor are not updated by this method, they + * <em>must</em> be updated by the caller</p> + */ + 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/math3/geometry/partitioning/utilities/OrderedTuple.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/OrderedTuple.java new file mode 100644 index 0000000..2dad2d7 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/OrderedTuple.java @@ -0,0 +1,431 @@ +/* + * 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.math3.geometry.partitioning.utilities; + +import java.util.Arrays; + +import org.apache.commons.math3.util.FastMath; + +/** This class implements an ordering operation for T-uples. + * + * <p>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 2<sup>k</sup> bit of component 0 is + * followed by the 2<sup>k</sup> bit of component 1, then by the + * 2<sup>k</sup> bit of component 2 up to the 2<sup>k</sup> bit of + * component {@code t}, which is followed by the 2<sup>k-1</sup> + * bit of component 0, followed by the 2<sup>k-1</sup> bit of + * component 1 ... The binary representations are extended as needed + * to handle numbers with different scales and a suitable + * 2<sup>p</sup> offset is added to the components in order to avoid + * negative numbers (this offset is adjusted as needed during the + * comparison operations).</p> + * + * <p>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 + * picture:</p> + * + * <img src="doc-files/OrderedTuple.png" /> + * + * <p>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.</p> + * + * <p>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.</p> + * + * <p>T-uples with negative infinite or positive infinite components + * are sorted logically.</p> + * + * <p>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 <em>regular</em> and consistent T-uples + * together.</p> + * <ul> + * <li>instances with different dimensions are sorted according to + * their dimension regardless of their components values</li> + * <li>instances with {@code Double.NaN} components are sorted + * after all other ones (even after instances with positive infinite + * components</li> + * <li>instances with both positive and negative infinite components + * are considered as if they had {@code Double.NaN} + * components</li> + * </ul> + * + * @since 3.0 + * @deprecated as of 3.4, this class is not used anymore and considered + * to be out of scope of Apache Commons Math + */ +@Deprecated +public class OrderedTuple implements Comparable<OrderedTuple> { + + /** Sign bit mask. */ + private static final long SIGN_MASK = 0x8000000000000000L; + + /** Exponent bits mask. */ + private static final long EXPONENT_MASK = 0x7ff0000000000000L; + + /** Mantissa bits mask. */ + private static final long MANTISSA_MASK = 0x000fffffffffffffL; + + /** Implicit MSB for normalized numbers. */ + private static final long IMPLICIT_ONE = 0x0010000000000000L; + + /** Double components of the T-uple. */ + private double[] components; + + /** Offset scale. */ + private int offset; + + /** Least Significant Bit scale. */ + private int lsb; + + /** Ordering encoding of the double components. */ + private long[] encoding; + + /** Positive infinity marker. */ + private boolean posInf; + + /** Negative infinity marker. */ + private boolean negInf; + + /** Not A Number marker. */ + private boolean nan; + + /** Build an ordered T-uple from its components. + * @param components double components of the T-uple + */ + public OrderedTuple(final double ... components) { + this.components = components.clone(); + int msb = Integer.MIN_VALUE; + lsb = Integer.MAX_VALUE; + posInf = false; + negInf = false; + nan = false; + for (int i = 0; i < components.length; ++i) { + if (Double.isInfinite(components[i])) { + if (components[i] < 0) { + negInf = true; + } else { + posInf = true; + } + } else if (Double.isNaN(components[i])) { + nan = true; + } else { + final long b = Double.doubleToLongBits(components[i]); + final long m = mantissa(b); + if (m != 0) { + final int e = exponent(b); + msb = FastMath.max(msb, e + computeMSB(m)); + lsb = FastMath.min(lsb, e + computeLSB(m)); + } + } + } + + if (posInf && negInf) { + // instance cannot be sorted logically + posInf = false; + negInf = false; + nan = true; + } + + if (lsb <= msb) { + // encode the T-upple with the specified offset + encode(msb + 16); + } else { + encoding = new long[] { + 0x0L + }; + } + + } + + /** Encode the T-uple with a given offset. + * @param minOffset minimal scale of the offset to add to all + * components (must be greater than the MSBs of all components) + */ + private void encode(final int minOffset) { + + // choose an offset with some margins + offset = minOffset + 31; + offset -= offset % 32; + + if ((encoding != null) && (encoding.length == 1) && (encoding[0] == 0x0L)) { + // the components are all zeroes + return; + } + + // allocate an integer array to encode the components (we use only + // 63 bits per element because there is no unsigned long in Java) + final int neededBits = offset + 1 - lsb; + final int neededLongs = (neededBits + 62) / 63; + encoding = new long[components.length * neededLongs]; + + // mix the bits from all components + int eIndex = 0; + int shift = 62; + long word = 0x0L; + for (int k = offset; eIndex < encoding.length; --k) { + for (int vIndex = 0; vIndex < components.length; ++vIndex) { + if (getBit(vIndex, k) != 0) { + word |= 0x1L << shift; + } + if (shift-- == 0) { + encoding[eIndex++] = word; + word = 0x0L; + shift = 62; + } + } + } + + } + + /** Compares this ordered T-uple with the specified object. + + * <p>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.</p> + + * <p>T-uples with negative infinite, positive infinite are sorted + * logically.</p> + + * <p>Some arbitrary choices have been made to handle specific + * cases. The rationale for these choices is to keep + * <em>normal</em> and consistent T-uples together.</p> + * <ul> + * <li>instances with different dimensions are sorted according to + * their dimension regardless of their components values</li> + * <li>instances with {@code Double.NaN} components are sorted + * after all other ones (evan after instances with positive infinite + * components</li> + * <li>instances with both positive and negative infinite components + * are considered as if they had {@code Double.NaN} + * components</li> + * </ul> + + * @param ot T-uple to compare instance with + * @return a negative integer if the instance is less than the + * object, zero if they are equal, or a positive integer if the + * instance is greater than the object + + */ + public int compareTo(final OrderedTuple ot) { + if (components.length == ot.components.length) { + if (nan) { + return +1; + } else if (ot.nan) { + return -1; + } else if (negInf || ot.posInf) { + return -1; + } else if (posInf || ot.negInf) { + return +1; + } else { + + if (offset < ot.offset) { + encode(ot.offset); + } else if (offset > ot.offset) { + ot.encode(offset); + } + + final int limit = FastMath.min(encoding.length, ot.encoding.length); + for (int i = 0; i < limit; ++i) { + if (encoding[i] < ot.encoding[i]) { + return -1; + } else if (encoding[i] > ot.encoding[i]) { + return +1; + } + } + + if (encoding.length < ot.encoding.length) { + return -1; + } else if (encoding.length > ot.encoding.length) { + return +1; + } else { + return 0; + } + + } + } + + return components.length - ot.components.length; + + } + + /** {@inheritDoc} */ + @Override + public boolean equals(final Object other) { + if (this == other) { + return true; + } else if (other instanceof OrderedTuple) { + return compareTo((OrderedTuple) other) == 0; + } else { + return false; + } + } + + /** {@inheritDoc} */ + @Override + public int hashCode() { + // the following constants are arbitrary small primes + final int multiplier = 37; + final int trueHash = 97; + final int falseHash = 71; + + // hash fields and combine them + // (we rely on the multiplier to have different combined weights + // for all int fields and all boolean fields) + int hash = Arrays.hashCode(components); + hash = hash * multiplier + offset; + hash = hash * multiplier + lsb; + hash = hash * multiplier + (posInf ? trueHash : falseHash); + hash = hash * multiplier + (negInf ? trueHash : falseHash); + hash = hash * multiplier + (nan ? trueHash : falseHash); + + return hash; + + } + + /** Get the components array. + * @return array containing the T-uple components + */ + public double[] getComponents() { + return components.clone(); + } + + /** Extract the sign from the bits of a double. + * @param bits binary representation of the double + * @return sign bit (zero if positive, non zero if negative) + */ + private static long sign(final long bits) { + return bits & SIGN_MASK; + } + + /** Extract the exponent from the bits of a double. + * @param bits binary representation of the double + * @return exponent + */ + private static int exponent(final long bits) { + return ((int) ((bits & EXPONENT_MASK) >> 52)) - 1075; + } + + /** Extract the mantissa from the bits of a double. + * @param bits binary representation of the double + * @return mantissa + */ + private static long mantissa(final long bits) { + return ((bits & EXPONENT_MASK) == 0) ? + ((bits & MANTISSA_MASK) << 1) : // subnormal number + (IMPLICIT_ONE | (bits & MANTISSA_MASK)); // normal number + } + + /** Compute the most significant bit of a long. + * @param l long from which the most significant bit is requested + * @return scale of the most significant bit of {@code l}, + * or 0 if {@code l} is zero + * @see #computeLSB + */ + private static int computeMSB(final long l) { + + long ll = l; + long mask = 0xffffffffL; + int scale = 32; + int msb = 0; + + while (scale != 0) { + if ((ll & mask) != ll) { + msb |= scale; + ll >>= scale; + } + scale >>= 1; + mask >>= scale; + } + + return msb; + + } + + /** Compute the least significant bit of a long. + * @param l long from which the least significant bit is requested + * @return scale of the least significant bit of {@code l}, + * or 63 if {@code l} is zero + * @see #computeMSB + */ + private static int computeLSB(final long l) { + + long ll = l; + long mask = 0xffffffff00000000L; + int scale = 32; + int lsb = 0; + + while (scale != 0) { + if ((ll & mask) == ll) { + lsb |= scale; + ll >>= scale; + } + scale >>= 1; + mask >>= scale; + } + + return lsb; + + } + + /** Get a bit from the mantissa of a double. + * @param i index of the component + * @param k scale of the requested bit + * @return the specified bit (either 0 or 1), after the offset has + * been added to the double + */ + private int getBit(final int i, final int k) { + final long bits = Double.doubleToLongBits(components[i]); + final int e = exponent(bits); + if ((k < e) || (k > offset)) { + return 0; + } else if (k == offset) { + return (sign(bits) == 0L) ? 1 : 0; + } else if (k > (e + 52)) { + return (sign(bits) == 0L) ? 0 : 1; + } else { + final long m = (sign(bits) == 0L) ? mantissa(bits) : -mantissa(bits); + return (int) ((m >> (k - e)) & 0x1L); + } + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/doc-files/OrderedTuple.png b/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/doc-files/OrderedTuple.png Binary files differnew file mode 100644 index 0000000..4eca233 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/doc-files/OrderedTuple.png diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/package-info.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/package-info.java new file mode 100644 index 0000000..31f57f1 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/package-info.java @@ -0,0 +1,24 @@ +/* + * 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. + */ +/** + * + * <p> + * This package provides multidimensional ordering features for partitioning. + * </p> + * + */ +package org.apache.commons.math3.geometry.partitioning.utilities; |
