From c9ee9ecc84881595b04f19280c93bcd587171270 Mon Sep 17 00:00:00 2001 From: dkazanc Date: Tue, 4 Dec 2018 16:13:38 +0000 Subject: GPU version, this completes implementation of nltv #68 --- Wrappers/Matlab/demos/demoMatlab_denoise.m | 10 +-- Wrappers/Python/ccpi/filters/regularisers.py | 8 ++- Wrappers/Python/conda-recipe/run_test.py | 2 +- Wrappers/Python/demos/demo_cpu_regularisers.py | 18 ++++-- .../Python/demos/demo_cpu_vs_gpu_regularisers.py | 55 ++++++++++++++++ Wrappers/Python/demos/demo_gpu_regularisers.py | 74 +++++++++++++++++++++- Wrappers/Python/setup-regularisers.py.in | 1 + Wrappers/Python/src/cpu_regularisers.pyx | 6 +- Wrappers/Python/src/gpu_regularisers.pyx | 32 ++++++++++ 9 files changed, 189 insertions(+), 17 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Matlab/demos/demoMatlab_denoise.m b/Wrappers/Matlab/demos/demoMatlab_denoise.m index 54b8bac..3506cca 100644 --- a/Wrappers/Matlab/demos/demoMatlab_denoise.m +++ b/Wrappers/Matlab/demos/demoMatlab_denoise.m @@ -138,17 +138,17 @@ figure; imshow(u_diff4, [0 1]); title('Diffusion 4thO denoised image (CPU)'); fprintf('Weights pre-calculation for Non-local TV (takes time on CPU) \n'); SearchingWindow = 7; PatchWindow = 2; -NeighboursNumber = 15; % the number of neibours to include +NeighboursNumber = 20; % the number of neibours to include h = 0.23; % edge related parameter for NLM -[H_i, H_j, Weights] = PatchSelect(single(u0), SearchingWindow, PatchWindow, NeighboursNumber, h); +tic; [H_i, H_j, Weights] = PatchSelect(single(u0), SearchingWindow, PatchWindow, NeighboursNumber, h); toc; %% fprintf('Denoise using Non-local Total Variation (CPU) \n'); -iter_nltv = 2; % number of nltv iterations -lambda_nltv = 0.085; % regularisation parameter for nltv +iter_nltv = 3; % number of nltv iterations +lambda_nltv = 0.05; % regularisation parameter for nltv tic; u_nltv = Nonlocal_TV(single(u0), H_i, H_j, 0, Weights, lambda_nltv, iter_nltv); toc; rmse_nltv = (RMSE(u_nltv(:),Im(:))); fprintf('%s %f \n', 'RMSE error for Non-local Total Variation is:', rmse_nltv); -figure; imshow(u_nltv, [0 1]); title('Non-local Total Variation denoised image (CPU)'); +figure; imagesc(u_nltv, [0 1]); colormap(gray); daspect([1 1 1]); title('Non-local Total Variation denoised image (CPU)'); %% %>>>>>>>>>>>>>> MULTI-CHANNEL priors <<<<<<<<<<<<<<< % diff --git a/Wrappers/Python/ccpi/filters/regularisers.py b/Wrappers/Python/ccpi/filters/regularisers.py index bf7e23c..0a65590 100644 --- a/Wrappers/Python/ccpi/filters/regularisers.py +++ b/Wrappers/Python/ccpi/filters/regularisers.py @@ -4,7 +4,7 @@ script which assigns a proper device core function based on a flag ('cpu' or 'gp 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, PATCHSEL_CPU, NLTV_CPU try: - 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.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU, Diff4th_GPU, TGV_GPU, LLT_ROF_GPU, PATCHSEL_GPU gpu_enabled = True except ImportError: gpu_enabled = False @@ -153,7 +153,11 @@ def PatchSelect(inputData, searchwindow, patchwindow, neighbours, edge_parameter neighbours, edge_parameter) elif device == 'gpu' and gpu_enabled: - return 1 + return PATCHSEL_GPU(inputData, + searchwindow, + patchwindow, + neighbours, + edge_parameter) else: if not gpu_enabled and device == 'gpu': raise ValueError ('GPU is not available') diff --git a/Wrappers/Python/conda-recipe/run_test.py b/Wrappers/Python/conda-recipe/run_test.py index 6ffaca1..499ae7f 100755 --- a/Wrappers/Python/conda-recipe/run_test.py +++ b/Wrappers/Python/conda-recipe/run_test.py @@ -2,7 +2,7 @@ import unittest 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, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, DIFF4th from PIL import Image class TiffReader(object): diff --git a/Wrappers/Python/demos/demo_cpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_regularisers.py index 31e4cad..78e9aff 100644 --- a/Wrappers/Python/demos/demo_cpu_regularisers.py +++ b/Wrappers/Python/demos/demo_cpu_regularisers.py @@ -400,20 +400,29 @@ plt.title('{}'.format('CPU results')) print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("___Nonlocal patches pre-calculation____") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +start_time = timeit.default_timer() # set parameters pars = {'algorithm' : PatchSelect, \ 'input' : u0,\ 'searchwindow': 7, \ 'patchwindow': 2,\ 'neighbours' : 15 ,\ - 'edge_parameter':0.23} + 'edge_parameter':0.18} H_i, H_j, Weights = PatchSelect(pars['input'], pars['searchwindow'], pars['patchwindow'], pars['neighbours'], pars['edge_parameter'],'cpu') - + +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +""" +plt.figure() +plt.imshow(Weights[0,:,:],cmap="gray",interpolation="nearest",vmin=0, vmax=1) +plt.show() +""" #%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("___Nonlocal Total Variation penalty____") @@ -431,10 +440,9 @@ pars2 = {'algorithm' : NLTV, \ 'H_j': H_j,\ 'H_k' : 0,\ 'Weights' : Weights,\ - 'regularisation_parameter': 0.085,\ - 'iterations': 2 + 'regularisation_parameter': 0.04,\ + 'iterations': 3 } -#%% start_time = timeit.default_timer() nltv_cpu = NLTV(pars2['input'], pars2['H_i'], diff --git a/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py index 3d6e92f..616eab0 100644 --- a/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py +++ b/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py @@ -13,6 +13,7 @@ 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 PatchSelect from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -728,6 +729,60 @@ 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 ("____Non-local regularisation bench_________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure() +plt.suptitle('Comparison of Nonlocal TV regulariser using CPU and GPU implementations') +a=fig.add_subplot(1,2,1) +a.set_title('Noisy Image') +imgplot = plt.imshow(u0,cmap="gray") + +pars = {'algorithm' : PatchSelect, \ + 'input' : u0,\ + 'searchwindow': 7, \ + 'patchwindow': 2,\ + 'neighbours' : 15 ,\ + 'edge_parameter':0.18} + +print ("############## Nonlocal Patches on CPU##################") +start_time = timeit.default_timer() +H_i, H_j, WeightsCPU = PatchSelect(pars['input'], + pars['searchwindow'], + pars['patchwindow'], + pars['neighbours'], + pars['edge_parameter'],'cpu') +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) + +print ("############## Nonlocal Patches on GPU##################") +start_time = timeit.default_timer() +start_time = timeit.default_timer() +H_i, H_j, WeightsGPU = PatchSelect(pars['input'], + pars['searchwindow'], + pars['patchwindow'], + pars['neighbours'], + pars['edge_parameter'],'gpu') +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) + +print ("--------Compare the results--------") +tolerance = 1e-05 +diff_im = np.zeros(np.shape(u0)) +diff_im = abs(WeightsCPU[0,:,:] - WeightsGPU[0,:,:]) +diff_im[diff_im > tolerance] = 1 +a=fig.add_subplot(1,2,2) +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: diff --git a/Wrappers/Python/demos/demo_gpu_regularisers.py b/Wrappers/Python/demos/demo_gpu_regularisers.py index de0cbde..2ada559 100644 --- a/Wrappers/Python/demos/demo_gpu_regularisers.py +++ b/Wrappers/Python/demos/demo_gpu_regularisers.py @@ -13,6 +13,7 @@ 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 PatchSelect, NLTV from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -84,7 +85,7 @@ pars = {'algorithm': ROF_TV, \ 'input' : u0,\ 'regularisation_parameter':0.04,\ 'number_of_iterations': 1200,\ - 'time_marching_parameter': 0.0025 + 'time_marching_parameter': 0.0025 } print ("##############ROF TV GPU##################") start_time = timeit.default_timer() @@ -392,6 +393,77 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(diff4_gpu, cmap="gray") plt.title('{}'.format('GPU results')) +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("___Nonlocal patches pre-calculation____") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +start_time = timeit.default_timer() +# set parameters +pars = {'algorithm' : PatchSelect, \ + 'input' : u0,\ + 'searchwindow': 7, \ + 'patchwindow': 2,\ + 'neighbours' : 15 ,\ + 'edge_parameter':0.18} + +H_i, H_j, Weights = PatchSelect(pars['input'], + pars['searchwindow'], + pars['patchwindow'], + pars['neighbours'], + pars['edge_parameter'],'gpu') + +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +""" +plt.figure() +plt.imshow(Weights[0,:,:],cmap="gray",interpolation="nearest",vmin=0, vmax=1) +plt.show() +""" +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("___Nonlocal Total Variation penalty____") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +## plot +fig = plt.figure() +plt.suptitle('Performance of NLTV regulariser using the CPU') +a=fig.add_subplot(1,2,1) +a.set_title('Noisy Image') +imgplot = plt.imshow(u0,cmap="gray") + +pars2 = {'algorithm' : NLTV, \ + 'input' : u0,\ + 'H_i': H_i, \ + 'H_j': H_j,\ + 'H_k' : 0,\ + 'Weights' : Weights,\ + 'regularisation_parameter': 0.02,\ + 'iterations': 3 + } +start_time = timeit.default_timer() +nltv_cpu = NLTV(pars2['input'], + pars2['H_i'], + pars2['H_j'], + pars2['H_k'], + pars2['Weights'], + pars2['regularisation_parameter'], + pars2['iterations']) + +rms = rmse(Im, nltv_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(nltv_cpu, cmap="gray") +plt.title('{}'.format('CPU results')) #%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("____________FGP-dTV bench___________________") diff --git a/Wrappers/Python/setup-regularisers.py.in b/Wrappers/Python/setup-regularisers.py.in index 542dcb4..462edda 100644 --- a/Wrappers/Python/setup-regularisers.py.in +++ b/Wrappers/Python/setup-regularisers.py.in @@ -45,6 +45,7 @@ extra_include_dirs += [os.path.join(".." , ".." , "Core"), os.path.join(".." , ".." , "Core", "regularisers_GPU" , "NDF" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "dTV_FGP" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "DIFF4th" ) , + os.path.join(".." , ".." , "Core", "regularisers_GPU" , "PatchSelect" ) , "."] if platform.system() == 'Windows': diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx index e51e6d8..4aa3251 100644 --- a/Wrappers/Python/src/cpu_regularisers.pyx +++ b/Wrappers/Python/src/cpu_regularisers.pyx @@ -456,7 +456,7 @@ def PATCHSEL_CPU(inputData, searchwindow, patchwindow, neighbours, edge_paramete if inputData.ndim == 2: return PatchSel_2D(inputData, searchwindow, patchwindow, neighbours, edge_parameter) elif inputData.ndim == 3: - return PatchSel_3D(inputData, searchwindow, patchwindow, neighbours, edge_parameter) + return 1 def PatchSel_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, int searchwindow, int patchwindow, @@ -480,7 +480,7 @@ def PatchSel_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, # Run patch-based weight selection function PatchSelect_CPU_main(&inputData[0,0], &H_j[0,0,0], &H_i[0,0,0], &H_i[0,0,0], &Weights[0,0,0], dims[2], dims[1], 0, searchwindow, patchwindow, neighbours, edge_parameter, 1) return H_i, H_j, Weights - +""" def PatchSel_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, int searchwindow, int patchwindow, @@ -507,7 +507,7 @@ def PatchSel_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, # Run patch-based weight selection function PatchSelect_CPU_main(&inputData[0,0,0], &H_i[0,0,0,0], &H_j[0,0,0,0], &H_k[0,0,0,0], &Weights[0,0,0,0], dims[2], dims[1], dims[0], searchwindow, patchwindow, neighbours, edge_parameter, 1) return H_i, H_j, H_k, Weights - +""" #****************************************************************# #***************Non-local Total Variation******************# diff --git a/Wrappers/Python/src/gpu_regularisers.pyx b/Wrappers/Python/src/gpu_regularisers.pyx index 82d3e01..302727e 100644 --- a/Wrappers/Python/src/gpu_regularisers.pyx +++ b/Wrappers/Python/src/gpu_regularisers.pyx @@ -26,6 +26,7 @@ cdef extern void LLT_ROF_GPU_main(float *Input, float *Output, float lambdaROF, 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); +cdef extern void PatchSelect_GPU_main(float *Input, unsigned short *H_i, unsigned short *H_j, float *Weights, int N, int M, int SearchWindow, int SimilarWin, int NumNeighb, float h); # Total-variation Rudin-Osher-Fatemi (ROF) def TV_ROF_GPU(inputData, @@ -542,3 +543,34 @@ def Diff4th_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 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 +#****************************************************************# +#************Patch-based weights pre-selection******************# +#****************************************************************# +def PATCHSEL_GPU(inputData, searchwindow, patchwindow, neighbours, edge_parameter): + if inputData.ndim == 2: + return PatchSel_2D(inputData, searchwindow, patchwindow, neighbours, edge_parameter) + elif inputData.ndim == 3: + return 1 +def PatchSel_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + int searchwindow, + int patchwindow, + int neighbours, + float edge_parameter): + cdef long dims[3] + dims[0] = neighbours + dims[1] = inputData.shape[0] + dims[2] = inputData.shape[1] + + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Weights = \ + np.zeros([dims[0], dims[1],dims[2]], dtype='float32') + + cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_i = \ + np.zeros([dims[0], dims[1],dims[2]], dtype='uint16') + + cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_j = \ + np.zeros([dims[0], dims[1],dims[2]], dtype='uint16') + + # Run patch-based weight selection function + PatchSelect_GPU_main(&inputData[0,0], &H_j[0,0,0], &H_i[0,0,0], &Weights[0,0,0], dims[2], dims[1], searchwindow, patchwindow, neighbours, edge_parameter) + + return H_i, H_j, Weights -- cgit v1.2.3