Cell record cleaning.

Read in cell records and do a basic cleaning:
  1. Remove duplicates.
     Only the first occurance is considered.
  2. Sort the records if not sorted.
  3. Ignore subsequent trees if more than one are specified.
  4. Renumber records if numbering is not contiguous.

src/swcio.hpp

@@ -2,8 +2,11 @@
 
 #include <exception>
 #include <iostream>
+#include <map>
 #include <sstream>
 #include <type_traits>
+#include <unordered_set>
+#include <vector>
 
 namespace nestmc
 {
@@ -20,9 +23,10 @@ static bool starts_with(const std::string &str, const std::string &prefix)
 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 +52,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()
@@ -123,12 +112,12 @@ public:
         return type_;
     }
 
-    int id() const
+    id_type id() const
     {
         return id_;
     }
 
-    int parent() const
+    id_type parent() const
     {
         return parent_id_;
     }
@@ -158,12 +147,47 @@ public:
         return 2*r_;
     }
 
+    void 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));
+    }
+
 private:
+    void 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");
+    }
+
     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
@@ -249,7 +273,7 @@ private:
 template<>
 cell_record::kind swc_parser::parse_value_strict(std::istream &is)
 {
-    int val;
+    cell_record::id_type val;
     check_parse_status(is >> val);
 
     // Let cell_record's constructor check for the type validity
@@ -296,5 +320,54 @@ std::ostream &operator<<(std::ostream &os, const cell_record &cell)
     return os;
 }
 
+//
+// 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::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

tests/test_swcio.cpp

+#include <array>
 #include <iostream>
 #include <fstream>
 #include <numeric>
@@ -191,7 +192,7 @@ TEST(swc_parser, valid_input)
                 ++nr_records;
             }
         } catch (std::exception &e) {
-            ADD_FAILURE();
+            ADD_FAILURE() << "unexpected exception thrown\n";
         }
     }
 }
@@ -216,3 +217,79 @@ TEST(swc_parser, from_allen_db)
     // 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;
+        }
+        
+    }
+}