View Javadoc
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.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  /** This class represents a 3D region: a set of polyhedrons.
54   */
55  public class PolyhedronsSet
56      extends AbstractRegion<Euclidean3D, Vector3D, Plane, SubPlane,
57                             Euclidean2D, Vector2D, org.hipparchus.geometry.euclidean.twod.Line, SubLine> {
58  
59      /** Build a polyhedrons set representing the whole real line.
60       * @param tolerance tolerance below which points are considered identical
61       */
62      public PolyhedronsSet(final double tolerance) {
63          super(tolerance);
64      }
65  
66      /** Build a polyhedrons set from a BSP tree.
67       * <p>The leaf nodes of the BSP tree <em>must</em> have a
68       * {@code Boolean} attribute representing the inside status of
69       * the corresponding cell (true for inside cells, false for outside
70       * cells). In order to avoid building too many small objects, it is
71       * recommended to use the predefined constants
72       * {@code Boolean.TRUE} and {@code Boolean.FALSE}</p>
73       * <p>
74       * This constructor is aimed at expert use, as building the tree may
75       * be a difficult task. It is not intended for general use and for
76       * performances reasons does not check thoroughly its input, as this would
77       * require walking the full tree each time. Failing to provide a tree with
78       * the proper attributes, <em>will</em> therefore generate problems like
79       * {@link NullPointerException} or {@link ClassCastException} only later on.
80       * This limitation is known and explains why this constructor is for expert
81       * use only. The caller does have the responsibility to provided correct arguments.
82       * </p>
83       * @param tree inside/outside BSP tree representing the region
84       * @param tolerance tolerance below which points are considered identical
85       */
86      public PolyhedronsSet(final BSPTree<Euclidean3D, Vector3D, Plane, SubPlane> tree, final double tolerance) {
87          super(tree, tolerance);
88      }
89  
90      /** Build a polyhedrons set from a Boundary REPresentation (B-rep) specified by sub-hyperplanes.
91       * <p>The boundary is provided as a collection of {@link
92       * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the
93       * interior part of the region on its minus side and the exterior on
94       * its plus side.</p>
95       * <p>The boundary elements can be in any order, and can form
96       * several non-connected sets (like for example polyhedrons with holes
97       * or a set of disjoint polyhedrons considered as a whole). In
98       * fact, the elements do not even need to be connected together
99       * (their topological connections are not used here). However, if the
100      * boundary does not really separate an inside open from an outside
101      * open (open having here its topological meaning), then subsequent
102      * calls to the {@link Region#checkPoint(Point) checkPoint} method will
103      * not be meaningful anymore.</p>
104      * <p>If the boundary is empty, the region will represent the whole
105      * space.</p>
106      * @param boundary collection of boundary elements, as a
107      * collection of {@link SubHyperplane SubHyperplane} objects
108      * @param tolerance tolerance below which points are considered identical
109      */
110     public PolyhedronsSet(final Collection<SubPlane> boundary,
111                           final double tolerance) {
112         super(boundary, tolerance);
113     }
114 
115     /** Build a polyhedrons set from a Boundary REPresentation (B-rep) specified by connected vertices.
116      * <p>
117      * The boundary is provided as a list of vertices and a list of facets.
118      * Each facet is specified as an integer array containing the arrays vertices
119      * indices in the vertices list. Each facet normal is oriented by right hand
120      * rule to the facet vertices list.
121      * </p>
122      * <p>
123      * Some basic sanity checks are performed but not everything is thoroughly
124      * assessed, so it remains under caller responsibility to ensure the vertices
125      * and facets are consistent and properly define a polyhedrons set.
126      * </p>
127      * @param vertices list of polyhedrons set vertices
128      * @param facets list of facets, as vertices indices in the vertices list
129      * @param tolerance tolerance below which points are considered identical
130      * @exception MathIllegalArgumentException if some basic sanity checks fail
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     /** Build a polyhedrons set from a Boundary REPresentation (B-rep) specified by connected vertices.
138      * <p>
139      * Some basic sanity checks are performed but not everything is thoroughly
140      * assessed, so it remains under caller responsibility to ensure the vertices
141      * and facets are consistent and properly define a polyhedrons set.
142      * </p>
143      * @param brep Boundary REPresentation of the polyhedron to build
144      * @param tolerance tolerance below which points are considered identical
145      * @exception MathIllegalArgumentException if some basic sanity checks fail
146      * @since 1.2
147      */
148     public PolyhedronsSet(final BRep brep, final double tolerance) {
149         super(buildBoundary(brep.getVertices(), brep.getFacets(), tolerance), tolerance);
150     }
151 
152     /** Build a parallellepipedic box.
153      * @param xMin low bound along the x direction
154      * @param xMax high bound along the x direction
155      * @param yMin low bound along the y direction
156      * @param yMax high bound along the y direction
157      * @param zMin low bound along the z direction
158      * @param zMax high bound along the z direction
159      * @param tolerance tolerance below which points are considered identical
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     /** Build a parallellepipedic box boundary.
169      * @param xMin low bound along the x direction
170      * @param xMax high bound along the x direction
171      * @param yMin low bound along the y direction
172      * @param yMax high bound along the y direction
173      * @param zMin low bound along the z direction
174      * @param zMax high bound along the z direction
175      * @param tolerance tolerance below which points are considered identical
176      * @return boundary tree
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             // too thin box, build an empty polygons set
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     /** Build boundary from vertices and facets.
198      * @param vertices list of polyhedrons set vertices
199      * @param facets list of facets, as vertices indices in the vertices list
200      * @param tolerance tolerance below which points are considered identical
201      * @return boundary as a list of sub-hyperplanes
202      * @exception MathIllegalArgumentException if some basic sanity checks fail
203      */
204     private static List<SubPlane> buildBoundary(final List<Vector3D> vertices, final List<int[]> facets, final double tolerance) {
205 
206         // check vertices distances
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         // find how vertices are referenced by facets
218         final int[][] references = findReferences(vertices, facets);
219 
220         // find how vertices are linked together by edges along the facets they belong to
221         final int[][] successors = successors(vertices, facets, references);
222 
223         // check edges orientations
224         for (int vA = 0; vA < vertices.size(); ++vA) {
225             for (final int vB : successors[vA]) {
226 
227                 if (vB >= 0) {
228                     // when facets are properly oriented, if vB is the successor of vA on facet f1,
229                     // then there must be an adjacent facet f2 where vA is the successor of vB
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             // define facet plane from the first 3 points
250             Plane plane = new Plane(vertices.get(facet[0]), vertices.get(facet[1]), vertices.get(facet[2]),
251                                     tolerance);
252 
253             // check all points are in the plane
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             // create the polygonal facet
265             boundary.add(new SubPlane(plane, new PolygonsSet(tolerance, two2Points)));
266 
267         }
268 
269         return boundary;
270 
271     }
272 
273     /** Find the facets that reference each edges.
274      * @param vertices list of polyhedrons set vertices
275      * @param facets list of facets, as vertices indices in the vertices list
276      * @return references array such that r[v][k] = f for some k if facet f contains vertex v
277      * @exception MathIllegalArgumentException if some facets have fewer than 3 vertices
278      */
279     private static int[][] findReferences(final List<Vector3D> vertices, final List<int[]> facets) {
280 
281         // find the maximum number of facets a vertex belongs to
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         // set up the references array
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                 // vertex v is referenced by facet f
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     /** Find the successors of all vertices among all facets they belong to.
315      * @param vertices list of polyhedrons set vertices
316      * @param facets list of facets, as vertices indices in the vertices list
317      * @param references facets references array
318      * @return indices of vertices that follow vertex v in some facet (the array
319      * may contain extra entries at the end, set to negative indices)
320      * @exception MathIllegalArgumentException if the same vertex appears more than
321      * once in the successors list (which means one facet orientation is wrong)
322      */
323     private static int[][] successors(final List<Vector3D> vertices, final List<int[]> facets,
324                                       final int[][] references) {
325 
326         // create an array large enough
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                 // look for vertex v
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                 // we have found vertex v, we deduce its successor on current facet
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     /** {@inheritDoc} */
362     @Override
363     public PolyhedronsSet buildNew(final BSPTree<Euclidean3D, Vector3D, Plane, SubPlane> tree) {
364         return new PolyhedronsSet(tree, getTolerance());
365     }
366 
367     /** {@inheritDoc} */
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     /** Get the boundary representation of the instance.
378      * <p>
379      * The boundary representation can be extracted <em>only</em> from
380      * bounded polyhedrons sets. If the polyhedrons set is unbounded,
381      * a {@link MathRuntimeException} will be thrown.
382      * </p>
383      * <p>
384      * The boundary representation extracted is not minimal, as for
385      * example canonical facets may be split into several smaller
386      * independent sub-facets sharing the same plane and connected by
387      * their edges.
388      * </p>
389      * <p>
390      * As the {@link BRep B-Rep} representation does not support
391      * facets with several boundary loops (for example facets with
392      * holes), an exception is triggered when attempting to extract
393      * B-Rep from such complex polyhedrons sets.
394      * </p>
395      * @return boundary representation of the instance
396      * @exception MathRuntimeException if polyhedrons is unbounded
397      * @since 1.2
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     /** Visitor extracting BRep. */
406     private static class BRepExtractor implements BSPTreeVisitor<Euclidean3D, Vector3D, Plane, SubPlane> {
407 
408         /** Tolerance for vertices identification. */
409         private final double tolerance;
410 
411         /** Extracted vertices. */
412         private final List<Vector3D> vertices;
413 
414         /** Extracted facets. */
415         private final List<int[]> facets;
416 
417         /** Simple constructor.
418          * @param tolerance tolerance for vertices identification
419          */
420         BRepExtractor(final double tolerance) {
421             this.tolerance = tolerance;
422             this.vertices  = new ArrayList<>();
423             this.facets    = new ArrayList<>();
424         }
425 
426         /** Get the BRep.
427          * @return extracted BRep
428          */
429         public BRep getBRep() {
430             return new BRep(vertices, facets);
431         }
432 
433         /** {@inheritDoc} */
434         @Override
435         public Order visitOrder(final BSPTree<Euclidean3D, Vector3D, Plane, SubPlane> node) {
436             return Order.MINUS_SUB_PLUS;
437         }
438 
439         /** {@inheritDoc} */
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         /** {@inheritDoc} */
454         @Override
455         public void visitLeafNode(final BSPTree<Euclidean3D, Vector3D, Plane, SubPlane> node) {
456         }
457 
458         /** Add the contribution of a boundary facet.
459          * @param facet boundary facet
460          * @param reversed if true, the facet has the inside on its plus side
461          * @exception MathRuntimeException if facet is unbounded
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         /** Get the index of a vertex.
488          * @param vertex vertex as a 3D point
489          * @return index of the vertex
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                     // the vertex is already known
496                     return i;
497                 }
498             }
499 
500             // the vertex is a new one, add it
501             vertices.add(vertex);
502             return vertices.size() - 1;
503 
504         }
505 
506     }
507 
508     /** {@inheritDoc} */
509     @Override
510     protected void computeGeometricalProperties() {
511 
512         // compute the contribution of all boundary facets
513         getTree(true).visit(new FacetsContributionVisitor());
514 
515         if (getSize() < 0) {
516             // the polyhedrons set as a finite outside
517             // surrounded by an infinite inside
518             setSize(Double.POSITIVE_INFINITY);
519             setBarycenter(Vector3D.NaN);
520         } else {
521             // the polyhedrons set is finite, apply the remaining scaling factors
522             setSize(getSize() / 3.0);
523             setBarycenter(new Vector3D(1.0 / (4 * getSize()), getBarycenter()));
524         }
525 
526     }
527 
528     /** Visitor computing geometrical properties. */
529     private class FacetsContributionVisitor implements BSPTreeVisitor<Euclidean3D, Vector3D, Plane, SubPlane> {
530 
531         /** Simple constructor. */
532         FacetsContributionVisitor() {
533             setSize(0);
534             setBarycenter(new Vector3D(0, 0, 0));
535         }
536 
537         /** {@inheritDoc} */
538         @Override
539         public Order visitOrder(final BSPTree<Euclidean3D, Vector3D, Plane, SubPlane> node) {
540             return Order.MINUS_SUB_PLUS;
541         }
542 
543         /** {@inheritDoc} */
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         /** {@inheritDoc} */
558         @Override
559         public void visitLeafNode(final BSPTree<Euclidean3D, Vector3D, Plane, SubPlane> node) {
560         }
561 
562         /** Add he contribution of a boundary facet.
563          * @param facet boundary facet
564          * @param reversed if true, the facet has the inside on its plus side
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     /** Get the first sub-hyperplane crossed by a semi-infinite line.
594      * @param point start point of the part of the line considered
595      * @param line line to consider (contains point)
596      * @return the first sub-hyperplane crossed by the line after the
597      * given point, or null if the line does not intersect any
598      * sub-hyperplane
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     /** Get the first sub-hyperplane crossed by a semi-infinite line.
605      * @param node current node
606      * @param point start point of the part of the line considered
607      * @param line line to consider (contains point)
608      * @return the first sub-hyperplane crossed by the line after the
609      * given point, or null if the line does not intersect any
610      * sub-hyperplane
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         // establish search order
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             // search in the cut hyperplane
637             final SubPlane facet = boundaryFacet(point, node);
638             if (facet != null) {
639                 return facet;
640             }
641         }
642 
643         // search in the near branch
644         final SubPlane crossed = recurseFirstIntersection(near, point, line);
645         if (crossed != null) {
646             return crossed;
647         }
648 
649         if (!in) {
650             // search in the cut hyperplane
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         // search in the far branch
661         return recurseFirstIntersection(far, point, line);
662 
663     }
664 
665     /** Check if a point belongs to the boundary part of a node.
666      * @param point point to check
667      * @param node node containing the boundary facet to check
668      * @return the boundary facet this points belongs to (or null if it
669      * does not belong to any boundary facet)
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     /** Rotate the region around the specified point.
688      * <p>The instance is not modified, a new instance is created.</p>
689      * @param center rotation center
690      * @param rotation vectorial rotation operator
691      * @return a new instance representing the rotated region
692      */
693     public PolyhedronsSet rotate(final Vector3D center, final Rotation rotation) {
694         return (PolyhedronsSet) applyTransform(new RotationTransform(center, rotation));
695     }
696 
697     /** 3D rotation as a Transform. */
698     private static class RotationTransform
699         implements Transform<Euclidean3D, Vector3D, Plane, SubPlane,
700                              Euclidean2D, Vector2D, org.hipparchus.geometry.euclidean.twod.Line, SubLine> {
701 
702         /** Center point of the rotation. */
703         private final Vector3D   center;
704 
705         /** Vectorial rotation. */
706         private final Rotation   rotation;
707 
708         /** Cached original hyperplane. */
709         private Plane cachedOriginal;
710 
711         /** Cached 2D transform valid inside the cached original hyperplane. */
712         private Transform<Euclidean2D, Vector2D, org.hipparchus.geometry.euclidean.twod.Line, SubLine,
713                           Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> cachedTransform;
714 
715         /** Build a rotation transform.
716          * @param center center point of the rotation
717          * @param rotation vectorial rotation
718          */
719         RotationTransform(final Vector3D center, final Rotation rotation) {
720             this.center   = center;
721             this.rotation = rotation;
722         }
723 
724         /** {@inheritDoc} */
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         /** {@inheritDoc} */
732         @Override
733         public Plane apply(final Plane hyperplane) {
734             return hyperplane.rotate(center, rotation);
735         }
736 
737         /** {@inheritDoc} */
738         @Override
739         public SubLine apply(final SubLine sub, final Plane original, final Plane transformed) {
740             if (original != cachedOriginal) {
741                 // we have changed hyperplane, reset the in-hyperplane transform
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     /** Translate the region by the specified amount.
765      * <p>The instance is not modified, a new instance is created.</p>
766      * @param translation translation to apply
767      * @return a new instance representing the translated region
768      */
769     public PolyhedronsSet translate(final Vector3D translation) {
770         return (PolyhedronsSet) applyTransform(new TranslationTransform(translation));
771     }
772 
773     /** 3D translation as a transform. */
774     private static class TranslationTransform
775         implements Transform<Euclidean3D, Vector3D, Plane, SubPlane,
776                              Euclidean2D, Vector2D, org.hipparchus.geometry.euclidean.twod.Line, SubLine> {
777 
778         /** Translation vector. */
779         private final Vector3D   translation;
780 
781         /** Cached original hyperplane. */
782         private Plane cachedOriginal;
783 
784         /** Cached 2D transform valid inside the cached original hyperplane. */
785         private Transform<Euclidean2D, Vector2D, org.hipparchus.geometry.euclidean.twod.Line, SubLine,
786                           Euclidean1D, Vector1D, OrientedPoint, SubOrientedPoint> cachedTransform;
787 
788         /** Build a translation transform.
789          * @param translation translation vector
790          */
791         TranslationTransform(final Vector3D translation) {
792             this.translation = translation;
793         }
794 
795         /** {@inheritDoc} */
796         @Override
797         public Vector3D apply(final Vector3D point) {
798             return new Vector3D(1.0, point, 1.0, translation);
799         }
800 
801         /** {@inheritDoc} */
802         @Override
803         public Plane apply(final Plane hyperplane) {
804             return hyperplane.translate(translation);
805         }
806 
807         /** {@inheritDoc} */
808         @Override
809         public SubLine apply(final SubLine sub, final Plane original, final Plane transformed) {
810             if (original != cachedOriginal) {
811                 // we have changed hyperplane, reset the in-hyperplane transform
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     /** Container for Boundary REPresentation (B-Rep).
829      * <p>
830      * The boundary is provided as a list of vertices and a list of facets.
831      * Each facet is specified as an integer array containing the arrays vertices
832      * indices in the vertices list. Each facet normal is oriented by right hand
833      * rule to the facet vertices list.
834      * </p>
835      * @see PolyhedronsSet#PolyhedronsSet(BSPTree, double)
836      * @see PolyhedronsSet#getBRep()
837      * @since 1.2
838      */
839     public static class BRep {
840 
841         /** List of polyhedrons set vertices. */
842         private final List<Vector3D> vertices;
843 
844         /** List of facets, as vertices indices in the vertices list. */
845         private final List<int[]> facets;
846 
847         /** Simple constructor.
848          * @param vertices list of polyhedrons set vertices
849          * @param facets list of facets, as vertices indices in the vertices list
850          */
851         public BRep(final List<Vector3D> vertices, final List<int[]> facets) {
852             this.vertices = vertices;
853             this.facets   = facets;
854         }
855 
856         /** Get the extracted vertices.
857          * @return extracted vertices
858          */
859         public List<Vector3D> getVertices() {
860             return vertices;
861         }
862 
863         /** Get the extracted facets.
864          * @return extracted facets
865          */
866         public List<int[]> getFacets() {
867             return facets;
868         }
869 
870     }
871 
872 }