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.enclosing;
018
019import java.util.ArrayList;
020import java.util.List;
021
022import org.apache.commons.math3.exception.MathInternalError;
023import org.apache.commons.math3.geometry.Point;
024import org.apache.commons.math3.geometry.Space;
025
026/** Class implementing Emo Welzl algorithm to find the smallest enclosing ball in linear time.
027 * <p>
028 * The class implements the algorithm described in paper <a
029 * href="http://www.inf.ethz.ch/personal/emo/PublFiles/SmallEnclDisk_LNCS555_91.pdf">Smallest
030 * Enclosing Disks (Balls and Ellipsoids)</a> by Emo Welzl, Lecture Notes in Computer Science
031 * 555 (1991) 359-370. The pivoting improvement published in the paper <a
032 * href="http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf">Fast and
033 * Robust Smallest Enclosing Balls</a>, by Bernd Gärtner and further modified in
034 * paper <a
035 * href=http://www.idt.mdh.se/kurser/ct3340/ht12/MINICONFERENCE/FinalPapers/ircse12_submission_30.pdf">
036 * Efficient Computation of Smallest Enclosing Balls in Three Dimensions</a> by Linus Källberg
037 * to avoid performing local copies of data have been included.
038 * </p>
039 * @param <S> Space type.
040 * @param <P> Point type.
041 * @since 3.3
042 */
043public class WelzlEncloser<S extends Space, P extends Point<S>> implements Encloser<S, P> {
044
045    /** Tolerance below which points are consider to be identical. */
046    private final double tolerance;
047
048    /** Generator for balls on support. */
049    private final SupportBallGenerator<S, P> generator;
050
051    /** Simple constructor.
052     * @param tolerance below which points are consider to be identical
053     * @param generator generator for balls on support
054     */
055    public WelzlEncloser(final double tolerance, final SupportBallGenerator<S, P> generator) {
056        this.tolerance = tolerance;
057        this.generator = generator;
058    }
059
060    /** {@inheritDoc} */
061    public EnclosingBall<S, P> enclose(final Iterable<P> points) {
062
063        if (points == null || !points.iterator().hasNext()) {
064            // return an empty ball
065            return generator.ballOnSupport(new ArrayList<P>());
066        }
067
068        // Emo Welzl algorithm with Bernd Gärtner and Linus Källberg improvements
069        return pivotingBall(points);
070
071    }
072
073    /** Compute enclosing ball using Gärtner's pivoting heuristic.
074     * @param points points to be enclosed
075     * @return enclosing ball
076     */
077    private EnclosingBall<S, P> pivotingBall(final Iterable<P> points) {
078
079        final P first = points.iterator().next();
080        final List<P> extreme = new ArrayList<P>(first.getSpace().getDimension() + 1);
081        final List<P> support = new ArrayList<P>(first.getSpace().getDimension() + 1);
082
083        // start with only first point selected as a candidate support
084        extreme.add(first);
085        EnclosingBall<S, P> ball = moveToFrontBall(extreme, extreme.size(), support);
086
087        while (true) {
088
089            // select the point farthest to current ball
090            final P farthest = selectFarthest(points, ball);
091
092            if (ball.contains(farthest, tolerance)) {
093                // we have found a ball containing all points
094                return ball;
095            }
096
097            // recurse search, restricted to the small subset containing support and farthest point
098            support.clear();
099            support.add(farthest);
100            EnclosingBall<S, P> savedBall = ball;
101            ball = moveToFrontBall(extreme, extreme.size(), support);
102            if (ball.getRadius() < savedBall.getRadius()) {
103                // this should never happen
104                throw new MathInternalError();
105            }
106
107            // it was an interesting point, move it to the front
108            // according to Gärtner's heuristic
109            extreme.add(0, farthest);
110
111            // prune the least interesting points
112            extreme.subList(ball.getSupportSize(), extreme.size()).clear();
113
114
115        }
116    }
117
118    /** Compute enclosing ball using Welzl's move to front heuristic.
119     * @param extreme subset of extreme points
120     * @param nbExtreme number of extreme points to consider
121     * @param support points that must belong to the ball support
122     * @return enclosing ball, for the extreme subset only
123     */
124    private EnclosingBall<S, P> moveToFrontBall(final List<P> extreme, final int nbExtreme,
125                                                final List<P> support) {
126
127        // create a new ball on the prescribed support
128        EnclosingBall<S, P> ball = generator.ballOnSupport(support);
129
130        if (ball.getSupportSize() <= ball.getCenter().getSpace().getDimension()) {
131
132            for (int i = 0; i < nbExtreme; ++i) {
133                final P pi = extreme.get(i);
134                if (!ball.contains(pi, tolerance)) {
135
136                    // we have found an outside point,
137                    // enlarge the ball by adding it to the support
138                    support.add(pi);
139                    ball = moveToFrontBall(extreme, i, support);
140                    support.remove(support.size() - 1);
141
142                    // it was an interesting point, move it to the front
143                    // according to Welzl's heuristic
144                    for (int j = i; j > 0; --j) {
145                        extreme.set(j, extreme.get(j - 1));
146                    }
147                    extreme.set(0, pi);
148
149                }
150            }
151
152        }
153
154        return ball;
155
156    }
157
158    /** Select the point farthest to the current ball.
159     * @param points points to be enclosed
160     * @param ball current ball
161     * @return farthest point
162     */
163    public P selectFarthest(final Iterable<P> points, final EnclosingBall<S, P> ball) {
164
165        final P center = ball.getCenter();
166        P farthest   = null;
167        double dMax  = -1.0;
168
169        for (final P point : points) {
170            final double d = point.distance(center);
171            if (d > dMax) {
172                farthest = point;
173                dMax     = d;
174            }
175        }
176
177        return farthest;
178
179    }
180
181}