diff --git a/src/shared_utils.cpp b/src/shared_utils.cpp index 399a7088..c51418f7 100644 --- a/src/shared_utils.cpp +++ b/src/shared_utils.cpp @@ -2,10 +2,13 @@ * * SPDX-License-Identifier: Apache-2.0 */ +#include // std::max #include -#include +#include // std::fabs +#include // std::numeric_limits #include "error_message_generation.h" +#include "morphio/vector_types.h" #include "shared_utils.hpp" #include @@ -70,6 +73,7 @@ bool is_regular_file(const std::string& path) { std::string join_path(const std::string& dirname, const std::string& filename) { return (ghc::filesystem::path(dirname) / filename).string(); } + namespace details { ThreePointSomaStatus checkNeuroMorphoSoma(const std::array& points, floatType radius) { // NeuroMorpho is the main provider of morphologies, but they @@ -95,14 +99,20 @@ ThreePointSomaStatus checkNeuroMorphoSoma(const std::array& points, fl // 2 1 x y (z + r) r 1 <- have `+` first // 3 1 x y (z - r) r 1 - auto withinEpsilon = [](floatType a, floatType b) { - return std::fabs(a - b) < morphio::epsilon; + auto withinTolerance = [](floatType a, floatType b) { + floatType diff = std::fabs(a - b); + // either we below the absolute epsilon... + return diff < morphio::epsilon || + // or within a reasonable tolerance... + // note: `(std::max)` is to work around `#define max` existing on some MSVC versions + (diff <= + (std::max)(std::fabs(a), std::fabs(b)) * std::numeric_limits::epsilon()); }; std::bitset<3> column_mask = {}; - for (uint8_t i = 0; i < 3; ++i) { - column_mask[i] = (withinEpsilon(points[0][i], points[1][i]) && - withinEpsilon(points[0][i], points[2][i])); + for (size_t i = 0; i < 3; ++i) { + column_mask[i] = (withinTolerance(points[0][i], points[1][i]) && + withinTolerance(points[0][i], points[2][i])); } if (column_mask.none()) { @@ -114,11 +124,12 @@ ThreePointSomaStatus checkNeuroMorphoSoma(const std::array& points, fl } const size_t col = !column_mask[0] ? 0 : !column_mask[1] ? 1 : 2; + std::cout << "asdf\n"; - if (!(withinEpsilon(points[0][col], points[1][col] - radius) && - withinEpsilon(points[0][col], points[2][col] + radius)) && - !(withinEpsilon(points[0][col], points[1][col] + radius) && - withinEpsilon(points[0][col], points[2][col] - radius))) { + if (!(withinTolerance(points[0][col], points[1][col] - radius) && + withinTolerance(points[0][col], points[2][col] + radius)) && + !(withinTolerance(points[0][col], points[1][col] + radius) && + withinTolerance(points[0][col], points[2][col] - radius))) { return NotRadiusOffset; } @@ -128,19 +139,19 @@ ThreePointSomaStatus checkNeuroMorphoSoma(const std::array& points, fl std::ostream& operator<<(std::ostream& os, ThreePointSomaStatus s) { switch (s) { case ZeroColumnsAreTheSame: - os << "None of the columns (ie: all the X, Y or Z values) are the same."; + os << "Three Point Soma: None of the columns (ie: all the X, Y or Z values) are the same."; break; case OneColumnIsTheSame: - os << "Only one column has the same coordinates."; + os << "Three Point Soma: Only one column has the same coordinates."; break; case ThreeColumnsAreTheSame: - os << "All three columns have the same coordinates."; + os << "Three Point Soma: All three columns have the same coordinates."; break; case NotRadiusOffset: - os << "The non-constant columns is not offset by +/- the radius from the initial sample."; + os << "Three Point Soma: The non-constant columns is not offset by +/- the radius from the initial sample."; break; case Conforms: - os << "Three point soma conforms"; + os << "Three Point Soma: conforms to specification"; break; } return os; diff --git a/tests/test_1_swc.py b/tests/test_1_swc.py index 77b339d1..f2934970 100644 --- a/tests/test_1_swc.py +++ b/tests/test_1_swc.py @@ -305,8 +305,10 @@ def test_soma_type_3_point(tmp_path): 3 1 0 0 0 3.0 1 # PID is 1''') assert (Morphology(content, extension='swc').soma_type == SomaType.SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS) - assert strip_color_codes(err.getvalue()).strip() == \ - '$STRING$:0:warning\nOnly one column has the same coordinates.' + assert strip_color_codes(err.getvalue()).strip() == ( + "$STRING$:0:warning\n" + "Three Point Soma: Only one column has the same coordinates." + ) with captured_output() as (_, err): with ostream_redirect(stdout=True, stderr=True): @@ -315,8 +317,10 @@ def test_soma_type_3_point(tmp_path): 3 1 0 0 0 3.0 1 # PID is 1''') assert (Morphology(content, extension='swc').soma_type == SomaType.SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS) - assert strip_color_codes(err.getvalue()).strip() == \ - '$STRING$:0:warning\nThe non-constant columns is not offset by +/- the radius from the initial sample.' + assert strip_color_codes(err.getvalue()).strip() == ( + "$STRING$:0:warning\n" + "Three Point Soma: The non-constant columns is not offset by +/- the radius from the initial sample." + ) # If this configuration is not respected -> SOMA_CYLINDERS content = ( '''1 1 0 0 0 3.0 -1 @@ -326,6 +330,47 @@ def test_soma_type_3_point(tmp_path): SomaType.SOMA_CYLINDERS) +def test_swc_threepoint_soma_tolerance(): + # from https://github.com/BlueBrain/MorphIO/issues/492 + contents = """ + 1 1 5789.674999999998 1322.85 2846.65 5.165118027 -1 + 2 1 5789.674999999998 1317.6848819729996 2846.65 5.165118027 1 + 3 1 5789.674999999998 1328.015118027 2846.65 5.165118027 1 + """ + warnings = morphio.WarningHandlerCollector() + Morphology(contents, extension="swc", warning_handler=warnings) + assert len(warnings.get_all()) == 0 + + contents = """ + 1 1 10000 10000 10000 100.999999999999999999 -1 + 2 1 10000 9899.000000000001 10000 100.999999999999999999 1 + 3 1 10000 10100.99999999999 10000 100.999999999999999999 1 + """ + warnings = morphio.WarningHandlerCollector() + Morphology(contents, extension="swc", warning_handler=warnings) + assert len(warnings.get_all()) == 0 + + # This still passes, even with the third point being off by 0.0000008 which seems reasonable + contents = """ + 1 1 10000 10000 10000 100.999999999999999999 -1 + 2 1 10000 9899.000000000001 10000 100.999999999999999999 1 + 3 1 10000 10100.99999911111 10000 100.999999999999999999 1 + """ + warnings = morphio.WarningHandlerCollector() + Morphology(contents, extension="swc", warning_handler=warnings) + assert len(warnings.get_all()) == 0 + + # The third point being off by 0.0088888 creates warning + contents = """ + 1 1 10000 10000 10000 100.999999999999999999 -1 + 2 1 10000 9899.000000000001 10000 100.999999999999999999 1 + 3 1 10000 10100.99111111111 10000 100.999999999999999999 1 + """ + warnings = morphio.WarningHandlerCollector() + Morphology(contents, extension="swc", warning_handler=warnings) + assert len(warnings.get_all()) == 1 + + def test_read_weird_ids(): '''The ordering of IDs is not required''' content = ('''10000 3 0 0 5 0.5 100 # neurite 4th point