View Javadoc

    1   /*
    2    * Licensed to the Apache Software Foundation (ASF) under one or more
    3    * contributor license agreements.  See the NOTICE file distributed with
    4    * this work for additional information regarding copyright ownership.
    5    * The ASF licenses this file to You under the Apache License, Version 2.0
    6    * (the "License"); you may not use this file except in compliance with
    7    * the License.  You may obtain a copy of the License at
    8    *
    9    *      http://www.apache.org/licenses/LICENSE-2.0
    10   *
    11   * Unless required by applicable law or agreed to in writing, software
    12   * distributed under the License is distributed on an "AS IS" BASIS,
    13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    14   * See the License for the specific language governing permissions and
    15   * limitations under the License.
    16   */
    17  package org.apache.commons.geometry.enclosing.euclidean.threed;
    18
    19  import java.util.Arrays;
    20  import java.util.Collections;
    21  import java.util.List;
    22
    23  import org.apache.commons.geometry.core.precision.DoublePrecisionContext;
    24  import org.apache.commons.geometry.enclosing.EnclosingBall;
    25  import org.apache.commons.geometry.enclosing.SupportBallGenerator;
    26  import org.apache.commons.geometry.enclosing.euclidean.twod.DiskGenerator;
    27  import org.apache.commons.geometry.euclidean.threed.EmbeddingPlane;
    28  import org.apache.commons.geometry.euclidean.threed.Planes;
    29  import org.apache.commons.geometry.euclidean.threed.Vector3D;
    30  import org.apache.commons.geometry.euclidean.twod.Vector2D;
    31  import org.apache.commons.numbers.fraction.BigFraction;
    32
    33  /** Class generating a sphere from its support points.
    34   */
    35  public class SphereGenerator implements SupportBallGenerator<Vector3D> {
    36
    37      /** Precision context used to compare floating point numbers. */
    38      private final DoublePrecisionContext precision;
    39
    40      /** Construct a new instance with the given precision context.
    41       * @param precision precision context used to compare floating point numbers
    42       */
    43      public SphereGenerator(final DoublePrecisionContext precision) {
    44          this.precision = precision;
    45      }
    46
    47      /** {@inheritDoc} */
    48      @Override
    49      public EnclosingBall<Vector3D> ballOnSupport(final List<Vector3D> support) {
    50          if (support.isEmpty()) {
    51              return new EnclosingBall<>(Vector3D.ZERO, Double.NEGATIVE_INFINITY, Collections.emptyList());
    52          }
    53          final Vector3D vA = support.get(0);
    54          if (support.size() < 2) {
    55              return new EnclosingBall<>(vA, 0, Collections.singletonList(vA));
    56          }
    57          final Vector3D vB = support.get(1);
    58          if (support.size() < 3) {
    59              return new EnclosingBall<>(Vector3D.linearCombination(0.5, vA, 0.5, vB),
    60                                         0.5 * vA.distance(vB),
    61                                         Arrays.asList(vA, vB));
    62          }
    63          final Vector3D vC = support.get(2);
    64          if (support.size() < 4) {
    65              final EmbeddingPlane p = Planes.fromPoints(vA, vB, vC, precision).getEmbedding();
    66              final EnclosingBall<Vector2D> disk =
    67                      new DiskGenerator().ballOnSupport(Arrays.asList(p.toSubspace(vA),
    68                                                                      p.toSubspace(vB),
    69                                                                      p.toSubspace(vC)));
    70
    71              // convert back to 3D
    72              return new EnclosingBall<>(p.toSpace(disk.getCenter()),
    73                                                              disk.getRadius(),
    74                                                              Arrays.asList(vA, vB, vC));
    75
    76          }
    77          final Vector3D vD = support.get(3);
    78          // a sphere is 3D can be defined as:
    79          // (1)   (x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = r^2
    80          // which can be written:
    81          // (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
    82          // or simply:
    83          // (3)   (x^2 + y^2 + z^2) + a x + b y + c z + d = 0
    84          // with sphere center coordinates -a/2, -b/2, -c/2
    85          // If the sphere exists, a b, c and d are a non zero solution to
    86          // [ (x^2  + y^2  + z^2)    x    y   z    1 ]   [ 1 ]   [ 0 ]
    87          // [ (xA^2 + yA^2 + zA^2)   xA   yA  zA   1 ]   [ a ]   [ 0 ]
    88          // [ (xB^2 + yB^2 + zB^2)   xB   yB  zB   1 ] * [ b ] = [ 0 ]
    89          // [ (xC^2 + yC^2 + zC^2)   xC   yC  zC   1 ]   [ c ]   [ 0 ]
    90          // [ (xD^2 + yD^2 + zD^2)   xD   yD  zD   1 ]   [ d ]   [ 0 ]
    91          // So the determinant of the matrix is zero. Computing this determinant
    92          // by expanding it using the minors m_ij of first row leads to
    93          // (4)   m_11 (x^2 + y^2 + z^2) - m_12 x + m_13 y - m_14 z + m_15 = 0
    94          // So by identifying equations (2) and (4) we get the coordinates
    95          // of center as:
    96          //      x_0 = +m_12 / (2 m_11)
    97          //      y_0 = -m_13 / (2 m_11)
    98          //      z_0 = +m_14 / (2 m_11)
    99          // Note that the minors m_11, m_12, m_13 and m_14 all have the last column
    100         // filled with 1.0, hence simplifying the computation
    101         final BigFraction[] c2 = {
    102             BigFraction.from(vA.getX()), BigFraction.from(vB.getX()),
    103             BigFraction.from(vC.getX()), BigFraction.from(vD.getX())
    104         };
    105         final BigFraction[] c3 = {
    106             BigFraction.from(vA.getY()), BigFraction.from(vB.getY()),
    107             BigFraction.from(vC.getY()), BigFraction.from(vD.getY())
    108         };
    109         final BigFraction[] c4 = {
    110             BigFraction.from(vA.getZ()), BigFraction.from(vB.getZ()),
    111             BigFraction.from(vC.getZ()), BigFraction.from(vD.getZ())
    112         };
    113         final BigFraction[] c1 = {
    114             c2[0].multiply(c2[0]).add(c3[0].multiply(c3[0])).add(c4[0].multiply(c4[0])),
    115             c2[1].multiply(c2[1]).add(c3[1].multiply(c3[1])).add(c4[1].multiply(c4[1])),
    116             c2[2].multiply(c2[2]).add(c3[2].multiply(c3[2])).add(c4[2].multiply(c4[2])),
    117             c2[3].multiply(c2[3]).add(c3[3].multiply(c3[3])).add(c4[3].multiply(c4[3]))
    118         };
    119         final BigFraction twoM11 = minor(c2, c3, c4).multiply(2);
    120         final BigFraction m12 = minor(c1, c3, c4);
    121         final BigFraction m13 = minor(c1, c2, c4);
    122         final BigFraction m14 = minor(c1, c2, c3);
    123         final BigFraction centerX = m12.divide(twoM11);
    124         final BigFraction centerY = m13.divide(twoM11).negate();
    125         final BigFraction centerZ = m14.divide(twoM11);
    126         final BigFraction dx = c2[0].subtract(centerX);
    127         final BigFraction dy = c3[0].subtract(centerY);
    128         final BigFraction dz = c4[0].subtract(centerZ);
    129         final BigFraction r2 = dx.multiply(dx).add(dy.multiply(dy)).add(dz.multiply(dz));
    130         return new EnclosingBall<>(Vector3D.of(centerX.doubleValue(),
    131                                                centerY.doubleValue(),
    132                                                centerZ.doubleValue()),
    133                                    Math.sqrt(r2.doubleValue()),
    134                                    Arrays.asList(vA, vB, vC, vD));
    135     }
    136
    137     /** Compute a dimension 4 minor, when 4<sup>th</sup> column is known to be filled with 1.0.
    138      * @param c1 first column
    139      * @param c2 second column
    140      * @param c3 third column
    141      * @return value of the minor computed has an exact fraction
    142      */
    143     private BigFraction minor(final BigFraction[] c1, final BigFraction[] c2, final BigFraction[] c3) {
    144         return c2[0].multiply(c3[1]).multiply(c1[2].subtract(c1[3])).
    145             add(c2[0].multiply(c3[2]).multiply(c1[3].subtract(c1[1]))).
    146             add(c2[0].multiply(c3[3]).multiply(c1[1].subtract(c1[2]))).
    147             add(c2[1].multiply(c3[0]).multiply(c1[3].subtract(c1[2]))).
    148             add(c2[1].multiply(c3[2]).multiply(c1[0].subtract(c1[3]))).
    149             add(c2[1].multiply(c3[3]).multiply(c1[2].subtract(c1[0]))).
    150             add(c2[2].multiply(c3[0]).multiply(c1[1].subtract(c1[3]))).
    151             add(c2[2].multiply(c3[1]).multiply(c1[3].subtract(c1[0]))).
    152             add(c2[2].multiply(c3[3]).multiply(c1[0].subtract(c1[1]))).
    153             add(c2[3].multiply(c3[0]).multiply(c1[2].subtract(c1[1]))).
    154             add(c2[3].multiply(c3[1]).multiply(c1[0].subtract(c1[2]))).
    155             add(c2[3].multiply(c3[2]).multiply(c1[1].subtract(c1[0])));
    156     }
    157 }