From 68178a0b9d22c0facb4f00659f7bb8a9a79283a6 Mon Sep 17 00:00:00 2001 From: Mike Gevaert Date: Mon, 25 Nov 2024 14:08:24 +0100 Subject: [PATCH] make sure SWC reader follows the MorphIO invariant * when writing, this extra point is removed --- src/mut/writer_swc.cpp | 50 +++++++++++++++++++++++++++++------ src/readers/morphologySWC.cpp | 10 ++++++- src/shared_utils.cpp | 4 ++- src/shared_utils.hpp | 2 -- tests/test_1_swc.py | 21 +++++++++++++-- tests/test_utilities.cpp | 1 + 6 files changed, 74 insertions(+), 14 deletions(-) diff --git a/src/mut/writer_swc.cpp b/src/mut/writer_swc.cpp index 4373c5a1..830b3a4a 100644 --- a/src/mut/writer_swc.cpp +++ b/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& 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& section) { +bool _sameDiameter(const std::shared_ptr& 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 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]); diff --git a/src/readers/morphologySWC.cpp b/src/readers/morphologySWC.cpp index 0efc594a..a42c3805 100644 --- a/src/readers/morphologySWC.cpp +++ b/src/readers/morphologySWC.cpp @@ -478,7 +478,7 @@ class SWCBuilder std::make_unique(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) { diff --git a/src/shared_utils.cpp b/src/shared_utils.cpp index c0f6a693..130d8f47 100644 --- a/src/shared_utils.cpp +++ b/src/shared_utils.cpp @@ -2,6 +2,8 @@ * * SPDX-License-Identifier: Apache-2.0 */ +#include "shared_utils.hpp" + #include // std::max #include #include // std::fabs @@ -9,7 +11,7 @@ #include "error_message_generation.h" #include "morphio/vector_types.h" -#include "shared_utils.hpp" +#include "point_utils.h" #include diff --git a/src/shared_utils.hpp b/src/shared_utils.hpp index 5835d6c8..b700e1b5 100644 --- a/src/shared_utils.hpp +++ b/src/shared_utils.hpp @@ -6,8 +6,6 @@ #include #include -#include "point_utils.h" - namespace morphio { diff --git a/tests/test_1_swc.py b/tests/test_1_swc.py index 413a6e6d..a6d775bb 100644 --- a/tests/test_1_swc.py +++ b/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 diff --git a/tests/test_utilities.cpp b/tests/test_utilities.cpp index cca04f7e..f925e4d0 100644 --- a/tests/test_utilities.cpp +++ b/tests/test_utilities.cpp @@ -9,6 +9,7 @@ #include #include +#include "../src/point_utils.h" #include "../src/shared_utils.hpp" namespace fs = std::filesystem;