Make MOT videos

generate.R

library(ggplot2)
library(animation)
library(circular)

# utility functions for trajectory generation -------------------
#   - these function came from motAnalysis package 

random.positions <- function(n, xlim = c(-10, +10), ylim = xlim,
                             dot.radius = 0.5, min.dist = 1) {
  pos <- list(
    n = n, xlim = xlim, ylim = ylim,
    x = runif(n, xlim[1], xlim[2]),
    y = runif(n, ylim[1], ylim[2]),
    radius = dot.radius
  )
  while (!valid.positions(pos, min.dist = min.dist)) {
    pos$x <- runif(n, xlim[1], xlim[2])
    pos$y <- runif(n, ylim[1], ylim[2])
  }
  class(pos) <- "positions"
  return(pos)
}

valid.positions <- function(pos, min.dist = 1) {
  distances <- dist(data.frame(pos$x, pos$y))
  distances <- as.numeric(distances)
  valid.dist <- all(distances >= min.dist)
  border.dist <-
    (pos$x >= pos$xlim[1] + min.dist) &
      (pos$x <= pos$xlim[2] - min.dist) &
      (pos$y >= pos$ylim[1] + min.dist) &
      (pos$y <= pos$ylim[2] + min.dist)
  border.dist <- all(border.dist)
  ok <- valid.dist & border.dist
  return(ok)
}

vonmises.trajectory <- function(pos, speed = 5, secs = 10, fps = 100,
                                frame.step = 10,
                                initial.dir = numeric(pos$n),
                                kappa = 50, sep = 1) {
  time <- seq(0, secs, 1 / fps)
  n.frames <- length(time)
  x <- matrix(0, nrow = n.frames, ncol = pos$n) * NA
  y <- matrix(0, nrow = n.frames, ncol = pos$n) * NA
  if (initial.dir[1] == "runif") {
    dir <- runif(pos$n, min = 0, max = 2 * pi)
  } else {
    dir <- rep(1, pos$n) * initial.dir
    # initial directions 0=up, pi/2=right (=CW)
  }
  # is there any other good default direction? => SEARCH
  step <- speed / fps
  for (f in 1:n.frames) {
    if (f == 1) { # first frame from input positions
      x[f, ] <- pos$x
      y[f, ] <- pos$y
      last.change <- 1
    } else {
      if (!is.na(frame.step) & (f > last.change + frame.step)) {
        # dir = runif(pos$n, min=0, max=2*pi)
        jitt <- rvonmises(n = pos$n, mu = circular(0), kappa = kappa)
        dir <- (dir + as.numeric(jitt)) %% (2 * pi)
        last.change <- f
      }
      dir <- bounce(x[f - 1, ], y[f - 1, ], dir, step, pos$xlim, pos$ylim)
      dir <- bounce.mutual(x[f - 1, ], y[f - 1, ], dir, step, sep = sep)
      x[f, ] <- x[f - 1, ] + sin(dir) * step
      y[f, ] <- y[f - 1, ] - cos(dir) * step
    }
  }
  track <- list(
    n = pos$n, xlim = pos$xlim, ylim = pos$ylim,
    time = time, x = x, y = y
  )
  class(track) <- "trajectory"
  return(track)
}

bounce <- function(x, y, dir, step, xlim = c(-10, +10), ylim = xlim) {
  dir2 <- dir
  x2 <- x + sin(dir) * step
  y2 <- y - cos(dir) * step
  too.R <- x2 > xlim[2]
  too.L <- x2 < xlim[1]
  too.U <- y2 < ylim[1]
  too.D <- y2 > ylim[2]
  corner <- (too.R | too.L) & (too.U | too.D)
  side <- xor(too.R | too.L, too.U | too.D)
  dir2[corner] <- (dir2[corner] + pi) %% (2 * pi)
  dir2[too.R | too.L] <- 2 * pi - dir2[too.R | too.L]
  dir2[too.U | too.D] <- (pi - dir2[too.U | too.D]) %% (2 * pi)
  # bounce.inspect(x,y,dir,dir2,xlim,ylim)
  return(dir2)
}

bounce.mutual <- function(x, y, dir1, step, sep = 1) {
  # je to tak? http://en.wikipedia.org/wiki/Elastic_collision
  n <- length(x)
  dir2 <- dir1
  next.x <- x + sin(dir1) * step
  next.y <- y - cos(dir1) * step
  ro <- matrix(rep(1:n, n), n, n)
  ro <- ro[lower.tri(ro)]
  co <- matrix(rep(1:n, each = n), n, n)
  co <- co[lower.tri(co)]
  di <- as.matrix(dist(data.frame(next.x, next.y)))
  di <- di[lower.tri(di)]
  # print(co); print(ro); print(di)
  colli <- which(di < sep)
  for (i in colli) {
    e1 <- ro[i]
    e2 <- co[i]
    dir2[c(e1, e2)] <- dir2[c(e2, e1)]
  }
  return(dir2)
}
# bounce.mutual(c(1,2,3), c(1,2,1), c(0,0,0), .1, sep=1.9)

snapshot.trajectory <- function(track, time, time.index = NA) {
  if (is.na(time.index)) {
    time.index <- which.min(abs(track$time - time))
  }
  pos <- list(
    n = track$n, xlim = track$xlim, ylim = track$ylim,
    x = track$x[time.index, ],
    y = track$y[time.index, ],
    radius = 0.5
  ) # TODO
  class(pos) <- "positions"
  return(pos)
}

# test generation code --------------------------------------------
set.seed(101)
xy <- random.positions(8)
xyt <- vonmises.trajectory(xy, initial.dir = 3.14 * (1:8) / 4)

plot(xy)
plot(xyt)


# display functions ---------------------------------------------
xplot.positions <- function(pos, targets = NA, labels = F,
                            legend = F, expand = c(0, 1)) {
  expand.add <- c(-1, +1) * expand[2]
  expand <- as.numeric(na.omit(c(expand, 0, 0))) # add one or two zeros
  width <- diff(pos$xlim)
  height <- diff(pos$ylim)
  margin.x <- width * expand[1]
  margin.y <- height * expand[1]
  xlim.new <- c(
    pos$xlim[1] - margin.x - expand[2],
    pos$xlim[2] + margin.x + expand[2]
  )
  ylim.new <- c(
    pos$ylim[1] - margin.y - expand[2],
    pos$ylim[2] + margin.y + expand[2]
  )
  n <- pos$n
  d <- data.frame(
    dot = factor(1:n), x = pos$x, y = pos$y,
    type = "dot", stringsAsFactors = F
  )
  d$type[d$type != "target"] <- "distractor"
  if (!any(is.na(targets))) {
    d$type[targets] <- "target"
    d$type[d$type != "target"] <- "distractor"
  }

  pp <- qplot(x, y,
    data = d,
    geom = "point", colour = type,
    asp = 1, size = I(10)
  ) +
    labs(
      x = "", y = "", title = "",
      colour = "", shape = ""
    ) +
    coord_cartesian(xlim = xlim.new, ylim = ylim.new) +
    scale_y_reverse() +
    scale_color_manual(values = c("#AAAAAA", "#00FF00")) +
    theme(
      rect = element_rect(fill = "white"), text = element_blank(),
      line = element_blank(),
      panel.background = element_rect(fill = "white"),
      plot.background = element_rect(fill = "white")
    )
  if (labels) {
    pp <- pp + geom_text(aes(label = dot),
      colour = I("black"), size = I(8)
    )
  }
  # better or _no_ legend
  pp <- pp + theme(legend.position = "none")
  return(pp)
}
# xplot.positions(xy, targets = 1:3)

plot.trajectory.x <- function(tr, legend = F, expand = c(1, 1)) {
  expand.add <- c(-1, +1) * expand[2]
  xlim.new <- tr$xlim * expand[1] + expand.add
  ylim.new <- rev(tr$ylim * expand[1] + expand.add)
  n <- tr$n
  d <- long.trajectory(tr)
  pp <- qplot(x, y,
    data = d,
    geom = "point", colour = factor(object), asp = 1
  ) +
    labs(x = "", y = "", title = "") +
    coord_cartesian(xlim = xlim.new, ylim = ylim.new) +
    scale_y_reverse()
  pp <- pp + theme(legend.position = "none")
  pp <- pp + theme()
  return(pp)
}
# plot.trajectory.x(xyt)

make.videox <- function(tt, fname, fps = 25,
                        outdir = getwd(), targets) {
  tmin <- min(tt$time)
  tmax <- max(tt$time)
  tlen <- tmax - tmin
  oopt <- ani.options(
    interval = 1 / fps, nmax = fps * tlen * 2,
    outdir = outdir, ani.width = 1920 / 2, ani.height = 1080 / 2
  )
  saveVideo({
    for (i in 1:(2 * fps)) {
      p <- snapshot.trajectory(tt, tmin)
      print(xplot.positions(p, targets = targets))
    }
    for (tim in seq(tmin, tmax, 1 / fps)) {
      p <- snapshot.trajectory(tt, tim)
      print(xplot.positions(p, targets = NA))
    }
    for (i in 1:(2 * fps)) {
      p <- snapshot.trajectory(tt, tmax)
      print(xplot.positions(p, targets = targets))
    }
  }, video.name = fname, other.opts = "-pix_fmt yuv420p -b 300k", clean = T)
  ani.options(oopt)
}

# making videos ------------------------------------------------------------

ani.options(ffmpeg = "CHANGE/Apps/ffmpeg/ffmpeg") # TODO

xy1 <- random.positions(8)
xyt1 <- vonmises.trajectory(xy1, speed = 5, initial.dir = "runif")
make.videox(xyt1, "video_s05_n1_01.mp4", targets = 1)

xy2 <- random.positions(8)
xyt2 <- vonmises.trajectory(xy2, speed = 10, initial.dir = "runif")
make.videox(xyt2, "video_s10_n1_01.mp4", targets = 1)

xy3 <- random.positions(8)
xyt3 <- vonmises.trajectory(xy3, speed = 5, initial.dir = "runif")
make.videox(xyt3, "video_s05_n4_01.mp4", targets = 1:4)

xy4 <- random.positions(8)
xyt4 <- vonmises.trajectory(xy4, speed = 10, initial.dir = "runif")
make.videox(xyt4, "video_s10_n4_01.mp4", targets = 1:4)

xy5 <- random.positions(8)
xyt5 <- vonmises.trajectory(xy5, speed = 10, initial.dir = "runif")
make.videox(xyt5, "video_s10_n4_02.mp4", targets = 1:4)

Removed dependency on motAnalysis

• utility functions from motAnalysis package copied to the beginning
• some cleanup