Prox for OperatorCompositionFunction

[epapoutsellis]

Description

Example Usage

from cil.optimisation.functions import OperatorCompositionFunction, L1Sparsity, L1Norm, LeastSquares, MixedL21Norm, IndicatorBox, TotalVariation
from cil.framework import ImageGeometry
from cil.optimisation.algorithms import PD3O, FISTA
from cil.optimisation.operators import WaveletOperator, IdentityOperator, GradientOperator
import numpy as np
from skimage.data import shepp_logan_phantom
from skimage.transform import rescale
from cil.framework import AcquisitionGeometry
from cil.plugins.astra import ProjectionOperator
import matplotlib.pyplot as plt

shepp = shepp_logan_phantom()
shepp = rescale(shepp, scale=0.4, mode='reflect', channel_axis=None)

angles = np.linspace(0.0, 180.0, 300, endpoint=False)
ag = AcquisitionGeometry.create_Parallel2D()\
    .set_angles(angles, angle_unit="degree")\
    .set_panel(shepp.shape[0])\
    .set_labels(['angle', 'horizontal'])
ig = ag.get_ImageGeometry()
shepp_cil = ig.allocate()
shepp_cil.fill(shepp)
A = ProjectionOperator(ig, ag, device="cpu")
sino = A.direct(shepp_cil)
plt.imshow(sino.array)

### fista proximal g is exact
alpha = 0.01
F = LeastSquares(A=A, b=sino, c=0.5)
step_size = 1./F.L
W = WaveletOperator(ig, level=0) # orthogonal_scalar=1. by default (kwargs)
f1 = alpha * L1Norm()
G_exact_1 = OperatorCompositionFunction(f1, W)

fista_exact_1 = FISTA(initial=ig.allocate(), f = F, g=G_exact_1,  update_objective_interval=50, step_size=step_size)
fista_exact_1.run(500,verbose=1)

# L1Sparsity
G_exact_2 = alpha * L1Sparsity(W)
fista_exact_2 = FISTA(initial=ig.allocate(), f = F, g=G_exact_2,  update_objective_interval=50, step_size=step_size)
fista_exact_2.run(500,verbose=1)

np.testing.assert_allclose(fista_exact_2.solution.array, fista_exact_1.solution.array, atol=1e-4)


# other problems TV + Wavelet Reg see "Efficient MR Image Reconstruction for Compressed MR Imaging"

# PD3O: Smooth term (fidelity), proximal term ( composition with wavelet), composite term (MixedL21Norm with Grad)
beta = 0.05
Grad = GradientOperator(ig)
pd3O_1 = PD3O(initial=ig.allocate(), f=F, g=G_exact_1, h = beta*MixedL21Norm(), operator=Grad, update_objective_interval=50)
pd3O_1.run(500,verbose=1)

# PD3O: Smooth term (fidelity), proximal term ( inexact TV), composite term (L1Norm with WaveletOperator)
TV = beta * TotalVariation(max_iteration=100, warm_start=True)
pd3O_2 = PD3O(initial=ig.allocate(), f=F, g=TV, h = f1, operator=W, update_objective_interval=50)
pd3O_2.run(500,verbose=1)

np.testing.assert_allclose(pd3O_1.solution.array, pd3O_2.solution.array, atol=1e-4)

❤️ Thanks for your contribution!