diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9bedee3de0d7a52497259c182e34dd1a2b16be1a..70135e7d77117e77113b0fd99f172cf57c7c0fdd 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -154,7 +154,7 @@ set(CMAKE_CXX_EXTENSIONS OFF)
# in the same CMakeLists.txt in which the target is defined.
# Interface library `arbor-config-defs` collects configure-time defines
-# for arbor, arborenv and arbornml, of the form ARB_HAVE_XXX. These
+# for arbor, arborenv, arborio and arbornml, of the form ARB_HAVE_XXX. These
# defines should _not_ be used in any installed public headers.
add_library(arbor-config-defs INTERFACE)
@@ -181,6 +181,13 @@ add_library(arbornml-private-deps INTERFACE)
target_link_libraries(arbornml-private-deps INTERFACE arbor-config-defs)
install(TARGETS arbornml-private-deps EXPORT arbor-targets)
+# Interface library `arborio-private-deps` collects dependencies, options etc.
+# for the arborio library.
+
+add_library(arborio-private-deps INTERFACE)
+target_link_libraries(arborio-private-deps INTERFACE arbor-config-defs)
+install(TARGETS arborio-private-deps EXPORT arbor-targets)
+
# Interface library `arbor-public-deps` collects requirements for the
# users of the arbor library (e.g. mpi) that will become part
# of arbor's PUBLIC interface.
@@ -404,6 +411,9 @@ add_subdirectory(arbor)
# arborenv, arborenv-public-headers:
add_subdirectory(arborenv)
+# arborio, arborio-public-headers:
+add_subdirectory(arborio)
+
# arbornml, arbornml-public-headers:
if(ARB_WITH_NEUROML)
add_subdirectory(arbornml)
@@ -459,6 +469,7 @@ set(arbor_override_import_lang)
set(arbor_add_import_libs)
set(arborenv_add_import_libs)
set(arbornml_add_import_libs)
+set(arborio_add_import_libs)
if(ARB_WITH_GPU)
set(arbor_override_import_lang CXX)
diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt
index 0368fc88613d5213bdc9582a473ad4ae6dbcf39f..5f135d92ab3d6816b2a02777fbbe8f656f851356 100644
--- a/arbor/CMakeLists.txt
+++ b/arbor/CMakeLists.txt
@@ -53,7 +53,6 @@ set(arbor_sources
schedule.cpp
spike_event_io.cpp
spike_source_cell_group.cpp
- swcio.cpp
s_expr.cpp
threading/threading.cpp
thread_private_spike_store.cpp
diff --git a/arbor/swcio.cpp b/arbor/swcio.cpp
deleted file mode 100644
index c17b375ecf14e226f81083a5308be2fc5ee75558..0000000000000000000000000000000000000000
--- a/arbor/swcio.cpp
+++ /dev/null
@@ -1,193 +0,0 @@
-#include <iostream>
-#include <limits>
-#include <string>
-#include <sstream>
-#include <unordered_map>
-#include <unordered_set>
-#include <vector>
-
-#include <arbor/morph/segment_tree.hpp>
-#include <arbor/swcio.hpp>
-
-#include "io/save_ios.hpp"
-#include "util/rangeutil.hpp"
-
-namespace arb {
-
-// SWC exceptions:
-
-swc_error::swc_error(const std::string& msg, int record_id):
- arbor_exception(msg+": sample id "+std::to_string(record_id)),
- record_id(record_id)
-{}
-
-swc_no_such_parent::swc_no_such_parent(int record_id):
- swc_error("missing SWC parent record", record_id)
-{}
-
-swc_record_precedes_parent::swc_record_precedes_parent(int record_id):
- swc_error("SWC parent id is not less than sample id", record_id)
-{}
-
-swc_duplicate_record_id::swc_duplicate_record_id(int record_id):
- swc_error("duplicate SWC sample id", record_id)
-{}
-
-swc_spherical_soma::swc_spherical_soma(int record_id):
- swc_error("SWC with spherical somata are not supported", record_id)
-{}
-
-bad_swc_data::bad_swc_data(int record_id):
- swc_error("Cannot interpret bad SWC data", record_id)
-{}
-
-// Record I/O
-
-std::ostream& operator<<(std::ostream& out, const swc_record& record) {
- io::save_ios_flags save(out);
-
- out.precision(std::numeric_limits<double>::digits10+2);
- return out << record.id << ' ' << record.tag << ' '
- << record.x << ' ' << record.y << ' ' << record.z << ' ' << record.r << ' '
- << record.parent_id << '\n';
-}
-
-std::istream& operator>>(std::istream& in, swc_record& record) {
- std::string line;
- if (!getline(in, line, '\n')) return in;
-
- swc_record r;
- std::istringstream s(line);
- s >> r.id >> r.tag >> r.x >> r.y >> r.z >> r.r >> r.parent_id;
- if (s) {
- record = r;
- }
- else {
- in.setstate(std::ios_base::failbit);
- }
-
- return in;
-}
-
-// Parse SWC format data (comments and sequence of SWC records).
-
-static std::vector<swc_record> sort_and_validate_swc(std::vector<swc_record> records, swc_mode mode) {
- if (records.empty()) return {};
-
- std::unordered_set<int> seen;
- std::size_t n_rec = records.size();
- int first_id = records[0].id;
- int first_tag = records[0].tag;
-
- if (records.size()<2) {
- throw swc_spherical_soma(first_id);
- }
-
- for (std::size_t i = 0; i<n_rec; ++i) {
- swc_record& r = records[i];
-
- if (r.parent_id>=r.id) {
- throw swc_record_precedes_parent(r.id);
- }
-
- if (!seen.insert(r.id).second) {
- throw swc_duplicate_record_id(r.id);
- }
- }
-
- util::sort_by(records, [](auto& r) { return r.id; });
- bool first_tag_match = false;
-
- for (std::size_t i = 0; i<n_rec; ++i) {
- const swc_record& r = records[i];
- first_tag_match |= r.parent_id==first_id && r.tag==first_tag;
-
- if ((i==0 && r.parent_id!=-1) || (i>0 && !seen.count(r.parent_id))) {
- throw swc_no_such_parent(r.id);
- }
- }
-
- if (mode==swc_mode::strict && !first_tag_match) {
- throw swc_spherical_soma(first_id);
- }
-
- return records;
-}
-
-swc_data parse_swc(std::istream& in, swc_mode mode) {
- // Collect any initial comments (lines beginning with '#').
-
- swc_data data;
- std::string line;
-
- while (in) {
- auto c = in.get();
- if (c=='#') {
- getline(in, line, '\n');
- auto from = line.find_first_not_of(" \t");
- if (from != std::string::npos) {
- data.metadata.append(line, from);
- }
- data.metadata += '\n';
- }
- else {
- in.unget();
- break;
- }
- }
-
- swc_record r;
- while (in && in >> r) {
- data.records.push_back(r);
- }
-
- data.records = sort_and_validate_swc(std::move(data.records), mode);
- return data;
-}
-
-swc_data parse_swc(const std::string& text, swc_mode mode) {
- std::istringstream is(text);
- return parse_swc(is, mode);
-}
-
-swc_data parse_swc(std::vector<swc_record> records, swc_mode mode) {
- swc_data data;
- data.records = sort_and_validate_swc(std::move(records), mode);
- return data;
-}
-
-segment_tree as_segment_tree(const std::vector<swc_record>& records) {
- if (records.empty()) return {};
- if (records.size()<2) throw bad_swc_data{records.front().id};
-
- segment_tree tree;
- std::size_t n_seg = records.size()-1;
- tree.reserve(n_seg);
-
- std::unordered_map<int, std::size_t> id_to_index;
- id_to_index[records[0].id] = 0;
-
- // ith segment is built from i+1th SWC record and its parent.
- for (std::size_t i = 1; i<n_seg+1; ++i) {
- const auto& dist = records[i];
-
- auto iter = id_to_index.find(dist.parent_id);
- if (iter==id_to_index.end()) throw bad_swc_data{dist.id};
- auto parent_idx = iter->second;
-
- const auto& prox = records[parent_idx];
- msize_t seg_parent = parent_idx? parent_idx-1: mnpos;
-
- tree.append(seg_parent,
- mpoint{prox.x, prox.y, prox.z, prox.r},
- mpoint{dist.x, dist.y, dist.z, dist.r},
- dist.tag);
-
- id_to_index[dist.id] = i;
- }
-
- return tree;
-}
-
-} // namespace arb
-
diff --git a/arborio/CMakeLists.txt b/arborio/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e07ce90d4ad50dcedece2352043487c5aab618d6
--- /dev/null
+++ b/arborio/CMakeLists.txt
@@ -0,0 +1,21 @@
+set(arborio-sources
+ swcio.cpp
+)
+
+add_library(arborio ${arborio-sources})
+
+add_library(arborio-public-headers INTERFACE)
+target_include_directories(arborio-public-headers INTERFACE
+ $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
+ $<INSTALL_INTERFACE:include>
+)
+
+target_link_libraries(arborio PUBLIC arbor arborio-public-headers)
+target_link_libraries(arborio PRIVATE arbor-config-defs arborio-private-deps)
+
+install(DIRECTORY include/arborio
+ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
+ FILES_MATCHING PATTERN "*.hpp")
+
+install(TARGETS arborio-public-headers EXPORT arbor-targets)
+install(TARGETS arborio EXPORT arbor-targets ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR})
diff --git a/arbor/include/arbor/swcio.hpp b/arborio/include/arborio/swcio.hpp
similarity index 60%
rename from arbor/include/arbor/swcio.hpp
rename to arborio/include/arborio/swcio.hpp
index db5b672f2cfdf39f576699ba39f734d4b9db9944..9fff6d9020d895cf60bebe6698380b23c1654cb9 100644
--- a/arbor/include/arbor/swcio.hpp
+++ b/arborio/include/arborio/swcio.hpp
@@ -7,12 +7,12 @@
#include <arbor/arbexcept.hpp>
#include <arbor/morph/segment_tree.hpp>
-namespace arb {
+namespace arborio {
// SWC exceptions are thrown by `parse_swc`, and correspond
// to inconsistent, or in `strict` mode, dubious SWC data.
-struct swc_error: public arbor_exception {
+struct swc_error: public arb::arbor_exception {
swc_error(const std::string& msg, int record_id);
int record_id;
};
@@ -38,9 +38,59 @@ struct swc_spherical_soma: swc_error {
explicit swc_spherical_soma(int record_id);
};
-// Bad or inconsistent SWC data was fed to an `swc_data` consumer.
-struct bad_swc_data: swc_error {
- explicit bad_swc_data(int record_id);
+// Smells like a non-spherical soma.
+struct swc_non_spherical_soma: swc_error {
+ explicit swc_non_spherical_soma(int record_id);
+};
+
+// Missing soma.
+struct swc_no_soma: swc_error {
+ explicit swc_no_soma(int record_id);
+};
+
+// Non-consecutive soma samples.
+struct swc_non_consecutive_soma: swc_error {
+ explicit swc_non_consecutive_soma(int record_id);
+};
+
+// Non-serial soma samples.
+struct swc_non_serial_soma: swc_error {
+ explicit swc_non_serial_soma(int record_id);
+};
+
+// Sample connecting to the middle of a soma causing an unsupported branch.
+struct swc_branchy_soma: swc_error {
+ explicit swc_branchy_soma(int record_id);
+};
+
+// Sample connecting to the middle of a soma causing an unsupported branch.
+struct swc_collocated_soma: swc_error {
+ explicit swc_collocated_soma(int record_id);
+};
+
+// Sample is not part of a segment
+struct swc_single_sample_segment: swc_error {
+ explicit swc_single_sample_segment(int record_id);
+};
+
+// Segment cannot have samples with different tags
+struct swc_mismatched_tags: swc_error {
+ explicit swc_mismatched_tags(int record_id);
+};
+
+// Only tags 1, 2, 3, 4 supported
+struct swc_unsupported_tag: swc_error {
+ explicit swc_unsupported_tag(int record_id);
+};
+
+// No gaps allowed
+struct swc_unsupported_gaps: swc_error {
+ explicit swc_unsupported_gaps(int record_id);
+};
+
+// Can't form a segment from a single sample
+struct swc_bad_description: swc_error {
+ explicit swc_bad_description(int record_id);
};
struct swc_record {
@@ -119,11 +169,31 @@ swc_data parse_swc(std::vector<swc_record>, swc_mode = swc_mode::strict);
// each SWC record after the first: this record defines the tag and distal point
// of the segment, while the proximal point is taken from the parent record.
-segment_tree as_segment_tree(const std::vector<swc_record>&);
+arb::segment_tree as_segment_tree(const std::vector<swc_record>&);
-inline segment_tree as_segment_tree(const swc_data& data) {
+inline arb::segment_tree as_segment_tree(const swc_data& data) {
return as_segment_tree(data.records);
}
+// As above, will convert a valid, ordered sequence of SWC records to a morphological
+// segment tree.
+//
+// Note that 'one-point soma' SWC files are supported here; the swc_data is expected
+// to abide by the restrictions of `relaxed` mode parsing as described above.
+//
+// These functions comply with inferred SWC rules from the Allen institute and Neuron.
+// These rules are explicitly listed in the docs.
+
+arb::segment_tree load_swc_neuron(const std::vector<swc_record>& records);
+
+inline arb::segment_tree load_swc_neuron(const swc_data& data) {
+ return load_swc_neuron(data.records);
+}
+
+arb::segment_tree load_swc_allen(std::vector<swc_record>& records, bool no_gaps=false);
+
+inline arb::segment_tree load_swc_allen(swc_data& data) {
+ return load_swc_allen(data.records);
+}
-} // namespace arb
+} // namespace arborio
diff --git a/arborio/swcio.cpp b/arborio/swcio.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..984dc61ad520df5f2885766ae452637a3e4dd4c6
--- /dev/null
+++ b/arborio/swcio.cpp
@@ -0,0 +1,503 @@
+#include <cmath>
+#include <ios>
+#include <iostream>
+#include <limits>
+#include <numeric>
+#include <set>
+#include <string>
+#include <sstream>
+#include <unordered_map>
+#include <unordered_set>
+#include <vector>
+
+#include <arbor/morph/segment_tree.hpp>
+
+#include <arborio/swcio.hpp>
+
+namespace arborio {
+
+// SWC exceptions:
+
+swc_error::swc_error(const std::string& msg, int record_id):
+ arbor_exception(msg+": sample id "+std::to_string(record_id)),
+ record_id(record_id)
+{}
+
+swc_no_such_parent::swc_no_such_parent(int record_id):
+ swc_error("Missing SWC parent record", record_id)
+{}
+
+swc_record_precedes_parent::swc_record_precedes_parent(int record_id):
+ swc_error("SWC parent id is not less than sample id", record_id)
+{}
+
+swc_duplicate_record_id::swc_duplicate_record_id(int record_id):
+ swc_error("duplicate SWC sample id", record_id)
+{}
+
+swc_spherical_soma::swc_spherical_soma(int record_id):
+ swc_error("SWC with spherical somata are not supported", record_id)
+{}
+
+swc_non_spherical_soma::swc_non_spherical_soma(int record_id):
+ swc_error("SWC with multi-sample somata are not supported", record_id)
+{}
+
+swc_no_soma::swc_no_soma(int record_id):
+ swc_error("No soma (tag 1) found at the root", record_id)
+{}
+
+swc_non_consecutive_soma::swc_non_consecutive_soma (int record_id):
+ swc_error("Soma samples (tag 1) are not listed consecutively", record_id)
+{}
+
+swc_non_serial_soma::swc_non_serial_soma (int record_id):
+ swc_error("Soma samples (tag 1) are not listed serially", record_id)
+{}
+
+swc_branchy_soma::swc_branchy_soma (int record_id):
+ swc_error("Non-soma sample (tag >= 1) connected to a non-distal sample of the soma", record_id)
+{}
+
+swc_collocated_soma::swc_collocated_soma(int record_id):
+ swc_error("The samples that make the soma (tag 1) are not allowed to be collocated", record_id)
+{}
+
+swc_single_sample_segment::swc_single_sample_segment(int record_id):
+ swc_error("Segments connected to the soma (tag 1) must have 2 samples with the same tag", record_id)
+{}
+
+swc_mismatched_tags::swc_mismatched_tags(int record_id):
+ swc_error("Every record not attached to the soma (tag 1) must have the same tag as its parent", record_id)
+{}
+
+swc_unsupported_tag::swc_unsupported_tag(int record_id):
+ swc_error("Every record must have a tag of 2, 3 or 4, except for the first which must have tag 1", record_id)
+{}
+
+swc_unsupported_gaps::swc_unsupported_gaps(int record_id):
+ swc_error("No gaps are allowed between the soma and any axons, dendrites or apical dendrites", record_id)
+{}
+
+swc_bad_description::swc_bad_description(int record_id):
+ swc_error("Need at least 2 samples to form a segment", record_id)
+{}
+
+
+// Record I/O:
+
+std::ostream& operator<<(std::ostream& out, const swc_record& record) {
+ std::ios_base::fmtflags flags(out.flags());
+
+ out.precision(std::numeric_limits<double>::digits10+2);
+ out << record.id << ' ' << record.tag << ' '
+ << record.x << ' ' << record.y << ' ' << record.z << ' ' << record.r << ' '
+ << record.parent_id << '\n';
+
+ out.flags(flags);
+
+ return out;
+}
+
+std::istream& operator>>(std::istream& in, swc_record& record) {
+ std::string line;
+ if (!getline(in, line, '\n')) return in;
+
+ swc_record r;
+ std::istringstream s(line);
+ s >> r.id >> r.tag >> r.x >> r.y >> r.z >> r.r >> r.parent_id;
+ if (s) {
+ record = r;
+ }
+ else {
+ in.setstate(std::ios_base::failbit);
+ }
+
+ return in;
+}
+
+// Parse SWC format data (comments and sequence of SWC records).
+
+static std::vector<swc_record> sort_and_validate_swc(std::vector<swc_record> records, swc_mode mode) {
+ if (records.empty()) return {};
+
+ std::unordered_set<int> seen;
+ std::size_t n_rec = records.size();
+ int first_id = records[0].id;
+ int first_tag = records[0].tag;
+
+ if (records.size()<2) {
+ throw swc_spherical_soma(first_id);
+ }
+
+ for (std::size_t i = 0; i<n_rec; ++i) {
+ swc_record& r = records[i];
+
+ if (r.parent_id>=r.id) {
+ throw swc_record_precedes_parent(r.id);
+ }
+
+ if (!seen.insert(r.id).second) {
+ throw swc_duplicate_record_id(r.id);
+ }
+ }
+
+ std::sort(records.begin(), records.end(), [](const auto& lhs, const auto& rhs) { return lhs.id < rhs.id; });
+ bool first_tag_match = false;
+
+ for (std::size_t i = 0; i<n_rec; ++i) {
+ const swc_record& r = records[i];
+ first_tag_match |= r.parent_id==first_id && r.tag==first_tag;
+
+ if ((i==0 && r.parent_id!=-1) || (i>0 && !seen.count(r.parent_id))) {
+ throw swc_no_such_parent(r.id);
+ }
+ }
+
+ if (mode==swc_mode::strict && !first_tag_match) {
+ throw swc_spherical_soma(first_id);
+ }
+
+ return records;
+}
+
+swc_data parse_swc(std::istream& in, swc_mode mode) {
+ // Collect any initial comments (lines beginning with '#').
+
+ swc_data data;
+ std::string line;
+
+ while (in) {
+ auto c = in.get();
+ if (c=='#') {
+ getline(in, line, '\n');
+ auto from = line.find_first_not_of(" \t");
+ if (from != std::string::npos) {
+ data.metadata.append(line, from);
+ }
+ data.metadata += '\n';
+ }
+ else {
+ in.unget();
+ break;
+ }
+ }
+
+ swc_record r;
+ while (in && in >> r) {
+ data.records.push_back(r);
+ }
+
+ data.records = sort_and_validate_swc(std::move(data.records), mode);
+ return data;
+}
+
+swc_data parse_swc(const std::string& text, swc_mode mode) {
+ std::istringstream is(text);
+ return parse_swc(is, mode);
+}
+
+swc_data parse_swc(std::vector<swc_record> records, swc_mode mode) {
+ swc_data data;
+ data.records = sort_and_validate_swc(std::move(records), mode);
+ return data;
+}
+
+arb::segment_tree as_segment_tree(const std::vector<swc_record>& records) {
+ if (records.empty()) return {};
+ if (records.size()<2) throw swc_bad_description{records.front().id};
+
+ arb::segment_tree tree;
+ std::size_t n_seg = records.size()-1;
+ tree.reserve(n_seg);
+
+ std::unordered_map<int, std::size_t> id_to_index;
+ id_to_index[records[0].id] = 0;
+
+ // ith segment is built from i+1th SWC record and its parent.
+ for (std::size_t i = 1; i<n_seg+1; ++i) {
+ const auto& dist = records[i];
+
+ auto iter = id_to_index.find(dist.parent_id);
+ if (iter==id_to_index.end()) throw swc_no_such_parent{dist.id};
+ auto parent_idx = iter->second;
+
+ const auto& prox = records[parent_idx];
+ arb::msize_t seg_parent = parent_idx? parent_idx-1: arb::mnpos;
+
+ tree.append(seg_parent,
+ arb::mpoint{prox.x, prox.y, prox.z, prox.r},
+ arb::mpoint{dist.x, dist.y, dist.z, dist.r},
+ dist.tag);
+
+ id_to_index[dist.id] = i;
+ }
+
+ return tree;
+}
+
+arb::segment_tree load_swc_neuron(const std::vector<swc_record>& records) {
+ if (records.empty()) return {};
+
+ const int soma_tag = 1;
+ auto soma_prox = records.front();
+
+ // Assert that root sample has tag 1.
+ if (soma_prox.tag != soma_tag) throw swc_no_soma{soma_prox.id};
+
+ // check for single soma cell
+ bool has_children = false;
+
+ // Map of SWC record id to index in `records`.
+ std::unordered_map<int, std::size_t> record_index;
+ record_index[soma_prox.id] = 0;
+
+ // Vector of records that make up the soma
+ std::vector<swc_record> soma_records = {soma_prox};
+ int prev_tag = soma_prox.tag;
+ int prev_id = soma_prox.id;
+
+ // Preliminary error checking and building the record_index
+ for (std::size_t i = 1; i < records.size(); ++i) {
+ const auto& r = records[i];
+ record_index[r.id] = i;
+
+ if (r.tag == soma_tag) {
+ if (prev_tag != soma_tag) throw swc_non_consecutive_soma{r.id};
+ if (prev_id != r.parent_id) throw swc_non_serial_soma{r.id};
+ soma_records.push_back(r);
+ }
+
+ // Find record index of the parent
+ auto parent_iter = record_index.find(r.parent_id);
+
+ if (parent_iter == record_index.end() || records[parent_iter->second].id == r.id) throw swc_no_such_parent{r.id};
+
+ if (r.tag != soma_tag && records[parent_iter->second].tag == soma_tag) {
+ if (r.parent_id != soma_records.back().id) throw swc_branchy_soma{r.id};
+ has_children = true;
+ }
+
+ prev_tag = r.tag;
+ prev_id = r.id;
+ }
+
+ arb::segment_tree tree;
+ tree.reserve(records.size());
+
+ // Map of SWC record id to index in `tree`.
+ std::unordered_map<int, arb::msize_t> tree_index;
+
+ // First, construct the soma
+ if (soma_records.size() == 1) {
+ if (!has_children) {
+ // Model the soma as a 1 cylinder with total length=2*radius, extended along the y axis
+ tree.append(arb::mnpos, {soma_prox.x, soma_prox.y - soma_prox.r, soma_prox.z, soma_prox.r},
+ {soma_prox.x, soma_prox.y + soma_prox.r, soma_prox.z, soma_prox.r}, soma_tag);
+ return tree;
+ }
+ // Model the soma as a 2 cylinders with total length=2*radius, extended along the y axis, centered at the sample
+ auto p = tree.append(arb::mnpos, {soma_prox.x, soma_prox.y - soma_prox.r, soma_prox.z, soma_prox.r},
+ {soma_prox.x, soma_prox.y, soma_prox.z, soma_prox.r}, soma_tag);
+ tree.append(p, {soma_prox.x, soma_prox.y, soma_prox.z, soma_prox.r},
+ {soma_prox.x, soma_prox.y + soma_prox.r, soma_prox.z, soma_prox.r}, soma_tag);
+ tree_index[soma_prox.id] = p;
+ }
+ else {
+ if (!has_children) {
+ // Don't split soma at the midpoint
+ arb::msize_t parent = arb::mnpos;
+ bool collocated_samples = true;
+ for (std::size_t i = 0; i < soma_records.size() - 1; ++i) {
+ const auto& p0 = soma_records[i];
+ const auto& p1 = soma_records[i + 1];
+ parent = tree.append(parent, {p0.x, p0.y, p0.z, p0.r}, {p1.x, p1.y, p1.z, p1.r}, 1);
+ collocated_samples &= ((p0.x == p1.x) && (p0.y == p1.y) && (p0.z == p1.z));
+ }
+ if (collocated_samples) {
+ throw swc_collocated_soma{records[0].id};
+ }
+ return tree;
+ }
+ // Calculate segment lengths
+ bool collocated_samples = true;
+ std::vector<double> soma_segment_lengths;
+ for (std::size_t i = 0; i < soma_records.size() - 1; ++i) {
+ const auto& p0 = soma_records[i];
+ const auto& p1 = soma_records[i + 1];
+ soma_segment_lengths.push_back(distance(arb::mpoint{p0.x, p0.y, p0.z, p0.r}, arb::mpoint{p1.x, p1.y, p1.z, p1.r}));
+ collocated_samples &= ((p0.x == p1.x) && (p0.y == p1.y) && (p0.z == p1.z));
+ }
+ if (collocated_samples) {
+ throw swc_collocated_soma{records[0].id};
+ }
+ double midlength = std::accumulate(soma_segment_lengths.begin(), soma_segment_lengths.end(), 0.) / 2;
+
+ std::size_t idx = 0;
+ for (; idx < soma_segment_lengths.size(); ++idx) {
+ auto l = soma_segment_lengths[idx];
+ if (midlength > l) {
+ midlength -= l;
+ continue;
+ }
+ break;
+ }
+
+ // Interpolate along the segment that contains the midpoint of the soma
+ double pos_on_segment = midlength / soma_segment_lengths[idx];
+
+ auto& r0 = soma_records[idx];
+ auto& r1 = soma_records[idx + 1];
+
+ auto x = r0.x + pos_on_segment * (r1.x - r0.x);
+ auto y = r0.y + pos_on_segment * (r1.y - r0.y);
+ auto z = r0.z + pos_on_segment * (r1.z - r0.z);
+ auto r = r0.r + pos_on_segment * (r1.r - r0.r);
+
+ arb::mpoint mid_soma = {x, y, z, r};
+
+ // Construct the soma
+ arb::msize_t parent = arb::mnpos;
+ for (std::size_t i = 0; i < idx; ++i) {
+ const auto& p0 = soma_records[i];
+ const auto& p1 = soma_records[i + 1];
+ parent = tree.append(parent, {p0.x, p0.y, p0.z, p0.r}, {p1.x, p1.y, p1.z, p1.r}, 1);
+ }
+ auto soma_seg = tree.append(parent, {r0.x, r0.y, r0.z, r0.r}, mid_soma, 1);
+
+ arb::mpoint r1_p{r1.x, r1.y, r1.z, r1.r};
+ parent = mid_soma != r1_p ? tree.append(soma_seg, mid_soma, r1_p, 1) : soma_seg;
+
+ for (std::size_t i = idx + 1; i < soma_records.size() - 1; ++i) {
+ const auto& p0 = soma_records[i];
+ const auto& p1 = soma_records[i + 1];
+ parent = tree.append(parent, {p0.x, p0.y, p0.z, p0.r}, {p1.x, p1.y, p1.z, p1.r}, 1);
+ }
+
+ tree_index[soma_records.back().id] = soma_seg;
+ }
+
+ // Build branches off soma.
+ std::set<int> unused_samples; // Samples that are not part of a segment
+ for (const auto& r: records) {
+ // Skip the soma samples
+ if (r.tag == soma_tag) continue;
+
+ const auto p = r.parent_id;
+
+ // Find parent segment of the record
+ auto pseg_iter = tree_index.find(p);
+ if (pseg_iter == tree_index.end()) throw swc_no_such_parent{r.id};
+
+ // Find parent record of the record
+ auto prec_iter = record_index.find(p);
+ if (prec_iter == record_index.end() || records[prec_iter->second].id == r.id) throw swc_no_such_parent{r.id};
+
+ // If the sample has a soma sample as its parent don't create a segment.
+ if (records[prec_iter->second].tag == soma_tag) {
+ // Map the sample id to the segment id of the soma (parent)
+ tree_index[r.id] = pseg_iter->second;
+ unused_samples.insert(r.id);
+ continue;
+ }
+
+ const auto& prox = records[prec_iter->second];
+ tree_index[r.id] = tree.append(pseg_iter->second, {prox.x, prox.y, prox.z, prox.r}, {r.x, r.y, r.z, r.r}, r.tag);
+ unused_samples.erase(prox.id);
+ }
+
+ if (!unused_samples.empty()) {
+ throw swc_single_sample_segment(*unused_samples.begin());
+ }
+ return tree;
+}
+
+arb::segment_tree load_swc_allen(std::vector<swc_record>& records, bool no_gaps) {
+ // Assert that the file contains at least one sample.
+ if (records.empty()) return {};
+
+ // Map of SWC record id to index in `records`.
+ std::unordered_map<int, std::size_t> record_index;
+
+ int soma_id = records[0].id;
+ record_index[soma_id] = 0;
+
+ // Assert that root sample has tag 1.
+ if (records[0].tag != 1) {
+ throw swc_no_soma{records[0].id};
+ }
+
+ for (std::size_t i = 1; i < records.size(); ++i) {
+ const auto& r = records[i];
+ record_index[r.id] = i;
+
+ // Find record index of the parent
+ auto p = r.parent_id;
+ auto parent_iter = record_index.find(p);
+
+ if (parent_iter == record_index.end() || records[parent_iter->second].id == r.id)
+ {
+ throw swc_no_such_parent{r.id};
+ }
+
+ // Assert that all samples have the same tag as their parent, except those attached to the soma.
+ if (p != soma_id && r.tag != records[parent_iter->second].tag) {
+ throw swc_mismatched_tags{r.id};
+ }
+
+ // Assert that all non-root samples have a tag of 2, 3, or 4.
+ if (r.tag < 2) throw swc_non_spherical_soma{r.id};
+ if (r.tag > 4) throw swc_unsupported_tag{r.id};
+ }
+
+ // Translate the morphology so that the soma is centered at the origin (0,0,0)
+ arb::mpoint sloc{records[0].x, records[0].y, records[0].z, records[0].r};
+ for (auto& r: records) {
+ r.x -= sloc.x;
+ r.y -= sloc.y;
+ r.z -= sloc.z;
+ }
+
+ arb::segment_tree tree;
+
+ // Model the spherical soma as a cylinder with length=2*radius.
+ // The cylinder is centred on the origin, and extended along the y axis.
+ double soma_rad = sloc.radius;
+ tree.append(arb::mnpos, {0, -soma_rad, 0, soma_rad}, {0, soma_rad, 0, soma_rad}, 1);
+
+ // Build branches off soma.
+ std::unordered_map<int, arb::msize_t> smap; // SWC record id -> segment id
+ std::set<int> unused_samples; // Samples that are not part of a segment
+ for (const auto& r: records) {
+ int id = r.id;
+ if (id == soma_id) continue;
+
+ // If sample i has the root as its parent don't create a segment.
+ if (r.parent_id == soma_id) {
+ if (no_gaps) {
+ // Assert that this branch starts on the "surface" of the spherical soma.
+ auto d = std::fabs(soma_rad - std::sqrt(r.x * r.x + r.y * r.y + r.z * r.z));
+ if (d > 1e-3) { // 1 nm tolerance
+ throw swc_unsupported_gaps{r.id};
+ }
+ }
+ // This maps axons and apical dendrites to soma.prox, and dendrites to soma.dist.
+ smap[id] = r.tag == 3 ? 0 : arb::mnpos;
+ unused_samples.insert(id);
+ continue;
+ }
+
+ const auto p = r.parent_id;
+ const auto& prox = records[record_index[p]];
+ smap[id] = tree.append(smap.at(p), {prox.x, prox.y, prox.z, prox.r}, {r.x, r.y, r.z, r.r}, r.tag);
+ unused_samples.erase(p);
+ }
+
+ if (!unused_samples.empty()) {
+ throw swc_single_sample_segment{*unused_samples.begin()};
+ }
+
+ return tree;
+}
+} // namespace arborio
+
diff --git a/doc/concepts/morphology.rst b/doc/concepts/morphology.rst
index 123562e33d2f17bc12bc5c311a96dd74eabad5f7..e5b03e06a4e683a86ab19f798322f9b6ee7416f1 100644
--- a/doc/concepts/morphology.rst
+++ b/doc/concepts/morphology.rst
@@ -482,6 +482,102 @@ because it has two children: the dendrites attached to its distal end.
:width: 900
:align: center
+.. _morph-formats:
+
+Supported file formats
+----------------------
+
+Arbor supports morphologies described using the SWC file format and the NeuroML file format.
+
+SWC
+~~~
+
+Arbor supports reading morphologies described using the
+`SWC <http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html>`_ file format.
+SWC files describe the morphology as a list of samples with an id, an `x,y,z` location is space, a radius, a tag
+and a parent id. Arbor parses these samples, performs some checks, then generates a
+:ref:`segment tree <morph-segment_tree>` describing the morphology according to one of three possible
+interpretations.
+
+The SWC file format specifications are not very detailed, which has lead different simulators to interpret
+SWC files in different ways, especially when it comes to the soma. Arbor has its own an interpretation that
+is powerful, and simple to understand at the same time. However, we have also developed functions that will
+interpret SWC files similarly to how the NEURON simulator would, and how the Allen Institute would.
+
+Despite the differences between the interpretations, there is a common set of checks that are always performed
+to check the validity of SWC files:
+ * Check that there are no duplicate ids.
+ * Check that the parent id of a sample is less than the id of the sample.
+ * Check that the parent id of a sample refers to an existing sample.
+
+In addition, all interpretations agree that a *segment* is (in the common case) constructed between a sample and
+its parent and inherits the tag of the sample; and if more than 1 sample have the same parent, the parent sample
+is interpreted as a fork point in the morphology, and acts as the proximal point to a new branch for each of its
+"child" samples. There a couple of exceptions to these rules which are listed below.
+
+Arbor interpretation:
+"""""""""""""""""""""
+In addition to the previously listed checks, the arbor interpretation explicitly disallows SWC files where the soma is
+described by a single sample. It constructs the soma from 2 or more samples that form 1 or more segments in the segment
+tree. A *segment* is always constructed between a sample and its parent. This means that there are no gaps in the
+resulting segment tree and morphology.
+
+Arbor has no magic rules or transformations for the soma. It can be a single branch or multiple branches; segments
+of a different tag can connect to its distal end, proximal end or anywhere in the middle. For example, to create a
+morphology with a single segment soma; a single segment axon connected to one end of the soma; and a single segment
+dendrite connected to the other end of the soma, the following swc file can be used:
+
+.. code:: Python
+
+ # id, tag, x, y, z, r, parent
+ 1, 1, 0, 0, 0, 1, -1
+ 2, 1, 2, 0, 0, 1, 1
+ 3, 2, -3, 0, 0, 0.7, 1
+ 4, 3, 20, 0, 0, 1, 2
+
+Samples 1 and 2 will form the soma; samples 1 and 3 will form the axon, connected to the soma at the proximal end;
+samples 2 and 4 will form the dendrite, connected to the soma at the distal end. The morphology will look something
+like this:
+
+.. figure:: ../gen-images/swc_morph.svg
+ :width: 400
+ :align: center
+
+
+Allen interpretation:
+"""""""""""""""""""""
+In addition to the previously mentioned checks, the Allen interpretation expects a single-sample soma to be the first
+sample of the file and to be interpreted as a spherical soma. Arbor represents the spherical soma as a cylinder with
+length and diameter equal to the diameter of the sample representing the sphere.
+
+This interpretation also expects that samples have the same tag as their parent samples, with the exception of samples
+that have the soma sample as a parent. In this case, when a sample's parent is the soma, no *segment* is created
+between the 2 samples; instead there is a gap in the morphology (represented electrically as a zero-resistance wire).
+Samples with the soma as a parent start new segments, that connect to the distal end of the soma if they are dendrites,
+or to the proximal end of the soma if they are axons or apical dendrites. Only axons, dendrites and apical dendrites
+(tags 2, 3 and 4 respectively) are allowed in this interpretation, in addition to the spherical soma.
+
+Finally the Allen institute interpretation of SWC files centers the morphology around the soma at the origin (0, 0, 0)
+and all samples are translated in space towards the origin.
+
+NEURON interpretation:
+""""""""""""""""""""""
+The NEURON interpretation was obtained by experimenting with the `Import3d_SWC_read` function. We came up with the
+following set of rules that govern NEURON's SWC behavior and enforced them in arbor's NEURON-complaint SWC
+interpreter:
+ * SWC files must contain a soma sample and it must to be the first sample.
+ * A soma is represented by a series of n≥1 unbranched, serially listed samples.
+ * A soma is constructed as a single cylinder with diameter equal to the piecewise average diameter of all the
+ segments forming the soma.
+ * A single-sample soma at is constructed as a cylinder with length=diameter.
+ * If a non-soma sample is to have a soma sample as its parent, it must have the most distal sample of the soma
+ as the parent.
+ * Every non-soma sample that has a soma sample as its parent, attaches to the created soma cylinder at its midpoint.
+ * If a non-soma sample has a soma sample as its parent, no segment is created between the sample and its parent,
+ instead that sample is the proximal point of a new segment, and there is a gap in the morphology (represented
+ electrically as a zero-resistance wire)
+ * To create a segment with a certain tag, that is to be attached to the soma, we need at least 2 samples with that
+ tag.
API
---
diff --git a/doc/cpp/cable_cell.rst b/doc/cpp/cable_cell.rst
index cab2b45a57a68fba92fdf4e2114c170610014c16..29320f6b5f58195970cffc733953ce29c2c864ae 100644
--- a/doc/cpp/cable_cell.rst
+++ b/doc/cpp/cable_cell.rst
@@ -145,6 +145,84 @@ by two stitches:
cell.paint("\"soma\"", "hh");
+Supported morphology formats:
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Arbor supports morphologies described using the SWC file format and the NeuroML file format.
+
+SWC
+"""
+
+Arbor supports reading morphologies described using the
+`SWC <http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html>`_ file format. And
+has three different interpretation of that format.
+
+A :cpp:func:`parse_swc()` function is used to parse the SWC file and generate a :cpp:type:`swc_data` object.
+This object contains a vector of :cpp:type:`swc_record` objects that represent the SWC samples, with a number of
+basic checks performed on them. The :cpp:type:`swc_data` object can then be used to generate a
+:cpp:type:`segment_tree` object using one of the following functions: (See the morphology concepts
+:ref:`page <morph-formats>` for more details).
+
+ * :cpp:func:`as_segment_tree`
+ * :cpp:func:`load_swc_allen`
+ * :cpp:func:`load_swc_neuron`
+
+.. cpp:class:: swc_record
+
+ .. cpp:member:: int id
+
+ ID of the record
+
+ .. cpp:member:: int tag
+
+ Structure identifier (tag).
+
+ .. cpp:member:: double x
+
+ x coordinate in space.
+
+ .. cpp:member:: double y
+
+ y coordinate in space.
+
+ .. cpp:member:: double z
+
+ z coordinate in space.
+
+ .. cpp:member:: double r
+
+ Sample radius.
+
+ .. cpp:member:: int parent_id
+
+ Record parent's sample ID.
+
+.. cpp:class:: swc_data
+
+ .. cpp:member:: std::string metadata
+
+ Contains the comments of an SWC file.
+
+ .. cpp:member:: std::vector<swc_record> records
+
+ Stored the list of samples from an SWC file, after performing some checks.
+
+.. cpp:function:: swc_data parse_swc(std::istream&)
+
+ Returns an `swc_data` object given an std::istream object.
+
+.. cpp:function:: segment_tree as_segment_tree(const swc_data& data)
+
+ Returns a segment tree constructed according to Arbor's SWC specifications.
+
+.. cpp:function:: segment_tree load_swc_allen(swc_data& data)
+
+ Returns a segment tree constructed according to the Allen Institute's SWC specifications.
+
+.. cpp:function:: segment_tree load_swc_neuron(swc_data& data)
+
+ Returns a segment tree constructed according to NEURON's SWC specifications.
+
.. _locsets-and-regions:
Identifying sites and subsets of the morphology
diff --git a/doc/python/morphology.rst b/doc/python/morphology.rst
index 046eef2df0b97fbe95eea160f42233c7dcd3d44b..f183f28064bca2ff017ba2d031703ff981d5ee9d 100644
--- a/doc/python/morphology.rst
+++ b/doc/python/morphology.rst
@@ -264,7 +264,8 @@ Cell morphology
.. py:function:: load_swc(filename)
- Loads the morphology in an SWC file as a :class:`segment_tree`.
+ Loads the morphology in an SWC file as a :class:`segment_tree` according to arbor's SWC specifications.
+ (See the morphology concepts :ref:`page <morph-formats>` for more details).
The samples in the SWC files are treated as the end points of segments, where a
sample and its parent form a segment.
@@ -291,6 +292,28 @@ Cell morphology
:param str filename: the name of the SWC file.
:rtype: segment_tree
+.. py:function:: load_swc_neuron(filename)
+
+ Loads the morphology in an SWC file as a :class:`segment_tree` according to NEURON's SWC specifications.
+ Specifically:
+ * The first sample must be a soma sample.
+ * The soma is represented by a series of n≥1 unbranched, serially listed samples.
+ * The soma is constructed as a single cylinder with diameter equal to the piecewise average diameter of all the
+ segments forming the soma.
+ * A single-sample soma at is constructed as a cylinder with length=diameter.
+ * If a non-soma sample is to have a soma sample as its parent, it must have the most distal sample of the soma
+ as the parent.
+ * Every non-soma sample that has a soma sample as its parent, attaches to the created soma cylinder at its midpoint.
+ * If a non-soma sample has a soma sample as its parent, no segment is created between the sample and its parent,
+ instead that sample is the proximal point of a new segment, and there is a gap in the morphology (represented
+ electrically as a zero-resistance wire)
+ * To create a segment with a certain tag, that is to be attached to the soma, we need at least 2 samples with that
+ tag.
+
+ :param str filename: the name of the SWC file.
+ :rtype: segment_tree
+
+
.. py:function:: load_swc_allen(filename, no_gaps=False)
Generate a segment tree from an SWC file following the rules prescribed by
diff --git a/doc/scripts/inputs.py b/doc/scripts/inputs.py
index d3cd78efc1b39735e0578c9b30c5ff253168658f..1a323c3ae012acfe155d41950d7719c3111b2a2e 100644
--- a/doc/scripts/inputs.py
+++ b/doc/scripts/inputs.py
@@ -84,6 +84,11 @@ tmp = [
[[Segment((-3.0, 0.0, 1.5), (-5.5, -0.2, 0.5), 2), Segment((-5.5, -0.2, 0.5), (-14.5, -0.1, 0.5), 2)]],]
ysoma_morph3 = representation.make_morph(tmp)
+tmp = [
+ [[Segment((-3.0, 0.0, 0.7), (0.0, 0.0, 1.0), 2)], [Segment((0.0, 0.0, 1.0), (2.0, 0.0, 1.0), 1)], [Segment((2.0, 0.0, 1.0), (20.0, 0.0, 1.0), 3)]],
+]
+swc_morph = representation.make_morph(tmp)
+
############# locsets
diff --git a/doc/scripts/make_images.py b/doc/scripts/make_images.py
index 0df83945c9063a7b80ce3fd821eae19b240c6b1b..b1fbb2de2a1833c893b26d29e5fc4dfb81932a9d 100644
--- a/doc/scripts/make_images.py
+++ b/doc/scripts/make_images.py
@@ -235,6 +235,7 @@ def generate(path=''):
morph_image([inputs.label_morph], ['segments'], path+'/label_seg.svg')
morph_image([inputs.detached_morph], ['segments'], path+'/detached_seg.svg')
morph_image([inputs.stacked_morph], ['segments'], path+'/stacked_seg.svg')
+ morph_image([inputs.swc_morph], ['segments'], path+'/swc_morph.svg')
morph_image([inputs.label_morph, inputs.label_morph], ['segments', 'branches'], path+'/label_morph.svg')
morph_image([inputs.detached_morph, inputs.detached_morph], ['segments', 'branches'], path+'/detached_morph.svg')
diff --git a/example/single/CMakeLists.txt b/example/single/CMakeLists.txt
index ad5257c6e70364152b03ba1b4637d75b729bd8ac..1a2c889e611be54bb888698cf19106e2b9058f6e 100644
--- a/example/single/CMakeLists.txt
+++ b/example/single/CMakeLists.txt
@@ -1,3 +1,3 @@
add_executable(single-cell EXCLUDE_FROM_ALL single.cpp)
add_dependencies(examples single-cell)
-target_link_libraries(single-cell PRIVATE arbor arborenv arbor-sup ext-tinyopt)
+target_link_libraries(single-cell PRIVATE arbor arborio arborenv arbor-sup ext-tinyopt)
diff --git a/example/single/single.cpp b/example/single/single.cpp
index f280096269998f2b8d977e96dfa152da77e446a3..f859009f24585996c62bc6fa8fafd8b95d8940d3 100644
--- a/example/single/single.cpp
+++ b/example/single/single.cpp
@@ -10,10 +10,11 @@
#include <arbor/cable_cell.hpp>
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/segment_tree.hpp>
-#include <arbor/swcio.hpp>
#include <arbor/simulation.hpp>
#include <arbor/simple_sampler.hpp>
+#include <arborio/swcio.hpp>
+
#include <tinyopt/tinyopt.h>
struct options {
@@ -157,5 +158,5 @@ arb::morphology read_swc(const std::string& path) {
std::ifstream f(path);
if (!f) throw std::runtime_error("unable to open SWC file: "+path);
- return arb::morphology(arb::as_segment_tree(arb::parse_swc(f)));
+ return arb::morphology(arborio::as_segment_tree(arborio::parse_swc(f)));
}
diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index 5344e1f771575c7cbbd8090c38bd9d60c4a1786e..de05eeb06c4bdf9c764692a1493d21d1a3a1113c 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -40,7 +40,7 @@ set(pyarb_source
# use by both the Python wrapper target (pyarb) and for the
# unit tests of the C++ components in the Python wrapper.
add_library(pyarb_obj OBJECT ${pyarb_source})
-target_link_libraries(pyarb_obj PRIVATE arbor pybind11::module)
+target_link_libraries(pyarb_obj PRIVATE arbor arborio pybind11::module)
# The Python library. MODULE will make a Python-exclusive model.
add_library(pyarb MODULE $<TARGET_OBJECTS:pyarb_obj>)
@@ -53,7 +53,7 @@ set_target_properties(pyarb PROPERTIES PREFIX "${PYTHON_MODULE_PREFIX}" SUFFIX "
# This dependency has to be spelt out again, despite being added to
# pyarb_obj because CMake.
-target_link_libraries(pyarb PRIVATE arbor pybind11::module)
+target_link_libraries(pyarb PRIVATE arbor arborio pybind11::module)
# Add support for mpi4py if available.
if (ARB_WITH_MPI)
diff --git a/python/morphology.cpp b/python/morphology.cpp
index ef44b2dad621947d3ce0fdeedf0489b69428a9a9..dd420097ce4d97f22ab3e4bd111bc580d6399bcd 100644
--- a/python/morphology.cpp
+++ b/python/morphology.cpp
@@ -6,13 +6,30 @@
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/primitives.hpp>
#include <arbor/morph/segment_tree.hpp>
-#include <arbor/swcio.hpp>
+
+#include <arborio/swcio.hpp>
#include "error.hpp"
#include "strprintf.hpp"
namespace pyarb {
+arb::segment_tree load_swc_neuron(const std::string& fname) {
+ std::ifstream fid{fname};
+ if (!fid.good()) {
+ throw pyarb_error(util::pprintf("can't open file '{}'", fname));
+ }
+ try {
+ auto records = arborio::parse_swc(fid, arborio::swc_mode::relaxed).records;
+ return arborio::load_swc_neuron(records);
+ }
+ catch (arborio::swc_error& e) {
+ // Try to produce helpful error messages for SWC parsing errors.
+ throw pyarb_error(
+ util::pprintf("NEURON SWC: error parsing {}: {}", fname, e.what()));
+ }
+}
+
arb::segment_tree load_swc_allen(const std::string& fname, bool no_gaps=false) {
std::ifstream fid{fname};
if (!fid.good()) {
@@ -20,91 +37,10 @@ arb::segment_tree load_swc_allen(const std::string& fname, bool no_gaps=false) {
}
try {
using namespace arb;
- auto records = parse_swc(fid, swc_mode::relaxed).records;
-
- // Assert that the file contains at least one sample.
- if (records.empty()) {
- throw pyarb_error("Allen SWC: empty file");
- }
-
- // Map of SWC record id to index in `records`.
- std::unordered_map<int, std::size_t> record_index;
-
- int soma_id = records[0].id;
- record_index[soma_id] = 0;
-
- // Assert that root sample has tag 1.
- if (records[0].tag!=1) {
- throw pyarb_error("Allen SWC: the soma record does not have tag 1");
- }
-
- for (std::size_t i = 1; i<records.size(); ++i) {
- const auto& r = records[i];
- record_index[r.id] = i;
-
- // Assert that all samples have the same tag as their parent, except those attached to the soma.
- if (auto p = r.parent_id; p!=soma_id && r.tag!=records[record_index[p]].tag) {
- throw pyarb_error(
- "Allen SWC: every record not attached to the soma must have the same tag as its parent");
- }
-
- // Assert that all non-root samples have a tag of 2, 3, or 4.
- if (r.tag<2 || r.tag>4) {
- throw pyarb_error(
- "Allen SWC: every record must have a tag of 2, 3 or 4, except for the first which must have tag 1");
- }
- }
-
- // Translate the morphology so that the soma is centered at the origin (0,0,0)
- mpoint sloc{records[0].x, records[0].y, records[0].z, records[0].r};
- for (auto& r: records) {
- r.x -= sloc.x;
- r.y -= sloc.y;
- r.z -= sloc.z;
- }
-
- segment_tree tree;
-
- // Model the spherical soma as a cylinder with length=2*radius.
- // The cylinder is centred on the origin, and extended along the y axis.
- double soma_rad = sloc.radius;
- tree.append(mnpos, {0, -soma_rad, 0, soma_rad}, {0, soma_rad, 0, soma_rad}, 1);
-
- // Build branches off soma.
- std::unordered_map<int, msize_t> smap; // SWC record id -> segment id
- std::set<int> unused_samples;
- for (const auto& r: records) {
- int id = r.id;
- if (id==soma_id) continue;
-
- // If sample i has the root as its parent don't create a segment.
- if (r.parent_id==soma_id) {
- if (no_gaps) {
- // Assert that this branch starts on the "surface" of the spherical soma.
- auto d = std::fabs(soma_rad - std::sqrt(r.x*r.x + r.y*r.y + r.z*r.z));
- if (d>1e-3) { // 1 nm tolerance
- throw pyarb_error("Allen SWC: no gaps are allowed between the soma and any axons, dendrites or apical dendrites");
- }
- }
- // This maps axons and apical dendrites to soma.prox, and dendrites to soma.dist.
- smap[id] = r.tag==3? 0: mnpos;
- unused_samples.insert(id);
- continue;
- }
-
- const auto p = r.parent_id;
- const auto& prox = records[record_index[p]];
- smap[id] = tree.append(smap.at(p), {prox.x, prox.y, prox.z, prox.r}, {r.x, r.y, r.z, r.r}, r.tag);
- unused_samples.erase(p);
- }
-
- if (!unused_samples.empty()) {
- throw pyarb_error("Allen SWC: Every branch must contain at least one segment");
- }
-
- return tree;
+ auto records = arborio::parse_swc(fid, arborio::swc_mode::relaxed).records;
+ return arborio::load_swc_allen(records, no_gaps);
}
- catch (arb::swc_error& e) {
+ catch (arborio::swc_error& e) {
// Try to produce helpful error messages for SWC parsing errors.
throw pyarb_error(
util::pprintf("Allen SWC: error parsing {}: {}", fname, e.what()));
@@ -239,10 +175,10 @@ void register_morphology(pybind11::module& m) {
throw pyarb_error(util::pprintf("can't open file '{}'", fname));
}
try {
- auto records = arb::parse_swc(fid).records;
- return arb::as_segment_tree(records);
+ auto records = arborio::parse_swc(fid).records;
+ return arborio::as_segment_tree(records);
}
- catch (arb::swc_error& e) {
+ catch (arborio::swc_error& e) {
// Try to produce helpful error messages for SWC parsing errors.
throw pyarb_error(
util::pprintf("error parsing {}: {}", fname, e.what()));
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 33b43cabd681e7986804a0feed8e82afccbd2908..3fa832387e893e792ef7304dd6de792ef73d9486 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -220,7 +220,7 @@ endif()
target_compile_options(unit PRIVATE ${ARB_CXXOPT_ARCH})
target_compile_definitions(unit PRIVATE "-DDATADIR=\"${CMAKE_CURRENT_SOURCE_DIR}/swc\"")
target_include_directories(unit PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
-target_link_libraries(unit PRIVATE gtest arbor arborenv arbor-private-headers arbor-sup)
+target_link_libraries(unit PRIVATE gtest arbor arborenv arborio arbor-private-headers arbor-sup)
if(ARB_WITH_NEUROML)
target_link_libraries(unit PRIVATE arbornml)
diff --git a/test/unit/test_morphology.cpp b/test/unit/test_morphology.cpp
index f35d8217a582049ff709a43b5537150f17a6b1f3..4f110981f03574841cb17857e01a5dc9fb41f8aa 100644
--- a/test/unit/test_morphology.cpp
+++ b/test/unit/test_morphology.cpp
@@ -10,7 +10,8 @@
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/primitives.hpp>
#include <arbor/morph/segment_tree.hpp>
-#include <arbor/swcio.hpp>
+
+#include <arborio/swcio.hpp>
#include "util/span.hpp"
@@ -318,10 +319,10 @@ TEST(morphology, swc) {
}
// Load swc samples from file.
- auto swc = arb::parse_swc(fid, arb::swc_mode::strict);
+ auto swc = arborio::parse_swc(fid, arborio::swc_mode::strict);
// Build a segmewnt_tree from swc samples.
- auto sm = arb::as_segment_tree(swc);
+ auto sm = arborio::as_segment_tree(swc);
EXPECT_EQ(5798u, sm.size()); // SWC data contains 5799 samples.
// Test that the morphology contains the expected number of branches.
diff --git a/test/unit/test_swcio.cpp b/test/unit/test_swcio.cpp
index debcb072109ed22912b6bd68b3ed9b0d497bed93..9de95c337c6d77d3ed36f1ba1d8069f8844d85da 100644
--- a/test/unit/test_swcio.cpp
+++ b/test/unit/test_swcio.cpp
@@ -6,7 +6,8 @@
#include <arbor/cable_cell.hpp>
#include <arbor/morph/primitives.hpp>
#include <arbor/morph/segment_tree.hpp>
-#include <arbor/swcio.hpp>
+
+#include <arborio/swcio.hpp>
#include "../gtest.h"
@@ -16,7 +17,10 @@
# define DATADIR "../data"
#endif
-using namespace arb;
+using namespace arborio;
+using arb::segment_tree;
+using arb::mpoint;
+using arb::mnpos;
TEST(swc_record, construction) {
swc_record record(1, 7, 1., 2., 3., 4., -1);
@@ -208,14 +212,14 @@ TEST(swc_parser, segment_tree) {
{1, 1, 0., 0., 0., 1., -1},
{5, 3, 1., 1., 1., 1., 2}
};
- EXPECT_THROW(as_segment_tree(swc), bad_swc_data);
+ EXPECT_THROW(as_segment_tree(swc), swc_no_such_parent);
}
{
// A single SWC record will throw.
std::vector<swc_record> swc{
{1, 1, 0., 0., 0., 1., -1}
};
- EXPECT_THROW(as_segment_tree(swc), bad_swc_data);
+ EXPECT_THROW(as_segment_tree(swc), swc_bad_description);
}
{
// Otherwise, ensure segment ends and tags correspond.
@@ -257,6 +261,666 @@ TEST(swc_parser, segment_tree) {
EXPECT_EQ(p4, tree.segments()[3].dist);
}
}
+TEST(swc_parser, allen_compliant) {
+ using namespace arborio;
+ {
+ // One-point soma; interpretted as 1 segment
+ mpoint p0{0, 0, 0, 10};
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1}
+ };
+ segment_tree tree = load_swc_allen(swc);
+
+ mpoint prox{p0.x, p0.y-p0.radius, p0.z, p0.radius};
+ mpoint dist{p0.x, p0.y+p0.radius, p0.z, p0.radius};
+
+ ASSERT_EQ(1u, tree.segments().size());
+
+ EXPECT_EQ(mnpos, tree.parents()[0]);
+ EXPECT_EQ(1, tree.segments()[0].tag);
+ EXPECT_EQ(prox, tree.segments()[0].prox);
+ EXPECT_EQ(dist, tree.segments()[0].dist);
+ }
+ {
+ // One-point soma, two-point dendrite
+ mpoint p0{0, 0, 0, 10};
+ mpoint p1{0, 0, 0, 5};
+ mpoint p2{0, 200, 0, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 3, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 3, p2.x, p2.y, p2.z, p2.radius, 2}
+ };
+ segment_tree tree = load_swc_allen(swc);
+
+ mpoint prox{0, -10, 0, 10};
+ mpoint dist{0, 10, 0, 10};
+
+ ASSERT_EQ(2u, tree.segments().size());
+
+ EXPECT_EQ(mnpos, tree.parents()[0]);
+ EXPECT_EQ(1, tree.segments()[0].tag);
+ EXPECT_EQ(prox, tree.segments()[0].prox);
+ EXPECT_EQ(dist, tree.segments()[0].dist);
+
+ EXPECT_EQ(0u, tree.parents()[1]);
+ EXPECT_EQ(3, tree.segments()[1].tag);
+ EXPECT_EQ(p1, tree.segments()[1].prox);
+ EXPECT_EQ(p2, tree.segments()[1].dist);
+ }
+ {
+ // 1-point soma, 2-point dendrite, 2-point axon
+ mpoint p0{0, 0, 0, 1};
+ mpoint p1{0, 0, 10, 10};
+ mpoint p2{0, 0, 20, 10};
+ mpoint p3{0, 0, 21, 10};
+ mpoint p4{0, 0, 30, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 3, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 3, p2.x, p2.y, p2.z, p2.radius, 2},
+ {4, 2, p3.x, p3.y, p3.z, p3.radius, 1},
+ {5, 2, p4.x, p4.y, p4.z, p4.radius, 4}
+ };
+ segment_tree tree = load_swc_allen(swc);
+
+ mpoint prox{0, -1, 0, 1};
+ mpoint dist{0, 1, 0, 1};
+
+ ASSERT_EQ(3u, tree.segments().size());
+
+ EXPECT_EQ(mnpos, tree.parents()[0]);
+ EXPECT_EQ(1, tree.segments()[0].tag);
+ EXPECT_EQ(prox, tree.segments()[0].prox);
+ EXPECT_EQ(dist, tree.segments()[0].dist);
+
+ EXPECT_EQ(0u, tree.parents()[1]);
+ EXPECT_EQ(3, tree.segments()[1].tag);
+ EXPECT_EQ(p1, tree.segments()[1].prox);
+ EXPECT_EQ(p2, tree.segments()[1].dist);
+
+ EXPECT_EQ(mnpos, tree.parents()[2]);
+ EXPECT_EQ(2, tree.segments()[2].tag);
+ EXPECT_EQ(p3, tree.segments()[2].prox);
+ EXPECT_EQ(p4, tree.segments()[2].dist);
+ }
+}
+
+TEST(swc_parser, not_allen_compliant) {
+ using namespace arborio;
+ {
+ // multi-point soma
+ mpoint p0{0, 0, -10, 10};
+ mpoint p1{0, 0, 0, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1}
+ };
+ EXPECT_THROW(load_swc_allen(swc), swc_non_spherical_soma);
+ }
+ {
+ // unsupported tag
+ mpoint p0{0, 0, 0, 1};
+ mpoint p1{0, 200, 20, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 5, p1.x, p1.y, p1.z, p1.radius, 1}
+ };
+ EXPECT_THROW(load_swc_allen(swc), swc_unsupported_tag);
+ }
+ {
+ // 1-point soma; 2-point dendrite; 1-point axon connected to the proximal end of the dendrite
+ mpoint p0{0, 0, -15, 10};
+ mpoint p1{0, 0, 0, 10};
+ mpoint p2{0, 0, 80, 10};
+ mpoint p3{0, 0, -80, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 3, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 3, p2.x, p2.y, p2.z, p2.radius, 2},
+ {4, 2, p3.x, p3.y, p3.z, p3.radius, 2}
+ };
+ EXPECT_THROW(load_swc_allen(swc), swc_mismatched_tags);
+ }
+ {
+ // 1-point soma and 1-point dendrite
+ mpoint p0{0, 0, 0, 10};
+ mpoint p1{0, 200, 0, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 3, p1.x, p1.y, p1.z, p1.radius, 1}
+ };
+ EXPECT_THROW(load_swc_allen(swc), swc_single_sample_segment);
+ }
+ {
+ // 2-point dendrite and 1-point soma at the end
+ mpoint p0{0, 0, 0, 1};
+ mpoint p1{0, 0, 10, 1};
+ mpoint p2{0, 200, 20, 10};
+
+ std::vector<swc_record> swc{
+ {1, 3, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 3, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 1, p2.x, p2.y, p2.z, p2.radius, 2}
+ };
+ EXPECT_THROW(load_swc_allen(swc), swc_no_soma);
+ }
+ {
+ // non-existent parent sample
+ mpoint p0{0, 0, 0, 1};
+ mpoint p1{0, 200, 20, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 3, p1.x, p1.y, p1.z, p1.radius, 4}
+ };
+ EXPECT_THROW(load_swc_allen(swc), swc_no_such_parent);
+ }
+ {
+ // parent sample is self
+ mpoint p0{0, 0, 0, 1};
+ mpoint p1{0, 200, 20, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 2}
+ };
+ EXPECT_THROW(load_swc_allen(swc), swc_no_such_parent);
+ }
+}
+
+TEST(swc_parser, neuron_compliant) {
+ using namespace arborio;
+ {
+ // One-point soma; interpretted as 1 segment
+
+ mpoint p0{0, 0, 0, 10};
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1}
+ };
+ segment_tree tree = load_swc_neuron(swc);
+
+ mpoint prox{p0.x, p0.y-p0.radius, p0.z, p0.radius};
+ mpoint dist{p0.x, p0.y+p0.radius, p0.z, p0.radius};
+
+ ASSERT_EQ(1u, tree.segments().size());
+
+ EXPECT_EQ(mnpos, tree.parents()[0]);
+ EXPECT_EQ(1, tree.segments()[0].tag);
+ EXPECT_EQ(prox, tree.segments()[0].prox);
+ EXPECT_EQ(dist, tree.segments()[0].dist);
+ }
+ {
+ // Two-point soma; interpretted as 1
+ mpoint p0{0, 0, -10, 10};
+ mpoint p1{0, 0, 0, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1}
+ };
+ segment_tree tree = load_swc_neuron(swc);
+
+ ASSERT_EQ(1u, tree.segments().size());
+
+ EXPECT_EQ(mnpos, tree.parents()[0]);
+ EXPECT_EQ(1, tree.segments()[0].tag);
+ EXPECT_EQ(p0, tree.segments()[0].prox);
+ EXPECT_EQ(p1, tree.segments()[0].dist);
+ }
+ {
+ // Three-point soma; interpretted as 2 segments
+ mpoint p0{0, 0, -10, 10};
+ mpoint p1{0, 0, 0, 10};
+ mpoint p2{0, 0, 10, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 1, p2.x, p2.y, p2.z, p2.radius, 2}
+ };
+ segment_tree tree = load_swc_neuron(swc);
+
+ ASSERT_EQ(2u, tree.segments().size());
+
+ EXPECT_EQ(mnpos, tree.parents()[0]);
+ EXPECT_EQ(1, tree.segments()[0].tag);
+ EXPECT_EQ(p0, tree.segments()[0].prox);
+ EXPECT_EQ(p1, tree.segments()[0].dist);
+
+ EXPECT_EQ(0u, tree.parents()[1]);
+ EXPECT_EQ(1, tree.segments()[1].tag);
+ EXPECT_EQ(p1, tree.segments()[1].prox);
+ EXPECT_EQ(p2, tree.segments()[1].dist);
+ }
+ {
+ // 6-point soma; interpretted as 5 segments
+ mpoint p0{0, 0, -5, 2};
+ mpoint p1{0, 0, 0, 5};
+ mpoint p2{0, 0, 2, 6};
+ mpoint p3{0, 0, 6, 1};
+ mpoint p4{0, 0, 10, 7};
+ mpoint p5{0, 0, 15, 2};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 1, p2.x, p2.y, p2.z, p2.radius, 2},
+ {4, 1, p3.x, p3.y, p3.z, p3.radius, 3},
+ {5, 1, p4.x, p4.y, p4.z, p4.radius, 4},
+ {6, 1, p5.x, p5.y, p5.z, p5.radius, 5}
+ };
+ segment_tree tree = load_swc_neuron(swc);
+
+
+ ASSERT_EQ(5u, tree.segments().size());
+
+ EXPECT_EQ(mnpos, tree.parents()[0]);
+ EXPECT_EQ(1, tree.segments()[0].tag);
+ EXPECT_EQ(p0, tree.segments()[0].prox);
+ EXPECT_EQ(p1, tree.segments()[0].dist);
+
+ EXPECT_EQ(0u, tree.parents()[1]);
+ EXPECT_EQ(1, tree.segments()[1].tag);
+ EXPECT_EQ(p1, tree.segments()[1].prox);
+ EXPECT_EQ(p2, tree.segments()[1].dist);
+
+ EXPECT_EQ(1u, tree.parents()[2]);
+ EXPECT_EQ(1, tree.segments()[2].tag);
+ EXPECT_EQ(p2, tree.segments()[2].prox);
+ EXPECT_EQ(p3, tree.segments()[2].dist);
+
+ EXPECT_EQ(2u, tree.parents()[3]);
+ EXPECT_EQ(1, tree.segments()[3].tag);
+ EXPECT_EQ(p3, tree.segments()[3].prox);
+ EXPECT_EQ(p4, tree.segments()[3].dist);
+
+ EXPECT_EQ(3u, tree.parents()[4]);
+ EXPECT_EQ(1, tree.segments()[4].tag);
+ EXPECT_EQ(p4, tree.segments()[4].prox);
+ EXPECT_EQ(p5, tree.segments()[4].dist);
+ }
+ {
+ // One-point soma, two-point dendrite
+ mpoint p0{0, 0, 0, 10};
+ mpoint p1{0, 0, 0, 5};
+ mpoint p2{0, 200, 0, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 3, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 3, p2.x, p2.y, p2.z, p2.radius, 2}
+ };
+ segment_tree tree = load_swc_neuron(swc);
+
+ mpoint prox{0, -10, 0, 10};
+ mpoint dist{0, 10, 0, 10};
+
+ ASSERT_EQ(3u, tree.segments().size());
+
+ EXPECT_EQ(mnpos, tree.parents()[0]);
+ EXPECT_EQ(1, tree.segments()[0].tag);
+ EXPECT_EQ(prox, tree.segments()[0].prox);
+ EXPECT_EQ(p0, tree.segments()[0].dist);
+
+ EXPECT_EQ(0u, tree.parents()[1]);
+ EXPECT_EQ(1, tree.segments()[1].tag);
+ EXPECT_EQ(p0, tree.segments()[1].prox);
+ EXPECT_EQ(dist, tree.segments()[1].dist);
+
+ EXPECT_EQ(0u, tree.parents()[2]);
+ EXPECT_EQ(3, tree.segments()[2].tag);
+ EXPECT_EQ(p1, tree.segments()[2].prox);
+ EXPECT_EQ(p2, tree.segments()[2].dist);
+ }
+ {
+ // 6-point soma, 2-point dendrite
+ mpoint p0{0, 0, -5, 2};
+ mpoint p1{0, 0, 0, 5};
+ mpoint p2{0, 0, 2, 6};
+ mpoint p3{0, 0, 6, 1};
+ mpoint p4{0, 0, 10, 7};
+ mpoint p5{0, 0, 15, 2};
+ mpoint p6{0, 0, 16, 1};
+ mpoint p7{0, 0, 55, 9};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 1, p2.x, p2.y, p2.z, p2.radius, 2},
+ {4, 1, p3.x, p3.y, p3.z, p3.radius, 3},
+ {5, 1, p4.x, p4.y, p4.z, p4.radius, 4},
+ {6, 1, p5.x, p5.y, p5.z, p5.radius, 5},
+ {7, 3, p6.x, p6.y, p6.z, p6.radius, 6},
+ {8, 3, p7.x, p7.y, p7.z, p7.radius, 7}
+ };
+ segment_tree tree = load_swc_neuron(swc);
+
+ mpoint mid {0, 0, 5, 2.25};
+
+ ASSERT_EQ(7u, tree.segments().size());
+
+ EXPECT_EQ(mnpos, tree.parents()[0]);
+ EXPECT_EQ(1, tree.segments()[0].tag);
+ EXPECT_EQ(p0, tree.segments()[0].prox);
+ EXPECT_EQ(p1, tree.segments()[0].dist);
+
+ EXPECT_EQ(0u, tree.parents()[1]);
+ EXPECT_EQ(1, tree.segments()[1].tag);
+ EXPECT_EQ(p1, tree.segments()[1].prox);
+ EXPECT_EQ(p2, tree.segments()[1].dist);
+
+ EXPECT_EQ(1u, tree.parents()[2]);
+ EXPECT_EQ(1, tree.segments()[2].tag);
+ EXPECT_EQ(p2, tree.segments()[2].prox);
+ EXPECT_EQ(mid, tree.segments()[2].dist);
+
+ EXPECT_EQ(2u, tree.parents()[3]);
+ EXPECT_EQ(1, tree.segments()[3].tag);
+ EXPECT_EQ(mid, tree.segments()[3].prox);
+ EXPECT_EQ(p3, tree.segments()[3].dist);
+
+ EXPECT_EQ(3u, tree.parents()[4]);
+ EXPECT_EQ(1, tree.segments()[4].tag);
+ EXPECT_EQ(p3, tree.segments()[4].prox);
+ EXPECT_EQ(p4, tree.segments()[4].dist);
+
+ EXPECT_EQ(4u, tree.parents()[5]);
+ EXPECT_EQ(1, tree.segments()[5].tag);
+ EXPECT_EQ(p4, tree.segments()[5].prox);
+ EXPECT_EQ(p5, tree.segments()[5].dist);
+
+ EXPECT_EQ(2u, tree.parents()[6]);
+ EXPECT_EQ(3, tree.segments()[6].tag);
+ EXPECT_EQ(p6, tree.segments()[6].prox);
+ EXPECT_EQ(p7, tree.segments()[6].dist);
+ }
+ {
+ // Two-point soma, two-point dendrite
+ mpoint p0{0, 0, -20, 10};
+ mpoint p1{0, 0, 0, 4};
+ mpoint p2{0, 0, 0, 10};
+ mpoint p3{0, 200, 0, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 3, p2.x, p2.y, p2.z, p2.radius, 2},
+ {4, 3, p3.x, p3.y, p3.z, p3.radius, 3}
+ };
+ segment_tree tree = load_swc_neuron(swc);
+
+ mpoint mid{0, 0, -10, 7};
+
+ ASSERT_EQ(3u, tree.segments().size());
+
+ EXPECT_EQ(mnpos, tree.parents()[0]);
+ EXPECT_EQ(1, tree.segments()[0].tag);
+ EXPECT_EQ(p0, tree.segments()[0].prox);
+ EXPECT_EQ(mid, tree.segments()[0].dist);
+
+ EXPECT_EQ(0u, tree.parents()[1]);
+ EXPECT_EQ(1, tree.segments()[1].tag);
+ EXPECT_EQ(mid, tree.segments()[1].prox);
+ EXPECT_EQ(p1, tree.segments()[1].dist);
+
+ EXPECT_EQ(0u, tree.parents()[2]);
+ EXPECT_EQ(3, tree.segments()[2].tag);
+ EXPECT_EQ(p2, tree.segments()[2].prox);
+ EXPECT_EQ(p3, tree.segments()[2].dist);
+ }
+ {
+ // 2-point soma; 2-point dendrite; 1-point axon connected to the proximal end of the dendrite
+ mpoint p0{0, 0, -15, 10};
+ mpoint p1{0, 0, 0, 3};
+ mpoint p2{0, 0, 0, 10};
+ mpoint p3{0, 0, 80, 10};
+ mpoint p4{0, 0, -80, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 3, p2.x, p2.y, p2.z, p2.radius, 2},
+ {4, 3, p3.x, p3.y, p3.z, p3.radius, 3},
+ {5, 2, p4.x, p4.y, p4.z, p4.radius, 3}
+ };
+ segment_tree tree = load_swc_neuron(swc);
+
+ mpoint mid{0, 0, -7.5, 6.5};
+
+ ASSERT_EQ(4u, tree.segments().size());
+
+ EXPECT_EQ(mnpos, tree.parents()[0]);
+ EXPECT_EQ(1, tree.segments()[0].tag);
+ EXPECT_EQ(p0, tree.segments()[0].prox);
+ EXPECT_EQ(mid, tree.segments()[0].dist);
+
+ EXPECT_EQ(0u, tree.parents()[1]);
+ EXPECT_EQ(1, tree.segments()[1].tag);
+ EXPECT_EQ(mid, tree.segments()[1].prox);
+ EXPECT_EQ(p1, tree.segments()[1].dist);
+
+ EXPECT_EQ(0u, tree.parents()[2]);
+ EXPECT_EQ(3, tree.segments()[2].tag);
+ EXPECT_EQ(p2, tree.segments()[2].prox);
+ EXPECT_EQ(p3, tree.segments()[2].dist);
+
+ EXPECT_EQ(0u, tree.parents()[3]);
+ EXPECT_EQ(2, tree.segments()[3].tag);
+ EXPECT_EQ(p2, tree.segments()[3].prox);
+ EXPECT_EQ(p4, tree.segments()[3].dist);
+ }
+ {
+ // 2-point soma, 2-point dendrite, 2-point axon
+ mpoint p0{0, 0, 0, 1};
+ mpoint p1{0, 0, 9, 2};
+ mpoint p2{0, 0, 10, 10};
+ mpoint p3{0, 0, 20, 10};
+ mpoint p4{0, 0, 21, 10};
+ mpoint p5{0, 0, 30, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 3, p2.x, p2.y, p2.z, p2.radius, 2},
+ {4, 3, p3.x, p3.y, p3.z, p3.radius, 3},
+ {5, 2, p4.x, p4.y, p4.z, p4.radius, 4},
+ {6, 2, p5.x, p5.y, p5.z, p5.radius, 5}
+ };
+ segment_tree tree = load_swc_neuron(swc);
+
+ mpoint mid{0, 0, 4.5, 1.5};
+
+ ASSERT_EQ(5u, tree.segments().size());
+
+ EXPECT_EQ(mnpos, tree.parents()[0]);
+ EXPECT_EQ(1, tree.segments()[0].tag);
+ EXPECT_EQ(p0, tree.segments()[0].prox);
+ EXPECT_EQ(mid, tree.segments()[0].dist);
+
+ EXPECT_EQ(0u, tree.parents()[1]);
+ EXPECT_EQ(1, tree.segments()[1].tag);
+ EXPECT_EQ(mid, tree.segments()[1].prox);
+ EXPECT_EQ(p1, tree.segments()[1].dist);
+
+ EXPECT_EQ(0u, tree.parents()[2]);
+ EXPECT_EQ(3, tree.segments()[2].tag);
+ EXPECT_EQ(p2, tree.segments()[2].prox);
+ EXPECT_EQ(p3, tree.segments()[2].dist);
+
+ EXPECT_EQ(2u, tree.parents()[3]);
+ EXPECT_EQ(2, tree.segments()[3].tag);
+ EXPECT_EQ(p3, tree.segments()[3].prox);
+ EXPECT_EQ(p4, tree.segments()[3].dist);
+
+ EXPECT_EQ(3u, tree.parents()[4]);
+ EXPECT_EQ(2, tree.segments()[4].tag);
+ EXPECT_EQ(p4, tree.segments()[4].prox);
+ EXPECT_EQ(p5, tree.segments()[4].dist);
+ }
+}
+
+TEST(swc_parser, not_neuron_compliant) {
+ using namespace arborio;
+ {
+ // Two-point collocated soma
+ mpoint p0{0, 0, 0, 5};
+ mpoint p1{0, 0, 0, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1}
+ };
+
+ EXPECT_THROW(load_swc_neuron(swc), swc_collocated_soma);
+ }
+ {
+ // 3-point soma joined in the middle (1-0-2)
+ mpoint p0{0, 0, 0, 10};
+ mpoint p1{0, 0, -10, 10};
+ mpoint p2{0, 0, 10, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 1, p2.x, p2.y, p2.z, p2.radius, 1}
+ };
+ EXPECT_THROW(load_swc_neuron(swc), swc_non_serial_soma);
+ }
+ {
+ // 4-point branching soma
+ mpoint p0{0, 0, 0, 10};
+ mpoint p1{0, 0, 10, 10};
+ mpoint p2{0, -5, 20, 10};
+ mpoint p3{0, 5, 20, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 1, p2.x, p2.y, p2.z, p2.radius, 2},
+ {4, 1, p3.x, p3.y, p3.z, p3.radius, 2}
+ };
+ EXPECT_THROW(load_swc_neuron(swc), swc_non_serial_soma);
+ }
+ {
+ // 1-point soma and 1-point dendrite
+ mpoint p0{0, 0, 0, 10};
+ mpoint p1{0, 200, 0, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 3, p1.x, p1.y, p1.z, p1.radius, 1}
+ };
+ EXPECT_THROW(load_swc_neuron(swc), swc_single_sample_segment);
+ }
+ {
+ // 2-point soma and 1-point dendrite
+ mpoint p0{0, 0, -10, 10};
+ mpoint p1{0, 0, 0, 10};
+ mpoint p2{0, 200, 0, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 3, p2.x, p2.y, p2.z, p2.radius, 2}
+ };
+ EXPECT_THROW(load_swc_neuron(swc), swc_single_sample_segment);
+ }
+ {
+ // 2-point soma and two 1-point dendrite
+ mpoint p0{0, 0, -20, 10};
+ mpoint p1{0, 0, 0, 10};
+ mpoint p2{0, -5, 80, 10};
+ mpoint p3{0, 5, -90, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 3, p2.x, p2.y, p2.z, p2.radius, 2},
+ {4, 3, p3.x, p3.y, p3.z, p3.radius, 2}
+ };
+ EXPECT_THROW(load_swc_neuron(swc), swc_single_sample_segment);
+ }
+ {
+ // 2-point dendrite and 1-point soma at the end
+ mpoint p0{0, 0, 0, 1};
+ mpoint p1{0, 0, 10, 1};
+ mpoint p2{0, 200, 20, 10};
+
+ std::vector<swc_record> swc{
+ {1, 3, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 3, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 1, p2.x, p2.y, p2.z, p2.radius, 2}
+ };
+ EXPECT_THROW(load_swc_neuron(swc), swc_no_soma);
+ }
+ {
+ // 3-point non-consecutive soma and a 2 point dendrite
+ mpoint p0{0, 0, 0, 1};
+ mpoint p1{0, 0, 10, 1};
+ mpoint p2{0, 0, 10, 10};
+ mpoint p3{0, 200, 20, 10};
+ mpoint p4{0, 0, 20, 1};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 3, p2.x, p2.y, p2.z, p2.radius, 2},
+ {4, 3, p3.x, p3.y, p3.z, p3.radius, 3},
+ {5, 1, p4.x, p4.y, p4.z, p4.radius, 2}
+ };
+ EXPECT_THROW(load_swc_neuron(swc), swc_non_consecutive_soma);
+ }
+ {
+ // 3-point soma and a 2 point dendrite connected in the middle of the soma
+ mpoint p0{0, 0, 0, 1};
+ mpoint p1{0, 0, 10, 1};
+ mpoint p2{0, 0, 20, 1};
+ mpoint p3{0, 0, 10, 10};
+ mpoint p4{0, 200, 20, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 1, p2.x, p2.y, p2.z, p2.radius, 2},
+ {4, 3, p3.x, p3.y, p3.z, p3.radius, 2},
+ {5, 3, p4.x, p4.y, p4.z, p4.radius, 4}
+ };
+ EXPECT_THROW(load_swc_neuron(swc), swc_branchy_soma);
+ }
+ {
+ // non-existent parent sample
+ mpoint p0{0, 0, 0, 1};
+ mpoint p1{0, 0, 10, 1};
+ mpoint p2{0, 200, 20, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 3, p2.x, p2.y, p2.z, p2.radius, 4}
+ };
+ EXPECT_THROW(load_swc_neuron(swc), swc_no_such_parent);
+ }
+ {
+ // parent sample is self
+ mpoint p0{0, 0, 0, 1};
+ mpoint p1{0, 0, 10, 1};
+ mpoint p2{0, 200, 20, 10};
+
+ std::vector<swc_record> swc{
+ {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+ {2, 1, p1.x, p1.y, p1.z, p1.radius, 1},
+ {3, 3, p2.x, p2.y, p2.z, p2.radius, 3}
+ };
+ EXPECT_THROW(load_swc_neuron(swc), swc_no_such_parent);
+ }
+}
// hipcc bug in reading DATADIR
#ifndef ARB_HIP