MORE COMPLETE CODE WITH DATA AT https://github.com/schwilklab/bark-contour. R script to read digitized bark contours and calculate bark thickness distributions. See Schwilk, D.W., M. Gaetani and H.M. Poulos. 2013. Oak bark allometry and fire survival strategies in the Chihuahuan Desert Sky Islands, Texas, USA. PLoS One: 8(11) e79285. DOI: 10.137…

BarkContour.R

## R code to calculate bark thickness distributions from a digitized contour
## D.W. Schwilk 2013
## (See Adams and Jackson 1995 for description of the method)

SSPAR = 0.7

library(ggplot2)
library(plyr)
library(grid)

themeopts <-  theme(axis.title.y = element_text(size = textsize, angle = 90, vjust=0.3),
                   axis.title.x = element_text(size = textsize, vjust=-0.3),
                   panel.background = element_rect(size = 1.6, fill = NA),
                   panel.border = element_rect(size = 1.6, fill=NA),
                   axis.text.x  = element_text(size=14),
                   axis.text.y  = element_text(size=14),
                   panel.grid.minor = element_line(colour = NA),
                   panel.grid.major = element_line(colour = NA))

############################################################
## Helper functions
###########################################################

## Euclidian distance
dist <- function(x1,y1,x2,y2){
  return(sqrt( (x1-x2)**2 + (y1-y2)**2))
}

## Intersection of 2 circles
cintersect <- function(x1,y1,x2,y2,r){
 e <- x2-x1
 f <- y2-y1
 p <- sqrt(e^2 + f^2)
 k <- (p^2) / (2*p)
 x3 <- x1 + ((e*k) /p) + ((f/p) * sqrt(r^2 - k^2))
 y3 <- y1 + ((f*k) /p) - ((e/p) * sqrt(r^2 - k^2))
 x4 <- x1 + ((e*k) /p) - ((f/p) * sqrt(r^2 - k^2))
 y4 <- y1 + ((f*k) /p) + ((e/p) * sqrt(r^2 - k^2))
 ## only return the upper (higher y) intersection
 ##print(c(x1,y1,x2,y2,x3,y3,x4,y4))
 if (is.nan(y3)) {
   x <- NA
   y <- NA
 } else if (y3 > y1) {
   x<-x3
   y<-y3
 } else {
    x<-x4
    y<-y4
  }
   return(c(x,y))
}

# x, y: the x and y coordinates of the hull points
# n: the number of points in the curve.
bez <- function(x, y, t)
	{
	outx <- 0
	outy <- 0
	n <- length(x)-1
	for (i in 0:n)
		{
		outx <- outx + choose(n, i)*((1-t)^(n-i))*t^i*x[i+1]
		outy <- outy + choose(n, i)*((1-t)^(n-i))*t^i*y[i+1]
		}
 
	return (list(x=outx, y=outy))
	}
bezierCurve <- function(x, y, n=10)
	{
	outx <- NULL
	outy <- NULL
 
	i <- 1
	for (t in seq(0, 1, length.out=n))
		{
		b <- bez(x, y, t)
		outx[i] <- b$x
		outy[i] <- b$y
 
		i <- i+1
		}
 
	return (list(x=outx, y=outy))
	}
  
# Example usage
##  <- c(4,6,4,5,6,7)
## y <- 1:6
## plot(x, y, "o", pch=20)
## points(bezierCurve(x,y,20), type="l", col="red")


### Center of tree bole find circle intersections fo each point on tree bole
### and average as best estimate of center location. x, y and r are vectors.
### One value for each bore depth
bole.center <- function(x,y,r){
  cx <- c()
  cy <- c()
  for(p1 in 1:(length(x)-1)) {
    for(p2 in (p1+1):length(x)) {
      origin <- cintersect(x[p1],y[p1],x[p2],y[p2],r)
      cx <- append(cx,origin[1])
      cy <-append(cy,origin[2])
    }
  }
  return(c(mean(cx),mean(cy)))
}


cart2polar <- function(df) {
   return(data.frame( r = sqrt(df$x^2 + df$y^2), theta = atan2(df$y,df$x))  )
 }

polar2cart <- function(df) {
  return(data.frame(x = df$r * cos(df$theta), y = df$r * sin(df$theta)))
}

local.maxima <- function(x,y,n=5) {
  out.x <- c()
  out.y <- c()
  s <- seq(min(x),max(x),length.out=n+1)
  for(i in seq(1,n) ) {
    j <- which.max(y[x>s[i] & x<s[i+1]])
    out.x <- append(out.x,  x[x>s[i] & x<s[i+1]][j]  )
    out.y <- append(out.y,  y[x>s[i] & x<s[i+1]][j])
  }
  return(list(x = out.x, y = out.y))
}


#################################################################################
###  FUNCTION: get.curves Returns data frame with inner (cambial curve
###             approximation, and outer contour), both truncated to estimable
###             arc
 # Note: we do any unit conversions prior. For now we are assuming all coords and
  # lengths are in same units.  
  # Arguments:
  # - bore.x, bore.y and blengths are three vectors that contain the x,y and
  #   depth information. These obviously must all be in the correct order as the
  #   values really should be lines in a data frame
  # - diam: the diameter of the tree at contour height
  # - contour.x/y: long vector describing the actual outer contour from the
  # - carpenters tool

get.curves <- function(bore.x, bore.y, blengths, diam, contour.x, contour.y, method="bezierPolar") {
  # get tree radius (including bark)
  r <- diam/2.0
  #print(r)
  # trim contour to arc needed based on bores
  rowsneeded = contour.x < max(bore.x) & contour.x > min(bore.x)
  contour.x <- contour.x[rowsneeded]
  contour.y <- contour.y[rowsneeded]
    
  # solve for center of bole -- first rough pass
  #print(bore.x)
  #print(bore.y)
  #plot(bore.x,bore.y)
  origin <- bole.center(bore.x,bore.y,r)
  center.x <- origin[1]
  center.y <- origin[2]
  ## print( center.x)
  ## print( center.y)
  
  ## convert all coordinates to new origin (center of bole)
  bore.x <- bore.x-center.x
  bore.y <- bore.y- center.y
  contour.x <- contour.x - center.x
  contour.y <- contour.y - center.y

     # convert all pts, contour and bores pts to polar coords with origin center
  # of bole
  bore.R <- sqrt(bore.x^2 + bore.y^2)
  bore.theta <- atan2(bore.y,bore.x)
  bore.inner.R <- bore.R - blengths
  contour.R <- sqrt(contour.x^2 + contour.y^2)
  contour.theta <-atan2(contour.y,contour.x)

  #print(bore.R)
  #print(bore.theta)

  ### now get more accurate center:
  cmax <- local.maxima(contour.theta, contour.R , 6)
   # new origin (in cart coords)

  tx <- cmax$y * cos(cmax$x)
  ty <- cmax$y * sin(cmax$x)
  ## print(tx)
  ## print(ty)
  ## print(bore.x)
  ## print(bore.y)
  new.origin <- bole.center(cmax$y * cos(cmax$x), cmax$y * sin(cmax$x), r )
  
  center.x <- new.origin[1]
  center.y <- new.origin[2]
  #print(center.x)
  #print(center.y)

  # transform to new origin:
  bore.x <- bore.x-center.x
  bore.y <- bore.y- center.y
  contour.x <- contour.x - center.x
  contour.y <- contour.y - center.y

    # convert all pts, contour and bores pts to polar coords with origin center
  # of bole
  bore.R <- sqrt(bore.x^2 + bore.y^2)
  bore.theta <- atan2(bore.y,bore.x)
  bore.inner.R <- bore.R - blengths
  contour.R <- sqrt(contour.x^2 + contour.y^2)
  contour.theta <-atan2(contour.y,contour.x)

  ## for graphing:
  bore.inner.x <- bore.inner.R *cos(bore.theta)
  bore.inner.y <- bore.inner.R * sin(bore.theta)


  ## smoothing methods each produces a smooth innner bole based on bore pts and
  ## outputs result in cartesian and polar coords methods smooth on polar
  ## coordinates ath x and y -- seems to work best to enfroce overal
  ## circularity
  if(method == "bezier") {
  t <- bezierCurve(bore.theta,bore.inner.R, length(contour.theta))
  inner.theta <- t$x
  inner.r <- t$y
   # and cartesian version
  inner.x <- inner.r * cos(inner.theta)
  inner.y <- inner.r * sin(inner.theta)  
  } else if (method=="circle") {
  ### method 1: circular bole
  ## ## average inner radius (for  inner as circle rather than loess)
   ## then get inner contour of bores as cartesian coordinates (lx,ly)
  radius.ave <- mean(bore.R-blengths)
  inner.theta = contour.theta
  inner.r = radius.ave
   # and cartesian version
  inner.x <- radius.ave * cos(contour.theta)
  inner.y <- radius.ave * sin(contour.theta)
} else if (method=="spline") {
    fspline <- splinefun(x=bore.theta,y=bore.inner.R)
    inner.theta <- contour.theta
    inner.r <- fspline(inner.theta)
    inner.x <- inner.r * cos(inner.theta)
    inner.y <- inner.r * sin(inner.theta)      
} else if (method=="smoothSpline") {
    fspline <- smooth.spline(x=bore.theta,y=bore.inner.R,spar = SSPAR, all.knots=TRUE)
    inner.theta <- contour.theta
    inner.r <- predict(fspline, inner.theta)$y
    inner.x <- inner.r * cos(inner.theta)
    inner.y <- inner.r * sin(inner.theta)    
  } else {
    print("Unknown method")

  }
  
  
  ndata <- data.frame(outer.theta = contour.theta, outer.r = contour.R,  inner.r = inner.r, outer.x = contour.x, outer.y=  contour.y, inner.x=inner.x, inner.y=inner.y, bt = contour.R-inner.r)

  bdata <- data.frame(bore.inner.x=bore.inner.x,bore.inner.y=bore.inner.y)
  
  return( list(ndata, bdata, data.frame(x=bore.x, y = bore.y)) )
}

  
read.contour <- function(fname) {
    t <- read.table(fname, header=FALSE, col.names = c("X","Y"))
    t$Y <- 10000 - t$Y
    t <- t[order(t$X),]
    t <- ddply(t, .(X), summarize, Y = mean(Y))
    return(t)
}

read.pts <- function(fname) {
    t <- read.table(fname, header=FALSE, col.names = c("X","Y","depth"))
    t$Y <- 10000 - t$Y
    t <- t[order(t$X),]
    return(t)
}


######## trial script for actual directory layout

corrected.thickness <- function(df){
  theta.out <- seq(min(df$outer.theta), max(df$outer.theta), length.out=500)
  thick.fun <- approxfun(df$outer.theta, df$bt, rule=2)
  return(thick.fun(theta.out))
}




#### main script
DATA_FOLDER <- "."
##PX_PER_MM <- 11.7

### read master summary file
bark.sum <- read.csv(paste(DATA_FOLDER, "BarkData.csv", sep="/"))

DATA_FOLDER <- "./lumped"
#results.df <- bark.sum  # for output
results.df <- data.frame(TreeID = NULL, bark.thicknesses = NULL)



for( id in bark.sum$TreeID) {

  diam <- 10 * bark.sum$D60[bark.sum$TreeID==id]  # convert to mm from cm
  res <- bark.sum$Pix.mm[bark.sum$TreeID==id]

  pfile = paste(DATA_FOLDER, "/", id, "pts.txt",sep="")
  cfile = paste(DATA_FOLDER, "/", id, "xy.txt",sep="")

  if (file.exists(pfile) & file.exists(cfile)) {
    print(id)
    points <- read.pts(pfile)
    contour <- read.contour(cfile)
    ## correct units
    cx <- contour$X / res
    cy <- contour$Y / res
    px <- points$X / res
    py <- points$Y / res

      
    t <- get.curves(px,py, points$depth, diam,cx,cy, method="smoothSpline")
    thick.df <- t[[1]]
    bore.df <- t[[2]]
 
## make graph
    #pdf(paste(id,"bark-plot.pdf",sep="_"), width=10,height=5)
    #print( qplot(outer.x, outer.y, data=thick.df) + geom_line(aes(inner.x, inner.y, data=thick.df),  size=3, color= "blue") + geom_line(aes(bore.inner.x, bore.inner.y),data=bore.df,  size=4, color= "green") ) 
    #dev.off()


    # get thicknesses
       
    thicknesses <- corrected.thickness(thick.df)
    t.df<-data.frame(TreeID = id, bark.thicknesses = thicknesses)
    results.df <- rbind(results.df,t.df)

    
  } else {  # end if file exists, otehrwise skip to next file name

    print(paste("No file found for tree ID:", id))
  }
}



### now write results.df to a csv file if you want to save the info
write.csv(results.df, file = "barkthick.csv", row.names = FALSE)



### code to make example figure of a contour
ma <- function(x,n=3){filter(x,rep(1/n,n), sides=2)}

makeContour <- function(id) {
  pfile = paste(DATA_FOLDER, "/", id, "pts.txt",sep="")
  cfile = paste(DATA_FOLDER, "/", id, "xy.txt",sep="")
  diam <- 10 * bark.sum$D60[bark.sum$TreeID==id]  # convert to mm from cm
  res <- bark.sum$Pix.mm[bark.sum$TreeID==id]

  pfile = paste(DATA_FOLDER, "/", id, "pts.txt",sep="")
  cfile = paste(DATA_FOLDER, "/", id, "xy.txt",sep="")

    points <- read.pts(pfile)
    contour <- read.contour(cfile)
    ## correct units
    cx <- contour$X / res
    cy <- contour$Y / res
    px <- points$X / res
    py <- points$Y / res

      
    t <- get.curves(px,py, points$depth, diam,cx,cy, method="smoothSpline")
    thick.df <- t[[1]]
    bore.df <- t[[2]]
    outpoints <- t[[3]]
    thick.df$outer.y2 <- ma(thick.df$outer.y)
## make graph
      boresegs <- data.frame(xo = outpoints$x,
                             yo = outpoints$y,
                             xi = bore.df$bore.inner.x,
                             yi = bore.df$bore.inner.y)
      bark.polygon <- data.frame(x = c(thick.df$outer.x, thick.df$inner.x), y = c(thick.df$outer.y, thick.df$inner.y))
  
      ggplot(thick.df, aes(x=outer.x, y=outer.y2)) +
 #     geom_ribbon(aes(ymin=outer.y2,ymax = inner.y), color = "gray", alpha = 0.5) + # filled area of bark
 #     geom_polygon(aes(x=x,y=y), data = bark.polygon, color = gray, alpha = 0.5) +   
      geom_segment(aes(x=xo,y=yo,xend=xi,yend=yi), data = boresegs,
                   size = 1.25,
                   alpha = 0.9,
                   arrow = arrow(length = unit(0.4,"cm"))) +
      geom_line(aes(inner.x, inner.y), size=1.5, linetype="dashed") + # smoothed inner bole
      geom_line(size=1.5)  +  # outer bark contour line
#      geom_point(aes(bore.inner.x, bore.inner.y),data=bore.df,  size=3) +
      scale_x_continuous("X dimension (cm)") +
      scale_y_continuous("Y dimension (cm)") +
      coord_fixed(ratio=1) + themeopts
#      geom_point(aes(x,y), data = outpoints, color = "blue", size = 3)
  
  }
  
  
#
DATA_FOLDER <- "./lumped"
makeContour("QEMO515")
ggsave("example-contour-QEMO515.pdf")


makeContour("QEMO510")