From c65291e6b987283e4767a8ad2bd2d2433ca3782e Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Thu, 28 Nov 2019 23:01:03 +0000 Subject: all work completed on gpu version of pdtv --- demos/Matlab_demos/demoMatlab_3Ddenoise.m | 19 ++++++ demos/Matlab_demos/demoMatlab_denoise.m | 16 +++++ demos/demo_cpu_regularisers.py | 4 +- demos/demo_cpu_regularisers3D.py | 10 +-- demos/demo_cpu_vs_gpu_regularisers.py | 104 +++++++++++++++++++++++++++++- demos/demo_gpu_regularisers.py | 51 ++++++++++++++- demos/demo_gpu_regularisers3D.py | 49 +++++++++++++- 7 files changed, 244 insertions(+), 9 deletions(-) (limited to 'demos') diff --git a/demos/Matlab_demos/demoMatlab_3Ddenoise.m b/demos/Matlab_demos/demoMatlab_3Ddenoise.m index f018327..b7f92cb 100644 --- a/demos/Matlab_demos/demoMatlab_3Ddenoise.m +++ b/demos/Matlab_demos/demoMatlab_3Ddenoise.m @@ -62,6 +62,25 @@ fprintf('Denoise a volume using the FGP-TV model (GPU) \n'); % fprintf('%s %f \n', 'RMSE error for FGP-TV is:', rmse_fgpG); % figure; imshow(u_fgpG(:,:,7), [0 1]); title('FGP-TV denoised volume (GPU)'); %% +fprintf('Denoise a volume using the PD-TV model (CPU) \n'); +lambda_reg = 0.03; % regularsation parameter for all methods +iter_pd = 300; % number of FGP iterations +epsil_tol = 0.0; % tolerance +tic; [u_pd,infovec] = PD_TV(single(vol3D), lambda_reg, iter_pd, epsil_tol); toc; +energyfunc_val_fgp = TV_energy(single(u_pd),single(vol3D),lambda_reg, 1); % get energy function value +rmse_pd = (RMSE(Ideal3D(:),u_pd(:))); +fprintf('%s %f \n', 'RMSE error for PD-TV is:', rmse_pd); +figure; imshow(u_pd(:,:,7), [0 1]); title('PD-TV denoised volume (CPU)'); +%% +% fprintf('Denoise a volume using the PD-TV model (GPU) \n'); +% lambda_reg = 0.03; % regularsation parameter for all methods +% iter_pd = 300; % number of FGP iterations +% epsil_tol = 0.0; % tolerance +% tic; u_pdG = PD_TV_GPU(single(vol3D), lambda_reg, iter_pd, epsil_tol); toc; +% rmse_pdG = (RMSE(Ideal3D(:),u_pdG(:))); +% fprintf('%s %f \n', 'RMSE error for PD-TV is:', rmse_pdG); +% figure; imshow(u_pdG(:,:,7), [0 1]); title('PD-TV denoised volume (GPU)'); +%% fprintf('Denoise a volume using the SB-TV model (CPU) \n'); iter_sb = 150; % number of SB iterations epsil_tol = 0.0; % tolerance diff --git a/demos/Matlab_demos/demoMatlab_denoise.m b/demos/Matlab_demos/demoMatlab_denoise.m index b50eaf5..3d93cb6 100644 --- a/demos/Matlab_demos/demoMatlab_denoise.m +++ b/demos/Matlab_demos/demoMatlab_denoise.m @@ -46,6 +46,22 @@ figure; imshow(u_fgp, [0 1]); title('FGP-TV denoised image (CPU)'); % tic; u_fgpG = FGP_TV_GPU(single(u0), lambda_reg, iter_fgp, epsil_tol); toc; % figure; imshow(u_fgpG, [0 1]); title('FGP-TV denoised image (GPU)'); %% +fprintf('Denoise using the PD-TV model (CPU) \n'); +lambda_reg = 0.03; +iter_pd = 500; % number of FGP iterations +epsil_tol = 0.0; % tolerance +tic; [u_pd,infovec] = PD_TV(single(u0), lambda_reg, iter_pd, epsil_tol); toc; +energyfunc_val_pd = TV_energy(single(u_pd),single(u0),lambda_reg, 1); % get energy function value +rmsePD = (RMSE(u_pd(:),Im(:))); +fprintf('%s %f \n', 'RMSE error for PD-TV is:', rmsePD); +[ssimval] = ssim(u_pd*255,single(Im)*255); +fprintf('%s %f \n', 'MSSIM error for PD-TV is:', ssimval); +figure; imshow(u_pd, [0 1]); title('PD-TV denoised image (CPU)'); +%% +% fprintf('Denoise using the PD-TV model (GPU) \n'); +% tic; u_pdG = PD_TV_GPU(single(u0), lambda_reg, iter_pd, epsil_tol); toc; +% figure; imshow(u_pdG, [0 1]); title('PD-TV denoised image (GPU)'); +%% fprintf('Denoise using the SB-TV model (CPU) \n'); lambda_reg = 0.03; iter_sb = 200; % number of SB iterations diff --git a/demos/demo_cpu_regularisers.py b/demos/demo_cpu_regularisers.py index 09781d5..7d66b7f 100644 --- a/demos/demo_cpu_regularisers.py +++ b/demos/demo_cpu_regularisers.py @@ -130,7 +130,7 @@ pars = {'algorithm' : FGP_TV, \ 'input' : u0,\ 'regularisation_parameter':0.02, \ 'number_of_iterations' :1500 ,\ - 'tolerance_constant':1e-06,\ + 'tolerance_constant':1e-08,\ 'methodTV': 0 ,\ 'nonneg': 0} @@ -176,7 +176,7 @@ pars = {'algorithm' : PD_TV, \ 'input' : u0,\ 'regularisation_parameter':0.02, \ 'number_of_iterations' :1500 ,\ - 'tolerance_constant':1e-06,\ + 'tolerance_constant':1e-08,\ 'methodTV': 0 ,\ 'nonneg': 1, 'lipschitz_const' : 8, diff --git a/demos/demo_cpu_regularisers3D.py b/demos/demo_cpu_regularisers3D.py index 349300f..cfdd2d4 100644 --- a/demos/demo_cpu_regularisers3D.py +++ b/demos/demo_cpu_regularisers3D.py @@ -185,10 +185,11 @@ pars = {'algorithm' : PD_TV, \ 'input' : noisyVol,\ 'regularisation_parameter':0.02, \ 'number_of_iterations' :1000 ,\ - 'tolerance_constant': 1e-06,\ + 'tolerance_constant':1e-06,\ 'methodTV': 0 ,\ - 'lipschitz_const' : 12,\ - 'nonneg': 0} + 'nonneg': 0, + 'lipschitz_const' : 8, + 'tau' : 0.0025} print ("#############FGP TV GPU####################") start_time = timeit.default_timer() @@ -198,7 +199,8 @@ start_time = timeit.default_timer() pars['tolerance_constant'], pars['methodTV'], pars['nonneg'], - pars['lipschitz_const'], 'cpu') + pars['lipschitz_const'], + pars['tau'],'cpu') Qtools = QualityTools(idealVol, pd_cpu3D) pars['rmse'] = Qtools.rmse() diff --git a/demos/demo_cpu_vs_gpu_regularisers.py b/demos/demo_cpu_vs_gpu_regularisers.py index 21e3899..015dfc6 100644 --- a/demos/demo_cpu_vs_gpu_regularisers.py +++ b/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, LLT_ROF, FGP_dTV, NDF, Diff4th +from ccpi.filters.regularisers import ROF_TV, FGP_TV, PD_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th from ccpi.filters.regularisers import PatchSelect from ccpi.supp.qualitymetrics import QualityTools ############################################################################### @@ -220,6 +220,108 @@ if (diff_im.sum() > 1): else: print ("Arrays match") + +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("____________PD-TV bench___________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure() +plt.suptitle('Comparison of PD-TV 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' : PD_TV, \ + 'input' : u0,\ + 'regularisation_parameter':0.02, \ + 'number_of_iterations' :1500 ,\ + 'tolerance_constant':0.0,\ + 'methodTV': 0 ,\ + 'nonneg': 0, + 'lipschitz_const' : 8, + 'tau' : 0.0025} + +print ("#############PD TV CPU####################") +start_time = timeit.default_timer() +(pd_cpu,info_vec_cpu) = PD_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['nonneg'], + pars['lipschitz_const'], + pars['tau'],'cpu') + +Qtools = QualityTools(Im, pd_cpu) +pars['rmse'] = Qtools.rmse() + +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(pd_cpu, cmap="gray") +plt.title('{}'.format('CPU results')) + +# set parameters +pars = {'algorithm' : PD_TV, \ + 'input' : u0,\ + 'regularisation_parameter':0.02, \ + 'number_of_iterations' :1500 ,\ + 'tolerance_constant':0.0,\ + 'methodTV': 0 ,\ + 'nonneg': 0, + 'lipschitz_const' : 8, + 'tau' : 0.0025} + +print ("#############PD TV CPU####################") +start_time = timeit.default_timer() +(pd_gpu,info_vec_gpu) = PD_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['nonneg'], + pars['lipschitz_const'], + pars['tau'],'gpu') + +Qtools = QualityTools(Im, pd_gpu) +pars['rmse'] = Qtools.rmse() + +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(pd_gpu, cmap="gray") +plt.title('{}'.format('GPU results')) + + +print ("--------Compare the results--------") +tolerance = 1e-05 +diff_im = np.zeros(np.shape(pd_cpu)) +diff_im = abs(pd_cpu - pd_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 ("____________SB-TV bench___________________") diff --git a/demos/demo_gpu_regularisers.py b/demos/demo_gpu_regularisers.py index 3efcfce..5131847 100644 --- a/demos/demo_gpu_regularisers.py +++ b/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, LLT_ROF, FGP_dTV, NDF, Diff4th +from ccpi.filters.regularisers import ROF_TV, FGP_TV, PD_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th from ccpi.filters.regularisers import PatchSelect, NLTV from ccpi.supp.qualitymetrics import QualityTools ############################################################################### @@ -158,6 +158,55 @@ imgplot = plt.imshow(fgp_gpu, cmap="gray") plt.title('{}'.format('GPU results')) #%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("_______________PD-TV (2D)__________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure() +plt.suptitle('Performance of PD-TV 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' : PD_TV, \ + 'input' : u0,\ + 'regularisation_parameter':0.02, \ + 'number_of_iterations' :1500 ,\ + 'tolerance_constant':1e-06,\ + 'methodTV': 0 ,\ + 'nonneg': 1, + 'lipschitz_const' : 8, + 'tau' : 0.0025} + +print ("#############PD TV CPU####################") +start_time = timeit.default_timer() +(pd_gpu,info_vec_gpu) = PD_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['nonneg'], + pars['lipschitz_const'], + pars['tau'],'gpu') + +Qtools = QualityTools(Im, pd_gpu) +pars['rmse'] = Qtools.rmse() + +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(pd_gpu, cmap="gray") +plt.title('{}'.format('GPU results')) +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("____________SB-TV regulariser______________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") diff --git a/demos/demo_gpu_regularisers3D.py b/demos/demo_gpu_regularisers3D.py index ccf9694..2c25d01 100644 --- a/demos/demo_gpu_regularisers3D.py +++ b/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, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th +from ccpi.filters.regularisers import ROF_TV, FGP_TV, PD_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th from ccpi.supp.qualitymetrics import QualityTools ############################################################################### def printParametersToString(pars): @@ -172,7 +172,54 @@ 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 ("_______________PD-TV (3D)__________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure() +plt.suptitle('Performance of PD-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' : PD_TV, \ + 'input' : noisyVol,\ + 'regularisation_parameter':0.02, \ + 'number_of_iterations' :1000 ,\ + 'tolerance_constant':1e-06,\ + 'methodTV': 0 ,\ + 'nonneg': 0, + 'lipschitz_const' : 8, + 'tau' : 0.0025} + +print ("#############PD TV GPU####################") +start_time = timeit.default_timer() +(pd_gpu3D, info_vec_gpu) = PD_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['nonneg'], + pars['lipschitz_const'], + pars['tau'],'gpu') + +Qtools = QualityTools(idealVol, pd_gpu3D) +pars['rmse'] = Qtools.rmse() +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(pd_gpu3D[10,:,:], cmap="gray") +plt.title('{}'.format('Recovered volume on the GPU using PD-TV')) #%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________SB-TV (3D)__________________") -- cgit v1.2.3