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}