[Insight-users] Problem implementing GradientAnisotropicDiffusionImageFilter

[O'Connor, Michael]

 Hi,
   I am new to ITK.  I have just started changing some of my post-processing tools from MATLAB to ITK (Version 3.12).  So far I have only had one issue that I cannot solve and that is implementing GradientAnisotropicDiffusionImageFilter.  I can workaround this issue but would rather get to the bottom of my problem as I suspect that I am missing something fundamental about data type in ITK.  I have read pertinent sections of ITK Software Guide 2nd Edition (dated 11/21/05) and template guide on www.itk.org.
    I am processing data obtained on a prototype CT system and am processing in order to eventually segment between one tissue type (approx linear attenuation coefficient 0.024 per mm) from another (approx attenuation coefficient 0.031 per mm).  Setting parameters of 5 iterations, timestep of 0.125 and conductance 0.003, there is minimal filtering when using following ITK program (i.e. not consistent with conductance parameter).   Using same values and same test input file in MATLAB, I do get reasonable output.   If I take another test input file scaled up by 10000 (both data values and conductance - i.e. approx data coefficients 240 per mm and 310 per mm and conductance 30), I do obtain reasonable results with the same ITK program (as well as MATLAB program with only change in conductance).    
   As pixel type is float in either case, I cannot understand why the GradientAnisotropicDiffusionImageFilter does not work with both range of data values in ITK as it does in MATLAB.  While I can perform scaling to workaround my ITK issue, as stated, I suspect that I am doing something wrong in ITK and would rather correct my misunderstanding.   I would greatly appreciate any comments or suggestions.

Regards,
   Michael O'Connor

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

  Program:   ITKFilterTest
  Module:    ITKFilterTest.cxx
  Language:  C++
  Date:      24 Nov 2009
  Version:   Revision: 1.0   Mike O'Connor

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

// NOTE: much of this code and description is taken directly from ITK documentation. Module initially based
//  on the Filter Example modules in ITK source code.
//  The module runs one ITK filter.
//
//  Filter 1: is ITK GradientAnisotropicDiffusionImageFilter. This filter  implements an
//    N-dimensional version of the classic Perona-Malik anisotropic diffusion
//    equation for scalar-valued images (Perona 1990).
//
//    The conductance term for this implementation is chosen as a function of the
//    gradient magnitude of the image at each point, reducing the strength of
//    diffusion at edge pixels.
//
//    The numerical implementation of this equation is similar to that described
//    in the Perona-Malik paper, but uses a more robust technique
//    for gradient magnitude estimation and has been generalized to N-dimensions.
//  Input Parameters for this implementation are passed using a 'config file'.  Objects and methods for the 
//  config file are based on CBSS.  There are booleans to determine if one or both of the filters are run.
//  The input parameters for each filter as well as parameters for file reader.


#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkRawImageIO.h"
#include "itkGradientAnisotropicDiffusionImageFilter.h"
#include "FilterSystem.h"
#include "Filter.h"
#include "Files.h"
#include "Message.h"



const unsigned int Dimension = 2;

int main( int argc, char * argv[] )
{
 // This is structure for our config file (used for parameter passing)
  FilterSystem f_system;
  WhatFunction function = Filter;  
  if(parseArgs(argc,argv) == UtilError) return EXIT_FAILURE;
  if(f_system.initialize(function, usage) == UtilError)closeAndExit(EXIT_FAILURE);

  typedef    float    InputPixelType;
  typedef    float    OutputPixelType;

//
// Instatiate using input & output image types
//
  typedef itk::Image< 
			InputPixelType,  
			Dimension >   		InputImageType;
  typedef itk::Image< 
			OutputPixelType, 
			Dimension >   		OutputImageType;
  typedef itk::ImageFileReader< 
			InputImageType > 	ReaderType;
  typedef itk::GradientAnisotropicDiffusionImageFilter<
               		InputImageType, 
			OutputImageType >  	ADF_FilterType;
  typedef itk::RawImageIO<
			InputPixelType, 
			Dimension	>	ImageIOType;
  typedef itk::ImageFileWriter< 
			OutputImageType >  	WriterType;                                                                                                                                                  
//
// Filter object, Image object and Reader & Writer object are  created
//
  ImageIOType::Pointer 	rawIO 		= ImageIOType::New();
  ReaderType::Pointer 	reader 		= ReaderType::New();
  ADF_FilterType::Pointer ADfilter 	= ADF_FilterType::New();
  WriterType::Pointer 	writer2		= WriterType::New();
	
// Get parameters from config file
  rawIO->SetHeaderSize(f_system.m_cfg.header);    
  rawIO->SetDimensions(0, f_system.m_cfg.file_dim_x);
  rawIO->SetDimensions(1, f_system.m_cfg.file_dim_y);
  rawIO->SetByteOrderToLittleEndian();

  // Define Input file
  writeMsg(DEBUG1, "Input File: %s \n",files.input);
  reader->SetFileName( files.input );
  reader->SetImageIO( rawIO );
//  reader->Update();		//Read the file
  //
  //Get other parameters from Config file
  const unsigned int 	numberOfIterations =  f_system.m_cfg.iter;
  const double	  	timeStep = (double) f_system.m_cfg.timestep; 
  const double       	conductance = (double) f_system.m_cfg.conductance;
  writeMsg(DEBUG1,"Version 2: Iterations = %i, Timestep = %f,  Conductance = %f \n", numberOfIterations, timeStep, conductance);

  //Set Parameters in ADF Filter & invoke filter
  ADfilter->SetInput( reader->GetOutput() );
  ADfilter->SetNumberOfIterations( numberOfIterations );
  ADfilter->SetTimeStep( timeStep );
  ADfilter->SetConductanceParameter( conductance );
  ADfilter->Update();

  //
  // Write the output file
  //
  writer2->SetImageIO (rawIO);
  writer2->SetFileName( files.output );
  writer2->SetInput( ADfilter->GetOutput() );
  writeMsg(DEBUG1, "Writing output of the AnisoDiff to output File: %s \n",files.output);
  writer2->Update();
  return EXIT_SUCCESS;

}

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091201/23a86a44/attachment.htm>