diff --git a/arbor/include/arbor/morph/isometry.hpp b/arbor/include/arbor/morph/isometry.hpp new file mode 100644 index 0000000000..cb06ee28dd --- /dev/null +++ b/arbor/include/arbor/morph/isometry.hpp @@ -0,0 +1,62 @@ +#pragma once + +#include +#include + +namespace arb { +// Represent a 3-d isometry as a rotation (via quaterion q_) +// and a subsequent translation (tx_, ty_, tz_). +struct ARB_ARBOR_API isometry { + + isometry() = default; + + // Composition: rotations are interpreted as applied to intrinsic + // coordinates (composed by right multiplication), and translations + // as applied to absolute coordinates (composed by addition). + + friend isometry operator*(const isometry& a, const isometry& b) { + return isometry(b.q_*a.q_, a.tx_+b.tx_, a.ty_+b.ty_, a.tz_+b.tz_); + } + + template + PointLike apply(PointLike p) const { + auto w = quaternion(p.x, p.y, p.z)^q_; + p.x = w.x+tx_; + p.y = w.y+ty_; + p.z = w.z+tz_; + return p; + } + +private: + using quaternion = arb::math::quaternion; + quaternion q_{1, 0, 0, 0}; + double tx_ = 0, ty_ = 0, tz_ = 0; + + isometry(quaternion q, double tx, double ty, double tz): + q_(std::move(q)), tx_(tx), ty_(ty), tz_(tz) {} + +public: + static isometry translate(double x, double y, double z) { + return isometry(quaternion{1, 0, 0, 0}, x, y, z); + } + + template + static isometry translate(const PointLike& p) { + return translate(p.x, p.y, p.z); + } + + // Rotate about axis in direction (ax, ay, az) by theta radians. + static isometry rotate(double theta, double ax, double ay, double az) { + double norm = std::sqrt(ax*ax+ay*ay+az*az); + double vscale = std::sin(theta/2)/norm; + + return isometry(quaternion{std::cos(theta/2), ax*vscale, ay*vscale, az*vscale}, 0, 0, 0); + } + + template + static isometry rotate(double theta, const PointLike& p) { + return rotate(theta, p.x, p.y, p.z); + } +}; + +} // arb diff --git a/arbor/include/arbor/morph/place_pwlin.hpp b/arbor/include/arbor/morph/place_pwlin.hpp index 09cc179dc7..6088a8ec07 100644 --- a/arbor/include/arbor/morph/place_pwlin.hpp +++ b/arbor/include/arbor/morph/place_pwlin.hpp @@ -10,64 +10,10 @@ #include #include #include -#include +#include -namespace arb { - -struct isometry { - // Represent a 3-d isometry as a rotation (via quaterion q_) - // and a subsequent translation (tx_, ty_, tz_). - - isometry() = default; - - // Composition: rotations are interpreted as applied to intrinsic - // coordinates (composed by right multiplication), and translations - // as applied to absolute coordinates (composed by addition). - - friend isometry operator*(const isometry& a, const isometry& b) { - return isometry(b.q_*a.q_, a.tx_+b.tx_, a.ty_+b.ty_, a.tz_+b.tz_); - } - template - PointLike apply(PointLike p) const { - auto w = quaternion(p.x, p.y, p.z)^q_; - p.x = w.x+tx_; - p.y = w.y+ty_; - p.z = w.z+tz_; - return p; - } - -private: - using quaternion = arb::math::quaternion; - quaternion q_{1, 0, 0, 0}; - double tx_ = 0, ty_ = 0, tz_ = 0; - - isometry(quaternion q, double tx, double ty, double tz): - q_(std::move(q)), tx_(tx), ty_(ty), tz_(tz) {} - -public: - static isometry translate(double x, double y, double z) { - return isometry(quaternion{1, 0, 0, 0}, x, y, z); - } - - template - static isometry translate(const PointLike& p) { - return translate(p.x, p.y, p.z); - } - - // Rotate about axis in direction (ax, ay, az) by theta radians. - static isometry rotate(double theta, double ax, double ay, double az) { - double norm = std::sqrt(ax*ax+ay*ay+az*az); - double vscale = std::sin(theta/2)/norm; - - return isometry(quaternion{std::cos(theta/2), ax*vscale, ay*vscale, az*vscale}, 0, 0, 0); - } - - template - static isometry rotate(double theta, const PointLike& p) { - return rotate(theta, p.x, p.y, p.z); - } -}; +namespace arb { struct place_pwlin_data; diff --git a/arbor/include/arbor/morph/primitives.hpp b/arbor/include/arbor/morph/primitives.hpp index cdb6472acb..4be04670b1 100644 --- a/arbor/include/arbor/morph/primitives.hpp +++ b/arbor/include/arbor/morph/primitives.hpp @@ -22,17 +22,11 @@ constexpr msize_t mnpos = msize_t(-1); struct ARB_SYMBOL_VISIBLE mpoint { double x, y, z; // [µm] double radius; // [μm] - - friend bool operator==(const mpoint& l, const mpoint& r); friend std::ostream& operator<<(std::ostream&, const mpoint&); - friend bool operator==(const mpoint& a, const mpoint& b) { - return a.x==b.x && a.y==b.y && a.z==b.z && a.radius==b.radius; - } - friend bool operator!=(const mpoint& a, const mpoint& b) { - return !(a==b); - } }; +ARB_DEFINE_LEXICOGRAPHIC_ORDERING(mpoint, (a.x,a.y,a.z,a.radius), (b.x,b.y,b.z,b.radius)); + ARB_ARBOR_API mpoint lerp(const mpoint& a, const mpoint& b, double u); ARB_ARBOR_API bool is_collocated(const mpoint& a, const mpoint& b); ARB_ARBOR_API double distance(const mpoint& a, const mpoint& b); @@ -55,6 +49,7 @@ struct ARB_SYMBOL_VISIBLE msegment { friend std::ostream& operator<<(std::ostream&, const msegment&); }; +ARB_DEFINE_LEXICOGRAPHIC_ORDERING(msegment, (a.id,a.prox,a.dist,a.tag), (b.id,b.prox,b.dist,b.tag)); // Describe a specific location on a morpholology. struct ARB_SYMBOL_VISIBLE mlocation { diff --git a/arbor/include/arbor/morph/segment_tree.hpp b/arbor/include/arbor/morph/segment_tree.hpp index e934f102eb..7d59be59ae 100644 --- a/arbor/include/arbor/morph/segment_tree.hpp +++ b/arbor/include/arbor/morph/segment_tree.hpp @@ -7,6 +7,7 @@ #include #include +#include namespace arb { @@ -52,8 +53,36 @@ class ARB_ARBOR_API segment_tree { bool is_root(msize_t i) const; friend std::ostream& operator<<(std::ostream&, const segment_tree&); + + // compare two trees for _identity_, not _equivalence_ + friend bool operator==(const segment_tree& l, const segment_tree& r) { + return (l.size() == r.size()) && (l.parents() == r.parents()) && (l.segments() == r.segments()); + } + + // apply isometry by mapping over internal state + friend segment_tree apply(const segment_tree&, const isometry&); }; -} // namesapce arb +// Split a segment_tree T into two subtrees such that R is the subtree +// of T that starts at the given id and L is T without R. +ARB_ARBOR_API std::pair +split_at(const segment_tree&, msize_t); +// Join two subtrees L and R at a given id in L, such that `join_at` is inverse +// to `split_at` for a proper choice of id. +ARB_ARBOR_API segment_tree +join_at(const segment_tree&, msize_t, const segment_tree&); +// Trees are equivalent if +// 1. the current segments' prox and dist points and their tags are identical. +// 2. all sub-trees starting at the current segment are equivalent. +// Note that orderdoes *not* matter in opposition to ==. +ARB_ARBOR_API bool +equivalent(const segment_tree& a, + const segment_tree& b); + +// Apply isometry +ARB_ARBOR_API segment_tree +apply(const segment_tree&, const isometry&); + +} // namesapce arb diff --git a/arbor/morph/segment_tree.cpp b/arbor/morph/segment_tree.cpp index 70498e4709..3a10a015b5 100644 --- a/arbor/morph/segment_tree.cpp +++ b/arbor/morph/segment_tree.cpp @@ -1,6 +1,6 @@ #include #include -#include +#include #include #include @@ -12,6 +12,119 @@ namespace arb { +struct node_t { + msize_t parent; + msize_t id; +}; + +using node_p = std::function; + +node_p yes = [](const node_t&) { return true; }; + +// invert parent <*> child relation, returns a map of parent_id -> [children_id] +// For predictable ordering we sort the vectors. +std::map> tree_to_children(const segment_tree& tree) { + const auto& parents = tree.parents(); + std::map> result; + for (auto ix = 0; ix < tree.size(); ++ix) result[parents[ix]].push_back(ix); + for (auto& [k, v]: result) std::sort(v.begin(), v.end()); + return result; +} + +// Copy a segment tree into a new tree +// - tree to be (partially) copied +// - start={parent, id}: start by attaching segment=`id`` from `tree` to the +// output at `parent`, then its children to it recursively +// - predicate: if returning false for a given node, we prune that sub-tree starting +// at node (inclusive). Can be used to prune trees by parent or id. +// - init: initial tree to append to +// Note: this is basically a recursive function w/ an explicit stack. +segment_tree copy_if(const segment_tree& tree, + const node_t& start, + node_p predicate, + const segment_tree& init={}) { + auto children_of = tree_to_children(tree); + auto& segments = tree.segments(); + segment_tree result = init; + auto todo = std::vector{start}; + while (!todo.empty()) { + auto node = todo.back(); + todo.pop_back(); + if (!predicate(node)) continue; + const auto& segment = segments[node.id]; + auto current = result.append(node.parent, segment.prox, segment.dist, segment.tag); + for (auto child: children_of[node.id]) { + todo.push_back({current, child}); + } + } + return result; +} + +std::pair +split_at(const segment_tree& tree, msize_t at) { + if (at >= tree.size() || at == mnpos) throw invalid_segment_parent(at, tree.size()); + // span the sub-tree starting at the splitting node + segment_tree post = copy_if(tree, {mnpos, at}, yes); + // copy the original tree, but skip all nodes in the `post` tree + segment_tree pre = copy_if(tree, + {mnpos, 0}, + [=](auto& node) { return node.id != at; }); + return {pre, post}; +} + +segment_tree +join_at(const segment_tree& lhs, msize_t at, const segment_tree& rhs) { + if (at >= lhs.size() && at != mnpos) throw invalid_segment_parent(at, lhs.size()); + return copy_if(rhs, {at, 0}, yes, lhs); +} + +bool +equivalent(const segment_tree& a, + const segment_tree& b) { + if(a.size() != b.size()) return false; + + auto + a_children_of = tree_to_children(a), + b_children_of = tree_to_children(b); + + auto fetch_children = [&](auto cursor, const auto& segments, auto& children_of) { + std::vector segs; + for (auto ix: children_of[cursor]) segs.push_back(segments[ix]); + std::sort(segs.begin(), segs.end(), + [](auto l, auto r) { + l.id = r.id = 0; + return l < r; + }); + return segs; + }; + + std::vector> todo{{arb::mnpos, arb::mnpos}}; + while (!todo.empty()) { + const auto& [a_cursor, b_cursor] = todo.back(); + auto as = fetch_children(a_cursor, a.segments(), a_children_of); + auto bs = fetch_children(b_cursor, b.segments(), b_children_of); + todo.pop_back(); + if (as.size() != bs.size()) return false; + for (auto ix = 0; ix < as.size(); ++ix) { + if ((as[ix].prox != bs[ix].prox) || + (as[ix].dist != bs[ix].dist) || + (as[ix].tag != bs[ix].tag)) return false; + todo.emplace_back(as[ix].id, bs[ix].id); + } + } + return true; +} + +segment_tree +apply(const segment_tree& tree, const isometry& iso) { + auto result = tree; + for (auto& seg: result.segments_) { + seg.prox = iso.apply(seg.prox); + seg.dist = iso.apply(seg.dist); + } + return result; +} + void segment_tree::reserve(msize_t n) { segments_.reserve(n); parents_.reserve(n); diff --git a/arborio/include/arborio/neurolucida.hpp b/arborio/include/arborio/neurolucida.hpp index cf82961380..3fd4d1fc8f 100644 --- a/arborio/include/arborio/neurolucida.hpp +++ b/arborio/include/arborio/neurolucida.hpp @@ -33,6 +33,9 @@ struct ARB_SYMBOL_VISIBLE asc_unsupported: asc_exception { }; struct asc_morphology { + // Raw segment tree from ASC, identical to morphology. + arb::segment_tree segment_tree; + // Morphology constructed from asc description. arb::morphology morphology; diff --git a/arborio/include/arborio/neuroml.hpp b/arborio/include/arborio/neuroml.hpp index 27c2663b2e..3646548e6f 100644 --- a/arborio/include/arborio/neuroml.hpp +++ b/arborio/include/arborio/neuroml.hpp @@ -71,7 +71,7 @@ struct nml_morphology_data { // Morphology id. std::string id; - // Morphology constructed from a signle NeuroML element. + // Morphology constructed from a single NeuroML element. arb::morphology morphology; // One region expression for each segment id. diff --git a/arborio/include/arborio/swcio.hpp b/arborio/include/arborio/swcio.hpp index 26189571ba..38b59e0c5a 100644 --- a/arborio/include/arborio/swcio.hpp +++ b/arborio/include/arborio/swcio.hpp @@ -128,6 +128,7 @@ ARB_ARBORIO_API swc_data parse_swc(const std::string&); // parent record. ARB_ARBORIO_API arb::morphology load_swc_arbor(const swc_data& data); +ARB_ARBORIO_API arb::segment_tree load_swc_arbor_raw(const swc_data& data); // As above, will convert a valid, ordered sequence of SWC records into a morphology // @@ -136,5 +137,6 @@ ARB_ARBORIO_API arb::morphology load_swc_arbor(const swc_data& data); // Complies inferred SWC rules from NEURON, explicitly listed in the docs. ARB_ARBORIO_API arb::morphology load_swc_neuron(const swc_data& data); +ARB_ARBORIO_API arb::segment_tree load_swc_neuron_raw(const swc_data& data); } // namespace arborio diff --git a/arborio/neurolucida.cpp b/arborio/neurolucida.cpp index af033796b0..832a7d2451 100644 --- a/arborio/neurolucida.cpp +++ b/arborio/neurolucida.cpp @@ -778,7 +778,7 @@ asc_morphology parse_asc_string(const char* input) { labels.set("dend", arb::reg::tagged(3)); labels.set("apic", arb::reg::tagged(4)); - return {std::move(morphology), std::move(labels)}; + return {stree, std::move(morphology), std::move(labels)}; } ARB_ARBORIO_API asc_morphology load_asc(std::string filename) { diff --git a/arborio/swcio.cpp b/arborio/swcio.cpp index 9fc2f2613d..34a0fec657 100644 --- a/arborio/swcio.cpp +++ b/arborio/swcio.cpp @@ -160,7 +160,7 @@ ARB_ARBORIO_API swc_data parse_swc(const std::string& text) { return parse_swc(is); } -ARB_ARBORIO_API arb::morphology load_swc_arbor(const swc_data& data) { +ARB_ARBORIO_API arb::segment_tree load_swc_arbor_raw(const swc_data& data) { const auto& records = data.records(); if (records.empty()) return {}; @@ -202,10 +202,10 @@ ARB_ARBORIO_API arb::morphology load_swc_arbor(const swc_data& data) { throw swc_spherical_soma(first_id); } - return arb::morphology(tree); + return tree; } -ARB_ARBORIO_API arb::morphology load_swc_neuron(const swc_data& data) { +ARB_ARBORIO_API arb::segment_tree load_swc_neuron_raw(const swc_data& data) { constexpr int soma_tag = 1; const auto n_samples = data.records().size(); @@ -230,7 +230,7 @@ ARB_ARBORIO_API arb::morphology load_swc_neuron(const swc_data& data) { if (auto it=std::find_if(R.begin(), R.end(), [](const auto& r) {return r.tag<1 || r.tag>4;}); it!=R.end()) { throw swc_unsupported_tag(R[std::distance(R.begin(), it)].id); } - return load_swc_arbor(data); + return load_swc_arbor_raw(data); } // Make a copy of the records and canonicalise them. @@ -326,8 +326,11 @@ ARB_ARBORIO_API arb::morphology load_swc_neuron(const swc_data& data) { } } - return arb::morphology(tree); + return tree; } +ARB_ARBORIO_API arb::morphology load_swc_neuron(const swc_data& data) { return {load_swc_neuron_raw(data)}; } +ARB_ARBORIO_API arb::morphology load_swc_arbor(const swc_data& data) { return {load_swc_arbor_raw(data)}; } + } // namespace arborio diff --git a/doc/concepts/morphology.rst b/doc/concepts/morphology.rst index 3c2fb57cbc..bcc8b78f72 100644 --- a/doc/concepts/morphology.rst +++ b/doc/concepts/morphology.rst @@ -602,6 +602,36 @@ because it has two children: the dendrites attached to its distal end. :width: 900 :align: center +Editing morphologies +-------------------- + +While a reified morphology cannot be edited -- it is immutable by definition -- +the segment tree can be changed. If you need to make such modifications, first +consider whethnr they should be stored in a file as this is often easier for +tracking provenance and version history. + +For the remaining cases, Arbor offers a limited suite of tools. First, most +morphology loaders have a 'raw' variant (C++) or flag (Python), that loads a +segment tree instead of a morphology. Segment trees are similarly immutable, but +by traversing the existing tree and adding/pruning subtrees into a new tree +changes are possible. From these edited trees, new morphologies can be formed. + +Two common editing operations are provided + +- ``split_at(t, i) -> (l, r)`` Split a segment tree ``t`` into two subtrees + ``l`` and ``r`` at a given segment id ``i``. + - ``r`` is the subtree of ``t`` that has the segment ``i`` at its root + - ``l`` is ``t`` with the segment ``i`` and its children removed + + The id ``i`` must be a valid id in ``t`` or ``mnpos`` in which case ``l == t`` + and ``r = {}``. +- ``join_at(t, i, o) -> j`` Given two segment trees ``t`` and ``o``, attach ``o`` to + ``t`` such that the segment in ``t`` identified by ``i`` is the parent of + ``o``. The id ``i`` must be valid in ``t`` and not ``mnpos`` (else ``t`` would + have two root segments). + +Note that ``join_at`` and ``split_at`` are inverse to each other. + .. _morph-formats: Supported file formats diff --git a/doc/cpp/morphology.rst b/doc/cpp/morphology.rst index f2ea79e560..9a21ea2481 100644 --- a/doc/cpp/morphology.rst +++ b/doc/cpp/morphology.rst @@ -4,25 +4,102 @@ Cable cell morphology ===================== Cell morphologies are required to describe a :ref:`cppcablecell`. -Morphologies can be constructed directly, or read from a number of +Morphologies can be constructed from :cpp:type:`segment_trees`, or read from a number of file formats; see :ref:`cppcablecell-morphology-construction` for details. -Morphology API --------------- +Segment tree +------------ -.. todo:: +A ``segment_tree` is -- as the name implies -- a set of segments arranged in a +tree structure, ie each segment has exactly one parent and no child is the +parent of any of its ancestors. The tree starts at a *root* segment which has no +parent. Each segment comprises two points in 3d space together with the radii at +these points. The segment's endpoints are called proximal (at the parent's +distal end) and distal (farther from the root). - Describe morphology methods. +Segment trees are used to form morphologies which ignore the 3d information +encoded in the segments and just utilise the radii, length, and tree-structure. +Branches in the tree occur where a segment has more than one child. The tree is +constructed by *appending* segments to the tree. Segments are numbered starting +at ``0`` in the order that they are added, with the first segment getting id +``0``, the second segment id ``1``, and so forth. A segment can not be added +before its parent, hence the first segment is always at the root. In this +manner, a segment tree is always guaranteed to be in a correct state, with +consistent parent-child indexing, and with ``n`` segments numbered from ``0`` to +``n-1``. The first parent must be :data:`mnpos`, indicating 'no parent'. -.. _cppcablecell-morphology-construction: -Constructing cell morphologies ------------------------------- +.. cpp:class:: segment_tree + + + .. cpp:function:: segment_tree() + + Construct an empty segment tree. + + .. cpp:function:: msize_t append(msize_t parent, const mpoint& prox, const mpoint& dist, int tag) + + Append a segment to the tree. Returns the new parent's id. + + .. cpp:function:: msize_t append(msize_t parent, const mpoint& dist, int tag) + + Append a segment to the tree whose proximal end has the location and + radius of the distal end of the parent segment. Returns the new + parent's id. + + This version of append can't be used for a segment at the root of the + tree, that is, when ``parent`` is :data:`mnpos`, in which case both + proximal and distal ends of the segment must be specified. + + .. cpp:function:: bool empty() + + If the tree is empty (i.e. whether it has size 0) + + .. cpp:function:: msize_t size() + + The number of segments. + + .. cpp:function:: std::vector parents() + + A list of parent indices of the segments. + + .. cpp:function:: std::vector segments() + + A list of the segments. + +.. cpp:function:: std::pair split_at(const segment_tree& t, msize_t id) + + Split a segment_tree into a pair of subtrees at the given id, + such that one tree is the subtree rooted at id and the other is the + original tree without said subtree. + +.. cpp:function:: segment_tree join_at(const segment_tree& t, msize_t id, const segment_tree& o) + + Join two subtrees at a given id, such that said id becomes the parent + of the inserted sub-tree. + +.. cpp:function:: bool equivalent(const segment_tree& l, const segment_tree& r) + + Two trees are equivalent if + 1. the root segments' ``prox`` and ``dist`` points and their ``tags`` are identical. + 2. recursively: all sub-trees starting at the current segment are pairwise equivalent. + + Note that under 1 we do not consider the ``id`` field. + +.. cpp:function:: segment_tree apply(const segment_tree& t, const isometry& i) + + Apply an :cpp:type:`isometry` to the segment tree, returns the transformed tree as a copy. + Isometries are rotations around an arbritary axis and/or translations; they can + be instantiated using ``isometry::translate`` and ``isometry::rotate`` and combined + using the ``*`` operator. + +Morphology API +-------------- .. todo:: - Description of segment trees. + Describe morphology methods. +.. _cppcablecell-morphology-construction: The stitch-builder interface ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/doc/python/morphology.rst b/doc/python/morphology.rst index caa3f3cced..c6b78f64fd 100644 --- a/doc/python/morphology.rst +++ b/doc/python/morphology.rst @@ -242,6 +242,36 @@ Cable cell morphology :param float radius: distal radius (μm) :param int tag: tag meta data of segment + .. method:: split_at(id) + + Split a segment_tree ``T`` into a pair of subtrees ``(L, R)`` such that + ``R`` is the subtree of ``T`` that starts at the given id and L is ``T`` + without ``R``. Splitting above the root ``mnpos`` returns ``(T, {})``. + + .. method:: join_at(id, other) + + Join two subtrees ``L`` and ``R`` at a given ``id`` in ``L``, such that + ``join_at`` is inverse to ``split_at`` for a proper choice of ``id``. + The join point ``id`` must be in ``L``. + + .. method:: apply_isometry(iso) + + Apply an :type:`isometry` to the segment tree, returns the transformed tree as a copy. + Isometries are rotations around an arbritary axis and/or translations; they can + be instantiated using ``translate`` and ``rotate`` and combined + using the ``*`` operator. + + :return: new tree + :param iso: isometry + + .. method:: equivalent(other) + + Two trees are equivalent if + 1. the root segments' ``prox`` and ``dist`` points and their ``tags`` + are identical. + 2. recursively: all sub-trees starting at the current segment are + equivalent. + .. attribute:: empty :type: bool @@ -558,7 +588,7 @@ region. SWC --- -.. py:function:: load_swc_arbor(filename) +.. py:function:: load_swc_arbor(filename, raw=False) Loads the :class:`morphology` from an SWC file according to arbor's SWC specifications. (See the morphology concepts :ref:`page ` for more details). @@ -586,16 +616,18 @@ SWC :param str filename: the name of the SWC file. - :rtype: morphology + :param bool raw: return segment_tree instead of morphology? + :rtype: morphology or segment_tree -.. py:function:: load_swc_neuron(filename) +.. py:function:: load_swc_neuron(filename, raw=False) Loads the :class:`morphology` from an SWC file according to NEURON's ``Import3D`` interpretation of the SWC specification. See :ref:`the SWC file documention ` for more details. :param str filename: the name of the SWC file. - :rtype: morphology + :param bool raw: return segment_tree instead of morphology? + :rtype: morphology or segment_tree .. _pyneuroml: @@ -707,6 +739,10 @@ Neurolucida The cable cell morphology. + .. py:attribute:: segment_tree + + The raw segment tree. + .. py:attribute:: labels The labeled regions and locations extracted from the meta data. The four canonical regions are labeled diff --git a/python/morphology.cpp b/python/morphology.cpp index ba238754a0..52bf43bba5 100644 --- a/python/morphology.cpp +++ b/python/morphology.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include @@ -245,16 +246,33 @@ void register_morphology(py::module& m) { "A list with the parent index of each segment.") .def_property_readonly("segments", [](const arb::segment_tree& st){return st.segments();}, "A list of the segments.") + .def("apply_isometry", + [](const arb::segment_tree& t, const arb::isometry& i) { return arb::apply(t, i); }, + "Apply an isometry to all segments in the tree.") + .def("split_at", + [](const arb::segment_tree& t, arb::msize_t id) { return arb::split_at(t, id); }, + "Split into a pair of trees at the given id, such that one tree is the subtree rooted at id and the other is the original tree without said subtree.") + .def("join_at", + [](const arb::segment_tree& t, arb::msize_t id, const arb::segment_tree& o) { return arb::join_at(t, id, o); }, + "Join two subtrees at a given id, such that said id becomes the parent of the inserted sub-tree.") + .def("equivalent", + [](const arb::segment_tree& t, const arb::segment_tree& o) { return arb::equivalent(t, o); }, + "Two trees are equivalent, but not neccessarily identical, ie they have the same segments and structure.") .def("__str__", [](const arb::segment_tree& s) { return util::pprintf("", s);}); - // Function that creates a morphology from an swc file. + using morph_or_tree = std::variant; + + // Function that creates a morphology/segment_tree from an swc file. // Wraps calls to C++ functions arborio::parse_swc() and arborio::load_swc_arbor(). m.def("load_swc_arbor", - [](py::object fn) { + [](py::object fn, bool raw) -> morph_or_tree { try { auto contents = util::read_file_or_buffer(fn); auto data = arborio::parse_swc(contents); + if (raw) { + return arborio::load_swc_arbor_raw(data); + } return arborio::load_swc_arbor(data); } catch (arborio::swc_error& e) { @@ -263,28 +281,31 @@ void register_morphology(py::module& m) { } }, "filename_or_stream"_a, - "Generate a morphology from an SWC file following the rules prescribed by Arbor.\n" + pybind11::arg_v("raw", false, "Return a segment tree instead of a fully formed morphology"), + "Generate a morphology/segment_tree from an SWC file following the rules prescribed by Arbor.\n" "Specifically:\n" - "* Single-segment somas are disallowed.\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 morphology."); - + " * Single-segment somas are disallowed.\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 morphology."); m.def("load_swc_neuron", - [](py::object fn) { + [](py::object fn, bool raw) -> morph_or_tree { try { auto contents = util::read_file_or_buffer(fn); auto data = arborio::parse_swc(contents); + if (raw) { + return arborio::load_swc_neuron_raw(data); + } return arborio::load_swc_neuron(data); } catch (arborio::swc_error& e) { // Try to produce helpful error messages for SWC parsing errors. - throw pyarb_error( - util::pprintf("NEURON SWC: parse error: {}", e.what())); + throw pyarb_error(util::pprintf("NEURON SWC: parse error: {}", e.what())); } }, "filename_or_stream"_a, + pybind11::arg_v("raw", false, "Return a segment tree instead of a fully formed morphology"), "Generate a morphology from an SWC file following the rules prescribed by NEURON.\n" "See the documentation https://docs.arbor-sim.org/en/latest/fileformat/swc.html\n" "for a detailed description of the interpretation."); @@ -327,7 +348,10 @@ void register_morphology(py::module& m) { asc_morphology .def_readonly("morphology", &arborio::asc_morphology::morphology, - "The cable cell morphology") + "The cable cell morphology.") + .def_readonly("segment_tree", + &arborio::asc_morphology::segment_tree, + "The raw segment tree.") .def_property_readonly("labels", [](const arborio::asc_morphology& m) {return label_dict_proxy(m.labels);}, "The four canonical regions are labeled 'soma', 'axon', 'dend' and 'apic'."); diff --git a/test/unit/test_segment_tree.cpp b/test/unit/test_segment_tree.cpp index 4f5baff1d7..56909796ee 100644 --- a/test/unit/test_segment_tree.cpp +++ b/test/unit/test_segment_tree.cpp @@ -111,3 +111,226 @@ TEST(segment_tree, fuzz) { } } +::testing::AssertionResult trees_equivalent(const arb::segment_tree& a, + const arb::segment_tree& b) { + if (!arb::equivalent(a, b)) return ::testing::AssertionFailure() << "Trees are not equivalent:\n" + << a + << "\nand:\n" + << b; + return ::testing::AssertionSuccess(); +} + +TEST(segment_tree, split) { + // linear chain + { + arb::segment_tree tree; + tree.append(arb::mnpos, {0, 0, 0}, {0, 0, 1}, 0); + tree.append(0, {0, 0, 1}, {0, 0, 2}, 0); + tree.append(1, {0, 0, 2}, {0, 0, 3}, 0); + tree.append(2, {0, 0, 3}, {0, 0, 4}, 0); + tree.append(3, {0, 0, 4}, {0, 0, 5}, 0); + { + arb::segment_tree p, q; + q.append(arb::mnpos, {0, 0, 0}, {0, 0, 1}, 0); + q.append(0, {0, 0, 1}, {0, 0, 2}, 0); + q.append(1, {0, 0, 2}, {0, 0, 3}, 0); + q.append(2, {0, 0, 3}, {0, 0, 4}, 0); + q.append(3, {0, 0, 4}, {0, 0, 5}, 0); + const auto& [l, r] = arb::split_at(tree, 0); + EXPECT_TRUE(trees_equivalent(p, l)); + EXPECT_TRUE(trees_equivalent(q, r)); + } + { + arb::segment_tree p, q; + p.append(arb::mnpos, {0, 0, 0}, {0, 0, 1}, 0); + p.append(0, {0, 0, 1}, {0, 0, 2}, 0); + q.append(arb::mnpos, {0, 0, 2}, {0, 0, 3}, 0); + q.append(0, {0, 0, 3}, {0, 0, 4}, 0); + q.append(1, {0, 0, 4}, {0, 0, 5}, 0); + const auto& [l, r] = arb::split_at(tree, 2); + EXPECT_TRUE(trees_equivalent(p, l)); + EXPECT_TRUE(trees_equivalent(q, r)); + } + { + arb::segment_tree p, q; + p.append(arb::mnpos, {0, 0, 0}, {0, 0, 1}, 0); + p.append(0, {0, 0, 1}, {0, 0, 2}, 0); + p.append(1, {0, 0, 2}, {0, 0, 3}, 0); + p.append(2, {0, 0, 3}, {0, 0, 4}, 0); + q.append(arb::mnpos, {0, 0, 4}, {0, 0, 5}, 0); + const auto& [l, r] = arb::split_at(tree, 4); + EXPECT_TRUE(trees_equivalent(p, l)); + EXPECT_TRUE(trees_equivalent(q, r)); + } + } + // Error cases + { + arb::segment_tree t; + EXPECT_THROW(arb::split_at(t, arb::mnpos), arb::arbor_exception); + EXPECT_THROW(arb::split_at(t, 1), arb::arbor_exception); + } + // gnarly tree + // (npos) - 0 - 1 - 4 + // \ + // 2 - 3 + // \ + // 5 + // NB: Splitting _will_ re-order segments. So we have to be careful when + // building values to compare against. + { + arb::segment_tree tree; + tree.append(arb::mnpos, {0, 0, 0}, {0, 0, 1}, 0); // 0 + tree.append(0, {0, 0, 1}, {0, 0, 2}, 0); // 1 + tree.append(0, {0, 0, 2}, {0, 0, 3}, 0); // 2 + tree.append(2, {0, 0, 3}, {0, 0, 4}, 0); // 3 + tree.append(1, {0, 0, 4}, {0, 0, 5}, 0); // 4 + tree.append(2, {0, 0, 5}, {0, 0, 6}, 0); // 5 + { + arb::segment_tree p, q; + + q.append(arb::mnpos, {0, 0, 0}, {0, 0, 1}, 0); + q.append(0, {0, 0, 2}, {0, 0, 3}, 0); + q.append(1, {0, 0, 5}, {0, 0, 6}, 0); + q.append(1, {0, 0, 3}, {0, 0, 4}, 0); + q.append(0, {0, 0, 1}, {0, 0, 2}, 0); + q.append(4, {0, 0, 4}, {0, 0, 5}, 0); + + const auto& [l, r] = arb::split_at(tree, 0); + + EXPECT_TRUE(trees_equivalent(p, l)); + EXPECT_TRUE(trees_equivalent(q, r)); + } + { + arb::segment_tree p, q; + + p.append(arb::mnpos, {0, 0, 0}, {0, 0, 1}, 0); + p.append(0, {0, 0, 2}, {0, 0, 3}, 0); + p.append(1, {0, 0, 5}, {0, 0, 6}, 0); + p.append(1, {0, 0, 3}, {0, 0, 4}, 0); + + q.append(arb::mnpos, {0, 0, 1}, {0, 0, 2}, 0); + q.append(0, {0, 0, 4}, {0, 0, 5}, 0); + + const auto& [l, r] = arb::split_at(tree, 1); + EXPECT_TRUE(trees_equivalent(p, l)); + EXPECT_TRUE(trees_equivalent(q, r)); + } + { + arb::segment_tree p, q; + + p.append(arb::mnpos, {0, 0, 0}, {0, 0, 1}, 0); + p.append(0, {0, 0, 1}, {0, 0, 2}, 0); + p.append(1, {0, 0, 4}, {0, 0, 5}, 0); + + q.append(arb::mnpos, {0, 0, 2}, {0, 0, 3}, 0); + q.append(0, {0, 0, 5}, {0, 0, 6}, 0); + q.append(0, {0, 0, 3}, {0, 0, 4}, 0); + + const auto& [l, r] = arb::split_at(tree, 2); + EXPECT_TRUE(trees_equivalent(p, l)); + EXPECT_TRUE(trees_equivalent(q, r)); + } + } +} + +TEST(segment_tree, join) { + // linear chain + { + arb::segment_tree tree; + tree.append(arb::mnpos, {0, 0, 0}, {0, 0, 1}, 0); + tree.append(0, {0, 0, 1}, {0, 0, 2}, 0); + tree.append(1, {0, 0, 2}, {0, 0, 3}, 0); + tree.append(2, {0, 0, 3}, {0, 0, 4}, 0); + tree.append(3, {0, 0, 4}, {0, 0, 5}, 0); + { + arb::segment_tree p, q; + + q.append(arb::mnpos, {0, 0, 0}, {0, 0, 1}, 0); + q.append(0, {0, 0, 1}, {0, 0, 2}, 0); + q.append(1, {0, 0, 2}, {0, 0, 3}, 0); + q.append(2, {0, 0, 3}, {0, 0, 4}, 0); + q.append(3, {0, 0, 4}, {0, 0, 5}, 0); + + const auto& t = arb::join_at(p, arb::mnpos, q); + EXPECT_TRUE(trees_equivalent(tree, t)); + } + { + arb::segment_tree p, q; + + p.append(arb::mnpos, {0, 0, 0}, {0, 0, 1}, 0); + p.append(0, {0, 0, 1}, {0, 0, 2}, 0); + + q.append(arb::mnpos, {0, 0, 2}, {0, 0, 3}, 0); + q.append(0, {0, 0, 3}, {0, 0, 4}, 0); + q.append(1, {0, 0, 4}, {0, 0, 5}, 0); + + const auto& t = arb::join_at(p, 1, q); + EXPECT_TRUE(trees_equivalent(tree, t)); + } + } + + // Error cases + { + arb::segment_tree t; + EXPECT_THROW(arb::split_at(t, arb::mnpos), arb::arbor_exception); + EXPECT_THROW(arb::split_at(t, 1), arb::arbor_exception); + } + // gnarly tree + // (npos) - 0 - 1 - 4 + // \ + // 2 - 3 + // \ + // 5 + // NB: Joining _will_ re-order segments. So we have to be careful when + // building values to compare against. + { + arb::segment_tree tree; + tree.append(arb::mnpos, {0, 0, 0}, {0, 0, 1}, 0); // 0 + tree.append(0, {0, 0, 1}, {0, 0, 2}, 0); // 1 + tree.append(0, {0, 0, 2}, {0, 0, 3}, 0); // 2 + tree.append(2, {0, 0, 3}, {0, 0, 4}, 0); // 3 + tree.append(1, {0, 0, 4}, {0, 0, 5}, 0); // 4 + tree.append(2, {0, 0, 5}, {0, 0, 6}, 0); // 5 + { + arb::segment_tree p, q; + + q.append(arb::mnpos, {0, 0, 0}, {0, 0, 1}, 0); + q.append(0, {0, 0, 2}, {0, 0, 3}, 0); + q.append(1, {0, 0, 5}, {0, 0, 6}, 0); + q.append(1, {0, 0, 3}, {0, 0, 4}, 0); + q.append(0, {0, 0, 1}, {0, 0, 2}, 0); + q.append(4, {0, 0, 4}, {0, 0, 5}, 0); + + const auto& t = arb::join_at(p, arb::mnpos, q); + EXPECT_TRUE(trees_equivalent(tree, t)); + } + { + arb::segment_tree p, q; + + p.append(arb::mnpos, {0, 0, 0}, {0, 0, 1}, 0); + p.append(0, {0, 0, 2}, {0, 0, 3}, 0); + p.append(1, {0, 0, 5}, {0, 0, 6}, 0); + p.append(1, {0, 0, 3}, {0, 0, 4}, 0); + + q.append(arb::mnpos, {0, 0, 1}, {0, 0, 2}, 0); + q.append(0, {0, 0, 4}, {0, 0, 5}, 0); + + const auto& t = arb::join_at(p, 0, q); + EXPECT_TRUE(trees_equivalent(tree, t)); + } + { + arb::segment_tree p, q; + + p.append(arb::mnpos, {0, 0, 0}, {0, 0, 1}, 0); + p.append(0, {0, 0, 1}, {0, 0, 2}, 0); + p.append(1, {0, 0, 4}, {0, 0, 5}, 0); + + q.append(arb::mnpos, {0, 0, 2}, {0, 0, 3}, 0); + q.append(0, {0, 0, 5}, {0, 0, 6}, 0); + q.append(0, {0, 0, 3}, {0, 0, 4}, 0); + + const auto& t = arb::join_at(p, 0, q); + EXPECT_TRUE(trees_equivalent(tree, t)); + } + } +}