Section change (#510)

* A section within MorphIO is defined as a series of segments between
  bifurcation/multifurcation points.  However, SWC files allow section
  change without a branch. Normal parsing of this will raise an
  exception, but can be allowed using the option `UNIFURCATED_SECTION_CHANGE`
* Create a warning if the section type changes without a bifurcation
* Add option to allow for sections changes in SWC files without bifurcation

Author: mgeplf

binds/python/bind_enums.cpp

@@ -70,6 +70,8 @@ void bind_enums(py::module& m) {
         .value("soma_sphere", morphio::enums::Option::SOMA_SPHERE)
         .value("no_duplicates", morphio::enums::Option::NO_DUPLICATES)
         .value("nrn_order", morphio::enums::Option::NRN_ORDER)
+        .value("allow_unifurcated_section_change",
+               morphio::enums::Option::ALLOW_UNIFURCATED_SECTION_CHANGE)
         .export_values();
 
 
@@ -95,7 +97,8 @@ void bind_enums(py::module& m) {
         .value("write_empty_morphology", morphio::enums::WRITE_EMPTY_MORPHOLOGY)
         .value("zero_diameter", morphio::enums::Warning::ZERO_DIAMETER)
         .value("soma_non_contour", morphio::enums::Warning::SOMA_NON_CONTOUR)
-        .value("soma_non_cylinder_or_point", morphio::enums::Warning::SOMA_NON_CYLINDER_OR_POINT);
+        .value("soma_non_cylinder_or_point", morphio::enums::Warning::SOMA_NON_CYLINDER_OR_POINT)
+        .value("type_changed_within_section", morphio::enums::Warning::SECTION_TYPE_CHANGED);
 
     py::enum_<morphio::enums::SomaType>(m, "SomaType", py::arithmetic())
         .value("SOMA_UNDEFINED", morphio::enums::SomaType::SOMA_UNDEFINED)

binds/python/generated/docstrings.h

@@ -770,6 +770,18 @@ static const char *mkd_doc_morphio_SectionBuilderError = R"doc()doc";
 
 static const char *mkd_doc_morphio_SectionBuilderError_SectionBuilderError = R"doc()doc";
 
+static const char *mkd_doc_morphio_SectionTypeChanged = R"doc()doc";
+
+static const char *mkd_doc_morphio_SectionTypeChanged_SectionTypeChanged = R"doc()doc";
+
+static const char *mkd_doc_morphio_SectionTypeChanged_errorLevel = R"doc()doc";
+
+static const char *mkd_doc_morphio_SectionTypeChanged_lineNumber = R"doc()doc";
+
+static const char *mkd_doc_morphio_SectionTypeChanged_msg = R"doc()doc";
+
+static const char *mkd_doc_morphio_SectionTypeChanged_warning = R"doc()doc";
+
 static const char *mkd_doc_morphio_Section_Section = R"doc()doc";
 
 static const char *mkd_doc_morphio_Section_breadth_begin = R"doc(Breadth first iterator)doc";
@@ -1094,6 +1106,8 @@ static const char *mkd_doc_morphio_enums_Option =
 R"doc(The list of modifier flags that can be passed when loading a
 morphology See morphio::mut::modifiers for more information)doc";
 
+static const char *mkd_doc_morphio_enums_Option_ALLOW_UNIFURCATED_SECTION_CHANGE = R"doc(Allow section type to change without bifurcation)doc";
+
 static const char *mkd_doc_morphio_enums_Option_NO_DUPLICATES = R"doc(Skip duplicating points)doc";
 
 static const char *mkd_doc_morphio_enums_Option_NO_MODIFIER = R"doc(Read morphology as is without any modification)doc";
@@ -1211,6 +1225,8 @@ static const char *mkd_doc_morphio_enums_Warning_NO_SOMA_FOUND = R"doc(No soma f
 
 static const char *mkd_doc_morphio_enums_Warning_ONLY_CHILD = R"doc(Single child sections are not allowed in SWC format)doc";
 
+static const char *mkd_doc_morphio_enums_Warning_SECTION_TYPE_CHANGED = R"doc(In SWC, the type changed within a section, not post bifurcation)doc";
+
 static const char *mkd_doc_morphio_enums_Warning_SOMA_NON_CONFORM = R"doc(Soma does not conform the three point soma spec from NeuroMorpho.org)doc";
 
 static const char *mkd_doc_morphio_enums_Warning_SOMA_NON_CONTOUR = R"doc(Soma must be a contour for ASC and H5)doc";

doc/source/morphology.rst

@@ -172,6 +172,7 @@ The following flags are supported:
     each section is no longer the last point of the parent section.
 * ``morphio::NRN_ORDER``\: Neurite are reordered according to the
     `NEURON simulator ordering <https://github.com/neuronsimulator/nrn/blob/2dbf2ebf95f1f8e5a9f0565272c18b1c87b2e54c/share/lib/hoc/import3d/import3d_gui.hoc#L874>`_
+* ``morphio::UNIFURCATED_SECTION_CHANGE``\: Allow section type to change without bifurcation, emits warning
 
 Multiple flags can be passed by using the standard bit flag manipulation (works the same way in C++
 and Python):

include/morphio/enums.h

@@ -19,7 +19,8 @@ enum Option {
     TWO_POINTS_SECTIONS = 0x01,  //!< Read sections only with 2 or more points
     SOMA_SPHERE = 0x02,          //!< Interpret morphology soma as a sphere
     NO_DUPLICATES = 0x04,        //!< Skip duplicating points
-    NRN_ORDER = 0x08             //!< Order of neurites will be the same as in NEURON simulator
+    NRN_ORDER = 0x08,            //!< Order of neurites will be the same as in NEURON simulator
+    ALLOW_UNIFURCATED_SECTION_CHANGE = 0x10  //!< Allow section type to change without bifurcation
 };
 
 /**
@@ -42,6 +43,7 @@ enum Warning {
     ZERO_DIAMETER,               //!< Zero section diameter
     SOMA_NON_CONTOUR,            //!< Soma must be a contour for ASC and H5
     SOMA_NON_CYLINDER_OR_POINT,  //!< Soma must be stacked cylinders or a point
+    SECTION_TYPE_CHANGED,  //!< In SWC, the type changed within a section, not post bifurcation
 };
 
 enum AnnotationType {

include/morphio/warning_handling.h

@@ -54,6 +54,23 @@ struct ZeroDiameter: public WarningMessage {
     uint64_t lineNumber;
 };
 
+struct SectionTypeChanged: public WarningMessage {
+    SectionTypeChanged(std::string uri_, uint64_t lineNumber_)
+        : WarningMessage(std::move(uri_))
+        , lineNumber(lineNumber_) {}
+    morphio::enums::Warning warning() const final {
+        return Warning::SECTION_TYPE_CHANGED;
+    }
+    morphio::readers::ErrorLevel errorLevel = morphio::readers::ErrorLevel::WARNING;
+    std::string msg() const final {
+        static const char* description =
+            "Warning: Type changed within section, without bifurcation";
+        return "\n" + details::errorLink(uri, lineNumber, errorLevel) + description;
+    }
+
+    uint64_t lineNumber;
+};
+
 struct DisconnectedNeurite: public WarningMessage {
     DisconnectedNeurite(std::string uri_, uint64_t lineNumber_)
         : WarningMessage(std::move(uri_))

src/readers/morphologySWC.cpp

@@ -144,7 +144,6 @@ class SWCTokenizer
 struct SWCSample {
     enum : unsigned int { UNKNOWN_ID = 0xFFFFFFFE };
 
-    SWCSample() = default;  // XXX
     floatType diameter = -1.;
     Point point{};
     SectionType type = SECTION_UNDEFINED;
@@ -234,14 +233,15 @@ class SWCBuilder
     using Samples = std::vector<SWCSample>;
 
   public:
-    SWCBuilder(std::string path, WarningHandler* warning_handler)
+    SWCBuilder(std::string path, WarningHandler* warning_handler, unsigned int options)
         : path_(std::move(path))
-        , warning_handler_(warning_handler) {}
+        , warning_handler_(warning_handler)
+        , options_(options) {}
 
-    Property::Properties buildProperties(const std::string& contents, unsigned int options) {
+    Property::Properties buildProperties(const std::string& contents) {
         const Samples samples = readSamples(contents, path_);
         buildSWC(samples);
-        morph_.applyModifiers(options);
+        morph_.applyModifiers(options_);
         return morph_.buildReadOnly();
     }
 
@@ -305,10 +305,6 @@ class SWCBuilder
         for (const auto& s : soma_samples) {
             if (s.parentId == SWC_ROOT) {
                 parent_count++;
-            } else if (samples_.count(s.parentId) == 0) {
-                details::ErrorMessages err_(path_);
-                throw MissingParentError(
-                    err_.ERROR_MISSING_PARENT(s.id, static_cast<int>(s.parentId), s.lineNumber));
             } else if (samples_.at(s.parentId).type != SECTION_SOMA) {
                 details::ErrorMessages err_(path_);
                 throw SomaError(err_.ERROR_SOMA_WITH_NEURITE_PARENT(s.lineNumber));
@@ -477,7 +473,14 @@ class SWCBuilder
         while (children_count == 1) {
             sample = &samples_.at(id);
             if(sample->type != samples_.at(children_.at(id)[0]).type){
-                break;
+                if (options_ & ALLOW_UNIFURCATED_SECTION_CHANGE) {
+                    warning_handler_->emit(
+                        std::make_unique<SectionTypeChanged>(path_, sample->lineNumber));
+                    break;
+                }
+                throw RawDataError("Section type changed without a bifucation at line: " +
+                                   std::to_string(sample->lineNumber) +
+                                   ", consider using UNIFURCATED_SECTION_CHANGE option");
             }
             points.push_back(sample->point);
             diameters.push_back(sample->diameter);
@@ -518,9 +521,9 @@ class SWCBuilder
     mut::Morphology morph_;
     std::string path_;
     WarningHandler* warning_handler_;
+    unsigned int options_;
 };
 
-
 }  // namespace details
 
 namespace readers {
@@ -530,7 +533,7 @@ Property::Properties load(const std::string& path,
                           unsigned int options,
                           std::shared_ptr<WarningHandler>& warning_handler) {
     auto properties =
-        details::SWCBuilder(path, warning_handler.get()).buildProperties(contents, options);
+        details::SWCBuilder(path, warning_handler.get(), options).buildProperties(contents);
 
     properties._cellLevel._cellFamily = NEURON;
     properties._cellLevel._version = {"swc", 1, 0};

tests/test_1_swc.py

@@ -10,7 +10,7 @@
                    strip_color_codes)
 
 from morphio import (MorphioError, Morphology, RawDataError, SomaError,
-                     SomaType, Warning, ostream_redirect, set_raise_warnings)
+                     SomaType, Warning, ostream_redirect, set_raise_warnings, Option)
 
 
 DATA_DIR = Path(__file__).parent / "data"
@@ -568,13 +568,42 @@ def test_throw_on_negative_id():
         Morphology(content, extension='swc')
 
 
+def test_axon_carrying_dendrite():
+    contents =('''
+    1 1 0 0 1 1 -1
+    2 2 0 0 2 2 1
+    3 2 0 0 3 3 2
+
+    4 3 0 0 4 4 3  # dendrite splits off
+    5 3 0 0 5 5 4
+
+    6 2 0 0 6 6 3  # axon carries on
+    7 2 0 0 7 7 3
+               ''')
+    Morphology(contents, "swc")
+
+
 def test_multi_type_section():
+    """
+    A section within MorphIO is defined as a series of segments between
+    bifurcation/multifurcation points.  However, SWC files allow section
+    change without a branch. Normal parsing of this will raise an
+    exception, but can be allowed using the option `UNIFURCATED_SECTION_CHANGE`
+    """
     contents =('''1 1 0 4 0 3.0 -1
                   2 6 0 0 2 0.5 1        # <- type 6
                   3 7 0 0 3 0.5 2        # <- type 7
                   4 8 0 0 4 0.5 3        # <- type 8
                   5 9 0 0 5 0.5 4''')    # <- type 9
-    n = Morphology(contents, "swc")
+
+    with pytest.raises(RawDataError):
+        Morphology(contents, "swc")
+
+    warnings = morphio.WarningHandlerCollector()
+    n = Morphology(contents,
+                   "swc",
+                   warning_handler=warnings,
+                   options=Option.allow_unifurcated_section_change)
     assert_array_equal(n.soma.points, [[0, 4, 0]])
     assert_array_equal(n.soma.diameters, [6.0])
     assert len(n.root_sections) == 1
@@ -583,6 +612,10 @@ def test_multi_type_section():
                        np.array([[0, 0, 2], ]))
     assert len(n.sections) == 4
     assert_array_equal(n.section_offsets, [0, 1, 3, 5, 7])
+    warnings = [f.warning for f in warnings.get_all()]
+    assert len(warnings) == 3  # type 7, 8, and 9
+    for warning in warnings:
+        assert warning.warning() == Warning.type_changed_within_section
 
 
 def test_missing_parent():

tests/test_swc_reader.cpp

@@ -2,6 +2,7 @@
  *
  * SPDX-License-Identifier: Apache-2.0
  */
+#include <morphio/enums.h>
 #include <morphio/exceptions.h>
 #include <morphio/morphology.h>
 #include <morphio/section.h>
@@ -112,6 +113,27 @@ TEST_CASE("morphio::swc::errors") {
         )";
         CHECK_THROWS_AS(Morphology(multiple_soma, "swc"), SomaError);
     }
+    SECTION("multiple_soma") {
+        const auto* multiple_soma = R"(
+1 2 0 0 1 .5 -1
+2 1 0 0 1 .5 1
+        )";
+        CHECK_THROWS_AS(Morphology(multiple_soma, "swc"), SomaError);
+    }
+
+    SECTION("large_id") {
+        const auto* multiple_soma = R"(
+01234567890123456789 1 0 0 1 .5 -1
+        )";
+        CHECK_THROWS_AS(Morphology(multiple_soma, "swc"), RawDataError);
+    }
+
+    SECTION("large_parent_id") {
+        const auto* multiple_soma = R"(
+1 1 0 0 1 .5 01234567890123456789
+        )";
+        CHECK_THROWS_AS(Morphology(multiple_soma, "swc"), RawDataError);
+    }
 }
 
 TEST_CASE("morphio::swc::working") {
@@ -130,15 +152,29 @@ TEST_CASE("morphio::swc::working") {
     }
 
     SECTION("chimera-axon-on-dendrite") {
-        const auto* no_soma = R"(
+        const auto* aod = R"(
 1 1 0 0 1 1 -1
 2 2 0 0 2 2 1
 3 2 0 0 3 3 2
 4 3 0 0 4 4 3
-5 3 0 0 5 5 4
+5 3 0 0 5 5 3
         )";
-        const auto m = Morphology(no_soma, "swc");
-        REQUIRE(m.sections().size() == 2);
-        REQUIRE(m.diameters().size() == 5);
+        const auto m = Morphology(aod, "swc");
+        REQUIRE(m.sections().size() == 3);
+        REQUIRE(m.diameters().size() == 6);
+    }
+
+    SECTION("section_type_change") {
+        const auto* changes = R"(
+1 1 0 4 0 3.0 -1
+2 6 0 0 2 0.5 1        # <- type 6
+3 7 0 0 3 0.5 2        # <- type 7
+4 8 0 0 4 0.5 3        # <- type 8
+5 9 0 0 5 0.5 4        # <- type 9
+        )";
+        CHECK_THROWS_AS(Morphology(changes, "swc"), RawDataError);
+        const auto m =
+            Morphology(changes, "swc", morphio::Option::ALLOW_UNIFURCATED_SECTION_CHANGE);
+        REQUIRE(m.sections().size() == 4);
     }
 }