[Insight-users] Problem implementing GradientAnisotropicDiffusionImageFilter
O'Connor, Michael
Michael.OConnor at umassmed.edu
Tue Dec 1 17:17:30 EST 2009
Hi Bill,
I use the exact same file with ITK (the one attached to the previous email with MATLAB code).
Regards,
Michael
-----Original Message-----
From: Bill Lorensen [mailto:bill.lorensen at gmail.com]
Sent: Tue 12/1/2009 4:17 PM
To: O'Connor, Michael
Cc: insight-users at itk.org
Subject: Re: [Insight-users] Problem implementing GradientAnisotropicDiffusionImageFilter
Can you provide a png file or meta file, or whatever file you tried with itk?
On Tue, Dec 1, 2009 at 4:08 PM, O'Connor, Michael
<Michael.OConnor at umassmed.edu> wrote:
> Hi Bill,
> Thanks. Attached is simple MATLAB test code + the MATLAB filter
> function used (Dr. Guy Gilboa's code). I've also take the liberty of
> attaching a test image (specifically opened in MATLAB code). The test
> image is 400x400 float with background = 0.024 + additive Gaussian Noise and
> foreground (circle) = 0.031 + additive Gaussian Noise. As said previously
> MATLAB code seems to work as expected as shown in plot. In ITK there is no
> similar intraregion filtering (i.e. filtered image ~ input )
>
> Regards,
> Michael
>
> ******
> %%% simple test of anisotropic diffusion filter
> %%% uses Guy Gilboa's diffusion function (from his website)
> %%%
> clear all;
> close all;
> close all hidden;
> row = 400;
> col = 400;
> GAIN = 1;
> iter = 5;
> conductance = 0.003 * GAIN;
> timestep = 0.125;
>
>
> %%
> fid_raw = fopen('/Users/michaeloconnor/Desktop/phantom_ADF3.rec','r');
>
> if (fid_raw == -1 )
> fprintf('NO Input File(s) in Working Directory \n');
> break; %stop now - invalid file
> end
> %%
>
> fprintf('iterations = %i timestep = %4.3f conductance = %f Gain = %i
> \n',iter,timestep,conductance, GAIN)
>
> read_raw= fread(fid_raw, [row col],'single');
> I=read_raw*GAIN;
>
> J_adf=diffusion2(I, 'pm1', iter, conductance, timestep);
>
> %% Display Both input and Filtered Using Image Processing Toolkit
>
> imtool(J_adf, [])
> imtool(I, [])
>
> %% Plot 1D Profile
>
> p=plot([1:400],I(200,:),[1:400],J_adf(200,:));
> set(p,'LineWidth',2);
> title('1D Profile (row 200)','FontSize',14, 'FontWeight', 'bold');
> xlabel('Column (1:400)','FontSize',14, 'FontWeight', 'bold')
> ylabel('MU Value','FontSize',14, 'FontWeight', 'bold')
> legend('Pre-Filter', 'ADF Filter')
> fclose all;
>
> ******* filter function (Dr. Gilboa's code)
> function Jd = diffusion2(J,method,N,K,dt,theta,sigma2)
> % private function: diffusion (by Guy Gilboa):
> % Jd=diffusion(J,method,N,K)
> % Simulates N iterations of diffusion, parameters:
> % J = source image (2D gray-level matrix) for diffusio
> % method = 'lin': Linear diffusion (constant c=1).
> % 'pm1': perona-malik, c=exp{-(|grad(J)|/K)^2} [PM90]
> % 'pm2': perona-malik, c=1/{1+(|grad(J)|/K)^2} [PM90]
> % 'rmp': complex valued - ramp preserving [GSZ01]
> % K edge threshold parameter
> % N number of iterations
> % dt time increment (0 < dt <= 0.25, default 0.2)
> % sigma2 - if present, calculates gradients of diffusion coefficient
> % convolved with gaussian of var sigma2 (Catte et al [CLMC92])
>
> if ~exist('N','var')
> %fprintf('No iterations\n')
> N=1;
> end
> if ~exist('K','var')
> % fprintf('No edge threshold\n')
> K=1;
> end
> if ~exist('dt','var')
> % fprintf('No time increment\n')
> dt=0.2;
> end
> if ~exist('theta','var')
> % fprintf('No theta\n')
> theta=pi/30;
> end
>
> if ~exist('sigma2','var')
> %fprintf('No gradients of diffusion\n')
> sigma2=0;
> end
>
> [Ny,Nx]=size(J);
>
> if (nargin<3)
> error('not enough arguments (at least 3 should be given)');
> end
>
> for i=1:N;
> % gaussian filter with kernel 5x5 (Catte et al)
> if (sigma2>0)
> Jo = J; % save J original
> J=gauss(J,7,sigma2);
> end
>
> % calculate gradient in all directions (N,S,E,W)
> In=[J(1,:); J(1:Ny-1,:)]-J;
> Is=[J(2:Ny,:); J(Ny,:)]-J;
> Ie=[J(:,2:Nx) J(:,Nx)]-J;
> Iw=[J(:,1) J(:,1:Nx-1)]-J;
>
> % calculate diffusion coefficients in all directions according to
> method
> if (method == 'lin')
> Cn=K; Cs=K; Ce=K; Cw=K;
> elseif (method == 'pm1')
> Cn=exp(-(abs(In)/K).^2);
> Cs=exp(-(abs(Is)/K).^2);
> Ce=exp(-(abs(Ie)/K).^2);
> Cw=exp(-(abs(Iw)/K).^2);
> elseif (method == 'pm2')
> Cn=1./(1+(abs(In)/K).^2);
> Cs=1./(1+(abs(Is)/K).^2);
> Ce=1./(1+(abs(Ie)/K).^2);
> Cw=1./(1+(abs(Iw)/K).^2);
> elseif (method == 'rmp') % complex - ramp preserving
> k=K; j=sqrt(-1);
>
> Cn=exp(j*theta)./(1+(imag(In)/(k*theta)).^2);
> Cs=exp(j*theta)./(1+(imag(Is)/(k*theta)).^2);
> Ce=exp(j*theta)./(1+(imag(Ie)/(k*theta)).^2);
> Cw=exp(j*theta)./(1+(imag(Iw)/(k*theta)).^2);
>
> else
> error(['Unknown method "' method '"']);
> end
>
> if (sigma2>0) % calculate real gradiants (not smoothed)
> In=[Jo(1,:); Jo(1:Ny-1,:)]-Jo;
> Is=[Jo(2:Ny,:); Jo(Ny,:)]-Jo;
> Ie=[Jo(:,2:Nx) Jo(:,Nx)]-Jo;
> Iw=[Jo(:,1) Jo(:,1:Nx-1)]-Jo;
> J=Jo;
> end
>
> % Next Image J
> J=J+dt*(Cn.*In + Cs.*Is + Ce.*Ie + Cw.*Iw);
>
> end; % for i
>
> Jd = J;
>
>
> %%%% Refs:
> % [PM90] P. Perona, J. Malik, "Scale-space and edge detection using
> anisotropic diffusion", PAMI 12(7), pp. 629-639, 1990.
> % [CLMC92] F. Catte, P. L. Lions, J. M. Morel and T. Coll, "Image selective
> smoothing and edge detection by nonlinear diffusion", SIAM J. Num. Anal.,
> vol. 29, no. 1, pp. 182-193, 1992.
> % [GSZ01] G. Gilboa, N. Sochen, Y. Y. Zeevi, "Complex Diffusion Processes
> for Image Filtering", Scale-Space 2001, LNCS 2106, pp. 299-307,
> Springer-Verlag 2001.
>
>
>
> -----Original Message-----
> From: Bill Lorensen [mailto:bill.lorensen at gmail.com]
> Sent: Tue 12/1/2009 1:01 PM
> To: O'Connor, Michael
> Cc: insight-users at itk.org
> Subject: Re: [Insight-users] Problem implementing
> GradientAnisotropicDiffusionImageFilter
>
> 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
>>
>>
>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091201/bed4ee55/attachment-0001.htm>
More information about the Insight-users
mailing list