function [adj, prob, vertices]=tract2graph(SL, res)
% function tract2graph(tract, res)
%
% Construct 3D network graph out of fiber bundles using the epsilon-neighbor method. 
% The algorithm is first published in [1]. . The 4 possible cases follow the  
%
%
% SL        : collection of fiber bundles in the cell array format
% res       : resolution of seed. 6mm is recommanded
% adj       : adjacency matrix
% prob     : fiber concentration probability 
% vertices: the coordinates of the nodes
%
%
% (C) Moo K. Chung
%     University of Wisconsin-Madison
%     mkchung@wisc.edu
%
% Reference:
% [1] Chung, M.K. Adluru, N. Dalton, K. M., Alexander, A.L., Davidson, R.J. 2011. 
%      Scalable brain network construction on white matter fibers", Proc. SPIE 7962, 79624G; 
%      http://www.stat.wisc.edu/~mchung/papers/chung.2011.SPIE.pdf
%
% Update history: 
% October 5, 2009 
% July 8, 2010 
% December 25, 2011 Network constructio done on the sorted tract length. Case 3 bug fixed. 


%---------------------------------------
% obtain end points of tracts
endp=[];%zeros(number,2,3);
length=[]; %length of tracts
j=1;
for i=1:10000
    tract = SL{i};
    L =tract2length(tract);
    if L > 10 % tracts  with length 10mm is considered.
    %old code: 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
        length(j)=L;
        j=j+1;
    end;
end;

% sort the length of tracts in descending order and store the two end points depending on the length.   
[lengths, is]  = sort(length,'descend');
%figure; plot(lengths);
endp=endp(is,:,:);
%-------------------------------------------

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




% BUG: 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. 

%-----------------------------------------------
% construct a 3D graph
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
    n_node=size(nodelist,1); % number of nodes in the nodelist
    
    %----------------------------------------------------------------------
    % CHECK IF NODE1 is NEAR ANY NODE IN THE NODELIST
    node1=squeeze(endp(i,1,:))'; %get the 1st end point
    
    % compute the distance between node1 and the nodes in the nodelist
    repnode=repmat(node1,n_node,1);
    diff=sqrt(sum((nodelist-repnode).^2,2));
   
    % find nodes that are within the given resolution to node1
    ind1=find(diff<res); % index from nodelist 
   
    %----------------------------------------------------------------------
    % CHECK IF NODE2 IS NEAR ANY NODE IN THE NODELIST
    node2=squeeze(endp(i,2,:))';  %get the 2nd end point
    
    % compute the distance between node2 and the nodes in the nodelist
    repnode=repmat(node2,n_node,1);
    % checking if nodes are within 5mm radius
    diff=sqrt(sum((nodelist-repnode).^2,2));
    
    % find nodes that are within the given resolution to node2
    ind2=find(diff<res); % index from nodelist 
    
    %----------------------------------------------------------------------
    %THERE ARE 4 POSSIBILIES. SEE REFERENCE [1] FOR DETAIL
   

    % CASE 1. node1 and node2 are all close to nodes in the nodelist
    %         --> do nothing.
    %
    % 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 4. node1 and node2 are not close to any node in the nodelist
    %         --> add node1 and node2 to nodeliest and connect them in adj.

    if ~isempty(ind1) %node1 is in the nodelist: case 2 or 1
        if isempty(ind2) %node2 is not in the nodelist
            %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 not in the nodelist
        if isempty(ind2) %node2 not in the nodelist
        %CASE 4
            adj(n_node+1,n_node+2)=1;
            adj(n_node+2,n_node+1)=1;
            nodelist=[nodelist; node1; node2];
            %old code:
            %nodelist=[nodelist; node1];
            %nodelist=[nodelist; node2];
            prob(n_node+1)=1;
            prob(n_node+2)=1;
        else
            %CASE 3: node1 not in the nodelist. node2 in the nodelist.
            adj(ind2,n_node+1)=1;
            adj(n_node+1,ind2)=1;
            nodelist=[nodelist; node1];
            prob(ind2)=prob(ind2)+1;
            prob(n_node+1)=1;
        end;
    end;
end;

vertices=nodelist;

function arc_length =tract2length(tract)
%
% function [length  para]=parameterize_arclength(tract)
%
% Computes the arc length and perform the unit length parameterization of a
% curve.
%
% tract : 3 x n_vertex or 2 x n_vertex
% length: total arc length
% para  : number between 0 and 1 that maps the track to a unit interval.
%
%
% (C) Moo K. Chung & Nagesh Adluru 2008
%     mkchung@wisc.edu
%     University of Wisconsin-Madison
%
% Update history:
% 2008 created by M.K. Chung
% 2009 optimized by N. Adluru
% 2010 comments added

n_vertex=size(tract,2);
% n_vertex >=2 in order for this function to work.

p0=tract(:,1:(n_vertex-1));
p1=tract(:,2:n_vertex);
disp=p1-p0;
%January 22, 2009.
%Nagesh Adluru.
%Now computes directly: no function call overhead.
%L2=L2norm(disp');
L2=sqrt(sum(disp'.^2,2));
arc_length=sum(L2);
 