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