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 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1560115 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
5b71c17c81
commit
43ef84b729
|
|
@ -51,6 +51,15 @@ If the output is not quite correct, check for invisible trailing spaces!
|
|||
</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<double[]>, so one can iterate
|
||||
over the sub-intervals without building a full list containing
|
||||
|
|
|
|||
|
|
@ -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 AbstractRegion<Euclidean1D, Euclidean1D> imple
|
|||
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
|
||||
|
|
|
|||
|
|
@ -116,7 +116,16 @@ public class OrientedPoint implements Hyperplane<Euclidean1D> {
|
|||
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;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -252,7 +252,16 @@ public class Plane implements Hyperplane<Euclidean3D>, Embedding<Euclidean3D, Eu
|
|||
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;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -260,7 +260,16 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
|
|||
(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;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -273,6 +273,15 @@ public abstract class AbstractRegion<S extends Space, T extends Space> implement
|
|||
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
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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());
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -55,6 +55,13 @@ public interface Hyperplane<S extends Space> {
|
|||
*/
|
||||
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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -26,8 +26,10 @@ import org.apache.commons.math3.exception.MathIllegalArgumentException;
|
|||
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 AbstractRegion<Sphere1D, Sphere1D> implements Itera
|
|||
}
|
||||
}
|
||||
|
||||
/** {@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
|
||||
|
|
|
|||
|
|
@ -115,9 +115,12 @@ public class LimitAngle implements Hyperplane<Sphere1D> {
|
|||
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;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -134,9 +134,12 @@ public class Circle implements Hyperplane<Sphere2D>, Embedding<Sphere2D, Sphere1
|
|||
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;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -818,8 +818,15 @@ public class SphericalPolygonsSet extends AbstractRegion<Sphere2D, Sphere1D> {
|
|||
}
|
||||
|
||||
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;
|
||||
|
|
|
|||
|
|
@ -23,6 +23,7 @@ import org.apache.commons.math3.geometry.euclidean.oned.Interval;
|
|||
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;
|
||||
|
|
@ -99,6 +100,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[][] {
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
Loading…
Reference in New Issue