function diff_vol = fast_anisodiff3D(vol, num_iter, delta_t, kappa, option, voxel_spacing) %ANISODIFF2D Conventional anisotropic diffusion % DIFF_VOL = ANISODIFF3D(VOL, NUM_ITER, DELTA_T, KAPPA, OPTION, VOXEL_SPACING) perfoms % conventional anisotropic diffusion (Perona & Malik) upon a stack of gray scale images. % A 3D network structure of 26 neighboring nodes is considered for diffusion conduction. % % ARGUMENT DESCRIPTION: % VOL - gray scale volume data (MxNxP). % NUM_ITER - number of iterations. % DELTA_T - integration constant (0 <= delta_t <= 3/44). % Usually, due to numerical stability this % parameter is set to its maximum value. % KAPPA - gradient modulus threshold that controls the conduction. % OPTION - conduction coefficient functions proposed by Perona & Malik: % 1 - c(x,y,z,t) = exp(-(nablaI/kappa).^2), % privileges high-contrast edges over low-contrast ones. % 2 - c(x,y,z,t) = 1./(1 + (nablaI/kappa).^2), % privileges wide regions over smaller ones. % VOXEL_SPACING - 3x1 vector column with the x, y and z dimensions of % the voxel (milimeters). In particular, only cubic and % anisotropic voxels in the z-direction are considered. % When dealing with DICOM images, the voxel spacing % dimensions can be extracted using MATLAB's dicominfo(.). % % OUTPUT DESCRIPTION: % DIFF_VOL - (diffused) volume with the largest scale-space parameter. % % Example % ------------- % vol = randn(100,100,100); % num_iter = 4; % delta_t = 3/44; % kappa = 70; % option = 2; % voxel_spacing = ones(3,1); % diff_vol = anisodiff3D(vol, num_iter, delta_t, kappa, option, voxel_spacing); % figure, subplot 121, imshow(vol(:,:,50),[]), subplot 122, imshow(diff_vol(:,:,50),[]) % % See also anisodiff1D, anisodiff2D. % References: % P. Perona and J. Malik. % Scale-Space and Edge Detection Using Anisotropic Diffusion. % IEEE Transactions on Pattern Analysis and Machine Intelligence, % 12(7):629-639, July 1990. % % G. Grieg, O. Kubler, R. Kikinis, and F. A. Jolesz. % Nonlinear Anisotropic Filtering of MRI Data. % IEEE Transactions on Medical Imaging, % 11(2):221-232, June 1992. % % MATLAB implementation based on Peter Kovesi's anisodiff(.): % P. D. Kovesi. MATLAB and Octave Functions for Computer Vision and Image Processing. % School of Computer Science & Software Engineering, % The University of Western Australia. Available from: % . % % Credits: % Daniel Simoes Lopes % ICIST % Instituto Superior Tecnico - Universidade Tecnica de Lisboa % danlopes (at) civil ist utl pt % http://www.civil.ist.utl.pt/~danlopes % % May 2007 original version. % Convert input volume to double. vol = double(vol); % Useful variables. [rows cols pags] = size(vol); % PDE (partial differential equation) initial condition. diff_vol = vol; clear vol % Center voxel distances. x = voxel_spacing(1); y = voxel_spacing(2); z = voxel_spacing(3); dx = 1; dy = 1; dz = 1; % 3D convolution masks - finite differences. h1 = zeros(3,3,3); h1(2,2,2) = -1; h1(2,2,1) = 1; h2 = zeros(3,3,3); h2(2,2,2) = -1; h2(2,2,3) = 1; h3 = zeros(3,3,3); h3(2,2,2) = -1; h3(2,1,2) = 1; h4 = zeros(3,3,3); h4(2,2,2) = -1; h4(2,3,2) = 1; h5 = zeros(3,3,3); h5(2,2,2) = -1; h5(3,2,2) = 1; h6 = zeros(3,3,3); h6(2,2,2) = -1; h6(1,2,2) = 1; % Anisotropic diffusion. for t = 1:num_iter % Finite differences. [imfilter(.,.,'conv') can be replaced by convn(.,.,'same')] % Due to possible memory limitations, the diffusion % will be calculated at each page/slice of the volume. for p = 1:pags-2 diff3pp = diff_vol(:,:,p:p+2); aux = imfilter(diff3pp,h1,'conv'); nabla1 = aux(:,:,2); aux = imfilter(diff3pp,h2,'conv'); nabla2 = aux(:,:,2); aux = imfilter(diff3pp,h3,'conv'); nabla3 = aux(:,:,2); aux = imfilter(diff3pp,h4,'conv'); nabla4 = aux(:,:,2); aux = imfilter(diff3pp,h5,'conv'); nabla5 = aux(:,:,2); aux = imfilter(diff3pp,h6,'conv'); nabla6 = aux(:,:,2); % Diffusion function. if option == 1 c1 = exp(-(nabla1/kappa).^2); c2 = exp(-(nabla2/kappa).^2); c3 = exp(-(nabla3/kappa).^2); c4 = exp(-(nabla4/kappa).^2); c5 = exp(-(nabla5/kappa).^2); c6 = exp(-(nabla6/kappa).^2); elseif option == 2 c1 = 1./(1 + (nabla1/kappa).^2); c2 = 1./(1 + (nabla2/kappa).^2); c3 = 1./(1 + (nabla3/kappa).^2); c4 = 1./(1 + (nabla4/kappa).^2); c5 = 1./(1 + (nabla5/kappa).^2); c6 = 1./(1 + (nabla6/kappa).^2); end % Discrete PDE solution. diff_vol(:,:,p+1) = diff_vol(:,:,p+1) + ... delta_t*(... (1/(dz^2))*c1.*nabla1 + (1/(dz^2))*c2.*nabla2 + ... (1/(dx^2))*c3.*nabla3 + (1/(dx^2))*c4.*nabla4 + ... (1/(dy^2))*c5.*nabla5 + (1/(dy^2))*c6.*nabla6); end % Iteration warning. fprintf('\rIteration %d\n',t); end % Search at MATLAB Central's FileExchange for the beepbeep(.) function % ("beepbeep" as keyword). It is a useful application to alert the termination of % a program execution by emitting a ritmic sequence of OS beeps. % http://www.mathworks.com/matlabcentral/fileexchange/loadCategory.do