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:
%   <http://www.csse.uwa.edu.au/~pk/research/matlabfns/>.
% 
% 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