1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 package org.hipparchus.geometry.euclidean.threed;
23
24 import java.util.ArrayList;
25 import java.util.Arrays;
26 import java.util.Collection;
27 import java.util.List;
28
29 import org.hipparchus.exception.LocalizedCoreFormats;
30 import org.hipparchus.exception.MathIllegalArgumentException;
31 import org.hipparchus.exception.MathRuntimeException;
32 import org.hipparchus.geometry.LocalizedGeometryFormats;
33 import org.hipparchus.geometry.Point;
34 import org.hipparchus.geometry.euclidean.oned.Euclidean1D;
35 import org.hipparchus.geometry.euclidean.twod.Euclidean2D;
36 import org.hipparchus.geometry.euclidean.twod.PolygonsSet;
37 import org.hipparchus.geometry.euclidean.twod.SubLine;
38 import org.hipparchus.geometry.euclidean.twod.Vector2D;
39 import org.hipparchus.geometry.partitioning.AbstractRegion;
40 import org.hipparchus.geometry.partitioning.BSPTree;
41 import org.hipparchus.geometry.partitioning.BSPTreeVisitor;
42 import org.hipparchus.geometry.partitioning.BoundaryAttribute;
43 import org.hipparchus.geometry.partitioning.Hyperplane;
44 import org.hipparchus.geometry.partitioning.Region;
45 import org.hipparchus.geometry.partitioning.RegionFactory;
46 import org.hipparchus.geometry.partitioning.SubHyperplane;
47 import org.hipparchus.geometry.partitioning.Transform;
48 import org.hipparchus.util.FastMath;
49
50
51
52 public class PolyhedronsSet extends AbstractRegion<Euclidean3D, Euclidean2D> {
53
54
55
56
57 public PolyhedronsSet(final double tolerance) {
58 super(tolerance);
59 }
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81 public PolyhedronsSet(final BSPTree<Euclidean3D> tree, final double tolerance) {
82 super(tree, tolerance);
83 }
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105 public PolyhedronsSet(final Collection<SubHyperplane<Euclidean3D>> boundary,
106 final double tolerance) {
107 super(boundary, tolerance);
108 }
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127 public PolyhedronsSet(final List<Vector3D> vertices, final List<int[]> facets,
128 final double tolerance) {
129 super(buildBoundary(vertices, facets, tolerance), tolerance);
130 }
131
132
133
134
135
136
137
138
139
140
141
142
143 public PolyhedronsSet(final BRep brep, final double tolerance) {
144 super(buildBoundary(brep.getVertices(), brep.getFacets(), tolerance), tolerance);
145 }
146
147
148
149
150
151
152
153
154
155
156 public PolyhedronsSet(final double xMin, final double xMax,
157 final double yMin, final double yMax,
158 final double zMin, final double zMax,
159 final double tolerance) {
160 super(buildBoundary(xMin, xMax, yMin, yMax, zMin, zMax, tolerance), tolerance);
161 }
162
163
164
165
166
167
168
169
170
171
172
173 private static BSPTree<Euclidean3D> buildBoundary(final double xMin, final double xMax,
174 final double yMin, final double yMax,
175 final double zMin, final double zMax,
176 final double tolerance) {
177 if ((xMin >= xMax - tolerance) || (yMin >= yMax - tolerance) || (zMin >= zMax - tolerance)) {
178
179 return new BSPTree<Euclidean3D>(Boolean.FALSE);
180 }
181 final Plane pxMin = new Plane(new Vector3D(xMin, 0, 0), Vector3D.MINUS_I, tolerance);
182 final Plane pxMax = new Plane(new Vector3D(xMax, 0, 0), Vector3D.PLUS_I, tolerance);
183 final Plane pyMin = new Plane(new Vector3D(0, yMin, 0), Vector3D.MINUS_J, tolerance);
184 final Plane pyMax = new Plane(new Vector3D(0, yMax, 0), Vector3D.PLUS_J, tolerance);
185 final Plane pzMin = new Plane(new Vector3D(0, 0, zMin), Vector3D.MINUS_K, tolerance);
186 final Plane pzMax = new Plane(new Vector3D(0, 0, zMax), Vector3D.PLUS_K, tolerance);
187 final Region<Euclidean3D> boundary =
188 new RegionFactory<Euclidean3D>().buildConvex(pxMin, pxMax, pyMin, pyMax, pzMin, pzMax);
189 return boundary.getTree(false);
190 }
191
192
193
194
195
196
197
198
199 private static List<SubHyperplane<Euclidean3D>> buildBoundary(final List<Vector3D> vertices,
200 final List<int[]> facets,
201 final double tolerance) {
202
203
204 for (int i = 0; i < vertices.size() - 1; ++i) {
205 final Vector3D vi = vertices.get(i);
206 for (int j = i + 1; j < vertices.size(); ++j) {
207 if (Vector3D.distance(vi, vertices.get(j)) <= tolerance) {
208 throw new MathIllegalArgumentException(LocalizedGeometryFormats.CLOSE_VERTICES,
209 vi.getX(), vi.getY(), vi.getZ());
210 }
211 }
212 }
213
214
215 final int[][] references = findReferences(vertices, facets);
216
217
218 final int[][] successors = successors(vertices, facets, references);
219
220
221 for (int vA = 0; vA < vertices.size(); ++vA) {
222 for (final int vB : successors[vA]) {
223
224 if (vB >= 0) {
225
226
227 boolean found = false;
228 for (final int v : successors[vB]) {
229 found = found || (v == vA);
230 }
231 if (!found) {
232 final Vector3D start = vertices.get(vA);
233 final Vector3D end = vertices.get(vB);
234 throw new MathIllegalArgumentException(LocalizedGeometryFormats.EDGE_CONNECTED_TO_ONE_FACET,
235 start.getX(), start.getY(), start.getZ(),
236 end.getX(), end.getY(), end.getZ());
237 }
238 }
239 }
240 }
241
242 final List<SubHyperplane<Euclidean3D>> boundary = new ArrayList<>();
243
244 for (final int[] facet : facets) {
245
246
247 Plane plane = new Plane(vertices.get(facet[0]), vertices.get(facet[1]), vertices.get(facet[2]),
248 tolerance);
249
250
251 final Vector2D[] two2Points = new Vector2D[facet.length];
252 for (int i = 0 ; i < facet.length; ++i) {
253 final Vector3D v = vertices.get(facet[i]);
254 if (!plane.contains(v)) {
255 throw new MathIllegalArgumentException(LocalizedGeometryFormats.OUT_OF_PLANE,
256 v.getX(), v.getY(), v.getZ());
257 }
258 two2Points[i] = plane.toSubSpace(v);
259 }
260
261
262 boundary.add(new SubPlane(plane, new PolygonsSet(tolerance, two2Points)));
263
264 }
265
266 return boundary;
267
268 }
269
270
271
272
273
274
275
276 private static int[][] findReferences(final List<Vector3D> vertices, final List<int[]> facets) {
277
278
279 final int[] nbFacets = new int[vertices.size()];
280 int maxFacets = 0;
281 for (final int[] facet : facets) {
282 if (facet.length < 3) {
283 throw new MathIllegalArgumentException(LocalizedCoreFormats.WRONG_NUMBER_OF_POINTS,
284 3, facet.length, true);
285 }
286 for (final int index : facet) {
287 maxFacets = FastMath.max(maxFacets, ++nbFacets[index]);
288 }
289 }
290
291
292 final int[][] references = new int[vertices.size()][maxFacets];
293 for (int[] r : references) {
294 Arrays.fill(r, -1);
295 }
296 for (int f = 0; f < facets.size(); ++f) {
297 for (final int v : facets.get(f)) {
298
299 int k = 0;
300 while (k < maxFacets && references[v][k] >= 0) {
301 ++k;
302 }
303 references[v][k] = f;
304 }
305 }
306
307 return references;
308
309 }
310
311
312
313
314
315
316
317
318
319
320 private static int[][] successors(final List<Vector3D> vertices, final List<int[]> facets,
321 final int[][] references) {
322
323
324 final int[][] successors = new int[vertices.size()][references[0].length];
325 for (final int[] s : successors) {
326 Arrays.fill(s, -1);
327 }
328
329 for (int v = 0; v < vertices.size(); ++v) {
330 for (int k = 0; k < successors[v].length && references[v][k] >= 0; ++k) {
331
332
333 final int[] facet = facets.get(references[v][k]);
334 int i = 0;
335 while (i < facet.length && facet[i] != v) {
336 ++i;
337 }
338
339
340 successors[v][k] = facet[(i + 1) % facet.length];
341 for (int l = 0; l < k; ++l) {
342 if (successors[v][l] == successors[v][k]) {
343 final Vector3D start = vertices.get(v);
344 final Vector3D end = vertices.get(successors[v][k]);
345 throw new MathIllegalArgumentException(LocalizedGeometryFormats.FACET_ORIENTATION_MISMATCH,
346 start.getX(), start.getY(), start.getZ(),
347 end.getX(), end.getY(), end.getZ());
348 }
349 }
350
351 }
352 }
353
354 return successors;
355
356 }
357
358
359 @Override
360 public PolyhedronsSet buildNew(final BSPTree<Euclidean3D> tree) {
361 return new PolyhedronsSet(tree, getTolerance());
362 }
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386 public BRep getBRep() throws MathRuntimeException {
387 BRepExtractor extractor = new BRepExtractor(getTolerance());
388 getTree(true).visit(extractor);
389 return extractor.getBRep();
390 }
391
392
393 private static class BRepExtractor implements BSPTreeVisitor<Euclidean3D> {
394
395
396 private final double tolerance;
397
398
399 private final List<Vector3D> vertices;
400
401
402 private final List<int[]> facets;
403
404
405
406
407 BRepExtractor(final double tolerance) {
408 this.tolerance = tolerance;
409 this.vertices = new ArrayList<>();
410 this.facets = new ArrayList<>();
411 }
412
413
414
415
416 public BRep getBRep() {
417 return new BRep(vertices, facets);
418 }
419
420
421 @Override
422 public Order visitOrder(final BSPTree<Euclidean3D> node) {
423 return Order.MINUS_SUB_PLUS;
424 }
425
426
427 @Override
428 public void visitInternalNode(final BSPTree<Euclidean3D> node) {
429 @SuppressWarnings("unchecked")
430 final BoundaryAttribute<Euclidean3D> attribute =
431 (BoundaryAttribute<Euclidean3D>) node.getAttribute();
432 if (attribute.getPlusOutside() != null) {
433 addContribution(attribute.getPlusOutside(), false);
434 }
435 if (attribute.getPlusInside() != null) {
436 addContribution(attribute.getPlusInside(), true);
437 }
438 }
439
440
441 @Override
442 public void visitLeafNode(final BSPTree<Euclidean3D> node) {
443 }
444
445
446
447
448
449
450 private void addContribution(final SubHyperplane<Euclidean3D> facet, final boolean reversed)
451 throws MathRuntimeException {
452
453 final Plane plane = (Plane) facet.getHyperplane();
454 final PolygonsSet polygon = (PolygonsSet) ((SubPlane) facet).getRemainingRegion();
455 final Vector2D[][] loops2D = polygon.getVertices();
456 if (loops2D.length == 0) {
457 throw new MathRuntimeException(LocalizedGeometryFormats.OUTLINE_BOUNDARY_LOOP_OPEN);
458 } else if (loops2D.length > 1) {
459 throw new MathRuntimeException(LocalizedGeometryFormats.FACET_WITH_SEVERAL_BOUNDARY_LOOPS);
460 } else {
461 for (final Vector2D[] loop2D : polygon.getVertices()) {
462 final int[] loop3D = new int[loop2D.length];
463 for (int i = 0; i < loop2D.length ; ++i) {
464 if (loop2D[i] == null) {
465 throw new MathRuntimeException(LocalizedGeometryFormats.OUTLINE_BOUNDARY_LOOP_OPEN);
466 }
467 loop3D[reversed ? loop2D.length - 1 - i : i] = getVertexIndex(plane.toSpace(loop2D[i]));
468 }
469 facets.add(loop3D);
470 }
471 }
472
473 }
474
475
476
477
478
479 private int getVertexIndex(final Vector3D vertex) {
480
481 for (int i = 0; i < vertices.size(); ++i) {
482 if (Vector3D.distance(vertex, vertices.get(i)) <= tolerance) {
483
484 return i;
485 }
486 }
487
488
489 vertices.add(vertex);
490 return vertices.size() - 1;
491
492 }
493
494 }
495
496
497 @Override
498 protected void computeGeometricalProperties() {
499
500
501 getTree(true).visit(new FacetsContributionVisitor());
502
503 if (getSize() < 0) {
504
505
506 setSize(Double.POSITIVE_INFINITY);
507 setBarycenter((Point<Euclidean3D>) Vector3D.NaN);
508 } else {
509
510 setSize(getSize() / 3.0);
511 setBarycenter((Point<Euclidean3D>) new Vector3D(1.0 / (4 * getSize()), (Vector3D) getBarycenter()));
512 }
513
514 }
515
516
517 private class FacetsContributionVisitor implements BSPTreeVisitor<Euclidean3D> {
518
519
520 FacetsContributionVisitor() {
521 setSize(0);
522 setBarycenter((Point<Euclidean3D>) new Vector3D(0, 0, 0));
523 }
524
525
526 @Override
527 public Order visitOrder(final BSPTree<Euclidean3D> node) {
528 return Order.MINUS_SUB_PLUS;
529 }
530
531
532 @Override
533 public void visitInternalNode(final BSPTree<Euclidean3D> node) {
534 @SuppressWarnings("unchecked")
535 final BoundaryAttribute<Euclidean3D> attribute =
536 (BoundaryAttribute<Euclidean3D>) node.getAttribute();
537 if (attribute.getPlusOutside() != null) {
538 addContribution(attribute.getPlusOutside(), false);
539 }
540 if (attribute.getPlusInside() != null) {
541 addContribution(attribute.getPlusInside(), true);
542 }
543 }
544
545
546 @Override
547 public void visitLeafNode(final BSPTree<Euclidean3D> node) {
548 }
549
550
551
552
553
554 private void addContribution(final SubHyperplane<Euclidean3D> facet, final boolean reversed) {
555
556 final Region<Euclidean2D> polygon = ((SubPlane) facet).getRemainingRegion();
557 final double area = polygon.getSize();
558
559 if (Double.isInfinite(area)) {
560 setSize(Double.POSITIVE_INFINITY);
561 setBarycenter((Point<Euclidean3D>) Vector3D.NaN);
562 } else {
563
564 final Plane plane = (Plane) facet.getHyperplane();
565 final Vector3D facetB = plane.toSpace(polygon.getBarycenter());
566 double scaled = area * facetB.dotProduct(plane.getNormal());
567 if (reversed) {
568 scaled = -scaled;
569 }
570
571 setSize(getSize() + scaled);
572 setBarycenter((Point<Euclidean3D>) new Vector3D(1.0, (Vector3D) getBarycenter(), scaled, facetB));
573
574 }
575
576 }
577
578 }
579
580
581
582
583
584
585
586
587 public SubHyperplane<Euclidean3D> firstIntersection(final Vector3D point, final Line line) {
588 return recurseFirstIntersection(getTree(true), point, line);
589 }
590
591
592
593
594
595
596
597
598
599 private SubHyperplane<Euclidean3D> recurseFirstIntersection(final BSPTree<Euclidean3D> node,
600 final Vector3D point,
601 final Line line) {
602
603 final SubHyperplane<Euclidean3D> cut = node.getCut();
604 if (cut == null) {
605 return null;
606 }
607 final BSPTree<Euclidean3D> minus = node.getMinus();
608 final BSPTree<Euclidean3D> plus = node.getPlus();
609 final Plane plane = (Plane) cut.getHyperplane();
610
611
612 final double offset = plane.getOffset((Point<Euclidean3D>) point);
613 final boolean in = FastMath.abs(offset) < getTolerance();
614 final BSPTree<Euclidean3D> near;
615 final BSPTree<Euclidean3D> far;
616 if (offset < 0) {
617 near = minus;
618 far = plus;
619 } else {
620 near = plus;
621 far = minus;
622 }
623
624 if (in) {
625
626 final SubHyperplane<Euclidean3D> facet = boundaryFacet(point, node);
627 if (facet != null) {
628 return facet;
629 }
630 }
631
632
633 final SubHyperplane<Euclidean3D> crossed = recurseFirstIntersection(near, point, line);
634 if (crossed != null) {
635 return crossed;
636 }
637
638 if (!in) {
639
640 final Vector3D hit3D = plane.intersection(line);
641 if (hit3D != null && line.getAbscissa(hit3D) > line.getAbscissa(point)) {
642 final SubHyperplane<Euclidean3D> facet = boundaryFacet(hit3D, node);
643 if (facet != null) {
644 return facet;
645 }
646 }
647 }
648
649
650 return recurseFirstIntersection(far, point, line);
651
652 }
653
654
655
656
657
658
659
660 private SubHyperplane<Euclidean3D> boundaryFacet(final Vector3D point,
661 final BSPTree<Euclidean3D> node) {
662 final Vector2D point2D = ((Plane) node.getCut().getHyperplane()).toSubSpace((Point<Euclidean3D>) point);
663 @SuppressWarnings("unchecked")
664 final BoundaryAttribute<Euclidean3D> attribute =
665 (BoundaryAttribute<Euclidean3D>) node.getAttribute();
666 if ((attribute.getPlusOutside() != null) &&
667 (((SubPlane) attribute.getPlusOutside()).getRemainingRegion().checkPoint(point2D) != Location.OUTSIDE)) {
668 return attribute.getPlusOutside();
669 }
670 if ((attribute.getPlusInside() != null) &&
671 (((SubPlane) attribute.getPlusInside()).getRemainingRegion().checkPoint(point2D) != Location.OUTSIDE)) {
672 return attribute.getPlusInside();
673 }
674 return null;
675 }
676
677
678
679
680
681
682
683 public PolyhedronsSet rotate(final Vector3D center, final Rotation rotation) {
684 return (PolyhedronsSet) applyTransform(new RotationTransform(center, rotation));
685 }
686
687
688 private static class RotationTransform implements Transform<Euclidean3D, Euclidean2D> {
689
690
691 private final Vector3D center;
692
693
694 private final Rotation rotation;
695
696
697 private Plane cachedOriginal;
698
699
700 private Transform<Euclidean2D, Euclidean1D> cachedTransform;
701
702
703
704
705
706 RotationTransform(final Vector3D center, final Rotation rotation) {
707 this.center = center;
708 this.rotation = rotation;
709 }
710
711
712 @Override
713 public Vector3D apply(final Point<Euclidean3D> point) {
714 final Vector3D delta = ((Vector3D) point).subtract(center);
715 return new Vector3D(1.0, center, 1.0, rotation.applyTo(delta));
716 }
717
718
719 @Override
720 public Plane apply(final Hyperplane<Euclidean3D> hyperplane) {
721 return ((Plane) hyperplane).rotate(center, rotation);
722 }
723
724
725 @Override
726 public SubHyperplane<Euclidean2D> apply(final SubHyperplane<Euclidean2D> sub,
727 final Hyperplane<Euclidean3D> original,
728 final Hyperplane<Euclidean3D> transformed) {
729 if (original != cachedOriginal) {
730
731
732 final Plane oPlane = (Plane) original;
733 final Plane tPlane = (Plane) transformed;
734 final Vector3D p00 = oPlane.getOrigin();
735 final Vector3D p10 = oPlane.toSpace((Point<Euclidean2D>) new Vector2D(1.0, 0.0));
736 final Vector3D p01 = oPlane.toSpace((Point<Euclidean2D>) new Vector2D(0.0, 1.0));
737 final Vector2D tP00 = tPlane.toSubSpace((Point<Euclidean3D>) apply(p00));
738 final Vector2D tP10 = tPlane.toSubSpace((Point<Euclidean3D>) apply(p10));
739 final Vector2D tP01 = tPlane.toSubSpace((Point<Euclidean3D>) apply(p01));
740
741 cachedOriginal = (Plane) original;
742 cachedTransform =
743 org.hipparchus.geometry.euclidean.twod.Line.getTransform(tP10.getX() - tP00.getX(),
744 tP10.getY() - tP00.getY(),
745 tP01.getX() - tP00.getX(),
746 tP01.getY() - tP00.getY(),
747 tP00.getX(),
748 tP00.getY());
749
750 }
751 return ((SubLine) sub).applyTransform(cachedTransform);
752 }
753
754 }
755
756
757
758
759
760
761 public PolyhedronsSet translate(final Vector3D translation) {
762 return (PolyhedronsSet) applyTransform(new TranslationTransform(translation));
763 }
764
765
766 private static class TranslationTransform implements Transform<Euclidean3D, Euclidean2D> {
767
768
769 private final Vector3D translation;
770
771
772 private Plane cachedOriginal;
773
774
775 private Transform<Euclidean2D, Euclidean1D> cachedTransform;
776
777
778
779
780 TranslationTransform(final Vector3D translation) {
781 this.translation = translation;
782 }
783
784
785 @Override
786 public Vector3D apply(final Point<Euclidean3D> point) {
787 return new Vector3D(1.0, (Vector3D) point, 1.0, translation);
788 }
789
790
791 @Override
792 public Plane apply(final Hyperplane<Euclidean3D> hyperplane) {
793 return ((Plane) hyperplane).translate(translation);
794 }
795
796
797 @Override
798 public SubHyperplane<Euclidean2D> apply(final SubHyperplane<Euclidean2D> sub,
799 final Hyperplane<Euclidean3D> original,
800 final Hyperplane<Euclidean3D> transformed) {
801 if (original != cachedOriginal) {
802
803
804 final Plane oPlane = (Plane) original;
805 final Plane tPlane = (Plane) transformed;
806 final Vector2D shift = tPlane.toSubSpace((Point<Euclidean3D>) apply(oPlane.getOrigin()));
807
808 cachedOriginal = (Plane) original;
809 cachedTransform =
810 org.hipparchus.geometry.euclidean.twod.Line.getTransform(1, 0, 0, 1,
811 shift.getX(),
812 shift.getY());
813
814 }
815
816 return ((SubLine) sub).applyTransform(cachedTransform);
817
818 }
819
820 }
821
822
823
824
825
826
827
828
829
830
831
832
833 public static class BRep {
834
835
836 private final List<Vector3D> vertices;
837
838
839 private final List<int[]> facets;
840
841
842
843
844
845 public BRep(final List<Vector3D> vertices, final List<int[]> facets) {
846 this.vertices = vertices;
847 this.facets = facets;
848 }
849
850
851
852
853 public List<Vector3D> getVertices() {
854 return vertices;
855 }
856
857
858
859
860 public List<int[]> getFacets() {
861 return facets;
862 }
863
864 }
865
866 }