Another clip round the ear

In the last episode I presented some code to triangulate a polygon using the ear-clipping algorithm. In this instalment the code has been extended to deal with polygons that have holes. Here are two example results...


No, they're not the prettiest triangulations ever seen, but they could be improved with some post-processing.

Here is the Java code used to generate them. Note, you'll need the JTS library if you want to run this. Warning: this code isn't pretty, even by my standards, but it might be useful if you want to play with the algorithm.

import com.vividsolutions.jts.algorithm.CGAlgorithms;
import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Envelope;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryFactory;
import com.vividsolutions.jts.geom.LineString;
import com.vividsolutions.jts.geom.Polygon;
import com.vividsolutions.jts.io.WKTReader;
import java.awt.Color;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.Comparator;
import java.util.List;

/**
 * Demonstrates brute force approach to the ear clipping algorithm
 * to triangulate a polygon. 
 * 
 * This version attempts a general approach to holes.
 *
 * @author Michael Bedward
 */
public class EarClipping4 {
    private static final double EPS = 1.0E-4;

    private final GeometryFactory gf;
    private final Polygon inputPolygon;
    private Geometry ears;

    private int curRingIndex;
    private List<Coordinate> shellCoords;
    private List<Integer> ringIndices;

    /**
     * Demonstrate the ear-clipping algorithm
     * @throws Exception
     */
    public static void main(String[] args) throws Exception {
        WKTReader reader = new WKTReader();
        Polygon poly = null;

        /*
         * Square with one hole that can't reach any exterior vertex
         *
        poly = (Polygon) reader.read(
                "POLYGON ((0 0, 0 200, 200 200, 200 0, 0 0), " +
                "(20 20, 20 40, 40 40, 40 20, 20 20), " +
                "(20 160, 20 180, 40 180, 40 160, 20 160), " +
                "(160 160, 160 180, 180 180, 180 160, 160 160), " +
                "(160 20, 160 40, 180 40, 180 20, 160 20), " +
                "(90 90, 90 110, 110 110, 110 90, 90 90) )");
         * 
         */

        /*
         * Interlinked U-shaped holes
         *
        poly = (Polygon) reader.read(
                "POLYGON ((0 0, 0 200, 250 200, 250 0, 0 0), " +
                "(40 40, 40 170, 60 170, 60 60, 140 60, 140 120, 160 120, 160 40, 40 40), " +
                "(90 80, 90 170, 210 170, 210 40, 190 40, 190 150, 110 150, 110 80, 90 80))" );

        EarClipping4 clipper = new EarClipping4(poly);
        Geometry ears  = clipper.getResult();

        for (int i = 0; i < ears.getNumGeometries(); i++) {
            System.out.println(ears.getGeometryN(i));
        }
    }

    /**
     * Constructor
     *
     * @param inputPolygon the input polygon
     */
    public EarClipping4(Polygon inputPolygon) {
        gf = new GeometryFactory();
        this.inputPolygon = inputPolygon;
    }

    /**
     * Get the result triangular polygons.
     *
     * @return triangles as a GeometryCollection
     */
    public Geometry getResult() {
        if (ears == null) {
            ears = triangulate();
        }

        return ears;
    }

    /**
     * Perform the triangulation
     *
     * @return GeometryCollection of triangular polygons
     */
    private Geometry triangulate() {
        shellCoords = new ArrayList<Coordinate>();
        ringIndices = new ArrayList<Integer>();
        List<Polygon> earList = new ArrayList<Polygon>();
        connectHoles();

        removeColinearVertices();
        int N = shellCoords.size() - 1;

        boolean finished = false;
        boolean found = false;
        int k0 = 0;
        do {
            int k1 = (k0 + 1) % N;
            int k2 = (k1 + 1) % N;
            LineString ls = gf.createLineString(new Coordinate[] {shellCoords.get(k0), shellCoords.get(k2)});
            
            found = false;
            if (inputPolygon.covers(ls)) {
                Polygon ear = gf.createPolygon(gf.createLinearRing(
                        new Coordinate[]{
                            shellCoords.get(k0),
                            shellCoords.get(k1),
                            shellCoords.get(k2),
                            shellCoords.get(k0)}),
                        null);

                if (inputPolygon.covers(ear)) {
                    found = true;
                    for (Polygon p : earList) {
                        if (ear.overlaps(p)) {
                            // overlap between candidate ear and previous ear
                            found = false;
                            break;
                        }
                    }

                    if (found) {
                        earList.add(ear);
                        shellCoords.remove(k1);
                        ringIndices.remove(k1);
                        removeColinearVertices();
                        N = shellCoords.size() - 1;
                        k0 = 0;
                        if (N < 3) {
                            finished = true;
                        }
                    }
                }
            }

            if (!found) {
                k0++;

                if (k0 == N) {
                    finished = true;
                }
            }

        } while (!finished);

        return gf.createGeometryCollection(earList.toArray(new Geometry[0]));
    }

    /**
     * Remove co-linear vertices. TopologyPreservingSimplifier could be
     * used for this but that seems like over-kill.
     */
    private void removeColinearVertices() {
        int N = shellCoords.size() - 1;
        List<Coordinate> coordList = new ArrayList<Coordinate>();

        for (int j = 1; j <= N; ) {
            int i = (j - 1) % N;
            int k = (j + 1) % N;
            if (ringIndices.get(i) == ringIndices.get(j) &&
                ringIndices.get(k) == ringIndices.get(j) &&
                CGAlgorithms.isOnLine(shellCoords.get(j),
                    new Coordinate[]{shellCoords.get(i), shellCoords.get(k)})) {
                Coordinate c = shellCoords.remove(j);
                ringIndices.remove(j);
                N-- ;
            } else {
                j++ ;
            }
        }
    }

    /**
     * Connect any holes in the input polygon to the exterior ring to
     * form a self-intersecting shell. Holes are added from the lowest
     * upwards. As the resulting shell develops, a hole might be added
     * to what was originally another hole.
     *
     * @return a new polygon with holes (if any) connected to the exterior
     *         ring
     */
    private List<Coordinate> connectHoles() {
        // defensively copy the input polygon
        Polygon poly = (Polygon) inputPolygon.clone();
        poly.normalize();

        List<Geometry> orderedHoles = getOrderedHoles(poly);

        Coordinate[] coords = poly.getExteriorRing().getCoordinates();
        shellCoords.addAll(Arrays.asList(coords));

        curRingIndex = 0;
        for (int i = 0; i < shellCoords.size(); i++) {
            ringIndices.add(curRingIndex);
        }

        for (int i = 0; i < orderedHoles.size(); i++) {
            curRingIndex++ ;
            joinHoleToShell(orderedHoles.get(i));
        }

        return shellCoords;
    }

    /**
     * Returns a list of holes in the input polygon (if any) ordered
     * by y coordinate with ties broken using x coordinate.
     *
     * @param poly input polygon
     * @return a list of Geometry objects representing the ordered holes
     *         (may be empty)
     */
    private List<Geometry> getOrderedHoles(final Polygon poly) {
        List<Geometry> holes = new ArrayList<Geometry>();
        List<IndexedEnvelope> bounds = new ArrayList<IndexedEnvelope>();

        if (poly.getNumInteriorRing() > 0) {
            for (int i = 0; i < poly.getNumInteriorRing(); i++) {
                bounds.add( new IndexedEnvelope(i, poly.getInteriorRingN(i).getEnvelopeInternal()) );
            }

            Collections.sort(bounds, new IndexedEnvelopeComparator());

            for (int i = 0; i < bounds.size(); i++) {
                holes.add(poly.getInteriorRingN(bounds.get(i).index));
            }
        }

        return holes;
    }

    /**
     * Join a given hole to the current shell. The hole coordinates are
     * inserted into the list of shell coordinates.
     *
     * @param hole the hole to join
     */
    private void joinHoleToShell(Geometry hole) {
        double minD2 = Double.MAX_VALUE;
        int shellVertexIndex = -1;

        final int Ns = shellCoords.size() - 1;

        final int holeVertexIndex = getLowestVertex(hole);
        final Coordinate[] holeCoords = hole.getCoordinates();

        final Coordinate ch = holeCoords[holeVertexIndex];
        List<IndexedDouble> distanceList = new ArrayList<IndexedDouble>();

        /*
         * Note: it's important to scan the shell vertices in reverse so
         * that if a hole ends up being joined to what was originally
         * another hole, the previous hole's coordinates appear in the shell
         * before the new hole's coordinates (otherwise the triangulation
         * algorithm tends to get stuck).
         */
        for (int i = Ns - 1; i >= 0; i--) {
            Coordinate cs = shellCoords.get(i);
            double d2 = (ch.x - cs.x) * (ch.x - cs.x) + (ch.y - cs.y) * (ch.y - cs.y);
            if (d2 < minD2) {
                minD2 = d2;
                shellVertexIndex = i;
            }
            distanceList.add(new IndexedDouble(i, d2));
        }
        
        /*
         * Try a quick join: if the closest shell vertex is reachable without
         * crossing any holes.
         */
        LineString join = gf.createLineString(new Coordinate[]{ch, shellCoords.get(shellVertexIndex)});
        if (inputPolygon.covers(join)) {
            doJoinHole(shellVertexIndex, holeCoords, holeVertexIndex);
            return;
        }

        /*
         * Quick join didn't work. Sort the shell coords on distance to the
         * hole vertex nnd choose the closest reachable one.
         */
        Collections.sort(distanceList, new IndexedDoubleComparator());
        for (int i = 1; i < distanceList.size(); i++) {
            join = gf.createLineString(new Coordinate[] {ch, shellCoords.get(distanceList.get(i).index)});
            if (inputPolygon.covers(join)) {
                shellVertexIndex = distanceList.get(i).index;
                doJoinHole(shellVertexIndex, holeCoords, holeVertexIndex);
                return;
            }
        }

        throw new IllegalStateException("Failed to join hole to shell");
    }

    /**
     * Helper method for joinHoleToShell. Insert the hole coordinates into
     * the shell coordinate list.
     *
     * @param shellCoords list of current shell coordinates
     * @param shellVertexIndex insertion point in the shell coordinate list
     * @param holeCoords array of hole coordinates
     * @param holeVertexIndex attachment point of hole
     */
    private void doJoinHole(int shellVertexIndex, Coordinate[] holeCoords, int holeVertexIndex) {
        List<Coordinate> newCoords = new ArrayList<Coordinate>();
        List<Integer> newRingIndices = new ArrayList<Integer>();

        newCoords.add(new Coordinate(shellCoords.get(shellVertexIndex)));
        newRingIndices.add(ringIndices.get(shellVertexIndex));

        final int N = holeCoords.length - 1;
        int i = holeVertexIndex;
        do {
            newCoords.add(new Coordinate(holeCoords[i]));
            newRingIndices.add(curRingIndex);
            i = (i + 1) % N;
        } while (i != holeVertexIndex);

        newCoords.add(new Coordinate(holeCoords[holeVertexIndex]));
        newRingIndices.add(curRingIndex);

        shellCoords.addAll(shellVertexIndex, newCoords);
        ringIndices.addAll(shellVertexIndex, newRingIndices);
    }

    /**
     * Return the index of the lowest vertex
     *
     * @param geom input geometry
     * @return index of the first vertex found at lowest point
     *         of the geometry
     */
    private int getLowestVertex(Geometry geom) {
        Coordinate[] coords = geom.getCoordinates();
        double minY = geom.getEnvelopeInternal().getMinY();
        for (int i = 0; i < coords.length; i++) {
            if (Math.abs(coords[i].y - minY) < EPS) {
                return i;
            }
        }

        throw new IllegalStateException("Failed to find lowest vertex");
    }

    private static class IndexedEnvelope {
        int index;
        Envelope envelope;

        public IndexedEnvelope(int i, Envelope env) { index = i; envelope = env; }
    }

    private static class IndexedEnvelopeComparator implements Comparator<IndexedEnvelope> {
        public int compare(IndexedEnvelope o1, IndexedEnvelope o2) {
            double delta = o1.envelope.getMinY() - o2.envelope.getMinY();
            if (Math.abs(delta) < EPS) {
                delta = o1.envelope.getMinX() - o2.envelope.getMinX();
                if (Math.abs(delta) < EPS) {
                    return 0;
                }
            }
            return (delta > 0 ? 1 : -1);
        }
    }

    private static class IndexedDouble {
        int index;
        double value;

        public IndexedDouble(int i, double v) { index = i; value = v; }
    }

    private static class IndexedDoubleComparator implements Comparator<IndexedDouble> {
        public int compare(IndexedDouble o1, IndexedDouble o2) {
            double delta = o1.value - o2.value;
            if (Math.abs(delta) < EPS) {
                    return 0;
            }
            return (delta > 0 ? 1 : -1);
        }
    }

}