SingleSet_Measurements01.py

from ij.plugin.filter import GaussianBlur
from fiji.threshold import Auto_Local_Threshold as ALT
from ij.plugin.filter import ParticleAnalyzer as PA
from org.jfree.data.statistics import BoxAndWhiskerCalculator
from java.util import ArrayList, Arrays
import os

# size of juxtanuclear region. In pixels.
RIMSIZE = 15
# image background is expected to be black. 
Prefs.blackBackground = True
# verbose output
VERBOSE = False

root = '/Volumes/data/bio-it_centres_course/data/course'
filedapi = '--W00005--P00001--Z00000--T00000--dapi.tif'

nucpath = os.path.join(root, filedapi)

impN = IJ.openImage(nucpath)

class InstBWC(BoxAndWhiskerCalculator):
    def __init__(self):
        pass
        
def getOutlierBound(rt):
  """ Analyzes the results of the 1st partcile analysis.
  Since the dilation of nuclear perimeter often causes 
  overlap of neighboring neculeus 'terrirories', such nucleus 
  are discarded from the measurements. 

  Small nucelei are already removed, but since rejection of nuclei depends on 
  standard outlier detection method, outliers in both smaller and larger sizes
  are discarded. 
  """
  area = rt.getColumn(rt.getColumnIndex('Area'))
  circ = rt.getColumn(rt.getColumnIndex("Circ."))
  arealist = ArrayList(Arrays.asList(area.tolist()))
  circlist = ArrayList(Arrays.asList(circ.tolist()))
  bwc = InstBWC()
  ans = bwc.calculateBoxAndWhiskerStatistics(arealist)
  #anscirc = bwc.calculateBoxAndWhiskerStatistics(circlist)
  if (VERBOSE):
    print ans.toString()
    print ans.getOutliers()
  q1 = ans.getQ1()
  q3 = ans.getQ3()
  intrange = q3 - q1 
  outlier_offset = intrange * 1.5
  return q1, q3, outlier_offset
  
def backgroundSubtraction(imp):
  """ subtract background, Cihan's method.
    see Simpson(2007)
  """
  impstats = imp.getProcessor().getStatistics()
  imp.getProcessor().setThreshold(impstats.min, impstats.mean, ImageProcessor.RED_LUT)
  measOpt = ImageStatistics.MEAN + ImageStatistics.LIMIT
  impstats = ImageStatistics.getStatistics(imp.getProcessor(), measOpt, None)
  backlevel = impstats.mean
  imp.getProcessor().resetThreshold()
  imp.getProcessor().subtract(backlevel)
  print imp.getTitle(), " : background intensity - ", backlevel
  return backlevel

def nucleusSegmentation(imp2):
    """ Segmentation of nucleus image. 
    Nucleus are selected that:
    1. No overlapping with dilated regions
    2. close to circular shape. Deformed nuclei are rejected.
    Outputs a binary image.
    """
#Convert to 8bit
    ImageConverter(imp2).convertToGray8()
#blur slightly using Gaussian Blur 
    radius = 2.0
    accuracy = 0.01
    GaussianBlur().blurGaussian( imp2.getProcessor(), radius, radius, accuracy)
# Auto Local Thresholding
    imps = ALT().exec(imp2, "Bernsen", 15, 0, 0, True)
    imp2 = imps[0]


#ParticleAnalysis 0: prefiltering by size and circularity
    rt = ResultsTable()
    paOpt = PA.CLEAR_WORKSHEET +\
                    PA.SHOW_MASKS +\
                    PA.EXCLUDE_EDGE_PARTICLES +\
                    PA.INCLUDE_HOLES #+ \
#		PA.SHOW_RESULTS 
    measOpt = PA.AREA + PA.STD_DEV + PA.SHAPE_DESCRIPTORS + PA.INTEGRATED_DENSITY
    MINSIZE = 20
    MAXSIZE = 10000
    pa0 = PA(paOpt, measOpt, rt, MINSIZE, MAXSIZE, 0.8, 1.0)
    pa0.setHideOutputImage(True)
    pa0.analyze(imp2)
    imp2 = pa0.getOutputImage() # Overwrite 
    imp2.getProcessor().invertLut()
    impNuc = imp2.duplicate()	## for the ring. 
    #impNuc = Duplicator().run(imp2)

#Dilate the Nucleus Area
## this should be 40 pixels in Cihan's method, but should be smaller. 
    rf = RankFilters()
    rf.rank(imp2.getProcessor(), RIMSIZE, RankFilters.MAX)

#Particle Analysis 1: get distribution of sizes. 

    paOpt = PA.CLEAR_WORKSHEET +\
                    PA.SHOW_NONE +\
                    PA.EXCLUDE_EDGE_PARTICLES +\
                    PA.INCLUDE_HOLES #+ \
#		PA.SHOW_RESULTS 
    measOpt = PA.AREA + PA.STD_DEV + PA.SHAPE_DESCRIPTORS + PA.INTEGRATED_DENSITY
    rt1 = ResultsTable()
    MINSIZE = 20
    MAXSIZE = 10000
    pa = PA(paOpt, measOpt, rt1, MINSIZE, MAXSIZE)
    pa.analyze(imp2)
    #rt.show('after PA 1')
#particle Analysis 2: filter nucleus by size and circularity. 
    #print rt1.getHeadings()
    if (rt1.getColumnIndex('Area') > -1):
      q1, q3, outlier_offset = getOutlierBound(rt1)
    else:
      q1 = MINSIZE
      q3 = MAXSIZE
      outlier_offset = 0
      print imp2.getTitle(), ": no Nucleus segmented,probably too many overlaps"

    paOpt = PA.CLEAR_WORKSHEET +\
                    PA.SHOW_MASKS +\
                    PA.EXCLUDE_EDGE_PARTICLES +\
                    PA.INCLUDE_HOLES #+ \
#		PA.SHOW_RESULTS 
    rt2 = ResultsTable()
    pa = PA(paOpt, measOpt, rt2, q1-outlier_offset, q3+outlier_offset, 0.8, 1.0)
    pa.setHideOutputImage(True)
    pa.analyze(imp2)
    impDilatedNuc = pa.getOutputImage() 

#filter original nucleus

    filteredNuc = ImageCalculator().run("AND create", impDilatedNuc, impNuc)
    return filteredNuc    
      
imp2 = impN.duplicate()
impfilteredNuc = nucleusSegmentation(imp2)
impfilteredNuc.show()