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.util.ArrayList;
020
021import org.apache.commons.math3.geometry.Point;
022import org.apache.commons.math3.geometry.euclidean.twod.Euclidean2D;
023import org.apache.commons.math3.geometry.euclidean.twod.PolygonsSet;
024import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
025import org.apache.commons.math3.geometry.partitioning.AbstractSubHyperplane;
026import org.apache.commons.math3.geometry.partitioning.BSPTree;
027import org.apache.commons.math3.geometry.partitioning.BSPTreeVisitor;
028import org.apache.commons.math3.geometry.partitioning.BoundaryAttribute;
029import org.apache.commons.math3.geometry.partitioning.RegionFactory;
030import org.apache.commons.math3.geometry.partitioning.SubHyperplane;
031import org.apache.commons.math3.util.FastMath;
032
033/** Extractor for {@link PolygonsSet polyhedrons sets} outlines.
034 * <p>This class extracts the 2D outlines from {{@link PolygonsSet
035 * polyhedrons sets} in a specified projection plane.</p>
036 * @since 3.0
037 */
038public 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(polyhedronsSet.getTolerance());
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 neighbor in an array.
100     * <p>The neighbors 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 neighbors
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        /** Tolerance below which points are considered identical. */
128        private final double tolerance;
129
130        /** Simple constructor.
131         * @param tolerance tolerance below which points are considered identical
132         */
133        BoundaryProjector(final double tolerance) {
134            this.projected = new PolygonsSet(new BSPTree<Euclidean2D>(Boolean.FALSE), tolerance);
135            this.tolerance = tolerance;
136        }
137
138        /** {@inheritDoc} */
139        public Order visitOrder(final BSPTree<Euclidean3D> node) {
140            return Order.MINUS_SUB_PLUS;
141        }
142
143        /** {@inheritDoc} */
144        public void visitInternalNode(final BSPTree<Euclidean3D> node) {
145            @SuppressWarnings("unchecked")
146            final BoundaryAttribute<Euclidean3D> attribute =
147                (BoundaryAttribute<Euclidean3D>) node.getAttribute();
148            if (attribute.getPlusOutside() != null) {
149                addContribution(attribute.getPlusOutside(), false);
150            }
151            if (attribute.getPlusInside() != null) {
152                addContribution(attribute.getPlusInside(), true);
153            }
154        }
155
156        /** {@inheritDoc} */
157        public void visitLeafNode(final BSPTree<Euclidean3D> node) {
158        }
159
160        /** Add he contribution of a boundary facet.
161         * @param facet boundary facet
162         * @param reversed if true, the facet has the inside on its plus side
163         */
164        private void addContribution(final SubHyperplane<Euclidean3D> facet, final boolean reversed) {
165
166            // extract the vertices of the facet
167            @SuppressWarnings("unchecked")
168            final AbstractSubHyperplane<Euclidean3D, Euclidean2D> absFacet =
169                (AbstractSubHyperplane<Euclidean3D, Euclidean2D>) facet;
170            final Plane plane    = (Plane) facet.getHyperplane();
171
172            final double scal = plane.getNormal().dotProduct(w);
173            if (FastMath.abs(scal) > 1.0e-3) {
174                Vector2D[][] vertices =
175                    ((PolygonsSet) absFacet.getRemainingRegion()).getVertices();
176
177                if ((scal < 0) ^ reversed) {
178                    // the facet is seen from the inside,
179                    // we need to invert its boundary orientation
180                    final Vector2D[][] newVertices = new Vector2D[vertices.length][];
181                    for (int i = 0; i < vertices.length; ++i) {
182                        final Vector2D[] loop = vertices[i];
183                        final Vector2D[] newLoop = new Vector2D[loop.length];
184                        if (loop[0] == null) {
185                            newLoop[0] = null;
186                            for (int j = 1; j < loop.length; ++j) {
187                                newLoop[j] = loop[loop.length - j];
188                            }
189                        } else {
190                            for (int j = 0; j < loop.length; ++j) {
191                                newLoop[j] = loop[loop.length - (j + 1)];
192                            }
193                        }
194                        newVertices[i] = newLoop;
195                    }
196
197                    // use the reverted vertices
198                    vertices = newVertices;
199
200                }
201
202                // compute the projection of the facet in the outline plane
203                final ArrayList<SubHyperplane<Euclidean2D>> edges = new ArrayList<SubHyperplane<Euclidean2D>>();
204                for (Vector2D[] loop : vertices) {
205                    final boolean closed = loop[0] != null;
206                    int previous         = closed ? (loop.length - 1) : 1;
207                    Vector3D previous3D  = plane.toSpace((Point<Euclidean2D>) loop[previous]);
208                    int current          = (previous + 1) % loop.length;
209                    Vector2D pPoint       = new Vector2D(previous3D.dotProduct(u),
210                                                         previous3D.dotProduct(v));
211                    while (current < loop.length) {
212
213                        final Vector3D current3D = plane.toSpace((Point<Euclidean2D>) loop[current]);
214                        final Vector2D  cPoint    = new Vector2D(current3D.dotProduct(u),
215                                                                 current3D.dotProduct(v));
216                        final org.apache.commons.math3.geometry.euclidean.twod.Line line =
217                            new org.apache.commons.math3.geometry.euclidean.twod.Line(pPoint, cPoint, tolerance);
218                        SubHyperplane<Euclidean2D> edge = line.wholeHyperplane();
219
220                        if (closed || (previous != 1)) {
221                            // the previous point is a real vertex
222                            // it defines one bounding point of the edge
223                            final double angle = line.getAngle() + 0.5 * FastMath.PI;
224                            final org.apache.commons.math3.geometry.euclidean.twod.Line l =
225                                new org.apache.commons.math3.geometry.euclidean.twod.Line(pPoint, angle, tolerance);
226                            edge = edge.split(l).getPlus();
227                        }
228
229                        if (closed || (current != (loop.length - 1))) {
230                            // the current point is a real vertex
231                            // it defines one bounding point of the edge
232                            final double angle = line.getAngle() + 0.5 * FastMath.PI;
233                            final org.apache.commons.math3.geometry.euclidean.twod.Line l =
234                                new org.apache.commons.math3.geometry.euclidean.twod.Line(cPoint, angle, tolerance);
235                            edge = edge.split(l).getMinus();
236                        }
237
238                        edges.add(edge);
239
240                        previous   = current++;
241                        previous3D = current3D;
242                        pPoint     = cPoint;
243
244                    }
245                }
246                final PolygonsSet projectedFacet = new PolygonsSet(edges, tolerance);
247
248                // add the contribution of the facet to the global outline
249                projected = (PolygonsSet) new RegionFactory<Euclidean2D>().union(projected, projectedFacet);
250
251            }
252        }
253
254        /** Get the projection of the polyhedrons set on the plane.
255         * @return projection of the polyhedrons set on the plane
256         */
257        public PolygonsSet getProjected() {
258            return projected;
259        }
260
261    }
262
263}