Source code for atlalign.zoo

"""Contains different way how to generate displacement field.

Notes
-----
Ideally you want the only positional argument to be the shape and all the others be keyword arguments with
reasonable defaults.
"""

"""
    The package atlalign is a tool for registration of 2D images.

    Copyright (C) 2021 EPFL/Blue Brain Project

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with this program.  If not, see <https://www.gnu.org/licenses/>.
"""

import math

import cv2
import numpy as np
from scipy.interpolate import Rbf, SmoothBivariateSpline, griddata
from skimage.filters import gaussian
from skimage.transform import AffineTransform, ProjectiveTransform, SimilarityTransform

from atlalign.utils import griddata_custom


[docs]def affine(shape, matrix=None): """Affine transformation encoded in a 2 x 3 matrix. Parameters ---------- shape : tuple Of the form (height, width). matrix : np.ndarray Transformation matrix of the shape 2 x 3. Raises ------ ValueError In case the transformation matrix has a wrong shape. Returns ------- delta_x : np.ndarray Displacement vector field of the x coordinates. delta_y : np.ndarray Displacement vector field of the y coordinates. """ if matrix is None: matrix = np.eye(3) if matrix.shape == (2, 3): matrix = np.vstack( (matrix, [0, 0, 1]) ) # just add the homogeneous coordinates parts if matrix.shape != (3, 3): raise ValueError( "The shape of affine transformation matrix is {}, correct is (3, 3)".format( matrix.shape ) ) tform = AffineTransform(matrix) x, y = np.meshgrid(range(shape[1]), range(shape[0])) coords = np.hstack((x.reshape(-1, 1), y.reshape(-1, 1))).astype(float) coords_after = tform(coords) coords_delta = coords_after - coords delta_x = np.reshape(coords_delta[:, 0], shape) delta_y = np.reshape(coords_delta[:, 1], shape) return delta_x, delta_y
[docs]def affine_simple( shape, scale_x=1, scale_y=1, rotation=0, translation_x=0, translation_y=0, shear=0, apply_centering=True, ): """Just a human version of affine mapping. Notes ----- Instead of specifying the whole matrix one can just specify all the understandable quantities. Parameters ---------- shape : tuple Of the form (height, width). scale_x : float Scale on the x axis. If scale_x < 1 then zoom out, if scale_x > 1 zoom in. scale_y : float Scale on the y axis. If scale_y < 1 then zoom out, if scale_y > 1 zoom in. rotation : float Rotation angle in counter-clockwise direction as radians. translation_x : float Translation in the x direction. If translation_x > 0 then to the right, else to the left. translation_y : float Translation in the y direction. If translation_y > 0 then down, else to the up. shear : float Shear angle in counter-clockwise direction as radians. apply_centering : bool If True then (h // 2 - 0.5, w // 2 - 0.5) is considered a center of the image. And before performing all the other operations the image is first shifted so that the center corresponds to (0, 0). Then the actual transformation is applied and after that the image is shifted into the original center. Returns ------- delta_x : np.ndarray Displacement vector field of the x coordinates. delta_y : np.ndarray Displacement vector field of the y coordinates. """ matrix = np.array( [ [ scale_x * math.cos(rotation), -scale_y * math.sin(rotation + shear), translation_x, ], [ scale_x * math.sin(rotation), scale_y * math.cos(rotation + shear), translation_y, ], [0, 0, 1], ] ) if not apply_centering: return affine(shape, matrix) center_rc = np.array([(shape[0] / 2) - 0.5, (shape[1] / 2) - 0.5]) center_xy = np.array([center_rc[1], center_rc[0]]) tform1 = SimilarityTransform(translation=center_xy) tform2 = SimilarityTransform(matrix=matrix) tform3 = SimilarityTransform(translation=-center_xy) tform = tform3 + tform2 + tform1 return affine(shape, tform.params)
[docs]def control_points( shape, points=None, values_delta_x=None, values_delta_y=None, anchor_corners=True, interpolation_method="griddata", interpolator_kwargs=None, ): """Simply interpolate given control points. Notes ----- We assume there are N control points. This function is used by others functions from the `zoo`. See below a complete list - `edge_stretching` - `single_frequency` Additionally, it is also used in the `eliminate_bb` method of the ``DisplacementField`` class. Parameters ---------- shape : tuple Of the form (height, width). points : np.ndarray, optional An array of shape (N, 2) where each row represents a (row, column) of a given control point. values_delta_x : np.ndarray, optional An array of shape (N, ) where each row represents a delta_x of the transformation at the corresponding control point. values_delta_y : np.ndarray, optional An array of shape (N, ) where each row represents a delta_y of the transformation at the corresponding control point. anchor_corners : bool, optional If True then each of the 4 images corners are automatically added to the control points and identity transformation is assumed. interpolation_method : {'griddata', 'bspline', 'rbf'}, optional Interpolation method to use. interpolator_kwargs : dict, optional Additional parameters passed to the interpolator. Returns ------- delta_x : np.ndarray Displacement vector field of the x coordinates of shape = `shape`. delta_y : np.ndarray Displacement vector field of the y coordinates of shape = `shape`. Raises ------ TypeError: When either of the `points`, `values_delta_x` or `values_delta_y` is not a np.ndarray. ValueError: Various inconsistencies in inputs (different len, zero len, wrong other dimensions or default without anchor). IndexError: Some of the points are outside of the image domain. """ height, width = shape interpolator_kwargs = interpolator_kwargs or {} in_default_mode = ( points is None and values_delta_x is None and values_delta_y is None ) # Create a default if in_default_mode: if anchor_corners: points = np.array( [(0, 0), (0, width - 1), (height - 1, width - 1), (height - 1, 0)] ) # clockwise values_delta_x = np.zeros(4) values_delta_y = np.zeros(4) else: raise ValueError("Cannot instantiate") # Checks if not ( isinstance(points, np.ndarray) and isinstance(values_delta_x, np.ndarray) and isinstance(values_delta_y, np.ndarray) ): raise TypeError("All inputs need to be a np.ndarray") if not (len(points) == len(values_delta_x) == len(values_delta_y)): raise ValueError("Arrays have different lengths") if len(points) == 0: raise ValueError("No points given.") N = len(points) if not (points.ndim == 2 and values_delta_x.ndim == 1 and values_delta_y.ndim == 1): raise ValueError("Wrong dimensions of input data") if not np.all([0 <= p[0] < height and 0 <= p[1] < width for p in points]): raise IndexError("Some control points are outside of image domain.") if not in_default_mode and anchor_corners: points_anchors = np.array( [(0, 0), (0, width - 1), (height - 1, width - 1), (height - 1, 0)] ) # clockwise values_delta_x_anchors = np.zeros((4, 1)) values_delta_y_anchors = np.zeros((4, 1)) points = np.vstack((points_anchors, points)) values_delta_x = np.vstack( (values_delta_x_anchors, values_delta_x.reshape((N, 1))) ) values_delta_y = np.vstack( (values_delta_y_anchors, values_delta_y.reshape((N, 1))) ) x, y = np.meshgrid(list(range(width)), list(range(height))) xi = (y, x) if interpolation_method == "griddata": values_grid_delta_x = griddata( points=points, values=values_delta_x, xi=xi, **interpolator_kwargs ) values_grid_delta_y = griddata( points=points, values=values_delta_y, xi=xi, **interpolator_kwargs ) delta_x = values_grid_delta_x.reshape(shape) delta_y = values_grid_delta_y.reshape(shape) elif interpolation_method == "griddata_custom": # triangulation performed only once values_grid_x, values_grid_y = griddata_custom( points, values_delta_x, values_delta_y, xi ) delta_x = values_grid_x.reshape(shape) delta_y = values_grid_y.reshape(shape) elif interpolation_method == "bspline": x_, y_ = points[:, 1], points[:, 0] ip_delta_x = SmoothBivariateSpline( x_, y_, values_delta_x, **interpolator_kwargs ) ip_delta_y = SmoothBivariateSpline( x_, y_, values_delta_y, **interpolator_kwargs ) delta_x = ip_delta_x(x.ravel(), y.ravel(), grid=False).reshape(shape) delta_y = ip_delta_y(x.ravel(), y.ravel(), grid=False).reshape(shape) elif interpolation_method == "rbf": x_, y_ = points[:, 1], points[:, 0] ip_delta_x = Rbf(x_, y_, values_delta_x, **interpolator_kwargs) ip_delta_y = Rbf(x_, y_, values_delta_y, **interpolator_kwargs) delta_x = ip_delta_x(x.ravel(), y.ravel()).reshape(shape) delta_y = ip_delta_y(x.ravel(), y.ravel()).reshape(shape) else: raise ValueError( "Unrecognized interpolation method {}".format(interpolation_method) ) return delta_x, delta_y
[docs]def edge_stretching( shape, edge_mask=None, n_perturbation_points=3, radius_max=30, interpolation_method="griddata_custom", interpolator_kwargs=None, ): """Pick points on the edges and using them to stretch the image. Parameters ---------- shape : tuple Of the form (height, width). edge_mask : np.ndarray An array of dtype=bool of shape `shape`. The True elements represent an edge. n_perturbation_points : int Number of points to pick among the edges on which the perturbation defined. radius_max : float, optional Maxim value of radius, the actual value is a sample from uniform [0, radius_max]. interpolation_method : {'griddata', 'griddata_custom', 'bspline', 'rbf'}, optional Interpolation method to use. interpolator_kwargs : dict, optional Additional parameters passed to the interpolator. Returns ------- delta_x : np.ndarray Displacement vector field of the x coordinates. delta_y : np.ndarray Displacement vector field of the y coordinates. """ if edge_mask is None: edge_mask = np.zeros(shape, dtype=bool) # no edges if shape != edge_mask.shape: raise ValueError("The edge_mask has an incorrect shape.") n_edge_pixels = edge_mask.sum() n_perturbation_points = min(n_perturbation_points, n_edge_pixels) if n_perturbation_points == 0: return np.zeros(shape), np.zeros( shape ) # no perturbation points -> identity mapping edge_points = np.argwhere(edge_mask) ixs = np.random.choice(n_edge_pixels, replace=False, size=n_perturbation_points) points = edge_points[ixs] perturbation_radii = np.random.uniform(0, radius_max, size=n_perturbation_points) perturbation_angles = np.random.uniform( 0, 2 * np.pi, size=n_perturbation_points ) # in radians pixel_displacements_rc = np.array( [ np.array([-np.sin(angle), np.cos(angle)]) * radius for radius, angle in zip(perturbation_radii, perturbation_angles) ] ) return control_points( shape=shape, points=points, values_delta_x=pixel_displacements_rc[:, 0], values_delta_y=pixel_displacements_rc[:, 1], interpolation_method=interpolation_method, interpolator_kwargs=interpolator_kwargs, )
[docs]def paper( shape, n_pixels=10, v_min=-20000, v_max=20000, kernel_sigma=25, p=None, random_state=None, ): """Algorithm proposed in the reference paper. Notes ----- This algorithm has 2 steps 1) Pick `n_pixels` in the image and randomly sample x and y displacement from interval [`v_min`, `v_max`] 2) Apply a gaussian kernel on the displacement field with a given `kernel_size` and `kernel_sigma` Parameters ---------- shape : tuple Of the form (height, width). n_pixels : int Number of pixels to choose in the first step. v_min : float Minimum value for x and y displacement sampling. v_max : float Maximum value for x and y displacement sampling. kernel_sigma : float Standard deviation of the kernel in both the x and y direction. p : np.array Pixelwise probability of selection, where p.shape=shape. Note that if None, then uniform. random_state : int If None, then results not reproducible. Returns ------- delta_x : np.ndarray Displacement vector field of the x coordinates. delta_y : np.ndarray Displacement vector field of the y coordinates. References ---------- [1] Sokooti H., de Vos B., Berendsen F., Lelieveldt B.P.F., IĆĄgum I., Staring M. (2017) Nonrigid Image Registration Using Multi-scale 3D Convolutional Neural Network """ # Checks if n_pixels <= 0: raise ValueError( "The n_pixels needs to be higher than 0, current value is {}".format( n_pixels ) ) height, width = shape delta_x = np.zeros( shape, dtype=np.float32 ) # makes GaussianBlur way quicker than float64 delta_y = np.zeros( shape, dtype=np.float32 ) # makes GaussianBlur way quicker than float64 # Step 1 np.random.seed(random_state) # set seed pixel_count = np.prod(shape) if p is not None: p_ = p.ravel() else: p_ = None selected_pixels = np.random.choice(pixel_count, n_pixels, replace=False, p=p_) selected_pixels_rc = np.array([(x // width, x % width) for x in selected_pixels]) xy_deltas = np.random.uniform(v_min, v_max, size=(n_pixels, 2)) delta_y[selected_pixels_rc[:, 0], selected_pixels_rc[:, 1]] = xy_deltas[:, 1] delta_x[selected_pixels_rc[:, 0], selected_pixels_rc[:, 1]] = xy_deltas[:, 0] # Step 2 delta_x = cv2.GaussianBlur( delta_x, ksize=(0, 0), sigmaX=kernel_sigma, sigmaY=kernel_sigma ) delta_y = cv2.GaussianBlur( delta_y, ksize=(0, 0), sigmaX=kernel_sigma, sigmaY=kernel_sigma ) return delta_x, delta_y
[docs]def paper_microsoft(shape, alpha=1000, sigma=10, random_state=None): """Generate artificial displacement based on a paper from Microsoft. Parameters ---------- shape : tuple Of the form (height, width). alpha : float Constant that the per pixel displacements are multiplied with. The higher the crazier displacements. If set to 0 then zero transformation. sigma : float Standard deviation of the gaussian kernel. The closer to 0 the crazier displacement. If close to inf then zero transformation. Returns ------- delta_x : np.ndarray Displacement vector field of the x coordinates. delta_y : np.ndarray Displacement vector field of the y coordinates. References ---------- Simard, P. Y., Steinkraus, D., & Platt, J. C. (2003, August). Best practices for convolutional neural networks applied to visual document analysis. In null (p. 958). IEEE. """ np.random.seed(random_state) random_delta_x_, random_delta_y_ = np.random.uniform(-1, 1, size=(2, *shape)) delta_x = gaussian(random_delta_x_, sigma=sigma, cval=0) * alpha delta_y = gaussian(random_delta_y_, sigma=sigma, cval=0) * alpha return delta_x, delta_y
[docs]def patch_shift( shape, ul=(10, 10), height=100, width=120, shift_size=30, shift_direction="D" ): """For a fixed patch in an image redefine it with another same-shaped patch elsewhere in the image. Parameters ---------- shape : tuple Of the form (height, width). ul : tuple Of the form (row of the UPPER LEFT corner of the patch, column of the UPPER LEFT corner of the patch). height : int Height of the patch. width : int Width of the patch. shift_size : int How many pixels to shift the patch. shift_direction : str, {'U', 'D', 'L', 'R'} The direction of the shift. 'U' = Up, 'D' = Down, 'L' = Left, 'R' = Right. Raises ------ IndexError: If the starting or the ending patch are not in the image. Returns ------- delta_x : np.ndarray Displacement vector field of the x coordinates. delta_y : np.ndarray Displacement vector field of the y coordinates. """ delta_x = np.zeros(shape) delta_y = np.zeros(shape) if shift_direction == "R": delta_x[ul[0] : ul[0] + height, ul[1] : ul[1] + width] = shift_size elif shift_direction == "L": delta_x[ul[0] : ul[0] + height, ul[1] : ul[1] + width] = -shift_size elif shift_direction == "U": delta_y[ul[0] : ul[0] + height, ul[1] : ul[1] + width] = -shift_size elif shift_direction == "D": delta_y[ul[0] : ul[0] + height, ul[1] : ul[1] + width] = shift_size else: raise ValueError("Unrecognized shift direction {}".format(shift_direction)) return delta_x, delta_y
[docs]def projective(shape, matrix=None): """Projective transformation encoded in a 3 x 3 matrix. Parameters ---------- shape : tuple Of the form (height, width). matrix : np.ndarray Transformation matrix of the shape 3 x 3. Raises ------ ValueError In case the transformation matrix has a wrong shape. Returns ------- delta_x : np.ndarray Displacement vector field of the x coordinates. delta_y : np.ndarray Displacement vector field of the y coordinates. """ if matrix is None: matrix = np.eye(3) if matrix.shape != (3, 3): raise ValueError( "The shape of projective transformation matrix is {}, correct is (3, 3)".format( matrix.shape ) ) tform = ProjectiveTransform(matrix) x, y = np.meshgrid(range(shape[1]), range(shape[0])) coords = np.hstack((x.reshape(-1, 1), y.reshape(-1, 1))).astype(float) coords_after = tform(coords) coords_delta = coords_after - coords delta_x = np.reshape(coords_delta[:, 0], shape) delta_y = np.reshape(coords_delta[:, 1], shape) return delta_x, delta_y
[docs]def single_frequency( shape, p=None, grid_spacing=5, n_perturbation_points=20, radius_mean=None, interpolation_method="griddata", interpolator_kwargs=None, ): """Single frequency artificial warping generator. Notes ----- 1) The reason why this approach is called single frequency is that the grid_spacing is constant over the entire image regions. One can therefore capture displacements of more or less the same size pixelwise. 2) Calls the `control_points` so in this sense it is slightly unusual. 3) The idea is to create a regular grid and then sample some of the points on it. For a fixed point we then randomly select and angle - Uniform[0, 2pi] and also the diameter ~ Exp(`radius_mean`) Parameters ---------- shape : tuple Of the form (height, width). p : np.array, optional Pixelwise probability of selection, where p.shape=shape. Note that if None, then uniform. grid_spacing : int, optional Grid spacing size in both the columns and rows. n_perturbation_points : int, optional Number of grid points to which random perturbation will be applied (without replacement). radius_mean : float, optional If None then set to `grid_spacing / 2`. interpolation_method : {'griddata', 'bspline', 'rbf'}, optional Interpolation method to use. interpolator_kwargs : dict, optional Additional parameters passed to the interpolator. Returns ------- delta_x : np.ndarray Displacement vector field of the x coordinates. delta_y : np.ndarray Displacement vector field of the y coordinates. References ---------- [1] https://github.com/hsokooti/RegNet (README.md) """ height, width = shape pixel_count = np.prod(shape) radius_mean = grid_spacing / 2 if radius_mean is None else radius_mean # Create grid mask (True = on grid, False = off grid) ixs_x, ixs_y = ( list(range(width))[::grid_spacing], list(range(height))[::grid_spacing], ) grid_mask = np.zeros(shape, dtype=bool) grid_mask[ixs_y, :] = True grid_mask[:, ixs_x] = True grid_pixel_count = np.sum(grid_mask) if grid_pixel_count < n_perturbation_points: raise ValueError("Cannot have more perturbation points than grid points.") if p is not None: if not p.sum() == 1: raise ValueError("The p needs to sum up to one or be None.") p_c = p.copy() sum_offgrid = np.sum(p_c[~grid_mask]) p_c[~grid_mask] = 0 # make it impossible to choose off grid elements p_c[grid_mask] += sum_offgrid / np.sum(grid_mask) else: p_c = np.zeros(shape) p_c[grid_mask] = 1 / np.sum(grid_mask) p_ = p_c.ravel() selected_pixels = np.random.choice( pixel_count, n_perturbation_points, replace=False, p=p_ ) selected_pixels_rc = np.array([(x // width, x % width) for x in selected_pixels]) perturbation_radii = np.random.exponential(radius_mean, n_perturbation_points) perturbation_angles = np.random.uniform( 0, 2 * np.pi, n_perturbation_points ) # in radians pixel_displacements_rc = np.array( [ np.array([-np.sin(angle), np.cos(angle)]) * radius for radius, angle in zip(perturbation_radii, perturbation_angles) ] ) return control_points( shape=shape, points=selected_pixels_rc, values_delta_x=pixel_displacements_rc[:, 0], values_delta_y=pixel_displacements_rc[:, 1], interpolation_method=interpolation_method, interpolator_kwargs=interpolator_kwargs, )