diff options
Diffstat (limited to 'Wrappers/Python')
| -rw-r--r-- | Wrappers/Python/ccpi/filters/regularisers.py | 22 | ||||
| -rw-r--r-- | Wrappers/Python/demos/demo_cpu_regularisers.py | 311 | ||||
| -rw-r--r-- | Wrappers/Python/demos/demo_cpu_regularisers3D.py | 366 | ||||
| -rw-r--r-- | Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py | 92 | ||||
| -rw-r--r-- | Wrappers/Python/demos/demo_gpu_regularisers.py | 285 | ||||
| -rw-r--r-- | Wrappers/Python/demos/demo_gpu_regularisers3D.py | 367 | ||||
| -rw-r--r-- | Wrappers/Python/setup-regularisers.py.in | 1 | ||||
| -rw-r--r-- | Wrappers/Python/src/cpu_regularisers.pyx | 43 | ||||
| -rw-r--r-- | Wrappers/Python/src/gpu_regularisers.pyx | 60 |
9 files changed, 1028 insertions, 519 deletions
diff --git a/Wrappers/Python/ccpi/filters/regularisers.py b/Wrappers/Python/ccpi/filters/regularisers.py index a07b39a..0b79dac 100644 --- a/Wrappers/Python/ccpi/filters/regularisers.py +++ b/Wrappers/Python/ccpi/filters/regularisers.py @@ -2,8 +2,9 @@ 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, NDF_INPAINT_CPU, NVM_INPAINT_CPU -from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_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 +from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU, Diff4th_GPU +from ccpi.filters.cpu_regularisers import NDF_INPAINT_CPU, NVM_INPAINT_CPU def ROF_TV(inputData, regularisation_parameter, iterations, time_marching_parameter,device='cpu'): @@ -110,6 +111,23 @@ def NDF(inputData, regularisation_parameter, edge_parameter, iterations, else: raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ .format(device)) +def DIFF4th(inputData, regularisation_parameter, edge_parameter, iterations, + time_marching_parameter, device='cpu'): + if device == 'cpu': + return Diff4th_CPU(inputData, + regularisation_parameter, + edge_parameter, + iterations, + time_marching_parameter) + elif device == 'gpu': + return Diff4th_GPU(inputData, + regularisation_parameter, + edge_parameter, + 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 986e3e9..51e7fb5 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, NDF +from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, FGP_dTV, TNV, NDF, DIFF4th from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -221,9 +221,9 @@ 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 ,\ + 'regularisation_parameter':0.025, \ + 'edge_parameter':0.015,\ + 'number_of_iterations' :500 ,\ 'time_marching_parameter':0.025,\ 'penalty_type':1 } @@ -255,11 +255,56 @@ plt.title('{}'.format('CPU results')) print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("_____________FGP-dTV (2D)__________________") +print ("___Anisotropic Diffusion 4th Order (2D)____") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot fig = plt.figure(5) +plt.suptitle('Performance of DIFF4th 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' : DIFF4th, \ + 'input' : u0,\ + 'regularisation_parameter':3.5, \ + 'edge_parameter':0.02,\ + 'number_of_iterations' :500 ,\ + 'time_marching_parameter':0.005 + } + +print ("#############DIFF4th CPU################") +start_time = timeit.default_timer() +diff4_cpu = DIFF4th(pars['input'], + pars['regularisation_parameter'], + pars['edge_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'cpu') + +rms = rmse(Im, diff4_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(diff4_cpu, cmap="gray") +plt.title('{}'.format('CPU results')) + + +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("_____________FGP-dTV (2D)__________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure(6) plt.suptitle('Performance of FGP-dTV regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -311,7 +356,7 @@ print ("__________Total nuclear Variation__________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(6) +fig = plt.figure(7) plt.suptitle('Performance of TNV regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -355,257 +400,3 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) imgplot = plt.imshow(tnv_cpu[3,:,:], cmap="gray") plt.title('{}'.format('CPU results')) - - -# Uncomment to test 3D regularisation performance -#%% -""" -slices = 20 -perc = 0.05 - -noisyVol = np.zeros((slices,N,M),dtype='float32') -noisyRef = np.zeros((slices,N,M),dtype='float32') -idealVol = np.zeros((slices,N,M),dtype='float32') - -for i in range (slices): - noisyVol[i,:,:] = Im + np.random.normal(loc = 0 , scale = perc * Im , size = np.shape(Im)) - 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(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') -imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray") - -# set parameters -pars = {'algorithm': ROF_TV, \ - 'input' : noisyVol,\ - 'regularisation_parameter':0.04,\ - 'number_of_iterations': 500,\ - 'time_marching_parameter': 0.0025 - } -print ("#############ROF TV CPU####################") -start_time = timeit.default_timer() -rof_cpu3D = ROF_TV(pars['input'], - pars['regularisation_parameter'], - pars['number_of_iterations'], - pars['time_marching_parameter'],'cpu') -rms = rmse(idealVol, rof_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(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(8) -plt.suptitle('Performance of FGP-TV 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' : FGP_TV, \ - 'input' : noisyVol,\ - 'regularisation_parameter':0.04, \ - 'number_of_iterations' :300 ,\ - 'tolerance_constant':0.00001,\ - 'methodTV': 0 ,\ - 'nonneg': 0 ,\ - 'printingOut': 0 - } - -print ("#############FGP TV CPU####################") -start_time = timeit.default_timer() -fgp_cpu3D = FGP_TV(pars['input'], - pars['regularisation_parameter'], - pars['number_of_iterations'], - pars['tolerance_constant'], - pars['methodTV'], - pars['nonneg'], - pars['printingOut'],'cpu') - - -rms = rmse(idealVol, fgp_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(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(9) -plt.suptitle('Performance of SB-TV 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' : SB_TV, \ - 'input' : noisyVol,\ - 'regularisation_parameter':0.04, \ - 'number_of_iterations' :150 ,\ - 'tolerance_constant':0.00001,\ - 'methodTV': 0 ,\ - 'printingOut': 0 - } - -print ("#############SB TV CPU####################") -start_time = timeit.default_timer() -sb_cpu3D = SB_TV(pars['input'], - pars['regularisation_parameter'], - pars['number_of_iterations'], - pars['tolerance_constant'], - pars['methodTV'], - pars['printingOut'],'cpu') - -rms = rmse(idealVol, sb_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(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(11) -plt.suptitle('Performance of FGP-dTV 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' : FGP_dTV,\ - 'input' : noisyVol,\ - 'refdata' : noisyRef,\ - 'regularisation_parameter':0.04, \ - 'number_of_iterations' :300 ,\ - 'tolerance_constant':0.00001,\ - 'eta_const':0.2,\ - 'methodTV': 0 ,\ - 'nonneg': 0 ,\ - 'printingOut': 0 - } - -print ("#############FGP dTV CPU####################") -start_time = timeit.default_timer() -fgp_dTV_cpu3D = FGP_dTV(pars['input'], - pars['refdata'], - pars['regularisation_parameter'], - pars['number_of_iterations'], - pars['tolerance_constant'], - pars['eta_const'], - pars['methodTV'], - pars['nonneg'], - pars['printingOut'],'cpu') - - -rms = rmse(idealVol, fgp_dTV_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(fgp_dTV_cpu3D[10,:,:], cmap="gray") -plt.title('{}'.format('Recovered volume on the CPU using FGP-dTV')) -""" -#%% diff --git a/Wrappers/Python/demos/demo_cpu_regularisers3D.py b/Wrappers/Python/demos/demo_cpu_regularisers3D.py new file mode 100644 index 0000000..0f47ea9 --- /dev/null +++ b/Wrappers/Python/demos/demo_cpu_regularisers3D.py @@ -0,0 +1,366 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Thu Feb 22 11:39:43 2018 + +Demonstration of 3D CPU regularisers + +@authors: Daniil Kazantsev, Edoardo Pasca +""" + +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, NDF, DIFF4th +from qualitymetrics import rmse +############################################################################### +def printParametersToString(pars): + txt = r'' + for key, value in pars.items(): + if key== 'algorithm' : + txt += "{0} = {1}".format(key, value.__name__) + elif key == 'input': + txt += "{0} = {1}".format(key, np.shape(value)) + elif key == 'refdata': + txt += "{0} = {1}".format(key, np.shape(value)) + else: + txt += "{0} = {1}".format(key, value) + txt += '\n' + return txt +############################################################################### +#%% +filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif") + +# read image +Im = plt.imread(filename) +Im = np.asarray(Im, dtype='float32') + +Im = Im/255 +perc = 0.05 +u0 = Im + np.random.normal(loc = 0 , + scale = perc * Im , + size = np.shape(Im)) +u_ref = Im + np.random.normal(loc = 0 , + scale = 0.01 * Im , + size = np.shape(Im)) +(N,M) = np.shape(u0) +# map the u0 u0->u0>0 +# f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1) +u0 = u0.astype('float32') +u_ref = u_ref.astype('float32') + +# change dims to check that modules work with non-squared images +""" +M = M-100 +u_ref2 = np.zeros([N,M],dtype='float32') +u_ref2[:,0:M] = u_ref[:,0:M] +u_ref = u_ref2 +del u_ref2 + +u02 = np.zeros([N,M],dtype='float32') +u02[:,0:M] = u0[:,0:M] +u0 = u02 +del u02 + +Im2 = np.zeros([N,M],dtype='float32') +Im2[:,0:M] = Im[:,0:M] +Im = Im2 +del Im2 +""" + +# Uncomment to test 3D regularisation performance +#%% +slices = 20 + +noisyVol = np.zeros((slices,N,M),dtype='float32') +noisyRef = np.zeros((slices,N,M),dtype='float32') +idealVol = np.zeros((slices,N,M),dtype='float32') + +for i in range (slices): + noisyVol[i,:,:] = Im + np.random.normal(loc = 0 , scale = perc * Im , size = np.shape(Im)) + 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) +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') +imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray") + +# set parameters +pars = {'algorithm': ROF_TV, \ + 'input' : noisyVol,\ + 'regularisation_parameter':0.04,\ + 'number_of_iterations': 500,\ + 'time_marching_parameter': 0.0025 + } +print ("#############ROF TV CPU####################") +start_time = timeit.default_timer() +rof_cpu3D = ROF_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'cpu') +rms = rmse(idealVol, rof_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(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) +plt.suptitle('Performance of FGP-TV 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' : FGP_TV, \ + 'input' : noisyVol,\ + 'regularisation_parameter':0.04, \ + 'number_of_iterations' :300 ,\ + 'tolerance_constant':0.00001,\ + 'methodTV': 0 ,\ + 'nonneg': 0 ,\ + 'printingOut': 0 + } + +print ("#############FGP TV CPU####################") +start_time = timeit.default_timer() +fgp_cpu3D = FGP_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['nonneg'], + pars['printingOut'],'cpu') + + +rms = rmse(idealVol, fgp_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(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) +plt.suptitle('Performance of SB-TV 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' : SB_TV, \ + 'input' : noisyVol,\ + 'regularisation_parameter':0.04, \ + 'number_of_iterations' :150 ,\ + 'tolerance_constant':0.00001,\ + 'methodTV': 0 ,\ + 'printingOut': 0 + } + +print ("#############SB TV CPU####################") +start_time = timeit.default_timer() +sb_cpu3D = SB_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['printingOut'],'cpu') + +rms = rmse(idealVol, sb_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(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(4) +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.025, \ + 'edge_parameter':0.015,\ + 'number_of_iterations' :500 ,\ + '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 ("___Anisotropic Diffusion 4th Order (2D)____") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure(5) +plt.suptitle('Performance of Diff4th 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' : DIFF4th, \ + 'input' : noisyVol,\ + 'regularisation_parameter':3.5, \ + 'edge_parameter':0.02,\ + 'number_of_iterations' :300 ,\ + 'time_marching_parameter':0.005 + } + +print ("#############DIFF4th CPU################") +start_time = timeit.default_timer() +diff4th_cpu3D = DIFF4th(pars['input'], + pars['regularisation_parameter'], + pars['edge_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter']) + +rms = rmse(idealVol, diff4th_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(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) +plt.suptitle('Performance of FGP-dTV 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' : FGP_dTV,\ + 'input' : noisyVol,\ + 'refdata' : noisyRef,\ + 'regularisation_parameter':0.04, \ + 'number_of_iterations' :300 ,\ + 'tolerance_constant':0.00001,\ + 'eta_const':0.2,\ + 'methodTV': 0 ,\ + 'nonneg': 0 ,\ + 'printingOut': 0 + } + +print ("#############FGP dTV CPU####################") +start_time = timeit.default_timer() +fgp_dTV_cpu3D = FGP_dTV(pars['input'], + pars['refdata'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['eta_const'], + pars['methodTV'], + pars['nonneg'], + pars['printingOut'],'cpu') + + +rms = rmse(idealVol, fgp_dTV_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(fgp_dTV_cpu3D[10,:,:], cmap="gray") +plt.title('{}'.format('Recovered volume on the CPU using FGP-dTV')) +#%% diff --git a/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py index 05db23e..2910c65 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, NDF +from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, FGP_dTV, NDF, DIFF4th from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -207,7 +207,7 @@ plt.title('{}'.format('GPU results')) print ("--------Compare the results--------") tolerance = 1e-05 -diff_im = np.zeros(np.shape(rof_cpu)) +diff_im = np.zeros(np.shape(fgp_cpu)) diff_im = abs(fgp_cpu - fgp_gpu) diff_im[diff_im > tolerance] = 1 a=fig.add_subplot(1,4,4) @@ -293,7 +293,7 @@ plt.title('{}'.format('GPU results')) print ("--------Compare the results--------") tolerance = 1e-05 -diff_im = np.zeros(np.shape(rof_cpu)) +diff_im = np.zeros(np.shape(sb_cpu)) diff_im = abs(sb_cpu - sb_gpu) diff_im[diff_im > tolerance] = 1 a=fig.add_subplot(1,4,4) @@ -379,7 +379,7 @@ plt.title('{}'.format('GPU results')) print ("--------Compare the results--------") tolerance = 1e-05 -diff_im = np.zeros(np.shape(rof_cpu)) +diff_im = np.zeros(np.shape(ndf_cpu)) diff_im = abs(ndf_cpu - ndf_gpu) diff_im[diff_im > tolerance] = 1 a=fig.add_subplot(1,4,4) @@ -391,13 +391,93 @@ else: print ("Arrays match") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("___Anisotropic Diffusion 4th Order (2D)____") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure(5) +plt.suptitle('Comparison of Diff4th 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' : DIFF4th, \ + 'input' : u0,\ + 'regularisation_parameter':3.5, \ + 'edge_parameter':0.02,\ + 'number_of_iterations' :500 ,\ + 'time_marching_parameter':0.005 + } + +print ("#############Diff4th CPU####################") +start_time = timeit.default_timer() +diff4th_cpu = Diff4th(pars['input'], + pars['regularisation_parameter'], + pars['edge_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'cpu') + +rms = rmse(Im, diff4th_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(diff4th_cpu, cmap="gray") +plt.title('{}'.format('CPU results')) + +print ("##############Diff4th GPU##################") +start_time = timeit.default_timer() +diff4th_gpu = Diff4th(pars['input'], + pars['regularisation_parameter'], + pars['edge_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'], 'gpu') + +rms = rmse(Im, diff4th_gpu) +pars['rmse'] = rms +pars['algorithm'] = Diff4th +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(diff4th_gpu, cmap="gray") +plt.title('{}'.format('GPU results')) + +print ("--------Compare the results--------") +tolerance = 1e-05 +diff_im = np.zeros(np.shape(diff4th_cpu)) +diff_im = abs(diff4th_cpu - diff4th_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) +fig = plt.figure(6) plt.suptitle('Comparison of FGP-dTV regulariser using CPU and GPU implementations') a=fig.add_subplot(1,4,1) a.set_title('Noisy Image') @@ -475,7 +555,7 @@ plt.title('{}'.format('GPU results')) print ("--------Compare the results--------") tolerance = 1e-05 -diff_im = np.zeros(np.shape(rof_cpu)) +diff_im = np.zeros(np.shape(fgp_dtv_cpu)) diff_im = abs(fgp_dtv_cpu - fgp_dtv_gpu) diff_im[diff_im > tolerance] = 1 a=fig.add_subplot(1,4,4) diff --git a/Wrappers/Python/demos/demo_gpu_regularisers.py b/Wrappers/Python/demos/demo_gpu_regularisers.py index f3ed50c..8432696 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, NDF +from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, FGP_dTV, NDF, DIFF4th from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -219,9 +219,9 @@ 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 ,\ + 'regularisation_parameter':0.025, \ + 'edge_parameter':0.015,\ + 'number_of_iterations' :500 ,\ 'time_marching_parameter':0.025,\ 'penalty_type': 1 } @@ -253,246 +253,34 @@ plt.title('{}'.format('GPU results')) print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("____________FGP-dTV bench___________________") +print ("___Anisotropic Diffusion 4th Order (2D)____") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot fig = plt.figure(5) -plt.suptitle('Performance of the FGP-dTV regulariser using the GPU') +plt.suptitle('Performance of DIFF4th 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' : FGP_dTV, \ +pars = {'algorithm' : DIFF4th, \ 'input' : u0,\ - 'refdata' : u_ref,\ - 'regularisation_parameter':0.04, \ - 'number_of_iterations' :2000 ,\ - 'tolerance_constant':1e-06,\ - 'eta_const':0.2,\ - 'methodTV': 0 ,\ - 'nonneg': 0 ,\ - 'printingOut': 0 - } - -print ("##############FGP dTV GPU##################") -start_time = timeit.default_timer() -fgp_dtv_gpu = FGP_dTV(pars['input'], - pars['refdata'], - pars['regularisation_parameter'], - pars['number_of_iterations'], - pars['tolerance_constant'], - pars['eta_const'], - pars['methodTV'], - pars['nonneg'], - pars['printingOut'],'gpu') - -rms = rmse(Im, fgp_dtv_gpu) -pars['rmse'] = rms -pars['algorithm'] = FGP_dTV -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(fgp_dtv_gpu, cmap="gray") -plt.title('{}'.format('GPU results')) - - -# Uncomment to test 3D regularisation performance -#%% -""" -N = 512 -slices = 20 - -filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif") -Im = plt.imread(filename) -Im = np.asarray(Im, dtype='float32') - -Im = Im/255 -perc = 0.05 - -noisyVol = np.zeros((slices,N,N),dtype='float32') -noisyRef = np.zeros((slices,N,N),dtype='float32') -idealVol = np.zeros((slices,N,N),dtype='float32') - -for i in range (slices): - noisyVol[i,:,:] = Im + np.random.normal(loc = 0 , scale = perc * Im , size = np.shape(Im)) - 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(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') -imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray") - -# set parameters -pars = {'algorithm': ROF_TV, \ - 'input' : noisyVol,\ - 'regularisation_parameter':0.04,\ - 'number_of_iterations': 500,\ - 'time_marching_parameter': 0.0025 + 'regularisation_parameter':3.5, \ + 'edge_parameter':0.02,\ + 'number_of_iterations' :500 ,\ + 'time_marching_parameter':0.005 } -print ("#############ROF TV GPU####################") + +print ("#############DIFF4th CPU################") start_time = timeit.default_timer() -rof_gpu3D = ROF_TV(pars['input'], - pars['regularisation_parameter'], - pars['number_of_iterations'], - pars['time_marching_parameter'],'gpu') -rms = rmse(idealVol, rof_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(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(7) -plt.suptitle('Performance of FGP-TV 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' : FGP_TV, \ - 'input' : noisyVol,\ - 'regularisation_parameter':0.04, \ - 'number_of_iterations' :300 ,\ - 'tolerance_constant':0.00001,\ - 'methodTV': 0 ,\ - 'nonneg': 0 ,\ - 'printingOut': 0 - } - -print ("#############FGP TV GPU####################") -start_time = timeit.default_timer() -fgp_gpu3D = FGP_TV(pars['input'], - pars['regularisation_parameter'], - pars['number_of_iterations'], - pars['tolerance_constant'], - pars['methodTV'], - pars['nonneg'], - pars['printingOut'],'gpu') - -rms = rmse(idealVol, fgp_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(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(8) -plt.suptitle('Performance of SB-TV 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' : SB_TV, \ - 'input' : noisyVol,\ - 'regularisation_parameter':0.04, \ - 'number_of_iterations' :100 ,\ - 'tolerance_constant':1e-05,\ - 'methodTV': 0 ,\ - 'printingOut': 0 - } - -print ("#############SB TV GPU####################") -start_time = timeit.default_timer() -sb_gpu3D = SB_TV(pars['input'], - pars['regularisation_parameter'], - pars['number_of_iterations'], - pars['tolerance_constant'], - pars['methodTV'], - pars['printingOut'],'gpu') - -rms = rmse(idealVol, sb_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(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'], +diff4_gpu = DIFF4th(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['time_marching_parameter'],'gpu') + +rms = rmse(Im, diff4_gpu) pars['rmse'] = rms txtstr = printParametersToString(pars) @@ -505,49 +293,48 @@ 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')) - +imgplot = plt.imshow(diff4_gpu, cmap="gray") +plt.title('{}'.format('GPU results')) print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("_______________FGP-dTV (3D)________________") +print ("____________FGP-dTV bench___________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(10) -plt.suptitle('Performance of FGP-dTV regulariser using the GPU') +fig = plt.figure(6) +plt.suptitle('Performance of the FGP-dTV regulariser using the GPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') -imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray") +imgplot = plt.imshow(u0,cmap="gray") # set parameters pars = {'algorithm' : FGP_dTV, \ - 'input' : noisyVol,\ - 'refdata' : noisyRef,\ + 'input' : u0,\ + 'refdata' : u_ref,\ 'regularisation_parameter':0.04, \ - 'number_of_iterations' :300 ,\ - 'tolerance_constant':0.00001,\ + 'number_of_iterations' :2000 ,\ + 'tolerance_constant':1e-06,\ 'eta_const':0.2,\ 'methodTV': 0 ,\ 'nonneg': 0 ,\ 'printingOut': 0 } -print ("#############FGP TV GPU####################") +print ("##############FGP dTV GPU##################") start_time = timeit.default_timer() -fgp_dTV_gpu3D = FGP_dTV(pars['input'], +fgp_dtv_gpu = FGP_dTV(pars['input'], pars['refdata'], pars['regularisation_parameter'], pars['number_of_iterations'], pars['tolerance_constant'], - pars['eta_const'], + pars['eta_const'], pars['methodTV'], pars['nonneg'], pars['printingOut'],'gpu') - -rms = rmse(idealVol, fgp_dTV_gpu3D) + +rms = rmse(Im, fgp_dtv_gpu) pars['rmse'] = rms - +pars['algorithm'] = FGP_dTV txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) @@ -558,7 +345,5 @@ 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(fgp_dTV_gpu3D[10,:,:], cmap="gray") -plt.title('{}'.format('Recovered volume on the GPU using FGP-dTV')) -""" -#%% +imgplot = plt.imshow(fgp_dtv_gpu, cmap="gray") +plt.title('{}'.format('GPU results')) diff --git a/Wrappers/Python/demos/demo_gpu_regularisers3D.py b/Wrappers/Python/demos/demo_gpu_regularisers3D.py new file mode 100644 index 0000000..022df95 --- /dev/null +++ b/Wrappers/Python/demos/demo_gpu_regularisers3D.py @@ -0,0 +1,367 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Thu Feb 22 11:39:43 2018 + +Demonstration of GPU regularisers + +@authors: Daniil Kazantsev, Edoardo Pasca +""" + +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 qualitymetrics import rmse +############################################################################### +def printParametersToString(pars): + txt = r'' + for key, value in pars.items(): + if key== 'algorithm' : + txt += "{0} = {1}".format(key, value.__name__) + elif key == 'input': + txt += "{0} = {1}".format(key, np.shape(value)) + elif key == 'refdata': + txt += "{0} = {1}".format(key, np.shape(value)) + else: + txt += "{0} = {1}".format(key, value) + txt += '\n' + return txt +############################################################################### +#%% +filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif") + +# read image +Im = plt.imread(filename) +Im = np.asarray(Im, dtype='float32') + +Im = Im/255 +perc = 0.05 +u0 = Im + np.random.normal(loc = 0 , + scale = perc * Im , + size = np.shape(Im)) +u_ref = Im + np.random.normal(loc = 0 , + scale = 0.01 * Im , + size = np.shape(Im)) +(N,M) = np.shape(u0) +# map the u0 u0->u0>0 +# f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1) +u0 = u0.astype('float32') +u_ref = u_ref.astype('float32') +""" +M = M-100 +u_ref2 = np.zeros([N,M],dtype='float32') +u_ref2[:,0:M] = u_ref[:,0:M] +u_ref = u_ref2 +del u_ref2 + +u02 = np.zeros([N,M],dtype='float32') +u02[:,0:M] = u0[:,0:M] +u0 = u02 +del u02 + +Im2 = np.zeros([N,M],dtype='float32') +Im2[:,0:M] = Im[:,0:M] +Im = Im2 +del Im2 +""" + +#%% +slices = 20 + +filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif") +Im = plt.imread(filename) +Im = np.asarray(Im, dtype='float32') + +Im = Im/255 +perc = 0.05 + +noisyVol = np.zeros((slices,N,N),dtype='float32') +noisyRef = np.zeros((slices,N,N),dtype='float32') +idealVol = np.zeros((slices,N,N),dtype='float32') + +for i in range (slices): + noisyVol[i,:,:] = Im + np.random.normal(loc = 0 , scale = perc * Im , size = np.shape(Im)) + 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) +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') +imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray") + +# set parameters +pars = {'algorithm': ROF_TV, \ + 'input' : noisyVol,\ + 'regularisation_parameter':0.04,\ + 'number_of_iterations': 500,\ + 'time_marching_parameter': 0.0025 + } +print ("#############ROF TV GPU####################") +start_time = timeit.default_timer() +rof_gpu3D = ROF_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'gpu') +rms = rmse(idealVol, rof_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(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) +plt.suptitle('Performance of FGP-TV 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' : FGP_TV, \ + 'input' : noisyVol,\ + 'regularisation_parameter':0.04, \ + 'number_of_iterations' :300 ,\ + 'tolerance_constant':0.00001,\ + 'methodTV': 0 ,\ + 'nonneg': 0 ,\ + 'printingOut': 0 + } + +print ("#############FGP TV GPU####################") +start_time = timeit.default_timer() +fgp_gpu3D = FGP_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['nonneg'], + pars['printingOut'],'gpu') + +rms = rmse(idealVol, fgp_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(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) +plt.suptitle('Performance of SB-TV 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' : SB_TV, \ + 'input' : noisyVol,\ + 'regularisation_parameter':0.04, \ + 'number_of_iterations' :100 ,\ + 'tolerance_constant':1e-05,\ + 'methodTV': 0 ,\ + 'printingOut': 0 + } + +print ("#############SB TV GPU####################") +start_time = timeit.default_timer() +sb_gpu3D = SB_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['printingOut'],'gpu') + +rms = rmse(idealVol, sb_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(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(4) +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.025, \ + 'edge_parameter':0.015,\ + 'number_of_iterations' :500 ,\ + '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 ("___Anisotropic Diffusion 4th Order (3D)____") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure(5) +plt.suptitle('Performance of DIFF4th 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' : DIFF4th, \ + 'input' : noisyVol,\ + 'regularisation_parameter':3.5, \ + 'edge_parameter':0.02,\ + 'number_of_iterations' :300 ,\ + 'time_marching_parameter':0.005 + } + +print ("#############DIFF4th CPU################") +start_time = timeit.default_timer() +diff4_gpu3D = DIFF4th(pars['input'], + pars['regularisation_parameter'], + pars['edge_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'gpu') + +rms = rmse(idealVol, diff4_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(diff4_gpu3D[10,:,:], cmap="gray") +plt.title('{}'.format('GPU results')) + + +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("_______________FGP-dTV (3D)________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure(6) +plt.suptitle('Performance of FGP-dTV 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' : FGP_dTV, \ + 'input' : noisyVol,\ + 'refdata' : noisyRef,\ + 'regularisation_parameter':0.04, \ + 'number_of_iterations' :300 ,\ + 'tolerance_constant':0.00001,\ + 'eta_const':0.2,\ + 'methodTV': 0 ,\ + 'nonneg': 0 ,\ + 'printingOut': 0 + } + +print ("#############FGP TV GPU####################") +start_time = timeit.default_timer() +fgp_dTV_gpu3D = FGP_dTV(pars['input'], + pars['refdata'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['eta_const'], + pars['methodTV'], + pars['nonneg'], + pars['printingOut'],'gpu') + +rms = rmse(idealVol, fgp_dTV_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(fgp_dTV_gpu3D[10,:,:], cmap="gray") +plt.title('{}'.format('Recovered volume on the GPU using FGP-dTV')) +#%% diff --git a/Wrappers/Python/setup-regularisers.py.in b/Wrappers/Python/setup-regularisers.py.in index f55c6fe..76dfecf 100644 --- a/Wrappers/Python/setup-regularisers.py.in +++ b/Wrappers/Python/setup-regularisers.py.in @@ -40,6 +40,7 @@ extra_include_dirs += [os.path.join(".." , ".." , "Core"), os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_SB" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "NDF" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "dTV_FGP" ) , + os.path.join(".." , ".." , "Core", "regularisers_GPU" , "DIFF4th" ) , "."] if platform.system() == 'Windows': diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx index c934f1d..7dc3396 100644 --- a/Wrappers/Python/src/cpu_regularisers.pyx +++ b/Wrappers/Python/src/cpu_regularisers.pyx @@ -22,6 +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 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); 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); @@ -322,6 +323,48 @@ def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, return outputData +#****************************************************************# +#*************Anisotropic Fourth-Order diffusion*****************# +#****************************************************************# +def Diff4th_CPU(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter): + if inputData.ndim == 2: + return Diff4th_2D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type) + elif inputData.ndim == 3: + return Diff4th_3D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type) + +def Diff4th_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularisation_parameter, + float edge_parameter, + int iterationsNumb, + 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 Anisotropic Fourth-Order diffusion for 2D data + Diffus4th_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[1], dims[0], 1) + return outputData + +def Diff4th_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularisation_parameter, + float edge_parameter, + int iterationsNumb, + 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 Anisotropic Fourth-Order diffusion for 3D data + Diffus4th_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[2], dims[1], dims[0]) + + return outputData #*********************Inpainting WITH****************************# #***************Nonlinear (Isotropic) Diffusion******************# #****************************************************************# diff --git a/Wrappers/Python/src/gpu_regularisers.pyx b/Wrappers/Python/src/gpu_regularisers.pyx index 7eab5d5..b67e62b 100644 --- a/Wrappers/Python/src/gpu_regularisers.pyx +++ b/Wrappers/Python/src/gpu_regularisers.pyx @@ -23,6 +23,7 @@ cdef extern void TV_FGP_GPU_main(float *Input, float *Output, float lambdaPar, i 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); +cdef extern void Diffus4th_GPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int N, int M, int Z); # Total-variation Rudin-Osher-Fatemi (ROF) def TV_ROF_GPU(inputData, @@ -135,7 +136,26 @@ def NDF_GPU(inputData, edge_parameter, iterations, time_marching_parameter, - penalty_type) + penalty_type) +# Anisotropic Fourth-Order diffusion +def Diff4th_GPU(inputData, + regularisation_parameter, + edge_parameter, + iterations, + time_marching_parameter): + if inputData.ndim == 2: + return Diff4th_2D(inputData, + regularisation_parameter, + edge_parameter, + iterations, + time_marching_parameter) + elif inputData.ndim == 3: + return Diff4th_3D(inputData, + regularisation_parameter, + edge_parameter, + iterations, + time_marching_parameter) + #****************************************************************# #********************** Total-variation ROF *********************# #****************************************************************# @@ -403,3 +423,41 @@ def NDF_GPU_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 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 +#****************************************************************# +#************Anisotropic Fourth-Order diffusion******************# +#****************************************************************# +def Diff4th_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularisation_parameter, + float edge_parameter, + int iterationsNumb, + 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 Anisotropic Fourth-Order diffusion for 2D data + # Running CUDA code here + Diffus4th_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[1], dims[0], 1) + return outputData + +def Diff4th_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularisation_parameter, + float edge_parameter, + int iterationsNumb, + 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 Anisotropic Fourth-Order diffusion for 3D data + # Running CUDA code here + Diffus4th_GPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[2], dims[1], dims[0]) + + return outputData |