您可以使用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实施{。{3}}。

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