From 7a119f2335ed0d8009637dc4c0f3f6a25d00983b Mon Sep 17 00:00:00 2001
From: Dimitri RODARIE <d.rodarie@gmail.com>
Date: Tue, 19 Mar 2024 11:09:25 +0100
Subject: [PATCH] Add check to the region filter in case the input dataset is
 all equal to 0 for the region of interest. (#74)

---
 atlas_densities/densities/utils.py |  9 ++++-----
 tests/densities/test_utils.py      | 14 ++++++++++++++
 2 files changed, 18 insertions(+), 5 deletions(-)

diff --git a/atlas_densities/densities/utils.py b/atlas_densities/densities/utils.py
index bc6430b..cd917f9 100644
--- a/atlas_densities/densities/utils.py
+++ b/atlas_densities/densities/utils.py
@@ -68,11 +68,11 @@ def normalize_intensity(
         3D array, the normalized marker intensity array.
     """
 
-    outside_mean = np.mean(
-        marker_intensity[np.logical_and((annotation == region_id), (marker_intensity > 0.0))]
-    )
     output_intensity = copy_array(marker_intensity, copy=copy)
-    output_intensity -= outside_mean * threshold_scale_factor
+    mask_low_filter = np.logical_and((annotation == region_id), (marker_intensity > 0.0))
+    if np.any(mask_low_filter):
+        outside_mean = np.mean(marker_intensity[mask_low_filter])
+        output_intensity -= outside_mean * threshold_scale_factor
     output_intensity[output_intensity < 0.0] = 0.0
     max_intensity = np.max(output_intensity)
     if max_intensity > 0:
@@ -405,7 +405,6 @@ def materialize(node):
 
 
 def get_group_ids(region_map: "RegionMap", config: dict, _skip_check=False) -> dict[str, set[int]]:
-
     """
     Get AIBS structure ids for several region groups of interest.
 
diff --git a/tests/densities/test_utils.py b/tests/densities/test_utils.py
index 51b9ba8..5680a7f 100644
--- a/tests/densities/test_utils.py
+++ b/tests/densities/test_utils.py
@@ -59,6 +59,20 @@ def test_normalize():
     # In place modification
     actual = tested.normalize_intensity(marker, annotation_raw, copy=False)
     npt.assert_array_almost_equal(marker, expected)
+    # Test outside with no intensity
+    marker = np.array(
+        [
+            [
+                [0.0, 0.0, 0.0, 0.0],
+                [0.0, 0.0, 1.0, 0.0],
+                [0.0, 3.0, 5.0, 0.0],
+                [0.0, 0.0, 0.0, 0.0],
+            ]
+        ],
+        dtype=float,
+    )
+    actual = tested.normalize_intensity(marker, annotation_raw, copy=True)
+    npt.assert_array_almost_equal(marker / np.max(marker), actual)
 
 
 def test_compensate_cell_overlap():