From b712fde7a335ba29eedbf3429d618fa35e466962 Mon Sep 17 00:00:00 2001 From: Param Date: Wed, 24 Jun 2026 12:03:29 +0530 Subject: [PATCH 1/8] feat: Introduce Fast_OS_SART and Adaptive-Weighted TV to ART algorithms --- Python/tigre/algorithms/__init__.py | 6 + .../tigre/algorithms/art_family_algorithms.py | 111 +++++++++++++++++- 2 files changed, 116 insertions(+), 1 deletion(-) diff --git a/Python/tigre/algorithms/__init__.py b/Python/tigre/algorithms/__init__.py index c3f1e926..198be46d 100644 --- a/Python/tigre/algorithms/__init__.py +++ b/Python/tigre/algorithms/__init__.py @@ -7,6 +7,9 @@ from .art_family_algorithms import ossart from .art_family_algorithms import sart_tv from .art_family_algorithms import ossart_tv +from .art_family_algorithms import fast_os_sart +from .art_family_algorithms import aw_sart_tv +from .art_family_algorithms import aw_ossart_tv from .ista_algorithms import fista from .ista_algorithms import ista from .iterative_recon_alg import iterativereconalg @@ -38,6 +41,9 @@ "ossart", "sart_tv", "ossart_tv", + "fast_os_sart", + "aw_sart_tv", + "aw_ossart_tv", "iterativereconalg", "FDK", "asd_pocs", diff --git a/Python/tigre/algorithms/art_family_algorithms.py b/Python/tigre/algorithms/art_family_algorithms.py index 18f0fa33..38651536 100644 --- a/Python/tigre/algorithms/art_family_algorithms.py +++ b/Python/tigre/algorithms/art_family_algorithms.py @@ -1,5 +1,5 @@ import copy - +import numpy as np from tigre.algorithms.iterative_recon_alg import IterativeReconAlg from tigre.algorithms.iterative_recon_alg import decorator from tigre.utilities.im_3d_denoise import im3ddenoise @@ -150,3 +150,112 @@ def run_main_iter(self): self.error_measurement(res_prev, i) ossart_tv = decorator(OSSART_TV, name="ossart_tv") + +class Fast_OS_SART(IterativeReconAlg): + __doc__ = ( + "Fast_OS_SART solves Cone Beam CT image reconstruction using Nesterov accelerated\n" + "Oriented Subsets Simultaneous Algebraic Reconstruction Technique algorithm\n" + "Fast_OS_SART(PROJ,GEO,ALPHA,NITER,BLOCKSIZE=20) solves the reconstruction problem\n" + "using the projection data PROJ taken over ALPHA angles, corresponding\n" + "to the geometry described in GEO, using NITER iterations.\n" + ) + IterativeReconAlg.__doc__ + + def __init__(self, proj, geo, angles, niter, **kwargs): + self.blocksize = 20 if 'blocksize' not in kwargs else kwargs["blocksize"] + IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) + self.__t__ = 1.0 + + def run_main_iter(self): + Quameasopts = self.Quameasopts + t = self.__t__ + y_rec = copy.deepcopy(self.res) + + for i in range(self.niter): + res_prev = copy.deepcopy(self.res) if Quameasopts is not None else None + if self.verbose: + self._estimate_time_until_completion(i) + + x_rec_old = copy.deepcopy(self.res) + + self.res = copy.deepcopy(y_rec) + getattr(self, self.dataminimizing)() + + t_old = t + t = (1.0 + np.sqrt(1.0 + 4.0 * t ** 2)) / 2.0 + y_rec = self.res + (t_old - 1.0) / t * (self.res - x_rec_old) + + if Quameasopts is not None: + self.error_measurement(res_prev, i) + +fast_os_sart = decorator(Fast_OS_SART, name="fast_os_sart") + +class AwSART_TV(IterativeReconAlg): + __doc__ = ( + "AwSART_TV solves Cone Beam CT image reconstruction using Simultaneous \n" + "Algebraic Reconstruction Technique with Adaptive-Weighted TV regularization algorithm\n" + "AwSART_TV(PROJ,GEO,ALPHA,NITER,TVLAMBDA=50,TVITER=50,DELTA=-0.005) solves the reconstruction\n" + "problem using the projection data PROJ taken over ALPHA angles\n" + "corresponding to the geometry described in GEO, using NITER iterations. \n" + ) + IterativeReconAlg.__doc__ + + def __init__(self, proj, geo, angles, niter, **kwargs): + if "blocksize" in kwargs and kwargs['blocksize']>1: + print('Warning: blocksize is set to 1, please use an OS version of the algorithm for blocksize > 1') + kwargs.update(dict(blocksize=1)) + self.tvlambda = 50 if 'tvlambda' not in kwargs else kwargs['tvlambda'] + self.tviter = 50 if 'tviter' not in kwargs else kwargs['tviter'] + self.delta = np.float32(-0.005) if "delta" not in kwargs else kwargs["delta"] + + IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) + self.numiter_tv = self.tviter + + def run_main_iter(self): + Quameasopts = self.Quameasopts + for i in range(self.niter): + res_prev = None + if Quameasopts is not None: + res_prev = copy.deepcopy(self.res) + if self.verbose: + self._estimate_time_until_completion(i) + + getattr(self, self.dataminimizing)() + self.res = self.minimizeAwTV(self.res, self.tvlambda) + if Quameasopts is not None: + self.error_measurement(res_prev, i) + +aw_sart_tv = decorator(AwSART_TV, name="aw_sart_tv") + +class AwOSSART_TV(IterativeReconAlg): + __doc__ = ( + "AwOSSART_TV solves Cone Beam CT image reconstruction using Oriented Subsets\n" + "Simultaneous Algebraic Reconstruction Technique with Adaptive-Weighted TV regularization\n" + "AwOSSART_TV(PROJ,GEO,ALPHA,NITER,BLOCKSIZE=20,TVLAMBDA=50,TVITER=50,DELTA=-0.005) \n" + "solves the reconstruction problem using the projection data PROJ taken\n" + "over ALPHA angles, corresponding to the geometry described in GEO,\n" + "using NITER iterations.\n" + ) + IterativeReconAlg.__doc__ + + def __init__(self, proj, geo, angles, niter, **kwargs): + self.blocksize = 20 if 'blocksize' not in kwargs else kwargs['blocksize'] + self.tvlambda = 50 if 'tvlambda' not in kwargs else kwargs['tvlambda'] + self.tviter = 50 if 'tviter' not in kwargs else kwargs['tviter'] + self.delta = np.float32(-0.005) if "delta" not in kwargs else kwargs["delta"] + + IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) + self.numiter_tv = self.tviter + + def run_main_iter(self): + Quameasopts = self.Quameasopts + for i in range(self.niter): + res_prev = None + if Quameasopts is not None: + res_prev = copy.deepcopy(self.res) + if self.verbose: + self._estimate_time_until_completion(i) + + getattr(self, self.dataminimizing)() + self.res = self.minimizeAwTV(self.res, self.tvlambda) + if Quameasopts is not None: + self.error_measurement(res_prev, i) + +aw_ossart_tv = decorator(AwOSSART_TV, name="aw_ossart_tv") From 471dc2945b754241b913c965be66db729eb68e94 Mon Sep 17 00:00:00 2001 From: Param Date: Wed, 24 Jun 2026 12:14:59 +0530 Subject: [PATCH 2/8] test: add convergence and quality benchmarking script --- generate_benchmarks.py | 86 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) create mode 100644 generate_benchmarks.py diff --git a/generate_benchmarks.py b/generate_benchmarks.py new file mode 100644 index 00000000..c88256e7 --- /dev/null +++ b/generate_benchmarks.py @@ -0,0 +1,86 @@ +import sys +import time +import numpy as np +import matplotlib.pyplot as plt + +# Attempt to import TIGRE +try: + import tigre + import tigre.algorithms as algs + from tigre.demos.Test_data import data_loader + from tigre.utilities.Measure_Quality import Measure_Quality +except ImportError: + print("ERROR: TIGRE is not properly installed or compiled.") + print("Please run this script in an environment with TIGRE's C++/CUDA backend compiled.") + sys.exit(1) + +def run_benchmarks(): + print("--- Setting up TIGRE Geometry & Phantom ---") + # 1. Setup geometry and phantom + geo = tigre.geometry_default(high_resolution=False) + geo.nVoxel = np.array([64, 64, 64]) # Use small voxel size for fast benchmarking + + # Generate angles + angles = np.linspace(0, 2 * np.pi, 100) + + # Load phantom + head = data_loader.load_head_phantom(geo.nVoxel) + + # Generate projection data + print("Generating forward projections...") + proj = tigre.Ax(head, geo, angles) + + niter = 30 + blocksize = 20 + + print("\n--- Benchmark 1: Convergence Speed & Time (OS_SART vs Fast_OS_SART) ---") + + # Standard OS_SART + print("Running standard OS_SART...") + start_time = time.time() + res_os_sart, err_os_sart = algs.ossart(proj, geo, angles, niter=niter, blocksize=blocksize, computel2=True) + time_os_sart = time.time() - start_time + + # Fast OS_SART + print("Running Fast_OS_SART...") + start_time = time.time() + res_fast, err_fast = algs.fast_os_sart(proj, geo, angles, niter=niter, blocksize=blocksize, computel2=True) + time_fast = time.time() - start_time + + print(f"OS_SART Total Time: {time_os_sart:.2f}s ({time_os_sart/niter:.3f}s per iteration)") + print(f"Fast_OS_SART Total Time: {time_fast:.2f}s ({time_fast/niter:.3f}s per iteration)") + print(f"Final L2 Error -> OS_SART: {err_os_sart[0][-1]:.4f} | Fast_OS_SART: {err_fast[0][-1]:.4f}") + + # Plot convergence + plt.figure(figsize=(8, 5)) + plt.plot(err_os_sart[0], label="OS_SART", linewidth=2) + plt.plot(err_fast[0], label="Fast_OS_SART", linewidth=2) + plt.title("Convergence Speed: OS_SART vs Fast_OS_SART") + plt.xlabel("Iteration") + plt.ylabel("L2 Error") + plt.legend() + plt.grid(True) + plt.savefig("convergence_benchmark.png") + print("Saved convergence plot to convergence_benchmark.png") + + print("\n--- Benchmark 2: Image Quality (OSSART_TV vs AwOSSART_TV) ---") + + # Standard OSSART_TV + print("Running OSSART_TV...") + res_tv = algs.ossart_tv(proj, geo, angles, niter=15, blocksize=blocksize, tviter=20, tvlambda=50) + + # AwOSSART_TV + print("Running AwOSSART_TV...") + res_awtv = algs.aw_ossart_tv(proj, geo, angles, niter=15, blocksize=blocksize, tviter=20, tvlambda=50, delta=-0.005) + + # Calculate Metrics (RMSE and SSIM) + rmse_tv = Measure_Quality(res_tv, head, ['RMSE']) + rmse_awtv = Measure_Quality(res_awtv, head, ['RMSE']) + + print(f"RMSE -> OSSART_TV: {rmse_tv[0]:.6f} | AwOSSART_TV: {rmse_awtv[0]:.6f}") + + print("\n--- Benchmarks Complete! ---") + print("You can copy these results into your GitHub PR.") + +if __name__ == '__main__': + run_benchmarks() From 5bc87295d897d9e6c41c511dc560a96919873eea Mon Sep 17 00:00:00 2001 From: Param Date: Wed, 24 Jun 2026 12:31:25 +0530 Subject: [PATCH 3/8] fix: correct import path for data_loader --- generate_benchmarks.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/generate_benchmarks.py b/generate_benchmarks.py index c88256e7..9e0d8164 100644 --- a/generate_benchmarks.py +++ b/generate_benchmarks.py @@ -7,7 +7,7 @@ try: import tigre import tigre.algorithms as algs - from tigre.demos.Test_data import data_loader + from tigre.utilities.sample_loader import load_head_phantom from tigre.utilities.Measure_Quality import Measure_Quality except ImportError: print("ERROR: TIGRE is not properly installed or compiled.") @@ -24,7 +24,7 @@ def run_benchmarks(): angles = np.linspace(0, 2 * np.pi, 100) # Load phantom - head = data_loader.load_head_phantom(geo.nVoxel) + head = load_head_phantom(geo.nVoxel) # Generate projection data print("Generating forward projections...") From 3026f20262ececa50cf1ed49d1709f80dfb5210d Mon Sep 17 00:00:00 2001 From: Param Date: Wed, 24 Jun 2026 12:33:51 +0530 Subject: [PATCH 4/8] fix: explicit float32 cast to avoid float64 type mismatch in C++ backend --- .../tigre/algorithms/art_family_algorithms.py | 31 ++++++++++--------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/Python/tigre/algorithms/art_family_algorithms.py b/Python/tigre/algorithms/art_family_algorithms.py index 38651536..5c0246a9 100644 --- a/Python/tigre/algorithms/art_family_algorithms.py +++ b/Python/tigre/algorithms/art_family_algorithms.py @@ -174,21 +174,22 @@ def run_main_iter(self): res_prev = copy.deepcopy(self.res) if Quameasopts is not None else None if self.verbose: self._estimate_time_until_completion(i) - - x_rec_old = copy.deepcopy(self.res) - - self.res = copy.deepcopy(y_rec) - getattr(self, self.dataminimizing)() - - t_old = t - t = (1.0 + np.sqrt(1.0 + 4.0 * t ** 2)) / 2.0 - y_rec = self.res + (t_old - 1.0) / t * (self.res - x_rec_old) - - if Quameasopts is not None: - self.error_measurement(res_prev, i) - -fast_os_sart = decorator(Fast_OS_SART, name="fast_os_sart") - + + x_rec_old = copy.deepcopy(self.res) + + self.res = copy.deepcopy(y_rec) + getattr(self, self.dataminimizing)() + + t_old = t + t = (1.0 + np.sqrt(1.0 + 4.0 * t ** 2)) / 2.0 + y_rec = self.res + (t_old - 1.0) / t * (self.res - x_rec_old) + y_rec = np.float32(y_rec) + + if Quameasopts is not None: + self.error_measurement(res_prev, i) + +fast_os_sart = decorator(Fast_OS_SART, name="fast_os_sart") + class AwSART_TV(IterativeReconAlg): __doc__ = ( "AwSART_TV solves Cone Beam CT image reconstruction using Simultaneous \n" From 2cf1c92d6341fe44f5dfa9e9b3d8ac4747351907 Mon Sep 17 00:00:00 2001 From: Param Date: Wed, 24 Jun 2026 12:45:33 +0530 Subject: [PATCH 5/8] fix: cast RMSE to float to avoid index error --- generate_benchmarks.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/generate_benchmarks.py b/generate_benchmarks.py index 9e0d8164..4e36bf32 100644 --- a/generate_benchmarks.py +++ b/generate_benchmarks.py @@ -77,7 +77,7 @@ def run_benchmarks(): rmse_tv = Measure_Quality(res_tv, head, ['RMSE']) rmse_awtv = Measure_Quality(res_awtv, head, ['RMSE']) - print(f"RMSE -> OSSART_TV: {rmse_tv[0]:.6f} | AwOSSART_TV: {rmse_awtv[0]:.6f}") + print(f"RMSE -> OSSART_TV: {float(rmse_tv):.6f} | AwOSSART_TV: {float(rmse_awtv):.6f}") print("\n--- Benchmarks Complete! ---") print("You can copy these results into your GitHub PR.") From e0501e34300c770385915fd81a462e7d39d1ccc7 Mon Sep 17 00:00:00 2001 From: Param Date: Wed, 24 Jun 2026 15:01:21 +0530 Subject: [PATCH 6/8] fix: remove redundant AwSART_TV and AwOSSART_TV as they exist in pocs_algorithms --- Python/tigre/algorithms/__init__.py | 4 -- .../tigre/algorithms/art_family_algorithms.py | 71 ------------------- generate_benchmarks.py | 16 ----- 3 files changed, 91 deletions(-) diff --git a/Python/tigre/algorithms/__init__.py b/Python/tigre/algorithms/__init__.py index 198be46d..557ab968 100644 --- a/Python/tigre/algorithms/__init__.py +++ b/Python/tigre/algorithms/__init__.py @@ -8,8 +8,6 @@ from .art_family_algorithms import sart_tv from .art_family_algorithms import ossart_tv from .art_family_algorithms import fast_os_sart -from .art_family_algorithms import aw_sart_tv -from .art_family_algorithms import aw_ossart_tv from .ista_algorithms import fista from .ista_algorithms import ista from .iterative_recon_alg import iterativereconalg @@ -42,8 +40,6 @@ "sart_tv", "ossart_tv", "fast_os_sart", - "aw_sart_tv", - "aw_ossart_tv", "iterativereconalg", "FDK", "asd_pocs", diff --git a/Python/tigre/algorithms/art_family_algorithms.py b/Python/tigre/algorithms/art_family_algorithms.py index 5c0246a9..ce6f6465 100644 --- a/Python/tigre/algorithms/art_family_algorithms.py +++ b/Python/tigre/algorithms/art_family_algorithms.py @@ -189,74 +189,3 @@ def run_main_iter(self): self.error_measurement(res_prev, i) fast_os_sart = decorator(Fast_OS_SART, name="fast_os_sart") - -class AwSART_TV(IterativeReconAlg): - __doc__ = ( - "AwSART_TV solves Cone Beam CT image reconstruction using Simultaneous \n" - "Algebraic Reconstruction Technique with Adaptive-Weighted TV regularization algorithm\n" - "AwSART_TV(PROJ,GEO,ALPHA,NITER,TVLAMBDA=50,TVITER=50,DELTA=-0.005) solves the reconstruction\n" - "problem using the projection data PROJ taken over ALPHA angles\n" - "corresponding to the geometry described in GEO, using NITER iterations. \n" - ) + IterativeReconAlg.__doc__ - - def __init__(self, proj, geo, angles, niter, **kwargs): - if "blocksize" in kwargs and kwargs['blocksize']>1: - print('Warning: blocksize is set to 1, please use an OS version of the algorithm for blocksize > 1') - kwargs.update(dict(blocksize=1)) - self.tvlambda = 50 if 'tvlambda' not in kwargs else kwargs['tvlambda'] - self.tviter = 50 if 'tviter' not in kwargs else kwargs['tviter'] - self.delta = np.float32(-0.005) if "delta" not in kwargs else kwargs["delta"] - - IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) - self.numiter_tv = self.tviter - - def run_main_iter(self): - Quameasopts = self.Quameasopts - for i in range(self.niter): - res_prev = None - if Quameasopts is not None: - res_prev = copy.deepcopy(self.res) - if self.verbose: - self._estimate_time_until_completion(i) - - getattr(self, self.dataminimizing)() - self.res = self.minimizeAwTV(self.res, self.tvlambda) - if Quameasopts is not None: - self.error_measurement(res_prev, i) - -aw_sart_tv = decorator(AwSART_TV, name="aw_sart_tv") - -class AwOSSART_TV(IterativeReconAlg): - __doc__ = ( - "AwOSSART_TV solves Cone Beam CT image reconstruction using Oriented Subsets\n" - "Simultaneous Algebraic Reconstruction Technique with Adaptive-Weighted TV regularization\n" - "AwOSSART_TV(PROJ,GEO,ALPHA,NITER,BLOCKSIZE=20,TVLAMBDA=50,TVITER=50,DELTA=-0.005) \n" - "solves the reconstruction problem using the projection data PROJ taken\n" - "over ALPHA angles, corresponding to the geometry described in GEO,\n" - "using NITER iterations.\n" - ) + IterativeReconAlg.__doc__ - - def __init__(self, proj, geo, angles, niter, **kwargs): - self.blocksize = 20 if 'blocksize' not in kwargs else kwargs['blocksize'] - self.tvlambda = 50 if 'tvlambda' not in kwargs else kwargs['tvlambda'] - self.tviter = 50 if 'tviter' not in kwargs else kwargs['tviter'] - self.delta = np.float32(-0.005) if "delta" not in kwargs else kwargs["delta"] - - IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) - self.numiter_tv = self.tviter - - def run_main_iter(self): - Quameasopts = self.Quameasopts - for i in range(self.niter): - res_prev = None - if Quameasopts is not None: - res_prev = copy.deepcopy(self.res) - if self.verbose: - self._estimate_time_until_completion(i) - - getattr(self, self.dataminimizing)() - self.res = self.minimizeAwTV(self.res, self.tvlambda) - if Quameasopts is not None: - self.error_measurement(res_prev, i) - -aw_ossart_tv = decorator(AwOSSART_TV, name="aw_ossart_tv") diff --git a/generate_benchmarks.py b/generate_benchmarks.py index 4e36bf32..2c7a8a2c 100644 --- a/generate_benchmarks.py +++ b/generate_benchmarks.py @@ -63,22 +63,6 @@ def run_benchmarks(): plt.savefig("convergence_benchmark.png") print("Saved convergence plot to convergence_benchmark.png") - print("\n--- Benchmark 2: Image Quality (OSSART_TV vs AwOSSART_TV) ---") - - # Standard OSSART_TV - print("Running OSSART_TV...") - res_tv = algs.ossart_tv(proj, geo, angles, niter=15, blocksize=blocksize, tviter=20, tvlambda=50) - - # AwOSSART_TV - print("Running AwOSSART_TV...") - res_awtv = algs.aw_ossart_tv(proj, geo, angles, niter=15, blocksize=blocksize, tviter=20, tvlambda=50, delta=-0.005) - - # Calculate Metrics (RMSE and SSIM) - rmse_tv = Measure_Quality(res_tv, head, ['RMSE']) - rmse_awtv = Measure_Quality(res_awtv, head, ['RMSE']) - - print(f"RMSE -> OSSART_TV: {float(rmse_tv):.6f} | AwOSSART_TV: {float(rmse_awtv):.6f}") - print("\n--- Benchmarks Complete! ---") print("You can copy these results into your GitHub PR.") From 1ee2265a9c462ec882c32d11cec730ef2b40d78b Mon Sep 17 00:00:00 2001 From: Param Date: Tue, 1 Sep 2026 21:36:47 +0530 Subject: [PATCH 7/8] Refactor: Integrate Nesterov acceleration into base IterativeReconAlg (addresses maintainer review) --- Python/tigre/algorithms/__init__.py | 2 - .../tigre/algorithms/art_family_algorithms.py | 289 +++++++++++------- .../tigre/algorithms/iterative_recon_alg.py | 18 ++ generate_benchmarks.py | 14 +- 4 files changed, 207 insertions(+), 116 deletions(-) diff --git a/Python/tigre/algorithms/__init__.py b/Python/tigre/algorithms/__init__.py index 557ab968..c3f1e926 100644 --- a/Python/tigre/algorithms/__init__.py +++ b/Python/tigre/algorithms/__init__.py @@ -7,7 +7,6 @@ from .art_family_algorithms import ossart from .art_family_algorithms import sart_tv from .art_family_algorithms import ossart_tv -from .art_family_algorithms import fast_os_sart from .ista_algorithms import fista from .ista_algorithms import ista from .iterative_recon_alg import iterativereconalg @@ -39,7 +38,6 @@ "ossart", "sart_tv", "ossart_tv", - "fast_os_sart", "iterativereconalg", "FDK", "asd_pocs", diff --git a/Python/tigre/algorithms/art_family_algorithms.py b/Python/tigre/algorithms/art_family_algorithms.py index ce6f6465..3254b3b9 100644 --- a/Python/tigre/algorithms/art_family_algorithms.py +++ b/Python/tigre/algorithms/art_family_algorithms.py @@ -76,116 +76,191 @@ def __init__(self, proj, geo, angles, niter, **kwargs): if "blocksize" in kwargs and kwargs['blocksize']>1: print('Warning: blocksize is set to 1, please use an OS version of the algorithm for blocksize > 1') kwargs.update(dict(blocksize=1)) - self.tvlambda = 50 if 'tvlambda' not in kwargs else kwargs['tvlambda'] - self.tviter = 50 if 'tviter' not in kwargs else kwargs['tviter'] - # these two settings work well for nVoxel=[254,254,254] - - IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) - - # Override - def run_main_iter(self): - """ - Goes through the main iteration for the given configuration. - :return: None - """ - Quameasopts = self.Quameasopts - - for i in range(self.niter): - - res_prev = None - if Quameasopts is not None: - res_prev = copy.deepcopy(self.res) - if self.verbose: - self._estimate_time_until_completion(i) - - getattr(self, self.dataminimizing)() - # print("run_main_iter: gpuids = {}", self.gpuids) - self.res = im3ddenoise(self.res, self.tviter, self.tvlambda, self.gpuids) - if Quameasopts is not None: - self.error_measurement(res_prev, i) - - -sart_tv = decorator(SART_TV, name="sart_tv") - - -class OSSART_TV(IterativeReconAlg): - __doc__ = ( - "OSSART_TV solves Cone Beam CT image reconstruction using Oriented Subsets\n" - "Simultaneous Algebraic Reconstruction Technique with TV regularization algorithm\n" - "OSSART_TV(PROJ,GEO,ALPHA,NITER,BLOCKSIZE=20,TVLAMBDA=50,TVITER=50) \n" - "solves the reconstruction problem using the projection data PROJ taken\n" - "over ALPHA angles, corresponding to the geometry described in GEO,\n" - "using NITER iterations.\n" - ) + IterativeReconAlg.__doc__ - - def __init__(self, proj, geo, angles, niter, **kwargs): - - self.blocksize = 20 if 'blocksize' not in kwargs else kwargs['blocksize'] - self.tvlambda = 50 if 'tvlambda' not in kwargs else kwargs['tvlambda'] - self.tviter = 50 if 'tviter' not in kwargs else kwargs['tviter'] - # these two settings work well for nVoxel=[254,254,254] - - IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) - - # Override - def run_main_iter(self): - """ - Goes through the main iteration for the given configuration. - :return: None - """ - Quameasopts = self.Quameasopts - - for i in range(self.niter): - - res_prev = None - if Quameasopts is not None: - res_prev = copy.deepcopy(self.res) - if self.verbose: - self._estimate_time_until_completion(i) - - getattr(self, self.dataminimizing)() - # print("run_main_iter: gpuids = {}", self.gpuids) - self.res = im3ddenoise(self.res, self.tviter, self.tvlambda, self.gpuids) - if Quameasopts is not None: - self.error_measurement(res_prev, i) - -ossart_tv = decorator(OSSART_TV, name="ossart_tv") - -class Fast_OS_SART(IterativeReconAlg): - __doc__ = ( - "Fast_OS_SART solves Cone Beam CT image reconstruction using Nesterov accelerated\n" - "Oriented Subsets Simultaneous Algebraic Reconstruction Technique algorithm\n" - "Fast_OS_SART(PROJ,GEO,ALPHA,NITER,BLOCKSIZE=20) solves the reconstruction problem\n" - "using the projection data PROJ taken over ALPHA angles, corresponding\n" - "to the geometry described in GEO, using NITER iterations.\n" - ) + IterativeReconAlg.__doc__ - - def __init__(self, proj, geo, angles, niter, **kwargs): - self.blocksize = 20 if 'blocksize' not in kwargs else kwargs["blocksize"] - IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) - self.__t__ = 1.0 - - def run_main_iter(self): - Quameasopts = self.Quameasopts - t = self.__t__ - y_rec = copy.deepcopy(self.res) - - for i in range(self.niter): - res_prev = copy.deepcopy(self.res) if Quameasopts is not None else None - if self.verbose: - self._estimate_time_until_completion(i) - - x_rec_old = copy.deepcopy(self.res) - - self.res = copy.deepcopy(y_rec) +import copy +import numpy as np +from tigre.algorithms.iterative_recon_alg import IterativeReconAlg +from tigre.algorithms.iterative_recon_alg import decorator +from tigre.utilities.im_3d_denoise import im3ddenoise + + + +class SART(IterativeReconAlg): + __doc__ = ( + "SART solves Cone Beam CT image reconstruction using \n" + "Simultaneous Algebraic Reconstruction Technique algorithm\n" + "SART(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem\n" + "using the projection data PROJ taken over ALPHA angles, corresponding\n" + "to the geometry described in GEO, using NITER iterations. \n" + ) + IterativeReconAlg.__doc__ + + def __init__(self, proj, geo, angles, niter, **kwargs): + if "blocksize" in kwargs and kwargs['blocksize']>1: + print('Warning: blocksize is set to 1, please use an OS version of the algorithm for blocksize > 1') + kwargs.update(dict(blocksize=1)) + IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) + + +sart = decorator(SART, name="sart") + + +class SIRT(IterativeReconAlg): + __doc__ = ( + "SIRT solves Cone Beam CT image reconstruction using \n" + "Simultaneous Iterative Reconstructive Technique algorithm\n" + "SIRT(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem\n" + "using the projection data PROJ taken over ALPHA angles, corresponding\n" + "to the geometry described in GEO, using NITER iterations.\n" + ) + IterativeReconAlg.__doc__ + + def __init__(self, proj, geo, angles, niter, **kwargs): + if "blocksize" in kwargs and kwargs['blocksize']>1: + print('Warning: blocksize is set to {}, please do not specify blocksize for this algorithm'.format(angles.shape[0])) + kwargs.update(dict(blocksize=angles.shape[0])) + IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) + + +sirt = decorator(SIRT, name="sirt") + + +class OS_SART(IterativeReconAlg): + __doc__ = ( + "OS_SART solves Cone Beam CT image reconstruction using Oriented Subsets\n" + "Simultaneous Algebraic Reconstruction Technique algorithm\n" + "OS_SART(PROJ,GEO,ALPHA,NITER,BLOCKSIZE=20) solves the reconstruction problem\n" + "using the projection data PROJ taken over ALPHA angles, corresponding\n" + "to the geometry described in GEO, using NITER iterations.\n" + ) + IterativeReconAlg.__doc__ + + def __init__(self, proj, geo, angles, niter, **kwargs): + + self.blocksize = 20 if 'blocksize' not in kwargs else kwargs["blocksize"] + IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) + + +ossart = decorator(OS_SART, name="ossart") + + +class SART_TV(IterativeReconAlg): + __doc__ = ( + "SART_TV solves Cone Beam CT image reconstruction using Simultaneous \n" + "Algebraic Reconstruction Technique with TV regularization algorithm\n" + "SART_TV(PROJ,GEO,ALPHA,NITER,TVLAMBDA=50,TVITER=50) solves the reconstruction\n" + "problem using the projection data PROJ taken over ALPHA angles\n" + "corresponding to the geometry described in GEO, using NITER iterations. \n" + ) + IterativeReconAlg.__doc__ + + def __init__(self, proj, geo, angles, niter, **kwargs): + + if "blocksize" in kwargs and kwargs['blocksize']>1: + print('Warning: blocksize is set to 1, please use an OS version of the algorithm for blocksize > 1') + kwargs.update(dict(blocksize=1)) + self.tvlambda = 50 if 'tvlambda' not in kwargs else kwargs['tvlambda'] + self.tviter = 50 if 'tviter' not in kwargs else kwargs['tviter'] + # these two settings work well for nVoxel=[254,254,254] + + IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) + + # Override + def run_main_iter(self): + """ + Goes through the main iteration for the given configuration. + :return: None + """ + Quameasopts = self.Quameasopts + + nesterov = False + if isinstance(self.lmbda, str) and self.lmbda.lower() == "nesterov": + nesterov = True + self.lmbda = 1.0 + t = 1.0 + y_rec = copy.deepcopy(self.res) + + for i in range(self.niter): + + res_prev = None + if Quameasopts is not None: + res_prev = copy.deepcopy(self.res) + if self.verbose: + self._estimate_time_until_completion(i) + + if nesterov: + x_rec_old = copy.deepcopy(self.res) + self.res = copy.deepcopy(y_rec) + getattr(self, self.dataminimizing)() + # print("run_main_iter: gpuids = {}", self.gpuids) + self.res = im3ddenoise(self.res, self.tviter, self.tvlambda, self.gpuids) - t_old = t - t = (1.0 + np.sqrt(1.0 + 4.0 * t ** 2)) / 2.0 - y_rec = self.res + (t_old - 1.0) / t * (self.res - x_rec_old) - y_rec = np.float32(y_rec) + if nesterov: + t_old = t + t = (1.0 + np.sqrt(1.0 + 4.0 * t ** 2)) / 2.0 + y_rec = self.res + (t_old - 1.0) / t * (self.res - x_rec_old) + y_rec = np.float32(y_rec) + + if Quameasopts is not None: + self.error_measurement(res_prev, i) + + +sart_tv = decorator(SART_TV, name="sart_tv") + + +class OSSART_TV(IterativeReconAlg): + __doc__ = ( + "OSSART_TV solves Cone Beam CT image reconstruction using Oriented Subsets\n" + "Simultaneous Algebraic Reconstruction Technique with TV regularization algorithm\n" + "OSSART_TV(PROJ,GEO,ALPHA,NITER,BLOCKSIZE=20,TVLAMBDA=50,TVITER=50) \n" + "solves the reconstruction problem using the projection data PROJ taken\n" + "over ALPHA angles, corresponding to the geometry described in GEO,\n" + "using NITER iterations.\n" + ) + IterativeReconAlg.__doc__ + + def __init__(self, proj, geo, angles, niter, **kwargs): + + self.blocksize = 20 if 'blocksize' not in kwargs else kwargs['blocksize'] + self.tvlambda = 50 if 'tvlambda' not in kwargs else kwargs['tvlambda'] + self.tviter = 50 if 'tviter' not in kwargs else kwargs['tviter'] + # these two settings work well for nVoxel=[254,254,254] + + IterativeReconAlg.__init__(self, proj, geo, angles, niter, **kwargs) + + # Override + def run_main_iter(self): + """ + Goes through the main iteration for the given configuration. + :return: None + """ + Quameasopts = self.Quameasopts + + nesterov = False + if isinstance(self.lmbda, str) and self.lmbda.lower() == "nesterov": + nesterov = True + self.lmbda = 1.0 + t = 1.0 + y_rec = copy.deepcopy(self.res) + + for i in range(self.niter): + + res_prev = None + if Quameasopts is not None: + res_prev = copy.deepcopy(self.res) + if self.verbose: + self._estimate_time_until_completion(i) + if nesterov: + x_rec_old = copy.deepcopy(self.res) + self.res = copy.deepcopy(y_rec) + + getattr(self, self.dataminimizing)() + # print("run_main_iter: gpuids = {}", self.gpuids) + self.res = im3ddenoise(self.res, self.tviter, self.tvlambda, self.gpuids) + + if nesterov: + t_old = t + t = (1.0 + np.sqrt(1.0 + 4.0 * t ** 2)) / 2.0 + y_rec = self.res + (t_old - 1.0) / t * (self.res - x_rec_old) + y_rec = np.float32(y_rec) + if Quameasopts is not None: self.error_measurement(res_prev, i) -fast_os_sart = decorator(Fast_OS_SART, name="fast_os_sart") +ossart_tv = decorator(OSSART_TV, name="ossart_tv") diff --git a/Python/tigre/algorithms/iterative_recon_alg.py b/Python/tigre/algorithms/iterative_recon_alg.py index c6e6333d..860c22a9 100644 --- a/Python/tigre/algorithms/iterative_recon_alg.py +++ b/Python/tigre/algorithms/iterative_recon_alg.py @@ -312,6 +312,13 @@ def run_main_iter(self): """ Quameasopts = self.Quameasopts + nesterov = False + if isinstance(self.lmbda, str) and self.lmbda.lower() == "nesterov": + nesterov = True + self.lmbda = 1.0 + t = 1.0 + y_rec = copy.deepcopy(self.res) + for i in range(self.niter): res_prev = None @@ -320,7 +327,18 @@ def run_main_iter(self): if self.verbose: self._estimate_time_until_completion(i) + if nesterov: + x_rec_old = copy.deepcopy(self.res) + self.res = copy.deepcopy(y_rec) + getattr(self, self.dataminimizing)() + + if nesterov: + t_old = t + t = (1.0 + np.sqrt(1.0 + 4.0 * t ** 2)) / 2.0 + y_rec = self.res + (t_old - 1.0) / t * (self.res - x_rec_old) + y_rec = np.float32(y_rec) + self.error_measurement(res_prev, i) def art_data_minimizing(self): diff --git a/generate_benchmarks.py b/generate_benchmarks.py index 2c7a8a2c..aa89c356 100644 --- a/generate_benchmarks.py +++ b/generate_benchmarks.py @@ -33,7 +33,7 @@ def run_benchmarks(): niter = 30 blocksize = 20 - print("\n--- Benchmark 1: Convergence Speed & Time (OS_SART vs Fast_OS_SART) ---") + print("\n--- Benchmark 1: Convergence Speed & Time (OS_SART vs OS_SART with Nesterov) ---") # Standard OS_SART print("Running standard OS_SART...") @@ -42,20 +42,20 @@ def run_benchmarks(): time_os_sart = time.time() - start_time # Fast OS_SART - print("Running Fast_OS_SART...") + print("Running OS_SART with Nesterov acceleration...") start_time = time.time() - res_fast, err_fast = algs.fast_os_sart(proj, geo, angles, niter=niter, blocksize=blocksize, computel2=True) + res_fast, err_fast = algs.ossart(proj, geo, angles, niter=niter, blocksize=blocksize, lmbda='nesterov', computel2=True) time_fast = time.time() - start_time print(f"OS_SART Total Time: {time_os_sart:.2f}s ({time_os_sart/niter:.3f}s per iteration)") - print(f"Fast_OS_SART Total Time: {time_fast:.2f}s ({time_fast/niter:.3f}s per iteration)") - print(f"Final L2 Error -> OS_SART: {err_os_sart[0][-1]:.4f} | Fast_OS_SART: {err_fast[0][-1]:.4f}") + print(f"OS_SART (Nesterov) Total Time: {time_fast:.2f}s ({time_fast/niter:.3f}s per iteration)") + print(f"Final L2 Error -> OS_SART: {err_os_sart[0][-1]:.4f} | OS_SART (Nesterov): {err_fast[0][-1]:.4f}") # Plot convergence plt.figure(figsize=(8, 5)) plt.plot(err_os_sart[0], label="OS_SART", linewidth=2) - plt.plot(err_fast[0], label="Fast_OS_SART", linewidth=2) - plt.title("Convergence Speed: OS_SART vs Fast_OS_SART") + plt.plot(err_fast[0], label="OS_SART (Nesterov)", linewidth=2) + plt.title("Convergence Speed: OS_SART vs OS_SART (Nesterov)") plt.xlabel("Iteration") plt.ylabel("L2 Error") plt.legend() From 758beccefad7a515089219217d9e923ed8dcd021 Mon Sep 17 00:00:00 2001 From: Param Date: Tue, 1 Sep 2026 21:46:47 +0530 Subject: [PATCH 8/8] Refactor: Optimize Nesterov memory allocation by casting scalar to float32 --- Python/tigre/algorithms/art_family_algorithms.py | 8 ++++---- Python/tigre/algorithms/iterative_recon_alg.py | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/Python/tigre/algorithms/art_family_algorithms.py b/Python/tigre/algorithms/art_family_algorithms.py index 3254b3b9..359a732f 100644 --- a/Python/tigre/algorithms/art_family_algorithms.py +++ b/Python/tigre/algorithms/art_family_algorithms.py @@ -194,8 +194,8 @@ def run_main_iter(self): if nesterov: t_old = t t = (1.0 + np.sqrt(1.0 + 4.0 * t ** 2)) / 2.0 - y_rec = self.res + (t_old - 1.0) / t * (self.res - x_rec_old) - y_rec = np.float32(y_rec) + gamma = np.float32((t_old - 1.0) / t) + y_rec = self.res + gamma * (self.res - x_rec_old) if Quameasopts is not None: self.error_measurement(res_prev, i) @@ -257,8 +257,8 @@ def run_main_iter(self): if nesterov: t_old = t t = (1.0 + np.sqrt(1.0 + 4.0 * t ** 2)) / 2.0 - y_rec = self.res + (t_old - 1.0) / t * (self.res - x_rec_old) - y_rec = np.float32(y_rec) + gamma = np.float32((t_old - 1.0) / t) + y_rec = self.res + gamma * (self.res - x_rec_old) if Quameasopts is not None: self.error_measurement(res_prev, i) diff --git a/Python/tigre/algorithms/iterative_recon_alg.py b/Python/tigre/algorithms/iterative_recon_alg.py index 860c22a9..5de242b4 100644 --- a/Python/tigre/algorithms/iterative_recon_alg.py +++ b/Python/tigre/algorithms/iterative_recon_alg.py @@ -336,8 +336,8 @@ def run_main_iter(self): if nesterov: t_old = t t = (1.0 + np.sqrt(1.0 + 4.0 * t ** 2)) / 2.0 - y_rec = self.res + (t_old - 1.0) / t * (self.res - x_rec_old) - y_rec = np.float32(y_rec) + gamma = np.float32((t_old - 1.0) / t) + y_rec = self.res + gamma * (self.res - x_rec_old) self.error_measurement(res_prev, i)