Problems with a watershed-based touching particle separator

[John Minter]

I have been trying to develop a plug-in to perform a watershed
algorithm on cryo-TEM images of colloidal particles where some touch.
I want to generate the Voronoi polyhedra and write them as while lines
into the input image (the image I am attaching is a part of a larger
image that has been shading corrected and median filtered to reduce
noise.) The idea is then analyze the separated image and analyze the
blobs. I have two problems that I have been unable to solve:

1. This one blows my mind. I have build the attached plug-in on three
computers - all with Sun Java 1.6.0.29. The plugin writes the white
lines into the test image on two systems: my MacBook Pro and my work
Linux box with an AMD X86_64 processor running Xubuntu x86_64 Oneiric
Ocelot (11.10). I also have a home Linux box with an Intel X86_64
processor Xubuntu x86_64 Oneiric Ocelot set up identically with my
work box and this one does NOT write the white separator lines. This
example is just 8 bits/px (some day I would like to work w 16 bit/px),
so I don't see where byte order would be important in this case. I
hope someone has an idea what is wrong here. I am stumped...

2. Issues with the watershed. I don't seem to be able to tune the
watershed to actually separate the touching particles. it gets the
easy non-touching ones, but can't separate the touching ones. I would
also like to thicken the separator lines to 2 pixels, but have not
figured out how to do so.

the source ("Do_Separator.java") is pasted below. I'm attaching two
files: the input
image ("orig.jpg") and a processed version that worked ("proc.jpg"). I
would appreciate any and all tips.

Best regards,
John Minter

The source:


/**
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * John Minter
 * <[hidden email]>
 *
 * Date: 2012/02/28
 *
 * write lines from Voronoi to original
 *
 */

import ij.*;
import ij.process.*;
import ij.measure.*;
// import ij.gui.*;
import ij.ImagePlus;
// import ij.plugin.Thresholder;
import ij.plugin.filter.*;
import ij.plugin.filter.PlugInFilter;
// import ij.plugin.filter.BackgroundSubtracter;
import ij.process.ImageProcessor;
// import ij.WindowManager.*;

// import java.awt.*;

public class Do_Separator implements PlugInFilter {
        // our ImagePlus object
        private ImagePlus m_objImpInp; // for the input image
        private ImagePlus m_objImpCop; // for the smoothed image
        // A debug info variable. Set as appropriate
        private boolean m_bVerbose = false;

        public int setup(String arg, ImagePlus imp) {
    m_objImpInp = imp;
                return DOES_8G;
        }

        public void run(ImageProcessor objIpIn) {

    int iW = objIpIn.getWidth();
    int iH = objIpIn.getHeight();
    String strTitle = "sep-" + m_objImpInp.getTitle();
    Calibration objCalib = m_objImpInp.getCalibration();
    // first save a copy of the original
    m_objImpCop = m_objImpInp.duplicate();
    m_objImpCop.setCalibration(objCalib);
    m_objImpCop.show();

    GaussianBlur fltGaussianBlur = new GaussianBlur();
                fltGaussianBlur.blur(objIpIn, 5.0);
    m_objImpInp.updateAndDraw();
    if(m_bVerbose) IJ.showMessage("blur");


    int[] Hist = objIpIn.getHistogram();
    AutoThresholder objAthr = new AutoThresholder();
    int iThr = objAthr.getThreshold("Otsu", Hist);

    int iVal;
    byte byWhite = (byte) 255;
    byte byBlack = (byte)   0;
    byte[] pix_arr = (byte[]) objIpIn.getPixels();
    for(int y=0; y<iH; y++){
      for(int x=0; x<iW; x++){
        iVal = (int)  (pix_arr[y*iW+x] & 0xff );
        if(iVal > iThr){
          pix_arr[y*iW+x] = byWhite;
        } else {
          pix_arr[y*iW+x] = byBlack;
        }
      }
    }
    objIpIn.setPixels(pix_arr);
    m_objImpInp.updateAndDraw();


    objIpIn.setAutoThreshold(ImageProcessor.ISODATA2,
ImageProcessor.NO_LUT_UPDATE);
                double minThreshold = objIpIn.getMinThreshold();
                double maxThreshold = objIpIn.getMaxThreshold();
    int[] lut = new int[256];
    for (int j=0; j<256; j++) {
      if (j>=minThreshold && j<=maxThreshold)
        lut[j] = byWhite;
      else
        lut[j] = byBlack;
    }
    objIpIn.applyTable(lut);
    objIpIn.setThreshold(255, 255, ImageProcessor.NO_LUT_UPDATE);
    m_objImpInp.updateAndDraw();


    EDM filtEDM = new EDM();
    filtEDM.setup("voronoi", m_objImpInp);
    filtEDM.run(objIpIn);
    filtEDM.setup("final", m_objImpInp);
    m_objImpInp.updateAndDraw();

    if(m_bVerbose) IJ.showMessage("voronoi");

    byte[] pix_work  = (byte[]) objIpIn.getPixels();
    byte[] pix_copy = (byte[]) m_objImpCop.getProcessor().getPixels();
    int iVorVal;
    for(int y=0; y<iH; y++){
      for(int x=0; x<iW; x++){
        iVorVal = (int)  (pix_work[y*iW+x] & 0xff );
        if(iVorVal > 0){
          // this is a separator line
          // so set the value in the copy
          // of the input image to white
          pix_copy[y*iW+x] = byWhite;
        }
      }
    }
    // so now stuff the corrected pixels
    // from the copy into the main image
    objIpIn.setPixels(pix_copy);
    objIpIn.resetMinAndMax();
    m_objImpInp.updateAndDraw();
    for (int j=0; j<256; j++) {
        lut[j] = (byte) j;
    }
    objIpIn.applyTable(lut);
    m_objImpInp.setTitle(strTitle);
    m_objImpInp.updateAndDraw();

    if(!m_bVerbose) m_objImpCop.close();
        }
}

[Michael Schmid]

Hi John,

without looking at your plugin, I think that this can be done by ImageJ
without plugins:
- threshold the image (default settings are fine) and 'Apply'
- Gaussian Blur filter of the result (radius about 4 pixels)
- threshold again
- Process>Binary>Watershed
- Process>Binary>Voronoi and threshold from 1 to maximum, if desired

Michael
_______________________________________________________________________

On Fri, March 2, 2012 04:55, John Minter wrote:

> I have been trying to develop a plug-in to perform a watershed
> algorithm on cryo-TEM images of colloidal particles where some touch.
> I want to generate the Voronoi polyhedra and write them as while lines
> into the input image (the image I am attaching is a part of a larger
> image that has been shading corrected and median filtered to reduce
> noise.) The idea is then analyze the separated image and analyze the
> blobs. I have two problems that I have been unable to solve:
>
> 1. This one blows my mind. I have build the attached plug-in on three
> computers - all with Sun Java 1.6.0.29. The plugin writes the white
> lines into the test image on two systems: my MacBook Pro and my work
> Linux box with an AMD X86_64 processor running Xubuntu x86_64 Oneiric
> Ocelot (11.10). I also have a home Linux box with an Intel X86_64
> processor Xubuntu x86_64 Oneiric Ocelot set up identically with my
> work box and this one does NOT write the white separator lines. This
> example is just 8 bits/px (some day I would like to work w 16 bit/px),
> so I don't see where byte order would be important in this case. I
> hope someone has an idea what is wrong here. I am stumped...
>
> 2. Issues with the watershed. I don't seem to be able to tune the
> watershed to actually separate the touching particles. it gets the
> easy non-touching ones, but can't separate the touching ones. I would
> also like to thicken the separator lines to 2 pixels, but have not
> figured out how to do so.
>
> the source ("Do_Separator.java") is pasted below. I'm attaching two
> files: the input
> image ("orig.jpg") and a processed version that worked ("proc.jpg"). I
> would appreciate any and all tips.
>
> Best regards,
> John Minter
>

[John Minter-3]

You are correct. I prototyped the plug-in that way. I was trying to
improve my skills by developing the plug-in. I thought it might be
faster to have this as a java file for later processing . I still
don't why the same code runs differently on the amd x_86_64 processor
and on the intel x86_64 procressor...

Regards,
John.

On Fri, Mar 2, 2012 at 2:58 AM, Michael Schmid <[hidden email]> wrote:

> Hi John,
>
> without looking at your plugin, I think that this can be done by ImageJ
> without plugins:
> - threshold the image (default settings are fine) and 'Apply'
> - Gaussian Blur filter of the result (radius about 4 pixels)
> - threshold again
> - Process>Binary>Watershed
> - Process>Binary>Voronoi and threshold from 1 to maximum, if desired
>
> Michael
>

[Michael Schmid]

Hi John,

usually, there is hardly any speed gain when converting a macro to a
plugin unless the macro contains loops with many repetitions, and the loop
contains rather elementary operations such as single-pixel calculations.

Anyhow, I have had a short look at your plugin and could not find an
obvious fault, but it is rather complicated compared to the simplicity of
the task.  I think that the use of a second ImagePlus object is
unnecessary and might be the reason for problems.  Please also note that
the GaussianBlur.blur method is deprecated (kept for compatibility with
very old versions of ImageJ only); its 'radius' variable does not
correspond to the 'sigma' of the Gaussian Blur command.

I have written a simple plugin that does essentially the same job, based
on the 'segmented particles' feature of 'Find Maxima'.  I'll paste it to
the end.

Michael
________________________________________________________________________


On Sat, March 3, 2012 01:44, John Minter wrote:

> You are correct. I prototyped the plug-in that way. I was trying to
> improve my skills by developing the plug-in. I thought it might be
> faster to have this as a java file for later processing . I still
> don't why the same code runs differently on the amd x_86_64 processor
> and on the intel x86_64 procressor...
>
> Regards,
> John.
>
> On Fri, Mar 2, 2012 at 2:58 AM, Michael Schmid <[hidden email]>
> wrote:
>> Hi John,
>>
>> without looking at your plugin, I think that this can be done by ImageJ
>> without plugins:
>> - threshold the image (default settings are fine) and 'Apply'
>> - Gaussian Blur filter of the result (radius about 4 pixels)
>> - threshold again
>> - Process>Binary>Watershed
>> - Process>Binary>Voronoi and threshold from 1 to maximum, if desired
>>
>> Michael
>>
>

________________________________________________________________________

import ij.*;
import ij.process.*;
import ij.ImagePlus;
import ij.plugin.filter.*;

public class Do_Separator1 implements PlugInFilter {

    public int setup(String arg, ImagePlus imp) {
        //Ask to process stacks.
        //SNAPSHOT: We always need a snapshot, also if not needed for
        //undo (e.g. in case of a stack).
        return IJ.setupDialog(imp, DOES_8G | DOES_16 | SNAPSHOT);
    }

    public void run(ImageProcessor ip) {
        int width = ip.getWidth();
        int height = ip.getHeight();
        //don't use old deprecated 'GaussianBlur.blur' method
        new GaussianBlur().blurGaussian(ip, 5.0, 5.0, 0.02);
        //for images with foreground = low values:
        ip.invert();
        //segmentation lines get value=0 in 'byteLines'
        ByteProcessor byteLines = new MaximumFinder().findMaxima(ip, 0.0,
                ImageProcessor.NO_THRESHOLD,
                MaximumFinder.SEGMENTED, /*excludeOnEdges=*/ false, false);
        //undo Gaussian Blur&invert
        //(we have a snapshot due to the SNAPSHOT flag)
        ip.reset();
        int lineValueInt = 128;
        if (ip instanceof ShortProcessor)
            lineValueInt = 32768;
        //or whatever value you want for the lines
        byte[] bPixels = (byte[])byteLines.getPixels();
        for (int y=0, p=0; y<height; y++)   //'p' is pointer in pixels array
            for (int x=0; x<width; x++, p++)
                if (bPixels[p]==0)
                    ip.set(x,y,lineValueInt);
    }

}

[John Minter-3]

I'd like to publicly thank Michael Schmid for rewriting my problematic
plug-in. His version works on all platforms, solving the issues I had.
I did make one minor change - I find the maximum intensity in the
input image and write the Voronoi polyhedra lines with this value.
This makes threshold adjustment and particle analysis easier.

Regards,
John Minter

On Sun, Mar 4, 2012 at 2:39 PM, Michael Schmid <[hidden email]> wrote:

> Hi John,
>
> usually, there is hardly any speed gain when converting a macro to a
> plugin unless the macro contains loops with many repetitions, and the loop
> contains rather elementary operations such as single-pixel calculations.
>
> Anyhow, I have had a short look at your plugin and could not find an
> obvious fault, but it is rather complicated compared to the simplicity of
> the task.  I think that the use of a second ImagePlus object is
> unnecessary and might be the reason for problems.  Please also note that
> the GaussianBlur.blur method is deprecated (kept for compatibility with
> very old versions of ImageJ only); its 'radius' variable does not
> correspond to the 'sigma' of the Gaussian Blur command.
>
> I have written a simple plugin that does essentially the same job, based
> on the 'segmented particles' feature of 'Find Maxima'.  I'll paste it to
> the end.
>
> Michael
> ________________________________________________________________________
>
>
> On Sat, March 3, 2012 01:44, John Minter wrote:
>> You are correct. I prototyped the plug-in that way. I was trying to
>> improve my skills by developing the plug-in. I thought it might be
>> faster to have this as a java file for later processing . I still
>> don't why the same code runs differently on the amd x_86_64 processor
>> and on the intel x86_64 procressor...
>>
>> Regards,
>> John.
>>
>> On Fri, Mar 2, 2012 at 2:58 AM, Michael Schmid <[hidden email]>
>> wrote:
>>> Hi John,
>>>
>>> without looking at your plugin, I think that this can be done by ImageJ
>>> without plugins:
>>> - threshold the image (default settings are fine) and 'Apply'
>>> - Gaussian Blur filter of the result (radius about 4 pixels)
>>> - threshold again
>>> - Process>Binary>Watershed
>>> - Process>Binary>Voronoi and threshold from 1 to maximum, if desired
>>>
>>> Michael
>>>
>>
>
> ________________________________________________________________________
>
> import ij.*;
> import ij.process.*;
> import ij.ImagePlus;
> import ij.plugin.filter.*;
>
> public class Do_Separator1 implements PlugInFilter {
>
>    public int setup(String arg, ImagePlus imp) {
>        //Ask to process stacks.
>        //SNAPSHOT: We always need a snapshot, also if not needed for
>        //undo (e.g. in case of a stack).
>        return IJ.setupDialog(imp, DOES_8G | DOES_16 | SNAPSHOT);
>    }
>
>    public void run(ImageProcessor ip) {
>        int width = ip.getWidth();
>        int height = ip.getHeight();
>        //don't use old deprecated 'GaussianBlur.blur' method
>        new GaussianBlur().blurGaussian(ip, 5.0, 5.0, 0.02);
>        //for images with foreground = low values:
>        ip.invert();
>        //segmentation lines get value=0 in 'byteLines'
>        ByteProcessor byteLines = new MaximumFinder().findMaxima(ip, 0.0,
>                ImageProcessor.NO_THRESHOLD,
>                MaximumFinder.SEGMENTED, /*excludeOnEdges=*/ false, false);
>        //undo Gaussian Blur&invert
>        //(we have a snapshot due to the SNAPSHOT flag)
>        ip.reset();
>        int lineValueInt = 128;
>        if (ip instanceof ShortProcessor)
>            lineValueInt = 32768;
>        //or whatever value you want for the lines
>        byte[] bPixels = (byte[])byteLines.getPixels();
>        for (int y=0, p=0; y<height; y++)   //'p' is pointer in pixels array
>            for (int x=0; x<width; x++, p++)
>                if (bPixels[p]==0)
>                    ip.set(x,y,lineValueInt);
>    }
>
> }