Interpolation for curved detectors

[JevgenijaAksjonova]

Hi!

There is a possibility to specify a curved detector in ODL now, but no back-end to support it. So, the best thing to do now is to use interpolation to flat detector. What if we do the interpolation somewhere ray_trafo.py?

[ozanoktem]

To clarify, freely available software packages for computing ray transform and corresponding backprojection (like ASTRA, RTK, CASToR) do not support 3D helical CT scanning with curved detectors.

ODL has data structures to represent 3D helical scanning geometries and also data on curved detectors (#1478), but there is no back-end that handles the combination of both. For this reason one needs to interpolate from curved to flat detector. The question is where to put this interpolation code. Is it as part of the ray transform operator, part of the geometry object or in the binding code?

[kohr-h]

Detector curvature is already part of the geometry, as @JevgenijaAksjonova wrote, so that part is clear. What is not quite clear to me is where the interpolation code should live.

• It could be in the different ..Impl classes, but that would mean a lot of duplication and would go against the idea that those wrappers should be as thin as posible.
• Another option would be to place the interpolation somewhere in the RayTransform._call method. But that would mean to go from the current nice return self.get_impl(self.use_cache).call_forward(x, out, **kwargs) to something more complicated that involves interpolation at that point. It seems odd to do it there.
• The third option is to make it a separate operator and (if needed) a convenience function to create a ray transform for curved detector. We do pay a performance penalty with this, but I don't think it would be detrimental.
• The forth, but unfortunately currently unavailable option would be to use ASTRA's hypothetical native backend for curved detectors :wink:

[ozanoktem]

Adopting the ODL operator centric viewpoint, the interpolation from curved to flat detector is (similar to rebinning) a mapping between two geometry objects. So why not have it as an explicit operator instead of hiding it?

[kohr-h]

@ozanoktem Agreed, I prefer that option, too.

[ozanoktem]

@JevgenijaAksjonova , here is a suggestion on how to proceed:

1. Make sure ODL code for ray transform includes a check whether the back-end supports curved detectors.
2. Define interpolation as a mapping between two tomographic geometry objects.

[JevgenijaAksjonova]

I see interpolation from curved to flat detector primarily as a mapping between the data, because it is what take the most effort. Right now mapping between the data is done as following

# Create a *temporary* ray transform (we need its range)
spc = odl.uniform_discr([-1] * 3, [1] * 3, [32] * 3)
ray_trafo = odl.tomo.RayTransform(spc, geometry, interp='linear')

# convert coordinates
theta, up, vp = ray_trafo.range.grid.meshgrid
d = src_radius + det_radius
u = d * np.arctan(up / d)
v = d / np.sqrt(d**2 + up**2) * vp

# Calculate projection data in rectangular coordinates since we have no
# backend that supports cylindrical
interpolator = linear_interpolator(data_array, (theta, u, v))
proj_data = ray_trafo.range.element(interpolator)

Mapping between geometries can be done as following. Instead of creating geometry with

dpart = unifrom_partition(-a,a, n)
det_curvature_rad = r

one can create a geometry with

b = r*tg(a/r)
dpart = uniform_partition(-b,b, n)
det_curvature_rad = 0