fest: Add randomization option for correlation calculation (#421)

Co-authored-by: Sourcery AI <>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: sourcery-ai[bot] <58596630+sourcery-ai[bot]@users.noreply.github.com>

package/PartSeg/common_gui/algorithms_description.py

@@ -170,11 +170,11 @@ def _get_field_from_value_type(cls, ap: AlgorithmProperty):
             res = QCheckBox()
         elif issubclass(ap.value_type, (float, int)):
             res = cls._get_numeric_field(ap)
-        elif issubclass(ap.value_type, str):
-            res = QLineEdit()
         elif issubclass(ap.value_type, Enum):
+            # noinspection PyTypeChecker
             res = QEnumComboBox(enum_class=ap.value_type)
-            # noinspection PyUnresolvedReferences
+        elif issubclass(ap.value_type, str):
+            res = QLineEdit()
         elif issubclass(ap.value_type, ROIExtractionProfile):
             res = ProfileSelect()
         elif issubclass(ap.value_type, list):

package/PartSegCore/analysis/measurement_calculation.py

@@ -42,10 +42,15 @@
 
 NO_COMPONENT = -1
 
-PEARSON_CORRELATION = "Pearson correlation coefficient"
-MANDERS_COEFIICIENT = "Mander's overlap coefficient"
-INTENSITY_CORRELATION = "Intensity correlation quotient"
-SPEARMAN_CORRELATION = "Spearman rank correlation"
+
+class CorrelationEnum(str, Enum):
+    pearson = "Pearson correlation coefficient"
+    manders = "Mander's overlap coefficient"
+    intensity = "Intensity correlation quotient"
+    spearman = "Spearman rank correlation"
+
+    def __str__(self):
+        return self.value
 
 
 class ProhibitedDivision(Exception):
@@ -1614,43 +1619,57 @@ class ColocalizationMeasurement(MeasurementMethodBase):
     def get_fields(cls):
         return [
             AlgorithmProperty("channel_fst", "Channel 1", 0, value_type=Channel),
-            AlgorithmProperty("channel_scd", "Channel 2", 0, value_type=Channel),
+            AlgorithmProperty("channel_scd", "Channel 2", 1, value_type=Channel),
             AlgorithmProperty(
                 "colocalization",
                 "Colocalization",
-                PEARSON_CORRELATION,
-                possible_values=[PEARSON_CORRELATION, MANDERS_COEFIICIENT, INTENSITY_CORRELATION, SPEARMAN_CORRELATION],
+                CorrelationEnum.pearson,
+            ),
+            AlgorithmProperty(
+                "randomize", "Randomize channel", False, help_text="If randomize orders of pixels in one channel"
+            ),
+            AlgorithmProperty(
+                "randomize_repeat", "Randomize num", 10, help_text="Number of repetitions for mean_calculate"
             ),
         ]
 
     @staticmethod
-    def calculate_property(area_array, colocalization, channel_fst=0, channel_scd=1, **kwargs):  # pylint: disable=W0221
-        mask_binary = area_array > 0
-        data_1 = kwargs[f"channel_{channel_fst}"][mask_binary].astype(float)
-        data_2 = kwargs[f"channel_{channel_scd}"][mask_binary].astype(float)
-        # data_max = max(data_1.max(), data_2.max())
-        # data_1 = data_1 / data_max
-        # data_2 = data_2 / data_max
-        if colocalization == SPEARMAN_CORRELATION:
+    def _calculate_masked(data_1, data_2, colocalization):
+        if colocalization == CorrelationEnum.spearman:
             data_1 = data_1.argsort().argsort().astype(float)
             data_2 = data_2.argsort().argsort().astype(float)
-            colocalization = PEARSON_CORRELATION
-        if colocalization == PEARSON_CORRELATION:
+            colocalization = CorrelationEnum.pearson
+        if colocalization == CorrelationEnum.pearson:
             data_1_mean = np.mean(data_1)
             data_2_mean = np.mean(data_2)
             nominator = np.sum((data_1 - data_1_mean) * (data_2 - data_2_mean))
             numerator = np.sqrt(np.sum((data_1 - data_1_mean) ** 2) * np.sum((data_2 - data_2_mean) ** 2))
             return nominator / numerator
-        if colocalization == MANDERS_COEFIICIENT:
+        if colocalization == CorrelationEnum.manders:
             nominator = np.sum(data_1 * data_2)
             numerator = np.sqrt(np.sum(data_1 ** 2) * np.sum(data_2 ** 2))
             return nominator / numerator
-        if colocalization == INTENSITY_CORRELATION:
+        if colocalization == CorrelationEnum.intensity:
             data_1_mean = np.mean(data_1)
             data_2_mean = np.mean(data_2)
             return np.sum((data_1 > data_1_mean) == (data_2 > data_2_mean)) / data_1.size - 0.5
 
-        raise RuntimeError("Not supported colocalization method")  # pragma: no cover
+        raise RuntimeError(f"Not supported colocalization method {colocalization}")  # pragma: no cover
+
+    @classmethod
+    def calculate_property(
+        cls, area_array, colocalization, randomize=False, randomize_repeat=10, channel_fst=0, channel_scd=1, **kwargs
+    ):  # pylint: disable=W0221
+        mask_binary = area_array > 0
+        data_1 = kwargs[f"channel_{channel_fst}"][mask_binary].astype(float)
+        data_2 = kwargs[f"channel_{channel_scd}"][mask_binary].astype(float)
+        if not randomize:
+            return cls._calculate_masked(data_1, data_2, colocalization)
+        res_list = []
+        for _ in range(randomize_repeat):
+            rand_data2 = np.random.permutation(data_2)
+            res_list.append(cls._calculate_masked(data_1, rand_data2, colocalization))
+        return np.mean(res_list)
 
     @classmethod
     def get_units(cls, ndim) -> symbols:

package/tests/test_PartSegCore/test_measurements.py

@@ -15,12 +15,11 @@
 from PartSegCore.analysis.measurement_base import AreaType, Leaf, MeasurementEntry, Node, PerComponent
 from PartSegCore.analysis.measurement_calculation import (
     HARALIC_FEATURES,
-    INTENSITY_CORRELATION,
-    MANDERS_COEFIICIENT,
     MEASUREMENT_DICT,
     ColocalizationMeasurement,
     ComponentsInfo,
     ComponentsNumber,
+    CorrelationEnum,
     Diameter,
     DistanceMaskROI,
     DistancePoint,
@@ -183,9 +182,11 @@ def test_square(self):
         mask1 = image.get_channel(0)[0] > 40
         mask2 = image.get_channel(0)[0] > 60
         mask3 = mask1 * ~mask2
-        assert PixelBrightnessSum.calculate_property(mask1, image.get_channel(0)) == 60 * 60 * 50 + 40 * 40 * 20
-        assert PixelBrightnessSum.calculate_property(mask2, image.get_channel(0)) == 40 * 40 * 70
-        assert PixelBrightnessSum.calculate_property(mask3, image.get_channel(0)) == (60 * 60 - 40 * 40) * 50
+        assert PixelBrightnessSum.calculate_property(mask1, image.get_channel(0)) == 60 ** 2 * 50 + 40 * 40 * 20
+
+        assert PixelBrightnessSum.calculate_property(mask2, image.get_channel(0)) == 40 ** 2 * 70
+
+        assert PixelBrightnessSum.calculate_property(mask3, image.get_channel(0)) == (60 ** 2 - 40 ** 2) * 50
 
     def test_empty(self):
         image = get_cube_image()
@@ -266,7 +267,7 @@ def test_square(self):
         mask3 = mask1 * ~mask2
         assert Voxels.calculate_property(mask1) == 60 * 60
         assert Voxels.calculate_property(mask2) == 40 * 40
-        assert Voxels.calculate_property(mask3) == 60 * 60 - 40 * 40
+        assert Voxels.calculate_property(mask3) == 60 ** 2 - 40 ** 2
 
     def test_empty(self):
         image = get_cube_image()
@@ -1236,7 +1237,7 @@ def test_square_equal_radius(self):
                 voxel_size=image.voxel_size,
                 result_scalar=1,
             )
-            == (60 * 60 - 40 * 40) * result_scale
+            == (60 ** 2 - 40 ** 2) * result_scale
         )
 
         assert (
@@ -1249,7 +1250,7 @@ def test_square_equal_radius(self):
                 voxel_size=image.voxel_size,
                 result_scalar=1,
             )
-            == (60 * 60 - 30 * 30) * result_scale
+            == (60 * 60 - 30 ** 2) * result_scale
         )
 
         assert (
@@ -1309,7 +1310,7 @@ def test_square_equal_volume(self):
                 voxel_size=image.voxel_size,
                 result_scalar=1,
             )
-            == (60 * 60 - 44 * 44) * result_scale
+            == (60 ** 2 - 44 * 44) * result_scale
         )
 
         assert (
@@ -1434,10 +1435,10 @@ def test_cube_equal_volume(self, nr, sum_val, diff_array):
         "nr, sum_val, diff_array, equal_volume",
         [
             (3, (60 * 60 - 40 * 40) * 50, False, False),
-            (2, (60 * 60 - 40 * 40) * 50 + (40 * 40 - 30 * 30) * 70, False, False),
+            (2, (60 ** 2 - 40 * 40) * 50 + (40 * 40 - 30 * 30) * 70, False, False),
             (3, 0, True, False),
-            (2, (40 * 40 - 30 * 30) * 70, True, False),
-            (3, (60 * 60 - 50 * 50) * 50, False, True),
+            (2, (40 ** 2 - 30 ** 2) * 70, True, False),
+            (3, (60 * 60 - 50 ** 2) * 50, False, True),
             (2, (60 * 60 - 44 * 44) * 50, False, True),
             (3, 0, True, True),
             (2, 0, True, True),
@@ -2155,6 +2156,7 @@ def test_all_methods(method, dtype):
         channel=data,
         channel_num=0,
         channel_0=data,
+        channel_1=data,
         voxel_size=(1, 1, 1),
         result_scalar=1,
         roi_alternative={},
@@ -2172,13 +2174,13 @@ def test_all_methods(method, dtype):
 )
 @pytest.mark.parametrize("area", [AreaType.ROI, AreaType.Mask])
 def test_per_component(method, area):
-    data = np.zeros((10, 20, 20), dtype=np.uint8)
+    data = np.zeros((10, 20, 20, 2), dtype=np.uint8)
     data[1:-1, 3:-3, 3:-3] = 2
     data[1:-1, 4:-4, 4:-4] = 3
     data[1:-1, 6, 6] = 5
-    roi = (data > 2).astype(np.uint8)
-    mask = (data > 0).astype(np.uint8)
-    image = Image(data, image_spacing=(10 ** -8,) * 3, axes_order="ZYX")
+    roi = (data[..., 0] > 2).astype(np.uint8)
+    mask = (data[..., 0] > 0).astype(np.uint8)
+    image = Image(data, image_spacing=(10 ** -8,) * 3, axes_order="ZYXC")
     image.set_mask(mask, axes="ZYX")
 
     statistics = [
@@ -2207,39 +2209,44 @@ def test_per_component(method, area):
     assert result["Measurement per component"][0][0] == result["Measurement"][0]
 
 
-@pytest.mark.parametrize("method", ColocalizationMeasurement.get_fields()[-1].possible_values)
-def test_colocalization(method):
+@pytest.mark.parametrize("method", CorrelationEnum.__members__.values())
+@pytest.mark.parametrize("randomize", [True, False])
+def test_colocalization(method, randomize):
     area_array = np.ones((10, 10))
     data = np.random.rand(10, 10)
-    factor = 0.5 if method == INTENSITY_CORRELATION else 1
+    factor = 0.5 if method == CorrelationEnum.intensity else 1
     value = ColocalizationMeasurement.calculate_property(
         area_array=area_array,
         channel_0=data,
         channel_1=data,
         colocalization=method,
+        randomize=randomize,
     )
-    assert value == factor
+    assert value == factor or randomize
     value = ColocalizationMeasurement.calculate_property(
         area_array=area_array,
         channel_0=data,
         channel_1=data * 100,
         colocalization=method,
+        randomize=randomize,
     )
-    assert isclose(value, factor)
+    assert isclose(value, factor) or randomize
 
     value = ColocalizationMeasurement.calculate_property(
         area_array=area_array,
         channel_0=data,
         channel_1=data + 100,
         colocalization=method,
+        randomize=randomize,
     )
 
-    assert isclose(value, factor) or (method == MANDERS_COEFIICIENT and value < 1)
+    assert isclose(value, factor) or (method == CorrelationEnum.manders and value < 1) or randomize
 
     value = ColocalizationMeasurement.calculate_property(
         area_array=area_array,
         channel_0=data,
         channel_1=-data,
         colocalization=method,
+        randomize=randomize,
     )
-    assert value == -factor
+    assert value == -factor or randomize