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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.geometry.euclidean.oned;
23
24 import java.util.ArrayList;
25 import java.util.Collection;
26 import java.util.Iterator;
27 import java.util.List;
28 import java.util.NoSuchElementException;
29
30 import org.hipparchus.geometry.partitioning.AbstractRegion;
31 import org.hipparchus.geometry.partitioning.BSPTree;
32 import org.hipparchus.geometry.partitioning.BoundaryProjection;
33 import org.hipparchus.util.Precision;
34
35 /** This class represents a 1D region: a set of intervals.
36 */
37 public class IntervalsSet
38 extends AbstractRegion<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint,
39 Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint>
40 implements Iterable<double[]> {
41
42 /** Build an intervals set representing the whole real line.
43 * @param tolerance tolerance below which points are considered identical.
44 */
45 public IntervalsSet(final double tolerance) {
46 super(tolerance);
47 }
48
49 /** Build an intervals set corresponding to a single interval.
50 * @param lower lower bound of the interval, must be lesser or equal
51 * to {@code upper} (may be {@code Double.NEGATIVE_INFINITY})
52 * @param upper upper bound of the interval, must be greater or equal
53 * to {@code lower} (may be {@code Double.POSITIVE_INFINITY})
54 * @param tolerance tolerance below which points are considered identical.
55 */
56 public IntervalsSet(final double lower, final double upper, final double tolerance) {
57 super(buildTree(lower, upper, tolerance), tolerance);
58 }
59
60 /** Build an intervals set from an inside/outside BSP tree.
61 * <p>The leaf nodes of the BSP tree <em>must</em> have a
62 * {@code Boolean} attribute representing the inside status of
63 * the corresponding cell (true for inside cells, false for outside
64 * cells). In order to avoid building too many small objects, it is
65 * recommended to use the predefined constants
66 * {@code Boolean.TRUE} and {@code Boolean.FALSE}</p>
67 * @param tree inside/outside BSP tree representing the intervals set
68 * @param tolerance tolerance below which points are considered identical.
69 */
70 public IntervalsSet(final BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> tree,
71 final double tolerance) {
72 super(tree, tolerance);
73 }
74
75 /** Build an intervals set from a Boundary REPresentation (B-rep).
76 * <p>The boundary is provided as a collection of {@link
77 * org.hipparchus.geometry.partitioning.SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the
78 * interior part of the region on its minus side and the exterior on
79 * its plus side.</p>
80 * <p>The boundary elements can be in any order, and can form
81 * several non-connected sets (like for example polygons with holes
82 * or a set of disjoints polyhedrons considered as a whole). In
83 * fact, the elements do not even need to be connected together
84 * (their topological connections are not used here). However, if the
85 * boundary does not really separate an inside open from an outside
86 * open (open having here its topological meaning), then subsequent
87 * calls to the {@link
88 * org.hipparchus.geometry.partitioning.Region#checkPoint(org.hipparchus.geometry.Point)
89 * checkPoint} method will not be meaningful anymore.</p>
90 * <p>If the boundary is empty, the region will represent the whole
91 * space.</p>
92 * @param boundary collection of boundary elements
93 * @param tolerance tolerance below which points are considered identical.
94 */
95 public IntervalsSet(final Collection<SubOrientedPoint> boundary, final double tolerance) {
96 super(boundary, tolerance);
97 }
98
99 /** Build an inside/outside tree representing a single interval.
100 * @param lower lower bound of the interval, must be lesser or equal
101 * to {@code upper} (may be {@code Double.NEGATIVE_INFINITY})
102 * @param upper upper bound of the interval, must be greater or equal
103 * to {@code lower} (may be {@code Double.POSITIVE_INFINITY})
104 * @param tolerance tolerance below which points are considered identical.
105 * @return the built tree
106 */
107 private static BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint>
108 buildTree(final double lower, final double upper, final double tolerance) {
109 if (Double.isInfinite(lower) && (lower < 0)) {
110 if (Double.isInfinite(upper) && (upper > 0)) {
111 // the tree must cover the whole real line
112 return new BSPTree<>(Boolean.TRUE);
113 }
114 // the tree must be open on the negative infinity side
115 final SubOrientedPoint upperCut = new OrientedPoint(new Vector1D(upper), true, tolerance).wholeHyperplane();
116 return new BSPTree<>(upperCut, new BSPTree<>(Boolean.FALSE), new BSPTree<>(Boolean.TRUE), null);
117 }
118 final SubOrientedPoint lowerCut = new OrientedPoint(new Vector1D(lower), false, tolerance).wholeHyperplane();
119 if (Double.isInfinite(upper) && (upper > 0)) {
120 // the tree must be open on the positive infinity side
121 return new BSPTree<>(lowerCut, new BSPTree<>(Boolean.FALSE), new BSPTree<>(Boolean.TRUE), null);
122 }
123
124 // the tree must be bounded on the two sides
125 final SubOrientedPoint upperCut = new OrientedPoint(new Vector1D(upper), true, tolerance).wholeHyperplane();
126 return new BSPTree<>(lowerCut,
127 new BSPTree<>(Boolean.FALSE),
128 new BSPTree<>(upperCut, new BSPTree<>(Boolean.FALSE), new BSPTree<>(Boolean.TRUE), null),
129 null);
130
131 }
132
133 /** {@inheritDoc} */
134 @Override
135 public IntervalsSet buildNew(final BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> tree) {
136 return new IntervalsSet(tree, getTolerance());
137 }
138
139 /** {@inheritDoc} */
140 @Override
141 public Vector1D getInteriorPoint() {
142
143 // look for the midpoint of the longest interval
144 // or some finite point if interval extends to ±∞
145 double selectedPoint = Double.NaN;
146 double selectedLength = 0;
147 for (final double[] a : this) {
148 final double length = a[1] - a[0];
149 if (length > selectedLength) {
150 // this interval is longer than the selected one, change selection
151 if (Double.isInfinite(a[0])) {
152 // interval starts at -∞, it may extend to +∞ as well
153 selectedPoint = Double.isInfinite(a[1]) ? 0 : a[1] - 1.0e9 * getTolerance();
154 } else if (Double.isInfinite(a[1])) {
155 // interval ends at +∞
156 selectedPoint = a[0] + 1.0e9 * getTolerance();
157 } else {
158 selectedPoint = 0.5 * (a[0] + a[1]);
159 }
160 selectedLength = length;
161 }
162 }
163
164 return Double.isNaN(selectedPoint) ? null : new Vector1D(selectedPoint);
165
166 }
167
168 /** {@inheritDoc} */
169 @Override
170 protected void computeGeometricalProperties() {
171 if (getTree(false).getCut() == null) {
172 setBarycenter(Vector1D.NaN);
173 setSize(((Boolean) getTree(false).getAttribute()) ? Double.POSITIVE_INFINITY : 0);
174 } else {
175 double size = 0.0;
176 double sum = 0.0;
177 for (final Interval interval : asList()) {
178 size += interval.getSize();
179 sum += interval.getSize() * interval.getBarycenter();
180 }
181 setSize(size);
182 if (Double.isInfinite(size)) {
183 setBarycenter(Vector1D.NaN);
184 } else if (size >= Precision.SAFE_MIN) {
185 setBarycenter(new Vector1D(sum / size));
186 } else {
187 setBarycenter(getTree(false).getCut().getHyperplane().getLocation());
188 }
189 }
190 }
191
192 /** Get the lowest value belonging to the instance.
193 * @return lowest value belonging to the instance
194 * ({@code Double.NEGATIVE_INFINITY} if the instance doesn't
195 * have any low bound, {@code Double.POSITIVE_INFINITY} if the
196 * instance is empty)
197 */
198 public double getInf() {
199 BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> node = getTree(false);
200 double inf = Double.POSITIVE_INFINITY;
201 while (node.getCut() != null) {
202 final OrientedPoint op = node.getCut().getHyperplane();
203 inf = op.getLocation().getX();
204 node = op.isDirect() ? node.getMinus() : node.getPlus();
205 }
206 return ((Boolean) node.getAttribute()) ? Double.NEGATIVE_INFINITY : inf;
207 }
208
209 /** Get the highest value belonging to the instance.
210 * @return highest value belonging to the instance
211 * ({@code Double.POSITIVE_INFINITY} if the instance doesn't
212 * have any high bound, {@code Double.NEGATIVE_INFINITY} if the
213 * instance is empty)
214 */
215 public double getSup() {
216 BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> node = getTree(false);
217 double sup = Double.NEGATIVE_INFINITY;
218 while (node.getCut() != null) {
219 final OrientedPoint op = node.getCut().getHyperplane();
220 sup = op.getLocation().getX();
221 node = op.isDirect() ? node.getPlus() : node.getMinus();
222 }
223 return ((Boolean) node.getAttribute()) ? Double.POSITIVE_INFINITY : sup;
224 }
225
226 /** {@inheritDoc}
227 */
228 @Override
229 public BoundaryProjection<Euclidean1D, Vector1D> projectToBoundary(final Vector1D point) {
230
231 // get position of test point
232 final double x = point.getX();
233
234 double previous = Double.NEGATIVE_INFINITY;
235 for (final double[] a : this) {
236 if (x < a[0]) {
237 // the test point lies between the previous and the current intervals
238 // offset will be positive
239 final double previousOffset = x - previous;
240 final double currentOffset = a[0] - x;
241 if (previousOffset < currentOffset) {
242 return new BoundaryProjection<>(point, finiteOrNullPoint(previous), previousOffset);
243 } else {
244 return new BoundaryProjection<>(point, finiteOrNullPoint(a[0]), currentOffset);
245 }
246 } else if (x <= a[1]) {
247 // the test point lies within the current interval
248 // offset will be negative
249 final double offset0 = a[0] - x;
250 final double offset1 = x - a[1];
251 if (offset0 < offset1) {
252 return new BoundaryProjection<>(point, finiteOrNullPoint(a[1]), offset1);
253 } else {
254 return new BoundaryProjection<>(point, finiteOrNullPoint(a[0]), offset0);
255 }
256 }
257 previous = a[1];
258 }
259
260 // the test point if past the last sub-interval
261 return new BoundaryProjection<>(point, finiteOrNullPoint(previous), x - previous);
262
263 }
264
265 /** Build a finite point.
266 * @param x abscissa of the point
267 * @return a new point for finite abscissa, null otherwise
268 */
269 private Vector1D finiteOrNullPoint(final double x) {
270 return Double.isInfinite(x) ? null : new Vector1D(x);
271 }
272
273 /** Build an ordered list of intervals representing the instance.
274 * <p>This method builds this intervals set as an ordered list of
275 * {@link Interval Interval} elements. If the intervals set has no
276 * lower limit, the first interval will have its low bound equal to
277 * {@code Double.NEGATIVE_INFINITY}. If the intervals set has
278 * no upper limit, the last interval will have its upper bound equal
279 * to {@code Double.POSITIVE_INFINITY}. An empty tree will
280 * build an empty list while a tree representing the whole real line
281 * will build a one element list with both bounds being
282 * infinite.</p>
283 * @return a new ordered list containing {@link Interval Interval}
284 * elements
285 */
286 public List<Interval> asList() {
287 final List<Interval> list = new ArrayList<>();
288 for (final double[] a : this) {
289 list.add(new Interval(a[0], a[1]));
290 }
291 return list;
292 }
293
294 /** Get the first leaf node of a tree.
295 * @param root tree root
296 * @return first leaf node
297 */
298 private BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint>
299 getFirstLeaf(final BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> root) {
300
301 if (root.getCut() == null) {
302 return root;
303 }
304
305 // find the smallest internal node
306 BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> smallest = null;
307 for (BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> n = root; n != null; n = previousInternalNode(n)) {
308 smallest = n;
309 }
310
311 return leafBefore(smallest);
312
313 }
314
315 /** Get the node corresponding to the first interval boundary.
316 * @return smallest internal node,
317 * or null if there are no internal nodes (i.e. the set is either empty or covers the real line)
318 */
319 private BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> getFirstIntervalBoundary() {
320
321 // start search at the tree root
322 BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> node = getTree(false);
323 if (node.getCut() == null) {
324 return null;
325 }
326
327 // walk tree until we find the smallest internal node
328 node = getFirstLeaf(node).getParent();
329
330 // walk tree until we find an interval boundary
331 while (node != null && !(isIntervalStart(node) || isIntervalEnd(node))) {
332 node = nextInternalNode(node);
333 }
334
335 return node;
336
337 }
338
339 /** Check if an internal node corresponds to the start abscissa of an interval.
340 * @param node internal node to check
341 * @return true if the node corresponds to the start abscissa of an interval
342 */
343 private boolean isIntervalStart(final BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> node) {
344
345 if ((Boolean) leafBefore(node).getAttribute()) {
346 // it has an inside cell before it, it may end an interval but not start it
347 return false;
348 }
349
350 if (!(Boolean) leafAfter(node).getAttribute()) {
351 // it has an outside cell after it, it is a dummy cut away from real intervals
352 return false;
353 }
354
355 // the cell has an outside before and an inside after it,
356 // it is the start of an interval
357 return true;
358
359 }
360
361 /** Check if an internal node corresponds to the end abscissa of an interval.
362 * @param node internal node to check
363 * @return true if the node corresponds to the end abscissa of an interval
364 */
365 private boolean isIntervalEnd(final BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> node) {
366
367 if (!(Boolean) leafBefore(node).getAttribute()) {
368 // it has an outside cell before it, it may start an interval but not end it
369 return false;
370 }
371
372 if ((Boolean) leafAfter(node).getAttribute()) {
373 // it has an inside cell after it, it is a dummy cut in the middle of an interval
374 return false;
375 }
376
377 // the cell has an inside before and an outside after it,
378 // it is the end of an interval
379 return true;
380
381 }
382
383 /** Get the next internal node.
384 * @param node current internal node
385 * @return next internal node in ascending order, or null
386 * if this is the last internal node
387 */
388 private BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint>
389 nextInternalNode(BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> node) {
390
391 if (childAfter(node).getCut() != null) {
392 // the next node is in the subtree
393 return leafAfter(node).getParent();
394 }
395
396 // there is nothing left deeper in the tree, we backtrack
397 while (isAfterParent(node)) {
398 node = node.getParent();
399 }
400 return node.getParent();
401
402 }
403
404 /** Get the previous internal node.
405 * @param node current internal node
406 * @return previous internal node in ascending order, or null
407 * if this is the first internal node
408 */
409 private BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint>
410 previousInternalNode(BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> node) {
411
412 if (childBefore(node).getCut() != null) {
413 // the next node is in the subtree
414 return leafBefore(node).getParent();
415 }
416
417 // there is nothing left deeper in the tree, we backtrack
418 while (isBeforeParent(node)) {
419 node = node.getParent();
420 }
421 return node.getParent();
422
423 }
424
425 /** Find the leaf node just before an internal node.
426 * @param node internal node at which the subtree starts
427 * @return leaf node just before the internal node
428 */
429 private BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint>
430 leafBefore(BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> node) {
431
432 node = childBefore(node);
433 while (node.getCut() != null) {
434 node = childAfter(node);
435 }
436
437 return node;
438
439 }
440
441 /** Find the leaf node just after an internal node.
442 * @param node internal node at which the subtree starts
443 * @return leaf node just after the internal node
444 */
445 private BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint>
446 leafAfter(BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> node) {
447
448 node = childAfter(node);
449 while (node.getCut() != null) {
450 node = childBefore(node);
451 }
452
453 return node;
454
455 }
456
457 /** Check if a node is the child before its parent in ascending order.
458 * @param node child node considered
459 * @return true is the node has a parent end is before it in ascending order
460 */
461 private boolean isBeforeParent(final BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> node) {
462 final BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> parent = node.getParent();
463 if (parent == null) {
464 return false;
465 } else {
466 return node == childBefore(parent);
467 }
468 }
469
470 /** Check if a node is the child after its parent in ascending order.
471 * @param node child node considered
472 * @return true is the node has a parent end is after it in ascending order
473 */
474 private boolean isAfterParent(final BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> node) {
475 final BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> parent = node.getParent();
476 if (parent == null) {
477 return false;
478 } else {
479 return node == childAfter(parent);
480 }
481 }
482
483 /** Find the child node just before an internal node.
484 * @param node internal node at which the subtree starts
485 * @return child node just before the internal node
486 */
487 private BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint>
488 childBefore(BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> node) {
489 if (isDirect(node)) {
490 // smaller abscissas are on minus side, larger abscissas are on plus side
491 return node.getMinus();
492 } else {
493 // smaller abscissas are on plus side, larger abscissas are on minus side
494 return node.getPlus();
495 }
496 }
497
498 /** Find the child node just after an internal node.
499 * @param node internal node at which the subtree starts
500 * @return child node just after the internal node
501 */
502 private BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint>
503 childAfter(BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> node) {
504 if (isDirect(node)) {
505 // smaller abscissas are on minus side, larger abscissas are on plus side
506 return node.getPlus();
507 } else {
508 // smaller abscissas are on plus side, larger abscissas are on minus side
509 return node.getMinus();
510 }
511 }
512
513 /** Check if an internal node has a direct oriented point.
514 * @param node internal node to check
515 * @return true if the oriented point is direct
516 */
517 private boolean isDirect(final BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> node) {
518 return node.getCut().getHyperplane().isDirect();
519 }
520
521 /** Get the abscissa of an internal node.
522 * @param node internal node to check
523 * @return abscissa
524 */
525 private double getAngle(final BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> node) {
526 return node.getCut().getHyperplane().getLocation().getX();
527 }
528
529 /** {@inheritDoc}
530 * <p>
531 * The iterator returns the limit values of sub-intervals in ascending order.
532 * </p>
533 * <p>
534 * The iterator does <em>not</em> support the optional {@code remove} operation.
535 * </p>
536 */
537 @Override
538 public Iterator<double[]> iterator() {
539 return new SubIntervalsIterator();
540 }
541
542 /** Local iterator for sub-intervals. */
543 private class SubIntervalsIterator implements Iterator<double[]> {
544
545 /** Current node. */
546 private BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> current;
547
548 /** Sub-interval no yet returned. */
549 private double[] pending;
550
551 /** Simple constructor.
552 */
553 SubIntervalsIterator() {
554
555 current = getFirstIntervalBoundary();
556
557 if (current == null) {
558 // all the leaf tree nodes share the same inside/outside status
559 if ((Boolean) getFirstLeaf(getTree(false)).getAttribute()) {
560 // it is an inside node, it represents the full real line
561 pending = new double[] {
562 Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY
563 };
564 } else {
565 pending = null;
566 }
567 } else if (isIntervalEnd(current)) {
568 // the first boundary is an interval end,
569 // so the first interval starts at infinity
570 pending = new double[] {
571 Double.NEGATIVE_INFINITY, getAngle(current)
572 };
573 } else {
574 selectPending();
575 }
576 }
577
578 /** Walk the tree to select the pending sub-interval.
579 */
580 private void selectPending() {
581
582 // look for the start of the interval
583 BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> start = current;
584 while (start != null && !isIntervalStart(start)) {
585 start = nextInternalNode(start);
586 }
587
588 if (start == null) {
589 // we have exhausted the iterator
590 current = null;
591 pending = null;
592 return;
593 }
594
595 // look for the end of the interval
596 BSPTree<Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> end = start;
597 while (end != null && !isIntervalEnd(end)) {
598 end = nextInternalNode(end);
599 }
600
601 if (end != null) {
602
603 // we have identified the interval
604 pending = new double[] {
605 getAngle(start), getAngle(end)
606 };
607
608 // prepare search for next interval
609 current = end;
610
611 } else {
612
613 // the final interval is open toward infinity
614 pending = new double[] {
615 getAngle(start), Double.POSITIVE_INFINITY
616 };
617
618 // there won't be any other intervals
619 current = null;
620
621 }
622
623 }
624
625 /** {@inheritDoc} */
626 @Override
627 public boolean hasNext() {
628 return pending != null;
629 }
630
631 /** {@inheritDoc} */
632 @Override
633 public double[] next() {
634 if (pending == null) {
635 throw new NoSuchElementException();
636 }
637 final double[] next = pending;
638 selectPending();
639 return next;
640 }
641
642 /** {@inheritDoc} */
643 @Override
644 public void remove() {
645 throw new UnsupportedOperationException();
646 }
647
648 }
649
650 }