diff --git a/CHANGELOG.md b/CHANGELOG.md index 8523ed462..60e436b37 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,9 @@ # Change Log +## v.2026.3.1.0 + +Stripe Weighted Least Squares (SWLS) [method](http://ieeexplore.ieee.org/document/7967846/) added to the list of data fidelities. Demo with dendrites (RealData.py) demonstrates the benefit of using such method. + ## v.2026.3.0.0 We apologise for the abrupt change, but starting from this version, iterative reconstruction methods will no longer be accessible through the `RecToolsIR` class. Only `RecToolsIRCuPy` will be supported going forward. @@ -22,7 +26,6 @@ Changes: * Upgrade of dependencies to CuPy 14.* and numpy 2.4. - ## v.2026.2.0.0 The ADMM implementation has been replaced with a new linearised variant that avoids the use of SciPy linear algebra solvers and instead relies directly on gradient-descent–type iterations. With appropriate parameter selection, this approach converges quickly and exhibits improved numerical stability. diff --git a/Demos/RealData.py b/Demos/RealData.py index 8fe1befe6..40335f93b 100644 --- a/Demos/RealData.py +++ b/Demos/RealData.py @@ -12,19 +12,21 @@ import timeit import numpy as np import matplotlib.pyplot as plt -import scipy.io import cupy as cp from tomobar.supp.suppTools import normaliser from tomobar.methodsIR_CuPy import RecToolsIRCuPy from tomobar.methodsDIR_CuPy import RecToolsDIRCuPy # load raw X-ray data of dendrites -datadict = scipy.io.loadmat("../data/DendrRawData.mat") +data_npz = np.load( + "../tests/test_data/DendrRawData.npz" +) # you need to be in ToMoBAR/Demos folder + # extract data (print(datadict.keys())) -dataRaw = datadict["data_raw3D"] -angles = datadict["angles"] -flats = datadict["flats_ar"] -darks = datadict["darks_ar"] +dataRaw = data_npz["arr_0"] +angles = data_npz["arr_1"] +flats = data_npz["arr_2"] +darks = data_npz["arr_3"] # normalise the data data_norm = normaliser( @@ -39,6 +41,7 @@ N_size = 1000 angles_rad = np.linspace(0, np.pi, 360) +del data_npz # %% print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print("%%%%Reconstructing with Log-Polar Fourier method %%%%%%%%%%%") diff --git a/Demos/RealData3_XRF.py b/Demos/RealData3_XRF.py new file mode 100644 index 000000000..ba8c8f5eb --- /dev/null +++ b/Demos/RealData3_XRF.py @@ -0,0 +1,99 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# +# In this demo XRF tomographic data from I14 beamline at Diamond Light Source (UK synchrotron) is reconstructed + + +import numpy as np +import cupy as cp +import matplotlib.pyplot as plt +from tomobar.methodsIR_CuPy import RecToolsIRCuPy +from tomobar.methodsDIR_CuPy import RecToolsDIRCuPy + +data_npz = np.load( + "../tests/test_data/Ga-Ka_aligned.npz" +) # you need to be in ToMoBAR/Demos folder + +data = data_npz["arr_0"] +angles = data_npz["arr_1"] + +angles_tot, detectorVert, detectorHoriz = np.shape(data) +data_labels3D = ["angles", "detY", "detX"] +data_cp = cp.asarray(np.float32(data), order="C") +# %% +print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print("%%%%%%%%%%%%Reconstructing with FBP method %%%%%%%%%%%%%%%%%") +print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +RectoolsDIR_cp = RecToolsDIRCuPy( + DetectorsDimH=detectorHoriz, # Horizontal detector dimension + DetectorsDimH_pad=45, # Padding size of horizontal detector + DetectorsDimV=detectorVert, # Vertical detector dimension (3D case) + CenterRotOffset=-0.5, # Centre of Rotation scalar + AnglesVec=np.deg2rad(angles), # A vector of projection angles in radians + ObjSize=detectorHoriz, # Reconstructed object dimensions (scalar) + device_projector=0, +) + +FBPrec_cupy = RectoolsDIR_cp.FBP(data_cp, data_axes_labels_order=data_labels3D) + +FBPrec_np = cp.asnumpy(FBPrec_cupy) + +fig = plt.figure() +plt.subplot(121) +plt.imshow(FBPrec_np[20, :, :], cmap="gray") +plt.title("FBP axial") +plt.subplot(122) +plt.imshow(FBPrec_np[:, 45, :], cmap="gray") +plt.title("FBP saggital") + +# vmin=-0.0001, vmax=0.002, +# %% +print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print("%%%%%%%%%%%%Reconstructing with OSEM-TV method %%%%%%%%%%%%%") +print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +RectoolsIR_cp = RecToolsIRCuPy( + DetectorsDimH=detectorHoriz, # Horizontal detector dimension + DetectorsDimH_pad=45, # Padding size of horizontal detector + DetectorsDimV=detectorVert, # Vertical detector dimension (3D case) + CenterRotOffset=-0.5, # Centre of Rotation scalar + AnglesVec=np.deg2rad(angles), # A vector of projection angles in radians + ObjSize=detectorHoriz, # Reconstructed object dimensions (scalar) + device_projector=0, + OS_number=12, # The number of ordered subsets +) + +####################### Creating the data dictionary: ####################### +_data_ = { + "projection_data": data_cp, + "data_axes_labels_order": data_labels3D, +} + +####################### Creating the algorithm dictionary: ####################### +_algorithm_ = { + "iterations": 20, + "recon_mask_radius": 2.0, + "nonnegativity": True, +} # The number of iterations + +##### creating regularisation dictionary: ##### +_regularisation_ = { + "method": "PD_TV", # Selected regularisation method + "regul_param": 1.5, # Regularisation parameter + "iterations": 35, # The number of regularisation iterations + "half_precision": True, # enabling half-precision calculation +} + + +# RUN MLEM/OSEM METHOD: +Rec_OSEM = RectoolsIR_cp.OSEM(_data_, _algorithm_, _regularisation_) + + +Rec_OSEM_np = cp.asnumpy(Rec_OSEM) + +fig = plt.figure() +plt.subplot(121) +plt.imshow(Rec_OSEM_np[20, :, :], cmap="gray") +plt.title("OSEM-TV axial") +plt.subplot(122) +plt.imshow(Rec_OSEM_np[:, 45, :], cmap="gray") +plt.title("OSEM-TV saggital") diff --git a/data/DendrRawData.mat b/data/DendrRawData.mat deleted file mode 100644 index 2de7f0cd4..000000000 Binary files a/data/DendrRawData.mat and /dev/null differ diff --git a/data/RealDataDend.mat b/data/RealDataDend.mat deleted file mode 100644 index c4876c356..000000000 Binary files a/data/RealDataDend.mat and /dev/null differ diff --git a/data/data_icecream.h5 b/data/data_icecream.h5 deleted file mode 100755 index ac13d124e..000000000 Binary files a/data/data_icecream.h5 and /dev/null differ diff --git a/data/exportDemoRD2Data.m b/data/exportDemoRD2Data.m deleted file mode 100644 index 2c0791dc4..000000000 --- a/data/exportDemoRD2Data.m +++ /dev/null @@ -1,31 +0,0 @@ -function [Sino3D, Weights3D, angles_rad, size_det] = exportDemoRD2Data() - -load('DendrRawData.mat') % load raw data of 3D dendritic set - -angles_rad = angles*(pi/180); % conversion to radians -size_det = size(data_raw3D,1); % detectors dim -angSize = size(data_raw3D, 2); % angles dim -slices_tot = size(data_raw3D, 3); % no of slices - - -Sino3D = zeros(size_det, angSize, slices_tot, 'single'); % log-corrected sino -% normalizing the data -for jj = 1:slices_tot - sino = data_raw3D(:,:,jj); - for ii = 1:angSize - Sino3D(:,ii,jj) = log((flats_ar(:,jj)-darks_ar(:,jj))./(single(sino(:,ii)) - darks_ar(:,jj))); - end -end - -Sino3D = Sino3D.*1000; -Weights3D = single(data_raw3D); % weights for PWLS model -clear data_raw3D - -% hdf5write('DendrData.h5', '/Weights3D', Weights3D) -% hdf5write('DendrData.h5', '/Sino3D', Sino3D, 'WriteMode', 'append') -% hdf5write('DendrData.h5', '/angles_rad', angles_rad, 'WriteMode', 'append') -% hdf5write('DendrData.h5', '/size_det', size_det, 'WriteMode', 'append') -% hdf5write('DendrData.h5', '/angSize', angSize, 'WriteMode', 'append') -% hdf5write('DendrData.h5', '/slices_tot', slices_tot, 'WriteMode', 'append') - -return \ No newline at end of file diff --git a/tests/conftest.py b/tests/conftest.py index b400c0881..4b2d6471e 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -89,11 +89,27 @@ def data_raw_file(test_data_path): return np.load(in_file) +@pytest.fixture(scope="session") +def data_XRFfile(test_data_path): + in_file = os.path.join(test_data_path, "Ga-Ka_aligned.npz") + return np.load(in_file) + + @pytest.fixture def raw_data(data_raw_file): return np.float32(np.copy(data_raw_file["data"])) +@pytest.fixture +def raw_data_Xrf(data_XRFfile): + return np.float32(np.copy(data_XRFfile["arr_0"])) + + +@pytest.fixture +def angles_data_Xrf(data_XRFfile): + return np.float32(np.copy(data_XRFfile["arr_1"])) + + @pytest.fixture def raw_data_cupy(raw_data): return cp.asarray(raw_data) diff --git a/tests/test_RecToolsIRCuPy.py b/tests/test_RecToolsIRCuPy.py index de90c3cb0..c6f1f9a6b 100644 --- a/tests/test_RecToolsIRCuPy.py +++ b/tests/test_RecToolsIRCuPy.py @@ -1125,3 +1125,127 @@ def test_ADMM_SWLS_cp_3D( assert_allclose(cp.max(Iter_rec).get(), expected_max, rtol=0, atol=eps) assert Iter_rec.dtype == np.float32 assert Iter_rec.shape == (128, 160, 160) + + +def test_MLEM_cp_3D(raw_data_Xrf, angles_data_Xrf, ensure_clean_memory): + _, detY, detX = np.shape(raw_data_Xrf) + data_cp = cp.asarray(np.float32(raw_data_Xrf), order="C") + angles_rad = np.deg2rad(angles_data_Xrf) + data_labels3D = ["angles", "detY", "detX"] + + N_size = detX + RecTools = RecToolsIRCuPy( + DetectorsDimH=detX, # Horizontal detector dimension + DetectorsDimH_pad=0, # Padding size of horizontal detector + DetectorsDimV=detY, # Vertical detector dimension (3D case) + CenterRotOffset=-0.5, # Center of Rotation scalar or a vector + AnglesVec=angles_rad, # A vector of projection angles in radians + ObjSize=N_size, # Reconstructed object dimensions (scalar) + device_projector=0, + OS_number=1, # The number of ordered subsets + ) + + _data_ = { + "projection_data": data_cp, + "data_axes_labels_order": data_labels3D, + } + + ####################### Creating the algorithm dictionary: ####################### + _algorithm_ = { + "iterations": 5, + "recon_mask_radius": 2.0, + "nonnegativity": True, + } # The number of iterations + + Iter_rec = RecTools.OSEM(_data_, _algorithm_) + + Iter_rec = Iter_rec.get() + assert_allclose(np.min(Iter_rec), 0.06285925, rtol=1e-04) + assert_allclose(np.max(Iter_rec), 836700.8, rtol=1e-04) + assert Iter_rec.dtype == np.float32 + assert Iter_rec.shape == (41, 90, 90) + + +def test_OSEM_cp_3D(raw_data_Xrf, angles_data_Xrf, ensure_clean_memory): + _, detY, detX = np.shape(raw_data_Xrf) + data_cp = cp.asarray(np.float32(raw_data_Xrf), order="C") + angles_rad = np.deg2rad(angles_data_Xrf) + data_labels3D = ["angles", "detY", "detX"] + + N_size = detX + RecTools = RecToolsIRCuPy( + DetectorsDimH=detX, # Horizontal detector dimension + DetectorsDimH_pad=0, # Padding size of horizontal detector + DetectorsDimV=detY, # Vertical detector dimension (3D case) + CenterRotOffset=-0.5, # Center of Rotation scalar or a vector + AnglesVec=angles_rad, # A vector of projection angles in radians + ObjSize=N_size, # Reconstructed object dimensions (scalar) + device_projector=0, + OS_number=12, # The number of ordered subsets + ) + + _data_ = { + "projection_data": data_cp, + "data_axes_labels_order": data_labels3D, + } + + ####################### Creating the algorithm dictionary: ####################### + _algorithm_ = { + "iterations": 5, + "recon_mask_radius": 2.0, + "nonnegativity": True, + } # The number of iterations + + Iter_rec = RecTools.OSEM(_data_, _algorithm_) + + Iter_rec = Iter_rec.get() + assert_allclose(np.min(Iter_rec), 0.0, rtol=1e-04) + assert_allclose(np.max(Iter_rec), 8849.165, rtol=1e-04) + assert Iter_rec.dtype == np.float32 + assert Iter_rec.shape == (41, 90, 90) + + +def test_OSEM_PDTV_cp_3D(raw_data_Xrf, angles_data_Xrf, ensure_clean_memory): + _, detY, detX = np.shape(raw_data_Xrf) + data_cp = cp.asarray(np.float32(raw_data_Xrf), order="C") + angles_rad = np.deg2rad(angles_data_Xrf) + data_labels3D = ["angles", "detY", "detX"] + + N_size = detX + RecTools = RecToolsIRCuPy( + DetectorsDimH=detX, # Horizontal detector dimension + DetectorsDimH_pad=0, # Padding size of horizontal detector + DetectorsDimV=detY, # Vertical detector dimension (3D case) + CenterRotOffset=-0.5, # Center of Rotation scalar or a vector + AnglesVec=angles_rad, # A vector of projection angles in radians + ObjSize=N_size, # Reconstructed object dimensions (scalar) + device_projector=0, + OS_number=12, # The number of ordered subsets + ) + + _data_ = { + "projection_data": data_cp, + "data_axes_labels_order": data_labels3D, + } + + ####################### Creating the algorithm dictionary: ####################### + _algorithm_ = { + "iterations": 3, + "recon_mask_radius": 2.0, + "nonnegativity": True, + } # The number of iterations + + _regularisation_ = { + "method": "PD_TV", # Selected regularisation method + "regul_param": 1.5, # Regularisation parameter + "iterations": 5, # The number of regularisation iterations + "half_precision": True, # enabling half-precision calculation + } + + Iter_rec = RecTools.OSEM(_data_, _algorithm_, _regularisation_) + + Iter_rec = Iter_rec.get() + assert_allclose(np.min(Iter_rec), 3.2815723e-07, rtol=1e-04) + assert_allclose(np.max(Iter_rec), 8659.158, rtol=1e-04) + assert Iter_rec.dtype == np.float32 + assert Iter_rec.shape == (41, 90, 90) diff --git a/tests/test_data/DendrRawData.npz b/tests/test_data/DendrRawData.npz new file mode 100644 index 000000000..dda79854d Binary files /dev/null and b/tests/test_data/DendrRawData.npz differ diff --git a/tests/test_data/Ga-Ka_aligned.npz b/tests/test_data/Ga-Ka_aligned.npz new file mode 100644 index 000000000..293f07e3d Binary files /dev/null and b/tests/test_data/Ga-Ka_aligned.npz differ