pyproject.toml

@@ -20,7 +20,7 @@ dev_template = "{tag}"
 
 [project]
 name = "tomobar"
-version = "2024.11"
+version = "2025.03"
 description = "TOmographic MOdel-BAsed Reconstruction (ToMoBAR) software"
 readme = "Readme.md"
 license = {text = "GPLv3"}
@@ -34,10 +34,10 @@ classifiers = [
     "Programming Language :: Python :: 3.10",
     "Environment :: GPU :: NVIDIA CUDA"
 ]
-requires-python = ">=3.9"
+requires-python = ">=3.10"
 dependencies = [
     "numpy",
-    "astra-toolbox",
+    "astra-toolbox>=2.3.0",
     "scipy",
     "pillow",
     "scikit-image"

tests/test_RecToolsDIRCuPy.py

@@ -32,8 +32,8 @@ def test_Fourier3D_inv(data_cupy, angles, ensure_clean_memory):
             data_cupy, data_axes_labels_order=["angles", "detY", "detX"]
         )
     recon_data = Fourier_rec_cupy.get()
-    assert_allclose(np.min(recon_data), -0.04679, rtol=1e-05)
-    assert_allclose(np.max(recon_data), 0.105128, rtol=1e-05)
+    assert_allclose(np.min(recon_data), -0.037292, rtol=1e-05)
+    assert_allclose(np.max(recon_data), 0.103775, rtol=1e-05)
     assert recon_data.dtype == np.float32
     assert recon_data.shape == (128, 160, 160)
 

tomobar/astra_wrappers/astra_base.py

@@ -8,6 +8,7 @@
 from typing import Union
 from tomobar.supp.funcs import _vec_geom_init2D, _vec_geom_init3D
 from astra.experimental import direct_BP3D, direct_FP3D
+from astra.pythonutils import GPULink
 
 try:
     import cupy as xp
@@ -503,7 +504,7 @@ def runAstraBackproj3DCuPy(
             xp.ndarray: A CuPy array containing back-projected volume.
         """
         # set ASTRA configuration for 3D reconstructor using CuPy arrays
-        proj_link = astra.data3d.GPULink(
+        proj_link = GPULink(
             proj_data.data.ptr, *proj_data.shape[::-1], 4 * proj_data.shape[2]
         )
         if self.ordsub_number != 1 and os_index is not None:
@@ -519,7 +520,7 @@ def runAstraBackproj3DCuPy(
         # create a CuPy array with ASTRA link to it
         recon_volume = xp.empty(astra.geom_size(self.vol_geom), dtype=xp.float32)
 
-        rec_link = astra.data3d.GPULink(
+        rec_link = GPULink(
             recon_volume.data.ptr, *recon_volume.shape[::-1], 4 * recon_volume.shape[2]
         )
 
@@ -543,7 +544,7 @@ def runAstraProj3DCuPy(
             xp.ndarray: projected volume array as a cupy array
         """
         # Enable GPUlink to the volume
-        volume_link = astra.data3d.GPULink(
+        volume_link = GPULink(
             volume_data.data.ptr, *volume_data.shape[::-1], 4 * volume_data.shape[2]
         )
         volume_id = astra.data3d.link("-vol", self.vol_geom, volume_link)
@@ -553,7 +554,7 @@ def runAstraProj3DCuPy(
             proj_volume = xp.empty(
                 astra.geom_size(self.proj_geom_OS[os_index]), dtype=xp.float32
             )
-            gpu_link_sino = astra.data3d.GPULink(
+            gpu_link_sino = GPULink(
                 proj_volume.data.ptr, *proj_volume.shape[::-1], 4 * proj_volume.shape[2]
             )
             projector_id = astra.create_projector(
@@ -563,7 +564,7 @@ def runAstraProj3DCuPy(
             # traditional full data parallel beam projection geometry
             # Enabling GPUlink to the created empty CuPy array
             proj_volume = xp.empty(astra.geom_size(self.proj_geom), dtype=xp.float32)
-            gpu_link_sino = astra.data3d.GPULink(
+            gpu_link_sino = GPULink(
                 proj_volume.data.ptr, *proj_volume.shape[::-1], 4 * proj_volume.shape[2]
             )
             projector_id = astra.create_projector(

tomobar/methodsDIR_CuPy.py

@@ -131,7 +131,7 @@ def FBP(self, data: xp.ndarray, **kwargs) -> xp.ndarray:
         data = _filtersinc3D_cupy(data, cutoff=cutoff_freq)
         data = xp.ascontiguousarray(xp.swapaxes(data, 0, 1))
         cache = xp.fft.config.get_plan_cache()
-        cache.clear()  # flush FFT cache here before backprojection        
+        cache.clear()  # flush FFT cache here before backprojection
         xp._default_memory_pool.free_all_blocks()  # free everything related to the filtering before starting Astra
         reconstruction = self.Atools._backprojCuPy(data)  # 3d backprojecting
         xp._default_memory_pool.free_all_blocks()
@@ -157,12 +157,12 @@ def FOURIER_INV(self, data: xp.ndarray, **kwargs) -> xp.ndarray:
             xp.ndarray: The NUFFT reconstructed volume as a CuPy array.
         """
 
-        center_size = 768
-        block_dim = [64, 4]
-        block_dim_center = [32, 4]
+        cutoff_freq = 1.0  # default value
+        filter_type = "shepp"  # default filter
 
-        cutoff_freq = 0.35  # default value
-        filter_type = "ramp"  # default filter
+        center_size = 512
+        block_dim = [32, 32]
+        block_dim_center = [32, 8]
 
         for key, value in kwargs.items():
             if key == "data_axes_labels_order" and value is not None:
@@ -173,9 +173,9 @@ def FOURIER_INV(self, data: xp.ndarray, **kwargs) -> xp.ndarray:
                 block_dim = value
             elif key == "block_dim_center" and value is not None:
                 block_dim_center = value
-            elif key == "cutoff_freq" and value is not None:
+            if key == "cutoff_freq" and value is not None:
                 cutoff_freq = value
-            elif key == "filter_type" and value is not None:
+            if key == "filter_type" and value is not None:
                 if value not in [
                     "none",
                     "ramp",
@@ -187,7 +187,7 @@ def FOURIER_INV(self, data: xp.ndarray, **kwargs) -> xp.ndarray:
                     "parzen",
                 ]:
                     print(
-                        "Unknown filter name, please use: none, ramp, shepp, cosine, cosine2, hamming, hann or parzen. Set to shepp"
+                        "Unknown filter name, please use: none, ramp, shepp, cosine, cosine2, hamming, hann or parzen. Set to shepp filter"
                     )
                 cutoff_freq = value
 
@@ -225,7 +225,7 @@ def FOURIER_INV(self, data: xp.ndarray, **kwargs) -> xp.ndarray:
             odd_vert = True
             del data_p
 
-        rotation_axis = self.Atools.centre_of_rotation - 0.5
+        rotation_axis = self.Atools.centre_of_rotation + 0.5
         if odd_horiz:
             rotation_axis -= 1
         theta = xp.array(-self.Atools.angles_vec, dtype=xp.float32)
@@ -305,28 +305,23 @@ def FOURIER_INV(self, data: xp.ndarray, **kwargs) -> xp.ndarray:
         # STEP2: interpolation (gathering) in the frequency domain
         if center_size > 0:
 
-            tic = timeit.default_timer()
-
-            gather_kernel_partial(
-                (int(xp.ceil(n / block_dim[0])), int(xp.ceil(nproj / block_dim[1])), nz // 2),
-                (block_dim[0], block_dim[1], 1),
-                (
-                    datac,
-                    fde,
-                    theta,
-                    np.int32(m),
-                    np.float32(mu),
-                    np.int32(center_size),
-                    np.int32(n),
-                    np.int32(nproj),
-                    np.int32(nz // 2),
-                ),
-            )
-
-            # Run_time = (timeit.default_timer() - tic)
-            # print("gather_kernel_partial in {} seconds".format(Run_time))
-
-            tic = timeit.default_timer()
+            if center_size != (n * 2 + m * 2):
+
+                gather_kernel_partial(
+                    (int(xp.ceil(n / block_dim[0])), int(xp.ceil(nproj / block_dim[1])), nz // 2),
+                    (block_dim[0], block_dim[1], 1),
+                    (
+                        datac,
+                        fde,
+                        theta,
+                        np.int32(m),
+                        np.float32(mu),
+                        np.int32(center_size),
+                        np.int32(n),
+                        np.int32(nproj),
+                        np.int32(nz // 2),
+                    ),
+                )
 
             gather_kernel_center_prune(
                 (1, int(xp.ceil(center_size / 4)), center_size),
@@ -355,11 +350,6 @@ def FOURIER_INV(self, data: xp.ndarray, **kwargs) -> xp.ndarray:
             # plt.title("Angle max")
             # plt.show()
 
-            # Run_time = (timeit.default_timer() - tic)
-            # print("gather_kernel_center_prune in {} seconds".format(Run_time))
-
-            tic = timeit.default_timer()
-
             gather_kernel_center(
                 (
                     int(xp.ceil(center_size / block_dim_center[0])),
@@ -381,12 +371,7 @@ def FOURIER_INV(self, data: xp.ndarray, **kwargs) -> xp.ndarray:
                 ),
             )
 
-            # Run_time = (timeit.default_timer() - tic)
-            # print("gather_kernel_center in {} seconds".format(Run_time))
-
         else:
-            tic = timeit.default_timer()
-
             gather_kernel(
                 (int(xp.ceil(n / block_dim[0])), int(xp.ceil(nproj / block_dim[1])), nz // 2),
                 (block_dim[0], block_dim[1], 1),
@@ -402,9 +387,6 @@ def FOURIER_INV(self, data: xp.ndarray, **kwargs) -> xp.ndarray:
                 ),
             )
 
-            # Run_time = (timeit.default_timer() - tic)
-            # print("gather_kernel in {} seconds".format(Run_time))
-
         wrap_kernel(
             (
                 int(np.ceil((2 * n + 2 * m) / 32)),