functions for 'mdf' files which stores digitized kinematic data from MTrackJ — Motor control of ultrafast, ballistic movements (by K. Kagaya and S. N. Patek)

Kagaya_Patek_kinem_control_hsv.R

## functions for 'mdf' files which stores digitized kinematic data from MTrackJ
## The paper's title: Motor control of ultrafast, ballistic movements (by K. Kagaya and S. N. Patek)
## This script is written by Katsushi Kagaya (k.kagaya@me.com)
## Gist URL: https://gist.github.com/kagaya/540596f6d1191c157cb7


library(ggplot2)
library(plyr)


read.one.mdf <- function(file) {
    ## mdf — output file of MTrackJ
    ## MTrackJ — a plugin of ImageJ useful for digitizing
    ## return 'track' — format for digitized points
    track.df <- read.table(file, skip = 4, fill = T)
    d <- data.frame()
    for (i in 1:dim(track.df)[1]) {
        line <- track.df[i,]
        if (line[1] == "Track") {
            track.id <- as.numeric(as.character(line[1, 2]))
            next
        }
        if (line[1] == "Point") {
            point.id <- as.numeric(as.character(line[1, 2]))
            point.x  <- as.numeric(as.character(line[1, 3]))
            point.y  <- as.numeric(as.character(line[1, 4]))
            frame.id <- as.numeric(as.character(line[1, 6]))
            d <- rbind(
                d, data.frame(
                    track.id = track.id,
                    point.id = point.id,
                    point.x = point.x, point.y = point.y,
                    frame.id = frame.id
                )
            )
        }
    }
    return(d)
}


get.rotation <- function(track, target.v = c(1,2), ref.v = c(1,3)) {
    ## current method to obtain angle made by two lines
    ## The the information for the two lines must be specified by four points.
    ## The four points are consisted of two points for one target line ("target.v"), and
    ## the other two points for one reference line ("ref.v").
    focus.id <- c(target.v, ref.v)
    trackt <- data.frame()
    for (i in focus.id) {
        trackt <- rbind(trackt, subset(track, track.id == i))
    }
    dmin <-
        ddply(trackt, .(track.id), subset, frame.id == min(frame.id))
    dmax <-
        ddply(trackt, .(track.id), subset, frame.id == max(frame.id))
    start.frame <- max(dmin$frame.id)
    end.frame <- min(dmax$frame.id)
    d <-
        subset(track, frame.id >= start.frame & frame.id <= end.frame)
    
    p1 <-
        subset(d, track.id == target.v[1], select = c("point.x", "point.y"))
    p2 <-
        subset(d, track.id == target.v[2], select = c("point.x", "point.y"))
    p3 <-
        subset(d, track.id == ref.v[1], select = c("point.x", "point.y"))
    p4 <-
        subset(d, track.id == ref.v[2], select = c("point.x", "point.y"))
    
    v1 <- p2 - p1
    v0 <- p4 - p3
    
    theta <- atan2(v1$point.y, v1$point.x) * 180 / pi + 180
    theta0 <- atan2(v0$point.y, v0$point.x) * 180 / pi + 180
    
    ## for(i in 1:length(theta)){
    ##     if(theta[i] > 270 & theta[i] < 360){
    ##         theta[i] <- theta[i] - 360
    ##     }
    ## }
    
    ## for(i in 1:length(theta0)){
    ##     if(theta0[i] > 270 & theta0[i] < 360){
    ##         theta0[i] <- theta0[i] - 360
    ##     }
    ## }
    
    rot <- theta - theta0
    rot <- rot - rot[1]
    
    return(rot)
}

get.rotation2 <- function(track, arm.length = 1, focal.p = 1) {
    # previous method to obtain rotational displacement
    ## nb504 case: 'a_track_of_nb504', 212.18, 1 # 212.18 = nb504$prop.len * 1000 / nb504$mmperpix
    # return propodus rotation
    # theta = l / r
    trackt <- subset(track, track.id == focal.p)
    theta = sqrt(diff(trackt$point.x) ^ 2 + diff(trackt$point.y) ^ 2) / arm.length
    phi = theta * 180 / pi
    rot2 <- phi[1]
    for (i in 2:length(phi)) {
        rot2 <- c(rot2, rot2[i - 1] + phi[i])
    }
    return(rot2)
}


check.rotation <- function(animal = nb504, strike = 10,
                           arm.length = 212.18, focal.p = 1,
                           arm.length.m = 80.63, focal.p.m = 3,
                           vlines = c(1, 21, 36, 46, 55)) {
    ## compare the two methods of digitizing (Fig.9B)
    ## nb504 case, arm.length=212.18 (pix)
    ## focal.p: focal moving (rotating) point, nb504 case, focal.p=1
    track <- animal$strike$raw[[strike]]
    fps <- animal$fps.strike
    r1 <-
        animal$strike$output.propodus.rotation[[strike]] # propodus current
    r.m1 <-
        -animal$strike$output.meralV.rotation[[strike]] # meral-V current
    ## propodus previous
    r2 <-
        c(0, get.rotation2(
            track = track, arm.length = arm.length, focal.p = focal.p
        ))
    ## meral-V previous
    r.m2 <-
        c(0,-get.rotation2(
            track = track, arm.length = arm.length.m, focal.p = focal.p.m
        ))
    time.stamp <-
        seq(0, 1 / fps * (length(r1) - 1), 1 / fps) * 1000 # s to ms
    d <-
        data.frame(
            calc.method = c(
                rep("propodus1", length(r1)), # prpodus current
                rep("propodus2", length(r2)), # propodus previous
                rep("meralV1", length(r.m1)), # meralV current
                rep("meralV2", length(r.m2))
            ), # meralV previous
            rotation = c(r1, r2, r.m1, r.m2),
            time = c(time.stamp, time.stamp, time.stamp, time.stamp)
        )
    
    vtimes <- 1 / fps * (vlines - 1) * 1000
    
    ggplot(d) +
        geom_point(aes(time, rotation,
                       col = calc.method, shape = calc.method), size = 2) +
        #            geom_line(aes(time, rotation, lty=calc.method)) +
        geom_hline(aes(yintercept = 20), lty = 2, lwd = 0.2) +
        geom_hline(aes(yintercept = 0), lty = 2, lwd = 0.2) +
        geom_vline(aes(xintercept = vtimes), data = data.frame(vtimes =
                                                                   vtimes), lwd = 0.2) +
        xlab("time (ms)") +
        ylab("rotation (degrees)") +
        scale_colour_manual(values = c(2, 2, 1, 1)) +
        scale_linetype_manual(values = c(1,1,2,2)) +
        scale_shape_manual(values = c(1,3,1,3)) +
        my.theme()
}


plot.track <- function(track) {
    ## to show digitized point
    ggplot(track) +
        geom_point(aes(point.x, point.y,
                       color = factor(track.id))) + theme_bw() + coord_equal()
}


my.theme <- function() {
    ## to custamize the ggplot for publication
    ## usage: q + my.theme()
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()
    ) +
        theme(
            axis.line = element_line(),
            axis.text = element_text(colour = "black"),
            axis.ticks = element_line(colour = "black")
        )
}