mehrhardt
commented
5 years ago @adler-j @aringh @kohr-h Anyone of you have an idea what is going wrong here and how to fix it?
Open mehrhardt opened 5 years ago
mehrhardt
commented
5 years ago @adler-j @aringh @kohr-h Anyone of you have an idea what is going wrong here and how to fix it?
kohr-h
commented
5 years ago 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
commented
5 years ago 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
commented
5 years ago 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
commented
5 years ago 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
commented
5 years ago @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])
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 https://github.com/odlgroup/odl/issues/1328 https://github.com/odlgroup/odl/pull/1333 https://github.com/odlgroup/odl/issues/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:
which has a wrong adjoint
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.pyOn the first glance, two potential problems are that it treats images as
realand thedata_fitresults 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.