Merge pull request #224 from clEsperanto/orthogonal_interpolation

Interpolation for deskewing

pyclesperanto_prototype/__init__.py

@@ -10,5 +10,5 @@
 from ._tier10 import *
 from ._tier11 import *
 
-__version__ = "0.19.4"
+__version__ = "0.20.0"
 __common_alias__ = "cle"

pyclesperanto_prototype/_tier8/__init__.py

@@ -7,5 +7,3 @@
 from ._rotate import rotate
 from ._scale import scale
 from ._translate import translate
-
-

pyclesperanto_prototype/_tier8/_affine_transform.py

@@ -8,8 +8,9 @@
 from skimage.transform import AffineTransform
 import numpy as np
 
+
 @plugin_function(output_creator=create_none)
-def affine_transform(source : Image, destination : Image = None, transform : Union[np.ndarray, AffineTransform3D, AffineTransform] = None, linear_interpolation : bool = False, auto_size:bool = False) -> Image:
+def affine_transform(source: Image, destination: Image = None, transform: Union[np.ndarray, AffineTransform3D, AffineTransform] = None, linear_interpolation: bool = False, auto_size: bool = False) -> Image:
     """
     Applies an affine transform to an image.
 
@@ -48,7 +49,8 @@ def affine_transform(source : Image, destination : Image = None, transform : Uni
 
     # 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
@@ -80,10 +82,10 @@ def affine_transform(source : Image, destination : Image = None, transform : Uni
         # 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]],
+            [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]]
+            [matrix[2, 0], matrix[2, 1], 0, matrix[2, 2]]
         ])
         transform_matrix = np.linalg.inv(matrix)
     else:
@@ -147,15 +149,18 @@ def _determine_translation_and_bounding_box(source: Image, affine_transformation
     else:
         nx, ny, nz = source.shape
 
-    original_bounding_box = [list(x) + [1] for x in product((0, nz), (0, ny), (0, nx))]
+    original_bounding_box = [list(x) + [1]
+                             for x in product((0, nz), (0, ny), (0, nx))]
     # transform the corners using the given affine transform
-    transformed_bounding_box = np.asarray(list(map(lambda x: affine_transformation._matrix @ x, original_bounding_box)))
+    transformed_bounding_box = np.asarray(
+        list(map(lambda x: affine_transformation._matrix @ x, original_bounding_box)))
 
     # the min and max coordinates tell us from where to where the image ranges (bounding box)
     min_coordinate = transformed_bounding_box.min(axis=0)
     max_coordinate = transformed_bounding_box.max(axis=0)
     # determine the size of the transformed bounding box
-    new_shape = np.around((max_coordinate - min_coordinate)[0:3]).astype(int).tolist()[::-1]
+    new_shape = np.around((max_coordinate - min_coordinate)
+                          [0:3]).astype(int).tolist()[::-1]
 
     # we make a copy to not modify the original transform
     new_affine_transform = AffineTransform3D()

pyclesperanto_prototype/_tier8/_affine_transform_deskew_3d.py

@@ -0,0 +1,108 @@
+from typing import Union
+
+from .._tier0 import plugin_function
+from .._tier0 import Image
+from .._tier0 import push
+from .._tier0 import create, create_like, create_none
+from ._AffineTransform3D import AffineTransform3D
+from skimage.transform import AffineTransform
+from ._affine_transform import _determine_translation_and_bounding_box
+import numpy as np
+from enum import Enum
+
+
+class DeskewDirection(Enum):
+    X = 1
+    Y = 2
+
+@plugin_function(output_creator=create_none)
+def affine_transform_deskew_3d(source: Image, destination: Image = None,
+                               transform: Union[np.ndarray, AffineTransform3D, AffineTransform] = None,
+                               deskewing_angle_in_degrees: float = 30,
+                               voxel_size_x: float = 0.1449922,
+                               voxel_size_y: float = 0.1449922,
+                               voxel_size_z: float = 0.3,
+                               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)
+
+    Parameters
+    ----------
+    source : Image
+        image to be transformed
+    destination : Image, optional
+        image where the transformed image should be written to
+    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
+    voxel_size_x: float, optional
+        Pixel size in X axis in microns
+    voxel_size_y: float, optional
+        Pixel size in Y axis in microns
+    voxel_size_z: float, optional
+        Step size between image planes along coverslip in microns; Voxel size in Z in microns
+    deskew_direction: str, optional
+        Direction of skew, dependent on microscope configuration
+    auto_size:bool, optional
+        If true, modifies the transform and the destination image size will be determined automatically, depending on the provided transform.
+        the transform might be modified so that all voxels of the result image have positions x>=0, y>=0, z>=0 and sit
+        tight to the coordinate origin. No voxels will cropped, the result image will fit in the returned destination.
+        Hence, the applied transform may have an additional translation vector that was not explicitly provided. This
+        also means that any given translation vector will be neglected.
+        If false, the destination image will have the same size as the input image.
+        Note: The value of auto-size is ignored if: destination is not None or transform is not an instance of
+        AffineTransform3D.
+
+    Returns
+    -------
+    destination
+
+    """
+    import numpy as np
+    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)}"
+
+    # handle output creation
+    new_size, transform, _ = _determine_translation_and_bounding_box(source, transform)
+    if destination is None:
+        destination = create(new_size)
+
+    # we invert the transform because we go from the target image to the source image to read pixels
+    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))
+    sintheta = np.float32(np.sin(deskewing_angle_in_degrees * np.pi/180))
+    costheta = np.float32(np.cos(deskewing_angle_in_degrees * np.pi/180))
+
+    gpu_transform_matrix = push(transform_matrix)
+
+    if deskew_direction == DeskewDirection.Y:
+        kernel_suffix = 'deskew_y_'
+        # change step size from physical space (nm) to camera space (pixels)
+        pixel_step = np.float32(voxel_size_z/voxel_size_y)
+    else:
+        kernel_suffix = 'deskew_x_'
+        # change step size from physical space (nm) to camera space (pixels)
+        pixel_step = np.float32(voxel_size_z/voxel_size_x)
+
+    # pass the shape of the final image, pixel step and precalculated trig
+
+    parameters = {
+        "input": source,
+        "output": destination,
+        "mat": gpu_transform_matrix,
+        "pixel_step": pixel_step,
+        "tantheta": tantheta,
+        "costheta": costheta,
+        "sintheta": sintheta
+    }
+
+    execute(__file__, './affine_transform_' + kernel_suffix + str(len(destination.shape)) + 'd' + '_x.cl',
+            'affine_transform_'  + kernel_suffix + str(len(destination.shape)) + 'd', destination.shape, parameters)
+
+    return destination

pyclesperanto_prototype/_tier8/_deskew_x.py

@@ -1,4 +1,4 @@
-from numpy import angle
+import warnings
 from .._tier0 import plugin_function
 from .._tier0 import Image
 from .._tier0 import create_none
@@ -12,10 +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 = False,
+             linear_interpolation: bool = True
              ) -> Image:
     """
-    Deskew an image stack as acquired with oblique plane light-sheet microscopy.
+    Deskew an image stack as acquired with oblique plane light-sheet microscopy, with skew in the X direction.
 
     Parameters
     ----------
@@ -28,28 +28,41 @@ def deskew_x(input_image: Image,
     voxel_size_x: float, optional
     voxel_size_y: float, optional
     voxel_size_z: float, optional
-         default: 1 micron
-         Voxel size, typically provided in microns
+        default: 1 micron
+        Voxel size, typically provided in microns
     scale_factor: float, optional
         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: bool, optional
-        If true, bi-/tri-linear interpolation will be applied, if hardware supports it.
-        If false, nearest-neighbor interpolation wille be applied.
+        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 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
-
-    # shear in the X plane towards Y
+    from ._affine_transform_deskew_3d import affine_transform_deskew_3d, DeskewDirection
+
+    # define affine transformation
     transform = AffineTransform3D()
     transform._deskew_x(angle_in_degrees=angle_in_degrees, voxel_size_x=voxel_size_x, voxel_size_y=voxel_size_y,
                         voxel_size_z=voxel_size_z, scale_factor=scale_factor)
 
     # apply transform
-    return affine_transform(input_image, output_image, transform=transform, auto_size=True,linear_interpolation=linear_interpolation)
+    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)

pyclesperanto_prototype/_tier8/_deskew_y.py

@@ -1,4 +1,4 @@
-from numpy import angle
+import warnings
 from .._tier0 import plugin_function
 from .._tier0 import Image
 from .._tier0 import create_none
@@ -12,10 +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 = False,
+             linear_interpolation: bool = True
              ) -> Image:
     """
-    Deskew an image stack as acquired with oblique plane light-sheet microscopy.
+    Deskew an image stack as acquired with oblique plane light-sheet microscopy with skew in the Y direction.
 
     Parameters
     ----------
@@ -28,35 +28,41 @@ def deskew_y(input_image: Image,
     voxel_size_x: float, optional
     voxel_size_y: float, optional
     voxel_size_z: float, optional
-         default: 1 micron
-         Voxel size, typically provided in microns
+        default: 1 micron
+        Voxel size, typically provided in microns
     scale_factor: float, optional
         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.
+        The isotropic voxel size of the output image will then be voxel_size_y / scaling_factor.
     linear_interpolation: bool, optional
-        If true, bi-/tri-linear interpolation will be applied, if hardware supports it.
-        If false, nearest-neighbor interpolation wille be applied.
+        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
     """
-
-    # this is a workaround step, see https://github.com/clEsperanto/pyclesperanto_prototype/issues/172
-    # resolved by defining shear factor based on voxel in _deskew_y()
-    #from ._scale import scale
-    #scaled_image = scale(input_image, factor_z=voxel_size_z / 0.3, auto_size=True)
-    #voxel_size_z = 0.3
+    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
-
-    # shear in the X plane towards Y
+    from ._affine_transform_deskew_3d import affine_transform_deskew_3d, DeskewDirection
+
+    # define affine transformation matrix
     transform = AffineTransform3D()
     transform._deskew_y(angle_in_degrees=angle_in_degrees, voxel_size_x=voxel_size_x, voxel_size_y=voxel_size_y,
                         voxel_size_z=voxel_size_z, scale_factor=scale_factor)
 
-    # apply transform
-    return affine_transform(input_image, output_image, transform=transform, auto_size=True,linear_interpolation=linear_interpolation)
-
+    # apply transform using special affine transform method for deskewing in the Y direction
+    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)