Update microstructure characterization module to v0.91

Author: Usseglio-Viretta

src/Microstructure_characterization/Function_Connectivity_Algorithm.m

@@ -16,7 +16,7 @@
 unique_cluster = unique(Clusters_tmp); % Unique cluster
 unique_cluster(unique_cluster==0)=[]; % Remove background (0) if any      
 number_cluster = length(unique_cluster); % Number of cluster
-cluster_size = zeros(number_cluster,12); % Initialise
+cluster_size = zeros(number_cluster,17); % Initialise
 % Column
 % 1 : cluster number of voxel
 % 2 : cluster volume in um3
@@ -29,6 +29,13 @@
 % 9 : 1= connected cluster, 2=unknown cluster, 3=isolated cluster (face2face connection direction 1)
 % 10: 1= connected cluster, 2=unknown cluster, 3=isolated cluster (face2face connection direction 2)
 % 11: 1= connected cluster, 2=unknown cluster, 3=isolated cluster (face2face connection direction 3)
+% 12 : 1= connected cluster, 2=unknown cluster, 3=isolated cluster (face connection 1)
+% 13: 1= connected cluster, 2=unknown cluster, 3=isolated cluster (face connection 1)
+% 14: 1= connected cluster, 2=unknown cluster, 3=isolated cluster (face connection 1)
+% 15 : 1= connected cluster, 2=unknown cluster, 3=isolated cluster (face2face connection direction 1)
+% 16: 1= connected cluster, 2=unknown cluster, 3=isolated cluster (face2face connection direction 2)
+% 17: 1= connected cluster, 2=unknown cluster, 3=isolated cluster (face2face connection direction 3)
+
 Connectivity_structure.connected_id=connected_id;
 Connectivity_structure.unknown_id=unknown_id;
 Connectivity_structure.isolated_id=isolated_id;
@@ -79,6 +86,8 @@
 % Cluster does not connect the faces, and do not touch the domain's boundaryies: =3
 for current_direction=1:1:3 % Loop on every direction
     Clusters_TransportFace2Face.direction(current_direction).array = zeros(Domain_size,'uint8'); % Initialise
+    Clusters_TransportFromFace1.direction(current_direction).array = zeros(Domain_size,'uint8'); % Initialise
+    Clusters_TransportFromFace2.direction(current_direction).array = zeros(Domain_size,'uint8'); % Initialise
 end
 
 Clusters_LargestIsolatedUnknown.array(Clusters_==1)=connected_id; % Main (largest) cluster
@@ -139,6 +148,77 @@
     end
     Clusters_TransportFace2Face.direction(3).array(I) = connectivity_id_3;
     cluster_size(k,11)=connectivity_id_3;
+
+
+    % % Check if cluster connect to one face
+    % Along direction 1
+    if x_min==1
+        connectivity_id_1 = connected_id; % Connect
+    elseif y_min==1 || z_min==1 || y_max==Domain_size(2) || z_max==Domain_size(3)
+        connectivity_id_1 = unknown_id; % Does not connect, but connection actually unknown due to limited field of view
+    else
+        connectivity_id_1 = isolated_id; % Does not connect (true isolated cluster)
+    end
+    Clusters_TransportFromFace1.direction(1).array(I) = connectivity_id_1;
+    cluster_size(k,12)=connectivity_id_1;
+    
+    % Along direction 2
+    if y_min==1
+        connectivity_id_2 = 1; % Connect
+    elseif x_min==1 || z_min==1 || x_max==Domain_size(1) || z_max==Domain_size(3)
+        connectivity_id_2 = unknown_id; % Does not connect, but connection actually unknown due to limited field of view
+    else
+        connectivity_id_2 = isolated_id; % Does not connect (true isolated cluster)
+    end
+    Clusters_TransportFromFace1.direction(2).array(I) = connectivity_id_2;
+    cluster_size(k,13)=connectivity_id_2;
+    
+    % Along direction 3
+    if z_min==1
+        connectivity_id_3 = 1; % Does not connect (true isolated cluster)
+    elseif y_min==1 || x_min==1 || y_max==Domain_size(2) || x_max==Domain_size(1)
+        connectivity_id_3 = unknown_id; % Does not connect (true isolated cluster)
+    else
+        connectivity_id_3 = isolated_id; % Does not connect (true isolated cluster)
+    end
+    Clusters_TransportFromFace1.direction(3).array(I) = connectivity_id_3;
+    cluster_size(k,14)=connectivity_id_3;
+
+    % % Check if cluster connect to one face
+    % Along direction 1
+    if x_max==Domain_size(1)
+        connectivity_id_1 = connected_id; % Connect
+    elseif y_min==1 || z_min==1 || y_max==Domain_size(2) || z_max==Domain_size(3)
+        connectivity_id_1 = unknown_id; % Does not connect, but connection actually unknown due to limited field of view
+    else
+        connectivity_id_1 = isolated_id; % Does not connect (true isolated cluster)
+    end
+    Clusters_TransportFromFace2.direction(1).array(I) = connectivity_id_1;
+    cluster_size(k,15)=connectivity_id_1;
+    
+    % Along direction 2
+    if y_max==Domain_size(2)
+        connectivity_id_2 = 1; % Connect
+    elseif x_min==1 || z_min==1 || x_max==Domain_size(1) || z_max==Domain_size(3)
+        connectivity_id_2 = unknown_id; % Does not connect, but connection actually unknown due to limited field of view
+    else
+        connectivity_id_2 = isolated_id; % Does not connect (true isolated cluster)
+    end
+    Clusters_TransportFromFace2.direction(2).array(I) = connectivity_id_2;
+    cluster_size(k,16)=connectivity_id_2;
+    
+    % Along direction 3
+    if z_max==Domain_size(3)
+        connectivity_id_3 = 1; % Does not connect (true isolated cluster)
+    elseif y_min==1 || x_min==1 || y_max==Domain_size(2) || x_max==Domain_size(1)
+        connectivity_id_3 = unknown_id; % Does not connect (true isolated cluster)
+    else
+        connectivity_id_3 = isolated_id; % Does not connect (true isolated cluster)
+    end
+    Clusters_TransportFromFace2.direction(3).array(I) = connectivity_id_3;
+    cluster_size(k,17)=connectivity_id_3;
+
+
 end
 
 % Save results in structure
@@ -156,10 +236,20 @@
     Clusters_TransportFace2Face.direction(current_direction).connected_cluster_phasefraction = 100 * sum(sum(sum( Clusters_TransportFace2Face.direction(current_direction).array==connected_id))) / Numbervoxel_phase; 
     Clusters_TransportFace2Face.direction(current_direction).unknown_cluster_phasefraction   = 100 * sum(sum(sum( Clusters_TransportFace2Face.direction(current_direction).array==unknown_id))) / Numbervoxel_phase; 
     Clusters_TransportFace2Face.direction(current_direction).isolated_cluster_phasefraction  = 100 * sum(sum(sum( Clusters_TransportFace2Face.direction(current_direction).array==isolated_id))) / Numbervoxel_phase; 
+
+    Clusters_TransportFromFace1.direction(current_direction).connected_cluster_phasefraction = 100 * sum(sum(sum( Clusters_TransportFromFace1.direction(current_direction).array==connected_id))) / Numbervoxel_phase; 
+    Clusters_TransportFromFace1.direction(current_direction).unknown_cluster_phasefraction   = 100 * sum(sum(sum( Clusters_TransportFromFace1.direction(current_direction).array==unknown_id))) / Numbervoxel_phase; 
+    Clusters_TransportFromFace1.direction(current_direction).isolated_cluster_phasefraction  = 100 * sum(sum(sum( Clusters_TransportFromFace1.direction(current_direction).array==isolated_id))) / Numbervoxel_phase; 
+
+    Clusters_TransportFromFace2.direction(current_direction).connected_cluster_phasefraction = 100 * sum(sum(sum( Clusters_TransportFromFace2.direction(current_direction).array==connected_id))) / Numbervoxel_phase; 
+    Clusters_TransportFromFace2.direction(current_direction).unknown_cluster_phasefraction   = 100 * sum(sum(sum( Clusters_TransportFromFace2.direction(current_direction).array==unknown_id))) / Numbervoxel_phase; 
+    Clusters_TransportFromFace2.direction(current_direction).isolated_cluster_phasefraction  = 100 * sum(sum(sum( Clusters_TransportFromFace2.direction(current_direction).array==isolated_id))) / Numbervoxel_phase; 
+
 end
 % Save results in structure
 Connectivity_structure.Clusters_LargestIsolatedUnknown = Clusters_LargestIsolatedUnknown;
 Connectivity_structure.Clusters_TransportFace2Face = Clusters_TransportFace2Face;
-
+Connectivity_structure.Clusters_TransportFromFace1 = Clusters_TransportFromFace1;
+Connectivity_structure.Clusters_TransportFromFace2 = Clusters_TransportFromFace2;
 
 end

src/Microstructure_characterization/Function_Discrete_particle_size_PCRF.m

@@ -144,6 +144,10 @@
 density_fct_parameters.round_value = 3;
 density_fct_parameters.smooth_cumulative_fct = true;
 
+% Edge detection
+background = 0;
+edgewithbackground = false;
+
 for current_phase=1:1:number_phase % Loop over all phases
     fprintf ('Current phase calculated: %s\n',INFO.phase(current_phase).name)
     code_tmp = INFO.phase(current_phase).code;
@@ -155,7 +159,7 @@
 
     % Edge, watershed
     [index_border_phase,~,~,~] = Function_identify_edges(binary_phase);
-    [index_border_label,~,~,~] = Function_identify_labelsedges(Particle_label_tmp, 0); 
+    [index_border_label,~,~,~] = Function_identify_labelsedges(Particle_label_tmp, background, edgewithbackground); 
     index_(current_phase).border_phase = index_border_phase;
     index_(current_phase).border_label = index_border_label;
 

src/Microstructure_characterization/Function_Discrete_particle_size_watershed_immersion.m

@@ -142,6 +142,10 @@
 density_fct_parameters.round_value = 3;
 density_fct_parameters.smooth_cumulative_fct = true;
 
+% Edge detection
+background = 0;
+edgewithbackground = false;
+
 for current_phase=1:1:number_phase % Loop over all phases
     fprintf ('Current phase calculated: %s\n',INFO.phase(current_phase).name)
     code_tmp = INFO.phase(current_phase).code;
@@ -153,7 +157,7 @@
 
     % Edge, watershed
     [index_border_phase,~,~,~] = Function_identify_edges(binary_phase);
-    [index_border_label,~,~,~] = Function_identify_labelsedges(Particle_label_tmp, 0); 
+    [index_border_label,~,~,~] = Function_identify_labelsedges(Particle_label_tmp, background, edgewithbackground); 
     index_(current_phase).border_phase = index_border_phase;
     index_(current_phase).border_label = index_border_label;
 

src/Microstructure_characterization/Function_Discrete_particle_size_watershed_immersion_algorithm.m

@@ -329,6 +329,9 @@
 if visualize_2D
     close(video_handle) % Close video   
 
+    background = 0;
+    edgewithbackground = false;
+    
     sub_axes_illustration(3).s=subplot(2,2,3,'Parent',Fig_illustration); % Create axes  
     hold(sub_axes_illustration(3).s,'on'); % Active subplot
     slice_image = image(slice_color,'parent',sub_axes_illustration(3).s); % Display the slice
@@ -343,7 +346,7 @@
     hold(sub_axes_illustration(3).s,'off'); % Active subplot    
 
     % Find watershedlines
-    [index_border_label,~,~,~] = Function_identify_labelsedges(Label_lake, 0);
+    [index_border_label,~,~,~] = Function_identify_labelsedges(Label_lake, background, edgewithbackground);
     sub_axes_illustration(4).s=subplot(2,2,4,'Parent',Fig_illustration); % Create axes  
     hold(sub_axes_illustration(4).s,'on'); % Active subplot
     tmp = slice_grey; tmp(index_border_label) = 1;    

src/Microstructure_characterization/Function_along_direction.m

@@ -27,32 +27,39 @@
             for k=1:1:length(p.ignore_value)
                 different_values(different_values==p.ignore_value(k)) = [];
             end
-                     
-            % Minimum
-            Alongdirection.phase(current_phase).direction(current_direction).val(position,2) = min(different_values);
-                        
-            % Maximum
-            Alongdirection.phase(current_phase).direction(current_direction).val(position,4) = max(different_values);
 
-            % Weighed values
-            % Initialisation
-            values=zeros(length(different_values),2);
-            for current_val=1:1:length(different_values)
-                values(current_val,1)=different_values(current_val);
-                values(current_val,2)=sum(sum( slice_== different_values(current_val)));
+            if ~isempty(different_values)
+                % Minimum
+                Alongdirection.phase(current_phase).direction(current_direction).val(position,2) = min(different_values);
+                
+                % Maximum
+                Alongdirection.phase(current_phase).direction(current_direction).val(position,4) = max(different_values);
+                
+                % Weighed values
+                % Initialisation
+                values=zeros(length(different_values),2);
+                for current_val=1:1:length(different_values)
+                    values(current_val,1)=different_values(current_val);
+                    values(current_val,2)=sum(sum( slice_== different_values(current_val)));
+                end
+                % Mean
+                Alongdirection.phase(current_phase).direction(current_direction).val(position,3) = sum(values(:,1).*values(:,2))/sum(values(:,2));
+                
+                % Standard deviation
+                % The weighted starndard deviation formula is
+                % sqrt( sum(wi((xi-<x>)^2)  / ( (n-1)/n * sum wi ) )
+                % With wi the weight of the xi, and <x> the weighted mean (mean_size)
+                wi = values(:,2);
+                xi = values(:,1);
+                n = length(xi);
+                mean_ = Alongdirection.phase(current_phase).direction(current_direction).val(position,3);
+                Alongdirection.phase(current_phase).direction(current_direction).val(position,5) = sqrt( sum( wi.*((xi-mean_).^2)) / ( (n-1)/n * sum(wi)  ));
+            else
+                Alongdirection.phase(current_phase).direction(current_direction).val(position,2) = NaN;
+                Alongdirection.phase(current_phase).direction(current_direction).val(position,3) = NaN;
+                Alongdirection.phase(current_phase).direction(current_direction).val(position,4) = NaN;
+                Alongdirection.phase(current_phase).direction(current_direction).val(position,5) = NaN;
             end
-            % Mean
-            Alongdirection.phase(current_phase).direction(current_direction).val(position,3) = sum(values(:,1).*values(:,2))/sum(values(:,2));             
-            
-            % Standard deviation
-            % The weighted starndard deviation formula is
-            % sqrt( sum(wi((xi-<x>)^2)  / ( (n-1)/n * sum wi ) )
-            % With wi the weight of the xi, and <x> the weighted mean (mean_size)
-            wi = values(:,2);
-            xi = values(:,1);
-            n = length(xi);
-            mean_ = Alongdirection.phase(current_phase).direction(current_direction).val(position,3);
-            Alongdirection.phase(current_phase).direction(current_direction).val(position,5) = sqrt( sum( wi.*((xi-mean_).^2)) / ( (n-1)/n * sum(wi)  ));
 
         end
     end