Skip to content

Commit

Permalink
Make ThreePointSoma use relative tolerance (#494)
Browse files Browse the repository at this point in the history 
* Make ThreePointSoma use relative tolerance
* add test
  • Loading branch information
@mgeplf
mgeplf authored Apr 23, 2024
1 parent cb725a1 commit 6f3d653
Show file tree
Hide file tree
Showing 2 changed files with 75 additions and 19 deletions.
41 changes: 26 additions & 15 deletions src/shared_utils.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down Expand Up @@ -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
Expand All @@ -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()) {
Expand All @@ -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;
}

Expand All @@ -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;
Expand Down
53 changes: 49 additions & 4 deletions tests/test_1_swc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand All @@ -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
Expand All @@ -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
Expand Down

0 comments on commit 6f3d653

Please sign in to comment.