Merge branch 'master' of github.com:eth-cscs/cell_algorithms

main.cpp

 #include <iostream>
 #include <fstream>
 #include <numeric>
+#include <vector>
 
 #include "gtest/gtest.h"
 
 #include "cell_tree.hpp"
+#include "swcio.hpp"
 #include "json/src/json.hpp"
 
 using json = nlohmann::json;
@@ -250,6 +252,179 @@ TEST(cell_tree, json_load) {
     }
 }
 
+// SWC tests
+void expect_cell_equals(const neuron::io::cell_record &expected,
+                        const neuron::io::cell_record &actual)
+{
+    EXPECT_EQ(expected.id(), actual.id());
+    EXPECT_EQ(expected.type(), actual.type());
+    EXPECT_FLOAT_EQ(expected.x(), actual.x());
+    EXPECT_FLOAT_EQ(expected.y(), actual.y());
+    EXPECT_FLOAT_EQ(expected.z(), actual.z());
+    EXPECT_FLOAT_EQ(expected.radius(), actual.radius());
+    EXPECT_EQ(expected.parent(), actual.parent());
+}
+
+TEST(cell_record, construction)
+{
+    using namespace neuron::io;
+
+    {
+        // force an invalid type
+        cell_record::kind invalid_type = static_cast<cell_record::kind>(100);
+        EXPECT_THROW(cell_record cell(invalid_type, 7, 1., 1., 1., 1., 5),
+                     std::invalid_argument);
+    }
+
+    {
+        // invalid id
+        EXPECT_THROW(cell_record cell(
+                         cell_record::custom, -3, 1., 1., 1., 1., 5),
+                     std::invalid_argument);
+    }
+
+    {
+        // invalid parent id
+        EXPECT_THROW(cell_record cell(
+                         cell_record::custom, 0, 1., 1., 1., 1., -5),
+                     std::invalid_argument);
+    }
+
+    {
+        // invalid radius
+        EXPECT_THROW(cell_record cell(
+                         cell_record::custom, 0, 1., 1., 1., -1., -1),
+                     std::invalid_argument);
+    }
+
+    {
+        // parent_id > id
+        EXPECT_THROW(cell_record cell(
+                         cell_record::custom, 0, 1., 1., 1., 1., 2),
+                     std::invalid_argument);
+    }
+
+    {
+        // parent_id == id
+        EXPECT_THROW(cell_record cell(
+                         cell_record::custom, 0, 1., 1., 1., 1., 0),
+                     std::invalid_argument);
+    }
+
+    {
+        // check standard construction by value
+        cell_record cell(cell_record::custom, 0, 1., 1., 1., 1., -1);
+        EXPECT_EQ(cell.id(), 0);
+        EXPECT_EQ(cell.type(), cell_record::custom);
+        EXPECT_EQ(cell.x(), 1.);
+        EXPECT_EQ(cell.y(), 1.);
+        EXPECT_EQ(cell.z(), 1.);
+        EXPECT_EQ(cell.radius(), 1.);
+        EXPECT_EQ(cell.diameter(), 2*1.);
+        EXPECT_EQ(cell.parent(), -1);
+    }
+
+    {
+        // check copy constructor
+        cell_record cell_orig(cell_record::custom, 0, 1., 1., 1., 1., -1);
+        cell_record cell(cell_orig);
+        expect_cell_equals(cell_orig, cell);
+    }
+}
+
+TEST(swc_parser, invalid_input)
+{
+    using namespace neuron::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);
+    }
+
+    {
+        // Check non-parsable values
+        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);
+    }
+
+    {
+        // Check invalid cell type
+        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);
+    }
+}
+
+
+TEST(swc_parser, valid_input)
+{
+    using namespace neuron::io;
+
+    {
+        // check empty file; no cell may be parsed
+        cell_record cell, cell_orig;
+        std::istringstream is("");
+        EXPECT_NO_THROW(is >> cell);
+        expect_cell_equals(cell_orig, cell);
+    }
+
+    {
+        // check comment-only file not ending with a newline;
+        // no cell may be parsed
+        cell_record cell, cell_orig;
+        std::istringstream is("#comment\n#comment");
+        EXPECT_NO_THROW(is >> cell);
+        expect_cell_equals(cell_orig, cell);
+    }
+
+
+    {
+        // check last line case (no newline at the end)
+        std::istringstream is("1 1 14.566132 34.873772 7.857000 0.717830 -1");
+        cell_record cell;
+        EXPECT_NO_THROW(is >> cell);
+        EXPECT_EQ(0, cell.id());    // zero-based indexing
+        EXPECT_EQ(cell_record::soma, cell.type());
+        EXPECT_FLOAT_EQ(14.566132, cell.x());
+        EXPECT_FLOAT_EQ(34.873772, cell.y());
+        EXPECT_FLOAT_EQ( 7.857000, cell.z());
+        EXPECT_FLOAT_EQ( 0.717830, cell.radius());
+        EXPECT_FLOAT_EQ( -1, cell.parent());
+    }
+
+    {
+        // check valid input with a series of records
+        std::vector<cell_record> cells_orig = {
+            cell_record(cell_record::soma, 0,
+                        14.566132, 34.873772, 7.857000, 0.717830, -1),
+            cell_record(cell_record::dendrite, 1,
+                        14.566132+1, 34.873772+1, 7.857000+1, 0.717830+1, -1)
+        };
+
+        std::stringstream swc_input;
+        swc_input << "# this is a comment\n";
+        swc_input << "# this is a comment\n";
+        for (auto c : cells_orig)
+            swc_input << c << "\n";
+
+        swc_input << "# this is a final comment\n";
+        try {
+            std::size_t nr_records = 0;
+            cell_record cell;
+            while ( !(swc_input >> cell).eof()) {
+                ASSERT_LT(nr_records, cells_orig.size());
+                expect_cell_equals(cells_orig[nr_records], cell);
+                ++nr_records;
+            }
+        } catch (std::exception &e) {
+            ADD_FAILURE();
+        }
+    }
+}
+
 int main(int argc, char **argv) {
     ::testing::InitGoogleTest(&argc, argv);
     return RUN_ALL_TESTS();

makefile

 CC=clang++
-FLAGS=-std=c++11 -g
+FLAGS=-std=c++11 -g -pedantic
 
 test.exe : main.cpp *.hpp makefile gtest.o
 	$(CC) $(FLAGS) main.cpp -o test.exe gtest.o -pthread

swcio.hpp

+#pragma once
+
+#include <exception>
+#include <iostream>
+#include <sstream>
+#include <type_traits>
+
+namespace neuron
+{
+
+namespace io
+{
+
+
+static bool starts_with(const std::string &str, const std::string &prefix)
+{
+    return (str.find(prefix) == 0);
+}
+
+class cell_record 
+{
+public:
+
+    // 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 {
+        undefined = 0,
+        soma,
+        axon,
+        dendrite,
+        apical_dendrite,
+        fork_point,
+        end_point,
+        custom
+    };
+
+    // cell records assume zero-based indexing; root's parent remains -1
+    cell_record(kind type, int id, 
+                float x, float y, float z, float r,
+                int parent_id)
+        : type_(type)
+        , id_(id)
+        , x_(x)
+        , y_(y)
+        , z_(z)
+        , 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");
+    }
+    
+    cell_record()
+        : type_(cell_record::undefined)
+        , id_(0)
+        , x_(0)
+        , y_(0)
+        , z_(0)
+        , r_(0)
+        , parent_id_(-1)
+    { }
+
+    cell_record(const cell_record &other) = default;
+    cell_record &operator=(const cell_record &other) = default;
+
+    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);
+    }
+
+    friend std::ostream &operator<<(std::ostream &os, const cell_record &cell);
+
+    kind type() const
+    {
+        return type_;
+    }
+
+    int id() const
+    {
+        return id_;
+    }
+
+    int parent() const
+    {
+        return parent_id_;
+    }
+
+    float x() const
+    {
+        return x_;
+    }
+
+    float y() const
+    {
+        return y_;
+    }
+
+    float z() const
+    {
+        return z_;
+    }
+
+    float radius() const
+    {
+        return r_;
+    }
+
+    float diameter() const
+    {
+        return 2*r_;
+    }
+
+private:
+    kind type_;         // cell type
+    int id_;            // cell id
+    float x_, y_, z_;   // cell coordinates
+    float r_;           // cell radius
+    int parent_id_;     // cell parent's id
+};
+
+class swc_parser
+{
+public:
+    swc_parser(const std::string &delim,
+               std::string comment_prefix)
+        : delim_(delim)
+        , comment_prefix_(comment_prefix)
+    { }
+
+    swc_parser()
+        : delim_(" ")
+        , 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;
+    }
+
+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);
+
+    std::string delim_;
+    std::string comment_prefix_;
+    std::string linebuff_;
+};
+
+
+// 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)
+{
+    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