function [adj, prob, vertices]=points2graph(endp, res, choice)
% function [adj, prob, vertices]=points2graph(tract, res)
%
% Construct 3D network graph out of cloud points data
% using the epsilon-neighbor method
%
% endp     : coordinates of cloud points
% res     : resolution of seed. 6mm is recommanded
% adj     : adjacency matrix
% prob    : fiber concentration probability at node
% vertices: the coordinates of the nodes
% choice  : either 'rips' or 'epsilon'
%
%
% (C) Moo K. Chung
%     Department of Biostatistics and Medical Informatics 
%     University of Wisconsin-Madison
%     mkchung@wisc.edu
%
% October 5, 2009; July 8, 2010; October 1, 2010
% Code is based on tract2graph.m


% the graph starts with the first point.
nodelist=endp(1,:); % node list of graph
adj=0;  % adjacent matrix

ep_nodelist=endp(1,:);
ep_adj=0;


switch choice
    case 'rips'
        for i=2:size(endp,1);

            % the graph starts with the first vertex.

            point=endp(i,:); %current node
            n_node=size(nodelist,1); % number of nodes in the nodelist

            %-------------------------------------
            % check if node is near any node in the nodelist

            repnode=repmat(point,n_node,1);
            diff=sqrt(sum((nodelist-repnode).^2,2));

            ind=find(diff<=res); % index from nodelist

            if ~isempty(ind)
                adj(ind,n_node+1)=1;
                adj(n_node+1,ind)=1; %adj(n_node+1,ind)=-1;
                nodelist=[nodelist; point];
            else
                adj(:,n_node+1)=0;
                adj(n_node+1,:)=0;
                nodelist=[nodelist; point];
            end;

        end;

        prob=[];
        vertices=nodelist;
        
    case 'epsilon'
        
        % run rips complex to obtain adj and vertices of the full graph
        [adj, prob, vertices]=points2graph(endp, res, 'rips');
        
        n=size(endp,1); % total number of nodes in the rips complex
        
        for i=2:n
            point=endp(i,:); %current node
            n_node=size(nodelist,1); % number of nodes in the nodelist
            repnode=repmat(point,n_node,1);
            diff=sqrt(sum((nodelist-repnode).^2,2));
            child=find(diff<=res); % index from nodelist. it contains only children

            if ~isempty(child)
                % need to find the parents of children
                for j=1:length(child)
                    parent=find(adj(:,child(j))==1);
                    ep_adj(parent,i)=1;
                    ep_adj(i,parent)=1;
                end;
            else
                ep_nodelist=[ep_nodelist; point];
            end;

            nodelist=[nodelist; point];
        end;
        
        adj=ep_adj;
        prob=[];
        vertices=ep_nodelist;
end;

 