Matlab: code optimization for calculating statistics in multiscale circular neighborhoods

Question

I have an image stack ( imgstack ) that I’d like to calculate some statistics ( eg mean, std, median ) in multiscale circular surroundings. For each image in the stack ( currscale ), the size of the applied circular neighborhood is previously stored in the vector imgstack_radii(ss) . The size of the circular neighborhoods changes according to the images on the stack. Some exemplary values ​​for the radius of the circular filters are 1.4142,1.6818,2.0000,2.3784 .

The code below does the job, since the size of my stack is quite large ( 1200x1200x25 ), the calculation time is very expensive. I was wondering if there is a way to optimize / vectorize the code? Any suggestions are welcome!

 [rows, cols, scls] = size(imgstack); imgstack_feats = cell(scls,1); [XX, YY] = meshgrid(1:cols, 1:rows); for rr=1:rows for cc=1:cols distance = sqrt( (YY-rr).^2 + (XX-cc).^2  ); for ss=1:scls % imgstack_radii contains the radius associated to a given scale, ie: radii = scale * sqrt(2) mask_feat_radii = (distance <= imgstack_radii(ss)); currscale    = imgstack(:,:,ss); responses    = currscale(mask_feat_radii); imgstack_feats{ss}(rr,cc,:) = [mean(responses), std(responses), median(responses)]; end end end 

After feedback from @Shai and @Jonas, the final code is as follows. Thanks guys!

 function disk = diskfilter(radius) height = 2*ceil(radius)+1; width = 2*ceil(radius)+1; [XX,YY] = meshgrid(1:width,1:height); dist = sqrt((XX-ceil(width/2)).^2+(YY-ceil(height/2)).^2); circfilter = strel(dist <= radius); end for ss=1:scls fprintf('\tProcessing scale: %d radii: %1.4f\n', ss, imgstack_radii(ss)); disk = diskfilter(imgstack_radii(ss)); neigh = getnhood( disk ); imgstack_feats{ss}(:,:,1) = imfilter( imgstack(:,:,ss), double(neigh)/sum(neigh(:)), 'symmetric' ); tmp = imfilter( imgstack(:,:,ss).^2, double(neigh)/sum(neigh(:)), 'symmetric' ); imgstack_feats{ss}(:,:,2) = sqrt( tmp - imgstack_feats{ss}(:,:,1) ); imgstack_feats{ss}(:,:,3) = ordfilt2( imgstack(:,:,ss), ceil( sum(neigh(:))/2 ), neigh ); end 

Answer 1

You can replace all operations with filters, which should be much faster.
For each imagestack_radii first create a circular mask:

 n = getnhood( strel('disk', imagestack_radii(s), 0 ) ); 

• means: use imfilter with double(n)/sum(n(:)) as a filter

   imgstack_feats{ss}(:,:,1) = imfilter( imgstack(:,:,ss), double(n)/sum(n(:)), 'symmetric' ); 

• std: as soon as you have the average value, you can calculate the second moment using

   tmp = imfilter( imgstack(:,:,ss).^2, double(n)/sum(n(:)), 'symmetric' ); imgstack_feats{ss}(:,:,2) = sqrt( tmp - imgstack_feats{ss}(:,:,1) ); 

• median: use ordfilt2

   imgstack_feats{ss}(:,:,3) = ordfilt2( imgstack(:,:,ss), ceil( sum(n(:))/2 ), n );