1. pna2015
  2. pending
In [2]:
# - compatibility with Python 3
from __future__ import print_function  # print('me') instead of print 'me'
from __future__ import division  # 1/2 == 0.5, not 0
In [3]:
# - show figures inside the notebook
%matplotlib inline
In [4]:
# - import common modules
import numpy as np  # the Python array package
import matplotlib.pyplot as plt  # the Python plotting package
In [5]:
# - set gray colormap and nearest neighbor interpolation by default
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'
In [6]:
import nibabel as nib
img = nib.load('ds114_sub009_t2r1.nii')
data = img.get_data()
In [7]:
from scipy.optimize import fmin_powell
from scipy.ndimage import affine_transform
In [8]:
fixed_data = data[..., 1:]
In [11]:
run_mean = fixed_data.mean(axis=-1)
run_mean.shape
Out[11]:
(64, 64, 30)
In [10]:
plt.imshow(run_mean[..., 17])
Out[10]:
<matplotlib.image.AxesImage at 0x106be4bd0>
In [33]:
def cost_function(trans_rots):
    trans = trans_rots[:3]
    M = nib.eulerangles.euler2mat(*trans_rots[3:])
    out_vol = affine_transform(run_mean, M, trans, order=1)
    return mismatch(out_vol, single_vol)
In [34]:
def mismatch(X, Y):
    return -np.corrcoef(X.ravel(), Y.ravel())[0, 1]
In [35]:
def my_callback(params):
   print("Trying parameters " + str(params))
In [36]:
single_vol = fixed_data[..., 0]
best_params = fmin_powell(cost_function, np.zeros(6), callback=my_callback)
best_params

Trying parameters [ 0.          0.00450497  0.          0.          0.          0.        ]
Optimization terminated successfully.
         Current function value: -0.999397
         Iterations: 1
         Function evaluations: 250
Out[36]:
array([ 0.        ,  0.00450497,  0.        ,  0.        ,  0.        ,  0.        ])
In [37]:
single_vol = fixed_data[..., -1]
best_params = fmin_powell(cost_function, np.zeros(6), callback=my_callback)
best_params

Trying parameters [  1.47981169e-02   1.12356064e-02   0.00000000e+00   5.27367624e-05
   0.00000000e+00   0.00000000e+00]
Optimization terminated successfully.
         Current function value: -0.999278
         Iterations: 1
         Function evaluations: 195
Out[37]:
array([  1.47981169e-02,   1.12356064e-02,   0.00000000e+00,
         5.27367624e-05,   0.00000000e+00,   0.00000000e+00])
In [38]:
from nipy.algorithms.registration import HistogramRegistration
In [40]:
to_img = nib.Nifti1Image(run_mean, img.affine, img.header)
from_img = nib.Nifti1Image(data[..., -1], img.affine, img.header)
R = HistogramRegistration(from_img, to_img, similarity='cc', interp='tri')
In [41]:
res = R.optimize('rigid', optimizer='powell')

Initial guess...
translation : [ 0.  0.  0.]
rotation    : [ 0.  0.  0.]

Optimizing using fmin_powell
translation : [-0.04603319 -0.02067415 -0.29859494]
rotation    : [  2.61482986e-04  -1.21159129e-03  -5.02499789e-05]

cc = 0.998871198277

Optimization terminated successfully.
         Current function value: -0.998871
         Iterations: 1
         Function evaluations: 53

/Users/mb312/usr/local/lib/python2.7/site-packages/nipy/algorithms/registration/affine.py:210: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.
  if array == None:
In [42]:
res
Out[42]:
<nipy.algorithms.registration.affine.Rigid at 0x109460ed0>
In [43]:
res.as_affine()
Out[43]:
array([[  9.99999265e-01,   5.00915608e-05,  -1.21159754e-03,
         -4.60331939e-02],
       [ -5.04083712e-05,   9.99999965e-01,  -2.61452478e-04,
         -2.06741531e-02],
       [  1.21158441e-03,   2.61513360e-04,   9.99999232e-01,
         -2.98594938e-01],
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          1.00000000e+00]])
In [44]:
import nipy
In [45]:
nipy.__version__
Out[45]:
'0.4.0.dev'
In [63]:
img = nib.load('ds107_sub012_t1r2.nii')
In [64]:
from nipy.algorithms.registration import SpaceTimeRealign
In [65]:
R = SpaceTimeRealign(img, tr=2.5, slice_times='asc_alt_2', slice_info=2)
In [69]:
R.estimate(refscan=None)
In [70]:
corr_run = R.resample()
In [71]:
from nipy.algorithms.utils.pca import pca_image
In [72]:
nipy_img = nipy.load_image('ds107_sub012_t1r2.nii')
pca1 = pca_image(nipy_img, axis=-1)
In [73]:
pca2 = pca_image(corr_run[0], axis=-1)
In [74]:
fig, axes = plt.subplots(10, 1)
for i, ax in enumerate(axes):
    ax.plot(pca1['basis_vectors'][:, i])
In [75]:
fig, axes = plt.subplots(10, 1)
for i, ax in enumerate(axes):
    ax.plot(pca2['basis_vectors'][:, i])
In [77]:
!mcflirt -in ds107_sub012_t1r2.nii -refvol 1
In [78]:
flirted_img = nipy.load_image('ds107_sub012_t1r2_mcf.nii.gz')
In [80]:
pca3 = pca_image(flirted_img, axis=-1)
In [81]:
fig, axes = plt.subplots(10, 1)
for i, ax in enumerate(axes):
    ax.plot(pca2['basis_vectors'][:, i])
In [82]:
SpaceTimeRealign?
In [86]:
from nipy.algorithms.registration.groupwise_registration import SpaceRealign
In [89]:
realigner = SpaceRealign(img)
In [91]:
realigner.estimate()
In [93]:
corr_run, = realigner.resample()
In [94]:
pca4 = pca_image(corr_run, axis=-1)
In [95]:
fig, axes = plt.subplots(10, 1)
for i, ax in enumerate(axes):
    ax.plot(pca4['basis_vectors'][:, i])
In [96]:
fig, axes = plt.subplots(10, 1)
for i, ax in enumerate(axes):
    ax.plot(pca1['basis_vectors'][:, i])
In [110]:
fig, axes = plt.subplots(2, 5, figsize=(15, 10))
for i, ax in enumerate(axes.ravel()):
    ax.imshow(pca1['basis_projections'][:, :, 17, i])
In [115]:
aff = realigner._transforms[0][-1].as_affine()
In [116]:
import transforms3d

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-116-c800605f6aba> in <module>()
----> 1 import transforms3d

ImportError: No module named transforms3d
In [117]:
rzs, trans = nib.affines.to_matvec(aff)
In [118]:
rzs
Out[118]:
array([[  9.99996110e-01,   9.96587701e-04,   2.60495795e-03],
       [ -9.28735877e-04,   9.99663592e-01,  -2.59198720e-02],
       [ -2.62991304e-03,   2.59173519e-02,   9.99660630e-01]])
In [119]:
nib.eulerangles.mat2euler(rzs)
Out[119]:
(-0.0009965912472167965, 0.002604960893285254, 0.02592286319244751)
In [120]:
trans
Out[120]:
array([ 0.0805639 ,  0.2055159 , -0.38129098])
In [121]:
best_params
Out[121]:
array([  1.47981169e-02,   1.12356064e-02,   0.00000000e+00,
         5.27367624e-05,   0.00000000e+00,   0.00000000e+00])
In [125]:
!mcflirt -in ds107_sub012_t1r2.nii -plots -refvol 0
In [129]:
params = np.loadtxt('ds107_sub012_t1r2_mcf.par')
params[-1]
Out[129]:
array([ 0.0218999 , -0.00631685, -0.00120705,  0.0840234 ,  0.24327   ,
        0.901101  ])
In [130]:
params[0]
Out[130]:
array([ 0., -0.,  0.,  0.,  0.,  0.])
In [ ]: