New demo for MLEM/OSEM using real XRF data

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.

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 %%%%%%%%%%%")

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")

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)

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)