1""" Affine image registration module consisting of the following classes: 2 3 AffineMap: encapsulates the necessary information to perform affine 4 transforms between two domains, defined by a `static` and a `moving` 5 image. The `domain` of the transform is the set of points in the 6 `static` image's grid, and the `codomain` is the set of points in 7 the `moving` image. When we call the `transform` method, `AffineMap` 8 maps each point `x` of the domain (`static` grid) to the codomain 9 (`moving` grid) and interpolates the `moving` image at that point 10 to obtain the intensity value to be placed at `x` in the resulting 11 grid. The `transform_inverse` method performs the opposite operation 12 mapping points in the codomain to points in the domain. 13 14 ParzenJointHistogram: computes the marginal and joint distributions of 15 intensities of a pair of images, using Parzen windows [Parzen62] 16 with a cubic spline kernel, as proposed by Mattes et al. [Mattes03]. 17 It also computes the gradient of the joint histogram w.r.t. the 18 parameters of a given transform. 19 20 MutualInformationMetric: computes the value and gradient of the mutual 21 information metric the way `Optimizer` needs them. That is, given 22 a set of transform parameters, it will use `ParzenJointHistogram` 23 to compute the value and gradient of the joint intensity histogram 24 evaluated at the given parameters, and evaluate the the value and 25 gradient of the histogram's mutual information. 26 27 AffineRegistration: it runs the multi-resolution registration, putting 28 all the pieces together. It needs to create the scale space of the 29 images and run the multi-resolution registration by using the Metric 30 and the Optimizer at each level of the Gaussian pyramid. At each 31 level, it will setup the metric to compute value and gradient of the 32 metric with the input images with different levels of smoothing. 33 34 References 35 ---------- 36 [Parzen62] E. Parzen. On the estimation of a probability density 37 function and the mode. Annals of Mathematical Statistics, 38 33(3), 1065-1076, 1962. 39 [Mattes03] Mattes, D., Haynor, D. R., Vesselle, H., Lewellen, T. K., 40 & Eubank, W. PET-CT image registration in the chest using 41 free-form deformations. IEEE Transactions on Medical 42 Imaging, 22(1), 120-8, 2003. 43 44""" 45 46import numpy as np 47import numpy.linalg as npl 48import scipy.ndimage as ndimage 49from dipy.core.optimize import Optimizer 50from dipy.core.interpolation import (interpolate_scalar_2d, 51 interpolate_scalar_3d) 52from dipy.align import vector_fields as vf 53from dipy.align import VerbosityLevels 54from dipy.align.parzenhist import (ParzenJointHistogram, 55 sample_domain_regular, 56 compute_parzen_mi) 57from dipy.align.imwarp import (get_direction_and_spacings, ScaleSpace) 58from dipy.align.scalespace import IsotropicScaleSpace 59from dipy.utils.deprecator import deprecated_params 60 61_interp_options = ['nearest', 'linear'] 62_transform_method = {} 63_transform_method[(2, 'nearest')] = vf.transform_2d_affine_nn 64_transform_method[(3, 'nearest')] = vf.transform_3d_affine_nn 65_transform_method[(2, 'linear')] = vf.transform_2d_affine 66_transform_method[(3, 'linear')] = vf.transform_3d_affine 67_number_dim_affine_matrix = 2 68 69 70class AffineInversionError(Exception): 71 pass 72 73 74class AffineInvalidValuesError(Exception): 75 pass 76 77 78class AffineMap(object): 79 80 def __init__(self, affine, domain_grid_shape=None, domain_grid2world=None, 81 codomain_grid_shape=None, codomain_grid2world=None): 82 """ AffineMap 83 84 Implements an affine transformation whose domain is given by 85 `domain_grid` and `domain_grid2world`, and whose co-domain is 86 given by `codomain_grid` and `codomain_grid2world`. 87 88 The actual transform is represented by the `affine` matrix, which 89 operate in world coordinates. Therefore, to transform a moving image 90 towards a static image, we first map each voxel (i,j,k) of the static 91 image to world coordinates (x,y,z) by applying `domain_grid2world`. 92 Then we apply the `affine` transform to (x,y,z) obtaining (x', y', z') 93 in moving image's world coordinates. Finally, (x', y', z') is mapped 94 to voxel coordinates (i', j', k') in the moving image by multiplying 95 (x', y', z') by the inverse of `codomain_grid2world`. The 96 `codomain_grid_shape` is used analogously to transform the static 97 image towards the moving image when calling `transform_inverse`. 98 99 If the domain/co-domain information is not provided (None) then the 100 sampling information needs to be specified each time the `transform` 101 or `transform_inverse` is called to transform images. Note that such 102 sampling information is not necessary to transform points defined in 103 physical space, such as stream lines. 104 105 Parameters 106 ---------- 107 affine : array, shape (dim + 1, dim + 1) 108 the matrix defining the affine transform, where `dim` is the 109 dimension of the space this map operates in (2 for 2D images, 110 3 for 3D images). If None, then `self` represents the identity 111 transformation. 112 domain_grid_shape : sequence, shape (dim,), optional 113 the shape of the default domain sampling grid. When `transform` 114 is called to transform an image, the resulting image will have 115 this shape, unless a different sampling information is provided. 116 If None, then the sampling grid shape must be specified each time 117 the `transform` method is called. 118 domain_grid2world : array, shape (dim + 1, dim + 1), optional 119 the grid-to-world transform associated with the domain grid. 120 If None (the default), then the grid-to-world transform is assumed 121 to be the identity. 122 codomain_grid_shape : sequence of integers, shape (dim,) 123 the shape of the default co-domain sampling grid. When 124 `transform_inverse` is called to transform an image, the resulting 125 image will have this shape, unless a different sampling 126 information is provided. If None (the default), then the sampling 127 grid shape must be specified each time the `transform_inverse` 128 method is called. 129 codomain_grid2world : array, shape (dim + 1, dim + 1) 130 the grid-to-world transform associated with the co-domain grid. 131 If None (the default), then the grid-to-world transform is assumed 132 to be the identity. 133 134 """ 135 self.set_affine(affine) 136 self.domain_shape = domain_grid_shape 137 self.domain_grid2world = domain_grid2world 138 self.codomain_shape = codomain_grid_shape 139 self.codomain_grid2world = codomain_grid2world 140 141 def get_affine(self): 142 """Return the value of the transformation, not a reference. 143 144 Returns 145 ------- 146 affine : ndarray 147 Copy of the transform, not a reference. 148 149 """ 150 151 # returning a copy to insulate it from changes outside object 152 return self.affine.copy() 153 154 def set_affine(self, affine): 155 """Set the affine transform (operating in physical space). 156 157 Also sets `self.affine_inv` - the inverse of `affine`, or None if 158 there is no inverse. 159 160 Parameters 161 ---------- 162 affine : array, shape (dim + 1, dim + 1) 163 the matrix representing the affine transform operating in 164 physical space. The domain and co-domain information 165 remains unchanged. If None, then `self` represents the identity 166 transformation. 167 168 """ 169 170 if affine is None: 171 self.affine = None 172 self.affine_inv = None 173 return 174 175 try: 176 affine = np.array(affine) 177 except Exception: 178 raise TypeError("Input must be type ndarray, or be convertible" 179 " to one.") 180 181 if len(affine.shape) != _number_dim_affine_matrix: 182 raise AffineInversionError('Affine transform must be 2D') 183 184 if not affine.shape[0] == affine.shape[1]: 185 raise AffineInversionError("Affine transform must be a square " 186 "matrix") 187 188 if not np.all(np.isfinite(affine)): 189 raise AffineInvalidValuesError("Affine transform contains invalid" 190 " elements") 191 192 # checking on proper augmentation 193 # First n-1 columns in last row in matrix contain non-zeros 194 if not np.all(affine[-1, :-1] == 0.0): 195 raise AffineInvalidValuesError("First {n_1} columns in last row" 196 " in matrix contain non-zeros!" 197 .format(n_1=affine.shape[0] - 1)) 198 199 # Last row, last column in matrix must be 1.0! 200 if affine[-1, -1] != 1.0: 201 raise AffineInvalidValuesError("Last row, last column in matrix" 202 " is not 1.0!") 203 204 # making a copy to insulate it from changes outside object 205 self.affine = affine.copy() 206 207 try: 208 self.affine_inv = npl.inv(affine) 209 except npl.LinAlgError: 210 raise AffineInversionError('Affine cannot be inverted') 211 212 def __str__(self): 213 """Printable format - relies on ndarray's implementation.""" 214 215 return str(self.affine) 216 217 def __repr__(self): 218 """Relodable representation - also relies on ndarray's implementation.""" 219 220 return self.affine.__repr__() 221 222 def __format__(self, format_spec): 223 """Implementation various formatting options""" 224 225 if format_spec is None or self.affine is None: 226 return str(self.affine) 227 elif isinstance(format_spec, str): 228 format_spec = format_spec.lower() 229 if format_spec in ['', ' ', 'f', 'full']: 230 return str(self.affine) 231 # rotation part only (initial 3x3) 232 elif format_spec in ['r', 'rotation']: 233 return str(self.affine[:-1, :-1]) 234 # translation part only (4th col) 235 elif format_spec in ['t', 'translation']: 236 # notice unusual indexing to make it a column vector 237 # i.e. rows from 0 to n-1, cols from n to n 238 return str(self.affine[:-1, -1:]) 239 else: 240 allowed_formats_print_map = ['full', 'f', 241 'rotation', 'r', 242 'translation', 't'] 243 raise NotImplementedError("Format {} not recognized or" 244 "implemented.\nTry one of {}" 245 .format(format_spec, 246 allowed_formats_print_map)) 247 248 @deprecated_params('interp', 'interpolation', since='1.13', until='1.15') 249 def _apply_transform(self, image, interpolation='linear', 250 image_grid2world=None, sampling_grid_shape=None, 251 sampling_grid2world=None, resample_only=False, 252 apply_inverse=False): 253 """Transform the input image applying this affine transform. 254 255 This is a generic function to transform images using either this 256 (direct) transform or its inverse. 257 258 If applying the direct transform (`apply_inverse=False`): 259 by default, the transformed image is sampled at a grid defined by 260 `self.domain_shape` and `self.domain_grid2world`. 261 If applying the inverse transform (`apply_inverse=True`): 262 by default, the transformed image is sampled at a grid defined by 263 `self.codomain_shape` and `self.codomain_grid2world`. 264 265 If the sampling information was not provided at initialization of this 266 transform then `sampling_grid_shape` is mandatory. 267 268 Parameters 269 ---------- 270 image : 2D or 3D array 271 the image to be transformed 272 interpolation : string, either 'linear' or 'nearest' 273 the type of interpolation to be used, either 'linear' 274 (for k-linear interpolation) or 'nearest' for nearest neighbor 275 image_grid2world : array, shape (dim + 1, dim + 1), optional 276 the grid-to-world transform associated with `image`. 277 If None (the default), then the grid-to-world transform is assumed 278 to be the identity. 279 sampling_grid_shape : sequence, shape (dim,), optional 280 the shape of the grid where the transformed image must be sampled. 281 If None (the default), then `self.domain_shape` is used instead 282 (which must have been set at initialization, otherwise an exception 283 will be raised). 284 sampling_grid2world : array, shape (dim + 1, dim + 1), optional 285 the grid-to-world transform associated with the sampling grid 286 (specified by `sampling_grid_shape`, or by default 287 `self.domain_shape`). If None (the default), then the 288 grid-to-world transform is assumed to be the identity. 289 resample_only : Boolean, optional 290 If False (the default) the affine transform is applied normally. 291 If True, then the affine transform is not applied, and the input 292 image is just re-sampled on the domain grid of this transform. 293 apply_inverse : Boolean, optional 294 If False (the default) the image is transformed from the codomain 295 of this transform to its domain using the (direct) affine 296 transform. Otherwise, the image is transformed from the domain 297 of this transform to its codomain using the (inverse) affine 298 transform. 299 300 Returns 301 ------- 302 transformed : array, shape `sampling_grid_shape` or `self.domain_shape` 303 the transformed image, sampled at the requested grid 304 305 """ 306 # Verify valid interpolation requested 307 if interpolation not in _interp_options: 308 msg = 'Unknown interpolation method: %s' % (interpolation,) 309 raise ValueError(msg) 310 311 # Obtain sampling grid 312 if sampling_grid_shape is None: 313 if apply_inverse: 314 sampling_grid_shape = self.codomain_shape 315 else: 316 sampling_grid_shape = self.domain_shape 317 if sampling_grid_shape is None: 318 msg = 'Unknown sampling info. Provide a valid sampling_grid_shape' 319 raise ValueError(msg) 320 321 dim = len(sampling_grid_shape) 322 shape = np.array(sampling_grid_shape, dtype=np.int32) 323 324 # Verify valid image dimension 325 img_dim = len(image.shape) 326 if img_dim < 2 or img_dim > 3: 327 raise ValueError('Undefined transform for dim: %d' % (img_dim,)) 328 329 # Obtain grid-to-world transform for sampling grid 330 if sampling_grid2world is None: 331 if apply_inverse: 332 sampling_grid2world = self.codomain_grid2world 333 else: 334 sampling_grid2world = self.domain_grid2world 335 if sampling_grid2world is None: 336 sampling_grid2world = np.eye(dim + 1) 337 338 # Obtain world-to-grid transform for input image 339 if image_grid2world is None: 340 if apply_inverse: 341 image_grid2world = self.domain_grid2world 342 else: 343 image_grid2world = self.codomain_grid2world 344 if image_grid2world is None: 345 image_grid2world = np.eye(dim + 1) 346 image_world2grid = npl.inv(image_grid2world) 347 348 # Compute the transform from sampling grid to input image grid 349 if apply_inverse: 350 aff = self.affine_inv 351 else: 352 aff = self.affine 353 354 if (aff is None) or resample_only: 355 comp = image_world2grid.dot(sampling_grid2world) 356 else: 357 comp = image_world2grid.dot(aff.dot(sampling_grid2world)) 358 359 # Transform the input image 360 if interpolation == 'linear': 361 image = image.astype(np.float64) 362 transformed = _transform_method[(dim, interpolation)](image, shape, 363 comp) 364 return transformed 365 366 @deprecated_params('interp', 'interpolation', since='1.13', until='1.15') 367 def transform(self, image, interpolation='linear', image_grid2world=None, 368 sampling_grid_shape=None, sampling_grid2world=None, 369 resample_only=False): 370 """Transform the input image from co-domain to domain space. 371 372 By default, the transformed image is sampled at a grid defined by 373 `self.domain_shape` and `self.domain_grid2world`. If such 374 information was not provided then `sampling_grid_shape` is mandatory. 375 376 Parameters 377 ---------- 378 image : 2D or 3D array 379 the image to be transformed 380 interpolation : string, either 'linear' or 'nearest' 381 the type of interpolation to be used, either 'linear' 382 (for k-linear interpolation) or 'nearest' for nearest neighbor 383 image_grid2world : array, shape (dim + 1, dim + 1), optional 384 the grid-to-world transform associated with `image`. 385 If None (the default), then the grid-to-world transform is assumed 386 to be the identity. 387 sampling_grid_shape : sequence, shape (dim,), optional 388 the shape of the grid where the transformed image must be sampled. 389 If None (the default), then `self.codomain_shape` is used instead 390 (which must have been set at initialization, otherwise an exception 391 will be raised). 392 sampling_grid2world : array, shape (dim + 1, dim + 1), optional 393 the grid-to-world transform associated with the sampling grid 394 (specified by `sampling_grid_shape`, or by default 395 `self.codomain_shape`). If None (the default), then the 396 grid-to-world transform is assumed to be the identity. 397 resample_only : Boolean, optional 398 If False (the default) the affine transform is applied normally. 399 If True, then the affine transform is not applied, and the input 400 image is just re-sampled on the domain grid of this transform. 401 402 Returns 403 ------- 404 transformed : array, shape `sampling_grid_shape` or 405 `self.codomain_shape` 406 the transformed image, sampled at the requested grid 407 408 """ 409 transformed = self._apply_transform(image, interpolation, 410 image_grid2world, 411 sampling_grid_shape, 412 sampling_grid2world, 413 resample_only, 414 apply_inverse=False) 415 return np.array(transformed) 416 417 @deprecated_params('interp', 'interpolation', since='1.13', until='1.15') 418 def transform_inverse(self, image, interpolation='linear', 419 image_grid2world=None, sampling_grid_shape=None, 420 sampling_grid2world=None, resample_only=False): 421 """Transform the input image from domain to co-domain space. 422 423 By default, the transformed image is sampled at a grid defined by 424 `self.codomain_shape` and `self.codomain_grid2world`. If such 425 information was not provided then `sampling_grid_shape` is mandatory. 426 427 Parameters 428 ---------- 429 image : 2D or 3D array 430 the image to be transformed 431 interpolation : string, either 'linear' or 'nearest' 432 the type of interpolation to be used, either 'linear' 433 (for k-linear interpolation) or 'nearest' for nearest neighbor 434 image_grid2world : array, shape (dim + 1, dim + 1), optional 435 the grid-to-world transform associated with `image`. 436 If None (the default), then the grid-to-world transform is assumed 437 to be the identity. 438 sampling_grid_shape : sequence, shape (dim,), optional 439 the shape of the grid where the transformed image must be sampled. 440 If None (the default), then `self.codomain_shape` is used instead 441 (which must have been set at initialization, otherwise an exception 442 will be raised). 443 sampling_grid2world : array, shape (dim + 1, dim + 1), optional 444 the grid-to-world transform associated with the sampling grid 445 (specified by `sampling_grid_shape`, or by default 446 `self.codomain_shape`). If None (the default), then the 447 grid-to-world transform is assumed to be the identity. 448 resample_only : Boolean, optional 449 If False (the default) the affine transform is applied normally. 450 If True, then the affine transform is not applied, and the input 451 image is just re-sampled on the domain grid of this transform. 452 453 Returns 454 ------- 455 transformed : array, shape `sampling_grid_shape` or 456 `self.codomain_shape` 457 the transformed image, sampled at the requested grid 458 459 """ 460 transformed = self._apply_transform(image, interpolation, 461 image_grid2world, 462 sampling_grid_shape, 463 sampling_grid2world, 464 resample_only, 465 apply_inverse=True) 466 return np.array(transformed) 467 468 469class MutualInformationMetric(object): 470 471 def __init__(self, nbins=32, sampling_proportion=None): 472 r"""Initialize an instance of the Mutual Information metric. 473 474 This class implements the methods required by Optimizer to drive the 475 registration process. 476 477 Parameters 478 ---------- 479 nbins : int, optional 480 the number of bins to be used for computing the intensity 481 histograms. The default is 32. 482 sampling_proportion : None or float in interval (0, 1], optional 483 There are two types of sampling: dense and sparse. Dense sampling 484 uses all voxels for estimating the (joint and marginal) intensity 485 histograms, while sparse sampling uses a subset of them. If 486 `sampling_proportion` is None, then dense sampling is 487 used. If `sampling_proportion` is a floating point value in (0,1] 488 then sparse sampling is used, where `sampling_proportion` 489 specifies the proportion of voxels to be used. The default is 490 None. 491 492 Notes 493 ----- 494 Since we use linear interpolation, images are not, in general, 495 differentiable at exact voxel coordinates, but they are differentiable 496 between voxel coordinates. When using sparse sampling, selected voxels 497 are slightly moved by adding a small random displacement within one 498 voxel to prevent sampling points from being located exactly at voxel 499 coordinates. When using dense sampling, this random displacement is 500 not applied. 501 502 """ 503 self.histogram = ParzenJointHistogram(nbins) 504 self.sampling_proportion = sampling_proportion 505 self.metric_val = None 506 self.metric_grad = None 507 508 def setup(self, transform, static, moving, static_grid2world=None, 509 moving_grid2world=None, starting_affine=None): 510 r"""Prepare the metric to compute intensity densities and gradients. 511 512 The histograms will be setup to compute probability densities of 513 intensities within the minimum and maximum values of `static` and 514 `moving` 515 516 Parameters 517 ---------- 518 transform: instance of Transform 519 the transformation with respect to whose parameters the gradient 520 must be computed 521 static : array, shape (S, R, C) or (R, C) 522 static image 523 moving : array, shape (S', R', C') or (R', C') 524 moving image. The dimensions of the static (S, R, C) and moving 525 (S', R', C') images do not need to be the same. 526 static_grid2world : array (dim+1, dim+1), optional 527 the grid-to-space transform of the static image. The default is 528 None, implying the transform is the identity. 529 moving_grid2world : array (dim+1, dim+1) 530 the grid-to-space transform of the moving image. The default is 531 None, implying the spacing along all axes is 1. 532 starting_affine : array, shape (dim+1, dim+1), optional 533 the pre-aligning matrix (an affine transform) that roughly aligns 534 the moving image towards the static image. If None, no 535 pre-alignment is performed. If a pre-alignment matrix is available, 536 it is recommended to provide this matrix as `starting_affine` 537 instead of manually transforming the moving image to reduce 538 interpolation artifacts. The default is None, implying no 539 pre-alignment is performed. 540 541 """ 542 n = transform.get_number_of_parameters() 543 self.metric_grad = np.zeros(n, dtype=np.float64) 544 self.dim = len(static.shape) 545 if moving_grid2world is None: 546 moving_grid2world = np.eye(self.dim + 1) 547 if static_grid2world is None: 548 static_grid2world = np.eye(self.dim + 1) 549 self.transform = transform 550 self.static = np.array(static).astype(np.float64) 551 self.moving = np.array(moving).astype(np.float64) 552 self.static_grid2world = static_grid2world 553 self.static_world2grid = npl.inv(static_grid2world) 554 self.moving_grid2world = moving_grid2world 555 self.moving_world2grid = npl.inv(moving_grid2world) 556 self.static_direction, self.static_spacing = \ 557 get_direction_and_spacings(static_grid2world, self.dim) 558 self.moving_direction, self.moving_spacing = \ 559 get_direction_and_spacings(moving_grid2world, self.dim) 560 self.starting_affine = starting_affine 561 562 P = np.eye(self.dim + 1) 563 if self.starting_affine is not None: 564 P = self.starting_affine 565 566 self.affine_map = AffineMap(P, static.shape, static_grid2world, 567 moving.shape, moving_grid2world) 568 569 if self.dim == 2: 570 self.interp_method = interpolate_scalar_2d 571 else: 572 self.interp_method = interpolate_scalar_3d 573 574 if self.sampling_proportion is None: 575 self.samples = None 576 self.ns = 0 577 else: 578 k = int(np.ceil(1.0 / self.sampling_proportion)) 579 shape = np.array(static.shape, dtype=np.int32) 580 self.samples = sample_domain_regular(k, shape, static_grid2world) 581 self.samples = np.array(self.samples) 582 self.ns = self.samples.shape[0] 583 # Add a column of ones (homogeneous coordinates) 584 self.samples = np.hstack((self.samples, np.ones(self.ns)[:, None])) 585 if self.starting_affine is None: 586 self.samples_prealigned = self.samples 587 else: 588 self.samples_prealigned = \ 589 self.starting_affine.dot(self.samples.T).T 590 # Sample the static image 591 static_p = self.static_world2grid.dot(self.samples.T).T 592 static_p = static_p[..., :self.dim] 593 self.static_vals, inside = self.interp_method(static, static_p) 594 self.static_vals = np.array(self.static_vals, dtype=np.float64) 595 self.histogram.setup(self.static, self.moving) 596 597 def _update_histogram(self): 598 r"""Update the histogram according to the current affine transform. 599 600 The current affine transform is given by `self.affine_map`, which 601 must be set before calling this method. 602 603 Returns 604 ------- 605 static_values: array, shape(n,) if sparse sampling is being used, 606 array, shape(S, R, C) or (R, C) if dense sampling 607 the intensity values corresponding to the static image used to 608 update the histogram. If sparse sampling is being used, then 609 it is simply a sequence of scalars, obtained by sampling the static 610 image at the `n` sampling points. If dense sampling is being used, 611 then the intensities are given directly by the static image, 612 whose shape is (S, R, C) in the 3D case or (R, C) in the 2D case. 613 moving_values: array, shape(n,) if sparse sampling is being used, 614 array, shape(S, R, C) or (R, C) if dense sampling 615 the intensity values corresponding to the moving image used to 616 update the histogram. If sparse sampling is being used, then 617 it is simply a sequence of scalars, obtained by sampling the moving 618 image at the `n` sampling points (mapped to the moving space by the 619 current affine transform). If dense sampling is being used, 620 then the intensities are given by the moving imaged linearly 621 transformed towards the static image by the current affine, which 622 results in an image of the same shape as the static image. 623 624 """ 625 if self.sampling_proportion is None: # Dense case 626 static_values = self.static 627 moving_values = self.affine_map.transform(self.moving) 628 self.histogram.update_pdfs_dense(static_values, moving_values) 629 else: # Sparse case 630 sp_to_moving = self.moving_world2grid.dot(self.affine_map.affine) 631 pts = sp_to_moving.dot(self.samples.T).T # Points on moving grid 632 pts = pts[..., :self.dim] 633 self.moving_vals, inside = self.interp_method(self.moving, pts) 634 self.moving_vals = np.array(self.moving_vals) 635 static_values = self.static_vals 636 moving_values = self.moving_vals 637 self.histogram.update_pdfs_sparse(static_values, moving_values) 638 return static_values, moving_values 639 640 def _update_mutual_information(self, params, update_gradient=True): 641 r"""Update marginal and joint distributions and the joint gradient. 642 643 The distributions are updated according to the static and transformed 644 images. The transformed image is precisely the moving image after 645 transforming it by the transform defined by the `params` parameters. 646 647 The gradient of the joint PDF is computed only if update_gradient 648 is True. 649 650 Parameters 651 ---------- 652 params : array, shape (n,) 653 the parameter vector of the transform currently used by the metric 654 (the transform name is provided when self.setup is called), n is 655 the number of parameters of the transform 656 update_gradient : Boolean, optional 657 if True, the gradient of the joint PDF will also be computed, 658 otherwise, only the marginal and joint PDFs will be computed. 659 The default is True. 660 661 """ 662 # Get the matrix associated with the `params` parameter vector 663 current_affine = self.transform.param_to_matrix(params) 664 # Get the static-to-prealigned matrix (only needed for the MI gradient) 665 static2prealigned = self.static_grid2world 666 if self.starting_affine is not None: 667 current_affine = current_affine.dot(self.starting_affine) 668 static2prealigned = self.starting_affine.dot(static2prealigned) 669 self.affine_map.set_affine(current_affine) 670 671 # Update the histogram with the current joint intensities 672 static_values, moving_values = self._update_histogram() 673 674 H = self.histogram # Shortcut to `self.histogram` 675 grad = None # Buffer to write the MI gradient into (if needed) 676 if update_gradient: 677 grad = self.metric_grad 678 # Compute the gradient of the joint PDF w.r.t. parameters 679 if self.sampling_proportion is None: # Dense case 680 # Compute the gradient of moving img. at physical points 681 # associated with the >>static image's grid<< cells 682 # The image gradient must be eval. at current moved points 683 grid_to_world = current_affine.dot(self.static_grid2world) 684 mgrad, inside = vf.gradient(self.moving, 685 self.moving_world2grid, 686 self.moving_spacing, 687 self.static.shape, 688 grid_to_world) 689 # The Jacobian must be evaluated at the pre-aligned points 690 H.update_gradient_dense( 691 params, 692 self.transform, 693 static_values, 694 moving_values, 695 static2prealigned, 696 mgrad) 697 else: # Sparse case 698 # Compute the gradient of moving at the sampling points 699 # which are already given in physical space coordinates 700 pts = current_affine.dot(self.samples.T).T # Moved points 701 mgrad, inside = vf.sparse_gradient(self.moving, 702 self.moving_world2grid, 703 self.moving_spacing, 704 pts) 705 # The Jacobian must be evaluated at the pre-aligned points 706 pts = self.samples_prealigned[..., :self.dim] 707 H.update_gradient_sparse(params, self.transform, static_values, 708 moving_values, pts, mgrad) 709 710 # Call the cythonized MI computation with self.histogram fields 711 self.metric_val = compute_parzen_mi(H.joint, H.joint_grad, 712 H.smarginal, H.mmarginal, 713 grad) 714 715 def distance(self, params): 716 r"""Numeric value of the negative Mutual Information. 717 718 We need to change the sign so we can use standard minimization 719 algorithms. 720 721 Parameters 722 ---------- 723 params : array, shape (n,) 724 the parameter vector of the transform currently used by the metric 725 (the transform name is provided when self.setup is called), n is 726 the number of parameters of the transform 727 728 Returns 729 ------- 730 neg_mi : float 731 the negative mutual information of the input images after 732 transforming the moving image by the currently set transform 733 with `params` parameters 734 735 """ 736 try: 737 self._update_mutual_information(params, False) 738 except (AffineInversionError, AffineInvalidValuesError): 739 return np.inf 740 return -1 * self.metric_val 741 742 def gradient(self, params): 743 r"""Numeric value of the metric's gradient at the given parameters. 744 745 Parameters 746 ---------- 747 params : array, shape (n,) 748 the parameter vector of the transform currently used by the metric 749 (the transform name is provided when self.setup is called), n is 750 the number of parameters of the transform 751 752 Returns 753 ------- 754 grad : array, shape (n,) 755 the gradient of the negative Mutual Information 756 757 """ 758 try: 759 self._update_mutual_information(params, True) 760 except (AffineInversionError, AffineInvalidValuesError): 761 return 0 * self.metric_grad 762 return -1 * self.metric_grad 763 764 def distance_and_gradient(self, params): 765 r"""Numeric value of the metric and its gradient at given parameters. 766 767 Parameters 768 ---------- 769 params : array, shape (n,) 770 the parameter vector of the transform currently used by the metric 771 (the transform name is provided when self.setup is called), n is 772 the number of parameters of the transform 773 774 Returns 775 ------- 776 neg_mi : float 777 the negative mutual information of the input images after 778 transforming the moving image by the currently set transform 779 with `params` parameters 780 neg_mi_grad : array, shape (n,) 781 the gradient of the negative Mutual Information 782 783 """ 784 try: 785 self._update_mutual_information(params, True) 786 except (AffineInversionError, AffineInvalidValuesError): 787 return np.inf, 0 * self.metric_grad 788 return -1 * self.metric_val, -1 * self.metric_grad 789 790 791class AffineRegistration(object): 792 793 def __init__(self, 794 metric=None, 795 level_iters=None, 796 sigmas=None, 797 factors=None, 798 method='L-BFGS-B', 799 ss_sigma_factor=None, 800 options=None, 801 verbosity=VerbosityLevels.STATUS): 802 """Initialize an instance of the AffineRegistration class. 803 804 Parameters 805 ---------- 806 metric : None or object, optional 807 an instance of a metric. The default is None, implying 808 the Mutual Information metric with default settings. 809 level_iters : sequence, optional 810 the number of iterations at each scale of the scale space. 811 `level_iters[0]` corresponds to the coarsest scale, 812 `level_iters[-1]` the finest, where n is the length of the 813 sequence. By default, a 3-level scale space with iterations 814 sequence equal to [10000, 1000, 100] will be used. 815 sigmas : sequence of floats, optional 816 custom smoothing parameter to build the scale space (one parameter 817 for each scale). By default, the sequence of sigmas will be 818 [3, 1, 0]. 819 factors : sequence of floats, optional 820 custom scale factors to build the scale space (one factor for each 821 scale). By default, the sequence of factors will be [4, 2, 1]. 822 method : string, optional 823 optimization method to be used. If Scipy version < 0.12, then 824 only L-BFGS-B is available. Otherwise, `method` can be any 825 gradient-based method available in `dipy.core.Optimize`: CG, BFGS, 826 Newton-CG, dogleg or trust-ncg. 827 The default is 'L-BFGS-B'. 828 ss_sigma_factor : float, optional 829 If None, this parameter is not used and an isotropic scale 830 space with the given `factors` and `sigmas` will be built. 831 If not None, an anisotropic scale space will be used by 832 automatically selecting the smoothing sigmas along each axis 833 according to the voxel dimensions of the given image. 834 The `ss_sigma_factor` is used to scale the automatically computed 835 sigmas. For example, in the isotropic case, the sigma of the 836 kernel will be $factor * (2 ^ i)$ where 837 $i = 1, 2, ..., n_scales - 1$ is the scale (the finest resolution 838 image $i=0$ is never smoothed). The default is None. 839 options : dict, optional 840 extra optimization options. The default is None, implying 841 no extra options are passed to the optimizer. 842 843 """ 844 self.metric = metric 845 846 if self.metric is None: 847 self.metric = MutualInformationMetric() 848 849 if level_iters is None: 850 level_iters = [10000, 1000, 100] 851 self.level_iters = level_iters 852 self.levels = len(level_iters) 853 if self.levels == 0: 854 raise ValueError('The iterations sequence cannot be empty') 855 856 self.options = options 857 self.method = method 858 if ss_sigma_factor is not None: 859 self.use_isotropic = False 860 self.ss_sigma_factor = ss_sigma_factor 861 else: 862 self.use_isotropic = True 863 if factors is None: 864 factors = [4, 2, 1] 865 if sigmas is None: 866 sigmas = [3, 1, 0] 867 self.factors = factors 868 self.sigmas = sigmas 869 870 self.verbosity = verbosity 871 872 # Separately add a string that tells about the verbosity kwarg. This needs 873 # to be separate, because it is set as a module-wide option in __init__: 874 docstring_addendum = \ 875 """verbosity: int (one of {0, 1, 2, 3}), optional 876 Set the verbosity level of the algorithm: 877 0 : do not print anything 878 1 : print information about the current status of the algorithm 879 2 : print high level information of the components involved in 880 the registration that can be used to detect a failing 881 component. 882 3 : print as much information as possible to isolate the cause 883 of a bug. 884 Default: % s 885 """ % VerbosityLevels.STATUS 886 887 __init__.__doc__ = __init__.__doc__ + docstring_addendum 888 889 def _init_optimizer(self, static, moving, transform, params0, 890 static_grid2world, moving_grid2world, 891 starting_affine): 892 r"""Initialize the registration optimizer. 893 894 Initializes the optimizer by computing the scale space of the input 895 images 896 897 Parameters 898 ---------- 899 static : array, shape (S, R, C) or (R, C) 900 the image to be used as reference during optimization. 901 moving : array, shape (S', R', C') or (R', C') 902 the image to be used as "moving" during optimization. The 903 dimensions of the static (S, R, C) and moving (S', R', C') images 904 do not need to be the same. 905 transform : instance of Transform 906 the transformation with respect to whose parameters the gradient 907 must be computed 908 params0 : array, shape (n,) 909 parameters from which to start the optimization. If None, the 910 optimization will start at the identity transform. n is the 911 number of parameters of the specified transformation. 912 static_grid2world : array, shape (dim+1, dim+1) 913 the voxel-to-space transformation associated with the static image 914 moving_grid2world : array, shape (dim+1, dim+1) 915 the voxel-to-space transformation associated with the moving image 916 starting_affine : string, or matrix, or None 917 If string: 918 'mass': align centers of gravity 919 'voxel-origin': align physical coordinates of voxel (0,0,0) 920 'centers': align physical coordinates of central voxels 921 If matrix: 922 array, shape (dim+1, dim+1) 923 If None: 924 Start from identity 925 926 """ 927 self.dim = len(static.shape) 928 self.transform = transform 929 n = transform.get_number_of_parameters() 930 self.nparams = n 931 932 if params0 is None: 933 params0 = self.transform.get_identity_parameters() 934 self.params0 = params0 935 if starting_affine is None: 936 self.starting_affine = np.eye(self.dim + 1) 937 elif isinstance(starting_affine, str): 938 if starting_affine == 'mass': 939 affine_map = transform_centers_of_mass(static, 940 static_grid2world, 941 moving, 942 moving_grid2world) 943 self.starting_affine = affine_map.affine 944 elif starting_affine == 'voxel-origin': 945 affine_map = transform_origins(static, static_grid2world, 946 moving, moving_grid2world) 947 self.starting_affine = affine_map.affine 948 elif starting_affine == 'centers': 949 affine_map = transform_geometric_centers(static, 950 static_grid2world, 951 moving, 952 moving_grid2world) 953 self.starting_affine = affine_map.affine 954 else: 955 raise ValueError('Invalid starting_affine strategy') 956 elif (isinstance(starting_affine, np.ndarray) and 957 starting_affine.shape >= (self.dim, self.dim + 1)): 958 self.starting_affine = starting_affine 959 else: 960 raise ValueError('Invalid starting_affine matrix') 961 # Extract information from affine matrices to create the scale space 962 static_direction, static_spacing = \ 963 get_direction_and_spacings(static_grid2world, self.dim) 964 moving_direction, moving_spacing = \ 965 get_direction_and_spacings(moving_grid2world, self.dim) 966 967 static = ((static.astype(np.float64) - static.min()) / 968 (static.max() - static.min())) 969 moving = ((moving.astype(np.float64) - moving.min()) / 970 (moving.max() - moving.min())) 971 972 # Build the scale space of the input images 973 if self.use_isotropic: 974 self.moving_ss = IsotropicScaleSpace(moving, self.factors, 975 self.sigmas, 976 moving_grid2world, 977 moving_spacing, False) 978 979 self.static_ss = IsotropicScaleSpace(static, self.factors, 980 self.sigmas, 981 static_grid2world, 982 static_spacing, False) 983 else: 984 self.moving_ss = ScaleSpace(moving, self.levels, moving_grid2world, 985 moving_spacing, self.ss_sigma_factor, 986 False) 987 988 self.static_ss = ScaleSpace(static, self.levels, static_grid2world, 989 static_spacing, self.ss_sigma_factor, 990 False) 991 992 def optimize(self, static, moving, transform, params0, 993 static_grid2world=None, moving_grid2world=None, 994 starting_affine=None, ret_metric=False): 995 r""" Start the optimization process. 996 997 Parameters 998 ---------- 999 static : 2D or 3D array 1000 the image to be used as reference during optimization. 1001 moving : 2D or 3D array 1002 the image to be used as "moving" during optimization. It is 1003 necessary to pre-align the moving image to ensure its domain 1004 lies inside the domain of the deformation fields. This is assumed 1005 to be accomplished by "pre-aligning" the moving image towards the 1006 static using an affine transformation given by the 1007 'starting_affine' matrix 1008 transform : instance of Transform 1009 the transformation with respect to whose parameters the gradient 1010 must be computed 1011 params0 : array, shape (n,) 1012 parameters from which to start the optimization. If None, the 1013 optimization will start at the identity transform. n is the 1014 number of parameters of the specified transformation. 1015 static_grid2world : array, shape (dim+1, dim+1), optional 1016 the voxel-to-space transformation associated with the static 1017 image. The default is None, implying the transform is the 1018 identity. 1019 moving_grid2world : array, shape (dim+1, dim+1), optional 1020 the voxel-to-space transformation associated with the moving 1021 image. The default is None, implying the transform is the 1022 identity. 1023 starting_affine : string, or matrix, or None, optional 1024 If string: 1025 'mass': align centers of gravity 1026 'voxel-origin': align physical coordinates of voxel (0,0,0) 1027 'centers': align physical coordinates of central voxels 1028 If matrix: 1029 array, shape (dim+1, dim+1). 1030 If None: 1031 Start from identity. 1032 The default is None. 1033 ret_metric : boolean, optional 1034 if True, it returns the parameters for measuring the 1035 similarity between the images (default 'False'). 1036 The metric containing optimal parameters and 1037 the distance between the images. 1038 1039 Returns 1040 ------- 1041 affine_map : instance of AffineMap 1042 the affine resulting affine transformation 1043 xopt : optimal parameters 1044 the optimal parameters (translation, rotation shear etc.) 1045 fopt : Similarity metric 1046 the value of the function at the optimal parameters. 1047 1048 """ 1049 self._init_optimizer(static, moving, transform, params0, 1050 static_grid2world, moving_grid2world, 1051 starting_affine) 1052 del starting_affine # Now we must refer to self.starting_affine 1053 1054 # Multi-resolution iterations 1055 original_static_shape = self.static_ss.get_image(0).shape 1056 original_static_grid2world = self.static_ss.get_affine(0) 1057 original_moving_shape = self.moving_ss.get_image(0).shape 1058 original_moving_grid2world = self.moving_ss.get_affine(0) 1059 affine_map = AffineMap(None, 1060 original_static_shape, 1061 original_static_grid2world, 1062 original_moving_shape, 1063 original_moving_grid2world) 1064 1065 for level in range(self.levels - 1, -1, -1): 1066 self.current_level = level 1067 max_iter = self.level_iters[-1 - level] 1068 if self.verbosity >= VerbosityLevels.STATUS: 1069 print('Optimizing level %d [max iter: %d]' % (level, max_iter)) 1070 1071 # Resample the smooth static image to the shape of this level 1072 smooth_static = self.static_ss.get_image(level) 1073 current_static_shape = self.static_ss.get_domain_shape(level) 1074 current_static_grid2world = self.static_ss.get_affine(level) 1075 current_affine_map = AffineMap(None, 1076 current_static_shape, 1077 current_static_grid2world, 1078 original_static_shape, 1079 original_static_grid2world) 1080 current_static = current_affine_map.transform(smooth_static) 1081 1082 # The moving image is full resolution 1083 current_moving_grid2world = original_moving_grid2world 1084 1085 current_moving = self.moving_ss.get_image(level) 1086 # Prepare the metric for iterations at this resolution 1087 self.metric.setup(transform, current_static, current_moving, 1088 current_static_grid2world, 1089 current_moving_grid2world, self.starting_affine) 1090 1091 # Optimize this level 1092 if self.options is None: 1093 self.options = {'gtol': 1e-4, 1094 'disp': False} 1095 1096 if self.method == 'L-BFGS-B': 1097 self.options['maxfun'] = max_iter 1098 else: 1099 self.options['maxiter'] = max_iter 1100 1101 opt = Optimizer(self.metric.distance_and_gradient, 1102 self.params0, 1103 method=self.method, jac=True, 1104 options=self.options) 1105 params = opt.xopt 1106 1107 # Update starting_affine matrix with optimal parameters 1108 T = self.transform.param_to_matrix(params) 1109 self.starting_affine = T.dot(self.starting_affine) 1110 1111 # Start next iteration at identity 1112 self.params0 = self.transform.get_identity_parameters() 1113 1114 affine_map.set_affine(self.starting_affine) 1115 if ret_metric: 1116 return affine_map, opt.xopt, opt.fopt 1117 return affine_map 1118 1119 1120def transform_centers_of_mass(static, static_grid2world, 1121 moving, moving_grid2world): 1122 r""" Transformation to align the center of mass of the input images. 1123 1124 Parameters 1125 ---------- 1126 static : array, shape (S, R, C) 1127 static image 1128 static_grid2world : array, shape (dim+1, dim+1) 1129 the voxel-to-space transformation of the static image 1130 moving : array, shape (S, R, C) 1131 moving image 1132 moving_grid2world : array, shape (dim+1, dim+1) 1133 the voxel-to-space transformation of the moving image 1134 1135 Returns 1136 ------- 1137 affine_map : instance of AffineMap 1138 the affine transformation (translation only, in this case) aligning 1139 the center of mass of the moving image towards the one of the static 1140 image 1141 1142 """ 1143 dim = len(static.shape) 1144 if static_grid2world is None: 1145 static_grid2world = np.eye(dim + 1) 1146 if moving_grid2world is None: 1147 moving_grid2world = np.eye(dim + 1) 1148 c_static = ndimage.measurements.center_of_mass(np.array(static)) 1149 c_static = static_grid2world.dot(c_static + (1,)) 1150 c_moving = ndimage.measurements.center_of_mass(np.array(moving)) 1151 c_moving = moving_grid2world.dot(c_moving + (1,)) 1152 transform = np.eye(dim + 1) 1153 transform[:dim, dim] = (c_moving - c_static)[:dim] 1154 affine_map = AffineMap(transform, 1155 static.shape, static_grid2world, 1156 moving.shape, moving_grid2world) 1157 return affine_map 1158 1159 1160def transform_geometric_centers(static, static_grid2world, 1161 moving, moving_grid2world): 1162 r""" Transformation to align the geometric center of the input images. 1163 1164 With "geometric center" of a volume we mean the physical coordinates of 1165 its central voxel 1166 1167 Parameters 1168 ---------- 1169 static : array, shape (S, R, C) 1170 static image 1171 static_grid2world : array, shape (dim+1, dim+1) 1172 the voxel-to-space transformation of the static image 1173 moving : array, shape (S, R, C) 1174 moving image 1175 moving_grid2world : array, shape (dim+1, dim+1) 1176 the voxel-to-space transformation of the moving image 1177 1178 Returns 1179 ------- 1180 affine_map : instance of AffineMap 1181 the affine transformation (translation only, in this case) aligning 1182 the geometric center of the moving image towards the one of the static 1183 image 1184 1185 """ 1186 dim = len(static.shape) 1187 if static_grid2world is None: 1188 static_grid2world = np.eye(dim + 1) 1189 if moving_grid2world is None: 1190 moving_grid2world = np.eye(dim + 1) 1191 c_static = tuple((np.array(static.shape, dtype=np.float64)) * 0.5) 1192 c_static = static_grid2world.dot(c_static + (1,)) 1193 c_moving = tuple((np.array(moving.shape, dtype=np.float64)) * 0.5) 1194 c_moving = moving_grid2world.dot(c_moving + (1,)) 1195 transform = np.eye(dim + 1) 1196 transform[:dim, dim] = (c_moving - c_static)[:dim] 1197 affine_map = AffineMap(transform, 1198 static.shape, static_grid2world, 1199 moving.shape, moving_grid2world) 1200 return affine_map 1201 1202 1203def transform_origins(static, static_grid2world, 1204 moving, moving_grid2world): 1205 r""" Transformation to align the origins of the input images. 1206 1207 With "origin" of a volume we mean the physical coordinates of 1208 voxel (0,0,0) 1209 1210 Parameters 1211 ---------- 1212 static : array, shape (S, R, C) 1213 static image 1214 static_grid2world : array, shape (dim+1, dim+1) 1215 the voxel-to-space transformation of the static image 1216 moving : array, shape (S, R, C) 1217 moving image 1218 moving_grid2world : array, shape (dim+1, dim+1) 1219 the voxel-to-space transformation of the moving image 1220 1221 Returns 1222 ------- 1223 affine_map : instance of AffineMap 1224 the affine transformation (translation only, in this case) aligning 1225 the origin of the moving image towards the one of the static 1226 image 1227 1228 """ 1229 dim = len(static.shape) 1230 if static_grid2world is None: 1231 static_grid2world = np.eye(dim + 1) 1232 if moving_grid2world is None: 1233 moving_grid2world = np.eye(dim + 1) 1234 c_static = static_grid2world[:dim, dim] 1235 c_moving = moving_grid2world[:dim, dim] 1236 transform = np.eye(dim + 1) 1237 transform[:dim, dim] = (c_moving - c_static)[:dim] 1238 affine_map = AffineMap(transform, 1239 static.shape, static_grid2world, 1240 moving.shape, moving_grid2world) 1241 return affine_map 1242