Return-Path: X-Original-To: apmail-commons-commits-archive@minotaur.apache.org Delivered-To: apmail-commons-commits-archive@minotaur.apache.org Received: from mail.apache.org (hermes.apache.org [140.211.11.3]) by minotaur.apache.org (Postfix) with SMTP id 6CFDB10FB1 for ; Sun, 2 Feb 2014 20:53:19 +0000 (UTC) Received: (qmail 54875 invoked by uid 500); 2 Feb 2014 20:53:17 -0000 Delivered-To: apmail-commons-commits-archive@commons.apache.org Received: (qmail 54727 invoked by uid 500); 2 Feb 2014 20:53:16 -0000 Mailing-List: contact commits-help@commons.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: dev@commons.apache.org Delivered-To: mailing list commits@commons.apache.org Received: (qmail 54720 invoked by uid 99); 2 Feb 2014 20:53:16 -0000 Received: from nike.apache.org (HELO nike.apache.org) (192.87.106.230) by apache.org (qpsmtpd/0.29) with ESMTP; Sun, 02 Feb 2014 20:53:16 +0000 X-ASF-Spam-Status: No, hits=-2000.0 required=5.0 tests=ALL_TRUSTED X-Spam-Check-By: apache.org Received: from [140.211.11.4] (HELO eris.apache.org) (140.211.11.4) by apache.org (qpsmtpd/0.29) with ESMTP; Sun, 02 Feb 2014 20:53:12 +0000 Received: from eris.apache.org (localhost [127.0.0.1]) by eris.apache.org (Postfix) with ESMTP id BF40D23888D7; Sun, 2 Feb 2014 20:52:49 +0000 (UTC) Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Subject: svn commit: r1563712 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math3/geometry/enclosing/ main/java/org/apache/commons/math3/geometry/euclidean/threed/ main/java/org/apache/commons/math3/geometry/euclidean/twod/ test/java/org/apa... Date: Sun, 02 Feb 2014 20:52:49 -0000 To: commits@commons.apache.org From: luc@apache.org X-Mailer: svnmailer-1.0.9 Message-Id: <20140202205249.BF40D23888D7@eris.apache.org> X-Virus-Checked: Checked by ClamAV on apache.org Author: luc Date: Sun Feb 2 20:52:49 2014 New Revision: 1563712 URL: http://svn.apache.org/r1563712 Log: 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 Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser.java commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGenerator.java commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/DiskGenerator.java commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser3DTest.java commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGeneratorTest.java Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser.java?rev=1563712&r1=1563711&r2=1563712&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser.java Sun Feb 2 20:52:49 2014 @@ -89,6 +89,7 @@ public class WelzlEncloser 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(); } Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGenerator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGenerator.java?rev=1563712&r1=1563711&r2=1563712&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGenerator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGenerator.java Sun Feb 2 20:52:49 2014 @@ -19,12 +19,13 @@ package org.apache.commons.math3.geometr 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 // z_0 = +m_14 / (2 m_11) // Note that the minors m_11, m_12, m_13 and m_14 all have the last column // filled with 1.0, hence simplifying the computation - final double[] c1 = new double[] { - vA.getNormSq(), vB.getNormSq(), vC.getNormSq(), vD.getNormSq() + final BigFraction[] c2 = new BigFraction[] { + new BigFraction(vA.getX()), new BigFraction(vB.getX()), + new BigFraction(vC.getX()), new BigFraction(vD.getX()) }; - final double[] c2 = new double[] { - vA.getX(), vB.getX(), vC.getX(), vD.getX() + final BigFraction[] c3 = new BigFraction[] { + new BigFraction(vA.getY()), new BigFraction(vB.getY()), + new BigFraction(vC.getY()), new BigFraction(vD.getY()) }; - final double[] c3 = new double[] { - vA.getY(), vB.getY(), vC.getY(), vD.getY() + final BigFraction[] c4 = new BigFraction[] { + new BigFraction(vA.getZ()), new BigFraction(vB.getZ()), + new BigFraction(vC.getZ()), new BigFraction(vD.getZ()) }; - final double[] c4 = new double[] { - vA.getZ(), vB.getZ(), vC.getZ(), vD.getZ() + final BigFraction[] c1 = new BigFraction[] { + c2[0].multiply(c2[0]).add(c3[0].multiply(c3[0])).add(c4[0].multiply(c4[0])), + c2[1].multiply(c2[1]).add(c3[1].multiply(c3[1])).add(c4[1].multiply(c4[1])), + c2[2].multiply(c2[2]).add(c3[2].multiply(c3[2])).add(c4[2].multiply(c4[2])), + c2[3].multiply(c2[3]).add(c3[3].multiply(c3[3])).add(c4[3].multiply(c4[3])) }; - final double m11 = minor(c2, c3, c4); - final double m12 = minor(c1, c3, c4); - final double m13 = minor(c1, c2, c4); - final double m14 = minor(c1, c2, c3); - final Vector3D center = new Vector3D(0.5 * m12 / m11, -0.5 * m13 / m11, 0.5 * m14 / m11); - return new EnclosingBall(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 } /** Compute a dimension 4 minor, when 4th 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]))); } } Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/DiskGenerator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/DiskGenerator.java?rev=1563712&r1=1563711&r2=1563712&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/DiskGenerator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/DiskGenerator.java Sun Feb 2 20:52:49 2014 @@ -18,9 +18,10 @@ package org.apache.commons.math3.geometr 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 Su // 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 double[] c1 = new double[] { - vA.getNormSq(), vB.getNormSq(), vC.getNormSq() + final BigFraction[] c2 = new BigFraction[] { + new BigFraction(vA.getX()), new BigFraction(vB.getX()), new BigFraction(vC.getX()) }; - final double[] c2 = new double[] { - vA.getX(), vB.getX(), vC.getX() + final BigFraction[] c3 = new BigFraction[] { + new BigFraction(vA.getY()), new BigFraction(vB.getY()), new BigFraction(vC.getY()) }; - final double[] c3 = new double[] { - vA.getY(), vB.getY(), 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 double m11 = minor(c2, c3); - final double m12 = minor(c1, c3); - final double m13 = minor(c1, c2); - final Vector2D center = new Vector2D(0.5 * m12 / m11, -0.5 * m13 / m11); - return new EnclosingBall(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]))); } } Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser3DTest.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser3DTest.java?rev=1563712&r1=1563711&r2=1563712&view=diff ============================================================================== --- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser3DTest.java (original) +++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser3DTest.java Sun Feb 2 20:52:49 2014 @@ -28,7 +28,6 @@ import org.apache.commons.math3.random.R 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(); Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGeneratorTest.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGeneratorTest.java?rev=1563712&r1=1563711&r2=1563712&view=diff ============================================================================== --- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGeneratorTest.java (original) +++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/threed/SphereGeneratorTest.java Sun Feb 2 20:52:49 2014 @@ -24,6 +24,7 @@ import org.apache.commons.math3.geometry 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)); + } + + } + }