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.geometry.spherical.oned;
018
019import java.util.ArrayList;
020import java.util.Collections;
021import java.util.Comparator;
022import java.util.List;
023import java.util.Objects;
024
025import org.apache.commons.geometry.core.Transform;
026import org.apache.commons.geometry.core.partitioning.Hyperplane;
027import org.apache.commons.geometry.core.partitioning.HyperplaneLocation;
028import org.apache.commons.geometry.core.partitioning.HyperplaneSubset;
029import org.apache.commons.geometry.core.partitioning.Split;
030import org.apache.commons.geometry.core.partitioning.bsp.AbstractBSPTree;
031import org.apache.commons.geometry.core.partitioning.bsp.AbstractRegionBSPTree;
032import org.apache.commons.geometry.euclidean.twod.Vector2D;
033import org.apache.commons.numbers.angle.Angle;
034import org.apache.commons.numbers.core.Precision;
035
036/** BSP tree representing regions in 1D spherical space.
037 */
038public class RegionBSPTree1S extends AbstractRegionBSPTree<Point1S, RegionBSPTree1S.RegionNode1S> {
039    /** Comparator used to sort BoundaryPairs by ascending azimuth.  */
040    private static final Comparator<BoundaryPair> BOUNDARY_PAIR_COMPARATOR =
041            Comparator.comparingDouble(BoundaryPair::getMinValue);
042
043    /** Create a new, empty instance.
044     */
045    public RegionBSPTree1S() {
046        this(false);
047    }
048
049    /** Create a new region. If {@code full} is true, then the region will
050     * represent the entire circle. Otherwise, it will be empty.
051     * @param full whether or not the region should contain the entire
052     *      circle or be empty
053     */
054    public RegionBSPTree1S(final boolean full) {
055        super(full);
056    }
057
058    /** Return a deep copy of this instance.
059     * @return a deep copy of this instance.
060     * @see #copy(org.apache.commons.geometry.core.partitioning.bsp.BSPTree)
061     */
062    public RegionBSPTree1S copy() {
063        final RegionBSPTree1S result = RegionBSPTree1S.empty();
064        result.copy(this);
065
066        return result;
067    }
068
069    /** Add an interval to this region. The resulting region will be the
070     * union of the interval and the region represented by this instance.
071     * @param interval the interval to add
072     */
073    public void add(final AngularInterval interval) {
074        union(fromInterval(interval));
075    }
076
077    /** {@inheritDoc} */
078    @Override
079    public Point1S project(final Point1S pt) {
080        final BoundaryProjector1S projector = new BoundaryProjector1S(pt);
081        accept(projector);
082
083        return projector.getProjected();
084    }
085
086    /** {@inheritDoc}
087     *
088     * <p>Each interval of the region is transformed individually and the
089     * results are unioned. If the size of any transformed interval is greater
090     * than or equal to 2pi, then the region is set to the full space.</p>
091     */
092    @Override
093    public void transform(final Transform<Point1S> transform) {
094        if (!isFull() && !isEmpty()) {
095            // transform each interval individually to handle wrap-around
096            final List<AngularInterval> intervals = toIntervals();
097
098            setEmpty();
099
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}