Geometry `dtype` mistmatch causes error when computing `norm` of ASTRA projection operator

[tommheik]

Description

If the ImageGeometry and AcquisitionGeometry have different data types (due to data for example), then computing the norm of the ProjectionOperator does not work (at least with ASTRA and cpu). The operator itself works fine so likely the in-place computations are the issue (although in-place computations also work...).

Here is a minimal works example:

import numpy as np
import cil, sys
print(cil.version.version, cil.version.commit_hash, sys.version, sys.platform)
from cil.plugins.astra import ProjectionOperator
from cil.framework import AcquisitionGeometry, AcquisitionData

print(f"Using CIL version {cil.__version__}")

# Create acquisition geometry
ag = AcquisitionGeometry.create_Cone2D([0,-100], [0,200], units='mm',detector_direction_x=[1, 0])
ag.set_panel(120, pixel_size=0.1, origin='bottom-left')
ag.set_angles(np.linspace(0., 360., 120, endpoint=False))
ag.set_labels(('angle','horizontal'))

ig = ag.get_ImageGeometry()

A1 = ProjectionOperator(ig, ag, 'cpu')
print(f"This works fine, the operator norm is {A1.norm()}")
print(f"By default: {ig.dtype = }, {ag.dtype = }")

sinogram64 = np.zeros((120, 120), dtype='float64')

# Setup float64 data (this changes `dtype` of ag!)
data = AcquisitionData(sinogram64, deep_copy=False, geometry=ag)

print(f"Now: {ig.dtype = }, {ag.dtype = }")

# Creating operator works fine
A2 = ProjectionOperator(ig, ag, 'cpu')

# Using the operator works fine
bp1 = A1.adjoint(data)
bp2 = A2.adjoint(data)

# Using the operator in place also works
out1 = data.copy()
A1.adjoint(out1, out=out1)
out2 = data.copy()
A2.adjoint(out2, out=out2)

# This fails
print(f"This no longer works: the operator norm is {A2.norm()}")

Environment

import cil, sys
print(cil.version.version, cil.version.commit_hash, sys.version, sys.platform)

Output: 24.2.0 None 3.12.5 | packaged by conda-forge | (main, Aug 8 2024, 18:36:51) [GCC 12.4.0] linux

[MargaretDuff]

Thanks @tommheik

Some quick investigation leads to a more minimal example:

import numpy as np
import cil, sys
print(cil.version.version, cil.version.commit_hash, sys.version, sys.platform)
from cil.plugins.astra import ProjectionOperator
from cil.framework import AcquisitionGeometry, AcquisitionData

print(f"Using CIL version {cil.__version__}")

# Create acquisition geometry
ag = AcquisitionGeometry.create_Cone2D([0,-100], [0,200], units='mm',detector_direction_x=[1, 0])
ag.set_panel(120, pixel_size=0.1, origin='bottom-left')
ag.set_angles(np.linspace(0., 360., 120, endpoint=False))
ag.set_labels(('angle','horizontal'))

ig = ag.get_ImageGeometry()

A1 = ProjectionOperator(ig, ag, 'cpu')

sinogram64 = np.zeros((120, 120), dtype='float64')

# Setup float64 data (this changes `dtype` of ag!)
data = AcquisitionData(sinogram64, deep_copy=False, geometry=ag)

#%% This now fails 
x=A1.domain_geometry().allocate('random')
y_tmp = A1.range_geometry().allocate()
A1.direct(x, out=y_tmp)

There are two issues, the first that making the acquisition data here data = AcquisitionData(sinogram64, deep_copy=False, geometry=ag) silently changes the acquisition geometry dtype. This could be checked in the Acquisition data init, similar to where we check the shape: https://github.com/TomographicImaging/CIL/blob/69391bb87a05a85ea69ee8c3e43501b70036d302/Wrappers/Python/cil/framework/acquisition_data.py#L71-L72

The part that causes the fail is here https://github.com/TomographicImaging/CIL/blob/69391bb87a05a85ea69ee8c3e43501b70036d302/Wrappers/Python/cil/framework/processors.py#L131-L135

Will chat about this with the CIL team for ideas

[paskino]

Here issue and PR where we introduced this https://github.com/TomographicImaging/CIL/issues/1019