Extend place_pwlin interface for multivalued. (#1124)

* Add `equal_range()` method to `pw_elements`.
* Simplify `embed_pwlin`, `place_pwlin` implementations, in particular the handling of zero-length branches.
* Scale segment positions in `embed_pwlin` by divding by branch_length, instead of multiplying by the reciprocal, as the latter is a false economy after all the caveats that need to be checked to preserve ordering and bounds.
* Add new method `place_pwlin::all_at()`, which gives all the points that correspond to a single `mlocation`.
* Add new method `place_pwlin::all_segments()`, providing the set of segments and/or partial-segments that correspond to an extent.
* Add new method `place_pwlin::segments()`, which provides a minimal set of segments and/or partial-segments that map onto an extent.
* Extend unit tests to suit.

Fixes #1116 and #1068.

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_;
 };

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

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

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

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

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

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