From 787b534643d5b4cad4e6f8d9c4b524b52d804348 Mon Sep 17 00:00:00 2001 From: dkazanc Date: Mon, 18 Feb 2019 14:48:29 +0000 Subject: TGV3D_GPU added, demos updated fpr Python/Matlab --- Wrappers/Python/conda-recipe/run_test.py | 2 +- Wrappers/Python/demos/demo_cpu_regularisers.py | 4 +- Wrappers/Python/demos/demo_cpu_regularisers3D.py | 54 ++++++++++++++++++++-- .../Python/demos/demo_cpu_vs_gpu_regularisers.py | 4 +- Wrappers/Python/demos/demo_gpu_regularisers.py | 4 +- Wrappers/Python/demos/demo_gpu_regularisers3D.py | 51 +++++++++++++++++++- Wrappers/Python/src/cpu_regularisers.pyx | 35 ++++++++++---- Wrappers/Python/src/gpu_regularisers.pyx | 39 ++++++++++++---- 8 files changed, 164 insertions(+), 29 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/conda-recipe/run_test.py b/Wrappers/Python/conda-recipe/run_test.py index 37b9dcc..21f3216 100755 --- a/Wrappers/Python/conda-recipe/run_test.py +++ b/Wrappers/Python/conda-recipe/run_test.py @@ -303,7 +303,7 @@ class TestRegularisers(unittest.TestCase): 'input' : u0,\ 'regularisation_parameter':0.04, \ 'alpha1':1.0,\ - 'alpha0':0.7,\ + 'alpha0':2.0,\ 'number_of_iterations' :250 ,\ 'LipshitzConstant' :12 ,\ } diff --git a/Wrappers/Python/demos/demo_cpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_regularisers.py index 859b633..e6befa9 100644 --- a/Wrappers/Python/demos/demo_cpu_regularisers.py +++ b/Wrappers/Python/demos/demo_cpu_regularisers.py @@ -225,8 +225,8 @@ pars = {'algorithm' : TGV, \ 'input' : u0,\ 'regularisation_parameter':0.04, \ 'alpha1':1.0,\ - 'alpha0':0.7,\ - 'number_of_iterations' :250 ,\ + 'alpha0':2.0,\ + 'number_of_iterations' :1350 ,\ 'LipshitzConstant' :12 ,\ } diff --git a/Wrappers/Python/demos/demo_cpu_regularisers3D.py b/Wrappers/Python/demos/demo_cpu_regularisers3D.py index c42c37b..2d2fc22 100644 --- a/Wrappers/Python/demos/demo_cpu_regularisers3D.py +++ b/Wrappers/Python/demos/demo_cpu_regularisers3D.py @@ -12,7 +12,7 @@ import matplotlib.pyplot as plt import numpy as np import os import timeit -from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, LLT_ROF, FGP_dTV, NDF, Diff4th +from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -68,7 +68,7 @@ Im2[:,0:M] = Im[:,0:M] Im = Im2 del Im2 """ -slices = 20 +slices = 15 noisyVol = np.zeros((slices,N,M),dtype='float32') noisyRef = np.zeros((slices,N,M),dtype='float32') @@ -96,7 +96,7 @@ pars = {'algorithm': ROF_TV, \ 'input' : noisyVol,\ 'regularisation_parameter':0.04,\ 'number_of_iterations': 500,\ - 'time_marching_parameter': 0.0025 + 'time_marching_parameter': 0.0025 } print ("#############ROF TV CPU####################") start_time = timeit.default_timer() @@ -262,6 +262,54 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(lltrof_cpu3D[10,:,:], cmap="gray") plt.title('{}'.format('Recovered volume on the CPU using LLT-ROF')) +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("_______________TGV (3D)_________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure() +plt.suptitle('Performance of TGV regulariser using the CPU') +a=fig.add_subplot(1,2,1) +a.set_title('Noisy Image') +imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray") + +# set parameters +pars = {'algorithm' : TGV, \ + 'input' : noisyVol,\ + 'regularisation_parameter':0.04, \ + 'alpha1':1.0,\ + 'alpha0':2.0,\ + 'number_of_iterations' :250 ,\ + 'LipshitzConstant' :12 ,\ + } + +print ("#############TGV CPU####################") +start_time = timeit.default_timer() +tgv_cpu3D = TGV(pars['input'], + pars['regularisation_parameter'], + pars['alpha1'], + pars['alpha0'], + pars['number_of_iterations'], + pars['LipshitzConstant'],'cpu') + + +rms = rmse(idealVol, tgv_cpu3D) +pars['rmse'] = rms + +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(1,2,2) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) +# place a text box in upper left in axes coords +a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(tgv_cpu3D[10,:,:], cmap="gray") +plt.title('{}'.format('Recovered volume on the CPU using TGV')) + #%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("________________NDF (3D)___________________") diff --git a/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py index 275e844..230a761 100644 --- a/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py +++ b/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py @@ -323,8 +323,8 @@ pars = {'algorithm' : TGV, \ 'input' : u0,\ 'regularisation_parameter':0.04, \ 'alpha1':1.0,\ - 'alpha0':0.7,\ - 'number_of_iterations' :250 ,\ + 'alpha0':2.0,\ + 'number_of_iterations' :400 ,\ 'LipshitzConstant' :12 ,\ } diff --git a/Wrappers/Python/demos/demo_gpu_regularisers.py b/Wrappers/Python/demos/demo_gpu_regularisers.py index 9115494..e1c6575 100644 --- a/Wrappers/Python/demos/demo_gpu_regularisers.py +++ b/Wrappers/Python/demos/demo_gpu_regularisers.py @@ -223,8 +223,8 @@ pars = {'algorithm' : TGV, \ 'input' : u0,\ 'regularisation_parameter':0.04, \ 'alpha1':1.0,\ - 'alpha0':0.7,\ - 'number_of_iterations' :250 ,\ + 'alpha0':2.0,\ + 'number_of_iterations' :1250 ,\ 'LipshitzConstant' :12 ,\ } diff --git a/Wrappers/Python/demos/demo_gpu_regularisers3D.py b/Wrappers/Python/demos/demo_gpu_regularisers3D.py index cda2847..b6058d2 100644 --- a/Wrappers/Python/demos/demo_gpu_regularisers3D.py +++ b/Wrappers/Python/demos/demo_gpu_regularisers3D.py @@ -12,7 +12,7 @@ import matplotlib.pyplot as plt import numpy as np import os import timeit -from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, LLT_ROF, FGP_dTV, NDF, Diff4th +from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -67,7 +67,7 @@ Im = Im2 del Im2 """ -#%% + slices = 20 filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif") @@ -266,6 +266,53 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(lltrof_gpu3D[10,:,:], cmap="gray") plt.title('{}'.format('Recovered volume on the GPU using LLT-ROF')) +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("_______________TGV (3D)_________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure() +plt.suptitle('Performance of TGV regulariser using the GPU') +a=fig.add_subplot(1,2,1) +a.set_title('Noisy Image') +imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray") + +# set parameters +pars = {'algorithm' : TGV, \ + 'input' : noisyVol,\ + 'regularisation_parameter':0.04, \ + 'alpha1':1.0,\ + 'alpha0':2.0,\ + 'number_of_iterations' :600 ,\ + 'LipshitzConstant' :12 ,\ + } + +print ("#############TGV GPU####################") +start_time = timeit.default_timer() +tgv_gpu3D = TGV(pars['input'], + pars['regularisation_parameter'], + pars['alpha1'], + pars['alpha0'], + pars['number_of_iterations'], + pars['LipshitzConstant'],'gpu') + + +rms = rmse(idealVol, tgv_gpu3D) +pars['rmse'] = rms + +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(1,2,2) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) +# place a text box in upper left in axes coords +a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(tgv_gpu3D[10,:,:], cmap="gray") +plt.title('{}'.format('Recovered volume on the GPU using TGV')) #%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________NDF-TV (3D)_________________") diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx index 7d57ed1..11a0617 100644 --- a/Wrappers/Python/src/cpu_regularisers.pyx +++ b/Wrappers/Python/src/cpu_regularisers.pyx @@ -22,7 +22,7 @@ cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ); cdef extern float SB_TV_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int printM, int dimX, int dimY, int dimZ); cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, int dimX, int dimY, int dimZ); -cdef extern float TGV_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY); +cdef extern float TGV_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY, int dimZ); cdef extern float Diffusion_CPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ); cdef extern float Diffus4th_CPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ); cdef extern float TNV_CPU_main(float *Input, float *u, float lambdaPar, int maxIter, float tol, int dimX, int dimY, int dimZ); @@ -202,12 +202,8 @@ def TGV_CPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, Lip return TGV_2D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst) elif inputData.ndim == 3: - shape = inputData.shape - out = inputData.copy() - for i in range(shape[0]): - out[i,:,:] = TGV_2D(inputData[i,:,:], regularisation_parameter, - alpha1, alpha0, iterations, LipshitzConst) - return out + return TGV_3D(inputData, regularisation_parameter, alpha1, alpha0, + iterations, LipshitzConst) def TGV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameter, @@ -229,7 +225,30 @@ def TGV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, alpha0, iterationsNumb, LipshitzConst, - dims[1],dims[0]) + dims[1],dims[0],1) + return outputData +def TGV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularisation_parameter, + float alpha1, + float alpha0, + int iterationsNumb, + float LipshitzConst): + + cdef long dims[3] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + dims[2] = inputData.shape[2] + + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \ + np.zeros([dims[0], dims[1], dims[2]], dtype='float32') + + #/* Run TGV iterations for 3D data */ + TGV_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, + alpha1, + alpha0, + iterationsNumb, + LipshitzConst, + dims[2], dims[1], dims[0]) return outputData #***************************************************************# diff --git a/Wrappers/Python/src/gpu_regularisers.pyx b/Wrappers/Python/src/gpu_regularisers.pyx index 47a6149..b52f669 100644 --- a/Wrappers/Python/src/gpu_regularisers.pyx +++ b/Wrappers/Python/src/gpu_regularisers.pyx @@ -23,7 +23,7 @@ CUDAErrorMessage = 'CUDA error' cdef extern int TV_ROF_GPU_main(float* Input, float* Output, float lambdaPar, int iter, float tau, int N, int M, int Z); cdef extern int TV_FGP_GPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z); cdef extern int TV_SB_GPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int printM, int N, int M, int Z); -cdef extern int TGV_GPU_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY); +cdef extern int TGV_GPU_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY, int dimZ); cdef extern int LLT_ROF_GPU_main(float *Input, float *Output, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, int N, int M, int Z); cdef extern int NonlDiff_GPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int N, int M, int Z); cdef extern int dTV_FGP_GPU_main(float *Input, float *InputRef, float *Output, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int printM, int N, int M, int Z); @@ -102,12 +102,7 @@ def TGV_GPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, Lip if inputData.ndim == 2: return TGV2D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst) elif inputData.ndim == 3: - shape = inputData.shape - out = inputData.copy() - for i in range(shape[0]): - out[i,:,:] = TGV2D(inputData[i,:,:], regularisation_parameter, - alpha1, alpha0, iterations, LipshitzConst) - return out + return TGV3D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst) # Directional Total-variation Fast-Gradient-Projection (FGP) def dTV_FGP_GPU(inputData, refdata, @@ -393,7 +388,6 @@ def LLT_ROF_GPU3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, raise ValueError(CUDAErrorMessage); - #***************************************************************# #***************** Total Generalised Variation *****************# #***************************************************************# @@ -417,11 +411,38 @@ def TGV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, alpha0, iterationsNumb, LipshitzConst, - dims[1],dims[0])==0): + dims[1],dims[0], 1)==0): return outputData else: raise ValueError(CUDAErrorMessage); +def TGV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularisation_parameter, + float alpha1, + float alpha0, + int iterationsNumb, + float LipshitzConst): + + cdef long dims[3] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + dims[2] = inputData.shape[2] + + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \ + np.zeros([dims[0],dims[1],dims[2]], dtype='float32') + + # Running CUDA code here + if (TGV_GPU_main( + &inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, + alpha1, + alpha0, + iterationsNumb, + LipshitzConst, + dims[2], dims[1], dims[0])==0): + return outputData; + else: + raise ValueError(CUDAErrorMessage); + #****************************************************************# #**************Directional Total-variation FGP ******************# -- cgit v1.2.3