From 3537d17b50514c78b1cca101a87352124a6270b7 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 15 Nov 2023 17:07:47 +0000 Subject: [PATCH 1/9] 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 b4e1bf52e9a..c487a0c6f39 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 4c0bce1441d9e1f5d9f095c6bc668984a0ed938b Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 29 Nov 2023 12:24:08 +0000 Subject: [PATCH 2/9] Return a np.array for TGV case --- mantidimaging/core/reconstruct/cil_recon.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/mantidimaging/core/reconstruct/cil_recon.py b/mantidimaging/core/reconstruct/cil_recon.py index c487a0c6f39..69f34a1ad7a 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 @@ -384,8 +389,11 @@ def full(images: ImageStack, f': Objective {algo.get_last_objective():.2f}', 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('Reconstructed 3D volume with shape: {0}'.format(volume.shape)) t1 = time.perf_counter() LOG.info(f"full reconstruction time: {t1-t0}s for shape {images.data.shape}") From 8638ff2d0ec18e11338dd452f8932bac7ab9cd10 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 29 Nov 2023 12:40:22 +0000 Subject: [PATCH 3/9] 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 69f34a1ad7a..3457fc4bf6f 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 34fbf47ff65b322615ac0448314f2af57150ab30 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 30 Nov 2023 09:21:43 +0000 Subject: [PATCH 4/9] 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 3457fc4bf6f..9d16f3c6dbf 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 07e1783a3d1b7ad40afa4d331e6063608fc5c78c Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 30 Nov 2023 09:23:45 +0000 Subject: [PATCH 5/9] 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 9d16f3c6dbf..3694d577655 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 af1843ca8f9e11de59f2117b1d637a950a58a679 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 1 Dec 2023 13:42:31 +0000 Subject: [PATCH 6/9] 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 3694d577655..a1e18ee01f8 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 3bd0d21ef93cb3da7563c8bdc7986ea0260107d5 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 5 Dec 2023 14:04:12 +0000 Subject: [PATCH 7/9] 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 a1e18ee01f8..6e88cf5bde2 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 ea6556bc6bc..8a892e83abb 100644 --- a/mantidimaging/core/utility/data_containers.py +++ b/mantidimaging/core/utility/data_containers.py @@ -104,12 +104,14 @@ class ReconstructionParameters: tilt: Optional[Degrees] = 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: Optional[List[float]] = None stochastic: bool = False projections_per_subset: int = 50 regularisation_percent: int = 30 + regulariser: str = "" def to_dict(self) -> dict: return { @@ -120,9 +122,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 55a08561b14..8daef645c52 100644 --- a/mantidimaging/gui/windows/recon/view.py +++ b/mantidimaging/gui/windows/recon/view.py @@ -410,6 +410,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() @@ -426,6 +430,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) -> Optional[List[float]]: if not self.lbhc_enabled.isChecked(): @@ -446,11 +454,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 0a39716c43b09c578d80f4a0fe2df0fc5ac8b855 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 5 Dec 2023 14:29:33 +0000 Subject: [PATCH 8/9] 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 6e88cf5bde2..6f41365bb61 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 a4c5b5e9753f21edb40a31e5837c1fba3fd407f1 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 5 Dec 2023 14:35:10 +0000 Subject: [PATCH 9/9] 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 6f41365bb61..408c200817b 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