[Insight-users] Problem implementing GradientAnisotropicDiffusionImageFilter

[Bill Lorensen]

Can you post the matlab source code for the algorithm? Perhaps there
is a difference in parameters? Or implementation?

On Tue, Dec 1, 2009 at 8:48 AM, O'Connor, Michael
<Michael.OConnor at umassmed.edu> wrote:
>  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;
>
> }
>
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>
>