Can this code be improved?

I have written the following code in MATLAB to process large images of the order of 3000X2500 pixels. Currently the operation takes more than half hour to complete. Is there any scope to improve the code to consume less time? I heard parallel processing can make things faster but I have no idea on how to implement it. Can someone help me with suggestions:

function dirvar(subfn)
[fn,pn] = uigetfile({'*.TIF; *.tiff; *.tif; *.TIFF; *.jpg; *.bmp; *.JPG; *.png'}, ...
            'Select an image', '~/');
I = double(imread(fullfile(pn,fn)));
ld = input('Enter the lag distance = '); % prompt for lag distance
fh = eval(['@' subfn]); % Function handles
I2 = uint8(nlfilter(I, [7 7], fh));
imshow(I2); % Texture Layer Image
imwrite(I2,'result_mat.tif');

% Zero Degree Variogram
function [gamma] = ewvar(I)
    c = (size(I)+1)/2; % Finds the central pixel of moving window
    EW = I(c(1),c(2):end); % Determines the values from central pixel to margin of window
    h = length(EW) - ld; % Number of lags
    gamma = 1/(2 * h) * sum((EW(1:ld:end-1) - EW(2:ld:end)).^2);
end

The input lag distance is usually 1.

[Matt Fig]

Have you profiled the code with the profiler to see where the bottleneck is?

[Chethan S]

I have not done that. In fact I had no knowledge about it as I started using MATLAB very recently.

Answer 1

[Ben Mitch]

This way of doing processing, unfortunately because it's the easiest to program, is the slowest way to get things done in Matlab. Below, I've provided an implementation of your code using matrix operations instead of loops. It runs nearly 300 times faster on my machine than the looping code. With a bit of imagination, it could be as general as your original code, too.

This is a really nice illustration of how you need to think in Matlab if you are after performance - it's got to be all about the matrices.

OUTPUT:

 Elapsed time is 11.141832 seconds.
 Elapsed time is 0.039854 seconds.
 max error: 0

SCRIPT:

function dirvar()

% test image
I = double(imread('test.png'));

% seems to require greyscale?
I = mean(I, 3);

% pars
subfn = 'ewvar';
ld = 1;
L = 3;

% slow sequential implementation
tic
fh = eval(['@' subfn]); % Function handles
I2_slow = uint8(nlfilter(I, (2*L+1)*[1 1], fh));
toc

% fast parallel implementation
tic
I2 = zeros(size(I)); % output image
I = [I zeros(size(I, 1), L)]; % augment image
D = diff(I, 1, 2); % take derivative
D2 = D.^2; % take square
for o = 0:ld:L-1 % for each lag
  I2 = I2 + D2(:, (1:size(I2, 2)) + o);
end
I2 = uint8(I2);
toc

% compare results
E = double(I2) - double(I2_slow);
disp(['max error: ' num2str(max(abs(E(:))))]);

% show results
figure(1)
imshow(I2_slow); % Texture Layer Image
figure(2)
imshow(I2); % Texture Layer Image

% Zero Degree Variogram
  function [gamma] = ewvar(I)
    c = (size(I)+1)/2; % Finds the central pixel of moving window
    EW = I(c(1),c(2):end); % Determines the values from central pixel to margin of window
    h = length(EW) - ld; % Number of lags
    gamma = 1/(2 * h) * sum((EW(1:ld:end-1) - EW(2:ld:end)).^2);
  end

end