A small nipy registration example
%pylab inline
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
import numpy as np
import matplotlib.pyplot as plt
import nibabel as nib
We're going to use the example dataset. That should be here:
import os
DATA_PATH = os.path.expanduser(os.path.join('~', 'data', 'ds105'))
os.listdir(DATA_PATH)
['study_key.txt', 'release_history.txt', 'task_key.txt', 'models', 'references.txt', 'sub005', 'license.txt', 'sub004', 'sub003', 'scan_key.txt', 'sub001', 'sub006', 'README', 'sub002']
SUB01 = os.path.join(DATA_PATH, 'sub001')
anat_fname = os.path.join(SUB01, 'anatomy', 'highres001.nii.gz')
anat_img = nib.load(anat_fname)
anat_img.shape
(124, 256, 256)
func_fname = os.path.join(SUB01, 'BOLD', 'task001_run001', 'bold.nii.gz')
func_img = nib.load(func_fname)
func_img.shape
(40, 64, 64, 121)
We make a mean functional image to register against:
func_data = func_img.get_data()
func_mean_data = np.mean(func_data, axis=3) # Mean over time
func_mean_img = nib.Nifti1Image(func_mean_data, func_img.get_affine(), func_img.get_header())
Now for the registration:
import nipy.algorithms.registration as nar
Investigate nar.HistogramRegistration
register_obj = nar.HistogramRegistration(from_img = anat_img, to_img = func_mean_img, similarity='nmi')
transform = register_obj.optimize('rigid', optimizer='powell')
Initial guess... translation : [ 0. 0. 0.] rotation : [ 0. 0. 0.] Optimizing using fmin_powell translation : [ 1. 1. 1.] rotation : [ 0.00381966 -0.00618034 -0.0005701 ] nmi = 0.12714613287 translation : [ 0.38196603 1.381966 1.381966 ] rotation : [ -2.36067975e-03 3.81966025e-03 -9.41613423e-06] nmi = 0.12953704528 translation : [ 0.56019456 1.46248928 1.36699688] rotation : [ 0.00145898 -0.00236068 -0.00019591] nmi = 0.129906538601 Optimization terminated successfully. Current function value: -0.129907 Iterations: 3 Function evaluations: 114
print(transform)
translation : [ 0.56019456 1.46248928 1.36699688] rotation : [ 0.00145898 -0.00236068 -0.00019591]
Investigate nar.resample
resampled = nar.resample(anat_img, transform, func_mean_img)
resampled_anat_data = resampled.get_data()
resampled_anat_data.shape
(40, 64, 64)
fig, axes = plt.subplots(1, 2)
fig.set_size_inches(10.5,18.5)
axes[0].imshow(resampled_anat_data[:, :, 38], cmap="gray")
axes[1].imshow(func_mean_data[:, :, 38], cmap="gray")
<matplotlib.image.AxesImage at 0x599c4d0>