Remove covered zero-regions from region intersections (#934)

Stop region intersection from returning extra zero length cables if the location it describes is in the cover of another region in the `cable_list`.

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

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