diff options
Diffstat (limited to 'Wrappers/Python')
| -rw-r--r-- | Wrappers/Python/ccpi/filters/regularisers.py | 23 | ||||
| -rw-r--r-- | Wrappers/Python/demos/demo_cpu_regularisers.py | 106 | ||||
| -rw-r--r-- | Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py | 91 | ||||
| -rw-r--r-- | Wrappers/Python/demos/demo_gpu_regularisers.py | 114 | ||||
| -rw-r--r-- | Wrappers/Python/setup-regularisers.py.in | 1 | ||||
| -rw-r--r-- | Wrappers/Python/src/cpu_regularisers.pyx | 46 | ||||
| -rw-r--r-- | Wrappers/Python/src/gpu_regularisers.pyx | 67 |
7 files changed, 427 insertions, 21 deletions
diff --git a/Wrappers/Python/ccpi/filters/regularisers.py b/Wrappers/Python/ccpi/filters/regularisers.py index e6814e8..eec8c4d 100644 --- a/Wrappers/Python/ccpi/filters/regularisers.py +++ b/Wrappers/Python/ccpi/filters/regularisers.py @@ -2,8 +2,8 @@ script which assigns a proper device core function based on a flag ('cpu' or 'gpu') """ -from ccpi.filters.cpu_regularisers_cython import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU -from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU +from ccpi.filters.cpu_regularisers_cython import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU +from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU def ROF_TV(inputData, regularisation_parameter, iterations, time_marching_parameter,device='cpu'): @@ -91,3 +91,22 @@ def TNV(inputData, regularisation_parameter, iterations, tolerance_param): regularisation_parameter, iterations, tolerance_param) +def NDF(inputData, regularisation_parameter, edge_parameter, iterations, + time_marching_parameter, penalty_type, device='cpu'): + if device == 'cpu': + return NDF_CPU(inputData, + regularisation_parameter, + edge_parameter, + iterations, + time_marching_parameter, + penalty_type) + elif device == 'gpu': + return NDF_GPU(inputData, + regularisation_parameter, + edge_parameter, + iterations, + time_marching_parameter, + penalty_type) + else: + raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ + .format(device)) diff --git a/Wrappers/Python/demos/demo_cpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_regularisers.py index 7443b83..3567f91 100644 --- a/Wrappers/Python/demos/demo_cpu_regularisers.py +++ b/Wrappers/Python/demos/demo_cpu_regularisers.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, FGP_dTV, TNV +from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, FGP_dTV, TNV, NDF from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -190,11 +190,58 @@ plt.title('{}'.format('CPU results')) print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("_____________FGP-dTV (2D)__________________") +print ("________________NDF (2D)___________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot fig = plt.figure(4) +plt.suptitle('Performance of NDF regulariser using the CPU') +a=fig.add_subplot(1,2,1) +a.set_title('Noisy Image') +imgplot = plt.imshow(u0,cmap="gray") + +# set parameters +pars = {'algorithm' : NDF, \ + 'input' : u0,\ + 'regularisation_parameter':0.06, \ + 'edge_parameter':0.04,\ + 'number_of_iterations' :1000 ,\ + 'time_marching_parameter':0.025,\ + 'penalty_type':1 + } + +print ("#############NDF CPU################") +start_time = timeit.default_timer() +ndf_cpu = NDF(pars['input'], + pars['regularisation_parameter'], + pars['edge_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'], + pars['penalty_type'],'cpu') + +rms = rmse(Im, ndf_cpu) +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(ndf_cpu, cmap="gray") +plt.title('{}'.format('CPU results')) + + +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("_____________FGP-dTV (2D)__________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure(5) plt.suptitle('Performance of FGP-dTV regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -247,7 +294,7 @@ print ("__________Total nuclear Variation__________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(5) +fig = plt.figure(6) plt.suptitle('Performance of TNV regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -321,7 +368,7 @@ print ("_______________ROF-TV (3D)_________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(6) +fig = plt.figure(7) plt.suptitle('Performance of ROF-TV regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy 15th slice of a volume') @@ -361,7 +408,7 @@ print ("_______________FGP-TV (3D)__________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(7) +fig = plt.figure(8) plt.suptitle('Performance of FGP-TV regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -410,7 +457,7 @@ print ("_______________SB-TV (3D)_________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(8) +fig = plt.figure(9) plt.suptitle('Performance of SB-TV regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -451,13 +498,58 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(sb_cpu3D[10,:,:], cmap="gray") plt.title('{}'.format('Recovered volume on the CPU using SB-TV')) +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("________________NDF (3D)___________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure(10) +plt.suptitle('Performance of NDF regulariser using the CPU') +a=fig.add_subplot(1,2,1) +a.set_title('Noisy volume') +imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray") + +# set parameters +pars = {'algorithm' : NDF, \ + 'input' : noisyVol,\ + 'regularisation_parameter':0.06, \ + 'edge_parameter':0.04,\ + 'number_of_iterations' :1000 ,\ + 'time_marching_parameter':0.025,\ + 'penalty_type': 1 + } + +print ("#############NDF CPU################") +start_time = timeit.default_timer() +ndf_cpu3D = NDF(pars['input'], + pars['regularisation_parameter'], + pars['edge_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'], + pars['penalty_type']) + +rms = rmse(idealVol, ndf_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(ndf_cpu3D[10,:,:], cmap="gray") +plt.title('{}'.format('Recovered volume on the CPU using NDF iterations')) print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________FGP-dTV (3D)__________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(9) +fig = plt.figure(11) plt.suptitle('Performance of FGP-dTV regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') diff --git a/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py index d8e2da7..05db23e 100644 --- a/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py +++ b/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.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, FGP_dTV +from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, FGP_dTV, NDF from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -306,11 +306,98 @@ else: print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("____________FGP-dTV bench___________________") +print ("_______________NDF bench___________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot fig = plt.figure(4) +plt.suptitle('Comparison of NDF regulariser using CPU and GPU implementations') +a=fig.add_subplot(1,4,1) +a.set_title('Noisy Image') +imgplot = plt.imshow(u0,cmap="gray") + +# set parameters +pars = {'algorithm' : NDF, \ + 'input' : u0,\ + 'regularisation_parameter':0.06, \ + 'edge_parameter':0.04,\ + 'number_of_iterations' :1000 ,\ + 'time_marching_parameter':0.025,\ + 'penalty_type': 1 + } + +print ("#############NDF CPU####################") +start_time = timeit.default_timer() +ndf_cpu = NDF(pars['input'], + pars['regularisation_parameter'], + pars['edge_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'], + pars['penalty_type'],'cpu') + +rms = rmse(Im, ndf_cpu) +pars['rmse'] = rms + +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(1,4,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(ndf_cpu, cmap="gray") +plt.title('{}'.format('CPU results')) + + +print ("##############NDF GPU##################") +start_time = timeit.default_timer() +ndf_gpu = NDF(pars['input'], + pars['regularisation_parameter'], + pars['edge_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'], + pars['penalty_type'],'gpu') + +rms = rmse(Im, ndf_gpu) +pars['rmse'] = rms +pars['algorithm'] = NDF +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(1,4,3) + +# 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(ndf_gpu, cmap="gray") +plt.title('{}'.format('GPU results')) + +print ("--------Compare the results--------") +tolerance = 1e-05 +diff_im = np.zeros(np.shape(rof_cpu)) +diff_im = abs(ndf_cpu - ndf_gpu) +diff_im[diff_im > tolerance] = 1 +a=fig.add_subplot(1,4,4) +imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray") +plt.title('{}'.format('Pixels larger threshold difference')) +if (diff_im.sum() > 1): + print ("Arrays do not match!") +else: + print ("Arrays match") + + + +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("____________FGP-dTV bench___________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure(5) plt.suptitle('Comparison of FGP-dTV regulariser using CPU and GPU implementations') a=fig.add_subplot(1,4,1) a.set_title('Noisy Image') diff --git a/Wrappers/Python/demos/demo_gpu_regularisers.py b/Wrappers/Python/demos/demo_gpu_regularisers.py index 25d8d85..b873700 100644 --- a/Wrappers/Python/demos/demo_gpu_regularisers.py +++ b/Wrappers/Python/demos/demo_gpu_regularisers.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, FGP_dTV +from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, FGP_dTV, NDF from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -50,7 +50,7 @@ u0 = u0.astype('float32') u_ref = u_ref.astype('float32') print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("____________ROF-TV bench___________________") +print ("____________ROF-TV regulariser_____________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot @@ -92,7 +92,7 @@ plt.title('{}'.format('GPU results')) print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("____________FGP-TV bench___________________") +print ("____________FGP-TV regulariser_____________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot @@ -141,7 +141,7 @@ plt.title('{}'.format('GPU results')) print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("____________SB-TV bench___________________") +print ("____________SB-TV regulariser______________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot @@ -186,12 +186,60 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(sb_gpu, cmap="gray") plt.title('{}'.format('GPU results')) + print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("____________FGP-dTV bench___________________") +print ("_______________NDF regulariser_____________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot fig = plt.figure(4) +plt.suptitle('Performance of the NDF regulariser using the GPU') +a=fig.add_subplot(1,2,1) +a.set_title('Noisy Image') +imgplot = plt.imshow(u0,cmap="gray") + +# set parameters +pars = {'algorithm' : NDF, \ + 'input' : u0,\ + 'regularisation_parameter':0.06, \ + 'edge_parameter':0.04,\ + 'number_of_iterations' :1000 ,\ + 'time_marching_parameter':0.025,\ + 'penalty_type': 1 + } + +print ("##############NDF GPU##################") +start_time = timeit.default_timer() +ndf_gpu = NDF(pars['input'], + pars['regularisation_parameter'], + pars['edge_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'], + pars['penalty_type'],'gpu') + +rms = rmse(Im, ndf_gpu) +pars['rmse'] = rms +pars['algorithm'] = NDF +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(ndf_gpu, cmap="gray") +plt.title('{}'.format('GPU results')) + + +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("____________FGP-dTV bench___________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure(5) plt.suptitle('Performance of the FGP-dTV regulariser using the GPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -266,7 +314,7 @@ print ("_______________ROF-TV (3D)_________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(5) +fig = plt.figure(6) plt.suptitle('Performance of ROF-TV regulariser using the GPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy 15th slice of a volume') @@ -306,7 +354,7 @@ print ("_______________FGP-TV (3D)__________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(6) +fig = plt.figure(7) plt.suptitle('Performance of FGP-TV regulariser using the GPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -354,7 +402,7 @@ print ("_______________SB-TV (3D)__________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(7) +fig = plt.figure(8) plt.suptitle('Performance of SB-TV regulariser using the GPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -395,12 +443,60 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(sb_gpu3D[10,:,:], cmap="gray") plt.title('{}'.format('Recovered volume on the GPU using SB-TV')) + +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("_______________NDF-TV (3D)_________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure(9) +plt.suptitle('Performance of NDF 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' : NDF, \ + 'input' : noisyVol,\ + 'regularisation_parameter':0.06, \ + 'edge_parameter':0.04,\ + 'number_of_iterations' :1000 ,\ + 'time_marching_parameter':0.025,\ + 'penalty_type': 1 + } + +print ("#############NDF GPU####################") +start_time = timeit.default_timer() +ndf_gpu3D = NDF(pars['input'], + pars['regularisation_parameter'], + pars['edge_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'], + pars['penalty_type'],'gpu') + +rms = rmse(idealVol, ndf_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(ndf_gpu3D[10,:,:], cmap="gray") +plt.title('{}'.format('Recovered volume on the GPU using NDF')) + + print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________FGP-dTV (3D)________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(8) +fig = plt.figure(10) plt.suptitle('Performance of FGP-dTV regulariser using the GPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') diff --git a/Wrappers/Python/setup-regularisers.py.in b/Wrappers/Python/setup-regularisers.py.in index 0681cc4..b900efe 100644 --- a/Wrappers/Python/setup-regularisers.py.in +++ b/Wrappers/Python/setup-regularisers.py.in @@ -37,6 +37,7 @@ extra_include_dirs += [os.path.join(".." , ".." , "Core"), os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_FGP" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_ROF" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_SB" ) , + os.path.join(".." , ".." , "Core", "regularisers_GPU" , "NDF" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "dTV_FGP" ) , "."] diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx index abbf3b0..7ed8fa1 100644 --- a/Wrappers/Python/src/cpu_regularisers.pyx +++ b/Wrappers/Python/src/cpu_regularisers.pyx @@ -21,10 +21,10 @@ cimport numpy as np cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ); 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 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 TNV_CPU_main(float *Input, float *u, float lambdaPar, int maxIter, float tol, int dimX, int dimY, int dimZ); cdef extern float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ); - #****************************************************************# #********************** Total-variation ROF *********************# #****************************************************************# @@ -275,3 +275,47 @@ def TNV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, # Run TNV iterations for 3D (X,Y,Channels) data TNV_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, iterationsNumb, tolerance_param, dims[2], dims[1], dims[0]) return outputData +#****************************************************************# +#***************Nonlinear (Isotropic) Diffusion******************# +#****************************************************************# +def NDF_CPU(inputData, regularisation_parameter, edge_parameter, iterationsNumb,time_marching_parameter, penalty_type): + if inputData.ndim == 2: + return NDF_2D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type) + elif inputData.ndim == 3: + return NDF_3D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type) + +def NDF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularisation_parameter, + float edge_parameter, + int iterationsNumb, + float time_marching_parameter, + int penalty_type): + cdef long dims[2] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \ + np.zeros([dims[0],dims[1]], dtype='float32') + + # Run Nonlinear Diffusion iterations for 2D data + Diffusion_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[0], dims[1], 1) + return outputData + +def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularisation_parameter, + float edge_parameter, + int iterationsNumb, + float time_marching_parameter, + int penalty_type): + 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 Nonlinear Diffusion iterations for 3D data + Diffusion_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, 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 36eec95..b0775054 100644 --- a/Wrappers/Python/src/gpu_regularisers.pyx +++ b/Wrappers/Python/src/gpu_regularisers.pyx @@ -21,6 +21,7 @@ cimport numpy as np cdef extern void TV_ROF_GPU_main(float* Input, float* Output, float lambdaPar, int iter, float tau, int N, int M, int Z); cdef extern void 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 void 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 void 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 void 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); # Total-variation Rudin-Osher-Fatemi (ROF) @@ -114,6 +115,27 @@ def dTV_FGP_GPU(inputData, methodTV, nonneg, printM) +# Nonlocal Isotropic Diffusion (NDF) +def NDF_GPU(inputData, + regularisation_parameter, + edge_parameter, + iterations, + time_marching_parameter, + penalty_type): + if inputData.ndim == 2: + return NDF_GPU_2D(inputData, + regularisation_parameter, + edge_parameter, + iterations, + time_marching_parameter, + penalty_type) + elif inputData.ndim == 3: + return NDF_GPU_3D(inputData, + regularisation_parameter, + edge_parameter, + iterations, + time_marching_parameter, + penalty_type) #****************************************************************# #********************** Total-variation ROF *********************# #****************************************************************# @@ -336,3 +358,48 @@ def FGPdTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, printM, dims[2], dims[1], dims[0]); return outputData + +#****************************************************************# +#***************Nonlinear (Isotropic) Diffusion******************# +#****************************************************************# +def NDF_GPU_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularisation_parameter, + float edge_parameter, + int iterationsNumb, + float time_marching_parameter, + int penalty_type): + cdef long dims[2] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \ + np.zeros([dims[0],dims[1]], dtype='float32') + + #rangecheck = penalty_type < 1 and penalty_type > 3 + #if not rangecheck: +# raise ValueError('Choose penalty type as 1 for Huber, 2 - Perona-Malik, 3 - Tukey Biweight') + + # Run Nonlinear Diffusion iterations for 2D data + # Running CUDA code here + NonlDiff_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[0], dims[1], 1) + return outputData + +def NDF_GPU_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularisation_parameter, + float edge_parameter, + int iterationsNumb, + float time_marching_parameter, + int penalty_type): + 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 Nonlinear Diffusion iterations for 3D data + # Running CUDA code here + NonlDiff_GPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0]) + + return outputData |