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.spherical.oned;
    18
    19  import java.util.ArrayList;
    20  import java.util.Collections;
    21  import java.util.Comparator;
    22  import java.util.List;
    23  import java.util.Objects;
    24
    25  import org.apache.commons.geometry.core.Transform;
    26  import org.apache.commons.geometry.core.partitioning.Hyperplane;
    27  import org.apache.commons.geometry.core.partitioning.HyperplaneLocation;
    28  import org.apache.commons.geometry.core.partitioning.HyperplaneSubset;
    29  import org.apache.commons.geometry.core.partitioning.Split;
    30  import org.apache.commons.geometry.core.partitioning.bsp.AbstractBSPTree;
    31  import org.apache.commons.geometry.core.partitioning.bsp.AbstractRegionBSPTree;
    32  import org.apache.commons.geometry.euclidean.twod.Vector2D;
    33  import org.apache.commons.numbers.angle.Angle;
    34  import org.apache.commons.numbers.core.Precision;
    35
    36  /** BSP tree representing regions in 1D spherical space.
    37   */
    38  public class RegionBSPTree1S extends AbstractRegionBSPTree<Point1S, RegionBSPTree1S.RegionNode1S> {
    39      /** Comparator used to sort BoundaryPairs by ascending azimuth.  */
    40      private static final Comparator<BoundaryPair> BOUNDARY_PAIR_COMPARATOR =
    41              Comparator.comparingDouble(BoundaryPair::getMinValue);
    42
    43      /** Create a new, empty instance.
    44       */
    45      public RegionBSPTree1S() {
    46          this(false);
    47      }
    48
    49      /** Create a new region. If {@code full} is true, then the region will
    50       * represent the entire circle. Otherwise, it will be empty.
    51       * @param full whether or not the region should contain the entire
    52       *      circle or be empty
    53       */
    54      public RegionBSPTree1S(final boolean full) {
    55          super(full);
    56      }
    57
    58      /** Return a deep copy of this instance.
    59       * @return a deep copy of this instance.
    60       * @see #copy(org.apache.commons.geometry.core.partitioning.bsp.BSPTree)
    61       */
    62      public RegionBSPTree1S copy() {
    63          final RegionBSPTree1S result = RegionBSPTree1S.empty();
    64          result.copy(this);
    65
    66          return result;
    67      }
    68
    69      /** Add an interval to this region. The resulting region will be the
    70       * union of the interval and the region represented by this instance.
    71       * @param interval the interval to add
    72       */
    73      public void add(final AngularInterval interval) {
    74          union(fromInterval(interval));
    75      }
    76
    77      /** {@inheritDoc} */
    78      @Override
    79      public Point1S project(final Point1S pt) {
    80          final BoundaryProjector1S projector = new BoundaryProjector1S(pt);
    81          accept(projector);
    82
    83          return projector.getProjected();
    84      }
    85
    86      /** {@inheritDoc}
    87       *
    88       * <p>Each interval of the region is transformed individually and the
    89       * results are unioned. If the size of any transformed interval is greater
    90       * than or equal to 2pi, then the region is set to the full space.</p>
    91       */
    92      @Override
    93      public void transform(final Transform<Point1S> transform) {
    94          if (!isFull() && !isEmpty()) {
    95              // transform each interval individually to handle wrap-around
    96              final List<AngularInterval> intervals = toIntervals();
    97
    98              setEmpty();
    99
    100             for (final AngularInterval interval : intervals) {
    101                 union(interval.transform(transform).toTree());
    102             }
    103         }
    104     }
    105
    106     /** {@inheritDoc}
    107      *
    108      * <p>It is important to note that split operations occur according to the rules of the
    109      * {@link CutAngle} hyperplane class. In this class, the continuous circle is viewed
    110      * as a non-circular segment of the number line in the range {@code [0, 2pi)}. Hyperplanes
    111      * are placed along this line and partition it into the segments {@code [0, x]}
    112      * and {@code [x, 2pi)}, where {@code x} is the location of the hyperplane. For example,
    113      * a positive-facing {@link CutAngle} instance with an azimuth of {@code 0.5pi} has
    114      * a minus side consisting of the angles {@code [0, 0.5pi]} and a plus side consisting of
    115      * the angles {@code [0.5pi, 2pi)}. Similarly, a positive-facing {@link CutAngle} with
    116      * an azimuth of {@code 0pi} has a plus side of {@code [0, 2pi)} (the full space) and
    117      * a minus side that is completely empty (since no points exist in our domain that are
    118      * less than zero). These rules can result in somewhat non-intuitive behavior at times.
    119      * For example, splitting a non-empty region with a hyperplane at {@code 0pi} is
    120      * essentially a no-op, since the region will either lie entirely on the plus or minus
    121      * side of the hyperplane (depending on the hyperplane's orientation) regardless of the actual
    122      * content of the region. In these situations, a copy of the tree is returned on the
    123      * appropriate side of the split.</p>
    124      *
    125      * @see CutAngle
    126      * @see #splitDiameter(CutAngle)
    127      */
    128     @Override
    129     public Split<RegionBSPTree1S> split(final Hyperplane<Point1S> splitter) {
    130         // Handle the special case where the cut is on the azimuth equivalent to zero.
    131         // In this case, it is not possible for any points to lie between it and zero.
    132         if (!isEmpty() && splitter.classify(Point1S.ZERO) == HyperplaneLocation.ON) {
    133             final CutAngle cut = (CutAngle) splitter;
    134             if (cut.isPositiveFacing()) {
    135                 return new Split<>(null, copy());
    136             } else {
    137                 return new Split<>(copy(), null);
    138             }
    139         }
    140
    141         return split(splitter, RegionBSPTree1S.empty(), RegionBSPTree1S.empty());
    142     }
    143
    144     /** Split the instance along a circle diameter.The diameter is defined by the given
    145      * split point and its reversed antipodal point.
    146      * @param splitter split point defining one side of the split diameter
    147      * @return result of the split operation
    148      */
    149     public Split<RegionBSPTree1S> splitDiameter(final CutAngle splitter) {
    150
    151         final CutAngle opposite = CutAngles.fromPointAndDirection(
    152                 splitter.getPoint().antipodal(),
    153                 !splitter.isPositiveFacing(),
    154                 splitter.getPrecision());
    155
    156         final double plusPoleOffset = splitter.isPositiveFacing() ?
    157                 +Angle.PI_OVER_TWO :
    158                 -Angle.PI_OVER_TWO;
    159         final Point1S plusPole = Point1S.of(splitter.getAzimuth() + plusPoleOffset);
    160
    161         final boolean zeroOnPlusSide = splitter.getPrecision()
    162                 .lte(plusPole.distance(Point1S.ZERO), Angle.PI_OVER_TWO);
    163
    164         final Split<RegionBSPTree1S> firstSplit = split(splitter);
    165         final Split<RegionBSPTree1S> secondSplit = split(opposite);
    166
    167         RegionBSPTree1S minus = RegionBSPTree1S.empty();
    168         RegionBSPTree1S plus = RegionBSPTree1S.empty();
    169
    170         if (zeroOnPlusSide) {
    171             // zero wrap-around needs to be handled on the plus side of the split
    172             safeUnion(plus, firstSplit.getPlus());
    173             safeUnion(plus, secondSplit.getPlus());
    174
    175             minus = firstSplit.getMinus();
    176             if (minus != null) {
    177                 minus = minus.split(opposite).getMinus();
    178             }
    179         } else {
    180             // zero wrap-around needs to be handled on the minus side of the split
    181             safeUnion(minus, firstSplit.getMinus());
    182             safeUnion(minus, secondSplit.getMinus());
    183
    184             plus = firstSplit.getPlus();
    185             if (plus != null) {
    186                 plus = plus.split(opposite).getPlus();
    187             }
    188         }
    189
    190         return new Split<>(
    191                 (minus != null && !minus.isEmpty()) ? minus : null,
    192                 (plus != null && !plus.isEmpty()) ? plus : null);
    193     }
    194
    195
    196     /** Convert the region represented by this tree into a list of separate
    197      * {@link AngularInterval}s, arranged in order of ascending min value.
    198      * @return list of {@link AngularInterval}s representing this region in order of
    199      *      ascending min value
    200      */
    201     public List<AngularInterval> toIntervals() {
    202         if (isFull()) {
    203             return Collections.singletonList(AngularInterval.full());
    204         }
    205
    206         final List<BoundaryPair> insideBoundaryPairs = new ArrayList<>();
    207         for (final RegionNode1S node : nodes()) {
    208             if (node.isInside()) {
    209                 insideBoundaryPairs.add(getNodeBoundaryPair(node));
    210             }
    211         }
    212
    213         insideBoundaryPairs.sort(BOUNDARY_PAIR_COMPARATOR);
    214
    215         final int boundaryPairCount = insideBoundaryPairs.size();
    216
    217         // Get the start point for merging intervals together.
    218         final int startOffset = getIntervalStartIndex(insideBoundaryPairs);
    219
    220         // Go through the pairs starting at the start offset and create intervals
    221         // for each set of adjacent pairs.
    222         final List<AngularInterval> intervals = new ArrayList<>();
    223
    224         BoundaryPair start = null;
    225         BoundaryPair end = null;
    226         BoundaryPair current;
    227
    228         for (int i = 0; i < boundaryPairCount; ++i) {
    229             current = insideBoundaryPairs.get((i + startOffset) % boundaryPairCount);
    230
    231             if (start == null) {
    232                 start = current;
    233                 end = current;
    234             } else if (Objects.equals(end.getMax(), current.getMin())) {
    235                 // these intervals should be merged
    236                 end = current;
    237             } else {
    238                 // these intervals should be separate
    239                 intervals.add(createInterval(start, end));
    240
    241                 // queue up the next pair
    242                 start = current;
    243                 end = current;
    244             }
    245         }
    246
    247         if (start != null && end != null) {
    248             intervals.add(createInterval(start, end));
    249         }
    250
    251         return intervals;
    252     }
    253
    254     /** Get the index that should be used as the starting point for combining adjacent boundary pairs
    255      * into contiguous intervals. This is computed as the first boundary pair found that is not connected
    256      * to the pair before it, or {@code 0} if no such entry exists.
    257      * @param boundaryPairs list of boundary pairs to search; must be ordered by increasing azimuth
    258      * @return the index to use as a starting place for combining adjacent boundary pairs
    259      *      into contiguous intervals
    260      */
    261     private int getIntervalStartIndex(final List<BoundaryPair> boundaryPairs) {
    262         final int size = boundaryPairs.size();
    263
    264         if (size > 0) {
    265             BoundaryPair current;
    266             BoundaryPair previous = boundaryPairs.get(size - 1);
    267
    268             for (int i = 0; i < size; ++i, previous = current) {
    269                 current = boundaryPairs.get(i);
    270
    271                 if (!Objects.equals(current.getMin(), previous.getMax())) {
    272                     return i;
    273                 }
    274             }
    275         }
    276
    277         return 0;
    278     }
    279
    280     /** Create an interval instance from the min boundary from the start boundary pair and
    281      * the max boundary from the end boundary pair. The hyperplane directions are adjusted
    282      * as needed.
    283      * @param start starting boundary pair
    284      * @param end ending boundary pair
    285      * @return an interval created from the min boundary of the given start pair and the
    286      *      max boundary from the given end pair
    287      */
    288     private AngularInterval createInterval(final BoundaryPair start, final BoundaryPair end) {
    289         CutAngle min = start.getMin();
    290         CutAngle max = end.getMax();
    291
    292         final Precision.DoubleEquivalence precision = (min != null) ? min.getPrecision() : max.getPrecision();
    293
    294         // flip the hyperplanes if needed since there's no
    295         // guarantee that the inside will be on the minus side
    296         // of the hyperplane (for example, if the region is complemented)
    297
    298         if (min != null) {
    299             if (min.isPositiveFacing()) {
    300                 min = min.reverse();
    301             }
    302         } else {
    303             min = CutAngles.createNegativeFacing(0.0, precision);
    304         }
    305
    306         if (max != null) {
    307             if (!max.isPositiveFacing()) {
    308                 max = max.reverse();
    309             }
    310         } else {
    311             max = CutAngles.createPositiveFacing(Angle.TWO_PI, precision);
    312         }
    313
    314         return AngularInterval.of(min, max);
    315     }
    316
    317     /** Return the min/max boundary pair for the convex region represented by the given node.
    318      * @param node the node to compute the interval for
    319      * @return the min/max boundary pair for the convex region represented by the given node
    320      */
    321     private BoundaryPair getNodeBoundaryPair(final RegionNode1S node) {
    322         CutAngle min = null;
    323         CutAngle max = null;
    324
    325         CutAngle pt;
    326         RegionNode1S child = node;
    327         RegionNode1S parent;
    328
    329         while ((min == null || max == null) && (parent = child.getParent()) != null) {
    330             pt = (CutAngle) parent.getCutHyperplane();
    331
    332             if ((pt.isPositiveFacing() && child.isMinus()) ||
    333                     (!pt.isPositiveFacing() && child.isPlus())) {
    334
    335                 if (max == null) {
    336                     max = pt;
    337                 }
    338             } else if (min == null) {
    339                 min = pt;
    340             }
    341
    342             child = parent;
    343         }
    344
    345         return new BoundaryPair(min, max);
    346     }
    347
    348     /** {@inheritDoc} */
    349     @Override
    350     protected RegionSizeProperties<Point1S> computeRegionSizeProperties() {
    351         if (isFull()) {
    352             return new RegionSizeProperties<>(Angle.TWO_PI, null);
    353         } else if (isEmpty()) {
    354             return new RegionSizeProperties<>(0, null);
    355         }
    356
    357         double size = 0;
    358         Vector2D scaledCentroidSum = Vector2D.ZERO;
    359
    360         double intervalSize;
    361
    362         for (final AngularInterval interval : toIntervals()) {
    363             intervalSize = interval.getSize();
    364
    365             size += intervalSize;
    366             scaledCentroidSum = scaledCentroidSum.add(interval.getCentroid().getVector().withNorm(intervalSize));
    367         }
    368
    369         final Precision.DoubleEquivalence precision = ((CutAngle) getRoot().getCutHyperplane()).getPrecision();
    370
    371         final Point1S centroid = scaledCentroidSum.eq(Vector2D.ZERO, precision) ?
    372                  null :
    373                  Point1S.from(scaledCentroidSum);
    374
    375         return new RegionSizeProperties<>(size, centroid);
    376     }
    377
    378     /** {@inheritDoc} */
    379     @Override
    380     protected RegionNode1S createNode() {
    381         return new RegionNode1S(this);
    382     }
    383
    384     /** Return a new, empty BSP tree.
    385      * @return a new, empty BSP tree.
    386      */
    387     public static RegionBSPTree1S empty() {
    388         return new RegionBSPTree1S(false);
    389     }
    390
    391     /** Return a new, full BSP tree. The returned tree represents the
    392      * full space.
    393      * @return a new, full BSP tree.
    394      */
    395     public static RegionBSPTree1S full() {
    396         return new RegionBSPTree1S(true);
    397     }
    398
    399     /** Return a new BSP tree representing the same region as the given angular interval.
    400      * @param interval the input interval
    401      * @return a new BSP tree representing the same region as the given angular interval
    402      */
    403     public static RegionBSPTree1S fromInterval(final AngularInterval interval) {
    404         final CutAngle minBoundary = interval.getMinBoundary();
    405         final CutAngle maxBoundary = interval.getMaxBoundary();
    406
    407         final RegionBSPTree1S tree = full();
    408
    409         if (minBoundary != null) {
    410             tree.insert(minBoundary.span());
    411         }
    412
    413         if (maxBoundary != null) {
    414             tree.insert(maxBoundary.span());
    415         }
    416
    417         return tree;
    418     }
    419
    420     /** Perform a union operation with {@code target} and {@code input}, storing the result
    421      * in {@code target}; does nothing if {@code input} is null.
    422      * @param target target tree
    423      * @param input input tree
    424      */
    425     private static void safeUnion(final RegionBSPTree1S target, final RegionBSPTree1S input) {
    426         if (input != null) {
    427             target.union(input);
    428         }
    429     }
    430
    431     /** BSP tree node for one dimensional spherical space.
    432      */
    433     public static final class RegionNode1S extends AbstractRegionBSPTree.AbstractRegionNode<Point1S, RegionNode1S> {
    434         /** Simple constructor.
    435          * @param tree the owning tree instance
    436          */
    437         private RegionNode1S(final AbstractBSPTree<Point1S, RegionNode1S> tree) {
    438             super(tree);
    439         }
    440
    441         /** {@inheritDoc} */
    442         @Override
    443         protected RegionNode1S getSelf() {
    444             return this;
    445         }
    446     }
    447
    448     /** Internal class containing pairs of interval boundaries.
    449      */
    450     private static final class BoundaryPair {
    451
    452         /** The min boundary. */
    453         private final CutAngle min;
    454
    455         /** The max boundary. */
    456         private final CutAngle max;
    457
    458         /** Simple constructor.
    459          * @param min min boundary hyperplane
    460          * @param max max boundary hyperplane
    461          */
    462         BoundaryPair(final CutAngle min, final CutAngle max) {
    463             this.min = min;
    464             this.max = max;
    465         }
    466
    467         /** Get the minimum boundary hyperplane.
    468          * @return the minimum boundary hyperplane.
    469          */
    470         public CutAngle getMin() {
    471             return min;
    472         }
    473
    474         /** Get the maximum boundary hyperplane.
    475          * @return the maximum boundary hyperplane.
    476          */
    477         public CutAngle getMax() {
    478             return max;
    479         }
    480
    481         /** Get the minimum value of the interval or zero if no minimum value exists.
    482          * @return the minimum value of the interval or zero
    483          *      if no minimum value exists.
    484          */
    485         public double getMinValue() {
    486             return (min != null) ? min.getNormalizedAzimuth() : 0;
    487         }
    488     }
    489
    490     /** Class used to project points onto the region boundary.
    491      */
    492     private static final class BoundaryProjector1S extends BoundaryProjector<Point1S, RegionNode1S> {
    493         /** Simple constructor.
    494          * @param point the point to project onto the region's boundary
    495          */
    496         BoundaryProjector1S(final Point1S point) {
    497             super(point);
    498         }
    499
    500         /** {@inheritDoc} */
    501         @Override
    502         protected boolean isPossibleClosestCut(final HyperplaneSubset<Point1S> cut, final Point1S target,
    503                 final double minDist) {
    504             // since the space wraps around, consider any cut as possibly being the closest
    505             return true;
    506         }
    507
    508         /** {@inheritDoc} */
    509         @Override
    510         protected Point1S disambiguateClosestPoint(final Point1S target, final Point1S a, final Point1S b) {
    511             // prefer the point with the smaller normalize azimuth value
    512             return a.getNormalizedAzimuth() < b.getNormalizedAzimuth() ? a : b;
    513         }
    514     }
    515 }