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.twod;
018
019import java.util.ArrayList;
020import java.util.Collection;
021import java.util.List;
022
023import org.apache.commons.math3.geometry.Point;
024import org.apache.commons.math3.geometry.euclidean.oned.Euclidean1D;
025import org.apache.commons.math3.geometry.euclidean.oned.Interval;
026import org.apache.commons.math3.geometry.euclidean.oned.IntervalsSet;
027import org.apache.commons.math3.geometry.euclidean.oned.Vector1D;
028import org.apache.commons.math3.geometry.partitioning.AbstractRegion;
029import org.apache.commons.math3.geometry.partitioning.AbstractSubHyperplane;
030import org.apache.commons.math3.geometry.partitioning.BSPTree;
031import org.apache.commons.math3.geometry.partitioning.BSPTreeVisitor;
032import org.apache.commons.math3.geometry.partitioning.BoundaryAttribute;
033import org.apache.commons.math3.geometry.partitioning.Hyperplane;
034import org.apache.commons.math3.geometry.partitioning.Side;
035import org.apache.commons.math3.geometry.partitioning.SubHyperplane;
036import org.apache.commons.math3.util.FastMath;
037import org.apache.commons.math3.util.Precision;
038
039/** This class represents a 2D region: a set of polygons.
040 * @since 3.0
041 */
042public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> {
043
044    /** Default value for tolerance. */
045    private static final double DEFAULT_TOLERANCE = 1.0e-10;
046
047    /** Vertices organized as boundary loops. */
048    private Vector2D[][] vertices;
049
050    /** Build a polygons set representing the whole plane.
051     * @param tolerance tolerance below which points are considered identical
052     * @since 3.3
053     */
054    public PolygonsSet(final double tolerance) {
055        super(tolerance);
056    }
057
058    /** Build a polygons set from a BSP tree.
059     * <p>The leaf nodes of the BSP tree <em>must</em> have a
060     * {@code Boolean} attribute representing the inside status of
061     * the corresponding cell (true for inside cells, false for outside
062     * cells). In order to avoid building too many small objects, it is
063     * recommended to use the predefined constants
064     * {@code Boolean.TRUE} and {@code Boolean.FALSE}</p>
065     * <p>
066     * This constructor is aimed at expert use, as building the tree may
067     * be a difficult task. It is not intended for general use and for
068     * performances reasons does not check thoroughly its input, as this would
069     * require walking the full tree each time. Failing to provide a tree with
070     * the proper attributes, <em>will</em> therefore generate problems like
071     * {@link NullPointerException} or {@link ClassCastException} only later on.
072     * This limitation is known and explains why this constructor is for expert
073     * use only. The caller does have the responsibility to provided correct arguments.
074     * </p>
075     * @param tree inside/outside BSP tree representing the region
076     * @param tolerance tolerance below which points are considered identical
077     * @since 3.3
078     */
079    public PolygonsSet(final BSPTree<Euclidean2D> tree, final double tolerance) {
080        super(tree, tolerance);
081    }
082
083    /** Build a polygons set from a Boundary REPresentation (B-rep).
084     * <p>The boundary is provided as a collection of {@link
085     * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the
086     * interior part of the region on its minus side and the exterior on
087     * its plus side.</p>
088     * <p>The boundary elements can be in any order, and can form
089     * several non-connected sets (like for example polygons with holes
090     * or a set of disjoint polygons considered as a whole). In
091     * fact, the elements do not even need to be connected together
092     * (their topological connections are not used here). However, if the
093     * boundary does not really separate an inside open from an outside
094     * open (open having here its topological meaning), then subsequent
095     * calls to the {@link
096     * org.apache.commons.math3.geometry.partitioning.Region#checkPoint(org.apache.commons.math3.geometry.Point)
097     * checkPoint} method will not be meaningful anymore.</p>
098     * <p>If the boundary is empty, the region will represent the whole
099     * space.</p>
100     * @param boundary collection of boundary elements, as a
101     * collection of {@link SubHyperplane SubHyperplane} objects
102     * @param tolerance tolerance below which points are considered identical
103     * @since 3.3
104     */
105    public PolygonsSet(final Collection<SubHyperplane<Euclidean2D>> boundary, final double tolerance) {
106        super(boundary, tolerance);
107    }
108
109    /** Build a parallellepipedic box.
110     * @param xMin low bound along the x direction
111     * @param xMax high bound along the x direction
112     * @param yMin low bound along the y direction
113     * @param yMax high bound along the y direction
114     * @param tolerance tolerance below which points are considered identical
115     * @since 3.3
116     */
117    public PolygonsSet(final double xMin, final double xMax,
118                       final double yMin, final double yMax,
119                       final double tolerance) {
120        super(boxBoundary(xMin, xMax, yMin, yMax, tolerance), tolerance);
121    }
122
123    /** Build a polygon from a simple list of vertices.
124     * <p>The boundary is provided as a list of points considering to
125     * represent the vertices of a simple loop. The interior part of the
126     * region is on the left side of this path and the exterior is on its
127     * right side.</p>
128     * <p>This constructor does not handle polygons with a boundary
129     * forming several disconnected paths (such as polygons with holes).</p>
130     * <p>For cases where this simple constructor applies, it is expected to
131     * be numerically more robust than the {@link #PolygonsSet(Collection) general
132     * constructor} using {@link SubHyperplane subhyperplanes}.</p>
133     * <p>If the list is empty, the region will represent the whole
134     * space.</p>
135     * <p>
136     * Polygons with thin pikes or dents are inherently difficult to handle because
137     * they involve lines with almost opposite directions at some vertices. Polygons
138     * whose vertices come from some physical measurement with noise are also
139     * difficult because an edge that should be straight may be broken in lots of
140     * different pieces with almost equal directions. In both cases, computing the
141     * lines intersections is not numerically robust due to the almost 0 or almost
142     * &pi; angle. Such cases need to carefully adjust the {@code hyperplaneThickness}
143     * parameter. A too small value would often lead to completely wrong polygons
144     * with large area wrongly identified as inside or outside. Large values are
145     * often much safer. As a rule of thumb, a value slightly below the size of the
146     * most accurate detail needed is a good value for the {@code hyperplaneThickness}
147     * parameter.
148     * </p>
149     * @param hyperplaneThickness tolerance below which points are considered to
150     * belong to the hyperplane (which is therefore more a slab)
151     * @param vertices vertices of the simple loop boundary
152     */
153    public PolygonsSet(final double hyperplaneThickness, final Vector2D ... vertices) {
154        super(verticesToTree(hyperplaneThickness, vertices), hyperplaneThickness);
155    }
156
157    /** Build a polygons set representing the whole real line.
158     * @deprecated as of 3.3, replaced with {@link #PolygonsSet(double)}
159     */
160    @Deprecated
161    public PolygonsSet() {
162        this(DEFAULT_TOLERANCE);
163    }
164
165    /** Build a polygons set from a BSP tree.
166     * <p>The leaf nodes of the BSP tree <em>must</em> have a
167     * {@code Boolean} attribute representing the inside status of
168     * the corresponding cell (true for inside cells, false for outside
169     * cells). In order to avoid building too many small objects, it is
170     * recommended to use the predefined constants
171     * {@code Boolean.TRUE} and {@code Boolean.FALSE}</p>
172     * @param tree inside/outside BSP tree representing the region
173     * @deprecated as of 3.3, replaced with {@link #PolygonsSet(BSPTree, double)}
174     */
175    @Deprecated
176    public PolygonsSet(final BSPTree<Euclidean2D> tree) {
177        this(tree, DEFAULT_TOLERANCE);
178    }
179
180    /** Build a polygons set from a Boundary REPresentation (B-rep).
181     * <p>The boundary is provided as a collection of {@link
182     * SubHyperplane sub-hyperplanes}. Each sub-hyperplane has the
183     * interior part of the region on its minus side and the exterior on
184     * its plus side.</p>
185     * <p>The boundary elements can be in any order, and can form
186     * several non-connected sets (like for example polygons with holes
187     * or a set of disjoint polygons considered as a whole). In
188     * fact, the elements do not even need to be connected together
189     * (their topological connections are not used here). However, if the
190     * boundary does not really separate an inside open from an outside
191     * open (open having here its topological meaning), then subsequent
192     * calls to the {@link
193     * org.apache.commons.math3.geometry.partitioning.Region#checkPoint(org.apache.commons.math3.geometry.Point)
194     * checkPoint} method will not be meaningful anymore.</p>
195     * <p>If the boundary is empty, the region will represent the whole
196     * space.</p>
197     * @param boundary collection of boundary elements, as a
198     * collection of {@link SubHyperplane SubHyperplane} objects
199     * @deprecated as of 3.3, replaced with {@link #PolygonsSet(Collection, double)}
200     */
201    @Deprecated
202    public PolygonsSet(final Collection<SubHyperplane<Euclidean2D>> boundary) {
203        this(boundary, DEFAULT_TOLERANCE);
204    }
205
206    /** Build a parallellepipedic box.
207     * @param xMin low bound along the x direction
208     * @param xMax high bound along the x direction
209     * @param yMin low bound along the y direction
210     * @param yMax high bound along the y direction
211     * @deprecated as of 3.3, replaced with {@link #PolygonsSet(double, double, double, double, double)}
212     */
213    @Deprecated
214    public PolygonsSet(final double xMin, final double xMax,
215                       final double yMin, final double yMax) {
216        this(xMin, xMax, yMin, yMax, DEFAULT_TOLERANCE);
217    }
218
219    /** Create a list of hyperplanes representing the boundary of a box.
220     * @param xMin low bound along the x direction
221     * @param xMax high bound along the x direction
222     * @param yMin low bound along the y direction
223     * @param yMax high bound along the y direction
224     * @param tolerance tolerance below which points are considered identical
225     * @return boundary of the box
226     */
227    private static Line[] boxBoundary(final double xMin, final double xMax,
228                                      final double yMin, final double yMax,
229                                      final double tolerance) {
230        if ((xMin >= xMax - tolerance) || (yMin >= yMax - tolerance)) {
231            // too thin box, build an empty polygons set
232            return null;
233        }
234        final Vector2D minMin = new Vector2D(xMin, yMin);
235        final Vector2D minMax = new Vector2D(xMin, yMax);
236        final Vector2D maxMin = new Vector2D(xMax, yMin);
237        final Vector2D maxMax = new Vector2D(xMax, yMax);
238        return new Line[] {
239            new Line(minMin, maxMin, tolerance),
240            new Line(maxMin, maxMax, tolerance),
241            new Line(maxMax, minMax, tolerance),
242            new Line(minMax, minMin, tolerance)
243        };
244    }
245
246    /** Build the BSP tree of a polygons set from a simple list of vertices.
247     * <p>The boundary is provided as a list of points considering to
248     * represent the vertices of a simple loop. The interior part of the
249     * region is on the left side of this path and the exterior is on its
250     * right side.</p>
251     * <p>This constructor does not handle polygons with a boundary
252     * forming several disconnected paths (such as polygons with holes).</p>
253     * <p>For cases where this simple constructor applies, it is expected to
254     * be numerically more robust than the {@link #PolygonsSet(Collection) general
255     * constructor} using {@link SubHyperplane subhyperplanes}.</p>
256     * @param hyperplaneThickness tolerance below which points are consider to
257     * belong to the hyperplane (which is therefore more a slab)
258     * @param vertices vertices of the simple loop boundary
259     * @return the BSP tree of the input vertices
260     */
261    private static BSPTree<Euclidean2D> verticesToTree(final double hyperplaneThickness,
262                                                       final Vector2D ... vertices) {
263
264        final int n = vertices.length;
265        if (n == 0) {
266            // the tree represents the whole space
267            return new BSPTree<Euclidean2D>(Boolean.TRUE);
268        }
269
270        // build the vertices
271        final Vertex[] vArray = new Vertex[n];
272        for (int i = 0; i < n; ++i) {
273            vArray[i] = new Vertex(vertices[i]);
274        }
275
276        // build the edges
277        List<Edge> edges = new ArrayList<Edge>(n);
278        for (int i = 0; i < n; ++i) {
279
280            // get the endpoints of the edge
281            final Vertex start = vArray[i];
282            final Vertex end   = vArray[(i + 1) % n];
283
284            // get the line supporting the edge, taking care not to recreate it
285            // if it was already created earlier due to another edge being aligned
286            // with the current one
287            Line line = start.sharedLineWith(end);
288            if (line == null) {
289                line = new Line(start.getLocation(), end.getLocation(), hyperplaneThickness);
290            }
291
292            // create the edge and store it
293            edges.add(new Edge(start, end, line));
294
295            // check if another vertex also happens to be on this line
296            for (final Vertex vertex : vArray) {
297                if (vertex != start && vertex != end &&
298                    FastMath.abs(line.getOffset((Point<Euclidean2D>) vertex.getLocation())) <= hyperplaneThickness) {
299                    vertex.bindWith(line);
300                }
301            }
302
303        }
304
305        // build the tree top-down
306        final BSPTree<Euclidean2D> tree = new BSPTree<Euclidean2D>();
307        insertEdges(hyperplaneThickness, tree, edges);
308
309        return tree;
310
311    }
312
313    /** Recursively build a tree by inserting cut sub-hyperplanes.
314     * @param hyperplaneThickness tolerance below which points are consider to
315     * belong to the hyperplane (which is therefore more a slab)
316     * @param node current tree node (it is a leaf node at the beginning
317     * of the call)
318     * @param edges list of edges to insert in the cell defined by this node
319     * (excluding edges not belonging to the cell defined by this node)
320     */
321    private static void insertEdges(final double hyperplaneThickness,
322                                    final BSPTree<Euclidean2D> node,
323                                    final List<Edge> edges) {
324
325        // find an edge with an hyperplane that can be inserted in the node
326        int index = 0;
327        Edge inserted =null;
328        while (inserted == null && index < edges.size()) {
329            inserted = edges.get(index++);
330            if (inserted.getNode() == null) {
331                if (node.insertCut(inserted.getLine())) {
332                    inserted.setNode(node);
333                } else {
334                    inserted = null;
335                }
336            } else {
337                inserted = null;
338            }
339        }
340
341        if (inserted == null) {
342            // no suitable edge was found, the node remains a leaf node
343            // we need to set its inside/outside boolean indicator
344            final BSPTree<Euclidean2D> parent = node.getParent();
345            if (parent == null || node == parent.getMinus()) {
346                node.setAttribute(Boolean.TRUE);
347            } else {
348                node.setAttribute(Boolean.FALSE);
349            }
350            return;
351        }
352
353        // we have split the node by inserting an edge as a cut sub-hyperplane
354        // distribute the remaining edges in the two sub-trees
355        final List<Edge> plusList  = new ArrayList<Edge>();
356        final List<Edge> minusList = new ArrayList<Edge>();
357        for (final Edge edge : edges) {
358            if (edge != inserted) {
359                final double startOffset = inserted.getLine().getOffset((Point<Euclidean2D>) edge.getStart().getLocation());
360                final double endOffset   = inserted.getLine().getOffset((Point<Euclidean2D>) edge.getEnd().getLocation());
361                Side startSide = (FastMath.abs(startOffset) <= hyperplaneThickness) ?
362                                 Side.HYPER : ((startOffset < 0) ? Side.MINUS : Side.PLUS);
363                Side endSide   = (FastMath.abs(endOffset) <= hyperplaneThickness) ?
364                                 Side.HYPER : ((endOffset < 0) ? Side.MINUS : Side.PLUS);
365                switch (startSide) {
366                    case PLUS:
367                        if (endSide == Side.MINUS) {
368                            // we need to insert a split point on the hyperplane
369                            final Vertex splitPoint = edge.split(inserted.getLine());
370                            minusList.add(splitPoint.getOutgoing());
371                            plusList.add(splitPoint.getIncoming());
372                        } else {
373                            plusList.add(edge);
374                        }
375                        break;
376                    case MINUS:
377                        if (endSide == Side.PLUS) {
378                            // we need to insert a split point on the hyperplane
379                            final Vertex splitPoint = edge.split(inserted.getLine());
380                            minusList.add(splitPoint.getIncoming());
381                            plusList.add(splitPoint.getOutgoing());
382                        } else {
383                            minusList.add(edge);
384                        }
385                        break;
386                    default:
387                        if (endSide == Side.PLUS) {
388                            plusList.add(edge);
389                        } else if (endSide == Side.MINUS) {
390                            minusList.add(edge);
391                        }
392                        break;
393                }
394            }
395        }
396
397        // recurse through lower levels
398        if (!plusList.isEmpty()) {
399            insertEdges(hyperplaneThickness, node.getPlus(),  plusList);
400        } else {
401            node.getPlus().setAttribute(Boolean.FALSE);
402        }
403        if (!minusList.isEmpty()) {
404            insertEdges(hyperplaneThickness, node.getMinus(), minusList);
405        } else {
406            node.getMinus().setAttribute(Boolean.TRUE);
407        }
408
409    }
410
411    /** Internal class for holding vertices while they are processed to build a BSP tree. */
412    private static class Vertex {
413
414        /** Vertex location. */
415        private final Vector2D location;
416
417        /** Incoming edge. */
418        private Edge incoming;
419
420        /** Outgoing edge. */
421        private Edge outgoing;
422
423        /** Lines bound with this vertex. */
424        private final List<Line> lines;
425
426        /** Build a non-processed vertex not owned by any node yet.
427         * @param location vertex location
428         */
429        public Vertex(final Vector2D location) {
430            this.location = location;
431            this.incoming = null;
432            this.outgoing = null;
433            this.lines    = new ArrayList<Line>();
434        }
435
436        /** Get Vertex location.
437         * @return vertex location
438         */
439        public Vector2D getLocation() {
440            return location;
441        }
442
443        /** Bind a line considered to contain this vertex.
444         * @param line line to bind with this vertex
445         */
446        public void bindWith(final Line line) {
447            lines.add(line);
448        }
449
450        /** Get the common line bound with both the instance and another vertex, if any.
451         * <p>
452         * When two vertices are both bound to the same line, this means they are
453         * already handled by node associated with this line, so there is no need
454         * to create a cut hyperplane for them.
455         * </p>
456         * @param vertex other vertex to check instance against
457         * @return line bound with both the instance and another vertex, or null if the
458         * two vertices do not share a line yet
459         */
460        public Line sharedLineWith(final Vertex vertex) {
461            for (final Line line1 : lines) {
462                for (final Line line2 : vertex.lines) {
463                    if (line1 == line2) {
464                        return line1;
465                    }
466                }
467            }
468            return null;
469        }
470
471        /** Set incoming edge.
472         * <p>
473         * The line supporting the incoming edge is automatically bound
474         * with the instance.
475         * </p>
476         * @param incoming incoming edge
477         */
478        public void setIncoming(final Edge incoming) {
479            this.incoming = incoming;
480            bindWith(incoming.getLine());
481        }
482
483        /** Get incoming edge.
484         * @return incoming edge
485         */
486        public Edge getIncoming() {
487            return incoming;
488        }
489
490        /** Set outgoing edge.
491         * <p>
492         * The line supporting the outgoing edge is automatically bound
493         * with the instance.
494         * </p>
495         * @param outgoing outgoing edge
496         */
497        public void setOutgoing(final Edge outgoing) {
498            this.outgoing = outgoing;
499            bindWith(outgoing.getLine());
500        }
501
502        /** Get outgoing edge.
503         * @return outgoing edge
504         */
505        public Edge getOutgoing() {
506            return outgoing;
507        }
508
509    }
510
511    /** Internal class for holding edges while they are processed to build a BSP tree. */
512    private static class Edge {
513
514        /** Start vertex. */
515        private final Vertex start;
516
517        /** End vertex. */
518        private final Vertex end;
519
520        /** Line supporting the edge. */
521        private final Line line;
522
523        /** Node whose cut hyperplane contains this edge. */
524        private BSPTree<Euclidean2D> node;
525
526        /** Build an edge not contained in any node yet.
527         * @param start start vertex
528         * @param end end vertex
529         * @param line line supporting the edge
530         */
531        public Edge(final Vertex start, final Vertex end, final Line line) {
532
533            this.start = start;
534            this.end   = end;
535            this.line  = line;
536            this.node  = null;
537
538            // connect the vertices back to the edge
539            start.setOutgoing(this);
540            end.setIncoming(this);
541
542        }
543
544        /** Get start vertex.
545         * @return start vertex
546         */
547        public Vertex getStart() {
548            return start;
549        }
550
551        /** Get end vertex.
552         * @return end vertex
553         */
554        public Vertex getEnd() {
555            return end;
556        }
557
558        /** Get the line supporting this edge.
559         * @return line supporting this edge
560         */
561        public Line getLine() {
562            return line;
563        }
564
565        /** Set the node whose cut hyperplane contains this edge.
566         * @param node node whose cut hyperplane contains this edge
567         */
568        public void setNode(final BSPTree<Euclidean2D> node) {
569            this.node = node;
570        }
571
572        /** Get the node whose cut hyperplane contains this edge.
573         * @return node whose cut hyperplane contains this edge
574         * (null if edge has not yet been inserted into the BSP tree)
575         */
576        public BSPTree<Euclidean2D> getNode() {
577            return node;
578        }
579
580        /** Split the edge.
581         * <p>
582         * Once split, this edge is not referenced anymore by the vertices,
583         * it is replaced by the two half-edges and an intermediate splitting
584         * vertex is introduced to connect these two halves.
585         * </p>
586         * @param splitLine line splitting the edge in two halves
587         * @return split vertex (its incoming and outgoing edges are the two halves)
588         */
589        public Vertex split(final Line splitLine) {
590            final Vertex splitVertex = new Vertex(line.intersection(splitLine));
591            splitVertex.bindWith(splitLine);
592            final Edge startHalf = new Edge(start, splitVertex, line);
593            final Edge endHalf   = new Edge(splitVertex, end, line);
594            startHalf.node = node;
595            endHalf.node   = node;
596            return splitVertex;
597        }
598
599    }
600
601    /** {@inheritDoc} */
602    @Override
603    public PolygonsSet buildNew(final BSPTree<Euclidean2D> tree) {
604        return new PolygonsSet(tree, getTolerance());
605    }
606
607    /** {@inheritDoc} */
608    @Override
609    protected void computeGeometricalProperties() {
610
611        final Vector2D[][] v = getVertices();
612
613        if (v.length == 0) {
614            final BSPTree<Euclidean2D> tree = getTree(false);
615            if (tree.getCut() == null && (Boolean) tree.getAttribute()) {
616                // the instance covers the whole space
617                setSize(Double.POSITIVE_INFINITY);
618                setBarycenter((Point<Euclidean2D>) Vector2D.NaN);
619            } else {
620                setSize(0);
621                setBarycenter((Point<Euclidean2D>) new Vector2D(0, 0));
622            }
623        } else if (v[0][0] == null) {
624            // there is at least one open-loop: the polygon is infinite
625            setSize(Double.POSITIVE_INFINITY);
626            setBarycenter((Point<Euclidean2D>) Vector2D.NaN);
627        } else {
628            // all loops are closed, we compute some integrals around the shape
629
630            double sum  = 0;
631            double sumX = 0;
632            double sumY = 0;
633
634            for (Vector2D[] loop : v) {
635                double x1 = loop[loop.length - 1].getX();
636                double y1 = loop[loop.length - 1].getY();
637                for (final Vector2D point : loop) {
638                    final double x0 = x1;
639                    final double y0 = y1;
640                    x1 = point.getX();
641                    y1 = point.getY();
642                    final double factor = x0 * y1 - y0 * x1;
643                    sum  += factor;
644                    sumX += factor * (x0 + x1);
645                    sumY += factor * (y0 + y1);
646                }
647            }
648
649            if (sum < 0) {
650                // the polygon as a finite outside surrounded by an infinite inside
651                setSize(Double.POSITIVE_INFINITY);
652                setBarycenter((Point<Euclidean2D>) Vector2D.NaN);
653            } else {
654                setSize(sum / 2);
655                setBarycenter((Point<Euclidean2D>) new Vector2D(sumX / (3 * sum), sumY / (3 * sum)));
656            }
657
658        }
659
660    }
661
662    /** Get the vertices of the polygon.
663     * <p>The polygon boundary can be represented as an array of loops,
664     * each loop being itself an array of vertices.</p>
665     * <p>In order to identify open loops which start and end by
666     * infinite edges, the open loops arrays start with a null point. In
667     * this case, the first non null point and the last point of the
668     * array do not represent real vertices, they are dummy points
669     * intended only to get the direction of the first and last edge. An
670     * open loop consisting of a single infinite line will therefore be
671     * represented by a three elements array with one null point
672     * followed by two dummy points. The open loops are always the first
673     * ones in the loops array.</p>
674     * <p>If the polygon has no boundary at all, a zero length loop
675     * array will be returned.</p>
676     * <p>All line segments in the various loops have the inside of the
677     * region on their left side and the outside on their right side
678     * when moving in the underlying line direction. This means that
679     * closed loops surrounding finite areas obey the direct
680     * trigonometric orientation.</p>
681     * @return vertices of the polygon, organized as oriented boundary
682     * loops with the open loops first (the returned value is guaranteed
683     * to be non-null)
684     */
685    public Vector2D[][] getVertices() {
686        if (vertices == null) {
687            if (getTree(false).getCut() == null) {
688                vertices = new Vector2D[0][];
689            } else {
690
691                // build the unconnected segments
692                final SegmentsBuilder visitor = new SegmentsBuilder(getTolerance());
693                getTree(true).visit(visitor);
694                final List<ConnectableSegment> segments = visitor.getSegments();
695
696                // connect all segments, using topological criteria first
697                // and using Euclidean distance only as a last resort
698                int pending = segments.size();
699                pending -= naturalFollowerConnections(segments);
700                if (pending > 0) {
701                    pending -= splitEdgeConnections(segments);
702                }
703                if (pending > 0) {
704                    pending -= closeVerticesConnections(segments);
705                }
706
707                // create the segment loops
708                final ArrayList<List<Segment>> loops = new ArrayList<List<Segment>>();
709                for (ConnectableSegment s = getUnprocessed(segments); s != null; s = getUnprocessed(segments)) {
710                    final List<Segment> loop = followLoop(s);
711                    if (loop != null) {
712                        if (loop.get(0).getStart() == null) {
713                            // this is an open loop, we put it on the front
714                            loops.add(0, loop);
715                        } else {
716                            // this is a closed loop, we put it on the back
717                            loops.add(loop);
718                        }
719                    }
720                }
721
722                // transform the loops in an array of arrays of points
723                vertices = new Vector2D[loops.size()][];
724                int i = 0;
725
726                for (final List<Segment> loop : loops) {
727                    if (loop.size() < 2 ||
728                        (loop.size() == 2 && loop.get(0).getStart() == null && loop.get(1).getEnd() == null)) {
729                        // single infinite line
730                        final Line line = loop.get(0).getLine();
731                        vertices[i++] = new Vector2D[] {
732                            null,
733                            line.toSpace((Point<Euclidean1D>) new Vector1D(-Float.MAX_VALUE)),
734                            line.toSpace((Point<Euclidean1D>) new Vector1D(+Float.MAX_VALUE))
735                        };
736                    } else if (loop.get(0).getStart() == null) {
737                        // open loop with at least one real point
738                        final Vector2D[] array = new Vector2D[loop.size() + 2];
739                        int j = 0;
740                        for (Segment segment : loop) {
741
742                            if (j == 0) {
743                                // null point and first dummy point
744                                double x = segment.getLine().toSubSpace((Point<Euclidean2D>) segment.getEnd()).getX();
745                                x -= FastMath.max(1.0, FastMath.abs(x / 2));
746                                array[j++] = null;
747                                array[j++] = segment.getLine().toSpace((Point<Euclidean1D>) new Vector1D(x));
748                            }
749
750                            if (j < (array.length - 1)) {
751                                // current point
752                                array[j++] = segment.getEnd();
753                            }
754
755                            if (j == (array.length - 1)) {
756                                // last dummy point
757                                double x = segment.getLine().toSubSpace((Point<Euclidean2D>) segment.getStart()).getX();
758                                x += FastMath.max(1.0, FastMath.abs(x / 2));
759                                array[j++] = segment.getLine().toSpace((Point<Euclidean1D>) new Vector1D(x));
760                            }
761
762                        }
763                        vertices[i++] = array;
764                    } else {
765                        final Vector2D[] array = new Vector2D[loop.size()];
766                        int j = 0;
767                        for (Segment segment : loop) {
768                            array[j++] = segment.getStart();
769                        }
770                        vertices[i++] = array;
771                    }
772                }
773
774            }
775        }
776
777        return vertices.clone();
778
779    }
780
781    /** Connect the segments using only natural follower information.
782     * @param segments segments complete segments list
783     * @return number of connections performed
784     */
785    private int naturalFollowerConnections(final List<ConnectableSegment> segments) {
786        int connected = 0;
787        for (final ConnectableSegment segment : segments) {
788            if (segment.getNext() == null) {
789                final BSPTree<Euclidean2D> node = segment.getNode();
790                final BSPTree<Euclidean2D> end  = segment.getEndNode();
791                for (final ConnectableSegment candidateNext : segments) {
792                    if (candidateNext.getPrevious()  == null &&
793                        candidateNext.getNode()      == end &&
794                        candidateNext.getStartNode() == node) {
795                        // connect the two segments
796                        segment.setNext(candidateNext);
797                        candidateNext.setPrevious(segment);
798                        ++connected;
799                        break;
800                    }
801                }
802            }
803        }
804        return connected;
805    }
806
807    /** Connect the segments resulting from a line splitting a straight edge.
808     * @param segments segments complete segments list
809     * @return number of connections performed
810     */
811    private int splitEdgeConnections(final List<ConnectableSegment> segments) {
812        int connected = 0;
813        for (final ConnectableSegment segment : segments) {
814            if (segment.getNext() == null) {
815                final Hyperplane<Euclidean2D> hyperplane = segment.getNode().getCut().getHyperplane();
816                final BSPTree<Euclidean2D> end  = segment.getEndNode();
817                for (final ConnectableSegment candidateNext : segments) {
818                    if (candidateNext.getPrevious()                      == null &&
819                        candidateNext.getNode().getCut().getHyperplane() == hyperplane &&
820                        candidateNext.getStartNode()                     == end) {
821                        // connect the two segments
822                        segment.setNext(candidateNext);
823                        candidateNext.setPrevious(segment);
824                        ++connected;
825                        break;
826                    }
827                }
828            }
829        }
830        return connected;
831    }
832
833    /** Connect the segments using Euclidean distance.
834     * <p>
835     * This connection heuristic should be used last, as it relies
836     * only on a fuzzy distance criterion.
837     * </p>
838     * @param segments segments complete segments list
839     * @return number of connections performed
840     */
841    private int closeVerticesConnections(final List<ConnectableSegment> segments) {
842        int connected = 0;
843        for (final ConnectableSegment segment : segments) {
844            if (segment.getNext() == null && segment.getEnd() != null) {
845                final Vector2D end = segment.getEnd();
846                ConnectableSegment selectedNext = null;
847                double min = Double.POSITIVE_INFINITY;
848                for (final ConnectableSegment candidateNext : segments) {
849                    if (candidateNext.getPrevious() == null && candidateNext.getStart() != null) {
850                        final double distance = Vector2D.distance(end, candidateNext.getStart());
851                        if (distance < min) {
852                            selectedNext = candidateNext;
853                            min          = distance;
854                        }
855                    }
856                }
857                if (min <= getTolerance()) {
858                    // connect the two segments
859                    segment.setNext(selectedNext);
860                    selectedNext.setPrevious(segment);
861                    ++connected;
862                }
863            }
864        }
865        return connected;
866    }
867
868    /** Get first unprocessed segment from a list.
869     * @param segments segments list
870     * @return first segment that has not been processed yet
871     * or null if all segments have been processed
872     */
873    private ConnectableSegment getUnprocessed(final List<ConnectableSegment> segments) {
874        for (final ConnectableSegment segment : segments) {
875            if (!segment.isProcessed()) {
876                return segment;
877            }
878        }
879        return null;
880    }
881
882    /** Build the loop containing a segment.
883     * <p>
884     * The segment put in the loop will be marked as processed.
885     * </p>
886     * @param defining segment used to define the loop
887     * @return loop containing the segment (may be null if the loop is a
888     * degenerated infinitely thin 2 points loop
889     */
890    private List<Segment> followLoop(final ConnectableSegment defining) {
891
892        final List<Segment> loop = new ArrayList<Segment>();
893        loop.add(defining);
894        defining.setProcessed(true);
895
896        // add segments in connection order
897        ConnectableSegment next = defining.getNext();
898        while (next != defining && next != null) {
899            loop.add(next);
900            next.setProcessed(true);
901            next = next.getNext();
902        }
903
904        if (next == null) {
905            // the loop is open and we have found its end,
906            // we need to find its start too
907            ConnectableSegment previous = defining.getPrevious();
908            while (previous != null) {
909                loop.add(0, previous);
910                previous.setProcessed(true);
911                previous = previous.getPrevious();
912            }
913        }
914
915        // filter out spurious vertices
916        filterSpuriousVertices(loop);
917
918        if (loop.size() == 2 && loop.get(0).getStart() != null) {
919            // this is a degenerated infinitely thin closed loop, we simply ignore it
920            return null;
921        } else {
922            return loop;
923        }
924
925    }
926
927    /** Filter out spurious vertices on straight lines (at machine precision).
928     * @param loop segments loop to filter (will be modified in-place)
929     */
930    private void filterSpuriousVertices(final List<Segment> loop) {
931        for (int i = 0; i < loop.size(); ++i) {
932            final Segment previous = loop.get(i);
933            int j = (i + 1) % loop.size();
934            final Segment next = loop.get(j);
935            if (next != null &&
936                Precision.equals(previous.getLine().getAngle(), next.getLine().getAngle(), Precision.EPSILON)) {
937                // the vertex between the two edges is a spurious one
938                // replace the two segments by a single one
939                loop.set(j, new Segment(previous.getStart(), next.getEnd(), previous.getLine()));
940                loop.remove(i--);
941            }
942        }
943    }
944
945    /** Private extension of Segment allowing connection. */
946    private static class ConnectableSegment extends Segment {
947
948        /** Node containing segment. */
949        private final BSPTree<Euclidean2D> node;
950
951        /** Node whose intersection with current node defines start point. */
952        private final BSPTree<Euclidean2D> startNode;
953
954        /** Node whose intersection with current node defines end point. */
955        private final BSPTree<Euclidean2D> endNode;
956
957        /** Previous segment. */
958        private ConnectableSegment previous;
959
960        /** Next segment. */
961        private ConnectableSegment next;
962
963        /** Indicator for completely processed segments. */
964        private boolean processed;
965
966        /** Build a segment.
967         * @param start start point of the segment
968         * @param end end point of the segment
969         * @param line line containing the segment
970         * @param node node containing the segment
971         * @param startNode node whose intersection with current node defines start point
972         * @param endNode node whose intersection with current node defines end point
973         */
974        public ConnectableSegment(final Vector2D start, final Vector2D end, final Line line,
975                                  final BSPTree<Euclidean2D> node,
976                                  final BSPTree<Euclidean2D> startNode,
977                                  final BSPTree<Euclidean2D> endNode) {
978            super(start, end, line);
979            this.node      = node;
980            this.startNode = startNode;
981            this.endNode   = endNode;
982            this.previous  = null;
983            this.next      = null;
984            this.processed = false;
985        }
986
987        /** Get the node containing segment.
988         * @return node containing segment
989         */
990        public BSPTree<Euclidean2D> getNode() {
991            return node;
992        }
993
994        /** Get the node whose intersection with current node defines start point.
995         * @return node whose intersection with current node defines start point
996         */
997        public BSPTree<Euclidean2D> getStartNode() {
998            return startNode;
999        }
1000
1001        /** Get the node whose intersection with current node defines end point.
1002         * @return node whose intersection with current node defines end point
1003         */
1004        public BSPTree<Euclidean2D> getEndNode() {
1005            return endNode;
1006        }
1007
1008        /** Get the previous segment.
1009         * @return previous segment
1010         */
1011        public ConnectableSegment getPrevious() {
1012            return previous;
1013        }
1014
1015        /** Set the previous segment.
1016         * @param previous previous segment
1017         */
1018        public void setPrevious(final ConnectableSegment previous) {
1019            this.previous = previous;
1020        }
1021
1022        /** Get the next segment.
1023         * @return next segment
1024         */
1025        public ConnectableSegment getNext() {
1026            return next;
1027        }
1028
1029        /** Set the next segment.
1030         * @param next previous segment
1031         */
1032        public void setNext(final ConnectableSegment next) {
1033            this.next = next;
1034        }
1035
1036        /** Set the processed flag.
1037         * @param processed processed flag to set
1038         */
1039        public void setProcessed(final boolean processed) {
1040            this.processed = processed;
1041        }
1042
1043        /** Check if the segment has been processed.
1044         * @return true if the segment has been processed
1045         */
1046        public boolean isProcessed() {
1047            return processed;
1048        }
1049
1050    }
1051
1052    /** Visitor building segments. */
1053    private static class SegmentsBuilder implements BSPTreeVisitor<Euclidean2D> {
1054
1055        /** Tolerance for close nodes connection. */
1056        private final double tolerance;
1057
1058        /** Built segments. */
1059        private final List<ConnectableSegment> segments;
1060
1061        /** Simple constructor.
1062         * @param tolerance tolerance for close nodes connection
1063         */
1064        public SegmentsBuilder(final double tolerance) {
1065            this.tolerance = tolerance;
1066            this.segments  = new ArrayList<ConnectableSegment>();
1067        }
1068
1069        /** {@inheritDoc} */
1070        public Order visitOrder(final BSPTree<Euclidean2D> node) {
1071            return Order.MINUS_SUB_PLUS;
1072        }
1073
1074        /** {@inheritDoc} */
1075        public void visitInternalNode(final BSPTree<Euclidean2D> node) {
1076            @SuppressWarnings("unchecked")
1077            final BoundaryAttribute<Euclidean2D> attribute = (BoundaryAttribute<Euclidean2D>) node.getAttribute();
1078            final Iterable<BSPTree<Euclidean2D>> splitters = attribute.getSplitters();
1079            if (attribute.getPlusOutside() != null) {
1080                addContribution(attribute.getPlusOutside(), node, splitters, false);
1081            }
1082            if (attribute.getPlusInside() != null) {
1083                addContribution(attribute.getPlusInside(), node, splitters, true);
1084            }
1085        }
1086
1087        /** {@inheritDoc} */
1088        public void visitLeafNode(final BSPTree<Euclidean2D> node) {
1089        }
1090
1091        /** Add the contribution of a boundary facet.
1092         * @param sub boundary facet
1093         * @param node node containing segment
1094         * @param splitters splitters for the boundary facet
1095         * @param reversed if true, the facet has the inside on its plus side
1096         */
1097        private void addContribution(final SubHyperplane<Euclidean2D> sub,
1098                                     final BSPTree<Euclidean2D> node,
1099                                     final Iterable<BSPTree<Euclidean2D>> splitters,
1100                                     final boolean reversed) {
1101            @SuppressWarnings("unchecked")
1102            final AbstractSubHyperplane<Euclidean2D, Euclidean1D> absSub =
1103                (AbstractSubHyperplane<Euclidean2D, Euclidean1D>) sub;
1104            final Line line      = (Line) sub.getHyperplane();
1105            final List<Interval> intervals = ((IntervalsSet) absSub.getRemainingRegion()).asList();
1106            for (final Interval i : intervals) {
1107
1108                // find the 2D points
1109                final Vector2D startV = Double.isInfinite(i.getInf()) ?
1110                                        null : (Vector2D) line.toSpace((Point<Euclidean1D>) new Vector1D(i.getInf()));
1111                final Vector2D endV   = Double.isInfinite(i.getSup()) ?
1112                                        null : (Vector2D) line.toSpace((Point<Euclidean1D>) new Vector1D(i.getSup()));
1113
1114                // recover the connectivity information
1115                final BSPTree<Euclidean2D> startN = selectClosest(startV, splitters);
1116                final BSPTree<Euclidean2D> endN   = selectClosest(endV, splitters);
1117
1118                if (reversed) {
1119                    segments.add(new ConnectableSegment(endV, startV, line.getReverse(),
1120                                                        node, endN, startN));
1121                } else {
1122                    segments.add(new ConnectableSegment(startV, endV, line,
1123                                                        node, startN, endN));
1124                }
1125
1126            }
1127        }
1128
1129        /** Select the node whose cut sub-hyperplane is closest to specified point.
1130         * @param point reference point
1131         * @param candidates candidate nodes
1132         * @return node closest to point, or null if no node is closer than tolerance
1133         */
1134        private BSPTree<Euclidean2D> selectClosest(final Vector2D point, final Iterable<BSPTree<Euclidean2D>> candidates) {
1135
1136            BSPTree<Euclidean2D> selected = null;
1137            double min = Double.POSITIVE_INFINITY;
1138
1139            for (final BSPTree<Euclidean2D> node : candidates) {
1140                final double distance = FastMath.abs(node.getCut().getHyperplane().getOffset(point));
1141                if (distance < min) {
1142                    selected = node;
1143                    min      = distance;
1144                }
1145            }
1146
1147            return min <= tolerance ? selected : null;
1148
1149        }
1150
1151        /** Get the segments.
1152         * @return built segments
1153         */
1154        public List<ConnectableSegment> getSegments() {
1155            return segments;
1156        }
1157
1158    }
1159
1160}