[Insight-users] Bug in watershed filter

[Richard Beare]

I'll attempt to check this out over the next couple of days.

I suspect that you have run into issues with discrete nature of
images, rather than a bug as such. Suppose that you are applying a
watershed to a flat control surface. If the seed points are an odd
number of voxels apart, and you chose the algorithm marking the
watershed line, then the split will appear nicely down the middle. If
you don't chose to mark the line then there will be a probably random
choice of which voxels along the midpoint get assigned to which
region. Conversely, if the seeds are an even number of voxels apart
and you chose to mark the watershed line, then the line should lie
between the pixels, but it can't. Hence it is likely to be jaggedly
jumping between the two sides of the plateau. Similar things can
happen with distance map control functions, as adjacent voxels can be
equidistant from seeds, producing a plateau.


On Thu, Feb 14, 2013 at 9:13 AM, Dženan Zukić <dzenanz at gmail.com> wrote:
> I think I have discovered a bug in WatershedImageFilter. Below is MWE which
> reproduces the problem. Input, output and expected image attached.
>
> The problem is that watershed lines are wrongly placed and rough (see
> attached error highlights). I first discovered it while working with 3D
> images, where the problem is more pronounced. I can also supply 3D float
> images (1.3MB). MorphologicalWatershedFromMarkersImageFilter suffers from
> the same problem.
>
> #include <itkImageFileReader.h>
> #include <itkImageFileWriter.h>
> #include <itkWatershedImageFilter.h>
> #include <itkCastImageFilter.h>
>
> typedef itk::Image<unsigned char, 2> VisualizingImageType;
>
> void main()
> {
>     typedef itk::ImageFileReader<VisualizingImageType> rType;
>     rType::Pointer reader=rType::New();
>     reader->SetFileName("image2D.png");
>     reader->Update();
>     VisualizingImageType::Pointer dm=reader->GetOutput();
>     typedef itk::WatershedImageFilter<VisualizingImageType> wsType;
>     wsType::Pointer ws2=wsType::New();
>     ws2->SetThreshold(0);
>     ws2->SetLevel(0);
>     ws2->SetInput(dm);
>     ws2->Update();
>     typedef itk::CastImageFilter<wsType::OutputImageType,
> VisualizingImageType> castType;
>     castType::Pointer caster=castType::New();
>     caster->SetInput(ws2->GetOutput());
>     caster->Update();
>     typedef itk::ImageFileWriter<VisualizingImageType> wType;
>     wType::Pointer writer=wType::New();
>     writer->SetFileName("imageWSint.png");
>     writer->SetInput(caster->GetOutput());
>     writer->Update();
> }