DiagonalOperator as ProductSpace for dynamic CT

[gesa-sa]

I'm trying to implement an operator with ODL.tomo for dynamic CT with parallel beam. At each time step the object should be scanned from different rotating angles distributed in [0,180°]. Since the geometry is different for each time step, also the range would be different for each RayTransform. In ODL 0.7.0 it worked to set the 'range' parameter for ray_trafo to a common_range, such that the resulting DiagonalOperator operated on a PowerSpace. In ODL 0.8.1 I can still set the 'range' parameter but I think it is ignored and overwritten by the range it gets from the geometry. Then I don't get a PoweSpace but only a ProductSpace. Is there a way to set the range parameter to a common range?

  detector_partition = odl.uniform_partition(-1.25 * space_range, 1.25 * space_range, det_shape)
  
  space = odl.uniform_discr([-space_range, -space_range], [space_range, space_range], (phantom.size(3), phantom.size(4)), dtype='float32', weighting=1.0)
  
  operators = []
  
  common_range = odl.uniform_discr([0, 0], [num_angles, det_shape], (num_angles, det_shape), dtype='float32')
  
  for t in range(time_steps):
      angles_all = np.linspace(0, np.pi, num_angles*time_steps, endpoint=False)
      angles = angles_all[t::time_steps]
  
      angle_partition = odl.nonuniform_partition(angles)
      geometry = odl.tomo.Parallel2dGeometry(angle_partition, detector_partition)
  
      ray_trafo = odl.tomo.RayTransform(space, geometry, impl='astra_cpu', range=common_range)
      operators.append(ray_trafo)
  
  
  fwd_op_odl = odl.DiagonalOperator(*operators)

[Emvlt]

Hello, Can you clarify why you want an operator with a common range?

Reading your code, it appears that you select a subset of angles_all for each t, resulting in time_steps different geometries. As such, when you project the volume, you will have a different number of angles for each t, but the domain of each operator will remain the same.

Also, a PowerSpace is a ProductSpace in which all spaces are identical. In your case, the spaces are not identical so i'm not sure that you need a PowerSpace.

[gesa-sa]

Hello,

thanks for your answer! I have different angles for each time-step but always the same number of angles for each t, so the dimension of all the spaces in the ProductSpace is the same. If the space is a PowerSpace, I can use np.asarry(result).astype(dtype, copy=False) to convert the result to a numpy array where result is the output of the forward operator. 'asarray' only works if space.is_power_space is True. At the moment I use a loop over the spaces in my ProductSpace such that I can use 'asarray' for each space seperately which is probably less efficient than using 'asarray' on the whole ProductSpace.

[leftaroundabout]

Ok... I see how extracting to a single array is useful, but going through a PowerSpace seems conceptually the wrong way of getting there, would you agree? Instead, ProductSpace should allow converting to a single array whenever all the subspaces yield arrays of the same space, regardless of whether they have the same mathematical meaning. Would that solve your problem?

[gesa-sa]

Yes, that's right! The possibility to convert to a single array would solve my problem. But maybe it would also be helpful if the range parameter couldn't be specified manually in the construction of a Ray_Trafo if it is chosen automatically anyway. Sorry for the unclear description...

[Emvlt]

I understand the issue better. What happens is that although the angles are different for each time step, there is in your case always the same number of them, so you'd expect to be able, let's say, to put all different sinograms in a array of size [timesteps, n_angles, det_shape]. Then I agree that @leftaroundabout 's solution is adapted.

As for the clarity, you'll notice that range is not a mandatory argument of the RayTransform. Although it inherits from the Operator class, the API expects a geometry argument and no range. Passing range as a kwarg will not do anything as the range of the operator is inferred by the geometry.

@leftaroundabout should we add a check to make sure that there are no remaining kwargs when instanciating a RayTransform ?