Monday, July 27, 2015

Graph-based circle packing

The previous two posts showed examples of a simple circle packing algorithm using the packcircles package (available from CRAN and GitHub). The algorithm involved iterative pair-repulsion to jiggle the circles until (hopefully) a non-overlapping arrangement emerged. In this post we'll look an alternative approach.

An algorithm to find an arrangement of circles satisfying a prior specification of circle sizes and tangencies was described by Collins and Stephenson in their 2003 paper in Computation Geometry Theory and Applications. A version of their algorithm was implemented in Python by David Eppstein as part of his PADS library (see CirclePack.py). I've now ported David's version to R/Rcpp and included it in the packcircles package.

In the figure below, the graph on the left represents the desired pattern of circle tangencies: e.g. circle 7 should touch all of, and only, circles 1, 2, 6 and 8. Circles 5, 7, 8 and 9 are internal, while the remaining circles are external. The plot on the right shows an arrangement of circles which conforms to the input graph.

The packcircles package provides an experimental version of the algorithm via the function circleGraphLayout (see also the underlying C++ sources on GitHub). To use it, first encode the graph of tangencies as a list of vectors:

library(packcircles)
library(ggplot2)

# List of tangencies. Vector elements are circle IDs.
# The first element in each vector is an internal circle
# and the subsequent elements are its neighbours.
internal <- list(
  c(9, 4, 5, 6, 10, 11),
  c(5, 4, 8, 6, 9),
  c(8, 3, 2, 7, 6, 5, 4),
  c(7, 8, 2, 1, 6)
)
Next, specify the sizes of all external circles:
# External circle radii (constant for this example)
external <- data.frame(id = c(1, 2, 3, 4, 6, 10, 11), radius=10.0)
Note that the sizes of the internal circles will be derived by the algorithm.

Now, pass these two objects to the circleGraphLayout function which will attempt to find a corresponding arrangement of circles which can then be plotted:

# Search for a layout. The returned value is a four-column 
# data.frame: id, x, y, radius.
layout <- circleGraphLayout(internal, external)

# Generate circle vertices from the layout for plotting
plotdat <- circlePlotData(layout, xyr.cols = 2:4, id.col = 1)

# Draw the circles annotated with their IDs.
ggplot() +
  geom_polygon(data=plotdat, aes(x, y, group=id), fill=NA, colour="brown") +
  geom_text(data=layout, aes(x, y, label=id)) +
  coord_equal() +
  theme_bw()
The resulting plot should look similar to the arrangement in the figure above.

Beware: the present implementation is a bit brittle. In particular, you have to specify all circle tangencies in a consistent order, either clockwise or anti-clockwise, otherwise the returned layout will either be incomplete, or contain overlapping circles.

The functions in the packcircles package are very basic but should suffice for simple applications. For some much more impressive circle work, both aesthetically and mathematically, visit Ken Stephenson's page, and take a look at this blog post by Danny Calegari.

Thursday, May 28, 2015

Static and moving circles

After the previous post on the packcircles package for R someone suggested it would be useful to be able to fix the position of selected circles. As a first attempt, I've added an optional weights argument to the circleLayout function. Weights can be in the range 0-1 inclusive, where a weight of 0 prevents a circle from moving, while a weight of 1 allows full movement. The updated code is at GitHub.

Here's an example where the largest of a set of initially overlapping circles is fixed in place:

And here is the code for the example:

library(packcircles)
library(ggplot2)
library(gridExtra)

# Generate some random overlapping circles
ncircles <- 200
limits <- c(-50, 50)
inset <- diff(limits) / 3
rmax <- 20

xyr <- data.frame(
  x = runif(ncircles, min(limits) + inset, max(limits) - inset),
  y = runif(ncircles, min(limits) + inset, max(limits) - inset),
  r = rbeta(ncircles, 1, 10) * rmax)

# Index of the largest circle
largest.id <- which(xyr$r == max(xyr$r))

## Generate plot data for the `before` layout
dat.before <- circlePlotData(xyr)

# Add a column to the plot data for the 'before' circles
# to indicate whether a circle is static of free to move
dat.before$state <- ifelse(dat.before$id == largest.id, "static", "free")

# Run the layout algorithm with a weights vector to fix the position
# of the largest circle
wts <- rep(1.0, nrow(xyr))
wts[ largest.id ] <- 0.0

res <- circleLayout(xyr, limits, limits, maxiter = 1000, weights=wts)

# A plot function to colour circles based on the state column
doPlot <- function(dat, title)
  ggplot(dat) + 
  geom_polygon(aes(x, y, group=id, fill=state), colour="brown1") +
  scale_fill_manual(values=c("NA", "brown4")) +
  coord_equal(xlim=limits, ylim=limits) +
  theme_bw() +
  theme(axis.text=element_blank(),
        axis.ticks=element_blank(),
        axis.title=element_blank(),
        legend.position="none") +
  labs(title=title)

g.before <- doPlot(dat.before, "before")

# Generate a plot for the 'after' circles
dat.after <- circlePlotData(res$layout)
dat.after$state <- ifelse(dat.after$id == largest.id, "static", "free")

g.after <- doPlot(dat.after, "after")

grid.arrange(g.before, g.after, nrow=1)

Wednesday, April 29, 2015

Circle packing in R (again)

Back in 2010 I posted some R code for circle packing. Now, just five years later, I've ported the code to Rcpp and created a little package which you can find at GitHub.

The main function is circleLayout which takes a set of overlapping circles and tries to find a non-overlapping arrangement for them. Here's an example:


And here's the code:

# Create some random circles, positioned within the central portion
# of a bounding square, with smaller circles being more common than
# larger ones.

ncircles <- 200
limits <- c(-50, 50)
inset <- diff(limits) / 3
rmax <- 20

xyr <- data.frame(
  x = runif(ncircles, min(limits) + inset, max(limits) - inset),
  y = runif(ncircles, min(limits) + inset, max(limits) - inset),
  r = rbeta(ncircles, 1, 10) * rmax)

# Next, we use the `circleLayout` function to try to find a non-overlapping
# arrangement, allowing the circles to occupy any part of the bounding square.
# The returned value is a list with elements for the layout and the number
# of iterations performed.

library(packcircles)

res <- circleLayout(xyr, limits, limits, maxiter = 1000)
cat(res$niter, "iterations performed")

# Now draw the before and after layouts with ggplot

library(ggplot2)
library(gridExtra)

## plot data for the `before` layout
dat.before <- circlePlotData(xyr)

## plot dta for the `after` layout returned by circleLayout
dat.after <- circlePlotData(res$layout)

doPlot <- function(dat, title)
  ggplot(dat) + 
  geom_polygon(aes(x, y, group=id), colour="brown", fill="burlywood", alpha=0.3) +
  coord_equal(xlim=limits, ylim=limits) +
  theme_bw() +
  theme(axis.text=element_blank(),
        axis.ticks=element_blank(),
        axis.title=element_blank()) +
  labs(title=title)

grid.arrange(
  doPlot(dat.before, "before"),
  doPlot(dat.after, "after"),
  nrow=1)

Sunday, September 30, 2012

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. 
  #   Proceedings of the 6th International Conference in Central Europe 
  #   on Computer Graphics and Visualization. WSCG '98, p. 125-132 
  #
  # Adapted from the original Matlab code by Michael Bedward (2010)
  # michael.bedward@gmail.com
  #
  # Subsequently improved by John Minter (2012)
  # 
  # Arguments: 
  # x, y - x and y coordinates of the data points.
  #        If a single arg is provided it is assumed to be a
  #        two column matrix.
  #
  # Returns a list with the following elements: 
  #
  # coef - coefficients of the ellipse as described by the general 
  #        quadratic:  ax^2 + bxy + cy^2 + dx + ey + f = 0 
  #
  # center - center x and y
  #
  # major - major semi-axis length
  #
  # minor - minor semi-axis length
  #
  EPS <- 1.0e-8 
  dat <- xy.coords(x, y) 
  
  D1 <- cbind(dat$x * dat$x, dat$x * dat$y, dat$y * dat$y) 
  D2 <- cbind(dat$x, dat$y, 1) 
  S1 <- t(D1) %*% D1 
  S2 <- t(D1) %*% D2 
  S3 <- t(D2) %*% D2 
  T <- -solve(S3) %*% t(S2) 
  M <- S1 + S2 %*% T 
  M <- rbind(M[3,] / 2, -M[2,], M[1,] / 2) 
  evec <- eigen(M)$vec 
  cond <- 4 * evec[1,] * evec[3,] - evec[2,]^2 
  a1 <- evec[, which(cond > 0)] 
  f <- c(a1, T %*% a1) 
  names(f) <- letters[1:6] 
  
  # calculate the center and lengths of the semi-axes 
  #
  # see http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2288654/
  # J. R. Minter
  # for the center, linear algebra to the rescue
  # center is the solution to the pair of equations
  # 2ax +  by + d = 0
  # bx  + 2cy + e = 0
  # or
  # | 2a   b |   |x|   |-d|
  # |  b  2c | * |y| = |-e|
  # or
  # A x = b
  # or
  # x = Ainv b
  # or
  # x = solve(A) %*% b
  A <- matrix(c(2*f[1], f[2], f[2], 2*f[3]), nrow=2, ncol=2, byrow=T )
  b <- matrix(c(-f[4], -f[5]), nrow=2, ncol=1, byrow=T)
  soln <- solve(A) %*% b

  b2 <- f[2]^2 / 4
  
  center <- c(soln[1], soln[2]) 
  names(center) <- c("x", "y") 
  
  num  <- 2 * (f[1] * f[5]^2 / 4 + f[3] * f[4]^2 / 4 + f[6] * b2 - f[2]*f[4]*f[5]/4 - f[1]*f[3]*f[6]) 
  den1 <- (b2 - f[1]*f[3]) 
  den2 <- sqrt((f[1] - f[3])^2 + 4*b2) 
  den3 <- f[1] + f[3] 
  
  semi.axes <- sqrt(c( num / (den1 * (den2 - den3)),  num / (den1 * (-den2 - den3)) )) 
  
  # calculate the angle of rotation 
  term <- (f[1] - f[3]) / f[2] 
  angle <- atan(1 / term) / 2 
  
  list(coef=f, center = center, major = max(semi.axes), minor = min(semi.axes), angle = unname(angle)) 
}

Next here is a utility function which takes a fitted ellipse and returns a matrix of vertices for plotting:

get.ellipse <- function( fit, n=360 ) 
{
  # Calculate points on an ellipse described by 
  # the fit argument as returned by fit.ellipse 
  # 
  # n is the number of points to render 
  
  tt <- seq(0, 2*pi, length=n) 
  sa <- sin(fit$angle) 
  ca <- cos(fit$angle) 
  ct <- cos(tt) 
  st <- sin(tt) 
  
  x <- fit$center[1] + fit$maj * ct * ca - fit$min * st * sa 
  y <- fit$center[2] + fit$maj * ct * sa + fit$min * st * ca 
  
  cbind(x=x, y=y) 
}

And finally, some demo code from John:

create.test.ellipse <- function(Rx=300,         # X-radius
                                Ry=200,         # Y-radius
                                Cx=250,         # X-center
                                Cy=150,         # Y-center
                                Rotation=0.4,   # Radians
                                NoiseLevel=0.5) # Gaussian Noise level
{
  set.seed(42)
  t <- seq(0, 100, by=1)
  x <- Rx * cos(t)
  y <- Ry * sin(t)
  nx <- x*cos(Rotation)-y*sin(Rotation) + Cx
  nx <- nx + rnorm(length(t))*NoiseLevel 
  ny <- x*sin(Rotation)+y*cos(Rotation) + Cy
  ny  <- ny + rnorm(length(t))*NoiseLevel
  cbind(x=nx, y=ny)
}

X <- create.test.ellipse()
efit <- fit.ellipse(X)
e <- get.ellipse(efit)
plot(X) 
lines(e, col="red") 

print(efit)

1 Halíř R., Flusser J.: Numerically stable direct least squares fitting of ellipses. In: Proceedings of the 6th International Conference in Central Europe on Computer Graphics and Visualization. WSCG '98. CZ, Plzeň 1998, pp. 125-132. (postscript file)

Monday, August 6, 2012

IDE errors with Maven projects after upgrading to NetBeans 7.2

A recent upgrade of NetBeans from version 7.1 to 7.2 seemed to cause problems with Maven multi-module projects. Normally, a newly created class in one Maven project module will be instantly visible to other modules open in the IDE, at least for editing purposes. But after the upgrade I started to see "cannot find symbol" errors in the IDE for such new classes. Despite this, a full build of the project would succeed.

Googling for an answer drew a blank. The usual kludge fixes such as deleting the NetBeans cache (for NetBeans 7.2 on OSX this seems to have moved from ~/.netbeans to ~/Library/Application Support/Netbeans) and re-indexing my local repo didn't make any difference. Despair loomed.

Happily, I stumbled across the cause of the problem after only a moderate amount of pulling of hair and gnashing of teeth. My Maven settings.xml file included this section...


    <profiles>
        <profile>
            <id>netbeans-ide</id>
            <activation>
                <activeByDefault>true</activeByDefault>
            </activation>
            <properties>
                <netbeans.installation>/Applications/NetBeans/NetBeans 7.1.app/Contents/Resources/NetBeans</netbeans.installation>
            </properties>
        </profile>
    </profiles>

After editing the path to point to 7.2.app the IDE errors disappeared.

Tuesday, April 10, 2012

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 extensive use of META-INF/services files which need to be merged properly. The Maven Shade plugin does this nicely and I'd previously used it when creating single, executable jars for GeoTools projects. For the current case I realised all I needed to do was set up a multi-module project with one module for the application and another, inside a profile, for the combined libraries jar.

Following is the top level POM. All of the application dependencies are declared here to be inherited by the child application and libraries uber-jar POMs. Note that we define a profile "libs" to enclose the module for the libraries jar:

<?xml version="1.0" encoding="UTF-8"?>
<project xmlns="http://maven.apache.org/POM/4.0.0"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 http://maven.apache.org/maven-v4_0_0.xsd">

<modelVersion>4.0.0</modelVersion>
<groupId>org.cafeanimal</groupId>
<artifactId>trapsimulation</artifactId>
<packaging>pom</packaging>

<version>0.1-SNAPSHOT</version>
<name>Fauna trapping simulation</name>

<properties>
<geotools.version>8.0-M3</geotools.version>
</properties>

<build>
<plugins>
<plugin>
<artifactId>maven-compiler-plugin</artifactId>
<version>2.0.2</version>
<configuration>
<source>1.6</source>
<target>1.6</target>
<encoding>UTF-8</encoding>
</configuration>
</plugin>
<plugin>
<artifactId>maven-resources-plugin</artifactId>
<version>2.2</version>
<configuration>
<encoding>UTF-8</encoding>
</configuration>
</plugin>
</plugins>
</build>

<dependencies>
<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-main</artifactId>
<version>${geotools.version}</version>
</dependency>
<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-referencing</artifactId>
<version>${geotools.version}</version>
</dependency>
<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-epsg-hsql</artifactId>
<version>${geotools.version}</version>
</dependency>

<!-- Lots of other dependencies omitted here for brevity -->

</dependencies>

<modules>
<module>app</module>
</modules>

<!-- specify -Dlibs to build the lib module with its uber-jar of dependencies -->
<profiles>
<profile>
<id>libs</id>
<activation>
<property>
<name>libs</name>
</property>
</activation>

<modules>
<module>lib</module>
</modules>

</profile>
</profiles>

</project>


The Maven module for the application resides in a sub-directory (app) of the parent project directory. Since all of the dependencies were declared in the parent POM there is very little required for the application itself:

<?xml version="1.0" encoding="UTF-8"?>
<project xmlns="http://maven.apache.org/POM/4.0.0"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 http://maven.apache.org/maven-v4_0_0.xsd">

<modelVersion>4.0.0</modelVersion>

<parent>
<groupId>org.cafeanimal</groupId>
<artifactId>trapsimulation</artifactId>
<version>0.1-SNAPSHOT</version>
</parent>

<artifactId>app</artifactId>
<name>Simulation app</name>
<packaging>jar</packaging>

</project>


The Maven module for the libraries uber-jar is in another sub-directory (lib) of the parent project directory. It references the Maven Shade plugin which will merge the GeoTools META-INF/services files and combine the individual jars for the many dependencies into a single jar:

<?xml version="1.0" encoding="UTF-8"?>
<project xmlns="http://maven.apache.org/POM/4.0.0"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 http://maven.apache.org/xsd/maven-4.0.0.xsd">

<modelVersion>4.0.0</modelVersion>

<parent>
<artifactId>trapsimulation</artifactId>
<groupId>org.cafeanimal</groupId>
<version>0.1-SNAPSHOT</version>
</parent>

<artifactId>lib</artifactId>
<packaging>jar</packaging>
<name>Library uber-jar</name>

<build>
<plugins>
<plugin>
<groupId>org.apache.maven.plugins</groupId>
<artifactId>maven-shade-plugin</artifactId>
<version>1.5</version>
<executions>
<execution>
<phase>package</phase>
<goals>
<goal>shade</goal>
</goals>
<configuration>
<transformers>
<!-- This transformer combines the META-INF/services files from the individual GeoTools modules -->
<transformer implementation="org.apache.maven.plugins.shade.resource.ServicesResourceTransformer"/>
</transformers>
</configuration>
</execution>
</executions>
</plugin>
</plugins>
</build>

</project>


Now to build both the application and the libraries jar from the command line, cd to the parent project directory and enter:
mvn clean install -Dlibs

To build the application alone, simply omit the -Dlibs profile variable.

Monday, December 6, 2010

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