图像处理:算法在MATLAB中耗时过长
我在MATLAB中处理两个512x512图像，域图像和范围图像。我想要实现的目标如下：
- 将域和范围图像分成8x8像素块
- 对于域图像中的每个8x8块，我必须对其应用线性变换，并将每个4096转换块与每个4096范围块进行比较。
- 计算转换块和范围图像块之间的每种情况下的误差，并找到最小误差。
- 最后，对于每个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分钟。我不知道为什么代码运行得这么慢。
感谢您阅读我的问题。
回答 2
Stack Overflow用户
发布于 2018-05-29 17:47:55
您可以使用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);*如果您没有图像处理工具箱，可以在here中找到更有效的im2col实现。
*整个计算过程不到一分钟。
Stack Overflow用户
发布于 2018-05-29 17:08:44
首先要做的是--你的代码什么也不做。但是你可能做了一些错误最小化的事情，只是忘了在这里粘贴，或者仍然需要编写代码。现在先别管了。
代码的一个大问题是计算64x64块的结果图像和源图像的转换。复杂操作的64^5次迭代肯定会很慢。相反，您应该一次计算所有转换并保存它们。
allTransMats = cell(64);
for i = 1 : 64
for j = 1 : 64
allTransMats{i,j} = getTransformation(DomainImage, i, j)
end
end
function D_trans = getTransformation(DomainImage, i,j)
D_trans = zeros(8);
for m = 1 : 8
for n = 1 : 8
[m_dash,n_dash] = ApplyTransformation(8*i-8+m,8*j-8+n);
D_trans(m,n) = DomainImage(m_dash,n_dash);
end
end
end这是用来获取allTransMat的，它在k，l循环之外。优选地，作为简单函数。
现在，你做你的k，l，i，j循环，根据需要比较所有的元素。比较也可以按块进行，而不是填充一个小的8x8矩阵，但由于某些原因，它是按元素进行的。
m = 1 : 8;
n = m;
for ...
R = RangeImage(...); % This will give 8x8 output as n and m are vectors.
D = allTransMats{i,j};
difference = sum(sum((R-D).^2));
if (difference < minDifference) ...
end尽管这是一个简单的无转换情况，但这大大加快了代码的速度。
最后，您确定需要将转换后的输出的每个块与源代码中的每个块进行比较吗？通常，在相同位置上将block1(a,b)与block2(a,b) -block(或像素)进行比较。
编辑: allTransMats也需要k和l。唉哟。没有办法在单次迭代中做到这一点，因为您需要64^5次调用ApplyTransformation (或该函数的矢量化，但即使这样也可能不快-我们必须在这里看到该函数的帮助)。因此，我将重复我的建议以生成所有转换，然后执行查找:生成allTransMats的答案的上半部分应该更改为具有所有4个循环并生成allTransMats{i,j,k,l};。它会很慢，就像我在编辑的上半部分提到的那样，没有办法解决这个问题。但是，这是你一次支付的成本，因为在保存allTransMats之后，所有进一步的图像分析都将能够简单地加载它，而不是再次生成它。
但是..。你到底是做什么的？依赖于源和目标块索引加上像素索引(总共=6个值)的转换听起来像是某个地方的错误，或者是优化的首选，而不是所有其他的。
https://stackoverflow.com/questions/50579013
复制相似问题