clear all;
close all;
%Load pre-processed network (Gaussian filtered, binarized, outliers
%removed)

%This script takes a z-stack that has been filtered and binarized and
%converts it to a 3D skeleton. This skeleton is then converted into a nodal
%network and analyzed to determine its size and relative degree of
%branching. A large section of this code relies on functions published by
%Phillip Kollmannsberger on the MATLAB File Exchange, as detailed in the
%main manuscript. Use of this code is not restricted as long as it does not
%violate Kollmannsberger's stipulations and includes citations to his work
%and the manuscript attached to this code. 

FileList = dir('*.tif');
N = size(FileList,1); %Get all files with extension in directory
Data = cell(N,6);
f = waitbar(0,'Is the wind still rising?');

for k = 1:N
    x = k/N;
    f = waitbar(x,f); %create waitbar to gauge progress- useful when analyzing a large number of files
    FileTif = FileList(k).name; %use TIFF library to import filtered, binarized skeletons saved as .tif files
    InfoImage=imfinfo(FileTif);
    mImage=InfoImage(1).Width;
    nImage=InfoImage(1).Height;
    NumberImages=length(InfoImage);
    FinalImage=zeros(nImage,mImage,NumberImages,'logical');%create an empty array the same size as the image
    TifLink = Tiff(FileTif, 'r');
    %Read each slice and convert to binary from uint8 file
    for i=1:NumberImages
        TifLink.setDirectory(i);
        FinalImage(:,:,i)= TifLink.read();
        FinalImage(:,:,i)= FinalImage(:,:,i) > 0; %copy each stack to a new "slice" in the zeros array
    end
    TifLink.close();
    
    %Skeletonize the image
    skel = Skeleton3D(FinalImage);
    
    %Note: the code below is slightly modified from that provided by
    %Phillip Kollmannsberger; please use with caution and refer to the
    %original to ensure adequate export quality
    figure();
    col=[.7 .7 .8];
    hiso = patch(isosurface(FinalImage,0),'FaceColor',col,'EdgeColor','none');
    hiso2 = patch(isocaps(FinalImage,0),'FaceColor',col,'EdgeColor','none');
    axis equal;axis off;
    lighting phong;
    isonormals(FinalImage,hiso);
    alpha(0.5);
    set(gca,'DataAspectRatio',[1 1 1])
    camlight;
    hold on;
    w=size(skel,1);
    l=size(skel,2);
    h=size(skel,3);
    [x,y,z]=ind2sub([w,l,h],find(skel(:)));
    plot3(y,x,z,'square','Markersize',4,'MarkerFaceColor','r','Color','r');
    set(gcf,'Color','white');
    view(140,80)
    
    
    %Analyze the Skeleton
    w = size(skel,1);
    l = size(skel,2);
    h = size(skel,3);
    
    % initial step: condense, convert to voxels and back, detect cells
    [~,node,link] = Skel2Graph3D(skel,0);
    
    % total length of network
    wl = sum(cellfun('length',{node.links}));
    
    skel2 = Graph2Skel3D(node,link,w,l,h);
    [~,node2,link2] = Skel2Graph3D(skel2,0);
    
    % calculate new total length of network
    wl_new = sum(cellfun('length',{node2.links}));
    
    % iterate the same steps until network length changed by less than 0.5%
    while(wl_new~=wl)
        
        wl = wl_new;
        
        skel2 = Graph2Skel3D(node2,link2,w,l,h);
        [A2,node2,link2] = Skel2Graph3D(skel2,0);
        
        wl_new = sum(cellfun('length',{node2.links}));
        
    end
    
    dist = 0;
    for i=1:length(node2)
        for j=1:length(node2(i).links)
            for b=1:length(link2(node2(i).links(j)).point)-1
                [x3,y3,z3]=ind2sub([w,l,h],link2(node2(i).links(j)).point(b));
                [x2,y2,z2]=ind2sub([w,l,h],link2(node2(i).links(j)).point(b+1));
                dist = dist + sqrt((x3-x2)^2+(y3-y2)^2+(z3-z2)^2);
            end
        end
    end
    
    figure();
    hold on;
    for i=1:length(node2)
        x1 = node2(i).comx;
        y1 = node2(i).comy;
        z1 = node2(i).comz;
        
        if(node2(i).ep==1)
            ncol = 'c';
        else
            ncol = 'y';
        end
        
        for j=1:length(node2(i).links)    % draw all connections of each node
            if(node2(node2(i).conn(j)).ep==1)
                col='k'; % branches are black
            else
                col='r'; % links are red
            end
            if(node2(i).ep==1)
                col='k';
            end
            
            
            % draw edges as lines using voxel positions
            for q=1:length(link2(node2(i).links(j)).point)-1
                [x3,y3,z3]=ind2sub([w,l,h],link2(node2(i).links(j)).point(q));
                [x2,y2,z2]=ind2sub([w,l,h],link2(node2(i).links(j)).point(q+1));
                line([y3 y2],[x3 x2],[z3 z2],'Color',col,'LineWidth',2);
            end
        end
        
        % draw all nodes as yellow circles
        plot3(y1,x1,z1,'o','Markersize',9,...
            'MarkerFaceColor',ncol,...
            'Color','k');
    end
    axis image;axis off;
    set(gcf,'Color','white');
    drawnow;

Data(k,1) = cellstr(FileList(k).name); %Filename
Data(k,2) = num2cell(length(node2)); %Total number of end points
runsum = 0;
for i = 1:length(node2)
    runsum = runsum + node2(i).ep;
end
Data(k,3) = num2cell(length(node2)-runsum); %Branch points
Data(k,4) = num2cell(length(link2));%Links (e.g., branches)
Data(k,5) = num2cell(dist); %Network length
Data(k,6) = num2cell(h); %Slice number (depth)

end







