function [ MLINES, LENGTHS ] = ras2ml(RAS, DTHRESH)
% [ MLINES, LENGTHS ] = ras2ml(RAS, DTHRESH)
%
% Connect the points in RAS to form a list of multi-line objects.
%
% RAS is a cell array.
% DTHRESH is a scalar -- if not specified, default is 2.
%
% MLINES is a cell array of cell arrays.
% LENGTHS is a cell array.
%
% RAS{i} is a matrix of size [NDIM(i) x NPTS(i)], where
%   NDIM(i) is the number of dimensions of each point
%   NPTS(i) is the number of points
%
% MLINES{i}{j} is a matrix of size [NDIM{i} x NPML(i,j)], where
%   NPML(i,j) is the number of points on multiline j in space number i.
%
% LENGTHS{i}(j) is the sum of the lengths of all of the line segments which
%   make up the multiline object MLINES{i}{j}.
%
% Points are connected based on an exhaustive search of Euclidean distances.
% Multiline is extended in two directions from the starting point by making
% connections to the closest point.
% 
% If distance to the closest point is further than DTHRESH, the line is broken off,
% and a new one is started.
%
% Computing ML using Euclidean distances may be kind of time-consuming; you may find 
% that otl2ml does a similar thing more quickly.
%

% If RAS is a matrix, embed it in a cell array
if ~iscell(RAS),
  RAS = { RAS };
end

% If no threshold is specified, choose a threshold of 2
if nargin<2, DTHRESH=2; end
D2thresh = DTHRESH^2;

% Look at every point space in RAS
NC = length(RAS);
NDIM = zeros(1,NC);
NPTS = zeros(1,NC);

for ispace=1:NC,
  MLINES{ispace} = {};
  LENGTHS{ispace} = [];
  
  % Get the number of dimensions and number of points
  ras = RAS{ispace};
  [NDIM(ispace), NPTS(ispace)] = size(ras);
  
  % First, compute the distances between every set of points
  D2 = zeros(NPTS(ispace), NPTS(ispace));
  for idim=1:NDIM(ispace),
    D2 = D2 + (repmat(ras(idim,:)',[1 NPTS(ispace)]) - repmat(ras(idim,:),[NPTS(ispace) 1])).^2;
  end
  
  % Now, do the following loop as long as there are any points which have not been used
  while all(size(D2) > 0),
    
    % Pick an arbitrary starting point -- the one with the smallest total distance to remaining points
    [ totdist, startpt ] = min(sum(D2));
    ml = ras(:,startpt);
    D2(startpt,:) = Inf;
    lml = 0;
    
    % Extend the multiline both forward and backward from the starting point
    for direction=1:2,
      
      % Continue extending while there is any point with D < D2thresh in this direction
      [ mind, nextpt ] = min(D2(:,startpt));
      while mind < D2thresh,
	% Extend the multiline, and add sqrt(mind) to its total length
	ml = [ ml, ras(:,nextpt) ];
	lml = lml + sqrt(mind);
	
	% Mark this row as used
	D2(nextpt, :) = Inf;

	% Find the next row
	[ mind, nextpt ] = min(D2(:,nextpt));
      end
      
      % Flip the multi-line, so that we can extend in the other direction
      ml = fliplr(ml);
    end
    
    % Append this ml to MLINES
    MLINES{ispace} = { MLINES{ispace}{:} ml };
    
    % Append the total length to LENGTHS
    LENGTHS{ispace} = [ LENGTHS{ispace} lml ];
    
    % Get rid of the rows and columns which have been used
    unused = find(D2(:,1) < Inf);
    D2 = D2(unused, unused);
    ras = ras(:,unused);
    
  end
end  


	
  
    


