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 * @version $Id: PolyhedronsSet.java 1590560 2014-04-28 06:39:01Z luc $
040 * @since 3.0
041 */
042public class PolyhedronsSet extends AbstractRegion<Euclidean3D, Euclidean2D> {
043
044    /** Default value for tolerance. */
045    private static final double DEFAULT_TOLERANCE = 1.0e-10;
046
047    /** Build a polyhedrons set representing the whole real line.
048     * @param tolerance tolerance below which points are considered identical
049     * @since 3.3
050     */
051    public PolyhedronsSet(final double tolerance) {
052        super(tolerance);
053    }
054
055    /** Build a polyhedrons set from a BSP tree.
056     * <p>The leaf nodes of the BSP tree <em>must</em> have a
057     * {@code Boolean} attribute representing the inside status of
058     * the corresponding cell (true for inside cells, false for outside
059     * cells). In order to avoid building too many small objects, it is
060     * recommended to use the predefined constants
061     * {@code Boolean.TRUE} and {@code Boolean.FALSE}</p>
062     * <p>
063     * This constructor is aimed at expert use, as building the tree may
064     * be a difficult task. It is not intended for general use and for
065     * performances reasons does not check thoroughly its input, as this would
066     * require walking the full tree each time. Failing to provide a tree with
067     * the proper attributes, <em>will</em> therefore generate problems like
068     * {@link NullPointerException} or {@link ClassCastException} only later on.
069     * This limitation is known and explains why this constructor is for expert
070     * use only. The caller does have the responsibility to provided correct arguments.
071     * </p>
072     * @param tree inside/outside BSP tree representing the region
073     * @param tolerance tolerance below which points are considered identical
074     * @since 3.3
075     */
076    public PolyhedronsSet(final BSPTree<Euclidean3D> tree, final double tolerance) {
077        super(tree, tolerance);
078    }
079
080    /** Build a polyhedrons set from a Boundary REPresentation (B-rep).
081     * <p>The boundary is provided as a collection of {@link
082     * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the
083     * interior part of the region on its minus side and the exterior on
084     * its plus side.</p>
085     * <p>The boundary elements can be in any order, and can form
086     * several non-connected sets (like for example polyhedrons with holes
087     * or a set of disjoint polyhedrons considered as a whole). In
088     * fact, the elements do not even need to be connected together
089     * (their topological connections are not used here). However, if the
090     * boundary does not really separate an inside open from an outside
091     * open (open having here its topological meaning), then subsequent
092     * calls to the {@link Region#checkPoint(Point) checkPoint} method will
093     * not be meaningful anymore.</p>
094     * <p>If the boundary is empty, the region will represent the whole
095     * space.</p>
096     * @param boundary collection of boundary elements, as a
097     * collection of {@link SubHyperplane SubHyperplane} objects
098     * @param tolerance tolerance below which points are considered identical
099     * @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}