Crash when trying forward projection with flipped head phantom

[pjueon]

I'm trying to run Python/demos/d03_generateData.py with the head phantom that is flipped upside down. If I run the original demo code, it works. But if I add this one line before the forward projection, the program crashes.

head = np.flip(head, axis=0)  # flip the phantom data upside down

I'm not sure but it might be an exception handling issue between python and Cuda C++.

Or am I missing something?

Expected Behavior

The program doesn't crash.

Actual Behavior

The program crashes.

Code to reproduce the problem (If applicable)

import tigre
import numpy as np
from tigre.utilities import sample_loader
from tigre.utilities import CTnoise

geo = tigre.geometry_default(high_resolution=False)

#%% Define angles of projection and load phatom image

# define projection angles (in radians)
angles = np.linspace(0, 2 * np.pi, 50)
# load phatnom image
head = sample_loader.load_head_phantom(geo.nVoxel)


head = np.flip(head, axis=0) # flip the phantom data upside down
print("flip done")

# Simulate forward projection.
projections = tigre.Ax(head, geo, angles)   # crashes here!

print("Ax done")   # not excuted

Specifications

• Python version: 3.9.2, with venv.
• OS: Windows 11
• CUDA version: Cuda 11.6

[AnderBiguri]

I wonder if this has to do with np.flip() returning a view and not a copy.

Try to do

head = np.flip(head, axis=0).copy() # flip the phantom data upside down

[pjueon]

@AnderBiguri Thanks! That does work.

One more thing:

I also tried this one

head = np.flip(head.copy(), axis=0) 

and it didn't work. Maybe flipped view causes the crash (could be a bug)?

[pjueon]

Here, only casting img.data to C++ pointer might be the problem. I'm not familiar with the numpy implementation but probably the pointer that img.data attribute returns is not pointing the position we desire and the memory order is not flipped even if img is flipped.

https://github.com/CERN/TIGRE/blob/54a854bec2bc5ebf4184c700d6298b31f00c1ec3/Python/tigre/utilities/cuda_interface/_Ax.pyx#L73-L83

[pjueon]

Ok, after googling about the memory order of numpy ndarray, I think I found how to fix this. I`ll make a PR when I get home.

[AnderBiguri]

@pjueon yes, it makes sense that it only works when the np.array is in the C order In your second case, you make ac opy, but what np.flip returns is still flipped in indices, but not in memory, while in my case, the copy enforces the indices+memory map to be "the right way".

I am not sure what the solution is, other than telling the users to make the copy. I don't want to enforce copies, because sometimes that array can be VERY large. Perhaps in Ax.py we can detect if the array is in the right order (i.e. are the indices and memory map the same order, or has it not been a flip/transpose) and if so, raise an error and tell the user to make a copy?

What do you think?

[pjueon]

I've just created the PR for this. (#400) The trick is using numpy.ascontiguousarray or checking numpy.ndarray.flags['C_CONTIGUOUS'].