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.Vector;
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 1416643 2012-12-03 19:37:14Z tn $
040 * @since 3.0
041 */
042public 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}