function pm = pm_smooth_phasemap(pm,angvar,vxs,fwhm)
% Performs a weighted (by 1/angvar) gaussian smoothing of a phasemap
% FORMAT pm = pm_smooth_phasemap(pm,angvar,vxs,fwhm)
%
% Input:
% pm         : Phase-map
% angvar     : Map of uncertainty of the angular estimate.
% vxs        : Voxel sizes (mm) in the three directions.
% fwhm       : FWHM (mm) of gaussian kernel for the three
%              directions (or scalar for isotropic kernel).
%__________________________________________________________________________

% Jesper Andersson
% Copyright (C) 2003-2022 Wellcome Centre for Human Neuroimaging

if nargin ~= 4 || nargout ~= 1
   help pm_smooth_phasemap
   pm = [];
   return;
end
if length(fwhm) == 1
   fwhm = fwhm * [1 1 1];
end

%
% Determine kernel. The code has been knicked
% (and changed a bit) from spm_smooth.
%
fwhm  = fwhm./vxs;                  % voxel anisotropy
fwhm  = max(fwhm,ones(size(fwhm))); % lower bound on FWHM
s     = fwhm/sqrt(8*log(2));        % FWHM -> Gaussian parameter

x  = round(3*s(1)); x = [-x:x];
y  = round(3*s(2)); y = [-y:y];
z  = round(3*s(3)); z = [-z:z];
x  = exp(-(x).^2/(2*(s(1)).^2));
y  = exp(-(y).^2/(2*(s(2)).^2));
z  = exp(-(z).^2/(2*(s(3)).^2));
x  = x/sum(x);
y  = y/sum(y);
z  = z/sum(z);

i  = (length(x) - 1)/2;
j  = (length(y) - 1)/2;
k  = (length(z) - 1)/2;

q1 = zeros(size(pm));
q2 = zeros(size(pm));
spm_conv_vol(pm./angvar,q1,x,y,z,-[i,j,k]);
spm_conv_vol( 1./angvar,q2,x,y,z,-[i,j,k]);
pm = q1./q2;
