On using a custom Transform
object in nipy registration
%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 nipy
import nipy.core.api as nca
/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/nose/plugins/manager.py:418: UserWarning: Module nipy was already imported from /Users/mb312/usr/local/lib/python2.7/site-packages/nipy/__init__.pyc, but /Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages is being added to sys.path import pkg_resources
fname = 'ds107_sub001_run01.nii.gz'
img4f = nipy.load_image(fname)
iterable_img = nca.iter_axis(img4f, 't')
vols = nca.ImageList(iterable_img)
vol0 = vols[0]
vol1 = vols[1]
vol0.shape
(64, 64, 35)
We can use histogram registration to register these
import nipy.algorithms.registration as reg
We are going to calculate the transform from the world space of the from
image, to the world space of the to
image.
reggie = reg.HistogramRegistration(from_img=vol0, to_img=vol1, interp='rand')
And run the registration thus:
rigid_estimated = reggie.optimize('rigid')
Initial guess... translation : [ 0. 0. 0.] rotation : [ 0. 0. 0.] Optimizing using fmin_powell translation : [ 4.69804057e-05 2.50401518e-05 5.49921424e-05] rotation : [ 3.04449687e-07 8.62186857e-09 -6.00736588e-07] crl1 = 0.936693315873 Optimization terminated successfully. Current function value: -0.936694 Iterations: 1 Function evaluations: 106
print(rigid_estimated)
translation : [ 4.69804057e-05 2.50401518e-05 5.49921424e-05] rotation : [ 3.04449687e-07 8.62186857e-09 -6.00736588e-07]
The input string 'rigid' is a shortcut for using a Rigid transform object:
from nipy.algorithms.registration import affine
transform = affine.Rigid()
transform.param
array([ 0., 0., 0., 0., 0., 0.])
rigid_reestimated = reggie.optimize(transform)
Initial guess... translation : [ 0. 0. 0.] rotation : [ 0. 0. 0.] Optimizing using fmin_powell translation : [ -1.17341540e-04 2.84777093e-04 -2.47568627e-05] rotation : [ -2.91789583e-06 0.00000000e+00 1.98685446e-07] crl1 = 0.936665573091 Optimization terminated successfully. Current function value: -0.936679 Iterations: 1 Function evaluations: 131
rigid_reestimated.param
array([ -1.17341540e-04, 2.84777093e-04, -2.47568627e-05, -2.91789583e-04, 0.00000000e+00, 1.98685446e-05])
The transform is random to some degree, because we used 'rand' as our interpolation method. We can show the string 'rigid' is the same as using the affine.Rigid()
object by seeding the random number generator.
import numpy as np
np.random.seed(42)
print(reggie.optimize('rigid').param)
Initial guess... translation : [ 0. 0. 0.] rotation : [ 0. 0. 0.] Optimizing using fmin_powell translation : [ 3.79059250e-05 1.65194246e-05 4.86191008e-05] rotation : [ -1.89194185e-07 3.00152800e-06 -2.96596799e-06] crl1 = 0.936678251118 Optimization terminated successfully. Current function value: -0.936695 Iterations: 1 Function evaluations: 98 [ 3.79059250e-05 1.65194246e-05 4.86191008e-05 -1.89194185e-05 3.00152800e-04 -2.96596799e-04]
np.random.seed(42)
print(reggie.optimize(affine.Rigid()).param)
Initial guess... translation : [ 0. 0. 0.] rotation : [ 0. 0. 0.] Optimizing using fmin_powell translation : [ 3.79059250e-05 1.65194246e-05 4.86191008e-05] rotation : [ -1.89194185e-07 3.00152800e-06 -2.96596799e-06] crl1 = 0.936678251118 Optimization terminated successfully. Current function value: -0.936695 Iterations: 1 Function evaluations: 98 [ 3.79059250e-05 1.65194246e-05 4.86191008e-05 -1.89194185e-05 3.00152800e-04 -2.96596799e-04]
print(reggie.optimize('rigid'))
Initial guess... translation : [ 0. 0. 0.] rotation : [ 0. 0. 0.] Optimizing using fmin_powell translation : [ -7.74728790e-06 1.07013893e-05 1.13706640e-05] rotation : [ 2.13755088e-06 -1.18215733e-06 1.63740027e-06] crl1 = 0.936666705588 Optimization terminated successfully. Current function value: -0.936682 Iterations: 1 Function evaluations: 99 translation : [ -7.74728790e-06 1.07013893e-05 1.13706640e-05] rotation : [ 2.13755088e-06 -1.18215733e-06 1.63740027e-06]
What kind of thing must a new "Transform" object do?
For the optimization, the transform must be able to provide two things:
param
property / attributeapply
methodThe param
property returns and accepts values that parameterize the transform, and that the optimizer can change to find the next step of the optimization.
The apply
method accepts a set of coordinates xyz
of shape (usually) (M, N, P, 3) where (M, N, P) will be the shape of the image being optimized, and the last dimension indexes the X, Y, Z coordinate values.
Lastly, the transform should be able to compose itself with other transforms. So it needs a compose
method.
The compose
method of our transform mytransform
needs to take another transform other
as input, and return a new transform, new
. new
needs only to provide an apply
method, that will apply the transform other
and then the transform mytransform
.
Here's a very simple example of a transform that does only translations (in mm).
from nipy.algorithms.registration.transform import Transform
class PointShifter(object):
def __init__(self, param=None):
if param is None:
param = np.zeros((3,))
self.param = np.array(param, dtype=np.float64)
def apply(self, xyz):
""" The params are scalars to add to X, Y, Z """
return xyz + self.param # This will broadcast
def __str__(self): # To see values during optimization
return "Pointshifter: " + str(self.param)
def compose(self, other):
""" Make a new transform that is the composition of this transform on another
"""
apply_func = lambda pts : self.apply(other.apply(pts))
return Transform(apply_func)
adder = PointShifter([0, 0, 0])
adder.param
array([ 0., 0., 0.])
translations = reggie.optimize(adder)
Initial guess... Pointshifter: [ 0. 0. 0.] Optimizing using fmin_powell Pointshifter: [ 3.57023227e-05 4.59146450e-05 -3.93127536e-06] crl1 = 0.936691223873 Optimization terminated successfully. Current function value: -0.936697 Iterations: 1 Function evaluations: 51
translations.param
array([ 3.57023227e-05, 4.59146450e-05, -3.93127536e-06])
We can optimize two sets of transforms at the same time
class Combiner(object):
def __init__(self, transforms):
self.transforms = transforms
@property
def param(self):
return np.hstack([t.param for t in self.transforms])
@param.setter
def param(self, param):
offset = 0
for t in self.transforms:
n_param = len(t.param)
t.param = param[offset:offset+n_param]
offset += n_param
def apply(self, xyz):
if len(self.transforms) == 0:
return xyz
composed = self.transforms[0]
for t in self.transforms[1:]:
composed = t.compose(composed)
return composed.apply(xyz)
def __str__(self):
return '\n'.join(str(t) for t in self.transforms)
def compose(self, other):
""" Make a new transform that is the composition of this transform on another
"""
apply_func = lambda pts : self.apply(other.apply(pts))
return Transform(apply_func)
both = Combiner((PointShifter(), affine.Rigid()))
best_both = reggie.optimize(both)
Initial guess... Pointshifter: [ 0. 0. 0.] translation : [ 0. 0. 0.] rotation : [ 0. 0. 0.] Optimizing using fmin_powell Pointshifter: [ -5.37797973e-05 -3.05024211e-05 -2.20785971e-05] translation : [ 1.91941445e-06 3.17421997e-05 2.08340892e-04] rotation : [ 3.60652070e-09 -4.54284910e-08 1.22395543e-05] crl1 = 0.93666397136 Optimization terminated successfully. Current function value: -0.936686 Iterations: 1 Function evaluations: 166
best_both.param
array([ -5.37797973e-05, -3.05024211e-05, -2.20785971e-05, 1.91941445e-06, 3.17421997e-05, 2.08340892e-04, 3.60652070e-07, -4.54284910e-06, 1.22395543e-03])
Let's translate an image up by 3 pixels and left by 10 pixels
vol1_data = vol1.get_data()
plt.imshow(vol1_data[:, :, 15])
<matplotlib.image.AxesImage at 0x114370bd0>
new_vol = np.zeros_like(vol1_data)
new_vol[:-10, 3:, :] = vol1_data[10:, :-3, :]
plt.imshow(new_vol[:, :, 15])
<matplotlib.image.AxesImage at 0x1150a00d0>
How many mm is that?
vol1.metadata['header'].get_zooms()
(3.0, 3.0, 3.0, 3.0)
np.multiply((-10, 3, 0), vol1.metadata['header'].get_zooms()[:3])
array([-30., 9., 0.])
30mm left, 9mm forwards
shifted_img = nca.Image(new_vol, vol1.coordmap)
fig, axes = plt.subplots(1, 2)
axes[0].imshow(vol0.get_data()[..., 15])
axes[1].imshow(shifted_img.get_data()[...,15])
<matplotlib.image.AxesImage at 0x1156c6650>
adder = PointShifter([0, 0, 0])
reggie = reg.HistogramRegistration(vol0, shifted_img, interp='rand')
added = reggie.optimize(adder)
Initial guess... Pointshifter: [ 0. 0. 0.] Optimizing using fmin_powell Pointshifter: [-50.00718377 7.66912205 0.381966 ] crl1 = 0.385950590763 Pointshifter: [-21.40660419 8.66912205 -0.23606797] crl1 = 0.591779265694 Pointshifter: [-29.87874028 9.05108805 0.14589803] crl1 = 0.924590553756 Pointshifter: [ -3.00051911e+01 8.97564693e+00 8.21099270e-03] crl1 = 0.941022108705 Pointshifter: [ -3.00004488e+01 9.01386211e+00 -1.35348482e-03] crl1 = 0.942871747399 Optimization terminated successfully. Current function value: -0.942753 Iterations: 5 Function evaluations: 117
added.param
array([ -3.00004488e+01, 9.01386211e+00, -1.35348482e-03])
fixed = reg.resample(shifted_img, added, vol0)
fig, axes = plt.subplots(1, 3)
axes[0].imshow(vol0.get_data()[..., 15])
axes[1].imshow(shifted_img.get_data()[...,15])
axes[2].imshow(fixed.get_data()[...,15])
<matplotlib.image.AxesImage at 0x106b36c90>
Of course here all the parameters are in mm space - after applying the voxel to world transforms of the images. We might want to apply transforms in the voxel space - like for the eddy current correction.
We can do that by optimizing the corresponding parameters in mm space.