From ff2018db89ff50605d168dd8ae2e7102c68f042f Mon Sep 17 00:00:00 2001
From: DEBREUVE Eric <eric.debreuve@inria.fr>
Date: Thu, 19 Sep 2019 11:46:38 +0200
Subject: [PATCH] refact+reorg+excluding already connected extensions within
 loop

---
 connection.py             |  52 +++++++++++++
 extension.py              |  33 ++++----
 feedback.py               |  56 ++++++++++++++
 frangi_3d.py              |   2 +
 glial_cmp.py              |  30 ++++---
 main.py                   | 159 +++++++++++++++-----------------------
 map_labeling.py           |   2 +-
 ref_result_nutrimorph.txt |  69 +++++++++++++++++
 soma.py                   |  57 ++++----------
 9 files changed, 293 insertions(+), 167 deletions(-)
 create mode 100644 connection.py
 create mode 100644 feedback.py
 create mode 100644 ref_result_nutrimorph.txt

diff --git a/connection.py b/connection.py
new file mode 100644
index 0000000..0f8794c
--- /dev/null
+++ b/connection.py
@@ -0,0 +1,52 @@
+import dijkstra_1_to_n as dk_
+from extension import extension_t
+from soma import soma_t
+from type import array_t, site_h, site_path_h
+
+import itertools as it_
+from typing import Callable, Sequence, Tuple
+
+import numpy as np_
+
+
+def CandidateConnections(
+    somas: Sequence[soma_t],
+    influence_map: array_t,
+    dist_to_closest: array_t,
+    extensions: Sequence[extension_t],
+    max_straight_sq_dist: float = np_.inf,
+) -> list:
+    #
+    candidate_conn_nfo = []  # conn=connection
+    for soma, extension in it_.product(somas, extensions):
+        new_candidates = extension.EndPointsForSoma(soma.uid, influence_map)
+        candidate_conn_nfo.extend(
+            (ep_idx, soma, extension, end_point)
+            for ep_idx, end_point in enumerate(new_candidates)
+            if dist_to_closest[end_point] <= max_straight_sq_dist
+        )
+    candidate_conn_nfo.sort(key=lambda elm: dist_to_closest[elm[3]])
+
+    return candidate_conn_nfo
+
+
+def ShortestPathTo(
+    point: site_h,
+    costs: array_t,
+    candidate_points_fct: Callable,
+    max_straight_sq_dist: float = np_.inf,
+) -> Tuple[site_path_h, float]:
+    #
+    candidate_points, candidate_indexing = candidate_points_fct(
+        point, max_straight_sq_dist
+    )
+    if candidate_points is None:
+        return (), np_.inf
+
+    costs[point] = 0.0
+    costs[candidate_indexing] = 0.0
+    path, length = dk_.DijkstraShortestPath(costs, point, candidate_points)
+    costs[point] = np_.inf
+    costs[candidate_indexing] = np_.inf
+
+    return path, length
diff --git a/extension.py b/extension.py
index c0b1271..8633e4b 100644
--- a/extension.py
+++ b/extension.py
@@ -5,7 +5,7 @@ from glial_cmp import glial_cmp_t
 import map_labeling as ml_
 from type import array_t, site_h
 
-from typing import Sequence, Tuple
+from typing import Optional, Sequence, Tuple
 
 import numpy as np_
 import skimage.filters as fl_
@@ -26,27 +26,27 @@ class extension_t(glial_cmp_t):
     # soma_uid: connected to a soma somewhere upstream (as opposed to downstream extensions)
     # extensions: downstream (as opposed to being upstreamed connected)
     #
-    __slots__ = ("scales", "end_points", "ep_closest_somas", "soma_uid", "__cache__")
+    __slots__ = ("scales", "end_points", "soma_uid", "__cache__")
 
     def __init__(self):
         #
         super().__init__()
         self.scales = None
         self.end_points = None
-        self.ep_closest_somas = None
         self.soma_uid = None
         self.__cache__ = None
 
     @classmethod
-    def FromMaps(
-        cls, lmp: array_t, end_point_lmp: array_t, scales: array_t, uid: int
-    ) -> extension_t:
+    def FromMap(cls, lmp: array_t, scales: array_t, uid: int) -> extension_t:
         #
         instance = cls()
 
-        instance.InitializeFromMaps(lmp, uid)
-        instance.scales = scales[instance.sites]
+        bmp = lmp == uid
+        instance.InitializeFromMap(bmp, uid)
+        end_point_map = extension_t.EndPointMap(bmp)
+        end_point_lmp = end_point_map * lmp
         instance.end_points = (end_point_lmp == uid).nonzero()
+        instance.scales = scales[instance.sites]
         instance.__cache__ = {}
 
         return instance
@@ -60,13 +60,11 @@ class extension_t(glial_cmp_t):
 
         return self.__cache__[pty_name]
 
-    def CaptureClosestSomas(self, soma_influence_map: array_t) -> None:
-        #
-        self.ep_closest_somas = soma_influence_map[self.end_points]
-
-    def EndPointsForSoma(self, soma_uid: int) -> Tuple[site_h, ...]:
+    def EndPointsForSoma(
+        self, soma_uid: int, influence_map: array_t
+    ) -> Tuple[site_h, ...]:
         #
-        ep_bmp = self.ep_closest_somas == soma_uid  # bmp=boolean map
+        ep_bmp = influence_map[self.end_points] == soma_uid  # bmp=boolean map
         if ep_bmp.any():
             end_point_idc = ep_bmp.nonzero()[0]
             end_points = self.end_points_as_array[:, end_point_idc]
@@ -93,13 +91,14 @@ class extension_t(glial_cmp_t):
             f"Ext.{self.uid}, "
             f"sites={self.sites[0].__len__()}, "
             f"endpoints={self.end_points[0].__len__()}, "
-            f"closest soma(s)={np_.unique(self.ep_closest_somas).__str__()[1:-1]}, "
             f"soma={self.soma_uid}, "
             f"extensions={n_extensions}"
         )
 
     @staticmethod
-    def ExtensionWithSite(extensions: Sequence[extension_t], site: site_h):
+    def ExtensionWithSite(
+        extensions: Sequence[extension_t], site: site_h
+    ) -> Optional[extension_t]:
         #
         for extension in extensions:
             if site in tuple(zip(*extension.sites)):
@@ -205,7 +204,7 @@ def __HysterisisImage__(image: array_t, low: float, high: float) -> array_t:
     # hight = (image_f <hight).astype(int)
 
     result = fl_.apply_hysteresis_threshold(image, low, high)
-    result = result.astype(np_.int8)
+    result = result.astype(np_.int8, copy=False)
 
     return result
 
diff --git a/feedback.py b/feedback.py
new file mode 100644
index 0000000..179844f
--- /dev/null
+++ b/feedback.py
@@ -0,0 +1,56 @@
+from extension import extension_t
+import plot as ot_
+from soma import soma_t
+
+from typing import Sequence, Tuple
+
+import matplotlib.pyplot as pl_
+
+
+def PrintConnectedExtensions(extensions: Sequence[extension_t]) -> None:
+    #
+    ext_uids = tuple(
+        extension.uid for extension in extensions if extension.soma_uid is not None
+    )
+    print(
+        f"    Connected Ext = {ext_uids.__len__()}"
+        f"/{extensions.__len__()}\n"
+        f"    {ext_uids}"
+    )
+
+
+def PlotSomas(somas: Sequence[soma_t], som_nfo: dict, axes: dict) -> None:
+    #
+    for soma in somas:
+        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_.matshow(som_nfo["dist_to_closest"].max(axis=0)), pl_.title("Soma Distances")
+
+
+def PlotExtensions(
+    extensions: Sequence[extension_t], ext_nfo: dict, img_shape: Tuple[int, ...]
+) -> None:
+    #
+    for extension in extensions:
+        _ = ot_.PlotExtensions(extension, img_shape)
+        pl_.title(f"Extension.{extension.uid}")
+    pl_.matshow(ext_nfo["map"].max(axis=0))
+    pl_.title("Extensions Extremities")
+
+
+def PlotSomasWithExtensions(somas: Sequence[soma_t], som_nfo: dict, which: str) -> None:
+    #
+    any_plot = False
+    for soma in filter(lambda elm: elm.has_extensions, somas):
+        if which == "with_ext_of_ext":
+            for extension in soma.extensions:
+                if not extension.has_extensions:
+                    continue
+        _ = ot_.PlotSomaWithExtensions(soma, som_nfo["lmp"])
+        pl_.title(f"Soma.{soma.uid} + Ext.{tuple(ext.uid for ext in soma.extensions)}")
+        any_plot = True
+
+    if any_plot:
+        pl_.matshow(som_nfo["soma_w_ext_lmp"].max(axis=0))
+        pl_.title("Somas + Extensions")
diff --git a/frangi_3d.py b/frangi_3d.py
index d416e32..31e160f 100644
--- a/frangi_3d.py
+++ b/frangi_3d.py
@@ -1,9 +1,11 @@
 from type import array_t
 
 import multiprocessing as pl_
+
 # from ctypes import c_bool as c_bool_t
 from ctypes import c_float as c_float_t
 from multiprocessing.sharedctypes import Array as shared_array_t
+
 # from multiprocessing.sharedctypes import Value as shared_value_t
 from typing import List, Tuple
 
diff --git a/glial_cmp.py b/glial_cmp.py
index bf512b4..740f5d0 100644
--- a/glial_cmp.py
+++ b/glial_cmp.py
@@ -1,6 +1,8 @@
 from __future__ import annotations
 
-from type import array_t, site_path_h
+from type import array_t, np_array_picker_h, py_array_picker_h, site_path_h
+
+from typing import Dict, List, Optional, Tuple
 
 import numpy as np_
 
@@ -11,20 +13,25 @@ class glial_cmp_t:
 
     def __init__(self):
         #
-        self.uid = None
-        self.sites = None
-        self.connection_path = None
-        self.extensions = None
-        self.img_shape = None
+        self.uid = None  # type: Optional[int]
+        self.sites = None  # type: Optional[np_array_picker_h]
+        self.connection_path = None  # type: Optional[Dict[int, py_array_picker_h]]
+        self.extensions = None  # type: Optional[List[glial_cmp_t]]
+        self.img_shape = None  # type: Optional[Tuple[int, ...]]
 
-    def InitializeFromMaps(self, lmp: array_t, uid: int) -> None:
+    def InitializeFromMap(self, bmp: array_t, uid: int) -> None:
         #
         self.uid = uid
         # sites: might contain voxels that could be removed w/o breaking connectivity
-        self.sites = (lmp == uid).nonzero()
+        self.sites = bmp.nonzero()
         self.connection_path = {}
         self.extensions = []
-        self.img_shape = lmp.shape
+        self.img_shape = bmp.shape
+
+    @property
+    def has_extensions(self) -> bool:
+        #
+        return self.extensions.__len__() > 0
 
     def ExtendWith(
         self, extension: glial_cmp_t, through: site_path_h, costs: array_t
@@ -38,15 +45,16 @@ class glial_cmp_t:
         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
+        # 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
 
-    def BackReferenceSoma(self, glial_cmp: glial_cmp_t):
+    def BackReferenceSoma(self, glial_cmp: glial_cmp_t) -> None:
         #
         raise NotImplementedError
diff --git a/main.py b/main.py
index 0742fae..b289d5d 100644
--- a/main.py
+++ b/main.py
@@ -1,8 +1,8 @@
-import extension as ext_
-import plot as ot_
-import soma as soma_
+import connection as cn_
+import extension as xt_
+import feedback as fb_
+import soma as sm_
 
-import itertools as it_
 import time as tm_
 
 import matplotlib.pyplot as pl_
@@ -17,24 +17,24 @@ print(f"STARTED: {tm_.strftime('%a, %b %d %Y @ %H:%M:%S')}")
 start_time = tm_.time()
 
 
-soma_t = soma_.soma_t
-extension_t = ext_.extension_t
+soma_t = sm_.soma_t
+extension_t = xt_.extension_t
 
 
 # --- Parameters
+data_path = "./DIO_6H_6_1.70bis_2.2_3.tif"
 run = ("soma", "extension", "som-ext", "ext-ext")
 with_plot = False
-data_path = "./DIO_6H_6_1.70bis_2.2_3.tif"
 
-soma_low = 0.15
-soma_high = 0.7126
-ext_low = 0.2  # ext_low = 9.0e-4
-ext_high = 0.6  # high_ext = 8.0e-3
+soma_low_c = 0.15
+soma_high_c = 0.7126
+ext_low_c = 0.2  # 0.02  # 0.2  # ext_low_c = 9.0e-4
+ext_high_c = 0.6  # 0.04  # 0.6  # high_ext = 8.0e-3
 
-soma_selem = mp_.disk(2)
+soma_selem_c = mp_.disk(2)
 
-max_straight_sq_dist = 30 ** 2
-max_weighted_length = 20.0
+max_straight_sq_dist_c = 30 ** 2
+max_weighted_length_c = 20.0
 
 min_area_c = 1000
 
@@ -44,8 +44,8 @@ image = io_.imread(data_path)
 image = cl_.rgb2gray(image)[:, 512:, 512:]
 img_shape = image.shape
 
-image_for_soma = soma_.NormalizedImage(image)
-image_for_ext = ext_.NormalizedImage(image)
+image_for_soma = sm_.NormalizedImage(image)
+image_for_ext = xt_.NormalizedImage(image)
 costs = 1.0 / (image + 1.0)
 
 
@@ -58,40 +58,27 @@ n_extensions = 0
 ext_nfo = {}  # ext=extension, nfo=info
 extensions = None  # Tuple of extension objects
 
-axes = None
+axes = {}
 
 
 # --- Somas
 if "soma" in run:
     print("\n--- Soma Detection")
 
-    som_nfo["map"] = soma_t.Map(image_for_soma, soma_low, soma_high, soma_selem)
+    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)
-    som_nfo["contour_map"] = soma_t.ContourMap(som_nfo["map"])
-
     som_nfo["lmp"], n_somas = ms_.label(som_nfo["map"], return_num=True)
-    som_nfo["contour_lmp"] = som_nfo["contour_map"] * som_nfo["lmp"]
-
     som_nfo["dist_to_closest"], som_nfo["influence_map"] = soma_t.InfluenceMaps(
         som_nfo["lmp"]
     )
 
     somas = tuple(
-        soma_t().FromMaps(som_nfo["lmp"], som_nfo["contour_lmp"], uid)
-        for uid in range(1, n_somas + 1)
+        soma_t().FromMap(som_nfo["lmp"], uid) for uid in range(1, n_somas + 1)
     )
 
     print(f"    n = {n_somas}")
-
     if with_plot:
-        axes = {}
-        for soma in somas:
-            axes[soma.uid] = ot_.PlotLMap(som_nfo["lmp"], labels=soma.uid)
-            pl_.title(f"Soma.{soma.uid}")
-        # pl_.matshow(som_nfo["map"].max(axis=0)), pl_.title("Somas")
-        # pl_.matshow(som_nfo["contour_map"].max(axis=0)), pl_.title("Soma Contours")
-        pl_.matshow(som_nfo["influence_map"].max(axis=0)), pl_.title("Soma Influencess")
-        pl_.matshow(som_nfo["dist_to_closest"].max(axis=0)), pl_.title("Soma Distances")
+        fb_.PlotSomas(somas, som_nfo, axes)
 
 
 # -- Extentions
@@ -99,33 +86,20 @@ if "extension" in run:
     print("\n--- Extension Detection")
 
     enhanced_ext, ext_scales = extension_t.EnhancedForDetection(image_for_ext)
-    ext_nfo["coarse_map"] = extension_t.CoarseMap(enhanced_ext, ext_low, ext_high)
+    ext_nfo["coarse_map"] = extension_t.CoarseMap(enhanced_ext, ext_low_c, ext_high_c)
     ext_nfo["coarse_map"] = extension_t.FilteredCoarseMap(ext_nfo["coarse_map"])
     ext_nfo["map"] = extension_t.FineMapFromCoarseMap(ext_nfo["coarse_map"])
     ext_nfo["map"][som_nfo["map"] > 0] = 0
-    ext_nfo["end_point_map"] = extension_t.EndPointMap(ext_nfo["map"])
-
     ext_nfo["lmp"], n_extensions = ms_.label(ext_nfo["map"], return_num=True)
-    ext_nfo["end_point_lmp"] = ext_nfo["end_point_map"] * ext_nfo["lmp"]
 
     extensions = tuple(
-        extension_t().FromMaps(
-            ext_nfo["lmp"], ext_nfo["end_point_lmp"], ext_scales, uid
-        )
+        extension_t().FromMap(ext_nfo["lmp"], ext_scales, uid)
         for uid in range(1, n_extensions + 1)
     )
-    for extension in extensions:
-        extension.CaptureClosestSomas(som_nfo["influence_map"])
 
     print(f"    n = {n_extensions}")
-
     if with_plot:
-        for extension in extensions:
-            _ = ot_.PlotExtensions(extension, img_shape)
-            pl_.title(f"Extension.{extension.uid}")
-        # pl_.matshow(ext_nfo["map"].max(axis=0)), pl_.title("Extensions")
-        pl_.matshow((10 * ext_nfo["end_point_map"] + ext_nfo["map"]).max(axis=0))
-        pl_.title("Extensions Extremities")
+        fb_.PlotExtensions(extensions, ext_nfo, img_shape)
 
 
 # -- Soma-Extention
@@ -137,51 +111,37 @@ if "som-ext" in run:
     # (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 = []  # conn=connection
-    for soma, extension in it_.product(somas, extensions):
-        new_candidates = extension.EndPointsForSoma(soma.uid)
-        candidate_conn_nfo.extend(
-            (soma, extension, end_point) for end_point in new_candidates
-        )
-    candidate_conn_nfo.sort(key=lambda elm: som_nfo["dist_to_closest"][elm[2]])
+    candidate_conn_nfo = cn_.CandidateConnections(
+        somas,
+        som_nfo["influence_map"],
+        som_nfo["dist_to_closest"],
+        extensions,
+        max_straight_sq_dist_c,
+    )
+
+    for ep_idx, soma, extension, end_point in candidate_conn_nfo:
+        if extension.soma_uid is not None:
+            continue
 
-    for c_idx, (soma, extension, end_point) in enumerate(candidate_conn_nfo):
-        print(f"    Soma.{soma.uid} <-?-> Ext.{extension.uid}({c_idx})", end="")
-        path, length = soma_t.ShortestPathFrom(
+        print(f"    Soma.{soma.uid} <-?-> Ext.{extension.uid}({ep_idx})", end="")
+        path, length = cn_.ShortestPathTo(
             end_point,
             costs,
             soma.ContourPointsCloseTo,
-            max_straight_sq_dist=max_straight_sq_dist,
+            max_straight_sq_dist=max_straight_sq_dist_c,
         )
-        if length <= max_weighted_length:
+        if length <= max_weighted_length_c:
             soma.ExtendWith(extension, path, costs)
             print(": Made")
         else:
             print("")
-    som_nfo["soma_w_ext_lmp"] = soma_t.SomasWithExtensionsLMap(somas)
 
     # for soma in somas:
     #     soma.Extend(extensions, som_nfo["dist_to_closest"], costs)
 
-    connected_ext_uids = tuple(
-        extension.uid for extension in extensions if extension.soma_uid is not None
-    )
-    print(
-        f"    Connected Ext = {connected_ext_uids.__len__()}"
-        f"/{extensions.__len__()}\n"
-        f"    {connected_ext_uids}"
-    )
-
+    fb_.PrintConnectedExtensions(extensions)
     if with_plot:
-        # Somas without extensions have already been plotted
-        for soma in filter(lambda elm: elm.has_extensions, somas):
-            _ = ot_.PlotSomaWithExtensions(soma, som_nfo["lmp"])
-            pl_.title(
-                f"Soma.{soma.uid} + Ext.{tuple(ext.uid for ext in soma.extensions)}"
-            )
-        pl_.matshow(som_nfo["soma_w_ext_lmp"].max(axis=0)), pl_.title(
-            "Somas + Extensions"
-        )
+        fb_.PlotSomasWithExtensions(somas, som_nfo, "all")
 
 
 # -- Extention-Extention
@@ -197,27 +157,28 @@ if "ext-ext" in run:
         unconnected_extensions = list(
             filter(lambda ext: ext.soma_uid is None, extensions)
         )
-        for extension in unconnected_extensions:
-            extension.CaptureClosestSomas(som_nfo["influence_map"])
-
-        candidate_conn_nfo = []  # conn=connection
-        for soma, extension in it_.product(somas, unconnected_extensions):
-            new_candidates = extension.EndPointsForSoma(soma.uid)
-            candidate_conn_nfo.extend(
-                (soma, extension, end_point) for end_point in new_candidates
-            )
-        candidate_conn_nfo.sort(key=lambda elm: som_nfo["dist_to_closest"][elm[2]])
+
+        candidate_conn_nfo = cn_.CandidateConnections(
+            somas,
+            som_nfo["influence_map"],
+            som_nfo["dist_to_closest"],
+            unconnected_extensions,
+            max_straight_sq_dist_c,
+        )
 
         should_look_for_connections = False
-        for c_idx, (soma, extension, end_point) in enumerate(candidate_conn_nfo):
-            print(f"    Soma.{soma.uid} <-?-> Ext.{extension.uid}({c_idx})", end="")
-            path, length = soma_t.ShortestPathFrom(
+        for ep_idx, soma, extension, end_point in candidate_conn_nfo:
+            if extension.soma_uid is not None:
+                continue
+
+            print(f"    Soma.{soma.uid} <-?-> Ext.{extension.uid}({ep_idx})", end="")
+            path, length = cn_.ShortestPathTo(
                 end_point,
                 costs,
                 soma.ExtensionPointsCloseTo,
-                max_straight_sq_dist=max_straight_sq_dist,
+                max_straight_sq_dist=max_straight_sq_dist_c,
             )
-            if length <= max_weighted_length:
+            if length <= max_weighted_length_c:
                 tgt_extenstion = extension_t.ExtensionWithSite(extensions, path[-1])
                 tgt_extenstion.ExtendWith(extension, path, costs)
                 should_look_for_connections = True
@@ -225,13 +186,17 @@ if "ext-ext" in run:
             else:
                 print("")
 
+    fb_.PrintConnectedExtensions(extensions)
+    if with_plot:
+        fb_.PlotSomasWithExtensions(somas, som_nfo, "with_ext_of_ext")
+
 
 # -- Summary
 print("\n")
 for soma in somas:
     print(soma)
-for extension in extensions:
-    print(extension)
+# for extension in extensions:
+#     print(extension)
 
 elapsed_time = tm_.gmtime(tm_.time() - start_time)
 print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")
diff --git a/map_labeling.py b/map_labeling.py
index 9e7cc5b..c75e730 100644
--- a/map_labeling.py
+++ b/map_labeling.py
@@ -22,7 +22,7 @@ def PartLMap(map_: array_t) -> array_t:
     # neighboring pixels that belong to the skeleton (as expected,
     # isolated pixels receive 0).
     #
-    result = map_.copy()
+    result = np_.array(map_, dtype=np_.int8)
     result[result > 0] = 1
     padded_sm = np_.pad(map_ > 0, 1, "constant")
 
diff --git a/ref_result_nutrimorph.txt b/ref_result_nutrimorph.txt
new file mode 100644
index 0000000..497a805
--- /dev/null
+++ b/ref_result_nutrimorph.txt
@@ -0,0 +1,69 @@
+/usr/bin/python3.7 /home/eric/Code/project/mno/nutrimorph/amina/main.py
+STARTED: Wed, Sep 18 2019 @ 17:28:28
+
+--- Soma Detection
+    n = 14
+
+--- Extension Detection
+/!\ Reading from precomputed data file
+    n = 64
+
+--- Soma <-> Extension
+    Soma.1 <-?-> Ext.9(0): Made
+    Soma.1 <-?-> Ext.21(0): Made
+    Soma.1 <-?-> Ext.31(0): Made
+    Soma.1 <-?-> Ext.32(0): Made
+    Soma.1 <-?-> Ext.35(1): Made
+    Soma.1 <-?-> Ext.38(0): Made
+    Soma.2 <-?-> Ext.1(1): Made
+    Soma.2 <-?-> Ext.10(1): Made
+    Soma.2 <-?-> Ext.18(1): Made
+    Soma.2 <-?-> Ext.19(1): Made
+    Soma.2 <-?-> Ext.23(0): Made
+    Soma.2 <-?-> Ext.25(1): Made
+    Soma.2 <-?-> Ext.27(1): Made
+    Soma.2 <-?-> Ext.30(0): Made
+    Soma.2 <-?-> Ext.34(0): Made
+    Soma.2 <-?-> Ext.36(1): Made
+    Soma.2 <-?-> Ext.37(0): Made
+    Soma.2 <-?-> Ext.41(1): Made
+    Soma.2 <-?-> Ext.43(1): Made
+    Soma.2 <-?-> Ext.44(0): Made
+    Soma.2 <-?-> Ext.46(0): Made
+    Soma.2 <-?-> Ext.47(0): Made
+    Soma.2 <-?-> Ext.50(0): Made
+    Soma.4 <-?-> Ext.2(0): Made
+    Soma.4 <-?-> Ext.4(0): Made
+    Soma.4 <-?-> Ext.5(0): Made
+    Soma.4 <-?-> Ext.6(0): Made
+    Soma.4 <-?-> Ext.7(0): Made
+    Soma.4 <-?-> Ext.11(0): Made
+    Soma.4 <-?-> Ext.12(1): Made
+    Soma.4 <-?-> Ext.15(0): Made
+    Soma.4 <-?-> Ext.17(0): Made
+    Soma.4 <-?-> Ext.20(1): Made
+    Soma.13 <-?-> Ext.64(0): Made
+    Soma.14 <-?-> Ext.33(1): Made
+    Soma.14 <-?-> Ext.48(1): Made
+    Soma.14 <-?-> Ext.58(0): Made
+    Soma.14 <-?-> Ext.59(0): Made
+    Soma.1 <-?-> Ext.26(0): Made
+    Soma.4 <-?-> Ext.14(1): Made
+    Soma.14 <-?-> Ext.39(1): Made
+    Soma.4 <-?-> Ext.13(0): Made
+    Soma.9 <-?-> Ext.57(0): Made
+    Soma.9 <-?-> Ext.55(1)
+    Soma.9 <-?-> Ext.55(0)
+    Connected Ext = 43/64
+    (1, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 15, 17, 18, 19, 20, 21, 23, 25, 26, 27, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 41, 43, 44, 46, 47, 48, 50, 57, 58, 59, 64)
+
+--- Extension <-> Extension
+    Soma.9 <-?-> Ext.55(1)
+    Soma.9 <-?-> Ext.55(0)
+
+
+
+Elapsed Time=00h 03m 53s
+DONE: Wed, Sep 18 2019 @ 17:32:22
+
+Process finished with exit code 0
diff --git a/soma.py b/soma.py
index 37c7b30..af24ece 100644
--- a/soma.py
+++ b/soma.py
@@ -1,14 +1,12 @@
 from __future__ import annotations
 
-import dijkstra_1_to_n as dk_
 from glial_cmp import glial_cmp_t
 import map_labeling as ml_
-from extension import extension_t
-from type import array_t, py_array_picker_h, site_h, site_path_h
+from type import array_t, py_array_picker_h, site_h
 
 from collections import namedtuple as namedtuple_t
 import sys as sy_
-from typing import Callable, Optional, Sequence, Tuple
+from typing import Optional, Sequence, Tuple
 
 import numpy as np_
 import scipy.ndimage as im_
@@ -34,21 +32,19 @@ class soma_t(glial_cmp_t):
         self.contour_points = None
 
     @classmethod
-    def FromMaps(cls, lmp: array_t, contour_lmp: array_t, uid: int) -> soma_t:
+    def FromMap(cls, lmp: array_t, uid: int) -> soma_t:
         #
         instance = cls()
 
-        instance.InitializeFromMaps(lmp, uid)
+        bmp = lmp == uid
+        instance.InitializeFromMap(bmp, uid)
+        contour_map = soma_t.ContourMap(bmp)
+        contour_lmp = contour_map * lmp
         instance.contour_points = tuple(zip(*(contour_lmp == uid).nonzero()))
 
         return instance
 
-    @property
-    def has_extensions(self) -> bool:
-        #
-        return self.extensions.__len__() > 0
-
-    def Extensions(self, max_level: int = sy_.maxsize) -> Tuple[extension_t]:
+    def Extensions(self, max_level: int = sy_.maxsize) -> Tuple[glial_cmp_t]:
         #
         # max_level=1: directly connected extensions
         # ...
@@ -88,6 +84,9 @@ class soma_t(glial_cmp_t):
     def ExtensionPointsCloseTo(
         self, point: site_h, max_distance: float
     ) -> Tuple[Optional[Tuple[site_h, ...]], Optional[py_array_picker_h]]:
+        #
+        # Leave connection paths apart because only detected extension pieces
+        # (as opposed to invented=connection paths) are considered valid connection targets
         #
         points = []
         for extension in self.Extensions():
@@ -126,7 +125,7 @@ class soma_t(glial_cmp_t):
     #     while should_look_for_connections:
     #         som_ext_path_nfos = []
     #         for ep_idx, (ext_end_point, extension) in enumerate(unconnected_candidates):
-    #             path, length = self.ShortestPathFrom(
+    #             path, length = self.ShortestPathTo(
     #                 extension.uid, ext_end_point, ep_idx, costs
     #             )
     #             if length <= max_weighted_dist:
@@ -146,31 +145,9 @@ class soma_t(glial_cmp_t):
     #         else:
     #             should_look_for_connections = False
 
-    @staticmethod
-    def ShortestPathFrom(
-        point: site_h,
-        costs: array_t,
-        soma_candidate_points_fct: Callable,
-        max_straight_sq_dist: float = np_.inf,
-    ) -> Tuple[site_path_h, float]:
-        #
-        soma_candidate_points, candidate_indexing = soma_candidate_points_fct(
-            point, max_straight_sq_dist
-        )
-        if soma_candidate_points is None:
-            return (), np_.inf
-
-        costs[point] = 0.0
-        costs[candidate_indexing] = 0.0
-        path, length = dk_.DijkstraShortestPath(costs, point, soma_candidate_points)
-        costs[point] = np_.inf
-        costs[candidate_indexing] = np_.inf
-
-        return path, length
-
-    def BackReferenceSoma(self, glial_cmp: glial_cmp_t):
+    def BackReferenceSoma(self, glial_cmp: glial_cmp_t) -> None:
         #
-        raise NotImplementedError
+        raise NotImplementedError("Soma do not need to back-reference a soma")
 
     def __str__(self) -> str:
         #
@@ -199,7 +176,7 @@ class soma_t(glial_cmp_t):
         low = low * (max_image - min_image) + min_image
         high = high * (max_image - min_image) + min_image
         result = fl_.apply_hysteresis_threshold(image, low, high)
-        result = result.astype(np_.int8)
+        result = result.astype(np_.int8, copy=False)
 
         for dep in range(image.shape[0]):
             result[dep, :, :] = mp_.closing(result[dep, :, :], selem)
@@ -250,9 +227,7 @@ class soma_t(glial_cmp_t):
         return dist_map, np_.array(map_[tuple(idx_map)])
 
     @staticmethod
-    def SomasWithExtensionsLMap(
-        somas: Sequence[soma_t]
-    ) -> array_t:
+    def SomasWithExtensionsLMap(somas: Sequence[soma_t]) -> array_t:
         #
         shape = somas[0].img_shape
         result = np_.zeros(shape, dtype=np_.int64)
-- 
GitLab