Vectorize 'for' loops that call other functions

Question

I have the following code snippet in matlab with two 'for' loops: 'I' is a binary image that has been preallocated.

 ...
    [x,y] = find(bwmorph(I,'endpoints'));
    n=numel(x);
    m=numel(x)-1;
    n=m+1;
    r=i+1;
    for i= 1:m
        for j = r:n
            I=linept(I, x(i), y(i), x(j), y(j));
        end;
    end;
    ...

The linept function is given below.Its from the Matlab File Exchange:

function result=linept(matrix, X1, Y1, X2, Y2)
result = matrix;
a=max(1, X1);b=sign(X2 - X1);c=max(1, X2);
for x=a:b:c
    y = round(f(x, X1, Y1, X2, Y2));
    if y > 0
        result(x, y) = 1;
    end
end
d=max(1, Y1);e=sign(Y2 - Y1);g=max(1, Y2);
for y=d:e:g
    x = round(f2(y, X1, Y1, X2, Y2));
    if x > 0
        result(x, y) = 1;
    end
end

function y=f(x, X1, Y1, X2, Y2)
a = (Y2 - Y1)/(X2 - X1);
b = Y1 - X1 * a;
y = a * x + b;

function x=f2(y, X1, Y1, X2, Y2)
if X1==X2
    x = X1;
else
    a = (Y2 - Y1)/(X2 - X1);
    b = Y1 - X1 * a;
    x = (y - b)/a;
end

Due to the many 'for' loops and function calls, this code is running very slowly.It runs fast for a simple image with few end points, but it takes a lot of time when the number of edges is more.It is slightly faster if the size of the image is reduced.I tried to vectorise it and have preallocated some variables, but there isn't much improvement.Can anyone help me regarding how to vectorise codes that call functions in loops.Thank you

Answer 1

This was one great problem!!

Brief discussion and solution codes

Well, listed next is a vectorized approach that heavily uses bsxfun at various places to take care of connecting all points with all other points by expansions which is basically how bsxfun operates . The code in it uses a sample input data for demo purposes. Take a look -

%// Create demo input data
img = false(20);
img(2,5:15) = 1;
img(12,5:15) = 1;
figure,imagesc(img), colormap gray, title('Starting Binary image')

%// Find endpoints
[X,Y] = find(bwmorph(img,'endpoints'));

%// Make a new binary image with only endpoints in it
I = false(size(img));
I(sub2ind(size(I),X,Y)) = 1;
figure,imagesc(I), colormap gray, title('Endpoints detected')

%-------- Vectorized code starts here ...
[nrows,ncols] = size(I);  %// Parameters
npts = numel(X);

twopts = nchoosek(1:npts,2);
slopes = (Y(twopts(:,1)) - Y(twopts(:,2)))./(X(twopts(:,1)) - X(twopts(:,2)));

%// Find the connecting indices with X-Y as they are, to work with 
%// slopes within [-1 1]
stage1 = abs(slopes)<=1;
[colID,rowID] = connecting_idx(X,Y,twopts(stage1,:),slopes(stage1));
valid = colID>0 & rowID>0 & colID<=ncols & rowID<=nrows;
I((colID(valid)-1)*nrows + rowID(valid))=1;

%// Find the connecting indices with X-Y interchanged, to work with 
%// slopes outside [-1 1]
[rowID,colID] = connecting_idx(Y,X,twopts(~stage1,:),1./slopes(~stage1));
valid = colID>0 & rowID>0 & colID<=ncols & rowID<=nrows;
I((colID(valid)-1)*nrows + rowID(valid))=1;
figure,imagesc(I),colormap gray,title('Final output : Endpoints connected')

Associated function code (most important part of codebase) -

function [y_all,x_all] = connecting_idx(X,Y,twopts,slopes)

%// Find XY indices that connects the X, Y anchor points given the two points
%// combinations and the corresponding slopes

X_twopts = reshape(X(twopts),size(twopts));
Y_twopts = reshape(Y(twopts),size(twopts));

[sortedX_pairs1,sorted_idx1] = sort(X_twopts,2);
X_starts1 = sortedX_pairs1(:,1);
Y_pairs = Y_twopts;
Y_starts = Y_pairs(sorted_idx1==1);
offsets = Y_starts - slopes.*X_starts1;
max_X_len = max(diff(sortedX_pairs1,[],2));

x_all = bsxfun(@plus,X_starts1,[0:max_X_len]);
x_lens = diff(sortedX_pairs1,[],2);
mask = bsxfun(@gt,x_lens+1,0:max_X_len);
y_all = round(bsxfun(@plus,bsxfun(@times,x_all,slopes),offsets));

y_all = y_all(mask);
x_all = x_all(mask);

return;

Debug images after code run -