clear all; close all; clc
% ----------------------------------------
samples_no = [2, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17];
sample_names_on_title = [""; "A 1"; ""; ""; "B 1"; "B 2"; ...
    "F 1"; "F 2"; "C 1"; "D 1";  "D 2"; ...
    "E 1"; "E 2"; "D 1 t."; "D 1 b."; ...
    "A 2"; "C 2"];

save_flag = 1;
save_multip_struct = 0; 
plots_and_videos = 0;
excel_file = 0;
windows   = 1; 
%// Flag for saving files in Windows or MacOS
% ---------------------------------------
load colours.mat
if save_multip_struct
    load multip_struct.mat
end

for sam = 1 : length(samples_no)
    sample_no = samples_no(sam);
    sample = sample_names(sample_no);
    sample_name_on_title = sample_names_on_title(sample_no);
    load scan_data.mat
    load scan_load.mat
    % ------------------------------------------
    % ------------- READING FILE ---------------
    % ------------------------------------------
    % Excel files to be read:
    file_to_read = sprintf('%s.xlsx', sample);
    [~, sheet] = xlsfinfo(file_to_read);
    singlets = numel(sheet);
    multiplets = NaN(32,numel(sheet));
    diffuse_multiplets = NaN(32, numel(sheet));
    
    % Set up the movie.
    if windows
        newSubFolder = fullfile('D:', 'OneDrive - University of Southampton' , ...
            'MATLAB', sample, '\');
    else
        newSubFolder = fullfile('Users', 'sebastianrosini', 'OneDrive', ...
            'OneDrive - University of Southampton' , 'MATLAB', sample, '/');
    end
    % Finally, create the folder if it doesn't exist already.
    if ~exist(newSubFolder, 'dir')
        mkdir(newSubFolder);
    end
    if windows
        results_file = sprintf('%s%s_matlab.xlsx', newSubFolder, sample);
    else
        results_file = sprintf('%s%s_matlab.xlsx', newSubFolder, sample);
    end
    
    %// To delete a pre-existing version of the Excel file that could be wrong:
    if excel_file == 1
        if exist(results_file, 'file') == 2
            delete(results_file);
        end
    end
    
    
    if plots_and_videos
        if windows
            writerObj = VideoWriter(sprintf('%s%s_breaks', newSubFolder, sample), 'MPEG-4');
        else
            writerObj = VideoWriter(sprintf('%s%s_breaks', newSubFolder, sample), 'MPEG-4');
        end
        writerObj.FrameRate = 0.5; % How many frames per second.
        open(writerObj);
    end
    
    
    breaks_table = cell( numel(sheet)+1, 16);
    data = cell(1, numel(sheet));
    for oo = 1 : numel(sheet)
        
        data{oo} = xlsread( file_to_read, sheet{oo}, 'A:D');

        % --------- PARAMETERS SCANS -------------

        voxel = scan_data{10,sample_no}; % (µm), resolution of the scans
        df = scan_data{14,sample_no} / voxel;
        planar_distance = 2 * df; % no. of voxels
        ineffective_length = 70/voxel; 
        % no. of voxels - Anna Scott's length
        x_scan = scan_data{11,sample_no}; % (voxels)
        y_scan = scan_data{12,sample_no}; % (voxels)
        z_scan = scan_data{13,sample_no}; % (voxels)
        vol_singlet = scan_data{18,sample_no}; % voxels
        vol_duplet = scan_data{19,sample_no};
        vol_triplet = scan_data{20,sample_no};
        vol_quadruplet = scan_data{21,sample_no};
        vol_quintuplet = scan_data{22,sample_no};

        % --- RESULTS FOLDER & FILE CREATION ------

        % Create a name for a subfolder within the matlab folder.
        if windows
            newSubFolder = fullfile('D:', 'OneDrive - University of Southampton' , ...
                'MATLAB', sample, '\');
        else
            newSubFolder = fullfile('Users', 'sebastianrosini', 'OneDrive', ...
                'OneDrive - University of Southampton' , 'MATLAB', sample, '/');
        end
        % Create the folder if it doesn't exist already:
        if ~exist(newSubFolder, 'dir')
            mkdir(newSubFolder);
        end
        
        if windows
            results_file = sprintf('%s%s_matlab_v02.xlsx', newSubFolder, sample);
        else
            results_file = sprintf('%s%s_matlab_v02.xlsx', newSubFolder, sample);
        end
        % ------ RENAMING THE COLLECTED DATA -----
        Centroid(:,:) = data{oo}(:,1:3);
        volume(:,1) = data{oo}(:,4);
        Centroid_planar = Centroid(:, 1:2);
        
        % ---- PLANAR CLUSTERS -----
        % breaks that are touching 
        % Au_planar will be a cell containing on the first column the id of the
        % planar cluster repeated many times:
        % duplet - one time
        % triplet - two times
        % this is to add to the diffuse cluster evaluation
        Au_planar = cell(length(Centroid),1);
        k = 1;
        sing_count = 0; singlets_only = zeros( length(volume), 3);
        p_2 = 0;
        p_3 = 0;
        p_4 = 0;
        p_5 = 0;
        p_6 = 0;
        for ii = 1 : length(volume)
            if volume(ii) < vol_singlet
                sing_count = sing_count + 1;
                singlets_only(sing_count, :) = Centroid(ii,1:3);
            end
            if (volume(ii) > vol_singlet) && (volume(ii) <= vol_duplet)
                Au_planar{k,1} = ii;
                Au_planar{k,2} = Centroid(ii,1:3);
                k = k+1;
                p_2 = p_2 + 1;
            end
            if (volume(ii) > vol_duplet) && (volume(ii) <= vol_triplet)
                Au_planar{k,1} = [ii ii];
                Au_planar{k,2} = Centroid(ii,1:3);
                k = k+1;
                p_3 = p_3 + 1;
            end
            if (volume(ii) > vol_triplet) && (volume(ii) <= vol_quadruplet)
                Au_planar{k,1} = [ii ii ii];
                Au_planar{k,2} = Centroid(ii,1:3);
                k = k+1;
                p_4 = p_4 + 1;
            end
            if (volume(ii) > vol_quadruplet) && (volume(ii) <= vol_quintuplet)
                Au_planar{k,1} = [ii ii ii ii];
                Au_planar{k,2} = Centroid(ii,1:3);
                k = k+1;
                p_5 = p_5 + 1;
            end
            if (volume(ii) > vol_quintuplet)
                Au_planar{k,1} = [ii ii ii ii ii];
                Au_planar{k,2} = Centroid(ii,1:3);
                k = k+1;
                p_6 = p_6 + 1;
            end
        end
        Au_planar(k:end,:) = [];
        % Remember Au_planar will be modified afterwards in the code. For this
        % reason, and to make a fair total breaks count, a copy of
        % Au_planar_for_breaks_count is created.
        Au_planar_for_breaks_count = Au_planar;
        
        %// Au_planar shows a duplet ID once. It will now show it twice:
        ii = 1;
        for k = 1 : length(Au_planar_for_breaks_count(:,1))
            if Au_planar_for_breaks_count{k} > 0
                Au_planar_for_breaks_count(ii, 1) = Au_planar_for_breaks_count(k,1);
                Au_planar_for_breaks_count{ii}(length(Au_planar_for_breaks_count{ii})+1) = ...
                Au_planar_for_breaks_count{ii}(length(Au_planar_for_breaks_count{ii}));
                ii =ii +1;
            end
        end
        Au_planar_for_breaks_count(ii:end,:) = [];
        
        scaled_Centroid = (Centroid(:,:).*voxel);
        % ------- DIFFUSE CLUSTERS -----
        % PLANAR EVALUATION
        % Points are considered without the z coord. as laying on a plane and their
        % distances are assessed, to check whether is inferior of 2*d_f, in which case
        % they are candidates to be part of a cluster as we see it.
        [idx1, dist1] = rangesearch(Centroid_planar,Centroid_planar,planar_distance);
        results1 = [idx1, dist1];
        R = zeros(1000, 6);
        k = 1;
        for i = 1 : length(Centroid(:,1))
            if length(dist1{i,1}) > 1
                R(k, 1) = Centroid(i, 1);
                R(k, 2) = Centroid(i, 2);
                R(k, 3) = Centroid(i, 3);
                R(k, 4) = idx1{i, 1} (1);
                R(k, 5) = idx1{i, 1} (2);
                R(k, 6) = dist1{i, 1} (2);
                k = k + 1;
            end
        end
        R(k : end, :) = [];
        if R > 0
            R_2 = R(:,1:3); % matrix containing candidates for clusters. 
            % It has original IDs of the fibres
            % (i.e. their position in Centroid(:,:) ).
        end
        
        % T = clusterdata(R_2(:,2:4),70);
        % scatter3(R_2(:,2),R_2(:,3),R_2(:,4),100,T,'filled')
        % rotate3d on
        
        % DIFFUSE CLUSTERS EVALUATION:
        % 1. need to remove singlets from 'results1':
        % meaning that rangesearch
        % gives back all the breaks and says which 
        % ones are close and which do
        % not have a twin. Now I extrapolate ONLY
        % the breaks that have a close break
        % as identified by rangesearch:
        
        results_planar = cell(size(results1(:,1)));
        k = 1;
        for ii = 1 : length(results1(:,1))
            if length(results1{ii,1}) > 1
                results_planar(k,1) = results1(ii,1);
                results_planar(k,2) = results1(ii,2);
                k = k+1;
            end
        end
        
        % 2. Evaluate z- distance to check diffuse clustering:
        
        results_diffuse = cell(length(results_planar),1);
        jj = 1;
        ii = 1;
        k = 1;
        while ii <= length(results_planar)
            while jj < length (results_planar{ii})
                m = jj + 1;
                ref_elem = results_planar{ii}(jj);
                while m <= length(results_planar{ii})
                    if (abs(Centroid(results_planar{ii}(m),3) - Centroid(ref_elem,3))) ...
                            < ineffective_length
                        results_diffuse{k}(1) = ref_elem;
                        results_diffuse{k}(2) = results_planar{ii}(m);
                        k = k + 1;
                    end
                    m = m + 1;
                end
                jj = jj + 1;
            end
            ii = ii + 1;
            jj = 1;
        end
        results_diffuse(k:end, :) = [];
        
        % 3. Sort the clusters arrays and mantain unique elements
        
        for i = 1 : length(results_diffuse(:,1))
            results_diffuse{i,1} = sort(results_diffuse{i,1},2);
        end
        
        [Au, idx_uniquecell1 ,idx_uniquecell2] = uniquecell(results_diffuse); 
        % ATT: CONTENT OF Au IS MODIFIED 
        % IN THE CLUSTERS EVALUATION LOOP - BE CAREFUL
        
        % ------------
        % CLUSTERS EVALUATION
        % But until now one break can be in two different clusters and I want
        % to unify two clusters with a connecting break. For this reason, with
        % the following loop I intersect the cells containing a common break so
        % that at the end in Au2 that I will save as Au_diffuse there will be
        % the final diffuse clusters.
        
        Au2 = cell(length(Au),1);
        
        jj = 1; % index for saving clusters data
        ii = 1; % index for going through uncombined list
        while ii <= length(Au) % starts from the first line of Au - first cluster
            k = ii+1; % used to go on the following line of Au and compare to following cluster
            ref_vec = Au{ii}; % as my Au{ii} is going to be my reference to be compared to all 
            % the following cluster it is set as reference_vector
            while k <= length(Au) % until the end of Au 
                dum = intersect(ref_vec, Au{k}); % I call dum the intersection of the ref_vec 
                % with the following cluster
                if numel(dum) == 0 % if this intersection is zero so there are no elements in common
                    k = k+1; % go to the next cluster so increase the k to move further
                else
                    ref_vec = unique([ref_vec,Au{k}]); % otherwise take only once the elements
                    % in common between the arrays (between the two clusters then)
                    Au(k,:) = []; % and delete that array in the cell
                    k = ii+1;
                end
            end
            Au2{jj} = ref_vec;   % save that reference vector as the complete cluster that goes in Au2
            jj =jj+1;
            ii = ii+1;
        end
        Au2(jj:end,:) = []; % delete from Au2 all the empty arrays. 
        % You had to assign a dimension for Au2 before, but all that cell is not needed.
        
        % CREATE A Au_diffuse TO CHECK ONLY THE DIFFUSE CLUSTERS WITHOUT
        % CONSIDERING THE PLANAR ONES:
        
        Au_diffuse(:,:) = Au2(:,:);
        
        Au_diffuse_for_breaks_count = Au_diffuse;
        
        [~,ncols] = cellfun(@size, Au_diffuse_for_breaks_count);
        diffuse_clusters = ncols;
        
        for i = 1 : length(diffuse_clusters)
            Au_diffuse_for_breaks_count{i,1} = [num2str(i+1,'%02d'),'-plets diffuse'];
            Au_diffuse_for_breaks_count{i,2} = sum(diffuse_clusters==i+1);
        end

        % --DIFFUSE + PLANAR CLUSTERS GENERATION ---

        % With these 'for': Au2 before had only diffuse clusters without considering
        % if fibre break 1 for example was a planar multi-plet.This is now
        % considered so the diffuse duplet [1 3] is actually a [1 1 3] triplet.
        
        for ii = 1 : length(Au2)
            for jj = 1 : length(Au2{ii})
                for k = 1 : length(Au_planar)
                    for m = 1 : length(Au_planar{k})
                        if Au2{ii}(jj) == Au_planar{k}(m)
                            Au2{ii}(length(Au2{ii})+1) = Au_planar{k}(m);
                            Au_planar{k}(m) = 0;
                        end
                    end
                end
            end
        end
        
        % However, planar multiplets that are not part of a diffuse multi-plet are
        % not considered in the cluster count. With this, a planar cluster not part
        % of a diffuse cluster (for example [2 2]) is considered in the final
        % duplet count
        
        Au2_planar = cell(length(Au_planar),1);
        ii =1;
        for k = 1 : length(Au_planar(:,1))
            if Au_planar{k} > 0
                Au2_planar(ii, 1) = Au_planar(k,1);
                Au2_planar{ii}(length(Au2_planar{ii})+1) = Au2_planar{ii}(length(Au2_planar{ii}));
                ii =ii +1;
            end
        end
        Au2_planar(ii:end,:) = [];
        Au_planar(ii:end,:) = [];
        
        % Then, as I add the id of the planar cluster at the end of the vector, I
        % have to sort the arrays in the cell containing all the breaks:
        
        for i = 1 : length(Au2(:,1))
            Au2{i,1} = sort(Au2{i,1},2);
        end
        
        % Now, I combine Au2 containing planar+diffuse clusters together and
        % Au2_planar containing PLANAR clusters only obtaining the Au2_total, with
        % all the clusters.
        
        Au2_total = cell((length(Au2))+(length(Au2_planar)),1);
        Au2_total(1:length(Au2)) = Au2;
        Au2_total(length(Au2)+1:end) = Au2_planar;
        
        [nrows,ncols] = cellfun(@size, Au2_total);
        clusters = ncols;
        clusters_length = length(clusters);
        if clusters_length > 0
            A_ALL = cell(clusters_length, 2);
        end
        
        for i = 1 : length(clusters)
            A_ALL{i,1} = [num2str(i+1,'%02d'),'-plets'];
            A_ALL{i,2} = sum(clusters==i+1);
        end
        % --- EXPORT DATA --
        % Coordinates of the clusters are in voxels
        % - need to be multiplied by 0.65
        clusters_coord = string( zeros(1000, 4) );
        
        if isempty(Au2_total) == 0
            k = 1;
            for ii = 1 : length(Au2_total)
                clusters_coord(k,1) = "Cluster no.:";
                clusters_coord(k,2) = ii;
                clusters_coord(k,3:4) = "";
                k = k + 1;
                for jj = 1 : length(Au2_total{ii})
                    clusters_coord(k,1) = Au2_total{ii}(jj);
                    clusters_coord(k,2:4) = Centroid(Au2_total{ii}(jj),1:3);
                    k = k + 1;
                end
            end
        end
        clusters_coord(k : end, :) = [];
        
        Centroid_labeled = zeros(length(Centroid(:,1)),4);
        for ii = 1 : length(Centroid(:,1))
            Centroid_labeled(ii,1) = ii;
            Centroid_labeled(ii,2:4) = Centroid(ii,1:3);
            Centroid_labeled(ii,5) = volume(ii,1);
        end
        
        % Creating a matrix containing the singlets only:
        
        if isempty(Au2_total) == 0
            Singlets_labeled = zeros(10000, 5);
            ll = 1;
            for ii = 1 : length(Centroid_labeled(:,1))
                if isempty(find([Au2_total{:}] == ii,1)) == 1
                    Singlets_labeled(ll, :) = Centroid_labeled(ii,:);
                    ll = ll + 1;
                end
            end
            Singlets_labeled(ll : end, :) = [];
        end
        % ---- STRUCTURE TO SAVE THE DATA --

        if save_multip_struct
            if isempty(Au_planar_for_breaks_count) == 0
                multip_struct.(num2str(sample)).planar.(['load_step_', num2str(oo)]).id = ...
                    Au_planar_for_breaks_count(:,1);
                multip_struct.(num2str(sample)).planar.(['load_step_', num2str(oo)]).coord = ...
                    cell2mat(Au_planar_for_breaks_count(:,2));
            end
            if isempty(Au2_total) == 0 %// if Au2 is not empty:
                multip_struct.(num2str(sample)).pl_diff.(['load_step_', num2str(oo)]) = Au2_total;
            end
            multip_struct.(num2str(sample)).centroids.(['load_step_', num2str(oo)]) = Centroid;
            
            %     % APRIL 2018: Adding the info on the singlets:
            if isempty(Au2_total) == 0
                multip_struct.(num2str(sample)).singlets.(['load_step_', num2str(oo)]) = Singlets_labeled;
            else
                multip_struct.(num2str(sample)).singlets.(['load_step_', num2str(oo)]) = Centroid_labeled;
            end
        end
        
        % A planar cluster is counted by the 
        % ImageJ code as one entity with one
        % centroid. Number of singlets making 
        % a cluster has to be considered in the
        % final singlets count:
        singlets_in_planar = sum(cellfun('length',Au_planar_for_breaks_count(:,1)));
        total_singlets = sing_count + singlets_in_planar;
        singlets(1, oo) = total_singlets;
        
        % % WRITE EXCEL SPREADSHEET:
        % warning( 'off', 'MATLAB:xlswrite:AddSheet' ) ; 
        % suppresses the warnings whena adding new worksheets
        if excel_file == 1
            xlwrite(results_file,{'Load step:'},sheet{oo},'A1');
            xlwrite(results_file,{sheet{oo}},sheet{oo},'B1');
            xlwrite(results_file,{'Singlets:'},sheet{oo},'A3');
            % Add Java POI Libs to matlab javapath
            if windows
                javaaddpath('poi_library\poi-3.8-20120326.jar');
                javaaddpath('poi_library\poi-ooxml-3.8-20120326.jar');
                javaaddpath('poi_library\poi-ooxml-schemas-3.8-20120326.jar');
                javaaddpath('poi_library\xmlbeans-2.3.0.jar');
                javaaddpath('poi_library\dom4j-1.6.1.jar');
                javaaddpath('poi_library\stax-api-1.0.1.jar');
            else
                javaaddpath('poi_library/poi-3.8-20120326.jar');
                javaaddpath('poi_library/poi-ooxml-3.8-20120326.jar');
                javaaddpath('poi_library/poi-ooxml-schemas-3.8-20120326.jar');
                javaaddpath('poi_library/xmlbeans-2.3.0.jar');
                javaaddpath('poi_library/dom4j-1.6.1.jar');
                javaaddpath('poi_library/stax-api-1.0.1.jar');
            end
            if isempty(Au2_total) == 0
                xlwrite(results_file,Singlets_labeled,sheet{oo},'A4');
            else
                xlwrite(results_file,Centroid_labeled,sheet{oo},'A4');
            end
            
            if isempty(Au2_total) == 0
                xlwrite(results_file,clusters_coord,sheet{oo},'G3');
            end
            xlwrite(results_file,{'Centroids:'},sheet{oo},'M3');
            xlwrite(results_file,length(Centroid(:,1)),sheet{oo},'N3');
            xlwrite(results_file,{'Extra singlets:'},sheet{oo},'M4');
            xlwrite(results_file,singlets_in_planar,sheet{oo},'N4');
            xlwrite(results_file,{'Tot singlets:'},sheet{oo},'M5');
            xlwrite(results_file,total_singlets,sheet{oo},'N5');
            xlwrite(results_file,{'Planar clusters:'},sheet{oo},'M6');
            xlwrite(results_file,length(Au2_planar(:,1)),sheet{oo},'N6');
            xlwrite(results_file,{'Total clusters:'},sheet{oo},'M7');
            xlwrite(results_file,length(clusters),sheet{oo},'N7');
            xlwrite(results_file,{'-PLETS (both diffuse and planar all together):'},sheet{oo},'L10');
            if clusters > 0
                for ii = 1 : length(A_ALL(:,1))
                    xlwrite(results_file,A_ALL{ii,1}(1),sheet{oo},['L' , int2str(ii+10)]);
                    xlwrite(results_file,A_ALL{ii,1}(2),sheet{oo},['M' , int2str(ii+10)]);
                    xlwrite(results_file,A_ALL{ii,2},sheet{oo},['N' , int2str(ii+10)]);
                end
            end
        end
        
        % To plot the clusters:
        if isempty(clusters) == 0
            for k = 1 : length(A_ALL(:,2))
                multiplets(k, oo) = A_ALL{k,2}(1);
            end
        end
        
        % To plot the diffuse clusters only:
        if isempty(diffuse_clusters) == 0
            for k = 1 : length(Au_diffuse_for_breaks_count(:,2))
                diffuse_multiplets(k, oo) = Au_diffuse_for_breaks_count{k,2}(1);
            end
        end
        % Just singlets evaluation:
        just_singlets_without_diffuse_and_planar = total_singlets - sum(clusters);
        just_singlets_without_planar = sing_count;
        % ----------------------
        % FEBRUARY 18:
        % Write down the number of planar clusters:
        
        if excel_file == 1
            xlwrite(results_file,{'Planar clusters SITES:'},sheet{oo},'X3');
            xlwrite(results_file,{'Planar 2-plets:'},sheet{oo},'X4'); % to write the 2-plets
            xlwrite(results_file,p_2,sheet{oo},'Y4');
            xlwrite(results_file,{'Planar 3-plets:'},sheet{oo},'X5'); % to write the 3-plets
            xlwrite(results_file,p_3,sheet{oo},'Y5');
            xlwrite(results_file,{'Planar 4-plets:'},sheet{oo},'X6'); % to write the 4-plets
            xlwrite(results_file,p_4,sheet{oo},'Y6');
            xlwrite(results_file,{'Planar 5-plets:'},sheet{oo},'X7'); % to write the 5-plets
            xlwrite(results_file,p_5,sheet{oo},'Y7');
            xlwrite(results_file,{'Planar 6-plets:'},sheet{oo},'X8'); % to write the 6-plets
            xlwrite(results_file,p_6,sheet{oo},'Y8');
            
            
            xlwrite(results_file,{'Singlets without breaks from diffuse & planar clusters:'},...
                sheet{oo},'X9'); % to write the smallest no. of singlets
            xlwrite(results_file,just_singlets_without_diffuse_and_planar,sheet{oo},'X10');
            
            
            xlwrite(results_file,{'Singlets without breaks from planar clusters only:'},sheet{oo},...
                'X11'); % to write the smallest no. of singlets
            xlwrite(results_file,just_singlets_without_planar,sheet{oo},'X12');
        end
        % Diffuse 2-plets and n-plets evaluation:
        diffuse_2_sites = 0;
        diffuse_n_sites = 0;
        diffuse_n_tot = 0;
        
        if isempty (Au_diffuse) == 0
            for k = 1 : length (Au_diffuse(:,1))
                if length (Au_diffuse{k}) == 2
                    diffuse_2_sites = diffuse_2_sites + 1;
                else
                    diffuse_n_tot = diffuse_n_tot + length(Au_diffuse{k});
                    diffuse_n_sites = diffuse_n_sites + 1;
                end
            end
        end
        
        if excel_file == 1
            xlwrite(results_file,{'Diffuse 2-plets sites:'},sheet{oo},'X13');
            xlwrite(results_file,diffuse_2_sites,sheet{oo},'X14');
            xlwrite(results_file,{'Diffuse n-plets sites:'},sheet{oo},'X15');
            xlwrite(results_file,diffuse_n_sites,sheet{oo},'X16');
            xlwrite(results_file,{'Diffuse n-plets total no. of breaks:'},sheet{oo},'X17');
            xlwrite(results_file,diffuse_n_tot,sheet{oo},'X18');
            
            xlwrite(results_file,{'Diffuse&plan 2-plets tot no.'},sheet{oo},'O10');
            if exist('A_ALL') == 1
                xlwrite(results_file,A_ALL{1,2}*2,sheet{oo},'O11');
            end
            % Write down the total no. of n-plets:
            xlwrite(results_file,{'Diffuse&plan n-plets tot no.'},sheet{oo},'O12');
        end
        tot_n_plets = 0;
        if exist('A_ALL') == 1
            for k = 2 : length(A_ALL(:,2))
                l = k + 1; tot_n_plets = tot_n_plets + A_ALL{k,2}*l;
            end
        end
        
        % The clusters counted as single blobs:
        plan_clusters_as_blobs = p_3 + p_4 + p_5 + p_6;
        plan_clusters_tot = p_3*3 + p_4*4 + p_5*5 + p_6*6;
        
        if excel_file == 1
            xlwrite(results_file,tot_n_plets,sheet{oo},'O13');
            xlwrite(results_file,{'Tot planar 2-plet'},sheet{oo},'Z4');
            xlwrite(results_file,p_2*2,sheet{oo},'AA4');
            xlwrite(results_file,{'Planar n-plet as blobs'},sheet{oo},'Z5');
            xlwrite(results_file,plan_clusters_as_blobs,sheet{oo},'AA5');
            xlwrite(results_file,{'Tot planar n-plet'},sheet{oo},'Z6');
            xlwrite(results_file,plan_clusters_tot,sheet{oo},'AA6');
            
            if isempty(Au_planar) == 0
                xlwrite(results_file,{'Planar (not in diffuse):'},sheet{oo},'Q3');
                for ii = 1 : length(Au2_planar(:,1))
                    xlwrite(results_file,Au2_planar{ii},sheet{oo},['Q' , int2str(ii+3)]);
                end
            end
        end
        
        % ---- CREATE TABLE WITH RESULTS FOR PLOTS ----
        breaks_table{1,1}  = 'singlets_abs';
        breaks_table{1,2}  = '2pl planar sites';
        breaks_table{1,3}  = '2pl planar tot';
        breaks_table{1,4}  = '2pl diff+plan site';
        breaks_table{1,5}  = '2pl diff+plan tot';
        breaks_table{1,6}  = 'n-pl planar sites';
        breaks_table{1,7}  = 'n-pl planar tot';
        breaks_table{1,8}  = 'n-pl diff+plan site';
        breaks_table{1,9}  = 'n-pl diff+plan tot';
        breaks_table{1,10}  = 'singlets density';
        breaks_table{1,11} = '2pl diff sites'; % added Mar-2018
        breaks_table{1,12} = 'n-pl diff sites'; % added Mar-2018
        
        breaks_table{1,13} = 'p_3'; % added Mar-2018
        breaks_table{1,14} = 'p_4'; % added Mar-2018
        breaks_table{1,15} = 'p_5'; % added Mar-2018
        breaks_table{1,16} = 'p_6'; % added Mar-2018
        
        breaks_table{oo+1,1}  = just_singlets_without_planar;
        breaks_table{oo+1,2}  = p_2;
        breaks_table{oo+1,3}  = p_2 * 2;
        if exist('A_ALL') == 1
            
            npl_diffplan_sites = 0;
            for ii = 2 : (length(A_ALL(:,2))-1)
                npl_diffplan_sites = npl_diffplan_sites + A_ALL{ii,2};
            end
            breaks_table{oo+1,8}  = npl_diffplan_sites;
            breaks_table{oo+1,4}  = A_ALL{1,2};
            breaks_table{oo+1,5}  = A_ALL{1,2} * 2;
            diffuse_2_sites = A_ALL{1,2} - p_2;
            breaks_table{oo+1, 11} = diffuse_2_sites;
            diffuse_n_sites = npl_diffplan_sites - plan_clusters_as_blobs;
            breaks_table{oo+1, 12} = diffuse_n_sites;
        else
            breaks_table{oo+1,4}  = 0;
            breaks_table{oo+1,5}  = 0;
            breaks_table{oo+1,8}  = 0;
            breaks_table{oo+1, 11} = 0;
            breaks_table{oo+1, 12} = 0;
        end
        
        
        breaks_table{oo+1,6}  = plan_clusters_as_blobs;
        breaks_table{oo+1,7}  = plan_clusters_tot;
        breaks_table{oo+1,9}  = tot_n_plets;
        
        breaks_table{oo+1,10} = just_singlets_without_planar / ...
            (x_scan * y_scan * z_scan * 0.65^3 / 10^9); %// Singlets density
        breaks_table{oo+1,11} = diffuse_2_sites;
        breaks_table{oo+1,12} = diffuse_n_sites;
        breaks_table{oo+1,13}  = p_3;
        breaks_table{oo+1,14}  = p_4;
        breaks_table{oo+1,15}  = p_5;
        breaks_table{oo+1,16}  = p_6;
        % ----- PLOTS -----

        % 1. Plot all the clusters (SCALED AXIS) with a different color:
        % PLOT ONLY LAST LOAD STEP:
        if plots_and_videos
            if oo == numel(sheet)
                % Mrkr = '+o*xsd^v><ph';
                Mrkr = '+oxs';
                figure('visible', 'off');
                for ii = 1 : length(Au2_total)
                    curr_col = rand(1,3);
                    curr_marker = Mrkr(mod(ii,length(Mrkr))+1); 
                    % modulus, gives a number that goes from 
                    % ii to the length of Mrkr and then starts 
                    % again so if length is 4, no.5 will be 1 again!
                    for jj = 1 : length(Au2_total{ii})
                        scatter3(scaled_Centroid(Au2_total{ii}(jj),1),...
                            scaled_Centroid(Au2_total{ii}(jj),2),...
                            scaled_Centroid(Au2_total{ii}(jj),3), 50, ...
                            'Marker',curr_marker, 'MarkerFaceColor', curr_col, 'MarkerEdgeColor', curr_col)
                        hold on
                        rotate3d on
                    end
                end
                view(110, 15); xlabel('(µm)'); ylabel('(µm)'); zlabel('(µm)');
                title(sprintf('Fibre break clusters before failure - %s', sample_name_on_title))
                xlim([0 x_scan*voxel]); ylim([0 y_scan*voxel]); zlim([0 z_scan*voxel]);
                % 2. Save the plots in a subfolder in Matlab
                % SAVE .PNG:
                clusters_plot = sprintf('%s%s_%s_clusters.png', newSubFolder, sample, sheet{oo});
                saveas(gcf, clusters_plot);
                % SAVE .FIG:
                clusters_plot = sprintf('%s%s_%s_clusters.fig', newSubFolder, sample, sheet{oo});
                saveas(gcf, clusters_plot);
                close
            end
            
            % VIDEO ACCUMULATION OF SINGLETS & n-PLETS
            if isempty(Au2_total) == 0
                % Counting the no. of elements in Au2_total
                l = 0;
                k = 0;
                for ii = 1 : length(Au2_total)
                    k = numel(Au2_total{ii});
                    l = l + k;
                    k = 0;
                end
                kk = 1;
                Clusters_labeled = zeros(l, 5);
                clear l
                ll = 1;
                for ii = 1 : length(Au2_total)
                    for jj = 1 : length(Au2_total{ii})
                        Clusters_labeled(kk, 1) = Au2_total{ii}(jj);
                        Clusters_labeled(kk, 2:4) = scaled_Centroid(Au2_total{ii}(jj),:);
                        Clusters_labeled(kk, 5) = ll;
                        kk = kk + 1;
                    end
                    ll = ll + 1;
                end
            else
                Singlets_labeled(:,:) = scaled_Centroid; 
                % Singlets_labeled is scaled as can be seen here(!)
            end
            
            figure%('visible', 'off');
            if isempty(Au2_total) == 1
                scatter3(Singlets_labeled(:,1),Singlets_labeled(:,2),Singlets_labeled(:,3),...
                    30, 'r', 'fill')
            else
                scatter3(Singlets_labeled(:,2)*voxel, Singlets_labeled(:,3)*voxel,...
                    Singlets_labeled(:,4)*voxel, 30, 'r', 'fill')
                hold on
                ii = 1;
                for jj = 1 : ll
                    Mrkr = '+oxs';
                    curr_marker = Mrkr(mod(jj,length(Mrkr))+1);
                    colours = colours(mod(jj,length(colours(:,1)))+1,:);
                    while (ii <= length(Clusters_labeled(:,1))) && (Clusters_labeled(ii, 5) == jj)
                        scatter3(Clusters_labeled(ii,2), Clusters_labeled(ii,3),...
                            Clusters_labeled(ii,4), 'b', 'Marker', curr_marker, 'MarkerEdgeColor', colours)
                        ii = ii + 1;
                    end
                end
            end
            xlim([0 x_scan*voxel]); ylim([0 y_scan*voxel]); zlim([0 z_scan*voxel]);
            rotate3d on
            view(90, 0); %// Perspective view with 110, 15
            xlabel('(µm)'); ylabel('(µm)'); zlabel('(µm)');
            width_c = scan_data{16,sample_no}; % (µm)
            thickness_c = scan_data{17,sample_no}; % (µm)
            area_min = width_c * thickness_c / 10^6;
            title(sprintf('Ply-level stress: %0.f MPa', scan_load{2,sample_no-1}(oo) / area_min))
            legend('Singlets', 'N-plets', 'Location', 'NorthWest')
            legend('boxoff')
            hold on
            frame = getframe(gcf);
            writeVideo(writerObj, frame);
        end
        
        clearvars -except scan_data excel_file plots_and_videos windows file_to_read sheet sample ...
            oo singlets breaks_table A_ALL multiplets sample_no newSubFolder ...
            writerObj writerObj.FrameRate sample_name_on_title colours breaks_table ...
            multip_struct diffuse_multiplets scan_load save_flag ...
            sam samples samples_no sample_names_on_title save_multip_struct
        
        warning( 'off', 'MATLAB:xlswrite:AddSheet' ) ; 
        % suppresses the warnings whena adding new worksheets
        
    end

    if plots_and_videos
        % Closing the movie:
        hold off
        close(writerObj); % Saves the movie.
    end
    
    % Saving the workspace
    if save_flag
        if save_multip_struct
            save ('D:\OneDrive - University of Southampton\multip_struct.mat', 'multip_struct')
        end
        workspace = strcat('workspace_', scan_data{1,sample_no}, '.mat');
        % This clear because accum_plot loads 
        % separately scan_data without taking it
        % from this workspace. Otherwise if a modification is made
        % on this script it will affect the accum_plot:
        clear multip_struct scan_load scan_data
        SaveFolder = "D:\OneDrive - University of Southampton\MATLAB\workspaces_break_accumulation\";
        save(sprintf('%s%s', SaveFolder, workspace))
    end
end