bulletr: Algorithms for Matching Bullet Lands

#' Estimate the twist in a bullet land
#' 
#' Estimation of the twist in a barrel follows roughly the process described by Chu et al (2010).
#' At the moment, twist is estimated from a single land - but the twist should be the same for the whole barrel. Therefore all lands of the same barrel should
#' have the same twist.
#' A note on timing: at the moment calculating the twist rate for a bullet land takes several minutes.
#' XXX TODO XXX make the different methods a parameter. Also, accept other input than the path - if we start with the flattened bulletland we get results much faster.
#' @param path to a file in x3p format
#' @param bullet data in x3p format as returned by function read_x3p
#' @param twistlimit Constraint the possible twist value
#' @param cutoff Use this for the quantile cutoff
#' @return numeric value estimating the twist
#' @export
#' @importFrom robustbase lmrob
#' @importFrom stats lm coef median 
#' @importFrom utils head tail
#' @importFrom x3ptools read_x3p
#' @examples
#' \dontrun{
#' # execution takes several minutes
#' load("data/b1.rda")
#' twist <- getTwist(path="barrel 1 bullet 1", bullet = b1, twistlimit=c(-2,0)*1.5625)
#' }
getTwist <- function(path, bullet = NULL, twistlimit = NULL, cutoff = .75) {
    x <- NULL
    r.squared <- NULL
    r.squared.robust <- NULL
    
  if (is.null(bullet)) bullet <- read_x3p(path)
  cat(path)
  cat("\n")
  
  gg115 <- processBullets(
    bullet, 
    x=bullet$header.info$profile_inc*c(0,bullet$header.info$num_profiles), 
    name=path)
  # qplot(data=gg115, x=y, y=x, fill=resid, geom="tile")+scale_fill_gradient2()
  
  xs <- unique(gg115$x)
  twist <- NULL
  ccf <- NULL
  aligned <- list()
  
  for (i in seq_along(xs)[-1]) {
    profiles <- subset(gg115, x %in% xs[c(i,i-1)])  
    profiles$bullet <- sprintf("%s-%s", path, profiles$x)
    aligned[[i]] <- bulletAlign(profiles, value="resid")
    twist <- c(twist, aligned[[i]]$lag)
    ccf <- c(ccf, aligned[[i]]$ccf)
  }
  #qplot(xs, c(0, twist)) + ylim(c(-10,10))
  #qplot(xs, c(1, ccf)) + geom_hline(yintercept=0.95, colour="red")
  #qplot(xs, cumsum(c(0, twist))) +ylim(c(700,1100))

#  if (!is.null(twistlimit)) {
    # is twistlimit a vector with two numbers?
    twistConstraint <- pmin(twist, max(twistlimit))
    twistConstraint <- pmax(twistConstraint, min(twistlimit))
#  } else twistConstraint = twist
  dframe <- data.frame(x = xs, twist=cumsum(c(0, twist)), 
                       twistConstraint = cumsum(c(0, twistConstraint)), 
                       ccf = c(0,ccf))
  #qplot(x, twist, data=subset(dframe, between(x, 220, 600))) +geom_smooth(method="lm")
  
  Rs <- data.frame(zoo::rollapply(data=dframe$twist, width=200, FUN=function(twist) {
    x <- 1:length(twist)
    m <- lm(twist~x)
    m2 <- try(robustbase::lmrob(twist~x, na.action=na.omit), silent=TRUE)
    if (class(m2) == "try-error") {
      cat(sprintf("NAs in robust estimation of twist in land %s\n", path))
 #     browser()
      r.squared.robust=NA
      twistRobust=NA
    } else {
      r.squared.robust=summary(m2)$r.squared
      twistRobust=coef(m2)[2]
    }
    
    data.frame(r.squared=summary(m)$r.squared, twist=coef(m)[2],
               r.squared.robust=r.squared.robust, 
               twistRobust=twistRobust)
  }, by=1))
#  browser()
  RConstraint <- data.frame(zoo::rollapply(data=dframe$twistConstraint, width=200, FUN=function(twist) {
    x <- 1:length(twist)
    m <- lm(twist~x)
    m2 <- try(robustbase::lmrob(twist~x, na.action=na.omit), silent=TRUE)
    if (class(m2) == "try-error") {
      cat(sprintf("NAs in robust estimation of twist in land %s\n", path))
 #     browser()
      r.squared.robust=NA
      twistRobust=NA
    } else {
      r.squared.robust=summary(m2)$r.squared
      twistRobust=coef(m2)[2]
    }
    data.frame(r.squared=summary(m)$r.squared, twist=coef(m)[2],
               r.squared.robust=r.squared.robust, 
               twistRobust=twistRobust)
  }, by=1))
  
    
  q75 <- quantile(Rs$r.squared, probs=cutoff)
  twist <- median(subset(Rs, r.squared > q75)$twist)
  q75r <- quantile(Rs$r.squared.robust, probs=cutoff, na.rm=TRUE)
  twistRobust <- median(subset(Rs, r.squared.robust > q75r)$twistRobust, na.rm=TRUE)
  
  q75C <- quantile(RConstraint$r.squared, probs=cutoff)
  twistC <- median(subset(RConstraint, r.squared > q75)$twist)
  q75rC <- quantile(RConstraint$r.squared.robust, probs=cutoff, na.rm=TRUE)
  twistRobustC <- median(subset(RConstraint, r.squared.robust > q75r)$twistRobust, na.rm=TRUE)
  
  
    
  data.frame(twist=twist, min.r.squared=q75, twistRobust=twistRobust, min.r.squared.robust=q75r,
             twistC=twistC, min.r.squaredC=q75C, twistRobustC=twistRobustC, min.r.squared.robustC=q75rC)
}

#' Plot a bullet land using plotly
#' 
#' @param path The path to the x3p file
#' @param bullet If not null, use this pre-loaded bullet
#' @param sample integer value. take every 1 in sample values from the surface matrix
#' @param ... parameters passed on to plot_ly call
#' @importFrom plotly plot_ly
#' @importFrom x3ptools read_x3p
#' @export
#' @examples 
#' data(br411)
#' plot_3d_land(bullet=br411, sample=2)
plot_3d_land <- function(path, bullet = NULL, sample=1, ...) {
    if (is.null(bullet)) bullet <- read_x3p(path)
    surfmat <- bullet$surface.matrix
    if (sample != 1) {
      xindx <- seq.int(from=1, to=dim(surfmat)[1], by=sample)
      yindx <- seq.int(from=1, to=dim(surfmat)[2], by=sample)
      surfmat <- surfmat[xindx, yindx]
    }
    
    plot_ly(z = surfmat, type = "surface", ...)
}

#' Process x3p file 
#' 
#' x3p file of a 3d topological bullet surface is processed at surface crosscut x, 
#' the bullet grooves in the crosscuts are identified and removed, and a loess smooth 
#' is used (see \code{?loess} for details) to remove the big structure. 
#' @param bullet file as returned from read_x3p
#' @param name name of the bullet
#' @param x (vector) of surface crosscuts to process (in meters). 
#' @param grooves The grooves to use as a two element vector, if desired
#' @param span The span for the loess fit
#' @param window The mean window around the ideal crosscut
#' @param ... Additional arguments, passed to the get_grooves function
#' @return data frame
#' @import dplyr
#' @importFrom zoo na.trim
#' @importFrom x3ptools read_x3p
#' @importFrom x3ptools x3p_to_df
#' @export
#' @examples
#' data(br411)
#' br411_processed <- processBullets(br411, name = "br411")
processBullets <- function(bullet, name = "", x = as.numeric(bullet$feature.info$Axes$CX$Increment[[1]]) * 100, grooves = NULL, span = 0.75, window = 0, ...) {
    y <- value <- NULL
    
    ## Convert x to be meters
    if (x > 1) {
        x <- x * 1e-6
    }
    
    if (!is.data.frame(bullet)) {
        crosscuts <- unique(x3p_to_df(bullet)$x)
        crosscuts <- crosscuts[crosscuts >= min(x)]
        crosscuts <- crosscuts[crosscuts <= max(x)]
    } else {
        crosscuts <- x
    }
    
    if (length(x) > 2) crosscuts <- crosscuts[crosscuts %in% x]
    
    list_of_fits <- lapply(crosscuts, function(myx) {
        if (!is.data.frame(bullet)) bullet <- x3p_to_df(bullet)
        
        br111 <- bullet %>%
            na.trim %>%
            filter(x >= myx - window, x <= myx + window) %>%
            group_by(y) %>%
            summarise(x = myx, value = mean(value, na.rm = TRUE)) %>%
            dplyr::select(x, y, value) %>%
            as.data.frame
        if (is.null(grooves)) {
            br111.groove <- get_grooves(br111, ...)
        } else {
            br111.groove <- list(groove = grooves)
        }
        
        myspan <- span
        if (myspan > 1) {
            ## Use the nist method
            myspan <- myspan / diff(br111.groove$groove)
        }
        fit_loess(br111, br111.groove, span = myspan)$resid$data
    })
    lof <- list_of_fits %>% bind_rows
    
    data.frame(lof, bullet = name, stringsAsFactors = FALSE)
}