Allow editing morphologies

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

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;
 

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 {

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 @@ 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 <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

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);

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;
 

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.

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

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) {