Merge pull request #70 from dkazanc/multigpu

FFT filtering 1d to 2D + CuPy Filtering and GPUlink to reconstruction in Astra

Demos/Python/DemoFBP_CuPy_3D.py

@@ -0,0 +1,141 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+GPLv3 license (ASTRA toolbox)
+
+Script to generate 3D analytical phantoms and their projection data with added
+noise and then reconstruct using 3D FBP and 3D FBP with filtering on a GPU using CuPy
+
+Dependencies:
+    * astra-toolkit, install conda install -c astra-toolbox astra-toolbox
+    * CCPi-RGL toolkit (for regularisation), install with
+    conda install ccpi-regulariser -c ccpi -c conda-forge
+    or https://github.com/vais-ral/CCPi-Regularisation-Toolkit
+    * TomoPhantom, https://github.com/dkazanc/TomoPhantom
+
+@author: Daniil Kazantsev
+"""
+import timeit
+import os
+import matplotlib.pyplot as plt
+import numpy as np
+import tomophantom
+from tomophantom import TomoP3D
+from tomobar.methodsDIR import RecToolsDIR
+
+print ("Building 3D phantom using TomoPhantom software")
+tic=timeit.default_timer()
+model = 13 # select a model number from the library
+N_size = 600 # Define phantom dimensions using a scalar value (cubic phantom)
+path = os.path.dirname(tomophantom.__file__)
+path_library3D = os.path.join(path, "Phantom3DLibrary.dat")
+
+# 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 = np.linspace(0.0,179.9,angles_num,dtype='float32') # in degrees
+angles_rad = angles*(np.pi/180.0)
+
+print ("Generate 3D analytical projection data with TomoPhantom")
+projData3D_analyt= TomoP3D.ModelSino(model, N_size, Horiz_det, Vert_det, angles, path_library3D)
+#%%
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print ("%%%Reconstructing with 3D FBP method (slice-by-slice ) %%%%%")
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+RectoolsDIR = RecToolsDIR(DetectorsDimH = Horiz_det,     # Horizontal detector dimension
+                    DetectorsDimV = Vert_det,            # Vertical detector dimension (3D case)
+                    CenterRotOffset = None,              # 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')
+
+tic=timeit.default_timer()
+FBPrec = RectoolsDIR.FBP(projData3D_analyt) #perform 3D FBP (slice-by-slice)
+toc=timeit.default_timer()
+Run_time = toc - tic
+print("FBP 3D reconstruction slice-by-slice in {} seconds".format(Run_time))
+
+sliceSel = int(0.5*N_size)
+max_val = 1
+plt.figure()
+plt.subplot(131)
+plt.imshow(FBPrec[sliceSel,:,:],vmin=0, vmax=max_val)
+plt.title('3D FBP Reconstruction, axial view')
+
+plt.subplot(132)
+plt.imshow(FBPrec[:,sliceSel,:],vmin=0, vmax=max_val)
+plt.title('3D FBP Reconstruction, coronal view')
+
+plt.subplot(133)
+plt.imshow(FBPrec[:,:,sliceSel],vmin=0, vmax=max_val)
+plt.title('3D FBP Reconstruction, sagittal view')
+plt.show()
+
+print("Min {} and Max {} of the volume".format(np.min(FBPrec), np.max(FBPrec)))
+del(FBPrec)
+#%%
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print ("%%%%%%%%%%%Reconstructing with 3D FBP method %%%%%%%%%%%%%%%")
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+RectoolsDIR = RecToolsDIR(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')
+
+tic=timeit.default_timer()
+FBPrec = RectoolsDIR.FBP(projData3D_analyt) #perform 3D FBP with filtering on a CPU
+toc=timeit.default_timer()
+Run_time = toc - tic
+print("FBP 3D reconstruction with FFT filtering on a CPU done in {} seconds".format(Run_time))
+
+sliceSel = int(0.5*N_size)
+max_val = 1
+plt.figure() 
+plt.subplot(131)
+plt.imshow(FBPrec[sliceSel,:,:],vmin=0, vmax=max_val)
+plt.title('3D FBP Reconstruction, axial view')
+
+plt.subplot(132)
+plt.imshow(FBPrec[:,sliceSel,:],vmin=0, vmax=max_val)
+plt.title('3D FBP Reconstruction, coronal view')
+
+plt.subplot(133)
+plt.imshow(FBPrec[:,:,sliceSel],vmin=0, vmax=max_val)
+plt.title('3D FBP Reconstruction, sagittal view')
+plt.show()
+
+print("Min {} and Max {} of the volume".format(np.min(FBPrec), np.max(FBPrec)))
+del(FBPrec)
+#%%
+# It is recommend to re-run twice in order to get the optimal time
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print ("%%%%%%%%%Reconstructing with 3D FBP-CuPy method %%%%%%%%%%%%")
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+tic=timeit.default_timer()
+FBPrec = RectoolsDIR.FBP3D_cupy(projData3D_analyt) #perform FBP
+toc=timeit.default_timer()
+Run_time = toc - tic
+print("FBP 3D reconstruction with FFT filtering using CuPy (GPU) in {} seconds".format(Run_time))
+
+sliceSel = int(0.5*N_size)
+max_val = 1
+plt.figure() 
+plt.subplot(131)
+plt.imshow(FBPrec[sliceSel,:,:],vmin=0, vmax=max_val)
+plt.title('3D FBP Reconstruction, axial view')
+
+plt.subplot(132)
+plt.imshow(FBPrec[:,sliceSel,:],vmin=0, vmax=max_val)
+plt.title('3D FBP Reconstruction, coronal view')
+
+plt.subplot(133)
+plt.imshow(FBPrec[:,:,sliceSel],vmin=0, vmax=max_val)
+plt.title('3D FBP Reconstruction, sagittal view')
+plt.show()
+
+print("Min {} and Max {} of the volume".format(np.min(FBPrec), np.max(FBPrec)))
+del(FBPrec)
+#%%

src/Python/tomobar/methodsDIR.py

@@ -1,7 +1,6 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
-"""
-A reconstruction class for direct reconstructon methods.
+"""A reconstruction class for direct reconstructon methods.
 
 -- Fourier Slice Theorem reconstruction (adopted from Tim Day's code)
 -- Forward/Backward projection (ASTRA)
@@ -14,9 +13,9 @@
 from tomobar.supp.astraOP import AstraTools
 from tomobar.supp.astraOP import AstraTools3D
 from tomobar.supp.astraOP import parse_device_argument
+import scipy.fftpack
 
 def filtersinc3D(projection3D):
-    import scipy.fftpack
     # applies filters to __3D projection data__ in order to achieve FBP
     # Data format [DetectorVert, Projections, DetectorHoriz]
     # adopted from Matlabs code by  Waqas Akram
@@ -26,7 +25,7 @@ def filtersinc3D(projection3D):
     #roll off at high frequencies.
     a = 1.1
     [DetectorsLengthV, projectionsNum, DetectorsLengthH] = np.shape(projection3D)
-    w =  np.linspace(-np.pi,np.pi-(2*np.pi)/DetectorsLengthH, DetectorsLengthH,dtype='float32')
+    w = np.linspace(-np.pi,np.pi-(2*np.pi)/DetectorsLengthH, DetectorsLengthH,dtype='float32')
 
     rn1 = np.abs(2.0/a*np.sin(a*w/2.0))
     rn2 = np.sin(a*w/2.0)
@@ -36,27 +35,53 @@ def filtersinc3D(projection3D):
     r = rn1*(np.dot(rn2, np.linalg.pinv(rd_c)))**2
     multiplier = (1.0/projectionsNum)
     f = scipy.fftpack.fftshift(r)
-    filtered = np.zeros(np.shape(projection3D))
+    # making a 2d filter for projection
+    f_2d = np.zeros((DetectorsLengthV,DetectorsLengthH), dtype='float32')
+    f_2d[0::,:] = f
+    filtered = np.zeros(np.shape(projection3D), dtype='float32')
+
+    for i in range(0,projectionsNum):
+        IMG = scipy.fftpack.fft2(projection3D[:,i,:])
+        fimg = IMG*f_2d
+        filtered[:,i,:] = np.real(scipy.fftpack.ifft2(fimg))
+    return multiplier*filtered
+
+def filtersinc3D_cupy(projection3D):
+    import cupy as cp
+    # applies filters to __3D projection data__ in order to achieve FBP (using CuPy)
+    # Input: projection data must be a CuPy object
+    # Output: Filtered GPU stored CuPy projection data array
+    a = 1.1
+    [DetectorsLengthV, projectionsNum, DetectorsLengthH] = np.shape(projection3D)
+    w = np.linspace(-np.pi,np.pi-(2*np.pi)/DetectorsLengthH, DetectorsLengthH,dtype='float32')
+
+    rn1 = np.abs(2.0/a*np.sin(a*w/2.0))
+    rn2 = np.sin(a*w/2.0)
+    rd = (a*w)/2.0
+    rd_c = np.zeros([1,DetectorsLengthH])
+    rd_c[0,:] = rd
+    r = rn1*(np.dot(rn2, np.linalg.pinv(rd_c)))**2
+    multiplier = (1.0/projectionsNum)
+    f = scipy.fftpack.fftshift(r)
+    # making a 2d filter for projection 
+    f_2d = np.zeros((DetectorsLengthV,DetectorsLengthH), dtype='float32')
+    f_2d[0::,:] = f
+    filter_gpu = cp.asarray(f_2d)
+
+    filtered = cp.zeros(cp.shape(projection3D), dtype='float32')
+
+    for i in range(0,projectionsNum):
+        IMG = cp.fft.fft2(projection3D[:,i,:])
+        fimg = IMG*filter_gpu
+        filtered[:,i,:] = cp.real(cp.fft.ifft2(fimg))
+    return multiplier*filtered
 
-    for j in range(0,DetectorsLengthV):
-        for i in range(0,projectionsNum):
-            IMG = scipy.fftpack.fft(projection3D[j,i,:])
-            fimg = IMG*f
-            filtered[j,i,:] = multiplier*np.real(scipy.fftpack.ifft(fimg))
-    return np.float32(filtered)
 
 def filtersinc2D(sinogram):
-    import scipy.fftpack
-    # applies filters toa sinogram in order to achieve FBP
-    # Data format [Projections, DetectorHoriz]
-    # adopted from Matlabs code by  Waqas Akram
-    #"a":	This parameter varies the filter magnitude response.
-    #When "a" is very small (a<<1), the response approximates |w|
-    #As "a" is increased, the filter response starts to
-    #roll off at high frequencies.
+    # applies filters to __2D projection data__ in order to achieve FBP
     a = 1.1
     [projectionsNum, DetectorsLengthH] = np.shape(sinogram)
-    w =  np.linspace(-np.pi,np.pi-(2*np.pi)/DetectorsLengthH, DetectorsLengthH,dtype='float32')
+    w = np.linspace(-np.pi,np.pi-(2*np.pi)/DetectorsLengthH, DetectorsLengthH,dtype='float32')
 
     rn1 = np.abs(2.0/a*np.sin(a*w/2.0))
     rn2 = np.sin(a*w/2.0)
@@ -92,7 +117,7 @@ def __init__(self,
         self.DetectorsDimV = DetectorsDimV
         self.AnglesVec = AnglesVec
         self.CenterRotOffset = CenterRotOffset
-        self.OS_number = 1       
+        self.OS_number = 1
         self.device_projector, self.GPUdevice_index = parse_device_argument(device_projector)
 
         if DetectorsDimV is None:
@@ -220,4 +245,21 @@ def FBP(self, sinogram):
             Atools = AstraTools3D(self.DetectorsDimH, self.DetectorsDimV, self.AnglesVec, self.CenterRotOffset, self.ObjSize, self.OS_number , self.device_projector, self.GPUdevice_index) # initiate 3D ASTRA class object
             filtered_sino = filtersinc3D(sinogram) # filtering sinogram
             FBP_rec = Atools.backproj(filtered_sino) # backproject
+            del(filtered_sino)
+        return FBP_rec
+    def FBP3D_cupy(self, projection3d):
+        """
+        3D reconstruction using CuPy API to keep the data in GPU memory for FFT filtering
+        and backprojection
+        """
+        Atools = AstraTools3D(self.DetectorsDimH, self.DetectorsDimV, self.AnglesVec, self.CenterRotOffset, self.ObjSize, self.OS_number , self.device_projector, self.GPUdevice_index) # initiate 3D ASTRA class object
+        import cupy as cp
+        projection3d_gpu = cp.asarray(projection3d) # move data to GPU
+        projection3d_filtered_gpu = filtersinc3D_cupy(projection3d_gpu) # filter data on a GPU and keep the result on gpu
+        del(projection3d_gpu) # delete the raw data
+        cp._default_memory_pool.free_all_blocks()
+        # Using GPULink Astra capability to pass a pointer to GPU memory
+        FBP_rec = Atools.backprojCuPy(projection3d_filtered_gpu) # backproject while keeping data on a GPU
+        del(projection3d_filtered_gpu)
+        cp._default_memory_pool.free_all_blocks()
         return FBP_rec

src/Python/tomobar/supp/astraOP.py

@@ -319,10 +319,10 @@ def __init__(self, DetectorsDimH, DetectorsDimV, AnglesVec, CenterRotOffset, Obj
     def runAstraRecon(self, proj_data, method, iterations, os_index):
         # set ASTRA configuration for 3D reconstructor
         if self.OS_number != 1:
-            # traditional full data parallel beam projection geometry
+            # ordered-subsets
             proj_id = astra.data3d.create("-sino", self.proj_geom_OS[os_index], proj_data)
         else:
-            # ordered-subsets
+            # traditional full data parallel beam projection geometry
             proj_id = astra.data3d.create("-sino", self.proj_geom, proj_data)
 
         # Create a data object for the reconstruction
@@ -345,18 +345,53 @@ def runAstraRecon(self, proj_data, method, iterations, os_index):
         astra.data3d.delete(rec_id)
         astra.data3d.delete(proj_id)
         return recon_volume
+
+    def runAstraReconCuPy(self, proj_data, method, iterations, os_index):
+        # set ASTRA configuration for 3D reconstructor using CuPy arrays
+        if self.OS_number != 1:
+            # ordered-subsets
+            projector_id = astra.create_projector('cuda3d', self.proj_geom_OS[os_index], self.vol_geom)
+        else:
+            # traditional full data parallel beam projection geometry
+            projector_id = astra.create_projector('cuda3d', self.proj_geom, self.vol_geom)
+
 
+        proj_link = astra.data3d.GPULink(proj_data.data.ptr, *proj_data.shape[::-1],4*proj_data.shape[2])
+        proj_id = astra.data3d.link('-proj3d', self.proj_geom, proj_link)
+
+        # Create a data object for the reconstruction
+        rec_id = astra.data3d.create('-vol', self.vol_geom)
+
+        # Create algorithm object
+        cfg = astra.astra_dict(method)
+        cfg['option'] = {'GPUindex': self.GPUdevice_index}
+        cfg['ProjectionDataId'] = proj_id
+        cfg['ReconstructionDataId'] = rec_id
+        cfg['ProjectorId'] = projector_id
+
+        # Create the algorithm object from the configuration structure
+        alg_id = astra.algorithm.create(cfg)
+        astra.algorithm.run(alg_id, iterations)
+
+        # Get the result (this will copy the data back to the host)
+        recon_volume = astra.data3d.get(rec_id)
+
+        astra.algorithm.delete(alg_id)
+        astra.data3d.delete(rec_id)
+        astra.data3d.delete(proj_id)
+        return recon_volume
+
     def runAstraProj(self, volume_data, os_index):
          # set ASTRA configuration for 3D projector
         if isinstance(volume_data, np.ndarray):
             volume_id = astra.data3d.link('-vol', self.vol_geom, volume_data)
         else:
             volume_id = volume_data
         if self.OS_number != 1:
-            # traditional full data parallel beam projection geometry
+            # ordered-subsets
             proj_id = astra.data3d.create('-sino', self.proj_geom_OS[os_index], 0)
         else:
-            # ordered-subsets
+            # traditional full data parallel beam projection geometry
             proj_id = astra.data3d.create('-sino', self.proj_geom, 0)
 
         # Create algorithm object
@@ -441,6 +476,8 @@ def forwproj(self, object3D):
         return Astra3D.runAstraProj(self, object3D, None) # 3D forward projection
     def backproj(self, proj_data):
         return Astra3D.runAstraRecon(self, proj_data, 'BP3D_CUDA', 1, None) # 3D backprojection
+    def backprojCuPy(self, proj_data):
+        return Astra3D.runAstraReconCuPy(self, proj_data, 'BP3D_CUDA', 1, None) # 3D backprojection using CuPy array
     def sirt3D(self, proj_data, iterations):
         return Astra3D.runAstraRecon(self, proj_data, 'SIRT3D_CUDA', iterations, None) #3D SIRT reconstruction
     def cgls3D(self, proj_data, iterations):