incircle of a polygon

Etienne Bellemare Racine <etiennebr at gmail.com> writes:

I'm looking for a way to find the largest circle fitting in an irregular
polygon. This is not the real definition of an incircle (as it has to be
tangent to all edges, but it's a more synthetic title).

The mathematical definition of the incircle of a polygon is: the largest circle fitting inside it.
See the Wikipedia entry for 'incircle'. 
If the polygon is convex, then the incircle will be tangent to at least three of the edges of the
polygon. In general it cannot be tangent to all edges of the polygon (that can only happen for regular polygons).
If the polygon is not convex, then the incircle does not need to be tangent to any edges: it might just touch the 
polygon at three concave vertices.

Now I would like to iteratively find the incircle in a polygon until a
certain radii (i.e. as if I was stuffing balls in my polygon, from the
biggest to the smallest). So I have to subtract that circle to my polygon.
Should I subtract it from my polygon by adding an inside loop (setting
hole to TRUE, is it that simple) or can I substract it directly from the
owin object ?

You can use disc() to create an owin object that represents a circle,
and intersect.owin and complement.owin to remove the disc from your polygon.

Here is a piece of code that will do what you wanted.

It's rather fun to try e.g.
          data(demopat)
          stuffballs(demopat$window, 100, TRUE)

Please note that the current version of spatstat (1.15-2) converts the polygon to a pixel mask
and uses the 8-connected distance transform to compute the incircle. Effectively this treats the 
circles as octagons. So, don't expect the circles to be perfectly tangent to the edges. 
You can see a few anomalies in the example above.

In the next version of spatstat, the distmap algorithm for polygons will be changed so that this
problem is vastly ameliorated.

Adrian Baddeley

----------------------
require(spatstat)

stuffballs <- function(W, rmin, plotit=TRUE) {
  stopifnot(is.owin(W))
  stopifnot(rmin > 0)
  if(plotit)
    plot(W)
  W <- as.mask(W)
  B <- as.rectangle(W)
  circles <- data.frame(x=numeric(0), y=numeric(0), r=numeric(0))
  ncircles <- 0
  while(TRUE) {
    circ <- incircle(W)
    if(circ$r <= rmin)
      break
    circles <- rbind(circles, as.data.frame(circ))
    D <- with(circ, disc(r, c(x,y)))
    if(plotit)
      plot(D, add=TRUE)
    W <- intersect.owin(W, complement.owin(D, B))
    if(area.owin(W) == 0)
      break
  }
  return(circles)
}