Bug using Tam-Danielsson window and FBP

[chongchenmath]

There seems to be a bug in how the Tam-Danielsson window is used in the adutils. Currently the code in the adutils reads (see lines 115-119)

if use_window:
    window = odl.tomo.tam_danielson_window(ray_trafo,
                                           smoothing_width=0.05,
                                           n_half_rot=3)
    ray_trafo = window * ray_trafo

which will cause ray_trafo to not be a RayTransform anymore. Thus the standard fbp in odl does not work. Minimal working example to reproduce the crash:

import numpy as np
import adutils

# Discretization
reco_space = adutils.get_discretization()

# Forward operator (in the form of a broadcast operator)
A = adutils.get_ray_trafo(reco_space, use_window=True)

# Define fbp
fbp = adutils.get_fbp(A)

This can be compared to how the Tam-Danielsson window is applied in the example filtered_backprojection_helical_3d in odl.

edit: @adler-j improved typesetting #

[adler-j]

Well, isn't the bug in get_fbp then? Should be a rather simple fix.

Also: why do we have the use_window option? To improve convergence?

[aringh]

Long story short: for other reasons @chongchenmath and I had a closer look at the code. We thus found the use_window option, and after feedback from @ozanoktem on Monday we "side-tracked" into having a quick look at the Tam-Danielsson window, to understand what it does and how it is used.

Since I do not (yet) understand when and how the Tam-Danielsson window is supposed to be used, I cannot answer whether it is reasonable to have it in get_ray_trafo or in get_fbp. All I know is that if one set use_window = True, the get_fbp does not work.

I can try to fix it, but in which of the two places? Because I guess one should only use one of the two windows, no?

[adler-j]

Ok, sounds reasonable!

Regarding fixing, I think the best would be to do the (somewhat hacky) thing and add a test inside get_fbp like (with a note that the names are from the top of my head):

if isinstance(ray_transform, OperatorVectorMultiplciation):
    ray_transform = ray_transform.operator

P.S. I was perhaps a bit harsh in my first reply and I should note that this was an exceptionally good bug report! Conciece, reproducable and clear!

[aringh]

That was the fix I was thinking of, but I did not know if you would like it ;-) I will try to have it done later today

[adler-j]

There is no one-to-one relationship between what I like and what should be done :) Anyway sounds good!

[davlars]

Caught up with this now. Assuming you're on it @aringh, simply do a pull request whenever you have it all compiled. Or update whatever you've found here and we'll try to fix it. Apart from that I agree with the rest: good bug-spotting!

[davlars]

With #33 I guess this is done.