feat: :sparkles: Iterative Reconstruction Methods (#12)

[WYVERN2742]

Implement iterative reconstruction methods available in CIL. Also allow for selecting arbitrary regularisation methods present in CCPi-Regularisation-Toolkit for iterative methods.

• GD
• [x] CGLS
• [x] SIRT
• ISTA Use FISTA instead
• [x] FISTA
• PDHG Not implemented, take a look at #37
• LADMM Not implemented, take a look at #37
• SPDHG Not implemented, take a look at #37
• Tikhoniv Block Framework
  • [x] Projection Operator
  • [x] Identity Operator
  • [x] Gradient Operator
• Constraint Functions
  • [x] IndicatorBox
  • [x] Total Variation
  • [x] [ccpi-reg] Fast Gradient Projection Total Variation (FGP_TV)
  • [x] [ccpi-reg] Total Generalised Variation (TGV)
  • [ccpi-reg] Fast Gradient Projection Discrete Total Variation (FGP_dTV) Requires a reference image, use FGP_TV
  • [ccpi-reg] Total Nuclear Variation (TNV) 2D only

Bugs:

• [x] Errors when combining block frameworks and constraint functions
  • Error was actually due to using GradientOperator with SIRT
  • [x] Error when using a float weight for LeastSquares
    • Has been reported to CIL and confirmed as a bug
    • Not a bug, c is the float weight paramater, weight is the paramater for DataContainers

[WYVERN2742]

Should note: support for arbitrary Tikhonov regularisers is present, currently three methods are supported:

• 'Projection' - just the standalone projection operator, with no regularisation
• 'Identity' - penalise the solution compared to the identity matrix based on its solution
• 'Gradient' - encourage similar values of neighbor values with a finite differences gradient

[WYVERN2742]

As noted in the previous commit, there is an issue when using both constraints and Regularisers together, I am unsure if this is something they support, although there are some papers of using both methods together.

[WYVERN2742]

Thanks to Jakob Sauer Jørgensen, SIRT turns out to have a non-negativity constraint inside the algorithm which caused it to fail with GradientOperator. Using IdentityOperator instead works fine, but does go to show that not all operation components will work well even if they support a group of operations.

There is a good paper on implementation from Jakob available here

There was a brief mention from multiple people that SIRT does not support TotalVariation, however it is in CIL's docs as a constraint example, so it may work in practical but not in theory? More information is needed.

[WYVERN2742]

recon-select.webm

[WYVERN2742]

Table of reconstruction and different methods all using 20 iterations for the reconstruction algorithm, along with 30 iterations for constraints/regularisers. All times are in seconds.

• X denotes failed/unsupported operations in groups.
• Blank boxes denotes non applicable groups

Method	Projection	Identity	Gradient		TotalVariation (TV)	IndicatorBox	FGP-TV	TGV		LeastSquares
FBP
FDK
CGLS	23.98	35.94	60.731
SIRT	41.058	53.569	❌		2126.372	41.415	78.672	1520.614
FISTA					2063.85	31.43	75.723	834.209		31.43

[WYVERN2742]

For future reference: PDHG/SPDHG and LADMM were not implemented due to time-complexity constraints:

export class PDHGParams implements IterativeReconstructionParams {
    readonly method = "PDHG" as const;
    quality: ReconQuality;
    iterations: number;
    f1: Proximal;
    beta: number;
    g: Proximal;
    operator: TikhonovRegulariser;
    sigma: number;
    tau: number;
}

Assumes that f = BlockOperator(f1, beta*L2NormSquared).

PDHG:
    f: BlockFunction(alpha * MixedL21Norm, beta * L2NormSquared(data))
    f: BlockFunction(alpha * MixedL21Norm, 1 * L2NormSquared(data))
    f: BlockFunction(alpha_tikhonov * L2NormSquared(), 0.5 * L2NormSquared(data))
    f: BlockFunction(L1Norm(data), alpha * MixedL21Norm, beta * MixedL21Norm)
    f: 0.5 * L2NormSquared(data)

    g: IndicatorBox(lower=0)
    g: IndicatorBox(lower=0)
    g: ZeroFunction
    g: ZeroFunction
    g: alpha * L1Norm
    g: IndicatorBox(lower=0)

    op: BlockOperator(ProjectionOperator, GradientOperator)
    op: BlockOperator(ProjectionOperator, GradientOperator)
    op: ProjectionOperator
    op: BlockOperator(ProjectionOperator, GradientOperator)

    op:
        Gradient
        Projection
        Identity

    g:
        IndicatorBox
        ZeroFunction

    f: alpha * MixedL21Norm | L2NormSquared
    f: beta * L2NormSquared

I'll open a feature issue to support PDHG