diff --git a/arborio/include/arborio/swcio.hpp b/arborio/include/arborio/swcio.hpp
index 2ec98a5e4ba792f7fe9bb6a78d5c356487346b71..2cb76b28693707562dea08013c48d1fb89a4d3ac 100644
--- a/arborio/include/arborio/swcio.hpp
+++ b/arborio/include/arborio/swcio.hpp
@@ -38,41 +38,6 @@ struct swc_spherical_soma: swc_error {
     explicit swc_spherical_soma(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);
@@ -83,16 +48,6 @@ 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 {
     int id = 0;          // sample number
     int tag = 0;         // structure identifier (tag)
@@ -177,10 +132,8 @@ arb::morphology load_swc_arbor(const swc_data& data);
 //
 // Note that 'one-point soma' SWC files are supported here
 //
-// These functions comply with inferred SWC rules from the Allen institute and Neuron.
-// These rules are explicitly listed in the docs.
+// Complies inferred SWC rules from NEURON, explicitly listed in the docs.
 
 arb::morphology load_swc_neuron(const swc_data& data);
-arb::morphology load_swc_allen(const swc_data& data, bool no_gaps=false);
 
 } // namespace arborio
diff --git a/arborio/neurolucida.cpp b/arborio/neurolucida.cpp
index f4b547d08bf5a4ca0457339dec18c29afc972689..dcb4cbadbe3c0bcecf4a25f38cd7de833ce9f574 100644
--- a/arborio/neurolucida.cpp
+++ b/arborio/neurolucida.cpp
@@ -1,5 +1,5 @@
 #include <cstring>
-#include <iostream>
+#include <ostream>
 #include <fstream>
 #include <numeric>
 
@@ -162,6 +162,7 @@ struct asc_color {
     uint8_t b;
 };
 
+[[maybe_unused]]
 std::ostream& operator<<(std::ostream& o, const asc_color& c) {
     return o << "(asc-color " << (int)c.r << " " << (int)c.g << " " << (int)c.b << ")";
 }
@@ -251,6 +252,7 @@ struct zsmear {
     double beta;
 };
 
+[[maybe_unused]]
 std::ostream& operator<<(std::ostream& o, const zsmear& z) {
     return o << "(zsmear " << z.alpha << " " << z.beta << ")";
 }
diff --git a/arborio/swcio.cpp b/arborio/swcio.cpp
index 69fa7833feb89dad57a60afce8bbafda56be76b5..1b6bec16f8f648d59abda2f78162e01ebb253d8a 100644
--- a/arborio/swcio.cpp
+++ b/arborio/swcio.cpp
@@ -1,8 +1,8 @@
 #include <cmath>
 #include <ios>
-#include <iostream>
 #include <limits>
 #include <numeric>
+#include <iostream>
 #include <set>
 #include <string>
 #include <sstream>
@@ -12,6 +12,8 @@
 
 #include <arbor/morph/segment_tree.hpp>
 
+#include "arbor/morph/primitives.hpp"
+
 #include <arborio/swcio.hpp>
 
 namespace arborio {
@@ -39,51 +41,14 @@ 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_error("Every record not attached to a soma sample 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_error("Only SWC record identifiers of 1, 2, 3 or 4 are supported.", 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) {
@@ -241,272 +206,128 @@ arb::morphology load_swc_arbor(const swc_data& data) {
 }
 
 arb::morphology load_swc_neuron(const swc_data& data) {
-    const auto& records = data.records();
-
-    // Assert that the file contains at least one sample.
-    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;
+    constexpr int soma_tag = 1;
 
-    // Map of SWC record id to index in `records`.
-    std::unordered_map<int, std::size_t> record_index;
-    record_index[soma_prox.id] = 0;
+    const auto n_samples = data.records().size();
 
-    // 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 (n_samples==0) {
+        return {};
+    }
 
-        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);
+    // The NEURON interpretation is only applied when the cell has a soma.
+    if (data.records()[0].tag != soma_tag) {
+        const auto& R = data.records();
+        // Search for other soma samples
+        if (auto it=std::find_if(R.begin(), R.end(), [](auto& r) {return r.tag==soma_tag;}); it!=R.end()) {
+            // The presence of a soma tag when there is a non-soma tag at the root
+            // violates the requirement that the parent of a soma sample is also a
+            // soma sample.
+            throw swc_mismatched_tags(it->id);
         }
 
-        // 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;
+        // No soma: fall back to Arbor interpretation.
+        // Check for tags unsupported by NEURON beforehand.
+        if (auto it=std::find_if(R.begin(), R.end(), [](const auto& r) {return r.tag<1 || r.tag>4;}); it!=R.end()) {
+            throw swc_unsupported_tag(R[std::distance(R.begin(), it)].id);
         }
-
-        prev_tag = r.tag;
-        prev_id = r.id;
+        return load_swc_arbor(data);
     }
 
-    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;
+    // Make a copy of the records and canonicalise them.
+    auto records = data.records();
+    std::unordered_map<int, int> record_index = {{-1, -1}};
+    std::vector<int> old_record_index(n_samples);
+
+    for (std::size_t i=0; i<n_samples; ++i) {
+        auto& r = records[i];
+        record_index[r.id] = i;
+        old_record_index[i] = r.id;
+        r.id = i;
+        if (!record_index.count(r.parent_id)) {
+            throw swc_no_such_parent(r.parent_id);
         }
-        // Model the soma as a 2 cylinders with total length=2*radius, extended along the y axis, centred 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;
+        r.parent_id = record_index[r.parent_id];
     }
-    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};
+    // Calculate meta-data
+    std::vector<int> child_count(n_samples);
+    std::size_t n_soma_samples = 0;
 
-        // 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);
+    for (std::size_t i=0; i<n_samples; ++i) {
+        auto& r = records[i];
+        // Only accept soma, axon, dend and apic samples.
+        if (!(r.tag>=0 && r.tag<=4)) {
+            throw swc_unsupported_tag(old_record_index[i]);
         }
-        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);
+        if (r.tag==soma_tag) {
+            ++n_soma_samples;
         }
-
-        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;
+        int pid = r.parent_id;
+        if (pid!=-1) {
+            ++child_count[pid];
+            const int ptag = records[pid].tag;
+            // Assert that sample has the same tag as its parent, or the parent is tagged soma.
+            if (r.tag!=ptag && ptag!=soma_tag) {
+                throw swc_mismatched_tags(old_record_index[i]);
+            }
         }
-
-        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 arb::morphology(tree);
-}
-
-arb::morphology load_swc_allen(const swc_data& data, bool no_gaps) {
-    auto records = data.records();
-
-    // 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;
+    const bool spherical_soma = n_soma_samples==1;
 
-    int soma_id = records[0].id;
-    record_index[soma_id] = 0;
+    // Construct the segment tree.
+    arb::segment_tree tree;
 
-    // Assert that root sample has tag 1.
-    if (records[0].tag != 1) {
-        throw swc_no_soma{records[0].id};
+    // It isn't possible to a-priori calculate the exact number of segments
+    // without more meta-data, but this will be accurate in the absence of
+    // single-sample sub-trees, which should be rare.
+    tree.reserve(n_samples+spherical_soma);
+
+    std::unordered_map<int, arb::msize_t> segmap;
+
+    // Construct a soma composed of two cylinders if is represented by a single sample.
+    if (spherical_soma) {
+        auto& sr = records[0];
+        arb::msize_t pid = arb::mnpos;
+        pid = tree.append(pid, {sr.x-sr.r, sr.y, sr.z, sr.r},
+                               {sr.x,      sr.y, sr.z, sr.r}, soma_tag);
+        // Children of the soma sample attach to the distal end of the first segment in the soma.
+        segmap[0] = pid;
+        pid = tree.append(pid, {sr.x,      sr.y, sr.z, sr.r},
+                               {sr.x+sr.r, sr.y, sr.z, sr.r}, soma_tag);
+    }
+    else {
+        segmap[0] = arb::mnpos;
     }
 
-    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);
+    // Build the tree.
+    for (std::size_t i=1; i<n_samples; ++i) {
+        auto& r = records[i];
+        const auto pid = r.parent_id;
+        auto& p = records[pid];
+        const int tag = r.tag;
 
-        if (parent_iter == record_index.end() || records[parent_iter->second].id == r.id)
-        {
-            throw swc_no_such_parent{r.id};
+        // Constructing a segment inside the soma or a sub-tree.
+        if (tag==p.tag) {
+            segmap[i] = tree.append(segmap.at(pid), {p.x, p.y, p.z, p.r}, {r.x, r.y, r.z, r.r}, tag);
         }
-
-        // 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};
+        // The start of a sub-tree.
+        else if (child_count[i]) {
+            // Do not create a segment, instead set up the segmap so that the
+            // first segment in the sub-tree will be connected to the soma with
+            // a "zero resistance cable".
+            segmap[i] = segmap.at(pid);
         }
-
-        // 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 centred 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;
+        // Sub-tree defined with a single sample.
+        else {
+            // The sub-tree is composed of a single segment connecting the soma
+            // to the sample, with constant radius defined by the sample.
+            segmap[i] = tree.append(segmap.at(pid), {p.x, p.y, p.z, r.r}, {r.x, r.y, r.z, r.r}, tag);
         }
-
-        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 arb::morphology(tree);
 }
+
 } // namespace arborio
 
diff --git a/doc/cpp/morphology.rst b/doc/cpp/morphology.rst
index f70352c7330aa49147fa0237ae937c5825f36f77..ae8fa4bca85afbb8b6930de7e2995f6c4de2cf4a 100644
--- a/doc/cpp/morphology.rst
+++ b/doc/cpp/morphology.rst
@@ -12,7 +12,7 @@ Morphology API
 
 .. todo::
 
-   TODO: Describe morphology methods.
+   Describe morphology methods.
 
 .. _cppcablecell-morphology-construction:
 
@@ -21,7 +21,7 @@ Constructing cell morphologies
 
 .. todo::
 
-   TODO: Description of segment trees.
+   Description of segment trees.
 
 
 The stitch-builder interface
@@ -354,7 +354,6 @@ basic checks performed on them. The :cpp:type:`swc_data` object can then be used
 :ref:`page <morph-formats>` for more details).
 
   * :cpp:func:`load_swc_arbor`
-  * :cpp:func:`load_swc_allen`
   * :cpp:func:`load_swc_neuron`
 
 .. cpp:class:: swc_record
@@ -403,18 +402,13 @@ basic checks performed on them. The :cpp:type:`swc_data` object can then be used
 
 .. cpp:function:: morphology load_swc_arbor(const swc_data& data)
 
-   Returns a :cpp:type:`morphology` constructed according to Arbor's SWC specifications.
-
-.. cpp:function:: morphology load_swc_allen(const swc_data& data, bool no_gaps=false)
-
-   Returns a :cpp:type:`morphology` constructed according to the Allen Institute's SWC
-   specifications. By default, gaps in the morphology are allowed, this can be toggled
-   using the ``no_gaps`` argument.
+   Returns a :cpp:type:`morphology` constructed according to Arbor's
+   :ref:`SWC specifications <formatswc-arbor>`.
 
 .. cpp:function:: morphology load_swc_neuron(const swc_data& data)
 
-   Returns a :cpp:type:`morphology` constructed according to NEURON's SWC specifications.
-
+   Returns a :cpp:type:`morphology` constructed according to NEURON's
+   :ref:`SWC specifications <formatswc-neuron>`.
 
 .. _cppasc:
 
@@ -424,8 +418,8 @@ Neurolucida ASCII
 Arbor supports reading morphologies described using the
 :ref:`Neurolucida ASCII file format <formatasc>`.
 
-The :cpp:func:`parse_asc()` function is used to parse the SWC file and generate a :cpp:type:`asc_morphology` object,
-which a simple struct with two members representing the morphology and a label dictionary with labeled
+The :cpp:func:`parse_asc()` function is used to parse the SWC file and generate a :cpp:type:`asc_morphology` object:
+a simple struct with two members representing the morphology and a label dictionary with labeled
 regions and locations.
 
 .. cpp:class:: asc_morphology
diff --git a/doc/fileformat/swc.rst b/doc/fileformat/swc.rst
index 953cc4e3a1f124f8a90bbc1e6e96270ba3e9bd63..65b19ade89ef80e234ded9fb0f6df388e987c53d 100644
--- a/doc/fileformat/swc.rst
+++ b/doc/fileformat/swc.rst
@@ -13,10 +13,11 @@ The description of the morphology is encoded as a list of samples with an id,
 an `x,y,z` location in space, a radius, a tag and a parent id. Arbor parses these samples, performs some checks,
 then generates a 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 interpretation that
+The SWC file format specifications does not describe how the file should be interpretted to reconstruct
+a morphology from SWC samples. This has lead different simulators to interpret SWC files in different
+ways, specifically the reconstruction of 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.
+interpret SWC files similarly to how the NEURON simulator would.
 
 Despite the differences between the interpretations, there is a common set of checks that are always performed
 to validate an SWC file:
@@ -30,9 +31,11 @@ its parent and inherits the tag of the sample; and if more than 1 sample have th
 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.
 
+.. _formatswc-arbor:
+
 Arbor interpretation
 """"""""""""""""""""
-In addition to the previously listed checks, the Arbor interpretation explicitly disallows SWC files where the soma is
+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, forming 1 or more segments. A *segment* is
 always constructed between a sample and its parent. This means that there are no gaps in the resulting morphology.
 
@@ -54,42 +57,60 @@ like this:
    :width: 400
    :align: center
 
+.. _formatswc-neuron:
 
-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.
+NEURON interpretation
+"""""""""""""""""""""
+Arbor provides support for interpreting SWC inputs in the same way as NEURON,
+to ease porting of cell models developed in NEURON to Arbor.
 
-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.
+The NEURON interpretations is based on the observed output of NEURON's ``Import3d_SWC_read``
+function, which is based on the `Neuromorpho approach <http://neuromorpho.org/SomaFormat.html>`_.
+However, there are differences and undocumented interpretations of the soma and how dendrites,
+axons and apical dendrites are attached to it that are not described explicitly by Neuromorpho.
 
-Finally the Allen institute interpretation of SWC files centres the morphology around the soma at the origin (0, 0, 0)
-and all samples are translated in space towards the origin.
+.. Warning::
+
+   The interpretation of SWC files by NEURON's import 3D method changed in NEURON
+   8 to address bugs in earlier versions. Arbor follows the NEURON 8 approach,
+   and can't guarantee compatibility with reconstructed SWC morphologies from NEURON 7.
+
+.. Note::
+
+    The rules below are applied to the morphology representation only when a soma
+    sample is present, otherwise the default
+    :ref:`Arbor interpretation <formatswc-arbor>` is applied.
+
+**Every sample must have the same SWC identifier (tag) as its parent, except for
+samples whose parent is tagged as soma**:
+This enforces that axons, dendrites and apical dendrites can only attach to the soma.
+Conversely, it isn't possible to attach an axon to a dendrite, for example.
+
+**The first sample is tagged as soma**:
+This requirement is a corollary of the previous rule.
+
+**Single-sample somas are permitted**:
+The `Neuromorpho guidelines <http://neuromorpho.org/SomaFormat.html>`_ regarding
+interpretation of a spherical soma described with a single soma sample can be summarised:
+
+* The soma is composed of two cylinders that have their proximal ends at the soma
+  center, extended first along the negative y-axis and then positive y-axis.
+
+Following the Neuromorpho specification, NEURON constructs the soma from two cylinders,
+joined at the soma center. It differs in two ways:
+
+* The soma is extended along the x-axis, not the y-axis.
+* The soma is constructed from three points, the first at ``x=x0-r``, the second with
+  ``x=x0`` and the third at ``x=x0+r``, to form a single section, with all dendrites, axons
+  and apical dendrites attached to the center of the soma with "zero resistance wires".
+
+**The axon, dendrite and apical sub-trees follow special rules for attachment to the soma**:
+By default, the sub-tree starts at the first sample with the dendrite, axon or apical tag, and not
+at the parent location on the soma, and the sub-tree is connected to its parent with a "zero resistance wire".
+**Except** when the sub tree is defined by a single child sample. In which case the sub-tree is
+composed of a single a segment from the parent location on the soma to the child sample,
+with constant radius of the child.
 
-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 enforce 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/python/morphology.rst b/doc/python/morphology.rst
index c8dfce3744e02093e8f4e8329159e56a22ebd45d..6bf92faa18acb575148fae08c89fcebb74176fdf 100644
--- a/doc/python/morphology.rst
+++ b/doc/python/morphology.rst
@@ -332,16 +332,16 @@ Cable cell morphology
 
     .. note::
         This method does not interpret the first sample, typically associated with the soma,
-        as a sphere. SWCs with single point somas are, unfortunately, reasonably common, for example
+        as a sphere. SWC files with single point somas are common, for example
         `SONATA <https://github.com/AllenInstitute/sonata/blob/master/docs/SONATA_DEVELOPER_GUIDE.md#representing-biophysical-neuron-morphologies>`_
         model descriptions.
 
-        Such representations are unfortunate because simulation tools like Arbor and NEURON require
-        the use of cylinders or fustrums to describe morphologies, and it is not possible to
-        infer how branches attached to the soma should be connected.
+        Such representations are challenging to consistently interpret in different
+        simulation tools because they require heuristics and, often undocumented, rules
+        for how to interpret the connectin of axons and dendrites to the soma.
 
-        The :func:`load_swc_allen` and :func:`load_swc_neuron` functions provide support for interpreting
-        such SWC files.
+        The :func:`load_swc_neuron` function provides support for loading
+        SWC files according to the interpretation used by NEURON.
 
 
     :param str filename: the name of the SWC file.
@@ -349,52 +349,13 @@ Cable cell morphology
 
 .. py:function:: load_swc_neuron(filename)
 
-    Loads the :class:`morphology` from an SWC file 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.
+    Loads the :class:`morphology` from an SWC file according to NEURON's ``Import3D``
+    interpretation of the SWC specification.
+    See :ref:`the SWC file documention <formatswc-neuron>` for more details.
 
     :param str filename: the name of the SWC file.
     :rtype: morphology
 
-
-.. py:function:: load_swc_allen(filename, no_gaps=False)
-
-    Loads the :class:`morphology` from an SWC file according to the AllenDB and Sonata's SWC specifications.
-    Specifically:
-
-        * The first sample (the root) is treated as the centre of the soma.
-        * The morphology is translated such that the soma is centred at (0,0,0).
-        * The first sample has tag 1 (soma).
-        * All other samples have tags 2, 3 or 4 (axon, apic and dend respectively)
-
-    SONATA prescribes that there should be no gaps, however some models in AllenDB
-    have gaps between the start of sections and the soma. The ``no_gaps`` flag can be
-    used to enforce this requirement.
-
-    Arbor does not support modelling the soma as a sphere, so a cylinder with length
-    equal to the soma diameter is used. The cylinder is centred on the origin, and
-    aligned along the z axis.
-    Axons and apical dendrites are attached to the proximal end of the cylinder, and
-    dendrites to the distal end, with a gap between the start of each branch and the
-    end of the soma cylinder to which it is attached.
-
-    :param str filename: the name of the SWC file.
-    :param bool no_gaps: enforce that distance between soma centre and branches attached to soma is the soma radius.
-    :rtype: morphology
-
 .. py:class:: place_pwlin
 
     A :class:`place_pwlin` object allows the querying of the 3-d location of locations and cables
diff --git a/python/morphology.cpp b/python/morphology.cpp
index 3d67b41ae2dde89c0669125336cfc2857a83ee86..24ffeb0235d3dbbed10759a395a251c319a07a77 100644
--- a/python/morphology.cpp
+++ b/python/morphology.cpp
@@ -260,49 +260,12 @@ void register_morphology(py::module& m) {
         "filename"_a,
         "Generate a morphology from an SWC file following the rules prescribed by Arbor.\n"
         "Specifically:\n"
-        "* Single-segment somas are disallowed. These are usually interpreted as spherical somas\n"
-        "  and are a special case. This behavior is not allowed using this SWC loader.\n"
+        "* Single-segment somas are disallowed.\n"
         "* There are no special rules related to somata. They can be one or multiple branches\n"
         "  and other segments can connect anywhere along them.\n"
         "* A segment is always created between a sample and its parent, meaning there\n"
         "  are no gaps in the resulting morphology.");
 
-    m.def("load_swc_allen",
-        [](std::string fname, bool no_gaps=false) {
-            std::ifstream fid{fname};
-            if (!fid.good()) {
-                throw pyarb_error(util::pprintf("can't open file '{}'", fname));
-            }
-            try {
-                auto data = arborio::parse_swc(fid);
-                check_trailing(fid, fname);
-                return arborio::load_swc_allen(data, no_gaps);
-
-            }
-            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()));
-            }
-        },
-        "filename"_a, "no_gaps"_a=false,
-        "Generate a morphology from an SWC file following the rules prescribed by AllenDB\n"
-        " and Sonata. Specifically:\n"
-        "* The first sample (the root) is treated as the centre of the soma.\n"
-        "* The first morphology is translated such that the soma is centred at (0,0,0).\n"
-        "* The first sample has tag 1 (soma).\n"
-        "* All other samples have tags 2, 3 or 4 (axon, apic and dend respectively)\n"
-        "SONATA prescribes that there should be no gaps, however the models in AllenDB\n"
-        "have gaps between the start of sections and the soma. The flag no_gaps can be\n"
-        "used to enforce this requirement.\n"
-        "\n"
-        "Arbor does not support modelling the soma as a sphere, so a cylinder with length\n"
-        "equal to the soma diameter is used. The cylinder is centred on the origin, and\n"
-        "aligned along the z axis.\n"
-        "Axons and apical dendrites are attached to the proximal end of the cylinder, and\n"
-        "dendrites to the distal end, with a gap between the start of each branch and the\n"
-        "end of the soma cylinder to which it is attached.");
-
     m.def("load_swc_neuron",
         [](std::string fname) {
             std::ifstream fid{fname};
@@ -322,23 +285,8 @@ void register_morphology(py::module& m) {
         },
         "filename"_a,
         "Generate a morphology from an SWC file following the rules prescribed by NEURON.\n"
-        " Specifically:\n"
-        "* The first sample must be a soma sample.\n"
-        "* The soma is represented by a series of n≥1 unbranched, serially listed samples.\n"
-        "* The soma is constructed as a single cylinder with diameter equal to the piecewise\n"
-        "  average diameter of all the segments forming the soma.\n"
-        "* A single-sample soma at is constructed as a cylinder with length=diameter.\n"
-        "* If a non-soma sample is to have a soma sample as its parent, it must have the\n"
-        "  most distal sample of the soma as the parent.\n"
-        "* Every non-soma sample that has a soma sample as its parent, attaches to the\n"
-        "  created soma cylinder at its midpoint.\n"
-        "* If a non-soma sample has a soma sample as its parent, no segment is created\n"
-        "  between the sample and its parent, instead that sample is the proximal point of\n"
-        "  a new segment, and there is a gap in the morphology (represented electrically as a\n"
-        "  zero-resistance wire)\n"
-        "* To create a segment with a certain tag, that is to be attached to the soma,\n"
-        "  we need at least 2 samples with that tag."
-        );
+        "See the documentation https://arbor.readthedocs.io/en/latest/fileformat/swc.html\n"
+        "for a detailed description of the interpretation.");
 
     // arb::morphology
 
diff --git a/test/unit/test_swcio.cpp b/test/unit/test_swcio.cpp
index ea36cc005c2a72148f771b2c8a7fc0a9c5639404..d73c61988c4488c14cd28ef52899eb35937273c7 100644
--- a/test/unit/test_swcio.cpp
+++ b/test/unit/test_swcio.cpp
@@ -4,8 +4,8 @@
 #include <sstream>
 
 #include <arbor/cable_cell.hpp>
+#include <arbor/morph/morphology.hpp>
 #include <arbor/morph/primitives.hpp>
-#include <arbor/morph/segment_tree.hpp>
 
 #include <arborio/swcio.hpp>
 
@@ -18,7 +18,6 @@
 #endif
 
 using namespace arborio;
-using arb::segment_tree;
 using arb::mpoint;
 using arb::mnpos;
 
@@ -254,7 +253,7 @@ TEST(swc_parser, stream_validity) {
 
 }
 
-TEST(swc_parser, arbor_complaint) {
+TEST(swc_parser, arbor_compliant) {
     {
         // Otherwise, ensure segment ends and tags correspond.
         mpoint p0{0.1, 0.2, 0.3, 0.4};
@@ -285,7 +284,7 @@ TEST(swc_parser, arbor_complaint) {
         EXPECT_EQ(1u, segs_0.size());
         EXPECT_EQ(2u, segs_1.size());
         EXPECT_EQ(1u, segs_2.size());
-        
+
         EXPECT_EQ(1,  segs_0[0].tag);
         EXPECT_EQ(p0, segs_0[0].prox);
         EXPECT_EQ(p1, segs_0[0].dist);
@@ -342,7 +341,7 @@ TEST(swc_parser, arbor_complaint) {
     }
 }
 
-TEST(swc_parser, not_arbor_complaint) {
+TEST(swc_parser, not_arbor_compliant) {
     {
         // Missing parent record will throw.
         std::vector<swc_record> swc{
@@ -369,194 +368,6 @@ TEST(swc_parser, not_arbor_complaint) {
     }
 }
 
-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}
-        };
-        auto morpho = 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, morpho.num_branches());
-        EXPECT_EQ(mnpos, morpho.branch_parent(0));
-
-        auto segs_0 = morpho.branch_segments(0);
-
-        EXPECT_EQ(1u, segs_0.size());
-
-        EXPECT_EQ(1,     segs_0[0].tag);
-        EXPECT_EQ(prox,  segs_0[0].prox);
-        EXPECT_EQ(dist,  segs_0[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}
-        };
-        auto morpho = load_swc_allen(swc);
-
-        mpoint prox{0, -10, 0, 10};
-        mpoint dist{0,  10, 0, 10};
-
-        ASSERT_EQ(1u, morpho.num_branches());
-
-        EXPECT_EQ(mnpos, morpho.branch_parent(0));
-
-        auto segs_0 = morpho.branch_segments(0);
-
-        EXPECT_EQ(2u, segs_0.size());
-
-        EXPECT_EQ(1,    segs_0[0].tag);
-        EXPECT_EQ(prox, segs_0[0].prox);
-        EXPECT_EQ(dist, segs_0[0].dist);
-
-        EXPECT_EQ(3,    segs_0[1].tag);
-        EXPECT_EQ(p1,   segs_0[1].prox);
-        EXPECT_EQ(p2,   segs_0[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}
-        };
-        auto morpho = load_swc_allen(swc);
-
-        mpoint prox{0, -1, 0, 1};
-        mpoint dist{0,  1, 0, 1};
-
-        ASSERT_EQ(2u, morpho.num_branches());
-
-        EXPECT_EQ(mnpos, morpho.branch_parent(0));
-        EXPECT_EQ(mnpos, morpho.branch_parent(1));
-
-        auto segs_0 = morpho.branch_segments(0);
-        auto segs_1 = morpho.branch_segments(1);
-
-        EXPECT_EQ(2u, segs_0.size());
-        EXPECT_EQ(1u, segs_1.size());
-
-        EXPECT_EQ(1,     segs_0[0].tag);
-        EXPECT_EQ(prox,  segs_0[0].prox);
-        EXPECT_EQ(dist,  segs_0[0].dist);
-
-        EXPECT_EQ(3,   segs_0[1].tag);
-        EXPECT_EQ(p1,  segs_0[1].prox);
-        EXPECT_EQ(p2,  segs_0[1].dist);
-
-        EXPECT_EQ(2,   segs_1[0].tag);
-        EXPECT_EQ(p3,  segs_1[0].prox);
-        EXPECT_EQ(p4,  segs_1[0].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_record_precedes_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_record_precedes_parent);
-    }
-}
-
 TEST(swc_parser, neuron_compliant) {
     using namespace arborio;
     {
@@ -568,20 +379,23 @@ TEST(swc_parser, neuron_compliant) {
         };
         auto morpho = 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};
+        mpoint prox{p0.x-p0.radius, p0.y, p0.z, p0.radius};
+        mpoint dist{p0.x+p0.radius, p0.y, p0.z, p0.radius};
 
         ASSERT_EQ(1u, morpho.num_branches());
 
         EXPECT_EQ(mnpos, morpho.branch_parent(0));
 
-        auto segs_0 = morpho.branch_segments(0);
+        auto segs = morpho.branch_segments(0);
 
-        EXPECT_EQ(1u, segs_0.size());
+        EXPECT_EQ(2u, segs.size());
 
-        EXPECT_EQ(1,     segs_0[0].tag);
-        EXPECT_EQ(prox,  segs_0[0].prox);
-        EXPECT_EQ(dist,  segs_0[0].dist);
+        EXPECT_EQ(1,     segs[0].tag);
+        EXPECT_EQ(prox,  segs[0].prox);
+        EXPECT_EQ(p0,    segs[0].dist);
+        EXPECT_EQ(1,     segs[1].tag);
+        EXPECT_EQ(p0,    segs[1].prox);
+        EXPECT_EQ(dist,  segs[1].dist);
     }
     {
         // Two-point soma; interpretted as 1 segment
@@ -645,12 +459,12 @@ TEST(swc_parser, neuron_compliant) {
         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}
+            { 2, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+            { 4, 1, p1.x, p1.y, p1.z, p1.radius,  2},
+            { 6, 1, p2.x, p2.y, p2.z, p2.radius,  4},
+            { 8, 1, p3.x, p3.y, p3.z, p3.radius,  6},
+            {10, 1, p4.x, p4.y, p4.z, p4.radius,  8},
+            {12, 1, p5.x, p5.y, p5.z, p5.radius, 10}
         };
         auto morpho = load_swc_neuron(swc);
 
@@ -685,8 +499,8 @@ TEST(swc_parser, neuron_compliant) {
     {
         // One-point soma, two-point dendrite
         mpoint p0{0,   0, 0, 10};
-        mpoint p1{0,   0, 0,  5};
-        mpoint p2{0, 200, 0, 10};
+        mpoint p1{0,  10, 0,  2};
+        mpoint p2{0, 200, 0,  1};
 
         std::vector<swc_record> swc{
             {1, 1, p0.x, p0.y, p0.z, p0.radius, -1},
@@ -695,8 +509,8 @@ TEST(swc_parser, neuron_compliant) {
         };
         auto morpho = load_swc_neuron(swc);
 
-        mpoint prox{0, -10, 0, 10};
-        mpoint dist{0,  10, 0, 10};
+        mpoint prox{-10, 0, 0, 10};
+        mpoint dist{ 10, 0, 0, 10};
 
         ASSERT_EQ(3u, morpho.num_branches());
 
@@ -725,29 +539,19 @@ TEST(swc_parser, neuron_compliant) {
         EXPECT_EQ(p2, segs_2[0].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};
+        // One-point soma, one-point dendrite
+        mpoint p0{0,   0, 0, 10};
+        mpoint p1{0, 200, 0,  1};
 
+        // Use some non-contiguous record ids to test normalisation.
         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}
+            {23, 1, p0.x, p0.y, p0.z, p0.radius, -1},
+            {83, 3, p1.x, p1.y, p1.z, p1.radius, 23}
         };
         auto morpho = load_swc_neuron(swc);
 
-        mpoint mid {0, 0, 5, 2.25};
+        mpoint prox{-10, 0, 0, 10};
+        mpoint dist{ 10, 0, 0, 10};
 
         ASSERT_EQ(3u, morpho.num_branches());
 
@@ -759,37 +563,21 @@ TEST(swc_parser, neuron_compliant) {
         auto segs_1 = morpho.branch_segments(1);
         auto segs_2 = morpho.branch_segments(2);
 
-        EXPECT_EQ(3u, segs_0.size());
-        EXPECT_EQ(3u, segs_1.size());
+        EXPECT_EQ(1u, segs_0.size());
+        EXPECT_EQ(1u, segs_1.size());
         EXPECT_EQ(1u, segs_2.size());
 
-        EXPECT_EQ(1,  segs_0[0].tag);
-        EXPECT_EQ(p0, segs_0[0].prox);
-        EXPECT_EQ(p1, segs_0[0].dist);
-
-        EXPECT_EQ(1,  segs_0[1].tag);
-        EXPECT_EQ(p1, segs_0[1].prox);
-        EXPECT_EQ(p2, segs_0[1].dist);
-
-        EXPECT_EQ(1,  segs_0[2].tag);
-        EXPECT_EQ(p2, segs_0[2].prox);
-        EXPECT_EQ(mid,segs_0[2].dist);
-
-        EXPECT_EQ(1,  segs_1[0].tag);
-        EXPECT_EQ(mid,segs_1[0].prox);
-        EXPECT_EQ(p3, segs_1[0].dist);
-
-        EXPECT_EQ(1,  segs_1[1].tag);
-        EXPECT_EQ(p3, segs_1[1].prox);
-        EXPECT_EQ(p4, segs_1[1].dist);
+        EXPECT_EQ(1,    segs_0[0].tag);
+        EXPECT_EQ(prox, segs_0[0].prox);
+        EXPECT_EQ(p0,   segs_0[0].dist);
 
-        EXPECT_EQ(1,  segs_1[2].tag);
-        EXPECT_EQ(p4, segs_1[2].prox);
-        EXPECT_EQ(p5, segs_1[2].dist);
+        EXPECT_EQ(1,    segs_1[0].tag);
+        EXPECT_EQ(p0,   segs_1[0].prox);
+        EXPECT_EQ(dist, segs_1[0].dist);
 
         EXPECT_EQ(3,  segs_2[0].tag);
-        EXPECT_EQ(p6, segs_2[0].prox);
-        EXPECT_EQ(p7, segs_2[0].dist);
+        EXPECT_EQ((mpoint{0, 0, 0, 1}), segs_2[0].prox);
+        EXPECT_EQ(p1, segs_2[0].dist);
     }
     {
         // Two-point soma, two-point dendrite
@@ -806,296 +594,192 @@ TEST(swc_parser, neuron_compliant) {
         };
         auto morpho = load_swc_neuron(swc);
 
-        mpoint mid{0, 0, -10, 7};
-
-        ASSERT_EQ(3u, morpho.num_branches());
+        ASSERT_EQ(1u, morpho.num_branches());
 
         EXPECT_EQ(mnpos, morpho.branch_parent(0));
-        EXPECT_EQ(0u,    morpho.branch_parent(1));
-        EXPECT_EQ(0u,    morpho.branch_parent(2));
 
         auto segs_0 = morpho.branch_segments(0);
-        auto segs_1 = morpho.branch_segments(1);
-        auto segs_2 = morpho.branch_segments(2);
 
-        EXPECT_EQ(1u, segs_0.size());
-        EXPECT_EQ(1u, segs_1.size());
-        EXPECT_EQ(1u, segs_2.size());
+        EXPECT_EQ(2u, segs_0.size());
 
         EXPECT_EQ(1,   segs_0[0].tag);
         EXPECT_EQ(p0,  segs_0[0].prox);
-        EXPECT_EQ(mid, segs_0[0].dist);
-
-        EXPECT_EQ(1,   segs_1[0].tag);
-        EXPECT_EQ(mid, segs_1[0].prox);
-        EXPECT_EQ(p1,  segs_1[0].dist);
+        EXPECT_EQ(p1,  segs_0[0].dist);
 
-        EXPECT_EQ(3,   segs_2[0].tag);
-        EXPECT_EQ(p2,  segs_2[0].prox);
-        EXPECT_EQ(p3,  segs_2[0].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}
+        EXPECT_EQ(3,   segs_0[1].tag);
+        EXPECT_EQ(p2,  segs_0[1].prox);
+        EXPECT_EQ(p3,  segs_0[1].dist);
+    }
+    {
+        // A more elaborate example with multiple dendrites attached to different
+        // parts of a 7-sample branchy soma with branches.
+        // Some of the dendrites are defined by a single sample, and others by 2
+        // samples, so this is a pretty thorough stress test for connections to the soma.
+        //
+        // Values taken from a test/pynrn/test_swc.py test case.
+
+        mpoint p0 {  0,  0, 0, 5};
+        mpoint p1 { 10,  0, 0, 3};
+        mpoint p2 { 15,  0, 0, 4};
+        mpoint p3 { 20,  0, 0, 2};
+        mpoint p4 { -8, -6, 0, 3};
+        mpoint p5 { -8,  6, 0, 4};
+        mpoint p6 {-16, 12, 0, 2};
+        mpoint p7 { 30,  0, 0, 1};
+        mpoint p8 {-16,-12, 0, 1};
+        mpoint p9 {-24, 18, 0, 1};
+        mpoint p10{  0,  5, 0, 1};
+        mpoint p11{  0, 15, 0, 1};
+        mpoint p12{ 10,  3, 0, 1};
+        mpoint p13{ 10, 15, 0, 1};
+        mpoint p14{ 15,  4, 0, 1};
+        mpoint p15{ 15, 15, 0, 1};
+        mpoint p16{ -8, 10, 0, 1};
+        mpoint p17{ -8, 15, 0, 1};
+
+        std::vector<swc_record> swc{
+            { 1, 1,  0,  0, 0, 5, -1},
+            { 2, 1, 10,  0, 0, 3,  1},
+            { 3, 1, 15,  0, 0, 4,  2},
+            { 4, 1, 20,  0, 0, 2,  3},
+            { 5, 1, -8, -6, 0, 3,  1},
+            { 6, 1, -8,  6, 0, 4,  1},
+            { 7, 1,-16, 12, 0, 2,  6},
+            { 8, 3, 30,  0, 0, 1,  4},
+            { 9, 3,-16,-12, 0, 1,  5},
+            {10, 3,-24, 18, 0, 1,  7},
+            {11, 3,  0,  5, 0, 1,  1},
+            {12, 3,  0, 15, 0, 1, 11},
+            {13, 3, 10,  3, 0, 1,  2},
+            {14, 3, 10, 15, 0, 1, 13},
+            {15, 3, 15,  4, 0, 1,  3},
+            {16, 3, 15, 15, 0, 1, 15},
+            {17, 3, -8, 10, 0, 1,  6},
+            {18, 3, -8, 15, 0, 1, 17}
         };
-        auto morpho = load_swc_neuron(swc);
 
-        mpoint mid{0, 0, -7.5, 6.5};
+        auto morpho = load_swc_neuron(swc);
 
-        ASSERT_EQ(4u, morpho.num_branches());
+        ASSERT_EQ(10u, morpho.num_branches());
 
         EXPECT_EQ(mnpos, morpho.branch_parent(0));
         EXPECT_EQ(0u,    morpho.branch_parent(1));
-        EXPECT_EQ(0u,    morpho.branch_parent(2));
-        EXPECT_EQ(0u,    morpho.branch_parent(3));
+        EXPECT_EQ(1u,    morpho.branch_parent(2));
+        EXPECT_EQ(mnpos, morpho.branch_parent(3));
+        EXPECT_EQ(mnpos, morpho.branch_parent(4));
+        EXPECT_EQ(4u,    morpho.branch_parent(5));
+        EXPECT_EQ(mnpos, morpho.branch_parent(6));
+        EXPECT_EQ(0u,    morpho.branch_parent(7));
+        EXPECT_EQ(1u,    morpho.branch_parent(8));
+        EXPECT_EQ(4u,    morpho.branch_parent(9));
 
         auto segs_0 = morpho.branch_segments(0);
         auto segs_1 = morpho.branch_segments(1);
         auto segs_2 = morpho.branch_segments(2);
         auto segs_3 = morpho.branch_segments(3);
-
-        EXPECT_EQ(1u, segs_0.size());
-        EXPECT_EQ(1u, segs_1.size());
-        EXPECT_EQ(1u, segs_2.size());
-        EXPECT_EQ(1u, segs_3.size());
+        auto segs_4 = morpho.branch_segments(4);
+        auto segs_5 = morpho.branch_segments(5);
+        auto segs_6 = morpho.branch_segments(6);
+        auto segs_7 = morpho.branch_segments(7);
+        auto segs_8 = morpho.branch_segments(8);
+        auto segs_9 = morpho.branch_segments(9);
+
+        for (auto i: {0, 1, 4, 6, 7, 8, 9})
+            EXPECT_EQ(1u, morpho.branch_segments(i).size());
+        for (auto i: {2, 3, 5})
+            EXPECT_EQ(2u, morpho.branch_segments(i).size());
 
         EXPECT_EQ(1,   segs_0[0].tag);
         EXPECT_EQ(p0,  segs_0[0].prox);
-        EXPECT_EQ(mid, segs_0[0].dist);
+        EXPECT_EQ(p1,  segs_0[0].dist);
 
         EXPECT_EQ(1,   segs_1[0].tag);
-        EXPECT_EQ(mid, segs_1[0].prox);
-        EXPECT_EQ(p1,  segs_1[0].dist);
+        EXPECT_EQ(p1,  segs_1[0].prox);
+        EXPECT_EQ(p2,  segs_1[0].dist);
 
-        EXPECT_EQ(3,   segs_2[0].tag);
+        EXPECT_EQ(1,   segs_2[0].tag);
         EXPECT_EQ(p2,  segs_2[0].prox);
         EXPECT_EQ(p3,  segs_2[0].dist);
+        EXPECT_EQ(3,   segs_2[1].tag);
+        EXPECT_EQ((mpoint{p3.x, p3.y, p3.z, p7.radius}),  segs_2[1].prox);
+        EXPECT_EQ(p7,  segs_2[1].dist);
 
-        EXPECT_EQ(2,   segs_3[0].tag);
-        EXPECT_EQ(p2,  segs_3[0].prox);
+        EXPECT_EQ(1,   segs_3[0].tag);
+        EXPECT_EQ(p0,  segs_3[0].prox);
         EXPECT_EQ(p4,  segs_3[0].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}
-        };
-        auto morpho = load_swc_neuron(swc);
-
-        mpoint mid{0, 0, 4.5, 1.5};
-
-        ASSERT_EQ(3u, morpho.num_branches());
-
-        EXPECT_EQ(mnpos, morpho.branch_parent(0));
-        EXPECT_EQ(0u,    morpho.branch_parent(1));
-        EXPECT_EQ(0u,    morpho.branch_parent(2));
-
-        auto segs_0 = morpho.branch_segments(0);
-        auto segs_1 = morpho.branch_segments(1);
-        auto segs_2 = morpho.branch_segments(2);
+        EXPECT_EQ(3,   segs_3[1].tag);
+        EXPECT_EQ((mpoint{p4.x, p4.y, p4.z, p8.radius}),  segs_3[1].prox);
+        EXPECT_EQ(p8,  segs_3[1].dist);
 
-        EXPECT_EQ(1u, segs_0.size());
-        EXPECT_EQ(1u, segs_1.size());
-        EXPECT_EQ(3u, segs_2.size());
+        EXPECT_EQ(1,   segs_4[0].tag);
+        EXPECT_EQ(p0,  segs_4[0].prox);
+        EXPECT_EQ(p5,  segs_4[0].dist);
 
-        EXPECT_EQ(1,   segs_0[0].tag);
-        EXPECT_EQ(p0,  segs_0[0].prox);
-        EXPECT_EQ(mid, segs_0[0].dist);
+        EXPECT_EQ(1,   segs_5[0].tag);
+        EXPECT_EQ(p5,  segs_5[0].prox);
+        EXPECT_EQ(p6,  segs_5[0].dist);
+        EXPECT_EQ(3,   segs_5[1].tag);
+        EXPECT_EQ((mpoint{p6.x, p6.y, p6.z, p9.radius}),  segs_5[1].prox);
+        EXPECT_EQ(p9,  segs_5[1].dist);
 
-        EXPECT_EQ(1,   segs_1[0].tag);
-        EXPECT_EQ(mid, segs_1[0].prox);
-        EXPECT_EQ(p1,  segs_1[0].dist);
+        EXPECT_EQ(3,   segs_6[0].tag);
+        EXPECT_EQ(p10, segs_6[0].prox);
+        EXPECT_EQ(p11, segs_6[0].dist);
 
-        EXPECT_EQ(3,   segs_2[0].tag);
-        EXPECT_EQ(p2,  segs_2[0].prox);
-        EXPECT_EQ(p3,  segs_2[0].dist);
+        EXPECT_EQ(3,   segs_7[0].tag);
+        EXPECT_EQ(p12, segs_7[0].prox);
+        EXPECT_EQ(p13, segs_7[0].dist);
 
-        EXPECT_EQ(2,   segs_2[1].tag);
-        EXPECT_EQ(p3,  segs_2[1].prox);
-        EXPECT_EQ(p4,  segs_2[1].dist);
+        EXPECT_EQ(3,   segs_8[0].tag);
+        EXPECT_EQ(p14, segs_8[0].prox);
+        EXPECT_EQ(p15, segs_8[0].dist);
 
-        EXPECT_EQ(2,   segs_2[2].tag);
-        EXPECT_EQ(p4,  segs_2[2].prox);
-        EXPECT_EQ(p5,  segs_2[2].dist);
+        EXPECT_EQ(3,   segs_9[0].tag);
+        EXPECT_EQ(p16, segs_9[0].prox);
+        EXPECT_EQ(p17, segs_9[0].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};
-
+        // The first soma sample is not the first sample.
         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}
+            {1, 3, 0, 0, 0, 2, -1}, // first sample is tagged as dendrite
+            {2, 1, 5, 0, 0, 2,  1}  // attached ot it is a soma sample (error)
         };
 
-        EXPECT_THROW(load_swc_neuron(swc), swc_collocated_soma);
+        EXPECT_THROW(load_swc_neuron(swc), swc_mismatched_tags);
+        try { load_swc_neuron(swc); }
+        catch (swc_mismatched_tags& e) {EXPECT_EQ(2, e.record_id);}
     }
     {
-        // 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};
-
+        // A dendrite attached to an axon instead of the soma.
         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}
+            {1, 1, 0, 0, 0, 2, -1}, // two sample soma
+            {2, 1, 5, 0, 0, 2,  1},
+            {3, 2, 5, 0, 0, 1,  2}, // with an axon attached
+            {4, 2,10, 0, 0, 1,  3},
+            {5, 3,20, 0, 0, 1,  4}, // and a dendrite attached to that (error)
         };
-        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);
+        EXPECT_THROW(load_swc_neuron(swc), swc_mismatched_tags);
+        try { load_swc_neuron(swc); }
+        catch (swc_mismatched_tags& e) {EXPECT_EQ(5, e.record_id);}
     }
     {
-        // 1-point soma and 1-point dendrite
-        mpoint p0{0,   0, 0, 10};
-        mpoint p1{0, 200, 0, 10};
-
+        // A tag other than 1,2,3,4
         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}
+            {1, 1, 0, 0, 0, 2, -1}, // two sample soma
+            {2, 1, 5, 0, 0, 2,  1},
+            {3, 2, 5, 0, 0, 1,  2}, // with an axon attached
+            {4, 2,10, 0, 0, 1,  3},
+            {5, 5, 5, 0, 0, 1,  2}, // and another branch with tag 5 attached (error)
+            {6, 5, 5, 5, 0, 1,  5},
         };
-        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_record_precedes_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_record_precedes_parent);
+        EXPECT_THROW(load_swc_neuron(swc), swc_unsupported_tag);
+        try { load_swc_neuron(swc); }
+        catch (swc_unsupported_tag& e) {EXPECT_EQ(5, e.record_id);}
     }
 }