<!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>RE: [Insight-users] Problem implementing GradientAnisotropicDiffusionImageFilter</TITLE>
</HEAD>
<BODY>
<!-- Converted from text/plain format -->
<P><FONT SIZE=2>Hi Bill,<BR>
I use the exact same file with ITK (the one attached to the previous email with MATLAB code).<BR>
<BR>
Regards,<BR>
Michael<BR>
<BR>
<BR>
-----Original Message-----<BR>
From: Bill Lorensen [<A HREF="mailto:bill.lorensen@gmail.com">mailto:bill.lorensen@gmail.com</A>]<BR>
Sent: Tue 12/1/2009 4:17 PM<BR>
To: O'Connor, Michael<BR>
Cc: insight-users@itk.org<BR>
Subject: Re: [Insight-users] Problem implementing GradientAnisotropicDiffusionImageFilter<BR>
<BR>
Can you provide a png file or meta file, or whatever file you tried with itk?<BR>
<BR>
On Tue, Dec 1, 2009 at 4:08 PM, O'Connor, Michael<BR>
<Michael.OConnor@umassmed.edu> wrote:<BR>
> Hi Bill,<BR>
> Thanks. Attached is simple MATLAB test code + the MATLAB filter<BR>
> function used (Dr. Guy Gilboa's code). I've also take the liberty of<BR>
> attaching a test image (specifically opened in MATLAB code). The test<BR>
> image is 400x400 float with background = 0.024 + additive Gaussian Noise and<BR>
> foreground (circle) = 0.031 + additive Gaussian Noise. As said previously<BR>
> MATLAB code seems to work as expected as shown in plot. In ITK there is no<BR>
> similar intraregion filtering (i.e. filtered image ~ input )<BR>
><BR>
> Regards,<BR>
> Michael<BR>
><BR>
> ******<BR>
> %%% simple test of anisotropic diffusion filter<BR>
> %%% uses Guy Gilboa's diffusion function (from his website)<BR>
> %%%<BR>
> clear all;<BR>
> close all;<BR>
> close all hidden;<BR>
> row = 400;<BR>
> col = 400;<BR>
> GAIN = 1;<BR>
> iter = 5;<BR>
> conductance = 0.003 * GAIN;<BR>
> timestep = 0.125;<BR>
><BR>
><BR>
> %%<BR>
> fid_raw = fopen('/Users/michaeloconnor/Desktop/phantom_ADF3.rec','r');<BR>
><BR>
> if (fid_raw == -1 )<BR>
> fprintf('NO Input File(s) in Working Directory \n');<BR>
> break; %stop now - invalid file<BR>
> end<BR>
> %%<BR>
><BR>
> fprintf('iterations = %i timestep = %4.3f conductance = %f Gain = %i<BR>
> \n',iter,timestep,conductance, GAIN)<BR>
><BR>
> read_raw= fread(fid_raw, [row col],'single');<BR>
> I=read_raw*GAIN;<BR>
><BR>
> J_adf=diffusion2(I, 'pm1', iter, conductance, timestep);<BR>
><BR>
> %% Display Both input and Filtered Using Image Processing Toolkit<BR>
><BR>
> imtool(J_adf, [])<BR>
> imtool(I, [])<BR>
><BR>
> %% Plot 1D Profile<BR>
><BR>
> p=plot([1:400],I(200,:),[1:400],J_adf(200,:));<BR>
> set(p,'LineWidth',2);<BR>
> title('1D Profile (row 200)','FontSize',14, 'FontWeight', 'bold');<BR>
> xlabel('Column (1:400)','FontSize',14, 'FontWeight', 'bold')<BR>
> ylabel('MU Value','FontSize',14, 'FontWeight', 'bold')<BR>
> legend('Pre-Filter', 'ADF Filter')<BR>
> fclose all;<BR>
><BR>
> ******* filter function (Dr. Gilboa's code)<BR>
> function Jd = diffusion2(J,method,N,K,dt,theta,sigma2)<BR>
> % private function: diffusion (by Guy Gilboa):<BR>
> % Jd=diffusion(J,method,N,K)<BR>
> % Simulates N iterations of diffusion, parameters:<BR>
> % J = source image (2D gray-level matrix) for diffusio<BR>
> % method = 'lin': Linear diffusion (constant c=1).<BR>
> % 'pm1': perona-malik, c=exp{-(|grad(J)|/K)^2} [PM90]<BR>
> % 'pm2': perona-malik, c=1/{1+(|grad(J)|/K)^2} [PM90]<BR>
> % 'rmp': complex valued - ramp preserving [GSZ01]<BR>
> % K edge threshold parameter<BR>
> % N number of iterations<BR>
> % dt time increment (0 < dt <= 0.25, default 0.2)<BR>
> % sigma2 - if present, calculates gradients of diffusion coefficient<BR>
> % convolved with gaussian of var sigma2 (Catte et al [CLMC92])<BR>
><BR>
> if ~exist('N','var')<BR>
> %fprintf('No iterations\n')<BR>
> N=1;<BR>
> end<BR>
> if ~exist('K','var')<BR>
> % fprintf('No edge threshold\n')<BR>
> K=1;<BR>
> end<BR>
> if ~exist('dt','var')<BR>
> % fprintf('No time increment\n')<BR>
> dt=0.2;<BR>
> end<BR>
> if ~exist('theta','var')<BR>
> % fprintf('No theta\n')<BR>
> theta=pi/30;<BR>
> end<BR>
><BR>
> if ~exist('sigma2','var')<BR>
> %fprintf('No gradients of diffusion\n')<BR>
> sigma2=0;<BR>
> end<BR>
><BR>
> [Ny,Nx]=size(J);<BR>
><BR>
> if (nargin<3)<BR>
> error('not enough arguments (at least 3 should be given)');<BR>
> end<BR>
><BR>
> for i=1:N;<BR>
> % gaussian filter with kernel 5x5 (Catte et al)<BR>
> if (sigma2>0)<BR>
> Jo = J; % save J original<BR>
> J=gauss(J,7,sigma2);<BR>
> end<BR>
><BR>
> % calculate gradient in all directions (N,S,E,W)<BR>
> In=[J(1,:); J(1:Ny-1,:)]-J;<BR>
> Is=[J(2:Ny,:); J(Ny,:)]-J;<BR>
> Ie=[J(:,2:Nx) J(:,Nx)]-J;<BR>
> Iw=[J(:,1) J(:,1:Nx-1)]-J;<BR>
><BR>
> % calculate diffusion coefficients in all directions according to<BR>
> method<BR>
> if (method == 'lin')<BR>
> Cn=K; Cs=K; Ce=K; Cw=K;<BR>
> elseif (method == 'pm1')<BR>
> Cn=exp(-(abs(In)/K).^2);<BR>
> Cs=exp(-(abs(Is)/K).^2);<BR>
> Ce=exp(-(abs(Ie)/K).^2);<BR>
> Cw=exp(-(abs(Iw)/K).^2);<BR>
> elseif (method == 'pm2')<BR>
> Cn=1./(1+(abs(In)/K).^2);<BR>
> Cs=1./(1+(abs(Is)/K).^2);<BR>
> Ce=1./(1+(abs(Ie)/K).^2);<BR>
> Cw=1./(1+(abs(Iw)/K).^2);<BR>
> elseif (method == 'rmp') % complex - ramp preserving<BR>
> k=K; j=sqrt(-1);<BR>
><BR>
> Cn=exp(j*theta)./(1+(imag(In)/(k*theta)).^2);<BR>
> Cs=exp(j*theta)./(1+(imag(Is)/(k*theta)).^2);<BR>
> Ce=exp(j*theta)./(1+(imag(Ie)/(k*theta)).^2);<BR>
> Cw=exp(j*theta)./(1+(imag(Iw)/(k*theta)).^2);<BR>
><BR>
> else<BR>
> error(['Unknown method "' method '"']);<BR>
> end<BR>
><BR>
> if (sigma2>0) % calculate real gradiants (not smoothed)<BR>
> In=[Jo(1,:); Jo(1:Ny-1,:)]-Jo;<BR>
> Is=[Jo(2:Ny,:); Jo(Ny,:)]-Jo;<BR>
> Ie=[Jo(:,2:Nx) Jo(:,Nx)]-Jo;<BR>
> Iw=[Jo(:,1) Jo(:,1:Nx-1)]-Jo;<BR>
> J=Jo;<BR>
> end<BR>
><BR>
> % Next Image J<BR>
> J=J+dt*(Cn.*In + Cs.*Is + Ce.*Ie + Cw.*Iw);<BR>
><BR>
> end; % for i<BR>
><BR>
> Jd = J;<BR>
><BR>
><BR>
> %%%% Refs:<BR>
> % [PM90] P. Perona, J. Malik, "Scale-space and edge detection using<BR>
> anisotropic diffusion", PAMI 12(7), pp. 629-639, 1990.<BR>
> % [CLMC92] F. Catte, P. L. Lions, J. M. Morel and T. Coll, "Image selective<BR>
> smoothing and edge detection by nonlinear diffusion", SIAM J. Num. Anal.,<BR>
> vol. 29, no. 1, pp. 182-193, 1992.<BR>
> % [GSZ01] G. Gilboa, N. Sochen, Y. Y. Zeevi, "Complex Diffusion Processes<BR>
> for Image Filtering", Scale-Space 2001, LNCS 2106, pp. 299-307,<BR>
> Springer-Verlag 2001.<BR>
><BR>
><BR>
><BR>
> -----Original Message-----<BR>
> From: Bill Lorensen [<A HREF="mailto:bill.lorensen@gmail.com">mailto:bill.lorensen@gmail.com</A>]<BR>
> Sent: Tue 12/1/2009 1:01 PM<BR>
> To: O'Connor, Michael<BR>
> Cc: insight-users@itk.org<BR>
> Subject: Re: [Insight-users] Problem implementing<BR>
> GradientAnisotropicDiffusionImageFilter<BR>
><BR>
> Can you post the matlab source code for the algorithm? Perhaps there<BR>
> is a difference in parameters? Or implementation?<BR>
><BR>
> On Tue, Dec 1, 2009 at 8:48 AM, O'Connor, Michael<BR>
> <Michael.OConnor@umassmed.edu> wrote:<BR>
>> Hi,<BR>
>> I am new to ITK. I have just started changing some of my<BR>
>> post-processing<BR>
>> tools from MATLAB to ITK (Version 3.12). So far I have only had one issue<BR>
>> that I cannot solve and that is implementing<BR>
>> GradientAnisotropicDiffusionImageFilter. I can workaround this issue but<BR>
>> would rather get to the bottom of my problem as I suspect that I am<BR>
>> missing<BR>
>> something fundamental about data type in ITK. I have read pertinent<BR>
>> sections of ITK Software Guide 2nd Edition (dated 11/21/05) and template<BR>
>> guide on www.itk.org.<BR>
>> I am processing data obtained on a prototype CT system and am<BR>
>> processing<BR>
>> in order to eventually segment between one tissue type (approx linear<BR>
>> attenuation coefficient 0.024 per mm) from another (approx attenuation<BR>
>> coefficient 0.031 per mm). Setting parameters of 5 iterations, timestep<BR>
>> of<BR>
>> 0.125 and conductance 0.003, there is minimal filtering when using<BR>
>> following<BR>
>> ITK program (i.e. not consistent with conductance parameter). Using same<BR>
>> values and same test input file in MATLAB, I do get reasonable output.<BR>
>> If<BR>
>> I take another test input file scaled up by 10000 (both data values and<BR>
>> conductance - i.e. approx data coefficients 240 per mm and 310 per mm and<BR>
>> conductance 30), I do obtain reasonable results with the same ITK program<BR>
>> (as well as MATLAB program with only change in conductance).<BR>
>> As pixel type is float in either case, I cannot understand why the<BR>
>> GradientAnisotropicDiffusionImageFilter does not work with both range of<BR>
>> data values in ITK as it does in MATLAB. While I can perform scaling to<BR>
>> workaround my ITK issue, as stated, I suspect that I am doing something<BR>
>> wrong in ITK and would rather correct my misunderstanding. I would<BR>
>> greatly<BR>
>> appreciate any comments or suggestions.<BR>
>><BR>
>> Regards,<BR>
>> Michael O'Connor<BR>
>><BR>
>><BR>
>> /*=========================================================================<BR>
>><BR>
>> Program: ITKFilterTest<BR>
>> Module: ITKFilterTest.cxx<BR>
>> Language: C++<BR>
>> Date: 24 Nov 2009<BR>
>> Version: Revision: 1.0 Mike O'Connor<BR>
>><BR>
>><BR>
>> =========================================================================*/<BR>
>><BR>
>> // NOTE: much of this code and description is taken directly from ITK<BR>
>> documentation. Module initially based<BR>
>> // on the Filter Example modules in ITK source code.<BR>
>> // The module runs one ITK filter.<BR>
>> //<BR>
>> // Filter 1: is ITK GradientAnisotropicDiffusionImageFilter. This filter<BR>
>> implements an<BR>
>> // N-dimensional version of the classic Perona-Malik anisotropic<BR>
>> diffusion<BR>
>> // equation for scalar-valued images (Perona 1990).<BR>
>> //<BR>
>> // The conductance term for this implementation is chosen as a function<BR>
>> of the<BR>
>> // gradient magnitude of the image at each point, reducing the strength<BR>
>> of<BR>
>> // diffusion at edge pixels.<BR>
>> //<BR>
>> // The numerical implementation of this equation is similar to that<BR>
>> described<BR>
>> // in the Perona-Malik paper, but uses a more robust technique<BR>
>> // for gradient magnitude estimation and has been generalized to<BR>
>> N-dimensions.<BR>
>> // Input Parameters for this implementation are passed using a 'config<BR>
>> file'. Objects and methods for the<BR>
>> // config file are based on CBSS. There are booleans to determine if one<BR>
>> or both of the filters are run.<BR>
>> // The input parameters for each filter as well as parameters for file<BR>
>> reader.<BR>
>><BR>
>><BR>
>> #include "itkImage.h"<BR>
>> #include "itkImageFileReader.h"<BR>
>> #include "itkImageFileWriter.h"<BR>
>> #include "itkRawImageIO.h"<BR>
>> #include "itkGradientAnisotropicDiffusionImageFilter.h"<BR>
>> #include "FilterSystem.h"<BR>
>> #include "Filter.h"<BR>
>> #include "Files.h"<BR>
>> #include "Message.h"<BR>
>><BR>
>><BR>
>><BR>
>> const unsigned int Dimension = 2;<BR>
>><BR>
>> int main( int argc, char * argv[] )<BR>
>> {<BR>
>> // This is structure for our config file (used for parameter passing)<BR>
>> FilterSystem f_system;<BR>
>> WhatFunction function = Filter;<BR>
>> if(parseArgs(argc,argv) == UtilError) return EXIT_FAILURE;<BR>
>> if(f_system.initialize(function, usage) ==<BR>
>> UtilError)closeAndExit(EXIT_FAILURE);<BR>
>><BR>
>> typedef float InputPixelType;<BR>
>> typedef float OutputPixelType;<BR>
>><BR>
>> //<BR>
>> // Instatiate using input & output image types<BR>
>> //<BR>
>> typedef itk::Image<<BR>
>> InputPixelType,<BR>
>> Dimension > InputImageType;<BR>
>> typedef itk::Image<<BR>
>> OutputPixelType,<BR>
>> Dimension > OutputImageType;<BR>
>> typedef itk::ImageFileReader<<BR>
>> InputImageType > ReaderType;<BR>
>> typedef itk::GradientAnisotropicDiffusionImageFilter<<BR>
>> InputImageType,<BR>
>> OutputImageType > ADF_FilterType;<BR>
>> typedef itk::RawImageIO<<BR>
>> InputPixelType,<BR>
>> Dimension > ImageIOType;<BR>
>> typedef itk::ImageFileWriter<<BR>
>> OutputImageType ><BR>
>> WriterType;<BR>
>> //<BR>
>> // Filter object, Image object and Reader & Writer object are created<BR>
>> //<BR>
>> ImageIOType::Pointer rawIO = ImageIOType::New();<BR>
>> ReaderType::Pointer reader = ReaderType::New();<BR>
>> ADF_FilterType::Pointer ADfilter = ADF_FilterType::New();<BR>
>> WriterType::Pointer writer2 = WriterType::New();<BR>
>><BR>
>> // Get parameters from config file<BR>
>> rawIO->SetHeaderSize(f_system.m_cfg.header);<BR>
>> rawIO->SetDimensions(0, f_system.m_cfg.file_dim_x);<BR>
>> rawIO->SetDimensions(1, f_system.m_cfg.file_dim_y);<BR>
>> rawIO->SetByteOrderToLittleEndian();<BR>
>><BR>
>> // Define Input file<BR>
>> writeMsg(DEBUG1, "Input File: %s \n",files.input);<BR>
>> reader->SetFileName( files.input );<BR>
>> reader->SetImageIO( rawIO );<BR>
>> // reader->Update(); //Read the file<BR>
>> //<BR>
>> //Get other parameters from Config file<BR>
>> const unsigned int numberOfIterations = f_system.m_cfg.iter;<BR>
>> const double timeStep = (double) f_system.m_cfg.timestep;<BR>
>> const double conductance = (double) f_system.m_cfg.conductance;<BR>
>> writeMsg(DEBUG1,"Version 2: Iterations = %i, Timestep = %f, Conductance<BR>
>> =<BR>
>> %f \n", numberOfIterations, timeStep, conductance);<BR>
>><BR>
>> //Set Parameters in ADF Filter & invoke filter<BR>
>> ADfilter->SetInput( reader->GetOutput() );<BR>
>> ADfilter->SetNumberOfIterations( numberOfIterations );<BR>
>> ADfilter->SetTimeStep( timeStep );<BR>
>> ADfilter->SetConductanceParameter( conductance );<BR>
>> ADfilter->Update();<BR>
>><BR>
>> //<BR>
>> // Write the output file<BR>
>> //<BR>
>> writer2->SetImageIO (rawIO);<BR>
>> writer2->SetFileName( files.output );<BR>
>> writer2->SetInput( ADfilter->GetOutput() );<BR>
>> writeMsg(DEBUG1, "Writing output of the AnisoDiff to output File: %s<BR>
>> \n",files.output);<BR>
>> writer2->Update();<BR>
>> return EXIT_SUCCESS;<BR>
>><BR>
>> }<BR>
>><BR>
>><BR>
>> _____________________________________<BR>
>> Powered by www.kitware.com<BR>
>><BR>
>> Visit other Kitware open-source projects at<BR>
>> <A HREF="http://www.kitware.com/opensource/opensource.html">http://www.kitware.com/opensource/opensource.html</A><BR>
>><BR>
>> Kitware offers ITK Training Courses, for more information visit:<BR>
>> <A HREF="http://www.kitware.com/products/protraining.html">http://www.kitware.com/products/protraining.html</A><BR>
>><BR>
>> Please keep messages on-topic and check the ITK FAQ at:<BR>
>> <A HREF="http://www.itk.org/Wiki/ITK_FAQ">http://www.itk.org/Wiki/ITK_FAQ</A><BR>
>><BR>
>> Follow this link to subscribe/unsubscribe:<BR>
>> <A HREF="http://www.itk.org/mailman/listinfo/insight-users">http://www.itk.org/mailman/listinfo/insight-users</A><BR>
>><BR>
>><BR>
><BR>
><BR>
><BR>
><BR>
<BR>
<BR>
</FONT>
</P>
</BODY>
</HTML>