dTV implementation

[epapoutsellis]

I am trying to add a unittest that compares the solution from FGP_dTV vs a solution coming from cxvpy.

The input data and the reference data are:

Case 1
The parameters for the FGD_dTV are:

alpha  = 10.0
eta = 0.01
pars = {'algorithm' : FGP_dTV, \
        'input' : data,\
        'refdata' : reference,\
        'regularisation_parameter': alpha, \
        'number_of_iterations' : 50000 ,\
        'tolerance_constant':0.0,\
        'eta_const':eta,\
        'methodTV': 0 ,\
        'nonneg': 0}

and the solutions pass the unittest up to two decimal points.

Case 2
The parameters for the FGD_dTV are:

alpha  = 10.0
eta = 1.0
pars = {'algorithm' : FGP_dTV, \
        'input' : data,\
        'refdata' : reference,\
        'regularisation_parameter': alpha, \
        'number_of_iterations' : 50000 ,\
        'tolerance_constant':0.0,\
        'eta_const':eta,\
        'methodTV': 0 ,\
        'nonneg': 0}

Now, I increased eta and the solution coming from cvxpy is:

This is "a constant" image and it should be as the reference image has no influence for this value of eta.

https://github.com/vais-ral/CCPi-Regularisation-Toolkit/blob/413c6001003c6f1272aeb43152654baaf0c8a423/src/Core/regularisers_CPU/FGP_dTV_core.c#L202-L205
Therefore B_x,B_y are small, and R1,R2 do not change
https://github.com/vais-ral/CCPi-Regularisation-Toolkit/blob/413c6001003c6f1272aeb43152654baaf0c8a423/src/Core/regularisers_CPU/FGP_dTV_core.c#L219-L220

In fact it behaves like a TV reconstruction and the solution should converge to the mean value of data, i.e., np.mean(data) = 0.171875

However, the solution coming from FGP_dTV is:

Am I doing something wrong when I set up the FGD_dTV parameters?

[epapoutsellis]

To continue, I implemented a dTV reconstruction using the PDHG algorithm from CIL:

ig = ImageGeometry(N,M)
data_cil = ig.allocate()
data_cil.fill(data)

g = L2NormSquared(b=data_cil)

ref_cil = ig.allocate()
ref_cil.fill(reference)


DY = FiniteDifferenceOperator(ig, direction=1)
DX = FiniteDifferenceOperator(ig, direction=0)

Grad = BlockOperator(DY, DX)
grad_ref = Grad.direct(ref_cil)
denom = (eta**2 + grad_ref.pnorm(2)**2).sqrt()
xi = grad_ref/denom

A1 = DY - CompositionOperator(DiagonalOperator(xi[0]**2),DY) - CompositionOperator(DiagonalOperator(xi[0]*xi[1]),DX)
A2 = DX - CompositionOperator(DiagonalOperator(xi[0]*xi[1]),DY) - CompositionOperator(DiagonalOperator(xi[1]**2),DX)

operator = BlockOperator(A1, A2)

f = alpha * MixedL21Norm()

pdhg = PDHG(f = f, g = g, operator = operator, 
            max_iteration=500, update_objective_interval = 100)
pdhg.run(verbose=2)

In the figure below, there are 3 reconstructions now: FGP_dTV from the regularisation toolkit, CVXpy and CIL with PDHG.

In order for the FGP_dTV implementation to match the results from CVXpy or PDHG with cil is to use

# set parameters
pars = {'algorithm' : FGP_dTV, \
        'input' : data_cil.array,\
        'refdata' : ref_cil.array,\
        'regularisation_parameter':alpha, \
        'number_of_iterations' :5000 ,\
        'tolerance_constant':0,\
        'eta_const':100*eta,\
        'methodTV': 0 ,\
        'nonneg': 0}

So now eta is not 1.0 but 100.0 and indeed converges to the mean value of the data.