    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     */
    017    package org.apache.commons.math3.geometry.euclidean.twod;
    018
    019    import java.util.ArrayList;
    020    import java.util.Collection;
    021    import java.util.List;
    022
    023    import org.apache.commons.math3.exception.MathInternalError;
    024    import org.apache.commons.math3.geometry.euclidean.oned.Euclidean1D;
    025    import org.apache.commons.math3.geometry.euclidean.oned.Interval;
    026    import org.apache.commons.math3.geometry.euclidean.oned.IntervalsSet;
    027    import org.apache.commons.math3.geometry.euclidean.oned.Vector1D;
    028    import org.apache.commons.math3.geometry.partitioning.AbstractRegion;
    029    import org.apache.commons.math3.geometry.partitioning.AbstractSubHyperplane;
    030    import org.apache.commons.math3.geometry.partitioning.BSPTree;
    031    import org.apache.commons.math3.geometry.partitioning.BSPTreeVisitor;
    032    import org.apache.commons.math3.geometry.partitioning.BoundaryAttribute;
    033    import org.apache.commons.math3.geometry.partitioning.Side;
    034    import org.apache.commons.math3.geometry.partitioning.SubHyperplane;
    035    import org.apache.commons.math3.geometry.partitioning.utilities.AVLTree;
    036    import org.apache.commons.math3.geometry.partitioning.utilities.OrderedTuple;
    037    import org.apache.commons.math3.util.FastMath;
    038
    039    /** This class represents a 2D region: a set of polygons.
    040     * @version $Id: PolygonsSet.java 1422195 2012-12-15 06:45:18Z psteitz $
    041     * @since 3.0
    042     */
    043    public 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>();
    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    }