From fffd81d5b5990ed2d3a4548f538399ff6f4072e7 Mon Sep 17 00:00:00 2001 From: Param Date: Wed, 24 Jun 2026 12:03:29 +0530 Subject: [PATCH 01/11] 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 9faca5c7abe93bf06eade19220d0ccd016f19a1d Mon Sep 17 00:00:00 2001 From: Param Date: Wed, 24 Jun 2026 12:14:59 +0530 Subject: [PATCH 02/11] 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 fa350b9a98e9b6498d44c7385c734bf0237ddcc9 Mon Sep 17 00:00:00 2001 From: Param Date: Wed, 24 Jun 2026 12:31:25 +0530 Subject: [PATCH 03/11] 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 a3a66a741a55212e5c222c05a95f28ab61ccdeef Mon Sep 17 00:00:00 2001 From: Param Date: Wed, 24 Jun 2026 12:33:51 +0530 Subject: [PATCH 04/11] 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 912ca2059381638d6bb789c927bc265972dec757 Mon Sep 17 00:00:00 2001 From: Param Date: Wed, 24 Jun 2026 12:45:33 +0530 Subject: [PATCH 05/11] 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 435e4aed0a3ebba461ed9c358648cd4e26de786d Mon Sep 17 00:00:00 2001 From: Param Date: Wed, 24 Jun 2026 15:01:21 +0530 Subject: [PATCH 06/11] 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 5b1eb01f8a4b808b50122a0ab28da67b19a94088 Mon Sep 17 00:00:00 2001 From: Param Date: Tue, 1 Sep 2026 21:36:47 +0530 Subject: [PATCH 07/11] 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 d9bc31f5..51a35773 100644 --- a/Python/tigre/algorithms/iterative_recon_alg.py +++ b/Python/tigre/algorithms/iterative_recon_alg.py @@ -335,6 +335,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 @@ -343,7 +350,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 1e380f72b50cebb466baa9d4475a9b33dfc14f9d Mon Sep 17 00:00:00 2001 From: Param Date: Tue, 1 Sep 2026 21:46:47 +0530 Subject: [PATCH 08/11] 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 51a35773..788936a5 100644 --- a/Python/tigre/algorithms/iterative_recon_alg.py +++ b/Python/tigre/algorithms/iterative_recon_alg.py @@ -359,8 +359,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) From 69e30102141112b4a705587685ca2622f0c373ab Mon Sep 17 00:00:00 2001 From: Param Date: Mon, 28 Sep 2026 12:13:42 +0530 Subject: [PATCH 09/11] fix(build): import setuptools before Cython to fix Python 3.12 wheel build --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 7da395e5..02743908 100644 --- a/setup.py +++ b/setup.py @@ -8,9 +8,9 @@ import shutil import re +from setuptools import setup, find_packages, Extension from Cython.Distutils import build_ext import numpy -from setuptools import setup, find_packages, Extension # Macros for compilation From a8fcc597f41ceed919b848ca6b98582dea391bcc Mon Sep 17 00:00:00 2001 From: CircuitPilot Bot Date: Mon, 28 Sep 2026 12:28:14 +0530 Subject: [PATCH 10/11] fix(build): remove uv from cibuildwheel to fix Python 3.12 wheel build dependency resolution --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index daf32d44..3cd737cb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -47,7 +47,7 @@ where = ["Python"] "tigre.utilities.cuda_interface" = ["*.pyx", "*.pxd", "*.cpp"] [tool.cibuildwheel] -build-frontend = "build[uv]" +# build-frontend = "build[uv]" skip = "pp* *musl*" archs = ["native"] From 2d0b5cfc9417d647a62f688c66b65fe79377f496 Mon Sep 17 00:00:00 2001 From: CircuitPilot Bot Date: Mon, 28 Sep 2026 12:43:36 +0530 Subject: [PATCH 11/11] fix(build): pin Cython<3.0.0 and numpy<2.0.0 to fix Python 3.12 wheel build --- pyproject.toml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 3cd737cb..3447fe9b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,8 +5,8 @@ line-length = 100 requires = [ "setuptools", "wheel", - "Cython", - "numpy", + "Cython<3.0.0", + "numpy<2.0.0", ] build-backend = "setuptools.build_meta" @@ -15,7 +15,7 @@ name="tigre" version="3.1.3" dependencies = [ "matplotlib", - "numpy>=1.20.0", + "numpy>=1.20.0,<2.0.0", "scipy", "h5py", "tqdm"