2

I'm working on an algorithm, which requires filtering of a 3D matrix (non-sparse, 512^3) to find edges. I only want to find edges in each slice, so I have been doing the following:

% 2D loop appaoch    
[x,y]=ndgrid(floor(-3*sigma):ceil(3*sigma),floor(-3*sigma):ceil(3*sigma));    
DGauss=-(x./(2*pi*sigma^4)).*exp(-(x.^2+y.^2)/(2*sigma^2));    
filteredVolume = zeros(size(vol))    
for n = 1:size(vol,3)      
    filteredVolume(:,:,n) = imfilter(vol(:,:,n),DGauss,'conv','symmetric');    
end

I also tried to do the same by calling imfilter on the entire volume:

% 3D matrix approach    
filteredVolume = imfilter(vol,DGauss,'conv','symmetric');

I compared the performance of both of these approaches, but the loop version is significantly faster (6.5 seconds to 20 seconds). Should this behavior be expected? If so, why?

klurie
  • 1,060
  • 8
  • 16

1 Answers1

2

The reason it takes longer with the 3D version is because imfilter decides that the filter is non-separable. The function imfilter>isSeparable says the following:

function separable = isSeparable(a, h)

% check for filter separability only if the kernel has at least
% 289 elements [17x17] (non-double input) or 49 [7x7] (double input),
% both the image and the filter kernel are two-dimensional and the
% kernel is not a row or column vector, nor does it contain any NaNs of Infs

Since the input image is not 2D, the function returns false and a 2D filtering operation is done instead of two sequential 1D filters, which are faster.

On a side note, imfilter does not benefit from the JIT compiler. All the time is spent in the compiled function images\private\imfilter_mex.

chappjc
  • 30,359
  • 6
  • 75
  • 132