Improved accuracy of internal computation. The improvement is obtained by using high accuracy linear combinations and by preserving cos/sin as much as possible. This last part is important in the frequent case of lines along y axis, i.e. when angle is Ï/2. We want the cosine to remain 0 and not get about 10^-17.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/5f667c03 Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/5f667c03 Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/5f667c03 Branch: refs/heads/master Commit: 5f667c031cdf02adbee9ef83cb71990333c58ba3 Parents: c29e3d6 Author: Luc Maisonobe <l...@apache.org> Authored: Sun Nov 30 11:44:07 2014 +0100 Committer: Luc Maisonobe <l...@apache.org> Committed: Tue Dec 2 15:24:31 2014 +0100 ---------------------------------------------------------------------- .../math3/geometry/euclidean/twod/Line.java | 54 ++++++++++---------- 1 file changed, 28 insertions(+), 26 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-math/blob/5f667c03/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java ---------------------------------------------------------------------- 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 index 70ab1ac..b40f0fa 100644 --- 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 @@ -31,6 +31,7 @@ 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. @@ -147,8 +148,8 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc */ public Line(final Line line) { angle = MathUtils.normalizeAngle(line.angle, FastMath.PI); - cos = FastMath.cos(angle); - sin = FastMath.sin(angle); + cos = line.cos; + sin = line.sin; originOffset = line.originOffset; tolerance = line.tolerance; } @@ -174,9 +175,9 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc originOffset = p1.getY(); } else { angle = FastMath.PI + FastMath.atan2(-dy, -dx); - cos = FastMath.cos(angle); - sin = FastMath.sin(angle); - originOffset = (p2.getX() * p1.getY() - p1.getX() * p2.getY()) / d; + cos = dx / d; + sin = dy / d; + originOffset = MathArrays.linearCombination(p2.getX(), p1.getY(), -p1.getX(), p2.getY()) / d; } } @@ -188,7 +189,7 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc this.angle = MathUtils.normalizeAngle(alpha, FastMath.PI); cos = FastMath.cos(this.angle); sin = FastMath.sin(this.angle); - originOffset = cos * p.getY() - sin * p.getX(); + originOffset = MathArrays.linearCombination(cos, p.getY(), -sin, p.getX()); } /** Revert the instance. @@ -235,14 +236,14 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc /** {@inheritDoc} */ public Vector1D toSubSpace(final Point<Euclidean2D> point) { Vector2D p2 = (Vector2D) point; - return new Vector1D(cos * p2.getX() + sin * p2.getY()); + 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(abscissa * cos - originOffset * sin, - abscissa * sin + originOffset * cos); + 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. @@ -251,12 +252,12 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc * or null if there are no intersection points */ public Vector2D intersection(final Line other) { - final double d = sin * other.cos - other.sin * cos; + final double d = MathArrays.linearCombination(sin, other.cos, -other.sin, cos); if (FastMath.abs(d) < tolerance) { return null; } - return new Vector2D((cos * other.originOffset - other.cos * originOffset) / d, - (sin * other.originOffset - other.sin * originOffset) / d); + return new Vector2D(MathArrays.linearCombination(cos, other.originOffset, -other.cos, originOffset) / d, + MathArrays.linearCombination(sin, other.originOffset, -other.sin, originOffset) / d); } /** {@inheritDoc} @@ -298,7 +299,7 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc */ public double getOffset(final Line line) { return originOffset + - ((cos * line.cos + sin * line.sin > 0) ? -line.originOffset : line.originOffset); + (MathArrays.linearCombination(cos, line.cos, sin, line.sin) > 0 ? -line.originOffset : line.originOffset); } /** Get the offset (oriented distance) of a vector. @@ -312,13 +313,13 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc /** {@inheritDoc} */ public double getOffset(final Point<Euclidean2D> point) { Vector2D p2 = (Vector2D) point; - return sin * p2.getX() - cos * p2.getY() + originOffset; + 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 (sin * otherL.sin + cos * otherL.cos) >= 0.0; + return MathArrays.linearCombination(sin, otherL.sin, cos, otherL.cos) >= 0.0; } /** Get one point from the plane. @@ -330,7 +331,8 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc public Vector2D getPointAt(final Vector1D abscissa, final double offset) { final double x = abscissa.getX(); final double dOffset = offset - originOffset; - return new Vector2D(x * cos + dOffset * sin, x * sin - dOffset * cos); + return new Vector2D(MathArrays.linearCombination(x, cos, dOffset, sin), + MathArrays.linearCombination(x, sin, -dOffset, cos)); } /** Check if the line contains a point. @@ -360,14 +362,14 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc * (they can have either the same or opposite orientations) */ public boolean isParallelTo(final Line line) { - return FastMath.abs(sin * line.cos - cos * line.sin) < tolerance; + 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 = cos * p.getY() - sin * p.getX(); + originOffset = MathArrays.linearCombination(cos, p.getY(), -sin, p.getX()); } /** Get the angle of the line. @@ -457,9 +459,9 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc cYY = m[3]; cY1 = m[5]; - c1Y = cXY * cY1 - cYY * cX1; - c1X = cXX * cY1 - cYX * cX1; - c11 = cXX * cYY - cYX * cXY; + 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); @@ -472,16 +474,16 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc final Vector2D p2D = (Vector2D) point; final double x = p2D.getX(); final double y = p2D.getY(); - return new Vector2D(cXX * x + cXY * y + cX1, - cYX * x + cYY * y + cY1); + 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 = c1X * line.cos + c1Y * line.sin + c11 * line.originOffset; - final double rCos = cXX * line.cos + cXY * line.sin; - final double rSin = cYX * line.cos + cYY * line.sin; + final double 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,
