Merge branch 'orthogonal_interpolation' of https://github.com/clEsper…

…anto/pyclesperanto_prototype into orthogonal_interpolation
clEsperanto · Dec 2, 2022 · 8959d82 · 8959d82
2 parents 0493060 + 654b296
commit 8959d82
Show file tree

Hide file tree

Showing 8 changed files with 875 additions and 74 deletions.
diff --git a/demo/transforms/deskew_artifact_interpolation-0.19.4.ipynb b/demo/transforms/deskew_artifact_interpolation-0.19.4.ipynb
diff --git a/demo/transforms/deskew_artifact_interpolation.ipynb b/demo/transforms/deskew_artifact_interpolation.ipynb
diff --git a/pyclesperanto_prototype/_tier8/_affine_transform_deskew_3d.py b/pyclesperanto_prototype/_tier8/_affine_transform_deskew_3d.py
@@ -22,8 +22,7 @@ def affine_transform_deskew_3d(source: Image, destination: Image = None,
                                voxel_size_x: float = 0.1449922,
                                voxel_size_y: float = 0.1449922,
                                voxel_size_z: float = 0.3,
-                               deskew_direction: DeskewDirection = DeskewDirection.Y,
-                               auto_size: bool = False) -> Image:
+                               deskew_direction: DeskewDirection = DeskewDirection.Y) -> Image:
     """
     Applies an affine transform to deskew an image. 
     Uses orthogonal interpolation (Sapoznik et al. (2020)  https://doi.org/10.7554/eLife.57681)
@@ -34,7 +33,7 @@ def affine_transform_deskew_3d(source: Image, destination: Image = None,
         image to be transformed
     destination : Image, optional
         image where the transformed image should be written to
-    transform : 4x4 numpy array or AffineTransform3D object or skimage.transform.AffineTransform object or str, optional
+    transform : AffineTransform3D object, optional
         transform matrix or object or string describing the transformation
     deskewing_angle_in_degrees: float, optional
         Oblique plane or deskewing acquisition angle
@@ -65,38 +64,15 @@ def affine_transform_deskew_3d(source: Image, destination: Image = None,
     from .._tier0 import execute
     from .._tier0 import create
 
-    assert len(
-        source.shape) == 3, f"Image needs to be 3D, got shape of {len(source.shape)}"
+    assert len(source.shape) == 3, f"Image needs to be 3D, got shape of {len(source.shape)}"
 
     # handle output creation
-    if auto_size and isinstance(transform, AffineTransform3D):
-        new_size, transform, _ = _determine_translation_and_bounding_box(
-            source, transform)
+    new_size, transform, _ = _determine_translation_and_bounding_box(source, transform)
     if destination is None:
-        if auto_size and isinstance(transform, AffineTransform3D):
-            # This modifies the given transform
-            destination = create(new_size)
-        else:
-            destination = create_like(source)
-
-    if isinstance(transform, str):
-        transform = AffineTransform3D(transform, source)
-
+        destination = create(new_size)
+
     # we invert the transform because we go from the target image to the source image to read pixels
-    if isinstance(transform, AffineTransform3D):
-        transform_matrix = np.asarray(transform.copy().inverse())
-    elif isinstance(transform, AffineTransform):
-        # Question: Don't we have to invert this one as well? haesleinhuepf
-        matrix = np.asarray(transform.params)
-        matrix = np.asarray([
-            [matrix[0, 0], matrix[0, 1], 0, matrix[0, 2]],
-            [matrix[1, 0], matrix[1, 1], 0, matrix[1, 2]],
-            [0, 0, 1, 0],
-            [matrix[2, 0], matrix[2, 1], 0, matrix[2, 2]]
-        ])
-        transform_matrix = np.linalg.inv(matrix)
-    else:
-        transform_matrix = np.linalg.inv(transform)
+    transform_matrix = np.asarray(transform.copy().inverse())
 
     # precalculate these functions that are dependent on deskewing angle
     tantheta = np.float32(np.tan(deskewing_angle_in_degrees * np.pi/180))
@@ -120,9 +96,6 @@ def affine_transform_deskew_3d(source: Image, destination: Image = None,
         "input": source,
         "output": destination,
         "mat": gpu_transform_matrix,
-        "deskewed_Nx": destination.shape[2],
-        "deskewed_Ny": destination.shape[1],
-        "deskewed_Nz": destination.shape[0],
         "pixel_step": pixel_step,
         "tantheta": tantheta,
         "costheta": costheta,

diff --git a/pyclesperanto_prototype/_tier8/_deskew_x.py b/pyclesperanto_prototype/_tier8/_deskew_x.py
@@ -12,11 +12,10 @@ def deskew_x(input_image: Image,
              voxel_size_y: float = 1,
              voxel_size_z: float = 1,
              scale_factor: float = 1,
-             linear_interpolation: bool = None
+             linear_interpolation: bool = True
              ) -> Image:
     """
     Deskew an image stack as acquired with oblique plane light-sheet microscopy, with skew in the X direction.
-    Uses orthogonal interpolation (Sapoznik et al. (2020)  https://doi.org/10.7554/eLife.57681)
 
     Parameters
     ----------
@@ -35,16 +34,19 @@ def deskew_x(input_image: Image,
         default: 1
         If the resulting image becomes too huge, it is possible to reduce output image size by this factor.
         The isotropic voxel size of the output image will then be voxel_size_x / scaling_factor.
-    linear_interpolation: obsolete
+    linear_interpolation: bool, optional
+        If True (default), uses orthogonal interpolation (Sapoznik et al. (2020)  https://doi.org/10.7554/eLife.57681)
+        If False, nearest-neighbor interpolation wille be applied.
 
     Returns
     -------
     output_image
     """
-    if linear_interpolation is not None:
-        warnings.warn("In function deskew_ the linear_interpolation parameters is obsolete. Orthogonal interpolation will always be used. Please remove it from your code.", DeprecationWarning)
+    if not linear_interpolation:
+        warnings.warn("linear_interpolation = False is deprecated due to deskewing artifacts. The linear_interpolation parameter will be removed in a future version.")
 
     from ._AffineTransform3D import AffineTransform3D
+    from ._affine_transform import affine_transform
     from ._affine_transform_deskew_3d import affine_transform_deskew_3d, DeskewDirection
 
     # define affine transformation
@@ -53,11 +55,14 @@ def deskew_x(input_image: Image,
                         voxel_size_z=voxel_size_z, scale_factor=scale_factor)
 
     # apply transform
-    return affine_transform_deskew_3d(source=input_image, destination=output_image,
-                                      transform=transform,
-                                      deskewing_angle_in_degrees=angle_in_degrees,
-                                      voxel_size_x=voxel_size_x,
-                                      voxel_size_y=voxel_size_y,
-                                      voxel_size_z=voxel_size_z,
-                                      deskew_direction=DeskewDirection.X,
-                                      auto_size=True)
+    if linear_interpolation:
+        return affine_transform_deskew_3d(source=input_image, destination=output_image,
+                                          transform=transform,
+                                          deskewing_angle_in_degrees=angle_in_degrees,
+                                          voxel_size_x=voxel_size_x,
+                                          voxel_size_y=voxel_size_y,
+                                          voxel_size_z=voxel_size_z,
+                                          deskew_direction=DeskewDirection.X)
+    else:
+        return affine_transform(input_image, output_image, transform=transform, auto_size=True,
+                                linear_interpolation=False)
diff --git a/pyclesperanto_prototype/_tier8/_deskew_y.py b/pyclesperanto_prototype/_tier8/_deskew_y.py
@@ -12,11 +12,10 @@ def deskew_y(input_image: Image,
              voxel_size_y: float = 1,
              voxel_size_z: float = 1,
              scale_factor: float = 1,
-             linear_interpolation: bool = None
+             linear_interpolation: bool = True
              ) -> Image:
     """
     Deskew an image stack as acquired with oblique plane light-sheet microscopy with skew in the Y direction.
-    Uses orthogonal interpolation (Sapoznik et al. (2020)  https://doi.org/10.7554/eLife.57681)
 
     Parameters
     ----------
@@ -35,16 +34,19 @@ def deskew_y(input_image: Image,
         default: 1
         If the resulting image becomes too huge, it is possible to reduce output image size by this factor.
         The isotropic voxel size of the output image will then be voxel_size_y / scaling_factor.
-    linear_interpolation: obsolete
+    linear_interpolation: bool, optional
+        If True (default), uses orthogonal interpolation (Sapoznik et al. (2020)  https://doi.org/10.7554/eLife.57681)
+        If False, nearest-neighbor interpolation wille be applied.
 
     Returns
     -------
     output_image
     """
-    if linear_interpolation is not None:
-        warnings.warn("In function deskew_ the linear_interpolation parameters is obsolete. Orthogonal interpolation will always be used. Please remove it from your code.", DeprecationWarning)
+    if not linear_interpolation:
+        warnings.warn("linear_interpolation = False is deprecated due to deskewing artifacts. The linear_interpolation parameter will be removed in a future version.")
 
     from ._AffineTransform3D import AffineTransform3D
+    from ._affine_transform import affine_transform
     from ._affine_transform_deskew_3d import affine_transform_deskew_3d, DeskewDirection
 
     # define affine transformation matrix
@@ -53,11 +55,14 @@ def deskew_y(input_image: Image,
                         voxel_size_z=voxel_size_z, scale_factor=scale_factor)
 
     # apply transform using special affine transform method for deskewing in the Y direction
-    return affine_transform_deskew_3d(source=input_image, destination=output_image,
-                                      transform=transform,
-                                      deskewing_angle_in_degrees=angle_in_degrees,
-                                      voxel_size_x=voxel_size_x,
-                                      voxel_size_y=voxel_size_y,
-                                      voxel_size_z=voxel_size_z,
-                                      deskew_direction=DeskewDirection.Y,
-                                      auto_size=True)
+    if linear_interpolation:
+        return affine_transform_deskew_3d(source=input_image, destination=output_image,
+                                          transform=transform,
+                                          deskewing_angle_in_degrees=angle_in_degrees,
+                                          voxel_size_x=voxel_size_x,
+                                          voxel_size_y=voxel_size_y,
+                                          voxel_size_z=voxel_size_z,
+                                          deskew_direction=DeskewDirection.Y)
+    else:
+        return affine_transform(input_image, output_image, transform=transform, auto_size=True,
+                                linear_interpolation=False)
diff --git a/pyclesperanto_prototype/_tier8/affine_transform_deskew_x_3d_x.cl b/pyclesperanto_prototype/_tier8/affine_transform_deskew_x_3d_x.cl
@@ -40,9 +40,7 @@
 
 __kernel void
 affine_transform_deskew_x_3d(IMAGE_input_TYPE input, IMAGE_output_TYPE output,
-                             IMAGE_mat_TYPE mat,
-                             const int deskewed_Nx, const int deskewed_Ny,
-                             const int deskewed_Nz, float pixel_step,
+                             IMAGE_mat_TYPE mat, float pixel_step,
                              float tantheta, float costheta, float sintheta) {
 
   const sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | SAMPLER_ADDRESS;
@@ -66,8 +64,8 @@ affine_transform_deskew_x_3d(IMAGE_input_TYPE input, IMAGE_output_TYPE output,
 
   // ensure within bounds of final image/deskewed image
 
-  if (x >= 0 && y >= 0 && z >= 0 && x < deskewed_Nx && y < deskewed_Ny &&
-      z < deskewed_Nz) {
+  if (x >= 0 && y >= 0 && z >= 0 && x < GET_IMAGE_WIDTH(output) && y < GET_IMAGE_HEIGHT(output) &&
+      z < GET_IMAGE_DEPTH(output)) {
 
     float virtual_plane = (x - z / tantheta);
     // get plane before
@@ -127,7 +125,7 @@ affine_transform_deskew_x_3d(IMAGE_input_TYPE input, IMAGE_output_TYPE output,
     }
   }
 
-  int4 pos = (int4){x, y, z, 0};
+  int4 pos = (int4){x, y, (GET_IMAGE_DEPTH(output) - 1 - z), 0};
 
   WRITE_output_IMAGE(output, pos, CONVERT_output_PIXEL_TYPE(pix));
 }
diff --git a/pyclesperanto_prototype/_tier8/affine_transform_deskew_y_3d_x.cl b/pyclesperanto_prototype/_tier8/affine_transform_deskew_y_3d_x.cl
@@ -41,9 +41,7 @@
 
 __kernel void
 affine_transform_deskew_y_3d(IMAGE_input_TYPE input, IMAGE_output_TYPE output,
-                             IMAGE_mat_TYPE mat,
-                             const int deskewed_Nx, const int deskewed_Ny,
-                             const int deskewed_Nz, float pixel_step,
+                             IMAGE_mat_TYPE mat, float pixel_step,
                              float tantheta, float costheta, float sintheta) {
 
   const sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | SAMPLER_ADDRESS;
@@ -72,8 +70,8 @@ affine_transform_deskew_y_3d(IMAGE_input_TYPE input, IMAGE_output_TYPE output,
 
   // ensure within bounds of final image/deskewed image
 
-  if (x >= 0 && y >= 0 && z >= 0 && x < deskewed_Nx && y < deskewed_Ny &&
-      z < deskewed_Nz) {
+  if (x >= 0 && y >= 0 && z >= 0 && x < GET_IMAGE_WIDTH(output) && y < GET_IMAGE_HEIGHT(output) &&
+      z < GET_IMAGE_DEPTH(output)) {
     //printf("%d\n,%d\n,%d\n", z,y,x);
     float virtual_plane = (y - z / tantheta);
     // get plane before
@@ -133,8 +131,7 @@ affine_transform_deskew_y_3d(IMAGE_input_TYPE input, IMAGE_output_TYPE output,
     }
   }
 
-  // if rotate coverslip, apply flipping on Z axis
-  int4 pos = (int4){x, y, z, 0};
+  int4 pos = (int4){x, y, (GET_IMAGE_DEPTH(output) - 1 - z), 0};
 
   WRITE_output_IMAGE(output, pos, CONVERT_output_PIXEL_TYPE(pix));
 }
diff --git a/tests/test_deskew.py b/tests/test_deskew.py
@@ -8,7 +8,7 @@
 
 reference = np.zeros((5, 19, 10),dtype=np.float32)
 #will this change with device used?
-reference[1, 2, 1] = 0.16987294
+reference[3, 2, 1] = 0.16987294
 
 
 def test_deskew_y():