diff --git a/src/algorithms.hpp b/src/algorithms.hpp
index 39df831c26a8c8a0321128a19396ea6a5d092b30..76ac6856eb8bdb5132789d4110d4c785e6e52ec0 100644
--- a/src/algorithms.hpp
+++ b/src/algorithms.hpp
@@ -171,7 +171,7 @@ namespace algorithms{
}
auto num_child = child_count(parent_index);
- std::vector<typename C::value_type> branch_index(
+ std::vector<typename C::value_type> branch_runs(
parent_index.size(), 0
);
@@ -182,16 +182,16 @@ namespace algorithms{
++num_branches;
}
- branch_index[i] = num_branches;
+ branch_runs[i] = num_branches;
}
- return branch_index;
+ return branch_runs;
}
template<typename C>
std::vector<typename C::value_type> branches_fast(const C &parent_index)
{
- return branches<C, false>(parent_index);
+ return branches<C,false>(parent_index);
}
} // namespace algorithms
diff --git a/src/swcio.cpp b/src/swcio.cpp
index 796972a60e4f8b4ab5fe3595386c91732481295a..f96676449d9526154d4f7d0cd3531da62bba715e 100644
--- a/src/swcio.cpp
+++ b/src/swcio.cpp
@@ -5,7 +5,9 @@
#include <unordered_set>
#include <algorithms.hpp>
+#include <point.hpp>
#include <swcio.hpp>
+#include <util.hpp>
namespace nest {
namespace mc {
@@ -160,10 +162,10 @@ swc_record swc_parser::parse_record(std::istringstream &is)
{
auto id = parse_value_strict<int>(is, *this);
auto type = parse_value_strict<swc_record::kind>(is, *this);
- auto x = parse_value_strict<float>(is, *this);
- auto y = parse_value_strict<float>(is, *this);
- auto z = parse_value_strict<float>(is, *this);
- auto r = parse_value_strict<float>(is, *this);
+ auto x = parse_value_strict<swc_record::coord_type>(is, *this);
+ auto y = parse_value_strict<swc_record::coord_type>(is, *this);
+ auto z = parse_value_strict<swc_record::coord_type>(is, *this);
+ auto r = parse_value_strict<swc_record::coord_type>(is, *this);
auto parent_id = parse_value_strict<swc_record::id_type>(is, *this);
// Convert to zero-based, leaving parent_id as-is if -1
@@ -223,11 +225,98 @@ swc_record_range_clean::swc_record_range_clean(std::istream &is)
parent_list.push_back(records_[i].parent());
}
- if (!nest::mc::algorithms::is_contiguously_numbered(parent_list)) {
+ if (!nest::mc::algorithms::has_contiguous_segments(parent_list)) {
throw swc_parse_error("branches are not contiguously numbered", 0);
}
}
+//
+// 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)
+{
+ cell newcell;
+ std::vector<swc_record::id_type> parent_list;
+ 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());
+ }
+
+ // The parent of soma must be 0
+ if (!parent_list.empty()) {
+ parent_list[0] = 0;
+ }
+
+ 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]]);
+ }
+ }
+
+ branch_run.push_back(swc_records[i]);
+ }
+
+ if (!branch_run.empty()) {
+ make_cable(newcell, branch_index, branch_run);
+ }
+
+ return newcell;
+}
+
} // namespace io
} // namespace mc
} // namespace nest
diff --git a/src/swcio.hpp b/src/swcio.hpp
index 6a4af5f8192354f44f769dc74ee21507a8a3f096..a58bbac9aee9a32c03ae9bb33eb9289782a537b2 100644
--- a/src/swcio.hpp
+++ b/src/swcio.hpp
@@ -6,6 +6,9 @@
#include <string>
#include <vector>
+#include <cell.hpp>
+#include <point.hpp>
+
namespace nest {
namespace mc {
namespace io {
@@ -14,6 +17,7 @@ class swc_record
{
public:
using id_type = int;
+ using coord_type = double;
// More on SWC files: http://research.mssm.edu/cnic/swc.html
enum class kind {
@@ -29,7 +33,7 @@ public:
// swc records assume zero-based indexing; root's parent remains -1
swc_record(swc_record::kind type, int id,
- float x, float y, float z, float r,
+ coord_type x, coord_type y, coord_type z, coord_type r,
int parent_id)
: type_(type)
, id_(id)
@@ -119,41 +123,46 @@ public:
return parent_id_;
}
- float x() const
+ coord_type x() const
{
return x_;
}
- float y() const
+ coord_type y() const
{
return y_;
}
- float z() const
+ coord_type z() const
{
return z_;
}
- float radius() const
+ coord_type radius() const
{
return r_;
}
- float diameter() const
+ coord_type diameter() const
{
return 2*r_;
}
+ nest::mc::point<coord_type> coord() const
+ {
+ return nest::mc::point<coord_type>(x_, y_, z_);
+ }
+
void renumber(id_type new_id, std::map<id_type, id_type> &idmap);
private:
void check_consistency() const;
- kind type_; // record type
- id_type id_; // record id
- float x_, y_, z_; // record coordinates
- float r_; // record radius
- id_type parent_id_; // record parent's id
+ kind type_; // record type
+ id_type id_; // record id
+ coord_type x_, y_, z_; // record coordinates
+ coord_type r_; // record radius
+ id_type parent_id_; // record parent's id
};
@@ -398,11 +407,13 @@ struct swc_io_clean
};
template<typename T = swc_io_clean>
- typename T::record_range_type swc_get_records(std::istream &is)
+typename T::record_range_type swc_get_records(std::istream &is)
{
return typename T::record_range_type(is);
}
+cell swc_read_cell(std::istream &is);
+
} // namespace io
} // namespace mc
} // namespace nest
diff --git a/tests/test_algorithms.cpp b/tests/test_algorithms.cpp
index 05da8c0c41f8aeb39c29f35f65c814542cd1355b..7804bb60e77e80f1a318aca46aad12d977bbdd6b 100644
--- a/tests/test_algorithms.cpp
+++ b/tests/test_algorithms.cpp
@@ -195,9 +195,9 @@ TEST(algorithms, has_contiguous_segments)
// 1
// |
// 2
- // /|\
+ // /|\.
// 3 6 5
- // / \
+ // / \.
// 4 7
//
EXPECT_FALSE(
@@ -212,9 +212,9 @@ TEST(algorithms, has_contiguous_segments)
// 1
// |
// 2
- // /|\
+ // /|\.
// 3 7 5
- // / \
+ // / \.
// 4 6
//
EXPECT_TRUE(
@@ -227,11 +227,11 @@ TEST(algorithms, has_contiguous_segments)
// 0
// |
// 1
- // / \
+ // / \.
// 2 7
- // / \
+ // / \.
// 3 5
- // / \
+ // / \.
// 4 6
//
EXPECT_TRUE(
@@ -260,17 +260,17 @@ TEST(algorithms, child_count)
{
//
// 0
- // /|\
+ // /|\.
// 1 4 6
- // / | \
+ // / | \.
// 2 5 7
- // / \
+ // / \.
// 3 8
- // / \
+ // / \.
// 9 11
- // / \
+ // / \.
// 10 12
- // \
+ // \.
// 13
//
std::vector<int> parent_index =
@@ -292,17 +292,17 @@ TEST(algorithms, branches)
{
//
// 0
- // /|\
+ // /|\.
// 1 4 6
- // / | \
+ // / | \.
// 2 5 7
- // / \
+ // / \.
// 3 8
- // / \
+ // / \.
// 9 11
- // / \
+ // / \.
// 10 12
- // \
+ // \.
// 13
//
std::vector<int> parent_index =
@@ -340,9 +340,9 @@ TEST(algorithms, branches)
// 1
// |
// 2
- // / \
+ // / \.
// 3 4
- // \
+ // \.
// 5
//
std::vector<int> parent_index =
diff --git a/tests/test_swcio.cpp b/tests/test_swcio.cpp
index 635ccb442666b3662745a1f02e064bdbf81b5168..17a836c2065b0a56954d26fa3e1297cbfc16de45 100644
--- a/tests/test_swcio.cpp
+++ b/tests/test_swcio.cpp
@@ -3,10 +3,12 @@
#include <iostream>
#include <fstream>
#include <numeric>
+#include <type_traits>
#include <vector>
#include "gtest.h"
+#include <cell.hpp>
#include <swcio.hpp>
// SWC tests
@@ -250,9 +252,9 @@ TEST(swc_parser, input_cleaning)
// Check duplicates
std::stringstream is;
is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
- is << "2 1 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "2 1 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "2 1 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
EXPECT_EQ(2u, swc_get_records<swc_io_clean>(is).size());
}
@@ -261,9 +263,9 @@ TEST(swc_parser, input_cleaning)
// Check multiple trees
std::stringstream is;
is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
- is << "2 1 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
is << "3 1 14.566132 34.873772 7.857000 0.717830 -1\n";
- is << "4 1 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "4 2 14.566132 34.873772 7.857000 0.717830 1\n";
auto records = swc_get_records<swc_io_clean>(is);
EXPECT_EQ(2u, records.size());
@@ -272,9 +274,9 @@ TEST(swc_parser, input_cleaning)
{
// Check unsorted input
std::stringstream is;
- is << "3 1 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "2 1 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "4 1 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "3 2 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "4 2 14.566132 34.873772 7.857000 0.717830 1\n";
is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
std::array<swc_record::id_type, 4> expected_id_list = {{ 0, 1, 2, 3 }};
@@ -293,11 +295,11 @@ TEST(swc_parser, input_cleaning)
// Check holes in numbering
std::stringstream is;
is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
- is << "21 1 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "31 1 14.566132 34.873772 7.857000 0.717830 21\n";
- is << "41 1 14.566132 34.873772 7.857000 0.717830 21\n";
- is << "51 1 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "61 1 14.566132 34.873772 7.857000 0.717830 51\n";
+ is << "21 2 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "31 2 14.566132 34.873772 7.857000 0.717830 21\n";
+ is << "41 2 14.566132 34.873772 7.857000 0.717830 21\n";
+ is << "51 2 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "61 2 14.566132 34.873772 7.857000 0.717830 51\n";
std::array<swc_record::id_type, 6> expected_id_list =
{{ 0, 1, 2, 3, 4, 5 }};
@@ -327,9 +329,9 @@ TEST(swc_record_ranges, raw)
// Check valid usage
std::stringstream is;
is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
- is << "2 1 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "3 1 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "4 1 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "3 2 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "4 2 14.566132 34.873772 7.857000 0.717830 1\n";
std::vector<swc_record> records;
for (auto c : swc_get_records<swc_io_raw>(is)) {
@@ -381,9 +383,9 @@ TEST(swc_record_ranges, raw)
// Check parse error context
std::stringstream is;
is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
- is << "2 1 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
is << "3 10 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "4 1 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "4 2 14.566132 34.873772 7.857000 0.717830 1\n";
std::vector<swc_record> records;
try {
@@ -404,3 +406,84 @@ TEST(swc_record_ranges, raw)
EXPECT_TRUE(swc_get_records<swc_io_clean>(is).empty());
}
}
+
+TEST(swc_io, cell_construction)
+{
+ using namespace nest::mc;
+
+ {
+ //
+ // 0
+ // |
+ // 1
+ // |
+ // 2
+ // / \.
+ // 3 4
+ // \.
+ // 5
+ //
+
+ std::stringstream is;
+ is << "1 1 0 0 0 2.1 -1\n";
+ is << "2 2 0.1 1.2 1.2 1.3 1\n";
+ is << "3 2 1.0 2.0 2.2 1.1 2\n";
+ is << "4 2 1.5 3.3 1.3 2.2 3\n";
+ is << "5 2 2.5 5.3 2.5 0.7 3\n";
+ is << "6 2 3.5 2.3 3.7 3.4 5\n";
+
+ using point_type = point<double>;
+ std::vector<point_type> points = {
+ { 0.0, 0.0, 0.0 },
+ { 0.1, 1.2, 1.2 },
+ { 1.0, 2.0, 2.2 },
+ { 1.5, 3.3, 1.3 },
+ { 2.5, 5.3, 2.5 },
+ { 3.5, 2.3, 3.7 },
+ };
+
+ cell cell = io::swc_read_cell(is);
+ EXPECT_TRUE(cell.has_soma());
+ EXPECT_EQ(4, cell.num_segments());
+
+ EXPECT_EQ(norm(points[1]-points[2]), cell.cable(1)->length());
+ EXPECT_EQ(norm(points[2]-points[3]), cell.cable(2)->length());
+ EXPECT_EQ(norm(points[2]-points[4]) + norm(points[4]-points[5]),
+ cell.cable(3)->length());
+
+
+ // Check each segment separately
+ EXPECT_TRUE(cell.segment(0)->is_soma());
+ EXPECT_EQ(2.1, cell.soma()->radius());
+ EXPECT_EQ(point_type(0, 0, 0), cell.soma()->center());
+
+ for (auto i = 1; i < cell.num_segments(); ++i) {
+ EXPECT_TRUE(cell.segment(i)->is_dendrite());
+ }
+
+ EXPECT_EQ(1, cell.cable(1)->num_sub_segments());
+ EXPECT_EQ(1, cell.cable(2)->num_sub_segments());
+ EXPECT_EQ(2, cell.cable(3)->num_sub_segments());
+
+
+ // Check the radii
+ EXPECT_EQ(1.3, cell.cable(1)->radius(0));
+ EXPECT_EQ(1.1, cell.cable(1)->radius(1));
+
+ EXPECT_EQ(1.1, cell.cable(2)->radius(0));
+ EXPECT_EQ(2.2, cell.cable(2)->radius(1));
+
+ EXPECT_EQ(1.1, cell.cable(3)->radius(0));
+ EXPECT_EQ(3.4, cell.cable(3)->radius(1));
+
+ auto len_ratio = norm(points[2]-points[4]) / cell.cable(3)->length();
+ EXPECT_NEAR(.7, cell.cable(3)->radius(len_ratio), 1e-15);
+
+ // Double-check radii at joins are equal
+ EXPECT_EQ(cell.cable(1)->radius(1),
+ cell.cable(2)->radius(0));
+
+ EXPECT_EQ(cell.cable(1)->radius(1),
+ cell.cable(3)->radius(0));
+ }
+}
diff --git a/tests/test_tree.cpp b/tests/test_tree.cpp
index 9bb320d1e4cf5f1b0e8b1f5e320384f32935d728..a9b923170e14a42edab6a1deb070cb380e8d2378 100644
--- a/tests/test_tree.cpp
+++ b/tests/test_tree.cpp
@@ -30,8 +30,17 @@ TEST(cell_tree, from_parent_index) {
EXPECT_EQ(tree.num_segments(), 1u);
EXPECT_EQ(tree.num_children(0), 0u);
}
- // tree with two segments off the root node
+
{
+ //
+ // 0 0
+ // / \ / \
+ // 1 4 => 1 2
+ // / \
+ // 2 5
+ // /
+ // 3
+ //
std::vector<int> parent_index =
{0, 0, 1, 2, 0, 4};
cell_tree tree(parent_index);
@@ -43,9 +52,18 @@ TEST(cell_tree, from_parent_index) {
EXPECT_EQ(tree.num_children(2), 0u);
}
{
- // tree with three segments off the root node
+ //
+ // 0 0
+ // /|\ /|\
+ // 1 4 6 => 1 2 3
+ // / | \
+ // 2 5 7
+ // / \
+ // 3 8
+ //
std::vector<int> parent_index =
{0, 0, 1, 2, 0, 4, 0, 6, 7, 8};
+
cell_tree tree(parent_index);
EXPECT_EQ(tree.num_segments(), 4u);
// the root has 3 children
@@ -54,9 +72,29 @@ TEST(cell_tree, from_parent_index) {
EXPECT_EQ(tree.num_children(1), 0u);
EXPECT_EQ(tree.num_children(2), 0u);
EXPECT_EQ(tree.num_children(3), 0u);
+
+ // Check new structure
+ EXPECT_EQ(-1, tree.parent(0));
+ EXPECT_EQ(0, tree.parent(1));
+ EXPECT_EQ(0, tree.parent(2));
+ EXPECT_EQ(0, tree.parent(3));
}
{
- // tree with three segments off the root node, and another 2 segments off of the third branch from the root node
+ //
+ // 0 0
+ // /|\ /|\
+ // 1 4 6 => 1 2 3
+ // / | \ / \
+ // 2 5 7 4 5
+ // / \
+ // 3 8
+ // / \
+ // 9 11
+ // / \
+ // 10 12
+ // \
+ // 13
+ //
std::vector<int> parent_index =
{0, 0, 1, 2, 0, 4, 0, 6, 7, 8, 9, 8, 11, 12};
cell_tree tree(parent_index);
@@ -70,6 +108,14 @@ TEST(cell_tree, from_parent_index) {
EXPECT_EQ(tree.num_children(2), 0u);
EXPECT_EQ(tree.num_children(4), 0u);
EXPECT_EQ(tree.num_children(5), 0u);
+
+ // Check new structure
+ EXPECT_EQ(-1, tree.parent(0));
+ EXPECT_EQ(0, tree.parent(1));
+ EXPECT_EQ(0, tree.parent(2));
+ EXPECT_EQ(0, tree.parent(3));
+ EXPECT_EQ(3, tree.parent(4));
+ EXPECT_EQ(3, tree.parent(5));
}
{
//