[Insight-users] Driver for class MRFImageFilter

Fri, 12 Dec 2003 10:54:33 -0500

Hello,

I've recently been able to spend some time on the MRF code contained
in ITK, specifically with the MRFImageFilter class.  I've come to some
conclusions and have included some code here so that others might 
save some time.

1.   Currently, the only driver for MRF segmentation is GibbsPriorFilter1.cxx
which is a driver for the RGBGibbsPriorFilter class.  This class is derived
from MRFImageFilter, so it seems to be the logical place to start.   However,
there are numerous problems with both the "Gibbs" driver and class code,
and in my opinion the entirely wrong place to start.    In fact, in my opinion
the "Gibbs" code ought to be removed entirely until it can be rewritten.
(I recognize that ITK is a work in progress, but I'll outline my reasoning below.)
Wanting to be part of the solution and not the problem, I've provided a 
driver for the MRFImageFilter class at the end of this message which I 
think is a good place to start for testing MRF segmentation.

2.   This message is also a follow-on to an October 23 message that I 
sent to the Insight
users list concerning EstimateGaussianModelPrameters contained 
in itkImageGaussModelEstimator.txx.   There were no comments, but perhaps
with a driver for the MRFImageFilter the number of users will increase.  In
short, I modified a few lines of code so that negative pixels are ignored
and n  classes are labeled 0 through n-1.   The original code had a behavior
such that class zero was always computed with no pixels and negative
pixels were silently converted to large positive classes.  In another
part of the documentation it is recommend that ignored pixels should
be the maximum pixel value, but this was not implemented as far a I
can see.   The
problem here is with C array, which like things to start at zero.  I didn't
want to modify everything, so I just assumed on the input that 0 implied
ignore, and thus subtracted 1 from a converted short value before
sending the image to the estimator. The driver
below is written to work with the changes to EstimateGaussianModelPrameters
I outlined in the October 23rd message, but can easily be modified for other schemes.

3.   The reason I believe that the "Gibbs" code should be removed (or
at least deprecated) is three-fold:  the algorithm is incorrect, the code
has numeric and memory inconsistencies, and it not a proper ITK 
derived object.

        a.   My first big clue ;-) about the algorithm being incorrect was the following line
in itkRGBGibbsPriorFilter.

#define RGBGibbsPriorFilterNeedsDebugging 1

This removes about a third of the active code, and reduces it more-or-less to MRFImageFilter.   (There
are still some differences, to  be sure.)   I haven't tried undoing this yet, as the other two
problems still need to be addressed.

        b.    There are a number of numeric and memory errors. In one
section the input type data is converted to labeled type and
then converted back to input type.  (GreyScalarBoundary).    If the input data were floats between 0 and
1 the algorithm would be completely wrong.  Similarly the memory problem stems from
generating an array m_LabelStatus via new[] which is never deleted []. 
Finally, the labeling method assumes
that there exists only one class for segmentation (index 1) and that no other classes exist.
This assumption is evident in RegionEraser(), were there is the line

labelledImageIt.Set(1-label);

Clearly, label is assumed to be either 0 or 1.

        c.  Finally, the code itself is not written in ITK object oriented style.   Here are just
a few examples.   First, in itkGibbsPriorFilter.h, the following lines appear in the class
definition of RGBGibbsPriorFilter:

InputImageConstPointer m_InputImage; /** the input */
TrainingImageType m_TrainingImage; /** image to train the filter. */
LabelledImageType m_LabelledImage; /** output */
unsigned int m_NumberOfClasses; /** the number of class need to be classified. */ 
unsigned int m_MaximumNumberOfIterations; /** number of the iteration. */

The first three are entirely unnecessary:  this is a ProcessObject, and as such image inputs
and outputs are already defined.   By doing it this way all of the nice pipeline properties
no longer valid.   The last two already exist in the hierarchy.

In another area of incompatibility with ITK/C++, the neighborhood
iterators which are used in MRFImageFilter are not used at all, and the "Gibbs" code takes
on a much more "C" flavor, with the same dangers that are associated with standard C
programming.  

4.  At the end of this message I've placed the MRFImageFilter driver.
 I've put an example input image, training
image, and output file on www.mypacs.net, case number 254930.
Unfortunately, there must be a small
amount of compression on the site because I cannot download the training image properly without
creating some changes in the values.   This is not so critical for viewing, but is critical
to describe the classes properly.  It is useful for describing the output, so I'll just leave
it there for now.     I borrowed as much code for the driver
as I could from RGBGibbsPrior1.cxx, but left out, for instance, the section where the
training image size was hardcoded in. ;-)   Probably it needs a few more comments too.

Greg Cook

MRFImageFilter1.cxx:

/*=========================================================================

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/

#ifdef _MSC_VER
#pragma warning ( disable : 4786 )
#endif

#include "itkMRFImageFilter.h"

// classes help the Gibbs filter to segment the image
#include "itkImageClassifierBase.h"
#include "itkImageGaussianModelEstimator.h"
#include "itkMahalanobisDistanceMembershipFunction.h"
#include "itkMinimumDecisionRule.h"

#include "itkScalarToArrayCastImageFilter.h"
#include "itkShiftScaleImageFilter.h"

#include "itkMinimumMaximumImageCalculator.h"

// image storage and I/O classes
#include "itkSize.h"
#include "itkImage.h"
#include "itkVector.h"
#include "itkImageRegionIterator.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

// itk console output routines
#include "itkMacro.h"
#include "itkOutputWindow.h"

int main( int argc, char *argv[] )
{
  if( argc != 4 )
    {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " inputImage trainingImage outputImage" << std::endl;
    return 1;
  }

  std::cout << "MRFImageFilter Test: " << std::endl;

  //  The input is a single channel image; the channel number is  
  //  \code{NUMBANDS} = 1, and \code{NDIMENSION} is set to 3.
  //

  const unsigned short NUMBANDS = 1;
  const unsigned short NDIMENSION = 3;

  typedef itk::Image<itk::Vector<unsigned short,NUMBANDS>,NDIMENSION> VectorImageType; 

  typedef itk::Image< short, NDIMENSION > ClassImageType;

  typedef itk::Image<short, NDIMENSION > InputImageType;
  typedef itk::Image<unsigned char, NDIMENSION > OutputImageType;

  // We instantiate reader and imageWriter types
  //
  typedef  itk::ImageFileReader< InputImageType > ReaderType;
  typedef  itk::ImageFileWriter<  OutputImageType  > WriterType;

  typedef itk::ScalarToArrayCastImageFilter<InputImageType,VectorImageType>
           ScalarToVectorFilterType;

  ReaderType::Pointer inputImageReader = ReaderType::New();
  ReaderType::Pointer trainingImageReader = ReaderType::New();
  WriterType::Pointer imageWriter = WriterType::New();

  ScalarToVectorFilterType::Pointer applyScalarToVectorFilter
   = ScalarToVectorFilterType::New();

  applyScalarToVectorFilter->SetInput(inputImageReader->GetOutput());

  inputImageReader->SetFileName( argv[1] );
  trainingImageReader->SetFileName( argv[2] );
  imageWriter->SetFileName( argv[3] );

  typedef itk::ShiftScaleImageFilter<ClassImageType,ClassImageType>
   ShiftFilterType;

  ShiftFilterType::Pointer applyShiftFilter = ShiftFilterType::New();

  applyShiftFilter->SetInput(trainingImageReader->GetOutput());
  applyShiftFilter->SetShift(-1);

  typedef itk::MinimumMaximumImageCalculator<ClassImageType> GetMaxTrainingValueType;
  GetMaxTrainingValueType::Pointer getMaxTrainingValue = GetMaxTrainingValueType::New();

  // MinimumMaximumImageCalculator is not a pipeline type, so need to update
  applyShiftFilter->Update();
  getMaxTrainingValue->SetImage(applyShiftFilter->GetOutput());
  getMaxTrainingValue->ComputeMaximum();

  std::cout << "Number of classes = "
         << getMaxTrainingValue->GetMaximum() + 1
   << std::endl;

  //----------------------------------------------------------------------
  //Set membership function (Using the statistics objects)
  //----------------------------------------------------------------------

  namespace stat = itk::Statistics;

  typedef VectorImageType::PixelType VectorImagePixelType;
  typedef stat::MahalanobisDistanceMembershipFunction< VectorImagePixelType > 
    MembershipFunctionType ;
  typedef MembershipFunctionType::Pointer MembershipFunctionPointer ;

  typedef std::vector< MembershipFunctionPointer > 
    MembershipFunctionPointerVector;  

  //----------------------------------------------------------------------
  // Set the image model estimator (train the class models)
  //----------------------------------------------------------------------

  typedef itk::ImageGaussianModelEstimator<VectorImageType,
    MembershipFunctionType, ClassImageType> 
    ImageGaussianModelEstimatorType;

  ImageGaussianModelEstimatorType::Pointer 
    applyEstimateModel = ImageGaussianModelEstimatorType::New();  

  applyEstimateModel->SetNumberOfModels(getMaxTrainingValue->GetMaximum()+1);
  applyEstimateModel->SetInputImage(applyScalarToVectorFilter->GetOutput());
  applyEstimateModel->SetTrainingImage(applyShiftFilter->GetOutput());  

// shouldn't need to do this here, but unfortunately the
// pipeline model isn't followed in ImageGaussianModelEstimator
  applyScalarToVectorFilter->Update();
  applyShiftFilter->Update();

  //Run the gaussian classifier algorithm
  applyEstimateModel->Update();

  std::cout << "Model Estimated" << std::endl;

  applyEstimateModel->Print(std::cout); 

  MembershipFunctionPointerVector membershipFunctions = 
    applyEstimateModel->GetMembershipFunctions();

  //----------------------------------------------------------------------
  //Set the decision rule 
  //----------------------------------------------------------------------  
  typedef itk::DecisionRuleBase::Pointer DecisionRuleBasePointer;

  typedef itk::MinimumDecisionRule DecisionRuleType;
  DecisionRuleType::Pointer  myDecisionRule = DecisionRuleType::New();

  // Set the classifier to be used and assign the parameters for the 
  // supervised classifier algorithm

  typedef VectorImagePixelType MeasurementVectorType;

  typedef itk::ImageClassifierBase< VectorImageType, ClassImageType > ClassifierType;
  typedef itk::ClassifierBase<VectorImageType>::Pointer ClassifierBasePointer;

  typedef ClassifierType::Pointer ClassifierPointer;
  ClassifierPointer myClassifier = ClassifierType::New();

  // Set the Classifier parameters
  myClassifier->SetNumberOfClasses(applyEstimateModel->GetNumberOfModels());

  // Set the decison rule 
  myClassifier->SetDecisionRule((DecisionRuleBasePointer) myDecisionRule );

  //Add the membership functions
  for( unsigned int i=0; i<myClassifier->GetNumberOfClasses(); i++ )
    {
    myClassifier->AddMembershipFunction( membershipFunctions[i] );
    }

  typedef itk::MRFImageFilter<VectorImageType,ClassImageType> MRFImageFilterType;

  MRFImageFilterType::Pointer applyMRFImageFilter
 = MRFImageFilterType::New();

  applyMRFImageFilter->SetInput(applyScalarToVectorFilter->GetOutput() );
  // shouldn't have to do this because the classifier *already* knows
  // the number of classes; need to fix this in itkMRFImageFilter
  applyMRFImageFilter->SetNumberOfClasses( myClassifier->GetNumberOfClasses());

  applyMRFImageFilter->SetClassifier( myClassifier );
  applyMRFImageFilter->SetNeighborhoodRadius( 1 );
  applyMRFImageFilter->SetErrorTolerance(0.05);  // default 0.2

// On output, put the shift back in if you want the exact class values
// that you started with in the training image.  SetShift(0) if you 
// don't.

  typedef itk::ShiftScaleImageFilter<ClassImageType,OutputImageType>
   OutputShiftFilterType;

  OutputShiftFilterType::Pointer applyOutputShiftFilter
   = OutputShiftFilterType::New();

  applyOutputShiftFilter->SetInput(applyMRFImageFilter->GetOutput());
  applyOutputShiftFilter->SetShift(+1);
  imageWriter->SetInput( applyOutputShiftFilter->GetOutput());

  // do it
  imageWriter->Update();

  std::cout << "MRF (ICM) Segmentation Completed" << std::endl;

  applyMRFImageFilter->Print(std::cout);

  return 0;
}