%
% nl2m.m
%
% Neurolucida v4.3 to Matlab Converter
%
% [As,indicies,list,len,lpoint] = convert(filename)
%
% filename    is a .asc file produced by Neurolucida v4.3
%
% As 	    is the surface area of the soma
% indicies  is the tree structure of the branches in base 3 numbering
% list      is a cell array of radii for all the dendrites
% len 	    is a vector of branch lengths 
% lpoint    is a cell array of accumulated lengths along branches
%

function [As,indicies,list,len,lpoint] = convert(filename);

fid = fopen(filename,'r'); 
xfid = fopen('xdat','w');
yfid = fopen('ydat','w');
afid = fopen('adat','w');
ifid = fopen('idat','w');
Dfid = fopen('Ddat','w');

word = ''; 		% word on a line
number = [];		% number on a line
isNumber = 0;		% 0 if not number, 1 if number
num = 0;
index = 1;		%index of dendrites

x = [];
y = [];

storeword = 0; % switch to store word

line = fgetl(fid);
word = getword(line);	%returns any words found on a line

% loop that runs through lines until CellBody data is found

while ~strcmp(word,'CellBody')
   
%   disp('Finding Cell body');
   line = fgetl(fid);
   word = getword(line);
   
end

disp(' ')
%disp('Cell body found, Processing');

line = fgetl(fid);
word = getword(line);

while ~strcmp(word,'Endofcontour'), % loops through cellbody data
   
   num = 0;
   number = '';
   
   for j=1:length(line),	% goes through each character in a line
      
      if line(j) ~= ' ' & line(j) ~= '(' & num <= 2, 
         % the character is not a space or '(' then it is a number 
         % and stored in number
     	 number = [number line(j)];
         isNumber = 1; 
      elseif isNumber,				
         % if a number has just been read in and the character 
         % is not a number, then reset
         num = num + 1;
         isNumber = 0;
         storeword = 1;
      end
         
      if num == 1 & storeword,  %store x coordinates in x vector, y in y vector
         x = [x str2num(number)];
         number = '';
         storeword = 0;
      elseif num == 2 & storeword,
         y = [y str2num(number)];
         number = '';
         storeword = 0;
      end      

   end	% end for j
   
   line = fgetl(fid);
   word = getword(line);

end  % end while
   
As = abs(area(x,y));	% calculates area of soma
disp(['The cross-sectional area of the soma is ' num2str(As)]);
disp(' ')
   
line = fgetl(fid);
word = getword(line);	% returns any words found on a line
   
% loop that runs through lines until branch data is found

while ~strcmp(word,'Dendrite')
   line = fgetl(fid);
   word = getword(line);
%   disp('Finding 1st Dendrite');
end
   
%disp('1st Dendrite Found');
line = fgetl(fid);
word = getword(line);
num = 0;
number = '';
isNumber = 0;
x = [];
y = [];
r = [];
indicies = [];
len = [];
index = 1;

D = 0;			% number of dendrites
   
while ~feof(fid),
       
   dend = getdendnum(line);

   if dend ~= -1,

      num = 0;
      number = '';
         
      for j=1:length(line),	% goes through each character in a line

      	  if line(j) ~= ' ' & line(j) ~= '(' & line(j) ~= ')' & num <= 4, 

             % the character is not a space or '(' then it is a number
             % and stored in number

             number = [number line(j)];
             isNumber = 1; 

  	  elseif isNumber,

             % if a number has just been read in and the character is 
             % not a number, then reset

             num = num + 1;
             isNumber = 0;
             storeword = 1;

          end
         
      	  if num == 1 & storeword,  

	      %store x coordinates in x vector, y in y vector

              x = [x str2num(number)];
              number = '';
              storeword = 0;

          elseif num == 2 & storeword,

              y = [y str2num(number)];
              number = '';
              storeword = 0;

          elseif num == 3 & storeword,

              number = '';
              storeword = 0;

          elseif num == 4 & storeword,

              r = [r ' ' number];
              number = '';
              storeword = 0;

          end

      end   % end for

    end   % end if dend
    
    line = fgetl(fid);
    word = getword(line);
    
    if dend ~= getdendnum(line) & dend > -1,
        
        disp(['Processing branch ' num2str(dend)]);  
        drawnow 
        indicies = [indicies dend];
        if dend == 0,
           D = D + 1;
        end
        arc = getlength(x,y);
        %text(mean(x),mean(y),num2str(dend),'color','r');
        len = [len floor(arc(length(arc)))];
        fwrite(xfid,x,'double');
        fwrite(yfid,y,'double');
        fwrite(afid,arc,'double');
        fwrite(ifid,dend*ones(1,length(x)),'double');
        fwrite(Dfid,D*ones(1,length(x)),'double');
        lpoint(index,1) = {num2str(arc)}; 
        list(index,1) = {r};
        index = index + 1;
        r = [];
        x = [];
        y = [];

    end
      
end   % end while
   
% sort elements of indicies, list, len and lpoint by dendrite 

%disp(' ')
%disp('Prior to sorting')
%disp(' ')
%disp('The branch indicies are')
%disp(indicies)
%disp('Their respective lengths are')
%disp(len)
%disp(' ')

rb = find(indicies==0);                 % roots of dendrites
D = length(rb);                         % number of dendrites

for d=1:D,

    if d < D
       ind = rb(d):rb(d+1)-1;
    else
       ind = rb(d):length(indicies);
    end

    [indicies(ind),I] = sort(indicies(ind));
    listtemp = list(ind);
    list(ind) = listtemp(I);
    lentemp = len(ind);
    len(ind) = lentemp(I);
    lpointtemp = lpoint(ind);
    lpoint(ind) = lpointtemp(I);

end

hold off
disp(' ')
disp('Finished Conversion')
disp(' ')
disp('The branch indicies are')
disp(indicies)
disp('Their respective lengths are')
disp(len)
disp(' ')
   
return

function arc = getlength(x,y)
arc = 0;
sum = 0;
for i=2:length(x),
   sum = sum + sqrt((x(i)-x(i-1))^2 + (y(i)-y(i-1))^2);
   plot([x(i) x(i-1)],[y(i) y(i-1)])
   hold on
   arc = [arc sum];
end
return

function word = getword(line)
word = '';
for j=1:length(line),
   if line(j) == '"',
      word = '';
      return;
   end
   if isletter(line(j)),
      word = [word line(j)];
   end
end
return

function dend = getdendnum(line)

button = 0; 		% 0 = wrong number, 1 = right number
number = ''; 		% initialize number to nothing
word = getword(line);
sw = length(word);
% disp(['in getdendnum: word = ' word])

if strcmp(word,'Root'),

   dend = 0;

elseif strcmp(word,'Incomplete') | ...
       strcmp(word,'Endofsplit') | ...
       strcmp(word,'Endoftree') | ...
       strcmp(word,'Dendrite') | ...
       strcmp(word,'Apical') | ...
       (sw>4 & strcmp(word(1:5),'Color')) | ...
       strcmp(word,'') | strcmp(word,'Normal'),

   dend = -1;

else 

   for j=1:length(line),

      if line(j) == ',',
         button = 0;
      end
      
      if line(j) ~= '-' & button,
         number = [number line(j)];
      end

      if line(j) == 'R',
         button = 1;
         number = 0;
      end

   end

   %convert unique branch ID from base 3 to base 10

   dend = base2dec(number,3);

end

if dend==[],
   disp(['dend is empty and line = ' line])
end

return

function  a = area(x,y)

% AREA(X,Y)  Calculates the area of the 2-dimensional polygon
%      formed by vertices with coordinates vectors  X and Y.
%   Kirill K. Pankratov, kirill@plume.mit.edu.
%   April 20, 1994

x = [x(:); x(1)]; y = [y(:); y(1)];
plot(x,y,'k')
hold on
lx = length(x);
a = abs( (x(2:lx)-x(1:lx-1))' * (y(1:lx-1)+y(2:lx)) )/2;

return