From ca63370300d02fe16246807ba1edc0e1561199ca Mon Sep 17 00:00:00 2001
From: Mike Gevaert <michael.gevaert@epfl.ch>
Date: Tue, 23 Apr 2024 11:30:30 +0200
Subject: [PATCH] add more tests to show pass/fail

---
 src/shared_utils.cpp | 24 +++++++++++++++---------
 tests/test_1_swc.py  | 29 +++++++++++++++++++++++++++++
 2 files changed, 44 insertions(+), 9 deletions(-)

diff --git a/src/shared_utils.cpp b/src/shared_utils.cpp
index b71e5bbd..c51418f7 100644
--- a/src/shared_utils.cpp
+++ b/src/shared_utils.cpp
@@ -8,6 +8,7 @@
 #include <limits>  // std::numeric_limits
 
 #include "error_message_generation.h"
+#include "morphio/vector_types.h"
 #include "shared_utils.hpp"
 
 #include <ghc/filesystem.hpp>
@@ -98,16 +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 withinRelativeEpsilon = [](floatType a, floatType b) {
+    auto withinTolerance = [](floatType a, floatType b) {
         floatType diff = std::fabs(a - b);
-        return diff <=
-               (std::max)(std::fabs(a), std::fabs(b)) * std::numeric_limits<floatType>::epsilon();
+        // 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 (size_t i = 0; i < 3; ++i) {
-        column_mask[i] = (withinRelativeEpsilon(points[0][i], points[1][i]) &&
-                          withinRelativeEpsilon(points[0][i], points[2][i]));
+        column_mask[i] = (withinTolerance(points[0][i], points[1][i]) &&
+                          withinTolerance(points[0][i], points[2][i]));
     }
 
     if (column_mask.none()) {
@@ -119,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 (!(withinRelativeEpsilon(points[0][col], points[1][col] - radius) &&
-          withinRelativeEpsilon(points[0][col], points[2][col] + radius)) &&
-        !(withinRelativeEpsilon(points[0][col], points[1][col] + radius) &&
-          withinRelativeEpsilon(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;
     }
 
diff --git a/tests/test_1_swc.py b/tests/test_1_swc.py
index bc21dd55..f2934970 100644
--- a/tests/test_1_swc.py
+++ b/tests/test_1_swc.py
@@ -341,6 +341,35 @@ def test_swc_threepoint_soma_tolerance():
     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'''