WIP filtering

brainglobe_registration/elastix/register.py

@@ -5,6 +5,8 @@
 import numpy as np
 from brainglobe_atlasapi import BrainGlobeAtlas
 
+from brainglobe_registration.utils.utils import filter_image
+
 
 def get_atlas_by_name(atlas_name: str) -> BrainGlobeAtlas:
     """
@@ -60,8 +62,14 @@ def run_registration(
     npt.NDArray
         The inverse moving image.
     """
+    # Filter images
+    filtered_atlas_image = filter_image(atlas_image)
+    # filtered_moving_image = filter_image(moving_image)
+
     # convert to ITK, view only
-    atlas_image_elastix = itk.GetImageViewFromArray(atlas_image).astype(itk.F)
+    atlas_image_elastix = itk.GetImageViewFromArray(
+        filtered_atlas_image
+    ).astype(itk.F)
     moving_image_elastix = itk.GetImageViewFromArray(moving_image).astype(
         itk.F
     )
@@ -134,7 +142,7 @@ def run_registration(
     )
 
     inverse_moving = itk.transformix_filter(
-        moving_image_elastix,
+        itk.GetImageViewFromArray(moving_image).astype(itk.F),
         inverse_transform_parameters,
     )
 

brainglobe_registration/registration_widget.py

@@ -28,6 +28,7 @@
     QPushButton,
     QTabWidget,
 )
+from scipy.ndimage import gaussian_filter
 from skimage.segmentation import find_boundaries
 from skimage.transform import rescale
 
@@ -407,16 +408,22 @@ def _on_scale_moving_image(self, x: float, y: float):
         if self._moving_image_data_backup is None:
             self._moving_image_data_backup = self._moving_image.data.copy()
 
+        filtered_image = gaussian_filter(
+            self._moving_image_data_backup, sigma=20
+        )
+
         x_factor = x / self._atlas.resolution[0]
         y_factor = y / self._atlas.resolution[1]
 
         self._moving_image.data = rescale(
-            self._moving_image_data_backup,
+            filtered_image,
             (y_factor, x_factor),
             mode="constant",
             preserve_range=True,
             anti_aliasing=True,
         )
+        #
+        # self._moving_image.data = filtered_image
 
     def _on_adjust_atlas_rotation(self, pitch: float, yaw: float, roll: float):
         if not (

brainglobe_registration/utils/utils.py

@@ -5,6 +5,8 @@
 import numpy as np
 import numpy.typing as npt
 from pytransform3d.rotations import active_matrix_from_angle
+from scipy.ndimage import gaussian_filter
+from skimage import morphology
 
 
 def adjust_napari_image_layer(
@@ -145,12 +147,87 @@ def calculate_rotated_bounding_box(
         ]
     )
 
-    transformed_corners = np.dot(rotation_matrix, corners.T)
-    min_corner = np.min(transformed_corners, axis=1)
-    max_corner = np.max(transformed_corners, axis=1)
+    transformed_corners = np.dot(corners, rotation_matrix.T)
+    min_corner = np.min(transformed_corners, axis=0)
+    max_corner = np.max(transformed_corners, axis=0)
 
     return (
-        int(np.round(max_corner[0] - min_corner[0])),
-        int(np.round(max_corner[1] - min_corner[1])),
-        int(np.round(max_corner[2] - min_corner[2])),
+        int(np.ceil(max_corner[0] - min_corner[0])),
+        int(np.ceil(max_corner[1] - min_corner[1])),
+        int(np.ceil(max_corner[2] - min_corner[2])),
     )
+
+
+def filter_image(image, preprocessing_args=None) -> npt.NDArray:
+    """
+    Filter a 2D image prior to registration
+
+    Parameters
+    ----------
+    image : npt.NDArray
+        The image to filter
+    preprocessing_args : Dict, optional
+        The preprocessing arguments, by default None
+
+    Returns
+    -------
+    npt.NDArray
+        The filtered image
+    """
+    if preprocessing_args and preprocessing_args.preprocessing == "skip":
+        pass
+    else:  # default pre-processing
+        image = filter_plane(image)
+
+    return image
+
+
+def filter_plane(img_plane):
+    """
+    Apply a set of filter to the plane (typically to avoid overfitting details
+    in the image during registration)
+    The filter is composed of a despeckle filter using opening and a pseudo
+    flatfield filter
+
+    Parameters
+    ----------
+    img_plane : npt.NDArray
+        A 2D array to filter
+
+    Returns
+    -------
+    npt.NDArray
+        The filtered array
+    """
+    # img_plane = despeckle_by_opening(img_plane)
+    img_plane = gaussian_filter(img_plane, 2)
+    return img_plane
+
+
+def despeckle_by_opening(img_plane, radius=2):
+    """
+    Despeckle the image plane using a grayscale opening operation
+
+    :param np.array img_plane:
+    :param int radius: The radius of the opening kernel
+    :return: The despeckled image
+    :rtype: np.array
+    """
+    kernel = morphology.disk(radius)
+    morphology.opening(img_plane, out=img_plane, footprint=kernel)
+    return img_plane
+
+
+def pseudo_flatfield(img_plane, sigma=5):
+    """
+    Pseudo flat field filter implementation using a de-trending by a
+    heavily gaussian filtered copy of the image.
+
+    :param np.array img_plane: The image to filter
+    :param int sigma: The sigma of the gaussian filter applied to the
+        image used for de-trending
+    :return: The pseudo flat field filtered image
+    :rtype: np.array
+    """
+    filtered_img = gaussian_filter(img_plane, sigma)
+    return img_plane / (filtered_img + 1)