matlab - 如何减少for循环消耗的时间？

Question

我正在尝试实现一个简单的像素级中心环绕图像增强。中心环绕技术利用窗口中心像素与周围邻域之间的统计数据作为决定需要进行何种增强的手段。在下面给出的代码中，我将中心像素与周围信息的平均值进行了比较，并基于此在两种情况之间切换以增强对比度。我写的代码如下:

im = normalize8(im,1);     %to set the range of pixel from 0-255
s1 = floor(K1/2);          %K1 is the size of the window for surround
M = 1000;                  %is a constant value
out1 = padarray(im,[s1,s1],'symmetric');
out1 = CE(out1,s1,M);
out = (out1(s1+1:end-s1,s1+1:end-s1));
out = normalize8(out,0);   %to set the range of pixel from 0-1



function [out] = CE(out,s,M)
  B = 255;
  out1 = out;
  for i = s+1 : size(out,1) - s
    for j = s+1 : size(out,2) - s
        temp = out(i-s:i+s,j-s:j+s);
        Yij = out1(i,j);
        Sij = (1/(2*s+1)^2)*sum(sum(temp));
          if (Yij>=Sij)
            Aij = A(Yij-Sij,M);
            out1(i,j) = ((B + Aij)*Yij)/(Aij+Yij);
          else
            Aij = A(Sij-Yij,M);
            out1(i,j) = (Aij*Yij)/(Aij+B-Yij);
          end
      end
   end
out = out1;
function [Ax] = A(x,M)
   if x == 0
       Ax = M;
   else 
       Ax = M/x;
    end

代码做了以下事情:
1) 将图像归一化到 0-255 范围内，并用额外的元素填充它以执行窗口操作。
2) 调用函数CE。
3) 在CE函数中获取窗口图像(temp)。
4) 求窗口的平均值(Sij)。
5) 比较窗口中心(Yij)与平均值(Sij)。
6) 根据比较结果进行两种增强操作之一。
7) 最后将范围设置回 0-1。

我必须针对多个窗口大小(K1、K2、K3 等)运行此程序，图像的大小为 1728*2034。当窗口大小选择为100时，消耗的时间非常高。

我可以在某个阶段使用矢量化来减少循环时间吗？

分析器结果(窗口大小为 21)如下:

分析器结果(窗口大小为 100)如下:

我已经更改了我的函数的代码，并在没有子函数的情况下编写了它。代码如下:

function [out] = CE(out,s,M)
B = 255;
Aij = zeros(1,2);
out1 = out;
n_factor = (1/(2*s+1)^2);
for i = s+1 : size(out,1) - s
    for j = s+1 : size(out,2) - s
        temp = out(i-s:i+s,j-s:j+s);
        Yij = out1(i,j);
        Sij = n_factor*sum(sum(temp));
        if Yij-Sij == 0
            Aij(1) = M;
            Aij(2) = M;
        else
            Aij(1) = M/(Yij-Sij);
            Aij(2) = M/(Sij-Yij);
        end
        if (Yij>=Sij)
            out1(i,j) = ((B + Aij(1))*Yij)/(Aij(1)+Yij);
        else
            out1(i,j) = (Aij(2)*Yij)/(Aij(2)+B-Yij);
        end
    end
end
out = out1;

从 93 秒到 88 秒，速度略有提高。欢迎对我的代码提出任何其他改进建议。

我尝试结合给出的建议，用卷积替换滑动窗口，然后对其余部分进行矢量化。下面的代码是我的实现，我没有得到预期的结果。

function [out_im] = CE_conv(im,s,M)
B = 255;
temp = ones(2*s,2*s);
temp = temp ./ numel(temp);
out1 = conv2(im,temp,'same');
out_im = im;
Aij = im-out1;                             %same as Yij-Sij
Aij1 = out1-im;                            %same as Sij-Yij
Mij = Aij;
Mij(Aij>0) = M./Aij(Aij>0);                % if Yij>Sij  Mij = M/Yij-Sij;
Mij(Aij<0) = M./Aij1(Aij<0);               % if Yij<Sij  Mij = M/Sij-Yij;
Mij(Aij==0) = M;                           % if Yij-Sij == 0 Mij = M;
out_im(Aij>=0) = ((B + Mij(Aij>=0)).*im(Aij>=0))./(Mij(Aij>=0)+im(Aij>=0));
out_im(Aij<0) = (Mij(Aij<0).*im(Aij<0))./ (Mij(Aij<0)+B-im(Aij<0));

我不知道哪里出了问题。

以下论文详细解释了我要实现的内容:
Vonikakis、Vassilios 和 Ioannis Andreadis。 “多尺度图像对比度增强。” In Control, Automation, Robotics and Vision, 2008. ICARCV 2008. 第 10 届国际 session ，第 856-861 页。 IEEE, 2008.

Answer 1

我试着看看我是否可以通过使用 colfilt 和 nlfilter 来处理这些时间，因为两者通常比滑动窗口的 for 循环快得多图像处理。 两者都适用于相对较小的 window 。对于 2048x2048 像素的图像和 10x10 的窗口，使用 colfilt 的解决方案大约需要 5 秒(在我的个人电脑上)。使用 21x21 的窗口，时间跳到了 27 秒，但这仍然是问题上显示时间的相对改进。不幸的是，我没有足够的内存来使用 100x100 的窗口进行 colfilt，但是使用 nlfilter 的解决方案有效，尽管需要大约 120 秒。

这里是代码

使用 colfilt 的解决方案:

function outval = enhancematrix(inputmatrix,M,B)
%Inputmatrix is a 2D matrix or column vector, outval is a 1D row vector.

% If inputmatrix is made of integers...
inputmatrix = double(inputmatrix);

%1. Compute S and Y
normFactor = 1 / (size(inputmatrix,1) + 1).^2; %Size of column.
S = normFactor*sum(inputmatrix,1); % Sum over the columns. 
Y = inputmatrix(ceil(size(inputmatrix,1)/2),:); % Center row.
% So far we have all S and Y, one value per column.

%2. Compute A(abs(Y-S)) 
A = Afunc(abs(S-Y),M);
% And all A: one value per column.

%3. The tricky part. If Y(i)-S(i) > 0 do something.
doPositive = (Y > S);
doNegative = ~doPositive;

outval = zeros(1,size(inputmatrix,2));

outval(doPositive) = (B + A(doPositive) .* Y(doPositive)) ./ (A(doPositive) + Y(doPositive));
outval(doNegative) = (A(doNegative) .* Y(doNegative)) ./ (A(doNegative) + B - Y(doNegative));

end

function out = Afunc(x,M)
% Input x is a row vector. Output is another row vector.
    out = x;
    out(x == 0) =  M;
    out(x ~= 0) = M./x(x ~= 0);
end

And to call it, simply do:

M = 1000; B = 255; enhancenow = @(x) enhancematrix(x,M,B);
w = 21 % windowsize
result = colfilt(inputImage,[w w],'sliding',enhancenow);

使用 nlfilter 的解决方案:

function outval = enhanceimagecontrast(neighbourhood,M,B)

%1. Compute S and Y
normFactor = 1 / (length(neighbourhood) + 1).^2;
S = normFactor*sum(neighbourhood(:));
Y = neighbourhood(ceil(size(neighbourhood,1)/2),ceil(size(neighbourhood,2)/2));


%2. Compute A(abs(Y-S))
test = (Y>=S);
A = Afunc(abs(Y-S),M);

%3. Return outval
if test
    outval = ((B + A) * Y) / (A + Y);
else
    outval = (A * Y) / (A + B - Y);
end


function aval = Afunc(x,M)
if (x == 0)
    aval = M;
else
    aval = M/x;
end

And to call it, simply do:

M = 1000; B = 255; enhancenow = @(x) enhanceimagecontrast(x,M,B);
w = 21 % windowsize
result = nlfilter(inputImage,[w w], enhancenow);

我没有花太多时间检查一切是否 100% 正确，但我确实看到了一些不错的对比度增强(头发看起来特别漂亮)。