matlab: optimizing code for computing statistics in multi-scale circular neighborhoods

Question

I've a stack of images (imgstack) over which I would like to compute some statistics (e.g. mean, std, median) in multi-scale circular neighborhoods. For each image on the stack (currscale), the size of the circular neighborhood to be applied is pre-stored in a vector imgstack_radii(ss). The size of the circular neighboorhoods change across the images on the stack. Some example values for the radius of the circular filters are 1.4142,1.6818,2.0000,2.3784.

The code below does the work, however, since the size of my stack is pretty big (1200x1200x25) the computational time is very expensive. I was wondering if there's a way to optimize/vectorize the code? Any suggestions are appreciated!

[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, i.e.: 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 the feedback of @Shai and @Jonas, the final code looks like below. 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 using filters, which should be significantly faster.
For each imagestack_radii, first create a circular mask:

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

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

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

• std: once you have the mean, you can compute 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 );