@@ -0,0 +1,2350 @@
+%22%22%22Total variation denoising of complex data using the Chambolle-Pock solver.%0A%0ASolves the optimization problem%0A%0A    min_%7Bx%7D  1/2 %7C%7Cx - g%7C%7C_2%5E2 + lam %7C%7C %7Cgrad(x)%7C %7C%7C_1%0A%0AWhere %60%60grad%60%60 the spatial gradient and %60%60g%60%60 is given noisy data.%0A%0AFor further details and a description of the solution method used, see%0A:ref:%60chambolle_pock%60 in the ODL documentation.%0A%22%22%22%0A%0Aimport numpy as np%0Aimport scipy%0Aimport odl%0A%0A# Read test image: use only every second pixel, convert integer to float,%0A# and rotate to get the image upright%0Adtype = np.complex64%0Aimage = np.rot90(scipy.misc.ascent()%5B::1, ::1%5D, 3).astype(dtype)%0Aimage = image + 1j*image.T%0Ashape = image.shape%0A%0A# Rescale max to 1%0Aimage /= image.real.max()%0A%0A# Discretized spaces%0Aspace = odl.uniform_discr(%5B0, 0%5D, shape, shape, dtype=dtype)%0A%0A# Original image%0Aorig = space.element(image)%0A%0A# Add noise%0Aimage += (np.random.normal(0, 0.05, shape) +%0A          1j*np.random.normal(0, 0.05, shape))%0A%0A# Data of noisy image%0Anoisy = space.element(image)%0A%0A# Gradient operator%0Agradient = odl.Gradient(space, method='forward')%0A%0A# Matrix of operators%0Aop = odl.BroadcastOperator(odl.IdentityOperator(space), gradient)%0A%0A# Set up the functionals%0A%0A# l2-squared data matching%0Al2_norm = odl.solvers.L2NormSquared(space).translated(noisy)%0A%0A# Isotropic TV-regularization: l1-norm of grad(x)%0Al1_norm = 0.15 * odl.solvers.L1Norm(gradient.range)%0A%0A# Make separable sum of functionals, order must correspond to the operator K%0Af = odl.solvers.SeparableSum(l2_norm, l1_norm)%0A%0Ag = odl.solvers.ZeroFunctional(op.domain)%0A%0A# --- Select solver parameters and solve using Chambolle-Pock --- #%0A%0A%0A# Estimated operator norm, add 10 percent to ensure %7C%7CK%7C%7C_2%5E2 * sigma * tau %3C 1%0Aop_norm = 1.1 * odl.power_method_opnorm(op, xstart=noisy)%0A%0Aniter = 200  # Number of iterations%0Atau = 1.0 / op_norm  # Step size for the primal variable%0Asigma = 1.0 / op_norm  # Step size for the dual variable%0A%0A# Optional: pass callback objects to solver%0Acallback = (odl.solvers.CallbackPrintIteration() &%0A            odl.solvers.CallbackShow(step=5))%0A%0A# Starting point%0Ax = op.domain.zero()%0A%0A# Run algorithm (and display intermediates)%0Aodl.solvers.chambolle_pock_solver(%0A    x, f, g, op, tau=tau, sigma=sigma, niter=niter, callback=callback)%0A%0A# Display images%0Aorig.show(title='original image')%0Anoisy.show(title='noisy image')%0Ax.show(title='reconstruction', force_show=True)%0A

add example denoising a complex-valued image