adjoint in odl.tomo.RayTransform is incorrect

[hongtao-argmin]

The forward and adjoint operators are incorrect by using the following commands:

ray_trafo = odl.tomo.RayTransform(reco_space, geometry) proj_data = ray_trafo(phantom) backproj = ray_trafo.adjoint(proj_data)

check: y = torch.randn_like(proj_data) torch.sum(y proj_data) \neq torch.sum( ray_trafo.adjoint(y) phantom)

for different number of views.

[mehrhardt]

ODL's adjoint is not with respect to the inner product given by "torch.sum(y*x)" but rather y.inner(x).

[hongtao-argmin]

Many thanks. For my understanding, the results of inner product should give us a scalar? So y.inner(x)'' will not give us a scalar so that looks weird if we define the inner product like this? If we usey=y.flatten()'' and ``y.inner(x.flatten)'', then it is equivalent to torch.sum(y*x).

Definition of inner product: https://en.wikipedia.org/wiki/Inner_product_space.

[mehrhardt]

Sorry, I can't follow your argument regarding scale.

I don't know from the top of my head what y=y.flatten() gives you. Will this be a numpy array? Is it even defined?

In any case, if you want to have that ODL computes the adjoint with respect to the inner product of torch, you could just compute constant scaling factors alpha and beta and correct for it:

ray_trafo = odl.tomo.RayTransform(reco_space, geometry)
proj_data = ray_trafo(phantom)
backproj = ray_trafo.adjoint(proj_data)

alpha = backproj.inner(phantom) / torch.sum(backproj.to_array() * phantom.to_array())
beta = proj_data.inner(backproj) / torch.sum(proj_data.to_array() * backproj.to_array())

Op = beta * ray_trafo
Adj = alpha * ray_trafo.adjoint

Note, that I haven't tested the code above but hopefully you get the idea.

[hongtao-argmin]

Thanks. y=y.flatten()'' (defined in PyTorch) will lety'' become a vector since images are represented in multi-dimension.

For my understanding, for any forward model A'', we should have<A^Tx,y> = <y,Ax>'' that A^T is the adjoint of A and the inner product <>'' will give us a scalar. However, x.inner(y) will not give us a scalar since images are represented in multi-dimension. thatalpha'' or ``beta'' showed in the template code is not a scalar.

[mehrhardt]

with "scale", you mean scalar https://en.wikipedia.org/wiki/Scalar_(mathematics)?

[hongtao-argmin]

yeah, sorry about the typo error. I corrected the word.

[JevgenijaAksjonova]

Hi,

The definition of the inner product in L2 function space is <f, g> = \int f(x) g(x) dx. Discretizing the above yields \sum f_i g_i cell_volume. Hence, for the ray transform it is <Ax, y > A.range.cell_volume = <x, A^T y> A.domain.cell_volume, where <a,b> = sum a_i b_i

[hongtao-argmin]

Thank you so much.

Suppose we worked on the 2D case and the image size is [m, n] and the size of the measurements is [nviews,detector_size]. Then, we should have <Ax,y>(nviewsdetector_size)= <x,A^T y>(mn). -- x represents the image and y denotes the measurements.

If my understanding is correct, then looks the ``='' is still not held. Or is it possible for you to show me one example that <Ax, y > A.range.cell_volume = <x, A^T y> A.domain.cell_volume is held for a given geometry but different image size, nviews, detector_size? I think this would help me a lot. Thank you so much.

[JevgenijaAksjonova]

If A = odl.tomo.RayTransform(reco_space, geometry), then use A.range.cell_volume and A.domain.cell_volume. Here, domain and range are input and output spaces of your operator. When you define the ray transform, you provide the domain (reco_space), range is computed automatically. space.cell_volume = product of space.cell_sides. Space.cell_sides = ( space.max_pt - space.min_pt) / space.shape

[hongtao-argmin]

It works now. Thank you so much. This really helps. :)