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