diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt
index 48715d7958f54b8bf0245bfb45f497772bb8469a..42a66fc7693b13594690bbcb0a4c97f2e7b0d257 100644
--- a/arbor/CMakeLists.txt
+++ b/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
diff --git a/arbor/include/arbor/morph/place_pwlin.hpp b/arbor/include/arbor/morph/place_pwlin.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5be0c0a11d3163c9e96f5e685ecc07aed88f84ca
--- /dev/null
+++ b/arbor/include/arbor/morph/place_pwlin.hpp
@@ -0,0 +1,82 @@
+#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
+
+
diff --git a/arbor/include/arbor/morph/primitives.hpp b/arbor/include/arbor/morph/primitives.hpp
index a77864fa643242b08655dc25004ade29852b8149..7fc8e87b256ebd9c0e507d0f3cc979f0c764eab8 100644
--- a/arbor/include/arbor/morph/primitives.hpp
+++ b/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);
diff --git a/arbor/morph/embed_pwlin.cpp b/arbor/morph/embed_pwlin.cpp
index f20bfb272f6917d63581eed8bba5974748fe5406..64971a03e62d5f7b6780ca4cc7e4eeade5b3aae3 100644
--- a/arbor/morph/embed_pwlin.cpp
+++ b/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.
//
diff --git a/arbor/morph/place_pwlin.cpp b/arbor/morph/place_pwlin.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4f7a78ef55ea7f2daf752afd26f6142c5f4475ec
--- /dev/null
+++ b/arbor/morph/place_pwlin.cpp
@@ -0,0 +1,144 @@
+#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
diff --git a/arbor/morph/pwlin_common.hpp b/arbor/morph/pwlin_common.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..54d1d983f25b6e8c3b549d79ca37cea342dec056
--- /dev/null
+++ b/arbor/morph/pwlin_common.hpp
@@ -0,0 +1,38 @@
+#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
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 3ae498dfa8b32a8a2cea72d4decc3f4ac1b11ea4..f914964ad32669b19af14639e333dca9cf6f2997 100644
--- a/test/unit/CMakeLists.txt
+++ b/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
diff --git a/test/unit/test_morph_place.cpp b/test/unit/test_morph_place.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7ba2308efc8a54301f58b375b3ceb3ca3c94f502
--- /dev/null
+++ b/test/unit/test_morph_place.cpp
@@ -0,0 +1,293 @@
+#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));
+}