make sure SWC reader follows the MorphIO invariant

src/mut/writer_swc.cpp

@@ -24,7 +24,23 @@
 #include "../shared_utils.hpp"
 #include "writer_utils.h"
 
+
 namespace {
+bool anySamplesMatch(const morphio::Point& point,
+                     const morphio::floatType diameter,
+                     const morphio::Points& points,
+                     const std::vector<morphio::floatType>& diameters) {
+    for (size_t i = 0; i < points.size(); ++i) {
+        if (std::fabs(diameters[i] - diameter) < morphio::epsilon &&
+            std::fabs(points[0][0] - point[0]) < morphio::epsilon &&
+            std::fabs(points[0][1] - point[1]) < morphio::epsilon &&
+            std::fabs(points[0][2] - point[2]) < morphio::epsilon) {
+            return true;
+        }
+    }
+    return false;
+}
+
 void writeLine(std::ofstream& myfile,
                int id,
                int parentId,
@@ -95,8 +111,7 @@ int writeSoma(std::ofstream& myfile,
     return startIdOnDisk;
 }
 
-/* Only skip duplicate if it has the same diameter */
-bool _skipDuplicate(const std::shared_ptr<morphio::mut::Section>& section) {
+bool _sameDiameter(const std::shared_ptr<morphio::mut::Section>& section) {
     return std::fabs(section->diameters().front() - section->parent()->diameters().back()) <
            morphio::epsilon;
 }
@@ -154,6 +169,7 @@ void swc(const Morphology& morph,
 
     std::ofstream myfile(filename);
     writeHeader(myfile);
+
     int segmentIdOnDisk = writeSoma(myfile, soma, handler);
 
     std::unordered_map<uint32_t, int32_t> newIds;
@@ -173,13 +189,31 @@ void swc(const Morphology& morph,
             }
         }
 
-        // skips duplicate point for non-root sections
-        unsigned int firstPoint = ((isRootSection || !_skipDuplicate(section)) ? 0 : 1);
+        unsigned int firstPoint = 0;
+        // When there is only a single point in a section on read, MorphIO adds a point to keep the
+        // invariant that each section has a segment, which in turn has 2 points. When writing, that
+        // point is dropped
+        if (isRootSection &&
+            anySamplesMatch(points[0], diameters[0], soma->points(), soma->diameters())) {
+            firstPoint = 1;
+        } else if (!isRootSection && _sameDiameter(section)) {
+            // skips duplicate point for non-root sections, if they have the same diameter
+            firstPoint = 1;
+        }
+
         for (unsigned int i = firstPoint; i < points.size(); ++i) {
-            int parentIdOnDisk = (i > firstPoint)
-                                     ? segmentIdOnDisk - 1
-                                     : (isRootSection ? (soma->points().empty() ? -1 : 1)
-                                                      : newIds[section->parent()->id()]);
+            int parentIdOnDisk = 0;
+            if (i > firstPoint) {
+                parentIdOnDisk = segmentIdOnDisk - 1;
+            } else if (isRootSection) {
+                if (soma->points().empty()) {
+                    parentIdOnDisk = -1;
+                } else {
+                    parentIdOnDisk = 1;
+                }
+            } else {
+                parentIdOnDisk = newIds[section->parent()->id()];
+            }
 
             writeLine(
                 myfile, segmentIdOnDisk, parentIdOnDisk, section->type(), points[i], diameters[i]);

src/readers/morphologySWC.cpp

@@ -478,7 +478,7 @@ class SWCBuilder
                         std::make_unique<SectionTypeChanged>(path_, sample->lineNumber));
                     break;
                 }
-                throw RawDataError("Section type changed without a bifucation at line: " +
+                throw RawDataError("Section type changed without a bifurcation at line: " +
                                    std::to_string(sample->lineNumber) +
                                    ", consider using UNIFURCATED_SECTION_CHANGE option");
             }
@@ -491,6 +491,14 @@ class SWCBuilder
         sample = &samples_.at(id);
         points.push_back(sample->point);
         diameters.push_back(sample->diameter);
+
+        if (points.size() == 1) {
+            // To maintain the invariant that each section has a segment, which in turn has 2 points
+            // a point is added
+            diameters.insert(diameters.begin(), start_diameter);
+            points.insert(points.begin(), start_point);
+        }
+
         appendSection(DeclaredID(id), parent_id, sample->type);
 
         if (children_count == 0) {

src/shared_utils.cpp

@@ -2,14 +2,16 @@
  *
  * SPDX-License-Identifier: Apache-2.0
  */
+#include "shared_utils.hpp"
+
 #include <algorithm>  // std::max
 #include <bitset>
 #include <cmath>   // std::fabs
 #include <limits>  // std::numeric_limits
 
 #include "error_message_generation.h"
 #include "morphio/vector_types.h"
-#include "shared_utils.hpp"
+#include "point_utils.h"
 
 #include <ghc/filesystem.hpp>
 

src/shared_utils.hpp

@@ -6,8 +6,6 @@
 #include <morphio/errorMessages.h>
 #include <morphio/types.h>
 
-#include "point_utils.h"
-
 
 namespace morphio {
 

tests/test_1_swc.py

@@ -609,15 +609,32 @@ def test_multi_type_section():
     assert len(n.root_sections) == 1
     assert n.root_sections[0].id == 0
     assert_array_equal(n.root_sections[0].points,
-                       np.array([[0, 0, 2], ]))
+                       np.array([[0, 4, 0], [0, 0, 2]]))
     assert len(n.sections) == 4
-    assert_array_equal(n.section_offsets, [0, 1, 3, 5, 7])
+    assert_array_equal(n.section_offsets, [0, 2, 4, 6, 8])
     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_single_point_section(tmp_path):
+    '''SWC allows for single point sections, but MorphIO defines a section as having one segment, and a segment has two points'''
+    contents =('''\
+1 1 1 2 1 1 -1
+2 3 1 2 2 1 1       # <- single point section, bifurcates right away
+3 3 1 2 3 1 2
+4 3 1 2 4 1 2
+5 3 1 2 6 1 4
+6 3 1 2 7 1 4''')
+    m = Morphology(contents, "swc")
+    path = tmp_path / 'test_single_point_section.swc'
+    m.as_mutable().write(path)
+
+    with open(path) as fd:
+        # added MorphIO header, and column names
+        assert len(fd.readlines()) == len(contents.splitlines()) + 2
+
 def test_missing_parent():
     contents =('''
 1 1  0  0 0 10 -1

tests/test_utilities.cpp

@@ -9,6 +9,7 @@
 #include <filesystem>
 #include <fstream>
 
+#include "../src/point_utils.h"
 #include "../src/shared_utils.hpp"
 
 namespace fs = std::filesystem;