From c5b7da0451df0464102c62d01c69c7d9e682dd63 Mon Sep 17 00:00:00 2001 From: Nora Abi Akar <nora.abiakar@gmail.com> Date: Tue, 17 Nov 2020 18:07:55 +0100 Subject: [PATCH] Return morphology from `load_swc_xxx` (#1244) * Return `arb::morphology` from `load_swc_xxx` methods for consistency with the NeuroML interface. * Fix unit tests. * Fix documentation --- arborio/include/arborio/swcio.hpp | 20 +- arborio/swcio.cpp | 12 +- doc/concepts/morphology.rst | 8 +- doc/cpp/cable_cell.rst | 18 +- doc/python/morphology.rst | 122 +++---- example/ring/branch_cell.hpp | 2 +- example/single/single.cpp | 2 +- python/example/single_cell_swc.py | 4 +- python/morphology.cpp | 16 +- test/ubench/fvm_discretize.cpp | 2 +- test/unit/test_morphology.cpp | 6 +- test/unit/test_swcio.cpp | 515 +++++++++++++++++------------- 12 files changed, 402 insertions(+), 325 deletions(-) diff --git a/arborio/include/arborio/swcio.hpp b/arborio/include/arborio/swcio.hpp index cd7de765..2ec98a5e 100644 --- a/arborio/include/arborio/swcio.hpp +++ b/arborio/include/arborio/swcio.hpp @@ -5,7 +5,7 @@ #include <vector> #include <arbor/arbexcept.hpp> -#include <arbor/morph/segment_tree.hpp> +#include <arbor/morph/morphology.hpp> namespace arborio { @@ -162,25 +162,25 @@ public: swc_data parse_swc(std::istream&); swc_data parse_swc(const std::string&); -// Convert a valid, ordered sequence of SWC records to a morphological segment tree. +// Convert a valid, ordered sequence of SWC records into a morphology. // // Note that 'one-point soma' SWC files are explicitly not supported. // -// The generated segment tree will be contiguous. There will be one segment for -// each SWC record after the first: this record defines the tag and distal point -// of the segment, while the proximal point is taken from the parent record. +// The segments of the generated morphology will be contiguous. There will be +// one segment for each SWC record after the first: this record defines the tag +// and distal point of the segment, while the proximal point is taken from the +// parent record. -arb::segment_tree load_swc_arbor(const swc_data& data); +arb::morphology load_swc_arbor(const swc_data& data); -// As above, will convert a valid, ordered sequence of SWC records to a morphological -// segment tree. +// As above, will convert a valid, ordered sequence of SWC records into a morphology // // Note that 'one-point soma' SWC files are supported here // // These functions comply with inferred SWC rules from the Allen institute and Neuron. // These rules are explicitly listed in the docs. -arb::segment_tree load_swc_neuron(const swc_data& data); -arb::segment_tree load_swc_allen(const swc_data& data, bool no_gaps=false); +arb::morphology load_swc_neuron(const swc_data& data); +arb::morphology load_swc_allen(const swc_data& data, bool no_gaps=false); } // namespace arborio diff --git a/arborio/swcio.cpp b/arborio/swcio.cpp index 89c5ea31..abb72967 100644 --- a/arborio/swcio.cpp +++ b/arborio/swcio.cpp @@ -195,7 +195,7 @@ swc_data parse_swc(const std::string& text) { return parse_swc(is); } -arb::segment_tree load_swc_arbor(const swc_data& data) { +arb::morphology load_swc_arbor(const swc_data& data) { const auto& records = data.records(); if (records.empty()) return {}; @@ -237,10 +237,10 @@ arb::segment_tree load_swc_arbor(const swc_data& data) { throw swc_spherical_soma(first_id); } - return tree; + return arb::morphology(tree); } -arb::segment_tree load_swc_neuron(const swc_data& data) { +arb::morphology load_swc_neuron(const swc_data& data) { const auto& records = data.records(); // Assert that the file contains at least one sample. @@ -416,10 +416,10 @@ arb::segment_tree load_swc_neuron(const swc_data& data) { if (!unused_samples.empty()) { throw swc_single_sample_segment(*unused_samples.begin()); } - return tree; + return arb::morphology(tree); } -arb::segment_tree load_swc_allen(const swc_data& data, bool no_gaps) { +arb::morphology load_swc_allen(const swc_data& data, bool no_gaps) { auto records = data.records(); // Assert that the file contains at least one sample. @@ -506,7 +506,7 @@ arb::segment_tree load_swc_allen(const swc_data& data, bool no_gaps) { throw swc_single_sample_segment{*unused_samples.begin()}; } - return tree; + return arb::morphology(tree); } } // namespace arborio diff --git a/doc/concepts/morphology.rst b/doc/concepts/morphology.rst index 6bcad30a..6cddcc0f 100644 --- a/doc/concepts/morphology.rst +++ b/doc/concepts/morphology.rst @@ -495,8 +495,7 @@ Arbor supports reading morphologies described using the SWC files may contain comments, which are stored as metadata. A blank line anywhere in the file is interpreted as end of data. The description of the morphology is encoded as a list of samples with an id, an `x,y,z` location in space, a radius, a tag and a parent id. Arbor parses these samples, performs some checks, -then generates a :ref:`segment tree <morph-segment_tree>` describing the morphology according to one of three -possible interpretations. +then generates a morphology according to one of three possible interpretations. The SWC file format specifications are not very detailed, which has lead different simulators to interpret SWC files in different ways, especially when it comes to the soma. Arbor has its own an interpretation that @@ -517,9 +516,8 @@ is interpreted as a fork point in the morphology, and acts as the proximal point Arbor interpretation: """"""""""""""""""""" In addition to the previously listed checks, the arbor interpretation explicitly disallows SWC files where the soma is -described by a single sample. It constructs the soma from 2 or more samples that form 1 or more segments in the segment -tree. A *segment* is always constructed between a sample and its parent. This means that there are no gaps in the -resulting segment tree and morphology. +described by a single sample. It constructs the soma from 2 or more samples, forming 1 or more segments. A *segment* is +always constructed between a sample and its parent. This means that there are no gaps in the resulting morphology. Arbor has no magic rules or transformations for the soma. It can be a single branch or multiple branches; segments of a different tag can connect to its distal end, proximal end or anywhere in the middle. For example, to create a diff --git a/doc/cpp/cable_cell.rst b/doc/cpp/cable_cell.rst index 6e7d146f..6df7f504 100644 --- a/doc/cpp/cable_cell.rst +++ b/doc/cpp/cable_cell.rst @@ -160,7 +160,7 @@ has three different interpretation of that format. A :cpp:func:`parse_swc()` function is used to parse the SWC file and generate a :cpp:type:`swc_data` object. This object contains a vector of :cpp:type:`swc_record` objects that represent the SWC samples, with a number of basic checks performed on them. The :cpp:type:`swc_data` object can then be used to generate a -:cpp:type:`segment_tree` object using one of the following functions: (See the morphology concepts +:cpp:type:`morphology` object using one of the following functions: (See the morphology concepts :ref:`page <morph-formats>` for more details). * :cpp:func:`load_swc_arbor` @@ -211,19 +211,19 @@ basic checks performed on them. The :cpp:type:`swc_data` object can then be used Returns an :cpp:type:`swc_data` object given an std::istream object. -.. cpp:function:: segment_tree load_swc_arbor(const swc_data& data) +.. cpp:function:: morphology load_swc_arbor(const swc_data& data) - Returns a segment tree constructed according to Arbor's SWC specifications. + Returns a :cpp:type:`morphology` constructed according to Arbor's SWC specifications. -.. cpp:function:: segment_tree load_swc_allen(const swc_data& data, bool no_gaps=false) +.. cpp:function:: morphology load_swc_allen(const swc_data& data, bool no_gaps=false) - Returns a segment tree constructed according to the Allen Institute's SWC specifications. - By default, gaps in the segment tree are allowed, this can be toggled using the ``no_gaps`` - argument. + Returns a :cpp:type:`morphology` constructed according to the Allen Institute's SWC + specifications. By default, gaps in the morphology are allowed, this can be toggled + using the ``no_gaps`` argument. -.. cpp:function:: segment_tree load_swc_neuron(const swc_data& data) +.. cpp:function:: morphology load_swc_neuron(const swc_data& data) - Returns a segment tree constructed according to NEURON's SWC specifications. + Returns a :cpp:type:`morphology` constructed according to NEURON's SWC specifications. .. _locsets-and-regions: diff --git a/doc/python/morphology.rst b/doc/python/morphology.rst index 84a4d566..84252e93 100644 --- a/doc/python/morphology.rst +++ b/doc/python/morphology.rst @@ -262,9 +262,62 @@ Cell morphology A list of the segments. +.. py:class:: morphology + + A *morphology* describes the geometry of a cell as unbranched cables + with variable radius and their associated tree structure. + + .. note:: + A morphology takes a segment tree and construct the cable branches. + Meta data about branches and their properties that may be expensive to calculate + is stored for fast look up during later stages of model building, and + querying by users. + + For this reason, morphologies are read only. To change a morphology, a new + morphology should be created using a new segment tree. + + There is one *constructor* for a morphology: + + .. function:: morphology(segment_tree) + + Construct from a segment tree. + + The morphology provides an interface for querying morphology properties: + + .. attribute:: empty + :type: bool + + Indicates if the morphology is empty. + + .. attribute:: num_branches + :type: int + + The number of branches in the morphology. + + .. method:: branch_parent(i) + + The parent branch of a branch. + + :param int i: branch index + :rtype: int + + .. method:: branch_children(i) + + The child branches of a branch. + + :param int i: branch index + :rtype: list + + .. method:: branch_segments(i) + + A list of the segments in a branch, ordered from proximal to distal. + + :param int i: branch index + :rtype: list + .. py:function:: load_swc(filename) - Loads the morphology in an SWC file as a :class:`segment_tree` according to arbor's SWC specifications. + Loads the :class:`morphology` from an SWC file according to arbor's SWC specifications. (See the morphology concepts :ref:`page <morph-formats>` for more details). The samples in the SWC files are treated as the end points of segments, where a @@ -290,12 +343,13 @@ Cell morphology :param str filename: the name of the SWC file. - :rtype: segment_tree + :rtype: morphology .. py:function:: load_swc_neuron(filename) - Loads the morphology in an SWC file as a :class:`segment_tree` according to NEURON's SWC specifications. + Loads the :class:`morphology` from an SWC file according to NEURON's SWC specifications. Specifically: + * The first sample must be a soma sample. * The soma is represented by a series of n≥1 unbranched, serially listed samples. * The soma is constructed as a single cylinder with diameter equal to the piecewise average diameter of all the @@ -311,13 +365,13 @@ Cell morphology tag. :param str filename: the name of the SWC file. - :rtype: segment_tree + :rtype: morphology .. py:function:: load_swc_allen(filename, no_gaps=False) - Generate a segment tree from an SWC file following the rules prescribed by - AllenDB and Sonata. Specifically: + Loads the :class:`morphology` from an SWC file according to the AllenDB and Sonata's SWC specifications. + Specifically: * The first sample (the root) is treated as the center of the soma. * The morphology is translated such that the soma is centered at (0,0,0). @@ -337,58 +391,4 @@ Cell morphology :param str filename: the name of the SWC file. :param bool no_gaps: enforce that distance between soma center and branches attached to soma is the soma radius. - :rtype: segment_tree - -.. py:class:: morphology - - A *morphology* describes the geometry of a cell as unbranched cables - with variable radius and their associated tree structure. - - .. note:: - A morphology takes a segment tree and construct the cable branches. - Meta data about branches and their properties that may be expensive to calculate - is stored for fast look up during later stages of model building, and - querying by users. - - For this reason, morphologies are read only. To change a morphology, a new - morphology should be created using a new segment tree. - - There is one *constructor* for a morphology: - - .. function:: morphology(segment_tree) - - Construct from a segment tree. - - The morphology provides an interface for querying morphology properties: - - .. attribute:: empty - :type: bool - - Indicates if the morphology is empty. - - .. attribute:: num_branches - :type: int - - The number of branches in the morphology. - - .. method:: branch_parent(i) - - The parent branch of a branch. - - :param int i: branch index - :rtype: int - - .. method:: branch_children(i) - - The child branches of a branch. - - :param int i: branch index - :rtype: list - - .. method:: branch_segments(i) - - A list of the segments in a branch, ordered from proximal to distal. - - :param int i: branch index - :rtype: list - + :rtype: morphology \ No newline at end of file diff --git a/example/ring/branch_cell.hpp b/example/ring/branch_cell.hpp index 837a6569..ec8d0075 100644 --- a/example/ring/branch_cell.hpp +++ b/example/ring/branch_cell.hpp @@ -105,7 +105,7 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param arb::label_dict d; using arb::reg::tagged; - d.set("soma", tagged(stag)); + d.set("soma", tagged(stag)); d.set("dend", tagged(dtag)); arb::cable_cell cell(arb::morphology(tree), d); diff --git a/example/single/single.cpp b/example/single/single.cpp index 6d8e2d46..2510ed46 100644 --- a/example/single/single.cpp +++ b/example/single/single.cpp @@ -158,5 +158,5 @@ arb::morphology read_swc(const std::string& path) { std::ifstream f(path); if (!f) throw std::runtime_error("unable to open SWC file: "+path); - return arb::morphology(arborio::load_swc_arbor(arborio::parse_swc(f))); + return arborio::load_swc_arbor(arborio::parse_swc(f)); } diff --git a/python/example/single_cell_swc.py b/python/example/single_cell_swc.py index bdfb718a..05cd7aa8 100755 --- a/python/example/single_cell_swc.py +++ b/python/example/single_cell_swc.py @@ -23,7 +23,7 @@ if len(sys.argv) < 2: sys.exit(0) filename = sys.argv[1] -tree = arbor.load_swc(filename) +morpho = arbor.load_swc(filename) # Define the regions and locsets in the model. defs = {'soma': '(tag 1)', # soma has tag 1 in swc files. @@ -35,7 +35,7 @@ defs = {'soma': '(tag 1)', # soma has tag 1 in swc files. labels = arbor.label_dict(defs) # Combine morphology with region and locset definitions to make a cable cell. -cell = arbor.cable_cell(tree, labels) +cell = arbor.cable_cell(morpho, labels) print(cell.locations('"axon_end"')) diff --git a/python/morphology.cpp b/python/morphology.cpp index c7de23ac..6f6425e9 100644 --- a/python/morphology.cpp +++ b/python/morphology.cpp @@ -139,7 +139,7 @@ void register_morphology(pybind11::module& m) { .def("__str__", [](const arb::segment_tree& s) { return util::pprintf("<arbor.segment_tree:\n{}>", s);}); - // Function that creates a segment_tree from an swc file. + // Function that creates a morphology from an swc file. // Wraps calls to C++ functions arborio::parse_swc() and arborio::load_swc_arbor(). m.def("load_swc", [](std::string fname) { @@ -158,14 +158,14 @@ void register_morphology(pybind11::module& m) { } }, "filename"_a, - "Generate a segment tree from an SWC file following the rules prescribed by\n" - "Arbor. Specifically:\n" + "Generate a morphology from an SWC file following the rules prescribed by Arbor.\n" + "Specifically:\n" "* Single-segment somas are disallowed. These are usually interpreted as spherical somas\n" " and are a special case. This behavior is not allowed using this SWC loader.\n" "* There are no special rules related to somata. They can be one or multiple branches\n" " and other segments can connect anywhere along them.\n" "* A segment is always created between a sample and its parent, meaning there\n" - " are no gaps in the resulting segment tree."); + " are no gaps in the resulting morphology."); m.def("load_swc_allen", [](std::string fname, bool no_gaps=false) { @@ -186,8 +186,8 @@ void register_morphology(pybind11::module& m) { } }, "filename"_a, "no_gaps"_a=false, - "Generate a segment tree from an SWC file following the rules prescribed by\n" - "AllenDB and Sonata. Specifically:\n" + "Generate a morphology from an SWC file following the rules prescribed by AllenDB\n" + " and Sonata. Specifically:\n" "* The first sample (the root) is treated as the center of the soma.\n" "* The first morphology is translated such that the soma is centered at (0,0,0).\n" "* The first sample has tag 1 (soma).\n" @@ -221,8 +221,8 @@ void register_morphology(pybind11::module& m) { } }, "filename"_a, - "Generate a segment tree from an SWC file following the rules prescribed by\n" - "NEURON. Specifically:\n" + "Generate a morphology from an SWC file following the rules prescribed by NEURON.\n" + " Specifically:\n" "* The first sample must be a soma sample.\n" "* The soma is represented by a series of n≥1 unbranched, serially listed samples.\n" "* The soma is constructed as a single cylinder with diameter equal to the piecewise\n" diff --git a/test/ubench/fvm_discretize.cpp b/test/ubench/fvm_discretize.cpp index dbfa258c..6f554244 100644 --- a/test/ubench/fvm_discretize.cpp +++ b/test/ubench/fvm_discretize.cpp @@ -28,7 +28,7 @@ arb::morphology from_swc(const std::string& path) { std::ifstream in(path); if (!in) throw std::runtime_error("could not open "+path); - return morphology(arborio::load_swc_arbor(arborio::parse_swc(in))); + return arborio::load_swc_arbor(arborio::parse_swc(in)); } void run_cv_geom(benchmark::State& state) { diff --git a/test/unit/test_morphology.cpp b/test/unit/test_morphology.cpp index b4648c6e..4822413e 100644 --- a/test/unit/test_morphology.cpp +++ b/test/unit/test_morphology.cpp @@ -322,11 +322,7 @@ TEST(morphology, swc) { auto swc = arborio::parse_swc(fid); // Build a segmewnt_tree from swc samples. - auto sm = arborio::load_swc_arbor(swc); - EXPECT_EQ(5798u, sm.size()); // SWC data contains 5799 samples. - - // Test that the morphology contains the expected number of branches. - auto m = arb::morphology(sm); + auto m = arborio::load_swc_arbor(swc); EXPECT_EQ(221u, m.num_branches()); // 219 branches + 2 from divided soma. } #endif diff --git a/test/unit/test_swcio.cpp b/test/unit/test_swcio.cpp index 06c6397b..ea36cc00 100644 --- a/test/unit/test_swcio.cpp +++ b/test/unit/test_swcio.cpp @@ -271,28 +271,36 @@ TEST(swc_parser, arbor_complaint) { {7, 3, p4.x, p4.y, p4.z, p4.radius, 4} }; - segment_tree tree = load_swc_arbor(swc); - ASSERT_EQ(4u, tree.segments().size()); - - EXPECT_EQ(mnpos, tree.parents()[0]); - EXPECT_EQ(1, tree.segments()[0].tag); - EXPECT_EQ(p0, tree.segments()[0].prox); - EXPECT_EQ(p1, tree.segments()[0].dist); - - EXPECT_EQ(0u, tree.parents()[1]); - EXPECT_EQ(3, tree.segments()[1].tag); - EXPECT_EQ(p1, tree.segments()[1].prox); - EXPECT_EQ(p2, tree.segments()[1].dist); - - EXPECT_EQ(0u, tree.parents()[2]); - EXPECT_EQ(2, tree.segments()[2].tag); - EXPECT_EQ(p1, tree.segments()[2].prox); - EXPECT_EQ(p3, tree.segments()[2].dist); - - EXPECT_EQ(1u, tree.parents()[3]); - EXPECT_EQ(3, tree.segments()[3].tag); - EXPECT_EQ(p2, tree.segments()[3].prox); - EXPECT_EQ(p4, tree.segments()[3].dist); + auto morpho = load_swc_arbor(swc); + ASSERT_EQ(3u, morpho.num_branches()); + + EXPECT_EQ(mnpos, morpho.branch_parent(0)); + EXPECT_EQ(0u, morpho.branch_parent(1)); + EXPECT_EQ(0u, morpho.branch_parent(2)); + + auto segs_0 = morpho.branch_segments(0); + auto segs_1 = morpho.branch_segments(1); + auto segs_2 = morpho.branch_segments(2); + + EXPECT_EQ(1u, segs_0.size()); + EXPECT_EQ(2u, segs_1.size()); + EXPECT_EQ(1u, segs_2.size()); + + EXPECT_EQ(1, segs_0[0].tag); + EXPECT_EQ(p0, segs_0[0].prox); + EXPECT_EQ(p1, segs_0[0].dist); + + EXPECT_EQ(3, segs_1[0].tag); + EXPECT_EQ(p1, segs_1[0].prox); + EXPECT_EQ(p2, segs_1[0].dist); + + EXPECT_EQ(3, segs_1[1].tag); + EXPECT_EQ(p2, segs_1[1].prox); + EXPECT_EQ(p4, segs_1[1].dist); + + EXPECT_EQ(2, segs_2[0].tag); + EXPECT_EQ(p1, segs_2[0].prox); + EXPECT_EQ(p3, segs_2[0].dist); } { // Otherwise, ensure segment ends and tags correspond. @@ -308,23 +316,29 @@ TEST(swc_parser, arbor_complaint) { {4, 3, p3.x, p3.y, p3.z, p3.radius, 2}, }; - segment_tree tree = load_swc_arbor(swc); - ASSERT_EQ(3u, tree.segments().size()); + auto morpho = load_swc_arbor(swc); + ASSERT_EQ(2u, morpho.num_branches()); - EXPECT_EQ(mnpos, tree.parents()[0]); - EXPECT_EQ(1, tree.segments()[0].tag); - EXPECT_EQ(p0, tree.segments()[0].prox); - EXPECT_EQ(p1, tree.segments()[0].dist); + EXPECT_EQ(mnpos, morpho.branch_parent(0)); + EXPECT_EQ(mnpos, morpho.branch_parent(1)); - EXPECT_EQ(mnpos, tree.parents()[1]); - EXPECT_EQ(2, tree.segments()[1].tag); - EXPECT_EQ(p0, tree.segments()[1].prox); - EXPECT_EQ(p2, tree.segments()[1].dist); + auto segs_0 = morpho.branch_segments(0); + auto segs_1 = morpho.branch_segments(1); - EXPECT_EQ(0u, tree.parents()[2]); - EXPECT_EQ(3, tree.segments()[2].tag); - EXPECT_EQ(p1, tree.segments()[2].prox); - EXPECT_EQ(p3, tree.segments()[2].dist); + EXPECT_EQ(2u, segs_0.size()); + EXPECT_EQ(1u, segs_1.size()); + + EXPECT_EQ(1, segs_0[0].tag); + EXPECT_EQ(p0, segs_0[0].prox); + EXPECT_EQ(p1, segs_0[0].dist); + + EXPECT_EQ(3, segs_0[1].tag); + EXPECT_EQ(p1, segs_0[1].prox); + EXPECT_EQ(p3, segs_0[1].dist); + + EXPECT_EQ(2, segs_1[0].tag); + EXPECT_EQ(p0, segs_1[0].prox); + EXPECT_EQ(p2, segs_1[0].dist); } } @@ -363,17 +377,21 @@ TEST(swc_parser, allen_compliant) { std::vector<swc_record> swc{ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1} }; - segment_tree tree = load_swc_allen(swc); + auto morpho = load_swc_allen(swc); mpoint prox{p0.x, p0.y-p0.radius, p0.z, p0.radius}; mpoint dist{p0.x, p0.y+p0.radius, p0.z, p0.radius}; - ASSERT_EQ(1u, tree.segments().size()); + ASSERT_EQ(1u, morpho.num_branches()); + EXPECT_EQ(mnpos, morpho.branch_parent(0)); - EXPECT_EQ(mnpos, tree.parents()[0]); - EXPECT_EQ(1, tree.segments()[0].tag); - EXPECT_EQ(prox, tree.segments()[0].prox); - EXPECT_EQ(dist, tree.segments()[0].dist); + auto segs_0 = morpho.branch_segments(0); + + EXPECT_EQ(1u, segs_0.size()); + + EXPECT_EQ(1, segs_0[0].tag); + EXPECT_EQ(prox, segs_0[0].prox); + EXPECT_EQ(dist, segs_0[0].dist); } { // One-point soma, two-point dendrite @@ -386,22 +404,26 @@ TEST(swc_parser, allen_compliant) { {2, 3, p1.x, p1.y, p1.z, p1.radius, 1}, {3, 3, p2.x, p2.y, p2.z, p2.radius, 2} }; - segment_tree tree = load_swc_allen(swc); + auto morpho = load_swc_allen(swc); mpoint prox{0, -10, 0, 10}; mpoint dist{0, 10, 0, 10}; - ASSERT_EQ(2u, tree.segments().size()); + ASSERT_EQ(1u, morpho.num_branches()); + + EXPECT_EQ(mnpos, morpho.branch_parent(0)); - EXPECT_EQ(mnpos, tree.parents()[0]); - EXPECT_EQ(1, tree.segments()[0].tag); - EXPECT_EQ(prox, tree.segments()[0].prox); - EXPECT_EQ(dist, tree.segments()[0].dist); + auto segs_0 = morpho.branch_segments(0); - EXPECT_EQ(0u, tree.parents()[1]); - EXPECT_EQ(3, tree.segments()[1].tag); - EXPECT_EQ(p1, tree.segments()[1].prox); - EXPECT_EQ(p2, tree.segments()[1].dist); + EXPECT_EQ(2u, segs_0.size()); + + EXPECT_EQ(1, segs_0[0].tag); + EXPECT_EQ(prox, segs_0[0].prox); + EXPECT_EQ(dist, segs_0[0].dist); + + EXPECT_EQ(3, segs_0[1].tag); + EXPECT_EQ(p1, segs_0[1].prox); + EXPECT_EQ(p2, segs_0[1].dist); } { // 1-point soma, 2-point dendrite, 2-point axon @@ -418,27 +440,33 @@ TEST(swc_parser, allen_compliant) { {4, 2, p3.x, p3.y, p3.z, p3.radius, 1}, {5, 2, p4.x, p4.y, p4.z, p4.radius, 4} }; - segment_tree tree = load_swc_allen(swc); + auto morpho = load_swc_allen(swc); mpoint prox{0, -1, 0, 1}; mpoint dist{0, 1, 0, 1}; - ASSERT_EQ(3u, tree.segments().size()); + ASSERT_EQ(2u, morpho.num_branches()); + + EXPECT_EQ(mnpos, morpho.branch_parent(0)); + EXPECT_EQ(mnpos, morpho.branch_parent(1)); + + auto segs_0 = morpho.branch_segments(0); + auto segs_1 = morpho.branch_segments(1); + + EXPECT_EQ(2u, segs_0.size()); + EXPECT_EQ(1u, segs_1.size()); - EXPECT_EQ(mnpos, tree.parents()[0]); - EXPECT_EQ(1, tree.segments()[0].tag); - EXPECT_EQ(prox, tree.segments()[0].prox); - EXPECT_EQ(dist, tree.segments()[0].dist); + EXPECT_EQ(1, segs_0[0].tag); + EXPECT_EQ(prox, segs_0[0].prox); + EXPECT_EQ(dist, segs_0[0].dist); - EXPECT_EQ(0u, tree.parents()[1]); - EXPECT_EQ(3, tree.segments()[1].tag); - EXPECT_EQ(p1, tree.segments()[1].prox); - EXPECT_EQ(p2, tree.segments()[1].dist); + EXPECT_EQ(3, segs_0[1].tag); + EXPECT_EQ(p1, segs_0[1].prox); + EXPECT_EQ(p2, segs_0[1].dist); - EXPECT_EQ(mnpos, tree.parents()[2]); - EXPECT_EQ(2, tree.segments()[2].tag); - EXPECT_EQ(p3, tree.segments()[2].prox); - EXPECT_EQ(p4, tree.segments()[2].dist); + EXPECT_EQ(2, segs_1[0].tag); + EXPECT_EQ(p3, segs_1[0].prox); + EXPECT_EQ(p4, segs_1[0].dist); } } @@ -538,20 +566,25 @@ TEST(swc_parser, neuron_compliant) { std::vector<swc_record> swc{ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1} }; - segment_tree tree = load_swc_neuron(swc); + auto morpho = load_swc_neuron(swc); mpoint prox{p0.x, p0.y-p0.radius, p0.z, p0.radius}; mpoint dist{p0.x, p0.y+p0.radius, p0.z, p0.radius}; - ASSERT_EQ(1u, tree.segments().size()); + ASSERT_EQ(1u, morpho.num_branches()); - EXPECT_EQ(mnpos, tree.parents()[0]); - EXPECT_EQ(1, tree.segments()[0].tag); - EXPECT_EQ(prox, tree.segments()[0].prox); - EXPECT_EQ(dist, tree.segments()[0].dist); + EXPECT_EQ(mnpos, morpho.branch_parent(0)); + + auto segs_0 = morpho.branch_segments(0); + + EXPECT_EQ(1u, segs_0.size()); + + EXPECT_EQ(1, segs_0[0].tag); + EXPECT_EQ(prox, segs_0[0].prox); + EXPECT_EQ(dist, segs_0[0].dist); } { - // Two-point soma; interpretted as 1 + // Two-point soma; interpretted as 1 segment mpoint p0{0, 0, -10, 10}; mpoint p1{0, 0, 0, 10}; @@ -559,14 +592,19 @@ TEST(swc_parser, neuron_compliant) { {1, 1, p0.x, p0.y, p0.z, p0.radius, -1}, {2, 1, p1.x, p1.y, p1.z, p1.radius, 1} }; - segment_tree tree = load_swc_neuron(swc); + auto morpho = load_swc_neuron(swc); + + ASSERT_EQ(1u, morpho.num_branches()); + + EXPECT_EQ(mnpos, morpho.branch_parent(0)); + + auto segs_0 = morpho.branch_segments(0); - ASSERT_EQ(1u, tree.segments().size()); + EXPECT_EQ(1u, segs_0.size()); - EXPECT_EQ(mnpos, tree.parents()[0]); - EXPECT_EQ(1, tree.segments()[0].tag); - EXPECT_EQ(p0, tree.segments()[0].prox); - EXPECT_EQ(p1, tree.segments()[0].dist); + EXPECT_EQ(1, segs_0[0].tag); + EXPECT_EQ(p0, segs_0[0].prox); + EXPECT_EQ(p1, segs_0[0].dist); } { // Three-point soma; interpretted as 2 segments @@ -579,19 +617,23 @@ TEST(swc_parser, neuron_compliant) { {2, 1, p1.x, p1.y, p1.z, p1.radius, 1}, {3, 1, p2.x, p2.y, p2.z, p2.radius, 2} }; - segment_tree tree = load_swc_neuron(swc); + auto morpho = load_swc_neuron(swc); - ASSERT_EQ(2u, tree.segments().size()); + ASSERT_EQ(1u, morpho.num_branches()); - EXPECT_EQ(mnpos, tree.parents()[0]); - EXPECT_EQ(1, tree.segments()[0].tag); - EXPECT_EQ(p0, tree.segments()[0].prox); - EXPECT_EQ(p1, tree.segments()[0].dist); + EXPECT_EQ(mnpos, morpho.branch_parent(0)); - EXPECT_EQ(0u, tree.parents()[1]); - EXPECT_EQ(1, tree.segments()[1].tag); - EXPECT_EQ(p1, tree.segments()[1].prox); - EXPECT_EQ(p2, tree.segments()[1].dist); + auto segs_0 = morpho.branch_segments(0); + + EXPECT_EQ(2u, segs_0.size()); + + EXPECT_EQ(1, segs_0[0].tag); + EXPECT_EQ(p0, segs_0[0].prox); + EXPECT_EQ(p1, segs_0[0].dist); + + EXPECT_EQ(1, segs_0[1].tag); + EXPECT_EQ(p1, segs_0[1].prox); + EXPECT_EQ(p2, segs_0[1].dist); } { // 6-point soma; interpretted as 5 segments @@ -610,35 +652,35 @@ TEST(swc_parser, neuron_compliant) { {5, 1, p4.x, p4.y, p4.z, p4.radius, 4}, {6, 1, p5.x, p5.y, p5.z, p5.radius, 5} }; - segment_tree tree = load_swc_neuron(swc); + auto morpho = load_swc_neuron(swc); + + ASSERT_EQ(1u, morpho.num_branches()); + + EXPECT_EQ(mnpos, morpho.branch_parent(0)); + auto segs_0 = morpho.branch_segments(0); - ASSERT_EQ(5u, tree.segments().size()); + EXPECT_EQ(5u, segs_0.size()); - EXPECT_EQ(mnpos, tree.parents()[0]); - EXPECT_EQ(1, tree.segments()[0].tag); - EXPECT_EQ(p0, tree.segments()[0].prox); - EXPECT_EQ(p1, tree.segments()[0].dist); + EXPECT_EQ(1, segs_0[0].tag); + EXPECT_EQ(p0, segs_0[0].prox); + EXPECT_EQ(p1, segs_0[0].dist); - EXPECT_EQ(0u, tree.parents()[1]); - EXPECT_EQ(1, tree.segments()[1].tag); - EXPECT_EQ(p1, tree.segments()[1].prox); - EXPECT_EQ(p2, tree.segments()[1].dist); + EXPECT_EQ(1, segs_0[1].tag); + EXPECT_EQ(p1, segs_0[1].prox); + EXPECT_EQ(p2, segs_0[1].dist); - EXPECT_EQ(1u, tree.parents()[2]); - EXPECT_EQ(1, tree.segments()[2].tag); - EXPECT_EQ(p2, tree.segments()[2].prox); - EXPECT_EQ(p3, tree.segments()[2].dist); + EXPECT_EQ(1, segs_0[2].tag); + EXPECT_EQ(p2, segs_0[2].prox); + EXPECT_EQ(p3, segs_0[2].dist); - EXPECT_EQ(2u, tree.parents()[3]); - EXPECT_EQ(1, tree.segments()[3].tag); - EXPECT_EQ(p3, tree.segments()[3].prox); - EXPECT_EQ(p4, tree.segments()[3].dist); + EXPECT_EQ(1, segs_0[3].tag); + EXPECT_EQ(p3, segs_0[3].prox); + EXPECT_EQ(p4, segs_0[3].dist); - EXPECT_EQ(3u, tree.parents()[4]); - EXPECT_EQ(1, tree.segments()[4].tag); - EXPECT_EQ(p4, tree.segments()[4].prox); - EXPECT_EQ(p5, tree.segments()[4].dist); + EXPECT_EQ(1, segs_0[4].tag); + EXPECT_EQ(p4, segs_0[4].prox); + EXPECT_EQ(p5, segs_0[4].dist); } { // One-point soma, two-point dendrite @@ -651,27 +693,36 @@ TEST(swc_parser, neuron_compliant) { {2, 3, p1.x, p1.y, p1.z, p1.radius, 1}, {3, 3, p2.x, p2.y, p2.z, p2.radius, 2} }; - segment_tree tree = load_swc_neuron(swc); + auto morpho = load_swc_neuron(swc); mpoint prox{0, -10, 0, 10}; mpoint dist{0, 10, 0, 10}; - ASSERT_EQ(3u, tree.segments().size()); + ASSERT_EQ(3u, morpho.num_branches()); - EXPECT_EQ(mnpos, tree.parents()[0]); - EXPECT_EQ(1, tree.segments()[0].tag); - EXPECT_EQ(prox, tree.segments()[0].prox); - EXPECT_EQ(p0, tree.segments()[0].dist); + EXPECT_EQ(mnpos, morpho.branch_parent(0)); + EXPECT_EQ(0u, morpho.branch_parent(1)); + EXPECT_EQ(0u, morpho.branch_parent(2)); - EXPECT_EQ(0u, tree.parents()[1]); - EXPECT_EQ(1, tree.segments()[1].tag); - EXPECT_EQ(p0, tree.segments()[1].prox); - EXPECT_EQ(dist, tree.segments()[1].dist); + auto segs_0 = morpho.branch_segments(0); + auto segs_1 = morpho.branch_segments(1); + auto segs_2 = morpho.branch_segments(2); - EXPECT_EQ(0u, tree.parents()[2]); - EXPECT_EQ(3, tree.segments()[2].tag); - EXPECT_EQ(p1, tree.segments()[2].prox); - EXPECT_EQ(p2, tree.segments()[2].dist); + EXPECT_EQ(1u, segs_0.size()); + EXPECT_EQ(1u, segs_1.size()); + EXPECT_EQ(1u, segs_2.size()); + + EXPECT_EQ(1, segs_0[0].tag); + EXPECT_EQ(prox, segs_0[0].prox); + EXPECT_EQ(p0, segs_0[0].dist); + + EXPECT_EQ(1, segs_1[0].tag); + EXPECT_EQ(p0, segs_1[0].prox); + EXPECT_EQ(dist, segs_1[0].dist); + + EXPECT_EQ(3, segs_2[0].tag); + EXPECT_EQ(p1, segs_2[0].prox); + EXPECT_EQ(p2, segs_2[0].dist); } { // 6-point soma, 2-point dendrite @@ -694,46 +745,51 @@ TEST(swc_parser, neuron_compliant) { {7, 3, p6.x, p6.y, p6.z, p6.radius, 6}, {8, 3, p7.x, p7.y, p7.z, p7.radius, 7} }; - segment_tree tree = load_swc_neuron(swc); + auto morpho = load_swc_neuron(swc); mpoint mid {0, 0, 5, 2.25}; - ASSERT_EQ(7u, tree.segments().size()); - - EXPECT_EQ(mnpos, tree.parents()[0]); - EXPECT_EQ(1, tree.segments()[0].tag); - EXPECT_EQ(p0, tree.segments()[0].prox); - EXPECT_EQ(p1, tree.segments()[0].dist); - - EXPECT_EQ(0u, tree.parents()[1]); - EXPECT_EQ(1, tree.segments()[1].tag); - EXPECT_EQ(p1, tree.segments()[1].prox); - EXPECT_EQ(p2, tree.segments()[1].dist); - - EXPECT_EQ(1u, tree.parents()[2]); - EXPECT_EQ(1, tree.segments()[2].tag); - EXPECT_EQ(p2, tree.segments()[2].prox); - EXPECT_EQ(mid, tree.segments()[2].dist); - - EXPECT_EQ(2u, tree.parents()[3]); - EXPECT_EQ(1, tree.segments()[3].tag); - EXPECT_EQ(mid, tree.segments()[3].prox); - EXPECT_EQ(p3, tree.segments()[3].dist); - - EXPECT_EQ(3u, tree.parents()[4]); - EXPECT_EQ(1, tree.segments()[4].tag); - EXPECT_EQ(p3, tree.segments()[4].prox); - EXPECT_EQ(p4, tree.segments()[4].dist); - - EXPECT_EQ(4u, tree.parents()[5]); - EXPECT_EQ(1, tree.segments()[5].tag); - EXPECT_EQ(p4, tree.segments()[5].prox); - EXPECT_EQ(p5, tree.segments()[5].dist); - - EXPECT_EQ(2u, tree.parents()[6]); - EXPECT_EQ(3, tree.segments()[6].tag); - EXPECT_EQ(p6, tree.segments()[6].prox); - EXPECT_EQ(p7, tree.segments()[6].dist); + ASSERT_EQ(3u, morpho.num_branches()); + + EXPECT_EQ(mnpos, morpho.branch_parent(0)); + EXPECT_EQ(0u, morpho.branch_parent(1)); + EXPECT_EQ(0u, morpho.branch_parent(2)); + + auto segs_0 = morpho.branch_segments(0); + auto segs_1 = morpho.branch_segments(1); + auto segs_2 = morpho.branch_segments(2); + + EXPECT_EQ(3u, segs_0.size()); + EXPECT_EQ(3u, segs_1.size()); + EXPECT_EQ(1u, segs_2.size()); + + EXPECT_EQ(1, segs_0[0].tag); + EXPECT_EQ(p0, segs_0[0].prox); + EXPECT_EQ(p1, segs_0[0].dist); + + EXPECT_EQ(1, segs_0[1].tag); + EXPECT_EQ(p1, segs_0[1].prox); + EXPECT_EQ(p2, segs_0[1].dist); + + EXPECT_EQ(1, segs_0[2].tag); + EXPECT_EQ(p2, segs_0[2].prox); + EXPECT_EQ(mid,segs_0[2].dist); + + EXPECT_EQ(1, segs_1[0].tag); + EXPECT_EQ(mid,segs_1[0].prox); + EXPECT_EQ(p3, segs_1[0].dist); + + EXPECT_EQ(1, segs_1[1].tag); + EXPECT_EQ(p3, segs_1[1].prox); + EXPECT_EQ(p4, segs_1[1].dist); + + EXPECT_EQ(1, segs_1[2].tag); + EXPECT_EQ(p4, segs_1[2].prox); + EXPECT_EQ(p5, segs_1[2].dist); + + EXPECT_EQ(3, segs_2[0].tag); + EXPECT_EQ(p6, segs_2[0].prox); + EXPECT_EQ(p7, segs_2[0].dist); } { // Two-point soma, two-point dendrite @@ -748,26 +804,35 @@ TEST(swc_parser, neuron_compliant) { {3, 3, p2.x, p2.y, p2.z, p2.radius, 2}, {4, 3, p3.x, p3.y, p3.z, p3.radius, 3} }; - segment_tree tree = load_swc_neuron(swc); + auto morpho = load_swc_neuron(swc); mpoint mid{0, 0, -10, 7}; - ASSERT_EQ(3u, tree.segments().size()); + ASSERT_EQ(3u, morpho.num_branches()); + + EXPECT_EQ(mnpos, morpho.branch_parent(0)); + EXPECT_EQ(0u, morpho.branch_parent(1)); + EXPECT_EQ(0u, morpho.branch_parent(2)); - EXPECT_EQ(mnpos, tree.parents()[0]); - EXPECT_EQ(1, tree.segments()[0].tag); - EXPECT_EQ(p0, tree.segments()[0].prox); - EXPECT_EQ(mid, tree.segments()[0].dist); + auto segs_0 = morpho.branch_segments(0); + auto segs_1 = morpho.branch_segments(1); + auto segs_2 = morpho.branch_segments(2); - EXPECT_EQ(0u, tree.parents()[1]); - EXPECT_EQ(1, tree.segments()[1].tag); - EXPECT_EQ(mid, tree.segments()[1].prox); - EXPECT_EQ(p1, tree.segments()[1].dist); + EXPECT_EQ(1u, segs_0.size()); + EXPECT_EQ(1u, segs_1.size()); + EXPECT_EQ(1u, segs_2.size()); - EXPECT_EQ(0u, tree.parents()[2]); - EXPECT_EQ(3, tree.segments()[2].tag); - EXPECT_EQ(p2, tree.segments()[2].prox); - EXPECT_EQ(p3, tree.segments()[2].dist); + EXPECT_EQ(1, segs_0[0].tag); + EXPECT_EQ(p0, segs_0[0].prox); + EXPECT_EQ(mid, segs_0[0].dist); + + EXPECT_EQ(1, segs_1[0].tag); + EXPECT_EQ(mid, segs_1[0].prox); + EXPECT_EQ(p1, segs_1[0].dist); + + EXPECT_EQ(3, segs_2[0].tag); + EXPECT_EQ(p2, segs_2[0].prox); + EXPECT_EQ(p3, segs_2[0].dist); } { // 2-point soma; 2-point dendrite; 1-point axon connected to the proximal end of the dendrite @@ -784,31 +849,42 @@ TEST(swc_parser, neuron_compliant) { {4, 3, p3.x, p3.y, p3.z, p3.radius, 3}, {5, 2, p4.x, p4.y, p4.z, p4.radius, 3} }; - segment_tree tree = load_swc_neuron(swc); + auto morpho = load_swc_neuron(swc); mpoint mid{0, 0, -7.5, 6.5}; - ASSERT_EQ(4u, tree.segments().size()); + ASSERT_EQ(4u, morpho.num_branches()); + + EXPECT_EQ(mnpos, morpho.branch_parent(0)); + EXPECT_EQ(0u, morpho.branch_parent(1)); + EXPECT_EQ(0u, morpho.branch_parent(2)); + EXPECT_EQ(0u, morpho.branch_parent(3)); - EXPECT_EQ(mnpos, tree.parents()[0]); - EXPECT_EQ(1, tree.segments()[0].tag); - EXPECT_EQ(p0, tree.segments()[0].prox); - EXPECT_EQ(mid, tree.segments()[0].dist); + auto segs_0 = morpho.branch_segments(0); + auto segs_1 = morpho.branch_segments(1); + auto segs_2 = morpho.branch_segments(2); + auto segs_3 = morpho.branch_segments(3); - EXPECT_EQ(0u, tree.parents()[1]); - EXPECT_EQ(1, tree.segments()[1].tag); - EXPECT_EQ(mid, tree.segments()[1].prox); - EXPECT_EQ(p1, tree.segments()[1].dist); + EXPECT_EQ(1u, segs_0.size()); + EXPECT_EQ(1u, segs_1.size()); + EXPECT_EQ(1u, segs_2.size()); + EXPECT_EQ(1u, segs_3.size()); - EXPECT_EQ(0u, tree.parents()[2]); - EXPECT_EQ(3, tree.segments()[2].tag); - EXPECT_EQ(p2, tree.segments()[2].prox); - EXPECT_EQ(p3, tree.segments()[2].dist); + EXPECT_EQ(1, segs_0[0].tag); + EXPECT_EQ(p0, segs_0[0].prox); + EXPECT_EQ(mid, segs_0[0].dist); - EXPECT_EQ(0u, tree.parents()[3]); - EXPECT_EQ(2, tree.segments()[3].tag); - EXPECT_EQ(p2, tree.segments()[3].prox); - EXPECT_EQ(p4, tree.segments()[3].dist); + EXPECT_EQ(1, segs_1[0].tag); + EXPECT_EQ(mid, segs_1[0].prox); + EXPECT_EQ(p1, segs_1[0].dist); + + EXPECT_EQ(3, segs_2[0].tag); + EXPECT_EQ(p2, segs_2[0].prox); + EXPECT_EQ(p3, segs_2[0].dist); + + EXPECT_EQ(2, segs_3[0].tag); + EXPECT_EQ(p2, segs_3[0].prox); + EXPECT_EQ(p4, segs_3[0].dist); } { // 2-point soma, 2-point dendrite, 2-point axon @@ -827,36 +903,43 @@ TEST(swc_parser, neuron_compliant) { {5, 2, p4.x, p4.y, p4.z, p4.radius, 4}, {6, 2, p5.x, p5.y, p5.z, p5.radius, 5} }; - segment_tree tree = load_swc_neuron(swc); + auto morpho = load_swc_neuron(swc); mpoint mid{0, 0, 4.5, 1.5}; - ASSERT_EQ(5u, tree.segments().size()); + ASSERT_EQ(3u, morpho.num_branches()); + + EXPECT_EQ(mnpos, morpho.branch_parent(0)); + EXPECT_EQ(0u, morpho.branch_parent(1)); + EXPECT_EQ(0u, morpho.branch_parent(2)); + + auto segs_0 = morpho.branch_segments(0); + auto segs_1 = morpho.branch_segments(1); + auto segs_2 = morpho.branch_segments(2); + + EXPECT_EQ(1u, segs_0.size()); + EXPECT_EQ(1u, segs_1.size()); + EXPECT_EQ(3u, segs_2.size()); - EXPECT_EQ(mnpos, tree.parents()[0]); - EXPECT_EQ(1, tree.segments()[0].tag); - EXPECT_EQ(p0, tree.segments()[0].prox); - EXPECT_EQ(mid, tree.segments()[0].dist); + EXPECT_EQ(1, segs_0[0].tag); + EXPECT_EQ(p0, segs_0[0].prox); + EXPECT_EQ(mid, segs_0[0].dist); - EXPECT_EQ(0u, tree.parents()[1]); - EXPECT_EQ(1, tree.segments()[1].tag); - EXPECT_EQ(mid, tree.segments()[1].prox); - EXPECT_EQ(p1, tree.segments()[1].dist); + EXPECT_EQ(1, segs_1[0].tag); + EXPECT_EQ(mid, segs_1[0].prox); + EXPECT_EQ(p1, segs_1[0].dist); - EXPECT_EQ(0u, tree.parents()[2]); - EXPECT_EQ(3, tree.segments()[2].tag); - EXPECT_EQ(p2, tree.segments()[2].prox); - EXPECT_EQ(p3, tree.segments()[2].dist); + EXPECT_EQ(3, segs_2[0].tag); + EXPECT_EQ(p2, segs_2[0].prox); + EXPECT_EQ(p3, segs_2[0].dist); - EXPECT_EQ(2u, tree.parents()[3]); - EXPECT_EQ(2, tree.segments()[3].tag); - EXPECT_EQ(p3, tree.segments()[3].prox); - EXPECT_EQ(p4, tree.segments()[3].dist); + EXPECT_EQ(2, segs_2[1].tag); + EXPECT_EQ(p3, segs_2[1].prox); + EXPECT_EQ(p4, segs_2[1].dist); - EXPECT_EQ(3u, tree.parents()[4]); - EXPECT_EQ(2, tree.segments()[4].tag); - EXPECT_EQ(p4, tree.segments()[4].prox); - EXPECT_EQ(p5, tree.segments()[4].dist); + EXPECT_EQ(2, segs_2[2].tag); + EXPECT_EQ(p4, segs_2[2].prox); + EXPECT_EQ(p5, segs_2[2].dist); } } -- GitLab