function out = FastNonLocalMeans3D( V, sigma, beta, rs, rc, ps, flag, block )
% A fast implementation of the non-local means based on distances in
% the features space. The full algorithm is discussed in detail in the
% following paper:
%
% A. Tristán-Vega, V. García-Pérez, S. Aja-Fernández, C.-F. Westin
% "Efficient and robust nonlocal means denoising of MR data based on
% salient features matching"
% Computer Methods and Programs in Biomedicine, vol. 105, pp. 131-144
% (2012)
%
% If you are willing to use this software for your research, please cite
% this work.
%
% NOTE: Some of the computational features described in the paper above
% cannot be exploited in the matlab implementation. If performance is an
% issue for you, we strongly encourage you use the C++/ITK implementation
% available at: http://www.nitrc.org/projects/unlmeans, for which both
% source code and pre-compiled executables can be downloaded.
%
% USAGE: out = FastNonLocalMeans( V, sigma [beta, rs, rc, ps, flag] )
%
% V: The input volume to be filtered (3D). - MANDATORY
% sigma: The noise powe in the input image. In the Gaussian case, this
% is the standard deviation of the Gaussian noise at each pixel.
% In the Rician case, it is the standard deviation of noise in
% the original, Gaussian distributed, real and imaginary parts of
% the signal (whose modulus is computed to get teh Rician
% variable). - MANDATORY
% beta: The filtering parameter. The larger its value, the more
% aggressive the filtering. The smaller its value, the better
% details are preserved. It should be in the range of 0.8 to 1.2
% for best performance (Default: 1.0).
% rs: A 3x1 vector with the search radii (Default: 2,2,2).
% rc: A 3x1 vector with the comparison radii (Default: 1,1,1).
% ps: The preselection threshold. All those pixels in the search
% window whose normalized distance to the center pixel is larger
% than this value are automatically removed from the weighted
% average (Default: 2.0).
% flag: Must be either 'gaussian' (the default) or 'rician'. In the
% latter case, the weighted average is performed over the squared
% pixels, and the filtered value is computed as
% sqrt(mu-2·sigma^2) so that the estimate becomes unbiased.
% block: This second flag tells the algorithm if the computation of the
% weights within the search window must be done with a loop (0,
% the default since it seems to be faster) or it must be done
% with vector operations (1).
%
% out: The filtered volume.
if( nargin<2 )
error('At least the input volume and the noise power must be provided');
end
if( nargin<3 )
beta = 1.0;
end
h = beta*sigma;
if( nargin<4 )
rs = [2;2;2];
else
if( length(rs)~=3 )
rs = rs(1).*ones(3,1);
end
end
if( nargin<5 )
rc = [1;1;1];
else
if( length(rc)~=3 )
rc = rc(1).*ones(3,1);
end
end
if( nargin<6 )
ps = 2.0;
end
if( nargin<7 )
FLAG = 0;
else
if( strcmpi('gaussian',flag) )
FLAG = 0;
elseif( strcmpi('rician',flag) )
FLAG = 1;
else
error(['Unknown filtering type: ',flag]);
end
end
if( nargin<8 )
block = 0;
end
% Compute the size of the image:
[Y,X,Z] = size(V);
% Compute the features map:
[mu,Gx,Gy,Gz,factors,hcorr] = ComputeLocalFeatures3D( V, rc );
% Compute the effective value of h as described in the paper:
h = hcorr*h;
% Initiallize the output:
out = zeros(Y,X,Z);
% Loop along the pixels:
for x=1:X
for y=1:Y
for z=1:Z
% We are filtering the pixel (x,y,z). First, create a
% neighborhood around this pixel checking for out-of-bound
% indices:
mx = max( x-rs(1), 1 );
MX = min( x+rs(1), X );
my = max( y-rs(2), 1 );
MY = min( y+rs(2), Y );
mz = max( z-rs(3), 1 );
MZ = min( z+rs(3), Z );
% Keep the center values:
mu0 = mu(y,x,z);
gx0 = Gx(y,x,z);
gy0 = Gy(y,x,z);
gz0 = Gz(y,x,z);
if( block==1 )
% VECTOR IMPLEMENTATION (SEEMS TO BE SLOWER):
% Get the valeus of the pixels in the whole search neihborhood:
vals = V(my:MY,mx:MX,mz:MZ);
% Get the mean values and gradients of the pixels in the whole
% search neighborhood:
mui = mu(my:MY,mx:MX,mz:MZ);
gxi = Gx(my:MY,mx:MX,mz:MZ);
gyi = Gy(my:MY,mx:MX,mz:MZ);
gzi = Gz(my:MY,mx:MX,mz:MZ);
% Compute the distances:
dists = (mui-mu0).*(mui-mu0) + ...
(gxi-gx0).*(gxi-gx0)*factors(1) + ...
(gyi-gy0).*(gyi-gy0)*factors(2) + ...
(gzi-gz0).*(gzi-gz0)*factors(3);
% Normalize the distances:
dists = dists./(h*h);
% Compute the weights:
wis = exp(-dists);
% Set to 0 the normalized distances above the threshold to
% execute pre-selection:
wis(dists>ps) = 0;
% Compute the normalization factor:
NORM = sum(wis(:));
% Filter the pixel; average the pixels or their squared values
% depending on the filtering type:
if( FLAG==0 ) % Gaussian
pixel = sum(wis(:).*vals(:));
else % Rician
pixel = sum(wis(:).*vals(:).*vals(:));
end
else
% LOOP IMPLEMENTATION (SEEMS TO BE FASTER):
pixel = 0.0;
NORM = 0.0;
for s=mx:MX
for t=my:MY
for u=mz:MZ
% Get the current features:
mui = mu(t,s,u);
gxi = Gx(t,s,u);
gyi = Gy(t,s,u);
gzi = Gz(t,s,u);
% Compute the distance and normalize:
dist = (mu0-mui)*(mu0-mui) + ...
(gx0-gxi)*(gx0-gxi)*factors(1) + ...
(gy0-gyi)*(gy0-gyi)*factors(2) + ...
(gz0-gzi)*(gz0-gzi)*factors(3);
dist = dist/(h*h);
% Compute the weight in case the distance is below
% the pre-selection threshold, otherwise set to 0:
if( dist