deconvolution.py

import numpy as np
import scipy as sp
import odl
from scipy import signal


class Convolve(odl.Operator):
    def __init__(self, space, kernel):
        super(Convolve, self).__init__(domain=space, range=space,
                                       linear=True)
        self.kernel = kernel

    def _call(self, x):
        return signal.fftconvolve(x, kernel, mode='same')

    @property
    def adjoint(self):
        return self


space = odl.uniform_discr([-1, -1], [1, 1], [100, 100])
kernel = space.element(lambda x: np.exp(-(x[0]**2 + x[1]**2) / 0.1 **2))
kernel.show('kernel')

conv = Convolve(space, kernel)
x = odl.phantom.forbild(space)
x.show('x')

data = conv(x)
data.show('data')

grad = odl.Gradient(space)
l1_norm = 0.0001 * odl.solvers.L1Norm(grad.range)

data_discr = odl.solvers.L2NormSquared(space).translated(data)

x0 = space.zero()

n1 = grad.norm(estimate=True)
n2 = conv.norm(estimate=True)

odl.solvers.douglas_rachford_pd(x0, f=odl.solvers.ZeroFunctional(space),
                                g=[l1_norm, data_discr],
                                L=[grad, conv],
                                tau=1.0, sigma=[1/n1**2, 1/n2**2], niter=1000,
                                callback=odl.solvers.CallbackShow(step=10))