some corrections on features and functions

brick/processing/best_fit_ellipsoid.py

@@ -104,20 +104,21 @@ def polyToParams3D(vec, printMe):  # gets 3D parameters of an ellipsoid.
     return (center, axes, inve)
 
 
-def GetConvexHull3D(soma_sites: site_h) -> array_t:
+def GetConvexHull3D(soma_sites: site_h) -> tuple:
     #
-    volume = np.stack((soma_sites[0], soma_sites[1], soma_sites[2]), axis=-1)
-    hull_volume = ConvexHull(volume)
+    three_D = np.stack((soma_sites[0], soma_sites[1], soma_sites[2]), axis=-1)
+    hull_three_D = ConvexHull(three_D)
+    volume_of_CH = hull_three_D.volume
 
-    len_hull = len(hull_volume.vertices)
+    len_hull = len(hull_three_D.vertices)
     hull = np.zeros((len_hull, 3))
 
-    for i in range(len(hull_volume.vertices)):
-        hull[i] = volume[hull_volume.vertices[i]]
+    for i in range(len(hull_three_D.vertices)):
+        hull[i] = three_D[hull_three_D.vertices[i]]
 
     convex_hull = np.transpose(hull)
 
-    return convex_hull
+    return convex_hull, volume_of_CH
 
 
 def FindBestFittingEllipsoid3D(soma: soma_t) -> tuple:
@@ -127,7 +128,7 @@ def FindBestFittingEllipsoid3D(soma: soma_t) -> tuple:
     """
 
     # get convex hull
-    convex_hull = GetConvexHull3D(soma.sites)
+    convex_hull = GetConvexHull3D(soma.sites)[0]
 
     # fit ellipsoid on the convex hull
     # # get ellipsoid polynomial coefficients

nutrimorph.py

@@ -119,7 +119,7 @@ image = in_.ImageVerification(image, channel)
 # iv_.image_verification(image, channel)  # -> PySide2 user interface  # TODO: must return the modified image!
 # /!\ conflicts between some versions of PySide2 and Python3
 
-image = image[:, 512:, 512:]  # 512 # 562 # Just for development
+image = image[:, 800:, 800:]  # 512 # 562 # Just for development
 img_shape = image.shape
 
 #
@@ -173,7 +173,7 @@ print("\n--- Extension Detection")
 
 # Change the extensions parameters from micron to pixel dimensions
 ext_min_area_c = in_.ToPixel(ext_min_area_c, size_voxel_in_micron, dimension=(0, 1))
-ext_selem_pixel_c = mp_.disk(in_.ToPixel(ext_selem_micron_c, size_voxel_in_micron))
+ext_selem_pixel_c = mp_.disk(in_.ToPixel(ext_selem_micron_c, size_voxel_in_micron))  # TODO if 0 : do not do the Cleaning => None et tester dans Cleaning
 scale_range_pixel = []
 #
 for value in scale_range:
@@ -203,12 +203,14 @@ elapsed_time = tm_.gmtime(tm_.time() - start_time)
 print(f"Elapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}\n")
 
 # Creation of the enhanced maps
-ext_nfo["coarse_map"] = extension_t.CoarseMap(enhanced_ext, ext_low_c, ext_high_c, ext_selem_pixel_c)
-ext_nfo["coarse_map"], ext_lmp = extension_t.FilteredCoarseMap(ext_nfo["coarse_map"], ext_min_area_c)
-ext_nfo["map"] = extension_t.FineMapFromCoarseMap(ext_nfo["coarse_map"])
+ext_nfo["coarse_map"] = extension_t.CoarseMap(enhanced_ext, ext_low_c, ext_high_c, ext_selem_pixel_c) # seuillage
+ext_nfo["coarse_map"], ext_lmp = extension_t.FilteredCoarseMap(ext_nfo["coarse_map"], ext_min_area_c)  # min
+ext_nfo["map"] = extension_t.FineMapFromCoarseMap(ext_nfo["coarse_map"]) # skeleton
 ext_nfo["map"][som_nfo["map"] > 0] = 0
 ext_nfo["lmp"], n_extensions = ms_.label(ext_nfo["map"], return_num=True)
 
+# TODO = regarder le MIP ici ! Visionner les 15 images en difference !
+
 # Use relabel instead of label to optimize the algorithm. BUT PROBLEM WITH THE NUMBER OF EXTENSIONS DETECTED !
 # ext_nfo["lmp"] = relabel_sequential(ext_lmp)[0]
 # n_extensions = ext_nfo["lmp"].max()
@@ -360,11 +362,6 @@ print('\n--- Graph extraction')
 print('\n- Graph roots')
 
 # Create the graphs
-# # Initialization of overall somas measures
-all_ext_lengths = []
-P_ext_lengths = []
-S_ext_lengths = []
-
 for soma in somas:
     ext_map = skl_map_t.FromShapeMap(ext_nfo['lmp_soma'] == soma.uid, store_widths=True, skeletonize=False)  # do_post_thinning=True
     # to remove pixel that are not breaking connectivity - FineMap in FromShapeMap
@@ -404,16 +401,22 @@ for soma in somas:
     # Soma features
     print('***Soma***')
     # # Volume of the soma
-    volume_pixel_micron = round(mt_.prod(size_voxel_in_micron[axis] for axis in (0, 1, 2)), 4)
+    volume_pixel_micron = round(np_.prod(size_voxel_in_micron), 4)
     soma.volume_soma_micron = volume_pixel_micron * len(soma.sites[0])
-    print("Soma volume = ", soma.volume_soma_micron)
+    volume_convex_hull = volume_pixel_micron * bf_.GetConvexHull3D(soma.sites)[1]
+    Coef_V_soma__V_convex_hull = soma.volume_soma_micron / volume_convex_hull  ## TODO keep as prm
+    print(f"Volume soma = {soma.volume_soma_micron}\n"
+          f"Volume soma / Volume Convex Hull = {Coef_V_soma__V_convex_hull}"
+          )
 
     # # Axes of the best fitting ellipsoid
     soma.axes_ellipsoid = bf_.FindBestFittingEllipsoid3D(soma)[2]
-    print("Axes of best fitting ellipsoid : ", soma.axes_ellipsoid)
+    Coef_axes_ellips_b__a = soma.axes_ellipsoid[1] / soma.axes_ellipsoid[0]  ## TODO keep as prm (2)
+    Coef_axes_ellips_c__a = soma.axes_ellipsoid[2] / soma.axes_ellipsoid[0]
+
 
     # -- Extension features
-    # # Graph features
+    # # Graph features # TODO keep as prms (4)
     N_nodes = soma.skl_graph.n_nodes  # number of nodes
     N_ext = soma.skl_graph.n_edges - len(soma.graph_roots)  # number of edges except the constructed ones from node soma to the roots
     N_primary_ext = len(soma.graph_roots)  # number of primary edges = linked to the soma except the constructed ones from node soma to the roots
@@ -427,23 +430,33 @@ for soma in somas:
         f"N secondary extensions = {N_sec_ext}\n"
     )
 
-    if len(soma.graph_roots) > 0:
-        highest_degree = soma.skl_graph.max_degree  # highest degree of the nodes except the soma
-        highest_degree_w_node = soma.skl_graph.highest_degree_w_nodes  # highest degree of the nodes with the node coordinates except the soma
-        N_nodes_of_highest_degree = len(highest_degree_w_node[1])  # number of nodes of highest degree
-        # min, mean, median, max and standard deviation of the degrees of non-leaves nodes
-        min_degree = soma.skl_graph.min_degree_except_leaves
-        mean_degree = soma.skl_graph.mean_degree_except_leaves
-        median_degree = soma.skl_graph.median_degree_except_leaves
-        max_degree = soma.skl_graph.max_degree_except_leaves
-        std_degree = soma.skl_graph.std_degree_except_leaves
-        #
+    if N_primary_ext > 0:
+        if N_sec_ext == 0:
+            highest_degree = 1
+            highest_degree_w_node = soma.skl_graph.highest_degree_w_nodes(soma)  # highest degree of the nodes with the node coordinates except the soma
+            # min, mean, median, max and standard deviation of the degrees of non-leaves nodes
+            min_degree = 1
+            mean_degree = 1
+            median_degree = 1
+            max_degree = 1
+            std_degree = 0
+
+        elif N_sec_ext > 0:
+            highest_degree = soma.skl_graph.max_degree  # highest degree of the nodes except the soma
+            if highest_degree == 2:
+                highest_degree = 1
+            highest_degree_w_node = soma.skl_graph.highest_degree_w_nodes(soma)  # highest degree of the nodes with the node coordinates except the soma
+            # min, mean, median, max and standard deviation of the degrees of non-leaves nodes
+            min_degree = soma.skl_graph.min_degree_except_leaves_and_roots
+            mean_degree = soma.skl_graph.mean_degree_except_leaves_and_roots
+            median_degree = soma.skl_graph.median_degree_except_leaves_and_roots
+            max_degree = soma.skl_graph.max_degree_except_leaves_an_roots
+            std_degree = soma.skl_graph.std_degree_except_leaves_and_roots
+
         # Calculate the extensions lengths
         ext_lengths = soma.skl_graph.edge_lengths
         total_ext_length = soma.skl_graph.length
         #
-        all_ext_lengths += list(ext_lengths)
-        #
         # min, mean, median, max and standard deviation of the ALL extensions
         min_length = soma.skl_graph.min_length
         mean_length = soma.skl_graph.mean_length
@@ -455,7 +468,6 @@ for soma in somas:
         print(
             f"NODES DEGREES\n"
             f"Highest degree (except soma) = {highest_degree}/{highest_degree_w_node}\n"
-            f"N nodes with highest degree = {N_nodes_of_highest_degree}\n"
             f"Min/Mean/Median/Max degree (except soma & leaves) = {min_degree} / {mean_degree} / {median_degree} / {max_degree}\n"
             f"Standard deviation (except soma & leaves) = {std_degree}\n\n"
             #
@@ -463,20 +475,18 @@ for soma in somas:
             f"  Min/Mean/Median/Max Length = {min_length} / {mean_length} / {median_length} / {max_length}\n"
             f"  Standard Deviation = {std_length} / Entropy = {entropy_length}")
 
-        pl_.hist(ext_lengths, color='r')
-        pl_.title(f"Histogram of all the extensions lengths of soma {soma.uid}")
-        pl_.xlabel("Lengths")
-        pl_.ylabel('Number of extensions')
-        pl_.savefig(f"path/hist_all_ext_soma_{soma.uid}.png")
-        # pl_.show()
-        pl_.close()
+        # pl_.hist(ext_lengths, color='r')
+        # pl_.title(f"Histogram of all the extensions lengths of soma {soma.uid}")
+        # pl_.xlabel("Lengths")
+        # pl_.ylabel('Number of extensions')
+        # pl_.savefig(f"path/hist_all_ext_soma_{soma.uid}.png")
+        # # pl_.show()
+        # pl_.close()
 
         # PRIMARY extensions
         ext_lengths_P = soma.skl_graph.primary_edge_lengths(soma)
         total_ext_length_P = sum(ext_lengths_P)
         #
-        P_ext_lengths += list(ext_lengths_P)
-        #
         # min, mean, median, max and standard deviation of the PRIMARY extensions
         if total_ext_length_P > 0:
             min_length_P = min(ext_lengths_P)
@@ -491,20 +501,18 @@ for soma in somas:
                 f"  Min/Mean/Median/Max Length = {min_length_P} / {mean_length_P} / {median_length_P} / {max_length_P}\n"
                 f"  Standard Deviation = {std_length_P} / Entropy = {entropy_length_P}")
 
-            pl_.hist(ext_lengths_P, color='b')
-            pl_.title(f"Histogram of Primary extensions lengths of soma {soma.uid}")
-            pl_.xlabel("Lengths")
-            pl_.ylabel('Number of Primary extensions')
-            pl_.savefig(f"path/hist_P_ext_soma_{soma.uid}.png")
-            # pl_.show(block=True)
-            pl_.close()
+            # pl_.hist(ext_lengths_P, color='b')
+            # pl_.title(f"Histogram of Primary extensions lengths of soma {soma.uid}")
+            # pl_.xlabel("Lengths")
+            # pl_.ylabel('Number of Primary extensions')
+            # pl_.savefig(f"path/hist_P_ext_soma_{soma.uid}.png")
+            # # pl_.show(block=True)
+            # pl_.close()
 
         # SECONDARY extensions
         ext_lengths_S = soma.skl_graph.secondary_edge_lengths(soma)
         total_ext_length_S = sum(ext_lengths_S)
         #
-        S_ext_lengths += list(ext_lengths_S)
-        #
         # min, mean, median, max and standard deviation of the PRIMARY extensions
         if total_ext_length_S > 0:
             min_length_S = min(ext_lengths_S)
@@ -519,26 +527,13 @@ for soma in somas:
                 f"  Min/Mean/Median/Max Length = {min_length_S} / {mean_length_S} / {median_length_S} / {max_length_S}\n"
                 f"  Standard Deviation = {std_length_S} / Entropy = {entropy_length_S}")
 
-            pl_.hist(ext_lengths_S, color='g')
-            pl_.title(f"Histogram of Secondary extensions lengths of soma {soma.uid}")
-            pl_.xlabel("Lengths")
-            pl_.ylabel('Number of Secondary extensions')
-            pl_.savefig(f"path/hist_S_ext_soma_{soma.uid}.png")
-            # pl_.show(block=True)
-            pl_.close()
-
-all_ext_lengths = tuple(all_ext_lengths)
-P_ext_lengths = tuple(P_ext_lengths)
-S_ext_lengths = tuple(S_ext_lengths)
-
-hist_all = pl_.hist(all_ext_lengths, label="All extensions", alpha=0.7, color='r')
-hist_P = pl_.hist(P_ext_lengths, label="Primary extensions", alpha=0.5, color='b')
-pl_.legend()
-pl_.title("Histogram of extensions lengths")
-pl_.xlabel("Lengths")
-pl_.ylabel('Number of extensions')
-pl_.savefig("path/hist_ext_all_soma.png")
-pl_.show(block=True)
+            # pl_.hist(ext_lengths_S, color='g')
+            # pl_.title(f"Histogram of Secondary extensions lengths of soma {soma.uid}")
+            # pl_.xlabel("Lengths")
+            # pl_.ylabel('Number of Secondary extensions')
+            # pl_.savefig(f"path/hist_S_ext_soma_{soma.uid}.png")
+            # # pl_.show(block=True)
+            # pl_.close()
 
 
     # TODO FRANGI:WIDTH OF EXTENSIONS, CURVATURE EXTENSIONS

sklgraph/skl_fgraph.py

@@ -57,39 +57,39 @@ class skl_graph_t(skl_nfgraph_t):
         return np_.std(list(degree for node, degree in self.degree if "S" not in node and degree != 1))
 
     @property
-    def max_degree_except_leaves(self) -> int:
-        return max(degree for node, degree in self.degree if "S" not in node and degree != 1)
+    def max_degree_except_leaves_an_roots(self) -> int:
+        try:
+            return max(degree for node, degree in self.degree if "S" not in node and degree > 2)
+        except:
+            return 1
 
     @property
-    def min_degree_except_leaves(self) -> int:
-        return min(degree for node, degree in self.degree if "S" not in node and degree != 1)
+    def min_degree_except_leaves_and_roots(self) -> int:
+        try:
+            return min(degree for node, degree in self.degree if "S" not in node and degree > 2)
+        except:
+            return 1
 
     @property
-    def mean_degree_except_leaves(self) -> int:
-        return np_.mean(list(degree for node, degree in self.degree if "S" not in node and degree != 1))
+    def mean_degree_except_leaves_and_roots(self) -> int:
+        try:
+            return np_.mean(list(degree for node, degree in self.degree if "S" not in node and degree > 2))
+        except:
+            return 1
 
     @property
-    def median_degree_except_leaves(self) -> int:
-        return np_.median(list(degree for node, degree in self.degree if "S" not in node and degree != 1))
+    def median_degree_except_leaves_and_roots(self) -> int:
+        try:
+            return np_.median(list(degree for node, degree in self.degree if "S" not in node and degree > 2))
+        except:
+            return 1
 
     @property
-    def std_degree_except_leaves(self) -> int:
-        return np_.std(list(degree for node, degree in self.degree if "S" not in node and degree != 1))
-
-    @property
-    def highest_degree_w_nodes(self) -> Tuple[int, List[str]]:
-        #
-        max_degree = -1
-        at_nodes = None
-        for node, degree in self.degree:
-            if "S" not in node:
-                if degree > max_degree:
-                    max_degree = degree
-                    at_nodes = [node]
-                elif degree == max_degree:
-                    at_nodes.append(node)
-
-        return max_degree, at_nodes
+    def std_degree_except_leaves_and_roots(self) -> int:
+        try:
+            return np_.std(list(degree for node, degree in self.degree if "S" not in node and degree > 2))
+        except:
+            return 0
 
     @property
     def edge_lengths(self) -> Tuple[float, ...]:
@@ -174,3 +174,18 @@ class skl_graph_t(skl_nfgraph_t):
     def secondary_edge_lengths(self, soma) -> Tuple[float, ...]:
         #
         return tuple(edge.lengths.length for node1_id, node2_id, edge in self.edges.data("as_edge_t") if (node1_id not in soma.graph_roots.values()) and (node2_id not in soma.graph_roots.values()))
+
+    def highest_degree_w_nodes(self, soma) -> Tuple[int, List[str]]:
+        #
+        max_degree = 1
+        at_nodes = None
+        for node, degree in self.degree:
+            if "S" not in node:
+                if degree > 2:
+                    if degree > max_degree:
+                        max_degree = degree
+                        at_nodes = [node]
+                    elif degree == max_degree:
+                        at_nodes.append(node)
+
+        return max_degree, at_nodes
\ No newline at end of file