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