001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math3.geometry.euclidean.threed;
018
019import java.awt.geom.AffineTransform;
020import java.util.Collection;
021
022import org.apache.commons.math3.geometry.Point;
023import org.apache.commons.math3.geometry.euclidean.oned.Euclidean1D;
024import org.apache.commons.math3.geometry.euclidean.twod.Euclidean2D;
025import org.apache.commons.math3.geometry.euclidean.twod.SubLine;
026import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
027import org.apache.commons.math3.geometry.partitioning.AbstractRegion;
028import org.apache.commons.math3.geometry.partitioning.BSPTree;
029import org.apache.commons.math3.geometry.partitioning.BSPTreeVisitor;
030import org.apache.commons.math3.geometry.partitioning.BoundaryAttribute;
031import org.apache.commons.math3.geometry.partitioning.Hyperplane;
032import org.apache.commons.math3.geometry.partitioning.Region;
033import org.apache.commons.math3.geometry.partitioning.RegionFactory;
034import org.apache.commons.math3.geometry.partitioning.SubHyperplane;
035import org.apache.commons.math3.geometry.partitioning.Transform;
036import org.apache.commons.math3.util.FastMath;
037
038/** This class represents a 3D region: a set of polyhedrons.
039 * @since 3.0
040 */
041public class PolyhedronsSet extends AbstractRegion<Euclidean3D, Euclidean2D> {
042
043    /** Default value for tolerance. */
044    private static final double DEFAULT_TOLERANCE = 1.0e-10;
045
046    /** Build a polyhedrons set representing the whole real line.
047     * @param tolerance tolerance below which points are considered identical
048     * @since 3.3
049     */
050    public PolyhedronsSet(final double tolerance) {
051        super(tolerance);
052    }
053
054    /** Build a polyhedrons set from a BSP tree.
055     * <p>The leaf nodes of the BSP tree <em>must</em> have a
056     * {@code Boolean} attribute representing the inside status of
057     * the corresponding cell (true for inside cells, false for outside
058     * cells). In order to avoid building too many small objects, it is
059     * recommended to use the predefined constants
060     * {@code Boolean.TRUE} and {@code Boolean.FALSE}</p>
061     * <p>
062     * This constructor is aimed at expert use, as building the tree may
063     * be a difficult task. It is not intended for general use and for
064     * performances reasons does not check thoroughly its input, as this would
065     * require walking the full tree each time. Failing to provide a tree with
066     * the proper attributes, <em>will</em> therefore generate problems like
067     * {@link NullPointerException} or {@link ClassCastException} only later on.
068     * This limitation is known and explains why this constructor is for expert
069     * use only. The caller does have the responsibility to provided correct arguments.
070     * </p>
071     * @param tree inside/outside BSP tree representing the region
072     * @param tolerance tolerance below which points are considered identical
073     * @since 3.3
074     */
075    public PolyhedronsSet(final BSPTree<Euclidean3D> tree, final double tolerance) {
076        super(tree, tolerance);
077    }
078
079    /** Build a polyhedrons set from a Boundary REPresentation (B-rep).
080     * <p>The boundary is provided as a collection of {@link
081     * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the
082     * interior part of the region on its minus side and the exterior on
083     * its plus side.</p>
084     * <p>The boundary elements can be in any order, and can form
085     * several non-connected sets (like for example polyhedrons with holes
086     * or a set of disjoint polyhedrons considered as a whole). In
087     * fact, the elements do not even need to be connected together
088     * (their topological connections are not used here). However, if the
089     * boundary does not really separate an inside open from an outside
090     * open (open having here its topological meaning), then subsequent
091     * calls to the {@link Region#checkPoint(Point) checkPoint} method will
092     * not be meaningful anymore.</p>
093     * <p>If the boundary is empty, the region will represent the whole
094     * space.</p>
095     * @param boundary collection of boundary elements, as a
096     * collection of {@link SubHyperplane SubHyperplane} objects
097     * @param tolerance tolerance below which points are considered identical
098     * @since 3.3
099     */
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}