图像处理:MATLAB中的算法耗时太长

Question

我在MATLAB中处理两个512x512图像,域图像和范围图像.我想要完成的是以下内容:

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

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

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*将图像转换为列格式,以便每个块形成[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可以在此处找到更有效的实施方案.

*整个计算只需不到一分钟.