Commit 007e99c8 authored by Chris Evenhuis's avatar Chris Evenhuis
Browse files

Initial commit

parents
# Biofilm analysis
## Summary
This script analysis the height of a biofilm and correlated changes in the biofilm height to localisations found on another channel.
## Creation of the height map
### Filtering of the biofilm channel
1. Median filter of 20px (4.14 $\mu$m)
2. Gaussian Filter of 15px (3.11 $\mu$m)
3. Threshold method `Default` in ImageJ
4. Height map built up by copy each layer by copying (with 0 transparent) to a 2D projection.
![](media/map_3D_comp.png)
### Particle Localisation
For each imageL
1. Median filter of 2px (0.41 $\mu$m)
2. Gaussian Filter of 1px (0.20 $\mu$m)
3. Threshold method `Otsu` in ImageJ
4. Connected components labelling in 3D to identify particles
5. Centroid, Bounding box (for min and max Z), volume and mean Intensity calculated for each particle roi.
6. Results saved to `$filename_res.csv`.
This fails for the control as there's only background fluorescence.
Global threshold used by comparing the thresholds across the images
| file | threshold for particle |
|-----:|------------------------|
| 3hr / wild | 1277 |
| nAg | 1931 |
| gent | 1916 |
| 6hr / wild | 2158 |
| nAg | 1646 |
| gent | 1255 |
| Global | **1800** |
### Particle projection to create 2D RoIs
1. Particle connected components projected (max).
2. Rois extracted from label image
3. Average biofilm height (from 2D projection) calculated for each 2D RoI.
### Summary statistics for image
The height map is split to the regions that have a particle localisation (`roi`) and the regions that do not (`back`).
Summary statistics (mean and standard deviation) of the heigh map is analysed as follows:
* `tot_mean`, `tot_stdev ` : mean of the biofilm height (whole image) and it’s standard deviation
* `back_mean`, `back_stdev` : mean of the biofilm height for the parts that are NOT a detection on the purple channel and it’s standard deviation
* `roi_mean`, `roi_stdev` : mean of the biofilm height for the parts that are a detection on the purple channel and it’s
![](media/summary_stats_split.png)
{
"cells": [],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 4
}
This diff is collapsed.
import sys
sys.path.insert(0,"/Users/evenhuis/Dropbox/MIF/Workflows/biofilm_map/biofilm_map")
def add_2Dmap_overlay(ov,imp,rt):
from ij.gui import Overlay,EllipseRoi
from math import sqrt
[width, height, nChannels, nSlices, nFrames]=imp.getDimensions()
cal = imp.getCalibration()
pixelDepth = cal.pixelDepth
px2um = cal.pixelWidth
from java.awt import Color
reds='''255,255,178,15
254,217,118,50
254,178,76,100
253,141,60,150
240,59,32,200
189,0,38,255'''
colors=reds.split('\n')
rgbs = [ map(int,rgb.strip().split(',')) for rgb in colors ]
alpha=255
colors_f = [Color(rgb[0],rgb[1],rgb[2],rgb[3]) for rgb in rgbs ]
# work out if this is a 2D or 3D resuts table
is3D = ( 'r3' in rt.getHeadings() )
means =rt.getColumnAsDoubles(rt.getColumnIndex('Mean'))
max_mean=max(means)
from math import floor
for i in range(rt.getCounter()):
rel_mean = (rt.getValue('Mean',i)/max_mean)**(0.5)
imean = min(len(colors)-1,int(floor(rel_mean*len(colors))))
x = rt.getValue('X',i)/px2um
y = rt.getValue('Y',i)/px2um
z = rt.getValue('Z',i)/pixelDepth
if( is3D ):
rz = rt.getValue('r1',i)/pixelDepth
r2 = rt.getValue('r2',i)/px2um
r3 = rt.getValue('r3',i)/px2um
rxy = sqrt(r2*r2+r3*r3)
else:
r1 = rt.getValue('r1',i)/px2um
r2 = rt.getValue('r2',i)/px2um
rxy = sqrt(r1*r1+r2*r2)
#rxy=4
islice = int( z)
rad_p = rxy
roi = EllipseRoi(x-rad_p,y-rad_p,x+rad_p,y+rad_p,1)
print(rt.getValue('Mean',i),max_mean,imean)
print(colors_f[imean])
roi.setStrokeColor(colors_f[imean])
ov.add(roi)
return ov
def main():
from ij import IJ, ImagePlus
imp = IJ.getImage()
w,h,nchan,nslices = imp.getDimensions()[:4]
from ij.plugin import Duplicator
dup = Duplicator()
imp_film = dup.run(imp, 2,2, 1,nslices, 1,1)
stk = imp_film.getImageStack()
method='median2D'
if( method=='gauss3D' ):
from ij.plugin import GaussianBlur3D
GaussianBlur3D().blur( imp_film, 15,15,1)
if(method=='gauss2D'):
for i in range(1,nslices+1):
ip = stk.getProcessor(i)
ip.blurGaussian(15)
if(method=='median2D'):
from ij.plugin.filter import RankFilters
rf = RankFilters()
for i in range(1,nslices+1):
w2 = w//2 ; h2 = h//2
ip = stk.getProcessor(i).resize(w2,h2,True)
rf.rank(ip,10,rf.MEDIAN)
ipt =ip.resize(w,h,True)
ipt.blurGaussian(15)
stk.setProcessor(ipt,i)
IJ.setAutoThreshold(imp_film, "Default dark stack");
IJ.run(imp_film, "Convert to Mask", "method=Default background=Dark black");
import createHeightMap
ip_map = createHeightMap.createHeightMap(imp_film)
imp_map=ImagePlus("map",ip_map)
from ij.plugin import LutLoader
lut_dir = IJ.getDirectory("luts")
lut = LutLoader.openLut(lut_dir+"Cyan Hot.lut")
imp_map.setLut(lut)
imp_map.setDisplayRange(0,15)
from ij.gui import Overlay
dots=False
if( dots ):
# get the dots on the purlple channel
imp_nag = dup.run(imp, 1,1, 1,nslices, 1,1)
stk = imp_nag.getImageStack()
if( False ):
from ij.plugin.filter import RankFilters
rf = RankFilters()
for i in range(1,nslices+1):
ip = stk.getProcessor(i)
rf.rank(ip,2,rf.MEDIAN)
ip.blurGaussian(3)
import seg3Dmax
rt = seg3Dmax.spotPreview(imp_nag,70,mode='3D',ch=1)
ov = Overlay()
seg3Dmax.add_spot_overlay(ov, imp,rt)
imp.setOverlay(ov)
ov_map = Overlay()
add_2Dmap_overlay(ov_map,imp,rt)
else: # blob approach
import nAgLocalisation
imp_nag = dup.run(imp, 1,1, 1,nslices, 1,1)
rt_results,roim_p,ip_lab = nAgLocalisation.get_nAg_results_and_roim(imp_nag,thresh=1800)
from ij.plugin.frame import RoiManager
roim = RoiManager.getInstance() or RoiManager()
roim.reset()
from ij.gui import TextRoi
from java.awt import Font,Color
font = Font("Helvetica", Font.PLAIN, 8)
ov = Overlay()
cal = imp.getCalibration()
pixelWidth = cal.pixelWidth
pixelDepth = cal.pixelDepth
for i in range(roim_p.getCount()):
roi = roim_p.getRoi(i)
ip_map.setRoi(roi)
stats = ip_map.getStatistics()
rt_results.setValue("Ave film height",i,stats.mean*pixelDepth)
rt_results.setValue("CentX",i,stats.xCentroid*pixelWidth)
rt_results.setValue("CentY",i,stats.yCentroid*pixelWidth)
roim.addRoi( roi )
ov.add(roi)
x = rt_results.getValue("Centroid.X",i)
y = rt_results.getValue("Centroid.Y",i)
z = rt_results.getValue("Centroid.Z",i)
text = TextRoi(int(x/pixelWidth),int(y/pixelWidth),"{}".format(i+1))
text.setPosition(-1,-1,-1)
text.setStrokeColor(Color.RED)
text.setFont(font)
ov.add(text)
rt_results.show("Particle Results")
# stats on height map
stats = ip_map.getStatistics()
hist = stats.getHistogram()
from ij.measure import ResultsTable
rt_depth = ResultsTable.getResultsTable("Depth") or ResultsTable()
import depthStats
info=imp.getOriginalFileInfo()
if( info is None ):
name = "unknown/"+imp.getTitle()
else:
dirname=str(info.directory)
filename=str(info.fileName)
print(dirname,filename)
lastdir=dirname.split('/')[-1]
name = lastdir+'/'+filename
import os
root = os.path.splitext(filename)[0]
rt_results.save(dirname+"/"+root+"_res.csv")
print("saved results")
depthStats.add_depth_stats(name,ip_map,ip_lab,rt_depth)
rt_depth.show("Depth")
ImagePlus("lab",ip_lab).show()
imp_map.setCalibration(imp.getCalibration())
imp_map.setOverlay(ov)
imp_map.show()
if(__name__=='__main__'):
main()
\ No newline at end of file
from __future__ import print_function
from ij import IJ,ImagePlus
def remove_unattached_regions(imp):
'''removes regions from a binary Z-stack that don't touch the bottom slice'
INPUT a binary Z-stack
OUTPUT a binary Z-stack'''
# fill the bottom layer completely (set to 255)
# connected components labelling
# get the label of the bottom slice.
# keep only this one
pass
return
def createHeightMap(imp):
'''create a height / contour map from a binary Z-stack
INPUT: imp black and white (thresholded Z stack - 1 chan
OUTPUT: a 2D project with the maximum height for each pixel
'''
w,h,nchan,nslices = imp.getDimensions()[:4]
stk = imp.getImageStack()
ip_sum = stk.getProcessor(1).duplicate()
ip_sum.multiply(1/255.)
from ij.process import Blitter
for i in range(2,nslices+1):
ipt = stk.getProcessor(i).duplicate()
ipt.multiply(i/255.)
ip_sum.copyBits(ipt,0,0,Blitter.COPY_ZERO_TRANSPARENT)
return ip_sum
if( __name__=='__main__'):
imp = IJ.getImage()
ip_sum = createHeightMap(imp)
ImagePlus("D",ip_sum).show()
\ No newline at end of file
from __future__ import print_function
import sys
sys.path.insert(0,"/Users/evenhuis/Dropbox/MIF/Workflows/biofilm_map/biofilm_map")
def add_depth_stats( filename, imp_map, ip_lab, rt_depth,\
showSplit=False ):
from ij import ImagePlus
from ij.plugin.filter import ThresholdToSelection
t2s = ThresholdToSelection()
ip_map = imp_map.getProcessor()
ic = rt_depth.getCounter()
ip_map.setRoi(None)
stats = ip_map.getStatistics()
hist_tot = stats.getHistogram()
rt_depth.setValue("file",ic,filename)
print(stats.area,stats.mean,stats.stdDev)
rt_depth.setValue("tot_area", ic,stats.area)
rt_depth.setValue("tot_mean", ic,stats.mean)
rt_depth.setValue("tot_stdev",ic,stats.stdDev)
ip_lab.setThreshold(0,0,ip_lab.NO_LUT_UPDATE)
roi_back = t2s.convert(ip_lab)
ip_lab.setThreshold(1,ip_lab.getMax(),ip_lab.NO_LUT_UPDATE)
roi_part = t2s.convert(ip_lab)
ip_map.setRoi(roi_back)
stats = ip_map.getStatistics()
hist_back = stats.getHistogram()
rt_depth.setValue("back_area", ic,stats.area)
rt_depth.setValue("back_mean", ic,stats.mean)
rt_depth.setValue("back_stdev",ic,stats.stdDev)
print(stats.area,stats.mean,stats.stdDev)
if( showSplit ):
from java.awt import Color
from ij.gui import Overlay
ip_back = ip_map.duplicate()
ip_back.setRoi(roi_part)
ip_back.setValue(0)
ip_back.fill(roi_part.getMask())
imp_back = ImagePlus("back",ip_back)
ov = Overlay()
roi_back.setFillColor(Color.GRAY)
roi_part.setFillColor(Color.GRAY)
ov.add(roi_part)
imp_back.setOverlay(ov)
imp_back.show()
ip_map.setRoi(roi_part)
stats = ip_map.getStatistics()
hist_roi = stats.getHistogram()
rt_depth.setValue("roi_area", ic,stats.area)
rt_depth.setValue("roi_mean", ic,stats.mean)
rt_depth.setValue("roi_stdev",ic,stats.stdDev)
print(stats.area,stats.mean,stats.stdDev)
if( showSplit ):
ip_roi = ip_map.duplicate()
ip_roi.setRoi(roi_back)
ip_roi.setValue(0)
ip_roi.fill(roi_part.getMask())
imp_roi = ImagePlus("roi",ip_roi)
ov = Overlay()
ov.add(roi_back)
imp_roi.setOverlay(ov)
imp_roi.show()
#for i,(r,b,t) in enumerate(zip(hist_roi,hist_back,hist_tot[:20])):
# print(("{},"*4).format(i,float(r),float(b),float(t)))
def main():
from ij import IJ, ImagePlus, WindowManager
imp_map = WindowManager.getImage("map")
ip_map = imp_map.getProcessor()
ip_lab = WindowManager.getImage("lab").getProcessor()
from ij.measure import ResultsTable
rt_depth = ResultsTable()
add_depth_stats("R",imp_map,ip_lab,rt_depth,showSplit=True)
rt_depth.show("Depth")
if(__name__=='__main__'):
main()
\ No newline at end of file
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def get_adj_matrix( ip ):
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
from inra.ijpb.label import RegionAdjacencyGraph
from inra.ijpb.label import LabelImages
labels=LabelImages.findAllLabels(ip)
A_dict=dict( [ (i, set()) for i in labels ] )
adj=RegionAdjacencyGraph.computeAdjacencies(ip)
for kv in adj:
l1 = int(kv.label1)
l2 = int(kv.label2)
A_dict[l1].add(l2)
A_dict[l2].add(l1)
return A_dict
def getRoisFromIsolatedLabelMap( ip ):
from ij import ImagePlus
ip.setThreshold(1,ip.maxValue(),ip.NO_LUT_UPDATE)
ip_mask = ip.createMask()
from ij.measure import Measurements, ResultsTable
from ij.plugin.filter import ParticleAnalyzer
from ij.plugin.frame import RoiManager
from java.lang import Double
roim = RoiManager(False)
rt = ResultsTable()
measure_int = Measurements.MEDIAN + Measurements.CENTROID +Measurements.ELLIPSE + Measurements.AREA+Measurements.SHAPE_DESCRIPTORS
pa = ParticleAnalyzer(ParticleAnalyzer.SHOW_NONE+ParticleAnalyzer.INCLUDE_HOLES+ParticleAnalyzer.FOUR_CONNECTED, measure_int , rt, 0, Double.POSITIVE_INFINITY, 0.0, 1.0)
pa.setRoiManager(roim)
pa.analyze(ImagePlus("",ip_mask),ip)
roim_inst = RoiManager.getInstance()
if( not roim_inst): roim_inst=RoiManager()
roim_inst.reset()
rois = roim.getRoisAsArray()
indices=map(int,rt.getColumn(rt.getColumnIndex('Median')))
#isort = sorted(range(len(indices)),key=indices.__getitem__)
return rois,indices
def isolateLabelsAndReturnConnected(ip):
''' converts a label image to a list of rois'''
from inra.ijpb.label import LabelImages
labels=LabelImages.findAllLabels(ip)
adj = get_adj_matrix(ip)
conn = { i:len(adj[i]) for i in labels }
labels_to_add = [ l for l,c in zip(labels,conn) if c==0 ]
for l in labels_to_add:
adj.pop(l)
# 1. add all isolated labels
labels_removed=[]
while True:
#print("sie of adj",len(adj))
if( len(adj)==0): break
conn = { i:len(adj[i]) for i in adj.keys() }
# find the maximum connected label
i_max = max(conn, key=conn.get)
labels_to_add.append(i_max)
# remove it the neigbours from the adjecncy matrix
neighs = adj[i_max]
for n in neighs:
labels_removed.append(n)
if( n in adj.keys() ) :
adj.pop(n)
adj.pop(i_max)
#print(labels_to_add)
#print(labels_removed)
ip_to_add = LabelImages.keepLabels(ip,labels_to_add)
ip_removed = LabelImages.keepLabels(ip,labels_removed)
#print(type(ip_to_add))
return ip_to_add,labels_to_add, ip_removed,labels_removed
#rois_add,indices_add = getRoisFromIsolatedLabelMap(ip_to_add)
#return rois_add,indices_add, ip_removed,labels_removed
def label_to_rois( ip ):
#print("pass 1")
ip_add,labels,ip_rm, lab_rm = isolateLabelsAndReturnConnected(ip)
rois ,labels = getRoisFromIsolatedLabelMap(ip_add)
#rois,labels,ip_rm, lab_rm = isolateLabelsAndReturnConnected(ip)
from java.awt import Color
colors = [ getattr(Color,col) for col in 'RED GREEN YELLOW BLUE MAGENTA CYAN'.split()]
ia=0
col = colors[ia]
while( len(lab_rm)>0):
ip_add,lab_add,ip_rm,lab_rm = isolateLabelsAndReturnConnected(ip_rm)
rois_add, lab_add = getRoisFromIsolatedLabelMap(ip_add)
#rois_add,lab_add,ip_rm,lab_rm = isolateLabelsAndReturnConnected(ip_rm)
rois = rois+rois_add
labels = labels+lab_add
#print(len(lab_rm))
#for roi,ind in zip(rois_add,lab_add):
# roi.setStrokeColor(colors[ia])
# rois.append(roi)
ia=(ia+1)%len(colors)
#print(col)
# sort the labels by index
isort = sorted(range(len(labels)),key=labels.__getitem__)
return [ rois[i] for i in isort ]
def label_to_roim(ip,roim=None,pre=None):
rois = label_to_rois(ip)
#rois = getRoisFromIsolatedLabelMap(ip)
if( roim is None ):
from ij.plugin.frame import RoiManager
roim = RoiManager.getInstance() or RoiManager()
else:
pass
#print(roim.getCount())
#roim.reset()
for i,roi in enumerate(rois):
if( pre ):
roi.setName("{}-{}".format(pre,i+1))
else:
roi.setName("{}".format(i+1))
roim.addRoi(roi)
def label_to_rois_slow(ip):
from inra.ijpb.label import LabelImages
labels=LabelImages.findAllLabels(ip)
rois = []
from ij.plugin.filter import ThresholdToSelection
t2s = ThresholdToSelection()
for l in labels:
ip.setThreshold(l,l,ip.NO_LUT_UPDATE)
rois.append( t2s.convert(ip) )
return rois
def label_to_roim_slow(ip,roim=None,pre=None):
rois = label_to_rois_slow(ip)
#rois = getRoisFromIsolatedLabelMap(ip)
if( roim is None ):
from ij.plugin.frame import RoiManager
roim = RoiManager.getInstance() or RoiManager()
else:
pass
#print(roim.getCount())
#roim.reset()
for i,roi in enumerate(rois):
if( pre ):
roi.setName("{}-{}".format(pre,i+1))
else:
roi.setName("{}".format(i+1))
roim.addRoi(roi)
def roim_to_label( ip_ex,roim):
from ij.process import ShortProcessor
from ij.gui import PolygonRoi
ip_seg = ip_ex.convertToShortProcessor()
ip_seg.multiply(0)
for i,roi in enumerate(roim.getRoisAsArray()):
ip_seg.setValue(i+1)
if( roi.getType() in [roi.RECTANGLE] ):
ip_seg.fill()
elif( roi.getType() in [roi.FREELINE, roi.POLYLINE,roi.LINE ] ):
roi.drawPixels(ip_seg)
else:
ip_seg.setRoi(roi)
ip_seg.fill(roi.getMask())
#print(roi.LINE,roi.FREELINE)
return ip_seg
def main():
from ij import IJ,ImagePlus
imp = IJ.getImage()
ip = imp.getProcessor()
if( True ):