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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.geometry.euclidean.threed;
23
24 import java.util.Arrays;
25 import java.util.List;
26
27 import org.hipparchus.fraction.BigFraction;
28 import org.hipparchus.geometry.enclosing.EnclosingBall;
29 import org.hipparchus.geometry.enclosing.SupportBallGenerator;
30 import org.hipparchus.geometry.euclidean.twod.DiskGenerator;
31 import org.hipparchus.geometry.euclidean.twod.Euclidean2D;
32 import org.hipparchus.geometry.euclidean.twod.Vector2D;
33 import org.hipparchus.util.FastMath;
34
35 /** Class generating an enclosing ball from its support points.
36 */
37 public class SphereGenerator implements SupportBallGenerator<Euclidean3D, Vector3D> {
38
39 /** Empty constructor.
40 * <p>
41 * This constructor is not strictly necessary, but it prevents spurious
42 * javadoc warnings with JDK 18 and later.
43 * </p>
44 * @since 3.0
45 */
46 public SphereGenerator() { // NOPMD - unnecessary constructor added intentionally to make javadoc happy
47 // nothing to do
48 }
49
50 /** {@inheritDoc} */
51 @Override
52 public EnclosingBall<Euclidean3D, Vector3D> ballOnSupport(final List<Vector3D> support) {
53
54 if (support.isEmpty()) {
55 return new EnclosingBall<>(Vector3D.ZERO, Double.NEGATIVE_INFINITY);
56 } else {
57 final Vector3D vA = support.get(0);
58 if (support.size() < 2) {
59 return new EnclosingBall<>(vA, 0, vA);
60 } else {
61 final Vector3D vB = support.get(1);
62 if (support.size() < 3) {
63
64 final Vector3D center = new Vector3D(0.5, vA, 0.5, vB);
65
66 // we could have computed r directly from the vA and vB
67 // (it was done this way up to Hipparchus 1.0), but as center
68 // is approximated in the computation above, it is better to
69 // take the final value of center and compute r from the distances
70 // to center of all support points, using a max to ensure all support
71 // points belong to the ball
72 // see <https://github.com/Hipparchus-Math/hipparchus/issues/20>
73 final double r = FastMath.max(Vector3D.distance(vA, center),
74 Vector3D.distance(vB, center));
75 return new EnclosingBall<>(center, r, vA, vB);
76
77 } else {
78 final Vector3D vC = support.get(2);
79 if (support.size() < 4) {
80
81 // delegate to 2D disk generator
82 final Plane p = new Plane(vA, vB, vC,
83 1.0e-10 * (vA.getNorm1() + vB.getNorm1() + vC.getNorm1()));
84 final EnclosingBall<Euclidean2D, Vector2D> disk =
85 new DiskGenerator().ballOnSupport(Arrays.asList(p.toSubSpace(vA),
86 p.toSubSpace(vB),
87 p.toSubSpace(vC)));
88
89 // convert back to 3D
90 final Vector3D center = p.toSpace(disk.getCenter());
91
92 // we could have computed r directly from the vA and vB
93 // (it was done this way up to Hipparchus 1.0), but as center
94 // is approximated in the computation above, it is better to
95 // take the final value of center and compute r from the distances
96 // to center of all support points, using a max to ensure all support
97 // points belong to the ball
98 // see <https://github.com/Hipparchus-Math/hipparchus/issues/20>
99 final double r = FastMath.max(Vector3D.distance(vA, center),
100 FastMath.max(Vector3D.distance(vB, center),
101 Vector3D.distance(vC, center)));
102 return new EnclosingBall<>(center, r, vA, vB, vC);
103
104 } else {
105 final Vector3D vD = support.get(3);
106 // a sphere is 3D can be defined as:
107 // (1) (x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = r^2
108 // which can be written:
109 // (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
110 // or simply:
111 // (3) (x^2 + y^2 + z^2) + a x + b y + c z + d = 0
112 // with sphere center coordinates -a/2, -b/2, -c/2
113 // If the sphere exists, a b, c and d are a non zero solution to
114 // [ (x^2 + y^2 + z^2) x y z 1 ] [ 1 ] [ 0 ]
115 // [ (xA^2 + yA^2 + zA^2) xA yA zA 1 ] [ a ] [ 0 ]
116 // [ (xB^2 + yB^2 + zB^2) xB yB zB 1 ] * [ b ] = [ 0 ]
117 // [ (xC^2 + yC^2 + zC^2) xC yC zC 1 ] [ c ] [ 0 ]
118 // [ (xD^2 + yD^2 + zD^2) xD yD zD 1 ] [ d ] [ 0 ]
119 // So the determinant of the matrix is zero. Computing this determinant
120 // by expanding it using the minors m_ij of first row leads to
121 // (4) m_11 (x^2 + y^2 + z^2) - m_12 x + m_13 y - m_14 z + m_15 = 0
122 // So by identifying equations (2) and (4) we get the coordinates
123 // of center as:
124 // x_0 = +m_12 / (2 m_11)
125 // y_0 = -m_13 / (2 m_11)
126 // z_0 = +m_14 / (2 m_11)
127 // Note that the minors m_11, m_12, m_13 and m_14 all have the last column
128 // filled with 1.0, hence simplifying the computation
129 final BigFraction[] c2 = {
130 new BigFraction(vA.getX()), new BigFraction(vB.getX()),
131 new BigFraction(vC.getX()), new BigFraction(vD.getX())
132 };
133 final BigFraction[] c3 = {
134 new BigFraction(vA.getY()), new BigFraction(vB.getY()),
135 new BigFraction(vC.getY()), new BigFraction(vD.getY())
136 };
137 final BigFraction[] c4 = {
138 new BigFraction(vA.getZ()), new BigFraction(vB.getZ()),
139 new BigFraction(vC.getZ()), new BigFraction(vD.getZ())
140 };
141 final BigFraction[] c1 = {
142 c2[0].multiply(c2[0]).add(c3[0].multiply(c3[0])).add(c4[0].multiply(c4[0])),
143 c2[1].multiply(c2[1]).add(c3[1].multiply(c3[1])).add(c4[1].multiply(c4[1])),
144 c2[2].multiply(c2[2]).add(c3[2].multiply(c3[2])).add(c4[2].multiply(c4[2])),
145 c2[3].multiply(c2[3]).add(c3[3].multiply(c3[3])).add(c4[3].multiply(c4[3]))
146 };
147 final BigFraction twoM11 = minor(c2, c3, c4).multiply(2);
148 final BigFraction m12 = minor(c1, c3, c4);
149 final BigFraction m13 = minor(c1, c2, c4);
150 final BigFraction m14 = minor(c1, c2, c3);
151 final Vector3D center = new Vector3D( m12.divide(twoM11).doubleValue(),
152 -m13.divide(twoM11).doubleValue(),
153 m14.divide(twoM11).doubleValue());
154
155 // we could have computed r directly from the minors above
156 // (it was done this way up to Hipparchus 1.0), but as center
157 // is approximated in the computation above, it is better to
158 // take the final value of center and compute r from the distances
159 // to center of all support points, using a max to ensure all support
160 // points belong to the ball
161 // see <https://github.com/Hipparchus-Math/hipparchus/issues/20>
162 final double r = FastMath.max(Vector3D.distance(vA, center),
163 FastMath.max(Vector3D.distance(vB, center),
164 FastMath.max(Vector3D.distance(vC, center),
165 Vector3D.distance(vD, center))));
166 return new EnclosingBall<>(center, r, vA, vB, vC, vD);
167
168 }
169 }
170 }
171 }
172 }
173
174 /** Compute a dimension 4 minor, when 4<sup>th</sup> column is known to be filled with 1.0.
175 * @param c1 first column
176 * @param c2 second column
177 * @param c3 third column
178 * @return value of the minor computed has an exact fraction
179 */
180 private BigFraction minor(final BigFraction[] c1, final BigFraction[] c2, final BigFraction[] c3) {
181 return c2[0].multiply(c3[1]).multiply(c1[2].subtract(c1[3])).
182 add(c2[0].multiply(c3[2]).multiply(c1[3].subtract(c1[1]))).
183 add(c2[0].multiply(c3[3]).multiply(c1[1].subtract(c1[2]))).
184 add(c2[1].multiply(c3[0]).multiply(c1[3].subtract(c1[2]))).
185 add(c2[1].multiply(c3[2]).multiply(c1[0].subtract(c1[3]))).
186 add(c2[1].multiply(c3[3]).multiply(c1[2].subtract(c1[0]))).
187 add(c2[2].multiply(c3[0]).multiply(c1[1].subtract(c1[3]))).
188 add(c2[2].multiply(c3[1]).multiply(c1[3].subtract(c1[0]))).
189 add(c2[2].multiply(c3[3]).multiply(c1[0].subtract(c1[1]))).
190 add(c2[3].multiply(c3[0]).multiply(c1[2].subtract(c1[1]))).
191 add(c2[3].multiply(c3[1]).multiply(c1[0].subtract(c1[2]))).
192 add(c2[3].multiply(c3[2]).multiply(c1[1].subtract(c1[0])));
193 }
194
195 }