From 7ae3f39fdfe970cff70ee74f472cd08bd9f3cfa8 Mon Sep 17 00:00:00 2001
From: Ben Cumming <bcumming@cscs.ch>
Date: Fri, 19 Mar 2021 09:46:00 +0100
Subject: [PATCH] Fix bug attaching branches to the soma in ASCII file
 descriptions (#1440)

Fix bug attaching branches to the soma in ASCII file descriptions.

Branches in an ASCII description attached to the soma were incorrectly creating a segment
from the center of the soma to the first point in the branch.
---
 arborio/neurolucida.cpp  |  21 ++++--
 test/unit/CMakeLists.txt |   2 +-
 test/unit/test_asc.cpp   | 137 ++++++++++++++++++++++++++++-----------
 3 files changed, 115 insertions(+), 45 deletions(-)

diff --git a/arborio/neurolucida.cpp b/arborio/neurolucida.cpp
index cd231c3e..705082d5 100644
--- a/arborio/neurolucida.cpp
+++ b/arborio/neurolucida.cpp
@@ -643,12 +643,22 @@ asc_morphology parse_asc_string(const char* input) {
             tails.pop_back();
 
             auto parent = head.parent_id;
+            auto& branch = *head.child;
             auto prox_sample = head.sample;
-            auto& b = *head.child;
 
-            for (auto& dist_sample: b.samples) {
-                parent = stree.append(parent, prox_sample, dist_sample, tag);
-                prox_sample = dist_sample;
+            if (!branch.samples.empty()) { // Skip empty branches, which are permitted
+                auto it = branch.samples.begin();
+                // Don't connect the first sample to the distal end of the parent
+                // branch if the parent is the soma center.
+                if (parent==arb::mnpos) {
+                    prox_sample = *it;
+                    ++it;
+                }
+                do {
+                    parent = stree.append(parent, prox_sample, *it, tag);
+                    prox_sample = *it;
+                    ++it;
+                } while (it!=branch.samples.end());
             }
 
             // Push child branches to stack in reverse order.
@@ -656,7 +666,8 @@ asc_morphology parse_asc_string(const char* input) {
             // order they were described in the file, so that segments in the
             // segment tree were added in the same order they are described
             // in the file, to give deterministic branch numbering.
-            for (auto it=b.children.rbegin(); it!=b.children.rend(); ++it) {
+            const auto& kids = branch.children;
+            for (auto it=kids.rbegin(); it!=kids.rend(); ++it) {
                 tails.push_back({&(*it), parent, prox_sample});
             }
         }
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index a6d0f519..7333b9ae 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -216,7 +216,7 @@ if(${CMAKE_POSITION_INDEPENDENT_CODE})
     MECHS dummy
     ARBOR "${PROJECT_SOURCE_DIR}"
     STANDALONE ON
-    VERBOSE ON)
+    VERBOSE OFF)
   target_compile_definitions(unit PRIVATE USE_DYNAMIC_CATALOGUES)
   add_dependencies(unit dummy-catalogue)
 endif()
diff --git a/test/unit/test_asc.cpp b/test/unit/test_asc.cpp
index 05807fa7..e1ecd098 100644
--- a/test/unit/test_asc.cpp
+++ b/test/unit/test_asc.cpp
@@ -135,10 +135,8 @@ TEST(asc, sub_tree_meta) {
     }
 }
 
-TEST(asc, branching) {
-    {
-        // Soma composed of 2 branches, and stick and fork dendrite composed of 3 branches.
-        const char* input =
+// Soma composed of 2 branches, and stick and fork dendrite composed of 3 branches.
+const char *asc_ball_and_y_dendrite =
 "((CellBody)\
  (0 0 0 2)\
 )\
@@ -151,20 +149,9 @@ TEST(asc, branching) {
   (6 5 0 1)\
  )\
  )";
-        auto result = arborio::parse_asc_string(input);
-        const auto& m = result.morphology;
-        EXPECT_EQ(m.num_branches(), 5u);
-        EXPECT_EQ(m.branch_children(0).size(), 0u);
-        EXPECT_EQ(m.branch_children(1).size(), 0u);
-        EXPECT_EQ(m.branch_children(2).size(), 2u);
-        EXPECT_EQ(m.branch_children(2)[0], 3u);
-        EXPECT_EQ(m.branch_children(2)[1], 4u);
-        EXPECT_EQ(m.branch_children(3).size(), 0u);
-        EXPECT_EQ(m.branch_children(4).size(), 0u);
-    }
-    {
-        // Soma composed of 2 branches, and a dendrite with a bit more interesting branching.
-        const char* input =
+
+// Soma composed of 2 branches, and a dendrite with a bit more interesting branching.
+const char *asc_ball_and_fancy_dendrite=
 "((CellBody)\
  (0 0 0 2)\
 )\
@@ -182,25 +169,10 @@ TEST(asc, branching) {
   (6 5 0 1)\
  )\
  )";
-        auto result = arborio::parse_asc_string(input);
-        const auto& m = result.morphology;
-        EXPECT_EQ(m.num_branches(), 7u);
-        EXPECT_EQ(m.branch_children(0).size(), 0u);
-        EXPECT_EQ(m.branch_children(1).size(), 0u);
-        EXPECT_EQ(m.branch_children(2).size(), 2u);
-        EXPECT_EQ(m.branch_children(2)[0], 3u);
-        EXPECT_EQ(m.branch_children(2)[1], 6u);
-        EXPECT_EQ(m.branch_children(3).size(), 2u);
-        EXPECT_EQ(m.branch_children(3)[0], 4u);
-        EXPECT_EQ(m.branch_children(3)[1], 5u);
-        EXPECT_EQ(m.branch_children(4).size(), 0u);
-        EXPECT_EQ(m.branch_children(5).size(), 0u);
-        EXPECT_EQ(m.branch_children(6).size(), 0u);
-    }
-    {
-        // Soma composed of 2 branches, and stick and fork dendrite and axon
-        // composed of 3 branches each.
-        const char* input =
+
+// Soma composed of 2 branches, and stick and fork dendrite and axon
+// composed of 3 branches each.
+const char* asc_ball_and_y_dendrite_and_y_axon =
 "((CellBody)\
  (0 0 0 2)\
 )\
@@ -222,7 +194,38 @@ TEST(asc, branching) {
   (6 -5 0 1)\
  )\
 )";
-        auto result = arborio::parse_asc_string(input);
+
+TEST(asc, branching) {
+    {
+        auto result = arborio::parse_asc_string(asc_ball_and_y_dendrite);
+        const auto& m = result.morphology;
+        EXPECT_EQ(m.num_branches(), 5u);
+        EXPECT_EQ(m.branch_children(0).size(), 0u);
+        EXPECT_EQ(m.branch_children(1).size(), 0u);
+        EXPECT_EQ(m.branch_children(2).size(), 2u);
+        EXPECT_EQ(m.branch_children(2)[0], 3u);
+        EXPECT_EQ(m.branch_children(2)[1], 4u);
+        EXPECT_EQ(m.branch_children(3).size(), 0u);
+        EXPECT_EQ(m.branch_children(4).size(), 0u);
+    }
+    {
+        auto result = arborio::parse_asc_string(asc_ball_and_fancy_dendrite);
+        const auto& m = result.morphology;
+        EXPECT_EQ(m.num_branches(), 7u);
+        EXPECT_EQ(m.branch_children(0).size(), 0u);
+        EXPECT_EQ(m.branch_children(1).size(), 0u);
+        EXPECT_EQ(m.branch_children(2).size(), 2u);
+        EXPECT_EQ(m.branch_children(2)[0], 3u);
+        EXPECT_EQ(m.branch_children(2)[1], 6u);
+        EXPECT_EQ(m.branch_children(3).size(), 2u);
+        EXPECT_EQ(m.branch_children(3)[0], 4u);
+        EXPECT_EQ(m.branch_children(3)[1], 5u);
+        EXPECT_EQ(m.branch_children(4).size(), 0u);
+        EXPECT_EQ(m.branch_children(5).size(), 0u);
+        EXPECT_EQ(m.branch_children(6).size(), 0u);
+    }
+    {
+        auto result = arborio::parse_asc_string(asc_ball_and_y_dendrite_and_y_axon);
         const auto& m = result.morphology;
         EXPECT_EQ(m.num_branches(), 8u);
         EXPECT_EQ(m.branch_children(0).size(), 0u);
@@ -238,6 +241,62 @@ TEST(asc, branching) {
         EXPECT_EQ(m.branch_children(6).size(), 0u);
         EXPECT_EQ(m.branch_children(7).size(), 0u);
     }
-
 }
 
+// Test that
+//  * first segment in branches connected to the soma has as its
+//    proximal sample the first sample in the branch
+//  * all other non-soma branches have the distal end of the parent
+//    branch as their proximal sample.
+//  * the soma is composed of two branches, attached to the soma center,
+//    and extending along the z axis
+TEST(asc, soma_connection) {
+    {
+        auto result = arborio::parse_asc_string(asc_ball_and_y_dendrite);
+        const auto& m = result.morphology;
+        EXPECT_EQ(m.num_branches(), 5u);
+        // Test soma
+        EXPECT_EQ(m.branch_segments(0)[0].prox, (arb::mpoint{0, 0, 0, 2}));
+        EXPECT_EQ(m.branch_segments(0)[0].dist, (arb::mpoint{0, 0,-2, 2}));
+        EXPECT_EQ(m.branch_segments(1)[0].prox, (arb::mpoint{0, 0, 0, 2}));
+        EXPECT_EQ(m.branch_segments(1)[0].dist, (arb::mpoint{0, 0, 2, 2}));
+        // Test dendrite proximal ends
+        EXPECT_EQ(m.branch_segments(2)[0].prox, (arb::mpoint{0, 2, 0, 1}));
+        EXPECT_EQ(m.branch_segments(3)[0].prox, (arb::mpoint{0, 5, 0, 1}));
+        EXPECT_EQ(m.branch_segments(4)[0].prox, (arb::mpoint{0, 5, 0, 1}));
+    }
+    {
+        auto result = arborio::parse_asc_string(asc_ball_and_fancy_dendrite);
+        const auto& m = result.morphology;
+        EXPECT_EQ(m.num_branches(), 7u);
+        // Test soma
+        EXPECT_EQ(m.branch_segments(0)[0].prox, (arb::mpoint{0, 0, 0, 2}));
+        EXPECT_EQ(m.branch_segments(0)[0].dist, (arb::mpoint{0, 0,-2, 2}));
+        EXPECT_EQ(m.branch_segments(1)[0].prox, (arb::mpoint{0, 0, 0, 2}));
+        EXPECT_EQ(m.branch_segments(1)[0].dist, (arb::mpoint{0, 0, 2, 2}));
+        // Test dendrite proximal ends
+        EXPECT_EQ(m.branch_segments(2)[0].prox, (arb::mpoint{ 0, 2, 0, 1}));
+        EXPECT_EQ(m.branch_segments(3)[0].prox, (arb::mpoint{ 0, 5, 0, 1}));
+        EXPECT_EQ(m.branch_segments(4)[0].prox, (arb::mpoint{-5, 5, 0, 1}));
+        EXPECT_EQ(m.branch_segments(5)[0].prox, (arb::mpoint{-5, 5, 0, 1}));
+        EXPECT_EQ(m.branch_segments(6)[0].prox, (arb::mpoint{ 0, 5, 0, 1}));
+    }
+    {
+        auto result = arborio::parse_asc_string(asc_ball_and_y_dendrite_and_y_axon);
+        const auto& m = result.morphology;
+        EXPECT_EQ(m.num_branches(), 8u);
+        // Test soma
+        EXPECT_EQ(m.branch_segments(0)[0].prox, (arb::mpoint{0, 0, 0, 2}));
+        EXPECT_EQ(m.branch_segments(0)[0].dist, (arb::mpoint{0, 0,-2, 2}));
+        EXPECT_EQ(m.branch_segments(1)[0].prox, (arb::mpoint{0, 0, 0, 2}));
+        EXPECT_EQ(m.branch_segments(1)[0].dist, (arb::mpoint{0, 0, 2, 2}));
+        // Test dendrite proximal ends
+        EXPECT_EQ(m.branch_segments(2)[0].prox, (arb::mpoint{0, 2, 0, 1}));
+        EXPECT_EQ(m.branch_segments(3)[0].prox, (arb::mpoint{0, 5, 0, 1}));
+        EXPECT_EQ(m.branch_segments(4)[0].prox, (arb::mpoint{0, 5, 0, 1}));
+        // Test axon proximal ends
+        EXPECT_EQ(m.branch_segments(5)[0].prox, (arb::mpoint{0,-2, 0, 1}));
+        EXPECT_EQ(m.branch_segments(6)[0].prox, (arb::mpoint{0,-5, 0, 1}));
+        EXPECT_EQ(m.branch_segments(7)[0].prox, (arb::mpoint{0,-5, 0, 1}));
+    }
+}
-- 
GitLab