clear all
clc
close all

% extract diameter-related metrics

Diameters = readmatrix('Day0_AdjacencyMatrices/Tumour 1/AdjacencyMatrixDiameters.txt'); % change for other tumours and for largest connected components
DiametersNZ = Diameters(Diameters~=0);  % non-zero elements
OrderedDiametersNZ = sort(DiametersNZ); % sort in ascending order
OrderedDiametersNZUnique = OrderedDiametersNZ(1:2:end); % get rid of duplicate values
MeanDiameter = mean(OrderedDiametersNZUnique) % mean
StDevDiameter = std(OrderedDiametersNZUnique) % standard deviation
MinDiameter = min(OrderedDiametersNZUnique) % min
MaxDiameter = max(OrderedDiametersNZUnique) % max

% radii related metrics are commented out because they are not used anymore - RadiiNZ =  DiametersNZ/2; %OrderedRadiiNZ = sort(RadiiNZ); %OrderedRadiiNZUnique = OrderedRadiiNZ(1:2:end); %MeanRadius = mean(OrderedRadiiNZUnique)

VesselCount = length(OrderedDiametersNZUnique)  % vessel count

% extract node connectivity

AdjacencyMatrix = Diameters;
AdjacencyMatrix(AdjacencyMatrix ~= 0) = 1;  % create adjacency matrix
NumberOfNeighbours = sum(AdjacencyMatrix);   % vector with number of neighbours for each node
NodesWith1Neighbour = sum(NumberOfNeighbours(:,:)==1) % inlets/outlets/sprout tips
NodesWith2Neighbours = sum(NumberOfNeighbours(:,:)==2)  % vessels split into 2 segments (a few occasions)
NodesWith3Neighbours = sum(NumberOfNeighbours(:,:)==3)  % bifurcation nodes
NodesWith4Neighbours = sum(NumberOfNeighbours(:,:)==4)  % trifurcations

% extract length-related metrics

Lengths = readmatrix('Day0_AdjacencyMatrices/Tumour 1/AdjacencyMatrixLengths.txt'); % change for other tumours and for largest connected components
LengthsNZ = Lengths(Lengths~=0); % non-zero elements
OrderedLengthsNZ = sort(LengthsNZ); % sort in ascending order
OrderedLengthsNZUnique = OrderedLengthsNZ(1:2:end); % get rid of duplicate values
MeanLength = mean(OrderedLengthsNZUnique) % mean
StDevLength = std(OrderedLengthsNZUnique) % standard deviation

% extract resistance-related metrics

Resistances = LengthsNZ./(DiametersNZ.^4);   % geometric resistances
OrderedResistances = sort(Resistances);  % sort in ascending order
OrderedResistancesUnique = OrderedResistances(1:2:end); % odd-indexed elements to get rid of duplicates
MeanGeometricResistance = mean(OrderedResistancesUnique) % mean geometric resistance
StDevGeometricResistance = std(OrderedResistancesUnique) % standard deviation in geometric resistance

Mu45 = 6*exp(-0.085*OrderedDiametersNZUnique)+3.2-2.44*exp(-0.06*OrderedDiametersNZUnique.^(0.645)); % mu_{45}
MuRel = (1+(Mu45-1).*(OrderedDiametersNZUnique./(OrderedDiametersNZUnique-1.1)).^2).*(OrderedDiametersNZUnique./(OrderedDiametersNZUnique-1.1)).^2; % relative viscosity
ResistancesWithVisc = 128/pi*MuRel.*OrderedLengthsNZUnique./(OrderedDiametersNZUnique.^4); % "resistance with viscosity", which in turn is dependent on diameter
MeanResistanceWithVisc = mean(ResistancesWithVisc) % mean "resistance with viscosity"
StDevResistanceWithVisc = std(ResistancesWithVisc) % standard deviation in "resistance with viscosity"

