mismatch of RayTransform in astra_cpu and astra_cuda

[AceCoooool]

There is an "obvious" gap between RayTransform using impl='astra_cuda' and impl='astra_cpu. For example:

import numpy as np
import odl

size = 4
space = odl.uniform_discr([-size//2, -size//2], [size//2, size//2], [size, size], dtype='float32')

geometry = odl.tomo.parallel_beam_geometry(space, num_angles=2)
operator = odl.tomo.RayTransform(space, geometry, impl='astra_cpu')

operator_cuda = odl.tomo.RayTransform(space, geometry, impl='astra_cuda')

a = np.arange(0, size*size).reshape((size, size))

out = operator(a)
recs = operator.adjoint(out)

out_cuda = operator_cuda(a)
recs_cuda = operator_cuda.adjoint(out_cuda)

print(out)
print(out_cuda)

print(recs)
print(recs_cuda)

The recs and recs_cuda as follows:

recs

[[  62.04081 ,   37.632668,   38.44899 ,   66.122406],
 [  53.79593 ,   83.265305,   88.653046,   65.38776 ],
 [  57.061237,  104.81631 ,  110.20407 ,   68.65309 ],
 [  78.367294,   83.99999 ,   84.81635 ,   82.44899 ]]

recs_cuda

[[  54.845592,   42.754524,   44.19093 ,   59.154797],
 [  57.199715,   79.81814 ,   84.6332  ,   68.26626 ],
 [  62.94532 ,   99.07839 ,  103.893456,   74.011856],
 [  72.08241 ,   87.02065 ,   88.45705 ,   76.39162 ]]

This is due to the backend of astra-toolbox, may related question: astra-issue38

[AceCoooool]

after I test the astra create_sino. I think it's better to set the default RayTransform cpu version with proj_type='linear' rather than 'line' (However, I find it's strange of parallel_beam_geometry in deal with the "detector space", such as proj_space). The demo code as follows:

import astra
import numpy as np

size = 2
vol_geom = astra.create_vol_geom(size, size)
proj_geom = astra.create_proj_geom('parallel', 1.0, int(np.ceil(size * np.sqrt(2))), np.linspace(0, np.pi, 2, False))

proj_id_cuda = astra.create_projector('cuda', proj_geom, vol_geom)
proj_id = astra.create_projector('linear', proj_geom, vol_geom)

P = np.arange(0, size * size).reshape((size, size))

sinogram_id, sinogram = astra.create_sino(P, proj_id)
sinogram_id_cuda, sinogram_cuda = astra.create_sino(P, proj_id_cuda)

print(sinogram)
print(sinogram_cuda)

And when using "proj space" not equal to 1.0, the astra_cpu and astra_gpu also give different results. !!! (And the library using space not equal to 1.0 in most cases.)

[kohr-h]

Thanks for bringing up the issue, @AceCoooool. The root cause is likely ASTRA as you point out yourself. Regarding choice of projection kernel, we probably overthought this whole thing with correspondence of interpolation scheme and projector type a bit. If I had to do it again, I'd simply give users a way to choose for themselves, instead of doing it automatically, and just choose a sensible default.

[AceCoooool]

It's may convenient change parallel_beam_geometry parameters (det_shape is inconvenient to user to change "detector space" --- due to int form, the detector space can not be 1.0 in most cases)

Maybe, this advice is bad. haha :smile: