Automatic range inference of RayTransform breaks for unweighted space

[JevgenijaAksjonova]

I am trying to estimate the norm of the Radon transform:

xlim = 14
space = odl.uniform_discr([-xlim, -xlim], [xlim, xlim], [28, 28],
                          dtype='float32')

geometry = odl.tomo.parallel_beam_geometry(space, num_angles=5)
operator = odl.tomo.RayTransform(space, geometry)
print(operator.norm(estimate=True))

The result is 11.8375811049. However, when I am changing xlim to 14.0000000000001, the result is 9.22141281723, which does not correspond to my expectations that an operators norm should be continuous with respect to the input.

[adler-j]

That sounds very shady indeed. What impl are you using? Could you try "all of them"? Does this also happen for other geometries (e.g. fan beam?)

[JevgenijaAksjonova]

I have tried all implementations, the numbers are slightly different, but the inconsistency persists.

The same happens, while using Parallel2dGeometry and FanBeamGeometry.

[adler-j]

What happens with 14.0? Could this be an integer issue? Are you running python 2 or 3?

[JevgenijaAksjonova]

Python 3. Same problem with 14.0.

I discovered this problem when I was trying to remove the influence of space cell volume on the operator norm. The dependency correspond to my theoretical considerations apart from the case when operator.domain.cell_volume = 1, which seem to be an outlier.

[adler-j]

Oh god damnit, now I know what it is. I knew this would happen :D @kohr-h!

This is because we no longer have properly unweighted spaces, only spaces with "no weight" so to say. This is as mentioned here:

I don't see how the weightings of domain and range are connected for ray transforms. They're totally different spaces.

See e.g. our old conversation:

Where did this go? Overall NoWeighting is different from ConstWeighting (conceptually at least), For example, if we make a ray transform with domain odl.uniform_discr([-1, -1], [1, 1], [2, 2]), we still want the range to be weighted, even though the domain has weighting 1 and is thus in some sense equivalent to NoWeighting.

@adler-j

The only situation where weighting vs. no weighting makes a difference is when the cell sizes are changed from 1 to something else, or the other way around; plus the case when axes are removed or added. The only place where we change cell sizes is uniform_discr_fromdiscr, I'll have to check what we do with weightings there. Removing axes will necessarily destroy weightings, unless there's a way to regenerate the lower-dimensional versions of them. This possibility exists in DiscreteLp, but not in TensorSpace. Later on, when we implement per-axis weighting, that will be preserved, too.

@kohr-h

From: https://github.com/odlgroup/odl/pull/1088#discussion_r148410042

[kohr-h]

Yeah, that's one of those cases where earlier mistakes come back to bite us. I still think that "no weighting" is a broken concept and I don't regret removing it :wink:

But the immediate issue is that RayTransform tries to be intelligent about how to generate its range from its domain and the geometry. As it turns out, it tries to be a bit too clever. I would suggest to remove that cleverness and instead factor out the code for generating the ray transform range into a function ray_trafo_range that can take a weighting option. Since RayTransform already allows custom range, nothing else needs to be changed.