Output of Ray transform with cone beam geometry is asymmetric.

[JevgenijaAksjonova]

Hi,

I define a simple image with the central pixel equal to 1 and 0, otherwise. Then I project this image with a cone beam geometry. As a result, I expect the projected data be symmetric with respect to center of the detector, but this is not the case. I execute the following code

space = odl.uniform_discr([-1] * 3, [1] * 3, [3] * 3)
im = np.zeros((3,3,3))
im[1,1,1] = 1
phantom = space.element(im)

apart = odl.uniform_partition(0, np.pi, 2)
dpart = odl.uniform_partition([-1,-1] , [1,1] , [2,1])
geom = odl.tomo.ConeBeamGeometry(apart, dpart,1, 1)

op = odl.tomo.RayTransform(space, geom)
y = op(phantom)
y.asarray()

and get an output

array([[[ 0.3796193 ], [ 0.44643226]], [[ 0.44643226], [ 0.3796193 ]]])

In contrast, when I do the same thing in 2D with fan beam geometry I get

array([[ 0.44643226, 0.44643226], [ 0.44643226, 0.44643226]])

which corresponds to my expectations.

Is it a bug, or am I missing something?

[JevgenijaAksjonova]

I am using astra 1.9.9.dev2 and ODL from git

[adler-j]

Looks like a bug to me. With that said, you have a rather weird geometry, with the detector inside the object. I also did the math (might be wrong) and the correct value should be ~0.48686, so in a sense they are both rather wrong, just one more than the other.

[JevgenijaAksjonova]

With det_radius = src_radius = 2, in 3D I get

array([[[ 0.42841417], [ 0.44538102]], [[ 0.44538102], [ 0.42841414]]])

and 0.44538102 in 2D.

[adler-j]

To me this seems like "within the error margins" so to say. Astra does not use exact integration on GPUs, instead using an approximate numerical scheme. You could try a more finely sampled volume and see what happens, I'd expect it to converge to the correct value (~0.48, if my math is correct).

[JevgenijaAksjonova]

Unfortunatelly, it does not behave that way. I increased discretization to 30 and I get an output

array([[[ 0.22226636], [ 0.32298374]], [[ 0.28075498], [ 0.27698457]]])

[JevgenijaAksjonova]

def cub(d): s = d // 2 w = d//3 v = np.zeros(d, dtype=np.float32) v[(s- w//2):(s + w//2+1)] = 1.0 c = np.zeros([d]*3, dtype=np.float32) c[:,:, d//2] = np.outer(v, v) return c

plt.imshow(cub(d)[:,:,d//2])

this is the code I used to create an image

[adler-j]

I get the correct results:

import numpy as np
import odl

space = odl.uniform_discr([-1] * 3, [1] * 3, [300] * 3)
im = odl.phantom.cuboid(space, [-1/3] * 3, [1/3] * 3)
im.show()

apart = odl.uniform_partition(0, np.pi, 2)
dpart = odl.uniform_partition([-1,-1] , [1,1] , [2,1])
geom = odl.tomo.ConeBeamGeometry(apart, dpart,1, 1)

op = odl.tomo.RayTransform(space, geom)
y = op(im)
y.asarray()

array([[[ 0.48748174],
        [ 0.48648977]],
       [[ 0.4864898 ],
        [ 0.48748171]]])

This looks like it's correct to the analytical result down to rounding.

Colab here: https://colab.research.google.com/drive/1Q-z9dlN06zqsPd1SKMFw7_e5Siw5KgVu

[adler-j]

The reason is basically that ASTRA uses an approximate axis aligned algorithm. At the angles 45, 135, etc there is a breakpoint where it goes from integrating along x to integrating along y. Here it seems it choses two different methods which gives slightly different results. If you try apart = odl.uniform_partition(-0.000001, np.pi, 2) the result is more symmetric.

[JevgenijaAksjonova]

Ok, so it means that the error should correspond to discretization in the order of magnitude?

[adler-j]

Yes, I'd estimate it to be roughly on the order of "size of a voxel" or less. Although, I'm quite sure ASTRA has no actual guarantees.

[JevgenijaAksjonova]

Thank you for the clarification!