1# distutils: language = c 2# cython: wraparound=False, cdivision=True, boundscheck=False 3 4import cython 5import numpy as np 6from libc.math cimport sqrt 7from libc.stdlib cimport malloc, free 8 9cimport numpy as cnp 10 11from dipy.tracking import Streamlines 12 13 14cdef extern from "dpy_math.h" nogil: 15 bint dpy_isnan(double x) 16 17 18cdef double c_length(Streamline streamline) nogil: 19 cdef: 20 cnp.npy_intp i 21 double out = 0.0 22 double dn, sum_dn_sqr 23 24 for i in range(1, streamline.shape[0]): 25 sum_dn_sqr = 0.0 26 for j in range(streamline.shape[1]): 27 dn = streamline[i, j] - streamline[i-1, j] 28 sum_dn_sqr += dn*dn 29 30 out += sqrt(sum_dn_sqr) 31 32 return out 33 34 35cdef void c_arclengths_from_arraysequence(Streamline points, 36 cnp.npy_intp[:] offsets, 37 cnp.npy_intp[:] lengths, 38 double[:] arclengths) nogil: 39 cdef: 40 cnp.npy_intp i, j, k 41 cnp.npy_intp offset 42 double dn, sum_dn_sqr 43 44 for i in range(offsets.shape[0]): 45 offset = offsets[i] 46 47 arclengths[i] = 0 48 for j in range(1, lengths[i]): 49 sum_dn_sqr = 0.0 50 for k in range(points.shape[1]): 51 dn = points[offset+j, k] - points[offset+j-1, k] 52 sum_dn_sqr += dn*dn 53 54 arclengths[i] += sqrt(sum_dn_sqr) 55 56 57def length(streamlines): 58 """ Euclidean length of streamlines 59 60 Length is in mm only if streamlines are expressed in world coordinates. 61 62 Parameters 63 ------------ 64 streamlines : ndarray or a list or :class:`dipy.tracking.Streamlines` 65 If ndarray, must have shape (N,3) where N is the number of points 66 of the streamline. 67 If list, each item must be ndarray shape (Ni,3) where Ni is the number 68 of points of streamline i. 69 If :class:`dipy.tracking.Streamlines`, its `common_shape` must be 3. 70 71 Returns 72 --------- 73 lengths : scalar or ndarray shape (N,) 74 If there is only one streamline, a scalar representing the length of the 75 streamline. 76 If there are several streamlines, ndarray containing the length of every 77 streamline. 78 79 Examples 80 ---------- 81 >>> from dipy.tracking.streamline import length 82 >>> import numpy as np 83 >>> streamline = np.array([[1, 1, 1], [2, 3, 4], [0, 0, 0]]) 84 >>> expected_length = np.sqrt([1+2**2+3**2, 2**2+3**2+4**2]).sum() 85 >>> length(streamline) == expected_length 86 True 87 >>> streamlines = [streamline, np.vstack([streamline, streamline[::-1]])] 88 >>> expected_lengths = [expected_length, 2*expected_length] 89 >>> lengths = [length(streamlines[0]), length(streamlines[1])] 90 >>> np.allclose(lengths, expected_lengths) 91 True 92 >>> length([]) 93 0.0 94 >>> length(np.array([[1, 2, 3]])) 95 0.0 96 97 """ 98 if isinstance(streamlines, Streamlines): 99 if len(streamlines) == 0: 100 return 0.0 101 102 arclengths = np.zeros(len(streamlines), dtype=np.float64) 103 104 if streamlines._data.dtype == np.float32: 105 c_arclengths_from_arraysequence[float2d]( 106 streamlines._data, 107 streamlines._offsets.astype(np.intp), 108 streamlines._lengths.astype(np.intp), 109 arclengths) 110 else: 111 c_arclengths_from_arraysequence[double2d]( 112 streamlines._data, 113 streamlines._offsets.astype(np.intp), 114 streamlines._lengths.astype(np.intp), 115 arclengths) 116 117 return arclengths 118 119 only_one_streamlines = False 120 if type(streamlines) is cnp.ndarray: 121 only_one_streamlines = True 122 streamlines = [streamlines] 123 124 if len(streamlines) == 0: 125 return 0.0 126 127 dtype = streamlines[0].dtype 128 for streamline in streamlines: 129 if streamline.dtype != dtype: 130 dtype = None 131 break 132 133 # Allocate memory for each streamline length. 134 streamlines_length = np.empty(len(streamlines), dtype=np.float64) 135 cdef cnp.npy_intp i 136 137 if dtype is None: 138 # List of streamlines having different dtypes 139 for i in range(len(streamlines)): 140 dtype = streamlines[i].dtype 141 # HACK: To avoid memleaks we have to recast with astype(dtype). 142 streamline = streamlines[i].astype(dtype) 143 if dtype != np.float32 and dtype != np.float64: 144 is_integer = dtype == np.int64 or dtype == np.uint64 145 dtype = np.float64 if is_integer else np.float32 146 streamline = streamlines[i].astype(dtype) 147 148 if dtype == np.float32: 149 streamlines_length[i] = c_length[float2d](streamline) 150 else: 151 streamlines_length[i] = c_length[double2d](streamline) 152 153 elif dtype == np.float32: 154 # All streamlines have composed of float32 points 155 for i in range(len(streamlines)): 156 # HACK: To avoid memleaks we have to recast with astype(dtype). 157 streamline = streamlines[i].astype(dtype) 158 streamlines_length[i] = c_length[float2d](streamline) 159 elif dtype == np.float64: 160 # All streamlines are composed of float64 points 161 for i in range(len(streamlines)): 162 # HACK: To avoid memleaks we have to recast with astype(dtype). 163 streamline = streamlines[i].astype(dtype) 164 streamlines_length[i] = c_length[double2d](streamline) 165 elif dtype == np.int64 or dtype == np.uint64: 166 # All streamlines are composed of int64 or uint64 points so convert 167 # them in float64 one at the time. 168 for i in range(len(streamlines)): 169 streamline = streamlines[i].astype(np.float64) 170 streamlines_length[i] = c_length[double2d](streamline) 171 else: 172 # All streamlines are composed of points with a dtype fitting in 173 # 32 bits so convert them in float32 one at the time. 174 for i in range(len(streamlines)): 175 streamline = streamlines[i].astype(np.float32) 176 streamlines_length[i] = c_length[float2d](streamline) 177 178 if only_one_streamlines: 179 return streamlines_length[0] 180 else: 181 return streamlines_length 182 183 184cdef void c_arclengths(Streamline streamline, double* out) nogil: 185 cdef cnp.npy_intp i = 0 186 cdef double dn 187 188 out[0] = 0.0 189 for i in range(1, streamline.shape[0]): 190 out[i] = 0.0 191 for j in range(streamline.shape[1]): 192 dn = streamline[i, j] - streamline[i-1, j] 193 out[i] += dn*dn 194 195 out[i] = out[i-1] + sqrt(out[i]) 196 197 198cdef void c_set_number_of_points(Streamline streamline, Streamline out) nogil: 199 cdef: 200 cnp.npy_intp N = streamline.shape[0] 201 cnp.npy_intp D = streamline.shape[1] 202 cnp.npy_intp new_N = out.shape[0] 203 double ratio, step, next_point, delta 204 cnp.npy_intp i, j, k, dim 205 206 # Get arclength at each point. 207 arclengths = <double*> malloc(streamline.shape[0] * sizeof(double)) 208 c_arclengths(streamline, arclengths) 209 210 step = arclengths[N-1] / (new_N-1) 211 212 next_point = 0.0 213 i = 0 214 j = 0 215 k = 0 216 217 while next_point < arclengths[N-1]: 218 if next_point == arclengths[k]: 219 for dim in range(D): 220 out[i, dim] = streamline[j, dim] 221 222 next_point += step 223 i += 1 224 j += 1 225 k += 1 226 elif next_point < arclengths[k]: 227 ratio = 1 - ((arclengths[k]-next_point) / 228 (arclengths[k]-arclengths[k-1])) 229 230 for dim in range(D): 231 delta = streamline[j, dim] - streamline[j-1, dim] 232 out[i, dim] = streamline[j-1, dim] + ratio * delta 233 234 next_point += step 235 i += 1 236 else: 237 j += 1 238 k += 1 239 240 # Last resampled point always the one from original streamline. 241 for dim in range(D): 242 out[new_N-1, dim] = streamline[N-1, dim] 243 244 free(arclengths) 245 246 247cdef void c_set_number_of_points_from_arraysequence(Streamline points, 248 cnp.npy_intp[:] offsets, 249 cnp.npy_intp[:] lengths, 250 long nb_points, 251 Streamline out) nogil: 252 cdef: 253 cnp.npy_intp i, j, k 254 cnp.npy_intp offset, length 255 cnp.npy_intp offset_out = 0 256 double dn, sum_dn_sqr 257 258 for i in range(offsets.shape[0]): 259 offset = offsets[i] 260 length = lengths[i] 261 262 c_set_number_of_points(points[offset:offset+length, :], 263 out[offset_out:offset_out+nb_points, :]) 264 265 offset_out += nb_points 266 267 268def set_number_of_points(streamlines, nb_points=3): 269 """ Change the number of points of streamlines 270 (either by downsampling or upsampling) 271 272 Change the number of points of streamlines in order to obtain 273 `nb_points`-1 segments of equal length. Points of streamlines will be 274 modified along the curve. 275 276 Parameters 277 ---------- 278 streamlines : ndarray or a list or :class:`dipy.tracking.Streamlines` 279 If ndarray, must have shape (N,3) where N is the number of points 280 of the streamline. 281 If list, each item must be ndarray shape (Ni,3) where Ni is the number 282 of points of streamline i. 283 If :class:`dipy.tracking.Streamlines`, its `common_shape` must be 3. 284 285 nb_points : int 286 integer representing number of points wanted along the curve. 287 288 Returns 289 ------- 290 new_streamlines : ndarray or a list or :class:`dipy.tracking.Streamlines` 291 Results of the downsampling or upsampling process. 292 293 Examples 294 -------- 295 >>> from dipy.tracking.streamline import set_number_of_points 296 >>> import numpy as np 297 298 One streamline, a semi-circle: 299 300 >>> theta = np.pi*np.linspace(0, 1, 100) 301 >>> x = np.cos(theta) 302 >>> y = np.sin(theta) 303 >>> z = 0 * x 304 >>> streamline = np.vstack((x, y, z)).T 305 >>> modified_streamline = set_number_of_points(streamline, 3) 306 >>> len(modified_streamline) 307 3 308 309 Multiple streamlines: 310 311 >>> streamlines = [streamline, streamline[::2]] 312 >>> new_streamlines = set_number_of_points(streamlines, 10) 313 >>> [len(s) for s in streamlines] 314 [100, 50] 315 >>> [len(s) for s in new_streamlines] 316 [10, 10] 317 318 """ 319 if isinstance(streamlines, Streamlines): 320 if len(streamlines) == 0: 321 return Streamlines() 322 323 nb_streamlines = len(streamlines) 324 dtype = streamlines._data.dtype 325 new_streamlines = Streamlines() 326 new_streamlines._data = np.zeros((nb_streamlines * nb_points, 3), 327 dtype=dtype) 328 new_streamlines._offsets = nb_points * np.arange(nb_streamlines, 329 dtype=np.intp) 330 new_streamlines._lengths = nb_points * np.ones(nb_streamlines, 331 dtype=np.intp) 332 333 if dtype == np.float32: 334 c_set_number_of_points_from_arraysequence[float2d]( 335 streamlines._data, streamlines._offsets.astype(np.intp), 336 streamlines._lengths.astype(np.intp), nb_points, 337 new_streamlines._data) 338 else: 339 c_set_number_of_points_from_arraysequence[double2d]( 340 streamlines._data, streamlines._offsets.astype(np.intp), 341 streamlines._lengths.astype(np.intp), nb_points, 342 new_streamlines._data) 343 344 return new_streamlines 345 346 only_one_streamlines = False 347 if type(streamlines) is cnp.ndarray: 348 only_one_streamlines = True 349 streamlines = [streamlines] 350 351 if len(streamlines) == 0: 352 return [] 353 354 if nb_points < 2: 355 raise ValueError("nb_points must be at least 2") 356 357 dtype = streamlines[0].dtype 358 for streamline in streamlines: 359 if streamline.dtype != dtype: 360 dtype = None 361 362 if len(streamline) < 2: 363 raise ValueError("All streamlines must have at least 2 points.") 364 365 # Allocate memory for each modified streamline 366 new_streamlines = [] 367 cdef cnp.npy_intp i 368 369 if dtype is None: 370 # List of streamlines having different dtypes 371 for i in range(len(streamlines)): 372 dtype = streamlines[i].dtype 373 # HACK: To avoid memleaks we have to recast with astype(dtype). 374 streamline = streamlines[i].astype(dtype) 375 if dtype != np.float32 and dtype != np.float64: 376 dtype = np.float32 377 if dtype == np.int64 or dtype == np.uint64: 378 dtype = np.float64 379 380 streamline = streamline.astype(dtype) 381 382 new_streamline = np.empty((nb_points, streamline.shape[1]), 383 dtype=dtype) 384 if dtype == np.float32: 385 c_set_number_of_points[float2d](streamline, new_streamline) 386 else: 387 c_set_number_of_points[double2d](streamline, new_streamline) 388 389 # HACK: To avoid memleaks we have to recast with astype(dtype). 390 new_streamlines.append(new_streamline.astype(dtype)) 391 392 elif dtype == np.float32: 393 # All streamlines have composed of float32 points 394 for i in range(len(streamlines)): 395 streamline = streamlines[i].astype(dtype) 396 modified_streamline = np.empty((nb_points, streamline.shape[1]), 397 dtype=streamline.dtype) 398 c_set_number_of_points[float2d](streamline, modified_streamline) 399 # HACK: To avoid memleaks we have to recast with astype(dtype). 400 new_streamlines.append(modified_streamline.astype(dtype)) 401 elif dtype == np.float64: 402 # All streamlines are composed of float64 points 403 for i in range(len(streamlines)): 404 streamline = streamlines[i].astype(dtype) 405 modified_streamline = np.empty((nb_points, streamline.shape[1]), 406 dtype=streamline.dtype) 407 c_set_number_of_points[double2d](streamline, modified_streamline) 408 # HACK: To avoid memleaks we have to recast with astype(dtype). 409 new_streamlines.append(modified_streamline.astype(dtype)) 410 elif dtype == np.int64 or dtype == np.uint64: 411 # All streamlines are composed of int64 or uint64 points so convert 412 # them in float64 one at the time 413 for i in range(len(streamlines)): 414 streamline = streamlines[i].astype(np.float64) 415 modified_streamline = np.empty((nb_points, streamline.shape[1]), 416 dtype=streamline.dtype) 417 c_set_number_of_points[double2d](streamline, modified_streamline) 418 # HACK: To avoid memleaks we've to recast with astype(np.float64). 419 new_streamlines.append(modified_streamline.astype(np.float64)) 420 else: 421 # All streamlines are composed of points with a dtype fitting in 422 # 32bits so convert them in float32 one at the time 423 for i in range(len(streamlines)): 424 streamline = streamlines[i].astype(np.float32) 425 modified_streamline = np.empty((nb_points, streamline.shape[1]), 426 dtype=streamline.dtype) 427 c_set_number_of_points[float2d](streamline, modified_streamline) 428 # HACK: To avoid memleaks we've to recast with astype(np.float32). 429 new_streamlines.append(modified_streamline.astype(np.float32)) 430 431 if only_one_streamlines: 432 return new_streamlines[0] 433 else: 434 return new_streamlines 435 436 437cdef double c_norm_of_cross_product(double bx, double by, double bz, 438 double cx, double cy, double cz) nogil: 439 """ Computes the norm of the cross-product in 3D. """ 440 cdef double ax, ay, az 441 ax = by*cz - bz*cy 442 ay = bz*cx - bx*cz 443 az = bx*cy - by*cx 444 return sqrt(ax*ax + ay*ay + az*az) 445 446 447cdef double c_dist_to_line(Streamline streamline, cnp.npy_intp prev, 448 cnp.npy_intp next, cnp.npy_intp curr) nogil: 449 """ Computes the shortest Euclidean distance between a point `curr` and 450 the line passing through `prev` and `next`. """ 451 452 cdef: 453 double dn, norm1, norm2 454 cnp.npy_intp D = streamline.shape[1] 455 456 # Compute cross product of next-prev and curr-next 457 norm1 = c_norm_of_cross_product(streamline[next, 0]-streamline[prev, 0], 458 streamline[next, 1]-streamline[prev, 1], 459 streamline[next, 2]-streamline[prev, 2], 460 streamline[curr, 0]-streamline[next, 0], 461 streamline[curr, 1]-streamline[next, 1], 462 streamline[curr, 2]-streamline[next, 2]) 463 464 # Norm of next-prev 465 norm2 = 0.0 466 for d in range(D): 467 dn = streamline[next, d]-streamline[prev, d] 468 norm2 += dn*dn 469 norm2 = sqrt(norm2) 470 471 return norm1 / norm2 472 473 474cdef double c_segment_length(Streamline streamline, 475 cnp.npy_intp start, cnp.npy_intp end) nogil: 476 """ Computes the length of the segment going from `start` to `end`. """ 477 cdef: 478 cnp.npy_intp D = streamline.shape[1] 479 cnp.npy_intp d 480 double segment_length = 0.0 481 double dn 482 483 for d in range(D): 484 dn = streamline[end, d] - streamline[start, d] 485 segment_length += dn*dn 486 487 return sqrt(segment_length) 488 489 490cdef cnp.npy_intp c_compress_streamline(Streamline streamline, Streamline out, 491 double tol_error, double max_segment_length) nogil: 492 """ Compresses a streamline (see function `compress_streamlines`).""" 493 cdef: 494 cnp.npy_intp N = streamline.shape[0] 495 cnp.npy_intp D = streamline.shape[1] 496 cnp.npy_intp nb_points = 0 497 cnp.npy_intp d, prev, next, curr 498 double segment_length 499 500 # Copy first point since it is always kept. 501 for d in range(D): 502 out[0, d] = streamline[0, d] 503 504 nb_points = 1 505 prev = 0 506 507 # Loop through the points of the streamline checking if we can use the 508 # linearized segment: next-prev. We start with next=2 (third points) since 509 # we already added point 0 and segment between the two firsts is linear. 510 for next in range(2, N): 511 # Euclidean distance between last added point and current point. 512 if c_segment_length(streamline, prev, next) > max_segment_length: 513 for d in range(D): 514 out[nb_points, d] = streamline[next-1, d] 515 516 nb_points += 1 517 prev = next-1 518 continue 519 520 # Check that each point is not offset by more than `tol_error` mm. 521 for curr in range(prev+1, next): 522 dist = c_dist_to_line(streamline, prev, next, curr) 523 524 if dpy_isnan(dist) or dist > tol_error: 525 for d in range(D): 526 out[nb_points, d] = streamline[next-1, d] 527 528 nb_points += 1 529 prev = next-1 530 break 531 532 # Copy last point since it is always kept. 533 for d in range(D): 534 out[nb_points, d] = streamline[N-1, d] 535 536 nb_points += 1 537 return nb_points 538 539 540def compress_streamlines(streamlines, tol_error=0.01, max_segment_length=10): 541 """ Compress streamlines by linearization as in [Presseau15]_. 542 543 The compression consists in merging consecutive segments that are 544 nearly collinear. The merging is achieved by removing the point the two 545 segments have in common. 546 547 The linearization process [Presseau15]_ ensures that every point being 548 removed are within a certain margin (in mm) of the resulting streamline. 549 Recommendations for setting this margin can be found in [Presseau15]_ 550 (in which they called it tolerance error). 551 552 The compression also ensures that two consecutive points won't be too far 553 from each other (precisely less or equal than `max_segment_length`mm). 554 This is a tradeoff to speed up the linearization process [Rheault15]_. A low 555 value will result in a faster linearization but low compression, whereas 556 a high value will result in a slower linearization but high compression. 557 558 Parameters 559 ---------- 560 streamlines : one or a list of array-like of shape (N,3) 561 Array representing x,y,z of N points in a streamline. 562 tol_error : float (optional) 563 Tolerance error in mm (default: 0.01). A rule of thumb is to set it 564 to 0.01mm for deterministic streamlines and 0.1mm for probabilitic 565 streamlines. 566 max_segment_length : float (optional) 567 Maximum length in mm of any given segment produced by the compression. 568 The default is 10mm. (In [Presseau15]_, they used a value of `np.inf`). 569 570 Returns 571 ------- 572 compressed_streamlines : one or a list of array-like 573 Results of the linearization process. 574 575 Examples 576 -------- 577 >>> from dipy.tracking.streamline import compress_streamlines 578 >>> import numpy as np 579 >>> # One streamline: a wiggling line 580 >>> rng = np.random.RandomState(42) 581 >>> streamline = np.linspace(0, 10, 100*3).reshape((100, 3)) 582 >>> streamline += 0.2 * rng.rand(100, 3) 583 >>> c_streamline = compress_streamlines(streamline, tol_error=0.2) 584 >>> len(streamline) 585 100 586 >>> len(c_streamline) 587 10 588 >>> # Multiple streamlines 589 >>> streamlines = [streamline, streamline[::2]] 590 >>> c_streamlines = compress_streamlines(streamlines, tol_error=0.2) 591 >>> [len(s) for s in streamlines] 592 [100, 50] 593 >>> [len(s) for s in c_streamlines] 594 [10, 7] 595 596 597 Notes 598 ----- 599 Be aware that compressed streamlines have variable step sizes. One needs to 600 be careful when computing streamlines-based metrics [Houde15]_. 601 602 References 603 ---------- 604 .. [Presseau15] Presseau C. et al., A new compression format for fiber 605 tracking datasets, NeuroImage, no 109, 73-83, 2015. 606 .. [Rheault15] Rheault F. et al., Real Time Interaction with Millions of 607 Streamlines, ISMRM, 2015. 608 .. [Houde15] Houde J.-C. et al. How to Avoid Biased Streamlines-Based 609 Metrics for Streamlines with Variable Step Sizes, ISMRM, 2015. 610 """ 611 only_one_streamlines = False 612 if type(streamlines) is cnp.ndarray: 613 only_one_streamlines = True 614 streamlines = [streamlines] 615 616 if len(streamlines) == 0: 617 return [] 618 619 compressed_streamlines = [] 620 cdef cnp.npy_intp i 621 for i in range(len(streamlines)): 622 dtype = streamlines[i].dtype 623 # HACK: To avoid memleaks we have to recast with astype(dtype). 624 streamline = streamlines[i].astype(dtype) 625 shape = streamline.shape 626 627 if dtype != np.float32 and dtype != np.float64: 628 dtype = np.float64 if dtype == np.int64 or dtype == np.uint64 else np.float32 629 streamline = streamline.astype(dtype) 630 631 if shape[0] <= 2: 632 compressed_streamlines.append(streamline.copy()) 633 continue 634 635 compressed_streamline = np.empty(shape, dtype) 636 637 if dtype == np.float32: 638 nb_points = c_compress_streamline[float2d](streamline, compressed_streamline, 639 tol_error, max_segment_length) 640 else: 641 nb_points = c_compress_streamline[double2d](streamline, compressed_streamline, 642 tol_error, max_segment_length) 643 644 compressed_streamline.resize((nb_points, streamline.shape[1])) 645 # HACK: To avoid memleaks we have to recast with astype(dtype). 646 compressed_streamlines.append(compressed_streamline.astype(dtype)) 647 648 if only_one_streamlines: 649 return compressed_streamlines[0] 650 else: 651 return compressed_streamlines 652