R/egg.polygon.area.R

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

#' Polygon area
#'
#' Iteratively find the area of a polygon on the surface of the egg
#'
#' @param egg an object of class \code{egg.fit}
#' @param vertices points describing the vertices of the polygon in spherical coordinates, as obtained from \code{egg.coords}
#' @param N number of triangle subdivisions
#'
#' @return The area of the polygon in pixels^2
#'
#' @examples
#' stop()
#'
#' @export
egg.polygon.area = function(egg, vertices, N=6) {

  # Step 1: triangulate the input polygon
  p = sph2image(egg, vertices)
  if(all(p[nrow(p),] == p[1,])) {
    p = p[-nrow(p),]
  }
  t = triangulate.polygon(p)

  # Step 2: subdivide each triangle N times
  p = cbind(p, 0)
  for(i in 1:nrow(p)) {
    p[i,] = sph2cart(get.sph(p[i,1], p[i,2], egg))
  }
  if(N > 0) {
    for(i in 1:N) {
      new.t = matrix(0, nrow=nrow(t)*4, ncol=3)
      n = 1
      for(j in 1:nrow(t)) {
        v1 = p[t[j,1],]
        v2 = p[t[j,2],]
        v3 = p[t[j,3],]
        v12 = sph2cart(egg.geodesic(egg, rbind(cart2sph(v1),cart2sph(v2)), n=2, return.points=TRUE)$points[2,])
        v13 = sph2cart(egg.geodesic(egg, rbind(cart2sph(v1),cart2sph(v3)), n=2, return.points=TRUE)$points[2,])
        v23 = sph2cart(egg.geodesic(egg, rbind(cart2sph(v2),cart2sph(v3)), n=2, return.points=TRUE)$points[2,])
        p = rbind(p, v12, v13, v23)
        new.t[n,] = c(t[j,1], nrow(p)-2, nrow(p)-1)
        new.t[n+1,] = c(t[j,2], nrow(p)-2, nrow(p))
        new.t[n+2,] = c(t[j,3], nrow(p)-1, nrow(p))
        new.t[n+3,] = c(nrow(p), nrow(p)-1, nrow(p)-2)
        n = n + 4
      }
      t = new.t
    }
  }

  # Step 3: find sum area of triangles
  A = 0
  for(i in 1:nrow(t)) {
    #egg.geodesic(egg, rbind(cart2sph(p[t[i,1],]),cart2sph(p[t[i,2],])), n=10, plot=TRUE, return.points=FALSE, col="white", lwd=1)
    #egg.geodesic(egg, rbind(cart2sph(p[t[i,1],]),cart2sph(p[t[i,3],])), n=10, plot=TRUE, return.points=FALSE, col="white", lwd=1)
    #egg.geodesic(egg, rbind(cart2sph(p[t[i,2],]),cart2sph(p[t[i,3],])), n=10, plot=TRUE, return.points=FALSE, col="white", lwd=1)
    A = A + triangle.area(p[t[i,1],], p[t[i,2],], p[t[i,3],])
  }
  return (A)
}