function [adj, prob, vertices]=tract2graph(SL, res)
% Construct 3D graph out of DTI tracts
% (C) Moo K. Chung
% mkchung@wisc.edu
% October 5, 2009
% 
%function tract2graph(tract, res)
% SL: collection of tracts in cell data format
% res: resolution of seed
% adj: adjacency matrix
% prob: fiber concentration probability 



% obtain end points of tracts
endp=[];%zeros(number,2,3);
j=1;
for i=1:10000
    tract = SL{i};
    if size(tract,1)> 10 %removes short tract with less than 10 points
        endp(j,1,:)=tract(1,:);  % one end
        endp(j,2,:)=tract(end,:);% the other end
        j=j+1;
    end;
end;


%figure;
%plot3(sgraph(:,1,1),sgraph(:,1,2),sgraph(:,1,3),'.')
%hold on;
%plot3(sgraph(:,2,1),sgraph(:,2,2),sgraph(:,2,3),'.r')


% construct a 3D graph

%
% BUG 2: nodevalue computes the degree of a node but this does not seem to 
% match with the degree computed from the adjacency matrix.
% Note If you sum up column or row in the adjacency matrix, you get the adjaceny matrix. 
%
% Other than these not so serious bugs, 3D graph looks cool. Let's get some
% paper out of this.

adj=[]  % adjacent matrix

% the graph starts with two ends points of tract 1.
nodelist=[squeeze(endp(1,1,:))';squeeze(endp(1,2,:))'];

adj(1,2)=1;
adj(2,1)=1;

prob(1)=1;
prob(2)=1;


for i=2:size(endp,1);
    i
    node1=squeeze(endp(i,1,:))'; %get the 1st end point
    n_node=size(nodelist,1); % number of nodes in the nodelist
    
    %-------------------------------------
    % check if node1 is near any node in the nodelist
 
    repnode=repmat(node1,n_node,1);
    diff=sqrt(sum((nodelist-repnode).^2,2));
    
    % checking if nodes are within the given resolution
    ind1=find(diff<res); % index from nodelist 
   
    
    node2=squeeze(endp(i,2,:))';  %get the 2nd end point
    %-------------------------------------
    % check if node2 is near any node in the nodelist
 
    repnode=repmat(node2,n_node,1);
    % checking if nodes are within 5mm radius
    diff=sqrt(sum((nodelist-repnode).^2,2));
    
    ind2=find(diff<res); % index from nodelist 
    
    
    %there are 4 possibilies
   
    % case 4. node1 and node2 are not close to any node in the nodelist
    %         --> add them to nodeliest. connect node1 to all nodes in ind2
    %                                    connect node2 to all nodes in ind1
    %
    % case 2. node1 is a neighbor but not node2
    %         --> add node2 to nodelist, and connect node2 to all nodes in ind1 
    %
    % case 3. node2 is in the nodelist but not node 1
    %         --> add node1 to nodelist, and connect node1 to all nodes in ind2
    %
    % case 1. node1 and node2 are all close to nodes in the nodelist
    %         --> do nothing.

    if ~isempty(ind1)
        %node1 is in the nodelist: case 2 or 1
        if isempty(ind2)
            % case 2: add node2
            adj(ind1,n_node+1)=1;
            adj(n_node+1,ind1)=1;
            nodelist=[nodelist; node2];
            prob(ind1)=prob(ind1)+1;
            prob(n_node+1)=1;
        else
            %case 1
            prob(ind1)=prob(ind1)+1;
            prob(ind2)=prob(ind2)+1;
        end;
    else % node1 is not in the list
        if isempty(ind2)
        %case 4: node2 not in the list
        
            adj(n_node+1,n_node+2)=1;
            adj(n_node+2,n_node+1)=1;
            nodelist=[nodelist; node1];
            nodelist=[nodelist; node2];
            prob(n_node+1)=1;
            prob(n_node+2)=1;
        else
            %case 3
            adj(ind2,n_node+1)=1;
            adj(n_node+1,ind2)=1;
            nodelist=[nodelist; node2];
            prob(ind2)=prob(ind2)+1;
            prob(n_node+1)=1;
        end;
    end;
end;

vertices=nodelist;

 