Repository: BlueBrain/MorphIO (language: python)

Issue:
Make three point soma warning more clear
This throws some warnings, without mention of three points soma, so it is un-understandable from a user not familiar with the internals of the code: https://github.com/BlueBrain/MorphIO/blob/54ffbd9fb01c16c9ccc8e57d4dc28f152a2962b4/src/shared_utils.cpp#L128 
I would be useful to mention three point soma in the warning, or something like that.


Discussion / hints:
Yup, I agree, thanks.

Fix patch:
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 <algorithm>  // std::max
 #include <bitset>
-#include <cstdint>
+#include <cmath>  // std::fabs
+#include <limits>  // std::numeric_limits
 
 #include "error_message_generation.h"
+#include "morphio/vector_types.h"
 #include "shared_utils.hpp"
 
 #include <ghc/filesystem.hpp>
@@ -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<Point, 3>& points, floatType radius) {
     //  NeuroMorpho is the main provider of morphologies, but they
@@ -95,14 +99,20 @@ ThreePointSomaStatus checkNeuroMorphoSoma(const std::array<Point, 3>& 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<floatType>::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<Point, 3>& 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<Point, 3>& 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;


Test patch:
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