Astra curved detector

[NogginBops]

This PR is based on #1634 but remade for odl 1.0.0.

This PR uses cyl_cone_vec from ASTRA 2.4.0 to allow curved cone beam detectors.

From #1634 there was a note about a off-by-half bug related to pixel centers that I'm not sure if it still applies.

Known bug: there is a misinterpretion of the detector center if the detector partition isn't centered around 0, which is the case for quarter-pixel-shifted detectors. Fixing this needs to be coordinated with the Mayo CT data loader, which currently is assuming the presence of this bug.

There is also a FIXME in astra_cuda.py related to not being able to arbitrarily transposing pytorch tensors. Not sure what the best solution for that is, guidance is appreciated. :) This is fixed by using permute_dims.

@leftaroundabout could you take a look at and see if you can see anything that is obviously wrong that I could fix? You can do a full review when you have time.

[NogginBops]

This is the transpose issue. We store transpose info in transpose_tuple to be able to do the transpose at a later time. But for pytorch tensors a transpose can only swap the positions of two dimensions and not arbitrarily rearrange them like we need to do in this case. With two transposes we can get the correct format, but I'm not sure what the best way to store that information is.

I could make a transpose_tuple2 that only gets used for pytorch tensors and only if transpose_tuple2 != None but I'm not sure if this is a good solution, seems a little bit messy.

[leftaroundabout]

I think permute_dims is now the preferred tool to use.

[NogginBops]

I've managed to use permute_dims but as I'm unfamiliar with the standardized python array API I might be doing something weird.

[leftaroundabout]

At a glance, it looks reasonable. Though, better use proj_space.array_namespace and store it in a variable, instead of repeatedly looking up proj_data.__array_namespace__.

[NogginBops]

Do you think it matters enough that I should change it?
I don't know how expensive __array_namespace__ really is, but I expect that compared to the other operations in these functions I expect it to be negligible.

[leftaroundabout]

Performance-wise the difference is likely negligible, but it's just a bit more structured.

[NogginBops]

I think I've fixed it now.

[NogginBops]

There is a bug in this code that I will investigate tomorrow, but because creating pull request threads is down (https://www.githubstatus.com/) can't comment this on the affected line of code.

On line 253 in astra_cuda.py:
I have observed a crash here that I will investigate tomorrow:

Traceback (most recent call last):
  File "/workspace/reconstruction.py", line 182, in <module>
    main()
  File "/workspace/reconstruction.py", line 133, in main
    resid = sinogram - A(x)   # residual in data space
                       ^^^^
  File "/opt/odl/src/odl/core/operator/operator.py", line 663, in __call__
    out = self._call_out_of_place(x, **kwargs)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/odl/src/odl/applications/tomo/operators/ray_trafo.py", line 304, in _call
    return self.get_impl(self.use_cache).call_forward(x, out, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/odl/src/odl/applications/tomo/backends/util.py", line 47, in wrapper
    return fn(self, x, out, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/odl/src/odl/applications/tomo/backends/astra_cuda.py", line 175, in call_forward
    return self._call_forward_real(x, out, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/odl/src/odl/applications/tomo/backends/astra_cuda.py", line 253, in _call_forward_real
    return self.proj_space.element(proj_data)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/odl/src/odl/core/discr/discr_space.py", line 369, in element
    self, self.tspace.element(inp, copy=copy)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/odl/src/odl/core/space/base_tensors.py", line 686, in element
    return wrapped_array(arr)
           ^^^^^^^^^^^^^^^^^^
  File "/opt/odl/src/odl/core/space/base_tensors.py", line 625, in wrapped_array
    raise ValueError(
ValueError: shape of `inp` not equal to space shape: (1560, 230, 11066) != (11066, 1560, 230)

[leftaroundabout]

@NogginBops looks like permute_dims behaves in a subtly different way to the old NumPy transpose. Not very surprising... in my experience this kind of thing always requires a lot of trial and error to get right in Python.
Good work regardless, getting it this far!

[NogginBops]

I've seemingly fixed the issue. The issue was that the old permutation had order 2 which meant applying the same permutation again returned to the original shape, but the new permutation is order 3 which means applying the same permutation again does not return to the dimentions to the original order. The fix is to introduce explicitly specify the inverse permutation.

[leftaroundabout]

@NogginBops sounds reasonable. Thanks a lot for your efforts!

I don't see anything wrong with the code, but still need to do a bunch of testing before I can merge it.

[leftaroundabout]

Well, this backprojection seems to suggest that at least the axes are correctly mapped:

detector_partition = odl.uniform_partition([-np.pi/8, -160], [np.pi/8, 160], [512, 512])
geometry = odl.applications.tomo.ConeBeamGeometry(
    angle_partition, detector_partition, src_radius=200,
    det_radius=200, det_curvature_radius=(400, None),
    axis=[1, 0, 0])

(everything else the same as in ray_trafo_cone_3d).

What's strange is that I have to make the detector very "oversized" (so the projection of the actual phantom only takes up a small part of the data space, which would be rather inefficient). If not, then the backprojection comes out heavily cropped on the sides, even though the projection is still well within the detector area:

Not sure if this is because of some actual bug in the code, or I'm just interpreting the results wrong.