Skip to main content

Smooth as...

For a current project I needed a polygon smoothing algorithm in Java that:

  1. generated a curve passing through all polygon vertices

  2. had some way of controlling tightness of fit to the polygon

  3. fell within the domain of my very limited mathematical grasp

Luckily for me, Maxim Shemanarev had not only worked out a Bezier curve algorithm that satisfied all of the above, but he also published a beautifully clear explanation of how it worksnote 1. Maxim also provided some illustrative C code, but my Java implementation below is an independent effort (ie. all mistakes are my fault). It uses geometry objects from the JTS library.

Below is an example of the algorithm in action. The original polygon is in black. The smoothed polygon composed of fitted cubic Bezier segments is in blue. The red points and lines indicate positions of the Bezier control points calculated with Maxim's algorithm. There are a pair of control points for each vertex of the input polygon.

The tightness of fit of the smoothed boundary to the original boundary is controlled by a single parameter which varies between 0 and 1. In the figure above the tightness was set to 0. Below is the same polygon smoothed with tightness set to 0.6.

Notice that the higher tightness parameter has shifted each pair of control points closer to their respective polygon vertices.

Finally, here's the Java code. Share and enjoy...
/*
* Copyright 2010 Michael Bedward
*
* This file is part of jai-tools.
*
* jai-tools is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as
* published by the Free Software Foundation, either version 3 of the
* License, or (at your option) any later version.
*
* jai-tools is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with jai-tools. If not, see <http://www.gnu.org/licenses/>.
*/

package mxb.interpolation;

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.GeometryFactory;
import com.vividsolutions.jts.geom.LinearRing;
import com.vividsolutions.jts.geom.Polygon;


/**
* Polygon smoothing by interpolation with cubic Bazier curves.
* This code implements an algorithm which was devised by Maxim Shemanarev
* and described by him at:
* <a href="http://www.antigrain.com/research/bezier_interpolation/index.html">
* http://www.antigrain.com/research/bezier_interpolation/index.html</a>.
*
* Note: the code here is <b>not</b> that written by Maxim to accompany his
* algorithm description. Rather, it is an original implementation and any
* errors are my fault.
*/
public class PolygonSmoother {
private GeometryFactory geomFactory;

/**
* Calculates a pair of Bezier control points for each vertex in an
* array of {@code Coordinates}.
*
* @param coords vertex coordinates
* @param N number of coordinates in {@coords} to use
* @param alpha tightness of fit
*
* @return array of {@code Coordinates} for positions of control points
*/
private Coordinate[][] getControlPoints(Coordinate[] coords, int N, double alpha) {
if (alpha < 0.0 || alpha > 1.0) {
throw new IllegalArgumentException("alpha must be a value between 0 and 1 inclusive");
}

Coordinate[][] ctrl = new Coordinate[N][2];

Coordinate[] v = new Coordinate[3];

Coordinate[] mid = new Coordinate[2];
mid[0] = new Coordinate();
mid[1] = new Coordinate();

Coordinate anchor = new Coordinate();
double[] vdist = new double[2];
double mdist;

v[1] = coords[N - 1];
v[2] = coords[0];
mid[1].x = (v[1].x + v[2].x) / 2.0;
mid[1].y = (v[1].y + v[2].y) / 2.0;
vdist[1] = v[1].distance(v[2]);

for (int i = 0; i < N; i++) {
v[0] = v[1];
v[1] = v[2];
v[2] = coords[(i + 1) % N];

mid[0].x = mid[1].x;
mid[0].y = mid[1].y;
mid[1].x = (v[1].x + v[2].x) / 2.0;
mid[1].y = (v[1].y + v[2].y) / 2.0;

vdist[0] = vdist[1];
vdist[1] = v[1].distance(v[2]);

double p = vdist[0] / (vdist[0] + vdist[1]);
anchor.x = mid[0].x + p * (mid[1].x - mid[0].x);
anchor.y = mid[0].y + p * (mid[1].y - mid[0].y);

double xdelta = anchor.x - v[1].x;
double ydelta = anchor.y - v[1].y;

ctrl[i][0] = new Coordinate(
alpha*(v[1].x - mid[0].x + xdelta) + mid[0].x - xdelta,
alpha*(v[1].y - mid[0].y + ydelta) + mid[0].y - ydelta);

ctrl[i][1] = new Coordinate(
alpha*(v[1].x - mid[1].x + xdelta) + mid[1].x - xdelta,
alpha*(v[1].y - mid[1].y + ydelta) + mid[1].y - ydelta);
}

return ctrl;
}


/**
* Calculates vertices along a cubic Bazier curve given start point, end point
* and two control points.
*
* @param start start position
* @param end end position
* @param ctrl1 first control point
* @param ctrl2 second control point
* @param nv number of vertices including the start and end points
*
* @return vertices along the Bezier curve
*/
private Coordinate[] cubicBezier(final Coordinate start, final Coordinate end,
final Coordinate ctrl1, final Coordinate ctrl2, final int nv) {

final Coordinate[] curve = new Coordinate[nv];

final Coordinate[] buf = new Coordinate[3];
for (int i = 0; i < buf.length; i++) {
buf[i] = new Coordinate();
}

curve[0] = new Coordinate(start);
curve[nv - 1] = new Coordinate(end);

for (int i = 1; i < nv-1; i++) {
double t = (double) i / (nv - 1);
double tc = 1.0 - t;

double t0 = tc*tc*tc;
double t1 = 3.0*tc*tc*t;
double t2 = 3.0*tc*t*t;
double t3 = t*t*t;
double tsum = t0 + t1 + t2 + t3;

Coordinate c = new Coordinate();

c.x = t0*start.x + t1*ctrl1.x + t2*ctrl2.x + t3*end.x;
c.x /= tsum;
c.y = t0*start.y + t1*ctrl1.y + t2*ctrl2.y + t3*end.y;
c.y /= tsum;

curve[i] = c;
}

return curve;
}

/**
* Creates a new {@code Polygon} whose exterior shell is a smoothed
* version of the input {@code Polygon}.
* <p>
* Note: this method presently ignores holes.
*
* @param p the input {@code Polygon}
*
* @param alpha a value between 0 and 1 (inclusive) specifying the tightness
* of fit of the smoothed boundary (0 is loose)
*
* @param pointsPerSegment the number of vertices to use to represent each
* Bezier curve derived from an edge of the input {@code Polygon}
*
* @return the smoothed {@code Polygon}
*/
public Polygon smooth(Polygon p, double alpha, int pointsPerSegment) {

if (geomFactory == null) {
geomFactory = new GeometryFactory();
}

Coordinate[] coords = p.getExteriorRing().getCoordinates();
int Nvertices = coords.length - 1; // first coord == last coord

Coordinate[][] controlPoints = getControlPoints(coords, Nvertices, alpha);

Coordinate[] smoothCoords = new Coordinate[Nvertices * pointsPerSegment];
for (int i = 0, k = 0; i < Nvertices; i++) {
int next = (i + 1) % Nvertices;

Coordinate[] segment = cubicBezier(
coords[i], coords[next],
controlPoints[i][1], controlPoints[next][0],
pointsPerSegment);

for (int j = 0; j < pointsPerSegment; j++, k++) {
smoothCoords[k] = segment[j];
}
}

LinearRing shell = geomFactory.createLinearRing(smoothCoords);
return geomFactory.createPolygon(shell, null);
}
}


1 Maxim Shemanarev's algorithm http://www.antigrain.com/research/bezier_interpolation/index.html

Comments

  1. do you have any example of the usage of your code? I tried to use it with multipolygon(call smoooth() for each polygon in mutipolygon) - the result is rather strange:) No smoothing, some new straight lines. Also, I posted a question to mailing list with my current issue( http://osgeo-org.1560.n6.nabble.com/multipolygon-smooth-td4994497.html ). Thank you for any help.

    ReplyDelete
  2. Very late reply here... I think I've already replied to your post on the GeoTools user list.

    You can see MultiPolygon smoothing implemented with a version of this code in the GeoTools library class JTS:

    https://github.com/geotools/geotools/blob/8.1/modules/library/api/src/main/java/org/geotools/geometry/jts/JTS.java#L971

    (look for the method smoothMultiPolygon)

    Also Java classes that implement this algorithm for polygons and line strings, in a fairly user friendly way, are included in the JAITools library:
    http://jaitools.googlecode.com/git/utils/src/main/java/org/jaitools/jts/


    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…