commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
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 GMT
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<S extends Spa
 
             // 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<S extends Spa
             EnclosingBall<S, P> 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<Euclidean3D, Vector3D>(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<Euclidean3D, Vector3D>(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 4<sup>th</sup> column is known to be
filled with 1.0.
-     * <p>
-     * 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.
-     * </p>
      * @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<Euclidean2D, Vector2D>(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<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.
-     * <p>
-     * 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.
-     * </p>
      * @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<Vector3D> 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<Vector3D> points = new ArrayList<Vector3D>();
             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<Vector3D> 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<Euclidean3D, Vector3D> 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));
+        }
+
+    }
+
 }



Mime
View raw message