diff --git a/arbor/include/arbor/morph/place_pwlin.hpp b/arbor/include/arbor/morph/place_pwlin.hpp
index 5be0c0a11d3163c9e96f5e685ecc07aed88f84ca..da278fab7804ee74546da648dd34ff572d312fbc 100644
--- a/arbor/include/arbor/morph/place_pwlin.hpp
+++ b/arbor/include/arbor/morph/place_pwlin.hpp
@@ -71,8 +71,19 @@ struct place_pwlin_data;
struct place_pwlin {
explicit place_pwlin(const morphology& m, const isometry& iso = isometry{});
+
+ // Any point corresponding to the location loc.
mpoint at(mlocation loc) const;
+ // All points corresponding to the location loc.
+ std::vector<mpoint> all_at(mlocation loc) const;
+
+ // A minimal set of segments or part segments whose union is coterminous with extent.
+ std::vector<msegment> segments(const mextent& extent) const;
+
+ // Maximal set of segments or part segments whose union is coterminous with extent.
+ std::vector<msegment> all_segments(const mextent& extent) const;
+
private:
std::shared_ptr<place_pwlin_data> data_;
};
diff --git a/arbor/morph/embed_pwlin.cpp b/arbor/morph/embed_pwlin.cpp
index ef5ab6500abd236f0df373bd01b30d8201cdcf0b..b46afbe2a0991d2b773e0e5bb5df7c1361a90440 100644
--- a/arbor/morph/embed_pwlin.cpp
+++ b/arbor/morph/embed_pwlin.cpp
@@ -6,7 +6,6 @@
#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"
@@ -17,6 +16,34 @@ 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>>;
+
+// Special case handling required for degenerate branches of length zero:
+template <typename Elem>
+static bool is_degenerate(const util::pw_elements<Elem>& pw) {
+ return pw.bounds().second==0;
+}
+
+template <unsigned p, unsigned q>
+double interpolate(const branch_pw_ratpoly<p, q>& f, unsigned bid, double pos) {
+ const auto& pw = f.at(bid);
+ if (is_degenerate(pw)) pos = 0;
+
+ auto piece = pw(pos);
+ auto& bounds = piece.first; // TODO: C++17 structured binding.
+ auto& element = piece.second;
+
+ 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.
//
@@ -221,11 +248,11 @@ embed_pwlin::embed_pwlin(const arb::morphology& m) {
}
double branch_length = seg_pos.back();
- double length_scale = branch_length>0? 1./branch_length: 0;
- for (auto& d: seg_pos) {
- d *= length_scale;
+ if (branch_length!=0) {
+ for (auto& d: seg_pos) {
+ d /= branch_length;
+ }
}
- seg_pos.back() = 1; // Circumvent any rounding infelicities.
for (auto d: seg_pos) {
all_segment_ends_.push_back({bid, d});
@@ -243,55 +270,44 @@ embed_pwlin::embed_pwlin(const arb::morphology& m) {
segment_cables_[seg.id] = mcable{bid, pos0, pos1};
}
-
double length_0 = parent==mnpos? 0: data_->length[parent].back().second[1];
data_->length[bid].push_back(0., 1, rat_element<1, 0>(length_0, length_0+branch_length));
double area_0 = parent==mnpos? 0: data_->area[parent].back().second[2];
double ixa_0 = parent==mnpos? 0: data_->ixa[parent].back().second[2];
- if (length_scale==0) {
- // Zero-length branch? Weird, but make best show of it.
- auto& s = segments.front().prox;
- double r = s.radius;
- double z = s.z;
- data_->radius[bid].push_back(0., 1., rat_element<1, 0>(r, r));
- data_->directed_projection[bid].push_back(0., 1., rat_element<1, 0>(z-proj_shift, z-proj_shift));
- data_->area[bid].push_back(0., 1., rat_element<2, 0>(area_0, area_0, area_0));
- data_->ixa[bid].push_back(0., 1., rat_element<1, 1>(ixa_0, ixa_0, ixa_0));
+ for (auto i: util::count_along(segments)) {
+ auto prox = segments[i].prox;
+ auto dist = segments[i].dist;
+
+ double p0 = seg_pos[i];
+ double p1 = seg_pos[i+1];
+
+ double z0 = prox.z - proj_shift;
+ double z1 = dist.z - proj_shift;
+ data_->directed_projection[bid].push_back(p0, p1, rat_element<1, 0>(z0, z1));
+
+ double r0 = prox.radius;
+ double r1 = dist.radius;
+ data_->radius[bid].push_back(p0, p1, rat_element<1, 0>(r0, r1));
+
+ double dx = (p1-p0)*branch_length;
+ double dr = r1-r0;
+ double c = pi*std::sqrt(dr*dr+dx*dx);
+ double area_half = area_0 + (0.75*r0+0.25*r1)*c;
+ double area_1 = area_0 + (r0+r1)*c;
+ data_->area[bid].push_back(p0, p1, rat_element<2, 0>(area_0, area_half, area_1));
+ area_0 = area_1;
+
+ // (Test for positive dx explicitly in case r0 is zero.)
+ double ixa_half = ixa_0 + (dx>0? dx/(pi*r0*(r0+r1)): 0);
+ double ixa_1 = ixa_0 + (dx>0? dx/(pi*r0*r1): 0);
+ data_->ixa[bid].push_back(p0, p1, rat_element<1, 1>(ixa_0, ixa_half, ixa_1));
+ ixa_0 = ixa_1;
}
- else {
- for (auto i: util::count_along(segments)) {
- auto prox = segments[i].prox;
- auto dist = segments[i].dist;
-
- double p0 = seg_pos[i];
- double p1 = seg_pos[i+1];
- if (p0==p1) continue;
-
- double z0 = prox.z - proj_shift;
- double z1 = dist.z - proj_shift;
- data_->directed_projection[bid].push_back(p0, p1, rat_element<1, 0>(z0, z1));
-
- double r0 = prox.radius;
- double r1 = dist.radius;
- data_->radius[bid].push_back(p0, p1, rat_element<1, 0>(r0, r1));
-
- double dx = (p1-p0)*branch_length;
- double dr = r1-r0;
- double c = pi*std::sqrt(dr*dr+dx*dx);
- double area_half = area_0 + (0.75*r0+0.25*r1)*c;
- double area_1 = area_0 + (r0+r1)*c;
- data_->area[bid].push_back(p0, p1, rat_element<2, 0>(area_0, area_half, area_1));
- area_0 = area_1;
-
- double ixa_half = ixa_0 + dx/(pi*r0*(r0+r1));
- double ixa_1 = ixa_0 + dx/(pi*r0*r1);
- data_->ixa[bid].push_back(p0, p1, rat_element<1, 1>(ixa_0, ixa_half, ixa_1));
- ixa_0 = ixa_1;
- }
- arb_assert((data_->radius[bid].size()>0));
+ arb_assert((data_->radius[bid].size()>0));
+ if (branch_length!=0) {
arb_assert((data_->radius[bid].bounds()==std::pair<double, double>(0., 1.)));
arb_assert((data_->area[bid].bounds()==std::pair<double, double>(0., 1.)));
arb_assert((data_->ixa[bid].bounds()==std::pair<double, double>(0., 1.)));
diff --git a/arbor/morph/place_pwlin.cpp b/arbor/morph/place_pwlin.cpp
index 9fe9afdb1643e88886fc1e8820fe27cec00dcb4e..136fcf90dcec081c79b4877ec15c58f9f4830f06 100644
--- a/arbor/morph/place_pwlin.cpp
+++ b/arbor/morph/place_pwlin.cpp
@@ -6,7 +6,6 @@
#include <arbor/morph/place_pwlin.hpp>
#include <arbor/morph/primitives.hpp>
-#include "morph/pwlin_common.hpp"
#include "util/piecewise.hpp"
#include "util/rangeutil.hpp"
#include "util/ratelem.hpp"
@@ -17,18 +16,111 @@ namespace arb {
using util::rat_element;
struct place_pwlin_data {
- branch_pw_ratpoly<1, 0> x, y, z, r; // [µm]
+ // Piecewise-constant indices into segment data, by branch.
+ std::vector<util::pw_elements<std::size_t>> segment_index;
+
+ // Segments from segment tree, after isometry is applied.
+ std::vector<msegment> segments;
explicit place_pwlin_data(msize_t n_branch):
- x(n_branch), y(n_branch), z(n_branch), r(n_branch)
+ segment_index(n_branch)
{}
};
+static mpoint interpolate_segment(const std::pair<double, double>& bounds, const msegment& seg, double pos) {
+ if (bounds.first==bounds.second) {
+ return seg.prox;
+ }
+ else {
+ double u = (pos-bounds.first)/(bounds.second-bounds.first);
+
+ util::rat_element<1, 0> x{seg.prox.x, seg.dist.x};
+ util::rat_element<1, 0> y{seg.prox.y, seg.dist.y};
+ util::rat_element<1, 0> z{seg.prox.z, seg.dist.z};
+ util::rat_element<1, 0> r{seg.prox.radius, seg.dist.radius};
+
+ return {x(u), y(u), z(u), r(u)};
+ }
+}
+
+template <typename Elem>
+static bool is_degenerate(const util::pw_elements<Elem>& pw) {
+ return pw.bounds().second==0;
+}
+
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) };
+ const auto& index = data_->segment_index.at(loc.branch);
+ double pos = is_degenerate(index)? 0: loc.pos;
+
+ auto piece = index(pos); // Here and below, TODO: C++17 structured bindings for piece.first and second.
+ return interpolate_segment(piece.first, data_->segments.at(piece.second), pos);
+}
+
+std::vector<mpoint> place_pwlin::all_at(mlocation loc) const {
+ std::vector<mpoint> result;
+ const auto& index = data_->segment_index.at(loc.branch);
+ double pos = is_degenerate(index)? 0: loc.pos;
+
+ for (auto piece: util::make_range(index.equal_range(pos))) {
+ result.push_back(interpolate_segment(piece.first, data_->segments.at(piece.second), pos));
+ }
+ return result;
+}
+
+template <bool exclude_trivial>
+static std::vector<msegment> extent_segments_impl(const place_pwlin_data& data, const mextent& extent) {
+ std::vector<msegment> result;
+
+ for (mcable c: extent) {
+ const auto& index = data.segment_index.at(c.branch);
+ if (is_degenerate(index)) {
+ c.prox_pos = c.dist_pos = 0;
+ }
+
+ auto b = index.equal_range(c.prox_pos).first;
+ auto e = index.equal_range(c.dist_pos).second;
+
+ for (auto piece: util::make_range(b, e)) {
+ const auto& bounds = piece.first;
+ const msegment& seg = data.segments.at(piece.second);
+
+ auto partial_bounds = bounds;
+ msegment partial = seg;
+
+ if (c.prox_pos>bounds.first) {
+ arb_assert(c.prox_pos<=bounds.second);
+ partial.prox = interpolate_segment(bounds, seg, c.prox_pos);
+ partial_bounds.first = c.prox_pos;
+ }
+ if (c.dist_pos<bounds.second) {
+ arb_assert(c.dist_pos>=bounds.first);
+ partial.dist = interpolate_segment(bounds, seg, c.dist_pos);
+ partial_bounds.second = c.dist_pos;
+ }
+
+ // With exclude_trivial set, skip zero-length segments if cable is non-trivial.
+ if (exclude_trivial && partial_bounds.first==partial_bounds.second && c.prox_pos!=c.dist_pos) {
+ continue;
+ }
+
+ result.push_back(partial);
+
+ // With exclude_trivial set, keep only one zero-length (partial) segment if cable is trivial.
+ if (exclude_trivial && c.prox_pos==c.dist_pos) {
+ break;
+ }
+ }
+ }
+
+ return result;
+}
+
+std::vector<msegment> place_pwlin::all_segments(const mextent& extent) const {
+ return extent_segments_impl<false>(*data_, extent);
+}
+
+std::vector<msegment> place_pwlin::segments(const mextent& extent) const {
+ return extent_segments_impl<true>(*data_, extent);
}
place_pwlin::place_pwlin(const arb::morphology& m, const isometry& iso) {
@@ -37,8 +129,6 @@ place_pwlin::place_pwlin(const arb::morphology& m, const isometry& iso) {
if (!n_branch) return;
- std::vector<double> sample_pos_on_branch;
-
std::vector<double> seg_pos;
for (msize_t bid = 0; bid<n_branch; ++bid) {
auto& segments = m.branch_segments(bid);
@@ -51,35 +141,22 @@ place_pwlin::place_pwlin(const arb::morphology& m, const isometry& iso) {
}
double branch_length = seg_pos.back();
- double length_scale = branch_length>0? 1./branch_length: 0;
- for (auto& x: seg_pos) {
- x *= length_scale;
+ if (branch_length!=0) {
+ for (auto& x: seg_pos) {
+ x /= branch_length;
+ }
}
- if (length_scale==0) {
- // Zero-length branch case?
- mpoint p = iso.apply(segments[0].prox);
+ for (auto i: util::count_along(segments)) {
+ msegment seg = segments[i];
+ double p0 = seg_pos[i];
+ double p1 = seg_pos[i+1];
- 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(segments)) {
- auto& seg = segments[i];
- double p0 = seg_pos[i];
- double p1 = seg_pos[i+1];
- if (p0==p1) continue;
-
- mpoint u0 = iso.apply(seg.prox);
- mpoint u1 = iso.apply(seg.dist);
-
- 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));
- }
+ seg.prox = iso.apply(seg.prox);
+ seg.dist = iso.apply(seg.dist);
+
+ data_->segment_index[bid].push_back(p0, p1, data_->segments.size());
+ data_->segments.push_back(seg);
}
}
};
diff --git a/arbor/morph/pwlin_common.hpp b/arbor/morph/pwlin_common.hpp
deleted file mode 100644
index 54d1d983f25b6e8c3b549d79ca37cea342dec056..0000000000000000000000000000000000000000
--- a/arbor/morph/pwlin_common.hpp
+++ /dev/null
@@ -1,38 +0,0 @@
-#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/arbor/util/piecewise.hpp b/arbor/util/piecewise.hpp
index 4e4d42d1e7cd777a7abf27c8283524e67d7cdefa..ef067d62cdaa0e8fc64ef2b364bbdccf618f48d2 100644
--- a/arbor/util/piecewise.hpp
+++ b/arbor/util/piecewise.hpp
@@ -138,6 +138,7 @@ struct pw_elements {
value_type front() const { return (*this)[0]; }
value_type back() const { return (*this)[size()-1]; }
+ // Return index of right-most element whose corresponding closed interval contains x.
size_type index_of(double x) const {
if (empty()) return npos;
@@ -146,6 +147,17 @@ struct pw_elements {
else return partn.index(x);
}
+ // Return iterator pair spanning elements whose corresponding closed intervals contain x.
+ std::pair<iterator, iterator> equal_range(double x) const {
+ auto eq = std::equal_range(vertex_.begin(), vertex_.end(), x);
+
+ if (eq.first==vertex_.end()) return {end(), end()};
+ if (eq.first>vertex_.begin()) --eq.first;
+ if (eq.second==vertex_.end()) --eq.second;
+
+ return {begin()+(eq.first-vertex_.begin()), begin()+(eq.second-vertex_.begin())};
+ }
+
value_type operator()(double x) const {
size_type i = index_of(x);
if (i==npos) {
diff --git a/test/unit/test_morph_place.cpp b/test/unit/test_morph_place.cpp
index dc7bebb5eb3a743fc05519a9e1e75cbc68fe4490..6d02a26e73cf87ed60ce10d6c3b1218f3194cda2 100644
--- a/test/unit/test_morph_place.cpp
+++ b/test/unit/test_morph_place.cpp
@@ -7,9 +7,9 @@
#include <arbor/morph/primitives.hpp>
#include "util/piecewise.hpp"
+#include "util/rangeutil.hpp"
#include "../test/gtest.h"
-#include "common.hpp"
#include "common_cells.hpp"
using namespace arb;
@@ -291,3 +291,253 @@ TEST(place_pwlin, branched) {
EXPECT_TRUE(mpoint_almost_eq(x_1, p_1));
EXPECT_TRUE(mpoint_almost_eq(x_2, p_2));
}
+
+TEST(place_pwlin, all_at) {
+ // One branch, two discontinguous segments.
+ {
+ mpoint p0p{0, 0, 0, 1};
+ mpoint p0d{2, 3, 4, 5};
+ mpoint p1p{3, 4, 5, 7};
+ mpoint p1d{6, 6, 7, 8};
+
+ segment_tree tree;
+ msize_t s0 = tree.append(mnpos, p0p, p0d, 0);
+ msize_t s1 = tree.append(s0, p1p, p1d, 0);
+
+ morphology m(tree);
+ mprovider p(m, label_dict{});
+ place_pwlin place(m);
+
+ mlocation_list end_s0 = thingify(ls::most_distal(reg::segment(s0)), p);
+ ASSERT_EQ(1u, end_s0.size());
+ mlocation_list end_s1 = thingify(ls::most_distal(reg::segment(s1)), p);
+ ASSERT_EQ(1u, end_s1.size());
+
+ auto points_end_s1 = place.all_at(end_s1[0]);
+ ASSERT_EQ(1u, points_end_s1.size());
+ EXPECT_TRUE(mpoint_almost_eq(p1d, points_end_s1[0]));
+
+ auto points_end_s0 = place.all_at(end_s0[0]);
+ ASSERT_EQ(2u, points_end_s0.size());
+ EXPECT_TRUE(mpoint_almost_eq(p0d, points_end_s0[0]));
+ EXPECT_TRUE(mpoint_almost_eq(p1p, points_end_s0[1]));
+ }
+
+ // One branch, multiple zero-length segments at end.
+ {
+ mpoint p0p{0, 0, 0, 1};
+ mpoint p0d{2, 3, 4, 5};
+ mpoint p1p{3, 4, 5, 7};
+ mpoint p1d = p1p;
+ mpoint p2p{6, 6, 7, 8};
+ mpoint p2d = p2p;
+
+ segment_tree tree;
+ msize_t s0 = tree.append(mnpos, p0p, p0d, 0);
+ msize_t s1 = tree.append(s0, p1p, p1d, 0);
+ (void)tree.append(s1, p2p, p2d, 0);
+
+ morphology m(tree);
+ mprovider p(m, label_dict{});
+ place_pwlin place(m);
+
+ auto points_end_b0 = place.all_at(mlocation{0, 1});
+ ASSERT_EQ(3u, points_end_b0.size());
+ EXPECT_TRUE(mpoint_almost_eq(p0d, points_end_b0[0]));
+ EXPECT_TRUE(mpoint_almost_eq(p1d, points_end_b0[1]));
+ EXPECT_TRUE(mpoint_almost_eq(p2d, points_end_b0[2]));
+ }
+
+ // Zero length branch comprising multiple zero-length segments.
+ // (Please don't do this.)
+ {
+ mpoint p0p{2, 3, 4, 5};
+ mpoint p0d = p0p;
+ mpoint p1p{3, 4, 5, 7};
+ mpoint p1d = p1p;
+ mpoint p2p{6, 6, 7, 8};
+ mpoint p2d = p2p;
+
+ segment_tree tree;
+ msize_t s0 = tree.append(mnpos, p0p, p0d, 0);
+ msize_t s1 = tree.append(s0, p1p, p1d, 0);
+ (void)tree.append(s1, p2p, p2d, 0);
+
+ morphology m(tree);
+ mprovider p(m, label_dict{});
+ place_pwlin place(m);
+
+ auto points_begin_b0 = place.all_at(mlocation{0, 0});
+ ASSERT_EQ(3u, points_begin_b0.size());
+ EXPECT_TRUE(mpoint_almost_eq(p0d, points_begin_b0[0]));
+ EXPECT_TRUE(mpoint_almost_eq(p1d, points_begin_b0[1]));
+ EXPECT_TRUE(mpoint_almost_eq(p2d, points_begin_b0[2]));
+
+ auto points_mid_b0 = place.all_at(mlocation{0, 0.5});
+ ASSERT_EQ(3u, points_mid_b0.size());
+ EXPECT_TRUE(mpoint_almost_eq(p0d, points_mid_b0[0]));
+ EXPECT_TRUE(mpoint_almost_eq(p1d, points_mid_b0[1]));
+ EXPECT_TRUE(mpoint_almost_eq(p2d, points_mid_b0[2]));
+
+ auto points_end_b0 = place.all_at(mlocation{0, 1});
+ ASSERT_EQ(3u, points_end_b0.size());
+ EXPECT_TRUE(mpoint_almost_eq(p0d, points_end_b0[0]));
+ EXPECT_TRUE(mpoint_almost_eq(p1d, points_end_b0[1]));
+ EXPECT_TRUE(mpoint_almost_eq(p2d, points_end_b0[2]));
+ }
+}
+
+TEST(place_pwlin, segments) {
+ // Y-shaped morphology with some discontinuous
+ // and zero-length segments.
+
+ segment_tree tree;
+
+ // root branch
+
+ mpoint p0p{0, 0, 0, 1};
+ mpoint p0d{1, 0, 0, 1};
+ mpoint p1p = p0d;
+ mpoint p1d{2, 0, 0, 1};
+
+ msize_t s0 = tree.append(mnpos, p0p, p0d, 0);
+ msize_t s1 = tree.append(s0, p1p, p1d, 0);
+
+ // branch A (total length 2)
+
+ mpoint p2p{2, 0, 0, 1};
+ mpoint p2d{2, 1, 0, 1};
+ mpoint p3p{2, 1, 0, 0.5};
+ mpoint p3d{2, 2, 0, 0.5};
+ mpoint p4p{8, 9, 7, 1.5}; // some random zero-length segments on the end...
+ mpoint p4d = p4p;
+ mpoint p5p{3, 1, 3, 0.5};
+ mpoint p5d = p5p;
+
+ msize_t s2 = tree.append(s1, p2p, p2d, 0);
+ msize_t s3 = tree.append(s2, p3p, p3d, 0);
+ msize_t s4 = tree.append(s3, p4p, p4d, 0);
+ msize_t s5 = tree.append(s4, p5p, p5d, 0);
+ (void)s5;
+
+ // branch B (total length 4)
+
+ mpoint p6p{2, 0, 0, 1};
+ mpoint p6d{2, 0, 2, 1};
+ mpoint p7p{2, 0, 2, 0.5}; // a zero-length segment in the middle...
+ mpoint p7d = p7p;
+ mpoint p8p{2, 0, 3, 1};
+ mpoint p8d{2, 0, 5, 1};
+
+ msize_t s6 = tree.append(s1, p6p, p6d, 0);
+ msize_t s7 = tree.append(s6, p7p, p7d, 0);
+ msize_t s8 = tree.append(s7, p8p, p8d, 0);
+ (void)s8;
+
+ morphology m(tree);
+ mprovider p(m, label_dict{});
+ place_pwlin place(m);
+
+ // Thingify a segment expression to work out which branch id is A and
+ // which is B.
+
+ mextent s2_extent = thingify(reg::segment(2), p);
+ msize_t branch_A = s2_extent.front().branch;
+
+ mextent s6_extent = thingify(reg::segment(6), p);
+ msize_t branch_B = s6_extent.front().branch;
+
+ ASSERT_TRUE((branch_A==1 && branch_B==2) || (branch_A==2 && branch_B==1));
+
+ // Region 1: all of branch A, middle section of branch B.
+
+ region r1 = join(reg::branch(branch_A), reg::cable(branch_B, 0.25, 0.75));
+ mextent x1 = thingify(r1, p);
+
+ std::vector<msegment> x1min = place.segments(x1);
+ std::vector<msegment> x1all = place.all_segments(x1);
+
+ auto seg_id = [](const msegment& s) { return s.id; };
+
+ util::sort_by(x1min, seg_id);
+ std::vector<msize_t> x1min_seg_ids = util::assign_from(util::transform_view(x1min, seg_id));
+
+ util::sort_by(x1all, seg_id);
+ std::vector<msize_t> x1all_seg_ids = util::assign_from(util::transform_view(x1all, seg_id));
+
+ ASSERT_EQ((std::vector<msize_t>{2, 3, 6, 8}), x1min_seg_ids);
+ ASSERT_EQ((std::vector<msize_t>{2, 3, 4, 5, 6, 7, 8}), x1all_seg_ids);
+
+ EXPECT_TRUE(mpoint_almost_eq(p2p, x1all[0].prox));
+ EXPECT_TRUE(mpoint_almost_eq(p2d, x1all[0].dist));
+
+ EXPECT_TRUE(mpoint_almost_eq(p3p, x1all[1].prox));
+ EXPECT_TRUE(mpoint_almost_eq(p3d, x1all[1].dist));
+
+ EXPECT_TRUE(mpoint_almost_eq(p4p, x1all[2].prox));
+ EXPECT_TRUE(mpoint_almost_eq(p4d, x1all[2].dist));
+
+ EXPECT_TRUE(mpoint_almost_eq(p5p, x1all[3].prox));
+ EXPECT_TRUE(mpoint_almost_eq(p5d, x1all[3].dist));
+
+ EXPECT_FALSE(mpoint_almost_eq(p6p, x1all[4].prox));
+ EXPECT_TRUE(mpoint_almost_eq(p6d, x1all[4].dist));
+ EXPECT_TRUE(mpoint_almost_eq(mpoint{2, 0, 1, 1}, x1all[4].prox));
+
+ EXPECT_TRUE(mpoint_almost_eq(p7p, x1all[5].prox));
+ EXPECT_TRUE(mpoint_almost_eq(p7d, x1all[5].dist));
+
+ EXPECT_TRUE(mpoint_almost_eq(p8p, x1all[6].prox));
+ EXPECT_FALSE(mpoint_almost_eq(p8d, x1all[6].dist));
+ EXPECT_TRUE(mpoint_almost_eq(mpoint{2, 0, 4, 1}, x1all[6].dist));
+
+ // Region 2: first half of branch A. Test exclusion of zero-length partial
+ // segments.
+
+ region r2 = reg::cable(branch_A, 0., 0.5);
+ mextent x2 = thingify(r2, p);
+
+ std::vector<msegment> x2min = place.segments(x2);
+ std::vector<msegment> x2all = place.all_segments(x2);
+
+ util::sort_by(x2min, seg_id);
+ std::vector<msize_t> x2min_seg_ids = util::assign_from(util::transform_view(x2min, seg_id));
+
+ util::sort_by(x2all, seg_id);
+ std::vector<msize_t> x2all_seg_ids = util::assign_from(util::transform_view(x2all, seg_id));
+
+ ASSERT_EQ((std::vector<msize_t>{2}), x2min_seg_ids);
+ ASSERT_EQ((std::vector<msize_t>{2, 3}), x2all_seg_ids);
+
+ EXPECT_TRUE(mpoint_almost_eq(p3p, x2all[1].prox));
+ EXPECT_TRUE(mpoint_almost_eq(p3p, x2all[1].dist));
+
+ // Region 3: trivial region, midpont of branch B.
+
+ region r3 = reg::cable(branch_B, 0.5, 0.5);
+ mextent x3 = thingify(r3, p);
+
+ std::vector<msegment> x3min = place.segments(x3);
+ std::vector<msegment> x3all = place.all_segments(x3);
+
+ util::sort_by(x3min, seg_id);
+ std::vector<msize_t> x3min_seg_ids = util::assign_from(util::transform_view(x3min, seg_id));
+
+ util::sort_by(x3all, seg_id);
+ std::vector<msize_t> x3all_seg_ids = util::assign_from(util::transform_view(x3all, seg_id));
+
+ ASSERT_EQ(1u, x3min_seg_ids.size()); // Could be end of s6, all of s7, or beginning of s8
+ ASSERT_EQ((std::vector<msize_t>{6, 7, 8}), x3all_seg_ids);
+
+ EXPECT_TRUE(mpoint_almost_eq(x3min[0].prox, x3min[0].dist));
+ EXPECT_TRUE(mpoint_almost_eq(p6d, x3min[0].prox) ||
+ mpoint_almost_eq(p7d, x3min[0].prox) ||
+ mpoint_almost_eq(p8p, x3min[0].prox));
+
+ EXPECT_TRUE(mpoint_almost_eq(p6d, x3all[0].prox));
+ EXPECT_TRUE(mpoint_almost_eq(p6d, x3all[0].dist));
+ EXPECT_TRUE(mpoint_almost_eq(p7d, x3all[1].prox));
+ EXPECT_TRUE(mpoint_almost_eq(p7d, x3all[1].dist));
+ EXPECT_TRUE(mpoint_almost_eq(p8p, x3all[2].prox));
+ EXPECT_TRUE(mpoint_almost_eq(p8p, x3all[2].dist));
+}
diff --git a/test/unit/test_piecewise.cpp b/test/unit/test_piecewise.cpp
index 03047be6c591964c5bdc28c9d0ac9c799de24db9..65b74361eb74186dffab9ce942ba6310a04560ff 100644
--- a/test/unit/test_piecewise.cpp
+++ b/test/unit/test_piecewise.cpp
@@ -231,6 +231,67 @@ TEST(piecewise, index_of) {
EXPECT_EQ(pw_npos, v0.index_of(0.));
}
+TEST(piecewise, equal_range) {
+ {
+ pw_elements<int> p{{1, 2, 3, 4}, {10, 9, 8}};
+
+ auto er0 = p.equal_range(0.0);
+ ASSERT_EQ(er0.first, er0.second);
+
+ auto er1 = p.equal_range(1.0);
+ ASSERT_EQ(1u, er1.second-er1.first);
+ EXPECT_EQ(10, er1.first->second);
+
+ auto er2 = p.equal_range(2.0);
+ ASSERT_EQ(2u, er2.second-er2.first);
+ auto iter = er2.first;
+ EXPECT_EQ(10, iter++->second);
+ EXPECT_EQ(9, iter->second);
+
+ auto er3_5 = p.equal_range(3.5);
+ ASSERT_EQ(1u, er3_5.second-er3_5.first);
+ EXPECT_EQ(8, er3_5.first->second);
+
+ auto er4 = p.equal_range(4.0);
+ ASSERT_EQ(1u, er4.second-er4.first);
+ EXPECT_EQ(8, er4.first->second);
+
+ auto er5 = p.equal_range(5.0);
+ ASSERT_EQ(er5.first, er5.second);
+ }
+
+ {
+ pw_elements<int> p{{1, 1, 2, 2, 2, 3, 3}, {10, 11, 12, 13, 14, 15}};
+
+ auto er0 = p.equal_range(0.0);
+ ASSERT_EQ(er0.first, er0.second);
+
+ auto er1 = p.equal_range(1.0);
+ ASSERT_EQ(2u, er1.second-er1.first);
+ auto iter = er1.first;
+ EXPECT_EQ(10, iter++->second);
+ EXPECT_EQ(11, iter++->second);
+
+ auto er2 = p.equal_range(2.0);
+ ASSERT_EQ(4u, er2.second-er2.first);
+ iter = er2.first;
+ EXPECT_EQ(11, iter++->second);
+ EXPECT_EQ(12, iter++->second);
+ EXPECT_EQ(13, iter++->second);
+ EXPECT_EQ(14, iter++->second);
+
+ auto er3 = p.equal_range(3.0);
+ ASSERT_EQ(2u, er3.second-er3.first);
+ iter = er3.first;
+ EXPECT_EQ(14, iter++->second);
+ EXPECT_EQ(15, iter++->second);
+
+ auto er5 = p.equal_range(5.0);
+ ASSERT_EQ(er5.first, er5.second);
+ }
+
+}
+
TEST(piecewise, push) {
pw_elements<int> q;
using dp = std::pair<double, double>;