Code for Birkedal to do slice by slice analysis

A ImageJ Analysis.ijm

// taken from https://gist.github.com/kmader/10871528
// Load the data in
imageStackDirectory=getDirectory("Select the folder where the images should loaded from");
run("Image Sequence...", "open="+imageStackDirectory+" sort use");
// crop the image
makeRectangle(1089, 993, 675, 702);
run("Crop");
// clear the meaningless pixel sizes and set them to 0
run("Properties...", "unit=pixel pixel_width=1 pixel_height=1 voxel_depth=1 global");
// set which metrics to record for each analysis
run("Set Measurements...", "area mean min centroid center perimeter bounding fit shape redirect=None decimal=3");
// Apply a threshold
setThreshold(15316, 38138);
setOption("BlackBackground", false);
run("Convert to Mask", "method=Default background=Light 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))
  strsplit(last.name,"_")[[1]][parm.loc]
}
file.list<-Sys.glob(file.path("z:/Data10/alberto/previews","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"),
                                  sample=parm.from.filename(in.name,parm.loc=2),
                                  block=parm.from.filename(in.name,parm.loc=3),
                                  slice=parm.from.filename(in.name,parm.loc=4))
# 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
# scale the area to um2
all.results$AreaScaled<-0.325^2*all.results$Area
# factors for the groups they used
all.results$Group<-as.factor(floor(as.numeric(sapply(all.results$sample,function(x) gsub("[^0-9]","",x)))/100))

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=Area,color=block))+
  geom_density()+
  scale_x_log10()+
  facet_wrap(~sample)+
  theme_bw()

ggplot(all.results,aes(x=Area,color=sample))+
  geom_density()+
  scale_x_log10()+
  theme_bw()

# the largest lacuna we want to keep has a volume of 1000um3
# the radius of that is cuberoot(1000*3/(4*pi))
max.radius<-(1000*3/(4*pi))^(0.33)
# that makes a circle radius pi rš2
max.area<-pi*max.radius^2
# the smallest meaningful lacunae is 10 pixels
lacuna.results<-subset(all.results,AreaScaled<max.area & Area>10)
ggplot(lacuna.results,
       aes(x=AreaScaled,color=block))+
  geom_density()+
  facet_wrap(~sample)+
  labs(x="Area (um2)")+
theme_bw()

# just the different groups
ggplot(lacuna.results,
       aes(x=AreaScaled,color=Group))+
  geom_density()+
  labs(x="Area (um2)")+
  theme_bw()

# just the different groups
ggplot(lacuna.results,
       aes(x=AreaScaled,color=sample,linetype=block))+
  geom_density()+
  facet_wrap(~Group)
  labs(x="Area (um2)")+
  theme_bw()

# just the objects that are in the size range
ggplot(lacuna.results,
       aes(x=sample,y=AreaScaled,color=block))+
  geom_boxplot()+
  theme_bw()

results.table<-ddply(lacuna.results,.(sample,block),function(c.group) {
  data.frame(count=nrow(c.group),
             mean.area.um2=mean(c.group$AreaScaled),
             mean.perimeter.px=mean(c.group$Perim.))
})

write.csv(results.table,"results.csv")