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.enclosing;
23
24 import java.util.ArrayList;
25 import java.util.List;
26
27 import org.hipparchus.exception.MathRuntimeException;
28 import org.hipparchus.geometry.Point;
29 import org.hipparchus.geometry.Space;
30
31 /** Class implementing Emo Welzl algorithm to find the smallest enclosing ball in linear time.
32 * <p>
33 * The class implements the algorithm described in paper <a
34 * href="http://www.inf.ethz.ch/personal/emo/PublFiles/SmallEnclDisk_LNCS555_91.pdf">Smallest
35 * Enclosing Disks (Balls and Ellipsoids)</a> by Emo Welzl, Lecture Notes in Computer Science
36 * 555 (1991) 359-370. The pivoting improvement published in the paper <a
37 * href="http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf">Fast and
38 * Robust Smallest Enclosing Balls</a>, by Bernd Gärtner and further modified in
39 * paper <a href="http://www.idt.mdh.se/kurser/ct3340/ht12/MINICONFERENCE/FinalPapers/ircse12_submission_30.pdf">
40 * Efficient Computation of Smallest Enclosing Balls in Three Dimensions</a> by Linus Källberg
41 * to avoid performing local copies of data have been included.
42 * </p>
43 * @param <S> Space type.
44 * @param <P> Point type.
45 */
46 public class WelzlEncloser<S extends Space, P extends Point<S>> implements Encloser<S, P> {
47
48 /** Tolerance below which points are consider to be identical. */
49 private final double tolerance;
50
51 /** Generator for balls on support. */
52 private final SupportBallGenerator<S, P> generator;
53
54 /** Simple constructor.
55 * @param tolerance below which points are consider to be identical
56 * @param generator generator for balls on support
57 */
58 public WelzlEncloser(final double tolerance, final SupportBallGenerator<S, P> generator) {
59 this.tolerance = tolerance;
60 this.generator = generator;
61 }
62
63 /** {@inheritDoc} */
64 @Override
65 public EnclosingBall<S, P> enclose(final Iterable<P> points) {
66
67 if (points == null || !points.iterator().hasNext()) {
68 // return an empty ball
69 return generator.ballOnSupport(new ArrayList<P>());
70 }
71
72 // Emo Welzl algorithm with Bernd Gärtner and Linus Källberg improvements
73 return pivotingBall(points);
74
75 }
76
77 /** Compute enclosing ball using Gärtner's pivoting heuristic.
78 * @param points points to be enclosed
79 * @return enclosing ball
80 */
81 private EnclosingBall<S, P> pivotingBall(final Iterable<P> points) {
82
83 final P first = points.iterator().next();
84 final List<P> extreme = new ArrayList<>(first.getSpace().getDimension() + 1);
85 final List<P> support = new ArrayList<>(first.getSpace().getDimension() + 1);
86
87 // start with only first point selected as a candidate support
88 extreme.add(first);
89 EnclosingBall<S, P> ball = moveToFrontBall(extreme, extreme.size(), support);
90
91 while (true) {
92
93 // select the point farthest to current ball
94 final P farthest = selectFarthest(points, ball);
95
96 if (ball.contains(farthest, tolerance)) {
97 // we have found a ball containing all points
98 return ball;
99 }
100
101 // recurse search, restricted to the small subset containing support and farthest point
102 support.clear();
103 support.add(farthest);
104 EnclosingBall<S, P> savedBall = ball;
105 ball = moveToFrontBall(extreme, extreme.size(), support);
106 if (ball.getRadius() < savedBall.getRadius()) {
107 // this should never happen
108 throw MathRuntimeException.createInternalError();
109 }
110
111 // it was an interesting point, move it to the front
112 // according to Gärtner's heuristic
113 extreme.add(0, farthest);
114
115 // prune the least interesting points
116 extreme.subList(ball.getSupportSize(), extreme.size()).clear();
117
118
119 }
120 }
121
122 /** Compute enclosing ball using Welzl's move to front heuristic.
123 * @param extreme subset of extreme points
124 * @param nbExtreme number of extreme points to consider
125 * @param support points that must belong to the ball support
126 * @return enclosing ball, for the extreme subset only
127 */
128 private EnclosingBall<S, P> moveToFrontBall(final List<P> extreme, final int nbExtreme,
129 final List<P> support) {
130
131 // create a new ball on the prescribed support
132 EnclosingBall<S, P> ball = generator.ballOnSupport(support);
133
134 if (ball.getSupportSize() <= ball.getCenter().getSpace().getDimension()) {
135
136 for (int i = 0; i < nbExtreme; ++i) {
137 final P pi = extreme.get(i);
138 if (!ball.contains(pi, tolerance)) {
139
140 // we have found an outside point,
141 // enlarge the ball by adding it to the support
142 support.add(pi);
143 ball = moveToFrontBall(extreme, i, support);
144 support.remove(support.size() - 1);
145
146 // it was an interesting point, move it to the front
147 // according to Welzl's heuristic
148 for (int j = i; j > 0; --j) {
149 extreme.set(j, extreme.get(j - 1));
150 }
151 extreme.set(0, pi);
152
153 }
154 }
155
156 }
157
158 return ball;
159
160 }
161
162 /** Select the point farthest to the current ball.
163 * @param points points to be enclosed
164 * @param ball current ball
165 * @return farthest point
166 */
167 public P selectFarthest(final Iterable<P> points, final EnclosingBall<S, P> ball) {
168
169 final P center = ball.getCenter();
170 P farthest = null;
171 double dMax = -1.0;
172
173 for (final P point : points) {
174 final double d = point.distance(center);
175 if (d > dMax) {
176 farthest = point;
177 dMax = d;
178 }
179 }
180
181 return farthest;
182
183 }
184
185 }