[rows, cols, scls] = size(scalespace);
    scalespace_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
                % scalespace_radii contaings the radius associated to a given scale, i.e.: radii = scale * sqrt(2)
                mask_feat_radii = (distance <= scalespace_radii(ss));
                currscale       = scalespace(:,:,ss);
                responses       = currscale(mask_feat_radii);

                scalespace_feats{ss}(rr,cc,:) = [mean(responses), std(responses), median(responses)];
            end
            fprintf('\n');
        end
    end