summaryrefslogtreecommitdiff
path: root/src/main/java/org/apache/commons/math3/geometry/enclosing/WelzlEncloser.java
blob: 987e7d9f28057fe9441481532fd4b04e46f85b18 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.commons.math3.geometry.enclosing;

import java.util.ArrayList;
import java.util.List;

import org.apache.commons.math3.exception.MathInternalError;
import org.apache.commons.math3.geometry.Point;
import org.apache.commons.math3.geometry.Space;

/** Class implementing Emo Welzl algorithm to find the smallest enclosing ball in linear time.
 * <p>
 * The class implements the algorithm described in paper <a
 * href="http://www.inf.ethz.ch/personal/emo/PublFiles/SmallEnclDisk_LNCS555_91.pdf">Smallest
 * Enclosing Disks (Balls and Ellipsoids)</a> by Emo Welzl, Lecture Notes in Computer Science
 * 555 (1991) 359-370. The pivoting improvement published in the paper <a
 * href="http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf">Fast and
 * Robust Smallest Enclosing Balls</a>, by Bernd Gärtner and further modified in
 * paper <a
 * href=http://www.idt.mdh.se/kurser/ct3340/ht12/MINICONFERENCE/FinalPapers/ircse12_submission_30.pdf">
 * Efficient Computation of Smallest Enclosing Balls in Three Dimensions</a> by Linus Källberg
 * to avoid performing local copies of data have been included.
 * </p>
 * @param <S> Space type.
 * @param <P> Point type.
 * @since 3.3
 */
public class WelzlEncloser<S extends Space, P extends Point<S>> implements Encloser<S, P> {

    /** Tolerance below which points are consider to be identical. */
    private final double tolerance;

    /** Generator for balls on support. */
    private final SupportBallGenerator<S, P> generator;

    /** Simple constructor.
     * @param tolerance below which points are consider to be identical
     * @param generator generator for balls on support
     */
    public WelzlEncloser(final double tolerance, final SupportBallGenerator<S, P> generator) {
        this.tolerance = tolerance;
        this.generator = generator;
    }

    /** {@inheritDoc} */
    public EnclosingBall<S, P> enclose(final Iterable<P> points) {

        if (points == null || !points.iterator().hasNext()) {
            // return an empty ball
            return generator.ballOnSupport(new ArrayList<P>());
        }

        // Emo Welzl algorithm with Bernd Gärtner and Linus Källberg improvements
        return pivotingBall(points);

    }

    /** Compute enclosing ball using Gärtner's pivoting heuristic.
     * @param points points to be enclosed
     * @return enclosing ball
     */
    private EnclosingBall<S, P> pivotingBall(final Iterable<P> points) {

        final P first = points.iterator().next();
        final List<P> extreme = new ArrayList<P>(first.getSpace().getDimension() + 1);
        final List<P> support = new ArrayList<P>(first.getSpace().getDimension() + 1);

        // start with only first point selected as a candidate support
        extreme.add(first);
        EnclosingBall<S, P> ball = moveToFrontBall(extreme, extreme.size(), support);

        while (true) {

            // select the point farthest to current ball
            final P farthest = selectFarthest(points, ball);

            if (ball.contains(farthest, tolerance)) {
                // we have found a ball containing all points
                return ball;
            }

            // recurse search, restricted to the small subset containing support and farthest point
            support.clear();
            support.add(farthest);
            EnclosingBall<S, P> savedBall = ball;
            ball = moveToFrontBall(extreme, extreme.size(), support);
            if (ball.getRadius() < savedBall.getRadius()) {
                // this should never happen
                throw new MathInternalError();
            }

            // it was an interesting point, move it to the front
            // according to Gärtner's heuristic
            extreme.add(0, farthest);

            // prune the least interesting points
            extreme.subList(ball.getSupportSize(), extreme.size()).clear();


        }
    }

    /** Compute enclosing ball using Welzl's move to front heuristic.
     * @param extreme subset of extreme points
     * @param nbExtreme number of extreme points to consider
     * @param support points that must belong to the ball support
     * @return enclosing ball, for the extreme subset only
     */
    private EnclosingBall<S, P> moveToFrontBall(final List<P> extreme, final int nbExtreme,
                                                final List<P> support) {

        // create a new ball on the prescribed support
        EnclosingBall<S, P> ball = generator.ballOnSupport(support);

        if (ball.getSupportSize() <= ball.getCenter().getSpace().getDimension()) {

            for (int i = 0; i < nbExtreme; ++i) {
                final P pi = extreme.get(i);
                if (!ball.contains(pi, tolerance)) {

                    // we have found an outside point,
                    // enlarge the ball by adding it to the support
                    support.add(pi);
                    ball = moveToFrontBall(extreme, i, support);
                    support.remove(support.size() - 1);

                    // it was an interesting point, move it to the front
                    // according to Welzl's heuristic
                    for (int j = i; j > 0; --j) {
                        extreme.set(j, extreme.get(j - 1));
                    }
                    extreme.set(0, pi);

                }
            }

        }

        return ball;

    }

    /** Select the point farthest to the current ball.
     * @param points points to be enclosed
     * @param ball current ball
     * @return farthest point
     */
    public P selectFarthest(final Iterable<P> points, final EnclosingBall<S, P> ball) {

        final P center = ball.getCenter();
        P farthest   = null;
        double dMax  = -1.0;

        for (final P point : points) {
            final double d = point.distance(center);
            if (d > dMax) {
                farthest = point;
                dMax     = d;
            }
        }

        return farthest;

    }

}