1#!python 2#cython: boundscheck=False 3#cython: wraparound=False 4#cython: cdivision=True 5 6 7import numpy as np 8cimport cython 9cimport numpy as cnp 10from dipy.align.fused_types cimport floating, number 11cdef extern from "dpy_math.h" nogil: 12 int dpy_isinf(double) 13 double floor(double) 14 15cdef inline int ifloor(double x) nogil: 16 return int(floor(x)) 17 18def quantize_positive_2d(floating[:, :] v, int num_levels): 19 r"""Quantizes a 2D image to num_levels quantization levels 20 21 Quantizes the input image at num_levels intensity levels considering <=0 22 as a special value. Those input pixels <=0, and only those, will be 23 assigned a quantization level of 0. The positive values are divided into 24 the remaining num_levels-1 uniform quantization levels. 25 26 The following are undefined, and raise a ValueError: 27 * Quantizing at zero levels because at least one level must be assigned 28 * Quantizing at one level because positive values should be assigned a 29 level different from the secial level 0 (at least 2 levels are needed) 30 31 Parameters 32 ---------- 33 v : array, shape (R, C) 34 the image to be quantized 35 num_levels : int 36 the number of levels 37 38 Returns 39 ------- 40 out : array, shape (R, C), same shape as v 41 the quantized image 42 levels: array, shape (num_levels,) 43 the quantization values: levels[0]=0, and levels[i] is the mid-point 44 of the interval of intensities that are assigned to quantization 45 level i, i=1, ..., num_levels-1. 46 hist: array, shape (num_levels,) 47 histogram: the number of pixels that were assigned to each quantization 48 level 49 """ 50 ftype = np.asarray(v).dtype 51 cdef: 52 cnp.npy_intp nrows = v.shape[0] 53 cnp.npy_intp ncols = v.shape[1] 54 cnp.npy_intp npix = nrows * ncols 55 cnp.npy_intp i, j, l 56 double epsilon, delta 57 double min_val = -1 58 double max_val = -1 59 cnp.npy_int32[:] hist = np.zeros(shape=(num_levels,), dtype=np.int32) 60 cnp.npy_int32[:, :] out = np.zeros(shape=(nrows, ncols,), dtype=np.int32) 61 floating[:] levels = np.zeros(shape=(num_levels,), dtype=ftype) 62 63 #Quantizing at zero levels is undefined 64 #Quantizing at one level is not supported because we want to make sure the 65 #maximum level in the quantization is never greater than num_levels-1 66 if(num_levels < 2): 67 raise ValueError('Quantization levels must be at least 2') 68 if (num_levels >= 2**31): 69 raise ValueError('Quantization levels must be < 2**31') 70 71 num_levels -= 1 # zero is one of the levels 72 73 with nogil: 74 75 for i in range(nrows): 76 for j in range(ncols): 77 if(v[i, j] > 0): 78 if((min_val < 0) or (v[i, j] < min_val)): 79 min_val = v[i, j] 80 if(v[i, j] > max_val): 81 max_val = v[i, j] 82 epsilon = 1e-8 83 delta = (max_val - min_val + epsilon) / num_levels 84 # notice that we decreased num_levels, so levels[0..num_levels] are well 85 # defined 86 if((num_levels < 2) or (delta < epsilon)): 87 for i in range(nrows): 88 for j in range(ncols): 89 if(v[i, j] > 0): 90 out[i, j] = 1 91 else: 92 out[i, j] = 0 93 hist[0] += 1 94 levels[0] = 0 95 levels[1] = 0.5 * (min_val + max_val) 96 hist[1] = npix - hist[0] 97 with gil: 98 return out, levels, hist 99 100 levels[0] = 0 101 levels[1] = min_val + delta * 0.5 102 for i in range(2, 1 + num_levels): 103 levels[i] = levels[i - 1] + delta 104 for i in range(nrows): 105 for j in range(ncols): 106 if(v[i, j] > 0): 107 l = ifloor((v[i, j] - min_val) / delta) 108 out[i, j] = l + 1 109 hist[l + 1] += 1 110 else: 111 out[i, j] = 0 112 hist[0] += 1 113 114 return np.asarray(out), np.array(levels), np.array(hist) 115 116 117def quantize_positive_3d(floating[:, :, :] v, int num_levels): 118 r"""Quantizes a 3D volume to num_levels quantization levels 119 120 Quantizes the input volume at num_levels intensity levels considering <=0 121 as a special value. Those input voxels <=0, and only those, will be 122 assigned a quantization level of 0. The positive values are divided into 123 the remaining num_levels-1 uniform quantization levels. 124 125 The following are undefined, and raise a ValueError: 126 * Quantizing at zero levels because at least one level must be assigned 127 * Quantizing at one level because positive values should be assigned a 128 level different from the secial level 0 (at least 2 levels are needed) 129 130 Parameters 131 ---------- 132 v : array, shape (S, R, C) 133 the volume to be quantized 134 num_levels : int 135 the number of levels 136 137 Returns 138 ------- 139 out : array, shape (S, R, C), same shape as v 140 the quantized volume 141 levels: array, shape (num_levels,) 142 the quantization values: levels[0]=0, and levels[i] is the mid-point 143 of the interval of intensities that are assigned to quantization 144 level i, i=1, ..., num_levels-1. 145 hist: array, shape (num_levels,) 146 histogram: the number of voxels that were assigned to each quantization 147 level 148 """ 149 ftype = np.asarray(v).dtype 150 cdef: 151 cnp.npy_intp nslices = v.shape[0] 152 cnp.npy_intp nrows = v.shape[1] 153 cnp.npy_intp ncols = v.shape[2] 154 cnp.npy_intp nvox = nrows * ncols * nslices 155 cnp.npy_intp i, j, k, l 156 double epsilon, delta 157 double min_val = -1 158 double max_val = -1 159 int[:] hist = np.zeros(shape=(num_levels,), dtype=np.int32) 160 int[:, :, :] out = np.zeros(shape=(nslices, nrows, ncols), 161 dtype=np.int32) 162 floating[:] levels = np.zeros(shape=(num_levels,), dtype=ftype) 163 164 #Quantizing at zero levels is undefined 165 #Quantizing at one level is not supported because we want to make sure the 166 #maximum level in the quantization is never greater than num_levels-1 167 if(num_levels < 2): 168 raise ValueError('Quantization levels must be at least 2') 169 170 num_levels -= 1 # zero is one of the levels 171 172 with nogil: 173 174 for k in range(nslices): 175 for i in range(nrows): 176 for j in range(ncols): 177 if(v[k, i, j] > 0): 178 if((min_val < 0) or (v[k, i, j] < min_val)): 179 min_val = v[k, i, j] 180 if(v[k, i, j] > max_val): 181 max_val = v[k, i, j] 182 epsilon = 1e-8 183 delta = (max_val - min_val + epsilon) / num_levels 184 # notice that we decreased num_levels, so levels[0..num_levels] are well 185 # defined 186 if((num_levels < 2) or (delta < epsilon)): 187 for k in range(nslices): 188 for i in range(nrows): 189 for j in range(ncols): 190 if(v[k, i, j] > 0): 191 out[k, i, j] = 1 192 else: 193 out[k, i, j] = 0 194 hist[0] += 1 195 levels[0] = 0 196 levels[1] = 0.5 * (min_val + max_val) 197 hist[1] = nvox - hist[0] 198 with gil: 199 return out, levels, hist 200 levels[0] = 0 201 levels[1] = min_val + delta * 0.5 202 for i in range(2, 1 + num_levels): 203 levels[i] = levels[i - 1] + delta 204 for k in range(nslices): 205 for i in range(nrows): 206 for j in range(ncols): 207 if(v[k, i, j] > 0): 208 l = ifloor((v[k, i, j] - min_val) / delta) 209 out[k, i, j] = l + 1 210 hist[l + 1] += 1 211 else: 212 out[k, i, j] = 0 213 hist[0] += 1 214 return np.asarray(out), np.asarray(levels), np.asarray(hist) 215 216 217def compute_masked_class_stats_2d(int[:, :] mask, floating[:, :] v, 218 int num_labels, int[:, :] labels): 219 r"""Computes the mean and std. for each quantization level. 220 221 Computes the mean and standard deviation of the intensities in 'v' for 222 each corresponding label in 'labels'. In other words, for each label 223 L, it computes the mean and standard deviation of the intensities in 'v' 224 at pixels whose label in 'labels' is L. This is used by the EM metric 225 to compute statistics for each hidden variable represented by the labels. 226 227 Parameters 228 ---------- 229 mask : array, shape (R, C) 230 the mask of pixels that will be taken into account for computing the 231 statistics. All zero pixels in mask will be ignored 232 v : array, shape (R, C) 233 the image which the statistics will be computed from 234 num_labels : int 235 the number of different labels in 'labels' (equal to the 236 number of hidden variables in the EM metric) 237 labels : array, shape (R, C) 238 the label assigned to each pixel 239 240 Returns 241 ------- 242 means : array, shape (num_labels,) 243 means[i], 0<=i<num_labels will be the mean intensity in v of all 244 voxels labeled i, or 0 if no voxels are labeled i 245 variances : array, shape (num_labels,) 246 variances[i], 0<=i<num_labels will be the standard deviation of the 247 intensities in v of all voxels labeled i, or infinite if less than 2 248 voxels are labeled i. 249 """ 250 ftype=np.asarray(v).dtype 251 cdef: 252 cnp.npy_intp nrows = v.shape[0] 253 cnp.npy_intp ncols = v.shape[1] 254 cnp.npy_intp i, j 255 double INF64 = np.inf 256 int[:] counts = np.zeros(shape=(num_labels,), dtype=np.int32) 257 floating diff 258 floating[:] means = np.zeros(shape=(num_labels,), dtype=ftype) 259 floating[:] variances = np.zeros(shape=(num_labels, ), dtype=ftype) 260 261 with nogil: 262 for i in range(nrows): 263 for j in range(ncols): 264 if(mask[i, j] != 0): 265 means[labels[i, j]] += v[i, j] 266 counts[labels[i, j]] += 1 267 for i in range(num_labels): 268 if(counts[i] > 0): 269 means[i] /= counts[i] 270 for i in range(nrows): 271 for j in range(ncols): 272 if(mask[i, j] != 0): 273 diff = v[i, j] - means[labels[i, j]] 274 variances[labels[i, j]] += diff ** 2 275 276 for i in range(num_labels): 277 if(counts[i] > 1): 278 variances[i] /= counts[i] 279 else: 280 variances[i] = INF64 281 return np.asarray(means), np.asarray(variances) 282 283 284def compute_masked_class_stats_3d(int[:, :, :] mask, floating[:, :, :] v, 285 int num_labels, int[:, :, :] labels): 286 r"""Computes the mean and std. for each quantization level. 287 288 Computes the mean and standard deviation of the intensities in 'v' for 289 each corresponding label in 'labels'. In other words, for each label 290 L, it computes the mean and standard deviation of the intensities in 'v' 291 at voxels whose label in 'labels' is L. This is used by the EM metric 292 to compute statistics for each hidden variable represented by the labels. 293 294 Parameters 295 ---------- 296 mask : array, shape (S, R, C) 297 the mask of voxels that will be taken into account for computing the 298 statistics. All zero voxels in mask will be ignored 299 v : array, shape (S, R, C) 300 the volume which the statistics will be computed from 301 num_labels : int 302 the number of different labels in 'labels' (equal to the 303 number of hidden variables in the EM metric) 304 labels : array, shape (S, R, C) 305 the label assigned to each pixel 306 307 Returns 308 ------- 309 means : array, shape (num_labels,) 310 means[i], 0<=i<num_labels will be the mean intensity in v of all 311 voxels labeled i, or 0 if no voxels are labeled i 312 variances : array, shape (num_labels,) 313 variances[i], 0<=i<num_labels will be the standard deviation of the 314 intensities in v of all voxels labeled i, or infinite if less than 2 315 voxels are labeled i. 316 """ 317 ftype=np.asarray(v).dtype 318 cdef: 319 cnp.npy_intp nslices = v.shape[0] 320 cnp.npy_intp nrows = v.shape[1] 321 cnp.npy_intp ncols = v.shape[2] 322 cnp.npy_intp i, j, k 323 double INF64 = np.inf 324 floating diff 325 int[:] counts = np.zeros(shape=(num_labels,), dtype=np.int32) 326 floating[:] means = np.zeros(shape=(num_labels,), dtype=ftype) 327 floating[:] variances = np.zeros(shape=(num_labels, ), dtype=ftype) 328 329 with nogil: 330 for k in range(nslices): 331 for i in range(nrows): 332 for j in range(ncols): 333 if(mask[k, i, j] != 0): 334 means[labels[k, i, j]] += v[k, i, j] 335 counts[labels[k, i, j]] += 1 336 for i in range(num_labels): 337 if(counts[i] > 0): 338 means[i] /= counts[i] 339 for k in range(nslices): 340 for i in range(nrows): 341 for j in range(ncols): 342 if(mask[k, i, j] != 0): 343 diff = means[labels[k, i, j]] - v[k, i, j] 344 variances[labels[k, i, j]] += diff ** 2 345 for i in range(num_labels): 346 if(counts[i] > 1): 347 variances[i] /= counts[i] 348 else: 349 variances[i] = INF64 350 return np.asarray(means), np.asarray(variances) 351 352@cython.boundscheck(False) 353@cython.wraparound(False) 354@cython.cdivision(True) 355def compute_em_demons_step_2d(floating[:,:] delta_field, 356 floating[:,:] sigma_sq_field, 357 floating[:,:,:] gradient_moving, 358 double sigma_sq_x, 359 floating[:,:,:] out): 360 r"""Demons step for EM metric in 2D 361 362 Computes the demons step [Vercauteren09] for SSD-driven registration 363 ( eq. 4 in [Vercauteren09] ) using the EM algorithm [Arce14] to handle 364 multi-modality images. 365 366 In this case, $\sigma_i$ in eq. 4 of [Vercauteren] is estimated using the EM 367 algorithm, while in the original version of diffeomorphic demons it is 368 estimated by the difference between the image values at each pixel. 369 370 Parameters 371 ---------- 372 delta_field : array, shape (R, C) 373 contains, at each pixel, the difference between the moving image (warped 374 under the current deformation s(. , .) ) J and the static image I: 375 delta_field[i,j] = J(s(i,j)) - I(i,j). The order is important, changing 376 to delta_field[i,j] = I(i,j) - J(s(i,j)) yields the backward demons step 377 warping the static image towards the moving, which may not be the 378 intended behavior unless the 'gradient_moving' passed corresponds to 379 the gradient of the static image 380 sigma_sq_field : array, shape (R, C) 381 contains, at each pixel (i, j), the estimated variance (not std) of the 382 hidden variable associated to the intensity at static[i,j] (which must 383 have been previously quantized) 384 gradient_moving : array, shape (R, C, 2) 385 the gradient of the moving image 386 sigma_sq_x : float 387 parameter controlling the amount of regularization. It corresponds to 388 $\sigma_x^2$ in algorithm 1 of Vercauteren et al.[2] 389 out : array, shape (R, C, 2) 390 the resulting demons step will be written to this array 391 392 Returns 393 ------- 394 demons_step : array, shape (R, C, 2) 395 the demons step to be applied for updating the current displacement 396 field 397 energy : float 398 the current em energy (before applying the returned demons_step) 399 400 References 401 ---------- 402 [Arce14] Arce-santana, E., Campos-delgado, D. U., & Vigueras-g, F. (2014). 403 Non-rigid Multimodal Image Registration Based on the 404 Expectation-Maximization Algorithm, (168140), 36-47. 405 406 [Vercauteren09] Vercauteren, T., Pennec, X., Perchant, A., & Ayache, N. 407 (2009). Diffeomorphic demons: efficient non-parametric 408 image registration. NeuroImage, 45(1 Suppl), S61-72. 409 doi:10.1016/j.neuroimage.2008.10.040 410 """ 411 cdef: 412 cnp.npy_intp nr = delta_field.shape[0] 413 cnp.npy_intp nc = delta_field.shape[1] 414 cnp.npy_intp i, j 415 double delta, sigma_sq_i, nrm2, energy, den, prod 416 417 if out is None: 418 out = np.zeros((nr, nc, 2), dtype=np.asarray(delta_field).dtype) 419 420 with nogil: 421 422 energy = 0 423 for i in range(nr): 424 for j in range(nc): 425 sigma_sq_i = sigma_sq_field[i,j] 426 delta = delta_field[i,j] 427 energy += (delta**2) 428 if dpy_isinf(sigma_sq_i) != 0: 429 out[i, j, 0], out[i, j, 1] = 0, 0 430 else: 431 nrm2 = (gradient_moving[i, j, 0]**2 + 432 gradient_moving[i, j, 1]**2) 433 if(sigma_sq_i == 0): 434 if nrm2 == 0: 435 out[i, j, 0], out[i, j, 1] = 0, 0 436 else: 437 out[i, j, 0] = (delta * 438 gradient_moving[i, j, 0] / nrm2) 439 out[i, j, 1] = (delta * 440 gradient_moving[i, j, 1] / nrm2) 441 else: 442 den = (sigma_sq_x * nrm2 + sigma_sq_i) 443 prod = sigma_sq_x * delta 444 out[i, j, 0] = prod * gradient_moving[i, j, 0] / den 445 out[i, j, 1] = prod * gradient_moving[i, j, 1] / den 446 return np.asarray(out), energy 447 448@cython.boundscheck(False) 449@cython.wraparound(False) 450@cython.cdivision(True) 451def compute_em_demons_step_3d(floating[:,:,:] delta_field, 452 floating[:,:,:] sigma_sq_field, 453 floating[:,:,:,:] gradient_moving, 454 double sigma_sq_x, 455 floating[:,:,:,:] out): 456 r"""Demons step for EM metric in 3D 457 458 Computes the demons step [Vercauteren09] for SSD-driven registration 459 ( eq. 4 in [Vercauteren09] ) using the EM algorithm [Arce14] to handle 460 multi-modality images. 461 462 In this case, $\sigma_i$ in eq. 4 of [Vercauteren09] is estimated using 463 the EM algorithm, while in the original version of diffeomorphic demons 464 it is estimated by the difference between the image values at each pixel. 465 466 Parameters 467 ---------- 468 delta_field : array, shape (S, R, C) 469 contains, at each pixel, the difference between the moving image (warped 470 under the current deformation s ) J and the static image I: 471 delta_field[k,i,j] = J(s(k,i,j)) - I(k,i,j). The order is important, 472 changing to delta_field[k,i,j] = I(k,i,j) - J(s(k,i,j)) yields the 473 backward demons step warping the static image towards the moving, which 474 may not be the intended behavior unless the 'gradient_moving' passed 475 corresponds to the gradient of the static image 476 sigma_sq_field : array, shape (S, R, C) 477 contains, at each pixel (k, i, j), the estimated variance (not std) of 478 the hidden variable associated to the intensity at static[k,i,j] (which 479 must have been previously quantized) 480 gradient_moving : array, shape (S, R, C, 2) 481 the gradient of the moving image 482 sigma_sq_x : float 483 parameter controlling the amount of regularization. It corresponds to 484 $\sigma_x^2$ in algorithm 1 of Vercauteren et al.[2]. 485 out : array, shape (S, R, C, 2) 486 the resulting demons step will be written to this array 487 488 Returns 489 ------- 490 demons_step : array, shape (S, R, C, 3) 491 the demons step to be applied for updating the current displacement 492 field 493 energy : float 494 the current em energy (before applying the returned demons_step) 495 496 References 497 ---------- 498 [Arce14] Arce-santana, E., Campos-delgado, D. U., & Vigueras-g, F. (2014). 499 Non-rigid Multimodal Image Registration Based on the 500 Expectation-Maximization Algorithm, (168140), 36-47. 501 502 [Vercauteren09] Vercauteren, T., Pennec, X., Perchant, A., & Ayache, N. 503 (2009). Diffeomorphic demons: efficient non-parametric 504 image registration. NeuroImage, 45(1 Suppl), S61-72. 505 doi:10.1016/j.neuroimage.2008.10.040 506 """ 507 cdef: 508 cnp.npy_intp ns = delta_field.shape[0] 509 cnp.npy_intp nr = delta_field.shape[1] 510 cnp.npy_intp nc = delta_field.shape[2] 511 cnp.npy_intp i, j, k 512 double delta, sigma_sq_i, nrm2, energy, den 513 514 if out is None: 515 out = np.zeros((ns, nr, nc, 3), dtype=np.asarray(delta_field).dtype) 516 517 with nogil: 518 519 energy = 0 520 for k in range(ns): 521 for i in range(nr): 522 for j in range(nc): 523 sigma_sq_i = sigma_sq_field[k,i,j] 524 delta = delta_field[k,i,j] 525 energy += (delta**2) 526 if dpy_isinf(sigma_sq_i) != 0: 527 out[k, i, j, 0] = 0 528 out[k, i, j, 1] = 0 529 out[k, i, j, 2] = 0 530 else: 531 nrm2 = (gradient_moving[k, i, j, 0]**2 + 532 gradient_moving[k, i, j, 1]**2 + 533 gradient_moving[k, i, j, 2]**2) 534 if(sigma_sq_i == 0): 535 if nrm2 == 0: 536 out[k, i, j, 0] = 0 537 out[k, i, j, 1] = 0 538 out[k, i, j, 2] = 0 539 else: 540 out[k, i, j, 0] = (delta * 541 gradient_moving[k, i, j, 0] / nrm2) 542 out[k, i, j, 1] = (delta * 543 gradient_moving[k, i, j, 1] / nrm2) 544 out[k, i, j, 2] = (delta * 545 gradient_moving[k, i, j, 2] / nrm2) 546 else: 547 den = (sigma_sq_x * nrm2 + sigma_sq_i) 548 out[k, i, j, 0] = (sigma_sq_x * delta * 549 gradient_moving[k, i, j, 0] / den) 550 out[k, i, j, 1] = (sigma_sq_x * delta * 551 gradient_moving[k, i, j, 1] / den) 552 out[k, i, j, 2] = (sigma_sq_x * delta * 553 gradient_moving[k, i, j, 2] / den) 554 return np.asarray(out), energy 555