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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
|
.. _coordinate_map:
#############################
Basics of the Coordinate Map
#############################
When you load an image it will have an associated Coordinate Map
**Coordinate Map**
The Coordinate Map contains information defining the input (domain) and
output (range) Coordinate Systems of the image, and the mapping between the
two Coordinate systems.
The *input* or *domain* in an image are voxel coordinates in the image array.
The *output* or *range* are the millimetre coordinates in some space, that
correspond to the input (voxel) coordinates.
>>> import nipy
Get a filename for an example file:
>>> from nipy.testing import anatfile
Get the coordinate map for the image:
>>> anat_img = nipy.load_image(anatfile)
>>> coordmap = anat_img.coordmap
For more on Coordinate Systems and their properties
:mod:`nipy.core.reference.coordinate_system`
You can inspect a coordinate map::
>>> coordmap.function_domain.coord_names
>>> ('i', 'j', 'k')
>>> coordmap.function_range.coord_names
('aligned-x=L->R', 'aligned-y=P->A', 'aligned-z=I->S')
>>> coordmap.function_domain.name
'voxels'
>>> coordmap.function_range.name
'aligned'
A Coordinate Map has a mapping from the *input* Coordinate System to the
*output* Coordinate System
Here we can see we have a voxel to millimeter mapping from the voxel
space (i,j,k) to the millimeter space (x,y,z)
We can also get the name of the respective Coordinate Systems that our
Coordinate Map maps between.
A Coordinate Map is two Coordinate Systems with a mapping between
them. Formally the mapping is a function that takes points from the
input Coordinate System and returns points from the output Coordinate
System. This is the same as saying that the mapping takes points in the mapping
function *domain* and transforms them to points in the mapping function *range*.
Often this is simple as applying an Affine transform. In that case the
Coordinate System may well have an affine property which returns the
affine matrix corresponding to the transform.
>>> coordmap.affine
array([[ -2., 0., 0., 32.],
[ 0., 2., 0., -40.],
[ 0., 0., 2., -16.],
[ 0., 0., 0., 1.]])
If you call the Coordinate Map you will apply the mapping function
between the two Coordinate Systems. In this case from (i,j,k) to (x,y,z):
>>> coordmap([1,2,3])
array([ 30., -36., -10.])
It can also be used to get the inverse mapping, or in this example from (x,y,z)
back to (i,j,k):
>>> coordmap.inverse()([30.,-36.,-10.])
array([1., 2., 3.])
We can see how this works if we just apply the affine
ourselves using dot product.
.. Note::
Notice the affine is using homogeneous coordinates so we need to add a 1 to
our input. (And note how a direct call to the coordinate map does this work
for you)
>>> coordmap.affine
array([[ -2., 0., 0., 32.],
[ 0., 2., 0., -40.],
[ 0., 0., 2., -16.],
[ 0., 0., 0., 1.]])
>>> import numpy as np
>>> np.dot(coordmap.affine, np.transpose([1,2,3,1]))
array([ 30., -36., -10., 1.])
.. Note::
The answer is the same as above (except for the added 1)
.. _normalize-coordmap:
***************************************************
Use of the Coordinate Map for spatial normalization
***************************************************
The Coordinate Map can be used to describe the transformations needed to perform
spatial normalization. Suppose we have an anatomical Image from one subject
*subject_img* and we want to create an Image in a standard space like Tailarach
space. An affine registration algorithm will produce a 4-by-4 matrix
representing the affine transformation, *T*, that takes a point in the subject's
coordinates *subject_world* to a point in Tailarach space *tailarach_world*. The
subject's Image has its own Coordinate Map, *subject_cmap* and there is a
Coordinate Map for Tailarach space which we will call *tailarach_cmap*.
Having found the transformation matrix *T*, the next step in spatial
normalization is usually to resample the array of *subject_img* so that it has
the same shape as some atlas *atlas_img*. Note that because it is an atlas
Image, *tailarach_camp=atlas_img.coordmap*.
A resampling algorithm uses an interpolator which needs to know
which voxel of *subject_img* corresponds to which voxel of *atlas_img*.
This is therefore a function from *atlas_voxel* to *subject_voxel*.
This function, paired with the information that it is a map from atlas-voxel to
subject-voxel is another example of a Coordinate Map. The code to do this might
look something like the following:
>>> from nipy.testing import anatfile, funcfile
>>> from nipy.algorithms.registration import HistogramRegistration
>>> from nipy.algorithms.kernel_smooth import LinearFilter
We'll make a smoothed version of the anatomical example image, and pretend it's
the template
>>> smoother = LinearFilter(anat_img.coordmap, anat_img.shape)
>>> atlas_im = smoother.smooth(anat_img)
>>> subject_im = anat_img
We do an affine registration between the two.
>>> reggie = HistogramRegistration(subject_im, atlas_im)
>>> aff = reggie.optimize('affine').as_affine() #doctest: +ELLIPSIS
Initial guess...
...
Now we make a coordmap with this transformation
>>> from nipy.core.api import AffineTransform
>>> subject_cmap = subject_im.coordmap
>>> talairach_cmap = atlas_im.coordmap
>>> subject_world_to_talairach_world = AffineTransform(
... subject_cmap.function_range,
... talairach_cmap.function_range,
... aff)
...
We resample the 'subject' image to the 'atlas image
>>> from nipy.algorithms.resample import resample
>>> normalized_subject_im = resample(subject_im, talairach_cmap,
... subject_world_to_talairach_world,
... atlas_im.shape)
>>> normalized_subject_im.shape == atlas_im.shape
True
>>> normalized_subject_im.coordmap == atlas_im.coordmap
True
>>> # Normalized image now has atlas affine.
>>> assert np.all(normalized_subject_im.affine == atlas_im.affine)
***********************
Mathematical definition
***********************
For a more formal mathematical description of the coordinate map, see
:ref:`math-coordmap`.
|
