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
This commit is contained in:
Luc Maisonobe
parent
7897aa6a83
commit
19c1c3bb9b
|
|
@ -89,6 +89,7 @@ public class WelzlEncloser<S extends Space, P extends Point<S>> 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<S extends Space, P extends Point<S>> implements Enclo
|
|||
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();
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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<Euclidean3D, Vector
|
|||
// 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 SupportBallGenerator<Euclidean3D, Vector
|
|||
}
|
||||
|
||||
/** 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])));
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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<Euclidean2D, Vector2D
|
|||
// 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])));
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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<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();
|
||||
|
|
|
|||
|
|
@ -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<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));
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue