From 19c1c3bb9bc175f9c19b026e7cda1d0fe47ade90 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Sun, 2 Feb 2014 20:52:49 +0000 Subject: [PATCH] Fixed sphere generation in degenerated cases. In almost coplanar / almost colinear cases, sphere generation could diverge in such a way the support points did not belong to the sphere anymore! This new implementation uses exact arithmetic for the computation. It is much slower, but since very few balls should be generated during the Welzl algorithm execution, this is probably acceptable. JIRA: MATH-1096 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1563712 13f79535-47bb-0310-9956-ffa450edef68 --- .../geometry/enclosing/WelzlEncloser.java | 3 +- .../euclidean/threed/SphereGenerator.java | 91 +++++++++---------- .../euclidean/twod/DiskGenerator.java | 52 ++++++----- .../enclosing/WelzlEncloser3DTest.java | 6 +- .../euclidean/threed/SphereGeneratorTest.java | 35 ++++++- 5 files changed, 110 insertions(+), 77 deletions(-) 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)); + } + + } + }