Voxel-based analysis of images

Voxel-based DIC Analysis.R

library("plyr")
library("ggplot2")
library("grid")
library("png")

# change the path to read in the correct file
img.files<-Sys.glob(file.path("thirdflow","*.png"))
# get the parameters from the file name (parm.loc means 4th spot by _ in the name)
# functions for converting images back and forth
im.to.df<-function(in.img) {
    out.im<-expand.grid(x=1:nrow(in.img),y=1:ncol(in.img))
    out.im$val<-as.vector(in.img)
    out.im
}
parm.from.filename<-function(in.name,parm.loc=4) {
  last.name<-sub("^([^.]*).*", "\\1", basename(in.name))
  as.numeric(strsplit(last.name,"_")[[1]][parm.loc])
}
all.images<-ldply(img.files,
                  function(file.name) {data.frame(file.name=file.name,
                                                  time=parm.from.filename(file.name,1),
                                                  im.to.df(readPNG(file.name)))
                                       })
 
all.images$phase<-all.images$val>0
 
 
 

# show the two images

start.img<-subset(all.images,time==min(all.images$time))
next.time<-min(subset(all.images,time>start.img[1,"time"])$time)
next.img<-subset(all.images,time==next.time)
ggplot()+
  geom_raster(data=subset(start.img,phase),aes(x,y,fill=as.factor(time)),alpha=0.75)+
  geom_raster(data=subset(next.img,phase),aes(x,y,fill=as.factor(time)),alpha=0.75)+
  coord_equal()+
  labs(fill="time")+
  theme_bw(20)
  
  
#' Calculate the cross correlation 
#' @author Kevin Mader (kevin.mader@gmail.com)
#' Generates flow with given object count, frame count and randomness
#' the box and crop are introduced to allow for objects entering and 
#' leaving the field of view
#'
#' @param img.a is the starting or I_0 image
#' @param img.b is the destination or I_1 image
#' @param tr.x is the function transforming the x coordinate
#' @param 
#' 
cc.imfun<-function(img.a,img.b,tr.x=function(x,y) x,tr.y=function(x,y) y) {
  # get the positions in image a
  x.vals<-unique(img.a$x)
  y.vals<-unique(img.a$y)
  
  # transform image b
  tr.img.b<-img.b
  # round is used to put everything back on an integer lattice
  tr.img.b$x<-round(with(img.b,tr.x(x,y)))
  tr.img.b$y<-round(with(img.b,tr.y(x,y)))
  # count the overlapping pixels in the window to normalize
  tr.img.b<-subset(tr.img.b,
                   ((x %in% x.vals) & (y %in% y.vals))
                   )
  norm.f<-nrow(tr.img.b)
  if(norm.f<1) norm.f<-1
  
  # keep only the in phase objects
  tr.img.a<-subset(img.a,phase)
  tr.img.b<-subset(tr.img.b,phase)
  
  if (nrow(tr.img.a)>0 & nrow(tr.img.b)>0) {
    matches<-ddply(rbind(cbind(tr.img.a,label="A"),cbind(tr.img.b,label="B")),.(x,y),
                   function(c.pos) {
                     if(nrow(c.pos)>1) data.frame(e.val=1)
                     else data.frame(e.val=c())
                     })
    data.frame(e.val=sum(matches$e.val)/norm.f,count=nrow(matches),size=norm.f)
    } else {
      data.frame(e.val=0,count=0,size=norm.f)
    }
  
  }
 
cc.points<-expand.grid(vx=seq(-10,10,1),vy=seq(-10,10,1))
cc.img<-ddply(cc.points,.(vx,vy),function(c.pt) {
  tr.x<-function(x,y) (x+c.pt[1,"vx"])
  tr.y<-function(x,y) (y+c.pt[1,"vy"])
  cc.imfun(start.img,next.img,tr.x=tr.x,tr.y=tr.y)
},.parallel=T)
 
 
ggplot(cc.img,aes(vx,vy,fill=e.val))+
  geom_raster()+geom_density2d(data=subset(cc.img,e.val>0),aes(weight=e.val))+
  labs(x="u",y="v",fill="Correlation",title="Correlation vs R")+
  scale_fill_gradient2(high="red")+
  theme_bw(25)