diff --git a/neurom/apps/cli.py b/neurom/apps/cli.py index 845d66df0..34d7012d5 100644 --- a/neurom/apps/cli.py +++ b/neurom/apps/cli.py @@ -76,6 +76,8 @@ def view(input_file, is_3d, plane, backend, realistic_diameters): plot(load_morphology(input_file)) if is_matplotlib: + if not is_3d: + plt.axis('equal') plt.show() diff --git a/neurom/check/morphtree.py b/neurom/check/morphtree.py index 531715e7b..0be18c03c 100644 --- a/neurom/check/morphtree.py +++ b/neurom/check/morphtree.py @@ -53,7 +53,7 @@ def is_monotonic(neurite, tol): if sec[point_id + 1][COLS.R] > sec[point_id][COLS.R] + tol: return False # Check that section boundary points satisfy monotonicity - if(node.parent is not None and + if (node.parent is not None and sec[0][COLS.R] > node.parent.points[-1][COLS.R] + tol): return False diff --git a/neurom/features/morphology.py b/neurom/features/morphology.py index ddee3aa19..dd8909dad 100644 --- a/neurom/features/morphology.py +++ b/neurom/features/morphology.py @@ -50,10 +50,11 @@ import numpy as np from neurom import morphmath -from neurom.core.morphology import iter_neurites, iter_segments, Morphology +from neurom.core.morphology import iter_neurites, iter_segments, Morphology, iter_sections from neurom.core.types import tree_type_checker as is_type from neurom.core.dataformat import COLS from neurom.core.types import NeuriteType +from neurom.core.morphology import Section from neurom.exceptions import NeuroMError from neurom.features import feature, NameSpace, neurite as nf, section as sf from neurom.utils import str_to_plane @@ -457,7 +458,9 @@ def neurite_volume_density(morph, neurite_type=NeuriteType.all): @feature(shape=(...,)) -def sholl_crossings(morph, neurite_type=NeuriteType.all, center=None, radii=None): +def sholl_crossings( + morph, neurite_type=NeuriteType.all, center=None, radii=None, distance_type="radial" +): """Calculate crossings of neurites. Args: @@ -466,6 +469,7 @@ def sholl_crossings(morph, neurite_type=NeuriteType.all, center=None, radii=None center(Point): center point, if None then soma center is taken radii(iterable of floats): radii for which crossings will be counted, if None then soma radius is taken + distance_type(str): either `radial` or `path` as distance measure to use Returns: Array of same length as radii, with a count of the number of crossings @@ -480,14 +484,25 @@ def sholl_crossings(morph, neurite_type=NeuriteType.all, center=None, radii=None """ def _count_crossings(neurite, radius): """Used to count_crossings of segments in neurite with radius.""" - r2 = radius ** 2 count = 0 - for start, end in iter_segments(neurite): - start_dist2, end_dist2 = (morphmath.point_dist2(center, start), - morphmath.point_dist2(center, end)) - count += int(start_dist2 <= r2 <= end_dist2 or - end_dist2 <= r2 <= start_dist2) + if distance_type == 'radial': + r2 = radius ** 2 + for start, end in iter_segments(neurite): + start_dist2, end_dist2 = (morphmath.point_dist2(center, start), + morphmath.point_dist2(center, end)) + + count += int(start_dist2 <= r2 <= end_dist2 or + end_dist2 <= r2 <= start_dist2) + + if distance_type == 'path': + for section in iter_sections(neurite): + base_length = sf.section_path_length(section) - section.length + path_dists = base_length + sf.segment_lengths(section, prepend_zero=True) + + forward = np.logical_and(path_dists[:-1] <= radius, path_dists[1:] >= radius) + backward = np.logical_and(path_dists[:-1] >= radius, path_dists[1:] <= radius) + count += sum(np.logical_or(forward, backward)) return count @@ -505,7 +520,9 @@ def _count_crossings(neurite, radius): @feature(shape=(...,)) -def sholl_frequency(morph, neurite_type=NeuriteType.all, step_size=10, bins=None): +def sholl_frequency( + morph, neurite_type=NeuriteType.all, step_size=10, bins=None, distance_type='euclidean' +): """Perform Sholl frequency calculations on a morph. Args: @@ -514,6 +531,7 @@ def sholl_frequency(morph, neurite_type=NeuriteType.all, step_size=10, bins=None step_size(float): step size between Sholl radii bins(iterable of floats): custom binning to use for the Sholl radii. If None, it uses intervals of step_size between min and max radii of ``morphologies``. + distance_type(str): either `euclidean` or `path` as distance measure to use Note: Given a morphology, the soma center is used for the concentric circles, @@ -528,19 +546,34 @@ def sholl_frequency(morph, neurite_type=NeuriteType.all, step_size=10, bins=None neurite_filter = is_type(neurite_type) if bins is None: - min_soma_edge = morph.soma.radius - max_radius_per_neurite = [ - np.max(np.linalg.norm(n.points[:, COLS.XYZ] - morph.soma.center, axis=1)) - for n in morph.neurites if neurite_filter(n) - ] - - if not max_radius_per_neurite: - return [] - - bins = np.arange(min_soma_edge, min_soma_edge + max(max_radius_per_neurite), step_size) - - return sholl_crossings(morph, neurite_type, morph.soma.center, bins) + if distance_type == "euclidean": + max_radius_per_neurite = [ + np.max(np.linalg.norm(n.points[:, COLS.XYZ] - morph.soma.center, axis=1)) + for n in morph.neurites + if neurite_filter(n) + ] + if not max_radius_per_neurite: + return [] + min_distance = morph.soma.radius + max_distance = min_distance + max(max_radius_per_neurite) + + if distance_type == "path": + max_path_distance_per_neurite = [ + max(sf.section_path_length(section) for section in iter_sections(n, Section.ileaf)) + for n in morph.neurites + if neurite_filter(n) + ] + if not max_path_distance_per_neurite: + return [] + min_distance = 0 + max_distance = max(max_path_distance_per_neurite) + + bins = np.arange(min_distance, max_distance, step_size) + + return sholl_crossings( + morph, neurite_type, morph.soma.center, bins, distance_type=distance_type + ) def _extent_along_axis(morph, axis, neurite_type): diff --git a/setup.py b/setup.py index 8adcd2bc7..62c73421a 100644 --- a/setup.py +++ b/setup.py @@ -55,8 +55,6 @@ extras_require={ 'plotly': ['plotly>=3.6.0', 'psutil>=5.5.1'], # for plotly image saving 'docs': [ - 'Jinja2<3.1', # New release 3.1.0 breaks sphinx-bluebrain-theme - 'sphinx', 'sphinx-bluebrain-theme', 'sphinx-autorun', ], diff --git a/tests/check/test_morphology_checks.py b/tests/check/test_morphology_checks.py index 3a3d0186c..0322f0fd2 100644 --- a/tests/check/test_morphology_checks.py +++ b/tests/check/test_morphology_checks.py @@ -387,7 +387,7 @@ def test__bool__(): def test_has_multifurcation(): m = load_morphology(StringIO(u""" - ((CellBody) (0 0 0 2)) + ((CellBody) (-1 0 0 2) (1 0 0 2)) ( (Color Blue) (Axon) (0 5 0 2) @@ -418,11 +418,7 @@ def test_has_multifurcation(): def test_has_unifurcation(): m = load_morphology(StringIO(u""" -("CellBody" - (Color Red) - (CellBody) - (0 0 0 2) - ) +((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 2) diff --git a/tests/check/test_morphtree.py b/tests/check/test_morphtree.py index ea4076d17..7137f49b7 100644 --- a/tests/check/test_morphtree.py +++ b/tests/check/test_morphtree.py @@ -41,8 +41,7 @@ def _generate_back_track_tree(n, dev): points = np.array(dev) + np.array([1, 3 if n == 0 else -3, 0]) m = load_morphology(StringIO(u""" - ((CellBody) - (0 0 0 0.4)) + ((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 0.4) diff --git a/tests/core/test_iter.py b/tests/core/test_iter.py index 01f3bc405..d7c86d08c 100644 --- a/tests/core/test_iter.py +++ b/tests/core/test_iter.py @@ -200,8 +200,7 @@ def test_iter_segments_pop(): def test_iter_segments_section(): sec = load_morphology(StringIO(u""" - ((CellBody) - (0 0 0 2)) + ((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (1 2 3 8) diff --git a/tests/data/h5/v1/ordering/sample.h5 b/tests/data/h5/v1/ordering/sample.h5 index 58dc61ce5..8f770a9be 100644 Binary files a/tests/data/h5/v1/ordering/sample.h5 and b/tests/data/h5/v1/ordering/sample.h5 differ diff --git a/tests/data/h5/v1/ordering/sample_mixed_tree_sections.h5 b/tests/data/h5/v1/ordering/sample_mixed_tree_sections.h5 index de692d4d9..5b9faf0e9 100644 Binary files a/tests/data/h5/v1/ordering/sample_mixed_tree_sections.h5 and b/tests/data/h5/v1/ordering/sample_mixed_tree_sections.h5 differ diff --git a/tests/features/test_get_features.py b/tests/features/test_get_features.py index b0ac1be0f..3a29d02b0 100644 --- a/tests/features/test_get_features.py +++ b/tests/features/test_get_features.py @@ -402,7 +402,7 @@ def test_segment_meander_angles(): def test_segment_meander_angles_single_section(): - m = nm.load_morphology(StringIO(u"""((CellBody) (0 0 0 0)) + m = nm.load_morphology(StringIO(u"""((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 2) (1 0 0 2) @@ -736,6 +736,7 @@ def test_sholl_frequency(): # check that if there is no neurite of a specific type, an empty list is returned assert features.get('sholl_frequency', m, neurite_type=NeuriteType.axon) == [] + assert features.get('sholl_frequency', m, neurite_type=NeuriteType.axon, distance_type='path') == [] def test_bifurcation_partitions(): diff --git a/tests/features/test_morphology.py b/tests/features/test_morphology.py index 2337c839a..d31c04a17 100644 --- a/tests/features/test_morphology.py +++ b/tests/features/test_morphology.py @@ -506,6 +506,16 @@ def test_sholl_analysis_custom(): assert (list(morphology.sholl_crossings(morph_A, center=center, radii=radii)) == [2, 2, 2, 2, 2, 2, 2, 2]) + assert list(morphology.sholl_frequency(morph_A)) == [2, 2, 2, 2, 2, 2, 2, 2] + + assert ( + list(morphology.sholl_frequency(morph_A, distance_type='path')) == [2, 2, 2, 2, 2, 2, 2, 2] + ) + + assert list( + morphology.sholl_crossings(morph_A, center=center, radii=radii, distance_type='path') + ) == [2, 2, 2, 2, 2, 2, 2, 2] + morph_B = load_swc("""\ 1 1 0 0 0 1. -1 2 3 0 0 0 1. 1 @@ -522,8 +532,16 @@ def test_sholl_analysis_custom(): 13 4 -51 -5 0 1. 9 14 4 -51 -10 0 1. 9 """) - assert (list(morphology.sholl_crossings(morph_B, center=center, radii=radii)) == - [2, 2, 2, 10, 10, 0, 0, 0]) + assert (list( + morphology.sholl_crossings(morph_B, center=center, radii=radii) + ) == [2, 2, 2, 10, 10, 0, 0, 0]) + + assert (list( + morphology.sholl_crossings(morph_B, center=center, radii=radii, distance_type='path') + ) == [2, 2, 2, 10, 10, 0, 0, 0]) + + assert list(morphology.sholl_frequency(morph_B)) == [2, 2, 2, 2, 10, 10] + assert list(morphology.sholl_frequency(morph_B, distance_type='path')) == [2, 2, 2, 2, 10, 10] morph_C = load_swc("""\ 1 1 0 0 0 1. -1 @@ -544,6 +562,15 @@ def test_sholl_analysis_custom(): assert (list(morphology.sholl_crossings(morph_C, center=center, radii=radii)) == [2, 2, 2, 2, 2, 2, 10, 10]) + assert (list( + morphology.sholl_crossings(morph_C, center=center, radii=radii, distance_type='path') + ) == [2, 2, 2, 2, 2, 2, 10, 10]) + + assert list(morphology.sholl_frequency(morph_C)) == [2, 2, 2, 2, 2, 2, 2, 10, 10] + assert (list( + morphology.sholl_frequency(morph_C, distance_type='path') + ) == [2, 2, 2, 2, 2, 2, 2, 10, 10]) + def test_extent_along_axis(): morph = load_swc(""" diff --git a/tests/features/test_section.py b/tests/features/test_section.py index d549410a9..6afb6defe 100644 --- a/tests/features/test_section.py +++ b/tests/features/test_section.py @@ -76,7 +76,7 @@ def test_segment_taper_rates(): def test_section_area(): - sec = load_morphology(StringIO(u"""((CellBody) (0 0 0 2)) + sec = load_morphology(StringIO(u"""((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 2) (1 0 0 2))"""), reader='asc').sections[SECTION_ID] @@ -85,7 +85,7 @@ def test_section_area(): def test_segment_areas(): - sec = load_morphology(StringIO(u"""((CellBody) (0 0 0 2)) + sec = load_morphology(StringIO(u"""((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 4) (1 0 0 4) @@ -95,7 +95,7 @@ def test_segment_areas(): def test_segment_volumes(): - sec = load_morphology(StringIO(u"""((CellBody) (0 0 0 2)) + sec = load_morphology(StringIO(u"""((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 4) (1 0 0 4) @@ -105,7 +105,7 @@ def test_segment_volumes(): def test_segment_mean_radii(): - sec = load_morphology(StringIO(u"""((CellBody) (0 0 0 2)) + sec = load_morphology(StringIO(u"""((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 2) (1 0 0 4) @@ -115,7 +115,7 @@ def test_segment_mean_radii(): def test_segment_midpoints(): - sec = load_morphology(StringIO(u"""((CellBody) (0 0 0 2)) + sec = load_morphology(StringIO(u"""((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 2) (1 0 0 4) @@ -126,7 +126,7 @@ def test_segment_midpoints(): def test_section_tortuosity(): sec_a = load_morphology(StringIO(u""" - ((CellBody) (0 0 0 2)) + ((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 2) (1 0 0 2) @@ -134,7 +134,7 @@ def test_section_tortuosity(): (3 0 0 2))"""), reader='asc').sections[SECTION_ID] sec_b = load_morphology(StringIO(u""" - ((CellBody) (0 0 0 2)) + ((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 2) (1 0 0 2) @@ -149,7 +149,7 @@ def test_section_tortuosity(): morphmath.section_length(s.points) / morphmath.point_dist(s.points[0], s.points[-1])) def test_setion_tortuosity_single_point(): - sec = load_morphology(StringIO(u"""((CellBody) (0 0 0 2)) + sec = load_morphology(StringIO(u"""((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (1 2 3 2))"""), reader='asc').sections[SECTION_ID] assert section.section_tortuosity(sec) == 1.0 @@ -157,7 +157,7 @@ def test_setion_tortuosity_single_point(): def test_section_tortuosity_looping_section(): sec = load_morphology(StringIO(u""" - ((CellBody) (0 0 0 2)) + ((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 2) (1 0 0 2) @@ -169,7 +169,7 @@ def test_section_tortuosity_looping_section(): def test_section_meander_angles(): - s0 = load_morphology(StringIO(u"""((CellBody) (0 0 0 0)) + s0 = load_morphology(StringIO(u"""((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 2) (1 0 0 2) @@ -178,7 +178,7 @@ def test_section_meander_angles(): (4 0 0 2))"""), reader='asc').sections[SECTION_ID] assert section.section_meander_angles(s0) == [math.pi, math.pi, math.pi] - s1 = load_morphology(StringIO(u"""((CellBody) (0 0 0 0)) + s1 = load_morphology(StringIO(u"""((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 2) (1 0 0 2) @@ -187,7 +187,7 @@ def test_section_meander_angles(): (2 2 0 2))"""), reader='asc').sections[SECTION_ID] assert section.section_meander_angles(s1) == [math.pi / 2, math.pi / 2, math.pi / 2] - s2 = load_morphology(StringIO(u"""((CellBody) (0 0 0 0)) + s2 = load_morphology(StringIO(u"""((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 2) (0 0 1 2) @@ -197,7 +197,7 @@ def test_section_meander_angles(): def test_section_meander_angles_single_segment(): - s = load_morphology(StringIO(u"""((CellBody) (0 0 0 0)) + s = load_morphology(StringIO(u"""((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 2) (1 1 1 2))"""), reader='asc').sections[SECTION_ID] @@ -212,7 +212,7 @@ def test_strahler_order(): def test_locate_segment_position(): - s = load_morphology(StringIO(u"""((CellBody) (0 0 0 0)) + s = load_morphology(StringIO(u"""((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 0) (3 0 4 200) @@ -231,8 +231,7 @@ def test_locate_segment_position(): def test_mean_radius(): n = load_morphology(StringIO(u""" - ((CellBody) - (0 0 0 1)) + ((CellBody) (-1 0 0 2) (1 0 0 2)) ((Dendrite) (0 0 0 0)