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    *      http://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  package org.apache.commons.math3.geometry.euclidean.threed;
18  
19  import java.awt.geom.AffineTransform;
20  import java.util.Collection;
21  
22  import org.apache.commons.math3.geometry.Point;
23  import org.apache.commons.math3.geometry.euclidean.oned.Euclidean1D;
24  import org.apache.commons.math3.geometry.euclidean.twod.Euclidean2D;
25  import org.apache.commons.math3.geometry.euclidean.twod.SubLine;
26  import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
27  import org.apache.commons.math3.geometry.partitioning.AbstractRegion;
28  import org.apache.commons.math3.geometry.partitioning.BSPTree;
29  import org.apache.commons.math3.geometry.partitioning.BSPTreeVisitor;
30  import org.apache.commons.math3.geometry.partitioning.BoundaryAttribute;
31  import org.apache.commons.math3.geometry.partitioning.Hyperplane;
32  import org.apache.commons.math3.geometry.partitioning.Region;
33  import org.apache.commons.math3.geometry.partitioning.RegionFactory;
34  import org.apache.commons.math3.geometry.partitioning.SubHyperplane;
35  import org.apache.commons.math3.geometry.partitioning.Transform;
36  import org.apache.commons.math3.util.FastMath;
37  
38  /** This class represents a 3D region: a set of polyhedrons.
39   * @since 3.0
40   */
41  public class PolyhedronsSet extends AbstractRegion<Euclidean3D, Euclidean2D> {
42  
43      /** Default value for tolerance. */
44      private static final double DEFAULT_TOLERANCE = 1.0e-10;
45  
46      /** Build a polyhedrons set representing the whole real line.
47       * @param tolerance tolerance below which points are considered identical
48       * @since 3.3
49       */
50      public PolyhedronsSet(final double tolerance) {
51          super(tolerance);
52      }
53  
54      /** Build a polyhedrons set from a BSP tree.
55       * <p>The leaf nodes of the BSP tree <em>must</em> have a
56       * {@code Boolean} attribute representing the inside status of
57       * the corresponding cell (true for inside cells, false for outside
58       * cells). In order to avoid building too many small objects, it is
59       * recommended to use the predefined constants
60       * {@code Boolean.TRUE} and {@code Boolean.FALSE}</p>
61       * <p>
62       * This constructor is aimed at expert use, as building the tree may
63       * be a difficult task. It is not intended for general use and for
64       * performances reasons does not check thoroughly its input, as this would
65       * require walking the full tree each time. Failing to provide a tree with
66       * the proper attributes, <em>will</em> therefore generate problems like
67       * {@link NullPointerException} or {@link ClassCastException} only later on.
68       * This limitation is known and explains why this constructor is for expert
69       * use only. The caller does have the responsibility to provided correct arguments.
70       * </p>
71       * @param tree inside/outside BSP tree representing the region
72       * @param tolerance tolerance below which points are considered identical
73       * @since 3.3
74       */
75      public PolyhedronsSet(final BSPTree<Euclidean3D> tree, final double tolerance) {
76          super(tree, tolerance);
77      }
78  
79      /** Build a polyhedrons set from a Boundary REPresentation (B-rep).
80       * <p>The boundary is provided as a collection of {@link
81       * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the
82       * interior part of the region on its minus side and the exterior on
83       * its plus side.</p>
84       * <p>The boundary elements can be in any order, and can form
85       * several non-connected sets (like for example polyhedrons with holes
86       * or a set of disjoint polyhedrons considered as a whole). In
87       * fact, the elements do not even need to be connected together
88       * (their topological connections are not used here). However, if the
89       * boundary does not really separate an inside open from an outside
90       * open (open having here its topological meaning), then subsequent
91       * calls to the {@link Region#checkPoint(Point) checkPoint} method will
92       * not be meaningful anymore.</p>
93       * <p>If the boundary is empty, the region will represent the whole
94       * space.</p>
95       * @param boundary collection of boundary elements, as a
96       * collection of {@link SubHyperplane SubHyperplane} objects
97       * @param tolerance tolerance below which points are considered identical
98       * @since 3.3
99       */
100     public PolyhedronsSet(final Collection<SubHyperplane<Euclidean3D>> boundary,
101                           final double tolerance) {
102         super(boundary, tolerance);
103     }
104 
105     /** Build a parallellepipedic box.
106      * @param xMin low bound along the x direction
107      * @param xMax high bound along the x direction
108      * @param yMin low bound along the y direction
109      * @param yMax high bound along the y direction
110      * @param zMin low bound along the z direction
111      * @param zMax high bound along the z direction
112      * @param tolerance tolerance below which points are considered identical
113      * @since 3.3
114      */
115     public PolyhedronsSet(final double xMin, final double xMax,
116                           final double yMin, final double yMax,
117                           final double zMin, final double zMax,
118                           final double tolerance) {
119         super(buildBoundary(xMin, xMax, yMin, yMax, zMin, zMax, tolerance), tolerance);
120     }
121 
122     /** Build a polyhedrons set representing the whole real line.
123      * @deprecated as of 3.3, replaced with {@link #PolyhedronsSet(double)}
124      */
125     @Deprecated
126     public PolyhedronsSet() {
127         this(DEFAULT_TOLERANCE);
128     }
129 
130     /** Build a polyhedrons set from a BSP tree.
131      * <p>The leaf nodes of the BSP tree <em>must</em> have a
132      * {@code Boolean} attribute representing the inside status of
133      * the corresponding cell (true for inside cells, false for outside
134      * cells). In order to avoid building too many small objects, it is
135      * recommended to use the predefined constants
136      * {@code Boolean.TRUE} and {@code Boolean.FALSE}</p>
137      * @param tree inside/outside BSP tree representing the region
138      * @deprecated as of 3.3, replaced with {@link #PolyhedronsSet(BSPTree, double)}
139      */
140     @Deprecated
141     public PolyhedronsSet(final BSPTree<Euclidean3D> tree) {
142         this(tree, DEFAULT_TOLERANCE);
143     }
144 
145     /** Build a polyhedrons set from a Boundary REPresentation (B-rep).
146      * <p>The boundary is provided as a collection of {@link
147      * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the
148      * interior part of the region on its minus side and the exterior on
149      * its plus side.</p>
150      * <p>The boundary elements can be in any order, and can form
151      * several non-connected sets (like for example polyhedrons with holes
152      * or a set of disjoint polyhedrons considered as a whole). In
153      * fact, the elements do not even need to be connected together
154      * (their topological connections are not used here). However, if the
155      * boundary does not really separate an inside open from an outside
156      * open (open having here its topological meaning), then subsequent
157      * calls to the {@link Region#checkPoint(Point) checkPoint} method will
158      * not be meaningful anymore.</p>
159      * <p>If the boundary is empty, the region will represent the whole
160      * space.</p>
161      * @param boundary collection of boundary elements, as a
162      * collection of {@link SubHyperplane SubHyperplane} objects
163      * @deprecated as of 3.3, replaced with {@link #PolyhedronsSet(Collection, double)}
164      */
165     @Deprecated
166     public PolyhedronsSet(final Collection<SubHyperplane<Euclidean3D>> boundary) {
167         this(boundary, DEFAULT_TOLERANCE);
168     }
169 
170     /** Build a parallellepipedic box.
171      * @param xMin low bound along the x direction
172      * @param xMax high bound along the x direction
173      * @param yMin low bound along the y direction
174      * @param yMax high bound along the y direction
175      * @param zMin low bound along the z direction
176      * @param zMax high bound along the z direction
177      * @deprecated as of 3.3, replaced with {@link #PolyhedronsSet(double, double,
178      * double, double, double, double, double)}
179      */
180     @Deprecated
181     public PolyhedronsSet(final double xMin, final double xMax,
182                           final double yMin, final double yMax,
183                           final double zMin, final double zMax) {
184         this(xMin, xMax, yMin, yMax, zMin, zMax, DEFAULT_TOLERANCE);
185     }
186 
187     /** Build a parallellepipedic box boundary.
188      * @param xMin low bound along the x direction
189      * @param xMax high bound along the x direction
190      * @param yMin low bound along the y direction
191      * @param yMax high bound along the y direction
192      * @param zMin low bound along the z direction
193      * @param zMax high bound along the z direction
194      * @param tolerance tolerance below which points are considered identical
195      * @return boundary tree
196      * @since 3.3
197      */
198     private static BSPTree<Euclidean3D> buildBoundary(final double xMin, final double xMax,
199                                                       final double yMin, final double yMax,
200                                                       final double zMin, final double zMax,
201                                                       final double tolerance) {
202         if ((xMin >= xMax - tolerance) || (yMin >= yMax - tolerance) || (zMin >= zMax - tolerance)) {
203             // too thin box, build an empty polygons set
204             return new BSPTree<Euclidean3D>(Boolean.FALSE);
205         }
206         final Plane pxMin = new Plane(new Vector3D(xMin, 0,    0),   Vector3D.MINUS_I, tolerance);
207         final Plane pxMax = new Plane(new Vector3D(xMax, 0,    0),   Vector3D.PLUS_I,  tolerance);
208         final Plane pyMin = new Plane(new Vector3D(0,    yMin, 0),   Vector3D.MINUS_J, tolerance);
209         final Plane pyMax = new Plane(new Vector3D(0,    yMax, 0),   Vector3D.PLUS_J,  tolerance);
210         final Plane pzMin = new Plane(new Vector3D(0,    0,   zMin), Vector3D.MINUS_K, tolerance);
211         final Plane pzMax = new Plane(new Vector3D(0,    0,   zMax), Vector3D.PLUS_K,  tolerance);
212         @SuppressWarnings("unchecked")
213         final Region<Euclidean3D> boundary =
214         new RegionFactory<Euclidean3D>().buildConvex(pxMin, pxMax, pyMin, pyMax, pzMin, pzMax);
215         return boundary.getTree(false);
216     }
217 
218     /** {@inheritDoc} */
219     @Override
220     public PolyhedronsSet buildNew(final BSPTree<Euclidean3D> tree) {
221         return new PolyhedronsSet(tree, getTolerance());
222     }
223 
224     /** {@inheritDoc} */
225     @Override
226     protected void computeGeometricalProperties() {
227 
228         // compute the contribution of all boundary facets
229         getTree(true).visit(new FacetsContributionVisitor());
230 
231         if (getSize() < 0) {
232             // the polyhedrons set as a finite outside
233             // surrounded by an infinite inside
234             setSize(Double.POSITIVE_INFINITY);
235             setBarycenter((Point<Euclidean3D>) Vector3D.NaN);
236         } else {
237             // the polyhedrons set is finite, apply the remaining scaling factors
238             setSize(getSize() / 3.0);
239             setBarycenter((Point<Euclidean3D>) new Vector3D(1.0 / (4 * getSize()), (Vector3D) getBarycenter()));
240         }
241 
242     }
243 
244     /** Visitor computing geometrical properties. */
245     private class FacetsContributionVisitor implements BSPTreeVisitor<Euclidean3D> {
246 
247         /** Simple constructor. */
248         public FacetsContributionVisitor() {
249             setSize(0);
250             setBarycenter((Point<Euclidean3D>) new Vector3D(0, 0, 0));
251         }
252 
253         /** {@inheritDoc} */
254         public Order visitOrder(final BSPTree<Euclidean3D> node) {
255             return Order.MINUS_SUB_PLUS;
256         }
257 
258         /** {@inheritDoc} */
259         public void visitInternalNode(final BSPTree<Euclidean3D> node) {
260             @SuppressWarnings("unchecked")
261             final BoundaryAttribute<Euclidean3D> attribute =
262                 (BoundaryAttribute<Euclidean3D>) node.getAttribute();
263             if (attribute.getPlusOutside() != null) {
264                 addContribution(attribute.getPlusOutside(), false);
265             }
266             if (attribute.getPlusInside() != null) {
267                 addContribution(attribute.getPlusInside(), true);
268             }
269         }
270 
271         /** {@inheritDoc} */
272         public void visitLeafNode(final BSPTree<Euclidean3D> node) {
273         }
274 
275         /** Add he contribution of a boundary facet.
276          * @param facet boundary facet
277          * @param reversed if true, the facet has the inside on its plus side
278          */
279         private void addContribution(final SubHyperplane<Euclidean3D> facet, final boolean reversed) {
280 
281             final Region<Euclidean2D> polygon = ((SubPlane) facet).getRemainingRegion();
282             final double area    = polygon.getSize();
283 
284             if (Double.isInfinite(area)) {
285                 setSize(Double.POSITIVE_INFINITY);
286                 setBarycenter((Point<Euclidean3D>) Vector3D.NaN);
287             } else {
288 
289                 final Plane    plane  = (Plane) facet.getHyperplane();
290                 final Vector3D facetB = plane.toSpace(polygon.getBarycenter());
291                 double   scaled = area * facetB.dotProduct(plane.getNormal());
292                 if (reversed) {
293                     scaled = -scaled;
294                 }
295 
296                 setSize(getSize() + scaled);
297                 setBarycenter((Point<Euclidean3D>) new Vector3D(1.0, (Vector3D) getBarycenter(), scaled, facetB));
298 
299             }
300 
301         }
302 
303     }
304 
305     /** Get the first sub-hyperplane crossed by a semi-infinite line.
306      * @param point start point of the part of the line considered
307      * @param line line to consider (contains point)
308      * @return the first sub-hyperplaned crossed by the line after the
309      * given point, or null if the line does not intersect any
310      * sub-hyperplaned
311      */
312     public SubHyperplane<Euclidean3D> firstIntersection(final Vector3D point, final Line line) {
313         return recurseFirstIntersection(getTree(true), point, line);
314     }
315 
316     /** Get the first sub-hyperplane crossed by a semi-infinite line.
317      * @param node current node
318      * @param point start point of the part of the line considered
319      * @param line line to consider (contains point)
320      * @return the first sub-hyperplaned crossed by the line after the
321      * given point, or null if the line does not intersect any
322      * sub-hyperplaned
323      */
324     private SubHyperplane<Euclidean3D> recurseFirstIntersection(final BSPTree<Euclidean3D> node,
325                                                                 final Vector3D point,
326                                                                 final Line line) {
327 
328         final SubHyperplane<Euclidean3D> cut = node.getCut();
329         if (cut == null) {
330             return null;
331         }
332         final BSPTree<Euclidean3D> minus = node.getMinus();
333         final BSPTree<Euclidean3D> plus  = node.getPlus();
334         final Plane               plane = (Plane) cut.getHyperplane();
335 
336         // establish search order
337         final double offset = plane.getOffset((Point<Euclidean3D>) point);
338         final boolean in    = FastMath.abs(offset) < 1.0e-10;
339         final BSPTree<Euclidean3D> near;
340         final BSPTree<Euclidean3D> far;
341         if (offset < 0) {
342             near = minus;
343             far  = plus;
344         } else {
345             near = plus;
346             far  = minus;
347         }
348 
349         if (in) {
350             // search in the cut hyperplane
351             final SubHyperplane<Euclidean3D> facet = boundaryFacet(point, node);
352             if (facet != null) {
353                 return facet;
354             }
355         }
356 
357         // search in the near branch
358         final SubHyperplane<Euclidean3D> crossed = recurseFirstIntersection(near, point, line);
359         if (crossed != null) {
360             return crossed;
361         }
362 
363         if (!in) {
364             // search in the cut hyperplane
365             final Vector3D hit3D = plane.intersection(line);
366             if (hit3D != null) {
367                 final SubHyperplane<Euclidean3D> facet = boundaryFacet(hit3D, node);
368                 if (facet != null) {
369                     return facet;
370                 }
371             }
372         }
373 
374         // search in the far branch
375         return recurseFirstIntersection(far, point, line);
376 
377     }
378 
379     /** Check if a point belongs to the boundary part of a node.
380      * @param point point to check
381      * @param node node containing the boundary facet to check
382      * @return the boundary facet this points belongs to (or null if it
383      * does not belong to any boundary facet)
384      */
385     private SubHyperplane<Euclidean3D> boundaryFacet(final Vector3D point,
386                                                      final BSPTree<Euclidean3D> node) {
387         final Vector2D point2D = ((Plane) node.getCut().getHyperplane()).toSubSpace((Point<Euclidean3D>) point);
388         @SuppressWarnings("unchecked")
389         final BoundaryAttribute<Euclidean3D> attribute =
390             (BoundaryAttribute<Euclidean3D>) node.getAttribute();
391         if ((attribute.getPlusOutside() != null) &&
392             (((SubPlane) attribute.getPlusOutside()).getRemainingRegion().checkPoint(point2D) == Location.INSIDE)) {
393             return attribute.getPlusOutside();
394         }
395         if ((attribute.getPlusInside() != null) &&
396             (((SubPlane) attribute.getPlusInside()).getRemainingRegion().checkPoint(point2D) == Location.INSIDE)) {
397             return attribute.getPlusInside();
398         }
399         return null;
400     }
401 
402     /** Rotate the region around the specified point.
403      * <p>The instance is not modified, a new instance is created.</p>
404      * @param center rotation center
405      * @param rotation vectorial rotation operator
406      * @return a new instance representing the rotated region
407      */
408     public PolyhedronsSet rotate(final Vector3D center, final Rotation rotation) {
409         return (PolyhedronsSet) applyTransform(new RotationTransform(center, rotation));
410     }
411 
412     /** 3D rotation as a Transform. */
413     private static class RotationTransform implements Transform<Euclidean3D, Euclidean2D> {
414 
415         /** Center point of the rotation. */
416         private Vector3D   center;
417 
418         /** Vectorial rotation. */
419         private Rotation   rotation;
420 
421         /** Cached original hyperplane. */
422         private Plane cachedOriginal;
423 
424         /** Cached 2D transform valid inside the cached original hyperplane. */
425         private Transform<Euclidean2D, Euclidean1D>  cachedTransform;
426 
427         /** Build a rotation transform.
428          * @param center center point of the rotation
429          * @param rotation vectorial rotation
430          */
431         public RotationTransform(final Vector3D center, final Rotation rotation) {
432             this.center   = center;
433             this.rotation = rotation;
434         }
435 
436         /** {@inheritDoc} */
437         public Vector3D apply(final Point<Euclidean3D> point) {
438             final Vector3D delta = ((Vector3D) point).subtract(center);
439             return new Vector3D(1.0, center, 1.0, rotation.applyTo(delta));
440         }
441 
442         /** {@inheritDoc} */
443         public Plane apply(final Hyperplane<Euclidean3D> hyperplane) {
444             return ((Plane) hyperplane).rotate(center, rotation);
445         }
446 
447         /** {@inheritDoc} */
448         public SubHyperplane<Euclidean2D> apply(final SubHyperplane<Euclidean2D> sub,
449                                                 final Hyperplane<Euclidean3D> original,
450                                                 final Hyperplane<Euclidean3D> transformed) {
451             if (original != cachedOriginal) {
452                 // we have changed hyperplane, reset the in-hyperplane transform
453 
454                 final Plane    oPlane = (Plane) original;
455                 final Plane    tPlane = (Plane) transformed;
456                 final Vector3D p00    = oPlane.getOrigin();
457                 final Vector3D p10    = oPlane.toSpace((Point<Euclidean2D>) new Vector2D(1.0, 0.0));
458                 final Vector3D p01    = oPlane.toSpace((Point<Euclidean2D>) new Vector2D(0.0, 1.0));
459                 final Vector2D tP00   = tPlane.toSubSpace((Point<Euclidean3D>) apply(p00));
460                 final Vector2D tP10   = tPlane.toSubSpace((Point<Euclidean3D>) apply(p10));
461                 final Vector2D tP01   = tPlane.toSubSpace((Point<Euclidean3D>) apply(p01));
462                 final AffineTransform at =
463                     new AffineTransform(tP10.getX() - tP00.getX(), tP10.getY() - tP00.getY(),
464                                         tP01.getX() - tP00.getX(), tP01.getY() - tP00.getY(),
465                                         tP00.getX(), tP00.getY());
466 
467                 cachedOriginal  = (Plane) original;
468                 cachedTransform = org.apache.commons.math3.geometry.euclidean.twod.Line.getTransform(at);
469 
470             }
471             return ((SubLine) sub).applyTransform(cachedTransform);
472         }
473 
474     }
475 
476     /** Translate the region by the specified amount.
477      * <p>The instance is not modified, a new instance is created.</p>
478      * @param translation translation to apply
479      * @return a new instance representing the translated region
480      */
481     public PolyhedronsSet translate(final Vector3D translation) {
482         return (PolyhedronsSet) applyTransform(new TranslationTransform(translation));
483     }
484 
485     /** 3D translation as a transform. */
486     private static class TranslationTransform implements Transform<Euclidean3D, Euclidean2D> {
487 
488         /** Translation vector. */
489         private Vector3D   translation;
490 
491         /** Cached original hyperplane. */
492         private Plane cachedOriginal;
493 
494         /** Cached 2D transform valid inside the cached original hyperplane. */
495         private Transform<Euclidean2D, Euclidean1D>  cachedTransform;
496 
497         /** Build a translation transform.
498          * @param translation translation vector
499          */
500         public TranslationTransform(final Vector3D translation) {
501             this.translation = translation;
502         }
503 
504         /** {@inheritDoc} */
505         public Vector3D apply(final Point<Euclidean3D> point) {
506             return new Vector3D(1.0, (Vector3D) point, 1.0, translation);
507         }
508 
509         /** {@inheritDoc} */
510         public Plane apply(final Hyperplane<Euclidean3D> hyperplane) {
511             return ((Plane) hyperplane).translate(translation);
512         }
513 
514         /** {@inheritDoc} */
515         public SubHyperplane<Euclidean2D> apply(final SubHyperplane<Euclidean2D> sub,
516                                                 final Hyperplane<Euclidean3D> original,
517                                                 final Hyperplane<Euclidean3D> transformed) {
518             if (original != cachedOriginal) {
519                 // we have changed hyperplane, reset the in-hyperplane transform
520 
521                 final Plane   oPlane = (Plane) original;
522                 final Plane   tPlane = (Plane) transformed;
523                 final Vector2D shift  = tPlane.toSubSpace((Point<Euclidean3D>) apply(oPlane.getOrigin()));
524                 final AffineTransform at =
525                     AffineTransform.getTranslateInstance(shift.getX(), shift.getY());
526 
527                 cachedOriginal  = (Plane) original;
528                 cachedTransform =
529                         org.apache.commons.math3.geometry.euclidean.twod.Line.getTransform(at);
530 
531             }
532 
533             return ((SubLine) sub).applyTransform(cachedTransform);
534 
535         }
536 
537     }
538 
539 }