计算图像中已识别像素的局部均值和方差

Question

我想计算由窗口（虚线）创建的子区域的均值和标准微分，并以识别的像素红色为中心（称为局部均值和标准微分）。这是描述它的图

我们可以通过使用蒙版卷积图像来做到这一点。但是，这需要很长时间，因为我只关心服务器点的均值和标准推导，而卷积计算的是图像中的整个点。您能否有一种更快的方法来解决它，只计算已识别像素处的均值和标准导数？我是用matlab做的。这是我的卷积函数代码

I=[18    36    70    33    64    40    62    76    71    37     5
   82    49    86    45    96    29    74     7    60    56    45
   25    32    55    48    25    30    12    82    95    77     8
   24    18    78    74    19    57    67    59    16    46    78
   28     9    59     2    29    11     7    31    75    15    25
   83    26    96     8    82    26    85    12    11    28    19
   81    64    78    70    26    33    17    72    81    16    54
   75    39    78    34    59    31    77    31    61    81    89
   89    84    29    99    79    25    26    35    65    56    76
   93    90    45     7    61    13    34    24    11    34    92
   88    82    91    81   100     4    88    70    85     8    19];
identified_position=[30 36 84 90] %indices of pixel 78, 48,72 60
mask=1/9.*ones(3,3);
mean_all=imfilter(I,mask,'same');
%Mean of identified pixels
mean_all(identified_position)
% Compute the variance
std_all=stdfilt(I,ones(3));
%std of identified pixels
std_all(identified_position)

这是对比代码

function compare_mean(dimx,dimy)
I=randi(100,[dimx,dimy]);
rad=3;
identified_position=randi(max(I(:)),[1,5]);% Get 5 random position
function way1()

    mask=ones(rad,rad);
    mask=mask./sum(mask(:));
    mean_all=conv2(I,mask,'same');
    mean_out =mean_all(identified_position);
end
function way2()
    box_size = rad; %// Edit your window size here (an odd number is preferred)
    bxr = floor(box_size/2); %// box radius
    %// Get neighboring indices and those elements for all identified positions
    off1 = bsxfun(@plus,[-bxr:bxr]',[-bxr:bxr]*size(I,1)); %//'#neighborhood offsets
    idx = bsxfun(@plus,off1(:),identified_position); %// all absolute offsets
    I_selected_neigh = I(idx); %// all offsetted elements
    mean_out = mean(I_selected_neigh,1); %// mean output

end
way2()
time_way1=@()way1();timeit(time_way1)
time_way2=@()way2();timeit(time_way2)
end

有时way2有错误

Subscript indices must either be real positive integers or logicals.

Error in compare_mean/way2 (line 18)
    I_selected_neigh = I(idx); %// all offsetted elements

Error in compare_mean (line 22)
way2()

Answer 1

讨论和解决方案代码

假设I 作为输入图像，identified_position 作为选定点的线性索引，bxsz 作为窗口/框大小，接下来列出的方法一定非常有效 -

%// Get XY coordinates
[X,Y] =  ind2sub(size(I),identified_position);
pts = [X(:) Y(:)];

%// Parameters
bxr = (bxsz-1)/2;
Isz = size(I);

%// XY coordinates of neighboring elements
[offx,offy] = ndgrid(-bxr:bxr,-bxr:bxr);
x_idx = bsxfun(@plus,offx(:),pts(:,1)');  %//'
y_idx = bsxfun(@plus,offy(:),pts(:,2)');  %//'

%// Outside image boundary elements
invalids = x_idx>Isz(1) | x_idx<1 | y_idx>Isz(2) | y_idx<1;

%// All neighboring indices
all_idx = (y_idx-1)*size(I,1) + x_idx;
all_idx(invalids) = 1;

%// All neighboring elements
all_vals = I(all_idx);
all_vals(invalids) = 0;

mean_out = mean(all_vals,1); %// final mean output

stdfilts = stdfilt(all_vals,ones(bxsz^2,1))
std_out = stdfilts(ceil(size(stdfilts,1)/2),:)  %// final stdfilt output

基本上，它使用bsxfun 一次性获取所有已识别位置的所有相邻索引，从而获得所有这些相邻元素。然后使用这些选定的元素来获取 mean 和 stdfilt 输出。整个想法是 keep the memory requirement minimum 并同时在这些选定的元素中执行 everything in a vectorized fashion。希望这必须更快！

---

基准测试

基准代码

dx = 10000; %// x-dimension of input image
dy = 10000; %// y-dimension of input image
npts = 1000; %// number of points
I=randi(100,[dx,dy]); %// create input image of random intensities
identified_position=randi(max(I(:)),[1,npts]);
rad=5; %// blocksize (rad x rad)

%// Run the approaches fed with the inputs
func1 = @() way1(I,identified_position,rad); %// original approach 
time1 = timeit(func1);
clear func1

func2 = @() way2(I,identified_position,rad); %// proposed approach 
time2 = timeit(func2);
clear func2

disp(['Input size: ' num2str(dx) 'x' num2str(dy) ' & Points: ' num2str(npts)])
disp(['With Original Approach: Elapsed Time = ' num2str(time1) '(s)'])
disp(['With Proposed Approach: Elapsed Time = ' num2str(time2) '(s)'])
disp(['**Speedup w/ Proposed Approach : ' num2str(time1/time2) 'x!**'])

相关功能码

%// OP's stated approach
function mean_out = way1(I,identified_position,rad)
mask=ones(rad,rad);
mask=mask./sum(mask(:));
mean_all=conv2(I,mask,'same');
mean_out =mean_all(identified_position);
return;

function mean_out = way2(I,identified_position,rad)
%//.... code from proposed approach stated earlier until mean_out %//

运行时结果

Input size: 10000x10000 & Points: 1000
With Original Approach: Elapsed Time = 0.46394(s)
With Proposed Approach: Elapsed Time = 0.00049403(s)
**Speedup w/ Proposed Approach : 939.0778x!**