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