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.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  /** This class represents a 3D region: a set of polyhedrons.
51   */
52  public class PolyhedronsSet extends AbstractRegion<Euclidean3D, Euclidean2D> {
53  
54      /** Build a polyhedrons set representing the whole real line.
55       * @param tolerance tolerance below which points are considered identical
56       */
57      public PolyhedronsSet(final double tolerance) {
58          super(tolerance);
59      }
60  
61      /** Build a polyhedrons set from a BSP tree.
62       * <p>The leaf nodes of the BSP tree <em>must</em> have a
63       * {@code Boolean} attribute representing the inside status of
64       * the corresponding cell (true for inside cells, false for outside
65       * cells). In order to avoid building too many small objects, it is
66       * recommended to use the predefined constants
67       * {@code Boolean.TRUE} and {@code Boolean.FALSE}</p>
68       * <p>
69       * This constructor is aimed at expert use, as building the tree may
70       * be a difficult task. It is not intended for general use and for
71       * performances reasons does not check thoroughly its input, as this would
72       * require walking the full tree each time. Failing to provide a tree with
73       * the proper attributes, <em>will</em> therefore generate problems like
74       * {@link NullPointerException} or {@link ClassCastException} only later on.
75       * This limitation is known and explains why this constructor is for expert
76       * use only. The caller does have the responsibility to provided correct arguments.
77       * </p>
78       * @param tree inside/outside BSP tree representing the region
79       * @param tolerance tolerance below which points are considered identical
80       */
81      public PolyhedronsSet(final BSPTree<Euclidean3D> tree, final double tolerance) {
82          super(tree, tolerance);
83      }
84  
85      /** Build a polyhedrons set from a Boundary REPresentation (B-rep) specified by sub-hyperplanes.
86       * <p>The boundary is provided as a collection of {@link
87       * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the
88       * interior part of the region on its minus side and the exterior on
89       * its plus side.</p>
90       * <p>The boundary elements can be in any order, and can form
91       * several non-connected sets (like for example polyhedrons with holes
92       * or a set of disjoint polyhedrons considered as a whole). In
93       * fact, the elements do not even need to be connected together
94       * (their topological connections are not used here). However, if the
95       * boundary does not really separate an inside open from an outside
96       * open (open having here its topological meaning), then subsequent
97       * calls to the {@link Region#checkPoint(Point) checkPoint} method will
98       * not be meaningful anymore.</p>
99       * <p>If the boundary is empty, the region will represent the whole
100      * space.</p>
101      * @param boundary collection of boundary elements, as a
102      * collection of {@link SubHyperplane SubHyperplane} objects
103      * @param tolerance tolerance below which points are considered identical
104      */
105     public PolyhedronsSet(final Collection<SubHyperplane<Euclidean3D>> boundary,
106                           final double tolerance) {
107         super(boundary, tolerance);
108     }
109 
110     /** Build a polyhedrons set from a Boundary REPresentation (B-rep) specified by connected vertices.
111      * <p>
112      * The boundary is provided as a list of vertices and a list of facets.
113      * Each facet is specified as an integer array containing the arrays vertices
114      * indices in the vertices list. Each facet normal is oriented by right hand
115      * rule to the facet vertices list.
116      * </p>
117      * <p>
118      * Some basic sanity checks are performed but not everything is thoroughly
119      * assessed, so it remains under caller responsibility to ensure the vertices
120      * and facets are consistent and properly define a polyhedrons set.
121      * </p>
122      * @param vertices list of polyhedrons set vertices
123      * @param facets list of facets, as vertices indices in the vertices list
124      * @param tolerance tolerance below which points are considered identical
125      * @exception MathIllegalArgumentException if some basic sanity checks fail
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     /** Build a polyhedrons set from a Boundary REPresentation (B-rep) specified by connected vertices.
133      * <p>
134      * Some basic sanity checks are performed but not everything is thoroughly
135      * assessed, so it remains under caller responsibility to ensure the vertices
136      * and facets are consistent and properly define a polyhedrons set.
137      * </p>
138      * @param brep Boundary REPresentation of the polyhedron to build
139      * @param tolerance tolerance below which points are considered identical
140      * @exception MathIllegalArgumentException if some basic sanity checks fail
141      * @since 1.2
142      */
143     public PolyhedronsSet(final BRep brep, final double tolerance) {
144         super(buildBoundary(brep.getVertices(), brep.getFacets(), tolerance), tolerance);
145     }
146 
147     /** Build a parallellepipedic box.
148      * @param xMin low bound along the x direction
149      * @param xMax high bound along the x direction
150      * @param yMin low bound along the y direction
151      * @param yMax high bound along the y direction
152      * @param zMin low bound along the z direction
153      * @param zMax high bound along the z direction
154      * @param tolerance tolerance below which points are considered identical
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     /** Build a parallellepipedic box boundary.
164      * @param xMin low bound along the x direction
165      * @param xMax high bound along the x direction
166      * @param yMin low bound along the y direction
167      * @param yMax high bound along the y direction
168      * @param zMin low bound along the z direction
169      * @param zMax high bound along the z direction
170      * @param tolerance tolerance below which points are considered identical
171      * @return boundary tree
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             // too thin box, build an empty polygons set
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     /** Build boundary from vertices and facets.
193      * @param vertices list of polyhedrons set vertices
194      * @param facets list of facets, as vertices indices in the vertices list
195      * @param tolerance tolerance below which points are considered identical
196      * @return boundary as a list of sub-hyperplanes
197      * @exception MathIllegalArgumentException if some basic sanity checks fail
198      */
199     private static List<SubHyperplane<Euclidean3D>> buildBoundary(final List<Vector3D> vertices,
200                                                                   final List<int[]> facets,
201                                                                   final double tolerance) {
202 
203         // check vertices distances
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         // find how vertices are referenced by facets
215         final int[][] references = findReferences(vertices, facets);
216 
217         // find how vertices are linked together by edges along the facets they belong to
218         final int[][] successors = successors(vertices, facets, references);
219 
220         // check edges orientations
221         for (int vA = 0; vA < vertices.size(); ++vA) {
222             for (final int vB : successors[vA]) {
223 
224                 if (vB >= 0) {
225                     // when facets are properly oriented, if vB is the successor of vA on facet f1,
226                     // then there must be an adjacent facet f2 where vA is the successor of vB
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             // define facet plane from the first 3 points
247             Plane plane = new Plane(vertices.get(facet[0]), vertices.get(facet[1]), vertices.get(facet[2]),
248                                     tolerance);
249 
250             // check all points are in the plane
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             // create the polygonal facet
262             boundary.add(new SubPlane(plane, new PolygonsSet(tolerance, two2Points)));
263 
264         }
265 
266         return boundary;
267 
268     }
269 
270     /** Find the facets that reference each edges.
271      * @param vertices list of polyhedrons set vertices
272      * @param facets list of facets, as vertices indices in the vertices list
273      * @return references array such that r[v][k] = f for some k if facet f contains vertex v
274      * @exception MathIllegalArgumentException if some facets have fewer than 3 vertices
275      */
276     private static int[][] findReferences(final List<Vector3D> vertices, final List<int[]> facets) {
277 
278         // find the maximum number of facets a vertex belongs to
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         // set up the references array
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                 // vertex v is referenced by facet f
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     /** Find the successors of all vertices among all facets they belong to.
312      * @param vertices list of polyhedrons set vertices
313      * @param facets list of facets, as vertices indices in the vertices list
314      * @param references facets references array
315      * @return indices of vertices that follow vertex v in some facet (the array
316      * may contain extra entries at the end, set to negative indices)
317      * @exception MathIllegalArgumentException if the same vertex appears more than
318      * once in the successors list (which means one facet orientation is wrong)
319      */
320     private static int[][] successors(final List<Vector3D> vertices, final List<int[]> facets,
321                                       final int[][] references) {
322 
323         // create an array large enough
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                 // look for vertex v
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                 // we have found vertex v, we deduce its successor on current facet
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     /** {@inheritDoc} */
359     @Override
360     public PolyhedronsSet buildNew(final BSPTree<Euclidean3D> tree) {
361         return new PolyhedronsSet(tree, getTolerance());
362     }
363 
364     /** Get the boundary representation of the instance.
365      * <p>
366      * The boundary representation can be extracted <em>only</em> from
367      * bounded polyhedrons sets. If the polyhedrons set is unbounded,
368      * a {@link MathRuntimeException} will be thrown.
369      * </p>
370      * <p>
371      * The boundary representation extracted is not minimal, as for
372      * example canonical facets may be split into several smaller
373      * independent sub-facets sharing the same plane and connected by
374      * their edges.
375      * </p>
376      * <p>
377      * As the {@link BRep B-Rep} representation does not support
378      * facets with several boundary loops (for example facets with
379      * holes), an exception is triggered when attempting to extract
380      * B-Rep from such complex polyhedrons sets.
381      * </p>
382      * @return boundary representation of the instance
383      * @exception MathRuntimeException if polyhedrons is unbounded
384      * @since 1.2
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     /** Visitor extracting BRep. */
393     private static class BRepExtractor implements BSPTreeVisitor<Euclidean3D> {
394 
395         /** Tolerance for vertices identification. */
396         private final double tolerance;
397 
398         /** Extracted vertices. */
399         private final List<Vector3D> vertices;
400 
401         /** Extracted facets. */
402         private final List<int[]> facets;
403 
404         /** Simple constructor.
405          * @param tolerance tolerance for vertices identification
406          */
407         BRepExtractor(final double tolerance) {
408             this.tolerance = tolerance;
409             this.vertices  = new ArrayList<>();
410             this.facets    = new ArrayList<>();
411         }
412 
413         /** Get the BRep.
414          * @return extracted BRep
415          */
416         public BRep getBRep() {
417             return new BRep(vertices, facets);
418         }
419 
420         /** {@inheritDoc} */
421         @Override
422         public Order visitOrder(final BSPTree<Euclidean3D> node) {
423             return Order.MINUS_SUB_PLUS;
424         }
425 
426         /** {@inheritDoc} */
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         /** {@inheritDoc} */
441         @Override
442         public void visitLeafNode(final BSPTree<Euclidean3D> node) {
443         }
444 
445         /** Add he contribution of a boundary facet.
446          * @param facet boundary facet
447          * @param reversed if true, the facet has the inside on its plus side
448          * @exception MathRuntimeException if facet is unbounded
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         /** Get the index of a vertex.
476          * @param vertex vertex as a 3D point
477          * @return index of the vertex
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                     // the vertex is already known
484                     return i;
485                 }
486             }
487 
488             // the vertex is a new one, add it
489             vertices.add(vertex);
490             return vertices.size() - 1;
491 
492         }
493 
494     }
495 
496     /** {@inheritDoc} */
497     @Override
498     protected void computeGeometricalProperties() {
499 
500         // compute the contribution of all boundary facets
501         getTree(true).visit(new FacetsContributionVisitor());
502 
503         if (getSize() < 0) {
504             // the polyhedrons set as a finite outside
505             // surrounded by an infinite inside
506             setSize(Double.POSITIVE_INFINITY);
507             setBarycenter((Point<Euclidean3D>) Vector3D.NaN);
508         } else {
509             // the polyhedrons set is finite, apply the remaining scaling factors
510             setSize(getSize() / 3.0);
511             setBarycenter((Point<Euclidean3D>) new Vector3D(1.0 / (4 * getSize()), (Vector3D) getBarycenter()));
512         }
513 
514     }
515 
516     /** Visitor computing geometrical properties. */
517     private class FacetsContributionVisitor implements BSPTreeVisitor<Euclidean3D> {
518 
519         /** Simple constructor. */
520         FacetsContributionVisitor() {
521             setSize(0);
522             setBarycenter((Point<Euclidean3D>) new Vector3D(0, 0, 0));
523         }
524 
525         /** {@inheritDoc} */
526         @Override
527         public Order visitOrder(final BSPTree<Euclidean3D> node) {
528             return Order.MINUS_SUB_PLUS;
529         }
530 
531         /** {@inheritDoc} */
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         /** {@inheritDoc} */
546         @Override
547         public void visitLeafNode(final BSPTree<Euclidean3D> node) {
548         }
549 
550         /** Add he contribution of a boundary facet.
551          * @param facet boundary facet
552          * @param reversed if true, the facet has the inside on its plus side
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     /** Get the first sub-hyperplane crossed by a semi-infinite line.
581      * @param point start point of the part of the line considered
582      * @param line line to consider (contains point)
583      * @return the first sub-hyperplane crossed by the line after the
584      * given point, or null if the line does not intersect any
585      * sub-hyperplane
586      */
587     public SubHyperplane<Euclidean3D> firstIntersection(final Vector3D point, final Line line) {
588         return recurseFirstIntersection(getTree(true), point, line);
589     }
590 
591     /** Get the first sub-hyperplane crossed by a semi-infinite line.
592      * @param node current node
593      * @param point start point of the part of the line considered
594      * @param line line to consider (contains point)
595      * @return the first sub-hyperplane crossed by the line after the
596      * given point, or null if the line does not intersect any
597      * sub-hyperplane
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         // establish search order
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             // search in the cut hyperplane
626             final SubHyperplane<Euclidean3D> facet = boundaryFacet(point, node);
627             if (facet != null) {
628                 return facet;
629             }
630         }
631 
632         // search in the near branch
633         final SubHyperplane<Euclidean3D> crossed = recurseFirstIntersection(near, point, line);
634         if (crossed != null) {
635             return crossed;
636         }
637 
638         if (!in) {
639             // search in the cut hyperplane
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         // search in the far branch
650         return recurseFirstIntersection(far, point, line);
651 
652     }
653 
654     /** Check if a point belongs to the boundary part of a node.
655      * @param point point to check
656      * @param node node containing the boundary facet to check
657      * @return the boundary facet this points belongs to (or null if it
658      * does not belong to any boundary facet)
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     /** Rotate the region around the specified point.
678      * <p>The instance is not modified, a new instance is created.</p>
679      * @param center rotation center
680      * @param rotation vectorial rotation operator
681      * @return a new instance representing the rotated region
682      */
683     public PolyhedronsSet rotate(final Vector3D center, final Rotation rotation) {
684         return (PolyhedronsSet) applyTransform(new RotationTransform(center, rotation));
685     }
686 
687     /** 3D rotation as a Transform. */
688     private static class RotationTransform implements Transform<Euclidean3D, Euclidean2D> {
689 
690         /** Center point of the rotation. */
691         private final Vector3D   center;
692 
693         /** Vectorial rotation. */
694         private final Rotation   rotation;
695 
696         /** Cached original hyperplane. */
697         private Plane cachedOriginal;
698 
699         /** Cached 2D transform valid inside the cached original hyperplane. */
700         private Transform<Euclidean2D, Euclidean1D>  cachedTransform;
701 
702         /** Build a rotation transform.
703          * @param center center point of the rotation
704          * @param rotation vectorial rotation
705          */
706         RotationTransform(final Vector3D center, final Rotation rotation) {
707             this.center   = center;
708             this.rotation = rotation;
709         }
710 
711         /** {@inheritDoc} */
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         /** {@inheritDoc} */
719         @Override
720         public Plane apply(final Hyperplane<Euclidean3D> hyperplane) {
721             return ((Plane) hyperplane).rotate(center, rotation);
722         }
723 
724         /** {@inheritDoc} */
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                 // we have changed hyperplane, reset the in-hyperplane transform
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     /** Translate the region by the specified amount.
757      * <p>The instance is not modified, a new instance is created.</p>
758      * @param translation translation to apply
759      * @return a new instance representing the translated region
760      */
761     public PolyhedronsSet translate(final Vector3D translation) {
762         return (PolyhedronsSet) applyTransform(new TranslationTransform(translation));
763     }
764 
765     /** 3D translation as a transform. */
766     private static class TranslationTransform implements Transform<Euclidean3D, Euclidean2D> {
767 
768         /** Translation vector. */
769         private final Vector3D   translation;
770 
771         /** Cached original hyperplane. */
772         private Plane cachedOriginal;
773 
774         /** Cached 2D transform valid inside the cached original hyperplane. */
775         private Transform<Euclidean2D, Euclidean1D>  cachedTransform;
776 
777         /** Build a translation transform.
778          * @param translation translation vector
779          */
780         TranslationTransform(final Vector3D translation) {
781             this.translation = translation;
782         }
783 
784         /** {@inheritDoc} */
785         @Override
786         public Vector3D apply(final Point<Euclidean3D> point) {
787             return new Vector3D(1.0, (Vector3D) point, 1.0, translation);
788         }
789 
790         /** {@inheritDoc} */
791         @Override
792         public Plane apply(final Hyperplane<Euclidean3D> hyperplane) {
793             return ((Plane) hyperplane).translate(translation);
794         }
795 
796         /** {@inheritDoc} */
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                 // we have changed hyperplane, reset the in-hyperplane transform
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     /** Container for Boundary REPresentation (B-Rep).
823      * <p>
824      * The boundary is provided as a list of vertices and a list of facets.
825      * Each facet is specified as an integer array containing the arrays vertices
826      * indices in the vertices list. Each facet normal is oriented by right hand
827      * rule to the facet vertices list.
828      * </p>
829      * @see PolyhedronsSet#PolyhedronsSet(BSPTree, double)
830      * @see PolyhedronsSet#getBRep()
831      * @since 1.2
832      */
833     public static class BRep {
834 
835         /** List of polyhedrons set vertices. */
836         private final List<Vector3D> vertices;
837 
838         /** List of facets, as vertices indices in the vertices list. */
839         private final List<int[]> facets;
840 
841         /** Simple constructor.
842          * @param vertices list of polyhedrons set vertices
843          * @param facets list of facets, as vertices indices in the vertices list
844          */
845         public BRep(final List<Vector3D> vertices, final List<int[]> facets) {
846             this.vertices = vertices;
847             this.facets   = facets;
848         }
849 
850         /** Get the extracted vertices.
851          * @return extracted vertices
852          */
853         public List<Vector3D> getVertices() {
854             return vertices;
855         }
856 
857         /** Get the extracted facets.
858          * @return extracted facets
859          */
860         public List<int[]> getFacets() {
861             return facets;
862         }
863 
864     }
865 
866 }