Load Stack of Images and Perform Analysis on Each Slice

A analysis.ijm

// Load the data in
imageStackDirectory=getDirectory("Select the folder where the images should loaded from");
run("Image Sequence...", "open="+imageStackDirectory+" sort");
run("Set Measurements...", "area mean min centroid center perimeter bounding fit shape redirect=None decimal=3");
// Apply a threshold
setThreshold(129, 255);
setOption("BlackBackground", false);
run("Convert to Mask", "method=Default background=Dark black");

imgDirectory=getInfo("image.directory");
// gather information about the current image
getDimensions(width, height, channels, slices, frames);

for(i=1;i<=slices;i++) { 
  setSlice(i);
  sliceName=getInfo("slice.label");
  run("Analyze Particles...", "display clear record slice");
  saveAs("Results", imgDirectory+File.separator+"shape_"+sliceName+".txt");
}

B read in images.R

library("plyr")
library("ggplot2")
# get the parameters from the file name (parm.loc means 4th spot by _ in the name)
parm.from.filename<-function(in.name,parm.loc=4) {
  last.name<-sub("^([^.]*).*", "\\1", basename(in.name))
  as.numeric(strsplit(last.name,"_")[[1]][parm.loc])
}
file.list<-Sys.glob(file.path("/Users/mader/Dropbox/TeachingRelated/QuantBig/Exercises/Ex9/firstflow","shape*.txt"))
# function to read each file and add a column for threshold
read.fcn<-function(in.name) cbind(read.csv(in.name,sep="\t"),
          time=parm.from.filename(in.name,parm.loc=2))
# read in all of the files
all.results<-ldply(file.list,read.fcn)
# change the names to com.x and com.y to keep things clear
all.results$com.x<-all.results$X.1
all.results$com.y<-all.results$Y

C show some plots.R

# show all the points first 
# note X is called X.1 but Y is just Y
ggplot(all.results,aes(x=com.x,y=com.y,color=time))+
  geom_point()+
  scale_color_gradientn(colours=rainbow(10))+
  theme_bw()

D run the tracking.R

track.fun<-function(in.results) {
  # run the function on time values
  # through n-1 (since the last time cant be tracked)
  ddply(subset(in.results,time<max(in.results$time)),.(time),
        function(cur.frame) {
          # get and save the current time
          c.time<-cur.frame[1,"time"]
          # calculate the next frame time
          next.time<-min(subset(in.results,time>c.time)$time)
          next.frame<-subset(in.results,time==next.time)
          print(next.time)
          ddply(cur.frame,.(com.x,com.y),function(c.pos) {
            cx<-c.pos[1,"com.x"]
            cy<-c.pos[1,"com.y"]
            r.dist<-with(next.frame,sqrt((com.x-cx)^2+(com.y-cy)^2))
            min.idx<-which.min(r.dist)
            # copy all of the information from the matched
            # point to the current point
            matched.pt<-next.frame[min.idx,,drop=F]
            names(matched.pt)<-sapply(names(matched.pt),function(in.name) paste("M_",in.name,sep=""))
            cbind(c.pos,matched.pt,dist=min(r.dist))
            })
          
          })
  }
  
  
  track.stats<-function(in.tracks) {
  ddply(in.tracks,.(time),function(cur.frame) {
    data.frame(Vx.Mean=with(cur.frame,mean(M_com.x-com.x)),
               Vx.Sd=with(cur.frame,sd(M_com.x-com.x)),
               Vy.Mean=with(cur.frame,mean(M_com.y-com.y)),
               Vy.Sd=with(cur.frame,sd(M_com.y-com.y)),
               Area.Change=with(cur.frame,mean(abs(M_Area-Area)/Area))
               )
  })
}


  
  
all.tracks<-track.fun(all.results)

track.stats(all.tracks)
  

E show some plots of the tracking.R

library(grid)
ggplot(all.tracks,aes(x=com.x,y=com.y,color=time))+
  geom_point()+
  # add the arrows to the next point
  geom_segment(aes(xend=M_com.x,yend=M_com.y),
               arrow=arrow(length = unit(0.3,"cm")))+
  scale_color_gradientn(colours=rainbow(10))+
  facet_wrap(~time)+
  theme_bw()

F Full Report Code.Rmd

Tracking Example
========================================================

This is an R Markdown document. Markdown is a simple formatting syntax for authoring web pages (click the **MD** toolbar button for help on Markdown).

When you click the **Knit HTML** button a web page will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

```{r}
library("plyr")
library("ggplot2")
# get the parameters from the file name (parm.loc means 4th spot by _ in the name)
parm.from.filename<-function(in.name,parm.loc=4) {
  last.name<-sub("^([^.]*).*", "\\1", basename(in.name))
  as.numeric(strsplit(last.name,"_")[[1]][parm.loc])
}
file.list<-Sys.glob(file.path("~/Downloads/firstflow","shape*.txt"))
# function to read each file and add a column for threshold
read.fcn<-function(in.name) cbind(read.csv(in.name,sep="\t"),
          time=parm.from.filename(in.name,parm.loc=2))
# read in all of the files
all.results<-ldply(file.list,read.fcn)
# change the names to com.x and com.y to keep things clear
all.results$com.x<-all.results$X.1
all.results$com.y<-all.results$Y
```

# Show a plot of the bubbles

```{r fig.width=7, fig.height=6}
# show all the points first 
# note X is called X.1 but Y is just Y
ggplot(all.results,aes(x=com.x,y=com.y,color=time))+
  geom_point()+
  scale_color_gradientn(colours=rainbow(10))+
  theme_bw()
```

# The tracking function in a named block

```{r tracking_function}
track.fun<-function(in.results) {
  # run the function on time values
  # through n-1 (since the last time cant be tracked)
  ddply(subset(in.results,time<max(in.results$time)),.(time),
        function(cur.frame) {
          # get and save the current time
          c.time<-cur.frame[1,"time"]
          # calculate the next frame time
          next.time<-min(subset(in.results,time>c.time)$time)
          next.frame<-subset(in.results,time==next.time)
          print(next.time)
          ddply(cur.frame,.(com.x,com.y),function(c.pos) {
            cx<-c.pos[1,"com.x"]
            cy<-c.pos[1,"com.y"]
            r.dist<-with(next.frame,sqrt((com.x-cx)^2+(com.y-cy)^2))
            min.idx<-which.min(r.dist)
            # copy all of the information from the matched
            # point to the current point
            matched.pt<-next.frame[min.idx,,drop=F]
            names(matched.pt)<-sapply(names(matched.pt),function(in.name) paste("M_",in.name,sep=""))
            cbind(c.pos,matched.pt,dist=min(r.dist))
            })
          
          })
  }
  
  
  track.stats<-function(in.tracks) {
  ddply(in.tracks,.(time),function(cur.frame) {
    data.frame(Vx.Mean=with(cur.frame,mean(M_com.x-com.x)),
               Vx.Sd=with(cur.frame,sd(M_com.x-com.x)),
               Vy.Mean=with(cur.frame,mean(M_com.y-com.y)),
               Vy.Sd=with(cur.frame,sd(M_com.y-com.y)),
               Area.Change=with(cur.frame,mean(abs(M_Area-Area)/Area))
               )
  })
}

  
all.tracks<-track.fun(all.results)
 

```

## Show some statistics on the tracking (asis allows the table to be formatted correctly)
```{r, results='asis'}
kable(track.stats(all.tracks))
```

# The tracking results

```{r, fig.width=7,fig.height=7}
library(grid)
ggplot(all.tracks,aes(x=com.x,y=com.y,color=time))+
  geom_point()+
  # add the arrows to the next point
  geom_segment(aes(xend=M_com.x,yend=M_com.y),
               arrow=arrow(length = unit(0.3,"cm")))+
  scale_color_gradientn(colours=rainbow(10))+
  facet_wrap(~time)+
  theme_bw()
```

G Voxel Based.R

library(png)
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




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 ([email protected])
#' 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)