ray_trafo with impl=astra_cuda using Parallel2dGeometry

[davlars]

Re-visiting this repo after some time...

We're running some of the example scripts and are noticing some funny behaviour in the shepp_logan_russian_doll.py, specifically when combining a odl.tomo.Parlallel2dGeometry with a ray transform using impl=astra_cuda.

Modifying the example slightly to highlight the problem:

import numpy as np
import odl
import odl_multigrid as multigrid

# %%

# Basic discretization
coarse_discr = odl.uniform_discr([-10, -10], [10, 10], [100, 100])

fine_min = [-5, -5]
fine_max = [-2.5, -2.5]
fine_discr = odl.uniform_discr(fine_min, fine_max, [500, 500])

# Geometry
angle_partition = odl.uniform_partition(0, 2*np.pi, 360)
det_partition = odl.uniform_partition(-15, 15, 3000)

geometry = odl.tomo.Parallel2dGeometry(angle_partition, det_partition,
                                       det_pos_init=[20, 0])

# Mask
coarse_mask = multigrid.operators.MaskingOperator(coarse_discr,
                                                  fine_min, fine_max)
ray_trafo_coarse = odl.tomo.RayTransform(coarse_discr, geometry,
                                         impl='astra_cuda')
masked_ray_trafo_coarse = ray_trafo_coarse * coarse_mask

# Phantom
phantom_c = odl.phantom.shepp_logan(coarse_discr, modified=True)
phantom_f = odl.phantom.shepp_logan(fine_discr, modified=True)

# Ray trafo
ray_trafo_fine = odl.tomo.RayTransform(fine_discr, geometry,
                                       impl='astra_cuda')
pspace_ray_trafo = odl.ReductionOperator(masked_ray_trafo_coarse,
                                         ray_trafo_fine)
pspace = pspace_ray_trafo.domain

phantom = pspace.element([phantom_c, phantom_f])

data = pspace_ray_trafo([phantom_c, phantom_f])
data.show('data')

noisy_data = data + odl.phantom.white_noise(ray_trafo_coarse.range, stddev=0.1)
noisy_data.show('noisy data')

reco = pspace_ray_trafo.domain.zero()
multigrid.graphics.show_both(phantom_c, phantom_f)

callback = (odl.solvers.CallbackPrintIteration(step=2) &
            odl.solvers.CallbackShow(step=2))
odl.solvers.conjugate_gradient_normal(
    pspace_ray_trafo, reco, data, niter=20, callback=callback)
multigrid.graphics.show_both(*reco)

In the above, the insert finely discretized phantom is shifted off-center. However, in the generated data (forward projected product space) we clearly see remnants of a small phantom in the center. Similarly, the reconstruction gives a shaded fine-discretized-phantom in the middle (see below).

Trying to de-bug the issue it all seems to arise when calling data = pspace_ray_trafo([phantom_c, phantom_f]). Even more funny is that the issue does not arise if I switch to a ray-trafo with impl=astra_cpu or if using something else than a parallel beam geometry (e.g. odl.tomo.FanFlatGeometry).

That there are differences in implementation (and output) of astra_CPU vs astra_cuda seems to be a known issue, but this is slightly different.

Do you know what might be going on?

[kohr-h]

Hey David! Nice to see that something is happening on this front as well. You might want to look into the shepp_logan_russian_doll_2_datasets.py example, I've written it not so long ago, and it describes things much more clearly.

That said, I don't remember the core issue, but some of this will go away when parallel 2D vec geometries are supported by ASTRA. They are supported in the Git master version so you can build it yourself and fix this ODL code accordingly if you want. But exactly what is going on where is hard to tell from this example.

Maybe for diagnosis it's better to look into these test-like examples. You can introduce shifts and see what happens.

[davlars]

Yes! Moved on to the shepp_logan_russian_doll_2_datasets.py even before your suggestion, which hinted on the problem with parallel vs. fan/cone.

Going back to odl tests seems like a decent idea, so does the astra master. Since it's not hindering progression on other fronts (...) I'm not stuck, but it's definitely something I'll check out. Will keep the issue posted in case of any new findings.

Slightly off-topic but also: making a push now with this (finally!) so hopefully we'll progress somewhere on the multigrid front after a year of idle state... :-)