diff --git a/Demos/Python/DemoFBP_CuPy_3D.py b/Demos/Python/DemoFBP_CuPy_3D.py new file mode 100644 index 000000000..deac199ee --- /dev/null +++ b/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) +#%% \ No newline at end of file diff --git a/src/Python/tomobar/methodsDIR.py b/src/Python/tomobar/methodsDIR.py index 6c222a23d..77f4896d3 100644 --- a/src/Python/tomobar/methodsDIR.py +++ b/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 \ No newline at end of file diff --git a/src/Python/tomobar/supp/astraOP.py b/src/Python/tomobar/supp/astraOP.py index 3f477bbce..3226a294c 100644 --- a/src/Python/tomobar/supp/astraOP.py +++ b/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,7 +345,42 @@ 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): @@ -353,10 +388,10 @@ def runAstraProj(self, volume_data, os_index): 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):