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