recfactoring

brick/component/connection.py

@@ -31,13 +31,13 @@
 
 import brick.processing.dijkstra_1_to_n as dk_
 from brick.component.extension import extension_t
+from brick.component.glial_cmp import glial_cmp_t
 from brick.component.soma import soma_t
 from brick.general.type import array_t, site_h, site_path_h
 
 import itertools as it_
 from typing import Callable, Sequence, Tuple
 
-import cython as cy_
 import numpy as np_
 
 
@@ -46,7 +46,7 @@ def CandidateConnections(
     influence_map: array_t,
     dist_to_closest: array_t,
     extensions: Sequence[extension_t],
-    max_straight_sq_dist: cy_.double = np_.inf,
+    max_straight_sq_dist: float = np_.inf,
 ) -> list:
     #
     candidate_conn_nfo = []  # conn=connection
@@ -84,3 +84,22 @@ def ShortestPathFromToN(
     costs[candidate_indexing] = np_.inf
 
     return path, length
+
+
+def ValidateConnection(
+    glial_cmp: glial_cmp_t, extension: glial_cmp_t, dijkstra_path: site_path_h, costs: array_t
+) -> None:
+    #
+    connection_path = tuple(zip(*dijkstra_path[1:-1]))
+    if connection_path.__len__() == 0:
+        connection_path = None
+
+    glial_cmp.connection_path[extension.uid] = connection_path
+    glial_cmp.extensions.append(extension)
+    extension.BackReferenceSoma(glial_cmp)
+
+    # TODO: Ideally, these paths should be dilated + put this outside
+    # but in ext-ext connections, there must not be dilation around the current ext
+    # (current ext plays the role of a soma in soma-ext step)
+    if connection_path is not None:
+        costs[connection_path] = np_.inf

brick/component/extension.py

@@ -70,8 +70,7 @@ class extension_t(glial_cmp_t):
         bmp = lmp == uid
         instance.InitializeFromMap(bmp, uid)
         end_point_map = cls.EndPointMap(bmp)
-        end_point_lmp = end_point_map * lmp
-        instance.end_points = (end_point_lmp == uid).nonzero()
+        instance.end_points = end_point_map.nonzero()
         instance.scales = scales[instance.sites]
         instance.__cache__ = {}
 
@@ -85,7 +84,7 @@ class extension_t(glial_cmp_t):
     @property
     def end_points_as_array(self) -> array_t:
         #
-        pty_name = "end_points_as_array"
+        pty_name = extension_t.end_points_as_array.__name__
         if pty_name not in self.__cache__:
             self.__cache__[pty_name] = np_.array(self.end_points)
 

brick/component/glial_cmp.py

@@ -33,6 +33,7 @@ from __future__ import annotations
 
 from brick.general.type import array_t, np_array_picker_h, py_array_picker_h, site_path_h
 
+from abc import abstractmethod
 from typing import Dict, List, Optional, Tuple
 
 import numpy as np_
@@ -66,28 +67,7 @@ class glial_cmp_t:
         #
         return self.extensions.__len__() > 0
 
-    def ExtendWith(
-        self, extension: glial_cmp_t, through: site_path_h, costs: array_t
-    ) -> None:
-        #
-        connection_path = tuple(zip(*through[1:-1]))
-        if connection_path.__len__() == 0:
-            connection_path = None
-        extension_path = extension.sites
-
-        self.connection_path[extension.uid] = connection_path
-        self.extensions.append(extension)
-
-        # TODO: put this outside
-        extension.BackReferenceSoma(self)
-
-        # TODO: Ideally, these paths should be dilated + put this outside
-        # but in ext-ext connections, there must not be dilation around the current ext
-        # (current ext plays the role of a soma in soma-ext step)
-        if connection_path is not None:
-            costs[connection_path] = np_.inf
-        costs[extension_path] = np_.inf
-
+    @abstractmethod
     def BackReferenceSoma(self, glial_cmp: glial_cmp_t) -> None:
         #
-        raise NotImplementedError
+        pass

brick/component/soma.py

@@ -66,8 +66,7 @@ class soma_t(glial_cmp_t):
         bmp = lmp == uid
         instance.InitializeFromMap(bmp, uid)
         contour_map = cls.ContourMap(bmp)
-        contour_lmp = contour_map * lmp
-        instance.contour_points = tuple(zip(*(contour_lmp == uid).nonzero()))
+        instance.contour_points = tuple(zip(*contour_map.nonzero()))
 
         return instance
 
@@ -164,9 +163,9 @@ class soma_t(glial_cmp_t):
     def FilteredMap(map_: array_t, min_area: int) -> array_t:
 
         result = map_.copy()
-        lmp = ms_.label(map_)
+        lmp = ms_.label(map_)  # label the image region by measuring the connectivity
 
-        for region in ms_.regionprops(lmp):
+        for region in ms_.regionprops(lmp):  # Measure properties of labeled image regions
             if region.area <= min_area:
                 region_sites = (lmp == region.label).nonzero()
                 result[region_sites] = 0

brick/general/feedback.py

@@ -35,7 +35,9 @@ from brick.component.soma import soma_t
 
 from typing import Sequence, Tuple
 
+import matplotlib as mpl_
 import matplotlib.pyplot as pl_
+mpl_.use('TkAgg')
 
 
 def PrintConnectedExtensions(extensions: Sequence[extension_t]) -> None:
@@ -56,7 +58,9 @@ def PlotSomas(somas: Sequence[soma_t], som_nfo: dict, axes: dict) -> None:
         axes[soma.uid] = ot_.PlotLMap(som_nfo["lmp"], labels=soma.uid)
         pl_.title(f"Soma.{soma.uid}")
     pl_.matshow(som_nfo["influence_map"].max(axis=0)), pl_.title("Soma Influencess")
+    pl_.show(block=True)
     pl_.matshow(som_nfo["dist_to_closest"].max(axis=0)), pl_.title("Soma Distances")
+    pl_.show(block=True)
 
 
 def PlotExtensions(
@@ -68,6 +72,7 @@ def PlotExtensions(
         pl_.title(f"Extension.{extension.uid}")
     pl_.matshow(ext_nfo["map"].max(axis=0))
     pl_.title("Extensions Extremities")
+    pl_.show(block=True)
 
 
 def PlotSomasWithExtensions(somas: Sequence[soma_t], som_nfo: dict, which: str) -> None:
@@ -80,8 +85,11 @@ def PlotSomasWithExtensions(somas: Sequence[soma_t], som_nfo: dict, which: str)
                     continue
         _ = ot_.PlotSomaWithExtensions(soma, som_nfo["lmp"])
         pl_.title(f"Soma.{soma.uid} + Ext.{tuple(ext.uid for ext in soma.extensions)}")
+        pl_.show(block=True)
         any_plot = True
 
     if any_plot:
         pl_.matshow(som_nfo["soma_w_ext_lmp"].max(axis=0))
         pl_.title("Somas + Extensions")
+        pl_.show(block=True)
+

brick/processing/normalization.py → brick/processing/input.py

@@ -36,7 +36,9 @@ import numpy as np_
 
 def NormalizedImage(image: array_t) -> array_t:
     #
-    print("This normalization does not bring anything; left as is for now to avoid the need for changing prms")
+    print(
+        "This normalization does not bring anything; left as is for now to avoid the need for changing prms"
+    )
     nonextreme_values = image[np_.logical_and(image > 0.0, image < image.max())]
 
     if nonextreme_values.size > 0:
@@ -46,3 +48,15 @@ def NormalizedImage(image: array_t) -> array_t:
         result = image.astype(np_.float32)
 
     return result
+
+
+def DijkstraCosts(image: array_t, som_map: array_t, ext_map: array_t) -> array_t:
+    #
+    # TODO: Ideally, the extension part should be dilated
+    # but in ext-ext connections, there must not be dilation around the current or the other exts
+    # (current ext plays the role of a soma in soma-ext step)
+    #
+    dijkstra_costs = 1.0 / (image + 1.0)
+    dijkstra_costs[np_.logical_or(som_map > 0, ext_map > 0)] = np_.inf
+
+    return dijkstra_costs

brick/processing/map_labeling.py

@@ -42,7 +42,7 @@ shifts_of_26_neighbors_c = tuple(
     for i in (-1, 0, 1)
     for j in (-1, 0, 1)
     for k in (-1, 0, 1)
-    if i != 0 or j != 0 or k != 0
+    if i != 0 or j != 0 or k != 0   # take every direct neighbors except itself
 )
 
 
@@ -54,8 +54,8 @@ def PartLMap(map_: array_t) -> array_t:
     # isolated pixels receive 0).
     #
     result = np_.array(map_, dtype=np_.int8)
-    result[result > 0] = 1
-    padded_sm = np_.pad(map_ > 0, 1, "constant")
+    result[result > 0] = 1  # binary map
+    padded_sm = np_.pad(map_ > 0, 1, "constant")  # pad for voxels at the border --> map of booleans / binary mask
 
     for shifts in shifts_of_26_neighbors_c:
         result += np_.roll(padded_sm, shifts, axis=(0, 1, 2))[1:-1, 1:-1, 1:-1]

brick/processing/plot.py

@@ -88,7 +88,7 @@ def PlotConnection(
     if new_axes:
         _, axes = __FigAndAxes__(shape, depth_limit)
 
-    # This test is put here but could be move outside this function
+    # TODO This test is put here but could be move outside this function
     if connection is not None:
         axes.plot(
             depth_factor * np_.array(connection[0]),

documentation/notations.txt

+_t=type
+_h=
+_c=constant
+
 som=soma
 ext=extension
 

nutrimorph.py

@@ -41,7 +41,7 @@ import brick.component.connection as cn_
 import brick.component.extension as xt_
 import brick.component.soma as sm_
 import brick.general.feedback as fb_
-import brick.processing.normalization as nm_
+import brick.processing.input as in_
 # from main_prm import *
 
 import os as os_
@@ -54,6 +54,7 @@ import numpy as np_
 import skimage.io as io_
 import skimage.measure as ms_
 
+print(sy_.argv, sy_.argv.__len__())
 
 if sy_.argv.__len__() < 2:
     print("Missing parameter file argument")
@@ -63,7 +64,7 @@ if not (os_.path.isfile(sy_.argv[1]) and os_.access(sy_.argv[1], os_.R_OK)):
     sy_.exit(0)
 
 data_path = None
-with_plot = None
+with_plot = True
 soma_low_c = None
 soma_high_c = None
 ext_low_c = None
@@ -88,9 +89,16 @@ start_time = tm_.time()
 
 # --- Images
 image = io_.imread(data_path)
-assert image.ndim == 3 # TODO: replace assertion with proper exception handling
-image = image[:, 512:, 512:]#562
+
+image = image[:,:,:,1]  # TODO For the green channel - change in case of other channel or greyscale
+# pl_.imshow(image[10,:,:])
+# pl_.show(block=True)
+
+assert image.ndim == 3  # TODO: replace assertion with proper exception handling
+image = image[:, 512:, 512:]  # 562
 img_shape = image.shape
+pl_.matshow(image[image.shape[0] // 2,:,:])
+pl_.show()
 #
 # max_img = image.max()
 # image.fill(image.min())
@@ -98,9 +106,8 @@ img_shape = image.shape
 #
 print(f"IMAGE: s.{img_shape} t.{image.dtype} m.{image.min()} M.{image.max()}")
 
-image_for_soma = nm_.NormalizedImage(image)
-image_for_ext = nm_.NormalizedImage(image)
-costs = 1.0 / (image + 1.0)
+image_for_soma = in_.NormalizedImage(image)
+image_for_ext = in_.NormalizedImage(image)
 print(f"NRM-IMG: t.{image_for_soma.dtype} m.{image_for_soma.min():.2f} M.{image_for_soma.max():.2f}")
 
 
@@ -120,6 +127,7 @@ print("\n--- Soma Detection")
 
 som_nfo["map"] = soma_t.Map(image_for_soma, soma_low_c, soma_high_c, soma_selem_c)
 som_nfo["map"] = soma_t.FilteredMap(som_nfo["map"], min_area_c)
+#TODO: change label into relabel + max
 som_nfo["lmp"], n_somas = ms_.label(som_nfo["map"], return_num=True)
 som_nfo["dist_to_closest"], som_nfo["influence_map"] = soma_t.InfluenceMaps(
     som_nfo["lmp"]
@@ -158,14 +166,13 @@ if with_plot:
     fb_.PlotExtensions(extensions, ext_nfo, img_shape)
 
 
+# -- Preparation for Connections
+dijkstra_costs = in_.DijkstraCosts(image, som_nfo["map"], ext_nfo["map"])
+
+
 # -- Soma-Extention
 print("\n--- Soma <-> Extension")
 
-# TODO: Ideally, the extension part should be dilated
-# but in ext-ext connections, there must not be dilation around the current ext
-# (current ext plays the role of a soma in soma-ext step)
-costs[np_.logical_or(som_nfo["map"] > 0, ext_nfo["map"] > 0)] = np_.inf
-
 candidate_conn_nfo = cn_.CandidateConnections(
     somas,
     som_nfo["influence_map"],
@@ -179,18 +186,18 @@ for ep_idx, soma, extension, end_point in candidate_conn_nfo:
         print(f"    Soma.{soma.uid} <-?-> Ext.{extension.uid}({ep_idx})", end="")
         path, length = cn_.ShortestPathFromToN(
             end_point,
-            costs,
+            dijkstra_costs,
             soma.ContourPointsCloseTo,
             max_straight_sq_dist=max_straight_sq_dist_c,
         )
         if length <= max_weighted_length_c:
-            soma.ExtendWith(extension, path, costs)
+            cn_.ValidateConnection(soma, extension, path, dijkstra_costs)
             print(": Made")
         else:
             print("")
 
 # for soma in somas:
-#     soma.Extend(extensions, som_nfo["dist_to_closest"], costs)
+#     soma.Extend(extensions, som_nfo["dist_to_closest"], dijkstra_costs)
 
 fb_.PrintConnectedExtensions(extensions)
 elapsed_time = tm_.gmtime(tm_.time() - start_time)
@@ -223,13 +230,13 @@ while should_look_for_connections:
             print(f"    Soma.{soma.uid} <-?-> Ext.{extension.uid}({ep_idx})", end="")
             path, length = cn_.ShortestPathFromToN(
                 end_point,
-                costs,
+                dijkstra_costs,
                 soma.ExtensionPointsCloseTo,
                 max_straight_sq_dist=max_straight_sq_dist_c,
             )
             if length <= max_weighted_length_c:
                 tgt_extenstion = extension_t.ExtensionWithSite(extensions, path[-1])
-                tgt_extenstion.ExtendWith(extension, path, costs)
+                cn_.ValidateConnection(tgt_extenstion, extension, path, dijkstra_costs)
                 should_look_for_connections = True
                 print(": Made")
             else:

parameters.py

@@ -33,7 +33,7 @@ import skimage.morphology as mp_
 
 
 data_path = "./data/DIO_6H_6_1.70bis_2.2_3.tif"
-with_plot = False
+with_plot = True
 
 soma_low_c = 0.15
 soma_high_c = 0.7126