Transposed output for certain inputs

[ucabova]

I noticed that using the Python API can result in transposed output in some cases, causing issues in subsequent processing where the connected component output isn't aligned with the original input image. I suspect that it has to do with the way numpy arrays are stored in the memory. I think the issue is caused by Fortran memory ordering. With C ordered input arrays I didn't experience the transposed output issue. Even with Fortran ordering there need to be some specific conditions, as not all Fortran ordered inputs lead to the transposed output. A default C ordered numpy array can become "F" ordered by transposing it, but I loading an nrrd file using pynrrd can also result in "F" ordered numpy arrays. I'd recommend changing the code so that if a numpy array is passed as input, first it is checked and ensured to be "C" ordered . I managed to illustrate the issue with the following code, where without modifications the output is transposed compared to the input:

import numpy as np
import matplotlib.pyplot as plt
import cc3d

a = np.zeros((50, 50, 50), order="F") # the array can be "F" ordered if it's loaded using pynrrd
a.T[10:20, 20:30, 30:40] = 1 
a[:20, :30, :40] = 2  # this and the above line are both needed to reproduce the issue
lbls_out, N_blobs = cc3d.largest_k(a, k=1, connectivity=6, return_N=True)
largest_cc = lbls_out.T == 1 # if lbls_out is transposed here, then the largest connected component is aligned with the input on the visualization below
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
i = 5
ax[0].imshow(a[i])
ax[0].set_title("input")
ax[1].imshow(largest_cc[i])
ax[1].set_title("labels out")

The input and the output can be aligned by:

1. initializing the array with "C" ordering
2. removing any of the two lines where values of the array are set to 1 or 2
3. transposing lbls_out

[william-silversmith]

Hi,

I need to take a closer look at the output again, but I believe this issue is due to the exclusive use of Fortran order for writing the output for largest_k in set_run_voxels when the output type of largest_k can be C or F.

[william-silversmith]

Hi @ucabova,

I don't think I'm able to reproduce your issue. I see that the output appears to be transposed, but the code is intentionally transposing it on the last line.

Would you mind editing the below test code to produce a failing case?

def test_largest_k_transpose_bug():
  a = np.zeros((50, 50, 50), order="C") # the array can be "F" ordered if it's loaded using pynrrd
  a.T[10:20, 20:30, 30:40] = 1
  a[:20, :30, :40] = 2 # this and the above line are both needed to reproduce the issue

  lbls_out, N_blobs = cc3d.largest_k(a, k=1, connectivity=6, return_N=True)
  largest_cc = lbls_out.T == 1 
  
  gt = np.zeros((50,50,50), order="F")
  gt[:20, :30, :40] = 1

  assert np.all(largest_cc == gt.T)

[ucabova]

Hi @william-silversmith,

Many thanks for looking into this! I noticed that my original example was confusing. The way I provided aligns the (otherwise transposed) output with the input by transposing lbls_out when largest_cc is created:
largest_cc = lbls_out.T == 1 # if lbls_out is transposed here, then the largest connected component is aligned with the input on the visualization below
As I expect the connected component function to work, this shouldn't be needed. However, without transposing there the output is misaligned (transposed compared to the input).

Regarding your code: to reproduce the issue array a should be "F" ordered, that is key. The order of gt doesn't matter. If you change the order of a to "F" in your code you can reproduce the issue. Normally, transposing shouldn't be needed, neither when the largest_cc is created nor for gt. Please try the modified test below. Here gt is never transposed, and the lbls_out is transposed based on the transpose parameter of the test.

def test_largest_k_transpose_bug(order_a="F", order_gt="C", transpose=False):
    a = np.zeros((50, 50, 50), order=order_a) # the array can be "F" ordered if it's loaded using pynrrd
    a.T[10:20, 20:30, 30:40] = 1
    a[:20, :30, :40] = 2 # this and the above line are both needed to reproduce the issue
    lbls_out, N_blobs = cc3d.largest_k(a, k=1, connectivity=6, return_N=True)
    if transpose:
        largest_cc = lbls_out.T == 1
    else:
        largest_cc = lbls_out == 1 
    
    gt = np.zeros((50,50,50), order=order_gt)
    gt[:20, :30, :40] = 1
    
    assert np.all(largest_cc == gt)

Both test_largest_k_transpose_bug("F", "C", False) and test_largest_k_transpose_bug("F", "F", False) give AssertionError. Any other combinations of parameters are fine.

[william-silversmith]

Ok I managed to make the test fail. But it seems that the key to making the test fail is transpose, which naturally makes the labels different than expected.

automated_test.py::test_largest_k_transpose_bug[True-F-F] FAILED                                                         [ 12%]
automated_test.py::test_largest_k_transpose_bug[True-F-C] FAILED                                                         [ 25%]
automated_test.py::test_largest_k_transpose_bug[True-C-F] FAILED                                                         [ 37%]
automated_test.py::test_largest_k_transpose_bug[True-C-C] FAILED                                                         [ 50%]
automated_test.py::test_largest_k_transpose_bug[False-F-F] PASSED                                                        [ 62%]
automated_test.py::test_largest_k_transpose_bug[False-F-C] PASSED                                                        [ 75%]
automated_test.py::test_largest_k_transpose_bug[False-C-F] PASSED                                                        [ 87%]
automated_test.py::test_largest_k_transpose_bug[False-C-C] PASSED                                                        [100%]

@pytest.mark.parametrize("order_a", ["F", "C"])
@pytest.mark.parametrize("order_gt", ["F", "C"])
@pytest.mark.parametrize("transpose", [True, False])
def test_largest_k_transpose_bug(order_a, order_gt, transpose):
    a = np.zeros((50, 50, 50), order=order_a) # the array can be "F" ordered if it's loaded using pynrrd
    a.T[10:20, 20:30, 30:40] = 1
    a[:20, :30, :40] = 2 # this and the above line are both needed to reproduce the issue
    lbls_out, N_blobs = cc3d.largest_k(a, k=1, connectivity=6, return_N=True)
    if transpose:
        largest_cc = lbls_out.T == 1
    else:
        largest_cc = lbls_out == 1 
    
    gt = np.zeros((50,50,50), order=order_gt)
    gt[:20, :30, :40] = 1
    
    assert np.all(largest_cc == gt)

[ucabova]

Yes, the example is not the most straightforward; it doesn't represent a typical use case. I encountered the issue by using segmentation masks saved as NRRD and loading them using pynrrd. I think some specific numpy operations may have caused the issue in the preprocessing before saving the array to the nrrd file. I can't share the file, and this was a way I found to reproduce the issue. I realize it may be a niche problem not impacting many users, but including the following check in the cc3d functions is also not too expensive:

if input_array.flags['F_CONTIGUOUS']:
    input_array = np.ascontiguousarray(input_array)

[william-silversmith]

Hi @ucabova,

I believe you that something is going on, but it's pretty difficult for me to apply an appropriate fix without a reproducible test case. If there's memory corruption introduced somewhere else, that could be impacting cc3d. Or perhaps I made a bad binary? Have you tried installing from source?

As for the check, I do have something like that for non-contiguous data:

  if not data.flags.c_contiguous and not data.flags.f_contiguous:
    data = np.copy(data, order=order)

However, it's... not entirely true that np.ascontiguousarray(input_array) is cheap. Most of my data is in F order. Many people apply cc3d in a loop, and in that context, repeated application of np.ascontiguousarray can become a significant cost if their data is not in C order.