1#cython: cdivision=True 2#cython: boundscheck=False 3#cython: nonecheck=False 4#cython: wraparound=False 5 6import numpy as np 7from warnings import warn 8 9cimport numpy as cnp 10cnp.import_array() 11 12DTYPE = np.intp 13cdef DTYPE_t BG_NODE_NULL = -999 14 15cdef struct s_shpinfo 16 17ctypedef s_shpinfo shape_info 18ctypedef size_t (* fun_ravel)(size_t, size_t, size_t, shape_info *) nogil 19 20 21# For having stuff concerning background in one place 22ctypedef struct bginfo: 23 ## The value in the image (i.e. not the label!) that identifies 24 ## the background. 25 DTYPE_t background_val 26 DTYPE_t background_node 27 ## Identification of the background in the labelled image 28 DTYPE_t background_label 29 30 31cdef void get_bginfo(background_val, bginfo *ret) except *: 32 if background_val is None: 33 ret.background_val = 0 34 else: 35 ret.background_val = background_val 36 37 # The node -999 doesn't exist, it will get subsituted by a meaningful value 38 # upon the first background pixel occurrence 39 ret.background_node = BG_NODE_NULL 40 ret.background_label = 0 41 42 43# A pixel has neighbors that have already been scanned. 44# In the paper, the pixel is denoted by E and its neighbors: 45# 46# z=1 z=0 x 47# ----------------------> 48# | A B C F G H 49# | D E . I J K 50# | . . . L M N 51# | 52# y V 53# 54# D_ea represents offset of A from E etc. - see the definition of 55# get_shape_info 56cdef enum: 57 # the 0D neighbor 58 # D_ee, # We don't need D_ee 59 # the 1D neighbor 60 D_ed, 61 # 2D neighbors 62 D_ea, D_eb, D_ec, 63 # 3D neighbors 64 D_ef, D_eg, D_eh, D_ei, D_ej, D_ek, D_el, D_em, D_en, 65 D_COUNT 66 67 68# Structure for centralized access to shape data 69# Contains information related to the shape of the input array 70cdef struct s_shpinfo: 71 DTYPE_t x 72 DTYPE_t y 73 DTYPE_t z 74 75 # Number of elements 76 DTYPE_t numels 77 # Dimensions of of the input array 78 DTYPE_t ndim 79 80 # Offsets between elements recalculated to linear index increments 81 # DEX[D_ea] is offset between E and A (i.e. to the point to upper left) 82 # The name DEX is supposed to evoke DE., where . = A, B, C, D, F etc. 83 DTYPE_t DEX[D_COUNT] 84 85 # Function pointer to a function that recalculates multi-index to linear 86 # index. Heavily depends on dimensions of the input array. 87 fun_ravel ravel_index 88 89 90cdef void get_shape_info(inarr_shape, shape_info *res) except *: 91 """ 92 Precalculates all the needed data from the input array shape 93 and stores them in the shape_info struct. 94 """ 95 res.y = 1 96 res.z = 1 97 res.ravel_index = ravel_index2D 98 # A shape (3, 1, 4) would make the algorithm crash, but the corresponding 99 # good_shape (i.e. the array with axis swapped) (1, 3, 4) is OK. 100 # Having an axis length of 1 when an axis on the left is longer than 1 101 # (in this case, it has length of 3) is NOT ALLOWED. 102 good_shape = tuple(sorted(inarr_shape)) 103 104 res.ndim = len(inarr_shape) 105 if res.ndim == 1: 106 res.x = inarr_shape[0] 107 res.ravel_index = ravel_index1D 108 elif res.ndim == 2: 109 res.x = inarr_shape[1] 110 res.y = inarr_shape[0] 111 res.ravel_index = ravel_index2D 112 if res.x == 1 and res.y > 1: 113 # Should not occur, but better be safe than sorry 114 raise ValueError( 115 "Swap axis of your %s array so it has a %s shape" 116 % (inarr_shape, good_shape)) 117 elif res.ndim == 3: 118 res.x = inarr_shape[2] 119 res.y = inarr_shape[1] 120 res.z = inarr_shape[0] 121 res.ravel_index = ravel_index3D 122 if ((res.x == 1 and res.y > 1) 123 or res.y == 1 and res.z > 1): 124 # Should not occur, but better be safe than sorry 125 raise ValueError( 126 "Swap axes of your %s array so it has a %s shape" 127 % (inarr_shape, good_shape)) 128 else: 129 raise NotImplementedError( 130 "Only for images of dimension 1-3 are supported, got a %sD one" 131 % res.ndim) 132 133 res.numels = res.x * res.y * res.z 134 135 # When reading this for the first time, look at the diagram by the enum 136 # definition above (keyword D_ee) 137 # Difference between E and G is (x=0, y=-1, z=-1), E and A (-1, -1, 0) etc. 138 # Here, it is recalculated to linear (raveled) indices of flattened arrays 139 # with their last (=contiguous) dimension is x. 140 141 # So now the 1st (needed for 1D, 2D and 3D) part, y = 1, z = 1 142 res.DEX[D_ed] = -1 143 144 # Not needed, just for illustration 145 # res.DEX[D_ee] = 0 146 147 # So now the 2nd (needed for 2D and 3D) part, y = 0, z = 1 148 res.DEX[D_ea] = res.ravel_index(-1, -1, 0, res) 149 res.DEX[D_eb] = res.DEX[D_ea] + 1 150 res.DEX[D_ec] = res.DEX[D_eb] + 1 151 152 # And now the 3rd (needed only for 3D) part, z = 0 153 res.DEX[D_ef] = res.ravel_index(-1, -1, -1, res) 154 res.DEX[D_eg] = res.DEX[D_ef] + 1 155 res.DEX[D_eh] = res.DEX[D_ef] + 2 156 res.DEX[D_ei] = res.DEX[D_ef] - res.DEX[D_eb] # DEX[D_eb] = one row up, remember? 157 res.DEX[D_ej] = res.DEX[D_ei] + 1 158 res.DEX[D_ek] = res.DEX[D_ei] + 2 159 res.DEX[D_el] = res.DEX[D_ei] - res.DEX[D_eb] 160 res.DEX[D_em] = res.DEX[D_el] + 1 161 res.DEX[D_en] = res.DEX[D_el] + 2 162 163 164cdef inline void join_trees_wrapper(DTYPE_t *data_p, DTYPE_t *forest_p, 165 DTYPE_t rindex, DTYPE_t idxdiff) nogil: 166 if data_p[rindex] == data_p[rindex + idxdiff]: 167 join_trees(forest_p, rindex, rindex + idxdiff) 168 169 170cdef size_t ravel_index1D(size_t x, size_t y, size_t z, 171 shape_info *shapeinfo) nogil: 172 """ 173 Ravel index of a 1D array - trivial. y and z are ignored. 174 """ 175 return x 176 177 178cdef size_t ravel_index2D(size_t x, size_t y, size_t z, 179 shape_info *shapeinfo) nogil: 180 """ 181 Ravel index of a 2D array. z is ignored 182 """ 183 cdef size_t ret = x + y * shapeinfo.x 184 return ret 185 186 187cdef size_t ravel_index3D(size_t x, size_t y, size_t z, 188 shape_info *shapeinfo) nogil: 189 """ 190 Ravel index of a 3D array 191 """ 192 cdef size_t ret = x + y * shapeinfo.x + z * shapeinfo.y * shapeinfo.x 193 return ret 194 195 196# Tree operations implemented by an array as described in Wu et al. 197# The term "forest" is used to indicate an array that stores one or more trees 198 199# Consider a following tree: 200# 201# 5 ----> 3 ----> 2 ----> 1 <---- 6 <---- 7 202# | | 203# 4 >----/ \----< 8 <---- 9 204# 205# The vertices are a unique number, so the tree can be represented by an 206# array where a the tuple (index, array[index]) represents an edge, 207# so for our example, array[2] == 1, array[7] == 6 and array[1] == 1, because 208# 1 is the root. 209# Last but not least, one array can hold more than one tree as long as their 210# indices are different. It is the case in this algorithm, so for that reason 211# the array is referred to as the "forest" = multiple trees next to each 212# other. 213# 214# In this algorithm, there are as many indices as there are elements in the 215# array to label and array[x] == x for all x. As the labelling progresses, 216# equivalence between so-called provisional (i.e. not final) labels is 217# discovered and trees begin to surface. 218# When we found out that label 5 and 3 are the same, we assign array[5] = 3. 219 220cdef DTYPE_t find_root(DTYPE_t *forest, DTYPE_t n) nogil: 221 """Find the root of node n. 222 Given the example above, for any integer from 1 to 9, 1 is always returned 223 """ 224 cdef DTYPE_t root = n 225 while (forest[root] < root): 226 root = forest[root] 227 return root 228 229 230cdef inline void set_root(DTYPE_t *forest, DTYPE_t n, DTYPE_t root) nogil: 231 """ 232 Set all nodes on a path to point to new_root. 233 Given the example above, given n=9, root=6, it would "reconnect" the tree. 234 so forest[9] = 6 and forest[8] = 6 235 The ultimate goal is that all tree nodes point to the real root, 236 which is element 1 in this case. 237 """ 238 cdef DTYPE_t j 239 while (forest[n] < n): 240 j = forest[n] 241 forest[n] = root 242 n = j 243 244 forest[n] = root 245 246 247cdef inline void join_trees(DTYPE_t *forest, DTYPE_t n, DTYPE_t m) nogil: 248 """Join two trees containing nodes n and m. 249 If we imagine that in the example tree, the root 1 is not known, we 250 rather have two disjoint trees with roots 2 and 6. 251 Joining them would mean that all elements of both trees become connected 252 to the element 2, so forest[9] == 2, forest[6] == 2 etc. 253 However, when the relationship between 1 and 2 can still be discovered later. 254 """ 255 cdef DTYPE_t root 256 cdef DTYPE_t root_m 257 258 if (n != m): 259 root = find_root(forest, n) 260 root_m = find_root(forest, m) 261 262 if (root > root_m): 263 root = root_m 264 265 set_root(forest, n, root) 266 set_root(forest, m, root) 267 268 269def _get_swaps(shp): 270 """ 271 What axes to swap if we want to convert an illegal array shape 272 to a legal one. 273 274 Args: 275 shp (tuple-like): The array shape 276 277 Returns: 278 list: List of tuples to be passed to np.swapaxes 279 """ 280 shp = np.array(shp) 281 swaps = [] 282 283 # Dimensions where the array is "flat" 284 ones = np.where(shp == 1)[0][::-1] 285 # The other dimensions 286 bigs = np.where(shp > 1)[0] 287 # We now go from opposite directions and swap axes if an index of a flat 288 # axis is higher than the one of a thick axis 289 for one, big in zip(ones, bigs): 290 if one < big: 291 # We are ordered already 292 break 293 else: 294 swaps.append((one, big)) 295 # TODO: Add more swaps so the array is longer along x and shorter along y 296 return swaps 297 298 299def _apply_swaps(arr, swaps): 300 arr2 = arr 301 for one, two in swaps: 302 arr2 = arr.swapaxes(one, two) 303 return arr2 304 305 306def reshape_array(arr): 307 """ 308 "Rotates" the array so it gains a shape that the algorithm can work with. 309 An illegal shape is (3, 1, 4), and correct one is (1, 3, 4) or (1, 4, 3). 310 The point is to have all 1s of the shape at the beginning, not in the 311 middle of the shape tuple. 312 313 Note: The greater-than-one shape component should go from greatest to 314 lowest numbers since it is more friendly to the CPU cache (so (1, 3, 4) is 315 less optimal than (1, 4, 3). Keyword to this is "memory spatial locality" 316 317 Args: 318 arr (np.ndarray): The array we want to fix 319 320 Returns: 321 tuple (result, swaps): The result is the "fixed" array and a record 322 of what has been done with it. 323 """ 324 swaps = _get_swaps(arr.shape) 325 reshaped = _apply_swaps(arr, swaps) 326 return reshaped, swaps 327 328 329def undo_reshape_array(arr, swaps): 330 """ 331 Does the opposite of what :func:`reshape_array` does 332 333 Args: 334 arr (np.ndarray): The array to "revert" 335 swaps (list): The second result of :func:`reshape_array` 336 337 Returns: 338 np.ndarray: The array of the original shape 339 """ 340 # It is safer to undo axes swaps in the opposite order 341 # than the application order 342 reshaped = _apply_swaps(arr, swaps[::-1]) 343 return reshaped 344 345 346def label_cython(input_, background=None, return_num=False, 347 connectivity=None): 348 # Connected components search as described in Fiorio et al. 349 # We have to ensure that the shape of the input can be handled by the 350 # algorithm. The input is reshaped as needed for compatibility. 351 input_, swaps = reshape_array(input_) 352 shape = input_.shape 353 ndim = input_.ndim 354 355 cdef cnp.ndarray[DTYPE_t, ndim=1] forest 356 357 # Having data a 2D array slows down access considerably using linear 358 # indices even when using the data_p pointer :-( 359 360 # np.array makes a copy so it is safe to modify data in-place 361 data = np.array(input_, order='C', dtype=DTYPE) 362 forest = np.arange(data.size, dtype=DTYPE) 363 364 cdef DTYPE_t *forest_p = <DTYPE_t*>forest.data 365 cdef DTYPE_t *data_p = <DTYPE_t*>cnp.PyArray_DATA(data) 366 367 cdef shape_info shapeinfo 368 cdef bginfo bg 369 370 get_shape_info(shape, &shapeinfo) 371 get_bginfo(background, &bg) 372 373 if connectivity is None: 374 # use the full connectivity by default 375 connectivity = ndim 376 377 if not 1 <= connectivity <= ndim: 378 raise ValueError( 379 f'Connectivity for {input_.ndim}D image should ' 380 f'be in [1, ..., {input_.ndim}]. Got {connectivity}.' 381 ) 382 383 cdef DTYPE_t conn = connectivity 384 # Label output 385 cdef DTYPE_t ctr 386 with nogil: 387 scanBG(data_p, forest_p, &shapeinfo, &bg) 388 # the data are treated as degenerated 3D arrays if needed 389 # without any performance sacrifice 390 scan3D(data_p, forest_p, &shapeinfo, &bg, conn) 391 ctr = resolve_labels(data_p, forest_p, &shapeinfo, &bg) 392 393 # Work around a bug in ndimage's type checking on 32-bit platforms 394 if data.dtype == np.int32: 395 data = data.view(np.int32) 396 397 if swaps: 398 data = undo_reshape_array(data, swaps) 399 400 if return_num: 401 return data, ctr 402 else: 403 return data 404 405 406cdef DTYPE_t resolve_labels(DTYPE_t *data_p, DTYPE_t *forest_p, 407 shape_info *shapeinfo, bginfo *bg) nogil: 408 """ 409 We iterate through the provisional labels and assign final labels based on 410 our knowledge of prov. labels relationship. 411 We also track how many distinct final labels we have. 412 """ 413 cdef DTYPE_t counter = 1, i 414 415 for i in range(shapeinfo.numels): 416 if i == forest_p[i]: 417 # We have stumbled across a root which is something new to us (root 418 # is the LOWEST of all prov. labels that are equivalent to it) 419 420 # If the root happens to be the background, 421 # assign the background label instead of a 422 # new label from the counter 423 if i == bg.background_node: 424 # Also, if there is no background in the image, 425 # bg.background_node == BG_NODE_NULL < 0 and this never occurs. 426 data_p[i] = bg.background_label 427 else: 428 data_p[i] = counter 429 # The background label is basically hardcoded to 0, so no need 430 # to check that the new counter != bg.background_label 431 counter += 1 432 else: 433 data_p[i] = data_p[forest_p[i]] 434 return counter - 1 435 436 437cdef void scanBG(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, 438 bginfo *bg) nogil: 439 """ 440 Settle all background pixels now and don't bother with them later. 441 Since this only requires one linar sweep through the array, it is fast 442 and it makes sense to do it separately. 443 444 The purpose of this function is update of forest_p and bg parameter inplace. 445 """ 446 cdef DTYPE_t i, bgval = bg.background_val, firstbg = shapeinfo.numels 447 # We find the provisional label of the background, which is the index of 448 # the first background pixel 449 for i in range(shapeinfo.numels): 450 if data_p[i] == bgval: 451 firstbg = i 452 bg.background_node = firstbg 453 break 454 455 # There is no background, therefore the first background element 456 # is not defined. 457 # Since BG_NODE_NULL < 0, this is enough to ensure 458 # that resolve_labels doesn't worry about background. 459 if bg.background_node == BG_NODE_NULL: 460 return 461 462 # And then we apply this provisional label to the whole background 463 for i in range(firstbg, shapeinfo.numels): 464 if data_p[i] == bgval: 465 forest_p[i] = firstbg 466 467 468# Here, we work with flat arrays regardless whether the data is 1, 2 or 3D. 469# The lookup to the neighbor in a 2D array is achieved by precalculating an 470# offset and adding it to the index. 471# The forward scan mask looks like this (the center point is actually E): 472# (take a look at shape_info docs for more exhaustive info) 473# A B C 474# D E 475# 476# So if I am in the point E and want to take a look to A, I take the index of 477# E and add shapeinfo.DEX[D_ea] to it and get the index of A. 478# The 1D indices are "raveled" or "linear", that's where "rindex" comes from. 479 480 481cdef void scan1D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, 482 bginfo *bg, DTYPE_t connectivity, DTYPE_t y, DTYPE_t z) nogil: 483 """ 484 Perform forward scan on a 1D object, usually the first row of an image 485 """ 486 if shapeinfo.numels == 0: 487 return 488 # Initialize the first row 489 cdef DTYPE_t x, rindex, bgval = bg.background_val 490 cdef DTYPE_t *DEX = shapeinfo.DEX 491 rindex = shapeinfo.ravel_index(0, y, z, shapeinfo) 492 493 for x in range(1, shapeinfo.x): 494 rindex += 1 495 # Handle the first row 496 if data_p[rindex] == bgval: 497 # Nothing to do if we are background 498 continue 499 500 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) 501 502 503cdef void scan2D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, 504 bginfo *bg, DTYPE_t connectivity, DTYPE_t z) nogil: 505 """ 506 Perform forward scan on a 2D array. 507 """ 508 if shapeinfo.numels == 0: 509 return 510 cdef DTYPE_t x, y, rindex, bgval = bg.background_val 511 cdef DTYPE_t *DEX = shapeinfo.DEX 512 scan1D(data_p, forest_p, shapeinfo, bg, connectivity, 0, z) 513 for y in range(1, shapeinfo.y): 514 # BEGINNING of x = 0 515 rindex = shapeinfo.ravel_index(0, y, 0, shapeinfo) 516 # Handle the first column 517 if data_p[rindex] != bgval: 518 # Nothing to do if we are background 519 520 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) 521 522 if connectivity >= 2: 523 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) 524 # END of x = 0 525 526 for x in range(1, shapeinfo.x - 1): 527 # We have just moved to another column (of the same row) 528 # so we increment the raveled index. It will be reset when we get 529 # to another row, so we don't have to worry about altering it here. 530 rindex += 1 531 if data_p[rindex] == bgval: 532 # Nothing to do if we are background 533 continue 534 535 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) 536 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) 537 538 if connectivity >= 2: 539 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) 540 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) 541 542 # Finally, the last column 543 # BEGINNING of x = max 544 rindex += 1 545 if data_p[rindex] != bgval: 546 # Nothing to do if we are background 547 548 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) 549 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) 550 551 if connectivity >= 2: 552 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) 553 # END of x = max 554 555 556cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, 557 bginfo *bg, DTYPE_t connectivity) nogil: 558 """ 559 Perform forward scan on a 3D array. 560 561 """ 562 if shapeinfo.numels == 0: 563 return 564 cdef DTYPE_t x, y, z, rindex, bgval = bg.background_val 565 cdef DTYPE_t *DEX = shapeinfo.DEX 566 # Handle first plane 567 scan2D(data_p, forest_p, shapeinfo, bg, connectivity, 0) 568 for z in range(1, shapeinfo.z): 569 # Handle first row in 3D manner 570 # BEGINNING of y = 0 571 # BEGINNING of x = 0 572 rindex = shapeinfo.ravel_index(0, 0, z, shapeinfo) 573 if data_p[rindex] != bgval: 574 # Nothing to do if we are background 575 576 # Now we have pixels below 577 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) 578 579 if connectivity >= 2: 580 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) 581 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) 582 if connectivity >= 3: 583 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) 584 # END of x = 0 585 586 for x in range(1, shapeinfo.x - 1): 587 rindex += 1 588 # Handle the first row 589 if data_p[rindex] == bgval: 590 # Nothing to do if we are background 591 continue 592 593 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) 594 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) 595 596 if connectivity >= 2: 597 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) 598 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) 599 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) 600 if connectivity >= 3: 601 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) 602 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) 603 604 # BEGINNING of x = max 605 rindex += 1 606 # Handle the last element of the first row 607 if data_p[rindex] != bgval: 608 # Nothing to do if we are background 609 610 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) 611 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) 612 613 if connectivity >= 2: 614 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) 615 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) 616 if connectivity >= 3: 617 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) 618 # END of x = max 619 # END of y = 0 620 621 # BEGINNING of y = ... 622 for y in range(1, shapeinfo.y - 1): 623 # BEGINNING of x = 0 624 rindex = shapeinfo.ravel_index(0, y, z, shapeinfo) 625 # Handle the first column in 3D manner 626 if data_p[rindex] != bgval: 627 # Nothing to do if we are background 628 629 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) 630 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) 631 632 if connectivity >= 2: 633 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) 634 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) 635 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) 636 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) 637 if connectivity >= 3: 638 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) 639 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) 640 # END of x = 0 641 642 # Handle the rest of columns 643 for x in range(1, shapeinfo.x - 1): 644 rindex += 1 645 if data_p[rindex] == bgval: 646 # Nothing to do if we are background 647 continue 648 649 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) 650 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) 651 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) 652 653 if connectivity >= 2: 654 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) 655 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) 656 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) 657 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) 658 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) 659 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) 660 if connectivity >= 3: 661 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ef]) 662 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) 663 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) 664 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) 665 666 # BEGINNING of x = max 667 rindex += 1 668 if data_p[rindex] != bgval: 669 # Nothing to do if we are background 670 671 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) 672 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) 673 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) 674 675 if connectivity >= 2: 676 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) 677 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) 678 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) 679 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) 680 if connectivity >= 3: 681 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ef]) 682 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) 683 # END of x = max 684 # END of y = ... 685 686 # BEGINNING of y = max 687 # BEGINNING of x = 0 688 rindex = shapeinfo.ravel_index(0, shapeinfo.y - 1, z, shapeinfo) 689 # Handle the first column in 3D manner 690 if data_p[rindex] != bgval: 691 # Nothing to do if we are background 692 693 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) 694 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) 695 696 if connectivity >= 2: 697 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) 698 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) 699 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) 700 if connectivity >= 3: 701 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) 702 # END of x = 0 703 704 # Handle the rest of columns 705 for x in range(1, shapeinfo.x - 1): 706 rindex += 1 707 if data_p[rindex] == bgval: 708 # Nothing to do if we are background 709 continue 710 711 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) 712 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) 713 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) 714 715 if connectivity >= 2: 716 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) 717 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) 718 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) 719 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) 720 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) 721 if connectivity >= 3: 722 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ef]) 723 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) 724 725 # BEGINNING of x = max 726 rindex += 1 727 if data_p[rindex] != bgval: 728 # Nothing to do if we are background 729 730 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) 731 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) 732 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) 733 734 if connectivity >= 2: 735 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) 736 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) 737 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) 738 if connectivity >= 3: 739 join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ef]) 740 # END of x = max 741 # END of y = max 742