Skip to main content

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);
}
}

}

Comments

  1. HI,Michael,,,,,,,,,,your on packing is great ,,,,,,,i am also working in field of circle packing but in matlab.

    If u have supported matlab material then please mail me at

    er.sandeepkumarrao@gmail.com

    ReplyDelete

Post a Comment

Popular posts from this blog

Fitting an ellipse to point data

Some time ago I wrote an R function to fit an ellipse to point data, using an algorithm developed by Radim Halíř and Jan Flusser1 in Matlab, and posted it to the r-help list. The implementation was a bit hacky, returning odd results for some data. A couple of days ago, an email arrived from John Minter asking for a pointer to the original code. I replied with a link and mentioned that I'd be interested to know if John made any improvements to the code. About ten minutes later, John emailed again with a much improved version ! Not only is it more reliable, but also more efficient. So with many thanks to John, here is the improved code: fit.ellipse <- function (x, y = NULL) { # from: # http://r.789695.n4.nabble.com/Fitting-a-half-ellipse-curve-tp2719037p2720560.html # # Least squares fitting of an ellipse to point data # using the algorithm described in: # Radim Halir & Jan Flusser. 1998. # Numerically stable direct least squares fitting of ellipses. …

Build an application plus a separate library uber-jar using Maven

I've been working on a small Java application with a colleague to simulate animal movements and look at the efficiency of different survey methods. It uses the GeoTools library to support map projections and shapefile output. GeoTools is great but comes at a cost in terms of size: the jar for our little application alone is less than 50kb but bundling it with GeoTools and its dependencies blows that out to 20Mb.

The application code has been changing on a daily basis as we explore ideas, add features and fix bugs. Working with my colleague at a distance, over a fairly feeble internet connection, I wanted to package the static libraries and the volatile application into separate jars so that he only needed to download the former once (another option would have been for my colleague to set up a local Maven repository but for various reasons this was impractical).

A slight complication with bundling GeoTools modules into a single jar (aka uber-jar) is that individual modules make ext…

Circle packing with R

To visualize the results of a simulation model of woodland trees within R, I needed an algorithm that could arrange a large number of circles within a rectangle such that no two circles overlapped by more than a specified amount. A colleague had approached this problem earlier by sorting the circles in order of descending size, then randomly dropping each one into the rectangle repeatedly until it landed in a position with acceptable overlap.

I suspected a faster and more robust algorithm could be constructed using some kind of "jiggling the circles" approach. Luckily for me, I discovered that Sean McCullough had written a really nice example of circles packing into a cluster using the Processing language. Sean's program is based on an iterative pair-repulsion algorithm in which overlapping circles move away from each other. Based on this, and modifying the algorithm a little, I came up with an R function to produce constrained random layouts of a given set of circles. He…