Merge pull request #6 from eth-cscs/swc-io

Cell record cleaning + minor other improvements

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)
-

src/CMakeLists.txt

+set(HEADERS
+    swcio.hpp
+    )
+set(BASE_SOURCES
+    swcio.cpp
+)
+
+add_library(cellalgo ${BASE_SOURCES} ${HEADERS})

src/swcio.cpp

+#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

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

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"
+)

tests/test_swcio.cpp

+#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;
+        }
+        
+    }
+}

vector @ 9c86d0a8

-Subproject commit a8dfadd460262ebbc1bc22b159efe9e33ad1d359
+Subproject commit 9c86d0a84efed0dd739888503d275378df67fe71