cwatershed fails on large 3D images (>2**31 voxels)

[williamjshipman]

I was attempting to process a 2304x2304x2475 voxel tomogram using cwatershed. I'd keep getting an array of zeros out. I eventually found that just processing 100 slices worked fine, so I kept increasing the number of slices until it failed. It seems like the issue is that as soon as the no. of voxels exceeds 2**31, cwatershed fails (either giving an array of zeros, or once even causing Python to segfault). Looking at _morph.cpp, I think the problem is that int is used for indices, e.g. line 599, instead of int64_t or uint64_t.

[luispedro]

Wow, I thought 2³¹ voxels should be enough for everybody, but I guess I was wrong (seriously, I never even thought this could be an issue).

Obviously, the fix is easy, so I will just do it.

[luispedro]

I just committed what I think is a fix. Do you have the ability to test the github version and tell me if it works?

[williamjshipman]

Thanks, will test it today.

…

On Mon, 18 Mar 2019, 5:55 AM Luis Pedro Coelho, ***@***.***> wrote: I just committed what I think is a fix. Do you have the ability to test the github version and tell me if it works? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#102 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AJEqgF-Zx5Lv3cxxf5p5chzrmA9wLs8Dks5vXw4igaJpZM4b313A> .

[williamjshipman]

Bad news, didn't fix it. Going above 2**31 voxels now leads to a seg fault, no exception error message to give you unfortunately.

Here's my very simple test code, if it helps. Change the 512 slices to 511 slices and it seems to work fine. With 512 slices, it seg faults.

import numpy as np
import mahotas
from matplotlib import pyplot as plt

xx, yy, zz = np.meshgrid(np.arange(2048), np.arange(2048), np.arange(512))
dist = np.sqrt((xx - 1024) ** 2 + (yy - 1024) ** 2 + (zz - 256) ** 2)
del xx, yy, zz
seed_img = np.zeros(dist.shape, dtype=np.uint64)
dist = dist.astype(np.float32)
wsimg = mahotas.cwatershed(dist, seed_img, Bc=np.ones((3,3,3)))
wsimg[dist > 100] = 0
plt.imsave('ws_test_img.png', wsimg[:, :, 256])

[luispedro]

I was trying to test this, but it's incredibly slow to work with such very large matrices on the only machine to that I have access to that has enough memory for this to work (it's a pretty old and slow machine, but has a lot of memory).

I switched a few remaining uses of int (32 bit s) to npy_intp (64 bits), so this may now work (if you wait long enough, at least it hasn't crashed yet).

[williamjshipman]

Thanks, I'll try to test it. I don't have access to the computer I was using, but I'll see what I can do.

…

On Mon, 8 Apr 2019 at 05:33, Luis Pedro Coelho ***@***.***> wrote: I was trying to test this, but it's incredibly slow to work with such very large matrices on the only machine to that I have access to that has enough memory for this work is pretty old and slow. I switched a few remaining uses of int (32 bit s) to npy_intp (64 bits), so this may now work (if you wait long enough, at least it hasn't crashed yet). — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#102 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AJEqgHSAc-3HU3EGpFcaS1cHgeheyulAks5veriVgaJpZM4b313A> .