[Insight-users] Still struggling with implementing grayscale morphological operators

[Neal R. Harvey]

I am still struggling to implement grayscale morphological operators.
The basic functions have been written:
e.g. 
http://www.itk.org/Doxygen310/html/classitk_1_1GrayscaleFunctionDilateImageFilter.html
Unfortunately there are no nice papers or tutorials, or even test code 
to provide insight as to exactly
how to go about using them.

It appears that they use a neighborhood kernel and I haven't been able 
to figure out how to create a
grayscale kernel. I have gone through the book - the part that describes 
neighborhood iterators, but
so far that hasn't been a great deal of help. The examples they provide 
don't exactly show how to
assign "grayscale" values to your kernel neighborhood iterator. The 
examples are for sobel and gaussian
kernels, where there is code in the background that calculates the 
values for you. You never see the details
of how it is used. The other examples they show in the book are for 
binary morphology, but that doesn't
require graylevel values kernels, so also isn't much use.

Basically, I think I need to be able to read in an input image that is 
to be processed and two images that
define the grayscale structuring element (SE), and output a processed 
image. One of the SE images
defines the region of support of the SE (i.e. anywhere where the values 
are greater than 0 are in the
SE's region of support). The other SE image defines the grayscale values 
in the SE. I then have to use
these two SE images to create a suitable kernel.

Below is the code that I have written so far, in attempting to write my 
own Grayscale Morphological
Operator based on the examples provided in the book. As you can see, 
it's not working and I am not
sure how to solve that problem.

If anyone has any suggestions or ideas of how to get this working, or 
where to look to get information
that could lead in that direction, it would be very much appreciated.

Kindest regards

Harve

#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#ifdef __BORLANDC__
#define ITK_LEAN_AND_MEAN
#endif

#include <iostream>
#include <string>
#include <sstream>

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkTIFFImageIO.h"
#include "itkNumericTraits.h"
#include "itkMinimumMaximumImageFilter.h"
#include "itkConstNeighborhoodIterator.h"
#include "itkConstShapedNeighborhoodIterator.h"
#include "itkNeighborhoodAlgorithm.h"
#include "itkImageRegionIterator.h"

using namespace std;

string itos(int i)    // convert int to string
{
  stringstream s;
  s << i;
  return s.str();
}

int main( int argc, char * argv[] )
{

  // Verify the number of parameters in the command line
  if( argc < 5 )
    {
      std::cerr << "Usage: " << argv[0];
      std::cerr << " inputImage inputSEValuesImage inputSEExtentImage 
outputImage ";
      return EXIT_FAILURE;
    }


  //  Then, as usual, we select the pixel types and the image
  //  dimension. Remember, if the file format represents pixels with a
  //  particular type, C-style casting will be performed to convert the 
data.

  typedef unsigned char       InputPixelType;
  typedef unsigned int        InputSEValuesPixelType;
  typedef int                 InputSEExtentPixelType;
  typedef unsigned char       OutputPixelType;
  typedef float               FloatPixelType;
  const   unsigned int        Dimension = 2;

  typedef itk::Image< InputPixelType,  Dimension >         InputImageType;
  typedef itk::Image< InputSEValuesPixelType, Dimension >  
InputSEValuesImageType;
  typedef itk::Image< InputSEExtentPixelType, Dimension >  
InputSEExtentImageType;
  typedef itk::Image< OutputPixelType, Dimension >         OutputImageType;
  typedef itk::Image< FloatPixelType, Dimension >          FloatImageType;

  typedef itk::ConstShapedNeighborhoodIterator< InputImageType > 
ShapedNeighborhoodIteratorType;

  typedef itk::ImageRegionIterator< InputImageType > IteratorType;

  //  Below we instantiate the MinimumMaximumImageCalculator class

  typedef itk::MinimumMaximumImageFilter< InputImageType > 
MinimumMaximumInputImageFilterType;

  typedef itk::MinimumMaximumImageFilter< InputSEValuesImageType > 
MinimumMaximumSEValuesImageFilterType;

  typedef itk::MinimumMaximumImageFilter< InputSEExtentImageType > 
MinimumMaximumSEExtentImageFilterType;

  //  We can now instantiate the reader and writer. These two classes are
  //  parameterized over the image type. We instantiate the
  //  TIFFImageIO class as well. Note that the ImageIO objects are
  //  not templated.

  typedef itk::ImageFileReader< InputImageType  >          
InputImageReaderType;

  typedef itk::ImageFileReader< InputSEValuesImageType  >   
InputSEValuesImageReaderType;

  typedef itk::ImageFileReader< InputSEExtentImageType  >  
InputSEExtentImageReaderType;

  typedef itk::ImageFileWriter< OutputImageType >          
OutputImageWriterType;

  typedef itk::TIFFImageIO                                 ImageIOType;

  //  Then, we create one object of each type using the New() method and
  //  assigning the result to a SmartPointer.

  InputImageReaderType::Pointer InputImageReader = 
InputImageReaderType::New();

  InputSEValuesImageReaderType::Pointer InputSEValuesImageReader = 
InputSEValuesImageReaderType::New();

  InputSEExtentImageReaderType::Pointer InputSEExtentImageReader = 
InputSEExtentImageReaderType::New();

  //  A MinMaxImageFilter object is constructed

  MinimumMaximumInputImageFilterType::Pointer 
MinimumMaximumInputImageFilter = MinimumMaximumInputImageFilterType::New();

  MinimumMaximumSEValuesImageFilterType::Pointer 
MinimumMaximumSEValuesImageFilter = 
MinimumMaximumSEValuesImageFilterType::New();

  MinimumMaximumSEExtentImageFilterType::Pointer 
MinimumMaximumSEExtentImageFilter = 
MinimumMaximumSEExtentImageFilterType::New();

  OutputImageWriterType::Pointer OutputImageWriter = 
OutputImageWriterType::New();
  ImageIOType::Pointer tiffIO = ImageIOType::New();

  //
  // Here we recover the file names from the command line arguments
  //
  const char * inputImageFileName         = argv[1];
  const char * inputSEValuesImageFileName = argv[2];
  const char * inputSEExtentImageFileName = argv[3];
  const char * outputImageFileName        = argv[4];

  //  The name of the file to be read or written is passed with the
  //  SetFileName() method.

  InputImageReader->SetFileName( inputImageFileName );

  InputImageReader->Update();

  InputSEValuesImageReader->SetFileName( inputSEValuesImageFileName );

  InputSEValuesImageReader->Update();

  InputSEExtentImageReader->SetFileName( inputSEExtentImageFileName );

  InputSEExtentImageReader->Update();

  MinimumMaximumInputImageFilter->SetInput( InputImageReader->GetOutput() );

  MinimumMaximumInputImageFilter->Update();

  InputPixelType imageMinVal = MinimumMaximumInputImageFilter->GetMinimum();
  InputPixelType imageMaxVal = MinimumMaximumInputImageFilter->GetMaximum();

  std::cerr << "Maximum Value from Input Image = " << 
static_cast<int>(imageMaxVal) << std::endl;
  std::cerr << "Minimum Value from Input Image = " << 
static_cast<int>(imageMinVal) << std::endl;

  MinimumMaximumSEValuesImageFilter->SetInput( 
InputSEValuesImageReader->GetOutput() );

  MinimumMaximumSEValuesImageFilter->Update();

  InputSEValuesPixelType SEValuesMinVal = 
MinimumMaximumSEValuesImageFilter->GetMinimum();
  InputSEValuesPixelType SEValuesMaxVal = 
MinimumMaximumSEValuesImageFilter->GetMaximum();

  std::cerr << "Maximum Value from Input SE Values Image = " << 
static_cast<int>(SEValuesMaxVal) << std::endl;
  std::cerr << "Minimum Value from Input SE Values Image = " << 
static_cast<int>(SEValuesMinVal) << std::endl; 

  MinimumMaximumSEExtentImageFilter->SetInput( 
InputSEExtentImageReader->GetOutput() );

  MinimumMaximumSEExtentImageFilter->Update();

  InputSEExtentPixelType SEExtentMinVal = 
MinimumMaximumSEExtentImageFilter->GetMinimum();
  InputSEExtentPixelType SEExtentMaxVal = 
MinimumMaximumSEExtentImageFilter->GetMaximum();

  std::cerr << "Maximum Value from Input SE Extent Image = " << 
static_cast<int>(SEExtentMaxVal) << std::endl;
  std::cerr << "Minimum Value from Input SE Extent Image = " << 
static_cast<int>(SEExtentMinVal) << std::endl; 

  long int inputImageCols = 
InputImageReader->GetOutput()->GetLargestPossibleRegion().GetSize()[0];

  long int inputImageRows = 
InputImageReader->GetOutput()->GetLargestPossibleRegion().GetSize()[0];

  std::cout << "Number of Columns for Input Image = " << inputImageCols 
<< "; Number of Rows for Input Image = " << inputImageRows << std::endl;

  long int inputSEValuesImageCols = 
InputSEValuesImageReader->GetOutput()->GetLargestPossibleRegion().GetSize()[0];

  long int inputSEValuesImageRows = 
InputSEValuesImageReader->GetOutput()->GetLargestPossibleRegion().GetSize()[0];

  std::cout << "Number of Columns for Input SE Values Image = " << 
inputSEValuesImageCols << "; Number of Rows for Input SE Values Image = 
" << inputSEValuesImageRows << std::endl;

  long int inputSEExtentImageCols = 
InputSEExtentImageReader->GetOutput()->GetLargestPossibleRegion().GetSize()[0];

  long int inputSEExtentImageRows = 
InputSEExtentImageReader->GetOutput()->GetLargestPossibleRegion().GetSize()[0];

  std::cout << "Number of Columns for Input SE Extent Image = " << 
inputSEExtentImageCols << "; Number of Rows for Input SE Extent Image = 
" << inputSEExtentImageRows << std::endl;

  // Check that the dimensions of the inputSEValuesImage and 
inputSEExtentImage are the same

  if ((inputSEValuesImageCols != inputSEExtentImageCols) || 
(inputSEValuesImageRows != inputSEExtentImageRows)) {

    std::cout << "inputSEValuesImage and inputSEExtentImage do not have 
the same dimensions - EXITING!" << std::endl;
    return EXIT_FAILURE;

  }

  // OK, so now we have read in the input image, an image containing the 
greyscale values for the Structuring Element
  // and an image containing 1s and 0s defining the region of support of 
the S.E.

  // First - calculate the SE radius from the size of the Input SE image

  // Which is the larger of the two dimensions

  long int SE_size = ((inputSEExtentImageCols > inputSEExtentImageRows) 
? inputSEExtentImageCols : inputSEExtentImageRows);

  std::cout << "SE Size = " << SE_size << std::endl;

  // Is this dimension an odd or even number

  bool even = ((SE_size % 2) == 0) ? true : false;

  string response = (even == true) ? "even" : "odd";

  std::cout << "SE dimension is " << response << std::endl;

  long int SE_radius = (SE_size / static_cast<long int>(2));

  std::cout << "SE radius is " << SE_radius << std::endl;

  // Now we need to do stuff for the iterator

  ShapedNeighborhoodIteratorType::RadiusType radius;
  radius.Fill(SE_radius);

  OutputImageType::Pointer OutputImage = OutputImageType::New();
  OutputImage->SetRegions( 
InputImageReader->GetOutput()->GetRequestedRegion());
  OutputImage->Allocate();

  typedef itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< 
InputImageType > FaceCalculatorType;

  FaceCalculatorType faceCalculator;
  FaceCalculatorType::FaceListType faceList;
  FaceCalculatorType::FaceListType::iterator fit;
 
  faceList = faceCalculator( InputImageReader->GetOutput(),
                             OutputImage->GetRequestedRegion(),
                             radius );

  IteratorType out;

  InputSEExtentImageType::IndexType pixelIndex;

  InputSEExtentPixelType::IndexType SEExtentPixelValue;

  InputPixelType InputPixelValue;

  // Looping through each of the regions determined by faceCalculator
  for ( fit=faceList.begin(); fit != faceList.end(); fit++) {

    ShapedNeighborhoodIteratorType InputIterator( radius, 
InputImageReader->GetOutput(), *fit );

    out = IteratorType ( outputImage, *fit );

    // Defining what the offsets are for the InputIterator

    // However, how does one set the "grayscale" values of the iterator 
kernel?

    for (register int j = -SE_radius; j <= SE_radius; j++) {
      for (register int i = -SE_radius; i <= SE_radius; i++) {

    ShapedNeighborhoodIteratorType::OffsetType off;

    pixelIndex[0] = (i + SE_radius); // Set the pixelIndex x/column position
    pixelIndex[1] = (j + SE_radius); // Set the pixelIndex y/row position

    SEExtentPixelValue = 
InputSEExtentImageReader->GetOutput()->GetPixel( pixelIndex );

    if (SEExtentPixelValue == static_cast<InputSEExtentPixelType>(1)) {

      off[0] = static_cast<int>(i);
      off[1] = static_cast<int>(j);

      InputIterator.ActivateOffset(off);
    }

      }
    }

    // Implements GrayScale Erosion

    for (InputIterator.GoToBegin(), out.GoToBegin(); 
!InputIterator.IsAtEnd(); ++InputIterator, ++out) {

      ShapedNeighborhoodIteratorType::ConstIterator ci;

      OutputPixelType minValue = static_cast<OutputPixelType>(imageMaxVal);

      for (ci = InputIterator.Begin(); ci != InputIterator.End(); ci++) {
   
    InputPixelValue = ci.Get();

        // Here we want to get the "grayscale" value of the 
corresponding SE point - but how?

    // Then we want to subtract the SE grayscale value from the 
inputPixelValue -
    // but how, if we don't know how to set or retrieve the SE kernel 
values?

    // Then see if this is less than the minValue - if it is set it to 
the minValue

      }

      // Then set the output value to be this minValue

      out.Set(minValue);
    }

  }