Last Resort Software

Smooth as...

2010-12-06T18:21:00.000-08:00

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

generated a curve passing through all polygon vertices

had some way of controlling tightness of fit to the polygon

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 works^{note 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);
  }
}

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

Monte Carlo testing of classification groups

2010-09-01T00:23:00.000-07:00

This is another article on the theme of defining groups in a hierarchical classification. A previous article described homogeneity analysis to visualize how any well any number of groups, defined at the same level accounts for the variability in the dataset, as measured by within-group pairwise distances. Here we will look at testing whether splitting a particular group into two daughter groups is supported by the data.

Suppose we have a group G containing N objects which we are considering splitting into two groups G1 and G2 consisting of N1 and N2 objects respectively. A Monte Carlo approach to testing whether this is worth doing is as follows:

1. Calculate the mean within-group pairwise distance that would result from defining G1 and G2.

2. Randomly allocate objects to two groups of sizes N1 and N2 and calculate the resulting mean within-group pairwise distance.

3. Repeat step 2 a large number of times (e.g. 999) to create a vector of Monte Carlo means. Add the actual mean from step 1 to this vector.

4. Calculate an empirical probability from the position of the actual mean in the sorted vector from step 3.

If the groups G1 and G2 are meaningful, the probability from step 4 should be low, ie. the actual mean should be lower than most or all of the Monte Carlo means.

Below is the R code for this algorithm. It takes as input a distance matrix; a vector of group IDs for all objects; the IDs of the two groups being tested and the number of Monte Carlo means to sample.



testGroups <- function (d, grp, g1=1, g2=2, nsamples=999)
{
# Performs a Monte-Carlo test of the informativeness of two groups of
# objects by comparing the actual mean of within-group pairwise distances
# to a vector of means derived by randomly allocating objects to the
# two groups. If the two groups are informative, the actual mean should
# be less than most or all of the Monte-Carlo means.
#
# Arguments:
# d - distance matrix
# grp - vector of group ids for all objects in the distance matrix
#       (ids of objects outside the two groups being considered can
#        be set to arbitrary values)
# g1 - id of first group to be tested
# g2 - id of second group to be tested
# nsamples - number of random groupings to generate (if possible)
#
# The function attempts to generate unique allocations using the following rules:
#
# 1. If nsamples is greater than choose(N1+N2, N1) it is reduced accordingly 
#    and a warning message issued. All unique allocations are then 
#    generated systematically.
#
# 2. If nsamples is greater than choose(N1+N2, N1) / 20 allocations are generated
#    systematically to guarantee uniqueness. 
#
# 3. If nsamples is less than choose(N1+N2, N1) / 20, allocations are
#    generated randomly with no check for uniqueness to reduce running
#    time.
#
# The function invisibly returns a list with the following elements:
# actual.mean - the mean within-group pairwise distance for actual groupings
# prob - probability of observing a mean less than or equal to the actual
#        mean calculated by comparing the actual value to a vector of 
#        means from random allocations plus the actual mean
# mc.mean - mean within-group pairwise distance for each random allocation
#
# Michael Bedward, August 2010

  # Constants
  SYSTEMATIC <- 1
  RANDOM <- 2

  Nobj <- length(grp)
  if (length(d) != Nobj*(Nobj-1)/2) {
    stop("Number of objects differs between dist object and group vector")
  }
  
  N1 <- sum(grp == g1)
  N2 <- sum(grp == g2)
  if (N1 == 1 & N2 == 1) {
    stop("Can't test split between two singleton groups")
  }

  objects <- which(grp == g1 | grp == g2)
  
  Nc <- choose(N1+N2, N1)
  
  if (nsamples >= Nc/20) {
    if (nsamples > Nc) {
      warning(paste("Requested number of samples (", nsamples, ") exceeds possible rearrangements (", Nc, ")", sep=""))
      nsamples <- Nc
    }
    sample.method <- SYSTEMATIC
  } else {
    sample.method <- RANDOM
  }

  # Helper function: calculates total pairwise distance between
  # objects within a group. Returns total and number of 
  # pairwise distances
  calc.grp.dist <- function(grp.obj) {
    if (length(grp.obj) == 1) {
      return( c(0, 0) )
    }

    pairwise <- combn(sort(grp.obj), 2)
    di <- Nobj * (pairwise[1,] - 1) - pairwise[1,]*(pairwise[1,]-1)/2 + pairwise[2,]-pairwise[1,]
    c(sum(d[di]), length(di)) 
  }

  # Helper function: calculates mean within-group distance.
  # Arg indices holds the indices of elements in the objects vector for
  # one of the groups.
  calc.mean <- function(indices) {
    g1dist <- calc.grp.dist(objects[indices])
    g2dist <- calc.grp.dist(objects[-indices])
    (g1dist[1] + g2dist[1]) / (g1dist[2] + g2dist[2])
  }

  samples <- numeric(nsamples+1)

  # actual within-group mean distance
  samples[1] <- calc.mean(which(grp == g1))

  # randomly allocate objects to groups and calculate means
  Nmin <- min(N1, N2)
  k <- 2
  if (sample.method == SYSTEMATIC) {
    combns <- combn(N1+N2, Nmin)[, sample(Nc, nsamples)]
  } else {  # sample.method == RANDOM
    combns <- replicate(nsamples, sample(N1+N2, Nmin))
  }

  apply( combns, 2, function(x){samples[k] <<- calc.mean(x); k <<- k+1} )
  prob <- ecdf(samples)(samples[1])
  print(paste("Prob =", prob, "from", nsamples, "samples plus actual value"))
  invisible(list(actual.mean=samples[1], prob=prob, mc.means=samples[-1]))
}

Fun with the proto package: building an MCMC sampler for Bayesian regression

2010-08-12T18:35:00.000-07:00

The proto package is my latest favourite R goodie. It brings prototype-based programming to the R language - a style of programming that lets you do many of the things you can do with classes, but with a lot less up-front work. Louis Kates and Thomas Petzoldt provide an excellent introduction to using proto in the package vignette.

As a learning exercise I concocted the example below involving Bayesian logistic regression. It was inspired by an article on Matt Shotwell's blog about using R environments rather than lists to store the state of a Markov Chain Monte Carlo sampler. Here I use proto to create a parent class-like object (or trait in proto-ese) to contain the regression functions and create child objects to hold both data and results for individual analyses.

First here's an example session...


# Make up some data with a continuous predictor and binary response
nrec <- 500
x <- rnorm(nrec)
y <- rbinom(nrec, 1, plogis(2 - 4*x))

# Predictor matrix with a col of 1s for intercept
pred <- matrix(c(rep(1, nrec), x), ncol=2)
colnames(pred) <- c("intercept", "X")

# Load the proto package
library(proto)

# Use the Logistic parent object to create a child object which will 
# hold the data and run the regression (the $ operator references 
# functions and data within a proto object)
lr <- Logistic$new(pred, y)
lr$run(5000, 1000)

# lr now contains both data and results
str(lr)

proto object 
 $ cov      : num [1:2, 1:2] 0.05 -0.0667 -0.0667 0.1621 
  ..- attr(*, "dimnames")=List of 2 
 $ prior.cov: num [1:2, 1:2] 100 0 0 100 
 $ prior.mu : num [1:2] 0 0 
 $ beta     : num [1:5000, 1:2] 2.09 2.09 2.09 2.21 2.21 ... 
  ..- attr(*, "dimnames")=List of 2 
 $ adapt    : num 1000 
 $ y        : num [1:500] 0 1 1 1 1 1 1 1 1 1 ... 
 $ x        : num [1:500, 1:2] 1 1 1 1 1 1 1 1 1 1 ... 
  ..- attr(*, "dimnames")=List of 2 
 parent: proto object 

# Use the Logistic summary function to tabulate and plot results
lr$summary()

From 5000 samples after 1000 iterations burning in
   intercept           X         
 Min.   :1.420   Min.   :-5.296  
 1st Qu.:1.840   1st Qu.:-3.915  
 Median :2.000   Median :-3.668  
 Mean   :1.994   Mean   :-3.693  
 3rd Qu.:2.128   3rd Qu.:-3.455  
 Max.   :2.744   Max.   :-2.437

And here's the code for the Logistic trait...


Logistic <- proto()

Logistic$new <- function(., x, y) {
  # Creates a child object to hold data and results
  #
  # x - a design matrix (ie. predictors
  # y - a binary reponse vector

  proto(., x=x, y=y)
}

Logistic$run <- function(., niter, adapt=1000) {
  # Perform the regression by running the MCMC
  # sampler
  #
  # niter - number of iterations to sample
  # adapt - number of prior iterations to run
  #         for the 'burning in' period

  require(mvtnorm)

  # Set up variables used by the sampler
  .$adapt <- adapt
  total.iter <- niter + adapt  
  .$beta <- matrix(0, nrow=total.iter, ncol=ncol(.$x))
  .$prior.mu <- rep(0, ncol(.$x))
  .$prior.cov <- diag(100, ncol(.$x))
  .$cov <- diag(ncol(.$x))

  # Run the sampler
  b <- rep(0, ncol(.$x))
  for (i in 1:total.iter) {
    b <- .$update(i, b)
    .$beta[i,] <- b
  }
  
  # Trim the results matrix to remove the burn-in
  # period
  if (.$adapt > 0) {
    .$beta <- .$beta[(.$adapt + 1):total.iter,]
  }
}

Logistic$update <- function(., it, beta.old) {
  # Perform a single iteration of the MCMC sampler using
  # Metropolis-Hastings algorithm.
  # Adapted from code by Brian Neelon published at:
  # http://www.duke.edu/~neelo003/r/
  #
  # it -       iteration number
  # beta.old - vector of coefficient values from 
  #            the previous iteration

  # Update the coefficient covariance if we are far
  # enough through the sampling
  if (.$adapt > 0 & it > 2 * .$adapt) {
    .$cov <- cov(.$beta[(it - .$adapt):(it - 1),])
  }
  
  # generate proposed new coefficient values
  beta.new <- c(beta.old + rmvnorm(1, sigma=.$cov))
  
  # calculate prior and current probabilities and log-likelihood
  if (it == 1) {
    .$..log.prior.old <- dmvnorm(beta.old, .$prior.mu, .$prior.cov, log=TRUE)
    .$..probs.old <- plogis(.$x %*% beta.old)
    .$..LL.old <- sum(log(ifelse(.$y, .$..probs.old, 1 - .$..probs.old)))
  }
  log.prior.new <- dmvnorm(beta.new, .$prior.mu, .$prior.cov, log=TRUE)
  probs.new <- plogis(.$x %*% beta.new)  
  LL.new <- sum(log(ifelse(.$y, probs.new, 1-probs.new)))
  
  # Metropolis-Hastings acceptance ratio (log scale)
  ratio <- LL.new + log.prior.new - .$..LL.old - .$..log.prior.old
  
  if (log(runif(1)) < ratio) {
   .$..log.prior.old <- log.prior.new
   .$..probs.old <- probs.new
   .$..LL.old <- LL.new
    return(beta.new)
  } else {
    return(beta.old)
  }
}

Logistic$summary <- function(., show.plot=TRUE) {
  # Summarize the results

  cat("From", nrow(.$beta), "samples after", .$adapt, "iterations burning in\n")
  base::print(base::summary(.$beta))
  
  if (show.plot) {
    par(mfrow=c(1, ncol(.$beta)))
    for (i in 1:ncol(.$beta)) {
      plot(density(.$beta[,i]), main=colnames(.$beta)[i])
    }
  }
}

Now that's probably not the greatest design in the world, but it only took me a few minutes to put it together and it's incredibly easy to modify or extend. Try it !

Thanks to Brian Neelon for making his MCMC logistic regression code available (http://www.duke.edu/~neelo003/r/).

Homogeneity analysis of hierarchical classifications

2010-08-10T00:16:00.000-07:00

I've spent more years than I care to remember analysing vegetation survey data (typically species abundances in plots) using a variety of software including my own algorithms coded in FORTRAN and C++. A recent query on the r-help list, about how to determine the number of groups to define in a hierarchical classification produced with the hclust function, prompted me to unearth one of these algorithms, homogeneity analysis¹, which can help to visualize how different levels of grouping partition the variability in a distance matrix.

This algorithm is extremely simple. The classification is progressively divided into groups, with all groups being defined at the same dendrogram height. At each level of grouping, the average of within-group pairwise distances is calculated. Homogeneity is then defined as:


      H = 1 - Dav_within-group - Dav_total

where Dav_total is the average pairwise distance in the dataset as a whole.

For data were there exist well-defined clusters of values, a plot of homogeneity against number of groups will display an 'elbow' where the initial rapid increase in homogeneity turns to a more gradual increase. The example above shows a classification of the USArrests dataset and corresponding homogeneity plot which suggests defining 7 groups. It was generated as follows:


> d <- dist(USArrests)
> hc <- hclust(d, method="average")
> h <- hmg(d, hc)
> plot(h, type="s", main="USArrests homogeneity")
> abline(v=7, lty=2)
> plot(hc, labels=FALSE, main="USArrests dendrogram")
> rect.hclust(hc, k=7)

Here is the code for the hmg function:


function(d, hc) {

# R version of homogeneity analysis as described in:
#   Bedward, Pressey and Keith. 1992. Homogeneity analysis: assessing the 
#   utility of classifications and maps of natural resources
#   Australian Journal of Ecology 17:133-139.
#
# Arguments:
# d  - either an object produced by dist() or a vector of
#      pairwise dissimilarity values ordered in the manner of
#      a dist result
#
# hc - classification produced by hclust()
#
# Value:
# A two column matrix of number of groups and corresponding homogeneity value
#
# This code by Michael Bedward, 2010

  m <- hc$merge

  # check args for consistency
  Nobj <- nrow(m) + 1
  if (length(d) != Nobj * (Nobj-1) / 2) {
    dname <- dparse(substitute(d))
    hcname <- dparse(substitute(hc))
    stop(paste(dname, "and", hcname, "refer to different numbers of objects"))
  }

  distMean <- mean(d)
  grp <- matrix(c(1:Nobj, rep(0, Nobj)), ncol=2)

  # helper function - decodes values in the merge matrix 
  # from hclust
  getGrp <- function( idx ) {
    ifelse( idx < 0, grp[-idx], getGrp( m[idx, 1] ) )
  }

  grpDistTotal <- numeric(Nobj)
  grpDistCount  <- numeric(Nobj)
  hmg <- numeric(Nobj)
  hmg[1] <- 0
  hmg[Nobj] <- 1

  distTotal <- 0
  distCount <- 0
  for (step in 1:(Nobj-2)) {
    mrec <- m[step, ]
    grp1 <- getGrp(mrec[1])
    grp2 <- getGrp(mrec[2])

    distTotal <- distTotal - grpDistTotal[grp1] - grpDistTotal[grp2]
    distCount <- distCount - grpDistCount[grp1] - grpDistCount[grp2]

    grpDistTotal[grp1] <- grpDistTotal[grp1] + grpDistTotal[grp2]
    grpDistCount[grp1] <- grpDistCount[grp1] + grpDistCount[grp2]

    grp1Obj <- which(grp == grp1)
    grp2Obj <- which(grp == grp2)

    for (i in grp1Obj) {
      for (j in grp2Obj) {
        ii <- min(i,j)
        jj <- max(i,j)
        distIdx <- Nobj*(ii-1) - ii*(ii-1)/2 + jj-ii
        grpDistTotal[grp1] <- grpDistTotal[grp1] + d[distIdx]
        grpDistCount[grp1] <- grpDistCount[grp1] + 1
      }
    }

    # update group vector
    grp[grp == grp2] <- grp1

    distTotal <- distTotal + grpDistTotal[grp1]
    distCount <- distCount + grpDistCount[grp1]

    hmg[Nobj - step] <- 1 - distTotal / ( distCount * distMean )
  }

  out <- matrix(c(1:Nobj, hmg), ncol=2)
  colnames(out) <- c("ngroups", "homogeneity")
  out
}

Reference
¹M. Bedward, D. A. Keith, R. L. Pressey. 1992. Homogeneity analysis: Assessing the utility of classifications and maps of natural resources. Australian Journal of Ecology 17(2) 133-139.

Meeting in the middle; or fudging model II regression with nls

2010-08-03T21:15:00.000-07:00

My colleague Karen needed an equation to predict trunk diameter given tree height, which she hoped to base on measurements of trees in semi-arid Australian woodlands. This is the dark art of allometry and a quick google found a large number of formulae that have been used in different studies of tree dimensions. No problem: I started to play with a few of them and eventually settled on this one:

dbh = exp( b0 + b1 / (b2 + h) )
where: dbh is trunk diameter at breast height; h is tree height.

Karen also needed to do reverse predictions, ie. predict a tree's height given its trunk diameter. Again no problem, the inverse equation is simply:

height = b1 / (log( dbh ) - b0) - b2

But then, the pièce de résistance: the forward predictions had to agree with the reverse predictions, ie. if plugging height h into the forward equation gave trunk diameter d, then plugging d into the reverse equation should get you back to h. Karen pointed out that this seemed to be a Model II regression problem because as well as needing symmetric predictions, both variables were subject to measurement error.

R has at least two packages with Model II regression functions: lmodel2 and smatr, but both seemed to be restricted to a single predictor (I'd be interested to hear if I'm wrong about that). Meanwhile, the trusty nls function seemed to do a good job of fitting the forward and reverse equations while dealing with the pattern of variance displayed by the tree data (mu^3).

A solution arrived in the form of this post in the r-help list archive explaining how the coefficients from the forward and reverse fits could be combined by taking their geometric mean to arrive at a form of Model II regression. We tried this and it appeared to work nicely.
The plot above shows the separate fits, plus the combined fit which does indeed give symmetric predictions, for one of the tree species.

Circle packing with R

2010-07-26T00:54:00.000-07:00

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. Here are two layouts for the same input set...

The algorithm is wonderfully simple. Each circle in the input list is compared to those following it. If two circles overlap by more than the allowed amount, they are moved apart until the overlap is acceptable. The distance moved by each circle is proportional to the radius of the other to approximate inertia (very loosely); thus when a small circle is overlapped by a large circle, the small circle moves furthest. This process is repeated iteratively until no more movement takes place (acceptable layout) or a maximum number of iterations is reached (layout failure). To avoid edge effects, the bounding rectangle is treated as a toroid. For my purposes, I only require a circle's centre to be within the plotted rectangle.

You can see in the plots that the algorithm tends to produce clusters of small circles around big ones. For the woodland simulation model this is a nice property (saplings clustering around parent trees) but for other applications the algorithm could be further modified to lessen or avoid this effect.

The code for the basic function is below. The set of input circles are described by the config matrix argument. Although this function produces good results, it takes a long time to run when the number of circles is large. However, a second version of the function, with most of the calculations moved into C code, runs much faster (code not posted here but available on request).


pack.circles <- function(config, size=c(100, 100), max.iter=1000, overlap=0 ) {
  #
  # Simple circle packing algorithm based on inverse size weighted repulsion
  #
  # config   - matrix with two cols: radius, N
  # size     - width and height of bounding rectangle
  # max.iter - maximum number of iterations to try
  # overlap  - allowable overlap expressed as proportion of joint radii

  # ============================================================================
  #  Global constants
  # ============================================================================
  # round-off tolerance
  TOL <- 0.0001

  # convert overlap to proportion of radius
  if (overlap < 0 | overlap >= 1) {
    stop("overlap should be in the range [0, 1)")
  }
  PRADIUS <- 1 - overlap

  NCIRCLES <- sum(config[,2])

  # ============================================================================
  #  Helper function - Draw a circle
  # ============================================================================
  draw.circle <- function(x, y, r, col) { 
    lines( cos(seq(0, 2*pi, pi/180)) * r + x, sin(seq(0, 2*pi, pi/180)) * r + y , col=col )
  }


  # ============================================================================
  #  Helper function - Move two circles apart. The proportion of the required
  #  distance moved by each circle is proportional to the size of the other 
  #  circle. For example, If a c1 with radius r1 overlaps c2 with radius r2,
  #  and the movement distance required to separate them is ds, then c1 will
  #  move ds * r2 / (r1 + r2) while c2 will move ds * r1 / (r1 + r2). Thus,
  #  when a big circle overlaps a little one, the little one moves a lot while
  #  the big one moves a little.
  # ============================================================================
  repel <- function(xyr, c0, c1) {
    dx <- xyr[c1, 1] - xyr[c0, 1]
    dy <- xyr[c1, 2] - xyr[c0, 2]
    d <- sqrt(dx*dx + dy*dy)
    r <- xyr[c1, 3] + xyr[c0, 3]
    w0 <- xyr[c1, 3] / r
    w1 <- xyr[c0, 3] / r

    if (d < r - TOL) {
      p <- (r - d) / d
      xyr[c1, 1] <<- toroid(xyr[c1, 1] + p*dx*w1, 1)
      xyr[c1, 2] <<- toroid(xyr[c1, 2] + p*dy*w1, 2)
      xyr[c0, 1] <<- toroid(xyr[c0, 1] - p*dx*w0, 1)
      xyr[c0, 2] <<- toroid(xyr[c0, 2] - p*dy*w0, 2)

      return(TRUE)
    }

    return(FALSE)
  }


  # ============================================================================
  #  Helper function - Adjust a coordinate such that if it is distance d beyond
  #  an edge (ie. outside the area) it is moved to be distance d inside the 
  #  opposite edge. This has the effect of treating the area as a toroid.
  # ============================================================================
  toroid <- function(coord, axis) {
    tcoord <- coord

    if (coord < 0) {
      tcoord <- coord + size[axis]
    } else if (coord >= size[axis]) {
      tcoord <- coord - size[axis]
    }

    tcoord
  }


  # ============================================================================
  #  Main program
  # ============================================================================

  # ------------------------------------------
  # create a random initial layout
  # ------------------------------------------
  xyr <- matrix(0, NCIRCLES, 3)

  pos0 <- 1
  for (i in 1:nrow(config)) {
    pos1 <- pos0 + config[i,2] - 1
    xyr[pos0:pos1, 1] <- runif(config[i, 2], 0, size[1])
    xyr[pos0:pos1, 2] <- runif(config[i, 2], 0, size[2])
    xyr[pos0:pos1, 3] <- config[i, 1] * PRADIUS
    pos0 <- pos1 + 1
  }

  # ------------------------------------------
  # iteratively adjust the layout
  # ------------------------------------------
  for (iter in 1:max.iter) {
    moved <- FALSE
    for (i in 1:(NCIRCLES-1)) {
      for (j in (i+1):NCIRCLES) {
        if (repel(xyr, i, j)) {
          moved <- TRUE
        }
      }
    }
    if (!moved) break
  }

  cat(paste(iter, "iterations\n"));
  
  # ------------------------------------------
  # draw the layout
  # ------------------------------------------
  plot(0, type="n", xlab="x", xlim=c(0,size[1]), ylab="y", ylim=c(0, size[2]))

  xyr[, 3] <- xyr[, 3] / PRADIUS
  for (i in 1:nrow(xyr)) {
    draw.circle(xyr[i, 1], xyr[i, 2], xyr[i, 3], "gray")
  }

  # ------------------------------------------
  # return the layout
  # ------------------------------------------
  colnames(xyr) <- c("x", "y", "radius")
  invisible(xyr) 
}

Another clip round the ear

2010-04-27T06:59:00.000-07:00

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

}

A clip round the ear

2010-04-23T02:24:00.001-07:00

Inspired by a recent discussion on the JTS user list about decomposing simple polygons into triangles, here is a program that demonstrates the brute force, ear-clipping algorithm. To run this you will need the JTS library which you can get from the project web site or as a maven artifact (version 1.11) from the Maven central repository.


import com.vividsolutions.jts.algorithm.Angle;
import com.vividsolutions.jts.geom.Coordinate;
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.util.ArrayList;
import java.util.List;

/**
 * Demonstrates brute force approach to the ear clipping algorithm
 * to triangulate a polygon. This demo doesn't yet cope with
 * polygons that have holes.
 *
 * @author Michael Bedward
 */
public class EarClipping {
    private static final double EPS = 1.0E-4;

    GeometryFactory gf = new GeometryFactory();

    public static void main(String[] args) throws Exception {
        new EarClipping().demo();
    }

    /**
     * Demonstrate the ear-clipping algorithm
     * @throws Exception
     */
    private void demo() throws Exception {
        WKTReader reader = new WKTReader(gf);
        Polygon poly = (Polygon )reader.read("POLYGON((0 0, 5 0, 2 3, 5 5, 10 0, 10 -5, 5 -2, 0 -5, 0 0))");
        System.out.println("Input polygon:");
        System.out.println(poly);

        Geometry ears = triangulate(poly);
        final int n = ears.getNumGeometries();

        System.out.println();
        System.out.println("Found " + n + " ears:");
        for (int i = 0; i < n; i++) {
            System.out.println(ears.getGeometryN(i));
        }
    }

    /**
     * Brute force approach to ear clipping
     *
     * @param inputPoly input polygon
     * @return GeometryCollection of triangular polygons
     */
    private Geometry triangulate(Polygon inputPoly) {
        if (inputPoly.getNumInteriorRing() > 0) {
            throw new IllegalArgumentException("Can't deal with polgons that have holes yet");
        }

        List<Polygon> ears = new ArrayList<Polygon>();
        Polygon workingPoly = (Polygon) inputPoly.clone();

        Coordinate[] coords = workingPoly.getBoundary().getCoordinates();
        coords = removeColinearVertices(coords);
        int N = coords.length - 1;

        boolean finished = false;
        int k0 = 0;
        do {
            int k1 = (k0 + 1) % N;
            int k2 = (k0 + 2) % N;
            LineString ls = gf.createLineString(new Coordinate[] {coords[k0], coords[k2]});
            
            if (workingPoly.covers(ls)) {
                Polygon ear = gf.createPolygon(
                        gf.createLinearRing(new Coordinate[]{coords[k0], coords[k1], coords[k2], coords[k0]}),
                        null);
                ears.add(ear);
                workingPoly = (Polygon) workingPoly.difference(ear);
                coords = workingPoly.getBoundary().getCoordinates();
                coords = removeColinearVertices(coords);
                N = coords.length - 1;
                k0 = 0;
                if (N == 3) {  // triangle
                    ears.add(gf.createPolygon(gf.createLinearRing(coords), null));
                    finished = true;
                }
            } else {
                k0++ ;
            }
        } while (!finished);

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

    /**
     * Remove co-linear vertices. TopologyPreservingSimplifier could be
     * used for this but that seems like over-kill.
     *
     * @param coords polygon vertices
     * @return coordinates with any co-linear vertices removed
     */
    private Coordinate[] removeColinearVertices(Coordinate[] coords) {
        final int N = coords.length - 1;
        List<Coordinate> coordList = new ArrayList<Coordinate>();

        for (int j = 1; j <= N; j++) {
            int i = (j - 1) % N;
            int k = (j + 1) % N;
            if (Math.abs(Math.PI - Angle.angleBetween(coords[i], coords[j], coords[k])) > EPS) {
                coordList.add(coords[j]);
            }
        }

        coordList.add(new Coordinate(coordList.get(0)));
        return coordList.toArray(new Coordinate[0]);
    }

}

Cubic spline interpolation

2010-03-20T23:15:00.000-07:00

I spent a while today casting around for a simple, fast and just-good-enough cubic interpolation algorithm to use in the GeoTools library where it is to serve as a 'filter function' that can be used when rendering maps.

The ACM digitial library has a mind-boggling collection of references on interpolation methods, most of which were fascinating but also over the top for this application.

In the end I put together the following code based on a Wikipedia article about cubic Hermite splines. Here it is...



/**
 * Cubic hermite spline interpolation. This is adapted from the description of the
 * algorithm at: http://en.wikipedia.org/wiki/Cubic_Hermite_spline.
 * Tangent caculations are done with simple finite differencing in the interests
 * of speed.
 * <p>
 * The input arrays xi and yi contain the coordinates of four interpolation
 * points defining three segments with the middle segment containing the point
 * for which we seek an interpolated value.
 *
 * @param x x ordinate of the point for which we seek an interpolated value
 *          and which lies between xi[1] and xi[2]
 * @param xi interpolation point x ordinates
 * @param yi interpolation point y ordinates
 *
 * @return interpolated value for x
 */
public double cubic(double x, double[] xi, double[] yi) {
    double span01 = xi[1] - xi[0];
    double span12 = xi[2] - xi[1];
    double span23 = xi[3] - xi[2];

    double t = (x - xi[1]) / span12;
    double t2 = t*t;
    double t3 = t2 * t;

    double m1 = 0.5 * ((yi[2] - yi[1]) / span12 + (yi[1] - yi[0]) / span01);
    double m2 = 0.5 * ((yi[3] - yi[2]) / span23 + (yi[2] - yi[1]) / span12);

    double y = (2*t3 - 3*t2 + 1) * yi[1] +
               (t3 - 2*t2 + t) * span12 * m1 +
               (-2*t3 + 3*t2) * yi[2] +
               (t3 - t2) * span12 * m2;

    return y;
}

I make no claims about the elegance or efficiency of this code other than suspecting that it rates badly in both areas, but it does a perfectly acceptable job when I throw some test data at it...

K ?

2009-11-26T23:23:00.000-08:00

Whisper the terrible words type erasure to any moderately experienced Java programmer and you're likely to see:

fear
loathing
soul-crushing sense of defeat
all of the above

Now if you are a moderately experienced Java programmer, but you're thinking to yourself "what's type erasure ?" my advice is this:

Stop reading now !

Sometimes ignorance really is bliss.

My own rule of thumb, when confronted by the need to query the identity of K, or V, or T or any of their evil friends, in the innards of a class with generic type parameters, is to run away. Sadly, there are circumstances when that's just not an option. Although I deeply empathise with those who equate Java reflection with incest and folk dancing there are times when it's black magic or nothing. But where to find the magic ?

Enter Richard Gomes and his article: Using TypeTokens to retrieve generic parameters. This demonstrates the seemingly impossible feat of pulling erased rabbits out of Java's hat.

I'm off to try it now... "Nothing up my sleeve..."

PragProWriBlockMo

2009-11-17T16:33:00.000-08:00

PragProWriMo has been even more challenging than I'd feared thanks to my sub-gnat concentration span and memory like a... you know, round thing with holes. While others have been trotting out chapter outlines, finding their voice and defining their readership, I've been shuffling along the beach of irrelevancy and gazing out at the ocean of unfinishedness.

It's amazing / pathetic / pathological (select all that apply) just how intimidating this simple daily writing exercise has become. The main problem has been the notion of a book hovering over the writing. No matter how many times I've told myself to treat the writing as a pump-priming exercise rather than an examination, I haven't been able to shake the anxiety of not "seeing the book" in my mind's eye.

But I will not surrender ! Well, that's not quite true because I did a few days ago, giving up on the whole thing as too hard. Now though, I'm putting the white flag back in the cupboard and making a new start. Rather than imagining the goal as a book I'm going to structure it as a set of short web tutorials.

I wonder if I can count this post ?

PragProWriMo begins

2009-11-01T18:00:00.000-08:00

November is here and this not-so-young-anymore geek's heart has turned to writing. I'm taking up the challenge of PragProWriMo by committing (forcing, tricking...) myself to write every day. The target is eighty pages by the end of the month.

As someone who has routinely taken years to write an eight page scientific paper (only four pages really if you subtract the tables, figures and compulsory, gratuitous citations of the publications of those who you suspect might review your paper) eighty pages seems an imposing target.

Ah well grasshopper, the longest journey etc. etc.

My chosen geek book working title is Biological Models in Java. That's likely to change as the writing progresses of course... I can already see it wandering towards Mathematical Recreations in Java, Ruby, Groovy, Clojure, R and that other language (you know, the one with the clicking sounds).

Unleash your inner literary geek

2009-10-28T20:41:00.000-07:00

Those good fellows at Pragmatic Programmers have come up with the geeky equivalent of NaNoWriMo - the annual event that encourages people to try their hand at writing - which they have christened PragProWriMo.

This is the chance for you to get started on that brilliant computing book that's been festering in the back of your mind for the last few years. All you have to do is commit to writing a few pages every day throughout November. And yes, every day means every day.

I'm planning a tome on displacement behaviour myself...

Measures of confusion with JSR-275

2009-08-10T22:35:00.000-07:00

I love the idea behind JSR-275, a units of measurement API for Java (see here for a background article by Jean-Marie Dautelle, the lead author) but I wish there were more examples out there about how to use it properly.

I'm presently writing some code for a spatial simulation. I want the user to be able to specify constraints on spatial entities such as their maximum area, width, movement over time and so forth. Bringing JSR-275 to bear on this means that I can use a consistent set of units within the code while leaving the user free to express their input data in whatever units are convenient for them.

The theory is beautiful but my, alas, my code is ugly.

Take this snippet for example...

// get a constraint from the user
Constraint constraint = ...

// get the unit used internally for this constraint (e.g. "m")
Unit<?> stdUnit = ... 

// the value for the constraint (in the user's units)
double value = constraint.getValue();

// get the label for the user's units (e.g. "ft")
String unitName = constraint.getUnitLabel(); 

// match the label to a jsr-275 class
Unit<?> unit = Unit.valueOf(unitName);

if ( !unit.isCompatible(stdUnit) ) { 
 // e.g. mixing up length and area units
 // throw an exception
}

if (!unit.equals(stdUnit)) {
 /*
    * The user is exercising their freedom to use a different unit
    * to us (curse them) so we convert their value to one
    * expressed in our units.
    *
    * (This is where it gets ugly)
    */
   UnitConverter converter = null;  

   // is it an area constraint ?
   if (unit.isCompatible(Area.UNIT)) {
       converter = unit.asType(Area.class).getConverterTo( 
               stdUnit.asType(Area.class) );

   // is it a length constraint ?
   } else if (unit.isCompatible(Length.UNIT)) {
       converter = unit.asType(Length.class).getConverterTo( 
               stdUnit.asType(Length.class) );

   // bummer - it was something else I forgot to allow for
   } else {
       throw new RuntimeException("Constraint system doesn't allow for "
               + unit.toString());
   }

   value = converter.convert(value);

Now I'd like to avoid the if-else block that says Is it an area ? Is it a length ? and just say...

   converter = unit.getConverterTo(stdUnit);

...but the compiler isn't happy with this. It wants the Unit references to be fully parameterized.

If I was only working with length, for instance, things would be simple...

   Unit<Length> stdUnit = ... // our standard unit (e.g. m)
   Unit<Length> unit = ...  // the user's unit (e.g. cubits)
   UnitConverter converter = unit.getConverterTo( stdUnit );

That's concise and elegant. But since I want to allow constraints to be applied to different quantities (length and area at the moment; later possibly others) I seem stuck with the if-else alternative parameterizations above.

I'm sure I must be missing something here. So, dear reader, if you are a JSR 275 guru please enlighten me !