Thursday, September 26, 2013

Testing image registration to correct for brain motion

Today I am going to test out a method for correcting for image motion, based on "Efficient subpixel image registration algorithms," Opt. Lett. 33, 156-158 (2008). There is a MATLAB implementation that seems to be widely used. This has been translated into python here:
http://image-registration.readthedocs.org/en/latest/
github page here:
https://github.com/keflavich/image_registration
from this page I cloned the repository (actually, I just downloaded the zipped directory into my src directory, then installed it using the following commands).

cd image_registration/
sudo python setup.py install

Attempting to run it informed me that I was missing some dependencies, which I downloaded using the following commands (additional information on these is located here http://pytest.org/latest/contents.html and here http://www.stsci.edu/institute/software_hardware/pyfits and here http://stsdas.stsci.edu/astrolib/pywcs/).

pip install -U pytest
sudo pip install pyfits
sudo pip install pywcs

Now, I am up and running! I tested image_registration.chi2_shift by modifying one of the examples. Since I think it will be helpful to "stabilize" and image, I want to not only find out the amount one image is shifted relative to another, but be able to actually shift it back into register. This test snippet allowed me to test this:

import image_registration
import numpy as np
import matplotlib.pyplot as plt
plt.ion()

rr = ((np.indices([100,100]) - np.array([50.,50.])[:,None,None])**2).sum(axis=0)**0.5
exampleImage = np.exp(-rr**2/(3.**2*2.)) * 20
fig = plt.figure()
ax1 = fig.add_subplot(221)
ax1.imshow(exampleImage)
ax1.set_title('Original image')

ax2 = fig.add_subplot(222)
shiftedImage = np.roll(np.roll(exampleImage,12,0),5,1) + np.random.randn(100,100)
ax2.imshow(shiftedImage)
ax2.set_title('Noisy shifted image')

dx,dy,edx,edy = image_registration.chi2_shift(exampleImage, shiftedImage, upsample_factor='auto')
shiftedBackImage = image_registration.fft_tools.shift2d(shiftedImage,-dx,-dy)
ax3 = fig.add_subplot(224)
ax3.imshow(shiftedBackImage)
ax3.set_title('Noisy image shifted to match original')

Looks good!