View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math3.geometry.euclidean.threed;
18  
19  import java.util.ArrayList;
20  
21  import org.apache.commons.math3.geometry.Point;
22  import org.apache.commons.math3.geometry.euclidean.twod.Euclidean2D;
23  import org.apache.commons.math3.geometry.euclidean.twod.PolygonsSet;
24  import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
25  import org.apache.commons.math3.geometry.partitioning.AbstractSubHyperplane;
26  import org.apache.commons.math3.geometry.partitioning.BSPTree;
27  import org.apache.commons.math3.geometry.partitioning.BSPTreeVisitor;
28  import org.apache.commons.math3.geometry.partitioning.BoundaryAttribute;
29  import org.apache.commons.math3.geometry.partitioning.RegionFactory;
30  import org.apache.commons.math3.geometry.partitioning.SubHyperplane;
31  import org.apache.commons.math3.util.FastMath;
32  
33  /** Extractor for {@link PolygonsSet polyhedrons sets} outlines.
34   * <p>This class extracts the 2D outlines from {{@link PolygonsSet
35   * polyhedrons sets} in a specified projection plane.</p>
36   * @version $Id: OutlineExtractor.java 1555174 2014-01-03 18:06:20Z luc $
37   * @since 3.0
38   */
39  public class OutlineExtractor {
40  
41      /** Abscissa axis of the projection plane. */
42      private Vector3D u;
43  
44      /** Ordinate axis of the projection plane. */
45      private Vector3D v;
46  
47      /** Normal of the projection plane (viewing direction). */
48      private Vector3D w;
49  
50      /** Build an extractor for a specific projection plane.
51       * @param u abscissa axis of the projection point
52       * @param v ordinate axis of the projection point
53       */
54      public OutlineExtractor(final Vector3D u, final Vector3D v) {
55          this.u = u;
56          this.v = v;
57          w = Vector3D.crossProduct(u, v);
58      }
59  
60      /** Extract the outline of a polyhedrons set.
61       * @param polyhedronsSet polyhedrons set whose outline must be extracted
62       * @return an outline, as an array of loops.
63       */
64      public Vector2D[][] getOutline(final PolyhedronsSet polyhedronsSet) {
65  
66          // project all boundary facets into one polygons set
67          final BoundaryProjector projector = new BoundaryProjector(polyhedronsSet.getTolerance());
68          polyhedronsSet.getTree(true).visit(projector);
69          final PolygonsSet projected = projector.getProjected();
70  
71          // Remove the spurious intermediate vertices from the outline
72          final Vector2D[][] outline = projected.getVertices();
73          for (int i = 0; i < outline.length; ++i) {
74              final Vector2D[] rawLoop = outline[i];
75              int end = rawLoop.length;
76              int j = 0;
77              while (j < end) {
78                  if (pointIsBetween(rawLoop, end, j)) {
79                      // the point should be removed
80                      for (int k = j; k < (end - 1); ++k) {
81                          rawLoop[k] = rawLoop[k + 1];
82                      }
83                      --end;
84                  } else {
85                      // the point remains in the loop
86                      ++j;
87                  }
88              }
89              if (end != rawLoop.length) {
90                  // resize the array
91                  outline[i] = new Vector2D[end];
92                  System.arraycopy(rawLoop, 0, outline[i], 0, end);
93              }
94          }
95  
96          return outline;
97  
98      }
99  
100     /** Check if a point is geometrically between its neighbor in an array.
101      * <p>The neighbors are computed considering the array is a loop
102      * (i.e. point at index (n-1) is before point at index 0)</p>
103      * @param loop points array
104      * @param n number of points to consider in the array
105      * @param i index of the point to check (must be between 0 and n-1)
106      * @return true if the point is exactly between its neighbors
107      */
108     private boolean pointIsBetween(final Vector2D[] loop, final int n, final int i) {
109         final Vector2D previous = loop[(i + n - 1) % n];
110         final Vector2D current  = loop[i];
111         final Vector2D next     = loop[(i + 1) % n];
112         final double dx1       = current.getX() - previous.getX();
113         final double dy1       = current.getY() - previous.getY();
114         final double dx2       = next.getX()    - current.getX();
115         final double dy2       = next.getY()    - current.getY();
116         final double cross     = dx1 * dy2 - dx2 * dy1;
117         final double dot       = dx1 * dx2 + dy1 * dy2;
118         final double d1d2      = FastMath.sqrt((dx1 * dx1 + dy1 * dy1) * (dx2 * dx2 + dy2 * dy2));
119         return (FastMath.abs(cross) <= (1.0e-6 * d1d2)) && (dot >= 0.0);
120     }
121 
122     /** Visitor projecting the boundary facets on a plane. */
123     private class BoundaryProjector implements BSPTreeVisitor<Euclidean3D> {
124 
125         /** Projection of the polyhedrons set on the plane. */
126         private PolygonsSet projected;
127 
128         /** Tolerance below which points are considered identical. */
129         private final double tolerance;
130 
131         /** Simple constructor.
132          * @param tolerance tolerance below which points are considered identical
133          */
134         public BoundaryProjector(final double tolerance) {
135             this.projected = new PolygonsSet(new BSPTree<Euclidean2D>(Boolean.FALSE), tolerance);
136             this.tolerance = tolerance;
137         }
138 
139         /** {@inheritDoc} */
140         public Order visitOrder(final BSPTree<Euclidean3D> node) {
141             return Order.MINUS_SUB_PLUS;
142         }
143 
144         /** {@inheritDoc} */
145         public void visitInternalNode(final BSPTree<Euclidean3D> node) {
146             @SuppressWarnings("unchecked")
147             final BoundaryAttribute<Euclidean3D> attribute =
148                 (BoundaryAttribute<Euclidean3D>) node.getAttribute();
149             if (attribute.getPlusOutside() != null) {
150                 addContribution(attribute.getPlusOutside(), false);
151             }
152             if (attribute.getPlusInside() != null) {
153                 addContribution(attribute.getPlusInside(), true);
154             }
155         }
156 
157         /** {@inheritDoc} */
158         public void visitLeafNode(final BSPTree<Euclidean3D> node) {
159         }
160 
161         /** Add he contribution of a boundary facet.
162          * @param facet boundary facet
163          * @param reversed if true, the facet has the inside on its plus side
164          */
165         private void addContribution(final SubHyperplane<Euclidean3D> facet, final boolean reversed) {
166 
167             // extract the vertices of the facet
168             @SuppressWarnings("unchecked")
169             final AbstractSubHyperplane<Euclidean3D, Euclidean2D> absFacet =
170                 (AbstractSubHyperplane<Euclidean3D, Euclidean2D>) facet;
171             final Plane plane    = (Plane) facet.getHyperplane();
172 
173             final double scal = plane.getNormal().dotProduct(w);
174             if (FastMath.abs(scal) > 1.0e-3) {
175                 Vector2D[][] vertices =
176                     ((PolygonsSet) absFacet.getRemainingRegion()).getVertices();
177 
178                 if ((scal < 0) ^ reversed) {
179                     // the facet is seen from the inside,
180                     // we need to invert its boundary orientation
181                     final Vector2D[][] newVertices = new Vector2D[vertices.length][];
182                     for (int i = 0; i < vertices.length; ++i) {
183                         final Vector2D[] loop = vertices[i];
184                         final Vector2D[] newLoop = new Vector2D[loop.length];
185                         if (loop[0] == null) {
186                             newLoop[0] = null;
187                             for (int j = 1; j < loop.length; ++j) {
188                                 newLoop[j] = loop[loop.length - j];
189                             }
190                         } else {
191                             for (int j = 0; j < loop.length; ++j) {
192                                 newLoop[j] = loop[loop.length - (j + 1)];
193                             }
194                         }
195                         newVertices[i] = newLoop;
196                     }
197 
198                     // use the reverted vertices
199                     vertices = newVertices;
200 
201                 }
202 
203                 // compute the projection of the facet in the outline plane
204                 final ArrayList<SubHyperplane<Euclidean2D>> edges = new ArrayList<SubHyperplane<Euclidean2D>>();
205                 for (Vector2D[] loop : vertices) {
206                     final boolean closed = loop[0] != null;
207                     int previous         = closed ? (loop.length - 1) : 1;
208                     Vector3D previous3D  = plane.toSpace((Point<Euclidean2D>) loop[previous]);
209                     int current          = (previous + 1) % loop.length;
210                     Vector2D pPoint       = new Vector2D(previous3D.dotProduct(u),
211                                                          previous3D.dotProduct(v));
212                     while (current < loop.length) {
213 
214                         final Vector3D current3D = plane.toSpace((Point<Euclidean2D>) loop[current]);
215                         final Vector2D  cPoint    = new Vector2D(current3D.dotProduct(u),
216                                                                  current3D.dotProduct(v));
217                         final org.apache.commons.math3.geometry.euclidean.twod.Line line =
218                             new org.apache.commons.math3.geometry.euclidean.twod.Line(pPoint, cPoint, tolerance);
219                         SubHyperplane<Euclidean2D> edge = line.wholeHyperplane();
220 
221                         if (closed || (previous != 1)) {
222                             // the previous point is a real vertex
223                             // it defines one bounding point of the edge
224                             final double angle = line.getAngle() + 0.5 * FastMath.PI;
225                             final org.apache.commons.math3.geometry.euclidean.twod.Line l =
226                                 new org.apache.commons.math3.geometry.euclidean.twod.Line(pPoint, angle, tolerance);
227                             edge = edge.split(l).getPlus();
228                         }
229 
230                         if (closed || (current != (loop.length - 1))) {
231                             // the current point is a real vertex
232                             // it defines one bounding point of the edge
233                             final double angle = line.getAngle() + 0.5 * FastMath.PI;
234                             final org.apache.commons.math3.geometry.euclidean.twod.Line l =
235                                 new org.apache.commons.math3.geometry.euclidean.twod.Line(cPoint, angle, tolerance);
236                             edge = edge.split(l).getMinus();
237                         }
238 
239                         edges.add(edge);
240 
241                         previous   = current++;
242                         previous3D = current3D;
243                         pPoint     = cPoint;
244 
245                     }
246                 }
247                 final PolygonsSet projectedFacet = new PolygonsSet(edges, tolerance);
248 
249                 // add the contribution of the facet to the global outline
250                 projected = (PolygonsSet) new RegionFactory<Euclidean2D>().union(projected, projectedFacet);
251 
252             }
253         }
254 
255         /** Get the projection of the polyhedrons set on the plane.
256          * @return projection of the polyhedrons set on the plane
257          */
258         public PolygonsSet getProjected() {
259             return projected;
260         }
261 
262     }
263 
264 }