


SMOOTH3F Smooth 3D data (fast version). W = SMOOTH3F(V) smoothes input data V with a Gaussian kernel. The smoothed data is returned in W. W = SMOOTH3F(V, 'filter') Filter can be 'gaussian' or 'box' (default) and determines the convolution kernel. W = SMOOTH3F(V, 'filter', SIZE) sets the size of the convolution kernel (default is [3 3 3]). If SIZE is a scalar, the size is interpreted as [SIZE SIZE SIZE]. Each element of SIZE is required to be an odd integer. W = SMOOTH3F(V, 'filter', SIZE, ARG) sets an attribute of the convolution kernel. When filter is 'gaussian', ARG is the standard deviation (default is .65). If filter is 'box', ARG has no effect. SMOOTH3F is similar to Matlab's built-in SMOOTH3 but uses a more efficient algorithm. (The only difference in calling the two functions is that SMOOTH3F requires an odd SIZE argument). See also SMOOTH3.



0001 function smoothed = smooth3f(data, filt, sz, arg) 0002 %SMOOTH3F Smooth 3D data (fast version). 0003 % W = SMOOTH3F(V) smoothes input data V with a Gaussian kernel. The 0004 % smoothed data is returned in W. 0005 % 0006 % W = SMOOTH3F(V, 'filter') Filter can be 'gaussian' or 'box' (default) 0007 % and determines the convolution kernel. 0008 % 0009 % W = SMOOTH3F(V, 'filter', SIZE) sets the size of the convolution 0010 % kernel (default is [3 3 3]). If SIZE is a scalar, the size is 0011 % interpreted as [SIZE SIZE SIZE]. Each element of SIZE is required 0012 % to be an odd integer. 0013 % 0014 % W = SMOOTH3F(V, 'filter', SIZE, ARG) sets an attribute of the 0015 % convolution kernel. When filter is 'gaussian', ARG is the standard 0016 % deviation (default is .65). If filter is 'box', ARG has no effect. 0017 % 0018 % SMOOTH3F is similar to Matlab's built-in SMOOTH3 but uses a more 0019 % efficient algorithm. (The only difference in calling the two 0020 % functions is that SMOOTH3F requires an odd SIZE argument). 0021 % 0022 % See also SMOOTH3. 0023 0024 % Modified from TMW's SMOOTH3. 0025 % The following commented code tests SMOOTH3f against SMOOTH3: 0026 % data = randn([50,50,50]); 0027 % tic; orig = smooth3(data); t(1) = toc; 0028 % tic; modf = smooth3f(data); t(2) = toc; 0029 % mse = sqrt(mean(abs(orig(:)-modf(:)).^2)); % mean squared error 0030 % printf('SMOOTH3: %4.3f sec SMOOTH3F: %4.3f sec MSE: %5.3f', t(1), t(2), mse); 0031 0032 %%%%%%%%%%%%%%%%%%%%%%%%%%% Parse Inputs %%%%%%%%%%%%%%%%%%%%%%%%%% 0033 if (nargin==1) %smooth3(data) 0034 arg = .65; sz = 3; filt = 'b'; 0035 elseif (nargin==2) %smooth3(data, filter) 0036 arg = .65; sz = 3; 0037 elseif (nargin==3) %smooth3(data, filter, sz) 0038 arg = .65; 0039 elseif ((nargin>4) | (nargin==0)) 0040 error('Wrong number of input arguments.'); 0041 end 0042 0043 if (ndims(data)~=3), error('V must be a 3D array.'); end; 0044 0045 if (numel(sz)==1) 0046 sz = [sz sz sz]; 0047 elseif (prod(size(sz))~=3) 0048 error('SIZE must be a scalar or a 3 element vector.') 0049 end 0050 sz = sz(:)'; % force column vector 0051 0052 szHalf = (sz-1)/2; 0053 if (~isequal(szHalf,round(szHalf))), error('SIZE must contain only odd integers.'); end; 0054 0055 %%%%%%%%%%%%%%%%%%%%%%%%%% Make the kernel %%%%%%%%%%%%%%%%%%%%%%%% 0056 % Make three kernels so that the full convolution kernel is the 0057 % outer product of the three ... 0058 switch(lower(filt(1))), 0059 case 'g', % 3gaussian 0060 kernel{1} = gausskernel(szHalf(1),arg); 0061 kernel{2} = gausskernel(szHalf(2),arg); 0062 kernel{3} = gausskernel(szHalf(3),arg); 0063 case 'b', % box 0064 kernel{1} = ones(sz(1),1)./sz(1); 0065 kernel{2} = ones(sz(2),1)./sz(2); 0066 kernel{3} = ones(sz(3),1)./sz(3); 0067 otherwise, 0068 error('Unknown filter.'); 0069 end 0070 0071 %%%%%%%%%%%%%%%%%%%%%%%%%%% Do the Smooth %%%%%%%%%%%%%%%%%%%%%%%%% 0072 % Its grossly inefficient to do the full convolution since the kernel 0073 % is separable. We use CONVNSEP to do three 1-D convolutions. 0074 smoothed = convnsep(kernel{:}, padreplicate(data,(sz-1)/2), 'valid'); 0075 0076 0077 %%%%% TAKEN FROM TMW's SMOOTH3 rev 1.7 -- pads an array by replicating values. 0078 function b=padreplicate(a, padSize) 0079 numDims = length(padSize); 0080 idx = cell(numDims,1); 0081 for k = 1:numDims 0082 M = size(a,k); 0083 onesVector = ones(1,padSize(k)); 0084 idx{k} = [onesVector 1:M M*onesVector]; 0085 end 0086 b = a(idx{:});