I’ve implemented a function that performs kernels 1 and 2 of the CCL algorithm found here: Traversing a matrix on a thread per row basis
The data is a binary 3D matrix (719x701x900) containing many particles as seen in the image below (I’ve left some orange noise for clarity):
Kernel 1 returns the spans and compressed label matrices as seen in the following sample:
Spans Matrix: Rows 698 to 718, Slice 0
Compressed Labels Matrix: Rows 698 to 718, Slice 0
Kernel 2 updates labels by contiguity and requires several iterations. It’s confirmed to work with very simple data. Below is a simple structure made up of 2 (50x50x50) cubes and one (250x25x25) parallelepiped and its finished spans and compressed label matrices. As can be seen, the function correctly identifies the spans and labels the object with a single label.
Spans Matrix: Rows 39 to 59, Slice 0 | Compressed Labels Matrix: Rows 39 to 59, Slice 0
However, the kernel doesn’t seem to work for larger data. There are no errors. The kernel simply blocks. In fact, it isn’t possible to retrieve any information - e.g. I can’t return the spans nor compressed labels matrices even after only 1 iteration.
Does anyone have any idea why this happens? Could it be the kernel is still working and, if so, is there any way to know?
Here’s a portion of the code for kernel 2:
@cuda.jit
def kernel_2(s_matrix, cl_matrix, changed):
row, image_slice = cuda.grid(2) # the position of the row. All positions start at the 0th column
pos_s, pos_l = 0, 0 # the position of the row in the spans and label matrices
pre_s, pre_l = 0, 0 # the previous row
post_s, post_l = 0, 0 # the following row
down_s, down_l = 0, 0 # the downwards row
up_s, up_l = 0, 0 # the upwards row
if row < s_matrix.shape[0] and image_slice < s_matrix.shape[2]: # guard for rows and slices
for column in range(cl_matrix.shape[1]):
current_span_label = cl_matrix[row, pos_l, image_slice]
if current_span_label == -1: # guard for a valid label i.e. it's not a null label
break
current_span_start = s_matrix[row, pos_s, image_slice]
current_span_end = s_matrix[row, pos_s + 1, image_slice]
# previous row
if image_slice > 0: # makes no sense to check the previous row of the first slice
while -1 < s_matrix[row, pre_s, image_slice - 1] < current_span_end:
if current_span_start <= s_matrix[row, pre_s + 1, image_slice - 1] \
and s_matrix[row, pre_s, image_slice - 1] <= current_span_end:
interval_label = cl_matrix[row, pre_l, image_slice - 1]
if interval_label == current_span_label:
pre_s = pre_s + 2
pre_l = pre_l + 1
continue
if interval_label == -1:
break
min_label = min(current_span_label, interval_label)
cuda.atomic.min(cl_matrix, (row, pos_l, image_slice), min_label)
cuda.atomic.min(cl_matrix, (row, pre_l, image_slice - 1), min_label)
cuda.atomic.add(changed, 0, 1)
pre_s = pre_s + 2
pre_l = pre_l + 1
# ...omitted code for other rows
pos_l = pos_l + 1
pos_s = pos_s + 2
Can’t believe I wasted so much time on this. Problem found. I wasn’t considering the case where the intervals don’t overlap.
All that’s needed is to… remove a tab from the following lines:
pre_s = pre_s + 2 # one tab back
pre_l = pre_l + 1 # one tab back