nanshe.imp.registration module

The registration module helps bring image stacks into alignment.

Overview

The registration algorithm developed by Wenzhi Sun. The strategy uses an area-based method for registration. Namely, it takes a template image (the mean projection) and finds the translations required for each frame to overlap optimally with the template image. Then a number of adjustments are made to the shifts to ensure that unnecessary translations are not preformed.

API

nanshe.imp.registration.find_offsets(*args, **kwargs)[source]

Computes the convolution of the template with the frames by taking advantage of their FFTs for faster computation that an ordinary convolution ( O(N*lg(N)) vs O(N^2) ) < http://ccrma.stanford.edu/~jos/ReviewFourier/FFT_Convolution_vs_Direct.html >. Once computed the maximum of the convolution is found to determine the best overlap of each frame with the template, which provides the needed offset. Some corrections are performed to make reasonable offsets.

Notes

Adapted from code provided by Wenzhi Sun with speed improvements provided by Uri Dubin.

Parameters:
  • frames2reg (numpy.ndarray) – image stack to register (time is the first dimension uses C-order tyx or tzyx).
  • template_fft (numpy.ndarray) – what to register the image stack against (single frame using C-order yx or zyx).
Returns:

an array containing the

translations to apply to each frame.

Return type:

(numpy.ndarray)

Examples

>>> a = numpy.zeros((5, 3, 4)); a[:,0] = 1; a[2,0] = 0; a[2,2] = 1
>>> a
array([[[ 1.,  1.,  1.,  1.],
        [ 0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.]],
<BLANKLINE>
       [[ 1.,  1.,  1.,  1.],
        [ 0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.]],
<BLANKLINE>
       [[ 0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.],
        [ 1.,  1.,  1.,  1.]],
<BLANKLINE>
       [[ 1.,  1.,  1.,  1.],
        [ 0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.]],
<BLANKLINE>
       [[ 1.,  1.,  1.,  1.],
        [ 0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.]]])
>>> af = numpy.fft.fftn(a, axes=range(1, a.ndim))
>>> af  # doctest: +SKIP
array([[[ 4.+0.j        ,  0.+0.j        ,  0.+0.j        ,  0.+0.j        ],
        [ 4.+0.j        ,  0.+0.j        ,  0.+0.j        ,  0.+0.j        ],
        [ 4.+0.j        ,  0.+0.j        ,  0.+0.j        ,  0.+0.j        ]],
<BLANKLINE>
       [[ 4.+0.j        ,  0.+0.j        ,  0.+0.j        ,  0.+0.j        ],
        [ 4.+0.j        ,  0.+0.j        ,  0.+0.j        ,  0.+0.j        ],
        [ 4.+0.j        ,  0.+0.j        ,  0.+0.j        ,  0.+0.j        ]],
<BLANKLINE>
       [[ 4.+0.j        ,  0.+0.j        ,  0.+0.j        ,  0.+0.j        ],
        [-2.+3.46410162j,  0.+0.j        ,  0.+0.j        ,  0.+0.j        ],
        [-2.-3.46410162j,  0.+0.j        ,  0.+0.j        ,  0.+0.j        ]],
<BLANKLINE>
       [[ 4.+0.j        ,  0.+0.j        ,  0.+0.j        ,  0.+0.j        ],
        [ 4.+0.j        ,  0.+0.j        ,  0.+0.j        ,  0.+0.j        ],
        [ 4.+0.j        ,  0.+0.j        ,  0.+0.j        ,  0.+0.j        ]],
<BLANKLINE>
       [[ 4.+0.j        ,  0.+0.j        ,  0.+0.j        ,  0.+0.j        ],
        [ 4.+0.j        ,  0.+0.j        ,  0.+0.j        ,  0.+0.j        ],
        [ 4.+0.j        ,  0.+0.j        ,  0.+0.j        ,  0.+0.j        ]]])
>>> tf = numpy.fft.fftn(a.mean(axis=0))
>>> tf  # doctest: +SKIP
array([[ 4.0+0.j        ,  0.0+0.j        ,  0.0+0.j        ,  0.0+0.j        ],
       [ 2.8+0.69282032j,  0.0+0.j        ,  0.0+0.j        ,  0.0+0.j        ],
       [ 2.8-0.69282032j,  0.0+0.j        ,  0.0+0.j        ,  0.0+0.j        ]])
>>> find_offsets(af, tf)
array([[ 0,  0],
       [ 0,  0],
       [-2,  0],
       [ 0,  0],
       [ 0,  0]])
nanshe.imp.registration.generate_unit_phase_shifts(*args, **kwargs)[source]

Computes the complex phase shift’s angle due to a unit spatial shift.

This is meant to be a helper function for register_mean_offsets. It does this by computing a table of the angle of the phase of a unit shift in each dimension (with a factor of \(2\pi\)).

This allows arbitrary phase shifts to be made in each dimensions by multiplying these angles by the size of the shift and added to the existing angle to induce the proper phase shift in fourier space, which is equivalent to the spatial translation.

Parameters:
  • shape (tuple of ints) – shape of the data to be shifted.
  • float_type (real type) – phase type (default numpy.float64)
Returns:

an array containing the angle of the

complex phase shift to use for each dimension.

Return type:

(numpy.ndarray)

Examples

>>> generate_unit_phase_shifts((2,4))
array([[[-0.        , -0.        , -0.        , -0.        ],
        [-3.14159265, -3.14159265, -3.14159265, -3.14159265]],
<BLANKLINE>
       [[-0.        , -1.57079633, -3.14159265, -4.71238898],
        [-0.        , -1.57079633, -3.14159265, -4.71238898]]])
nanshe.imp.registration.register_mean_offsets(*args, **kwargs)[source]

This algorithm registers the given image stack against its mean projection. This is done by computing translations needed to put each frame in alignment. Then the translation is performed and new translations are computed. This is repeated until no further improvement can be made.

The code for translations can be found in find_mean_offsets.

Notes

Adapted from code provided by Wenzhi Sun with speed improvements provided by Uri Dubin.

Parameters:
  • frames2reg (numpy.ndarray) – Image stack to register (time is the first dimension uses C-order tyx or tzyx).
  • max_iters (int) – Number of iterations to allow before forcing termination if stable point is not found yet. Set to -1 if no limit. (Default -1)
  • block_frame_length (int) – Number of frames to work with at a time. By default all. (Default -1)
  • include_shift (bool) – Whether to return the shifts used, as well. (Default False)
  • to_truncate (bool) – Whether to truncate the frames to remove all masked portions. (Default False)
  • float_type (type) – Type of float to use for calculation. (Default numpy.float64).
Returns:

an array containing the

translations to apply to each frame.

Return type:

(numpy.ndarray)

Examples

>>> a = numpy.zeros((5, 3, 4)); a[:,0] = 1; a[2,0] = 0; a[2,2] = 1
>>> a
array([[[ 1.,  1.,  1.,  1.],
        [ 0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.]],
<BLANKLINE>
       [[ 1.,  1.,  1.,  1.],
        [ 0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.]],
<BLANKLINE>
       [[ 0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.],
        [ 1.,  1.,  1.,  1.]],
<BLANKLINE>
       [[ 1.,  1.,  1.,  1.],
        [ 0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.]],
<BLANKLINE>
       [[ 1.,  1.,  1.,  1.],
        [ 0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.]]])
>>> register_mean_offsets(a, include_shift=True)
(masked_array(data =
 [[[1.0 1.0 1.0 1.0]
  [0.0 0.0 0.0 0.0]
  [0.0 0.0 0.0 0.0]]
<BLANKLINE>
 [[1.0 1.0 1.0 1.0]
  [0.0 0.0 0.0 0.0]
  [0.0 0.0 0.0 0.0]]
<BLANKLINE>
 [[-- -- -- --]
  [0.0 0.0 0.0 0.0]
  [0.0 0.0 0.0 0.0]]
<BLANKLINE>
 [[1.0 1.0 1.0 1.0]
  [0.0 0.0 0.0 0.0]
  [0.0 0.0 0.0 0.0]]
<BLANKLINE>
 [[1.0 1.0 1.0 1.0]
  [0.0 0.0 0.0 0.0]
  [0.0 0.0 0.0 0.0]]],
             mask =
 [[[False False False False]
  [False False False False]
  [False False False False]]
<BLANKLINE>
 [[False False False False]
  [False False False False]
  [False False False False]]
<BLANKLINE>
 [[ True  True  True  True]
  [False False False False]
  [False False False False]]
<BLANKLINE>
 [[False False False False]
  [False False False False]
  [False False False False]]
<BLANKLINE>
 [[False False False False]
  [False False False False]
  [False False False False]]],
       fill_value = 0.0)
, array([[0, 0],
       [0, 0],
       [1, 0],
       [0, 0],
       [0, 0]]))
nanshe.imp.registration.translate_fourier(*args, **kwargs)[source]

Translates frame(s) of data in Fourier space using the shift(s) given.

Parameters:
  • frame_fft (complex array) – Either a single frame with C-order axes or multiple frames with time on the 0th axis.
  • shift (array of ints) – Either the shift for each dimension with C-ordered values or multiple frames with time on the 0th axis.
Returns:

The frame(s) shifted.

Return type:

(numpy.ndarray)

Examples

>>> a = numpy.arange(12).reshape(3,4).astype(float)
>>> a
array([[  0.,   1.,   2.,   3.],
       [  4.,   5.,   6.,   7.],
       [  8.,   9.,  10.,  11.]])
>>> af = fft.fftn(a, axes=tuple(iters.irange(a.ndim)))
>>> numpy.around(af, decimals=10)  # doctest: +SKIP
array([[ 66. +0.j        ,  -6. +6.j        ,  -6. +0.j        ,  -6. -6.j        ],
       [-24.+13.85640646j,   0. +0.j        ,   0. +0.j        ,   0. +0.j        ],
       [-24.-13.85640646j,   0. +0.j        ,   0. +0.j        ,   0. +0.j        ]])
>>> s = numpy. array([1, -1])
>>> atf = translate_fourier(af, s)
>>> numpy.around(atf, decimals=10)  # doctest: +SKIP
array([[ 66. +0.j        ,  -6. -6.j        ,   6. -0.j        ,  -6. +6.j        ],
       [ 24.+13.85640646j,   0. +0.j        ,   0. +0.j        ,  -0. +0.j        ],
       [ 24.-13.85640646j,   0. -0.j        ,   0. +0.j        ,   0. +0.j        ]])
>>> fft.ifftn(
...     atf, axes=tuple(iters.irange(a.ndim))
... ).real.round().astype(int).astype(float)
array([[  9.,  10.,  11.,   8.],
       [  1.,   2.,   3.,   0.],
       [  5.,   6.,   7.,   4.]])
>>> a = a[None]; af = af[None]; s = s[None]
>>> atf = translate_fourier(af, s)
>>> numpy.around(atf, decimals=10)  # doctest: +SKIP
array([[[ 66. +0.j        ,  -6. -6.j        ,   6. -0.j        ,  -6. +6.j        ],
        [ 24.+13.85640646j,   0. +0.j        ,   0. +0.j        ,  -0. +0.j        ],
        [ 24.-13.85640646j,   0. -0.j        ,   0. +0.j        ,   0. +0.j        ]]])
>>> fft.ifftn(
...     atf, axes=tuple(iters.irange(1, a.ndim))
... ).real.round().astype(int).astype(float)
array([[[  9.,  10.,  11.,   8.],
        [  1.,   2.,   3.,   0.],
        [  5.,   6.,   7.,   4.]]])