draft for fista cupy, gpu tests added

Demos/Python/Demo_CuPy_3D.py

@@ -26,7 +26,7 @@
 print("Building 3D phantom using TomoPhantom software")
 tic = timeit.default_timer()
 model = 13  # select a model number from the library
-N_size = 128  # Define phantom dimensions using a scalar value (cubic phantom)
+N_size = 350  # Define phantom dimensions using a scalar value (cubic phantom)
 path = os.path.dirname(tomophantom.__file__)
 path_library3D = os.path.join(path, "Phantom3DLibrary.dat")
 
@@ -35,7 +35,7 @@
 # Projection geometry related parameters:
 Horiz_det = int(np.sqrt(2) * N_size)  # detector column count (horizontal)
 Vert_det = N_size  # detector row count (vertical) (no reason for it to be > N)
-angles_num = int(0.25 * np.pi * N_size)  # angles number
+angles_num = int(0.3 * np.pi * N_size)  # angles number
 angles = np.linspace(0.0, 179.9, angles_num, dtype="float32")  # in degrees
 angles_rad = angles * (np.pi / 180.0)
 
@@ -61,7 +61,7 @@
 )
 
 tic = timeit.default_timer()
-FBPrec_cupy = RecToolsCP.FBP3D(projData3D_analyt_cupy)
+FBPrec_cupy = RecToolsCP.FBP3D(projData3D_analyt_cupy, recon_mask_radius = 0.9)
 toc = timeit.default_timer()
 Run_time = toc - tic
 print(
@@ -93,9 +93,76 @@
     "Min {} and Max {} of the volume".format(np.min(FBPrec_cupy), np.max(FBPrec_cupy))
 )
 del FBPrec_cupy
+# # %%
+# print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+# print("%%%%%%%%Reconstructing using Landweber algorithm %%%%%%%%%%%")
+# print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+# RecToolsCP = RecToolsIRCuPy(
+#     DetectorsDimH=Horiz_det,  # Horizontal detector dimension
+#     DetectorsDimV=Vert_det,  # Vertical detector dimension (3D case)
+#     CenterRotOffset=0.0,  # 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="gpu",
+# )
+
+# # prepare dictionaries with parameters:
+# _data_ = {"projection_norm_data": projData3D_analyt_cupy}  # data dictionary
+# LWrec_cupy = RecToolsCP.Landweber(_data_)
+
+# lwrec = cp.asnumpy(LWrec_cupy)
+
+# sliceSel = int(0.5 * N_size)
+# plt.figure()
+# plt.subplot(131)
+# plt.imshow(lwrec[sliceSel, :, :])
+# plt.title("3D Landweber Reconstruction, axial view")
+
+# plt.subplot(132)
+# plt.imshow(lwrec[:, sliceSel, :])
+# plt.title("3D Landweber Reconstruction, coronal view")
+
+# plt.subplot(133)
+# plt.imshow(lwrec[:, :, sliceSel])
+# plt.title("3D Landweber Reconstruction, sagittal view")
+# plt.show()
+# # %%
+# print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+# print("%%%%%%%%%% Reconstructing using SIRT algorithm %%%%%%%%%%%%%")
+# print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+# RecToolsCP = RecToolsIRCuPy(
+#     DetectorsDimH=Horiz_det,  # Horizontal detector dimension
+#     DetectorsDimV=Vert_det,  # Vertical detector dimension (3D case)
+#     CenterRotOffset=0.0,  # 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="gpu",
+# )
+
+# # prepare dictionaries with parameters:
+# _data_ = {"projection_norm_data": projData3D_analyt_cupy}  # data dictionary
+# _algorithm_ = {"iterations": 300, "nonnegativity": True}
+# SIRTrec_cupy = RecToolsCP.SIRT(_data_, _algorithm_)
+
+# sirt_rec = cp.asnumpy(SIRTrec_cupy)
+
+# sliceSel = int(0.5 * N_size)
+# plt.figure()
+# plt.subplot(131)
+# plt.imshow(sirt_rec[sliceSel, :, :])
+# plt.title("3D SIRT Reconstruction, axial view")
+
+# plt.subplot(132)
+# plt.imshow(sirt_rec[:, sliceSel, :])
+# plt.title("3D SIRT Reconstruction, coronal view")
+
+# plt.subplot(133)
+# plt.imshow(sirt_rec[:, :, sliceSel])
+# plt.title("3D SIRT Reconstruction, sagittal view")
+# plt.show()
 # %%
 print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print("%%%%%%%%Reconstructing using Landweber algorithm %%%%%%%%%%%")
+print("%%%%%%%%%% Reconstructing using CGLS algorithm %%%%%%%%%%%%%")
 print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 RecToolsCP = RecToolsIRCuPy(
     DetectorsDimH=Horiz_det,  # Horizontal detector dimension
@@ -108,28 +175,26 @@
 
 # prepare dictionaries with parameters:
 _data_ = {"projection_norm_data": projData3D_analyt_cupy}  # data dictionary
-LWrec_cupy = RecToolsCP.Landweber(_data_)
+_algorithm_ = {"iterations": 20, "nonnegativity": True}
+CGLSrec_cupy = RecToolsCP.CGLS(_data_, _algorithm_)
 
-lwrec = cp.asnumpy(LWrec_cupy)
+cgls_rec = cp.asnumpy(CGLSrec_cupy)
 
 sliceSel = int(0.5 * N_size)
 plt.figure()
 plt.subplot(131)
-plt.imshow(lwrec[sliceSel, :, :])
-plt.title("3D Landweber Reconstruction, axial view")
+plt.imshow(cgls_rec[sliceSel, :, :])
+plt.title("3D CGLS Reconstruction, axial view")
 
 plt.subplot(132)
-plt.imshow(lwrec[:, sliceSel, :])
-plt.title("3D Landweber Reconstruction, coronal view")
+plt.imshow(cgls_rec[:, sliceSel, :])
+plt.title("3D CGLS Reconstruction, coronal view")
 
 plt.subplot(133)
-plt.imshow(lwrec[:, :, sliceSel])
-plt.title("3D Landweber Reconstruction, sagittal view")
+plt.imshow(cgls_rec[:, :, sliceSel])
+plt.title("3D CGLS Reconstruction, sagittal view")
 plt.show()
 # %%
-print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print("%%%%%%%%%% Reconstructing using SIRT algorithm %%%%%%%%%%%%%")
-print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 RecToolsCP = RecToolsIRCuPy(
     DetectorsDimH=Horiz_det,  # Horizontal detector dimension
     DetectorsDimV=Vert_det,  # Vertical detector dimension (3D case)
@@ -139,59 +204,36 @@
     device_projector="gpu",
 )
 
-# prepare dictionaries with parameters:
-_data_ = {"projection_norm_data": projData3D_analyt_cupy}  # data dictionary
-_algorithm_ = {"iterations": 300, "nonnegativity": True}
-SIRTrec_cupy = RecToolsCP.SIRT(_data_, _algorithm_)
-
-sirt_rec = cp.asnumpy(SIRTrec_cupy)
-
-sliceSel = int(0.5 * N_size)
-plt.figure()
-plt.subplot(131)
-plt.imshow(sirt_rec[sliceSel, :, :])
-plt.title("3D SIRT Reconstruction, axial view")
-
-plt.subplot(132)
-plt.imshow(sirt_rec[:, sliceSel, :])
-plt.title("3D SIRT Reconstruction, coronal view")
-
-plt.subplot(133)
-plt.imshow(sirt_rec[:, :, sliceSel])
-plt.title("3D SIRT Reconstruction, sagittal view")
-plt.show()
-# %%
-print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print("%%%%%%%%%% Reconstructing using CGLS algorithm %%%%%%%%%%%%%")
-print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-RecToolsCP = RecToolsIRCuPy(
-    DetectorsDimH=Horiz_det,  # Horizontal detector dimension
-    DetectorsDimV=Vert_det,  # Vertical detector dimension (3D case)
-    CenterRotOffset=0.0,  # 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="gpu",
-)
 
 # prepare dictionaries with parameters:
 _data_ = {"projection_norm_data": projData3D_analyt_cupy}  # data dictionary
-_algorithm_ = {"iterations": 20, "nonnegativity": True}
-CGLSrec_cupy = RecToolsCP.CGLS(_data_, _algorithm_)
-
-cgls_rec = cp.asnumpy(CGLSrec_cupy)
+_algorithm_ = {"iterations": 200, "lipschitz_const": 50000.346}
+_regularisation_ = {
+    "method": "PD_TV",
+    "regul_param": 0.001,
+    "iterations": 100,
+    "device_regulariser": "gpu",
+}
+start_time = timeit.default_timer()
+RecFISTA = RecToolsCP.FISTA(_data_, _algorithm_, _regularisation_)
+txtstr = "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+
+fista_rec_np = cp.asnumpy(RecFISTA)
 
 sliceSel = int(0.5 * N_size)
 plt.figure()
 plt.subplot(131)
-plt.imshow(cgls_rec[sliceSel, :, :])
-plt.title("3D CGLS Reconstruction, axial view")
+plt.imshow(fista_rec_np[sliceSel, :, :])
+plt.title("3D FISTA Reconstruction, axial view")
 
 plt.subplot(132)
-plt.imshow(cgls_rec[:, sliceSel, :])
-plt.title("3D CGLS Reconstruction, coronal view")
+plt.imshow(fista_rec_np[:, sliceSel, :])
+plt.title("3D FISTA Reconstruction, coronal view")
 
 plt.subplot(133)
-plt.imshow(cgls_rec[:, :, sliceSel])
-plt.title("3D CGLS Reconstruction, sagittal view")
+plt.imshow(fista_rec_np[:, :, sliceSel])
+plt.title("3D FISTA Reconstruction, sagittal view")
 plt.show()
+del RecFISTA
 # %%

tests_gpu/conftest.py

@@ -0,0 +1,44 @@
+# Defines common fixtures and makes them available to all tests
+
+import os
+import numpy as np
+import pytest
+import cupy as cp
+
+CUR_DIR = os.path.abspath(os.path.dirname(__file__))
+
+@pytest.fixture(scope="session")
+def test_data_path():
+    return os.path.join(CUR_DIR, "test_data")
+
+# only load from disk once per session, and we use np.copy for the elements,
+# to ensure data in this loaded file stays as originally loaded
+@pytest.fixture(scope="session")
+def data_file(test_data_path):
+    in_file = os.path.join(test_data_path, "normalised_data.npz")
+    return np.load(in_file)
+
+@pytest.fixture
+def ensure_clean_memory():
+    cp.get_default_memory_pool().free_all_blocks()
+    cp.get_default_pinned_memory_pool().free_all_blocks()
+    yield None
+    cp.get_default_memory_pool().free_all_blocks()
+    cp.get_default_pinned_memory_pool().free_all_blocks()
+
+
+@pytest.fixture
+def data(data_file):
+    return np.copy(data_file["data_norm"])
+
+@pytest.fixture
+def data_cupy(data):
+    return cp.asarray(data)
+
+@pytest.fixture
+def angles(data_file):
+    return np.copy(data_file["angles"])
+
+@pytest.fixture
+def angles_cupy(angles):
+    return cp.asarray(angles)

tests_gpu/test_direct_recon_cupy.py

@@ -0,0 +1,49 @@
+from unittest import mock
+import cupy as cp
+from cupy.cuda import nvtx
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+
+from tomobar.methodsDIR_CuPy import RecToolsDIRCuPy
+
+eps = 1e-06
+@cp.testing.gpu
+def test_rec_FBPcupy(data_cupy, angles, ensure_clean_memory):
+    detX=cp.shape(data_cupy)[2]
+    detY=cp.shape(data_cupy)[1]
+    N_size = detX
+    RecToolsCP = RecToolsDIRCuPy(
+        DetectorsDimH=detX,  # Horizontal detector dimension
+        DetectorsDimV=detY,  # Vertical detector dimension (3D case)
+        CenterRotOffset=0.0,  # Center of Rotation scalar or a vector
+        AnglesVec=angles,  # A vector of projection angles in radians
+        ObjSize=N_size,  # Reconstructed object dimensions (scalar)
+        device_projector="gpu",
+    )
+    FBPrec_cupy = RecToolsCP.FBP3D(data_cupy)
+    recon_data = FBPrec_cupy.get()
+    assert_allclose(np.min(recon_data), -0.014693323, rtol=eps)
+    assert_allclose(np.max(recon_data), 0.0340156, rtol=eps)
+    assert recon_data.dtype == np.float32
+    assert recon_data.shape == (128, 160, 160)
+
+@cp.testing.gpu
+def test_rec_FBP_mask_cupy(data_cupy, angles, ensure_clean_memory):
+    detX=cp.shape(data_cupy)[2]
+    detY=cp.shape(data_cupy)[1]
+    N_size = detX
+    RecToolsCP = RecToolsDIRCuPy(
+        DetectorsDimH=detX,  # Horizontal detector dimension
+        DetectorsDimV=detY,  # Vertical detector dimension (3D case)
+        CenterRotOffset=0.0,  # Center of Rotation scalar or a vector
+        AnglesVec=angles,  # A vector of projection angles in radians
+        ObjSize=N_size,  # Reconstructed object dimensions (scalar)
+        device_projector="gpu",
+    )
+    FBPrec_cupy = RecToolsCP.FBP3D(data_cupy, recon_mask_radius = 0.7)
+    recon_data = FBPrec_cupy.get()
+    assert_allclose(np.min(recon_data), -0.0129751, rtol=eps)
+    assert_allclose(np.max(recon_data), 0.0340156, rtol=eps)
+    assert recon_data.dtype == np.float32
+    assert recon_data.shape == (128, 160, 160)    

tomobar/methodsDIR.py

@@ -253,4 +253,4 @@ def FBP(self, sinogram):
         if self.geom == "3D":
             filtered_sino = filtersinc3D(sinogram)  # filtering projection data
             FBP_rec = self.Atools.backproj(filtered_sino)  # 3d backproject
-        return FBP_rec
+        return FBP_rec

tomobar/methodsDIR_CuPy.py

@@ -16,6 +16,7 @@
 
 from tomobar.cuda_kernels import load_cuda_module
 from tomobar.methodsDIR import RecToolsDIR
+from tomobar.supp.suppTools import _check_kwargs, _apply_circular_mask
 
 
 def _filtersinc3D_cupy(projection3D: cp.ndarray) -> cp.ndarray:
@@ -60,6 +61,7 @@ def _filtersinc3D_cupy(projection3D: cp.ndarray) -> cp.ndarray:
 
     return cupyx.scipy.fft.irfft(proj_f, projection3D.shape[2], axis=-1, norm="forward", overwrite_x=True)
 
+
 class RecToolsDIRCuPy(RecToolsDIR):
     def __init__(
         self,
@@ -79,12 +81,16 @@ def __init__(
             device_projector,
         )
 
-    def FBP3D(self, data: cp.ndarray) -> cp.ndarray:
+    def FBP3D(self, data: cp.ndarray, **kwargs) -> cp.ndarray:
         """Filtered backprojection on a CuPy array using a custom built SINC filter
 
         Args:
             data : cp.ndarray
                 Projection data as a CuPy array.
+            kwargs: dict
+                Optional parameters:
+                1. recon_mask_radius - zero the values outside the circular mask
+                of a certain radius. To see the effect of the cropping, set the value in the range [0.7-1.0].
 
         Returns:
             cp.ndarray
@@ -95,5 +101,5 @@ def FBP3D(self, data: cp.ndarray) -> cp.ndarray:
         )  # filter the data on the GPU and keep the result there
         data = cp.ascontiguousarray(cp.swapaxes(data, 0, 1))
         reconstruction = self.Atools.backprojCuPy(data)  # 3d backprojecting
-        cp._default_memory_pool.free_all_blocks()
-        return reconstruction
+        cp._default_memory_pool.free_all_blocks()        
+        return _check_kwargs(reconstruction, **kwargs)