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

@@ -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.
 

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():