odlgroup / odl

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

Dynamic process with time derivatives #1261

Open MJLagerwerf opened 6 years ago

MJLagerwerf commented 6 years ago

How would you implement a 2D parallel beam CT dynamic process with time derivatives?

My first idea was just using the diagonal operator. This makes the forward operator easy to implement, but this means that the reconstructed object is an element of a product space. In which case you cannot easily define the partial time derivative.

# Number of time steps
T = 10
# Object `space`
obj_space = odl.uniform_discr([-20, -20], [20, 20], [128, 128])
# Data space
angles = odl.uniform_partition(0, np.pi, 180)
det_part = odl.uniform_partition(-40, 40, 256)
# Geometry
data_geom = odl.tomo.Parallel2dGeometry(angles, det_part)
# Radon transform and adjoint
FwP = odl.tomo.RayTransform(obj_space, data_geom)
BwP = odl.tomo.RayBackProjection(obj_space, data_geom)
# Radon transform for all time steps
K = odl.DiagonalOperator(FwP, T)

To make sure you can get the partial time derivative, you could define a dynamic space:

dyn_space = odl.uniform_discr([-20, -20, 0], [20, 20, T], [128, 128, T])

However, the forward operator K does not accept elements of dyn_space.

Is there example code for dynamic processes? If not, how would this best be handled?

kohr-h commented 6 years ago

AFAIK there's no "built-in" way to do it, but to me the dyn_space approach seems more correct in the sense that it allows you to also define time steps that are not 1, and to use all the machinery that exists in ODL, e.g. with regard to derivatives as you write above.

However, to make that scenario more user-friendly, there should likely be a conversion utility from spaces like rn((2, 3)) to rn(3) ** 2, since some operators should use the first representation, some the second.

I created a gist that implements a basic version of these conversions. Depending on your use case, you might want to wrap them into operators as well. Note that some some results are not as in the docstrings due to a bug in the handling of powers like rn(4) ** (2, 3), but that doesn't matter for your case.

And btw, for simplicity you might want to add the time dimension as first one, not last one.