How to reduce the time consumed by the for loop?

Question

I am trying to implement a simple pixel level center-surround image enhancement. Center-surround technique makes use of statistics between the center pixel of the window and the surrounding neighborhood as a means to decide what enhancement needs to be done. In the code given below I have compared the center pixel with average of the surrounding information and based on that I switch between two cases to enhance the contrast. The code that I have written is as follows:

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

The code does the following things:
1) Normalize the image to 0-255 range and pad it with additional elements to perform windowing operation.
2) Calls the function CE.
3) In the function CE obtain the windowed image(temp).
4) Find the average of the window (Sij).
5) Compare the center of the window (Yij) with the average value (Sij).
6) Based on the result of comparison perform one of the two enhancement operation.
7) Finally set the range back to 0-1.

I have to run this for multiple window size (K1,K2,K3, etc.) and the images are of size 1728*2034. When the window size is selected as 100, the time consumed is very high.

Can I use vectorization at some stage to reduce the time for loops?

The profiler result (for window size 21) is as follows:

The profiler result (for window size 100) is as follows:

I have changed the code of my function and have written it without the sub-function. The code is as follows:

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;

There is a slight improvement in the speed from 93 sec to 88 sec. Suggestions for any other improvements to my code are welcomed.

I have tried to incorporate the suggestions given to replace sliding window with convolution and then vectorize the rest of it. The code below is my implementation and I'm not getting the result expected.

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));

I am not able to figure out where I'm going wrong.

A detailed explanation of what I'm trying to implement is given in the following paper:
Vonikakis, Vassilios, and Ioannis Andreadis. "Multi-scale image contrast enhancement." In Control, Automation, Robotics and Vision, 2008. ICARCV 2008. 10th International Conference on, pp. 856-861. IEEE, 2008.

Answer 1

This answer is the implementation that was suggested by Peter. I debugged the implementation and presenting the final working version of the fast implementation.

function [out_im] = CE_conv(im,s,M)
B = 255;
im = ( im - min(im(:)) ) ./ ( max(im(:)) - min(im(:)) )*255;
h = ones(s,s)./(s*s);
out1 = imfilter(im,h,'conv');
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));
out_im = ( out_im - min(out_im(:)) ) ./ ( max(out_im(:)) - min(out_im(:)) );

To call this use the following code

I = imread('pout.tif');
w_size = 51;
M = 4000;
output = CE_conv(I(:,:,1),w_size,M);

The output for the 'pout.tif' image is given below

The execution time for Bigger image and with 100*100 block size is around 5 secs with this implementation.

Answer 2

I've tried to see if I could get those times down by processing with colfiltand nlfilter, since both are usually much faster than for-loops for sliding window image processing. Both worked fine for relatively small windows. For an image of 2048x2048 pixels and a window of 10x10, the solution with colfilt takes about 5 seconds (on my personal computer). With a window of 21x21 the time jumped to 27 seconds, but that is still a relative improvement on the times displayed on the question. Unfortunately I don't have enough memory to colfilt using windows of 100x100, but the solution with nlfilter works, though taking about 120 seconds.

Here the code

Solution with 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);

Solution with 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);

I didn't spend much time checking that everything is 100% correct, but I did see some nice contrast enhancement (hair looks particularly nice).