motion_correction.image_registration package

Subpackages

Submodules

motion_correction.image_registration.chi2_shifts module

Chi^2 shifts

Various tools for calculating shifts based on the chi^2 method

motion_correction.image_registration.chi2_shifts.chi2_shift(im1, im2, err=None, upsample_factor='auto', boundary='wrap', nthreads=1, use_numpy_fft=False, zeromean=False, nfitted=2, verbose=False, return_error=True, return_chi2array=False, max_auto_size=512, max_nsig=1.1)[source]

Find the offsets between image 1 and image 2 using the DFT upsampling method (http://www.mathworks.com/matlabcentral/fileexchange/18401-efficient-subpixel-image-registration-by-cross-correlation/content/html/efficient_subpixel_registration.html) combined with \(\chi^2\) to measure the errors on the fit

Equation 1 gives the \(\chi^2\) value as a function of shift, where Y is the model as a function of shift:

\[\chi^2(dx,dy) = \Sigma_{ij} \frac{(X_{ij}-Y_{ij}(dx,dy))^2}{\sigma_{ij}^2}\]

Equation 2-4:

\begin{align} \mathrm{Term~1:} & f(dx,dy) & = & \Sigma_{ij} \frac{X_{ij}^2}{\sigma_{ij}^2} \\ & f(dx,dy) & = & f(0,0) , \forall dx,dy \\ \mathrm{Term~2:} & g(dx,dy) & = & -2 \Sigma_{ij} \frac{X_{ij}Y_{ij}(dx,dy)}{\sigma_{ij}^2} = -2 \Sigma_{ij} \left(\frac{X_{ij}}{\sigma_{ij}^2}\right) Y_{ij}(dx,dy) \\ \mathrm{Term~3:} & h(dx,dy) & = & \Sigma_{ij} \frac{Y_{ij}(dx,dy)^2}{\sigma_{ij}^2} = \Sigma_{ij} \left(\frac{1}{\sigma_{ij}^2}\right) Y^2_{ij}(dx,dy) \end{align}

The cross-correlation can be computed with fourier transforms, and is defined

\[CC_{m,n}(x,y) = \Sigma_{ij} x^*_{ij} y_{(n+i)(m+j)}\]

which can then be applied to our problem, noting that the cross-correlation has the same form as term 2 and 3 in \(\chi^2\) (term 1 is a constant, with no dependence on the shift)

\begin{align} \mathrm{Term~2:} & CC(X/\sigma^2,Y)[dx,dy] & = & \Sigma_{ij} \left(\frac{X_{ij}}{\sigma_{ij}^2}\right)^* Y_{ij}(dx,dy) \\ \mathrm{Term~3:} & CC(\sigma^{-2},Y^2)[dx,dy] & = & \Sigma_{ij} \left(\frac{1}{\sigma_{ij}^2}\right)^* Y^2_{ij}(dx,dy) \end{align}

Technically, only terms 2 and 3 has any effect on the resulting image, since term 1 is the same for all shifts, and the quantity of interest is \(\Delta \chi^2\) when determining the best-fit shift and error.

The resulting shifts are limited to an accuracy of +/-0.5 pixels in the upsampled image frame. That is not a Gaussian uncertainty but a quantized limit: the best solution will always be +/-0.5 pixels offset because we’re zooming in on an even grid, so the “best” position is required to be a discrete pixel center. If you’re looking at an image with exactly zero shift, it will have exactly +/- 1/usfac/2 systematic error in the resulting solution.

Parameters

im1 : np.ndarray im2 : np.ndarray

The images to register.

errnp.ndarray

Per-pixel error in image 2

boundary‘wrap’,’constant’,’reflect’,’nearest’

Option to pass to map_coordinates for determining what to do with shifts outside of the boundaries.

upsample_factorint or ‘auto’

upsampling factor; governs accuracy of fit (1/usfac is best accuracy) (can be “automatically” determined based on chi^2 error)

return_errorbool

Returns the “fit error” (1-sigma in x and y) based on the delta-chi2 values

return_chi2_arraybool

Returns the x and y shifts and the chi2 as a function of those shifts in addition to other returned parameters. i.e., the last return from this function will be a tuple (x, y, chi2)

zeromeanbool

Subtract the mean from the images before cross-correlating? If no, you may get a 0,0 offset because the DC levels are strongly correlated.

verbosebool

Print error message if upsampling factor is inadequate to measure errors

use_numpy_fftbool

Force use numpy’s fft over fftw? (only matters if you have fftw installed)

nthreadsbool

Number of threads to use for fft (only matters if you have fftw installed)

nfittedint

number of degrees of freedom in the fit (used for chi^2 computations). Should probably always be 2.

max_auto_sizeint

Maximum zoom image size to create when using auto-upsampling

Returns

dx,dyfloat,float

Measures the amount im2 is offset from im1 (i.e., shift im2 by -1 * these #’s to match im1)

errx,erryfloat,float

optional, error in x and y directions

xvals,yvals,chi2n_upsampledndarray,ndarray,ndarray,

x,y positions (in original chi^2 coordinates) of the chi^2 values and their corresponding chi^2 value

Examples

Create a 2d array, shift it in both directions, then use chi2_shift to determine the shift

>>> import image_registration
>>> rr = ((np.indices([100,100]) - np.array([50.,50.])[:,None,None])**2).sum(axis=0)**0.5
>>> image = np.exp(-rr**2/(3.**2*2.)) * 20
>>> shifted = np.roll(np.roll(image,12,0),5,1) + np.random.randn(100,100)
>>> dx,dy,edx,edy = chi2_shift(image, shifted, upsample_factor='auto')
>>> shifted2 = image_registration.fft_tools.shift2d(image,3.665,-4.25) + np.random.randn(100,100)
>>> dx2,dy2,edx2,edy2 = chi2_shift(image, shifted2, upsample_factor='auto')
motion_correction.image_registration.chi2_shifts.chi2_shift_iterzoom(im1, im2, err=None, upsample_factor='auto', boundary='wrap', nthreads=1, use_numpy_fft=False, zeromean=False, verbose=False, return_error=True, return_chi2array=False, zoom_shape=[10, 10], rezoom_shape=[100, 100], rezoom_factor=5, mindiff=1, **kwargs)[source]

Find the offsets between image 1 and image 2 using an iterative DFT upsampling method combined with \(\chi^2\) to measure the errors on the fit

A simpler version of chi2_shift() that only computes the \(\chi^2\) array on the largest scales, then uses a fourier upsampling technique to zoom in.

Parameters

im1 : np.ndarray im2 : np.ndarray

The images to register.

errnp.ndarray

Per-pixel error in image 2

boundary‘wrap’,’constant’,’reflect’,’nearest’

Option to pass to map_coordinates for determining what to do with shifts outside of the boundaries.

upsample_factorint or ‘auto’

upsampling factor; governs accuracy of fit (1/usfac is best accuracy) (can be “automatically” determined based on chi^2 error)

zeromeanbool

Subtract the mean from the images before cross-correlating? If no, you may get a 0,0 offset because the DC levels are strongly correlated.

verbosebool

Print error message if upsampling factor is inadequate to measure errors

use_numpy_fftbool

Force use numpy’s fft over fftw? (only matters if you have fftw installed)

nthreadsbool

Number of threads to use for fft (only matters if you have fftw installed)

nfittedint

number of degrees of freedom in the fit (used for chi^2 computations). Should probably always be 2.

zoom_shape[int,int]

Shape of iterative zoom image

rezoom_shape[int,int]

Shape of the final output chi^2 map to use for determining the errors

rezoom_factorint

Amount to zoom above the last zoom factor. Should be <= rezoom_shape/zoom_shape

Other Parameters

return_errorbool

Returns the “fit error” (1-sigma in x and y) based on the delta-chi2 values

return_chi2_arraybool

Returns the x and y shifts and the chi2 as a function of those shifts in addition to other returned parameters. i.e., the last return from this function will be a tuple (x, y, chi2)

Returns

dx,dyfloat,float

Measures the amount im2 is offset from im1 (i.e., shift im2 by -1 * these #’s to match im1)

errx,erryfloat,float

optional, error in x and y directions

xvals,yvals,chi2n_upsampledndarray,ndarray,ndarray,

x,y positions (in original chi^2 coordinates) of the chi^2 values and their corresponding chi^2 value

Examples

Create a 2d array, shift it in both directions, then use chi2_shift_iterzoom to determine the shift

>>> import image_registration
>>> np.random.seed(42) # so the doctest will pass
>>> image = np.random.randn(50,55)
>>> shifted = np.roll(np.roll(image,12,0),5,1)
>>> dx,dy,edx,edy = chi2_shift_iterzoom(image, shifted, upsample_factor='auto')
>>> shifted2 = image_registration.fft_tools.shift2d(image,3.665,-4.25)
>>> dx2,dy2,edx2,edy2 = chi2_shift_iterzoom(image, shifted2, upsample_factor='auto')
motion_correction.image_registration.chi2_shifts.chi2n_map(im1, im2, err=None, boundary='wrap', nthreads=1, zeromean=False, use_numpy_fft=False, return_all=False, reduced=False)[source]

Parameters

im1 : np.ndarray im2 : np.ndarray

The images to register.

errnp.ndarray

Per-pixel error in image 2

boundary‘wrap’,’constant’,’reflect’,’nearest’

Option to pass to map_coordinates for determining what to do with shifts outside of the boundaries.

zeromeanbool

Subtract the mean from the images before cross-correlating? If no, you may get a 0,0 offset because the DC levels are strongly correlated.

nthreadsbool

Number of threads to use for fft (only matters if you have fftw installed)

reducedbool

Return the reduced \(\chi^2\) array, or unreduced? (assumes 2 degrees of freedom for the fit)

Returns

chi2nnp.ndarray

the \(\chi^2\) array

term1float

Scalar, term 1 in the \(\chi^2\) equation

term2np.ndarray

Term 2 in the equation, -2 * cross-correlation(x/sigma^2,y)

term3np.ndarray | float

If error is an array, returns an array, otherwise is a scalar float corresponding to term 3 in the equation

motion_correction.image_registration.conftest module

Configure Test Suite.

This file is used to configure the behavior of pytest when using the Astropy test infrastructure. It needs to live inside the package in order for it to get picked up when running the tests inside an interpreter using image_registration.test().

motion_correction.image_registration.conftest.pytest_configure(config)[source]

Configure Pytest with Astropy.

Parameters

config : pytest configuration

motion_correction.image_registration.cross_correlation_shifts module

motion_correction.image_registration.cross_correlation_shifts.cross_correlation_shifts(image1, image2, errim1=None, errim2=None, maxoff=None, verbose=False, gaussfit=False, return_error=False, zeromean=True, **kwargs)[source]

Use cross-correlation and a 2nd order taylor expansion to measure the offset between two images

Given two images, calculate the amount image2 is offset from image1 to sub-pixel accuracy using 2nd order taylor expansion.

Parameters

image1: np.ndarray

The reference image

image2: np.ndarray

The offset image. Must have the same shape as image1

errim1: np.ndarray [optional]

The pixel-by-pixel error on the reference image

errim2: np.ndarray [optional]

The pixel-by-pixel error on the offset image.

maxoff: int

Maximum allowed offset (in pixels). Useful for low s/n images that you know are reasonably well-aligned, but might find incorrect offsets due to edge noise

zeromeanbool

Subtract the mean from each image before performing cross-correlation?

verbose: bool

Print out extra messages?

gaussfitbool

Use a Gaussian fitter to fit the peak of the cross-correlation?

return_error: bool

Return an estimate of the error on the shifts. WARNING: I still don’t understand how to make these agree with simulations. The analytic estimate comes from http://adsabs.harvard.edu/abs/2003MNRAS.342.1291Z At high signal-to-noise, the analytic version overestimates the error by a factor of about 1.8, while the gaussian version overestimates error by about 1.15. At low s/n, they both UNDERestimate the error. The transition zone occurs at a total S/N ~ 1000 (i.e., the total signal in the map divided by the standard deviation of the map - it depends on how many pixels have signal)

**kwargs are passed to correlate2d, which in turn passes them to convolve. The available options include image padding for speed and ignoring NaNs.

References

From http://solarmuri.ssl.berkeley.edu/~welsch/public/software/cross_cor_taylor.pro

Examples

>>> import numpy as np
>>> im1 = np.zeros([10,10])
>>> im2 = np.zeros([10,10])
>>> im1[4,3] = 1
>>> im2[5,5] = 1
>>> import image_registration
>>> yoff,xoff = image_registration.cross_correlation_shifts(im1,im2)
>>> im1_aligned_to_im2 = np.roll(np.roll(im1,int(yoff),1),int(xoff),0)
>>> assert (im1_aligned_to_im2-im2).sum() == 0
motion_correction.image_registration.cross_correlation_shifts.cross_correlation_shifts_FITS(fitsfile1, fitsfile2, return_cropped_images=False, sigma_cut=False, verbose=False, register_method=<function cross_correlation_shifts>, **kwargs)[source]

Determine the shift between two FITS images using the cross-correlation technique. Requires reproject

Parameters

fitsfile1: str

Reference fits file name

fitsfile2: str

Offset fits file name

return_cropped_images: bool

Returns the images used for the analysis in addition to the measured offsets

sigma_cut: bool or int

Perform a sigma-cut before cross-correlating the images to minimize noise correlation?

motion_correction.image_registration.iterative_zoom module

motion_correction.image_registration.iterative_zoom.centers_to_edges(arr)[source]
motion_correction.image_registration.iterative_zoom.iterative_zoom(image, mindiff=1.0, zoomshape=[10, 10], return_zoomed=False, zoomstep=2, verbose=False, minmax=<function min>, ploteach=False, return_center=True)[source]

Iteratively zoom in on the minimum position in an image until the delta-peak value is below mindiff

Parameters

imagenp.ndarray

Two-dimensional image with a minimum to zoom in on (or maximum, if specified using minmax)

mindifffloat

Minimum difference that must be present in image before zooming is done

zoomshape[int,int]

Shape of the “mini” image to create. Smaller is faster, but a bit less accurate. [10,10] seems to work well in preliminary tests (though unit tests have not been written)

return_zoomedbool

Return the zoomed image in addition to the measured offset?

zoomstepint

Amount to increase the zoom factor by on each iteration. Probably best to stick with small integers (2-5ish).

verbosebool

Print out information about zoom factor, offset at each iteration

minmaxnp.min or np.max

Can zoom in on the minimum or maximum of the image

ploteachbool

Primarily a debug tool, and to be used with extreme caution! Will open a new figure at each iteration showing the next zoom level.

return_centerbool

Return the center position in original image coordinates? If False, will retern the offset from center instead (but beware the conventions associated with the concept of ‘center’ for even images).

Returns

The y,x offsets (following numpy convention) of the center position of the original image. If return_zoomed, returns (zoomed_image, zoom_factor, offsets) because you can’t interpret the zoomed image without the zoom factor.

motion_correction.image_registration.iterative_zoom.iterative_zoom_1d(data, mindiff=1.0, zoomshape=(10, ), return_zoomed=False, zoomstep=2, verbose=False, minmax=<function min>, return_center=True)[source]

Iteratively zoom in on the minimum position in a spectrum or timestream until the delta-peak value is below mindiff

Parameters

datanp.ndarray

One-dimensional array with a minimum (or maximum, as specified by minmax) to zoom in on

mindifffloat

Minimum difference that must be present in image before zooming is done

zoomshapeint

Shape of the “mini” image to create. Smaller is faster, but a bit less accurate. 10 seems to work well in preliminary tests (though unit tests have not been written)

return_zoomedbool

Return the zoomed image in addition to the measured offset?

zoomstepint

Amount to increase the zoom factor by on each iteration. Probably best to stick with small integers (2-5ish).

verbosebool

Print out information about zoom factor, offset at each iteration

minmaxnp.min or np.max

Can zoom in on the minimum or maximum of the image

return_centerbool

Return the center position in original image coordinates? If False, will retern the offset from center instead (but beware the conventions associated with the concept of ‘center’ for even images).

Returns

The x offsets of the center position of the original spectrum. If return_zoomed, returns (zoomed_image, zoom_factor, offsets) because you can’t interpret the zoomed spectrum without the zoom factor.

motion_correction.image_registration.register_images module

motion_correction.image_registration.register_images.register_images(im1, im2, usfac=1, return_registered=False, return_error=False, zeromean=True, DEBUG=False, maxoff=None, nthreads=1, use_numpy_fft=False)[source]

Sub-pixel image registration (see dftregistration for lots of details)

Parameters

im1 : np.ndarray im2 : np.ndarray

The images to register.

usfacint

upsampling factor; governs accuracy of fit (1/usfac is best accuracy)

return_registeredbool

Return the registered image as the last parameter

return_errorbool

Does nothing at the moment, but in principle should return the “fit error” (it does nothing because I don’t know how to compute the “fit error”)

zeromeanbool

Subtract the mean from the images before cross-correlating? If no, you may get a 0,0 offset because the DC levels are strongly correlated.

maxoffint

Maximum allowed offset to measure (setting this helps avoid spurious peaks)

DEBUGbool

Test code used during development. Should DEFINITELY be removed.

Returns

dx,dyfloat,float

REVERSE of dftregistration order (also, signs flipped) for consistency with other routines. Measures the amount im2 is offset from im1 (i.e., shift im2 by these #’s to match im1)

motion_correction.image_registration.version module

Module contents

This is an Astropy affiliated package.