1"""
2=========================
3Affine Registration in 3D
4=========================
5
6This example explains how to compute an affine transformation to register two
73D volumes by maximization of their Mutual Information [Mattes03]_. The
8optimization strategy is similar to that implemented in ANTS [Avants11]_.
9
10We will do this twice. The first part of this tutorial will walk through the
11details of the process with the object-oriented interface implemented in
12the ``dipy.align`` module. The second part will use a simplified functional
13interface.
14"""
15
16from os.path import join as pjoin
17import numpy as np
18from dipy.viz import regtools
19from dipy.data import fetch_stanford_hardi
20from dipy.data.fetcher import fetch_syn_data
21from dipy.io.image import load_nifti
22from dipy.align.imaffine import (transform_centers_of_mass,
23                                 AffineMap,
24                                 MutualInformationMetric,
25                                 AffineRegistration)
26from dipy.align.transforms import (TranslationTransform3D,
27                                   RigidTransform3D,
28                                   AffineTransform3D)
29
30"""
31Let's fetch two b0 volumes, the static image will be the b0 from the Stanford
32HARDI dataset
33"""
34
35files, folder = fetch_stanford_hardi()
36static_data, static_affine, static_img = load_nifti(
37                                            pjoin(folder, 'HARDI150.nii.gz'),
38                                            return_img=True)
39static = np.squeeze(static_data)[..., 0]
40static_grid2world = static_affine
41
42"""
43Now the moving image
44"""
45
46files, folder2 = fetch_syn_data()
47moving_data, moving_affine, moving_img = load_nifti(
48                                            pjoin(folder2, 'b0.nii.gz'),
49                                            return_img=True)
50moving = moving_data
51moving_grid2world = moving_affine
52
53"""
54We can see that the images are far from aligned by drawing one on top of
55the other. The images don't even have the same number of voxels, so in order
56to draw one on top of the other we need to resample the moving image on a grid
57of the same dimensions as the static image, we can do this by "transforming"
58the moving image using an identity transform
59"""
60
61identity = np.eye(4)
62affine_map = AffineMap(identity,
63                       static.shape, static_grid2world,
64                       moving.shape, moving_grid2world)
65resampled = affine_map.transform(moving)
66regtools.overlay_slices(static, resampled, None, 0,
67                        "Static", "Moving", "resampled_0.png")
68regtools.overlay_slices(static, resampled, None, 1,
69                        "Static", "Moving", "resampled_1.png")
70regtools.overlay_slices(static, resampled, None, 2,
71                        "Static", "Moving", "resampled_2.png")
72
73"""
74.. figure:: resampled_0.png
75   :align: center
76.. figure:: resampled_1.png
77   :align: center
78.. figure:: resampled_2.png
79   :align: center
80
81   Input images before alignment.
82"""
83
84"""
85We can obtain a very rough (and fast) registration by just aligning the centers
86of mass of the two images
87"""
88
89c_of_mass = transform_centers_of_mass(static, static_grid2world,
90                                      moving, moving_grid2world)
91
92"""
93We can now transform the moving image and draw it on top of the static image,
94registration is not likely to be good, but at least they will occupy roughly
95the same space
96"""
97
98transformed = c_of_mass.transform(moving)
99regtools.overlay_slices(static, transformed, None, 0,
100                        "Static", "Transformed", "transformed_com_0.png")
101regtools.overlay_slices(static, transformed, None, 1,
102                        "Static", "Transformed", "transformed_com_1.png")
103regtools.overlay_slices(static, transformed, None, 2,
104                        "Static", "Transformed", "transformed_com_2.png")
105
106"""
107.. figure:: transformed_com_0.png
108   :align: center
109.. figure:: transformed_com_1.png
110   :align: center
111.. figure:: transformed_com_2.png
112   :align: center
113
114   Registration result by aligning the centers of mass of the images.
115"""
116
117"""
118This was just a translation of the moving image towards the static image, now
119we will refine it by looking for an affine transform. We first create the
120similarity metric (Mutual Information) to be used. We need to specify the
121number of bins to be used to discretize the joint and marginal probability
122distribution functions (PDF), a typical value is 32. We also need to specify
123the percentage (an integer in (0, 100]) of voxels to be used for computing the
124PDFs, the most accurate registration will be obtained by using all voxels, but
125it is also the most time-consuming choice. We specify full sampling by passing
126None instead of an integer
127"""
128
129nbins = 32
130sampling_prop = None
131metric = MutualInformationMetric(nbins, sampling_prop)
132
133"""
134To avoid getting stuck at local optima, and to accelerate convergence, we use a
135multi-resolution strategy (similar to ANTS [Avants11]_) by building a Gaussian
136Pyramid. To have as much flexibility as possible, the user can specify how this
137Gaussian Pyramid is built. First of all, we need to specify how many
138resolutions we want to use. This is indirectly specified by just providing a
139list of the number of iterations we want to perform at each resolution. Here we
140will just specify 3 resolutions and a large number of iterations, 10000 at the
141coarsest resolution, 1000 at the medium resolution and 100 at the finest. These
142are the default settings
143"""
144
145level_iters = [10000, 1000, 100]
146
147"""
148To compute the Gaussian pyramid, the original image is first smoothed at each
149level of the pyramid using a Gaussian kernel with the requested sigma. A good
150initial choice is [3.0, 1.0, 0.0], this is the default
151"""
152
153sigmas = [3.0, 1.0, 0.0]
154
155"""
156Now we specify the sub-sampling factors. A good configuration is [4, 2, 1],
157which means that, if the original image shape was (nx, ny, nz) voxels, then the
158shape of the coarsest image will be about (nx//4, ny//4, nz//4), the shape in
159the middle resolution will be about (nx//2, ny//2, nz//2) and the image at the
160finest scale has the same size as the original image. This set of factors is
161the default
162"""
163
164factors = [4, 2, 1]
165
166"""
167Now we go ahead and instantiate the registration class with the configuration
168we just prepared
169"""
170
171affreg = AffineRegistration(metric=metric,
172                            level_iters=level_iters,
173                            sigmas=sigmas,
174                            factors=factors)
175
176"""
177Using AffineRegistration we can register our images in as many stages as we
178want, providing previous results as initialization for the next (the same logic
179as in ANTS). The reason why it is useful is that registration is a non-convex
180optimization problem (it may have more than one local optima), which means that
181it is very important to initialize as close to the solution as possible. For
182example, let's start with our (previously computed) rough transformation
183aligning the centers of mass of our images, and then refine it in three stages.
184First look for an optimal translation. The dictionary regtransforms contains
185all available transforms, we obtain one of them by providing its name and the
186dimension (either 2 or 3) of the image we are working with (since we are
187aligning volumes, the dimension is 3)
188"""
189
190transform = TranslationTransform3D()
191params0 = None
192starting_affine = c_of_mass.affine
193translation = affreg.optimize(static, moving, transform, params0,
194                              static_grid2world, moving_grid2world,
195                              starting_affine=starting_affine)
196
197"""
198If we look at the result, we can see that this translation is much better than
199simply aligning the centers of mass
200"""
201
202transformed = translation.transform(moving)
203regtools.overlay_slices(static, transformed, None, 0,
204                        "Static", "Transformed", "transformed_trans_0.png")
205regtools.overlay_slices(static, transformed, None, 1,
206                        "Static", "Transformed", "transformed_trans_1.png")
207regtools.overlay_slices(static, transformed, None, 2,
208                        "Static", "Transformed", "transformed_trans_2.png")
209
210"""
211.. figure:: transformed_trans_0.png
212   :align: center
213.. figure:: transformed_trans_1.png
214   :align: center
215.. figure:: transformed_trans_2.png
216   :align: center
217
218   Registration result by translating the moving image, using
219   Mutual Information.
220"""
221
222"""
223Now let's refine with a rigid transform (this may even modify our previously
224found optimal translation)
225"""
226
227transform = RigidTransform3D()
228params0 = None
229starting_affine = translation.affine
230rigid = affreg.optimize(static, moving, transform, params0,
231                        static_grid2world, moving_grid2world,
232                        starting_affine=starting_affine)
233
234"""
235This produces a slight rotation, and the images are now better aligned
236"""
237
238transformed = rigid.transform(moving)
239regtools.overlay_slices(static, transformed, None, 0,
240                        "Static", "Transformed", "transformed_rigid_0.png")
241regtools.overlay_slices(static, transformed, None, 1,
242                        "Static", "Transformed", "transformed_rigid_1.png")
243regtools.overlay_slices(static, transformed, None, 2,
244                        "Static", "Transformed", "transformed_rigid_2.png")
245
246"""
247.. figure:: transformed_rigid_0.png
248   :align: center
249.. figure:: transformed_rigid_1.png
250   :align: center
251.. figure:: transformed_rigid_2.png
252   :align: center
253
254   Registration result with a rigid transform, using Mutual Information.
255"""
256
257"""
258Finally, let's refine with a full affine transform (translation, rotation, scale
259and shear), it is safer to fit more degrees of freedom now since we must be
260very close to the optimal transform
261"""
262
263transform = AffineTransform3D()
264params0 = None
265starting_affine = rigid.affine
266affine = affreg.optimize(static, moving, transform, params0,
267                         static_grid2world, moving_grid2world,
268                         starting_affine=starting_affine)
269
270"""
271This results in a slight shear and scale
272"""
273
274transformed = affine.transform(moving)
275regtools.overlay_slices(static, transformed, None, 0,
276                        "Static", "Transformed", "transformed_affine_0.png")
277regtools.overlay_slices(static, transformed, None, 1,
278                        "Static", "Transformed", "transformed_affine_1.png")
279regtools.overlay_slices(static, transformed, None, 2,
280                        "Static", "Transformed", "transformed_affine_2.png")
281
282"""
283.. figure:: transformed_affine_0.png
284   :align: center
285.. figure:: transformed_affine_1.png
286   :align: center
287.. figure:: transformed_affine_2.png
288   :align: center
289
290   Registration result with an affine transform, using Mutual Information.
291
292"""
293
294"""
295Now, let's repeat this process with a simplified functional interface:
296"""
297
298from dipy.align import (affine_registration, center_of_mass, translation,
299                        rigid, affine, register_dwi_to_template)
300
301"""
302This interface constructs a pipeline of operations as a sequence of functions
303that each implement one of the transforms.
304"""
305
306pipeline = [center_of_mass, translation, rigid, affine]
307
308"""
309And then applies the functions in the pipeline on the input (from left to
310right) with a call to an `affine_registration` function, which takes optional
311settings for things like the iterations, sigmas and factors.
312"""
313
314xformed_img, reg_affine = affine_registration(
315    moving,
316    static,
317    moving_affine=moving_affine,
318    static_affine=static_affine,
319    nbins=32,
320    metric='MI',
321    pipeline=pipeline,
322    level_iters=level_iters,
323    sigmas=sigmas,
324    factors=factors)
325
326regtools.overlay_slices(static, xformed_img, None, 0,
327                        "Static", "Transformed", "xformed_affine_0.png")
328regtools.overlay_slices(static, xformed_img, None, 1,
329                        "Static", "Transformed", "xformed_affine_1.png")
330regtools.overlay_slices(static, xformed_img, None, 2,
331                        "Static", "Transformed", "xformed_affine_2.png")
332
333
334"""
335
336.. figure:: xformed_affine_0.png
337   :align: center
338.. figure:: xformed_affine_1.png
339   :align: center
340.. figure:: xformed_affine_2.png
341   :align: center
342
343   Registration result with an affine transform, using functional interface.
344
345"""
346
347"""
348Alternatively, you can also use the `register_dwi_to_template` function that
349needs to also know about the gradient table of the DWI data, provided as a
350tuple of (bvals_file, bvecs_file). In this case, we are going to move the
351diffusion data to the B0 image (the opposite of the previous examples), which
352reverses what is the "moving" image and what is "static".
353
354"""
355
356xformed_dwi, reg_affine = register_dwi_to_template(
357    dwi=static_img,
358    gtab=(pjoin(folder, 'HARDI150.bval'),
359          pjoin(folder, 'HARDI150.bvec')),
360    template=moving_img,
361    reg_method="aff",
362    nbins=32,
363    metric='MI',
364    pipeline=pipeline,
365    level_iters=level_iters,
366    sigmas=sigmas,
367    factors=factors)
368
369regtools.overlay_slices(moving, xformed_dwi, None, 0,
370                        "Static", "Transformed", "xformed_dwi_0.png")
371regtools.overlay_slices(moving, xformed_dwi, None, 1,
372                        "Static", "Transformed", "xformed_dwi_1.png")
373regtools.overlay_slices(moving, xformed_dwi, None, 2,
374                        "Static", "Transformed", "xformed_dwi_2.png")
375
376
377
378"""
379
380.. figure:: xformed_dwi_0.png
381   :align: center
382.. figure:: xformed_dwi_1.png
383   :align: center
384.. figure:: xformed_dwi_2.png
385   :align: center
386
387   Same again, using the `dwi_to_template` functional interface.
388
389
390.. [Mattes03] Mattes, D., Haynor, D. R., Vesselle, H., Lewellen, T. K.,
391              Eubank, W. (2003). PET-CT image registration in the chest using
392              free-form deformations. IEEE Transactions on Medical Imaging,
393              22(1), 120-8.
394.. [Avants11] Avants, B. B., Tustison, N., & Song, G. (2011). Advanced
395              Normalization Tools (ANTS), 1-35.
396
397.. include:: ../links_names.inc
398
399"""
400