clear all
close all
topLevelFolder = pwd;       %  the current folder in which the OCTA subfolders (subjects)exist.
subfolders = dir (topLevelFolder)    %  Number of OCTA subfolders/subjects in the main folder
subfolders = subfolders(~ismember({subfolders(:).name},{'.','..'})) % remove ./ and ../ from

% Process Starts
for j = 1:length(subfolders) % k=1 to the number of OCTA subfloders
    display(j, 'Processing the case No.: ')
    currentsubfolder = subfolders(j).name; % Get the current OCTA subfolder (subject) file names
    dr3=dir([currentsubfolder '/*.png']); % get only png format file
    f1={dr3.name}; % get only file names to cell
    for dd=1:length(f1) % for each MCV image in each OCTA subfolder/subject
       %close all; % this is in case the program is terminated at any internal step
       image1 =imread([currentsubfolder '/' f1{dd}]); % read original one MCV image
       display(dd, 'Processing the MCV No.: ')
       %figure,imshow(image1)
       [m n] = size(image1);    % size of MCV image
       checkSTDEV=std2(image1);  % Calculate teh STDev of the MCV
       if checkSTDEV<1   % this MCV image has no vessels (outside Macula)
         totalVAD(dd) = 0; % Set BVDI = 0 in this MCV
         totalVSD(dd) = 0; % Set BVSI = 0 in this MCV
         totalVT(dd) = 0;  % Set BVTI = 0 in this MCV
         continue
       end
       BWVessels=image1>25; % threshold the MCV to show blood vessels pixels shades > 25 (10% of the 256 shades)
       %figure,imshow(BWVessels);
       BWVessels = bwareaopen(BWVessels,5); % remove vessels less than 5 px in size
       %figure, imshow(BWVessels)
       BWSkeleton = bwskel(BWVessels);   % generate the skeleton of blood vessels
       %figure, imshow(BWSkeleteon)
       BWSkeleton = bwareaopen(BWSkeleton,5); % remove skeletons (of vessels) less than 5 px in size
       %figure, imshow(BWSkeleteon)
    
       % make a region for each Vesssel in the MCV
       [labeledImage1, numRegions1] = bwlabel(BWVessels);  % label all vessels in the MCV
       [labeledImage2, numRegions2] = bwlabel(BWSkeleton);  % label all skeletons in the MCV
 
       % Condition in case the MCV has no vessels
       if numRegions2 < 1    % no vessels
          totalVAD(dd) = 0;     % Set BVDI = 0 in this MCV
          totalVSD(dd) =  0;    % Set BVSI = 0 in this MCV
          totalVT(dd) = 0;      % Set BVTI = 0 in this MCV
          continue
        end
%   
      % start Calculating the OCTA Metrics
      % OCTA parmeter No.1 : calculate the VAD (i.e. BVDI)
      for k = 1 : numRegions1     % calculate the VAD for all vessels in the MCV
        thisRegion = ismember(labeledImage1, k);  % extract vessel no. "k"
        stats = regionprops('table', thisRegion, 'Area');
        VADtemp = stats.Area;
        VAD(k) = VADtemp;
      end
      totalVAD(dd) = sum(VAD(:))/(m * n); % or totalVAD(dd) = sum(VAD(:))/ numRegions1
      VAD(:) = 0; % set zero to the previous VAD calculation (to start calculating on new MCV)
  
     % OCTA parmeter No.2: calculate the VSD (i.e. BVSI)
      for k = 1 : numRegions2     % calculate the VSD for all vessels in the MCV
         thisRegion = ismember(labeledImage2, k);  % extract vessel no. "k"
         stats = regionprops('table', thisRegion, 'Area');
         VSDtemp = stats.Area;
         VSD(k) = VSDtemp;
      end
      totalVSD(dd) = sum(VSD(:))/(m * n); % or totalVSD(dd) = sum(VSD(:))/ numRegions2
      VSD(:) = 0; % set zero to the previous VSD calculation (to start calculating on new MCV)
    
     % OCTA parmeter No.3: calculate the VT
     for k = 1 : numRegions2     % calculate the average tortuosity for all vessels in the MCV
          thisRegion = ismember(labeledImage2, k);  % extract vessel no. "k"
          endpoints = bwmorph(thisRegion, 'endpoints');   % Find endpoints x,y coordinates
          % Get coordinate
         [r, c] = find(endpoints);    % endpoints of a vessel no. "k"
         if size(r)< 2 | size(c)<2  % no vessels
            VT(k)= 0;
            continue
         end
         euclideanDistance(k) = sqrt((c(2) - c(1))^2 + (r(2) - r(1))^2);  % calculate Euclidean distance for vessel no. "k"       
         area(k) = sum(thisRegion(:)); % calculate area for vessel no. "k" 
         VT(k) = area(k) / euclideanDistance(k); % tortuosity for  vessel no. K
     end
     totalVT(dd) = sum(VT(:))/(m*n);      % Average tortuosity for all skeleton vessels in MCV no. dd
     VT(:)=0;    % set zero to the previous calculation (to start calculating on new  MCV)
    end
    results1(j,:)=totalVAD;
    results2(j,:)=totalVSD;
    results3(j,:)=totalVT;
end
    % save results to excel files
    xlswrite('X1.xlsx',results1)  % send BVDI results to excell file
    xlswrite('X2.xlsx',results2)  % send BVSI results to excell file
    xlswrite('X3.xlsx',results3)  % send BVTI results to excell file
   

