1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
|
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
import numpy as np
from nipy.neurospin.registration import *
from nipy.neurospin.registration.grid_transform import *
from nipy.neurospin.image import *
from nipy.utils import example_data
from nipy.io.imageformats import load, save
### DEBUG
from numpy.testing import *
from os.path import join
import time
print('Scanning data directory...')
# Input images are provided with the nipy-data package
source = 'ammon'
target = 'anubis'
source_file = example_data.get_filename('neurospin','sulcal2000','nobias_'+source+'.nii.gz')
target_file = example_data.get_filename('neurospin','sulcal2000','nobias_'+target+'.nii.gz')
# Optional arguments
similarity = 'cr'
interp = 'pv'
optimizer = 'powell'
# Make registration instance
I = Image(load(source_file))
J = Image(load(target_file))
R = IconicRegistration(I, J)
R.set_source_fov(fixed_npoints=64**3)
# Make Gaussian spline transform instance
spacing = 16
slices = [slice(0,s.stop,s.step*spacing) for s in R._slices]
cp = np.mgrid[slices]
cp = np.rollaxis(cp, 0, 4)
# Start with an affine registration
A0 = Affine()
A = R.optimize(A0)
###A = Affine()
# Save affinely transformed target
Jt = J.transform(A, reference=I)
save(asNifti1Image(Jt), 'affine_anubis_to_ammon.nii')
# Then add control points...
T0 = SplineTransform(I, cp, sigma=20., grid_coords=True, affine=A)
"""
# Test 1
s = R.eval(T0)
sa = R.eval(T0.affine)
assert_almost_equal(s, sa)
# Test 2
T = SplineTransform(I, cp, sigma=5., grid_coords=True, affine=A)
T.param += 1.
s0 = R.eval(T0)
s = R.eval(T)
print(s-s0)
"""
# Optimize spline transform
T = R.optimize(T0, method='steepest')
###T = R.optimize(T0)
###T = T0
###T.param = np.load('spline_param.npy')
# Resample target image
t = np.asarray(T)
Jt = J.transform(t, reference=I)
save(asNifti1Image(Jt), 'deform_anubis_to_ammon.nii')
# Test 3
"""
ts = t[R._slices+[slice(0,3)]]
tts = T[R._slices]()
"""
|
