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