From 9437bef67ab23ff7d3fdc89b53c3803cd24260cb Mon Sep 17 00:00:00 2001
From: Vasileios Karakasis <karakasis@cscs.ch>
Date: Fri, 17 Jun 2016 11:27:34 +0200
Subject: [PATCH] Re-implementation of swc_read_cell() using the basic tree
algorithms.
---
src/algorithms.hpp | 3 +-
src/swcio.cpp | 111 ++++++++++++++++---------------------------
tests/test_swcio.cpp | 1 -
3 files changed, 44 insertions(+), 71 deletions(-)
diff --git a/src/algorithms.hpp b/src/algorithms.hpp
index 40b6e9cc..5d609b1e 100644
--- a/src/algorithms.hpp
+++ b/src/algorithms.hpp
@@ -244,8 +244,9 @@ std::vector<typename C::value_type> make_parent_index(
"integral type required"
);
- if (parent_index.empty() && branch_index.empty())
+ if (parent_index.empty() && branch_index.empty()) {
return {};
+ }
EXPECTS(parent_index.size() == branch_index.back());
EXPECTS(has_contiguous_segments(parent_index));
diff --git a/src/swcio.cpp b/src/swcio.cpp
index 792d0fc7..203f69eb 100644
--- a/src/swcio.cpp
+++ b/src/swcio.cpp
@@ -1,4 +1,5 @@
#include <algorithm>
+#include <functional>
#include <iomanip>
#include <map>
#include <sstream>
@@ -230,88 +231,60 @@ swc_record_range_clean::swc_record_range_clean(std::istream &is)
}
}
-//
-// Convenience functions for returning the radii and the coordinates of a series
-// of swc records
-//
-static std::vector<swc_record::coord_type>
-swc_radii(const std::vector<swc_record> &records)
-{
- std::vector<swc_record::coord_type> radii;
- for (const auto &r : records) {
- radii.push_back(r.radius());
- }
-
- return radii;
-}
-
-static std::vector<nest::mc::point<swc_record::coord_type> >
-swc_points(const std::vector<swc_record> &records)
-{
- std::vector<nest::mc::point<swc_record::coord_type> > points;
- for (const auto &r : records) {
- points.push_back(r.coord());
- }
-
- return points;
-}
-
-static void make_cable(cell &cell,
- const std::vector<swc_record::id_type> &branch_index,
- const std::vector<swc_record> &branch_run)
-{
- auto new_parent = branch_index[branch_run.back().id()] - 1;
- cell.add_cable(new_parent, nest::mc::segmentKind::dendrite,
- swc_radii(branch_run), swc_points(branch_run));
-}
-
cell swc_read_cell(std::istream &is)
{
+ using namespace nest::mc;
+
cell newcell;
- std::vector<swc_record::id_type> parent_list;
+ std::vector<swc_record::id_type> parent_index;
std::vector<swc_record> swc_records;
for (const auto &r : swc_get_records<swc_io_clean>(is)) {
swc_records.push_back(r);
- parent_list.push_back(r.parent());
+ parent_index.push_back(r.parent());
}
- // The parent of soma must be 0
- if (!parent_list.empty()) {
- parent_list[0] = 0;
+ if (parent_index.empty()) {
+ return newcell;
}
- auto branch_index = nest::mc::algorithms::branches(parent_list);
- std::vector<swc_record> branch_run;
-
- branch_run.reserve(parent_list.size());
- auto last_branch_point = branch_index[0];
- for (auto i = 0u; i < swc_records.size(); ++i) {
- if (branch_index[i] != last_branch_point) {
- // New branch encountered; add to cell the current one
- const auto &p = branch_run.back();
- if (p.parent() == -1) {
- // This is a soma
- newcell.add_soma(p.radius(), p.coord());
- last_branch_point = i;
- } else {
- last_branch_point = i - 1;
- make_cable(newcell, branch_index, branch_run);
- }
-
- // Reset the branch run
- branch_run.clear();
- if (p.parent() != -1) {
- // Add parent of the current cell to the branch,
- // if not branching from soma
- branch_run.push_back(swc_records[parent_list[i]]);
- }
+ // The parent of soma must be 0, while in SWC files is -1
+ parent_index[0] = 0;
+ auto branch_index = algorithms::branches(parent_index);
+ auto new_parent_index = algorithms::make_parent_index(parent_index,
+ branch_index);
+
+ // sanity check
+ EXPECTS(new_parent_index.size() == branch_index.size() - 1);
+
+ // Add the soma first; then the segments
+ newcell.add_soma(swc_records[0].radius(), swc_records[0].coord());
+ for (std::size_t i = 1; i < new_parent_index.size(); ++i) {
+ auto b_start = std::next(swc_records.begin(), branch_index[i]);
+ auto b_end = std::next(swc_records.begin(), branch_index[i+1]);
+
+ std::vector<swc_record::coord_type> radii;
+ std::vector<nest::mc::point<swc_record::coord_type>> points;
+ if (new_parent_index[i] != 0) {
+ // include the parent of current record if not branching from soma
+ auto p = parent_index[branch_index[i]];
+ radii.push_back(swc_records[p].radius());
+ points.push_back(swc_records[p].coord());
}
- branch_run.push_back(swc_records[i]);
- }
+ // extract the radii and the points
+ std::for_each(b_start, b_end,
+ [&radii](const swc_record& r) {
+ radii.push_back(r.radius());
+ });
+
+ std::for_each(b_start, b_end,
+ [&points](const swc_record& r) {
+ points.push_back(r.coord());
+ });
- if (!branch_run.empty()) {
- make_cable(newcell, branch_index, branch_run);
+ // add the new cable
+ newcell.add_cable(new_parent_index[i],
+ nest::mc::segmentKind::dendrite, radii, points);
}
return newcell;
diff --git a/tests/test_swcio.cpp b/tests/test_swcio.cpp
index 62c52b09..5ba0b991 100644
--- a/tests/test_swcio.cpp
+++ b/tests/test_swcio.cpp
@@ -521,4 +521,3 @@ TEST(swc_parser, from_file_ball_and_stick)
EXPECT_TRUE(nest::mc::cell_basic_equality(local_cell, cell));
}
-
--
GitLab