Skip to content
Snippets Groups Projects
Commit a37c0431 authored by Vasileios Karakasis's avatar Vasileios Karakasis
Browse files

Construction of a cell from an swc file. 

+ some other minor improvements.
parent 393b145e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/algorithms.hpp
+ 4
− 4
View file @ a37c0431
......@@ -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
......
src/swcio.cpp
+ 94
− 5
View file @ a37c0431
......@@ -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
src/swcio.hpp
+ 23
− 12
View file @ a37c0431
......@@ -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
tests/test_algorithms.cpp
+ 21
− 21
View file @ a37c0431
......@@ -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 =
......
tests/test_swcio.cpp
+ 101
− 18
View file @ a37c0431
......@@ -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));
}
}
tests/test_tree.cpp
+ 49
− 3
View file @ a37c0431
......@@ -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));
}
{
//
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment