#ver. 2017.2.16 TrapsKin <- function(filename, fps, dpi, lmlength, rmlength, lmass, rmass, muscle.mass, mag = 10, spar= 5e-24, mandible = "both", fileout = filename) { # Calculates displacement, velocity, and acceleration from x-y coordinates # of mandible tip tracking data # ### Arguments: # filename: character string of filename of .csv file with headers in the # format left.x, left.y, right.x, right.y # fps: framerate that the video was filmed at # dpi: image resolution. Can find in Photoshop # mag: video magnification when filmed with microscope. 10x is default # lmlength: left mandible length in meters # rmlength: right mandible length in meters # lmass: the mass of the left mandible in kg # mandible: Which mandible should be calculated. Right, Left or Both # fileout: The name of the output file. Default is input filename # muscle.mass: The estimated mass of the adductor muscle in kg. # spar: How much smoothing. The default 5e-24 fits the data less well, and is a # conservative underestimate of max velocity and acceleration # ### Returns: # Plots of displacement, velocity, and acceleration, and a csv file with # various kinematic parameter estimates require("pspline") # Read in data. Needs to be in a csv file with headers labeled "right.x", # "right.y", "left.x", "left.y" coords <- read.table(file = paste(filename, ".csv", sep = ""), header = T, sep = ",") ### ERROR HANDLING if (mandible == "both") { if (length(coords) < 4) stop("Tracking both L and R mandible requires input to have headers 'left.x' 'left.y' 'right.x' 'right.y'") } if (mandible == "left") { if ((any(colnames(coords) == "left.x") & any(colnames(coords) == "left.y")) != T) stop("To calculate for Left Mandible, headers must be labeled 'left.x' and 'left.y'") } if (mandible == "right") { if ((any(colnames(coords) == "right.x") & any(colnames(coords) == "right.y")) != T) stop("To calculate for Right Mandible, headers must be labeled 'right.x' and 'right.y'") } ### LEFT MANDIBLE FUNCTION if (mandible == "left") { ###calculate distance between xy coords in pixels (chords of mandible movement circle) chords <- sqrt(diff(coords$left.x, differences = 1)^2+ diff(coords$left.y, differences = 1)^2) chords[which(chords == "NaN")] <- 0 chords.m <- (chords/dpi/mag)*0.0254 theta.1 <- sapply(chords.m, function(xx) 2*asin(xx/(2*lmlength))) theta <- cumsum(c(rep(0,5), theta.1, rep(0,5))) time <- seq(from=0, length.out=length(theta), by=1/fps) tt <- seq(from=0, length(theta)/fps-(1/fps), length=150) left.Psp <- smooth.Pspline(time, theta, norder=3, method=1, spar=spar) left.velocity <- predict(left.Psp, tt, 1) left.acceleration <- predict(left.Psp, tt, 2) tt2 <- tt[1:(length(tt)-50)] left.velocity2 <- left.velocity[26:(length(left.velocity)-25)] left.acceleration2 <- left.acceleration[26:(length(left.acceleration)-25)] left.inertia <- (1/3)*lmass*lmlength^2 left.kinetic <- 0.5*left.inertia*(predict(left.Psp, tt, 1)^2) left.kinetic.psp <-smooth.Pspline(tt, left.kinetic, norder=3, method=3) left.velocity.max <- max(left.velocity) left.vmax.time <- tt[which.max(left.velocity)]-tt[min(which(left.velocity > 250))] left.kinetic.vel.max <- 0.5*left.inertia*left.velocity.max^2 left.ave.power <- left.kinetic.vel.max/left.vmax.time left.mass.specific.power <- left.ave.power/(muscle.mass/2) #max stats left.duration <- time[which.max(theta)]-time[min(which(theta != 0))] left.total.angle <- theta[length(theta)] left.lin.velocity.max <- max(left.velocity*lmlength) left.acceleration.max <- max(left.acceleration) left.lin.acceleration.max <- max(left.acceleration*lmlength) stats.table <- data.frame(left.duration, left.total.angle, left.velocity.max, left.lin.velocity.max, left.acceleration.max, left.lin.acceleration.max, left.mass.specific.power) row.names(stats.table) <- filename original.dir <- getwd() dir.create(paste(getwd(),"/", fileout, sep="")) setwd(paste(getwd(), "/", fileout, sep="")) postscript(file=paste(fileout,".graph.eps", sep=""), horizontal = T, paper = "letter", onefile = T) plot(time, theta, xlab = "Time (sec)", ylab = "Angular Displacement (radians)", pch = 16, col = "blue") lines(tt, predict(left.Psp, tt), lwd=2, col="blue") plot(tt, left.velocity, xlab="Time (sec)", ylab=expression("Angular Velocity" ~ (rad ~sec^{-1})), pch=NA) lines(tt, left.velocity, lwd=2, col="blue") plot(tt, left.acceleration, xlab = "Time (sec)", ylab = expression("Angular acceleration" ~ (rad ~ sec^{-2})), pch = NA) lines(tt, left.acceleration, lwd=2, col="blue") dev.off() write.csv(stats.table, file=paste(filename,"stats.csv")) #go back to original directory setwd(original.dir) } ### RIGHT MANDIBLE if (mandible == "right") { ###calculate distance between xy coords in pixels (chords of mandible movement circle) chords <- sqrt(diff(coords$right.x, differences = 1)^2+ diff(coords$right.y, differences = 1)^2) chords[which(chords == "NaN")] <- 0 chords.m <- (chords/dpi/mag)*0.0254 theta.1 <- sapply(chords.m, function(xx) 2*asin(xx/(2*lmlength))) theta <- cumsum(c(rep(0,5), theta.1, rep(0,5))) time <- seq(from=0, length.out=length(theta), by=1/fps) tt <- seq(0, length(theta)/fps-(1/fps), length=150) right.Psp <- smooth.Pspline(time, theta, norder=3, method=1, spar=spar) right.velocity <- predict(right.Psp, tt, 1) right.acceleration <- predict(right.Psp, tt, 2) tt2 <- tt[1:(length(tt)-50)] right.velocity2 <- right.velocity[25:(length(right.velocity)-25)] right.acceleration2 <- right.acceleration[25:(length(right.acceleration)-25)] right.inertia <- (1/3)*rmass*rmlength^2 right.velocity.max <- max(right.velocity) right.vmax.time <- tt[which.max(right.velocity)]-tt[min(which(right.velocity > 250))] right.kinetic.vel.max <- 0.5*right.inertia*right.velocity.max^2 right.ave.power <- right.kinetic.vel.max/right.vmax.time right.mass.specific.power <- right.ave.power/(muscle.mass/2) #max stats right.duration <- time[which.max(theta)]-time[min(which(theta != 0))] right.total.angle <- theta[length(theta)] right.lin.velocity.max <- max(right.velocity*lmlength) right.acceleration.max <- max(right.acceleration) right.lin.acceleration.max <- max(right.acceleration*lmlength) stats.table <- data.frame(right.duration, right.total.angle, right.velocity.max, right.lin.velocity.max, right.acceleration.max, right.lin.acceleration.max, right.mass.specific.power) row.names(stats.table) <- filename write.csv(stats.table, file=paste(filename,"stats.csv")) original.dir <- getwd() dir.create(paste(getwd(),"/", fileout, sep="")) setwd(paste(getwd(), "/", fileout, sep="")) postscript(file=paste(fileout,".graph.eps", sep=""), horizontal = T, paper = "letter", onefile = T) plot(time, theta, xlab = "Time (sec)", ylab = "Angular Displacement (radians)", pch = 16, col = "red") lines(tt, predict(right.Psp, tt), lwd=2, col="blue") plot(tt, right.velocity, xlab="Time (sec)", ylab=expression("Angular Velocity" ~ (rad ~sec^{-1})), pch=NA) lines(tt, right.velocity, lwd=2, col="red") plot(tt, right.acceleration, xlab="Time (sec)", ylab=expression("Angular acceleration" ~ (rad ~ sec^{-2})), pch=NA) lines(tt, right.acceleration, lwd=2, col="red") dev.off() #go back to original directory setwd(original.dir) } ### BOTH MANDIBLE CALCULATIONS if (mandible == "both") { ###calculate distance between xy coords in pixels (chords of mandible movement circle) #Left Mandible left.chords <- sqrt(diff(coords$left.x, differences = 1)^2+ diff(coords$left.y, differences = 1)^2) right.chords <- sqrt(diff(coords$right.x, differences = 1)^2+ diff(coords$right.y, differences = 1)^2) #replace all "NaN" with 0 left.chords[which(left.chords=="NaN")] <- 0 right.chords[which(right.chords=="NaN")] <- 0 ###combine data from left and right mandibles chords <- data.frame(left.chords, right.chords) #convert chord length from pixels to m chords.m <- (chords/dpi/mag)*0.0254 #Calculate central angle of differents between each xy coord theta.list.1 <- list() theta.list.1$left.chords <- lapply(chords.m$left.chords, function(xx) 2*asin(xx/(2*lmlength))) theta.list.1$right.chords <- lapply(chords.m$right.chords, function(xx) 2*asin(xx/(2*rmlength))) #add points to beginning and end of dataframe to improve derivation theta.list.1$left.chords <- c(rep(0,5), theta.list.1$left.chords, rep(0,5)) theta.list.1$right.chords <- c(rep(0,5), theta.list.1$right.chords, rep(0,5)) theta <- as.data.frame(lapply(theta.list.1, cumsum)) colnames(theta) <- c("left", "right") time <- seq(from=0, length.out=length(theta$left), by=1/fps) tt <- seq(0, length(theta$left)/fps-(1/fps), length=150) #P-spline of theta left.Psp <- smooth.Pspline(time, theta$left, norder=3, method=1, spar=spar) right.Psp <- smooth.Pspline(time, theta$right, norder=3, method=1, spar=spar) #numerical derivitave of displacement data to estimate velocity and acceleration left.velocity <- predict(left.Psp, tt, 1) left.acceleration <- predict(left.Psp, tt, 2) right.velocity <- predict(right.Psp, tt, 1) right.acceleration <- predict(right.Psp, tt, 2) #filter out the first and last 50 elements tt2 <- tt[1:(length(tt)-50)] left.velocity2 <- left.velocity[25:(length(left.velocity)-25)] left.acceleration2 <- left.acceleration[25:(length(left.acceleration)-25)] right.velocity2 <- right.velocity[25:(length(right.velocity)-25)] right.acceleration2 <- right.acceleration[25:(length(right.acceleration)-25)] left.inertia <- (1/3)*lmass*lmlength^2 right.inertia <- (1/3)*rmass*rmlength^2 left.velocity.max <- max(left.velocity) right.velocity.max <- max(right.velocity) left.vmax.time <- tt[which.max(left.velocity)]-tt[min(which(left.velocity > 250))] right.vmax.time <- tt[which.max(right.velocity)]-tt[min(which(right.velocity > 250))] right.kinetic.vel.max <- 0.5*left.inertia*right.velocity.max^2 right.ave.power <- right.kinetic.vel.max/right.vmax.time right.mass.specific.power <- right.ave.power/(muscle.mass/2) left.kinetic.vel.max <- 0.5*left.inertia*left.velocity.max^2 left.ave.power <- left.kinetic.vel.max/left.vmax.time left.mass.specific.power <- left.ave.power/(muscle.mass/2) #max stats #leftmandible left.duration <- time[which.max(theta$left)]-time[min(which(theta$left != 0))] left.total.angle <- theta$left[length(theta$left)] left.lin.velocity.max <- max(left.velocity*lmlength) left.acceleration.max <- max(left.acceleration) left.lin.acceleration.max <- max(left.acceleration*lmlength) #rightmandible right.duration <- time[which.max(theta$right)]-time[min(which(theta$right != 0))] right.total.angle <- theta$right[length(theta$right)] right.lin.velocity.max <- max(right.velocity2*rmlength) right.acceleration.max <- max(right.acceleration) right.lin.acceleration.max <- max(right.acceleration*rmlength) start.stop <- c(time[which.max(theta$left)],time[min(which(theta$left != 0))], time[which.max(theta$right)], time[min(which(theta$right != 0))]) strike.duration <- max(start.stop)-min(start.stop) stats.table <- data.frame(strike.duration, left.duration, left.total.angle, left.velocity.max, left.lin.velocity.max, left.acceleration.max, left.lin.acceleration.max, left.mass.specific.power, right.duration, right.total.angle, right.velocity.max, right.lin.velocity.max, right.acceleration.max, right.lin.acceleration.max, right.mass.specific.power, left.vmax.time, right.vmax.time) row.names(stats.table) <- filename #save graphs to file original.dir <- getwd() dir.create(paste(getwd(),"/", fileout, sep="")) setwd(paste(getwd(),"/", fileout, sep="")) postscript(file=paste(fileout,".graph.eps", sep=""), horizontal=T, paper="letter", onefile=T) plot(time, theta$left, xlab="Time (sec)", ylab="Angular Displacement (radians)", pch=16, col="blue", ylim=c(0,max(left.Psp$ysmth, right.Psp$ysmth))) lines(tt, predict(left.Psp, tt), lwd=2, col="blue") points(time, theta$right, pch=16, col="red") lines(tt, predict(right.Psp, tt), lwd=2, col="red") plot(tt, left.velocity, xlab="Time (sec)", ylab=expression("Angular Velocity" ~ (rad ~sec^{-1})), pch=NA, ylim=c(0,max(left.velocity, right.velocity))) lines(tt, left.velocity, lwd=2, col="blue") lines(tt, right.velocity, lwd=2, col="red") plot(tt, left.acceleration, xlab="Time (sec)", ylab=expression("Angular acceleration" ~ (rad ~ sec^{-2})), pch=NA, ylim=c(min(left.acceleration, right.acceleration),max(left.acceleration, right.acceleration))) lines(tt, left.acceleration, lwd=2, col="blue") lines(tt, right.acceleration, lwd=2, col="red") dev.off() write.csv(stats.table, file=paste(filename,"stats.csv")) #go back to original directory setwd(original.dir) } write.table(stats.table, file="summary.csv", col.names=F, append=TRUE, sep="," ) }