    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.util.ArrayList;
    020
    021    import org.apache.commons.math3.geometry.euclidean.twod.Euclidean2D;
    022    import org.apache.commons.math3.geometry.euclidean.twod.PolygonsSet;
    023    import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
    024    import org.apache.commons.math3.geometry.partitioning.AbstractSubHyperplane;
    025    import org.apache.commons.math3.geometry.partitioning.BSPTree;
    026    import org.apache.commons.math3.geometry.partitioning.BSPTreeVisitor;
    027    import org.apache.commons.math3.geometry.partitioning.BoundaryAttribute;
    028    import org.apache.commons.math3.geometry.partitioning.RegionFactory;
    029    import org.apache.commons.math3.geometry.partitioning.SubHyperplane;
    030    import org.apache.commons.math3.util.FastMath;
    031
    032    /** Extractor for {@link PolygonsSet polyhedrons sets} outlines.
    033     * <p>This class extracts the 2D outlines from {{@link PolygonsSet
    034     * polyhedrons sets} in a specified projection plane.</p>
    035     * @version $Id: OutlineExtractor.java 1416643 2012-12-03 19:37:14Z tn $
    036     * @since 3.0
    037     */
    038    public class OutlineExtractor {
    039
    040        /** Abscissa axis of the projection plane. */
    041        private Vector3D u;
    042
    043        /** Ordinate axis of the projection plane. */
    044        private Vector3D v;
    045
    046        /** Normal of the projection plane (viewing direction). */
    047        private Vector3D w;
    048
    049        /** Build an extractor for a specific projection plane.
    050         * @param u abscissa axis of the projection point
    051         * @param v ordinate axis of the projection point
    052         */
    053        public OutlineExtractor(final Vector3D u, final Vector3D v) {
    054            this.u = u;
    055            this.v = v;
    056            w = Vector3D.crossProduct(u, v);
    057        }
    058
    059        /** Extract the outline of a polyhedrons set.
    060         * @param polyhedronsSet polyhedrons set whose outline must be extracted
    061         * @return an outline, as an array of loops.
    062         */
    063        public Vector2D[][] getOutline(final PolyhedronsSet polyhedronsSet) {
    064
    065            // project all boundary facets into one polygons set
    066            final BoundaryProjector projector = new BoundaryProjector();
    067            polyhedronsSet.getTree(true).visit(projector);
    068            final PolygonsSet projected = projector.getProjected();
    069
    070            // Remove the spurious intermediate vertices from the outline
    071            final Vector2D[][] outline = projected.getVertices();
    072            for (int i = 0; i < outline.length; ++i) {
    073                final Vector2D[] rawLoop = outline[i];
    074                int end = rawLoop.length;
    075                int j = 0;
    076                while (j < end) {
    077                    if (pointIsBetween(rawLoop, end, j)) {
    078                        // the point should be removed
    079                        for (int k = j; k < (end - 1); ++k) {
    080                            rawLoop[k] = rawLoop[k + 1];
    081                        }
    082                        --end;
    083                    } else {
    084                        // the point remains in the loop
    085                        ++j;
    086                    }
    087                }
    088                if (end != rawLoop.length) {
    089                    // resize the array
    090                    outline[i] = new Vector2D[end];
    091                    System.arraycopy(rawLoop, 0, outline[i], 0, end);
    092                }
    093            }
    094
    095            return outline;
    096
    097        }
    098
    099        /** Check if a point is geometrically between its neighbour in an array.
    100         * <p>The neighbours are computed considering the array is a loop
    101         * (i.e. point at index (n-1) is before point at index 0)</p>
    102         * @param loop points array
    103         * @param n number of points to consider in the array
    104         * @param i index of the point to check (must be between 0 and n-1)
    105         * @return true if the point is exactly between its neighbours
    106         */
    107        private boolean pointIsBetween(final Vector2D[] loop, final int n, final int i) {
    108            final Vector2D previous = loop[(i + n - 1) % n];
    109            final Vector2D current  = loop[i];
    110            final Vector2D next     = loop[(i + 1) % n];
    111            final double dx1       = current.getX() - previous.getX();
    112            final double dy1       = current.getY() - previous.getY();
    113            final double dx2       = next.getX()    - current.getX();
    114            final double dy2       = next.getY()    - current.getY();
    115            final double cross     = dx1 * dy2 - dx2 * dy1;
    116            final double dot       = dx1 * dx2 + dy1 * dy2;
    117            final double d1d2      = FastMath.sqrt((dx1 * dx1 + dy1 * dy1) * (dx2 * dx2 + dy2 * dy2));
    118            return (FastMath.abs(cross) <= (1.0e-6 * d1d2)) && (dot >= 0.0);
    119        }
    120
    121        /** Visitor projecting the boundary facets on a plane. */
    122        private class BoundaryProjector implements BSPTreeVisitor<Euclidean3D> {
    123
    124            /** Projection of the polyhedrons set on the plane. */
    125            private PolygonsSet projected;
    126
    127            /** Simple constructor.
    128             */
    129            public BoundaryProjector() {
    130                projected = new PolygonsSet(new BSPTree<Euclidean2D>(Boolean.FALSE));
    131            }
    132
    133            /** {@inheritDoc} */
    134            public Order visitOrder(final BSPTree<Euclidean3D> node) {
    135                return Order.MINUS_SUB_PLUS;
    136            }
    137
    138            /** {@inheritDoc} */
    139            public void visitInternalNode(final BSPTree<Euclidean3D> node) {
    140                @SuppressWarnings("unchecked")
    141                final BoundaryAttribute<Euclidean3D> attribute =
    142                    (BoundaryAttribute<Euclidean3D>) node.getAttribute();
    143                if (attribute.getPlusOutside() != null) {
    144                    addContribution(attribute.getPlusOutside(), false);
    145                }
    146                if (attribute.getPlusInside() != null) {
    147                    addContribution(attribute.getPlusInside(), true);
    148                }
    149            }
    150
    151            /** {@inheritDoc} */
    152            public void visitLeafNode(final BSPTree<Euclidean3D> node) {
    153            }
    154
    155            /** Add he contribution of a boundary facet.
    156             * @param facet boundary facet
    157             * @param reversed if true, the facet has the inside on its plus side
    158             */
    159            private void addContribution(final SubHyperplane<Euclidean3D> facet, final boolean reversed) {
    160
    161                // extract the vertices of the facet
    162                @SuppressWarnings("unchecked")
    163                final AbstractSubHyperplane<Euclidean3D, Euclidean2D> absFacet =
    164                    (AbstractSubHyperplane<Euclidean3D, Euclidean2D>) facet;
    165                final Plane plane    = (Plane) facet.getHyperplane();
    166
    167                final double scal = plane.getNormal().dotProduct(w);
    168                if (FastMath.abs(scal) > 1.0e-3) {
    169                    Vector2D[][] vertices =
    170                        ((PolygonsSet) absFacet.getRemainingRegion()).getVertices();
    171
    172                    if ((scal < 0) ^ reversed) {
    173                        // the facet is seen from the inside,
    174                        // we need to invert its boundary orientation
    175                        final Vector2D[][] newVertices = new Vector2D[vertices.length][];
    176                        for (int i = 0; i < vertices.length; ++i) {
    177                            final Vector2D[] loop = vertices[i];
    178                            final Vector2D[] newLoop = new Vector2D[loop.length];
    179                            if (loop[0] == null) {
    180                                newLoop[0] = null;
    181                                for (int j = 1; j < loop.length; ++j) {
    182                                    newLoop[j] = loop[loop.length - j];
    183                                }
    184                            } else {
    185                                for (int j = 0; j < loop.length; ++j) {
    186                                    newLoop[j] = loop[loop.length - (j + 1)];
    187                                }
    188                            }
    189                            newVertices[i] = newLoop;
    190                        }
    191
    192                        // use the reverted vertices
    193                        vertices = newVertices;
    194
    195                    }
    196
    197                    // compute the projection of the facet in the outline plane
    198                    final ArrayList<SubHyperplane<Euclidean2D>> edges = new ArrayList<SubHyperplane<Euclidean2D>>();
    199                    for (Vector2D[] loop : vertices) {
    200                        final boolean closed = loop[0] != null;
    201                        int previous         = closed ? (loop.length - 1) : 1;
    202                        Vector3D previous3D  = plane.toSpace(loop[previous]);
    203                        int current          = (previous + 1) % loop.length;
    204                        Vector2D pPoint       = new Vector2D(previous3D.dotProduct(u),
    205                                                             previous3D.dotProduct(v));
    206                        while (current < loop.length) {
    207
    208                            final Vector3D current3D = plane.toSpace(loop[current]);
    209                            final Vector2D  cPoint    = new Vector2D(current3D.dotProduct(u),
    210                                                                     current3D.dotProduct(v));
    211                            final org.apache.commons.math3.geometry.euclidean.twod.Line line =
    212                                new org.apache.commons.math3.geometry.euclidean.twod.Line(pPoint, cPoint);
    213                            SubHyperplane<Euclidean2D> edge = line.wholeHyperplane();
    214
    215                            if (closed || (previous != 1)) {
    216                                // the previous point is a real vertex
    217                                // it defines one bounding point of the edge
    218                                final double angle = line.getAngle() + 0.5 * FastMath.PI;
    219                                final org.apache.commons.math3.geometry.euclidean.twod.Line l =
    220                                    new org.apache.commons.math3.geometry.euclidean.twod.Line(pPoint, angle);
    221                                edge = edge.split(l).getPlus();
    222                            }
    223
    224                            if (closed || (current != (loop.length - 1))) {
    225                                // the current point is a real vertex
    226                                // it defines one bounding point of the edge
    227                                final double angle = line.getAngle() + 0.5 * FastMath.PI;
    228                                final org.apache.commons.math3.geometry.euclidean.twod.Line l =
    229                                    new org.apache.commons.math3.geometry.euclidean.twod.Line(cPoint, angle);
    230                                edge = edge.split(l).getMinus();
    231                            }
    232
    233                            edges.add(edge);
    234
    235                            previous   = current++;
    236                            previous3D = current3D;
    237                            pPoint     = cPoint;
    238
    239                        }
    240                    }
    241                    final PolygonsSet projectedFacet = new PolygonsSet(edges);
    242
    243                    // add the contribution of the facet to the global outline
    244                    projected = (PolygonsSet) new RegionFactory<Euclidean2D>().union(projected, projectedFacet);
    245
    246                }
    247            }
    248
    249            /** Get the projection of the polyhedrons set on the plane.
    250             * @return projection of the polyhedrons set on the plane
    251             */
    252            public PolygonsSet getProjected() {
    253                return projected;
    254            }
    255
    256        }
    257
    258    }