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}