#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.