diff --git a/src/algorithms.hpp b/src/algorithms.hpp
index bc25b98b84518e0a0c444bf6163bd05b8ed37711..8ac62055f351212fb75c72657fa6157463450f45 100644
--- a/src/algorithms.hpp
+++ b/src/algorithms.hpp
@@ -5,6 +5,7 @@
#include <algorithm>
#include <numeric>
#include <type_traits>
+#include <vector>
/*
* Some simple wrappers around stl algorithms to improve readability of code
@@ -117,6 +118,89 @@ namespace algorithms{
return true;
}
+ template<typename C>
+ bool has_contiguous_segments(const C &parent_index)
+ {
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "integral type required"
+ );
+
+ if (!is_minimal_degree(parent_index)) {
+ return false;
+ }
+
+ std::vector<bool> is_leaf(parent_index.size(), false);
+
+ for (std::size_t i = 1; i < parent_index.size(); ++i) {
+ auto p = parent_index[i];
+ if (is_leaf[p]) {
+ return false;
+ }
+
+ if (p != i-1) {
+ // we have a branch and i-1 is a leaf node
+ is_leaf[i-1] = true;
+ }
+ }
+
+ return true;
+ }
+
+ template<typename C>
+ std::vector<typename C::value_type> child_count(const C &parent_index)
+ {
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "integral type required"
+ );
+
+ std::vector<typename C::value_type> count(parent_index.size(), 0);
+ for (std::size_t i = 1; i < parent_index.size(); ++i) {
+ ++count[parent_index[i]];
+ }
+
+ return count;
+ }
+
+ template<typename C, bool CheckStrict = true>
+ std::vector<typename C::value_type> branches(const C &parent_index)
+ {
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "integral type required"
+ );
+
+ if (CheckStrict && !has_contiguous_segments(parent_index)) {
+ throw std::invalid_argument(
+ "parent_index has not contiguous branch numbering"
+ );
+ }
+
+ auto num_child = child_count(parent_index);
+ std::vector<typename C::value_type> branch_runs(
+ parent_index.size(), 0
+ );
+
+ std::size_t num_branches = (num_child[0] == 1) ? 1 : 0;
+ for (std::size_t i = 1; i < parent_index.size(); ++i) {
+ auto p = parent_index[i];
+ if (num_child[p] > 1) {
+ ++num_branches;
+ }
+
+ branch_runs[i] = num_branches;
+ }
+
+ return branch_runs;
+ }
+
+ template<typename C>
+ std::vector<typename C::value_type> branches_fast(const C &parent_index)
+ {
+ return branches<C,false>(parent_index);
+ }
+
} // namespace algorithms
} // namespace mc
} // namespace nest
diff --git a/src/point.hpp b/src/point.hpp
index cf74f3238d705e4e28c9430bb42cc85aa2b1ee14..3b048bb3db7fb580443e019c98292821f715ffb2 100644
--- a/src/point.hpp
+++ b/src/point.hpp
@@ -73,6 +73,20 @@ dot(
return lhs.x*rhs.x + lhs.y*rhs.y + lhs.z*rhs.z;
}
+template<typename T>
+bool operator==(const point<T> &rhs, const point<T> &lhs)
+{
+ return (rhs.x == lhs.x) &&
+ (rhs.y == lhs.y) &&
+ (rhs.z == lhs.z);
+}
+
+template<typename T>
+bool operator!=(const point<T> &rhs, const point<T> &lhs)
+{
+ return !(rhs == lhs);
+}
+
} // namespace mc
} // namespace nest
@@ -81,5 +95,3 @@ std::ostream& operator << (std::ostream& o, nest::mc::point<T> const& p)
{
return o << "[" << p.x << ", " << p.y << ", " << p.z << "]";
}
-
-
diff --git a/src/swcio.cpp b/src/swcio.cpp
index 98db848e2b4768d59399a3340c0d3611104fb4e9..f96676449d9526154d4f7d0cd3531da62bba715e 100644
--- a/src/swcio.cpp
+++ b/src/swcio.cpp
@@ -4,7 +4,10 @@
#include <sstream>
#include <unordered_set>
+#include <algorithms.hpp>
+#include <point.hpp>
#include <swcio.hpp>
+#include <util.hpp>
namespace nest {
namespace mc {
@@ -32,7 +35,8 @@ void swc_record::check_consistency() const
{
// Check record type as well; enum's do not offer complete type safety,
// since you can cast anything that fits to its underlying type
- if (type_ < 0 || type_ > custom) {
+ if (static_cast<int>(type_) < 0 ||
+ static_cast<int>(type_) > static_cast<int>(kind::custom)) {
throw std::invalid_argument("unknown record type");
}
@@ -65,7 +69,7 @@ std::ostream &operator<<(std::ostream &os, const swc_record &record)
{
// output in one-based indexing
os << record.id_+1 << " "
- << record.type_ << " "
+ << static_cast<int>(record.type_) << " "
<< std::setprecision(7) << record.x_ << " "
<< std::setprecision(7) << record.y_ << " "
<< std::setprecision(7) << record.z_ << " "
@@ -158,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
@@ -177,26 +181,26 @@ swc_record_range_clean::swc_record_range_clean(std::istream &is)
{
std::unordered_set<swc_record::id_type> ids;
- std::size_t num_trees = 0;
+ std::size_t num_trees = 0;
swc_record::id_type last_id = -1;
- bool needsort = false;
+ bool needsort = false;
swc_record curr_record;
- for (auto c : swc_get_records<swc_io_raw>(is)) {
- if (c.parent() == -1 && ++num_trees > 1) {
+ for (auto r : swc_get_records<swc_io_raw>(is)) {
+ if (r.parent() == -1 && ++num_trees > 1) {
// only a single tree is allowed
break;
}
- auto inserted = ids.insert(c.id());
+ auto inserted = ids.insert(r.id());
if (inserted.second) {
// not a duplicate; insert record
- records_.push_back(c);
- if (!needsort && c.id() < last_id) {
+ records_.push_back(r);
+ if (!needsort && r.id() < last_id) {
needsort = true;
}
- last_id = c.id();
+ last_id = r.id();
}
}
@@ -207,13 +211,110 @@ swc_record_range_clean::swc_record_range_clean(std::istream &is)
// Renumber records if necessary
std::map<swc_record::id_type, swc_record::id_type> idmap;
swc_record::id_type next_id = 0;
- for (auto &c : records_) {
- if (c.id() != next_id) {
- c.renumber(next_id, idmap);
+ for (auto &r : records_) {
+ if (r.id() != next_id) {
+ r.renumber(next_id, idmap);
}
++next_id;
}
+
+ // Reject if branches are not contiguously numbered
+ std::vector<swc_record::id_type> parent_list = { 0 };
+ for (std::size_t i = 1; i < records_.size(); ++i) {
+ parent_list.push_back(records_[i].parent());
+ }
+
+ 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
diff --git a/src/swcio.hpp b/src/swcio.hpp
index 0486a2b511e7afcb13fbc13f7a5977fa4f4db6a4..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,12 +17,10 @@ class swc_record
{
public:
using id_type = int;
+ using coord_type = double;
- // FIXME: enum's are not completely type-safe, since they can accept
- // anything that can be casted to their underlying type.
- //
// More on SWC files: http://research.mssm.edu/cnic/swc.html
- enum kind {
+ enum class kind {
undefined = 0,
soma,
axon,
@@ -31,9 +32,9 @@ public:
};
// swc records assume zero-based indexing; root's parent remains -1
- swc_record(kind type, int id,
- float x, float y, float z, float r,
- int parent_id)
+ swc_record(swc_record::kind type, int id,
+ coord_type x, coord_type y, coord_type z, coord_type r,
+ int parent_id)
: type_(type)
, id_(id)
, x_(x)
@@ -46,7 +47,7 @@ public:
}
swc_record()
- : type_(swc_record::undefined)
+ : type_(swc_record::kind::undefined)
, id_(0)
, x_(0)
, y_(0)
@@ -122,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
};
@@ -401,12 +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 cfe058a3c8b92d930c739fbbde3f5a3bad6f634c..6503c07d9a413bf2e90f8cbf151b0df1e3e5949b 100644
--- a/tests/test_algorithms.cpp
+++ b/tests/test_algorithms.cpp
@@ -169,3 +169,196 @@ TEST(algorithms, is_positive)
)
);
}
+
+TEST(algorithms, has_contiguous_segments)
+{
+ //
+ // 0
+ // |
+ // 1
+ // |
+ // 2
+ // /|\.
+ // 3 7 4
+ // / \.
+ // 5 6
+ //
+ EXPECT_FALSE(
+ nest::mc::algorithms::has_contiguous_segments(
+ std::vector<int>{0, 0, 1, 2, 2, 3, 4, 2}
+ )
+ );
+
+ //
+ // 0
+ // |
+ // 1
+ // |
+ // 2
+ // /|\.
+ // 3 6 5
+ // / \.
+ // 4 7
+ //
+ EXPECT_FALSE(
+ nest::mc::algorithms::has_contiguous_segments(
+ std::vector<int>{0, 0, 1, 2, 3, 2, 2, 5}
+ )
+ );
+
+ //
+ // 0
+ // |
+ // 1
+ // |
+ // 2
+ // /|\.
+ // 3 7 5
+ // / \.
+ // 4 6
+ //
+ EXPECT_TRUE(
+ nest::mc::algorithms::has_contiguous_segments(
+ std::vector<int>{0, 0, 1, 2, 3, 2, 5, 2}
+ )
+ );
+
+ //
+ // 0
+ // |
+ // 1
+ // / \.
+ // 2 7
+ // / \.
+ // 3 5
+ // / \.
+ // 4 6
+ //
+ EXPECT_TRUE(
+ nest::mc::algorithms::has_contiguous_segments(
+ std::vector<int>{0, 0, 1, 2, 3, 2, 5, 1}
+ )
+ );
+
+ // Soma-only list
+ EXPECT_TRUE(
+ nest::mc::algorithms::has_contiguous_segments(
+ std::vector<int>{0}
+ )
+ );
+
+ // Empty list
+ EXPECT_TRUE(
+ nest::mc::algorithms::has_contiguous_segments(
+ std::vector<int>{}
+ )
+ );
+}
+
+TEST(algorithms, child_count)
+{
+ {
+ //
+ // 0
+ // /|\.
+ // 1 4 6
+ // / | \.
+ // 2 5 7
+ // / \.
+ // 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 };
+ std::vector<int> expected_child_count =
+ { 3, 1, 1, 0, 1, 0, 1, 1, 2, 1, 0, 1, 1, 0 };
+
+ // auto count = nest::mc::algorithms::child_count(parent_index);
+ EXPECT_EQ(expected_child_count,
+ nest::mc::algorithms::child_count(parent_index));
+ }
+
+}
+
+TEST(algorithms, branches)
+{
+ using namespace nest::mc;
+
+ {
+ //
+ // 0
+ // /|\.
+ // 1 4 6
+ // / | \.
+ // 2 5 7
+ // / \.
+ // 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 };
+ std::vector<int> expected_branches =
+ { 0, 1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 5, 5, 5 };
+
+ auto actual_branches = algorithms::branches_fast(parent_index);
+ EXPECT_EQ(expected_branches, actual_branches);
+ }
+
+ {
+ //
+ // 0
+ // |
+ // 1
+ // |
+ // 2
+ // |
+ // 3
+ //
+ std::vector<int> parent_index =
+ { 0, 0, 1, 2 };
+ std::vector<int> expected_branches =
+ { 0, 1, 1, 1 };
+
+ auto actual_branches = algorithms::branches_fast(parent_index);
+ EXPECT_EQ(expected_branches, actual_branches);
+ }
+
+ {
+ //
+ // 0
+ // |
+ // 1
+ // |
+ // 2
+ // / \.
+ // 3 4
+ // \.
+ // 5
+ //
+ std::vector<int> parent_index =
+ { 0, 0, 1, 2, 2, 4 };
+ std::vector<int> expected_branches =
+ { 0, 1, 1, 2, 3, 3 };
+
+ auto actual_branches = algorithms::branches_fast(parent_index);
+ EXPECT_EQ(expected_branches, actual_branches);
+ }
+
+ {
+ std::vector<int> parent_index = { 0 };
+ std::vector<int> expected_branches = { 0 };
+
+ auto actual_branches = algorithms::branches_fast(parent_index);
+ EXPECT_EQ(expected_branches, actual_branches);
+ }
+}
diff --git a/tests/test_point.cpp b/tests/test_point.cpp
index 9666a1387f895d408778f40afc6b07db5363bec6..b5bfeb9c1783a6c921ca97f3eb6ae4afd535c311 100644
--- a/tests/test_point.cpp
+++ b/tests/test_point.cpp
@@ -12,6 +12,9 @@ TEST(point, construction)
{
// default constructor
point<float> p;
+
+ EXPECT_FALSE(p.is_set());
+
// expect NaN, which returns false when comparing for equality
EXPECT_NE(p.x, p.x);
EXPECT_NE(p.y, p.y);
@@ -36,10 +39,7 @@ TEST(point, construction)
// copy constructor
auto p2 = p1;
-
- EXPECT_EQ(p1.x, p2.x);
- EXPECT_EQ(p1.y, p2.y);
- EXPECT_EQ(p1.z, p2.z);
+ EXPECT_EQ(p1, p2);
}
}
@@ -64,9 +64,7 @@ TEST(point, addition)
auto p = p1 + p2;
- EXPECT_EQ(p.x, 2);
- EXPECT_EQ(p.y, 4);
- EXPECT_EQ(p.z, 6);
+ EXPECT_EQ(point<double>(2, 4, 6), p);
}
TEST(point, subtraction)
@@ -77,9 +75,7 @@ TEST(point, subtraction)
auto p = p1 - p2;
- EXPECT_EQ(p.x, 0);
- EXPECT_EQ(p.y, 0);
- EXPECT_EQ(p.z, 0);
+ EXPECT_EQ(point<double>(0, 0, 0), p);
}
TEST(point, scalar_prod)
@@ -89,9 +85,7 @@ TEST(point, scalar_prod)
auto p = 0.5 * p1;
- EXPECT_EQ(p.x, 0.5);
- EXPECT_EQ(p.y, 1.0);
- EXPECT_EQ(p.z, 1.5);
+ EXPECT_EQ(point<double>(0.5, 1.0, 1.5), p);
}
TEST(point, norm)
@@ -114,4 +108,3 @@ TEST(point, dot)
EXPECT_EQ(dot(p1,p2), 2.);
}
-
diff --git a/tests/test_swcio.cpp b/tests/test_swcio.cpp
index 3de060ed4a69c446167d07c372ee0c872255ee30..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
@@ -36,43 +38,43 @@ TEST(swc_record, construction)
{
// invalid id
EXPECT_THROW(swc_record record(
- swc_record::custom, -3, 1., 1., 1., 1., 5),
+ swc_record::kind::custom, -3, 1., 1., 1., 1., 5),
std::invalid_argument);
}
{
// invalid parent id
EXPECT_THROW(swc_record record(
- swc_record::custom, 0, 1., 1., 1., 1., -5),
+ swc_record::kind::custom, 0, 1., 1., 1., 1., -5),
std::invalid_argument);
}
{
// invalid radius
EXPECT_THROW(swc_record record(
- swc_record::custom, 0, 1., 1., 1., -1., -1),
+ swc_record::kind::custom, 0, 1., 1., 1., -1., -1),
std::invalid_argument);
}
{
// parent_id > id
EXPECT_THROW(swc_record record(
- swc_record::custom, 0, 1., 1., 1., 1., 2),
+ swc_record::kind::custom, 0, 1., 1., 1., 1., 2),
std::invalid_argument);
}
{
// parent_id == id
EXPECT_THROW(swc_record record(
- swc_record::custom, 0, 1., 1., 1., 1., 0),
+ swc_record::kind::custom, 0, 1., 1., 1., 1., 0),
std::invalid_argument);
}
{
// check standard construction by value
- swc_record record(swc_record::custom, 0, 1., 1., 1., 1., -1);
+ swc_record record(swc_record::kind::custom, 0, 1., 1., 1., 1., -1);
EXPECT_EQ(record.id(), 0);
- EXPECT_EQ(record.type(), swc_record::custom);
+ EXPECT_EQ(record.type(), swc_record::kind::custom);
EXPECT_EQ(record.x(), 1.);
EXPECT_EQ(record.y(), 1.);
EXPECT_EQ(record.z(), 1.);
@@ -83,7 +85,7 @@ TEST(swc_record, construction)
{
// check copy constructor
- swc_record record_orig(swc_record::custom, 0, 1., 1., 1., 1., -1);
+ swc_record record_orig(swc_record::kind::custom, 0, 1., 1., 1., 1., -1);
swc_record record(record_orig);
expect_record_equals(record_orig, record);
}
@@ -95,9 +97,9 @@ TEST(swc_record, comparison)
{
// check comparison operators
- swc_record record0(swc_record::custom, 0, 1., 1., 1., 1., -1);
- swc_record record1(swc_record::custom, 0, 2., 3., 4., 5., -1);
- swc_record record2(swc_record::custom, 1, 2., 3., 4., 5., -1);
+ swc_record record0(swc_record::kind::custom, 0, 1., 1., 1., 1., -1);
+ swc_record record1(swc_record::kind::custom, 0, 2., 3., 4., 5., -1);
+ swc_record record2(swc_record::kind::custom, 1, 2., 3., 4., 5., -1);
EXPECT_EQ(record0, record1);
EXPECT_LT(record0, record2);
EXPECT_GT(record2, record1);
@@ -131,6 +133,31 @@ TEST(swc_parser, invalid_input)
swc_record record;
EXPECT_THROW(is >> record, swc_parse_error);
}
+
+ {
+ // Non-contiguous numbering in branches is considered invalid
+ // 1
+ // / \
+ // 2 3
+ // /
+ // 4
+ 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 2\n";
+
+ std::vector<swc_record> records;
+ try {
+ for (auto c : swc_get_records<swc_io_clean>(is)) {
+ records.push_back(std::move(c));
+ }
+
+ FAIL() << "expected swc_parse_error, none was thrown\n";
+ } catch (const swc_parse_error &e) {
+ SUCCEED();
+ }
+ }
}
@@ -162,7 +189,7 @@ TEST(swc_parser, valid_input)
swc_record record;
EXPECT_NO_THROW(is >> record);
EXPECT_EQ(0, record.id()); // zero-based indexing
- EXPECT_EQ(swc_record::soma, record.type());
+ EXPECT_EQ(swc_record::kind::soma, record.type());
EXPECT_FLOAT_EQ(14.566132, record.x());
EXPECT_FLOAT_EQ(34.873772, record.y());
EXPECT_FLOAT_EQ( 7.857000, record.z());
@@ -173,9 +200,9 @@ TEST(swc_parser, valid_input)
{
// check valid input with a series of records
std::vector<swc_record> records_orig = {
- swc_record(swc_record::soma, 0,
+ swc_record(swc_record::kind::soma, 0,
14.566132, 34.873772, 7.857000, 0.717830, -1),
- swc_record(swc_record::dendrite, 1,
+ swc_record(swc_record::kind::dendrite, 1,
14.566132+1, 34.873772+1, 7.857000+1, 0.717830+1, -1)
};
@@ -225,37 +252,37 @@ 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(is).size());
+ EXPECT_EQ(2u, swc_get_records<swc_io_clean>(is).size());
}
{
// 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(is);
+ auto records = swc_get_records<swc_io_clean>(is);
EXPECT_EQ(2u, records.size());
}
{
// 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 }};
auto expected_id = expected_id_list.cbegin();
- for (auto c : swc_get_records(is)) {
+ for (auto c : swc_get_records<swc_io_clean>(is)) {
EXPECT_EQ(*expected_id, c.id());
++expected_id;
}
@@ -268,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 }};
@@ -281,7 +308,7 @@ TEST(swc_parser, input_cleaning)
auto expected_id = expected_id_list.cbegin();
auto expected_parent = expected_parent_list.cbegin();
- for (auto c : swc_get_records(is)) {
+ for (auto c : swc_get_records<swc_io_clean>(is)) {
EXPECT_EQ(*expected_id, c.id());
EXPECT_EQ(*expected_parent, c.parent());
++expected_id;
@@ -302,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)) {
@@ -356,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 {
@@ -379,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));
}
{
//