<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2//EN">
<HTML>
<HEAD>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=iso-8859-1">
<META NAME="Generator" CONTENT="MS Exchange Server version 6.5.7654.12">
<TITLE>Problem implementing GradientAnisotropicDiffusionImageFilter</TITLE>
</HEAD>
<BODY>
<!-- Converted from text/plain format -->

<P><FONT SIZE=2>&nbsp;Hi,<BR>
&nbsp;&nbsp; I am new to ITK.&nbsp; I have just started changing some of my post-processing tools from MATLAB to ITK (Version 3.12).&nbsp; So far I have only had one issue that I cannot solve and that is implementing GradientAnisotropicDiffusionImageFilter.&nbsp; 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.&nbsp; I have read pertinent sections of ITK Software Guide 2nd Edition (dated 11/21/05) and template guide on www.itk.org.<BR>
&nbsp;&nbsp;&nbsp; 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).&nbsp; 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).&nbsp;&nbsp; Using same values and same test input file in MATLAB, I do get reasonable output.&nbsp;&nbsp; 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).&nbsp;&nbsp;&nbsp;<BR>
&nbsp;&nbsp; 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.&nbsp; 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.&nbsp;&nbsp; I would greatly appreciate any comments or suggestions.<BR>
<BR>
Regards,<BR>
&nbsp;&nbsp; Michael O'Connor<BR>
<BR>
/*=========================================================================<BR>
<BR>
&nbsp; Program:&nbsp;&nbsp; ITKFilterTest<BR>
&nbsp; Module:&nbsp;&nbsp;&nbsp; ITKFilterTest.cxx<BR>
&nbsp; Language:&nbsp; C++<BR>
&nbsp; Date:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 24 Nov 2009<BR>
&nbsp; Version:&nbsp;&nbsp; Revision: 1.0&nbsp;&nbsp; Mike O'Connor<BR>
<BR>
=========================================================================*/<BR>
<BR>
// NOTE: much of this code and description is taken directly from ITK documentation. Module initially based<BR>
//&nbsp; on the Filter Example modules in ITK source code.<BR>
//&nbsp; The module runs one ITK filter.<BR>
//<BR>
//&nbsp; Filter 1: is ITK GradientAnisotropicDiffusionImageFilter. This filter&nbsp; implements an<BR>
//&nbsp;&nbsp;&nbsp; N-dimensional version of the classic Perona-Malik anisotropic diffusion<BR>
//&nbsp;&nbsp;&nbsp; equation for scalar-valued images (Perona 1990).<BR>
//<BR>
//&nbsp;&nbsp;&nbsp; The conductance term for this implementation is chosen as a function of the<BR>
//&nbsp;&nbsp;&nbsp; gradient magnitude of the image at each point, reducing the strength of<BR>
//&nbsp;&nbsp;&nbsp; diffusion at edge pixels.<BR>
//<BR>
//&nbsp;&nbsp;&nbsp; The numerical implementation of this equation is similar to that described<BR>
//&nbsp;&nbsp;&nbsp; in the Perona-Malik paper, but uses a more robust technique<BR>
//&nbsp;&nbsp;&nbsp; for gradient magnitude estimation and has been generalized to N-dimensions.<BR>
//&nbsp; Input Parameters for this implementation are passed using a 'config file'.&nbsp; Objects and methods for the<BR>
//&nbsp; config file are based on CBSS.&nbsp; There are booleans to determine if one or both of the filters are run.<BR>
//&nbsp; The input parameters for each filter as well as parameters for file reader.<BR>
<BR>
<BR>
#include &quot;itkImage.h&quot;<BR>
#include &quot;itkImageFileReader.h&quot;<BR>
#include &quot;itkImageFileWriter.h&quot;<BR>
#include &quot;itkRawImageIO.h&quot;<BR>
#include &quot;itkGradientAnisotropicDiffusionImageFilter.h&quot;<BR>
#include &quot;FilterSystem.h&quot;<BR>
#include &quot;Filter.h&quot;<BR>
#include &quot;Files.h&quot;<BR>
#include &quot;Message.h&quot;<BR>
<BR>
<BR>
<BR>
const unsigned int Dimension = 2;<BR>
<BR>
int main( int argc, char * argv[] )<BR>
{<BR>
&nbsp;// This is structure for our config file (used for parameter passing)<BR>
&nbsp; FilterSystem f_system;<BR>
&nbsp; WhatFunction function = Filter;&nbsp;<BR>
&nbsp; if(parseArgs(argc,argv) == UtilError) return EXIT_FAILURE;<BR>
&nbsp; if(f_system.initialize(function, usage) == UtilError)closeAndExit(EXIT_FAILURE);<BR>
<BR>
&nbsp; typedef&nbsp;&nbsp;&nbsp; float&nbsp;&nbsp;&nbsp; InputPixelType;<BR>
&nbsp; typedef&nbsp;&nbsp;&nbsp; float&nbsp;&nbsp;&nbsp; OutputPixelType;<BR>
<BR>
//<BR>
// Instatiate using input &amp; output image types<BR>
//<BR>
&nbsp; typedef itk::Image&lt;<BR>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; InputPixelType,&nbsp;<BR>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Dimension &gt;&nbsp;&nbsp; &nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; InputImageType;<BR>
&nbsp; typedef itk::Image&lt;<BR>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OutputPixelType,<BR>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Dimension &gt;&nbsp;&nbsp; &nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OutputImageType;<BR>
&nbsp; typedef itk::ImageFileReader&lt;<BR>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; InputImageType &gt; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ReaderType;<BR>
&nbsp; typedef itk::GradientAnisotropicDiffusionImageFilter&lt;<BR>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; InputImageType,<BR>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OutputImageType &gt;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp; ADF_FilterType;<BR>
&nbsp; typedef itk::RawImageIO&lt;<BR>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; InputPixelType,<BR>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Dimension&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &gt;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ImageIOType;<BR>
&nbsp; typedef itk::ImageFileWriter&lt;<BR>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OutputImageType &gt;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp; WriterType;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<BR>
//<BR>
// Filter object, Image object and Reader &amp; Writer object are&nbsp; created<BR>
//<BR>
&nbsp; ImageIOType::Pointer&nbsp; rawIO &nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = ImageIOType::New();<BR>
&nbsp; ReaderType::Pointer &nbsp; reader&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = ReaderType::New();<BR>
&nbsp; ADF_FilterType::Pointer ADfilter &nbsp;&nbsp;&nbsp;&nbsp; = ADF_FilterType::New();<BR>
&nbsp; WriterType::Pointer &nbsp; writer2 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = WriterType::New();<BR>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<BR>
// Get parameters from config file<BR>
&nbsp; rawIO-&gt;SetHeaderSize(f_system.m_cfg.header);&nbsp;&nbsp;&nbsp;<BR>
&nbsp; rawIO-&gt;SetDimensions(0, f_system.m_cfg.file_dim_x);<BR>
&nbsp; rawIO-&gt;SetDimensions(1, f_system.m_cfg.file_dim_y);<BR>
&nbsp; rawIO-&gt;SetByteOrderToLittleEndian();<BR>
<BR>
&nbsp; // Define Input file<BR>
&nbsp; writeMsg(DEBUG1, &quot;Input File: %s \n&quot;,files.input);<BR>
&nbsp; reader-&gt;SetFileName( files.input );<BR>
&nbsp; reader-&gt;SetImageIO( rawIO );<BR>
//&nbsp; reader-&gt;Update();&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; //Read the file<BR>
&nbsp; //<BR>
&nbsp; //Get other parameters from Config file<BR>
&nbsp; const unsigned int &nbsp;&nbsp; numberOfIterations =&nbsp; f_system.m_cfg.iter;<BR>
&nbsp; const double&nbsp; &nbsp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp; timeStep = (double) f_system.m_cfg.timestep;<BR>
&nbsp; const double&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp; conductance = (double) f_system.m_cfg.conductance;<BR>
&nbsp; writeMsg(DEBUG1,&quot;Version 2: Iterations = %i, Timestep = %f,&nbsp; Conductance = %f \n&quot;, numberOfIterations, timeStep, conductance);<BR>
<BR>
&nbsp; //Set Parameters in ADF Filter &amp; invoke filter<BR>
&nbsp; ADfilter-&gt;SetInput( reader-&gt;GetOutput() );<BR>
&nbsp; ADfilter-&gt;SetNumberOfIterations( numberOfIterations );<BR>
&nbsp; ADfilter-&gt;SetTimeStep( timeStep );<BR>
&nbsp; ADfilter-&gt;SetConductanceParameter( conductance );<BR>
&nbsp; ADfilter-&gt;Update();<BR>
<BR>
&nbsp; //<BR>
&nbsp; // Write the output file<BR>
&nbsp; //<BR>
&nbsp; writer2-&gt;SetImageIO (rawIO);<BR>
&nbsp; writer2-&gt;SetFileName( files.output );<BR>
&nbsp; writer2-&gt;SetInput( ADfilter-&gt;GetOutput() );<BR>
&nbsp; writeMsg(DEBUG1, &quot;Writing output of the AnisoDiff to output File: %s \n&quot;,files.output);<BR>
&nbsp; writer2-&gt;Update();<BR>
&nbsp; return EXIT_SUCCESS;<BR>
<BR>
}<BR>
<BR>
</FONT>
</P>

</BODY>
</HTML>