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