从一千幅图像中提取许多感兴趣的区域 ROI）

Question

我有大量的显微图像，每张图像都有数百个点 (ROI)。这些点在空间中是固定的。我想从每个图像中提取每个点并保存到工作区中，以便进一步分析它们。

我自己编写了一个代码，它运行良好，但速度太慢。从每张图像中完全读出所有点大约需要 250 秒。

我的代码核心如下：

for s=1:NumberImages   
  im1=imread(fn(s,1).name);    
  im=im1-medfilt2(im1,[15,15]);    
  for i=1:length(p_g_x)    
    GreenROI(i,s)=double(sum(sum(im(round(p_g_y(i))+(-2:2), round(p_g_x(i))+(-2:2)))));
    RedROI(i,s)=double(sum(sum(im(round(p_r_y(i))+(-2:2), round(p_r_x(i))+(-2:2)))));        
  end
end

从代码中可以看出，我正在提取 5x5 区域。 p_g_x 的长度在 500-700 之间。

感谢您的意见。我使用配置文件查看器来确定究竟哪个功能需要更多时间。它是中值滤波器，需要花费大量时间（~90%）。

任何加快速度的建议将不胜感激。

谢谢

马希帕尔

Answer 1

使用 Matlab 的分析工具！

profile on % Starts the profiler
% Run some code now.
profile viewer % Shows you how often each function was called, and
               % where most time was spent. Try to start with the slowest part.
profile off  % Resets the Profiler, so you can measure again.

预分配

预分配输出，因为您知道大小，这样会更快。 （Matlab 已经告诉过你了！）

GreenROI = zeros(length(p_g_x), NumberImages); % And the same for RedROI.

使用卷积

阅读 Matlab 的 conv2 代码。

for s=1:NumberImages   
  im=imread(fn(s,1).name);    
  im=im-medfilt2(im,[15,15]);    
  % Pre-compute the sums first. This will only be faster for large p_g_x
  roi_image = conv2(im, ones(5,5));
  for i=1:length(p_g_x)    
    GreenROI(i,s)=roi_image(round(p_g_y(i)), round(p_g_x(i))); % You might have to offset the indices by 2, because of the convolution. Check that results are the same.
    RedROI(i,s)=roi_image(round(p_r_y(i)), round(p_r_x(i)));        
  end
end

Matlab 化代码

现在，您已经使用卷积获得了 5x5 窗口上的总和图像（或者您可以使用 @Shai 的 accumarray，同样的事情），您可以通过不遍历中的每个元素来进一步加快速度p_g_x 但直接将其用作向量。

我把它留给读者作为练习。 （使用sub2ind 将p_g_x 和p_g_y 转换为索引，作为提示）。

更新

我们的答案（包括我的答案）表明过早优化是一件坏事。在不知道的情况下，我假设您的循环将花费大部分时间，但是当您测量它时（谢谢！）事实证明这不是问题所在。瓶颈是medfilt2 中值滤波器，它占用了 90% 的时间。所以你应该先解决这个问题。 （请注意，在我的计算机上，您的原始代码对我的口味来说已经足够快了，但它仍然是中值滤波器占用的大部分时间。）

查看中值滤波操作的作用可能有助于我们弄清楚如何使其更快。这是一张图片。在左侧，您可以看到原始图像。中间是中值滤波器，右边是结果。

在我看来，结果与边缘检测结果非常相似。 （数学上这并不奇怪。）

我建议您开始尝试各种边缘检测。看看 Canny 和 Sobel。或者只使用conv2(image, kernel_x) where kernel_x = [1, 2, 1; 0, 0, 0; -1, -2, -1] 和相同但转换为kernel_y。您可以在此处找到各种边缘检测选项：edge(im, option)。我尝试了{'sobel', 'canny', 'roberts', 'prewitt'} 的所有选项。除了 Canny，它们都与您的中值滤波方法花费的时间大致相同。 Canny 慢了 4 倍，其余的（包括原版）需要 7.x 秒。所有这一切都没有 GPU。 imgradient 是 9 秒。

因此，我会说你不能再快了。如果你有一个 GPU 并且它可以与 Matlab 一起工作，你可以加快它的速度。将您的图像加载为gpuArrays。 medfilt2 documentation 上有一个例子。你仍然可以做一些小的加速，但它们只能提高 10% 的速度，所以几乎不值得。

Answer 2

你应该做的几件事

1. Pre-allocate 由Didac Perez 建议。

2. 使用profiler 来查看究竟在您的代码中花费了多少时间，是中值过滤器吗？是索引吗？

3. 假设所有图像的大小相同，您可以使用accumarray 和固定掩码subs 快速求和：

   subs_g = zeros( h, w ); %// allocate mask for green
   subs_r = zeros( h, w );  
   subs_g( sub2ind( [h w], round(p_g_y), round(p_g_x) ) = 1:numel(p_g_x); %//index each region
   subs_g = conv2( subs_g, ones(5), 'same' );
   subs_r( sub2ind( [h w], round(p_r_y), round(p_r_x) ) = 1:numel(p_r_x); %//index each region
   subs_r = conv2( subs_r, ones(5), 'same' );
   sel_g = subs_g > 0;
   sel_r = subs_r > 0;
   subs_g = subs_g(sel_g);
   subs_r = subs_r(sel_r);
   

   一旦这些蒙版固定，您就可以处理所有图像

   %// pre-allocation goes here - I'll leave it to you
   for s=1:NumberImages   
       im1=imread(fn(s,1).name);    
       im=double( im1-medfilt2(im1,[15,15]) ); 
       accumarray( subs_g, im( sel_g ) ); % summing all the green ROIs 
       accumarray( subs_r, im( sel_r ) ); % summing all the green ROIs 
   end
   

Answer 3

首先，预先分配您的 GreenROI 和 RedROI 结构，因为您之前知道最终大小。现在，您在每次迭代中一次又一次地调整它们的大小。

其次，我建议你使用“tic”和“toc”来调查问题出在哪里，它会给你有用的时机。

Answer 4

对每个图像进行操作的矢量化代码 -

%// Pre-compute green and red indices to be used across all the images
r1 = round(bsxfun(@plus,permute(p_g_y,[3 2 1]),[-2:2]'));
c1 = round(bsxfun(@plus,permute(p_g_x,[3 2 1]),[-2:2]));
green_ind = reshape(bsxfun(@plus,(c1-1)*size(im,1),r1),[],numel(p_g_x));

r2 = round(bsxfun(@plus,permute(p_r_y,[3 2 1]),[-2:2]'));
c2 = round(bsxfun(@plus,permute(p_r_x,[3 2 1]),[-2:2]));
red_ind = reshape(bsxfun(@plus,(c2-1)*size(im,1),r2),[],numel(p_g_x));

for s=1:NumberImages
    im1=imread(fn(s,1).name);
    im=double(im1-medfilt2(im1,[15,15]));

    GreenROI=sum(im(green_ind));
    RedROI =sum(im(red_ind));
end