    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     */
    017    package org.apache.commons.math3.geometry.euclidean.threed;
    018
    019    import java.awt.geom.AffineTransform;
    020    import java.util.Collection;
    021
    022    import org.apache.commons.math3.geometry.Vector;
    023    import org.apache.commons.math3.geometry.euclidean.oned.Euclidean1D;
    024    import org.apache.commons.math3.geometry.euclidean.twod.Euclidean2D;
    025    import org.apache.commons.math3.geometry.euclidean.twod.SubLine;
    026    import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
    027    import org.apache.commons.math3.geometry.partitioning.AbstractRegion;
    028    import org.apache.commons.math3.geometry.partitioning.BSPTree;
    029    import org.apache.commons.math3.geometry.partitioning.BSPTreeVisitor;
    030    import org.apache.commons.math3.geometry.partitioning.BoundaryAttribute;
    031    import org.apache.commons.math3.geometry.partitioning.Hyperplane;
    032    import org.apache.commons.math3.geometry.partitioning.Region;
    033    import org.apache.commons.math3.geometry.partitioning.RegionFactory;
    034    import org.apache.commons.math3.geometry.partitioning.SubHyperplane;
    035    import org.apache.commons.math3.geometry.partitioning.Transform;
    036    import org.apache.commons.math3.util.FastMath;
    037
    038    /** This class represents a 3D region: a set of polyhedrons.
    039     * @version $Id: PolyhedronsSet.java 1416643 2012-12-03 19:37:14Z tn $
    040     * @since 3.0
    041     */
    042    public class PolyhedronsSet extends AbstractRegion<Euclidean3D, Euclidean2D> {
    043
    044        /** Build a polyhedrons set representing the whole real line.
    045         */
    046        public PolyhedronsSet() {
    047            super();
    048        }
    049
    050        /** Build a polyhedrons set from a BSP tree.
    051         * <p>The leaf nodes of the BSP tree <em>must</em> have a
    052         * {@code Boolean} attribute representing the inside status of
    053         * the corresponding cell (true for inside cells, false for outside
    054         * cells). In order to avoid building too many small objects, it is
    055         * recommended to use the predefined constants
    056         * {@code Boolean.TRUE} and {@code Boolean.FALSE}</p>
    057         * @param tree inside/outside BSP tree representing the region
    058         */
    059        public PolyhedronsSet(final BSPTree<Euclidean3D> tree) {
    060            super(tree);
    061        }
    062
    063        /** Build a polyhedrons set from a Boundary REPresentation (B-rep).
    064         * <p>The boundary is provided as a collection of {@link
    065         * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the
    066         * interior part of the region on its minus side and the exterior on
    067         * its plus side.</p>
    068         * <p>The boundary elements can be in any order, and can form
    069         * several non-connected sets (like for example polyhedrons with holes
    070         * or a set of disjoint polyhedrons considered as a whole). In
    071         * fact, the elements do not even need to be connected together
    072         * (their topological connections are not used here). However, if the
    073         * boundary does not really separate an inside open from an outside
    074         * open (open having here its topological meaning), then subsequent
    075         * calls to the {@link Region#checkPoint(Vector) checkPoint} method will
    076         * not be meaningful anymore.</p>
    077         * <p>If the boundary is empty, the region will represent the whole
    078         * space.</p>
    079         * @param boundary collection of boundary elements, as a
    080         * collection of {@link SubHyperplane SubHyperplane} objects
    081         */
    082        public PolyhedronsSet(final Collection<SubHyperplane<Euclidean3D>> boundary) {
    083            super(boundary);
    084        }
    085
    086        /** Build a parallellepipedic box.
    087         * @param xMin low bound along the x direction
    088         * @param xMax high bound along the x direction
    089         * @param yMin low bound along the y direction
    090         * @param yMax high bound along the y direction
    091         * @param zMin low bound along the z direction
    092         * @param zMax high bound along the z direction
    093         */
    094        public PolyhedronsSet(final double xMin, final double xMax,
    095                              final double yMin, final double yMax,
    096                              final double zMin, final double zMax) {
    097            super(buildBoundary(xMin, xMax, yMin, yMax, zMin, zMax));
    098        }
    099
    100        /** Build a parallellepipedic box boundary.
    101         * @param xMin low bound along the x direction
    102         * @param xMax high bound along the x direction
    103         * @param yMin low bound along the y direction
    104         * @param yMax high bound along the y direction
    105         * @param zMin low bound along the z direction
    106         * @param zMax high bound along the z direction
    107         * @return boundary tree
    108         */
    109        private static BSPTree<Euclidean3D> buildBoundary(final double xMin, final double xMax,
    110                                                          final double yMin, final double yMax,
    111                                                          final double zMin, final double zMax) {
    112            final Plane pxMin = new Plane(new Vector3D(xMin, 0,    0),   Vector3D.MINUS_I);
    113            final Plane pxMax = new Plane(new Vector3D(xMax, 0,    0),   Vector3D.PLUS_I);
    114            final Plane pyMin = new Plane(new Vector3D(0,    yMin, 0),   Vector3D.MINUS_J);
    115            final Plane pyMax = new Plane(new Vector3D(0,    yMax, 0),   Vector3D.PLUS_J);
    116            final Plane pzMin = new Plane(new Vector3D(0,    0,   zMin), Vector3D.MINUS_K);
    117            final Plane pzMax = new Plane(new Vector3D(0,    0,   zMax), Vector3D.PLUS_K);
    118            @SuppressWarnings("unchecked")
    119            final Region<Euclidean3D> boundary =
    120            new RegionFactory<Euclidean3D>().buildConvex(pxMin, pxMax, pyMin, pyMax, pzMin, pzMax);
    121            return boundary.getTree(false);
    122        }
    123
    124        /** {@inheritDoc} */
    125        @Override
    126        public PolyhedronsSet buildNew(final BSPTree<Euclidean3D> tree) {
    127            return new PolyhedronsSet(tree);
    128        }
    129
    130        /** {@inheritDoc} */
    131        @Override
    132        protected void computeGeometricalProperties() {
    133
    134            // compute the contribution of all boundary facets
    135            getTree(true).visit(new FacetsContributionVisitor());
    136
    137            if (getSize() < 0) {
    138                // the polyhedrons set as a finite outside
    139                // surrounded by an infinite inside
    140                setSize(Double.POSITIVE_INFINITY);
    141                setBarycenter(Vector3D.NaN);
    142            } else {
    143                // the polyhedrons set is finite, apply the remaining scaling factors
    144                setSize(getSize() / 3.0);
    145                setBarycenter(new Vector3D(1.0 / (4 * getSize()), (Vector3D) getBarycenter()));
    146            }
    147
    148        }
    149
    150        /** Visitor computing geometrical properties. */
    151        private class FacetsContributionVisitor implements BSPTreeVisitor<Euclidean3D> {
    152
    153            /** Simple constructor. */
    154            public FacetsContributionVisitor() {
    155                setSize(0);
    156                setBarycenter(new Vector3D(0, 0, 0));
    157            }
    158
    159            /** {@inheritDoc} */
    160            public Order visitOrder(final BSPTree<Euclidean3D> node) {
    161                return Order.MINUS_SUB_PLUS;
    162            }
    163
    164            /** {@inheritDoc} */
    165            public void visitInternalNode(final BSPTree<Euclidean3D> node) {
    166                @SuppressWarnings("unchecked")
    167                final BoundaryAttribute<Euclidean3D> attribute =
    168                    (BoundaryAttribute<Euclidean3D>) node.getAttribute();
    169                if (attribute.getPlusOutside() != null) {
    170                    addContribution(attribute.getPlusOutside(), false);
    171                }
    172                if (attribute.getPlusInside() != null) {
    173                    addContribution(attribute.getPlusInside(), true);
    174                }
    175            }
    176
    177            /** {@inheritDoc} */
    178            public void visitLeafNode(final BSPTree<Euclidean3D> node) {
    179            }
    180
    181            /** Add he contribution of a boundary facet.
    182             * @param facet boundary facet
    183             * @param reversed if true, the facet has the inside on its plus side
    184             */
    185            private void addContribution(final SubHyperplane<Euclidean3D> facet, final boolean reversed) {
    186
    187                final Region<Euclidean2D> polygon = ((SubPlane) facet).getRemainingRegion();
    188                final double area    = polygon.getSize();
    189
    190                if (Double.isInfinite(area)) {
    191                    setSize(Double.POSITIVE_INFINITY);
    192                    setBarycenter(Vector3D.NaN);
    193                } else {
    194
    195                    final Plane    plane  = (Plane) facet.getHyperplane();
    196                    final Vector3D facetB = plane.toSpace(polygon.getBarycenter());
    197                    double   scaled = area * facetB.dotProduct(plane.getNormal());
    198                    if (reversed) {
    199                        scaled = -scaled;
    200                    }
    201
    202                    setSize(getSize() + scaled);
    203                    setBarycenter(new Vector3D(1.0, (Vector3D) getBarycenter(), scaled, facetB));
    204
    205                }
    206
    207            }
    208
    209        }
    210
    211        /** Get the first sub-hyperplane crossed by a semi-infinite line.
    212         * @param point start point of the part of the line considered
    213         * @param line line to consider (contains point)
    214         * @return the first sub-hyperplaned crossed by the line after the
    215         * given point, or null if the line does not intersect any
    216         * sub-hyperplaned
    217         */
    218        public SubHyperplane<Euclidean3D> firstIntersection(final Vector3D point, final Line line) {
    219            return recurseFirstIntersection(getTree(true), point, line);
    220        }
    221
    222        /** Get the first sub-hyperplane crossed by a semi-infinite line.
    223         * @param node current node
    224         * @param point start point of the part of the line considered
    225         * @param line line to consider (contains point)
    226         * @return the first sub-hyperplaned crossed by the line after the
    227         * given point, or null if the line does not intersect any
    228         * sub-hyperplaned
    229         */
    230        private SubHyperplane<Euclidean3D> recurseFirstIntersection(final BSPTree<Euclidean3D> node,
    231                                                                    final Vector3D point,
    232                                                                    final Line line) {
    233
    234            final SubHyperplane<Euclidean3D> cut = node.getCut();
    235            if (cut == null) {
    236                return null;
    237            }
    238            final BSPTree<Euclidean3D> minus = node.getMinus();
    239            final BSPTree<Euclidean3D> plus  = node.getPlus();
    240            final Plane               plane = (Plane) cut.getHyperplane();
    241
    242            // establish search order
    243            final double offset = plane.getOffset(point);
    244            final boolean in    = FastMath.abs(offset) < 1.0e-10;
    245            final BSPTree<Euclidean3D> near;
    246            final BSPTree<Euclidean3D> far;
    247            if (offset < 0) {
    248                near = minus;
    249                far  = plus;
    250            } else {
    251                near = plus;
    252                far  = minus;
    253            }
    254
    255            if (in) {
    256                // search in the cut hyperplane
    257                final SubHyperplane<Euclidean3D> facet = boundaryFacet(point, node);
    258                if (facet != null) {
    259                    return facet;
    260                }
    261            }
    262
    263            // search in the near branch
    264            final SubHyperplane<Euclidean3D> crossed = recurseFirstIntersection(near, point, line);
    265            if (crossed != null) {
    266                return crossed;
    267            }
    268
    269            if (!in) {
    270                // search in the cut hyperplane
    271                final Vector3D hit3D = plane.intersection(line);
    272                if (hit3D != null) {
    273                    final SubHyperplane<Euclidean3D> facet = boundaryFacet(hit3D, node);
    274                    if (facet != null) {
    275                        return facet;
    276                    }
    277                }
    278            }
    279
    280            // search in the far branch
    281            return recurseFirstIntersection(far, point, line);
    282
    283        }
    284
    285        /** Check if a point belongs to the boundary part of a node.
    286         * @param point point to check
    287         * @param node node containing the boundary facet to check
    288         * @return the boundary facet this points belongs to (or null if it
    289         * does not belong to any boundary facet)
    290         */
    291        private SubHyperplane<Euclidean3D> boundaryFacet(final Vector3D point,
    292                                                         final BSPTree<Euclidean3D> node) {
    293            final Vector2D point2D = ((Plane) node.getCut().getHyperplane()).toSubSpace(point);
    294            @SuppressWarnings("unchecked")
    295            final BoundaryAttribute<Euclidean3D> attribute =
    296                (BoundaryAttribute<Euclidean3D>) node.getAttribute();
    297            if ((attribute.getPlusOutside() != null) &&
    298                (((SubPlane) attribute.getPlusOutside()).getRemainingRegion().checkPoint(point2D) == Location.INSIDE)) {
    299                return attribute.getPlusOutside();
    300            }
    301            if ((attribute.getPlusInside() != null) &&
    302                (((SubPlane) attribute.getPlusInside()).getRemainingRegion().checkPoint(point2D) == Location.INSIDE)) {
    303                return attribute.getPlusInside();
    304            }
    305            return null;
    306        }
    307
    308        /** Rotate the region around the specified point.
    309         * <p>The instance is not modified, a new instance is created.</p>
    310         * @param center rotation center
    311         * @param rotation vectorial rotation operator
    312         * @return a new instance representing the rotated region
    313         */
    314        public PolyhedronsSet rotate(final Vector3D center, final Rotation rotation) {
    315            return (PolyhedronsSet) applyTransform(new RotationTransform(center, rotation));
    316        }
    317
    318        /** 3D rotation as a Transform. */
    319        private static class RotationTransform implements Transform<Euclidean3D, Euclidean2D> {
    320
    321            /** Center point of the rotation. */
    322            private Vector3D   center;
    323
    324            /** Vectorial rotation. */
    325            private Rotation   rotation;
    326
    327            /** Cached original hyperplane. */
    328            private Plane cachedOriginal;
    329
    330            /** Cached 2D transform valid inside the cached original hyperplane. */
    331            private Transform<Euclidean2D, Euclidean1D>  cachedTransform;
    332
    333            /** Build a rotation transform.
    334             * @param center center point of the rotation
    335             * @param rotation vectorial rotation
    336             */
    337            public RotationTransform(final Vector3D center, final Rotation rotation) {
    338                this.center   = center;
    339                this.rotation = rotation;
    340            }
    341
    342            /** {@inheritDoc} */
    343            public Vector3D apply(final Vector<Euclidean3D> point) {
    344                final Vector3D delta = ((Vector3D) point).subtract(center);
    345                return new Vector3D(1.0, center, 1.0, rotation.applyTo(delta));
    346            }
    347
    348            /** {@inheritDoc} */
    349            public Plane apply(final Hyperplane<Euclidean3D> hyperplane) {
    350                return ((Plane) hyperplane).rotate(center, rotation);
    351            }
    352
    353            /** {@inheritDoc} */
    354            public SubHyperplane<Euclidean2D> apply(final SubHyperplane<Euclidean2D> sub,
    355                                                    final Hyperplane<Euclidean3D> original,
    356                                                    final Hyperplane<Euclidean3D> transformed) {
    357                if (original != cachedOriginal) {
    358                    // we have changed hyperplane, reset the in-hyperplane transform
    359
    360                    final Plane    oPlane = (Plane) original;
    361                    final Plane    tPlane = (Plane) transformed;
    362                    final Vector3D p00    = oPlane.getOrigin();
    363                    final Vector3D p10    = oPlane.toSpace(new Vector2D(1.0, 0.0));
    364                    final Vector3D p01    = oPlane.toSpace(new Vector2D(0.0, 1.0));
    365                    final Vector2D tP00   = tPlane.toSubSpace(apply(p00));
    366                    final Vector2D tP10   = tPlane.toSubSpace(apply(p10));
    367                    final Vector2D tP01   = tPlane.toSubSpace(apply(p01));
    368                    final AffineTransform at =
    369                        new AffineTransform(tP10.getX() - tP00.getX(), tP10.getY() - tP00.getY(),
    370                                            tP01.getX() - tP00.getX(), tP01.getY() - tP00.getY(),
    371                                            tP00.getX(), tP00.getY());
    372
    373                    cachedOriginal  = (Plane) original;
    374                    cachedTransform = org.apache.commons.math3.geometry.euclidean.twod.Line.getTransform(at);
    375
    376                }
    377                return ((SubLine) sub).applyTransform(cachedTransform);
    378            }
    379
    380        }
    381
    382        /** Translate the region by the specified amount.
    383         * <p>The instance is not modified, a new instance is created.</p>
    384         * @param translation translation to apply
    385         * @return a new instance representing the translated region
    386         */
    387        public PolyhedronsSet translate(final Vector3D translation) {
    388            return (PolyhedronsSet) applyTransform(new TranslationTransform(translation));
    389        }
    390
    391        /** 3D translation as a transform. */
    392        private static class TranslationTransform implements Transform<Euclidean3D, Euclidean2D> {
    393
    394            /** Translation vector. */
    395            private Vector3D   translation;
    396
    397            /** Cached original hyperplane. */
    398            private Plane cachedOriginal;
    399
    400            /** Cached 2D transform valid inside the cached original hyperplane. */
    401            private Transform<Euclidean2D, Euclidean1D>  cachedTransform;
    402
    403            /** Build a translation transform.
    404             * @param translation translation vector
    405             */
    406            public TranslationTransform(final Vector3D translation) {
    407                this.translation = translation;
    408            }
    409
    410            /** {@inheritDoc} */
    411            public Vector3D apply(final Vector<Euclidean3D> point) {
    412                return new Vector3D(1.0, (Vector3D) point, 1.0, translation);
    413            }
    414
    415            /** {@inheritDoc} */
    416            public Plane apply(final Hyperplane<Euclidean3D> hyperplane) {
    417                return ((Plane) hyperplane).translate(translation);
    418            }
    419
    420            /** {@inheritDoc} */
    421            public SubHyperplane<Euclidean2D> apply(final SubHyperplane<Euclidean2D> sub,
    422                                                    final Hyperplane<Euclidean3D> original,
    423                                                    final Hyperplane<Euclidean3D> transformed) {
    424                if (original != cachedOriginal) {
    425                    // we have changed hyperplane, reset the in-hyperplane transform
    426
    427                    final Plane   oPlane = (Plane) original;
    428                    final Plane   tPlane = (Plane) transformed;
    429                    final Vector2D shift  = tPlane.toSubSpace(apply(oPlane.getOrigin()));
    430                    final AffineTransform at =
    431                        AffineTransform.getTranslateInstance(shift.getX(), shift.getY());
    432
    433                    cachedOriginal  = (Plane) original;
    434                    cachedTransform =
    435                            org.apache.commons.math3.geometry.euclidean.twod.Line.getTransform(at);
    436
    437                }
    438
    439                return ((SubLine) sub).applyTransform(cachedTransform);
    440
    441            }
    442
    443        }
    444
    445    }