[Insight-users] Problem implementing GradientAnisotropicDiffusionImageFilter
Bill Lorensen
bill.lorensen at gmail.com
Tue Dec 1 23:39:26 EST 2009
It's a raw file. Can you provide the config file that describes the image.
On Tue, Dec 1, 2009 at 4:17 PM, Bill Lorensen <bill.lorensen at gmail.com> wrote:
> 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
>>>
>>>
>>
>>
>>
>>
>
More information about the Insight-users
mailing list