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    *      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 }