From 29dbf5aedbdad3a6d6f8eceec1becf24e8ec5431 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 15 Nov 2023 17:07:47 +0000 Subject: [PATCH 01/11] added initial TGV explicit to test --- mantidimaging/core/reconstruct/cil_recon.py | 60 +++++++++++++++++++++ 1 file changed, 60 insertions(+) diff --git a/mantidimaging/core/reconstruct/cil_recon.py b/mantidimaging/core/reconstruct/cil_recon.py index 1ebf8c9886e..00c927bca58 100644 --- a/mantidimaging/core/reconstruct/cil_recon.py +++ b/mantidimaging/core/reconstruct/cil_recon.py @@ -13,6 +13,8 @@ from cil.framework import AcquisitionData, AcquisitionGeometry, DataOrder, ImageGeometry, BlockGeometry from cil.optimisation.algorithms import PDHG, SPDHG from cil.optimisation.operators import GradientOperator, BlockOperator +from cil.optimisation.operators import SymmetrisedGradientOperator, ZeroOperator, IdentityOperator + from cil.optimisation.functions import MixedL21Norm, L2NormSquared, BlockFunction, ZeroFunction, IndicatorBox, Function from cil.plugins.astra.operators import ProjectionOperator @@ -82,6 +84,64 @@ def set_up_TV_regularisation( return (K, F, G) + + @staticmethod + def set_up_TGV_regularisation( + image_geometry: ImageGeometry, acquisition_data: AcquisitionData, + recon_params: ReconstructionParameters) -> tuple[BlockOperator, BlockFunction, Function]: + + # Forward operator + A2d = ProjectionOperator(image_geometry, acquisition_data.geometry, 'gpu') + + if recon_params.stochastic: + for partition_geometry, partition_operator in zip(acquisition_data.geometry, A2d, strict=True): + CILRecon.set_approx_norm(partition_operator, partition_geometry, image_geometry) + else: + CILRecon.set_approx_norm(A2d, acquisition_data.geometry, image_geometry) + + # Define Gradient Operator and BlockOperator + alpha = recon_params.alpha + gamma = recon_params.gamma + beta = alpha / gamma + + f2 = alpha * MixedL21Norm() + f3 = beta * MixedL21Norm() + + if recon_params.stochastic: + + # now, A2d is a BlockOperator as acquisition_data is a BlockDataContainer + fs = [] + for i, _ in enumerate(acquisition_data.geometry): + fs.append(L2NormSquared(b=acquisition_data.get_item(i))) + + F = BlockFunction(*fs, f2, f3) + + + else: + # Define BlockFunction F using the MixedL21Norm() and the L2NormSquared() + f1 = 0.5 * L2NormSquared(b=acquisition_data) + + F = BlockFunction(f1, f2, f3) + + # Define BlockOperator K + + #%% + K11 = A2d + K21 = GradientOperator(image_geometry) + K32 = SymmetrisedGradientOperator(K21.range) + K12 = ZeroOperator(K32.domain, image_geometry) + K22 = IdentityOperator(K21.range) + K31 = ZeroOperator(image_geometry, K32.range) + K = BlockOperator(K11, K12, K21, -K22, K31, K32, shape=(3,2) ) + + if recon_params.non_negative: + G = IndicatorBox(lower=0) + else: + # Define Function G simply as zero + G = ZeroFunction() + + return (K, F, G) + @staticmethod def set_approx_norm(A2d: BlockOperator, acquisition_data: AcquisitionGeometry, image_geometry: ImageGeometry) -> None: From ba07f18c817d7b7b01a17130e67c3e04477648db Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 29 Nov 2023 12:24:08 +0000 Subject: [PATCH 02/11] Return a np.array for TGV case --- mantidimaging/core/reconstruct/cil_recon.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/mantidimaging/core/reconstruct/cil_recon.py b/mantidimaging/core/reconstruct/cil_recon.py index 00c927bca58..14628c3b60a 100644 --- a/mantidimaging/core/reconstruct/cil_recon.py +++ b/mantidimaging/core/reconstruct/cil_recon.py @@ -10,7 +10,7 @@ import numpy as np -from cil.framework import AcquisitionData, AcquisitionGeometry, DataOrder, ImageGeometry, BlockGeometry +from cil.framework import AcquisitionData, AcquisitionGeometry, DataOrder, ImageGeometry, BlockGeometry, BlockDataContainer from cil.optimisation.algorithms import PDHG, SPDHG from cil.optimisation.operators import GradientOperator, BlockOperator from cil.optimisation.operators import SymmetrisedGradientOperator, ZeroOperator, IdentityOperator @@ -270,6 +270,11 @@ def single_sino(sino: np.ndarray, progress.mark_complete() t1 = time.perf_counter() LOG.info(f"single_sino time: {t1-t0}s for shape {sino.shape}") + + if isinstance(algo.solution, BlockDataContainer): + # TGV case + return algo.solution[0].as_array() + return algo.solution.as_array() @staticmethod @@ -385,7 +390,11 @@ def full(images: ImageStack, force_continue=False) algo.next() - volume = algo.solution.as_array() + if isinstance(algo.solution, BlockDataContainer): + # TGV case + volume = algo.solution[0].as_array() + else: + volume = algo.solution.as_array() LOG.info(f'Reconstructed 3D volume with shape: {volume.shape}') t1 = time.perf_counter() LOG.info(f"full reconstruction time: {t1-t0}s for shape {images.data.shape}") From ed767113e5e3c56caedd8f2d422434eeffa5cbe9 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 29 Nov 2023 12:40:22 +0000 Subject: [PATCH 03/11] add selector between TV and TGV --- mantidimaging/core/reconstruct/cil_recon.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/mantidimaging/core/reconstruct/cil_recon.py b/mantidimaging/core/reconstruct/cil_recon.py index 14628c3b60a..4e8615ff457 100644 --- a/mantidimaging/core/reconstruct/cil_recon.py +++ b/mantidimaging/core/reconstruct/cil_recon.py @@ -229,7 +229,10 @@ def single_sino(sino: np.ndarray, ig = ag.get_ImageGeometry() - K, F, G = CILRecon.set_up_TV_regularisation(ig, data, recon_params) + if recon_params.regulariser == 'TV': + K, F, G = CILRecon.set_up_TV_regularisation(ig, data, recon_params) + elif recon_params.regulariser == 'TGV': + K, F, G = CILRecon.set_up_TGV_regularisation(ig, data, recon_params) max_iteration = 100000 # this should set to a sensible number as evaluating the objective is costly @@ -353,7 +356,11 @@ def full(images: ImageStack, num_subsets) ig = ag.get_ImageGeometry() - K, F, G = CILRecon.set_up_TV_regularisation(ig, data, recon_params) + + if recon_params.regulariser == 'TV': + K, F, G = CILRecon.set_up_TV_regularisation(ig, data, recon_params) + elif recon_params.regulariser == 'TGV': + K, F, G = CILRecon.set_up_TGV_regularisation(ig, data, recon_params) max_iteration = 100000 # this should set to a sensible number as evaluating the objective is costly From 27e5a86c528bbd439b25f7e6e3ad31ec5da019d9 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 30 Nov 2023 09:21:43 +0000 Subject: [PATCH 04/11] updated TGV and moved alpha to the operators --- mantidimaging/core/reconstruct/cil_recon.py | 25 ++++++++++++--------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/mantidimaging/core/reconstruct/cil_recon.py b/mantidimaging/core/reconstruct/cil_recon.py index 4e8615ff457..5f850bf1f41 100644 --- a/mantidimaging/core/reconstruct/cil_recon.py +++ b/mantidimaging/core/reconstruct/cil_recon.py @@ -104,8 +104,8 @@ def set_up_TGV_regularisation( gamma = recon_params.gamma beta = alpha / gamma - f2 = alpha * MixedL21Norm() - f3 = beta * MixedL21Norm() + f2 = MixedL21Norm() + f3 = MixedL21Norm() if recon_params.stochastic: @@ -119,20 +119,25 @@ def set_up_TGV_regularisation( else: # Define BlockFunction F using the MixedL21Norm() and the L2NormSquared() - f1 = 0.5 * L2NormSquared(b=acquisition_data) + # mathematicians like to multiply 1/2 in front of L2NormSquared. This is not necessary + # it will mean that the regularisation parameter alpha is doubled + f1 = L2NormSquared(b=acquisition_data) F = BlockFunction(f1, f2, f3) # Define BlockOperator K - #%% + # Set up the 3 operator A, Grad and Epsilon K11 = A2d - K21 = GradientOperator(image_geometry) - K32 = SymmetrisedGradientOperator(K21.range) - K12 = ZeroOperator(K32.domain, image_geometry) - K22 = IdentityOperator(K21.range) - K31 = ZeroOperator(image_geometry, K32.range) - K = BlockOperator(K11, K12, K21, -K22, K31, K32, shape=(3,2) ) + K21 = alpha * GradientOperator(K11.domain) + # https://tomographicimaging.github.io/CIL/nightly/optimisation.html#cil.optimisation.operators.SymmetrisedGradientOperator + K32 = beta * SymmetrisedGradientOperator(K21.range) + # these define the domain and range of the other operators + K12 = ZeroOperator(K32.domain, K11.range) + K22 = -alpha * IdentityOperator(domain_geometry=K21.range, range_geometry=K32.range) + K31 = ZeroOperator(K11.domain, K32.range) + + K = BlockOperator(K11, K12, K21, K22, K31, K32, shape=(3,2) ) if recon_params.non_negative: G = IndicatorBox(lower=0) From cbeece9eba2f1cdd979707ec3994708da94d242f Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 30 Nov 2023 09:23:45 +0000 Subject: [PATCH 05/11] removed unhelpful comment --- mantidimaging/core/reconstruct/cil_recon.py | 1 - 1 file changed, 1 deletion(-) diff --git a/mantidimaging/core/reconstruct/cil_recon.py b/mantidimaging/core/reconstruct/cil_recon.py index 5f850bf1f41..c346effd68b 100644 --- a/mantidimaging/core/reconstruct/cil_recon.py +++ b/mantidimaging/core/reconstruct/cil_recon.py @@ -130,7 +130,6 @@ def set_up_TGV_regularisation( # Set up the 3 operator A, Grad and Epsilon K11 = A2d K21 = alpha * GradientOperator(K11.domain) - # https://tomographicimaging.github.io/CIL/nightly/optimisation.html#cil.optimisation.operators.SymmetrisedGradientOperator K32 = beta * SymmetrisedGradientOperator(K21.range) # these define the domain and range of the other operators K12 = ZeroOperator(K32.domain, K11.range) From d8157a94a7455ba952e6f2071d4fac2b29eaf925 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 1 Dec 2023 13:42:31 +0000 Subject: [PATCH 06/11] gamma is multiplicative instead in a ratio --- mantidimaging/core/reconstruct/cil_recon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mantidimaging/core/reconstruct/cil_recon.py b/mantidimaging/core/reconstruct/cil_recon.py index c346effd68b..d093087f1e9 100644 --- a/mantidimaging/core/reconstruct/cil_recon.py +++ b/mantidimaging/core/reconstruct/cil_recon.py @@ -102,7 +102,7 @@ def set_up_TGV_regularisation( # Define Gradient Operator and BlockOperator alpha = recon_params.alpha gamma = recon_params.gamma - beta = alpha / gamma + beta = alpha * gamma f2 = MixedL21Norm() f3 = MixedL21Norm() From d62b92d23f6afca6e93e485908afa22b48be7a0d Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 5 Dec 2023 14:04:12 +0000 Subject: [PATCH 07/11] cherry pick from Sam --- mantidimaging/core/reconstruct/cil_recon.py | 20 +++++++++++-------- mantidimaging/core/utility/data_containers.py | 4 ++++ mantidimaging/gui/windows/recon/view.py | 10 ++++++++++ 3 files changed, 26 insertions(+), 8 deletions(-) diff --git a/mantidimaging/core/reconstruct/cil_recon.py b/mantidimaging/core/reconstruct/cil_recon.py index d093087f1e9..ce146b9f147 100644 --- a/mantidimaging/core/reconstruct/cil_recon.py +++ b/mantidimaging/core/reconstruct/cil_recon.py @@ -10,7 +10,8 @@ import numpy as np -from cil.framework import AcquisitionData, AcquisitionGeometry, DataOrder, ImageGeometry, BlockGeometry, BlockDataContainer +from cil.framework import (AcquisitionData, AcquisitionGeometry, DataOrder, ImageGeometry, BlockGeometry, + BlockDataContainer) from cil.optimisation.algorithms import PDHG, SPDHG from cil.optimisation.operators import GradientOperator, BlockOperator from cil.optimisation.operators import SymmetrisedGradientOperator, ZeroOperator, IdentityOperator @@ -105,17 +106,16 @@ def set_up_TGV_regularisation( beta = alpha * gamma f2 = MixedL21Norm() - f3 = MixedL21Norm() + f3 = MixedL21Norm() if recon_params.stochastic: - + # now, A2d is a BlockOperator as acquisition_data is a BlockDataContainer fs = [] for i, _ in enumerate(acquisition_data.geometry): fs.append(L2NormSquared(b=acquisition_data.get_item(i))) - - F = BlockFunction(*fs, f2, f3) + F = BlockFunction(*fs, f2, f3) else: # Define BlockFunction F using the MixedL21Norm() and the L2NormSquared() @@ -123,11 +123,11 @@ def set_up_TGV_regularisation( # it will mean that the regularisation parameter alpha is doubled f1 = L2NormSquared(b=acquisition_data) - F = BlockFunction(f1, f2, f3) + F = BlockFunction(f1, f2, f3) # Define BlockOperator K - # Set up the 3 operator A, Grad and Epsilon + # Set up the 3 operator A, Grad and Symmetrised Gradient K11 = A2d K21 = alpha * GradientOperator(K11.domain) K32 = beta * SymmetrisedGradientOperator(K21.range) @@ -136,7 +136,7 @@ def set_up_TGV_regularisation( K22 = -alpha * IdentityOperator(domain_geometry=K21.range, range_geometry=K32.range) K31 = ZeroOperator(K11.domain, K32.range) - K = BlockOperator(K11, K12, K21, K22, K31, K32, shape=(3,2) ) + K = BlockOperator(K11, K12, K21, K22, K31, K32, shape=(3, 2)) if recon_params.non_negative: G = IndicatorBox(lower=0) @@ -237,6 +237,8 @@ def single_sino(sino: np.ndarray, K, F, G = CILRecon.set_up_TV_regularisation(ig, data, recon_params) elif recon_params.regulariser == 'TGV': K, F, G = CILRecon.set_up_TGV_regularisation(ig, data, recon_params) + else: + raise ValueError(f"Regulariser must be one of 'TV', 'TGV'. Received '{recon_params.regulariser}'") max_iteration = 100000 # this should set to a sensible number as evaluating the objective is costly @@ -365,6 +367,8 @@ def full(images: ImageStack, K, F, G = CILRecon.set_up_TV_regularisation(ig, data, recon_params) elif recon_params.regulariser == 'TGV': K, F, G = CILRecon.set_up_TGV_regularisation(ig, data, recon_params) + else: + raise ValueError(f"Regulariser must be one of 'TV', 'TGV'. Received '{recon_params.regulariser}'") max_iteration = 100000 # this should set to a sensible number as evaluating the objective is costly diff --git a/mantidimaging/core/utility/data_containers.py b/mantidimaging/core/utility/data_containers.py index 95dd73ce633..64ef1566bea 100644 --- a/mantidimaging/core/utility/data_containers.py +++ b/mantidimaging/core/utility/data_containers.py @@ -103,12 +103,14 @@ class ReconstructionParameters: tilt: Degrees | None = None pixel_size: float = 0.0 alpha: float = 0.0 + gamma: float = 1.0 non_negative: bool = False max_projection_angle: float = 360.0 beam_hardening_coefs: list[float] | None = None stochastic: bool = False projections_per_subset: int = 50 regularisation_percent: int = 30 + regulariser: str = "" def to_dict(self) -> dict: return { @@ -119,9 +121,11 @@ def to_dict(self) -> dict: 'tilt': str(self.tilt), 'pixel_size': self.pixel_size, 'alpha': self.alpha, + 'gamma': self.gamma, 'stochastic': self.stochastic, 'projections_per_subset': self.projections_per_subset, 'regularisation_percent': self.regularisation_percent, + 'regulariser': self.regulariser, } diff --git a/mantidimaging/gui/windows/recon/view.py b/mantidimaging/gui/windows/recon/view.py index 7ac9195c8d9..5bacb18ce7b 100644 --- a/mantidimaging/gui/windows/recon/view.py +++ b/mantidimaging/gui/windows/recon/view.py @@ -411,6 +411,10 @@ def pixel_size(self, value: int) -> None: def alpha(self) -> float: return self.alphaSpinBox.value() + @property + def gamma(self) -> float: + return 1 + @property def non_negative(self) -> bool: return self.nonNegativeCheckBox.isChecked() @@ -427,6 +431,10 @@ def projections_per_subset(self) -> int: def regularisation_percent(self) -> int: return self.regPercentSpinBox.value() + @property + def regulariser(self) -> str: + return "TV" + @property def beam_hardening_coefs(self) -> list[float] | None: if not self.lbhc_enabled.isChecked(): @@ -447,11 +455,13 @@ def recon_params(self) -> ReconstructionParameters: tilt=Degrees(self.tilt), pixel_size=self.pixel_size, alpha=self.alpha, + gamma=self.gamma, non_negative=self.non_negative, stochastic=self.stochastic, projections_per_subset=self.projections_per_subset, max_projection_angle=self.max_proj_angle, regularisation_percent=self.regularisation_percent, + regulariser=self.regulariser, beam_hardening_coefs=self.beam_hardening_coefs) def set_table_point(self, idx, slice_idx, cor): From a1a4f871085b9615863779acf9b49ca6c653f8df Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 5 Dec 2023 14:29:33 +0000 Subject: [PATCH 08/11] Fix non-negativity constraint --- mantidimaging/core/reconstruct/cil_recon.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/mantidimaging/core/reconstruct/cil_recon.py b/mantidimaging/core/reconstruct/cil_recon.py index ce146b9f147..2e06718992a 100644 --- a/mantidimaging/core/reconstruct/cil_recon.py +++ b/mantidimaging/core/reconstruct/cil_recon.py @@ -139,7 +139,8 @@ def set_up_TGV_regularisation( K = BlockOperator(K11, K12, K21, K22, K31, K32, shape=(3, 2)) if recon_params.non_negative: - G = IndicatorBox(lower=0) + G = BlockFunction(IndicatorBox(lower=0, upper=None), ZeroFunction()) + else: # Define Function G simply as zero G = ZeroFunction() From 880fc24eda417434aef6fc1c27ca399d613ad991 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 5 Dec 2023 14:35:10 +0000 Subject: [PATCH 09/11] autoyapf complaints --- mantidimaging/core/reconstruct/cil_recon.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/mantidimaging/core/reconstruct/cil_recon.py b/mantidimaging/core/reconstruct/cil_recon.py index 2e06718992a..a6bca8eeeeb 100644 --- a/mantidimaging/core/reconstruct/cil_recon.py +++ b/mantidimaging/core/reconstruct/cil_recon.py @@ -10,7 +10,7 @@ import numpy as np -from cil.framework import (AcquisitionData, AcquisitionGeometry, DataOrder, ImageGeometry, BlockGeometry, +from cil.framework import (AcquisitionData, AcquisitionGeometry, DataOrder, ImageGeometry, BlockGeometry, BlockDataContainer) from cil.optimisation.algorithms import PDHG, SPDHG from cil.optimisation.operators import GradientOperator, BlockOperator @@ -85,7 +85,6 @@ def set_up_TV_regularisation( return (K, F, G) - @staticmethod def set_up_TGV_regularisation( image_geometry: ImageGeometry, acquisition_data: AcquisitionData, @@ -122,7 +121,7 @@ def set_up_TGV_regularisation( # mathematicians like to multiply 1/2 in front of L2NormSquared. This is not necessary # it will mean that the regularisation parameter alpha is doubled f1 = L2NormSquared(b=acquisition_data) - + F = BlockFunction(f1, f2, f3) # Define BlockOperator K From 412573a124f303428d53c3a7d285633b0e38566e Mon Sep 17 00:00:00 2001 From: Sam Tygier Date: Thu, 14 Dec 2023 11:59:14 +0000 Subject: [PATCH 10/11] Add release note --- docs/release_notes/next/dev-1985-tgv-backend | 1 + 1 file changed, 1 insertion(+) create mode 100644 docs/release_notes/next/dev-1985-tgv-backend diff --git a/docs/release_notes/next/dev-1985-tgv-backend b/docs/release_notes/next/dev-1985-tgv-backend new file mode 100644 index 00000000000..2edb8213492 --- /dev/null +++ b/docs/release_notes/next/dev-1985-tgv-backend @@ -0,0 +1 @@ +#1985 : Backend work for TGV reconstruction From 0d07c184f1ef37615aa1023dabac39c10cd6fd25 Mon Sep 17 00:00:00 2001 From: Sam Tygier Date: Thu, 25 Apr 2024 10:05:33 +0100 Subject: [PATCH 11/11] Add better message if TGV with stochastic is attempted --- mantidimaging/core/reconstruct/cil_recon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mantidimaging/core/reconstruct/cil_recon.py b/mantidimaging/core/reconstruct/cil_recon.py index a6bca8eeeeb..f8622f2e637 100644 --- a/mantidimaging/core/reconstruct/cil_recon.py +++ b/mantidimaging/core/reconstruct/cil_recon.py @@ -108,7 +108,7 @@ def set_up_TGV_regularisation( f3 = MixedL21Norm() if recon_params.stochastic: - + raise ValueError("TGV reconstruction does not yet support stochastic mode") # now, A2d is a BlockOperator as acquisition_data is a BlockDataContainer fs = [] for i, _ in enumerate(acquisition_data.geometry):