You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2014/01/21 18:49:13 UTC

svn commit: r1560115 - in /commons/proper/math/trunk/src: changes/ main/java/org/apache/commons/math3/geometry/euclidean/oned/ main/java/org/apache/commons/math3/geometry/euclidean/threed/ main/java/org/apache/commons/math3/geometry/euclidean/twod/ mai...

Author: luc
Date: Tue Jan 21 17:49:13 2014
New Revision: 1560115

URL: http://svn.apache.org/r1560115
Log:
BSP tree now provide an API to compute a global signed distance.

The distance is defined from a test point to the region. It is positive
if the point is outside of the region, negative if the point is inside,
and zero when the point is at the boundary. The distance is continuous
everywhere, so it can be used with a root solver to identify accurately
boundary crossings. This API is available for all BSP trees, in
Euclidean and spherical geometries, and in all dimensions.

Fixes MATH-1091

Added:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjection.java   (with props)
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjector.java   (with props)
Modified:
    commons/proper/math/trunk/src/changes/changes.xml
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/oned/IntervalsSet.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/oned/OrientedPoint.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/Plane.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/Hyperplane.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/Region.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/LimitAngle.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Circle.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSetTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSetTest.java

Modified: commons/proper/math/trunk/src/changes/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/changes/changes.xml?rev=1560115&r1=1560114&r2=1560115&view=diff
==============================================================================
--- commons/proper/math/trunk/src/changes/changes.xml (original)
+++ commons/proper/math/trunk/src/changes/changes.xml Tue Jan 21 17:49:13 2014
@@ -51,6 +51,15 @@ If the output is not quite correct, chec
   </properties>
   <body>
     <release version="3.3" date="TBD" description="TBD">
+      <action dev="luc" type="add" issue="MATH-1091">
+       BSP tree now provide an API to compute a global signed distance from
+       a test point to the region. The distance is positive if the point is
+       outside of the region, negative if the point is inside, and zero
+       when the point is at the boundary. The distance is continuous
+       everywhere, so it can be used with a root solver to identify accurately
+       boundary crossings. This API is available for all BSP trees, in
+       Euclidean and spherical geometries, and in all dimensions.
+      </action>
       <action dev="luc" type="add" issue="MATH-1090">
        IntervalsSet now implements Iterable&lt;double[]&gt;, so one can iterate
        over the sub-intervals without building a full list containing

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/oned/IntervalsSet.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/oned/IntervalsSet.java?rev=1560115&r1=1560114&r2=1560115&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/oned/IntervalsSet.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/oned/IntervalsSet.java Tue Jan 21 17:49:13 2014
@@ -25,6 +25,7 @@ import java.util.NoSuchElementException;
 import org.apache.commons.math3.geometry.Point;
 import org.apache.commons.math3.geometry.partitioning.AbstractRegion;
 import org.apache.commons.math3.geometry.partitioning.BSPTree;
+import org.apache.commons.math3.geometry.partitioning.BoundaryProjection;
 import org.apache.commons.math3.geometry.partitioning.SubHyperplane;
 import org.apache.commons.math3.util.Precision;
 
@@ -268,6 +269,53 @@ public class IntervalsSet extends Abstra
         return ((Boolean) node.getAttribute()) ? Double.POSITIVE_INFINITY : sup;
     }
 
+    /** {@inheritDoc}
+     * @since 3.3
+     */
+    public BoundaryProjection<Euclidean1D> projectToBoundary(final Point<Euclidean1D> point) {
+
+        // get position of test point
+        final double x = ((Vector1D) point).getX();
+
+        double previous = Double.NEGATIVE_INFINITY;
+        for (final double[] a : this) {
+            if (x < a[0]) {
+                // the test point lies between the previous and the current intervals
+                // offset will be positive
+                final double previousOffset = x - previous;
+                final double currentOffset  = a[0] - x;
+                if (previousOffset < currentOffset) {
+                    return new BoundaryProjection<Euclidean1D>(point, finiteOrNullPoint(previous), previousOffset);
+                } else {
+                    return new BoundaryProjection<Euclidean1D>(point, finiteOrNullPoint(a[0]), currentOffset);
+                }
+            } else if (x <= a[1]) {
+                // the test point lies within the current interval
+                // offset will be negative
+                final double offset0 = a[0] - x;
+                final double offset1 = x - a[1];
+                if (offset0 < offset1) {
+                    return new BoundaryProjection<Euclidean1D>(point, finiteOrNullPoint(a[1]), offset1);
+                } else {
+                    return new BoundaryProjection<Euclidean1D>(point, finiteOrNullPoint(a[0]), offset0);
+                }
+            }
+            previous = a[1];
+        }
+
+        // the test point if past the last sub-interval
+        return new BoundaryProjection<Euclidean1D>(point, finiteOrNullPoint(previous), x - previous);
+
+    }
+
+    /** Build a finite point.
+     * @param x abscissa of the point
+     * @return a new point for finite abscissa, null otherwise
+     */
+    private Vector1D finiteOrNullPoint(final double x) {
+        return Double.isInfinite(x) ? null : new Vector1D(x);
+    }
+
     /** Build an ordered list of intervals representing the instance.
      * <p>This method builds this intervals set as an ordered list of
      * {@link Interval Interval} elements. If the intervals set has no

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/oned/OrientedPoint.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/oned/OrientedPoint.java?rev=1560115&r1=1560114&r2=1560115&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/oned/OrientedPoint.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/oned/OrientedPoint.java Tue Jan 21 17:49:13 2014
@@ -116,7 +116,16 @@ public class OrientedPoint implements Hy
         return !(direct ^ ((OrientedPoint) other).direct);
     }
 
-    /** {@inheritDoc} */
+    /** {@inheritDoc}
+     * @since 3.3
+     */
+    public Point<Euclidean1D> project(Point<Euclidean1D> point) {
+        return location;
+    }
+
+    /** {@inheritDoc}
+     * @since 3.3
+     */
     public double getTolerance() {
         return tolerance;
     }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/Plane.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/Plane.java?rev=1560115&r1=1560114&r2=1560115&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/Plane.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/Plane.java Tue Jan 21 17:49:13 2014
@@ -252,7 +252,16 @@ public class Plane implements Hyperplane
         return v;
     }
 
-    /** {@inheritDoc} */
+    /** {@inheritDoc}
+     * @since 3.3
+     */
+    public Point<Euclidean3D> project(Point<Euclidean3D> point) {
+        return toSpace(toSubSpace(point));
+    }
+
+    /** {@inheritDoc}
+     * @since 3.3
+     */
     public double getTolerance() {
         return tolerance;
     }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java?rev=1560115&r1=1560114&r2=1560115&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java Tue Jan 21 17:49:13 2014
@@ -260,7 +260,16 @@ public class Line implements Hyperplane<
                             (sin * other.originOffset - other.sin * originOffset) / d);
     }
 
-    /** {@inheritDoc} */
+    /** {@inheritDoc}
+     * @since 3.3
+     */
+    public Point<Euclidean2D> project(Point<Euclidean2D> point) {
+        return toSpace(toSubSpace(point));
+    }
+
+    /** {@inheritDoc}
+     * @since 3.3
+     */
     public double getTolerance() {
         return tolerance;
     }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java?rev=1560115&r1=1560114&r2=1560115&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java Tue Jan 21 17:49:13 2014
@@ -273,6 +273,15 @@ public abstract class AbstractRegion<S e
         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

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java?rev=1560115&r1=1560114&r2=1560115&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BSPTree.java Tue Jan 21 17:49:13 2014
@@ -682,6 +682,7 @@ public class BSPTree<S extends Space> {
      * @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,

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjection.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjection.java?rev=1560115&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjection.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjection.java Tue Jan 21 17:49:13 2014
@@ -0,0 +1,84 @@
+/*
+ * 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.
+ * @version $Id$
+ * @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> get0riginal() {
+        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;
+    }
+
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjection.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjection.java
------------------------------------------------------------------------------
    svn:keywords = "Author Date Id Revision"

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjector.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjector.java?rev=1560115&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjector.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjector.java Tue Jan 21 17:49:13 2014
@@ -0,0 +1,201 @@
+/*
+ * 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.
+ * @version $Id$
+ * @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
+     */
+    public 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());
+
+    }
+
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjector.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryProjector.java
------------------------------------------------------------------------------
    svn:keywords = "Author Date Id Revision"

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/Hyperplane.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/Hyperplane.java?rev=1560115&r1=1560114&r2=1560115&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/Hyperplane.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/Hyperplane.java Tue Jan 21 17:49:13 2014
@@ -55,6 +55,13 @@ public interface Hyperplane<S extends Sp
      */
     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

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/Region.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/Region.java?rev=1560115&r1=1560114&r2=1560115&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/Region.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/partitioning/Region.java Tue Jan 21 17:49:13 2014
@@ -112,6 +112,13 @@ public interface Region<S extends Space>
      */
     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

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java?rev=1560115&r1=1560114&r2=1560115&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/ArcsSet.java Tue Jan 21 17:49:13 2014
@@ -26,8 +26,10 @@ import org.apache.commons.math3.exceptio
 import org.apache.commons.math3.exception.MathInternalError;
 import org.apache.commons.math3.exception.NumberIsTooLargeException;
 import org.apache.commons.math3.exception.util.LocalizedFormats;
+import org.apache.commons.math3.geometry.Point;
 import org.apache.commons.math3.geometry.partitioning.AbstractRegion;
 import org.apache.commons.math3.geometry.partitioning.BSPTree;
+import org.apache.commons.math3.geometry.partitioning.BoundaryProjection;
 import org.apache.commons.math3.geometry.partitioning.Side;
 import org.apache.commons.math3.geometry.partitioning.SubHyperplane;
 import org.apache.commons.math3.util.FastMath;
@@ -473,6 +475,88 @@ public class ArcsSet extends AbstractReg
         }
     }
 
+    /** {@inheritDoc}
+     * @since 3.3
+     */
+    public BoundaryProjection<Sphere1D> projectToBoundary(final Point<Sphere1D> point) {
+
+        // get position of test point
+        final double alpha = ((S1Point) point).getAlpha();
+
+        boolean wrapFirst = false;
+        double first      = Double.NaN;
+        double previous   = Double.NaN;
+        for (final double[] a : this) {
+
+            if (Double.isNaN(first)) {
+                // remember the first angle in case we need it later
+                first = a[0];
+            }
+
+            if (!wrapFirst) {
+                if (alpha < a[0]) {
+                    // the test point lies between the previous and the current arcs
+                    // offset will be positive
+                    if (Double.isNaN(previous)) {
+                        // we need to wrap around the circle
+                        wrapFirst = true;
+                    } else {
+                        final double previousOffset = alpha - previous;
+                        final double currentOffset  = a[0] - alpha;
+                        if (previousOffset < currentOffset) {
+                            return new BoundaryProjection<Sphere1D>(point, new S1Point(previous), previousOffset);
+                        } else {
+                            return new BoundaryProjection<Sphere1D>(point, new S1Point(a[0]), currentOffset);
+                        }
+                    }
+                } else if (alpha <= a[1]) {
+                    // the test point lies within the current arc
+                    // offset will be negative
+                    final double offset0 = a[0] - alpha;
+                    final double offset1 = alpha - a[1];
+                    if (offset0 < offset1) {
+                        return new BoundaryProjection<Sphere1D>(point, new S1Point(a[1]), offset1);
+                    } else {
+                        return new BoundaryProjection<Sphere1D>(point, new S1Point(a[0]), offset0);
+                    }
+                }
+            }
+            previous = a[1];
+        }
+
+        if (Double.isNaN(previous)) {
+
+            // there are no points at all in the arcs set
+            return new BoundaryProjection<Sphere1D>(point, null, MathUtils.TWO_PI);
+
+        } else {
+
+            // the test point if before first arc and after last arc,
+            // somewhere around the 0/2 \pi crossing
+            if (wrapFirst) {
+                // the test point is between 0 and first
+                final double previousOffset = alpha - (previous - MathUtils.TWO_PI);
+                final double currentOffset  = first - alpha;
+                if (previousOffset < currentOffset) {
+                    return new BoundaryProjection<Sphere1D>(point, new S1Point(previous), previousOffset);
+                } else {
+                    return new BoundaryProjection<Sphere1D>(point, new S1Point(first), currentOffset);
+                }
+            } else {
+                // the test point is between last and 2\pi
+                final double previousOffset = alpha - previous;
+                final double currentOffset  = first + MathUtils.TWO_PI - alpha;
+                if (previousOffset < currentOffset) {
+                    return new BoundaryProjection<Sphere1D>(point, new S1Point(previous), previousOffset);
+                } else {
+                    return new BoundaryProjection<Sphere1D>(point, new S1Point(first), currentOffset);
+                }
+            }
+
+        }
+
+    }
+
     /** Build an ordered list of arcs representing the instance.
      * <p>This method builds this arcs set as an ordered list of
      * {@link Arc Arc} elements. An empty tree will build an empty list

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/LimitAngle.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/LimitAngle.java?rev=1560115&r1=1560114&r2=1560115&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/LimitAngle.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/oned/LimitAngle.java Tue Jan 21 17:49:13 2014
@@ -115,9 +115,12 @@ public class LimitAngle implements Hyper
         return location;
     }
 
-    /** Get the tolerance below which angles are considered identical.
-     * @return tolerance below which angles are considered identical
-     */
+    /** {@inheritDoc} */
+    public Point<Sphere1D> project(Point<Sphere1D> point) {
+        return location;
+    }
+
+    /** {@inheritDoc} */
     public double getTolerance() {
         return tolerance;
     }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Circle.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Circle.java?rev=1560115&r1=1560114&r2=1560115&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Circle.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/Circle.java Tue Jan 21 17:49:13 2014
@@ -134,9 +134,12 @@ public class Circle implements Hyperplan
         return new Circle(pole.negate(), x, y.negate(), tolerance);
     }
 
-    /** Get the tolerance below which close sub-arcs are merged together.
-     * @return tolerance below which close sub-arcs are merged together
-     */
+    /** {@inheritDoc} */
+    public Point<Sphere2D> project(Point<Sphere2D> point) {
+        return toSpace(toSubSpace(point));
+    }
+
+    /** {@inheritDoc} */
     public double getTolerance() {
         return tolerance;
     }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java?rev=1560115&r1=1560114&r2=1560115&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSet.java Tue Jan 21 17:49:13 2014
@@ -818,8 +818,15 @@ public class SphericalPolygonsSet extend
             }
 
             if (following == null) {
+                final Vector3D previousStart = previous.getStart().getLocation().getVector();
+                if (Vector3D.angle(point.getVector(), previousStart) <= tolerance) {
+                    // the edge connects back to itself
+                    return previous;
+                }
+
                 // this should never happen
                 throw new MathIllegalStateException(LocalizedFormats.OUTLINE_BOUNDARY_LOOP_OPEN);
+
             }
 
             return following;

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSetTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSetTest.java?rev=1560115&r1=1560114&r2=1560115&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSetTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSetTest.java Tue Jan 21 17:49:13 2014
@@ -23,6 +23,7 @@ import org.apache.commons.math3.geometry
 import org.apache.commons.math3.geometry.euclidean.oned.IntervalsSet;
 import org.apache.commons.math3.geometry.euclidean.oned.Vector1D;
 import org.apache.commons.math3.geometry.partitioning.BSPTree;
+import org.apache.commons.math3.geometry.partitioning.BoundaryProjection;
 import org.apache.commons.math3.geometry.partitioning.Region;
 import org.apache.commons.math3.geometry.partitioning.Region.Location;
 import org.apache.commons.math3.geometry.partitioning.RegionFactory;
@@ -100,6 +101,38 @@ public class PolygonsSetTest {
     }
 
     @Test
+    public void testEmpty() {
+        PolygonsSet empty = (PolygonsSet) new RegionFactory<Euclidean2D>().getComplement(new PolygonsSet(1.0e-10));
+        Assert.assertTrue(empty.isEmpty());
+        Assert.assertEquals(0, empty.getVertices().length);
+        Assert.assertEquals(0.0, empty.getBoundarySize(), 1.0e-10);
+        Assert.assertEquals(0.0, empty.getSize(), 1.0e-10);
+        for (double y = -1; y < 1; y += 0.1) {
+            for (double x = -1; x < 1; x += 0.1) {
+                Assert.assertEquals(Double.POSITIVE_INFINITY,
+                                    empty.projectToBoundary(new Vector2D(x, y)).getOffset(),
+                                    1.0e-10);
+            }
+        }
+    }
+
+    @Test
+    public void testFull() {
+        PolygonsSet empty = new PolygonsSet(1.0e-10);
+        Assert.assertFalse(empty.isEmpty());
+        Assert.assertEquals(0, empty.getVertices().length);
+        Assert.assertEquals(0.0, empty.getBoundarySize(), 1.0e-10);
+        Assert.assertEquals(Double.POSITIVE_INFINITY, empty.getSize(), 1.0e-10);
+        for (double y = -1; y < 1; y += 0.1) {
+            for (double x = -1; x < 1; x += 0.1) {
+                Assert.assertEquals(Double.NEGATIVE_INFINITY,
+                                    empty.projectToBoundary(new Vector2D(x, y)).getOffset(),
+                                    1.0e-10);
+            }
+        }
+    }
+
+    @Test
     public void testHole() {
         Vector2D[][] vertices = new Vector2D[][] {
             new Vector2D[] {
@@ -141,6 +174,40 @@ public class PolygonsSetTest {
             new Vector2D(3.0, 3.0)
         });
         checkVertices(set.getVertices(), vertices);
+
+        for (double x = -0.999; x < 3.999; x += 0.11) {
+            Vector2D v = new Vector2D(x, x + 0.5);
+            BoundaryProjection<Euclidean2D> projection = set.projectToBoundary(v);
+            Assert.assertTrue(projection.get0riginal() == v);
+            Vector2D p = (Vector2D) projection.getProjected();
+            if (x < -0.5) {
+                Assert.assertEquals(0.0,      p.getX(), 1.0e-10);
+                Assert.assertEquals(0.0,      p.getY(), 1.0e-10);
+                Assert.assertEquals(+v.distance(Vector2D.ZERO), projection.getOffset(), 1.0e-10);
+            } else if (x < 0.5) {
+                Assert.assertEquals(0.0,      p.getX(), 1.0e-10);
+                Assert.assertEquals(v.getY(), p.getY(), 1.0e-10);
+                Assert.assertEquals(-v.getX(), projection.getOffset(), 1.0e-10);
+            } else if (x < 1.25) {
+                Assert.assertEquals(1.0,      p.getX(), 1.0e-10);
+                Assert.assertEquals(v.getY(), p.getY(), 1.0e-10);
+                Assert.assertEquals(v.getX() - 1.0, projection.getOffset(), 1.0e-10);
+            } else if (x < 2.0) {
+                Assert.assertEquals(v.getX(), p.getX(), 1.0e-10);
+                Assert.assertEquals(2.0,      p.getY(), 1.0e-10);
+                Assert.assertEquals(2.0 - v.getY(), projection.getOffset(), 1.0e-10);
+            } else if (x < 3.0) {
+                Assert.assertEquals(v.getX(), p.getX(), 1.0e-10);
+                Assert.assertEquals(3.0,      p.getY(), 1.0e-10);
+                Assert.assertEquals(v.getY() - 3.0, projection.getOffset(), 1.0e-10);
+            } else {
+                Assert.assertEquals(3.0,      p.getX(), 1.0e-10);
+                Assert.assertEquals(3.0,      p.getY(), 1.0e-10);
+                Assert.assertEquals(+v.distance(new Vector2D(3, 3)), projection.getOffset(), 1.0e-10);
+            }
+            
+        }
+
     }
 
     @Test

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSetTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSetTest.java?rev=1560115&r1=1560114&r2=1560115&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSetTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/spherical/twod/SphericalPolygonsSetTest.java Tue Jan 21 17:49:13 2014
@@ -46,6 +46,7 @@ public class SphericalPolygonsSetTest {
         }
         Assert.assertEquals(4 * FastMath.PI, new SphericalPolygonsSet(0.01, new S2Point[0]).getSize(), 1.0e-10);
         Assert.assertEquals(0, new SphericalPolygonsSet(0.01, new S2Point[0]).getBoundarySize(), 1.0e-10);
+        Assert.assertEquals(0, full.getBoundaryLoops().size());
     }
 
     @Test
@@ -65,6 +66,7 @@ public class SphericalPolygonsSetTest {
                 Assert.assertEquals(Location.BOUNDARY, south.checkPoint(new S2Point(v)));
             }
         }
+        Assert.assertEquals(1, south.getBoundaryLoops().size());
     }
 
     @Test
@@ -312,16 +314,16 @@ public class SphericalPolygonsSetTest {
         double tol = 0.01;
         double alpha = 0.7;
         S2Point center = new S2Point(new Vector3D(1, 1, 1));
-        SphericalPolygonsSet octant = new SphericalPolygonsSet(center.getVector(), Vector3D.PLUS_K, alpha, 6, tol);
-        SphericalPolygonsSet hole   = new SphericalPolygonsSet(tol,
-                                                               new S2Point(FastMath.PI / 6, FastMath.PI / 3),
-                                                               new S2Point(FastMath.PI / 3, FastMath.PI / 3),
-                                                               new S2Point(FastMath.PI / 4, FastMath.PI / 6));
-        SphericalPolygonsSet octantWithHole =
-                (SphericalPolygonsSet) new RegionFactory<Sphere2D>().difference(octant, hole);
+        SphericalPolygonsSet hexa = new SphericalPolygonsSet(center.getVector(), Vector3D.PLUS_K, alpha, 6, tol);
+        SphericalPolygonsSet hole  = new SphericalPolygonsSet(tol,
+                                                              new S2Point(FastMath.PI / 6, FastMath.PI / 3),
+                                                              new S2Point(FastMath.PI / 3, FastMath.PI / 3),
+                                                              new S2Point(FastMath.PI / 4, FastMath.PI / 6));
+        SphericalPolygonsSet hexaWithHole =
+                (SphericalPolygonsSet) new RegionFactory<Sphere2D>().difference(hexa, hole);
 
         for (double phi = center.getPhi() - alpha + 0.1; phi < center.getPhi() + alpha - 0.1; phi += 0.07) {
-            Location l = octantWithHole.checkPoint(new S2Point(FastMath.PI / 4, phi));
+            Location l = hexaWithHole.checkPoint(new S2Point(FastMath.PI / 4, phi));
             if (phi < FastMath.PI / 6 || phi > FastMath.PI / 3) {
                 Assert.assertEquals(Location.INSIDE,  l);
             } else {
@@ -330,10 +332,10 @@ public class SphericalPolygonsSetTest {
         }
 
         // there should be two separate boundary loops
-        Assert.assertEquals(2, octantWithHole.getBoundaryLoops().size());
+        Assert.assertEquals(2, hexaWithHole.getBoundaryLoops().size());
 
-        Assert.assertEquals(octant.getBoundarySize() + hole.getBoundarySize(), octantWithHole.getBoundarySize(), 1.0e-10);
-        Assert.assertEquals(octant.getSize() - hole.getSize(), octantWithHole.getSize(), 1.0e-10);
+        Assert.assertEquals(hexa.getBoundarySize() + hole.getBoundarySize(), hexaWithHole.getBoundarySize(), 1.0e-10);
+        Assert.assertEquals(hexa.getSize() - hole.getSize(), hexaWithHole.getSize(), 1.0e-10);
 
     }
 
@@ -372,6 +374,25 @@ public class SphericalPolygonsSetTest {
                             triOut.getSize()    - triIn.getSize(),
                             concentric.getSize(), 1.0e-10);
 
+        // we expect lots of sign changes as we traverse all concentric rings
+        double phi = new S2Point(center).getPhi();
+        Assert.assertEquals(+0.207, concentric.projectToBoundary(new S2Point(-0.60,  phi)).getOffset(), 0.01);
+        Assert.assertEquals(-0.048, concentric.projectToBoundary(new S2Point(-0.21,  phi)).getOffset(), 0.01);
+        Assert.assertEquals(+0.027, concentric.projectToBoundary(new S2Point(-0.10,  phi)).getOffset(), 0.01);
+        Assert.assertEquals(-0.041, concentric.projectToBoundary(new S2Point( 0.01,  phi)).getOffset(), 0.01);
+        Assert.assertEquals(+0.049, concentric.projectToBoundary(new S2Point( 0.16,  phi)).getOffset(), 0.01);
+        Assert.assertEquals(-0.038, concentric.projectToBoundary(new S2Point( 0.29,  phi)).getOffset(), 0.01);
+        Assert.assertEquals(+0.097, concentric.projectToBoundary(new S2Point( 0.48,  phi)).getOffset(), 0.01);
+        Assert.assertEquals(-0.022, concentric.projectToBoundary(new S2Point( 0.64,  phi)).getOffset(), 0.01);
+        Assert.assertEquals(+0.072, concentric.projectToBoundary(new S2Point( 0.79,  phi)).getOffset(), 0.01);
+        Assert.assertEquals(-0.022, concentric.projectToBoundary(new S2Point( 0.93,  phi)).getOffset(), 0.01);
+        Assert.assertEquals(+0.091, concentric.projectToBoundary(new S2Point( 1.08,  phi)).getOffset(), 0.01);
+        Assert.assertEquals(-0.037, concentric.projectToBoundary(new S2Point( 1.28,  phi)).getOffset(), 0.01);
+        Assert.assertEquals(+0.051, concentric.projectToBoundary(new S2Point( 1.40,  phi)).getOffset(), 0.01);
+        Assert.assertEquals(-0.041, concentric.projectToBoundary(new S2Point( 1.55,  phi)).getOffset(), 0.01);
+        Assert.assertEquals(+0.027, concentric.projectToBoundary(new S2Point( 1.67,  phi)).getOffset(), 0.01);
+        Assert.assertEquals(-0.044, concentric.projectToBoundary(new S2Point( 1.79,  phi)).getOffset(), 0.01);
+        Assert.assertEquals(+0.201, concentric.projectToBoundary(new S2Point( 2.16,  phi)).getOffset(), 0.01);
     }
 
     private SubCircle create(Vector3D pole, Vector3D x, Vector3D y,