Merge branch 'Create_graphs' into 'master'

Create graphs and change the root search

See merge request !2

brick/component/connection.py

@@ -102,6 +102,7 @@ def ValidateConnection(
     if type(glial_cmp) is soma_t:  # restrain the verification to the soma <-> ext step
         extension.ext_root.append([end_point, ep_idx])
         glial_cmp.ext_root.append([end_point, ep_idx])
+        # TODO: delete one of this info (ep_idx or end_point) according to the code following and its needs
 
     # TODO: Ideally, these paths should be dilated + put this outside
     # but in ext-ext connections, there must not be dilation around the current ext

brick/processing/graph_extraction.py

+# Copyright CNRS/Inria/UNS
+# Contributor(s): Eric Debreuve (since 2019), Morgane Nadal (2020)
+#
+# eric.debreuve@cnrs.fr
+#
+# This software is governed by the CeCILL  license under French law and
+# abiding by the rules of distribution of free software.  You can  use,
+# modify and/ or redistribute the software under the terms of the CeCILL
+# license as circulated by CEA, CNRS and INRIA at the following URL
+# "http://www.cecill.info".
+#
+# As a counterpart to the access to the source code and  rights to copy,
+# modify and redistribute granted by the license, users are provided only
+# with a limited warranty  and the software's author,  the holder of the
+# economic rights,  and the successive licensors  have only  limited
+# liability.
+#
+# In this respect, the user's attention is drawn to the risks associated
+# with loading,  using,  modifying and/or developing or reproducing the
+# software by the user in light of its specific status of free software,
+# that may mean  that it is complicated to manipulate,  and  that  also
+# therefore means  that it is reserved for developers  and  experienced
+# professionals having in-depth computer knowledge. Users are therefore
+# encouraged to load and test the software's suitability as regards their
+# requirements in conditions enabling the security of their systems and/or
+# data to be ensured and,  more generally, to use and operate it in the
+# same conditions as regards security.
+#
+# The fact that you are presently reading this means that you have had
+# knowledge of the CeCILL license and that you accept its terms.
+
+import re as re_
+import networkx as nx_
+import numpy as np_
+
+from brick.component.soma import soma_t
+from sklgraph.skl_graph import skl_graph_t
+from typing import Tuple
+
+
+def FindGraphsRoot(somas: Tuple[soma_t, ...], ext_skl_graph: skl_graph_t) -> skl_graph_t:
+    """
+    Find the roots of the {extension+connexion} graphs to be lined to the soma.
+    Add a key "root" (bool) in the dict of nodes attributes.
+    """
+
+    node_degree_bool = tuple(degree == 1 for _, degree in ext_skl_graph.degree)
+    node_coord = tuple(xyz for xyz, _ in ext_skl_graph.degree)
+
+    # get the coordinates of the nodes (x,y,z)
+    coordinates = GetNodesCoordinates(node_coord)
+
+    # get a list with elements = (soma_uid, extension_uid, root coordinates)
+    roots = GetListRoots(somas)
+
+    # for each node in the graph, search among the degree 1 nodes the nodes that are roots (linked to soma)
+    for node in range(len(coordinates)):
+        if node_degree_bool[node]:
+            # compare the coor with end points
+            for ext_root in roots:
+                if ext_root[2] == coordinates[node]:
+                    # add the attribute 'root' = True in the dict of nodes attributes
+                    ext_skl_graph._node[node_coord[node]]['root'] = True
+                    # add the id of soma and of the extension in the nodes attributes dict
+                    ext_skl_graph._node[node_coord[node]]['soma_uid'] = ext_root[0]
+                    ext_skl_graph._node[node_coord[node]]['ext_uid'] = ext_root[1]
+                else:
+                    try:
+                        # verify if the node has already been covered to avoid overwriting on a True value:
+                        if ext_skl_graph._node[node_coord[node]]['root'] is not True:
+                            ext_skl_graph._node[node_coord[node]]['root'] = False
+                    except:
+                        # add the attribute 'root' = False for non-root endpoints in the dict of nodes attributes
+                        ext_skl_graph._node[node_coord[node]]['root'] = False
+        else:
+            # all the nodes with degree > 1 are not roots
+            ext_skl_graph._node[node_coord[node]]['root'] = False
+
+    return ext_skl_graph
+
+    # TODO: add to every edge the id of the extension
+
+
+def GetNodesCoordinates(node_coord: tuple) -> list:
+    """
+    Gets a list of coordinates - tuple (x,y,z) - from the nodes attributes  - str 'x-y-z' - .
+    """
+    coord = []
+    for c in node_coord:
+        coord.append(c)
+
+    for node in range(len(node_coord)):
+        coord_node = coord[node]
+        pattern = '\d+'
+        coord_node = re_.findall(pattern, coord_node)
+        coor = []
+        for i in range(3):
+            coor.append(int(coord_node[i]))
+        coor = tuple(coor)
+        coord[node] = coor
+
+    return coord
+
+
+def GetListRoots(somas: soma_t) -> list:
+    """
+    Gives a list containing the following information for all somas: [soma id: int, extension id: int, root = (x,y,z): tuple]
+    """
+    roots = []
+    for soma in somas:
+        for ext_id, ext_root in enumerate(soma.ext_root):
+            roots.append((soma.uid, soma.extensions[ext_id].uid, ext_root[0]))
+    return roots

brick/processing/input.py

@@ -35,7 +35,7 @@ import math as mt_
 import sys as sy_
 from PIL import Image
 from PIL.ExifTags import TAGS
-import re
+import re as re_
 
 
 def ImageVerification(image: array_t, channel: str) -> array_t:
@@ -46,7 +46,7 @@ def ImageVerification(image: array_t, channel: str) -> array_t:
 
     # Verification of the dimension of the image and its coherence with the parameters (channel)
     elif image.ndim == 3:
-        print('The image has only one color channel.')
+        print('1 CHANNEL')
         if channel is not None:
             raise ValueError('The image has only 3 dimensions. However, a value for the "channel" parameter is '
                              'specified. Give the channel the value None')
@@ -54,8 +54,8 @@ def ImageVerification(image: array_t, channel: str) -> array_t:
     elif image.ndim == 4:
         if channel == 'R' or channel == 'G' or channel == 'B':
             # not changed into if --channel in 'RGB'-- because if channel='RG' => True.
-            print('The image has multiple color channels. The channel', channel,
-                  'is specified in the parameters.')
+            print('MULTIPLE CHANNELS: ', channel,
+                  ' specified in the parameters.')
             image = image[:, :, :, 'RGB'.find(channel)]
 
             # The obtained image must not be constant
@@ -101,6 +101,9 @@ def IntensityNormalizedImage(image: array_t) -> array_t:
 
 
 def FindVoxelDimensionInMicron(data_path: str, size_voxel_in_micron: list = None) -> array_t:
+    '''
+    Find Voxel dimension in micron from the image metadata.
+    '''
     #
     if size_voxel_in_micron is not None:
         return np_.array(size_voxel_in_micron)
@@ -120,7 +123,7 @@ def FindVoxelDimensionInMicron(data_path: str, size_voxel_in_micron: list = None
             voxel_size = []  # Initialize the list of voxel size in str
             for axe in 'XYZ':
                 pattern = 'Voxel' + axe + '.+\= (\d.+E.\d.)'  # Regular expression
-                voxel_size.append(re.findall(pattern, metadata)[0])
+                voxel_size.append(re_.findall(pattern, metadata)[0])
 
             voxel_size = np_.array(list(map(float, voxel_size)))
             voxel_size_micron = 1.0e6 * voxel_size  # Conversion meters in micron

brick/processing/plot.py

@@ -41,18 +41,12 @@ import mpl_toolkits.mplot3d as p3_
 import numpy as np_
 import skimage.measure as ms_
 
-import SimpleITK as sitk
 
-
-def MaximumIntensityProjectionZ(img: array_t, cmap: str ='tab20', output_image_file_name: str = None) -> None:
+def MaximumIntensityProjectionZ(img: array_t, cmap: str ='tab20', axis: int = 0, output_image_file_name: str = None) -> None:
     """ Maximum Image Projection on the Z axis. """
     #
-    image = sitk.GetImageFromArray(img)
-    # image_size = image.GetSize()
-    # print('Image Size = ', image_size)
-    projection = sitk.MaximumProjection(image, 2)
-    np_img = sitk.GetArrayFromImage(projection)
-    pl_.imshow(np_img[0, :, :], cmap=cmap)
+    xy = np_.amax(img, axis=axis)
+    pl_.imshow(xy, cmap=cmap)
     pl_.show(block=True)
     if output_image_file_name is not None:
         pl_.imsave(output_image_file_name, np_img[0, :, :], cmap=cmap)

nutrimorph.py

@@ -49,6 +49,8 @@ from sklgraph.skl_fgraph import skl_graph_t
 from sklgraph.skl_graph import plot_mode_e
 from sklgraph.skl_map import skl_map_t
 
+import brick.processing.graph_extraction as ge_
+
 import os as os_
 import sys as sy_
 import time as tm_
@@ -61,6 +63,7 @@ import skimage.io as io_
 import skimage.morphology as mp_
 import skimage.measure as ms_
 import networkx as nx_
+import re as re_
 
 
 print(sy_.argv, sy_.argv.__len__())
@@ -338,7 +341,8 @@ po_.MaximumIntensityProjectionZ(ext_nfo['lmp_ext'])
 
 
 # Extract the graph of each extension of a soma
-ext_map = skl_map_t.FromShapeMap(ext_nfo['lmp_ext'], store_widths=True, skeletonize=False, do_post_thinning=False)
+print('\n--- Graph extraction\n')
+ext_map = skl_map_t.FromShapeMap(ext_nfo['lmp_ext'], store_widths=True, skeletonize=False)  # , do_post_thinning=True
 ext_skl_graph = skl_graph_t.FromSkeleton(ext_map)
 
 # --- Some info about the skeleton graphs
@@ -351,12 +355,49 @@ print(
     f"Length={ext_skl_graph.length}<-{ext_skl_graph.edge_lengths}\n"
     f"Width=Hom.{ext_skl_graph.reduced_width()}/Het.{ext_skl_graph.heterogeneous_reduced_width()}<-{ext_skl_graph.edge_reduced_widths()}\n"
     f"Area as WxL={ext_skl_graph.reduced_width() * ext_skl_graph.length}\n"
-    f"Area as WW Length={ext_skl_graph.ww_length}<-{ext_skl_graph.edge_ww_lengths}\n")
+    f"Area as WW Length={ext_skl_graph.ww_length}<-{ext_skl_graph.edge_ww_lengths}\n\n")
+
+
+elapsed_time = tm_.gmtime(tm_.time() - start_time)
+print(f"Elapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")
 
 ext_skl_graph.Plot(mode=plot_mode_e.SKL_Curve, w_directions=True, should_block=False)
 ext_skl_graph.Plot(should_block=True)
 pl_.show()
 
+
 if with_plot:
     pl_.show()
 
+# Find the root of the {ext+conn} graphs.
+    # Roots are the nodes of degree 1 that are to be linked to the soma
+print("\n--- Finding {ext+conn} graphs' root")
+
+# nx_.draw_networkx(ext_skl_graph, block=True)
+# pl_.show()
+
+ext_skl_graph = ge_.FindGraphsRoot(somas, ext_skl_graph)
+
+elapsed_time = tm_.gmtime(tm_.time() - start_time)
+print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")
+print(f"DONE: {tm_.strftime('%a, %b %d %Y @ %H:%M:%S')}")
+
+
+# for node_id,node_id2, edge_nfo in ext_skl_graph.edges.data('as_edge_t'):
+# ...     if (ext_skl_graph.degree[node_id] == 1) or (ext_skl_graph.degree[node_id2] == 1):
+# ...         print(edge_nfo.sites)
+
+# # Pour un soma donne
+# # primary_extension_sites = tuple(set(zip(*extension.sites) for extension in soma.extensions)
+# primary_extension_uids = tuple(extension.uid for extension in soma.extensions)
+# for node_id,node_id2, edge_nfo in ext_skl_graph.edges.data('as_edge_t'):
+#     if (ext_skl_graph.degree[node_id] == 1) or (ext_skl_graph.degree[node_id2] == 1):
+#         ext_uid = np_.unique(ext_nfo['lmp'][edge_nfo.sites])[1]
+#         if ext_uid in primary_extension_uids:
+#             if ext_skl_graph.degree[node_id] == 1:
+#                 root_node = node_id
+#             else:
+#                 root_node = node_id2
+#         # print(edge_nfo.sites)
+#         # edge_sites = set(zip(*edge_nfo.sites))
+#         # intersections = map(in)
\ No newline at end of file