diff --git a/CMakeLists.txt b/CMakeLists.txt
index 717eb124c4eb1e08d455bdf73446fabec96a6265..ea981fe4e2f67149ce0f349175cc99d9e8346d17 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -11,11 +11,11 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -pthread -Wall")
set(CMAKE_EXPORT_COMPILE_COMMANDS "YES")
# generated .a and .so go into /lib
-#set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
-#set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
+set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
+set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
include_directories(${CMAKE_SOURCE_DIR}/src)
include_directories(${CMAKE_SOURCE_DIR})
+add_subdirectory(src)
add_subdirectory(tests)
-
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..90c9d8f8c4446615bfd3e045cd64a082a0d1d016
--- /dev/null
+++ b/src/CMakeLists.txt
@@ -0,0 +1,8 @@
+set(HEADERS
+ swcio.hpp
+ )
+set(BASE_SOURCES
+ swcio.cpp
+)
+
+add_library(cellalgo ${BASE_SOURCES} ${HEADERS})
diff --git a/src/swcio.cpp b/src/swcio.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..129f70afe97cf6a1ca78475205f3d4ac75c54e78
--- /dev/null
+++ b/src/swcio.cpp
@@ -0,0 +1,214 @@
+#include <iomanip>
+#include <map>
+#include <sstream>
+#include <unordered_set>
+#include <swcio.hpp>
+
+namespace nestmc
+{
+namespace io
+{
+
+//
+// cell_record implementation
+//
+void cell_record::renumber(id_type new_id, std::map<id_type, id_type> &idmap)
+{
+ auto old_id = id_;
+ id_ = new_id;
+
+ // Obtain parent_id from the map
+ auto new_parent_id = idmap.find(parent_id_);
+ if (new_parent_id != idmap.end()) {
+ parent_id_ = new_parent_id->second;
+ }
+
+ check_consistency();
+ idmap.insert(std::make_pair(old_id, new_id));
+}
+
+void cell_record::check_consistency() const
+{
+ // Check cell 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) {
+ throw std::invalid_argument("unknown cell type");
+ }
+
+ if (id_ < 0) {
+ throw std::invalid_argument("negative ids not allowed");
+ }
+
+ if (parent_id_ < -1) {
+ throw std::invalid_argument("parent_id < -1 not allowed");
+ }
+
+ if (parent_id_ >= id_) {
+ throw std::invalid_argument("parent_id >= id is not allowed");
+ }
+
+ if (r_ < 0) {
+ throw std::invalid_argument("negative radii are not allowed");
+ }
+}
+
+std::istream &operator>>(std::istream &is, cell_record &cell)
+{
+ swc_parser parser;
+ parser.parse_record(is, cell);
+ return is;
+}
+
+
+std::ostream &operator<<(std::ostream &os, const cell_record &cell)
+{
+ // output in one-based indexing
+ os << cell.id_+1 << " "
+ << cell.type_ << " "
+ << std::setprecision(7) << cell.x_ << " "
+ << std::setprecision(7) << cell.y_ << " "
+ << std::setprecision(7) << cell.z_ << " "
+ << std::setprecision(7) << cell.r_ << " "
+ << ((cell.parent_id_ == -1) ? cell.parent_id_ : cell.parent_id_+1);
+
+ return os;
+}
+
+
+//
+// Utility functions
+//
+
+bool starts_with(const std::string &str, const std::string &prefix)
+{
+ return (str.find(prefix) == 0);
+}
+
+void check_parse_status(const std::istream &is)
+{
+ if (is.fail())
+ // If we try to read past the eof; fail bit will also be set
+ throw swc_parse_error("could not parse value");
+}
+
+template<typename T>
+T parse_value_strict(std::istream &is)
+{
+ T val;
+ check_parse_status(is >> val);
+
+ // everything's fine
+ return val;
+}
+
+// specialize parsing for cell types
+template<>
+cell_record::kind parse_value_strict(std::istream &is)
+{
+ cell_record::id_type val;
+ check_parse_status(is >> val);
+
+ // Let cell_record's constructor check for the type validity
+ return static_cast<cell_record::kind>(val);
+}
+
+//
+// swc_parser implementation
+//
+
+std::istream &swc_parser::parse_record(std::istream &is, cell_record &cell)
+{
+ while (!is.eof() && !is.bad()) {
+ // consume empty and comment lines first
+ std::getline(is, linebuff_);
+ if (!linebuff_.empty() && !starts_with(linebuff_, comment_prefix_))
+ break;
+ }
+
+ if (is.bad()) {
+ // let the caller check for such events
+ return is;
+ }
+
+ if (is.eof() &&
+ (linebuff_.empty() || starts_with(linebuff_, comment_prefix_))) {
+ // last line is either empty or a comment; don't parse anything
+ return is;
+ }
+
+ if (is.fail()) {
+ throw swc_parse_error("too long line detected");
+ }
+
+ std::istringstream line(linebuff_);
+ cell = parse_record(line);
+ return is;
+}
+
+cell_record swc_parser::parse_record(std::istringstream &is)
+{
+ auto id = parse_value_strict<int>(is);
+ auto type = parse_value_strict<cell_record::kind>(is);
+ auto x = parse_value_strict<float>(is);
+ auto y = parse_value_strict<float>(is);
+ auto z = parse_value_strict<float>(is);
+ auto r = parse_value_strict<float>(is);
+ auto parent_id = parse_value_strict<cell_record::id_type>(is);
+
+ // Convert to zero-based, leaving parent_id as-is if -1
+ if (parent_id != -1) {
+ parent_id--;
+ }
+
+ return cell_record(type, id-1, x, y, z, r, parent_id);
+}
+
+
+std::vector<cell_record> swc_read_cells(std::istream &is)
+{
+ std::vector<cell_record> cells;
+ std::unordered_set<cell_record::id_type> ids;
+
+ std::size_t num_trees = 0;
+ cell_record::id_type last_id = -1;
+ bool needsort = false;
+
+ cell_record curr_cell;
+ while ( !(is >> curr_cell).eof()) {
+ if (curr_cell.parent() == -1 && ++num_trees > 1) {
+ // only a single tree is allowed
+ break;
+ }
+
+ auto inserted = ids.insert(curr_cell.id());
+ if (inserted.second) {
+ // not a duplicate; insert cell
+ cells.push_back(curr_cell);
+ if (!needsort && curr_cell.id() < last_id) {
+ needsort = true;
+ }
+
+ last_id = curr_cell.id();
+ }
+ }
+
+ if (needsort) {
+ std::sort(cells.begin(), cells.end());
+ }
+
+ // Renumber cells if necessary
+ std::map<cell_record::id_type, cell_record::id_type> idmap;
+ cell_record::id_type next_id = 0;
+ for (auto &c : cells) {
+ if (c.id() != next_id) {
+ c.renumber(next_id, idmap);
+ }
+
+ ++next_id;
+ }
+
+ return std::move(cells);
+}
+
+} // end of nestmc::io
+} // end of nestmc
diff --git a/src/swcio.hpp b/src/swcio.hpp
index d362109dd15b771ac897bb884380682214f113b0..ba52720a7f1852c88cb60ae51b0f1aee9c10f9b9 100644
--- a/src/swcio.hpp
+++ b/src/swcio.hpp
@@ -2,27 +2,23 @@
#include <exception>
#include <iostream>
-#include <sstream>
-#include <type_traits>
+#include <string>
+#include <vector>
-namespace neuron
+namespace nestmc
{
namespace io
{
-static bool starts_with(const std::string &str, const std::string &prefix)
-{
- return (str.find(prefix) == 0);
-}
-
class cell_record
{
public:
+ using id_type = int;
- // FIXME: enum's are not completely type-safe, since they can accept anything
- // that can be casted to their underlying type.
+ // 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 {
@@ -48,22 +44,7 @@ public:
, r_(r)
, parent_id_(parent_id)
{
- // Check cell 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)
- throw std::invalid_argument("unknown cell type");
-
- if (id_ < 0)
- throw std::invalid_argument("negative ids not allowed");
-
- if (parent_id_ < -1)
- throw std::invalid_argument("parent_id < -1 not allowed");
-
- if (parent_id_ >= id_)
- throw std::invalid_argument("parent_id >= id is not allowed");
-
- if (r_ < 0)
- throw std::invalid_argument("negative radii are not allowed");
+ check_consistency();
}
cell_record()
@@ -79,18 +60,43 @@ public:
cell_record(const cell_record &other) = default;
cell_record &operator=(const cell_record &other) = default;
+ // Equality and comparison operators
friend bool operator==(const cell_record &lhs,
const cell_record &rhs)
{
return lhs.id_ == rhs.id_;
}
+ friend bool operator<(const cell_record &lhs,
+ const cell_record &rhs)
+ {
+ return lhs.id_ < rhs.id_;
+ }
+
+ friend bool operator<=(const cell_record &lhs,
+ const cell_record &rhs)
+ {
+ return (lhs < rhs) || (lhs == rhs);
+ }
+
friend bool operator!=(const cell_record &lhs,
const cell_record &rhs)
{
return !(lhs == rhs);
}
+ friend bool operator>(const cell_record &lhs,
+ const cell_record &rhs)
+ {
+ return !(lhs < rhs) && (lhs != rhs);
+ }
+
+ friend bool operator>=(const cell_record &lhs,
+ const cell_record &rhs)
+ {
+ return !(lhs < rhs);
+ }
+
friend std::ostream &operator<<(std::ostream &os, const cell_record &cell);
kind type() const
@@ -98,12 +104,12 @@ public:
return type_;
}
- int id() const
+ id_type id() const
{
return id_;
}
- int parent() const
+ id_type parent() const
{
return parent_id_;
}
@@ -133,12 +139,28 @@ public:
return 2*r_;
}
+ void renumber(id_type new_id, std::map<id_type, id_type> &idmap);
+
private:
+ void check_consistency() const;
+
kind type_; // cell type
- int id_; // cell id
+ id_type id_; // cell id
float x_, y_, z_; // cell coordinates
float r_; // cell radius
- int parent_id_; // cell parent's id
+ id_type parent_id_; // cell parent's id
+};
+
+class swc_parse_error : public std::runtime_error
+{
+public:
+ explicit swc_parse_error(const char *msg)
+ : std::runtime_error(msg)
+ { }
+
+ explicit swc_parse_error(const std::string &msg)
+ : std::runtime_error(msg)
+ { }
};
class swc_parser
@@ -155,51 +177,9 @@ public:
, comment_prefix_("#")
{ }
- std::istream &parse_record(std::istream &is, cell_record &cell)
- {
- while (!is.eof() && !is.bad()) {
- // consume empty and comment lines first
- std::getline(is, linebuff_);
- if (!linebuff_.empty() && !starts_with(linebuff_, comment_prefix_))
- break;
- }
-
- if (is.bad())
- // let the caller check for such events
- return is;
-
- if (is.eof() &&
- (linebuff_.empty() || starts_with(linebuff_, comment_prefix_)))
- // last line is either empty or a comment; don't parse anything
- return is;
-
- if (is.fail())
- throw std::runtime_error("too long line detected");
-
- std::istringstream line(linebuff_);
- cell = parse_record(line);
- return is;
- }
+ std::istream &parse_record(std::istream &is, cell_record &cell);
private:
- void check_parse_status(const std::istream &is)
- {
- if (is.fail())
- // If we try to read past the eof; fail bit will also be set
- // FIXME: better throw a custom parse_error exception
- throw std::logic_error("could not parse value");
- }
-
- template<typename T>
- T parse_value_strict(std::istream &is)
- {
- T val;
- check_parse_status(is >> val);
-
- // everything's fine
- return val;
- }
-
// Read the record from a string stream; will be treated like a single line
cell_record parse_record(std::istringstream &is);
@@ -209,55 +189,16 @@ private:
};
-// specialize parsing for cell types
-template<>
-cell_record::kind swc_parser::parse_value_strict(std::istream &is)
-{
- int val;
- check_parse_status(is >> val);
-
- // Let cell_record's constructor check for the type validity
- return static_cast<cell_record::kind>(val);
-}
-
-cell_record swc_parser::parse_record(std::istringstream &is)
-{
- auto id = parse_value_strict<int>(is);
- auto type = parse_value_strict<cell_record::kind>(is);
- auto x = parse_value_strict<float>(is);
- auto y = parse_value_strict<float>(is);
- auto z = parse_value_strict<float>(is);
- auto r = parse_value_strict<float>(is);
- auto parent_id = parse_value_strict<int>(is);
-
- // Convert to zero-based, leaving parent_id as-is if -1
- if (parent_id != -1)
- parent_id--;
-
- return cell_record(type, id-1, x, y, z, r, parent_id);
-}
+std::istream &operator>>(std::istream &is, cell_record &cell);
+//
+// Reads cells from an input stream until an eof is encountered and returns a
+// cleaned sequence of cell records.
+//
+// For more information check here:
+// https://github.com/eth-cscs/cell_algorithms/wiki/SWC-file-parsing
+//
+std::vector<cell_record> swc_read_cells(std::istream &is);
-std::istream &operator>>(std::istream &is, cell_record &cell)
-{
- swc_parser parser;
- parser.parse_record(is, cell);
- return is;
-}
-
-std::ostream &operator<<(std::ostream &os, const cell_record &cell)
-{
- // output in one-based indexing
- os << cell.id_+1 << " "
- << cell.type_ << " "
- << std::setprecision(7) << cell.x_ << " "
- << std::setprecision(7) << cell.y_ << " "
- << std::setprecision(7) << cell.z_ << " "
- << std::setprecision(7) << cell.r_ << " "
- << ((cell.parent_id_ == -1) ? cell.parent_id_ : cell.parent_id_+1);
-
- return os;
-}
-
-} // end of neuron::io
-} // end of neuron
+} // end of nestmc::io
+} // end of nestmc
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index ba8b481ad51109f2a758ebd19c4281d5f3792d2f..7e787ba46b50a48d29ca3d5dd8d70f8116354d5d 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -18,10 +18,9 @@ set(TEST_SOURCES
add_executable(test.exe ${TEST_SOURCES} ${HEADERS})
-#target_link_libraries(test_compiler LINK_PUBLIC compiler gtest)
-
-#set_target_properties( test.exe
-# PROPERTIES
-# RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests"
-#)
+target_link_libraries(test.exe LINK_PUBLIC cellalgo)
+set_target_properties(test.exe
+ PROPERTIES
+ RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests"
+)
diff --git a/tests/test_swcio.cpp b/tests/test_swcio.cpp
index 5a11a5ed214e884d0e5da248ba2fff2abcaac4b2..a3a7336a8716877cb2a6bc084b4b866a9b58a005 100644
--- a/tests/test_swcio.cpp
+++ b/tests/test_swcio.cpp
@@ -1,3 +1,4 @@
+#include <array>
#include <iostream>
#include <fstream>
#include <numeric>
@@ -8,8 +9,8 @@
#include <swcio.hpp>
// SWC tests
-void expect_cell_equals(const neuron::io::cell_record &expected,
- const neuron::io::cell_record &actual)
+void expect_cell_equals(const nestmc::io::cell_record &expected,
+ const nestmc::io::cell_record &actual)
{
EXPECT_EQ(expected.id(), actual.id());
EXPECT_EQ(expected.type(), actual.type());
@@ -22,7 +23,7 @@ void expect_cell_equals(const neuron::io::cell_record &expected,
TEST(cell_record, construction)
{
- using namespace neuron::io;
+ using namespace nestmc::io;
{
// force an invalid type
@@ -87,27 +88,45 @@ TEST(cell_record, construction)
}
}
+TEST(cell_record, comparison)
+{
+ using namespace nestmc::io;
+
+ {
+ // check comparison operators
+ cell_record cell0(cell_record::custom, 0, 1., 1., 1., 1., -1);
+ cell_record cell1(cell_record::custom, 0, 2., 3., 4., 5., -1);
+ cell_record cell2(cell_record::custom, 1, 2., 3., 4., 5., -1);
+ EXPECT_EQ(cell0, cell1);
+ EXPECT_LT(cell0, cell2);
+ EXPECT_GT(cell2, cell1);
+ }
+
+}
+
TEST(swc_parser, invalid_input)
{
- using namespace neuron::io;
+ using namespace nestmc::io;
{
// check incomplete lines; missing parent
std::istringstream is("1 1 14.566132 34.873772 7.857000 0.717830\n");
cell_record cell;
- EXPECT_THROW(is >> cell, std::logic_error);
+ EXPECT_THROW(is >> cell, swc_parse_error);
}
{
// Check non-parsable values
- std::istringstream is("1a 1 14.566132 34.873772 7.857000 0.717830 -1\n");
+ std::istringstream is(
+ "1a 1 14.566132 34.873772 7.857000 0.717830 -1\n");
cell_record cell;
- EXPECT_THROW(is >> cell, std::logic_error);
+ EXPECT_THROW(is >> cell, swc_parse_error);
}
{
// Check invalid cell type
- std::istringstream is("1 10 14.566132 34.873772 7.857000 0.717830 -1\n");
+ std::istringstream is(
+ "1 10 14.566132 34.873772 7.857000 0.717830 -1\n");
cell_record cell;
EXPECT_THROW(is >> cell, std::invalid_argument);
}
@@ -116,7 +135,7 @@ TEST(swc_parser, invalid_input)
TEST(swc_parser, valid_input)
{
- using namespace neuron::io;
+ using namespace nestmc::io;
{
// check empty file; no cell may be parsed
@@ -175,14 +194,14 @@ TEST(swc_parser, valid_input)
++nr_records;
}
} catch (std::exception &e) {
- ADD_FAILURE();
+ ADD_FAILURE() << "unexpected exception thrown\n";
}
}
}
TEST(swc_parser, from_allen_db)
{
- using namespace neuron;
+ using namespace nestmc;
auto fname = "../data/example.swc";
std::ifstream fid(fname);
@@ -197,6 +216,83 @@ TEST(swc_parser, from_allen_db)
while( !(fid >> node).eof()) {
nodes.push_back(std::move(node));
}
+
// verify that the correct number of nodes was read
EXPECT_EQ(nodes.size(), 1058u);
}
+
+TEST(swc_parser, input_cleaning)
+{
+ using namespace nestmc::io;
+
+ {
+ // 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";
+
+ auto cells = swc_read_cells(is);
+ EXPECT_EQ(2, cells.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 << "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";
+
+ auto cells = swc_read_cells(is);
+ EXPECT_EQ(2, cells.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 << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
+
+ std::array<cell_record::id_type, 4> expected_id_list = {{ 0, 1, 2, 3 }};
+ auto cells = swc_read_cells(is);
+ ASSERT_EQ(4, cells.size());
+
+ auto expected_id = expected_id_list.cbegin();
+ for (const auto &c : cells) {
+ EXPECT_EQ(*expected_id, c.id());
+ ++expected_id;
+ }
+ }
+
+ {
+ // 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";
+
+ auto cells = swc_read_cells(is);
+ std::array<cell_record::id_type, 6> expected_id_list =
+ {{ 0, 1, 2, 3, 4, 5 }};
+ std::array<cell_record::id_type, 6> expected_parent_list =
+ {{ -1, 0, 1, 1, 0, 4 }};
+ ASSERT_EQ(6, cells.size());
+
+ auto expected_id = expected_id_list.cbegin();
+ auto expected_parent = expected_parent_list.cbegin();
+ for (const auto &c : cells) {
+ EXPECT_EQ(*expected_id, c.id());
+ EXPECT_EQ(*expected_parent, c.parent());
+ ++expected_id;
+ ++expected_parent;
+ }
+
+ }
+}
diff --git a/vector b/vector
index a8dfadd460262ebbc1bc22b159efe9e33ad1d359..9c86d0a84efed0dd739888503d275378df67fe71 160000
--- a/vector
+++ b/vector
@@ -1 +1 @@
-Subproject commit a8dfadd460262ebbc1bc22b159efe9e33ad1d359
+Subproject commit 9c86d0a84efed0dd739888503d275378df67fe71