diff options
Diffstat (limited to 'Wrappers/Python')
| -rw-r--r-- | Wrappers/Python/ccpi/filters/regularisers.py | 14 | ||||
| -rw-r--r-- | Wrappers/Python/demos/demo_cpu_regularisers.py | 50 | ||||
| -rw-r--r-- | Wrappers/Python/demos/demo_cpu_regularisers3D.py | 67 | ||||
| -rw-r--r-- | Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py | 86 | ||||
| -rw-r--r-- | Wrappers/Python/demos/demo_gpu_regularisers.py | 49 | ||||
| -rw-r--r-- | Wrappers/Python/demos/demo_gpu_regularisers3D.py | 66 | ||||
| -rw-r--r-- | Wrappers/Python/setup-regularisers.py.in | 1 | ||||
| -rw-r--r-- | Wrappers/Python/src/cpu_regularisers.pyx | 45 | ||||
| -rw-r--r-- | Wrappers/Python/src/gpu_regularisers.pyx | 47 |
9 files changed, 397 insertions, 28 deletions
diff --git a/Wrappers/Python/ccpi/filters/regularisers.py b/Wrappers/Python/ccpi/filters/regularisers.py index 0e435a6..52c7974 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 import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, Diff4th_CPU, TGV_CPU -from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU, Diff4th_GPU, TGV_GPU +from ccpi.filters.cpu_regularisers import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, Diff4th_CPU, TGV_CPU, LLT_ROF_CPU +from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU, Diff4th_GPU, TGV_GPU, LLT_ROF_GPU from ccpi.filters.cpu_regularisers import NDF_INPAINT_CPU, NVM_INPAINT_CPU def ROF_TV(inputData, regularisation_parameter, iterations, @@ -147,7 +147,15 @@ def TGV(inputData, regularisation_parameter, alpha1, alpha0, iterations, else: raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ .format(device)) - +def LLT_ROF(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, + time_marching_parameter, device='cpu'): + if device == 'cpu': + return LLT_ROF_CPU(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter) + elif device == 'gpu': + return LLT_ROF_GPU(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter) + else: + raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ + .format(device)) def NDF_INP(inputData, maskData, regularisation_parameter, edge_parameter, iterations, time_marching_parameter, penalty_type): return NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter, diff --git a/Wrappers/Python/demos/demo_cpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_regularisers.py index 5c20244..b94f11c 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, TGV, FGP_dTV, TNV, NDF, DIFF4th +from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, TNV, NDF, DIFF4th from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -256,6 +256,54 @@ imgplot = plt.imshow(tgv_cpu, cmap="gray") plt.title('{}'.format('CPU results')) #%% + +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("______________LLT- ROF (2D)________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure() +plt.suptitle('Performance of LLT-ROF 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' : LLT_ROF, \ + 'input' : u0,\ + 'regularisation_parameterROF':0.04, \ + 'regularisation_parameterLLT':0.01, \ + 'number_of_iterations' :500 ,\ + 'time_marching_parameter' :0.0025 ,\ + } + +print ("#############LLT- ROF CPU####################") +start_time = timeit.default_timer() +lltrof_cpu = LLT_ROF(pars['input'], + pars['regularisation_parameterROF'], + pars['regularisation_parameterLLT'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'cpu') + +rms = rmse(Im, lltrof_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(lltrof_cpu, cmap="gray") +plt.title('{}'.format('CPU results')) + +#%% + + print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("________________NDF (2D)___________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") diff --git a/Wrappers/Python/demos/demo_cpu_regularisers3D.py b/Wrappers/Python/demos/demo_cpu_regularisers3D.py index 8ee157e..9c28de1 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, FGP_dTV, NDF, DIFF4th +from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, LLT_ROF, FGP_dTV, NDF, DIFF4th from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -85,7 +85,7 @@ print ("_______________ROF-TV (3D)_________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(1) +fig = plt.figure() 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') @@ -120,13 +120,13 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(rof_cpu3D[10,:,:], cmap="gray") plt.title('{}'.format('Recovered volume on the CPU using ROF-TV')) - +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________FGP-TV (3D)__________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(2) +fig = plt.figure() plt.suptitle('Performance of FGP-TV regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -170,12 +170,13 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(fgp_cpu3D[10,:,:], cmap="gray") plt.title('{}'.format('Recovered volume on the CPU using FGP-TV')) +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________SB-TV (3D)_________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(3) +fig = plt.figure() plt.suptitle('Performance of SB-TV regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -216,12 +217,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 ("_______________LLT-ROF (3D)_________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure() +plt.suptitle('Performance of LLT-ROF 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' : LLT_ROF, \ + 'input' : noisyVol,\ + 'regularisation_parameterROF':0.04, \ + 'regularisation_parameterLLT':0.015, \ + 'number_of_iterations' :300 ,\ + 'time_marching_parameter' :0.0025 ,\ + } + +print ("#############LLT ROF CPU####################") +start_time = timeit.default_timer() +lltrof_cpu3D = LLT_ROF(pars['input'], + pars['regularisation_parameterROF'], + pars['regularisation_parameterLLT'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'cpu') + +rms = rmse(idealVol, lltrof_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(lltrof_cpu3D[10,:,:], cmap="gray") +plt.title('{}'.format('Recovered volume on the CPU using LLT-ROF')) + +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("________________NDF (3D)___________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(4) +fig = plt.figure() plt.suptitle('Performance of NDF regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy volume') @@ -262,13 +309,13 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(ndf_cpu3D[10,:,:], cmap="gray") plt.title('{}'.format('Recovered volume on the CPU using NDF iterations')) - +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("___Anisotropic Diffusion 4th Order (2D)____") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(5) +fig = plt.figure() plt.suptitle('Performance of Diff4th regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy volume') @@ -307,13 +354,13 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(diff4th_cpu3D[10,:,:], cmap="gray") plt.title('{}'.format('Recovered volume on the CPU using DIFF4th iterations')) - +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________FGP-dTV (3D)__________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(6) +fig = plt.figure() 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 46b8ffc..e45dc40 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, TGV, 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): @@ -352,8 +352,7 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(tgv_cpu, cmap="gray") plt.title('{}'.format('CPU results')) - -print ("##############SB TV GPU##################") +print ("##############TGV GPU##################") start_time = timeit.default_timer() tgv_gpu = TGV(pars['input'], pars['regularisation_parameter'], @@ -392,6 +391,87 @@ else: print ("Arrays match") #%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("____________LLT-ROF bench___________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure() +plt.suptitle('Comparison of LLT-ROF 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' : LLT_ROF, \ + 'input' : u0,\ + 'regularisation_parameterROF':0.04, \ + 'regularisation_parameterLLT':0.01, \ + 'number_of_iterations' :500 ,\ + 'time_marching_parameter' :0.0025 ,\ + } + +print ("#############LLT- ROF CPU####################") +start_time = timeit.default_timer() +lltrof_cpu = LLT_ROF(pars['input'], + pars['regularisation_parameterROF'], + pars['regularisation_parameterLLT'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'cpu') + +rms = rmse(Im, lltrof_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(lltrof_cpu, cmap="gray") +plt.title('{}'.format('CPU results')) + +print ("#############LLT- ROF GPU####################") +start_time = timeit.default_timer() +lltrof_gpu = LLT_ROF(pars['input'], + pars['regularisation_parameterROF'], + pars['regularisation_parameterLLT'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'gpu') + +rms = rmse(Im, lltrof_gpu) +pars['rmse'] = rms +pars['algorithm'] = LLT_ROF +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(lltrof_gpu, cmap="gray") +plt.title('{}'.format('GPU results')) + +print ("--------Compare the results--------") +tolerance = 1e-05 +diff_im = np.zeros(np.shape(lltrof_gpu)) +diff_im = abs(lltrof_cpu - lltrof_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 ("_______________NDF bench___________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") diff --git a/Wrappers/Python/demos/demo_gpu_regularisers.py b/Wrappers/Python/demos/demo_gpu_regularisers.py index 792a019..de0cbde 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, TGV, 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): @@ -254,6 +254,53 @@ imgplot = plt.imshow(tgv_gpu, cmap="gray") plt.title('{}'.format('GPU results')) #%% + +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("______________LLT- ROF (2D)________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure() +plt.suptitle('Performance of LLT-ROF 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' : LLT_ROF, \ + 'input' : u0,\ + 'regularisation_parameterROF':0.04, \ + 'regularisation_parameterLLT':0.01, \ + 'number_of_iterations' :500 ,\ + 'time_marching_parameter' :0.0025 ,\ + } + +print ("#############LLT- ROF GPU####################") +start_time = timeit.default_timer() +lltrof_gpu = LLT_ROF(pars['input'], + pars['regularisation_parameterROF'], + pars['regularisation_parameterLLT'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'gpu') + + +rms = rmse(Im, lltrof_gpu) +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(lltrof_gpu, cmap="gray") +plt.title('{}'.format('GPU results')) + +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________NDF regulariser_____________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") diff --git a/Wrappers/Python/demos/demo_gpu_regularisers3D.py b/Wrappers/Python/demos/demo_gpu_regularisers3D.py index 13c4e7b..d5f9a39 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, FGP_dTV, NDF, DIFF4th +from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, LLT_ROF, FGP_dTV, NDF, DIFF4th from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -86,12 +86,13 @@ for i in range (slices): noisyRef[i,:,:] = Im + np.random.normal(loc = 0 , scale = 0.01 * Im , size = np.shape(Im)) idealVol[i,:,:] = Im +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________ROF-TV (3D)_________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(1) +fig = plt.figure() 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') @@ -125,13 +126,13 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) imgplot = plt.imshow(rof_gpu3D[10,:,:], cmap="gray") plt.title('{}'.format('Recovered volume on the GPU using ROF-TV')) - +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________FGP-TV (3D)__________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(2) +fig = plt.figure() plt.suptitle('Performance of FGP-TV regulariser using the GPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -174,12 +175,13 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(fgp_gpu3D[10,:,:], cmap="gray") plt.title('{}'.format('Recovered volume on the GPU using FGP-TV')) +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________SB-TV (3D)__________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(3) +fig = plt.figure() plt.suptitle('Performance of SB-TV regulariser using the GPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -219,14 +221,58 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) imgplot = plt.imshow(sb_gpu3D[10,:,:], cmap="gray") plt.title('{}'.format('Recovered volume on the GPU using SB-TV')) +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("_______________LLT-ROF (3D)_________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure() +plt.suptitle('Performance of LLT-ROF 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' : LLT_ROF, \ + 'input' : noisyVol,\ + 'regularisation_parameterROF':0.04, \ + 'regularisation_parameterLLT':0.015, \ + 'number_of_iterations' :300 ,\ + 'time_marching_parameter' :0.0025 ,\ + } + +print ("#############LLT ROF CPU####################") +start_time = timeit.default_timer() +lltrof_gpu3D = LLT_ROF(pars['input'], + pars['regularisation_parameterROF'], + pars['regularisation_parameterLLT'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'gpu') +rms = rmse(idealVol, lltrof_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(lltrof_gpu3D[10,:,:], cmap="gray") +plt.title('{}'.format('Recovered volume on the GPU using LLT-ROF')) + +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________NDF-TV (3D)_________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(4) +fig = plt.figure() plt.suptitle('Performance of NDF regulariser using the GPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -267,13 +313,13 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(ndf_gpu3D[10,:,:], cmap="gray") plt.title('{}'.format('Recovered volume on the GPU using NDF')) - +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("___Anisotropic Diffusion 4th Order (3D)____") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(5) +fig = plt.figure() plt.suptitle('Performance of DIFF4th regulariser using the GPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -312,13 +358,13 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(diff4_gpu3D[10,:,:], cmap="gray") plt.title('{}'.format('GPU results')) - +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________FGP-dTV (3D)________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(6) +fig = plt.figure() 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 89ebaf9..7108683 100644 --- a/Wrappers/Python/setup-regularisers.py.in +++ b/Wrappers/Python/setup-regularisers.py.in @@ -39,6 +39,7 @@ extra_include_dirs += [os.path.join(".." , ".." , "Core"), os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_ROF" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_SB" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TGV" ) , + os.path.join(".." , ".." , "Core", "regularisers_GPU" , "LLTROF" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "NDF" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "dTV_FGP" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "DIFF4th" ) , diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx index cf81bec..bf9c861 100644 --- a/Wrappers/Python/src/cpu_regularisers.pyx +++ b/Wrappers/Python/src/cpu_regularisers.pyx @@ -21,6 +21,7 @@ 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 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 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); @@ -222,7 +223,51 @@ def TGV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, LipshitzConst, dims[1],dims[0]) return outputData + +#***************************************************************# +#******************* ROF - LLT regularisation ******************# +#***************************************************************# +def LLT_ROF_CPU(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter): + if inputData.ndim == 2: + return LLT_ROF_2D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter) + elif inputData.ndim == 3: + return LLT_ROF_3D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter) + +def LLT_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularisation_parameterROF, + float regularisation_parameterLLT, + int iterations, + float time_marching_parameter): + + 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 ROF-LLT iterations for 2D data */ + LLT_ROF_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, dims[1],dims[0],1) + return outputData + +def LLT_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularisation_parameterROF, + float regularisation_parameterLLT, + int iterations, + float time_marching_parameter): + + 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 ROF-LLT iterations for 3D data */ + LLT_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, dims[2], dims[1], dims[0]) + return outputData + #****************************************************************# #**************Directional Total-variation FGP ******************# #****************************************************************# diff --git a/Wrappers/Python/src/gpu_regularisers.pyx b/Wrappers/Python/src/gpu_regularisers.pyx index 4a202d7..82d3e01 100644 --- a/Wrappers/Python/src/gpu_regularisers.pyx +++ b/Wrappers/Python/src/gpu_regularisers.pyx @@ -22,6 +22,7 @@ cdef extern void TV_ROF_GPU_main(float* Input, float* Output, float lambdaPar, i 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 TGV_GPU_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY); +cdef extern void LLT_ROF_GPU_main(float *Input, float *Output, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, 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); cdef extern void Diffus4th_GPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int N, int M, int Z); @@ -87,6 +88,12 @@ def TV_SB_GPU(inputData, tolerance_param, methodTV, printM) +# LLT-ROF model +def LLT_ROF_GPU(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter): + if inputData.ndim == 2: + return LLT_ROF_GPU2D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter) + elif inputData.ndim == 3: + return LLT_ROF_GPU3D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter) # Total Generilised Variation (TGV) def TGV_GPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst): if inputData.ndim == 2: @@ -324,6 +331,46 @@ def SBTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, return outputData #***************************************************************# +#************************ LLT-ROF model ************************# +#***************************************************************# +#************Joint LLT-ROF model for higher order **************# +def LLT_ROF_GPU2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularisation_parameterROF, + float regularisation_parameterLLT, + int iterations, + float time_marching_parameter): + + 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') + + # Running CUDA code here + LLT_ROF_GPU_main(&inputData[0,0], &outputData[0,0],regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, dims[1],dims[0],1); + return outputData + +def LLT_ROF_GPU3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularisation_parameterROF, + float regularisation_parameterLLT, + int iterations, + float time_marching_parameter): + + 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 + LLT_ROF_GPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, dims[2], dims[1], dims[0]); + return outputData + + +#***************************************************************# #***************** Total Generalised Variation *****************# #***************************************************************# def TGV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, |