Add place_pwlin for morphology geometry queries. (#1062)

* Add new `isometry` class representing a direct Euclidean isometry comprising a rotation and translation. Translations are always taken with respect to absolute coordinates, but rotations are composed with respect to intrinsic coordinates.
* Add new `place_pwlin` class, that takes a morphology and an isometry, and answers location queries: where in space does a particular `mlocation` lie?
* Split out (some of) the common code between `place_pwlin` and `embed_pwlin`.
* Add equality operator to `mpoint`.

This functionality is added to aid the presentation of the local field potential example.

arbor/CMakeLists.txt

@@ -36,6 +36,7 @@ set(arbor_sources
     morph/morphexcept.cpp
     morph/morphology.cpp
     morph/mprovider.cpp
+    morph/place_pwlin.cpp
     morph/primitives.cpp
     morph/region.cpp
     morph/sample_tree.cpp

arbor/include/arbor/morph/place_pwlin.hpp

+#pragma once
+
+// 'Place' morphology in 3-d by applying an isometry to
+// sample points and interpolating linearly.
+
+#include <cmath>
+#include <utility>
+
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arbor/math.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);
+    }
+};
+
+struct place_pwlin_data;
+
+struct place_pwlin {
+    explicit place_pwlin(const morphology& m, const isometry& iso = isometry{});
+    mpoint at(mlocation loc) const;
+
+private:
+    std::shared_ptr<place_pwlin_data> data_;
+};
+
+} // namespace arb
+
+

arbor/include/arbor/morph/primitives.hpp

@@ -23,6 +23,12 @@ struct mpoint {
     double radius;   // [μm]
 
     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);
+    }
 };
 
 mpoint lerp(const mpoint& a, const mpoint& b, double u);

arbor/morph/embed_pwlin.cpp

@@ -6,6 +6,7 @@
 #include <arbor/morph/morphology.hpp>
 #include <arbor/morph/primitives.hpp>
 
+#include "morph/pwlin_common.hpp"
 #include "util/piecewise.hpp"
 #include "util/range.hpp"
 #include "util/rangeutil.hpp"
@@ -16,27 +17,6 @@ namespace arb {
 
 using util::rat_element;
 
-template <unsigned p, unsigned q>
-using pw_ratpoly = util::pw_elements<rat_element<p, q>>;
-
-template <unsigned p, unsigned q>
-using branch_pw_ratpoly = std::vector<pw_ratpoly<p, q>>;
-
-template <unsigned p, unsigned q>
-double interpolate(const branch_pw_ratpoly<p, q>& f, unsigned bid, double pos) {
-    const auto& pw = f.at(bid);
-    unsigned index = pw.index_of(pos);
-
-    const auto& element = pw.element(index);
-    std::pair<double, double> bounds = pw.interval(index);
-
-    if (bounds.first==bounds.second) return element[0];
-    else {
-        double x = (pos-bounds.first)/(bounds.second-bounds.first);
-        return element(x);
-    }
-}
-
 // Length, area, and ixa are polynomial or rational polynomial functions of branch position,
 // continuos and monotonically increasing with respect to distance from root.
 //

arbor/morph/place_pwlin.cpp

+#include <cmath>
+#include <memory>
+#include <vector>
+
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/place_pwlin.hpp>
+#include <arbor/morph/primitives.hpp>
+
+#include "morph/pwlin_common.hpp"
+#include "util/piecewise.hpp"
+#include "util/ratelem.hpp"
+#include "util/span.hpp"
+
+namespace arb {
+
+using util::rat_element;
+
+struct place_pwlin_data {
+    branch_pw_ratpoly<1, 0> x, y, z, r; // [µm]
+
+    explicit place_pwlin_data(msize_t n_branch):
+        x(n_branch), y(n_branch), z(n_branch), r(n_branch)
+    {}
+};
+
+mpoint place_pwlin::at(mlocation loc) const {
+    return { interpolate(data_->x, loc.branch, loc.pos),
+             interpolate(data_->y, loc.branch, loc.pos),
+             interpolate(data_->z, loc.branch, loc.pos),
+             interpolate(data_->r, loc.branch, loc.pos) };
+}
+
+place_pwlin::place_pwlin(const arb::morphology& m, const isometry& iso) {
+    msize_t n_branch = m.num_branches();
+    data_ = std::make_shared<place_pwlin_data>(n_branch);
+
+    if (!n_branch) return;
+
+    const auto& samples = m.samples();
+
+    std::vector<double> sample_pos_on_branch;
+
+    for (msize_t bid = 0; bid<n_branch; ++bid) {
+        auto sample_indices = util::make_range(m.branch_indexes(bid));
+        if (bid==0 && m.spherical_root()) {
+            arb_assert(sample_indices.size()==1);
+            arb_assert(sample_indices.front()==0);
+
+            // Use the next distinct sample, if it exists, to determine the
+            // proximal-distal axis for a spherical root.
+
+            mpoint c = samples[0].loc;
+            mpoint d = c;
+            for (msize_t i = 1; i<samples.size(); ++i) {
+                const auto p = samples[i].loc;
+                if (p.x!=c.x || p.y!=c.y || p.z!=c.z) {
+                    d = p;
+                    break;
+                }
+            }
+
+            mpoint u0 = c, u1 = c;
+
+            if (c.x!=d.x || c.y!=d.y || c.z!=d.z) {
+                double dx = d.x-c.x;
+                double dy = d.y-c.y;
+                double dz = d.z-c.z;
+                double scale = c.radius/std::sqrt(dx*dx+dy*dy+dz*dz);
+
+                u0.x -= dx*scale;
+                u0.y -= dy*scale;
+                u0.z -= dz*scale;
+
+                u1.x += dx*scale;
+                u1.y += dy*scale;
+                u1.z += dz*scale;
+            }
+            else {
+                // No luck, so pick distal as negative z-axis.
+                u0.z += c.radius;
+                u1.z -= c.radius;
+            }
+
+            u0 = iso.apply(u0);
+            u1 = iso.apply(u1);
+
+            data_->x[bid].push_back(0., 1., rat_element<1, 0>(u0.x, u1.x));
+            data_->y[bid].push_back(0., 1., rat_element<1, 0>(u0.y, u1.y));
+            data_->z[bid].push_back(0., 1., rat_element<1, 0>(u0.z, u1.z));
+            data_->r[bid].push_back(0., 1., rat_element<1, 0>(u0.radius, u1.radius));
+        }
+        else {
+            arb_assert(sample_indices.size()>1);
+
+            sample_pos_on_branch.reserve(samples.size());
+            sample_pos_on_branch = {0};
+
+            for (auto i: util::count_along(sample_indices)) {
+                if (!i) continue;
+
+                sample_pos_on_branch.push_back(
+                    sample_pos_on_branch.back()+
+                    distance(samples[sample_indices[i-1]], samples[sample_indices[i]]));
+            }
+
+            double branch_length = sample_pos_on_branch.back();
+            double length_scale = branch_length>0? 1./branch_length: 0;
+
+            for (auto& x: sample_pos_on_branch) {
+                x *= length_scale;
+            }
+
+            if (length_scale==0) {
+                // Zero-length branch case?
+
+                mpoint p = iso.apply(samples[sample_indices[0]].loc);
+
+                data_->x[bid].push_back(0., 1., rat_element<1, 0>(p.x, p.x));
+                data_->y[bid].push_back(0., 1., rat_element<1, 0>(p.y, p.y));
+                data_->z[bid].push_back(0., 1., rat_element<1, 0>(p.z, p.z));
+                data_->r[bid].push_back(0., 1., rat_element<1, 0>(p.radius, p.radius));
+            }
+            else {
+                for (auto i: util::count_along(sample_indices)) {
+                    if (!i) continue;
+
+                    double p0 = i>1? sample_pos_on_branch[i-1]: 0;
+                    double p1 = sample_pos_on_branch[i];
+                    if (p0==p1) continue;
+
+                    mpoint u0 = iso.apply(samples[sample_indices[i-1]].loc);
+                    mpoint u1 = iso.apply(samples[sample_indices[i]].loc);
+
+                    data_->x[bid].push_back(p0, p1, rat_element<1, 0>(u0.x, u1.x));
+                    data_->y[bid].push_back(p0, p1, rat_element<1, 0>(u0.y, u1.y));
+                    data_->z[bid].push_back(p0, p1, rat_element<1, 0>(u0.z, u1.z));
+                    data_->r[bid].push_back(p0, p1, rat_element<1, 0>(u0.radius, u1.radius));
+                }
+            }
+        }
+    }
+};
+
+} // namespace arb

arbor/morph/pwlin_common.hpp

+#pragma once
+
+// Per-branch piecewise rational polynomial definitions and interpolation.
+
+#include <utility>
+#include <vector>
+
+#include "util/piecewise.hpp"
+#include "util/ratelem.hpp"
+
+namespace arb {
+
+template <unsigned p, unsigned q>
+using pw_ratpoly = util::pw_elements<util::rat_element<p, q>>;
+
+template <unsigned p, unsigned q>
+using branch_pw_ratpoly = std::vector<pw_ratpoly<p, q>>;
+
+template <unsigned p, unsigned q>
+double interpolate(const pw_ratpoly<p, q>& f, double pos) {
+    unsigned index = f.index_of(pos);
+
+    const auto& element = f.element(index);
+    std::pair<double, double> bounds = f.interval(index);
+
+    if (bounds.first==bounds.second) return element[0];
+    else {
+        double x = (pos-bounds.first)/(bounds.second-bounds.first);
+        return element(x);
+    }
+}
+
+template <unsigned p, unsigned q>
+double interpolate(const branch_pw_ratpoly<p, q>& f, unsigned bid, double pos) {
+    return interpolate(f.at(bid), pos);
+}
+
+} // namespace arb

test/unit/CMakeLists.txt

@@ -119,6 +119,7 @@ set(unit_sources
     test_morphology.cpp
     test_morph_embedding.cpp
     test_morph_expr.cpp
+    test_morph_place.cpp
     test_morph_primitives.cpp
     test_multi_event_stream.cpp
     test_optional.cpp

test/unit/test_morph_place.cpp

+#include <cmath>
+#include <vector>
+
+#include <arbor/math.hpp>
+#include <arbor/morph/place_pwlin.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/sample_tree.hpp>
+
+#include "util/piecewise.hpp"
+
+#include "../test/gtest.h"
+#include "common.hpp"
+
+using namespace arb;
+
+namespace {
+
+struct v3 {
+    double x, y, z;
+    friend double dot(v3 p, v3 q) {
+        return p.x*q.x+p.y*q.y+p.z*q.z;
+    }
+    friend v3 cross(v3 p, v3 q) {
+        return {p.y*q.z-p.z*q.y, p.z*q.x-p.x*q.z, p.x*q.y-p.y*q.x};
+    }
+    friend v3 operator*(v3 p, double s) {
+        return {s*p.x, s*p.y, s*p.z};
+    }
+    friend v3 operator*(double s, v3 p) {
+        return p*s;
+    }
+    friend v3 operator/(v3 p, double s) {
+        return p*(1./s);
+    }
+    friend v3 operator+(v3 p, v3 q) {
+        return {p.x+q.x, p.y+q.y, p.z+q.z};
+    }
+    friend double length(v3 x) {
+        return std::sqrt(dot(x, x));
+    }
+};
+
+::testing::AssertionResult mpoint_almost_eq(mpoint a, mpoint b) {
+    using FP = testing::internal::FloatingPoint<double>;
+    if (FP(a.x).AlmostEquals(FP(b.x)) &&
+        FP(a.y).AlmostEquals(FP(b.y)) &&
+        FP(a.z).AlmostEquals(FP(b.z)) &&
+        FP(a.radius).AlmostEquals(FP(b.radius)))
+    {
+        return ::testing::AssertionSuccess();
+    }
+    else {
+        return ::testing::AssertionFailure() << "mpoint values "
+            << '(' << a.x << ',' << a.y << ',' << a.z << ';' << a.radius << ')'
+            << " and "
+            << '(' << b.x << ',' << b.y << ',' << b.z << ';' << b.radius << ')'
+            << " differ";
+    }
+}
+
+::testing::AssertionResult v3_almost_eq(v3 a, v3 b) {
+    using FP = testing::internal::FloatingPoint<double>;
+    if (FP(a.x).AlmostEquals(FP(b.x)) &&
+        FP(a.y).AlmostEquals(FP(b.y)) &&
+        FP(a.z).AlmostEquals(FP(b.z)))
+    {
+        return ::testing::AssertionSuccess();
+    }
+    else {
+        return ::testing::AssertionFailure() << "xyz values "
+            << '(' << a.x << ',' << a.y << ',' << a.z << ')'
+            << " and "
+            << '(' << b.x << ',' << b.y << ',' << b.z << ')'
+            << " differ";
+    }
+}
+
+} // anonymous namespace
+
+TEST(isometry, translate) {
+    mpoint p{1, 2, 3, 4};
+
+    {
+        isometry id;
+        mpoint q = id.apply(p);
+        EXPECT_EQ(p.x, q.x);
+        EXPECT_EQ(p.y, q.y);
+        EXPECT_EQ(p.z, q.z);
+        EXPECT_EQ(p.radius, q.radius);
+    }
+    {
+        isometry shift = isometry::translate(0.5, 1.5, 3.5);
+        mpoint q = shift.apply(p);
+
+        EXPECT_EQ(p.x+0.5, q.x);
+        EXPECT_EQ(p.y+1.5, q.y);
+        EXPECT_EQ(p.z+3.5, q.z);
+        EXPECT_EQ(p.radius, q.radius);
+
+        // Should work with anything with floating point x, y, z.
+        struct v3 {
+            double x, y, z;
+        };
+
+        isometry shift_bis = isometry::translate(v3{0.5, 1.5, 3.5});
+        mpoint q_bis = shift_bis.apply(p);
+        EXPECT_EQ(q, q_bis);
+    }
+}
+
+TEST(isometry, rotate) {
+    {
+        // Rotation about axis through p should do nothing.
+        mpoint p{1, 2, 3, 4};
+
+        isometry rot = isometry::rotate(1.23, p);
+        mpoint q = rot.apply(p);
+
+        EXPECT_TRUE(mpoint_almost_eq(p, q));
+    }
+    {
+        // Rotate about an arbitrary axis.
+        //
+        // Geometry says rotating x about (normalized) u by
+        // θ should give a vector u(u.x) + (u×x)×u cos θ + u×x sin 0
+
+        double theta = 0.234;
+        v3 axis{-1, 3, 2.2};
+
+        isometry rot = isometry::rotate(theta, axis);
+        v3 x = {1, 2, 3};
+        v3 q = rot.apply(x);
+
+        v3 u = axis/length(axis);
+        v3 p = u*dot(u, x) + cross(cross(u, x), u)*std::cos(theta) + cross(u, x)*std::sin(theta);
+
+        EXPECT_TRUE(v3_almost_eq(p, q));
+    }
+}
+
+TEST(isometry, compose) {
+    // For an isometry (r, t) (r the rotation, t the translation),
+    // we compose them by:
+    //   (r, t) ∘ (s, u) = (sr, t+u).
+    //
+    // On the other hand, for a vector x:
+    //   (r, t) ((s, u) x)
+    //      = (r, t) (sx + u)
+    //      = rsx + ru + t
+    //      = (rs, ru+t) x
+    //      = ((s, ru) ∘ (r, t)) x
+
+    double theta1 = -0.3;
+    v3 axis1{1, -3, 4};
+    v3 shift1{-2, -1, 3};
+
+    double theta2 = 0.6;
+    v3 axis2{-1, 2, 0};
+    v3 shift2{3, 0, -1};
+
+    v3 p{5., 6., 7.};
+
+    // Construct (r, t):
+    isometry r = isometry::rotate(theta1, axis1);
+    isometry r_t = r*isometry::translate(shift1);
+
+    // Construct (s, u):
+    isometry s = isometry::rotate(theta2, axis2);
+    isometry s_u = s*isometry::translate(shift2);
+
+    // Construct (s, ru):
+    isometry s_ru = s*isometry::translate(r.apply(shift2));
+
+    // Compare application of s_u and r_t against application
+    // of (s, ru) ∘ (r, t).
+
+    v3 u = r_t.apply(s_u.apply(p));
+    v3 v = (s_ru * r_t).apply(p);
+
+    EXPECT_TRUE(v3_almost_eq(u, v));
+}
+
+TEST(place_pwlin, cable) {
+    using pvec = std::vector<msize_t>;
+    using svec = std::vector<msample>;
+
+    // L-shaped simple cable.
+    // 0.25 of the length in the z-direction,
+    // 0.75 of the length in the x-direction.
+
+    pvec parents = {mnpos, 0, 1};
+    svec samples = {
+        {{  0,  0,  0,  2}, 1},
+        {{  0,  0,  1,  2}, 1},
+        {{  3,  0,  1,  2}, 1}
+    };
+
+    sample_tree sm(samples, parents);
+    morphology m(sm, false);
+
+    {
+        // With no transformation:
+        place_pwlin pl(m);
+
+        // Interpolated points.
+        mpoint p_0 = pl.at(mlocation{0, 0.});
+        mpoint p_1 = pl.at(mlocation{0, 0.125});
+        mpoint p_2 = pl.at(mlocation{0, 0.25});
+        mpoint p_3 = pl.at(mlocation{0, 0.5});
+
+        // Expected results.
+        mpoint x_0{0.0, 0.0, 0.0, 2.};
+        mpoint x_1{0.0, 0.0, 0.5, 2.};
+        mpoint x_2{0.0, 0.0, 1.0, 2.};
+        mpoint x_3{1.0, 0.0, 1.0, 2.};
+
+        EXPECT_TRUE(mpoint_almost_eq(x_0, p_0));
+        EXPECT_TRUE(mpoint_almost_eq(x_1, p_1));
+        EXPECT_TRUE(mpoint_almost_eq(x_2, p_2));
+        EXPECT_TRUE(mpoint_almost_eq(x_3, p_3));
+    }
+    {
+        // With a rotation and translation:
+
+        double theta = -0.3;
+        v3 axis{1, -3, 4};
+        v3 shift{-2, -1, 3};
+
+        isometry iso = isometry::rotate(theta, axis)*isometry::translate(shift);
+        place_pwlin pl(m, iso);
+
+        // Interpolated points.
+        mpoint p_0 = pl.at(mlocation{0, 0.});
+        mpoint p_1 = pl.at(mlocation{0, 0.125});
+        mpoint p_2 = pl.at(mlocation{0, 0.25});
+        mpoint p_3 = pl.at(mlocation{0, 0.5});
+
+        // Expected results.
+        mpoint x_0 = iso.apply(mpoint{0.0, 0.0, 0.0, 2.});
+        mpoint x_1 = iso.apply(mpoint{0.0, 0.0, 0.5, 2.});
+        mpoint x_2 = iso.apply(mpoint{0.0, 0.0, 1.0, 2.});
+        mpoint x_3 = iso.apply(mpoint{1.0, 0.0, 1.0, 2.});
+
+        EXPECT_TRUE(mpoint_almost_eq(x_0, p_0));
+        EXPECT_TRUE(mpoint_almost_eq(x_1, p_1));
+        EXPECT_TRUE(mpoint_almost_eq(x_2, p_2));
+        EXPECT_TRUE(mpoint_almost_eq(x_3, p_3));
+    }
+}
+
+TEST(place_pwlin, branched) {
+    using pvec = std::vector<msize_t>;
+    using svec = std::vector<msample>;
+
+    // Y-shaped branched morphology.
+    // Second branch (branch 1) tapers radius.
+
+    pvec parents = {mnpos, 0, 1, 2, 3, 4, 2, 6};
+    svec samples = {
+        // branch 0
+        {{  0,  0,  0,  2}, 1},
+        {{  0,  0,  1,  2}, 1},
+        {{  3,  0,  1,  2}, 1},
+        // branch 1
+        {{  3,  0,  1,  1.0}, 1},
+        {{  3,  1,  1,  1.0}, 1},
+        {{  3,  2,  1,  0.0}, 1},
+        // branch 2
+        {{  3,  0,  1,  2} , 1},
+        {{  3,  -1, 1,  2} , 1}
+    };
+
+    sample_tree sm(samples, parents);
+    morphology m(sm, false);
+
+    isometry iso = isometry::translate(2, 3, 4);
+    place_pwlin pl(m, iso);
+
+    // Examnine points on branch 1.
+
+    mpoint p_0 = pl.at(mlocation{1, 0.});
+    mpoint p_1 = pl.at(mlocation{1, 0.25});
+    mpoint p_2 = pl.at(mlocation{1, 0.75});
+
+    mpoint x_0 = iso.apply(mpoint{3, 0,   1, 1});
+    mpoint x_1 = iso.apply(mpoint{3, 0.5, 1, 1});
+    mpoint x_2 = iso.apply(mpoint{3, 1.5, 1, 0.5});
+
+    EXPECT_TRUE(mpoint_almost_eq(x_0, p_0));
+    EXPECT_TRUE(mpoint_almost_eq(x_1, p_1));
+    EXPECT_TRUE(mpoint_almost_eq(x_2, p_2));
+}