function [ MAP, MLINE ] = otl2ml( X, Y, SZ, NFILL )
% [ MAP, MLINE ] = otl2ml( X, Y, SZ, NFILL )
%
% Round off the points in X and Y to the nearest integers, and connect to form multilines.
%
% The returned MAP contains ones on every line segment in the multline, and zeros elsewhere.
%
% MLINE is a cell array of multiline objects. Each multiline object is an Nx2 array, which contains
% an ordered list of N points. The points are arranged in order along a multiline path, from some
% starting point to some other ending point.
%
% Points are connected on the basis of proximity, as measured using a pseudo-Hamming distance
% criterion. If you want to connect points on the basis of Euclidean distance, with no
% rounding of coordinates, look at ras2ml.m.
%
% Find the unique entries, so 'sparse' doesn't complain
YX = unique( max(1, round([min(SZ(1),Y(:)) min(SZ(2), X(:)) ])), 'rows');
[n,m] = size(YX);
% Create a sparse MAP, and an empty MLINE
MAP = sparse(YX(:,1), YX(:,2), ones(n,1), SZ(1), SZ(2));
MLINE = {};
% Return if zero length
if length(X)==0, return; end;
% Make a copy of the map, so we can keep track of used and unused points
UNUSED = sparse(MAP);
% While there are still nonzero points in UNUSED, do the following loop
while nnz(UNUSED)>0,
[Y, X] = find(UNUSED);
newpt = [Y(1) X(1)];
UNUSED(newpt(1),newpt(2)) = 0;
mline = newpt;
% We look for extensions in two directions, both starting from the first point
for direction=1:2,
% This loop looks for another extension, and stops when no extension is found
oldpt = [-1000 -1000];
while any(newpt ~= oldpt),
oldpt = newpt;
% This loop starts with zero Hamming distance, and expands as necessary to find an extension
hammingdist = 0;
while newpt == oldpt & hammingdist <= NFILL+1,
hammingdist = hammingdist+1;
upperleft = max(1, newpt-hammingdist);
lowerright = min(SZ, newpt+hammingdist);
[ey, ex] = find(UNUSED(upperleft(1):lowerright(1), upperleft(2):lowerright(2)));
% If ey contains entries, then an extension was found
if length(ey) > 0,
% Create all of the newpts, and find the one with shortest distance to oldpt
newpts = [upperleft(1)+ey-1, upperleft(2)+ex-1];
[ mindist2, inew ] = min((newpts(:,1)-oldpt(1)).^2 + (newpts(:,2)-oldpt(2)).^2);
newpt = newpts(inew(1),:);
% Add the newpt to the mline, and erase it from UNUSED
mline = [mline; newpt];
UNUSED(newpt(1),newpt(2)) = 0;
% Now fill in all pixels in a straight line between oldpt and newpt
mindist = sqrt(mindist2);
for i=1:mindist,
pt = round((i*newpt + (mindist-i)*oldpt)/mindist);
MAP(pt(1), pt(2)) = 1;
end;
end;
end;
end;
% Now we flip the mline, and start looking for extensions in the other direction
newpt = mline(1,:);
mline = flipud(mline);
end;
% Add the new mline to the cell array MLINE
MLINE = { MLINE{:} mline };
end;