macromax.utils.ft package

Submodules

macromax.utils.ft.ft_implementation module

The ft module provides direct access to numpy.fft, scipy.fftpack, pyfftw, or mkl_fft, depending on availability. The back-end that is deemed most efficient is automatically imported. Independently of the back-end, the ft module provides at least the following functions:

ft.fft(a: array_like, axis: int): The 1-dimensional fast Fourier transform.
ft.ifft(a: array_like, axis: int): The 1-dimensional fast Fourier transform.
ft.fftn(a: array_like, axes: Sequence): The n-dimensional fast Fourier transform.
ft.ifftn(a: array_like, axes: Sequence): The n-dimensional inverse fast Fourier transform.
ft.fftshift(a: array_like, axes: Sequence): This fftshifts the input array.
ft.ifftshift(a: array_like, axes: Sequence): This ifftshifts the input array (only different from ft.fftshift for odd-shapes).

All functions return a numpy.ndarray of the same shape as the input array or array_like, a. With the exception of ft.fft() and ft.ifft(), all functions take the axes argument to limit the action to specific axes of the numpy.ndarray.

Note that axis indices should be unique and non-negative. Negative or repeated axis indices are not compatible with all back-end implementations!

macromax.utils.ft.ft_implementation.fftshift(x, axes=None)

Shift the zero-frequency component to the center of the spectrum.

This function swaps half-spaces for all axes listed (defaults to all). Note that y[0] is the Nyquist component only if len(x) is even.

Parameters:

x (array_like) – Input array.
axes (int or shape tuple, optional) – Axes over which to shift. Default is None, which shifts all axes.

Returns:

y – The shifted array.

Return type:

ndarray

See also

ifftshift: The inverse of fftshift.

Examples

>>> freqs = np.fft.fftfreq(10, 0.1)
>>> freqs
array([ 0.,  1.,  2., ..., -3., -2., -1.])
>>> np.fft.fftshift(freqs)
array([-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.])

Shift the zero-frequency component only along the second axis:

>>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
>>> freqs
array([[ 0.,  1.,  2.],
       [ 3.,  4., -4.],
       [-3., -2., -1.]])
>>> np.fft.fftshift(freqs, axes=(1,))
array([[ 2.,  0.,  1.],
       [-4.,  3.,  4.],
       [-1., -3., -2.]])

macromax.utils.ft.ft_implementation.ifftshift(x, axes=None)

The inverse of fftshift. Although identical for even-length x, the functions differ by one sample for odd-length x.

Parameters:

x (array_like) – Input array.
axes (int or shape tuple, optional) – Axes over which to calculate. Defaults to None, which shifts all axes.

Returns:

y – The shifted array.

Return type:

ndarray

See also

fftshift: Shift zero-frequency component to the center of the spectrum.

Examples

>>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
>>> freqs
array([[ 0.,  1.,  2.],
       [ 3.,  4., -4.],
       [-3., -2., -1.]])
>>> np.fft.ifftshift(np.fft.fftshift(freqs))
array([[ 0.,  1.,  2.],
       [ 3.,  4., -4.],
       [-3., -2., -1.]])

macromax.utils.ft.ft_implementation.fft(x, n=None, axis=-1, overwrite_x=False)[source]

Return discrete Fourier transform of real or complex sequence.

The returned complex array contains y(0), y(1),..., y(n-1), where

y(j) = (x * exp(-2*pi*sqrt(-1)*j*np.arange(n)/n)).sum().

Parameters:

x (array_like) – Array to Fourier transform.
n (int, optional) – Length of the Fourier transform. If n < x.shape[axis], x is truncated. If n > x.shape[axis], x is zero-padded. The default results in n = x.shape[axis].
axis (int, optional) – Axis along which the fft’s are computed; the default is over the last axis (i.e., axis=-1).
overwrite_x (bool, optional) – If True, the contents of x can be destroyed; the default is False.

Returns:

z –

with the elements:

[y(0),y(1),..,y(n/2),y(1-n/2),...,y(-1)]        if n is even
[y(0),y(1),..,y((n-1)/2),y(-(n-1)/2),...,y(-1)]  if n is odd

where:

y(j) = sum[k=0..n-1] x[k] * exp(-sqrt(-1)*j*k* 2*pi/n), j = 0..n-1

Return type:

complex ndarray

See also

ifft: Inverse FFT
rfft: FFT of a real sequence

Notes

The packing of the result is “standard”: If A = fft(a, n), then A[0] contains the zero-frequency term, A[1:n/2] contains the positive-frequency terms, and A[n/2:] contains the negative-frequency terms, in order of decreasingly negative frequency. So ,for an 8-point transform, the frequencies of the result are [0, 1, 2, 3, -4, -3, -2, -1]. To rearrange the fft output so that the zero-frequency component is centered, like [-4, -3, -2, -1, 0, 1, 2, 3], use fftshift.

Both single and double precision routines are implemented. Half precision inputs will be converted to single precision. Non-floating-point inputs will be converted to double precision. Long-double precision inputs are not supported.

This function is most efficient when n is a power of two, and least efficient when n is prime.

Note that if x is real-valued, then A[j] == A[n-j].conjugate(). If x is real-valued and n is even, then A[n/2] is real.

If the data type of x is real, a “real FFT” algorithm is automatically used, which roughly halves the computation time. To increase efficiency a little further, use rfft, which does the same calculation, but only outputs half of the symmetrical spectrum. If the data is both real and symmetrical, the dct can again double the efficiency by generating half of the spectrum from half of the signal.

Examples

>>> from scipy.fftpack import fft, ifft
>>> x = np.arange(5)
>>> np.allclose(fft(ifft(x)), x, atol=1e-15)  # within numerical accuracy.
True

macromax.utils.ft.ft_implementation.ifft(x, n=None, axis=-1, overwrite_x=False)[source]

Return discrete inverse Fourier transform of real or complex sequence.

The returned complex array contains y(0), y(1),..., y(n-1), where

y(j) = (x * exp(2*pi*sqrt(-1)*j*np.arange(n)/n)).mean().

Parameters:

x (array_like) – Transformed data to invert.
n (int, optional) – Length of the inverse Fourier transform. If n < x.shape[axis], x is truncated. If n > x.shape[axis], x is zero-padded. The default results in n = x.shape[axis].
axis (int, optional) – Axis along which the ifft’s are computed; the default is over the last axis (i.e., axis=-1).
overwrite_x (bool, optional) – If True, the contents of x can be destroyed; the default is False.

Returns:

ifft – The inverse discrete Fourier transform.

Return type:

ndarray of floats

See also

fft: Forward FFT

Notes

Both single and double precision routines are implemented. Half precision inputs will be converted to single precision. Non-floating-point inputs will be converted to double precision. Long-double precision inputs are not supported.

This function is most efficient when n is a power of two, and least efficient when n is prime.

If the data type of x is real, a “real IFFT” algorithm is automatically used, which roughly halves the computation time.

Examples

>>> from scipy.fftpack import fft, ifft
>>> import numpy as np
>>> x = np.arange(5)
>>> np.allclose(ifft(fft(x)), x, atol=1e-15)  # within numerical accuracy.
True

macromax.utils.ft.ft_implementation.fftn(x, shape=None, axes=None, overwrite_x=False)[source]

Return multidimensional discrete Fourier transform.

The returned array contains:

y[j_1,..,j_d] = sum[k_1=0..n_1-1, ..., k_d=0..n_d-1]
   x[k_1,..,k_d] * prod[i=1..d] exp(-sqrt(-1)*2*pi/n_i * j_i * k_i)

where d = len(x.shape) and n = x.shape.

Parameters:

x (array_like) – The (N-D) array to transform.
shape (int or array_like of ints or None, optional) – The shape of the result. If both shape and axes (see below) are None, shape is x.shape; if shape is None but axes is not None, then shape is numpy.take(x.shape, axes, axis=0). If shape[i] > x.shape[i], the ith dimension is padded with zeros. If shape[i] < x.shape[i], the ith dimension is truncated to length shape[i]. If any element of shape is -1, the size of the corresponding dimension of x is used.
axes (int or array_like of ints or None, optional) – The axes of x (y if shape is not None) along which the transform is applied. The default is over all axes.
overwrite_x (bool, optional) – If True, the contents of x can be destroyed. Default is False.

Returns:

y – The (N-D) DFT of the input array.

Return type:

complex-valued N-D NumPy array

See also

ifftn

Notes

If x is real-valued, then y[..., j_i, ...] == y[..., n_i-j_i, ...].conjugate().

Both single and double precision routines are implemented. Half precision inputs will be converted to single precision. Non-floating-point inputs will be converted to double precision. Long-double precision inputs are not supported.

Examples

>>> from scipy.fftpack import fftn, ifftn
>>> y = (-np.arange(16), 8 - np.arange(16), np.arange(16))
>>> np.allclose(y, fftn(ifftn(y)))
True

macromax.utils.ft.ft_implementation.ifftn(x, shape=None, axes=None, overwrite_x=False)[source]

Return inverse multidimensional discrete Fourier transform.

The sequence can be of an arbitrary type.

The returned array contains:

y[j_1,..,j_d] = 1/p * sum[k_1=0..n_1-1, ..., k_d=0..n_d-1]
   x[k_1,..,k_d] * prod[i=1..d] exp(sqrt(-1)*2*pi/n_i * j_i * k_i)

where d = len(x.shape), n = x.shape, and p = prod[i=1..d] n_i.

For description of parameters see fftn.

See also

fftn: for detailed information.

Examples

>>> from scipy.fftpack import fftn, ifftn
>>> import numpy as np
>>> y = (-np.arange(16), 8 - np.arange(16), np.arange(16))
>>> np.allclose(y, ifftn(fftn(y)))
True

macromax.utils.ft.grid module

class macromax.utils.ft.grid.Grid(shape=None, step=None, extent=None, first=None, center=None, last=None, include_last=False, ndim=None, flat=False, origin_at_center=True, center_at_index=True)[source]

Bases: Sequence

A class representing an immutable uniformly-spaced plaid Cartesian grid and its Fourier Transform.

macromax.utils.ft.subpixel module

Classes and functions to register a subject array to a reference array with subpixel precision. This is based on the algorithm described in: Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup, “Efficient subpixel image registration algorithms,” Optics Letters. 33, 156-158 (2008).

macromax.utils.ft.subpixel.register(subject=None, reference=None, precision=None)[source]

Registers a subject array to a reference array with subpixel precision. This is based on the algorithm described in Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup, “Efficient subpixel image registration algorithms,” Optics Letters. 33, 156-158 (2008).

Parameters:

subject (Union[int, float, Sequence, ndarray, None]) – The subject image array (n-dimensional) .
reference (Union[int, float, Sequence, ndarray, None]) – The optional reference array (n-dimensional).
precision (Optional[float]) – The registration precision in units of fracional pixels (default).

Return type:

Registration

Returns:

A Registration instance describing the registered image, the shift, and the amplitude.

macromax.utils.ft.subpixel.roll(subject, shift, axis=None)[source]

Rolls (shifting with wrapping around) and nd-array with sub-pixel precision.

Parameters:

subject (ndarray) – The to-be-shifted nd-array. Default: all zeros except first element.
shift (Union[int, float, Sequence, ndarray, None]) – The (fractional) shift. Default: all zeros except first element.
axis (Union[int, Sequence, ndarray, None]) – Optional int or sequence of ints. The axis or axes along which elements are shifted. Unlike numpy’s roll function, by default, the left-most axes are used.

Returns:

The shifted nd-array.

macromax.utils.ft.subpixel.roll_ft(subject_ft, shift, axes=None)[source]

Rolls (shifting with wrapping around) and nd-array with sub-pixel precision. The input and output array are Fourier transformed.

Parameters:

subject_ft (ndarray) – The Fourier transform of the to-be-shifted nd-array.
shift (Union[int, float, Sequence, ndarray, None]) – The (fractional) shift.
axes (Union[int, Sequence, ndarray, None]) – Optional int or sequence of ints. The axis or axes along which elements are shifted. Unlike numpy’s roll function, by default, the left-most axes are used.

Returns:

The Fourier transform of the shifted nd-array.

class macromax.utils.ft.subpixel.Registration(shift, factor=1.0, error=0.0, original_ft=None, original=None, registered_ft=None, registered=None)[source]

Bases: object

A class to represent the result of a registration of reference class.

__init__(shift, factor=1.0, error=0.0, original_ft=None, original=None, registered_ft=None, registered=None)[source]

Constructs a registration result.

Parameters:

shift (Union[int, float, Sequence, ndarray]) – translation in pixels of the registered with respect to original reference
factor (complex) – scaling factor between registered and original reference image
error (float) – The root-mean-square difference after registration. See `register` for more details.
original_ft (Union[int, float, Sequence, ndarray, None]) – The Fourier transform of the original image, prior to registration.
original (Union[int, float, Sequence, ndarray, None]) – The original image, prior to registration.
registered_ft (Union[int, float, Sequence, ndarray, None]) – The Fourier transform of the registered image, identical to the original but with a sub-pixel shift.
registered (Union[int, float, Sequence, ndarray, None]) – The registered image, identical to the original but with a sub-pixel shift.

property shift: ndarray: Vector indicating subpixel shift between the original and registration image.

property ndim: int: The number of registration dimensions.

property factor: (Complex) scaling factor indicating the ratio between original and registered image.

property error: float: The rms difference between the registered and the original image including rescaling factor.

property image_ft: ndarray: Fourier transform of the registered and renormalized array.

property image: ndarray: Registered and renormalized array. It is shifted and scaled so that it is as close as possible to the reference. I.e. it minimizes the l2-norm of the difference.

__array__()[source]

Return the registered and renormalized image as an array.

Return type:: ndarray

property original_ft: ndarray: Fourier transform of the original reference array.

property original: ndarray: Original reference array.

class macromax.utils.ft.subpixel.Reference(reference_image=None, precision=None, axes=None, ndim=None, reference_image_ft=None)[source]

Bases: object

Represents a reference ‘image’ (which can be n-dimensional) to be used to register against. Multi-channel arrays should be handled iteratively or by averaging over the channels.

__init__(reference_image=None, precision=None, axes=None, ndim=None, reference_image_ft=None)[source]

Construct a reference ‘image’ object. If neither reference_image or its Fourier transform reference_image_ft are specified, a point-reference is assumed, where the point is in the first element of the nd-array (top-left corner).

Parameters:

reference_image (Union[int, float, Sequence, ndarray, None]) – The reference image nd-array. As an alternative, its unshifted Fourier transform can be specified as reference_image_ft.
precision (Optional[float]) – (optional) The default sub-pixel precision (default: 1/128).
axes (Union[int, float, Sequence, ndarray, None]) – (optional) The axes to operate on. If not specified, all dimensions of the reference image or its Fourier transform are used.
ndim (Optional[int]) – The number of dimensions is usually determined from . If neither is specified, ndim determines the number of dimensions to operate on.
reference_image_ft (Union[int, float, Sequence, ndarray, None]) –

property ndim: int

property shape

register(subject=None, subject_ft=None, precision=None)[source]

Register an nd-image to sub-pixel precision.

Algorithm based on the 2-d implementation of: Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup, “Efficient subpixel image registration algorithms,” Opt. Lett. 33, 156-158 (2008).

Parameters:

subject (Optional[array]) – The subject image as an nd-array. If this is not specified, its Fourier transform should be.
subject_ft (Optional[array]) – (optional) The (non-fftshifted) Fourier transform of the subject image.
precision (Optional[float]) – (optional) The sub-pixel precision in units of pixel (default: 1/128).

Return type:

Registration

Returns:

A Registration object representing the registered image as well as the shift, global phase change, and error.