convex optimisation with complex unknowns (e.g. compressed sensing MRI)

[mehrhardt]

Solving complex problems with variational regularization is tricky and some special care needs to be given for this to make sense. I know that this is highly related to other issues raised #1328 #1333 #590 but it seems to be not properly solved (yet).

The most standard approach is to identify C^n with R^2n and define everything on R^2n, see e.g. https://www.springer.com/gp/book/9783030014575. For this one needs to define an isomorphism between these spaces. I tried to accomplish this in ODL with the following:

import odl
X = odl.cn(2)
Re = odl.RealPart(X)
Im = odl.ImagPart(X)
j = odl.BroadcastOperator(Re, Im)

which has a wrong adjoint

j.range
Out[26]: ProductSpace(rn(2), 2)

j.adjoint.domain
Out[27]: ProductSpace(cn(2), 2)

Is there another way to accomplish this in ODL?

There is an example in ODL which seems to do exactly what I want: odl/examples/solvers/douglas_rachford_pd_mri.py

On the first glance, two potential problems are that it treats images as real and the data_fit results in complex values (with zero imaginary part though). I modified the example to have the phantom purely imaginary (see below), but this resulted in weird artefacts. Thus, I guess ODL does implicitly something different (not quite sure what though). With the identification above, there should be no difference if the image was purely real or purely imaginary.

Note also that the example uses the indicator with box constraints on a complex space which is not defined. With the identification above, this should act on real and imaginary part separately but the code results in reconstructions with negative imaginary part.

"""Total variation MRI inversion using the Douglas-Rachford solver.

Solves the optimization problem

    min_{0 <= x <= 1} ||Ax - g||_2^2 + lam || |grad(x)| ||_1

where ``A`` is a simplified MRI imaging operator, ``grad`` is the spatial
gradient and ``g`` the given noisy data.
"""

import numpy as np
import odl

# Parameters
n = 256
subsampling = 0.5  # propotion of data to use
lam = 0.003

# Create a space
space = odl.uniform_discr([0, 0], [n, n], [n, n], dtype='complex')

# Create MRI operator. First fourier transform, then subsample
ft = odl.trafos.FourierTransform(space)
sampling_points = np.random.rand(*ft.range.shape) < subsampling
sampling_mask = ft.range.element(sampling_points)
mri_op = sampling_mask * ft

# Create noisy MRI data
phantom = 1j * odl.phantom.shepp_logan(space, modified=True)
noisy_data = mri_op(phantom) + odl.phantom.white_noise(mri_op.range) * 0.1
phantom.show('Phantom')
noisy_data.show('Noisy MRI Data')

# Gradient for TV regularization
gradient = odl.Gradient(space)

# Assemble all operators
lin_ops = [mri_op, gradient]

# Create functionals as needed
g = [odl.solvers.L2Norm(mri_op.range).translated(noisy_data),
     lam * odl.solvers.L1Norm(gradient.range)]
f = odl.solvers.IndicatorBox(space, 0, 1)

# Solve
x = mri_op.domain.zero()
callback = (odl.solvers.CallbackShow(step=5, clim=[-.2, 1.2]) &
            odl.solvers.CallbackPrintIteration())
odl.solvers.douglas_rachford_pd(x, f, g, lin_ops,
                                tau=2.0, sigma=[1.0, 0.1],
                                niter=500, callback=callback)

x.show('Douglas-Rachford Result')
ft.inverse(noisy_data).show('Fourier Inversion Result', force_show=True)

[mehrhardt]

@adler-j @aringh @kohr-h Anyone of you have an idea what is going wrong here and how to fix it?

[kohr-h]

Hi @mehrhardt. I just tried the code on current master, and it seems to work fine, i.e., converge to a reconstruction with small real part and Shepp-Logan imaginary part, no weird artefacts.

What may have solved the issue is the recently merged #1513. Could you try the example again with that fix?

[kohr-h]

In general, I agree that we should try to adopt the "C=R^2" interpretation of complex-valued problems throughout. We already do that partly by considering some operators as linear which are only linear in the above sense.

Related issue: #1328
WIP PR: #1333

[mehrhardt]

Hi @mehrhardt. I just tried the code on current master, and it seems to work fine, i.e., converge to a reconstruction with small real part and Shepp-Logan imaginary part, no weird artefacts.

What may have solved the issue is the recently merged #1513. Could you try the example again with that fix?

I can confirm that the example now works fine. The artefacts came from the clipping for plotting "clim=[-.2, 1.2]".

Still it uses the indicator of [0, 1] which I believe is done on the real part only. This is not obvious (or "correct") from its definition f = odl.solvers.IndicatorBox(space, 0, 1).

[mehrhardt]

In general, I agree that we should try to adopt the "C=R^2" interpretation of complex-valued problems throughout. We already do that partly by considering some operators as linear which are only linear in the above sense.

Is there a way to make this connection explicit in ODL? I.e. how could one define the bijection j : C -> R^2 in ODL? (my attempt at the top certainly did not work ...)

[mehrhardt]

@kohr-h, I have found a solution for MRI that is satisfactory, see example below. It requires a new definition of the Fourier transform as an operator from R^2n to R^2n. I could only make this work for the Discrete Fourier transform, for the "normal" Fourier transform the adjoint was off by a non-constant factor. Some parts in there (e.g. the gradient definition) are quite cumbersome, perhaps you know a better way to do this in ODL?

"""Total variation MRI inversion using the Douglas-Rachford solver.

Solves the optimization problem

    min_{0 <= x <= 1} ||Ax - g||_2^2 + lam || |grad(x)| ||_1

where ``A`` is a simplified MRI imaging operator, ``grad`` is the spatial
gradient and ``g`` the given noisy data.
"""

import numpy as np
import odl


class RealFourierTransform(odl.Operator):
    
    def __init__(self, domain):
        """TBC
        
        Parameters
        ----------
        TBC
        
        Examples
        --------
        >>> import odl
        >>> import myOperators
        >>> X = odl.uniform_discr(0, 1, 10) ** 2
        >>> F = myOperators.RealFourierTransform(X)
        >>> x = X.one()
        >>> y = F(x)
        """
        domain_complex = domain[0].complex_space
        self.fourier = odl.trafos.DiscreteFourierTransform(domain_complex)
        
        range = self.fourier.range.real_space ** 2
        
        super(RealFourierTransform, self).__init__(
                domain=domain, range=range, linear=True)
    
    def _call(self, x, out):
        Fx = self.fourier(x[0].asarray() + 1j * x[1].asarray())
        out[0][:] = np.real(Fx)
        out[1][:] = np.imag(Fx)
        
        out *= self.domain[0].cell_volume
                            
    @property
    def adjoint(self):
        op = self
        
        class RealFourierTransformAdjoint(odl.Operator):
    
            def __init__(self, op):        
                """TBC
                
                Parameters
                ----------
                TBC
                
                Examples
                --------
                >>> import odl
                >>> import myOperators
                >>> X = odl.uniform_discr(0, 2, 10) ** 2
                >>> A = myOperators.RealFourierTransform(X)
                >>> x = odl.phantom.white_noise(A.domain)
                >>> y = odl.phantom.white_noise(A.range)
                >>> t1 = A(x).inner(y)
                >>> t2 = x.inner(A.adjoint(y))
                >>> t1 / t2
                
                >>> import odl
                >>> import myOperators
                >>> X = odl.uniform_discr([-1, -1], [2, 1], [10, 30]) ** 2
                >>> A = myOperators.RealFourierTransform(X)
                >>> x = odl.phantom.white_noise(A.domain)
                >>> y = odl.phantom.white_noise(A.range)
                >>> t1 = A(x).inner(y)
                >>> t2 = x.inner(A.adjoint(y))
                >>> t1 / t2
                """
                self.op = op
                
                super(RealFourierTransformAdjoint, self).__init__(
                        domain=op.range, range=op.domain, linear=True)
            
            def _call(self, x, out):
                y = x[0].asarray() + 1j * x[1].asarray()
                Fadjy = self.op.fourier.adjoint(y)
                out[0][:] = np.real(Fadjy)
                out[1][:] = np.imag(Fadjy)
                
                out *= self.op.fourier.domain.size
        
            @property
            def adjoint(self):
                return op
        
        return RealFourierTransformAdjoint(op)
    
    

# Parameters
n = 256
subsampling = 0.5  # propotion of data to use
lam = 0.4

#%%
# Create a space
complex_space = odl.uniform_discr([-1, -1], [1, 1], [n, n], dtype='complex')
space = complex_space.real_space ** 2

M = odl.PointwiseNorm(space)

# Create MRI operator. First fourier transform, then subsample
ft = RealFourierTransform(space)
sampling_points = np.random.rand(*ft.range.shape) < subsampling
sampling_mask = ft.range.element(sampling_points)
mri_op = 1/2 * sampling_mask * ft

#%%
# Create noisy MRI data
# Ground truth image
phase = complex_space.element(lambda x: np.exp(1j * 0.1 * (np.sin(3 * x[0] ** 2) + np.cos(5 * x[1] ** 3))))
gt = odl.phantom.shepp_logan(complex_space, modified=True) * phase
phantom = space.element()
phantom[0] = np.real(gt)
phantom[1] = np.imag(gt)

data = mri_op(phantom)
noisy_data = mri_op(phantom) + odl.phantom.white_noise(mri_op.range, stddev=0.0001)
#phantom.show('Phantom')
#noisy_data.show('Noisy MRI Data')

M(phantom).show('Phantom')
M(mri_op.adjoint(data)).show('Linear recon', clim=[0,1])

#%%
# Gradient for TV regularization
pd_basic = [odl.discr.diff_ops.PartialDerivative(space[0], i) for i in range(2)]
cp_basic = [odl.operator.ComponentProjection(space, i) for i in range(2)]
gradient = odl.BroadcastOperator(*[pd * cp for cp in cp_basic for pd in pd_basic])

# Assemble all operators
lin_ops = [mri_op, 1 / (np.sqrt(2) * n) * gradient]

# Create functionals as needed
g = [odl.solvers.L2Norm(mri_op.range).translated(noisy_data),
     lam * odl.solvers.GroupL1Norm(gradient.range)]
f = odl.solvers.ZeroFunctional(space)

# Solve
x = mri_op.domain.zero()
callback = (odl.solvers.CallbackShow(step=50) &
            odl.solvers.CallbackPrintIteration())
odl.solvers.douglas_rachford_pd(x, f, g, lin_ops,
                                tau=2.0, sigma=[1.0, 1.0],
                                niter=500, callback=callback)

M(x).show('Douglas-Rachford Result', clim=[0,1])