R/egg.fit.R

In eggs: Distances and areas on the surface of an egg

#' Estimate egg shape
#'
#' Estimates egg shape using the method of Troscianko, J. (2014) A simple tool for calculating egg shape, volume and surface area from digital images. Ibis 156, 874–878.
#'
#' @return an object of class \code{egg.fit} containing the egg shape parameters
#'
#' @examples
#' stop()
#'
#' @export
egg.fit = function() {
  options(locatorBell=FALSE)

  # Get pointed and blunt ends
  pointed.end = locator(n=1, type="p", pch=21, fg="red", bg="red", cex=0.8)
  blunt.end = locator(n=1, type="p", pch=21, fg="red", bg="red", cex=0.8)

  # Get edge points
  coords = locator(n=256, type="p", pch=21, fg="green", bg="green", cex=0.8)

  # Get rotation and translation
  v = c(pointed.end$x-blunt.end$x, pointed.end$y-blunt.end$y)
  angle = atan2(v[1], v[2]) - atan2(0, 1)
  R = rbind(c(cos(angle),-sin(angle)), c(sin(angle),cos(angle)))
  Rinv = solve(R)
  h = magnitude(c(pointed.end$x-blunt.end$x, pointed.end$y-blunt.end$y))
  t = cbind(-c(0,h/2))
  mid.point = c((pointed.end$x+blunt.end$x)/2, (pointed.end$y+blunt.end$y)/2)

  # Rotate and translate points so that the pointed end is uppermost
  pp = matrix(0, nrow=length(coords$x)+2, ncol=2)
  pp[1,] = t(R%*%cbind(c(pointed.end$x-blunt.end$x, pointed.end$y-blunt.end$y)) + t)
  pp[2,] = c(0,0) + t
  for(i in 1:length(coords$x)) {
    pp[i+2,] = t(R%*%cbind(c(coords$x[i]-blunt.end$x, coords$y[i]-blunt.end$y)) + t)
  }

  # Fit curve
  optim.fit = function(x, z.hat, d0) {
    a = x[1]
    b = x[2]
    c = x[3]
    d = (a*exp(((-z.hat^2)/(2*b^2)) + ((c*z.hat)/(b^2)) - ((c^2)/(2*b^2)))*sqrt(1 - z.hat)*sqrt(z.hat))/(pi*b)
    return (sum((d - abs(d0))^2))
  }
  z.hat = normalise(pp[,2], min(pp[,2]), max(pp[,2]), 0, 1)
  x0 = c(h*5, 2, 0)
  opt = optim(par=x0, fn=optim.fit, z.hat=z.hat, d0=pp[,1])

  # Return fitted parameters
  out = list(pointed.end=c(pointed.end$x,pointed.end$y), blunt.end=c(blunt.end$x,blunt.end$y), mid.point=mid.point, length=h, R=R, t=t, a=opt$par[1], b=opt$par[2], c=opt$par[3])
  class(out) = "egg.fit"
  return (out)
}