[Insight-users] Classification using Markov Random Field ITK software guide

[inparanoia]

i've done 2d abd it funcs,
in 3d the changes i should do for a float mhd and a segmentation by 
kdmeans mhd char:
put dimension 3
put float in InputImage
for a radius of neighbour 1 put 9 weights;
how do i put weights? whta order?(because i should put a matrix in a vector)
i've done some tries but using the means and same number of classes by 
kdmeans it stops at the first iteration


<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 ./ScalarImageMarkovRandomField1  original.mhd outputbykdmeans.mhd 
outit10.mhd 10 3 3 threeaverage by kdmeans

Number of Iterations : 1
Stop condition:
  (1) Maximum number of iterations
  (2) Error tolerance: 
2
<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

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



// Software Guide : BeginCodeSnippet
#include "itkImage.h"
#include "itkFixedArray.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkScalarToArrayCastImageFilter.h"
// Software Guide : EndCodeSnippet

#include "itkRescaleIntensityImageFilter.h"

// Software Guide : BeginLatex
//
// The following headers are related to the statistical classification 
classes.
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
#include "itkMRFImageFilter.h"
#include "itkDistanceToCentroidMembershipFunction.h"
#include "itkMinimumDecisionRule.h"
#include "itkImageClassifierBase.h"
// Software Guide : EndCodeSnippet

int main( int argc, char * argv [] )
{
  if( argc < 5 )
    {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0];
    std::cerr << " inputScalarImage inputLabeledImage";
    std::cerr << " outputLabeledImage numberOfIterations";
    std::cerr << " smoothingFactor numberOfClasses";
    std::cerr << " mean1 mean2 ... meanN " << std::endl;
    return EXIT_FAILURE;
    }

  const char * inputImageFileName      = argv[1];
  const char * inputLabelImageFileName = argv[2];
  const char * outputImageFileName     = argv[3];

  const unsigned int numberOfIterations = atoi( argv[4] );
  const double       smoothingFactor    = atof( argv[5] );
  const unsigned int numberOfClasses    = atoi( argv[6] );

  const unsigned int numberOfArgumentsBeforeMeans = 7;

  if( static_cast<unsigned int>(argc) <
      numberOfClasses + numberOfArgumentsBeforeMeans )
    {
    std::cerr << "Error: " << std::endl;
    std::cerr << numberOfClasses << " classes has been specified ";
    std::cerr << "but no enough means have been provided in the command ";
    std::cerr << "line arguments " << std::endl;
    return EXIT_FAILURE;
    }




  typedef float       PixelType;//----------------------
  const unsigned int          Dimension = 3;//---------------

  typedef itk::Image<PixelType, Dimension > ImageType;

  typedef itk::ImageFileReader< ImageType > ReaderType;
  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName( inputImageFileName );
// Software Guide : EndCodeSnippet




  typedef unsigned char       LabelPixelType;

  typedef itk::Image<LabelPixelType, Dimension > LabelImageType;

  typedef itk::ImageFileReader< LabelImageType > LabelReaderType;
  LabelReaderType::Pointer labelReader = LabelReaderType::New();
  labelReader->SetFileName( inputLabelImageFileName );
// Software Guide : EndCodeSnippet



  typedef itk::FixedArray<LabelPixelType,1>  ArrayPixelType;

  typedef itk::Image< ArrayPixelType, Dimension > ArrayImageType;

  typedef itk::ScalarToArrayCastImageFilter<
                     ImageType, ArrayImageType > ScalarToArrayFilterType;

  ScalarToArrayFilterType::Pointer
    scalarToArrayFilter = ScalarToArrayFilterType::New();
  scalarToArrayFilter->SetInput( reader->GetOutput() );
// Software Guide : EndCodeSnippet


  typedef itk::MRFImageFilter< ArrayImageType, LabelImageType > 
MRFFilterType;

  MRFFilterType::Pointer mrfFilter = MRFFilterType::New();

  mrfFilter->SetInput( scalarToArrayFilter->GetOutput() );
// Software Guide : EndCodeSnippet




  mrfFilter->SetNumberOfClasses( numberOfClasses );
  mrfFilter->SetMaximumNumberOfIterations( numberOfIterations );
  mrfFilter->SetErrorTolerance( 1e-7 );
// Software Guide : EndCodeSnippet



  mrfFilter->SetSmoothingFactor( smoothingFactor );
// Software Guide : EndCodeSnippet




  typedef itk::ImageClassifierBase<
                              ArrayImageType,
                              LabelImageType >   SupervisedClassifierType;

  SupervisedClassifierType::Pointer classifier =
                                         SupervisedClassifierType::New();
// Software Guide : EndCodeSnippet



  typedef itk::MinimumDecisionRule DecisionRuleType;

  DecisionRuleType::Pointer  classifierDecisionRule = 
DecisionRuleType::New();

  classifier->SetDecisionRule( classifierDecisionRule.GetPointer() );
// Software Guide : EndCodeSnippet


  typedef itk::Statistics::DistanceToCentroidMembershipFunction<
                                                    ArrayPixelType >
                                                       
MembershipFunctionType;

  typedef MembershipFunctionType::Pointer MembershipFunctionPointer;


  double meanDistance = 0;
  vnl_vector<double> centroid(1);
  for( unsigned int i=0; i < numberOfClasses; i++ )
    {
    MembershipFunctionPointer membershipFunction =
                                         MembershipFunctionType::New();

    centroid[0] = atof( argv[i+numberOfArgumentsBeforeMeans] );

    membershipFunction->SetCentroid( centroid );

    classifier->AddMembershipFunction( membershipFunction );
    meanDistance += static_cast< double > (centroid[0]);
    }
  meanDistance /= numberOfClasses;

  mrfFilter->SetSmoothingFactor( smoothingFactor );
// Software Guide : EndCodeSnippet





// Software Guide : BeginCodeSnippet
  mrfFilter->SetNeighborhoodRadius( 1 );
// Software Guide : EndCodeSnippet




//weights copied from .h don't know the order----------------------
// Software Guide : BeginCodeSnippet
  std::vector< double > weights;
  weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.5);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);

  weights.push_back(1.7);
  weights.push_back(1.7);
  weights.push_back(1.7);
  weights.push_back(1.7);
  weights.push_back(0.0);
  weights.push_back(1.7);
  weights.push_back(1.7);
  weights.push_back(1.7);
  weights.push_back(1.7);
 
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.5);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);
  
  /*weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(0.0);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);

  weights.push_back(1.3);
  weights.push_back(1.3);
  weights.push_back(1.3);
  weights.push_back(1.3);
  weights.push_back(1.3);
  weights.push_back(1.3);
  weights.push_back(1.3);
  weights.push_back(1.3);
  weights.push_back(1.3);
 
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);
   weights.push_back(1.3);*/


 

  double totalWeight = 0;
  for(std::vector< double >::const_iterator wcIt = weights.begin();
      wcIt != weights.end(); ++wcIt )
    {
    totalWeight += *wcIt;
    }
  for(std::vector< double >::iterator wIt = weights.begin();
      wIt != weights.end(); wIt++ )
    {
    *wIt = static_cast< double > ( (*wIt) * meanDistance / (3 * 
totalWeight));//two or three? in 2d 2,already tried 2 in 3d
    }

  mrfFilter->SetMRFNeighborhoodWeight( weights );
// Software Guide : EndCodeSnippet




// Software Guide : BeginLatex
//
// Finally, the classifier class is connected to the Markof Random 
Fields filter.
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
  mrfFilter->SetClassifier( classifier );
// Software Guide : EndCodeSnippet




  typedef MRFFilterType::OutputImageType  OutputImageType;
// Software Guide : EndCodeSnippet

  // Rescale outputs to the dynamic range of the display
  typedef itk::Image< unsigned char, Dimension > RescaledOutputImageType;
  typedef itk::RescaleIntensityImageFilter<
             OutputImageType, RescaledOutputImageType >   RescalerType;

  RescalerType::Pointer intensityRescaler = RescalerType::New();
  intensityRescaler->SetOutputMinimum(   0 );
  intensityRescaler->SetOutputMaximum( 255 );
  intensityRescaler->SetInput( mrfFilter->GetOutput() ); 

// Software Guide : BeginCodeSnippet
  typedef itk::ImageFileWriter< OutputImageType > WriterType;
  //typedef itk::ImageFileWriter< OutputImageType > WriterType;
  WriterType::Pointer writer = WriterType::New();
 
  //writer->SetInput( intensityRescaler->GetOutput() );
  writer->SetInput( mrfFilter->GetOutput() );
 
  writer->SetFileName( outputImageFileName );
// Software Guide : EndCodeSnippet





  try
    {
    writer->Update();
    }
  catch( itk::ExceptionObject & excp )
    {
    std::cerr << "Problem encountered while writing ";
    std::cerr << " image file : " << argv[2] << std::endl;
    std::cerr << excp << std::endl;
    return EXIT_FAILURE;
    }
// Software Guide : EndCodeSnippet

  std::cout << "Number of Iterations : ";
  std::cout << mrfFilter->GetNumberOfIterations() << std::endl;
  std::cout << "Stop condition: " << std::endl;
  std::cout << "  (1) Maximum number of iterations " << std::endl;
  std::cout << "  (2) Error tolerance:  "  << std::endl;
  std::cout << mrfFilter->GetStopCondition() << std::endl;

 

  return EXIT_SUCCESS;
 
}