图像处理 : Algorithm taking too long in MATLAB

Question

我正在使用 MATLAB 处理两张 512x512 图像，域图像和距离图像。我想要完成的是以下内容:

1. 将域和范围图像划分为 8x8 像素 block
2. 对于域图像中的每个 8x8 block ，我必须对其应用线性变换并将 4096 个变换 block 中的每一个与 4096 个范围 block 中的每一个进行比较。
3. 在每种情况下计算变换 block 和范围图像 block 之间的误差，并找到最小误差。
4. 最后，对于每个 8x8 范围 block ，我将获得误差最小的 8x8 域 block 的 ID(范围 block 和转换后的域 block 之间的误差)

为此，我编写了以下代码:

RangeImagecolor = imread('input.png'); %input is 512x512
DomainImagecolor = imread('input.png'); %Range and Domain images are identical

RangeImagetemp = rgb2gray(RangeImagecolor);
DomainImagetemp = rgb2gray(DomainImagecolor);

RangeImage = im2double(RangeImagetemp);

DomainImage = im2double(DomainImagetemp);
%For the (k,l)th 8x8 range image block
for k = 1:64
    for l = 1:64
        minerror = 9999;
        min_i = 0;
        min_j = 0;

        for i = 1:64
            for j = 1:64

                 %here I compute for the (i,j)th domain block, the transformed domain block stored in D_trans

                 error = 0;
                 D_trans = zeros(8,8);
                 R = zeros(8,8); %Contains the pixel values of the (k,l)th range block
                 for m = 1:8
                     for n = 1:8
                         R(m,n) = RangeImage(8*k-8+m,8*l-8+n); 
                         %ApplyTransformation can depend on (k,l) so I can't compute the transformation outside the k,l loop.
                         [m_dash,n_dash] = ApplyTransformation(8*i-8+m,8*j-8+n);
                         D_trans(m,n) = DomainImage(m_dash,n_dash);
                         error = error + (R(m,n)-D_trans(m,n))^2; 
                     end 
                 end
                 if(error < minerror)
                     minerror = error; 
                     min_i = i;
                     min_j = j;
                 end
            end
        end 
    end
end

例如 ApplyTransformation，可以使用身份转换:

 function [x_dash,y_dash] = Iden(x,y)
    x_dash = x;
    y_dash = y;
end

现在我面临的问题是高计算时间。上面代码中的计算顺序是64^5，也就是10^9的顺序。这个计算应该在最坏的几分钟或一个小时内完成。计算仅 50 次迭代大约需要 40 分钟。我不知道为什么代码运行这么慢。

感谢阅读我的问题。

Answer 1

您可以使用 im2col* 将图像转换为列格式，以便每个 block 形成 [64 * 4096] 矩阵的一列。然后对每一列应用转换并使用 bsxfun向量化误差计算。

DomainImage=rand(512);
RangeImage=rand(512);
DomainImage_col = im2col(DomainImage,[8 8],'distinct');
R =  im2col(RangeImage,[8 8],'distinct');
[x y]=ndgrid(1:8);

function [x_dash, y_dash] = ApplyTransformation(x,y)
    x_dash = x;
    y_dash = y;
end

[x_dash, y_dash] = ApplyTransformation(x,y);
idx = sub2ind([8 8],x_dash, y_dash);

D_trans = DomainImage_col(idx,:);       %transformation is reduced to matrix indexing
Error = 0;
for mn = 1:64
    Error = Error + bsxfun(@minus,R(mn,:),D_trans(mn,:).').^2;
end
[minerror ,min_ij]= min(Error,[],2);         % linear index of minimum of each block;
[min_i  min_j]=ind2sub([64 64],min_ij); % convert linear index to subscript

解释:

我们的目标是尽可能减少循环次数。为此，我们应该避免矩阵索引，而应该使用向量化。嵌套循环应转换为一个循环。作为第一步，我们可以创建一个更优化的循环，如下所示:

min_ij = zeros(4096,1);

for kl = 1:4096                 %%% => 1:size(D_trans,2)
    minerror = 9999;
    min_ij(kl) = 0;
    for ij = 1:4096             %%% => 1:size(R,2)
        Error = 0;
        for mn = 1:64
            Error = Error + (R(mn,kl) - D_trans(mn,ij)).^2;
        end
        if(Error < minerror)
            minerror = Error; 
            min_ij(kl) = ij;
        end
    end
end

我们可以重新安排循环，我们可以将最内层的循环作为外层循环，并将最小值的计算与误差的计算分开。

% Computation of the error
Error = zeros(4096,4096);
for mn = 1:64
    for kl = 1:4096
        for ij = 1:4096
            Error(kl,ij) = Error(kl,ij) + (R(mn,kl) - D_trans(mn,ij)).^2;
        end
    end
end

% Computation of the min
min_ij = zeros(4096,1);
for kl = 1:4096
    minerror = 9999;
    min_ij(kl) = 0;
    for ij = 1:4096
        if(Error(kl,ij) < minerror)
            minerror = Error(kl,ij); 
            min_ij(kl) = ij;
        end
    end
end

现在代码以最适合矢量化的方式排列:

Error = 0;
for mn = 1:64
    Error = Error + bsxfun(@minus,R(mn,:),D_trans(mn,:).').^2;
end

[minerror ,min_ij] = min(Error, [], 2);   
[min_i  ,min_j] = ind2sub([64 64], min_ij);

*如果您没有图像处理工具箱，可以找到更有效的 im2col 实现 here .

*整个计算过程不到一分钟。