Autothresholding problem

scripting

#1

Dear ImageJ users,

With the following script I want to determine viral DNA(Channel 2)/RNA(Channel 3) counts for each nucleus (channel1). Therefore I am using an autothresholding method (Yen dark- no rest). However, the code gives very different RNA counts for images that look very similar. I attached the link for the two images and a screenshot of the results to demonstrate the problem.

Images: https://drive.google.com/open?id=1UvVqTHStbQ3SgCONWYqlyANnRi10ArT9

//@File(label="Input directory", description="Select the directory with OIB images", style="directory") inputDir
//@Integer(label="Nuclei Channel", value=1) nucleiChannel
//@Integer(label="DNA Channel", value=2) dnaChannel
//@Integer(label="RNA Channel", value=3) rnaChannel
//Double(label="Min Particle size in um", value=0.1) minParticleSizeInUnit
//String(label="Threshold Method for Tissue Detection",choices={'Moments','Huang','Percentile','Triangle','IJ_IsoData','Percentile'}) thresholdMethod

import ij.IJ
import ij.ImagePlus
import ij.gui.GenericDialog
import ij.Prefs
import ij.measure.Measurements
import ij.measure.Calibration
import ij.WindowManager

import ij.process.ImageConverter
import ij.process.ImageProcessor
import ij.process.ByteProcessor
import ij.process.ImageStatistics 

import ij.plugin.Duplicator
import ij.plugin.ImageCalculator
import ij.plugin.RoiScaler
import ij.plugin.ZProjector 
import ij.plugin.Thresholder
import ij.process.ImageStatistics
import ij.measure.ResultsTable
import ij.plugin.filter.ParticleAnalyzer
import ij.plugin.frame.RoiManager

import ij.gui.Roi
import ij.gui.ShapeRoi

import groovy.io.FileType
import groovy.time.TimeCategory 
import groovy.time.TimeDuration

import loci.plugins.BF
import loci.formats.ChannelSeparator
import loci.formats.FormatException
import loci.formats.IFormatReader
import loci.formats.meta.IMetadata
import loci.formats.MetadataTools

import loci.common.Region
import loci.plugins.in.ImporterOptions

import loci.plugins.util.ImageProcessorReader
import loci.plugins.util.LociPrefs
import loci.common.Region

import java.awt.Rectangle
import java.io.File
import java.lang.Double
import java.nio.file.Paths
import java.util.Collections

import static groovy.io.FileType.FILES

rm = new RoiManager(true)
rm.reset()
rt = new ResultsTable()

//Clear the log window
IJ.log("\\Clear")

folderName=inputDir.getAbsolutePath()
IJ.log('Analyze Dir '+inputDir.name)

minCellSize=1100
maxCellSize=Double.POSITIVE_INFINITY


minParticleSizedna=0
maxParticleSize=Double.POSITIVE_INFINITY

minParticleSizerna=0
maxParticleSize=Double.POSITIVE_INFINITY

fileNameList=[]
cellNumberList=[]
numberOfDNAParticleList=[]
numberOfRNAParticleList=[]
meanIntensityDNAParticleList=[]
areaDNAParticleList=[]
meanIntensityRNAParticleList=[]
areaRNAParticleList=[]
  
IJ.log('Start Analysis')
inputDir.eachFile (FILES) {  
  if(it.name.endsWith('.oif')) {
    IJ.log(it.name)
    def baseName=it.getAbsolutePath()[0..-5]
    def baseFileName=it.name[0..-5]
    rm.reset()

    def options = new ImporterOptions();
    options.setId(Paths.get(it.getAbsolutePath()).toString());
    options.setAutoscale(true);
    //options.setCrop(true);
    //options.setCropRegion(0, new Region(x, y, w. h));
    options.setColorMode(ImporterOptions.COLOR_MODE_COMPOSITE);
    def imps = BF.openImagePlus(options);
    imp = imps[0]
    //imp.show()
    //lll
    //IJ.run("Bio-Formats Importer", "open=["+Paths.get(it.getAbsolutePath()).toString()+"] color_mode=Composite rois_import=[ROI manager] view=Hyperstack stack_order=XYCZT");
    //imp = IJ.getImage()
    //Z Sum Projection
    imp = ZProjector.run(imp,"sum");
    //Duplicate the nuclei channel
    imp_nuclei=new Duplicator().run(imp, nucleiChannel.toInteger(), nucleiChannel.toInteger())
    //Detect the Nuclei and store the matching ROIs
    IJ.run(imp_nuclei, "Gaussian Blur...", "sigma=10");
    IJ.setAutoThreshold(imp_nuclei, "Otsu dark no-reset");
    Prefs.blackBackground = false;
    IJ.run(imp_nuclei, "Convert to Mask", "");
    IJ.run(imp_nuclei, "Watershed", "");
 
    particleAnalysisOptions = ParticleAnalyzer.ADD_TO_MANAGER + ParticleAnalyzer.EXCLUDE_EDGE_PARTICLES; 
    p  = new ParticleAnalyzer(particleAnalysisOptions, 0, rt, minCellSize, maxCellSize, 0.0, 1.0)
    p.setRoiManager(rm)
    p.analyze(imp_nuclei)
    def cellRois=rm.getRoisAsArray()
    rm.runCommand('Save', Paths.get(inputDir.getAbsolutePath(),baseFileName+'_cell_rois.zip').toString())
    //Then for each cell we will count the number of particle and take their intensity
    //Channel 2 is DNA
    thresholdMethod="Yen dark no-reset"
    imp_DNA=new Duplicator().run(imp, dnaChannel.toInteger(), dnaChannel.toInteger())
    imp_DNA_binary = new Duplicator().run(imp, dnaChannel.toInteger(), dnaChannel.toInteger())
    IJ.setAutoThreshold(imp_DNA_binary, thresholdMethod);
    Prefs.blackBackground = false;
    IJ.run(imp_DNA_binary, "Convert to Mask", "");
    //imp_DNA_binary.show()
    //lll
    //Channel 3 is RNA
    imp_RNA=new Duplicator().run(imp, rnaChannel.toInteger(), rnaChannel.toInteger())
    imp_RNA_binary = new Duplicator().run(imp, rnaChannel.toInteger(), rnaChannel.toInteger())
    IJ.setAutoThreshold(imp_RNA_binary, thresholdMethod);
    Prefs.blackBackground = false;
    IJ.run(imp_RNA_binary, "Convert to Mask", "");
    for(int i=0;i<cellRois.size();i++)
    {
      fileNameList.add(baseFileName)
      cellNumberList.add((i+1))
      rm.reset()
      roiCell=cellRois[i].clone()
      roiCell.setLocation(0, 0)
      //Measure on the DNA Channel
      imp_DNA.setRoi(cellRois[i])
      impCell = imp_DNA.duplicate();
      imp_DNA_binary.setRoi(cellRois[i])
      impCellMask=imp_DNA_binary.duplicate();
      //impCellMask.show()
      //lll
      (numberOfDNAParticleList,meanIntensityDNAParticleList,areaDNAParticleList)=analyzeCellParticle(impCell,impCellMask,roiCell,minParticleSizedna,maxParticleSize,
                                                                                 numberOfDNAParticleList,meanIntensityDNAParticleList,areaDNAParticleList,rm,"Yen dark no-reset")
      //Measure on the RNA Channel
      imp_RNA.setRoi(cellRois[i])
      impCell = imp_RNA.duplicate();
      imp_RNA_binary.setRoi(cellRois[i])
      impCellMask=imp_RNA_binary.duplicate();
      (numberOfRNAParticleList,meanIntensityRNAParticleList,areaRNAParticleList)=analyzeCellParticle(impCell,impCellMask,roiCell,minParticleSizerna,maxParticleSize,
                                                                                 numberOfRNAParticleList,meanIntensityRNAParticleList,areaRNAParticleList,rm,"Yen dark no-reset")

      /*
      IJ.setBackgroundColor(0, 0, 0);
      IJ.run(impCell, "Clear Outside", "");
      IJ.setAutoThreshold(impCell, "Yen dark no-reset");
      Prefs.blackBackground = false;
      IJ.run(impCell, "Convert to Mask", "");
      particleAnalysisOptions = ParticleAnalyzer.ADD_TO_MANAGER + 
ParticleAnalyzer.EXCLUDE_EDGE_PARTICLES; 
      p  = new ParticleAnalyzer(particleAnalysisOptions, 0, rt, minParticleSize, maxParticleSize, 0.0, 1.0)
      p.setRoiManager(rm)
      p.analyze(impCell)
      def particleRois=rm.getRoisAsArray()
      def meanParticleIntensity=0
      numberOfDNAParticleList.add(particleRois.size())
      for(int p=0;p<particleRois.size();p++)
      {
        impCell.setRoi(particleRois[p])
        stats=impCell.getStatistics(Measurements.MEAN+Measurements.MEDIAN+Measurements.AREA+Measurements.STD_DEV )
        meanParticlesIntensity+=stats.@mean
        areaParticles+=stats.@area
      }
      meanParticlesIntensity=meanParticlesIntensity/particleRois.size()
      meanIntensityDNAParticleList.add(meanParticlesIntensity)
      areaParticles=areaParticles/particleRois.size()
      areaDNAParticleList.add(areaParticles)      
    }   
    */
    }
  }
}
rt=populateResultTable(rt,fileNameList,cellNumberList,
                        numberOfDNAParticleList,meanIntensityDNAParticleList,areaDNAParticleList,
                        numberOfRNAParticleList,meanIntensityRNAParticleList,areaRNAParticleList)
rt.show('Result')
rt.save(Paths.get(inputDir.getAbsolutePath(), 'results.xls').toString())
IJ.log('Result saved in '+Paths.get(inputDir.getAbsolutePath(), 'results.xls').toString())
IJ.log('Analysis Done!')

def populateResultTable(rt,fileNameList,cellNumberList,
                        numberOfDNAParticleList,meanIntensityDNAParticleList,areaDNAParticleList,
                        numberOfRNAParticleList,meanIntensityRNAParticleList,areaRNAParticleList)
{              
  labels = ['FileName','Cell Number List',
                    '# DNA Particle','DNA Particle Intensity Mean','DNA Particle Area Mean',
                    '# RNA Particle','RNA Particle Intensity Mean','RNA Particle Area Mean'
                    ];

  rt.reset()
  for(j=0;j<fileNameList.size();j++)
  {
    label=0;
    //Values
    rt.incrementCounter()
    rt.addValue(labels[label++],fileNameList[j])              //
    rt.addValue(labels[label++],cellNumberList[j])        //
    rt.addValue(labels[label++],numberOfDNAParticleList[j])
    rt.addValue(labels[label++],meanIntensityDNAParticleList[j])
    rt.addValue(labels[label++],areaDNAParticleList[j])
    rt.addValue(labels[label++],numberOfRNAParticleList[j])
    rt.addValue(labels[label++],meanIntensityRNAParticleList[j])
    rt.addValue(labels[label++],areaRNAParticleList[j])    
  }
  return rt
}


def analyzeCellParticle(impCell,impCellMask,roiCell,minParticleSize,maxParticleSize,numberOfParticleList,meanParticleIntensityList,areaParticleList,rm,thresholdMethod)
{
  rm.reset()
  impCellMask.setRoi(roiCell)
  IJ.setBackgroundColor(255, 255, 255);
  IJ.run(impCellMask, "Clear Outside", "");
  particleAnalysisOptions = ParticleAnalyzer.ADD_TO_MANAGER
  p  = new ParticleAnalyzer(particleAnalysisOptions, 0, rt, minParticleSize, maxParticleSize, 0.0, 1.0)
  p.setRoiManager(rm)
  p.analyze(impCellMask)
  def particleRois=rm.getRoisAsArray()
  def meanParticlesIntensity=0
  def areaParticles=0
  numberOfParticleList.add(particleRois.size())
  if(particleRois.size()>0)
  {
    for(int p=0;p<particleRois.size();p++)
    {
      impCell.setRoi(particleRois[p])
      stats=impCell.getStatistics(Measurements.MEAN+Measurements.MEDIAN+Measurements.AREA+Measurements.STD_DEV )
      meanParticlesIntensity+=stats.@mean
      areaParticles+=stats.@area
    }
    
    meanParticlesIntensity=meanParticlesIntensity/particleRois.size()
    meanParticleIntensityList.add(meanParticlesIntensity)
    areaParticles=areaParticles/particleRois.size()
    areaParticleList.add(areaParticles) 
    //impCellMask.show()
    //lll
  }
  else
  {
    meanParticleIntensityList.add(-1)
    areaParticleList.add(-1)
  }
  rm.reset()
  return[numberOfParticleList,meanParticleIntensityList,areaParticleList]
}


As the RNA counts seem to be correct in the second image I assume the problem is with the image and not with the script?

I would appreciate any suggestion!
Best,
Julie


#2

The Yen (and many other) autothresholding methods compute the threshold value from the image histogram.
I think there is no way to guarantee that you can detect the same object across different images which have different histograms, regardless of whether the images look similar or not.
You probably need a different segmentation approach.


#3

I agree with Gabriel. Intensity based segmentation can be tricky if there is variability over the samples. Either from imaging, staining or the biology.

It seems that you want to segment the spots within each nucleus.
So maybe rather go with a spot or maxima detection method:
https://imagej.nih.gov/ij/docs/guide/146-29.html#Flo:Maxima-outputs

Also other processing steps such as Difference of Gaussian (DoG) or Laplacian of Gaussian (LoG) are really good for spot detection and make the detection more robust to noise:
http://bigwww.epfl.ch/sage/soft/LoG3D/


#4

Depending on what you want to quantify you can add more complex segmentation also. Such as Seeded watershed. Look into the morphoLibJ plugin: https://imagej.net/MorphoLibJ
Your images kind of remind me anyways of the examples given in the wiki there. So could be something worth to investigate.


#5

*Hi Iā€™m new to imageJ and I could not be able to create a .obj 3D file without losing some data from the .raw file stack.
*Basically the .raw file contains number of images arranged in stack.
*These images are black and white images only. Each and every pixel in this image is very important.
*After importing the .raw file into the imagej it can be saved as a .obj generally. But while saving this .obj even after with the low threshold and re-sampling factor as 1, some of the pixel datas are missing in the 3d model.
*Can anyone help me with this issue please?