diff --git a/arbor/morph/region.cpp b/arbor/morph/region.cpp
index f0536d6110cfb6fefc38ef6305dd2e9e7245880b..abe9342771007eab6f6f5a4e267241bc1b42a0a7 100644
--- a/arbor/morph/region.cpp
+++ b/arbor/morph/region.cpp
@@ -132,6 +132,44 @@ mcable_list remove_cover(mcable_list cables, const morphology& m) {
return merge(cables);
}
+// Remove zero-length regions from a sorted cable_list that are in the cover of another region in the list.
+mcable_list remove_covered_points(mcable_list cables, const morphology& m) {
+ struct branch_index_pair {
+ msize_t bid;
+ unsigned lid;
+ };
+ std::vector<branch_index_pair> end_branches;
+ std::vector<unsigned> erase_indices;
+
+ for (unsigned i = 0; i < cables.size(); i++) {
+ auto c = cables[i];
+ // Save zero length cables at the distal end of a branch
+ // Or at the proximal end of the soma
+ if ((c.prox_pos==1 && c.dist_pos==1) ||
+ (c.prox_pos==0 && c.dist_pos==0)) {
+ end_branches.push_back({c.branch, i});
+ }
+ // Look for branches that are children of the cables saved in end_branches
+ else if (c.prox_pos==0) {
+ auto parent = m.branch_parent(c.branch);
+ if (parent==mnpos) parent = 0;
+
+ auto it = std::find_if(end_branches.begin(), end_branches.end(),
+ [&](branch_index_pair& p) { return p.bid==parent;});
+
+ if (it!=end_branches.end()) {
+ erase_indices.push_back((*it).lid);
+ end_branches.erase(it);
+ }
+ }
+ }
+ util::sort(erase_indices);
+ for (auto it = erase_indices.rbegin(); it != erase_indices.rend(); it++) {
+ cables.erase(cables.begin() + *it);
+ }
+
+ return cables;
+}
// Empty region.
@@ -324,7 +362,7 @@ mcable_list thingify_(const reg_and& P, const mprovider& p) {
at_end = it.first==end.first || it.second==end.second;
}
- return remove_cover(L, m);
+ return remove_covered_points(remove_cover(L, m), m);
}
std::ostream& operator<<(std::ostream& o, const reg_and& x) {
diff --git a/test/unit/test_morph_expr.cpp b/test/unit/test_morph_expr.cpp
index 6304baa47e765974f9f727cc2765865f85626c7f..bee279ab1c1fb47da5759aaa0c20c59d133848e6 100644
--- a/test/unit/test_morph_expr.cpp
+++ b/test/unit/test_morph_expr.cpp
@@ -17,6 +17,13 @@
using namespace arb;
using embedding = embed_pwlin;
+namespace arb {
+ namespace reg {
+ mcable_list remove_covered_points(mcable_list cables, const morphology& m);
+
+ }
+}
+
::testing::AssertionResult cable_eq(mcable a, mcable b) {
if (a.branch!=b.branch) {
return ::testing::AssertionFailure()
@@ -471,5 +478,121 @@ TEST(region, thingify) {
std::swap(lhs, rhs);
EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand);
EXPECT_EQ(thingify(join(lhs, rhs), mp), ror);
+
+ // b1
+ // 123456789
+ // |----- | lhs
+ // |----- | rhs
+ // |xxxxx | rand
+ // |xxxxx | ror
+ lhs = cable({1,0,.5});
+ rhs = cable({1,0,.5});
+ rand = cl{{1,0,.5}};
+ ror = cl{{1,0,.5}};
+ EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand);
+ EXPECT_EQ(thingify(join(lhs, rhs), mp), ror);
+
+ // b3
+ // 123456789
+ // |----- | lhs
+ // |----- | rhs
+ // |xxxxx | rand
+ // |xxxxx | ror
+ lhs = cable({3,0,.5});
+ rhs = cable({3,0,.5});
+ rand = cl{{3,0,.5}};
+ ror = cl{{3,0,.5}};
+ EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand);
+ EXPECT_EQ(thingify(join(lhs, rhs), mp), ror);
+
+ // b0 b1
+ // 123456789 123456789
+ // |xxxxx | | lhs
+ // | |xxxxx | rhs
+ // x | | rand
+ // |xxxxx |xxxxx | ror
+ lhs = cable({0,0,.5});
+ rhs = cable({1,0,.5});
+ rand = cl{{0,0,0}};
+ ror = cl{{0,0,.5},{1,0,.5}};
+ EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand);
+ EXPECT_EQ(thingify(join(lhs, rhs), mp), ror);
+
+ // Assert communtativity
+ std::swap(lhs, rhs);
+ EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand);
+ EXPECT_EQ(thingify(join(lhs, rhs), mp), ror);
+
+ // b2 b3
+ // 123456789 123456789
+ // |xxxxx | | lhs
+ // | |xxxxx | rhs
+ // x | | rand
+ // |xxxxx |xxxxx | ror
+ lhs = cable({2,0,.5});
+ rhs = cable({3,0,.5});
+ rand = cl{{1,1,1}};
+ ror = cl{{2,0,.5},{3,0,.5}};
+ EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand);
+ EXPECT_EQ(thingify(join(lhs, rhs), mp), ror);
+
+ // Assert communtativity
+ std::swap(lhs, rhs);
+ EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand);
+ EXPECT_EQ(thingify(join(lhs, rhs), mp), ror);
+
+ // b0 b3
+ // 123456789 123456789
+ // |xxxxx |xxxxx | lhs
+ // |xxxxxxx |xxx | rhs
+ // |xxxxx |xxx | rand
+ // |xxxxxxx |xxxxx | ror
+ lhs = join(cable({0,0,.5}), cable({3,0,.5}));
+ rhs = join(cable({0,0,.7}), cable({3,0,.3}));
+ rand = cl{{0,0,.5},{3,0,.3}};
+ ror = cl{{0,0,.7},{3,0,.5}};
+ EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand);
+ EXPECT_EQ(thingify(join(lhs, rhs), mp), ror);
+
+ // Assert communtativity
+ std::swap(lhs, rhs);
+ EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand);
+ EXPECT_EQ(thingify(join(lhs, rhs), mp), ror);
+
+ }
+ {
+ pvec parents = {mnpos, 0, 1, 0, 3, 4, 5, 5, 7, 7, 4, 10};
+ svec samples = {
+ {{ 0, 0, 0, 2}, 3}, //0
+ {{ 10, 0, 0, 2}, 3}, //1
+ {{100, 0, 0, 2}, 3}, //2
+ {{ 0, 10, 0, 2}, 3}, //3
+ {{ 0,100, 0, 2}, 3}, //4
+ {{100,100, 0, 2}, 3}, //5
+ {{200,100, 0, 2}, 3}, //6
+ {{100,200, 0, 2}, 3}, //7
+ {{200,200, 0, 2}, 3}, //8
+ {{100,300, 0, 2}, 3}, //9
+ {{ 0,200, 0, 2}, 3}, //10
+ {{ 0,300, 0, 2}, 3}, //11
+ };
+ sample_tree sm(samples, parents);
+ auto m = morphology(sm, false);
+ std::cout << m.branch_parent(7);
+
+ {
+ auto in = cl{{0,0,0},{1,0,0.5},{1,1,1},{2,0,1},{2,1,1},{3,1,1},{4,0,1},{5,1,1},{7,0,1}};
+ auto out = reg::remove_covered_points(in, m);
+
+ auto expected = cl{{1,0,0.5},{2,0,1},{3,1,1},{4,0,1},{5,1,1},{7,0,1}};
+ EXPECT_TRUE(cablelist_eq(out, expected));
+ }
+ {
+ auto in = cl{{0,0,0},{1,0,0.5},{1,1,1},{2,1,1},{3,1,1},{4,0,1},{5,1,1},{7,0,1}};
+ auto out = reg::remove_covered_points(in, m);
+
+ auto expected = cl{{1,0,0.5},{3,1,1},{4,0,1},{5,1,1},{7,0,1}};
+ EXPECT_TRUE(cablelist_eq(out, expected));
+ }
}
}