001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.math3.geometry.euclidean.threed; 018 019import java.util.Arrays; 020import java.util.List; 021 022import org.apache.commons.math3.fraction.BigFraction; 023import org.apache.commons.math3.geometry.enclosing.EnclosingBall; 024import org.apache.commons.math3.geometry.enclosing.SupportBallGenerator; 025import org.apache.commons.math3.geometry.euclidean.twod.DiskGenerator; 026import org.apache.commons.math3.geometry.euclidean.twod.Euclidean2D; 027import org.apache.commons.math3.geometry.euclidean.twod.Vector2D; 028import org.apache.commons.math3.util.FastMath; 029 030/** Class generating an enclosing ball from its support points. 031 * @since 3.3 032 */ 033public class SphereGenerator implements SupportBallGenerator<Euclidean3D, Vector3D> { 034 035 /** {@inheritDoc} */ 036 public EnclosingBall<Euclidean3D, Vector3D> ballOnSupport(final List<Vector3D> support) { 037 038 if (support.size() < 1) { 039 return new EnclosingBall<Euclidean3D, Vector3D>(Vector3D.ZERO, Double.NEGATIVE_INFINITY); 040 } else { 041 final Vector3D vA = support.get(0); 042 if (support.size() < 2) { 043 return new EnclosingBall<Euclidean3D, Vector3D>(vA, 0, vA); 044 } else { 045 final Vector3D vB = support.get(1); 046 if (support.size() < 3) { 047 return new EnclosingBall<Euclidean3D, Vector3D>(new Vector3D(0.5, vA, 0.5, vB), 048 0.5 * vA.distance(vB), 049 vA, vB); 050 } else { 051 final Vector3D vC = support.get(2); 052 if (support.size() < 4) { 053 054 // delegate to 2D disk generator 055 final Plane p = new Plane(vA, vB, vC, 056 1.0e-10 * (vA.getNorm1() + vB.getNorm1() + vC.getNorm1())); 057 final EnclosingBall<Euclidean2D, Vector2D> disk = 058 new DiskGenerator().ballOnSupport(Arrays.asList(p.toSubSpace(vA), 059 p.toSubSpace(vB), 060 p.toSubSpace(vC))); 061 062 // convert back to 3D 063 return new EnclosingBall<Euclidean3D, Vector3D>(p.toSpace(disk.getCenter()), 064 disk.getRadius(), vA, vB, vC); 065 066 } else { 067 final Vector3D vD = support.get(3); 068 // a sphere is 3D can be defined as: 069 // (1) (x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = r^2 070 // which can be written: 071 // (2) (x^2 + y^2 + z^2) - 2 x_0 x - 2 y_0 y - 2 z_0 z + (x_0^2 + y_0^2 + z_0^2 - r^2) = 0 072 // or simply: 073 // (3) (x^2 + y^2 + z^2) + a x + b y + c z + d = 0 074 // with sphere center coordinates -a/2, -b/2, -c/2 075 // If the sphere exists, a b, c and d are a non zero solution to 076 // [ (x^2 + y^2 + z^2) x y z 1 ] [ 1 ] [ 0 ] 077 // [ (xA^2 + yA^2 + zA^2) xA yA zA 1 ] [ a ] [ 0 ] 078 // [ (xB^2 + yB^2 + zB^2) xB yB zB 1 ] * [ b ] = [ 0 ] 079 // [ (xC^2 + yC^2 + zC^2) xC yC zC 1 ] [ c ] [ 0 ] 080 // [ (xD^2 + yD^2 + zD^2) xD yD zD 1 ] [ d ] [ 0 ] 081 // So the determinant of the matrix is zero. Computing this determinant 082 // by expanding it using the minors m_ij of first row leads to 083 // (4) m_11 (x^2 + y^2 + z^2) - m_12 x + m_13 y - m_14 z + m_15 = 0 084 // So by identifying equations (2) and (4) we get the coordinates 085 // of center as: 086 // x_0 = +m_12 / (2 m_11) 087 // y_0 = -m_13 / (2 m_11) 088 // z_0 = +m_14 / (2 m_11) 089 // Note that the minors m_11, m_12, m_13 and m_14 all have the last column 090 // filled with 1.0, hence simplifying the computation 091 final BigFraction[] c2 = new BigFraction[] { 092 new BigFraction(vA.getX()), new BigFraction(vB.getX()), 093 new BigFraction(vC.getX()), new BigFraction(vD.getX()) 094 }; 095 final BigFraction[] c3 = new BigFraction[] { 096 new BigFraction(vA.getY()), new BigFraction(vB.getY()), 097 new BigFraction(vC.getY()), new BigFraction(vD.getY()) 098 }; 099 final BigFraction[] c4 = new BigFraction[] { 100 new BigFraction(vA.getZ()), new BigFraction(vB.getZ()), 101 new BigFraction(vC.getZ()), new BigFraction(vD.getZ()) 102 }; 103 final BigFraction[] c1 = new BigFraction[] { 104 c2[0].multiply(c2[0]).add(c3[0].multiply(c3[0])).add(c4[0].multiply(c4[0])), 105 c2[1].multiply(c2[1]).add(c3[1].multiply(c3[1])).add(c4[1].multiply(c4[1])), 106 c2[2].multiply(c2[2]).add(c3[2].multiply(c3[2])).add(c4[2].multiply(c4[2])), 107 c2[3].multiply(c2[3]).add(c3[3].multiply(c3[3])).add(c4[3].multiply(c4[3])) 108 }; 109 final BigFraction twoM11 = minor(c2, c3, c4).multiply(2); 110 final BigFraction m12 = minor(c1, c3, c4); 111 final BigFraction m13 = minor(c1, c2, c4); 112 final BigFraction m14 = minor(c1, c2, c3); 113 final BigFraction centerX = m12.divide(twoM11); 114 final BigFraction centerY = m13.divide(twoM11).negate(); 115 final BigFraction centerZ = m14.divide(twoM11); 116 final BigFraction dx = c2[0].subtract(centerX); 117 final BigFraction dy = c3[0].subtract(centerY); 118 final BigFraction dz = c4[0].subtract(centerZ); 119 final BigFraction r2 = dx.multiply(dx).add(dy.multiply(dy)).add(dz.multiply(dz)); 120 return new EnclosingBall<Euclidean3D, Vector3D>(new Vector3D(centerX.doubleValue(), 121 centerY.doubleValue(), 122 centerZ.doubleValue()), 123 FastMath.sqrt(r2.doubleValue()), 124 vA, vB, vC, vD); 125 } 126 } 127 } 128 } 129 } 130 131 /** Compute a dimension 4 minor, when 4<sup>th</sup> column is known to be filled with 1.0. 132 * @param c1 first column 133 * @param c2 second column 134 * @param c3 third column 135 * @return value of the minor computed has an exact fraction 136 */ 137 private BigFraction minor(final BigFraction[] c1, final BigFraction[] c2, final BigFraction[] c3) { 138 return c2[0].multiply(c3[1]).multiply(c1[2].subtract(c1[3])). 139 add(c2[0].multiply(c3[2]).multiply(c1[3].subtract(c1[1]))). 140 add(c2[0].multiply(c3[3]).multiply(c1[1].subtract(c1[2]))). 141 add(c2[1].multiply(c3[0]).multiply(c1[3].subtract(c1[2]))). 142 add(c2[1].multiply(c3[2]).multiply(c1[0].subtract(c1[3]))). 143 add(c2[1].multiply(c3[3]).multiply(c1[2].subtract(c1[0]))). 144 add(c2[2].multiply(c3[0]).multiply(c1[1].subtract(c1[3]))). 145 add(c2[2].multiply(c3[1]).multiply(c1[3].subtract(c1[0]))). 146 add(c2[2].multiply(c3[3]).multiply(c1[0].subtract(c1[1]))). 147 add(c2[3].multiply(c3[0]).multiply(c1[2].subtract(c1[1]))). 148 add(c2[3].multiply(c3[1]).multiply(c1[0].subtract(c1[2]))). 149 add(c2[3].multiply(c3[2]).multiply(c1[1].subtract(c1[0]))); 150 } 151 152}