odlgroup / odl

Operator Discretization Library https://odlgroup.github.io/odl/
Mozilla Public License 2.0
374 stars 105 forks source link

The problem about indexing parallel2Dgeometry #1625

Open dlutgy opened 1 year ago

dlutgy commented 1 year ago

I meet an inconsistent in the following code.

import odl
import numpy as np
X = odl.uniform_discr(min_pt=[-1, -1], max_pt=[1, 1], shape=[150, 150])
x = X.one()
ang_num = 2
det_num = 3
angle_partition = odl.uniform_partition(0, np.pi, ang_num)
detector_partition = odl.uniform_partition(-3, 3, det_num)
geometry = odl.tomo.Parallel2dGeometry(angle_partition, detector_partition)
op = odl.tomo.RayTransform(X, geometry)
print(op(x))

### split the geometry
geo1 = [geometry[i,j] for i in range(ang_num) for j in range(det_num)]
ops = [odl.tomo.RayTransform(X, geo) for geo in geo1]
list_opx = [np.array(subop(x)) for subop in ops]
print(np.array(list_opx))

The outcome is

[[ 0.        ,  1.42327523,  1.42327273,  0.        ],
 [ 0.        ,  1.99999833,  1.99999833,  0.        ],
 [ 0.        ,  1.42328048,  1.42327571,  0.        ]]
[[[ 2.30940294]]

 [[ 2.30940294]]

 [[ 2.30940294]]

 [[ 2.30940294]]

 [[ 1.99999833]]

 [[ 1.99999833]]

 [[ 1.99999833]]

 [[ 1.99999833]]

 [[ 2.30940175]]

 [[ 2.30940175]]

 [[ 2.30940175]]

 [[ 2.30940175]]]

Why is the result of the sinograms different when you use the whole geometry (i.e. 'geometry' in the code) than when you use the separate geometry(i.e. 'geo1' in the code) ? Thanks for your attention, I am looking forward to hear some better suggestions from you.

dlutgy commented 1 year ago

Sorry, there is some typo. The value of ang_num and det_num should be 3 and 4, respectively. Namely,

ang_num = 3
det_num = 4