diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/geometry/euclidean/twod')
16 files changed, 3762 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/DiskGenerator.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/DiskGenerator.java new file mode 100644 index 0000000..332b1b7 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/DiskGenerator.java @@ -0,0 +1,108 @@ +/* + * 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.euclidean.twod; + +import java.util.List; + +import org.apache.commons.math3.fraction.BigFraction; +import org.apache.commons.math3.geometry.enclosing.EnclosingBall; +import org.apache.commons.math3.geometry.enclosing.SupportBallGenerator; +import org.apache.commons.math3.util.FastMath; + +/** Class generating an enclosing ball from its support points. + * @since 3.3 + */ +public class DiskGenerator implements SupportBallGenerator<Euclidean2D, Vector2D> { + + /** {@inheritDoc} */ + public EnclosingBall<Euclidean2D, Vector2D> ballOnSupport(final List<Vector2D> support) { + + if (support.size() < 1) { + return new EnclosingBall<Euclidean2D, Vector2D>(Vector2D.ZERO, Double.NEGATIVE_INFINITY); + } else { + final Vector2D vA = support.get(0); + if (support.size() < 2) { + return new EnclosingBall<Euclidean2D, Vector2D>(vA, 0, vA); + } else { + final Vector2D vB = support.get(1); + if (support.size() < 3) { + return new EnclosingBall<Euclidean2D, Vector2D>(new Vector2D(0.5, vA, 0.5, vB), + 0.5 * vA.distance(vB), + vA, vB); + } else { + final Vector2D vC = support.get(2); + // a disk is 2D can be defined as: + // (1) (x - x_0)^2 + (y - y_0)^2 = r^2 + // which can be written: + // (2) (x^2 + y^2) - 2 x_0 x - 2 y_0 y + (x_0^2 + y_0^2 - r^2) = 0 + // or simply: + // (3) (x^2 + y^2) + a x + b y + c = 0 + // with disk center coordinates -a/2, -b/2 + // If the disk exists, a, b and c are a non-zero solution to + // [ (x^2 + y^2 ) x y 1 ] [ 1 ] [ 0 ] + // [ (xA^2 + yA^2) xA yA 1 ] [ a ] [ 0 ] + // [ (xB^2 + yB^2) xB yB 1 ] * [ b ] = [ 0 ] + // [ (xC^2 + yC^2) xC yC 1 ] [ c ] [ 0 ] + // So the determinant of the matrix is zero. Computing this determinant + // by expanding it using the minors m_ij of first row leads to + // (4) m_11 (x^2 + y^2) - m_12 x + m_13 y - m_14 = 0 + // So by identifying equations (2) and (4) we get the coordinates + // of center as: + // x_0 = +m_12 / (2 m_11) + // y_0 = -m_13 / (2 m_11) + // Note that the minors m_11, m_12 and m_13 all have the last column + // filled with 1.0, hence simplifying the computation + final BigFraction[] c2 = new BigFraction[] { + new BigFraction(vA.getX()), new BigFraction(vB.getX()), new BigFraction(vC.getX()) + }; + final BigFraction[] c3 = new BigFraction[] { + new BigFraction(vA.getY()), new BigFraction(vB.getY()), new BigFraction(vC.getY()) + }; + final BigFraction[] c1 = new BigFraction[] { + c2[0].multiply(c2[0]).add(c3[0].multiply(c3[0])), + c2[1].multiply(c2[1]).add(c3[1].multiply(c3[1])), + c2[2].multiply(c2[2]).add(c3[2].multiply(c3[2])) + }; + final BigFraction twoM11 = minor(c2, c3).multiply(2); + final BigFraction m12 = minor(c1, c3); + final BigFraction m13 = minor(c1, c2); + final BigFraction centerX = m12.divide(twoM11); + final BigFraction centerY = m13.divide(twoM11).negate(); + final BigFraction dx = c2[0].subtract(centerX); + final BigFraction dy = c3[0].subtract(centerY); + final BigFraction r2 = dx.multiply(dx).add(dy.multiply(dy)); + return new EnclosingBall<Euclidean2D, Vector2D>(new Vector2D(centerX.doubleValue(), + centerY.doubleValue()), + FastMath.sqrt(r2.doubleValue()), + vA, vB, vC); + } + } + } + } + + /** Compute a dimension 3 minor, when 3<sup>d</sup> column is known to be filled with 1.0. + * @param c1 first column + * @param c2 second column + * @return value of the minor computed has an exact fraction + */ + private BigFraction minor(final BigFraction[] c1, final BigFraction[] c2) { + return c2[0].multiply(c1[2].subtract(c1[1])). + add(c2[1].multiply(c1[0].subtract(c1[2]))). + add(c2[2].multiply(c1[1].subtract(c1[0]))); + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Euclidean2D.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Euclidean2D.java new file mode 100644 index 0000000..af7630d --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Euclidean2D.java @@ -0,0 +1,74 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math3.geometry.euclidean.twod; + +import java.io.Serializable; + +import org.apache.commons.math3.geometry.Space; +import org.apache.commons.math3.geometry.euclidean.oned.Euclidean1D; + +/** + * This class implements a two-dimensional space. + * @since 3.0 + */ +public class Euclidean2D implements Serializable, Space { + + /** Serializable version identifier. */ + private static final long serialVersionUID = 4793432849757649566L; + + /** Private constructor for the singleton. + */ + private Euclidean2D() { + } + + /** Get the unique instance. + * @return the unique instance + */ + public static Euclidean2D getInstance() { + return LazyHolder.INSTANCE; + } + + /** {@inheritDoc} */ + public int getDimension() { + return 2; + } + + /** {@inheritDoc} */ + public Euclidean1D getSubSpace() { + return Euclidean1D.getInstance(); + } + + // CHECKSTYLE: stop HideUtilityClassConstructor + /** Holder for the instance. + * <p>We use here the Initialization On Demand Holder Idiom.</p> + */ + private static class LazyHolder { + /** Cached field instance. */ + private static final Euclidean2D INSTANCE = new Euclidean2D(); + } + // CHECKSTYLE: resume HideUtilityClassConstructor + + /** Handle deserialization of the singleton. + * @return the singleton instance + */ + private Object readResolve() { + // return the singleton instance + return LazyHolder.INSTANCE; + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java new file mode 100644 index 0000000..c300fa1 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java @@ -0,0 +1,587 @@ +/* + * 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.euclidean.twod; + +import java.awt.geom.AffineTransform; + +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.Vector; +import org.apache.commons.math3.geometry.euclidean.oned.Euclidean1D; +import org.apache.commons.math3.geometry.euclidean.oned.IntervalsSet; +import org.apache.commons.math3.geometry.euclidean.oned.OrientedPoint; +import org.apache.commons.math3.geometry.euclidean.oned.Vector1D; +import org.apache.commons.math3.geometry.partitioning.Embedding; +import org.apache.commons.math3.geometry.partitioning.Hyperplane; +import org.apache.commons.math3.geometry.partitioning.SubHyperplane; +import org.apache.commons.math3.geometry.partitioning.Transform; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathArrays; +import org.apache.commons.math3.util.MathUtils; + +/** This class represents an oriented line in the 2D plane. + + * <p>An oriented line can be defined either by prolongating a line + * segment between two points past these points, or by one point and + * an angular direction (in trigonometric orientation).</p> + + * <p>Since it is oriented the two half planes at its two sides are + * unambiguously identified as a left half plane and a right half + * plane. This can be used to identify the interior and the exterior + * in a simple way by local properties only when part of a line is + * used to define part of a polygon boundary.</p> + + * <p>A line can also be used to completely define a reference frame + * in the plane. It is sufficient to select one specific point in the + * line (the orthogonal projection of the original reference frame on + * the line) and to use the unit vector in the line direction and the + * orthogonal vector oriented from left half plane to right half + * plane. We define two coordinates by the process, the + * <em>abscissa</em> along the line, and the <em>offset</em> across + * the line. All points of the plane are uniquely identified by these + * two coordinates. The line is the set of points at zero offset, the + * left half plane is the set of points with negative offsets and the + * right half plane is the set of points with positive offsets.</p> + + * @since 3.0 + */ +public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euclidean1D> { + + /** Default value for tolerance. */ + private static final double DEFAULT_TOLERANCE = 1.0e-10; + + /** Angle with respect to the abscissa axis. */ + private double angle; + + /** Cosine of the line angle. */ + private double cos; + + /** Sine of the line angle. */ + private double sin; + + /** Offset of the frame origin. */ + private double originOffset; + + /** Tolerance below which points are considered identical. */ + private final double tolerance; + + /** Reverse line. */ + private Line reverse; + + /** Build a line from two points. + * <p>The line is oriented from p1 to p2</p> + * @param p1 first point + * @param p2 second point + * @param tolerance tolerance below which points are considered identical + * @since 3.3 + */ + public Line(final Vector2D p1, final Vector2D p2, final double tolerance) { + reset(p1, p2); + this.tolerance = tolerance; + } + + /** Build a line from a point and an angle. + * @param p point belonging to the line + * @param angle angle of the line with respect to abscissa axis + * @param tolerance tolerance below which points are considered identical + * @since 3.3 + */ + public Line(final Vector2D p, final double angle, final double tolerance) { + reset(p, angle); + this.tolerance = tolerance; + } + + /** Build a line from its internal characteristics. + * @param angle angle of the line with respect to abscissa axis + * @param cos cosine of the angle + * @param sin sine of the angle + * @param originOffset offset of the origin + * @param tolerance tolerance below which points are considered identical + * @since 3.3 + */ + private Line(final double angle, final double cos, final double sin, + final double originOffset, final double tolerance) { + this.angle = angle; + this.cos = cos; + this.sin = sin; + this.originOffset = originOffset; + this.tolerance = tolerance; + this.reverse = null; + } + + /** Build a line from two points. + * <p>The line is oriented from p1 to p2</p> + * @param p1 first point + * @param p2 second point + * @deprecated as of 3.3, replaced with {@link #Line(Vector2D, Vector2D, double)} + */ + @Deprecated + public Line(final Vector2D p1, final Vector2D p2) { + this(p1, p2, DEFAULT_TOLERANCE); + } + + /** Build a line from a point and an angle. + * @param p point belonging to the line + * @param angle angle of the line with respect to abscissa axis + * @deprecated as of 3.3, replaced with {@link #Line(Vector2D, double, double)} + */ + @Deprecated + public Line(final Vector2D p, final double angle) { + this(p, angle, DEFAULT_TOLERANCE); + } + + /** Copy constructor. + * <p>The created instance is completely independent from the + * original instance, it is a deep copy.</p> + * @param line line to copy + */ + public Line(final Line line) { + angle = MathUtils.normalizeAngle(line.angle, FastMath.PI); + cos = line.cos; + sin = line.sin; + originOffset = line.originOffset; + tolerance = line.tolerance; + reverse = null; + } + + /** {@inheritDoc} */ + public Line copySelf() { + return new Line(this); + } + + /** Reset the instance as if built from two points. + * <p>The line is oriented from p1 to p2</p> + * @param p1 first point + * @param p2 second point + */ + public void reset(final Vector2D p1, final Vector2D p2) { + unlinkReverse(); + final double dx = p2.getX() - p1.getX(); + final double dy = p2.getY() - p1.getY(); + final double d = FastMath.hypot(dx, dy); + if (d == 0.0) { + angle = 0.0; + cos = 1.0; + sin = 0.0; + originOffset = p1.getY(); + } else { + angle = FastMath.PI + FastMath.atan2(-dy, -dx); + cos = dx / d; + sin = dy / d; + originOffset = MathArrays.linearCombination(p2.getX(), p1.getY(), -p1.getX(), p2.getY()) / d; + } + } + + /** Reset the instance as if built from a line and an angle. + * @param p point belonging to the line + * @param alpha angle of the line with respect to abscissa axis + */ + public void reset(final Vector2D p, final double alpha) { + unlinkReverse(); + this.angle = MathUtils.normalizeAngle(alpha, FastMath.PI); + cos = FastMath.cos(this.angle); + sin = FastMath.sin(this.angle); + originOffset = MathArrays.linearCombination(cos, p.getY(), -sin, p.getX()); + } + + /** Revert the instance. + */ + public void revertSelf() { + unlinkReverse(); + if (angle < FastMath.PI) { + angle += FastMath.PI; + } else { + angle -= FastMath.PI; + } + cos = -cos; + sin = -sin; + originOffset = -originOffset; + } + + /** Unset the link between an instance and its reverse. + */ + private void unlinkReverse() { + if (reverse != null) { + reverse.reverse = null; + } + reverse = null; + } + + /** Get the reverse of the instance. + * <p>Get a line with reversed orientation with respect to the + * instance.</p> + * <p> + * As long as neither the instance nor its reverse are modified + * (i.e. as long as none of the {@link #reset(Vector2D, Vector2D)}, + * {@link #reset(Vector2D, double)}, {@link #revertSelf()}, + * {@link #setAngle(double)} or {@link #setOriginOffset(double)} + * methods are called), then the line and its reverse remain linked + * together so that {@code line.getReverse().getReverse() == line}. + * When one of the line is modified, the link is deleted as both + * instance becomes independent. + * </p> + * @return a new line, with orientation opposite to the instance orientation + */ + public Line getReverse() { + if (reverse == null) { + reverse = new Line((angle < FastMath.PI) ? (angle + FastMath.PI) : (angle - FastMath.PI), + -cos, -sin, -originOffset, tolerance); + reverse.reverse = this; + } + return reverse; + } + + /** Transform a space point into a sub-space point. + * @param vector n-dimension point of the space + * @return (n-1)-dimension point of the sub-space corresponding to + * the specified space point + */ + public Vector1D toSubSpace(Vector<Euclidean2D> vector) { + return toSubSpace((Point<Euclidean2D>) vector); + } + + /** Transform a sub-space point into a space point. + * @param vector (n-1)-dimension point of the sub-space + * @return n-dimension point of the space corresponding to the + * specified sub-space point + */ + public Vector2D toSpace(Vector<Euclidean1D> vector) { + return toSpace((Point<Euclidean1D>) vector); + } + + /** {@inheritDoc} */ + public Vector1D toSubSpace(final Point<Euclidean2D> point) { + Vector2D p2 = (Vector2D) point; + return new Vector1D(MathArrays.linearCombination(cos, p2.getX(), sin, p2.getY())); + } + + /** {@inheritDoc} */ + public Vector2D toSpace(final Point<Euclidean1D> point) { + final double abscissa = ((Vector1D) point).getX(); + return new Vector2D(MathArrays.linearCombination(abscissa, cos, -originOffset, sin), + MathArrays.linearCombination(abscissa, sin, originOffset, cos)); + } + + /** Get the intersection point of the instance and another line. + * @param other other line + * @return intersection point of the instance and the other line + * or null if there are no intersection points + */ + public Vector2D intersection(final Line other) { + final double d = MathArrays.linearCombination(sin, other.cos, -other.sin, cos); + if (FastMath.abs(d) < tolerance) { + return null; + } + return new Vector2D(MathArrays.linearCombination(cos, other.originOffset, -other.cos, originOffset) / d, + MathArrays.linearCombination(sin, other.originOffset, -other.sin, originOffset) / d); + } + + /** {@inheritDoc} + * @since 3.3 + */ + public Point<Euclidean2D> project(Point<Euclidean2D> point) { + return toSpace(toSubSpace(point)); + } + + /** {@inheritDoc} + * @since 3.3 + */ + public double getTolerance() { + return tolerance; + } + + /** {@inheritDoc} */ + public SubLine wholeHyperplane() { + return new SubLine(this, new IntervalsSet(tolerance)); + } + + /** Build a region covering the whole space. + * @return a region containing the instance (really a {@link + * PolygonsSet PolygonsSet} instance) + */ + public PolygonsSet wholeSpace() { + return new PolygonsSet(tolerance); + } + + /** Get the offset (oriented distance) of a parallel line. + * <p>This method should be called only for parallel lines otherwise + * the result is not meaningful.</p> + * <p>The offset is 0 if both lines are the same, it is + * positive if the line is on the right side of the instance and + * negative if it is on the left side, according to its natural + * orientation.</p> + * @param line line to check + * @return offset of the line + */ + public double getOffset(final Line line) { + return originOffset + + (MathArrays.linearCombination(cos, line.cos, sin, line.sin) > 0 ? -line.originOffset : line.originOffset); + } + + /** Get the offset (oriented distance) of a vector. + * @param vector vector to check + * @return offset of the vector + */ + public double getOffset(Vector<Euclidean2D> vector) { + return getOffset((Point<Euclidean2D>) vector); + } + + /** {@inheritDoc} */ + public double getOffset(final Point<Euclidean2D> point) { + Vector2D p2 = (Vector2D) point; + return MathArrays.linearCombination(sin, p2.getX(), -cos, p2.getY(), 1.0, originOffset); + } + + /** {@inheritDoc} */ + public boolean sameOrientationAs(final Hyperplane<Euclidean2D> other) { + final Line otherL = (Line) other; + return MathArrays.linearCombination(sin, otherL.sin, cos, otherL.cos) >= 0.0; + } + + /** Get one point from the plane. + * @param abscissa desired abscissa for the point + * @param offset desired offset for the point + * @return one point in the plane, with given abscissa and offset + * relative to the line + */ + public Vector2D getPointAt(final Vector1D abscissa, final double offset) { + final double x = abscissa.getX(); + final double dOffset = offset - originOffset; + return new Vector2D(MathArrays.linearCombination(x, cos, dOffset, sin), + MathArrays.linearCombination(x, sin, -dOffset, cos)); + } + + /** Check if the line contains a point. + * @param p point to check + * @return true if p belongs to the line + */ + public boolean contains(final Vector2D p) { + return FastMath.abs(getOffset(p)) < tolerance; + } + + /** Compute the distance between the instance and a point. + * <p>This is a shortcut for invoking FastMath.abs(getOffset(p)), + * and provides consistency with what is in the + * org.apache.commons.math3.geometry.euclidean.threed.Line class.</p> + * + * @param p to check + * @return distance between the instance and the point + * @since 3.1 + */ + public double distance(final Vector2D p) { + return FastMath.abs(getOffset(p)); + } + + /** Check the instance is parallel to another line. + * @param line other line to check + * @return true if the instance is parallel to the other line + * (they can have either the same or opposite orientations) + */ + public boolean isParallelTo(final Line line) { + return FastMath.abs(MathArrays.linearCombination(sin, line.cos, -cos, line.sin)) < tolerance; + } + + /** Translate the line to force it passing by a point. + * @param p point by which the line should pass + */ + public void translateToPoint(final Vector2D p) { + originOffset = MathArrays.linearCombination(cos, p.getY(), -sin, p.getX()); + } + + /** Get the angle of the line. + * @return the angle of the line with respect to the abscissa axis + */ + public double getAngle() { + return MathUtils.normalizeAngle(angle, FastMath.PI); + } + + /** Set the angle of the line. + * @param angle new angle of the line with respect to the abscissa axis + */ + public void setAngle(final double angle) { + unlinkReverse(); + this.angle = MathUtils.normalizeAngle(angle, FastMath.PI); + cos = FastMath.cos(this.angle); + sin = FastMath.sin(this.angle); + } + + /** Get the offset of the origin. + * @return the offset of the origin + */ + public double getOriginOffset() { + return originOffset; + } + + /** Set the offset of the origin. + * @param offset offset of the origin + */ + public void setOriginOffset(final double offset) { + unlinkReverse(); + originOffset = offset; + } + + /** Get a {@link org.apache.commons.math3.geometry.partitioning.Transform + * Transform} embedding an affine transform. + * @param transform affine transform to embed (must be inversible + * otherwise the {@link + * org.apache.commons.math3.geometry.partitioning.Transform#apply(Hyperplane) + * apply(Hyperplane)} method would work only for some lines, and + * fail for other ones) + * @return a new transform that can be applied to either {@link + * Vector2D Vector2D}, {@link Line Line} or {@link + * org.apache.commons.math3.geometry.partitioning.SubHyperplane + * SubHyperplane} instances + * @exception MathIllegalArgumentException if the transform is non invertible + * @deprecated as of 3.6, replaced with {@link #getTransform(double, double, double, double, double, double)} + */ + @Deprecated + public static Transform<Euclidean2D, Euclidean1D> getTransform(final AffineTransform transform) + throws MathIllegalArgumentException { + final double[] m = new double[6]; + transform.getMatrix(m); + return new LineTransform(m[0], m[1], m[2], m[3], m[4], m[5]); + } + + /** Get a {@link org.apache.commons.math3.geometry.partitioning.Transform + * Transform} embedding an affine transform. + * @param cXX transform factor between input abscissa and output abscissa + * @param cYX transform factor between input abscissa and output ordinate + * @param cXY transform factor between input ordinate and output abscissa + * @param cYY transform factor between input ordinate and output ordinate + * @param cX1 transform addendum for output abscissa + * @param cY1 transform addendum for output ordinate + * @return a new transform that can be applied to either {@link + * Vector2D Vector2D}, {@link Line Line} or {@link + * org.apache.commons.math3.geometry.partitioning.SubHyperplane + * SubHyperplane} instances + * @exception MathIllegalArgumentException if the transform is non invertible + * @since 3.6 + */ + public static Transform<Euclidean2D, Euclidean1D> getTransform(final double cXX, + final double cYX, + final double cXY, + final double cYY, + final double cX1, + final double cY1) + throws MathIllegalArgumentException { + return new LineTransform(cXX, cYX, cXY, cYY, cX1, cY1); + } + + /** Class embedding an affine transform. + * <p>This class is used in order to apply an affine transform to a + * line. Using a specific object allow to perform some computations + * on the transform only once even if the same transform is to be + * applied to a large number of lines (for example to a large + * polygon)./<p> + */ + private static class LineTransform implements Transform<Euclidean2D, Euclidean1D> { + + /** Transform factor between input abscissa and output abscissa. */ + private double cXX; + + /** Transform factor between input abscissa and output ordinate. */ + private double cYX; + + /** Transform factor between input ordinate and output abscissa. */ + private double cXY; + + /** Transform factor between input ordinate and output ordinate. */ + private double cYY; + + /** Transform addendum for output abscissa. */ + private double cX1; + + /** Transform addendum for output ordinate. */ + private double cY1; + + /** cXY * cY1 - cYY * cX1. */ + private double c1Y; + + /** cXX * cY1 - cYX * cX1. */ + private double c1X; + + /** cXX * cYY - cYX * cXY. */ + private double c11; + + /** Build an affine line transform from a n {@code AffineTransform}. + * @param cXX transform factor between input abscissa and output abscissa + * @param cYX transform factor between input abscissa and output ordinate + * @param cXY transform factor between input ordinate and output abscissa + * @param cYY transform factor between input ordinate and output ordinate + * @param cX1 transform addendum for output abscissa + * @param cY1 transform addendum for output ordinate + * @exception MathIllegalArgumentException if the transform is non invertible + * @since 3.6 + */ + LineTransform(final double cXX, final double cYX, final double cXY, + final double cYY, final double cX1, final double cY1) + throws MathIllegalArgumentException { + + this.cXX = cXX; + this.cYX = cYX; + this.cXY = cXY; + this.cYY = cYY; + this.cX1 = cX1; + this.cY1 = cY1; + + c1Y = MathArrays.linearCombination(cXY, cY1, -cYY, cX1); + c1X = MathArrays.linearCombination(cXX, cY1, -cYX, cX1); + c11 = MathArrays.linearCombination(cXX, cYY, -cYX, cXY); + + if (FastMath.abs(c11) < 1.0e-20) { + throw new MathIllegalArgumentException(LocalizedFormats.NON_INVERTIBLE_TRANSFORM); + } + + } + + /** {@inheritDoc} */ + public Vector2D apply(final Point<Euclidean2D> point) { + final Vector2D p2D = (Vector2D) point; + final double x = p2D.getX(); + final double y = p2D.getY(); + return new Vector2D(MathArrays.linearCombination(cXX, x, cXY, y, cX1, 1), + MathArrays.linearCombination(cYX, x, cYY, y, cY1, 1)); + } + + /** {@inheritDoc} */ + public Line apply(final Hyperplane<Euclidean2D> hyperplane) { + final Line line = (Line) hyperplane; + final double rOffset = MathArrays.linearCombination(c1X, line.cos, c1Y, line.sin, c11, line.originOffset); + final double rCos = MathArrays.linearCombination(cXX, line.cos, cXY, line.sin); + final double rSin = MathArrays.linearCombination(cYX, line.cos, cYY, line.sin); + final double inv = 1.0 / FastMath.sqrt(rSin * rSin + rCos * rCos); + return new Line(FastMath.PI + FastMath.atan2(-rSin, -rCos), + inv * rCos, inv * rSin, + inv * rOffset, line.tolerance); + } + + /** {@inheritDoc} */ + public SubHyperplane<Euclidean1D> apply(final SubHyperplane<Euclidean1D> sub, + final Hyperplane<Euclidean2D> original, + final Hyperplane<Euclidean2D> transformed) { + final OrientedPoint op = (OrientedPoint) sub.getHyperplane(); + final Line originalLine = (Line) original; + final Line transformedLine = (Line) transformed; + final Vector1D newLoc = + transformedLine.toSubSpace(apply(originalLine.toSpace(op.getLocation()))); + return new OrientedPoint(newLoc, op.isDirect(), originalLine.tolerance).wholeHyperplane(); + } + + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/NestedLoops.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/NestedLoops.java new file mode 100644 index 0000000..83928fa --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/NestedLoops.java @@ -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.euclidean.twod; + +import java.util.ArrayList; +import java.util.Iterator; +import java.util.List; + +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.euclidean.oned.IntervalsSet; +import org.apache.commons.math3.geometry.partitioning.Region; +import org.apache.commons.math3.geometry.partitioning.RegionFactory; +import org.apache.commons.math3.geometry.partitioning.SubHyperplane; + +/** This class represent a tree of nested 2D boundary loops. + + * <p>This class is used for piecewise polygons construction. + * Polygons are built using the outline edges as + * representative of boundaries, the orientation of these lines are + * meaningful. However, we want to allow the user to specify its + * outline loops without having to take care of this orientation. This + * class is devoted to correct mis-oriented loops.<p> + + * <p>Orientation is computed assuming the piecewise polygon is finite, + * i.e. the outermost loops have their exterior side facing points at + * infinity, and hence are oriented counter-clockwise. The orientation of + * internal loops is computed as the reverse of the orientation of + * their immediate surrounding loop.</p> + + * @since 3.0 + */ +class NestedLoops { + + /** Boundary loop. */ + private Vector2D[] loop; + + /** Surrounded loops. */ + private List<NestedLoops> surrounded; + + /** Polygon enclosing a finite region. */ + private Region<Euclidean2D> polygon; + + /** Indicator for original loop orientation. */ + private boolean originalIsClockwise; + + /** Tolerance below which points are considered identical. */ + private final double tolerance; + + /** Simple Constructor. + * <p>Build an empty tree of nested loops. This instance will become + * the root node of a complete tree, it is not associated with any + * loop by itself, the outermost loops are in the root tree child + * nodes.</p> + * @param tolerance tolerance below which points are considered identical + * @since 3.3 + */ + NestedLoops(final double tolerance) { + this.surrounded = new ArrayList<NestedLoops>(); + this.tolerance = tolerance; + } + + /** Constructor. + * <p>Build a tree node with neither parent nor children</p> + * @param loop boundary loop (will be reversed in place if needed) + * @param tolerance tolerance below which points are considered identical + * @exception MathIllegalArgumentException if an outline has an open boundary loop + * @since 3.3 + */ + private NestedLoops(final Vector2D[] loop, final double tolerance) + throws MathIllegalArgumentException { + + if (loop[0] == null) { + throw new MathIllegalArgumentException(LocalizedFormats.OUTLINE_BOUNDARY_LOOP_OPEN); + } + + this.loop = loop; + this.surrounded = new ArrayList<NestedLoops>(); + this.tolerance = tolerance; + + // build the polygon defined by the loop + final ArrayList<SubHyperplane<Euclidean2D>> edges = new ArrayList<SubHyperplane<Euclidean2D>>(); + Vector2D current = loop[loop.length - 1]; + for (int i = 0; i < loop.length; ++i) { + final Vector2D previous = current; + current = loop[i]; + final Line line = new Line(previous, current, tolerance); + final IntervalsSet region = + new IntervalsSet(line.toSubSpace((Point<Euclidean2D>) previous).getX(), + line.toSubSpace((Point<Euclidean2D>) current).getX(), + tolerance); + edges.add(new SubLine(line, region)); + } + polygon = new PolygonsSet(edges, tolerance); + + // ensure the polygon encloses a finite region of the plane + if (Double.isInfinite(polygon.getSize())) { + polygon = new RegionFactory<Euclidean2D>().getComplement(polygon); + originalIsClockwise = false; + } else { + originalIsClockwise = true; + } + + } + + /** Add a loop in a tree. + * @param bLoop boundary loop (will be reversed in place if needed) + * @exception MathIllegalArgumentException if an outline has crossing + * boundary loops or open boundary loops + */ + public void add(final Vector2D[] bLoop) throws MathIllegalArgumentException { + add(new NestedLoops(bLoop, tolerance)); + } + + /** Add a loop in a tree. + * @param node boundary loop (will be reversed in place if needed) + * @exception MathIllegalArgumentException if an outline has boundary + * loops that cross each other + */ + private void add(final NestedLoops node) throws MathIllegalArgumentException { + + // check if we can go deeper in the tree + for (final NestedLoops child : surrounded) { + if (child.polygon.contains(node.polygon)) { + child.add(node); + return; + } + } + + // check if we can absorb some of the instance children + for (final Iterator<NestedLoops> iterator = surrounded.iterator(); iterator.hasNext();) { + final NestedLoops child = iterator.next(); + if (node.polygon.contains(child.polygon)) { + node.surrounded.add(child); + iterator.remove(); + } + } + + // we should be separate from the remaining children + RegionFactory<Euclidean2D> factory = new RegionFactory<Euclidean2D>(); + for (final NestedLoops child : surrounded) { + if (!factory.intersection(node.polygon, child.polygon).isEmpty()) { + throw new MathIllegalArgumentException(LocalizedFormats.CROSSING_BOUNDARY_LOOPS); + } + } + + surrounded.add(node); + + } + + /** Correct the orientation of the loops contained in the tree. + * <p>This is this method that really inverts the loops that where + * provided through the {@link #add(Vector2D[]) add} method if + * they are mis-oriented</p> + */ + public void correctOrientation() { + for (NestedLoops child : surrounded) { + child.setClockWise(true); + } + } + + /** Set the loop orientation. + * @param clockwise if true, the loop should be set to clockwise + * orientation + */ + private void setClockWise(final boolean clockwise) { + + if (originalIsClockwise ^ clockwise) { + // we need to inverse the original loop + int min = -1; + int max = loop.length; + while (++min < --max) { + final Vector2D tmp = loop[min]; + loop[min] = loop[max]; + loop[max] = tmp; + } + } + + // go deeper in the tree + for (final NestedLoops child : surrounded) { + child.setClockWise(!clockwise); + } + + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSet.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSet.java new file mode 100644 index 0000000..61fae9f --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSet.java @@ -0,0 +1,1160 @@ +/* + * 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.euclidean.twod; + +import java.util.ArrayList; +import java.util.Collection; +import java.util.List; + +import org.apache.commons.math3.geometry.Point; +import org.apache.commons.math3.geometry.euclidean.oned.Euclidean1D; +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.AbstractRegion; +import org.apache.commons.math3.geometry.partitioning.AbstractSubHyperplane; +import org.apache.commons.math3.geometry.partitioning.BSPTree; +import org.apache.commons.math3.geometry.partitioning.BSPTreeVisitor; +import org.apache.commons.math3.geometry.partitioning.BoundaryAttribute; +import org.apache.commons.math3.geometry.partitioning.Hyperplane; +import org.apache.commons.math3.geometry.partitioning.Side; +import org.apache.commons.math3.geometry.partitioning.SubHyperplane; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.Precision; + +/** This class represents a 2D region: a set of polygons. + * @since 3.0 + */ +public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> { + + /** Default value for tolerance. */ + private static final double DEFAULT_TOLERANCE = 1.0e-10; + + /** Vertices organized as boundary loops. */ + private Vector2D[][] vertices; + + /** Build a polygons set representing the whole plane. + * @param tolerance tolerance below which points are considered identical + * @since 3.3 + */ + public PolygonsSet(final double tolerance) { + super(tolerance); + } + + /** Build a polygons set from a 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}</p> + * <p> + * This constructor is aimed at expert use, as building the tree may + * be a difficult task. It is not intended for general use and for + * performances reasons does not check thoroughly its input, as this would + * require walking the full tree each time. Failing to provide a tree with + * the proper attributes, <em>will</em> therefore generate problems like + * {@link NullPointerException} or {@link ClassCastException} only later on. + * This limitation is known and explains why this constructor is for expert + * use only. The caller does have the responsibility to provided correct arguments. + * </p> + * @param tree inside/outside BSP tree representing the region + * @param tolerance tolerance below which points are considered identical + * @since 3.3 + */ + public PolygonsSet(final BSPTree<Euclidean2D> tree, final double tolerance) { + super(tree, tolerance); + } + + /** Build a polygons set 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 disjoint polygons 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 + * org.apache.commons.math3.geometry.partitioning.Region#checkPoint(org.apache.commons.math3.geometry.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 + * @since 3.3 + */ + public PolygonsSet(final Collection<SubHyperplane<Euclidean2D>> boundary, final double tolerance) { + super(boundary, tolerance); + } + + /** Build a parallellepipedic box. + * @param xMin low bound along the x direction + * @param xMax high bound along the x direction + * @param yMin low bound along the y direction + * @param yMax high bound along the y direction + * @param tolerance tolerance below which points are considered identical + * @since 3.3 + */ + public PolygonsSet(final double xMin, final double xMax, + final double yMin, final double yMax, + final double tolerance) { + super(boxBoundary(xMin, xMax, yMin, yMax, tolerance), tolerance); + } + + /** Build a polygon from a simple list of vertices. + * <p>The boundary is provided as a list of points considering to + * represent the vertices of a simple loop. The interior part of the + * region is on the left side of this path and the exterior is on its + * right side.</p> + * <p>This constructor does not handle polygons with a boundary + * forming several disconnected paths (such as polygons with holes).</p> + * <p>For cases where this simple constructor applies, it is expected to + * be numerically more robust than the {@link #PolygonsSet(Collection) general + * constructor} using {@link SubHyperplane subhyperplanes}.</p> + * <p>If the list is empty, the region will represent the whole + * space.</p> + * <p> + * Polygons with thin pikes or dents are inherently difficult to handle because + * they involve lines with almost opposite directions at some vertices. Polygons + * whose vertices come from some physical measurement with noise are also + * difficult because an edge that should be straight may be broken in lots of + * different pieces with almost equal directions. In both cases, computing the + * lines intersections is not numerically robust due to the almost 0 or almost + * π angle. Such cases need to carefully adjust the {@code hyperplaneThickness} + * parameter. A too small value would often lead to completely wrong polygons + * with large area wrongly identified as inside or outside. Large values are + * often much safer. As a rule of thumb, a value slightly below the size of the + * most accurate detail needed is a good value for the {@code hyperplaneThickness} + * parameter. + * </p> + * @param hyperplaneThickness tolerance below which points are considered to + * belong to the hyperplane (which is therefore more a slab) + * @param vertices vertices of the simple loop boundary + */ + public PolygonsSet(final double hyperplaneThickness, final Vector2D ... vertices) { + super(verticesToTree(hyperplaneThickness, vertices), hyperplaneThickness); + } + + /** Build a polygons set representing the whole real line. + * @deprecated as of 3.3, replaced with {@link #PolygonsSet(double)} + */ + @Deprecated + public PolygonsSet() { + this(DEFAULT_TOLERANCE); + } + + /** Build a polygons set from a 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}</p> + * @param tree inside/outside BSP tree representing the region + * @deprecated as of 3.3, replaced with {@link #PolygonsSet(BSPTree, double)} + */ + @Deprecated + public PolygonsSet(final BSPTree<Euclidean2D> tree) { + this(tree, DEFAULT_TOLERANCE); + } + + /** Build a polygons set 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 disjoint polygons 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 + * org.apache.commons.math3.geometry.partitioning.Region#checkPoint(org.apache.commons.math3.geometry.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 + * @deprecated as of 3.3, replaced with {@link #PolygonsSet(Collection, double)} + */ + @Deprecated + public PolygonsSet(final Collection<SubHyperplane<Euclidean2D>> boundary) { + this(boundary, DEFAULT_TOLERANCE); + } + + /** Build a parallellepipedic box. + * @param xMin low bound along the x direction + * @param xMax high bound along the x direction + * @param yMin low bound along the y direction + * @param yMax high bound along the y direction + * @deprecated as of 3.3, replaced with {@link #PolygonsSet(double, double, double, double, double)} + */ + @Deprecated + public PolygonsSet(final double xMin, final double xMax, + final double yMin, final double yMax) { + this(xMin, xMax, yMin, yMax, DEFAULT_TOLERANCE); + } + + /** Create a list of hyperplanes representing the boundary of a box. + * @param xMin low bound along the x direction + * @param xMax high bound along the x direction + * @param yMin low bound along the y direction + * @param yMax high bound along the y direction + * @param tolerance tolerance below which points are considered identical + * @return boundary of the box + */ + private static Line[] boxBoundary(final double xMin, final double xMax, + final double yMin, final double yMax, + final double tolerance) { + if ((xMin >= xMax - tolerance) || (yMin >= yMax - tolerance)) { + // too thin box, build an empty polygons set + return null; + } + final Vector2D minMin = new Vector2D(xMin, yMin); + final Vector2D minMax = new Vector2D(xMin, yMax); + final Vector2D maxMin = new Vector2D(xMax, yMin); + final Vector2D maxMax = new Vector2D(xMax, yMax); + return new Line[] { + new Line(minMin, maxMin, tolerance), + new Line(maxMin, maxMax, tolerance), + new Line(maxMax, minMax, tolerance), + new Line(minMax, minMin, tolerance) + }; + } + + /** Build the BSP tree of a polygons set from a simple list of vertices. + * <p>The boundary is provided as a list of points considering to + * represent the vertices of a simple loop. The interior part of the + * region is on the left side of this path and the exterior is on its + * right side.</p> + * <p>This constructor does not handle polygons with a boundary + * forming several disconnected paths (such as polygons with holes).</p> + * <p>For cases where this simple constructor applies, it is expected to + * be numerically more robust than the {@link #PolygonsSet(Collection) general + * constructor} using {@link SubHyperplane subhyperplanes}.</p> + * @param hyperplaneThickness tolerance below which points are consider to + * belong to the hyperplane (which is therefore more a slab) + * @param vertices vertices of the simple loop boundary + * @return the BSP tree of the input vertices + */ + private static BSPTree<Euclidean2D> verticesToTree(final double hyperplaneThickness, + final Vector2D ... vertices) { + + final int n = vertices.length; + if (n == 0) { + // the tree represents the whole space + return new BSPTree<Euclidean2D>(Boolean.TRUE); + } + + // build the vertices + final Vertex[] vArray = new Vertex[n]; + for (int i = 0; i < n; ++i) { + vArray[i] = new Vertex(vertices[i]); + } + + // build the edges + List<Edge> edges = new ArrayList<Edge>(n); + for (int i = 0; i < n; ++i) { + + // get the endpoints of the edge + final Vertex start = vArray[i]; + final Vertex end = vArray[(i + 1) % n]; + + // get the line supporting the edge, taking care not to recreate it + // if it was already created earlier due to another edge being aligned + // with the current one + Line line = start.sharedLineWith(end); + if (line == null) { + line = new Line(start.getLocation(), end.getLocation(), hyperplaneThickness); + } + + // create the edge and store it + edges.add(new Edge(start, end, line)); + + // check if another vertex also happens to be on this line + for (final Vertex vertex : vArray) { + if (vertex != start && vertex != end && + FastMath.abs(line.getOffset((Point<Euclidean2D>) vertex.getLocation())) <= hyperplaneThickness) { + vertex.bindWith(line); + } + } + + } + + // build the tree top-down + final BSPTree<Euclidean2D> tree = new BSPTree<Euclidean2D>(); + insertEdges(hyperplaneThickness, tree, edges); + + return tree; + + } + + /** Recursively build a tree by inserting cut sub-hyperplanes. + * @param hyperplaneThickness tolerance below which points are consider to + * belong to the hyperplane (which is therefore more a slab) + * @param node current tree node (it is a leaf node at the beginning + * of the call) + * @param edges list of edges to insert in the cell defined by this node + * (excluding edges not belonging to the cell defined by this node) + */ + private static void insertEdges(final double hyperplaneThickness, + final BSPTree<Euclidean2D> node, + final List<Edge> edges) { + + // find an edge with an hyperplane that can be inserted in the node + int index = 0; + Edge inserted =null; + while (inserted == null && index < edges.size()) { + inserted = edges.get(index++); + if (inserted.getNode() == null) { + if (node.insertCut(inserted.getLine())) { + inserted.setNode(node); + } else { + inserted = null; + } + } else { + inserted = null; + } + } + + if (inserted == null) { + // no suitable edge was found, the node remains a leaf node + // we need to set its inside/outside boolean indicator + final BSPTree<Euclidean2D> parent = node.getParent(); + if (parent == null || node == parent.getMinus()) { + node.setAttribute(Boolean.TRUE); + } else { + node.setAttribute(Boolean.FALSE); + } + return; + } + + // we have split the node by inserting an edge as a cut sub-hyperplane + // distribute the remaining edges in the two sub-trees + final List<Edge> plusList = new ArrayList<Edge>(); + final List<Edge> minusList = new ArrayList<Edge>(); + for (final Edge edge : edges) { + if (edge != inserted) { + final double startOffset = inserted.getLine().getOffset((Point<Euclidean2D>) edge.getStart().getLocation()); + final double endOffset = inserted.getLine().getOffset((Point<Euclidean2D>) edge.getEnd().getLocation()); + Side startSide = (FastMath.abs(startOffset) <= hyperplaneThickness) ? + Side.HYPER : ((startOffset < 0) ? Side.MINUS : Side.PLUS); + Side endSide = (FastMath.abs(endOffset) <= hyperplaneThickness) ? + Side.HYPER : ((endOffset < 0) ? Side.MINUS : Side.PLUS); + switch (startSide) { + case PLUS: + if (endSide == Side.MINUS) { + // we need to insert a split point on the hyperplane + final Vertex splitPoint = edge.split(inserted.getLine()); + minusList.add(splitPoint.getOutgoing()); + plusList.add(splitPoint.getIncoming()); + } else { + plusList.add(edge); + } + break; + case MINUS: + if (endSide == Side.PLUS) { + // we need to insert a split point on the hyperplane + final Vertex splitPoint = edge.split(inserted.getLine()); + minusList.add(splitPoint.getIncoming()); + plusList.add(splitPoint.getOutgoing()); + } else { + minusList.add(edge); + } + break; + default: + if (endSide == Side.PLUS) { + plusList.add(edge); + } else if (endSide == Side.MINUS) { + minusList.add(edge); + } + break; + } + } + } + + // recurse through lower levels + if (!plusList.isEmpty()) { + insertEdges(hyperplaneThickness, node.getPlus(), plusList); + } else { + node.getPlus().setAttribute(Boolean.FALSE); + } + if (!minusList.isEmpty()) { + insertEdges(hyperplaneThickness, node.getMinus(), minusList); + } else { + node.getMinus().setAttribute(Boolean.TRUE); + } + + } + + /** Internal class for holding vertices while they are processed to build a BSP tree. */ + private static class Vertex { + + /** Vertex location. */ + private final Vector2D location; + + /** Incoming edge. */ + private Edge incoming; + + /** Outgoing edge. */ + private Edge outgoing; + + /** Lines bound with this vertex. */ + private final List<Line> lines; + + /** Build a non-processed vertex not owned by any node yet. + * @param location vertex location + */ + Vertex(final Vector2D location) { + this.location = location; + this.incoming = null; + this.outgoing = null; + this.lines = new ArrayList<Line>(); + } + + /** Get Vertex location. + * @return vertex location + */ + public Vector2D getLocation() { + return location; + } + + /** Bind a line considered to contain this vertex. + * @param line line to bind with this vertex + */ + public void bindWith(final Line line) { + lines.add(line); + } + + /** Get the common line bound with both the instance and another vertex, if any. + * <p> + * When two vertices are both bound to the same line, this means they are + * already handled by node associated with this line, so there is no need + * to create a cut hyperplane for them. + * </p> + * @param vertex other vertex to check instance against + * @return line bound with both the instance and another vertex, or null if the + * two vertices do not share a line yet + */ + public Line sharedLineWith(final Vertex vertex) { + for (final Line line1 : lines) { + for (final Line line2 : vertex.lines) { + if (line1 == line2) { + return line1; + } + } + } + return null; + } + + /** Set incoming edge. + * <p> + * The line supporting the incoming edge is automatically bound + * with the instance. + * </p> + * @param incoming incoming edge + */ + public void setIncoming(final Edge incoming) { + this.incoming = incoming; + bindWith(incoming.getLine()); + } + + /** Get incoming edge. + * @return incoming edge + */ + public Edge getIncoming() { + return incoming; + } + + /** Set outgoing edge. + * <p> + * The line supporting the outgoing edge is automatically bound + * with the instance. + * </p> + * @param outgoing outgoing edge + */ + public void setOutgoing(final Edge outgoing) { + this.outgoing = outgoing; + bindWith(outgoing.getLine()); + } + + /** Get outgoing edge. + * @return outgoing edge + */ + public Edge getOutgoing() { + return outgoing; + } + + } + + /** Internal class for holding edges while they are processed to build a BSP tree. */ + private static class Edge { + + /** Start vertex. */ + private final Vertex start; + + /** End vertex. */ + private final Vertex end; + + /** Line supporting the edge. */ + private final Line line; + + /** Node whose cut hyperplane contains this edge. */ + private BSPTree<Euclidean2D> node; + + /** Build an edge not contained in any node yet. + * @param start start vertex + * @param end end vertex + * @param line line supporting the edge + */ + Edge(final Vertex start, final Vertex end, final Line line) { + + this.start = start; + this.end = end; + this.line = line; + this.node = null; + + // connect the vertices back to the edge + start.setOutgoing(this); + end.setIncoming(this); + + } + + /** Get start vertex. + * @return start vertex + */ + public Vertex getStart() { + return start; + } + + /** Get end vertex. + * @return end vertex + */ + public Vertex getEnd() { + return end; + } + + /** Get the line supporting this edge. + * @return line supporting this edge + */ + public Line getLine() { + return line; + } + + /** Set the node whose cut hyperplane contains this edge. + * @param node node whose cut hyperplane contains this edge + */ + public void setNode(final BSPTree<Euclidean2D> node) { + this.node = node; + } + + /** Get the node whose cut hyperplane contains this edge. + * @return node whose cut hyperplane contains this edge + * (null if edge has not yet been inserted into the BSP tree) + */ + public BSPTree<Euclidean2D> getNode() { + return node; + } + + /** Split the edge. + * <p> + * Once split, this edge is not referenced anymore by the vertices, + * it is replaced by the two half-edges and an intermediate splitting + * vertex is introduced to connect these two halves. + * </p> + * @param splitLine line splitting the edge in two halves + * @return split vertex (its incoming and outgoing edges are the two halves) + */ + public Vertex split(final Line splitLine) { + final Vertex splitVertex = new Vertex(line.intersection(splitLine)); + splitVertex.bindWith(splitLine); + final Edge startHalf = new Edge(start, splitVertex, line); + final Edge endHalf = new Edge(splitVertex, end, line); + startHalf.node = node; + endHalf.node = node; + return splitVertex; + } + + } + + /** {@inheritDoc} */ + @Override + public PolygonsSet buildNew(final BSPTree<Euclidean2D> tree) { + return new PolygonsSet(tree, getTolerance()); + } + + /** {@inheritDoc} */ + @Override + protected void computeGeometricalProperties() { + + final Vector2D[][] v = getVertices(); + + if (v.length == 0) { + final BSPTree<Euclidean2D> tree = getTree(false); + if (tree.getCut() == null && (Boolean) tree.getAttribute()) { + // the instance covers the whole space + setSize(Double.POSITIVE_INFINITY); + setBarycenter((Point<Euclidean2D>) Vector2D.NaN); + } else { + setSize(0); + setBarycenter((Point<Euclidean2D>) new Vector2D(0, 0)); + } + } else if (v[0][0] == null) { + // there is at least one open-loop: the polygon is infinite + setSize(Double.POSITIVE_INFINITY); + setBarycenter((Point<Euclidean2D>) Vector2D.NaN); + } else { + // all loops are closed, we compute some integrals around the shape + + double sum = 0; + double sumX = 0; + double sumY = 0; + + for (Vector2D[] loop : v) { + double x1 = loop[loop.length - 1].getX(); + double y1 = loop[loop.length - 1].getY(); + for (final Vector2D point : loop) { + final double x0 = x1; + final double y0 = y1; + x1 = point.getX(); + y1 = point.getY(); + final double factor = x0 * y1 - y0 * x1; + sum += factor; + sumX += factor * (x0 + x1); + sumY += factor * (y0 + y1); + } + } + + if (sum < 0) { + // the polygon as a finite outside surrounded by an infinite inside + setSize(Double.POSITIVE_INFINITY); + setBarycenter((Point<Euclidean2D>) Vector2D.NaN); + } else { + setSize(sum / 2); + setBarycenter((Point<Euclidean2D>) new Vector2D(sumX / (3 * sum), sumY / (3 * sum))); + } + + } + + } + + /** Get the vertices of the polygon. + * <p>The polygon boundary can be represented as an array of loops, + * each loop being itself an array of vertices.</p> + * <p>In order to identify open loops which start and end by + * infinite edges, the open loops arrays start with a null point. In + * this case, the first non null point and the last point of the + * array do not represent real vertices, they are dummy points + * intended only to get the direction of the first and last edge. An + * open loop consisting of a single infinite line will therefore be + * represented by a three elements array with one null point + * followed by two dummy points. The open loops are always the first + * ones in the loops array.</p> + * <p>If the polygon has no boundary at all, a zero length loop + * array will be returned.</p> + * <p>All line segments in the various loops have the inside of the + * region on their left side and the outside on their right side + * when moving in the underlying line direction. This means that + * closed loops surrounding finite areas obey the direct + * trigonometric orientation.</p> + * @return vertices of the polygon, organized as oriented boundary + * loops with the open loops first (the returned value is guaranteed + * to be non-null) + */ + public Vector2D[][] getVertices() { + if (vertices == null) { + if (getTree(false).getCut() == null) { + vertices = new Vector2D[0][]; + } else { + + // build the unconnected segments + final SegmentsBuilder visitor = new SegmentsBuilder(getTolerance()); + getTree(true).visit(visitor); + final List<ConnectableSegment> segments = visitor.getSegments(); + + // connect all segments, using topological criteria first + // and using Euclidean distance only as a last resort + int pending = segments.size(); + pending -= naturalFollowerConnections(segments); + if (pending > 0) { + pending -= splitEdgeConnections(segments); + } + if (pending > 0) { + pending -= closeVerticesConnections(segments); + } + + // create the segment loops + final ArrayList<List<Segment>> loops = new ArrayList<List<Segment>>(); + for (ConnectableSegment s = getUnprocessed(segments); s != null; s = getUnprocessed(segments)) { + final List<Segment> loop = followLoop(s); + if (loop != null) { + if (loop.get(0).getStart() == null) { + // this is an open loop, we put it on the front + loops.add(0, loop); + } else { + // this is a closed loop, we put it on the back + loops.add(loop); + } + } + } + + // transform the loops in an array of arrays of points + vertices = new Vector2D[loops.size()][]; + int i = 0; + + for (final List<Segment> loop : loops) { + if (loop.size() < 2 || + (loop.size() == 2 && loop.get(0).getStart() == null && loop.get(1).getEnd() == null)) { + // single infinite line + final Line line = loop.get(0).getLine(); + vertices[i++] = new Vector2D[] { + null, + line.toSpace((Point<Euclidean1D>) new Vector1D(-Float.MAX_VALUE)), + line.toSpace((Point<Euclidean1D>) new Vector1D(+Float.MAX_VALUE)) + }; + } else if (loop.get(0).getStart() == null) { + // open loop with at least one real point + final Vector2D[] array = new Vector2D[loop.size() + 2]; + int j = 0; + for (Segment segment : loop) { + + if (j == 0) { + // null point and first dummy point + double x = segment.getLine().toSubSpace((Point<Euclidean2D>) segment.getEnd()).getX(); + x -= FastMath.max(1.0, FastMath.abs(x / 2)); + array[j++] = null; + array[j++] = segment.getLine().toSpace((Point<Euclidean1D>) new Vector1D(x)); + } + + if (j < (array.length - 1)) { + // current point + array[j++] = segment.getEnd(); + } + + if (j == (array.length - 1)) { + // last dummy point + double x = segment.getLine().toSubSpace((Point<Euclidean2D>) segment.getStart()).getX(); + x += FastMath.max(1.0, FastMath.abs(x / 2)); + array[j++] = segment.getLine().toSpace((Point<Euclidean1D>) new Vector1D(x)); + } + + } + vertices[i++] = array; + } else { + final Vector2D[] array = new Vector2D[loop.size()]; + int j = 0; + for (Segment segment : loop) { + array[j++] = segment.getStart(); + } + vertices[i++] = array; + } + } + + } + } + + return vertices.clone(); + + } + + /** Connect the segments using only natural follower information. + * @param segments segments complete segments list + * @return number of connections performed + */ + private int naturalFollowerConnections(final List<ConnectableSegment> segments) { + int connected = 0; + for (final ConnectableSegment segment : segments) { + if (segment.getNext() == null) { + final BSPTree<Euclidean2D> node = segment.getNode(); + final BSPTree<Euclidean2D> end = segment.getEndNode(); + for (final ConnectableSegment candidateNext : segments) { + if (candidateNext.getPrevious() == null && + candidateNext.getNode() == end && + candidateNext.getStartNode() == node) { + // connect the two segments + segment.setNext(candidateNext); + candidateNext.setPrevious(segment); + ++connected; + break; + } + } + } + } + return connected; + } + + /** Connect the segments resulting from a line splitting a straight edge. + * @param segments segments complete segments list + * @return number of connections performed + */ + private int splitEdgeConnections(final List<ConnectableSegment> segments) { + int connected = 0; + for (final ConnectableSegment segment : segments) { + if (segment.getNext() == null) { + final Hyperplane<Euclidean2D> hyperplane = segment.getNode().getCut().getHyperplane(); + final BSPTree<Euclidean2D> end = segment.getEndNode(); + for (final ConnectableSegment candidateNext : segments) { + if (candidateNext.getPrevious() == null && + candidateNext.getNode().getCut().getHyperplane() == hyperplane && + candidateNext.getStartNode() == end) { + // connect the two segments + segment.setNext(candidateNext); + candidateNext.setPrevious(segment); + ++connected; + break; + } + } + } + } + return connected; + } + + /** Connect the segments using Euclidean distance. + * <p> + * This connection heuristic should be used last, as it relies + * only on a fuzzy distance criterion. + * </p> + * @param segments segments complete segments list + * @return number of connections performed + */ + private int closeVerticesConnections(final List<ConnectableSegment> segments) { + int connected = 0; + for (final ConnectableSegment segment : segments) { + if (segment.getNext() == null && segment.getEnd() != null) { + final Vector2D end = segment.getEnd(); + ConnectableSegment selectedNext = null; + double min = Double.POSITIVE_INFINITY; + for (final ConnectableSegment candidateNext : segments) { + if (candidateNext.getPrevious() == null && candidateNext.getStart() != null) { + final double distance = Vector2D.distance(end, candidateNext.getStart()); + if (distance < min) { + selectedNext = candidateNext; + min = distance; + } + } + } + if (min <= getTolerance()) { + // connect the two segments + segment.setNext(selectedNext); + selectedNext.setPrevious(segment); + ++connected; + } + } + } + return connected; + } + + /** Get first unprocessed segment from a list. + * @param segments segments list + * @return first segment that has not been processed yet + * or null if all segments have been processed + */ + private ConnectableSegment getUnprocessed(final List<ConnectableSegment> segments) { + for (final ConnectableSegment segment : segments) { + if (!segment.isProcessed()) { + return segment; + } + } + return null; + } + + /** Build the loop containing a segment. + * <p> + * The segment put in the loop will be marked as processed. + * </p> + * @param defining segment used to define the loop + * @return loop containing the segment (may be null if the loop is a + * degenerated infinitely thin 2 points loop + */ + private List<Segment> followLoop(final ConnectableSegment defining) { + + final List<Segment> loop = new ArrayList<Segment>(); + loop.add(defining); + defining.setProcessed(true); + + // add segments in connection order + ConnectableSegment next = defining.getNext(); + while (next != defining && next != null) { + loop.add(next); + next.setProcessed(true); + next = next.getNext(); + } + + if (next == null) { + // the loop is open and we have found its end, + // we need to find its start too + ConnectableSegment previous = defining.getPrevious(); + while (previous != null) { + loop.add(0, previous); + previous.setProcessed(true); + previous = previous.getPrevious(); + } + } + + // filter out spurious vertices + filterSpuriousVertices(loop); + + if (loop.size() == 2 && loop.get(0).getStart() != null) { + // this is a degenerated infinitely thin closed loop, we simply ignore it + return null; + } else { + return loop; + } + + } + + /** Filter out spurious vertices on straight lines (at machine precision). + * @param loop segments loop to filter (will be modified in-place) + */ + private void filterSpuriousVertices(final List<Segment> loop) { + for (int i = 0; i < loop.size(); ++i) { + final Segment previous = loop.get(i); + int j = (i + 1) % loop.size(); + final Segment next = loop.get(j); + if (next != null && + Precision.equals(previous.getLine().getAngle(), next.getLine().getAngle(), Precision.EPSILON)) { + // the vertex between the two edges is a spurious one + // replace the two segments by a single one + loop.set(j, new Segment(previous.getStart(), next.getEnd(), previous.getLine())); + loop.remove(i--); + } + } + } + + /** Private extension of Segment allowing connection. */ + private static class ConnectableSegment extends Segment { + + /** Node containing segment. */ + private final BSPTree<Euclidean2D> node; + + /** Node whose intersection with current node defines start point. */ + private final BSPTree<Euclidean2D> startNode; + + /** Node whose intersection with current node defines end point. */ + private final BSPTree<Euclidean2D> endNode; + + /** Previous segment. */ + private ConnectableSegment previous; + + /** Next segment. */ + private ConnectableSegment next; + + /** Indicator for completely processed segments. */ + private boolean processed; + + /** Build a segment. + * @param start start point of the segment + * @param end end point of the segment + * @param line line containing the segment + * @param node node containing the segment + * @param startNode node whose intersection with current node defines start point + * @param endNode node whose intersection with current node defines end point + */ + ConnectableSegment(final Vector2D start, final Vector2D end, final Line line, + final BSPTree<Euclidean2D> node, + final BSPTree<Euclidean2D> startNode, + final BSPTree<Euclidean2D> endNode) { + super(start, end, line); + this.node = node; + this.startNode = startNode; + this.endNode = endNode; + this.previous = null; + this.next = null; + this.processed = false; + } + + /** Get the node containing segment. + * @return node containing segment + */ + public BSPTree<Euclidean2D> getNode() { + return node; + } + + /** Get the node whose intersection with current node defines start point. + * @return node whose intersection with current node defines start point + */ + public BSPTree<Euclidean2D> getStartNode() { + return startNode; + } + + /** Get the node whose intersection with current node defines end point. + * @return node whose intersection with current node defines end point + */ + public BSPTree<Euclidean2D> getEndNode() { + return endNode; + } + + /** Get the previous segment. + * @return previous segment + */ + public ConnectableSegment getPrevious() { + return previous; + } + + /** Set the previous segment. + * @param previous previous segment + */ + public void setPrevious(final ConnectableSegment previous) { + this.previous = previous; + } + + /** Get the next segment. + * @return next segment + */ + public ConnectableSegment getNext() { + return next; + } + + /** Set the next segment. + * @param next previous segment + */ + public void setNext(final ConnectableSegment next) { + this.next = next; + } + + /** Set the processed flag. + * @param processed processed flag to set + */ + public void setProcessed(final boolean processed) { + this.processed = processed; + } + + /** Check if the segment has been processed. + * @return true if the segment has been processed + */ + public boolean isProcessed() { + return processed; + } + + } + + /** Visitor building segments. */ + private static class SegmentsBuilder implements BSPTreeVisitor<Euclidean2D> { + + /** Tolerance for close nodes connection. */ + private final double tolerance; + + /** Built segments. */ + private final List<ConnectableSegment> segments; + + /** Simple constructor. + * @param tolerance tolerance for close nodes connection + */ + SegmentsBuilder(final double tolerance) { + this.tolerance = tolerance; + this.segments = new ArrayList<ConnectableSegment>(); + } + + /** {@inheritDoc} */ + public Order visitOrder(final BSPTree<Euclidean2D> node) { + return Order.MINUS_SUB_PLUS; + } + + /** {@inheritDoc} */ + public void visitInternalNode(final BSPTree<Euclidean2D> node) { + @SuppressWarnings("unchecked") + final BoundaryAttribute<Euclidean2D> attribute = (BoundaryAttribute<Euclidean2D>) node.getAttribute(); + final Iterable<BSPTree<Euclidean2D>> splitters = attribute.getSplitters(); + if (attribute.getPlusOutside() != null) { + addContribution(attribute.getPlusOutside(), node, splitters, false); + } + if (attribute.getPlusInside() != null) { + addContribution(attribute.getPlusInside(), node, splitters, true); + } + } + + /** {@inheritDoc} */ + public void visitLeafNode(final BSPTree<Euclidean2D> node) { + } + + /** Add the contribution of a boundary facet. + * @param sub boundary facet + * @param node node containing segment + * @param splitters splitters for the boundary facet + * @param reversed if true, the facet has the inside on its plus side + */ + private void addContribution(final SubHyperplane<Euclidean2D> sub, + final BSPTree<Euclidean2D> node, + final Iterable<BSPTree<Euclidean2D>> splitters, + final boolean reversed) { + @SuppressWarnings("unchecked") + final AbstractSubHyperplane<Euclidean2D, Euclidean1D> absSub = + (AbstractSubHyperplane<Euclidean2D, Euclidean1D>) sub; + final Line line = (Line) sub.getHyperplane(); + final List<Interval> intervals = ((IntervalsSet) absSub.getRemainingRegion()).asList(); + for (final Interval i : intervals) { + + // find the 2D points + final Vector2D startV = Double.isInfinite(i.getInf()) ? + null : (Vector2D) line.toSpace((Point<Euclidean1D>) new Vector1D(i.getInf())); + final Vector2D endV = Double.isInfinite(i.getSup()) ? + null : (Vector2D) line.toSpace((Point<Euclidean1D>) new Vector1D(i.getSup())); + + // recover the connectivity information + final BSPTree<Euclidean2D> startN = selectClosest(startV, splitters); + final BSPTree<Euclidean2D> endN = selectClosest(endV, splitters); + + if (reversed) { + segments.add(new ConnectableSegment(endV, startV, line.getReverse(), + node, endN, startN)); + } else { + segments.add(new ConnectableSegment(startV, endV, line, + node, startN, endN)); + } + + } + } + + /** Select the node whose cut sub-hyperplane is closest to specified point. + * @param point reference point + * @param candidates candidate nodes + * @return node closest to point, or null if no node is closer than tolerance + */ + private BSPTree<Euclidean2D> selectClosest(final Vector2D point, final Iterable<BSPTree<Euclidean2D>> candidates) { + + BSPTree<Euclidean2D> selected = null; + double min = Double.POSITIVE_INFINITY; + + for (final BSPTree<Euclidean2D> node : candidates) { + final double distance = FastMath.abs(node.getCut().getHyperplane().getOffset(point)); + if (distance < min) { + selected = node; + min = distance; + } + } + + return min <= tolerance ? selected : null; + + } + + /** Get the segments. + * @return built segments + */ + public List<ConnectableSegment> getSegments() { + return segments; + } + + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Segment.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Segment.java new file mode 100644 index 0000000..2ef7f4e --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Segment.java @@ -0,0 +1,112 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.geometry.euclidean.twod; + +import org.apache.commons.math3.geometry.Point; +import org.apache.commons.math3.util.FastMath; + +/** Simple container for a two-points segment. + * @since 3.0 + */ +public class Segment { + + /** Start point of the segment. */ + private final Vector2D start; + + /** End point of the segment. */ + private final Vector2D end; + + /** Line containing the segment. */ + private final Line line; + + /** Build a segment. + * @param start start point of the segment + * @param end end point of the segment + * @param line line containing the segment + */ + public Segment(final Vector2D start, final Vector2D end, final Line line) { + this.start = start; + this.end = end; + this.line = line; + } + + /** Get the start point of the segment. + * @return start point of the segment + */ + public Vector2D getStart() { + return start; + } + + /** Get the end point of the segment. + * @return end point of the segment + */ + public Vector2D getEnd() { + return end; + } + + /** Get the line containing the segment. + * @return line containing the segment + */ + public Line getLine() { + return line; + } + + /** Calculates the shortest distance from a point to this line segment. + * <p> + * If the perpendicular extension from the point to the line does not + * cross in the bounds of the line segment, the shortest distance to + * the two end points will be returned. + * </p> + * + * Algorithm adapted from: + * <a href="http://www.codeguru.com/forum/printthread.php?s=cc8cf0596231f9a7dba4da6e77c29db3&t=194400&pp=15&page=1"> + * Thread @ Codeguru</a> + * + * @param p to check + * @return distance between the instance and the point + * @since 3.1 + */ + public double distance(final Vector2D p) { + final double deltaX = end.getX() - start.getX(); + final double deltaY = end.getY() - start.getY(); + + final double r = ((p.getX() - start.getX()) * deltaX + (p.getY() - start.getY()) * deltaY) / + (deltaX * deltaX + deltaY * deltaY); + + // r == 0 => P = startPt + // r == 1 => P = endPt + // r < 0 => P is on the backward extension of the segment + // r > 1 => P is on the forward extension of the segment + // 0 < r < 1 => P is on the segment + + // if point isn't on the line segment, just return the shortest distance to the end points + if (r < 0 || r > 1) { + final double dist1 = getStart().distance((Point<Euclidean2D>) p); + final double dist2 = getEnd().distance((Point<Euclidean2D>) p); + + return FastMath.min(dist1, dist2); + } + else { + // find point on line and see if it is in the line segment + final double px = start.getX() + r * deltaX; + final double py = start.getY() + r * deltaY; + + final Vector2D interPt = new Vector2D(px, py); + return interPt.distance((Point<Euclidean2D>) p); + } + } +} diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/SubLine.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/SubLine.java new file mode 100644 index 0000000..d930b76 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/SubLine.java @@ -0,0 +1,214 @@ +/* + * 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.euclidean.twod; + +import java.util.ArrayList; +import java.util.List; + +import org.apache.commons.math3.geometry.Point; +import org.apache.commons.math3.geometry.euclidean.oned.Euclidean1D; +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.OrientedPoint; +import org.apache.commons.math3.geometry.euclidean.oned.Vector1D; +import org.apache.commons.math3.geometry.partitioning.AbstractSubHyperplane; +import org.apache.commons.math3.geometry.partitioning.BSPTree; +import org.apache.commons.math3.geometry.partitioning.Hyperplane; +import org.apache.commons.math3.geometry.partitioning.Region; +import org.apache.commons.math3.geometry.partitioning.Region.Location; +import org.apache.commons.math3.geometry.partitioning.SubHyperplane; +import org.apache.commons.math3.util.FastMath; + +/** This class represents a sub-hyperplane for {@link Line}. + * @since 3.0 + */ +public class SubLine extends AbstractSubHyperplane<Euclidean2D, Euclidean1D> { + + /** Default value for tolerance. */ + private static final double DEFAULT_TOLERANCE = 1.0e-10; + + /** Simple constructor. + * @param hyperplane underlying hyperplane + * @param remainingRegion remaining region of the hyperplane + */ + public SubLine(final Hyperplane<Euclidean2D> hyperplane, + final Region<Euclidean1D> remainingRegion) { + super(hyperplane, remainingRegion); + } + + /** Create a sub-line from two endpoints. + * @param start start point + * @param end end point + * @param tolerance tolerance below which points are considered identical + * @since 3.3 + */ + public SubLine(final Vector2D start, final Vector2D end, final double tolerance) { + super(new Line(start, end, tolerance), buildIntervalSet(start, end, tolerance)); + } + + /** Create a sub-line from two endpoints. + * @param start start point + * @param end end point + * @deprecated as of 3.3, replaced with {@link #SubLine(Vector2D, Vector2D, double)} + */ + @Deprecated + public SubLine(final Vector2D start, final Vector2D end) { + this(start, end, DEFAULT_TOLERANCE); + } + + /** Create a sub-line from a segment. + * @param segment single segment forming the sub-line + */ + public SubLine(final Segment segment) { + super(segment.getLine(), + buildIntervalSet(segment.getStart(), segment.getEnd(), segment.getLine().getTolerance())); + } + + /** Get the endpoints of the sub-line. + * <p> + * A subline may be any arbitrary number of disjoints segments, so the endpoints + * are provided as a list of endpoint pairs. Each element of the list represents + * one segment, and each segment contains a start point at index 0 and an end point + * at index 1. If the sub-line is unbounded in the negative infinity direction, + * the start point of the first segment will have infinite coordinates. If the + * sub-line is unbounded in the positive infinity direction, the end point of the + * last segment will have infinite coordinates. So a sub-line covering the whole + * line will contain just one row and both elements of this row will have infinite + * coordinates. If the sub-line is empty, the returned list will contain 0 segments. + * </p> + * @return list of segments endpoints + */ + public List<Segment> getSegments() { + + final Line line = (Line) getHyperplane(); + final List<Interval> list = ((IntervalsSet) getRemainingRegion()).asList(); + final List<Segment> segments = new ArrayList<Segment>(list.size()); + + for (final Interval interval : list) { + final Vector2D start = line.toSpace((Point<Euclidean1D>) new Vector1D(interval.getInf())); + final Vector2D end = line.toSpace((Point<Euclidean1D>) new Vector1D(interval.getSup())); + segments.add(new Segment(start, end, line)); + } + + return segments; + + } + + /** Get the intersection of the instance and another sub-line. + * <p> + * This method is related to the {@link Line#intersection(Line) + * intersection} method in the {@link Line Line} class, but in addition + * to compute the point along infinite lines, it also checks the point + * lies on both sub-line ranges. + * </p> + * @param subLine other sub-line which may intersect instance + * @param includeEndPoints if true, endpoints are considered to belong to + * instance (i.e. they are closed sets) and may be returned, otherwise endpoints + * are considered to not belong to instance (i.e. they are open sets) and intersection + * occurring on endpoints lead to null being returned + * @return the intersection point if there is one, null if the sub-lines don't intersect + */ + public Vector2D intersection(final SubLine subLine, final boolean includeEndPoints) { + + // retrieve the underlying lines + Line line1 = (Line) getHyperplane(); + Line line2 = (Line) subLine.getHyperplane(); + + // compute the intersection on infinite line + Vector2D v2D = line1.intersection(line2); + if (v2D == null) { + return null; + } + + // check location of point with respect to first sub-line + Location loc1 = getRemainingRegion().checkPoint(line1.toSubSpace((Point<Euclidean2D>) v2D)); + + // check location of point with respect to second sub-line + Location loc2 = subLine.getRemainingRegion().checkPoint(line2.toSubSpace((Point<Euclidean2D>) v2D)); + + if (includeEndPoints) { + return ((loc1 != Location.OUTSIDE) && (loc2 != Location.OUTSIDE)) ? v2D : null; + } else { + return ((loc1 == Location.INSIDE) && (loc2 == Location.INSIDE)) ? v2D : null; + } + + } + + /** Build an interval set from two points. + * @param start start point + * @param end end point + * @param tolerance tolerance below which points are considered identical + * @return an interval set + */ + private static IntervalsSet buildIntervalSet(final Vector2D start, final Vector2D end, final double tolerance) { + final Line line = new Line(start, end, tolerance); + return new IntervalsSet(line.toSubSpace((Point<Euclidean2D>) start).getX(), + line.toSubSpace((Point<Euclidean2D>) end).getX(), + tolerance); + } + + /** {@inheritDoc} */ + @Override + protected AbstractSubHyperplane<Euclidean2D, Euclidean1D> buildNew(final Hyperplane<Euclidean2D> hyperplane, + final Region<Euclidean1D> remainingRegion) { + return new SubLine(hyperplane, remainingRegion); + } + + /** {@inheritDoc} */ + @Override + public SplitSubHyperplane<Euclidean2D> split(final Hyperplane<Euclidean2D> hyperplane) { + + final Line thisLine = (Line) getHyperplane(); + final Line otherLine = (Line) hyperplane; + final Vector2D crossing = thisLine.intersection(otherLine); + final double tolerance = thisLine.getTolerance(); + + if (crossing == null) { + // the lines are parallel + final double global = otherLine.getOffset(thisLine); + if (global < -tolerance) { + return new SplitSubHyperplane<Euclidean2D>(null, this); + } else if (global > tolerance) { + return new SplitSubHyperplane<Euclidean2D>(this, null); + } else { + return new SplitSubHyperplane<Euclidean2D>(null, null); + } + } + + // the lines do intersect + final boolean direct = FastMath.sin(thisLine.getAngle() - otherLine.getAngle()) < 0; + final Vector1D x = thisLine.toSubSpace((Point<Euclidean2D>) crossing); + final SubHyperplane<Euclidean1D> subPlus = + new OrientedPoint(x, !direct, tolerance).wholeHyperplane(); + final SubHyperplane<Euclidean1D> subMinus = + new OrientedPoint(x, direct, tolerance).wholeHyperplane(); + + final BSPTree<Euclidean1D> splitTree = getRemainingRegion().getTree(false).split(subMinus); + final BSPTree<Euclidean1D> plusTree = getRemainingRegion().isEmpty(splitTree.getPlus()) ? + new BSPTree<Euclidean1D>(Boolean.FALSE) : + new BSPTree<Euclidean1D>(subPlus, new BSPTree<Euclidean1D>(Boolean.FALSE), + splitTree.getPlus(), null); + final BSPTree<Euclidean1D> minusTree = getRemainingRegion().isEmpty(splitTree.getMinus()) ? + new BSPTree<Euclidean1D>(Boolean.FALSE) : + new BSPTree<Euclidean1D>(subMinus, new BSPTree<Euclidean1D>(Boolean.FALSE), + splitTree.getMinus(), null); + return new SplitSubHyperplane<Euclidean2D>(new SubLine(thisLine.copySelf(), new IntervalsSet(plusTree, tolerance)), + new SubLine(thisLine.copySelf(), new IntervalsSet(minusTree, tolerance))); + + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Vector2D.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Vector2D.java new file mode 100644 index 0000000..191d225 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Vector2D.java @@ -0,0 +1,460 @@ +/* + * 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.euclidean.twod; + +import java.text.NumberFormat; + +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.MathArithmeticException; +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; +import org.apache.commons.math3.util.MathArrays; +import org.apache.commons.math3.util.MathUtils; + +/** This class represents a 2D vector. + * <p>Instances of this class are guaranteed to be immutable.</p> + * @since 3.0 + */ +public class Vector2D implements Vector<Euclidean2D> { + + /** Origin (coordinates: 0, 0). */ + public static final Vector2D ZERO = new Vector2D(0, 0); + + // CHECKSTYLE: stop ConstantName + /** A vector with all coordinates set to NaN. */ + public static final Vector2D NaN = new Vector2D(Double.NaN, Double.NaN); + // CHECKSTYLE: resume ConstantName + + /** A vector with all coordinates set to positive infinity. */ + public static final Vector2D POSITIVE_INFINITY = + new Vector2D(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY); + + /** A vector with all coordinates set to negative infinity. */ + public static final Vector2D NEGATIVE_INFINITY = + new Vector2D(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY); + + /** Serializable UID. */ + private static final long serialVersionUID = 266938651998679754L; + + /** Abscissa. */ + private final double x; + + /** Ordinate. */ + private final double y; + + /** Simple constructor. + * Build a vector from its coordinates + * @param x abscissa + * @param y ordinate + * @see #getX() + * @see #getY() + */ + public Vector2D(double x, double y) { + this.x = x; + this.y = y; + } + + /** Simple constructor. + * Build a vector from its coordinates + * @param v coordinates array + * @exception DimensionMismatchException if array does not have 2 elements + * @see #toArray() + */ + public Vector2D(double[] v) throws DimensionMismatchException { + if (v.length != 2) { + throw new DimensionMismatchException(v.length, 2); + } + this.x = v[0]; + this.y = v[1]; + } + + /** Multiplicative constructor + * Build a vector from another one and a scale factor. + * The vector built will be a * u + * @param a scale factor + * @param u base (unscaled) vector + */ + public Vector2D(double a, Vector2D u) { + this.x = a * u.x; + this.y = a * u.y; + } + + /** Linear constructor + * Build a vector from two other ones and corresponding scale factors. + * The vector built will be a1 * u1 + a2 * u2 + * @param a1 first scale factor + * @param u1 first base (unscaled) vector + * @param a2 second scale factor + * @param u2 second base (unscaled) vector + */ + public Vector2D(double a1, Vector2D u1, double a2, Vector2D u2) { + this.x = a1 * u1.x + a2 * u2.x; + this.y = a1 * u1.y + a2 * u2.y; + } + + /** Linear constructor + * Build a vector from three other ones and corresponding scale factors. + * The vector built will be a1 * u1 + a2 * u2 + a3 * u3 + * @param a1 first scale factor + * @param u1 first base (unscaled) vector + * @param a2 second scale factor + * @param u2 second base (unscaled) vector + * @param a3 third scale factor + * @param u3 third base (unscaled) vector + */ + public Vector2D(double a1, Vector2D u1, double a2, Vector2D u2, + double a3, Vector2D u3) { + this.x = a1 * u1.x + a2 * u2.x + a3 * u3.x; + this.y = a1 * u1.y + a2 * u2.y + a3 * u3.y; + } + + /** Linear constructor + * Build a vector from four other ones and corresponding scale factors. + * The vector built will be a1 * u1 + a2 * u2 + a3 * u3 + a4 * u4 + * @param a1 first scale factor + * @param u1 first base (unscaled) vector + * @param a2 second scale factor + * @param u2 second base (unscaled) vector + * @param a3 third scale factor + * @param u3 third base (unscaled) vector + * @param a4 fourth scale factor + * @param u4 fourth base (unscaled) vector + */ + public Vector2D(double a1, Vector2D u1, double a2, Vector2D u2, + double a3, Vector2D u3, double a4, Vector2D u4) { + this.x = a1 * u1.x + a2 * u2.x + a3 * u3.x + a4 * u4.x; + this.y = a1 * u1.y + a2 * u2.y + a3 * u3.y + a4 * u4.y; + } + + /** Get the abscissa of the vector. + * @return abscissa of the vector + * @see #Vector2D(double, double) + */ + public double getX() { + return x; + } + + /** Get the ordinate of the vector. + * @return ordinate of the vector + * @see #Vector2D(double, double) + */ + public double getY() { + return y; + } + + /** Get the vector coordinates as a dimension 2 array. + * @return vector coordinates + * @see #Vector2D(double[]) + */ + public double[] toArray() { + return new double[] { x, y }; + } + + /** {@inheritDoc} */ + public Space getSpace() { + return Euclidean2D.getInstance(); + } + + /** {@inheritDoc} */ + public Vector2D getZero() { + return ZERO; + } + + /** {@inheritDoc} */ + public double getNorm1() { + return FastMath.abs(x) + FastMath.abs(y); + } + + /** {@inheritDoc} */ + public double getNorm() { + return FastMath.sqrt (x * x + y * y); + } + + /** {@inheritDoc} */ + public double getNormSq() { + return x * x + y * y; + } + + /** {@inheritDoc} */ + public double getNormInf() { + return FastMath.max(FastMath.abs(x), FastMath.abs(y)); + } + + /** {@inheritDoc} */ + public Vector2D add(Vector<Euclidean2D> v) { + Vector2D v2 = (Vector2D) v; + return new Vector2D(x + v2.getX(), y + v2.getY()); + } + + /** {@inheritDoc} */ + public Vector2D add(double factor, Vector<Euclidean2D> v) { + Vector2D v2 = (Vector2D) v; + return new Vector2D(x + factor * v2.getX(), y + factor * v2.getY()); + } + + /** {@inheritDoc} */ + public Vector2D subtract(Vector<Euclidean2D> p) { + Vector2D p3 = (Vector2D) p; + return new Vector2D(x - p3.x, y - p3.y); + } + + /** {@inheritDoc} */ + public Vector2D subtract(double factor, Vector<Euclidean2D> v) { + Vector2D v2 = (Vector2D) v; + return new Vector2D(x - factor * v2.getX(), y - factor * v2.getY()); + } + + /** {@inheritDoc} */ + public Vector2D normalize() throws MathArithmeticException { + double s = getNorm(); + if (s == 0) { + throw new MathArithmeticException(LocalizedFormats.CANNOT_NORMALIZE_A_ZERO_NORM_VECTOR); + } + return scalarMultiply(1 / s); + } + + /** Compute the angular separation between two vectors. + * <p>This method computes the angular separation between two + * vectors using the dot product for well separated vectors and the + * cross product for almost aligned vectors. This allows to have a + * good accuracy in all cases, even for vectors very close to each + * other.</p> + * @param v1 first vector + * @param v2 second vector + * @return angular separation between v1 and v2 + * @exception MathArithmeticException if either vector has a null norm + */ + public static double angle(Vector2D v1, Vector2D v2) throws MathArithmeticException { + + double normProduct = v1.getNorm() * v2.getNorm(); + if (normProduct == 0) { + throw new MathArithmeticException(LocalizedFormats.ZERO_NORM); + } + + double dot = v1.dotProduct(v2); + double threshold = normProduct * 0.9999; + if ((dot < -threshold) || (dot > threshold)) { + // the vectors are almost aligned, compute using the sine + final double n = FastMath.abs(MathArrays.linearCombination(v1.x, v2.y, -v1.y, v2.x)); + if (dot >= 0) { + return FastMath.asin(n / normProduct); + } + return FastMath.PI - FastMath.asin(n / normProduct); + } + + // the vectors are sufficiently separated to use the cosine + return FastMath.acos(dot / normProduct); + + } + + /** {@inheritDoc} */ + public Vector2D negate() { + return new Vector2D(-x, -y); + } + + /** {@inheritDoc} */ + public Vector2D scalarMultiply(double a) { + return new Vector2D(a * x, a * y); + } + + /** {@inheritDoc} */ + public boolean isNaN() { + return Double.isNaN(x) || Double.isNaN(y); + } + + /** {@inheritDoc} */ + public boolean isInfinite() { + return !isNaN() && (Double.isInfinite(x) || Double.isInfinite(y)); + } + + /** {@inheritDoc} */ + public double distance1(Vector<Euclidean2D> p) { + Vector2D p3 = (Vector2D) p; + final double dx = FastMath.abs(p3.x - x); + final double dy = FastMath.abs(p3.y - y); + return dx + dy; + } + + /** {@inheritDoc} + */ + public double distance(Vector<Euclidean2D> p) { + return distance((Point<Euclidean2D>) p); + } + + /** {@inheritDoc} */ + public double distance(Point<Euclidean2D> p) { + Vector2D p3 = (Vector2D) p; + final double dx = p3.x - x; + final double dy = p3.y - y; + return FastMath.sqrt(dx * dx + dy * dy); + } + + /** {@inheritDoc} */ + public double distanceInf(Vector<Euclidean2D> p) { + Vector2D p3 = (Vector2D) p; + final double dx = FastMath.abs(p3.x - x); + final double dy = FastMath.abs(p3.y - y); + return FastMath.max(dx, dy); + } + + /** {@inheritDoc} */ + public double distanceSq(Vector<Euclidean2D> p) { + Vector2D p3 = (Vector2D) p; + final double dx = p3.x - x; + final double dy = p3.y - y; + return dx * dx + dy * dy; + } + + /** {@inheritDoc} */ + public double dotProduct(final Vector<Euclidean2D> v) { + final Vector2D v2 = (Vector2D) v; + return MathArrays.linearCombination(x, v2.x, y, v2.y); + } + + /** + * Compute the cross-product of the instance and the given points. + * <p> + * The cross product can be used to determine the location of a point + * with regard to the line formed by (p1, p2) and is calculated as: + * \[ + * P = (x_2 - x_1)(y_3 - y_1) - (y_2 - y_1)(x_3 - x_1) + * \] + * with \(p3 = (x_3, y_3)\) being this instance. + * <p> + * If the result is 0, the points are collinear, i.e. lie on a single straight line L; + * if it is positive, this point lies to the left, otherwise to the right of the line + * formed by (p1, p2). + * + * @param p1 first point of the line + * @param p2 second point of the line + * @return the cross-product + * + * @see <a href="http://en.wikipedia.org/wiki/Cross_product">Cross product (Wikipedia)</a> + */ + public double crossProduct(final Vector2D p1, final Vector2D p2) { + final double x1 = p2.getX() - p1.getX(); + final double y1 = getY() - p1.getY(); + final double x2 = getX() - p1.getX(); + final double y2 = p2.getY() - p1.getY(); + return MathArrays.linearCombination(x1, y1, -x2, y2); + } + + /** Compute the distance between two vectors according to the L<sub>2</sub> norm. + * <p>Calling this method is equivalent to calling: + * <code>p1.subtract(p2).getNorm()</code> except that no intermediate + * vector is built</p> + * @param p1 first vector + * @param p2 second vector + * @return the distance between p1 and p2 according to the L<sub>2</sub> norm + */ + public static double distance(Vector2D p1, Vector2D p2) { + return p1.distance(p2); + } + + /** Compute the distance between two vectors according to the L<sub>∞</sub> norm. + * <p>Calling this method is equivalent to calling: + * <code>p1.subtract(p2).getNormInf()</code> except that no intermediate + * vector is built</p> + * @param p1 first vector + * @param p2 second vector + * @return the distance between p1 and p2 according to the L<sub>∞</sub> norm + */ + public static double distanceInf(Vector2D p1, Vector2D p2) { + return p1.distanceInf(p2); + } + + /** Compute the square of the distance between two vectors. + * <p>Calling this method is equivalent to calling: + * <code>p1.subtract(p2).getNormSq()</code> except that no intermediate + * vector is built</p> + * @param p1 first vector + * @param p2 second vector + * @return the square of the distance between p1 and p2 + */ + public static double distanceSq(Vector2D p1, Vector2D p2) { + return p1.distanceSq(p2); + } + + /** + * Test for the equality of two 2D vectors. + * <p> + * If all coordinates of two 2D vectors are exactly the same, and none are + * <code>Double.NaN</code>, the two 2D vectors are considered to be equal. + * </p> + * <p> + * <code>NaN</code> coordinates are considered to affect globally the vector + * and be equals to each other - i.e, if either (or all) coordinates of the + * 2D vector are equal to <code>Double.NaN</code>, the 2D vector is equal to + * {@link #NaN}. + * </p> + * + * @param other Object to test for equality to this + * @return true if two 2D vector objects are equal, false if + * object is null, not an instance of Vector2D, or + * not equal to this Vector2D instance + * + */ + @Override + public boolean equals(Object other) { + + if (this == other) { + return true; + } + + if (other instanceof Vector2D) { + final Vector2D rhs = (Vector2D)other; + if (rhs.isNaN()) { + return this.isNaN(); + } + + return (x == rhs.x) && (y == rhs.y); + } + return false; + } + + /** + * Get a hashCode for the 2D vector. + * <p> + * All NaN values have the same hash code.</p> + * + * @return a hash code value for this object + */ + @Override + public int hashCode() { + if (isNaN()) { + return 542; + } + return 122 * (76 * MathUtils.hash(x) + MathUtils.hash(y)); + } + + /** Get a string representation of this vector. + * @return a string representation of this vector + */ + @Override + public String toString() { + return Vector2DFormat.getInstance().format(this); + } + + /** {@inheritDoc} */ + public String toString(final NumberFormat format) { + return new Vector2DFormat(format).format(this); + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Vector2DFormat.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Vector2DFormat.java new file mode 100644 index 0000000..21261c5 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Vector2DFormat.java @@ -0,0 +1,138 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math3.geometry.euclidean.twod; + +import java.text.FieldPosition; +import java.text.NumberFormat; +import java.text.ParsePosition; +import java.util.Locale; + +import org.apache.commons.math3.exception.MathParseException; +import org.apache.commons.math3.geometry.Vector; +import org.apache.commons.math3.geometry.VectorFormat; +import org.apache.commons.math3.util.CompositeFormat; + +/** + * Formats a 2D vector in components list format "{x; y}". + * <p>The prefix and suffix "{" and "}" and the separator "; " can be replaced by + * any user-defined strings. The number format for components can be configured.</p> + * <p>White space is ignored at parse time, even if it is in the prefix, suffix + * or separator specifications. So even if the default separator does include a space + * character that is used at format time, both input string "{1;1}" and + * " { 1 ; 1 } " will be parsed without error and the same vector will be + * returned. In the second case, however, the parse position after parsing will be + * just after the closing curly brace, i.e. just before the trailing space.</p> + * <p><b>Note:</b> using "," as a separator may interfere with the grouping separator + * of the default {@link NumberFormat} for the current locale. Thus it is advised + * to use a {@link NumberFormat} instance with disabled grouping in such a case.</p> + * + * @since 3.0 + */ +public class Vector2DFormat extends VectorFormat<Euclidean2D> { + + /** + * Create an instance with default settings. + * <p>The instance uses the default prefix, suffix and separator: + * "{", "}", and "; " and the default number format for components.</p> + */ + public Vector2DFormat() { + super(DEFAULT_PREFIX, DEFAULT_SUFFIX, DEFAULT_SEPARATOR, + CompositeFormat.getDefaultNumberFormat()); + } + + /** + * Create an instance with a custom number format for components. + * @param format the custom format for components. + */ + public Vector2DFormat(final NumberFormat format) { + super(DEFAULT_PREFIX, DEFAULT_SUFFIX, DEFAULT_SEPARATOR, format); + } + + /** + * Create an instance with custom prefix, suffix and separator. + * @param prefix prefix to use instead of the default "{" + * @param suffix suffix to use instead of the default "}" + * @param separator separator to use instead of the default "; " + */ + public Vector2DFormat(final String prefix, final String suffix, + final String separator) { + super(prefix, suffix, separator, CompositeFormat.getDefaultNumberFormat()); + } + + /** + * Create an instance with custom prefix, suffix, separator and format + * for components. + * @param prefix prefix to use instead of the default "{" + * @param suffix suffix to use instead of the default "}" + * @param separator separator to use instead of the default "; " + * @param format the custom format for components. + */ + public Vector2DFormat(final String prefix, final String suffix, + final String separator, final NumberFormat format) { + super(prefix, suffix, separator, format); + } + + /** + * Returns the default 2D vector format for the current locale. + * @return the default 2D vector format. + */ + public static Vector2DFormat getInstance() { + return getInstance(Locale.getDefault()); + } + + /** + * Returns the default 2D vector format for the given locale. + * @param locale the specific locale used by the format. + * @return the 2D vector format specific to the given locale. + */ + public static Vector2DFormat getInstance(final Locale locale) { + return new Vector2DFormat(CompositeFormat.getDefaultNumberFormat(locale)); + } + + /** {@inheritDoc} */ + @Override + public StringBuffer format(final Vector<Euclidean2D> vector, final StringBuffer toAppendTo, + final FieldPosition pos) { + final Vector2D p2 = (Vector2D) vector; + return format(toAppendTo, pos, p2.getX(), p2.getY()); + } + + /** {@inheritDoc} */ + @Override + public Vector2D parse(final String source) throws MathParseException { + ParsePosition parsePosition = new ParsePosition(0); + Vector2D result = parse(source, parsePosition); + if (parsePosition.getIndex() == 0) { + throw new MathParseException(source, + parsePosition.getErrorIndex(), + Vector2D.class); + } + return result; + } + + /** {@inheritDoc} */ + @Override + public Vector2D parse(final String source, final ParsePosition pos) { + final double[] coordinates = parseCoordinates(2, source, pos); + if (coordinates == null) { + return null; + } + return new Vector2D(coordinates[0], coordinates[1]); + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/AbstractConvexHullGenerator2D.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/AbstractConvexHullGenerator2D.java new file mode 100644 index 0000000..b234ad5 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/AbstractConvexHullGenerator2D.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.euclidean.twod.hull; + +import java.util.Collection; + +import org.apache.commons.math3.exception.ConvergenceException; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.geometry.euclidean.twod.Vector2D; +import org.apache.commons.math3.util.MathUtils; + +/** + * Abstract base class for convex hull generators in the two-dimensional euclidean space. + * + * @since 3.3 + */ +abstract class AbstractConvexHullGenerator2D implements ConvexHullGenerator2D { + + /** Default value for tolerance. */ + private static final double DEFAULT_TOLERANCE = 1e-10; + + /** Tolerance below which points are considered identical. */ + private final double tolerance; + + /** + * Indicates if collinear points on the hull shall be present in the output. + * If {@code false}, only the extreme points are added to the hull. + */ + private final boolean includeCollinearPoints; + + /** + * Simple constructor. + * <p> + * The default tolerance (1e-10) will be used to determine identical points. + * + * @param includeCollinearPoints indicates if collinear points on the hull shall be + * added as hull vertices + */ + protected AbstractConvexHullGenerator2D(final boolean includeCollinearPoints) { + this(includeCollinearPoints, DEFAULT_TOLERANCE); + } + + /** + * Simple constructor. + * + * @param includeCollinearPoints indicates if collinear points on the hull shall be + * added as hull vertices + * @param tolerance tolerance below which points are considered identical + */ + protected AbstractConvexHullGenerator2D(final boolean includeCollinearPoints, final double tolerance) { + this.includeCollinearPoints = includeCollinearPoints; + this.tolerance = tolerance; + } + + /** + * Get the tolerance below which points are considered identical. + * @return the tolerance below which points are considered identical + */ + public double getTolerance() { + return tolerance; + } + + /** + * Returns if collinear points on the hull will be added as hull vertices. + * @return {@code true} if collinear points are added as hull vertices, or {@code false} + * if only extreme points are present. + */ + public boolean isIncludeCollinearPoints() { + return includeCollinearPoints; + } + + /** {@inheritDoc} */ + public ConvexHull2D generate(final Collection<Vector2D> points) + throws NullArgumentException, ConvergenceException { + // check for null points + MathUtils.checkNotNull(points); + + Collection<Vector2D> hullVertices = null; + if (points.size() < 2) { + hullVertices = points; + } else { + hullVertices = findHullVertices(points); + } + + try { + return new ConvexHull2D(hullVertices.toArray(new Vector2D[hullVertices.size()]), + tolerance); + } catch (MathIllegalArgumentException e) { + // the hull vertices may not form a convex hull if the tolerance value is to large + throw new ConvergenceException(); + } + } + + /** + * Find the convex hull vertices from the set of input points. + * @param points the set of input points + * @return the convex hull vertices in CCW winding + */ + protected abstract Collection<Vector2D> findHullVertices(Collection<Vector2D> points); + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/AklToussaintHeuristic.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/AklToussaintHeuristic.java new file mode 100644 index 0000000..f5d1b84 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/AklToussaintHeuristic.java @@ -0,0 +1,153 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.geometry.euclidean.twod.hull; + +import java.util.ArrayList; +import java.util.Collection; +import java.util.List; + +import org.apache.commons.math3.geometry.euclidean.twod.Vector2D; + +/** + * A simple heuristic to improve the performance of convex hull algorithms. + * <p> + * The heuristic is based on the idea of a convex quadrilateral, which is formed by + * four points with the lowest and highest x / y coordinates. Any point that lies inside + * this quadrilateral can not be part of the convex hull and can thus be safely discarded + * before generating the convex hull itself. + * <p> + * The complexity of the operation is O(n), and may greatly improve the time it takes to + * construct the convex hull afterwards, depending on the point distribution. + * + * @see <a href="http://en.wikipedia.org/wiki/Convex_hull_algorithms#Akl-Toussaint_heuristic"> + * Akl-Toussaint heuristic (Wikipedia)</a> + * @since 3.3 + */ +public final class AklToussaintHeuristic { + + /** Hide utility constructor. */ + private AklToussaintHeuristic() { + } + + /** + * Returns a point set that is reduced by all points for which it is safe to assume + * that they are not part of the convex hull. + * + * @param points the original point set + * @return a reduced point set, useful as input for convex hull algorithms + */ + public static Collection<Vector2D> reducePoints(final Collection<Vector2D> points) { + + // find the leftmost point + int size = 0; + Vector2D minX = null; + Vector2D maxX = null; + Vector2D minY = null; + Vector2D maxY = null; + for (Vector2D p : points) { + if (minX == null || p.getX() < minX.getX()) { + minX = p; + } + if (maxX == null || p.getX() > maxX.getX()) { + maxX = p; + } + if (minY == null || p.getY() < minY.getY()) { + minY = p; + } + if (maxY == null || p.getY() > maxY.getY()) { + maxY = p; + } + size++; + } + + if (size < 4) { + return points; + } + + final List<Vector2D> quadrilateral = buildQuadrilateral(minY, maxX, maxY, minX); + // if the quadrilateral is not well formed, e.g. only 2 points, do not attempt to reduce + if (quadrilateral.size() < 3) { + return points; + } + + final List<Vector2D> reducedPoints = new ArrayList<Vector2D>(quadrilateral); + for (final Vector2D p : points) { + // check all points if they are within the quadrilateral + // in which case they can not be part of the convex hull + if (!insideQuadrilateral(p, quadrilateral)) { + reducedPoints.add(p); + } + } + + return reducedPoints; + } + + /** + * Build the convex quadrilateral with the found corner points (with min/max x/y coordinates). + * + * @param points the respective points with min/max x/y coordinate + * @return the quadrilateral + */ + private static List<Vector2D> buildQuadrilateral(final Vector2D... points) { + List<Vector2D> quadrilateral = new ArrayList<Vector2D>(); + for (Vector2D p : points) { + if (!quadrilateral.contains(p)) { + quadrilateral.add(p); + } + } + return quadrilateral; + } + + /** + * Checks if the given point is located within the convex quadrilateral. + * @param point the point to check + * @param quadrilateralPoints the convex quadrilateral, represented by 4 points + * @return {@code true} if the point is inside the quadrilateral, {@code false} otherwise + */ + private static boolean insideQuadrilateral(final Vector2D point, + final List<Vector2D> quadrilateralPoints) { + + Vector2D p1 = quadrilateralPoints.get(0); + Vector2D p2 = quadrilateralPoints.get(1); + + if (point.equals(p1) || point.equals(p2)) { + return true; + } + + // get the location of the point relative to the first two vertices + final double last = point.crossProduct(p1, p2); + final int size = quadrilateralPoints.size(); + // loop through the rest of the vertices + for (int i = 1; i < size; i++) { + p1 = p2; + p2 = quadrilateralPoints.get((i + 1) == size ? 0 : i + 1); + + if (point.equals(p1) || point.equals(p2)) { + return true; + } + + // do side of line test: multiply the last location with this location + // if they are the same sign then the operation will yield a positive result + // -x * -y = +xy, x * y = +xy, -x * y = -xy, x * -y = -xy + if (last * point.crossProduct(p1, p2) < 0) { + return false; + } + } + return true; + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/ConvexHull2D.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/ConvexHull2D.java new file mode 100644 index 0000000..5d9734b --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/ConvexHull2D.java @@ -0,0 +1,172 @@ +/* + * 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.euclidean.twod.hull; + +import java.io.Serializable; + +import org.apache.commons.math3.exception.InsufficientDataException; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.geometry.euclidean.twod.Euclidean2D; +import org.apache.commons.math3.geometry.euclidean.twod.Line; +import org.apache.commons.math3.geometry.euclidean.twod.Segment; +import org.apache.commons.math3.geometry.euclidean.twod.Vector2D; +import org.apache.commons.math3.geometry.hull.ConvexHull; +import org.apache.commons.math3.geometry.partitioning.Region; +import org.apache.commons.math3.geometry.partitioning.RegionFactory; +import org.apache.commons.math3.util.MathArrays; +import org.apache.commons.math3.util.Precision; + +/** + * This class represents a convex hull in an two-dimensional euclidean space. + * + * @since 3.3 + */ +public class ConvexHull2D implements ConvexHull<Euclidean2D, Vector2D>, Serializable { + + /** Serializable UID. */ + private static final long serialVersionUID = 20140129L; + + /** Vertices of the hull. */ + private final Vector2D[] vertices; + + /** Tolerance threshold used during creation of the hull vertices. */ + private final double tolerance; + + /** + * Line segments of the hull. + * The array is not serialized and will be created from the vertices on first access. + */ + private transient Segment[] lineSegments; + + /** + * Simple constructor. + * @param vertices the vertices of the convex hull, must be ordered + * @param tolerance tolerance below which points are considered identical + * @throws MathIllegalArgumentException if the vertices do not form a convex hull + */ + public ConvexHull2D(final Vector2D[] vertices, final double tolerance) + throws MathIllegalArgumentException { + + // assign tolerance as it will be used by the isConvex method + this.tolerance = tolerance; + + if (!isConvex(vertices)) { + throw new MathIllegalArgumentException(LocalizedFormats.NOT_CONVEX); + } + + this.vertices = vertices.clone(); + } + + /** + * Checks whether the given hull vertices form a convex hull. + * @param hullVertices the hull vertices + * @return {@code true} if the vertices form a convex hull, {@code false} otherwise + */ + private boolean isConvex(final Vector2D[] hullVertices) { + if (hullVertices.length < 3) { + return true; + } + + int sign = 0; + for (int i = 0; i < hullVertices.length; i++) { + final Vector2D p1 = hullVertices[i == 0 ? hullVertices.length - 1 : i - 1]; + final Vector2D p2 = hullVertices[i]; + final Vector2D p3 = hullVertices[i == hullVertices.length - 1 ? 0 : i + 1]; + + final Vector2D d1 = p2.subtract(p1); + final Vector2D d2 = p3.subtract(p2); + + final double crossProduct = MathArrays.linearCombination(d1.getX(), d2.getY(), -d1.getY(), d2.getX()); + final int cmp = Precision.compareTo(crossProduct, 0.0, tolerance); + // in case of collinear points the cross product will be zero + if (cmp != 0.0) { + if (sign != 0.0 && cmp != sign) { + return false; + } + sign = cmp; + } + } + + return true; + } + + /** {@inheritDoc} */ + public Vector2D[] getVertices() { + return vertices.clone(); + } + + /** + * Get the line segments of the convex hull, ordered. + * @return the line segments of the convex hull + */ + public Segment[] getLineSegments() { + return retrieveLineSegments().clone(); + } + + /** + * Retrieve the line segments from the cached array or create them if needed. + * + * @return the array of line segments + */ + private Segment[] retrieveLineSegments() { + if (lineSegments == null) { + // construct the line segments - handle special cases of 1 or 2 points + final int size = vertices.length; + if (size <= 1) { + this.lineSegments = new Segment[0]; + } else if (size == 2) { + this.lineSegments = new Segment[1]; + final Vector2D p1 = vertices[0]; + final Vector2D p2 = vertices[1]; + this.lineSegments[0] = new Segment(p1, p2, new Line(p1, p2, tolerance)); + } else { + this.lineSegments = new Segment[size]; + Vector2D firstPoint = null; + Vector2D lastPoint = null; + int index = 0; + for (Vector2D point : vertices) { + if (lastPoint == null) { + firstPoint = point; + lastPoint = point; + } else { + this.lineSegments[index++] = + new Segment(lastPoint, point, new Line(lastPoint, point, tolerance)); + lastPoint = point; + } + } + this.lineSegments[index] = + new Segment(lastPoint, firstPoint, new Line(lastPoint, firstPoint, tolerance)); + } + } + return lineSegments; + } + + /** {@inheritDoc} */ + public Region<Euclidean2D> createRegion() throws InsufficientDataException { + if (vertices.length < 3) { + throw new InsufficientDataException(); + } + final RegionFactory<Euclidean2D> factory = new RegionFactory<Euclidean2D>(); + final Segment[] segments = retrieveLineSegments(); + final Line[] lineArray = new Line[segments.length]; + for (int i = 0; i < segments.length; i++) { + lineArray[i] = segments[i].getLine(); + } + return factory.buildConvex(lineArray); + } +} diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/ConvexHullGenerator2D.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/ConvexHullGenerator2D.java new file mode 100644 index 0000000..3e13e1a --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/ConvexHullGenerator2D.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.euclidean.twod.hull; + +import java.util.Collection; + +import org.apache.commons.math3.exception.ConvergenceException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.geometry.euclidean.twod.Euclidean2D; +import org.apache.commons.math3.geometry.euclidean.twod.Vector2D; +import org.apache.commons.math3.geometry.hull.ConvexHullGenerator; + +/** + * Interface for convex hull generators in the two-dimensional euclidean space. + * + * @since 3.3 + */ +public interface ConvexHullGenerator2D extends ConvexHullGenerator<Euclidean2D, Vector2D> { + + /** {@inheritDoc} */ + ConvexHull2D generate(Collection<Vector2D> points) throws NullArgumentException, ConvergenceException; + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/MonotoneChain.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/MonotoneChain.java new file mode 100644 index 0000000..4421344 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/MonotoneChain.java @@ -0,0 +1,181 @@ +/* + * 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.euclidean.twod.hull; + +import java.util.ArrayList; +import java.util.Collection; +import java.util.Collections; +import java.util.Comparator; +import java.util.List; + +import org.apache.commons.math3.geometry.euclidean.twod.Line; +import org.apache.commons.math3.geometry.euclidean.twod.Vector2D; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.Precision; + +/** + * Implements Andrew's monotone chain method to generate the convex hull of a finite set of + * points in the two-dimensional euclidean space. + * <p> + * The runtime complexity is O(n log n), with n being the number of input points. If the + * point set is already sorted (by x-coordinate), the runtime complexity is O(n). + * <p> + * The implementation is not sensitive to collinear points on the hull. The parameter + * {@code includeCollinearPoints} allows to control the behavior with regard to collinear points. + * If {@code true}, all points on the boundary of the hull will be added to the hull vertices, + * otherwise only the extreme points will be present. By default, collinear points are not added + * as hull vertices. + * <p> + * The {@code tolerance} parameter (default: 1e-10) is used as epsilon criteria to determine + * identical and collinear points. + * + * @see <a href="http://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain"> + * Andrew's monotone chain algorithm (Wikibooks)</a> + * @since 3.3 + */ +public class MonotoneChain extends AbstractConvexHullGenerator2D { + + /** + * Create a new MonotoneChain instance. + */ + public MonotoneChain() { + this(false); + } + + /** + * Create a new MonotoneChain instance. + * @param includeCollinearPoints whether collinear points shall be added as hull vertices + */ + public MonotoneChain(final boolean includeCollinearPoints) { + super(includeCollinearPoints); + } + + /** + * Create a new MonotoneChain instance. + * @param includeCollinearPoints whether collinear points shall be added as hull vertices + * @param tolerance tolerance below which points are considered identical + */ + public MonotoneChain(final boolean includeCollinearPoints, final double tolerance) { + super(includeCollinearPoints, tolerance); + } + + /** {@inheritDoc} */ + @Override + public Collection<Vector2D> findHullVertices(final Collection<Vector2D> points) { + + final List<Vector2D> pointsSortedByXAxis = new ArrayList<Vector2D>(points); + + // sort the points in increasing order on the x-axis + Collections.sort(pointsSortedByXAxis, new Comparator<Vector2D>() { + /** {@inheritDoc} */ + public int compare(final Vector2D o1, final Vector2D o2) { + final double tolerance = getTolerance(); + // need to take the tolerance value into account, otherwise collinear points + // will not be handled correctly when building the upper/lower hull + final int diff = Precision.compareTo(o1.getX(), o2.getX(), tolerance); + if (diff == 0) { + return Precision.compareTo(o1.getY(), o2.getY(), tolerance); + } else { + return diff; + } + } + }); + + // build lower hull + final List<Vector2D> lowerHull = new ArrayList<Vector2D>(); + for (Vector2D p : pointsSortedByXAxis) { + updateHull(p, lowerHull); + } + + // build upper hull + final List<Vector2D> upperHull = new ArrayList<Vector2D>(); + for (int idx = pointsSortedByXAxis.size() - 1; idx >= 0; idx--) { + final Vector2D p = pointsSortedByXAxis.get(idx); + updateHull(p, upperHull); + } + + // concatenate the lower and upper hulls + // the last point of each list is omitted as it is repeated at the beginning of the other list + final List<Vector2D> hullVertices = new ArrayList<Vector2D>(lowerHull.size() + upperHull.size() - 2); + for (int idx = 0; idx < lowerHull.size() - 1; idx++) { + hullVertices.add(lowerHull.get(idx)); + } + for (int idx = 0; idx < upperHull.size() - 1; idx++) { + hullVertices.add(upperHull.get(idx)); + } + + // special case: if the lower and upper hull may contain only 1 point if all are identical + if (hullVertices.isEmpty() && ! lowerHull.isEmpty()) { + hullVertices.add(lowerHull.get(0)); + } + + return hullVertices; + } + + /** + * Update the partial hull with the current point. + * + * @param point the current point + * @param hull the partial hull + */ + private void updateHull(final Vector2D point, final List<Vector2D> hull) { + final double tolerance = getTolerance(); + + if (hull.size() == 1) { + // ensure that we do not add an identical point + final Vector2D p1 = hull.get(0); + if (p1.distance(point) < tolerance) { + return; + } + } + + while (hull.size() >= 2) { + final int size = hull.size(); + final Vector2D p1 = hull.get(size - 2); + final Vector2D p2 = hull.get(size - 1); + + final double offset = new Line(p1, p2, tolerance).getOffset(point); + if (FastMath.abs(offset) < tolerance) { + // the point is collinear to the line (p1, p2) + + final double distanceToCurrent = p1.distance(point); + if (distanceToCurrent < tolerance || p2.distance(point) < tolerance) { + // the point is assumed to be identical to either p1 or p2 + return; + } + + final double distanceToLast = p1.distance(p2); + if (isIncludeCollinearPoints()) { + final int index = distanceToCurrent < distanceToLast ? size - 1 : size; + hull.add(index, point); + } else { + if (distanceToCurrent > distanceToLast) { + hull.remove(size - 1); + hull.add(point); + } + } + return; + } else if (offset > 0) { + hull.remove(size - 1); + } else { + break; + } + } + hull.add(point); + } + +} diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/package-info.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/package-info.java new file mode 100644 index 0000000..d0469a4 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/hull/package-info.java @@ -0,0 +1,25 @@ +/* + * 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 algorithms to generate the convex hull + * for a set of points in an two-dimensional euclidean space. + * </p> + * + */ +package org.apache.commons.math3.geometry.euclidean.twod.hull; diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/package-info.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/package-info.java new file mode 100644 index 0000000..feb43b1 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/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 basic 2D geometry components. + * </p> + * + */ +package org.apache.commons.math3.geometry.euclidean.twod; |