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.core.precision.DoublePrecisionContext;
    33  import org.apache.commons.geometry.euclidean.twod.Vector2D;
    34  import org.apache.commons.numbers.angle.PlaneAngleRadians;
    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                 +PlaneAngleRadians.PI_OVER_TWO :
    158                 -PlaneAngleRadians.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), PlaneAngleRadians.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         // Find the index of the first boundary pair that is not connected to pair before it.
    218         // This will be our start point for merging intervals together.
    219         int startOffset = 0;
    220         if (boundaryPairCount > 1) {
    221             BoundaryPair current = null;
    222             BoundaryPair previous = insideBoundaryPairs.get(boundaryPairCount - 1);
    223
    224             for (int i = 0; i < boundaryPairCount; ++i, previous = current) {
    225                 current = insideBoundaryPairs.get(i);
    226
    227                 if (!Objects.equals(current.getMin(), previous.getMax())) {
    228                     startOffset = i;
    229                     break;
    230                 }
    231             }
    232         }
    233
    234         // Go through the pairs starting at the start offset and create intervals
    235         // for each set of adjacent pairs.
    236         final List<AngularInterval> intervals = new ArrayList<>();
    237
    238         BoundaryPair start = null;
    239         BoundaryPair end = null;
    240         BoundaryPair current = null;
    241
    242         for (int i = 0; i < boundaryPairCount; ++i) {
    243             current = insideBoundaryPairs.get((i + startOffset) % boundaryPairCount);
    244
    245             if (start == null) {
    246                 start = current;
    247                 end = current;
    248             } else if (Objects.equals(end.getMax(), current.getMin())) {
    249                 // these intervals should be merged
    250                 end = current;
    251             } else {
    252                 // these intervals should be separate
    253                 intervals.add(createInterval(start, end));
    254
    255                 // queue up the next pair
    256                 start = current;
    257                 end = current;
    258             }
    259         }
    260
    261         if (start != null && end != null) {
    262             intervals.add(createInterval(start, end));
    263         }
    264
    265         return intervals;
    266     }
    267
    268     /** Create an interval instance from the min boundary from the start boundary pair and
    269      * the max boundary from the end boundary pair. The hyperplane directions are adjusted
    270      * as needed.
    271      * @param start starting boundary pair
    272      * @param end ending boundary pair
    273      * @return an interval created from the min boundary of the given start pair and the
    274      *      max boundary from the given end pair
    275      */
    276     private AngularInterval createInterval(final BoundaryPair start, final BoundaryPair end) {
    277         CutAngle min = start.getMin();
    278         CutAngle max = end.getMax();
    279
    280         final DoublePrecisionContext precision = (min != null) ? min.getPrecision() : max.getPrecision();
    281
    282         // flip the hyperplanes if needed since there's no
    283         // guarantee that the inside will be on the minus side
    284         // of the hyperplane (for example, if the region is complemented)
    285
    286         if (min != null) {
    287             if (min.isPositiveFacing()) {
    288                 min = min.reverse();
    289             }
    290         } else {
    291             min = CutAngles.createNegativeFacing(0.0, precision);
    292         }
    293
    294         if (max != null) {
    295             if (!max.isPositiveFacing()) {
    296                 max = max.reverse();
    297             }
    298         } else {
    299             max = CutAngles.createPositiveFacing(PlaneAngleRadians.TWO_PI, precision);
    300         }
    301
    302         return AngularInterval.of(min, max);
    303     }
    304
    305     /** Return the min/max boundary pair for the convex region represented by the given node.
    306      * @param node the node to compute the interval for
    307      * @return the min/max boundary pair for the convex region represented by the given node
    308      */
    309     private BoundaryPair getNodeBoundaryPair(final RegionNode1S node) {
    310         CutAngle min = null;
    311         CutAngle max = null;
    312
    313         CutAngle pt;
    314         RegionNode1S child = node;
    315         RegionNode1S parent;
    316
    317         while ((min == null || max == null) && (parent = child.getParent()) != null) {
    318             pt = (CutAngle) parent.getCutHyperplane();
    319
    320             if ((pt.isPositiveFacing() && child.isMinus()) ||
    321                     (!pt.isPositiveFacing() && child.isPlus())) {
    322
    323                 if (max == null) {
    324                     max = pt;
    325                 }
    326             } else if (min == null) {
    327                 min = pt;
    328             }
    329
    330             child = parent;
    331         }
    332
    333         return new BoundaryPair(min, max);
    334     }
    335
    336     /** {@inheritDoc} */
    337     @Override
    338     protected RegionSizeProperties<Point1S> computeRegionSizeProperties() {
    339         if (isFull()) {
    340             return new RegionSizeProperties<>(PlaneAngleRadians.TWO_PI, null);
    341         } else if (isEmpty()) {
    342             return new RegionSizeProperties<>(0, null);
    343         }
    344
    345         double size = 0;
    346         Vector2D scaledCentroidSum = Vector2D.ZERO;
    347
    348         double intervalSize;
    349
    350         for (final AngularInterval interval : toIntervals()) {
    351             intervalSize = interval.getSize();
    352
    353             size += intervalSize;
    354             scaledCentroidSum = scaledCentroidSum.add(interval.getCentroid().getVector().withNorm(intervalSize));
    355         }
    356
    357         final DoublePrecisionContext precision = ((CutAngle) getRoot().getCutHyperplane()).getPrecision();
    358
    359         final Point1S centroid = scaledCentroidSum.eq(Vector2D.ZERO, precision) ?
    360                  null :
    361                  Point1S.from(scaledCentroidSum);
    362
    363         return new RegionSizeProperties<>(size, centroid);
    364     }
    365
    366     /** {@inheritDoc} */
    367     @Override
    368     protected RegionNode1S createNode() {
    369         return new RegionNode1S(this);
    370     }
    371
    372     /** Return a new, empty BSP tree.
    373      * @return a new, empty BSP tree.
    374      */
    375     public static RegionBSPTree1S empty() {
    376         return new RegionBSPTree1S(false);
    377     }
    378
    379     /** Return a new, full BSP tree. The returned tree represents the
    380      * full space.
    381      * @return a new, full BSP tree.
    382      */
    383     public static RegionBSPTree1S full() {
    384         return new RegionBSPTree1S(true);
    385     }
    386
    387     /** Return a new BSP tree representing the same region as the given angular interval.
    388      * @param interval the input interval
    389      * @return a new BSP tree representing the same region as the given angular interval
    390      */
    391     public static RegionBSPTree1S fromInterval(final AngularInterval interval) {
    392         final CutAngle minBoundary = interval.getMinBoundary();
    393         final CutAngle maxBoundary = interval.getMaxBoundary();
    394
    395         final RegionBSPTree1S tree = full();
    396
    397         if (minBoundary != null) {
    398             tree.insert(minBoundary.span());
    399         }
    400
    401         if (maxBoundary != null) {
    402             tree.insert(maxBoundary.span());
    403         }
    404
    405         return tree;
    406     }
    407
    408     /** Perform a union operation with {@code target} and {@code input}, storing the result
    409      * in {@code target}; does nothing if {@code input} is null.
    410      * @param target target tree
    411      * @param input input tree
    412      */
    413     private static void safeUnion(final RegionBSPTree1S target, final RegionBSPTree1S input) {
    414         if (input != null) {
    415             target.union(input);
    416         }
    417     }
    418
    419     /** BSP tree node for one dimensional spherical space.
    420      */
    421     public static final class RegionNode1S extends AbstractRegionBSPTree.AbstractRegionNode<Point1S, RegionNode1S> {
    422         /** Simple constructor.
    423          * @param tree the owning tree instance
    424          */
    425         private RegionNode1S(final AbstractBSPTree<Point1S, RegionNode1S> tree) {
    426             super(tree);
    427         }
    428
    429         /** {@inheritDoc} */
    430         @Override
    431         protected RegionNode1S getSelf() {
    432             return this;
    433         }
    434     }
    435
    436     /** Internal class containing pairs of interval boundaries.
    437      */
    438     private static final class BoundaryPair {
    439
    440         /** The min boundary. */
    441         private final CutAngle min;
    442
    443         /** The max boundary. */
    444         private final CutAngle max;
    445
    446         /** Simple constructor.
    447          * @param min min boundary hyperplane
    448          * @param max max boundary hyperplane
    449          */
    450         BoundaryPair(final CutAngle min, final CutAngle max) {
    451             this.min = min;
    452             this.max = max;
    453         }
    454
    455         /** Get the minimum boundary hyperplane.
    456          * @return the minimum boundary hyperplane.
    457          */
    458         public CutAngle getMin() {
    459             return min;
    460         }
    461
    462         /** Get the maximum boundary hyperplane.
    463          * @return the maximum boundary hyperplane.
    464          */
    465         public CutAngle getMax() {
    466             return max;
    467         }
    468
    469         /** Get the minimum value of the interval or zero if no minimum value exists.
    470          * @return the minimum value of the interval or zero
    471          *      if no minimum value exists.
    472          */
    473         public double getMinValue() {
    474             return (min != null) ? min.getNormalizedAzimuth() : 0;
    475         }
    476     }
    477
    478     /** Class used to project points onto the region boundary.
    479      */
    480     private static final class BoundaryProjector1S extends BoundaryProjector<Point1S, RegionNode1S> {
    481         /** Simple constructor.
    482          * @param point the point to project onto the region's boundary
    483          */
    484         BoundaryProjector1S(final Point1S point) {
    485             super(point);
    486         }
    487
    488         /** {@inheritDoc} */
    489         @Override
    490         protected boolean isPossibleClosestCut(final HyperplaneSubset<Point1S> cut, final Point1S target,
    491                 final double minDist) {
    492             // since the space wraps around, consider any cut as possibly being the closest
    493             return true;
    494         }
    495
    496         /** {@inheritDoc} */
    497         @Override
    498         protected Point1S disambiguateClosestPoint(final Point1S target, final Point1S a, final Point1S b) {
    499             // prefer the point with the smaller normalize azimuth value
    500             return a.getNormalizedAzimuth() < b.getNormalizedAzimuth() ? a : b;
    501         }
    502     }
    503 }