Example for PET analysis

requirements-dev.txt

@@ -5,3 +5,4 @@ flake8
 mypy
 meshio
 deepdiff
+pyamg

src/darsia/__init__.py

@@ -25,6 +25,8 @@
 from darsia.utils.linear_solvers.mg import *
 from darsia.utils.andersonacceleration import *
 from darsia.utils.dtype import *
+from darsia.utils.grid import *
+from darsia.utils.fv import *
 from darsia.corrections.basecorrection import *
 from darsia.corrections.shape.curvature import *
 from darsia.corrections.shape.affine import *
@@ -67,3 +69,5 @@
 from darsia.assistants.base_assistant import *
 from darsia.assistants.rotation_correction_assistant import *
 from darsia.assistants.subregion_assistant import *
+
+# from darsia.assistants.curvature_correction_assistant import *

src/darsia/image/arithmetics.py

@@ -32,9 +32,10 @@ def weight(img: darsia.Image, weight: Union[float, int, darsia.Image]) -> darsia
         weighted_img.img *= weight
 
     elif isinstance(weight, darsia.Image):
-        assert darsia.check_equal_coordinatesystems(
+        equal_coordinate_system, log = darsia.check_equal_coordinatesystems(
             img.coordinatesystem, weight.coordinatesystem, exclude_size=True
         )
+        assert equal_coordinate_system, f"{log}"
         space_dim = img.space_dim
         assert len(weight.img.shape) == space_dim
 

src/darsia/image/coordinatesystem.py

@@ -225,7 +225,7 @@ def check_equal_coordinatesystems(
     coordinatesystem1: CoordinateSystem,
     coordinatesystem2: CoordinateSystem,
     exclude_size: bool = False,
-) -> bool:
+) -> tuple[bool, dict]:
     """Check whether two coordinate systems are equivalent, i.e., they share basic
     attributes.
 
@@ -236,38 +236,64 @@ def check_equal_coordinatesystems(
 
     Returns:
         bool: True iff the two coordinate systems are equivalent.
+        dict: log of the failed checks.
 
     """
     success = True
-    success = success and coordinatesystem1.indexing == coordinatesystem2.indexing
-    success = success and coordinatesystem1.dim == coordinatesystem2.dim
+    failure_log = []
+
+    if not (coordinatesystem1.indexing == coordinatesystem2.indexing):
+        failure_log.append("indexing")
+
+    if not (coordinatesystem1.dim == coordinatesystem2.dim):
+        failure_log.append("dim")
+
     if not exclude_size:
-        success = success and np.allclose(
-            coordinatesystem1.shape, coordinatesystem2.shape
-        )
-    success = success and np.allclose(
-        coordinatesystem1.dimensions, coordinatesystem2.dimensions
-    )
-    success = success and coordinatesystem1.axes == coordinatesystem2.axes
+        if not (np.allclose(coordinatesystem1.shape, coordinatesystem2.shape)):
+            failure_log.append("shape")
+
+    if not (np.allclose(coordinatesystem1.dimensions, coordinatesystem2.dimensions)):
+        failure_log.append("dimensions")
+
+    if not (coordinatesystem1.axes == coordinatesystem2.axes):
+        failure_log.append("axes")
+
     if not exclude_size:
+        voxel_size_equal = True
         for axis in coordinatesystem1.axes:
-            success = (
-                success
+            voxel_size_equal = (
+                voxel_size_equal
                 and coordinatesystem1.voxel_size[axis]
                 == coordinatesystem2.voxel_size[axis]
             )
-    success = success and np.allclose(
-        coordinatesystem1._coordinate_of_origin_voxel,
-        coordinatesystem2._coordinate_of_origin_voxel,
-    )
-    success = success and np.allclose(
-        coordinatesystem1._coordinate_of_opposite_voxel,
-        coordinatesystem2._coordinate_of_opposite_voxel,
-    )
-    if not exclude_size:
-        success = success and np.allclose(
-            coordinatesystem1._voxel_of_origin_coordinate,
-            coordinatesystem2._voxel_of_origin_coordinate,
+        if not voxel_size_equal:
+            failure_log.append("voxel_size")
+
+    if not (
+        np.allclose(
+            coordinatesystem1._coordinate_of_origin_voxel,
+            coordinatesystem2._coordinate_of_origin_voxel,
+        )
+    ):
+        failure_log.append("coordinate_of_origin_voxel")
+
+    if not (
+        np.allclose(
+            coordinatesystem1._coordinate_of_opposite_voxel,
+            coordinatesystem2._coordinate_of_opposite_voxel,
         )
+    ):
+        failure_log.append("coordinate_of_opposite_voxel")
+
+    if not exclude_size:
+        if not (
+            np.allclose(
+                coordinatesystem1._voxel_of_origin_coordinate,
+                coordinatesystem2._voxel_of_origin_coordinate,
+            )
+        ):
+            failure_log.append("voxel_of_origin_coordinate")
+
+    success = len(failure_log) == 0
 
-    return success
+    return success, failure_log

src/darsia/measure/emd.py

@@ -46,6 +46,10 @@ def __call__(
             float or array: distance between img_1 and img_2.
 
         """
+        # Two-dimensional
+        if not (img_1.space_dim == 2 and img_2.space_dim == 2):
+            raise NotImplementedError("EMD only implemented for 2d.")
+
         # FIXME investigation required regarding resize preprocessing...
         # Preprocess images
         preprocessed_img_1 = self._preprocess(img_1)
@@ -116,13 +120,11 @@ def _compatibility_check(
         # Series
         assert img_1.time_num == img_2.time_num
 
-        # Two-dimensional
-        assert img_1.space_dim == 2 and img_2.space_dim == 2
-
         # Check whether the coordinate system is compatible
-        assert darsia.check_equal_coordinatesystems(
+        equal_coordinate_system, log = darsia.check_equal_coordinatesystems(
             img_1.coordinatesystem, img_2.coordinatesystem
         )
+        assert equal_coordinate_system, f"{log}"
         assert np.allclose(img_1.voxel_size, img_2.voxel_size)
 
         # Compatible distributions - comparing sums is sufficient since it is implicitly