用于循环的 Matlab 矢量化以操作 RGB 图像中的像素

Question

我仍在学习 matlab，并且正在尝试了解矢量化。我认为我的问题的根源是我不明白如何引用不同的矩阵等。我希望这个问题的答案能帮助我理解。

FI = imread(ForegroundImageName);
BI = imread(BackgroundImageName);
refRows =size(FI,1);
refCols =size(FI,2);
refChan =size(FI,3);
CommonRGB = mode(mode(FI));
BI = imresize(BI, [refRows refCols]);
swappedPixels = 0;

for row=1:refRows
    for col=1:refCols
        if(FI(row,col,:)==CommonRGB)
            FI(row,col,:)=BI(row,col,:);
            swappedPixels = swappedPixels + 1;
        end
    end
end

这个问题的背景是，如果前景像素与最常见的颜色匹配，我将用背景像素替换前景像素。 CommonRGB 是一个 1x1x3 矩阵，因为 mode(mode(FI)) 会吐出它。这些图像是 3D 彩色 RGB 图像。我选择这个作为我的问题的一个简单例子。这个 for 循环执行我想要它做的事情并且似乎工作。只是运行for循环需要很长时间。我执行mode(mode(FI)) 是运行上述for 循环所需时间的一小部分，我尝试自己实现mode(mode(FI))，与上面的像素交换相比，它变得相当复杂。创建直方图时我遇到了类似的问题。我希望你能帮助我更多地了解 matlab 并将我的编程知识扩展到这门语言中。我知道如果我说我想对整个矩阵执行这些操作会更容易，所以如果我们对其进行矢量化，我会假设我们不需要 refRows、refCols 和 refChan 变量在这种情况下.

我的一次失败尝试

if(FI(:,:,:)==CommonRGB)
    FI(:,:,:)=BI(:,:,:);
    swappedPixels = swappedPixels + 1;
end

从目前的答案来看，他们已经展示了向量化逻辑的方法。我得到的结论是，任何其他不涉及长 for 循环的方法都将比像我所做的那样使用大型 for 循环快得多。即使这意味着创建额外的数组并执行额外的过程，例如创建掩码数组并多次运行图片。据我所知，问题的根源在于 matlab 的 JIT Just In Time 编译器必须在 for 循环的每次迭代中重新解析命令。这种对for循环的解析和处理才是速度问题的真正根源。如果 matlab 可以“看到” for 循环并提前计划，那么它会运行得更快。因此，我对原始代码无能为力，我只需将 row 和 col 替换为其他内容并删除两个 for 循环。我将不得不设计其他不使用大 for 循环的方法。只有这样它才能以合理的速度运行。这告诉我，matlab 是一种真正的前端脚本语言，因为使用 for 循环运行脚本同样会遭受不利的性能下降。因此，我不确定如何将此问题标记为已回答、未回答或投票支持其他答案。因为如果我改变 for 循环内的逻辑，那么我将不得不改变我实现加速的方式。

所以，除非有人能告诉我如何用 1:refRows 之类的其他东西替换 for 循环中的变量 row 和 col，使用 parfor，使用简单的东西，只更改 for 循环内的代码并删除 for循环，或者可以只写一个答案来验证上面的段落，那么我现在不知道该怎么做。

正如@Divakar 所说，“你说得对，这里没有明确的矢量化技术来插入循环参数并从中获得矢量化解决方案。如果你有一个循环而不是这两个嵌套循环在循环情况下，您可以使用逻辑索引获得更直接的即插即用矢量化解决方案。但是，是的，矢量化在大多数情况下无法推广，需要根据具体情况进行处理。祝你好运矢量化探索！- Divakar"

话虽如此，我想帮助其他人，向他们展示我做了什么来应用我从这个问题中学到的东西。似乎我已经弄清楚了我认为的算法，尽管我仍然不确定我使用光阈值选择常见色调的第 3 阶段是否有效。当光阈值升高时，我没有看到它选择不同的色调。我只是看到第 3 阶段始终选择与第 2 阶段不同的色调。直到今天我才意识到这一点。对问题本身并不重要，但下面对问题很重要，因为我提到了 for 循环方法，然后实现了矢量化方法。在第 2-3 阶段比在第 1 阶段更是如此，我不需要做太多的改变。 ...好吧，它并没有让我粘贴那么多信息，所以我会尝试另一种方法。

My PDF of my lab report with code and pictures 我知道 mode(mode(FI)) 不正确，但我仍然很好奇第 2 阶段如何获得与第 3 阶段不同的色调。我认为第 2 阶段是正确的，我一次逐个像素地比较以验证第 2 阶段是否相同作为第 3 阶段，但 mode(FI(:)) 应该是正确的。如果您发现任何问题，请告诉我。我理解了大部分矢量化，除了使用掩码作为索引的工作原理。我知道他们称之为逻辑索引，但我猜只有 1 表示它将使用它，而零表示它不会查看该值。

Answer 1

矢量化有两个非常重要的工具 - bsxfun 和 logical indexing。现在，bsxfun 确实可以与逻辑运算符一起使用，因此只需这两个函数，您就可以实现大部分矢量化代码。

图像处理中的bsxfun 和logical indexing

如何以及何时使用logical indexing？

在某些情况下，在处理 RGB 图像时，特别是在需要索引到二维时（就像这里的情况），它可能会变得很棘手，因为你不能在那里直接使用逻辑索引。

在这种情况下可以使用的技巧之一是使用reshape，前提是其余维度所需的索引是统一结构的，这在此问题案例中再次得到满足。

如何以及何时使用bsxfun？

bsxfun 在各种问题中找到最优解。对于给定的问题，我们可以将其用于以矢量化方式进行的许多相等比较检查。

这是通过将CommonRGB 重新整形为一个大小来完成的，这样在重新整形的FI 中存在一个与row 和col 的组合维度相对应的单一维度。然后，bsxfun 处理重整后的CommonRGB 急需的"expansion"，并使用内置的函数句柄@eq 执行相等性检查。

让我们跳到使用前面讨论过的这两个函数和技术的矢量化代码 -

%// Reshape the inputs and save as new variables
FIr = reshape(FI,[],3);
BIr = reshape(BI,[],3);
CommonRGBr = reshape(CommonRGB,[],3);

%// Create 1D equality mask corresponding the equality satisfied across 
%// two dimensions - `row` and `col`
mask = all(bsxfun(@eq,FIr,CommonRGBr),2);
swappedPixels = sum(mask); %// Get the count of "swappings"

%// Perform the "swaps" and then reshape FI back to its original size
FIr(mask,:) = BIr(mask,:);
FI = reshape(FIr,size(FI));

---

基准测试

本节将讨论针对原始代码对提议的矢量化代码进行基准测试。首先，我们创建了 原始代码和矢量化代码的函数形式。

(1) 原代码-

function [FI,swappedPixels] = org_code(FI,BI,CommonRGB)

refRows =size(FI,1);
refCols =size(FI,2);

swappedPixels = 0;
for row=1:refRows
    for col=1:refCols
        if(FI(row,col,:)==CommonRGB)
            FI(row,col,:)=BI(row,col,:);
            swappedPixels = swappedPixels + 1;
        end
    end
end

return;

(2) 向量化代码-

function [FI,swappedPixels] = vectorized_code(FI,BI,CommonRGB)

FIr = reshape(FI,[],3);
BIr = reshape(BI,[],3);
CommonRGBr = reshape(CommonRGB,[],3);

mask = all(bsxfun(@eq,FIr,CommonRGBr),2);
swappedPixels = sum(mask);

FIr(mask,:) = BIr(mask,:);
FI = reshape(FIr,size(FI));

return;

(3) 最后是标杆代码，它还绘制了加速比 -

N_arr = [100 200 500 1000 2000 4000]; %// datasize array
timeall = zeros(2,numel(N_arr));
for iter = 1:numel(N_arr)
    N = N_arr(iter);
    FI = uint8(randi(255,N,N,3));
    BI = uint8(randi(255,N,N,3));
    
    CommonRGB = mode(mode(FI));
    BI = imresize(BI, [size(FI,1) size(FI,2)]); 
    
    f = @() org_code(FI,BI,CommonRGB);
    timeall(1,iter) = timeit(f);
    clear f
    
    f = @() vectorized_code(FI,BI,CommonRGB);
    timeall(2,iter) = timeit(f);
    clear f
end
figure,hold on,grid on
plot(N_arr,timeall(1,:)./timeall(2,:),'-bo')
legend('Speedup with Vectorized Code over Original one'),
xlabel('Datasize ->'),ylabel('Speedup Factor (x)')

在我的系统上得到的绘图结果是这样的 -

结论

这表明建议的矢量化解决方案的加速比接近170x！希望这会在一定程度上吸引您进入矢量化业务，尤其是bsxfun。

---

系统配置

MATLAB Version: 8.3.0.532 (R2014a)
Operating System: Windows 7
RAM: 3GB
CPU Model: Intel® Pentium® Processor E5400 (2M Cache, 2.70 GHz)

Answer 2

您正确地推测逻辑索引是要走的路。您可以在数组中选择满足条件的像素并将结果用于索引：

在FI 和BI 的情况下，灰度图像：

FI(FI==value)=BI(FI==value);

---

这里发生了什么： FI==value 生成一个大小为 FI 的二进制矩阵，其中满足条件的位置为 1（真），其他位置为 0（假）。然后可以使用此二进制矩阵访问FI 或任何其他相同大小的矩阵中的相应元素。这里FI中满足条件的元素与BI中相同位置的元素互换。

---

同样的方法可以用于 RGB 图像，只要您创建一个 3D 数组，沿第 3 维复制二进制掩码，例如使用 repmat。

Mask=(FI(:,:,1)==CommonRGB(1,1,1)).*(FI(:,:,2)==CommonRGB(1,1,2)).*(FI(:,:,3)==CommonRGB(1,1,3));
MaskRGB=repmat(Mask, 1, 1, 3);

FI(MaskRGB)=BI(MaskRBG);

注意.* 是像素乘法，它与逻辑与相同。 See also this。因此，第一行创建了一个二进制掩码，其中 RGB 图像的所有 3 层都满足条件。

如果您需要交换像素的数量，请使用：

swappedPixels = nnz(Mask);

Answer 3

据我了解，这是一种可以做到的方式。

您可以使用find 代替 if 语句来定位 FI 的每个通道中与该特定通道的最常见值相似的像素（因此是 commonRGB 的第 3 维）。一旦有了这些像素的行和列，就可以轻松快速地使用逻辑索引将像素从 BI 交换到 FI。请注意，我使用了 Matlab 附带的测试图像，因此输出可能不代表任何内容，但我认为它可以满足您的需求。如果没有，请告诉我！

这里是注释代码：

clear
clc
close all

FI1 = imread('peppers.png');
BI = imread('pears.png');

FI = FI1; %// Only used to display original image

[refRows,refCols,refChan] = size(FI); %// 1-liner for assignments.

CommonRGB = squeeze(mode(mode(FI)));

BI = imresize(BI, [refRows refCols]);

swappedPixels = zeros(1,3); %// We will use it differently than in your code.

%// Loop through each channel and use find to locate pixels in FI that are
%// similar to commonRGB
for CheckChannel = 1:refChan

    [row,col] = find(FI(:,:,CheckChannel) == CommonRGB(CheckChannel));

    FI(row,col,CheckChannel) = BI(row,col,CheckChannel); %// Index pixel value from BI into FI.

    swappedPixels(CheckChannel) = numel(row); %// Calculate the number of elements in the vector, i.e. number of swapped pixels.
end

 swappedPixels %// To see the number of swapped pixels

%// Display results.
figure;

subplot(1,3,1)
imshow(FI1);
title('Original foreground image','FontSize',16);

subplot(1,3,2)
imshow(BI);
title('Original background image','FontSize',16);

subplot(1,3,3)
imshow(FI);
title('Modified foreground','FontSize',16);

数组swappedPixels 如下所示：

swappedPixels =

       11175        4508       10753

其中每个数字对应于图像的 3 个维度之一。

输出图像如下所示：

Answer 4

我会这样解决这个问题。

  refRows =size(fi,1);
  refCols =size(fi,2);
  CommonRGB = mode(mode(FI));
  % finding where the 3 channels are equals
  [indicesX indicesY ] = find(fi(:,:,1)==CommonRGB(1) & 
                              fi(:,:,2)== CommonRGB(2) & 
                              fi(:,:,3)== CommonRGB(3))


  for c=1:3
       %find the linear index for 3 channels and replace
       idx = sub2ind([refRows ,refCols ,3], indicesX, indicesY,[c c]')
       fi(idx ) = bi(idx );
  end

  %counting swapped pixels
  swappedPixels = length(indicesX);