diff --git a/arbor/include/arbor/morph/isometry.hpp b/arbor/include/arbor/morph/isometry.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cb06ee28dd58145a51fbc3fba2f8c4fa122da05f
--- /dev/null
+++ b/arbor/include/arbor/morph/isometry.hpp
@@ -0,0 +1,62 @@
+#pragma once
+
+#include <arbor/export.hpp>
+#include <arbor/math.hpp>
+
+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 <typename PointLike>
+ 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 <typename PointLike>
+ 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 <typename PointLike>
+ 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 09cc179dc72537d0af1b77900595f767e45a84f4..6088a8ec072e59d864f2418e7548de40bd9dca9c 100644
--- a/arbor/include/arbor/morph/place_pwlin.hpp
+++ b/arbor/include/arbor/morph/place_pwlin.hpp
@@ -10,64 +10,10 @@
#include <arbor/export.hpp>
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/primitives.hpp>
-#include <arbor/math.hpp>
+#include <arbor/morph/isometry.hpp>
-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 <typename PointLike>
- 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 <typename PointLike>
- 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 <typename PointLike>
- 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 cdb6472acb8f3e61e638e6e22de8ee7bd0560704..4be04670b1785f5e09d75a45bacbb3e4bd84d842 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 e934f102eb1d709e3fbb3c20136e0658e55fba29..7d59be59aead032a89f90235ee7110e18a8fa145 100644
--- a/arbor/include/arbor/morph/segment_tree.hpp
+++ b/arbor/include/arbor/morph/segment_tree.hpp
@@ -7,6 +7,7 @@
#include <arbor/export.hpp>
#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/isometry.hpp>
namespace arb {
@@ -52,8 +53,36 @@ public:
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 <L, R> 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<segment_tree, segment_tree>
+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 70498e47095995c6cea8166758fc5b12a4488c13..3a10a015b516662c907f014246115d4e4013c8e5 100644
--- a/arbor/morph/segment_tree.cpp
+++ b/arbor/morph/segment_tree.cpp
@@ -1,6 +1,6 @@
#include <iostream>
#include <stdexcept>
-#include <unordered_map>
+#include <map>
#include <vector>
#include <arbor/morph/morphexcept.hpp>
@@ -12,6 +12,119 @@
namespace arb {
+struct node_t {
+ msize_t parent;
+ msize_t id;
+};
+
+using node_p = std::function<bool(const node_t&)>;
+
+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<msize_t, std::vector<msize_t>> tree_to_children(const segment_tree& tree) {
+ const auto& parents = tree.parents();
+ std::map<msize_t, std::vector<msize_t>> 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<node_t>{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<segment_tree, segment_tree>
+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<arb::msegment> 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<std::pair<arb::msize_t, arb::msize_t>> 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 cf82961380c093048e314d67513c9899138649eb..3fd4d1fc8f83e694c46747d46ea7fb88f7ead720 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 27c2663b2e95bccc6049f3199c9ec5ce45ed5896..3646548e6f773b106a3414994434041f406b4b90 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 <morphology> element.
+ // Morphology constructed from a single NeuroML <morphology> 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 26189571ba91a9c7cb21e81ce44fca2a80c339f3..38b59e0c5ae26fbb90595c744d2ee19de41e5f44 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 af033796b0c12d75907172928a69967b96829188..832a7d2451f3b3a2a55bc8966c2a4c265076c1a9 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 9fc2f2613d5bf8a268a8e6ffd390f0b0ca10d1e1..34a0fec657b90a66aa410b3340f9e41388aad338 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 3c2fb57cbc677a8af67db96b100c0c09c0573594..bcc8b78f72b86ce2e8226240ea3201f3c48ce070 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 f2ea79e56042342d31e1d96c04049747f68f3951..9a21ea2481a0c3c6c19904744db9393747928921 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<msize_t> parents()
+
+ A list of parent indices of the segments.
+
+ .. cpp:function:: std::vector<msegment> segments()
+
+ A list of the segments.
+
+.. cpp:function:: std::pair<segment_tree, segment_tree> 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 caa3f3ccedbe32c775e0d0d0f88e127fb2fed54a..c6b78f64fde94f2cc9ec84868e795d09d45a94d9 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 <morph-formats>` 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 <formatswc-neuron>` 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 ba238754a000fc6227f73ecb67eae1b502f99d05..52bf43bba503f55c650531088df0a9c5ef536fd9 100644
--- a/python/morphology.cpp
+++ b/python/morphology.cpp
@@ -1,5 +1,6 @@
#include <fstream>
#include <tuple>
+#include <variant>
#include <pybind11/operators.h>
#include <pybind11/pybind11.h>
@@ -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("<arbor.segment_tree:\n{}>", s);});
- // Function that creates a morphology from an swc file.
+ using morph_or_tree = std::variant<arb::segment_tree, arb::morphology>;
+
+ // 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 4f5baff1d7aba5163dabbd0f2fd293b138e78098..56909796eea6eb1729c07480e0ccd602fa288f1b 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));
+ }
+ }
+}