From 61eb7e5fd5b48e8aad0eacaafaca5160396f6bcf Mon Sep 17 00:00:00 2001 From: michaeldeistler Date: Wed, 10 Apr 2024 20:36:01 +0200 Subject: [PATCH] Automatically add groups for SWC reader --- jaxley/modules/cell.py | 38 ++++++++++++++++++++++++++---- jaxley/utils/swc.py | 4 ++++ tests/jaxley_identical/test_swc.py | 3 ++- 3 files changed, 40 insertions(+), 5 deletions(-) diff --git a/jaxley/modules/cell.py b/jaxley/modules/cell.py index 7f901bd4..8d94ffa2 100644 --- a/jaxley/modules/cell.py +++ b/jaxley/modules/cell.py @@ -266,9 +266,22 @@ def read_swc( nseg: int, max_branch_len: float = 300.0, min_radius: Optional[float] = None, + assign_groups: bool = False, ): - """Reads SWC file into a `jx.Cell`.""" - parents, pathlengths, radius_fns, _, coords_of_branches = swc_to_jaxley( + """Reads SWC file into a `jx.Cell`. + + Args: + fname: Path to the swc file. + nseg: The number of compartments per branch. + max_branch_len: If a branch is longer than this value it is split into two + branches. + min_radius: If the radius of a reconstruction is below this value it is clipped. + assign_groups: If True, then the identity of reconstructed points in the SWC + file will be used to generate groups `undefined`, `soma`, `axon`, `basal`, + `apical`, `custom`. See here: + http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html + """ + parents, pathlengths, radius_fns, types, coords_of_branches = swc_to_jaxley( fname, max_branch_len=max_branch_len, sort=True, num_lines=None ) nbranches = len(parents) @@ -276,8 +289,8 @@ def read_swc( non_split = 1 / nseg range_ = np.linspace(non_split / 2, 1 - non_split / 2, nseg) - comp = Compartment().initialize() - branch = Branch([comp for _ in range(nseg)]).initialize() + comp = Compartment() + branch = Branch([comp for _ in range(nseg)]) cell = Cell( [branch for _ in range(nbranches)], parents=parents, xyzr=coords_of_branches ) @@ -295,4 +308,21 @@ def read_swc( cell.set("length", lengths_each) cell.set("radius", radiuses_each) + + # Description of SWC file format: + # http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html + ind_name_lookup = { + 0: "undefined", + 1: "soma", + 2: "axon", + 3: "basal", + 4: "apical", + 5: "custom", + } + types = np.asarray(types).astype(int) + if assign_groups: + for type_ind, name in ind_name_lookup.items(): + indices = np.where(types == type_ind)[0].tolist() + if len(indices) > 0: + cell.branch(indices).add_to_group(name) return cell diff --git a/jaxley/utils/swc.py b/jaxley/utils/swc.py index 52504b20..b6fe9d87 100644 --- a/jaxley/utils/swc.py +++ b/jaxley/utils/swc.py @@ -43,6 +43,10 @@ def swc_to_jaxley( radius_fns = [lambda x: content[0, 5] * np.ones_like(x)] + radius_fns sorted_branches = [[0]] + sorted_branches + # Type of padded section is assumed to be of `custom` type: + # http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html + types = [5.0] + types + all_coords_of_branches = [] for i, branch in enumerate(sorted_branches): coords_of_branch = content[np.asarray(branch) - 1, 2:6] diff --git a/tests/jaxley_identical/test_swc.py b/tests/jaxley_identical/test_swc.py index f610fcc2..3e30b478 100644 --- a/tests/jaxley_identical/test_swc.py +++ b/tests/jaxley_identical/test_swc.py @@ -23,7 +23,8 @@ def test_swc_cell(): dirname = os.path.dirname(__file__) fname = os.path.join(dirname, "../morph.swc") - cell = jx.read_swc(fname, nseg=2, max_branch_len=300.0) + cell = jx.read_swc(fname, nseg=2, max_branch_len=300.0, assign_groups=True) + _ = cell.soma # Only to test whether the `soma` group was created. cell.insert(HH()) cell.branch(1).loc(0.0).record() cell.branch(1).loc(0.0).stimulate(current)