diff --git a/src/main/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser.java b/src/main/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser.java index 1bacb2549..12a645f43 100644 --- a/src/main/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser.java +++ b/src/main/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser.java @@ -89,6 +89,7 @@ public class WelzlEncloser> implements Enclo // select the point farthest to current ball final P farthest = selectFarthest(points, ball); + if (ball.contains(farthest, tolerance)) { // we have found a ball containing all points return ball; @@ -100,7 +101,7 @@ public class WelzlEncloser> implements Enclo EnclosingBall savedBall = ball; ball = moveToFrontBall(extreme, extreme.size(), support); if (ball.getRadius() < savedBall.getRadius()) { - // TODO: fix this, it should never happen but it does! + // this should never happen throw new MathInternalError(); } diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGenerator.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGenerator.java index 03e445089..f5a6b7c95 100644 --- a/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGenerator.java +++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGenerator.java @@ -19,12 +19,13 @@ package org.apache.commons.math3.geometry.euclidean.threed; import java.util.Arrays; 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.geometry.euclidean.twod.DiskGenerator; import org.apache.commons.math3.geometry.euclidean.twod.Euclidean2D; import org.apache.commons.math3.geometry.euclidean.twod.Vector2D; -import org.apache.commons.math3.util.MathArrays; +import org.apache.commons.math3.util.FastMath; /** Class generating an enclosing ball from its support points. * @version $Id$ @@ -88,24 +89,39 @@ public class SphereGenerator implements SupportBallGenerator(center, center.distance(vA), + final BigFraction twoM11 = minor(c2, c3, c4).multiply(2); + final BigFraction m12 = minor(c1, c3, c4); + final BigFraction m13 = minor(c1, c2, c4); + final BigFraction m14 = minor(c1, c2, c3); + final BigFraction centerX = m12.divide(twoM11); + final BigFraction centerY = m13.divide(twoM11).negate(); + final BigFraction centerZ = m14.divide(twoM11); + final BigFraction dx = c2[0].subtract(centerX); + final BigFraction dy = c3[0].subtract(centerY); + final BigFraction dz = c4[0].subtract(centerZ); + final BigFraction r2 = dx.multiply(dx).add(dy.multiply(dy)).add(dz.multiply(dz)); + return new EnclosingBall(new Vector3D(centerX.doubleValue(), + centerY.doubleValue(), + centerZ.doubleValue()), + FastMath.sqrt(r2.doubleValue()), vA, vB, vC, vD); } } @@ -114,41 +130,24 @@ public class SphereGenerator implements SupportBallGeneratorth column is known to be filled with 1.0. - *

- * The computation is performed using {@link MathArrays#linearCombination(double[], double[]) - * high accuracy sum of products}, trying to avoid cancellations effect. This should reduce - * risks in case of near co-planar points. - *

* @param c1 first column * @param c2 second column * @param c3 third column - * @return value of the minor computed to high accuracy + * @return value of the minor computed has an exact fraction */ - private double minor(final double[] c1, final double[] c2, final double[] c3) { - final double m01 = c2[0] * c3[1]; - final double m02 = c2[0] * c3[2]; - final double m03 = c2[0] * c3[3]; - final double m10 = c2[1] * c3[0]; - final double m12 = c2[1] * c3[2]; - final double m13 = c2[1] * c3[3]; - final double m20 = c2[2] * c3[0]; - final double m21 = c2[2] * c3[1]; - final double m23 = c2[2] * c3[3]; - final double m30 = c2[3] * c3[0]; - final double m31 = c2[3] * c3[1]; - final double m32 = c2[3] * c3[2]; - return MathArrays.linearCombination(new double[] { - c1[2], c1[1], c1[3], -c1[1], -c1[3], -c1[2], - c1[0], c1[3], c1[2], -c1[3], -c1[0], -c1[2], - c1[1], c1[0], c1[3], -c1[0], -c1[3], -c1[1], - c1[0], c1[2], c1[1], -c1[2], -c1[0], -c1[1] - }, - new double[] { - m13, m32, m21, m23, m12, m31, - m23, m02, m30, m20, m32, m03, - m03, m31, m10, m13, m01, m30, - m12, m01, m20, m10, m21, m02 - }); + private BigFraction minor(final BigFraction[] c1, final BigFraction[] c2, final BigFraction[] c3) { + return c2[0].multiply(c3[1]).multiply(c1[2].subtract(c1[3])). + add(c2[0].multiply(c3[2]).multiply(c1[3].subtract(c1[1]))). + add(c2[0].multiply(c3[3]).multiply(c1[1].subtract(c1[2]))). + add(c2[1].multiply(c3[0]).multiply(c1[3].subtract(c1[2]))). + add(c2[1].multiply(c3[2]).multiply(c1[0].subtract(c1[3]))). + add(c2[1].multiply(c3[3]).multiply(c1[2].subtract(c1[0]))). + add(c2[2].multiply(c3[0]).multiply(c1[1].subtract(c1[3]))). + add(c2[2].multiply(c3[1]).multiply(c1[3].subtract(c1[0]))). + add(c2[2].multiply(c3[3]).multiply(c1[0].subtract(c1[1]))). + add(c2[3].multiply(c3[0]).multiply(c1[2].subtract(c1[1]))). + add(c2[3].multiply(c3[1]).multiply(c1[0].subtract(c1[2]))). + add(c2[3].multiply(c3[2]).multiply(c1[1].subtract(c1[0]))); } } 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 index 514df1458..d06c013cd 100644 --- 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 @@ -18,9 +18,10 @@ 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.MathArrays; +import org.apache.commons.math3.util.FastMath; /** Class generating an enclosing ball from its support points. * @version $Id$ @@ -66,42 +67,43 @@ public class DiskGenerator implements SupportBallGenerator(center, center.distance(vA), vA, vB, vC); + 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(new Vector2D(centerX.doubleValue(), + centerY.doubleValue()), + FastMath.sqrt(r2.doubleValue()), + vA, vB, vC); } } } } /** Compute a dimension 3 minor, when 3d column is known to be filled with 1.0. - *

- * The computation is performed using {@link MathArrays#linearCombination(double[], double[]) - * high accuracy sum of products}, trying to avoid cancellations effect. This should reduce - * risks in case of near co-planar points. - *

* @param c1 first column * @param c2 second column - * @return value of the minor computed to high accuracy + * @return value of the minor computed has an exact fraction */ - private double minor(final double[] c1, final double[] c2) { - return MathArrays.linearCombination(new double[] { - c1[0], c1[2], c1[1], -c1[2], -c1[0], -c1[1] - }, - new double[] { - c2[1], c2[0], c2[2], c2[1], c2[2], c2[0] - }); + 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/test/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser3DTest.java b/src/test/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser3DTest.java index 068288877..b4bc2d682 100644 --- a/src/test/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser3DTest.java +++ b/src/test/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser3DTest.java @@ -28,7 +28,6 @@ import org.apache.commons.math3.random.RandomGenerator; import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator; import org.apache.commons.math3.random.Well1024a; import org.junit.Assert; -import org.junit.Ignore; import org.junit.Test; @@ -53,7 +52,6 @@ public class WelzlEncloser3DTest { } @Test - @Ignore public void testReducingBall() { List list = Arrays.asList(new Vector3D(-7.140397329568118, -16.571661242582177, 11.714458961735405), @@ -106,14 +104,14 @@ public class WelzlEncloser3DTest { public void testLargeSamples() throws IOException { RandomGenerator random = new Well1024a(0x35ddecfc78131e1dl); final UnitSphereRandomVectorGenerator sr = new UnitSphereRandomVectorGenerator(3, random); - for (int k = 0; k < 100; ++k) { + for (int k = 0; k < 50; ++k) { // define the reference sphere we want to compute double d = 25 * random.nextDouble(); double refRadius = 10 * random.nextDouble(); Vector3D refCenter = new Vector3D(d, new Vector3D(sr.nextVector())); // set up a large sample inside the reference sphere - int nbPoints = random.nextInt(10000); + int nbPoints = random.nextInt(1000); List points = new ArrayList(); for (int i = 0; i < nbPoints; ++i) { double r = refRadius * random.nextDouble(); diff --git a/src/test/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGeneratorTest.java b/src/test/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGeneratorTest.java index 877d7d415..1950a066c 100644 --- a/src/test/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGeneratorTest.java +++ b/src/test/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGeneratorTest.java @@ -24,6 +24,7 @@ import org.apache.commons.math3.geometry.enclosing.EnclosingBall; import org.apache.commons.math3.random.RandomGenerator; import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator; import org.apache.commons.math3.random.Well1024a; +import org.apache.commons.math3.util.FastMath; import org.junit.Assert; import org.junit.Test; @@ -135,7 +136,7 @@ public class SphereGeneratorTest { public void testRandom() { final RandomGenerator random = new Well1024a(0xd015982e9f31ee04l); final UnitSphereRandomVectorGenerator sr = new UnitSphereRandomVectorGenerator(3, random); - for (int i = 0; i < 500; ++i) { + for (int i = 0; i < 100; ++i) { double d = 25 * random.nextDouble(); double refRadius = 10 * random.nextDouble(); Vector3D refCenter = new Vector3D(d, new Vector3D(sr.nextVector())); @@ -150,4 +151,36 @@ public class SphereGeneratorTest { } + @Test + public void testDegeneratedCase() { + final List support = + Arrays.asList(new Vector3D(FastMath.scalb(-8039905610797991.0, -50), // -7.140870659936730 + FastMath.scalb(-4663475464714142.0, -48), // -16.567993074240455 + FastMath.scalb( 6592658872616184.0, -49)), // 11.710914678204503 + new Vector3D(FastMath.scalb(-8036658568968473.0, -50), // -7.137986707455888 + FastMath.scalb(-4664256346424880.0, -48), // -16.570767323375720 + FastMath.scalb( 6591357011730307.0, -49)), // 11.708602108715928) + new Vector3D(FastMath.scalb(-8037820142977230.0, -50), // -7.139018392423351 + FastMath.scalb(-4665280434237813.0, -48), // -16.574405614157020 + FastMath.scalb( 6592435966112099.0, -49)), // 11.710518716711425 + new Vector3D(FastMath.scalb(-8038007803611611.0, -50), // -7.139185068549035 + FastMath.scalb(-4664291215918380.0, -48), // -16.570891204702250 + FastMath.scalb( 6595270610894208.0, -49))); // 11.715554057357394 + EnclosingBall sphere = new SphereGenerator().ballOnSupport(support); + + // the following values have been computed using Emacs calc with exact arithmetic from the + // rational representation corresponding to the scalb calls (i.e. -8039905610797991/2^50, ...) + // The results were converted to decimal representation rounded to 1.0e-30 when writing the reference + // values in this test + Assert.assertEquals( 0.003616820213530053297575846168, sphere.getRadius(), 1.0e-20); + Assert.assertEquals( -7.139325643360503322823511839511, sphere.getCenter().getX(), 1.0e-20); + Assert.assertEquals(-16.571096474251747245361467833760, sphere.getCenter().getY(), 1.0e-20); + Assert.assertEquals( 11.711945804096960876521111630800, sphere.getCenter().getZ(), 1.0e-20); + + for (Vector3D v : support) { + Assert.assertTrue(sphere.contains(v, 1.0e-14)); + } + + } + }