diff options
Diffstat (limited to 'Wrappers')
-rw-r--r-- | Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m | 15 | ||||
-rw-r--r-- | Wrappers/Matlab/demos/demoMatlab_denoise.m | 16 | ||||
-rw-r--r-- | Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c | 29 | ||||
-rw-r--r-- | Wrappers/Python/conda-recipe/conda_build_config.yaml | 2 | ||||
-rwxr-xr-x | Wrappers/Python/conda-recipe/run_test.py | 12 | ||||
-rw-r--r-- | Wrappers/Python/demos/demo_cpu_regularisers.py | 14 | ||||
-rw-r--r-- | Wrappers/Python/demos/demo_cpu_regularisers3D.py | 60 | ||||
-rw-r--r-- | Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py | 14 | ||||
-rw-r--r-- | Wrappers/Python/demos/demo_gpu_regularisers.py | 12 | ||||
-rw-r--r-- | Wrappers/Python/demos/demo_gpu_regularisers3D.py | 55 | ||||
-rw-r--r-- | Wrappers/Python/src/cpu_regularisers.pyx | 35 | ||||
-rw-r--r-- | Wrappers/Python/src/gpu_regularisers.pyx | 39 |
12 files changed, 227 insertions, 76 deletions
diff --git a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m index 5cc47b3..0c331a4 100644 --- a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m +++ b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m @@ -2,11 +2,13 @@ clear; close all Path1 = sprintf(['..' filesep 'mex_compile' filesep 'installed'], 1i); Path2 = sprintf(['..' filesep '..' filesep '..' filesep 'data' filesep], 1i); +Path3 = sprintf(['..' filesep 'supp'], 1i); addpath(Path1); addpath(Path2); +addpath(Path3); N = 512; -slices = 15; +slices = 7; vol3D = zeros(N,N,slices, 'single'); Ideal3D = zeros(N,N,slices, 'single'); Im = double(imread('lena_gray_512.tif'))/255; % loading image @@ -131,7 +133,16 @@ figure; imshow(u_diff4(:,:,7), [0 1]); title('Diffusion 4thO denoised volume (CP % fprintf('%s %f \n', 'RMSE error for Anis.Diff of 4th order is:', rmse_diff4); % figure; imshow(u_diff4_g(:,:,7), [0 1]); title('Diffusion 4thO denoised volume (GPU)'); %% - +fprintf('Denoise using the TGV model (CPU) \n'); +lambda_TGV = 0.03; % regularisation parameter +alpha1 = 1.0; % parameter to control the first-order term +alpha0 = 2.0; % parameter to control the second-order term +iter_TGV = 500; % number of Primal-Dual iterations for TGV +tic; u_tgv = TGV(single(vol3D), lambda_TGV, alpha1, alpha0, iter_TGV); toc; +rmseTGV = RMSE(Ideal3D(:),u_tgv(:)); +fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV); +figure; imshow(u_tgv(:,:,3), [0 1]); title('TGV denoised volume (CPU)'); +%% %>>>>>>>>>>>>>> MULTI-CHANNEL priors <<<<<<<<<<<<<<< % fprintf('Denoise a volume using the FGP-dTV model (CPU) \n'); diff --git a/Wrappers/Matlab/demos/demoMatlab_denoise.m b/Wrappers/Matlab/demos/demoMatlab_denoise.m index 3506cca..14d3096 100644 --- a/Wrappers/Matlab/demos/demoMatlab_denoise.m +++ b/Wrappers/Matlab/demos/demoMatlab_denoise.m @@ -60,20 +60,20 @@ figure; imshow(u_sb, [0 1]); title('SB-TV denoised image (CPU)'); % figure; imshow(u_sbG, [0 1]); title('SB-TV denoised image (GPU)'); %% fprintf('Denoise using the TGV model (CPU) \n'); -lambda_TGV = 0.04; % regularisation parameter -alpha1 = 1; % parameter to control the first-order term -alpha0 = 0.7; % parameter to control the second-order term -iter_TGV = 500; % number of Primal-Dual iterations for TGV +lambda_TGV = 0.045; % regularisation parameter +alpha1 = 1.0; % parameter to control the first-order term +alpha0 = 2.0; % parameter to control the second-order term +iter_TGV = 2000; % number of Primal-Dual iterations for TGV tic; u_tgv = TGV(single(u0), lambda_TGV, alpha1, alpha0, iter_TGV); toc; rmseTGV = (RMSE(u_tgv(:),Im(:))); fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV); figure; imshow(u_tgv, [0 1]); title('TGV denoised image (CPU)'); %% % fprintf('Denoise using the TGV model (GPU) \n'); -% lambda_TGV = 0.04; % regularisation parameter -% alpha1 = 1; % parameter to control the first-order term -% alpha0 = 0.7; % parameter to control the second-order term -% iter_TGV = 500; % number of Primal-Dual iterations for TGV +% lambda_TGV = 0.045; % regularisation parameter +% alpha1 = 1.0; % parameter to control the first-order term +% alpha0 = 2.0; % parameter to control the second-order term +% iter_TGV = 2000; % number of Primal-Dual iterations for TGV % tic; u_tgv_gpu = TGV_GPU(single(u0), lambda_TGV, alpha1, alpha0, iter_TGV); toc; % rmseTGV_gpu = (RMSE(u_tgv_gpu(:),Im(:))); % fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV_gpu); diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c index 5459bf5..aa4eed4 100644 --- a/Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c +++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c @@ -21,14 +21,14 @@ limitations under the License. #include "TGV_core.h" /* C-OMP implementation of Primal-Dual denoising method for - * Total Generilized Variation (TGV)-L2 model [1] (2D case only) + * Total Generilized Variation (TGV)-L2 model [1] (2D/3D) * * Input Parameters: - * 1. Noisy image (2D) (required) - * 2. lambda - regularisation parameter (required) - * 3. parameter to control the first-order term (alpha1) (default - 1) - * 4. parameter to control the second-order term (alpha0) (default - 0.5) - * 5. Number of Chambolle-Pock (Primal-Dual) iterations (default is 300) + * 1. Noisy image/volume (2D/3D) + * 2. lambda - regularisation parameter + * 3. parameter to control the first-order term (alpha1) + * 4. parameter to control the second-order term (alpha0) + * 5. Number of Chambolle-Pock (Primal-Dual) iterations * 6. Lipshitz constant (default is 12) * * Output: @@ -44,7 +44,7 @@ void mexFunction( { int number_of_dims, iter; - mwSize dimX, dimY; + mwSize dimX, dimY, dimZ; const mwSize *dim_array; float *Input, *Output=NULL, lambda, alpha0, alpha1, L2; @@ -55,7 +55,7 @@ void mexFunction( /*Handling Matlab input data*/ if ((nrhs < 2) || (nrhs > 6)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D), Regularisation parameter, alpha0, alpha1, iterations number, Lipshitz Constant"); - Input = (float *) mxGetData(prhs[0]); /*noisy image (2D) */ + Input = (float *) mxGetData(prhs[0]); /*noisy image/volume */ lambda = (float) mxGetScalar(prhs[1]); /* regularisation parameter */ alpha1 = 1.0f; /* parameter to control the first-order term */ alpha0 = 0.5f; /* parameter to control the second-order term */ @@ -69,12 +69,15 @@ void mexFunction( if (nrhs == 6) L2 = (float) mxGetScalar(prhs[5]); /* Lipshitz constant */ /*Handling Matlab output data*/ - dimX = dim_array[0]; dimY = dim_array[1]; + dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; if (number_of_dims == 2) { - Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - /* running the function */ - TGV_main(Input, Output, lambda, alpha1, alpha0, iter, L2, dimX, dimY); + dimZ = 1; /*2D case*/ + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); } - if (number_of_dims == 3) {mexErrMsgTxt("Only 2D images accepted");} + if (number_of_dims == 3) { + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + } + /* running the function */ + TGV_main(Input, Output, lambda, alpha1, alpha0, iter, L2, dimX, dimY, dimZ); } diff --git a/Wrappers/Python/conda-recipe/conda_build_config.yaml b/Wrappers/Python/conda-recipe/conda_build_config.yaml index b7977f3..d6a3915 100644 --- a/Wrappers/Python/conda-recipe/conda_build_config.yaml +++ b/Wrappers/Python/conda-recipe/conda_build_config.yaml @@ -2,6 +2,8 @@ python: - 2.7 # [not win] - 3.5 - 3.6 + - 3.7 numpy: - 1.12 + - 1.14 - 1.15 diff --git a/Wrappers/Python/conda-recipe/run_test.py b/Wrappers/Python/conda-recipe/run_test.py index cfb3f53..21f3216 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):
@@ -303,7 +303,7 @@ class TestRegularisers(unittest.TestCase): 'input' : u0,\
'regularisation_parameter':0.04, \
'alpha1':1.0,\
- 'alpha0':0.7,\
+ 'alpha0':2.0,\
'number_of_iterations' :250 ,\
'LipshitzConstant' :12 ,\
}
@@ -530,7 +530,7 @@ class TestRegularisers(unittest.TestCase): print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
# set parameters
- pars = {'algorithm' : DIFF4th, \
+ pars = {'algorithm' : Diff4th, \
'input' : u0,\
'regularisation_parameter':3.5, \
'edge_parameter':0.02,\
@@ -540,7 +540,7 @@ class TestRegularisers(unittest.TestCase): print ("#############Diff4th CPU####################")
start_time = timeit.default_timer()
- diff4th_cpu = DIFF4th(pars['input'],
+ diff4th_cpu = Diff4th(pars['input'],
pars['regularisation_parameter'],
pars['edge_parameter'],
pars['number_of_iterations'],
@@ -555,7 +555,7 @@ class TestRegularisers(unittest.TestCase): print ("##############Diff4th GPU##################")
start_time = timeit.default_timer()
try:
- diff4th_gpu = DIFF4th(pars['input'],
+ diff4th_gpu = Diff4th(pars['input'],
pars['regularisation_parameter'],
pars['edge_parameter'],
pars['number_of_iterations'],
@@ -565,7 +565,7 @@ class TestRegularisers(unittest.TestCase): self.skipTest("Results not comparable. GPU computing error.")
rms = rmse(Im, diff4th_gpu)
pars['rmse'] = rms
- pars['algorithm'] = DIFF4th
+ pars['algorithm'] = Diff4th
txtstr = printParametersToString(pars)
txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
print (txtstr)
diff --git a/Wrappers/Python/demos/demo_cpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_regularisers.py index 78e9aff..e6befa9 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, LLT_ROF, 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 ccpi.filters.regularisers import PatchSelect, NLTV from qualitymetrics import rmse ############################################################################### @@ -225,8 +225,8 @@ pars = {'algorithm' : TGV, \ 'input' : u0,\ 'regularisation_parameter':0.04, \ 'alpha1':1.0,\ - 'alpha0':0.7,\ - 'number_of_iterations' :250 ,\ + 'alpha0':2.0,\ + 'number_of_iterations' :1350 ,\ 'LipshitzConstant' :12 ,\ } @@ -358,13 +358,13 @@ print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot fig = plt.figure() -plt.suptitle('Performance of DIFF4th regulariser using the CPU') +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, \ +pars = {'algorithm' : Diff4th, \ 'input' : u0,\ 'regularisation_parameter':3.5, \ 'edge_parameter':0.02,\ @@ -372,9 +372,9 @@ pars = {'algorithm' : DIFF4th, \ 'time_marching_parameter':0.0015 } -print ("#############DIFF4th CPU################") +print ("#############Diff4th CPU################") start_time = timeit.default_timer() -diff4_cpu = DIFF4th(pars['input'], +diff4_cpu = Diff4th(pars['input'], pars['regularisation_parameter'], pars['edge_parameter'], pars['number_of_iterations'], diff --git a/Wrappers/Python/demos/demo_cpu_regularisers3D.py b/Wrappers/Python/demos/demo_cpu_regularisers3D.py index 9c28de1..2d2fc22 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, 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 qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -68,7 +68,7 @@ Im2[:,0:M] = Im[:,0:M] Im = Im2 del Im2 """ -slices = 20 +slices = 15 noisyVol = np.zeros((slices,N,M),dtype='float32') noisyRef = np.zeros((slices,N,M),dtype='float32') @@ -96,7 +96,7 @@ pars = {'algorithm': ROF_TV, \ 'input' : noisyVol,\ 'regularisation_parameter':0.04,\ 'number_of_iterations': 500,\ - 'time_marching_parameter': 0.0025 + 'time_marching_parameter': 0.0025 } print ("#############ROF TV CPU####################") start_time = timeit.default_timer() @@ -264,6 +264,54 @@ plt.title('{}'.format('Recovered volume on the CPU using LLT-ROF')) #%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("_______________TGV (3D)_________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure() +plt.suptitle('Performance of TGV 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' : TGV, \ + 'input' : noisyVol,\ + 'regularisation_parameter':0.04, \ + 'alpha1':1.0,\ + 'alpha0':2.0,\ + 'number_of_iterations' :250 ,\ + 'LipshitzConstant' :12 ,\ + } + +print ("#############TGV CPU####################") +start_time = timeit.default_timer() +tgv_cpu3D = TGV(pars['input'], + pars['regularisation_parameter'], + pars['alpha1'], + pars['alpha0'], + pars['number_of_iterations'], + pars['LipshitzConstant'],'cpu') + + +rms = rmse(idealVol, tgv_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(tgv_cpu3D[10,:,:], cmap="gray") +plt.title('{}'.format('Recovered volume on the CPU using TGV')) + +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("________________NDF (3D)___________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") @@ -322,7 +370,7 @@ a.set_title('Noisy volume') imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray") # set parameters -pars = {'algorithm' : DIFF4th, \ +pars = {'algorithm' : Diff4th, \ 'input' : noisyVol,\ 'regularisation_parameter':3.5, \ 'edge_parameter':0.02,\ @@ -330,9 +378,9 @@ pars = {'algorithm' : DIFF4th, \ 'time_marching_parameter':0.0015 } -print ("#############DIFF4th CPU################") +print ("#############Diff4th CPU################") start_time = timeit.default_timer() -diff4th_cpu3D = DIFF4th(pars['input'], +diff4th_cpu3D = Diff4th(pars['input'], pars['regularisation_parameter'], pars['edge_parameter'], pars['number_of_iterations'], diff --git a/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py index 6529b5c..230a761 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, 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 ccpi.filters.regularisers import PatchSelect from qualitymetrics import rmse ############################################################################### @@ -323,8 +323,8 @@ pars = {'algorithm' : TGV, \ 'input' : u0,\ 'regularisation_parameter':0.04, \ 'alpha1':1.0,\ - 'alpha0':0.7,\ - 'number_of_iterations' :250 ,\ + 'alpha0':2.0,\ + 'number_of_iterations' :400 ,\ 'LipshitzConstant' :12 ,\ } @@ -570,7 +570,7 @@ a.set_title('Noisy Image') imgplot = plt.imshow(u0,cmap="gray") # set parameters -pars = {'algorithm' : DIFF4th, \ +pars = {'algorithm' : Diff4th, \ 'input' : u0,\ 'regularisation_parameter':3.5, \ 'edge_parameter':0.02,\ @@ -580,7 +580,7 @@ pars = {'algorithm' : DIFF4th, \ print ("#############Diff4th CPU####################") start_time = timeit.default_timer() -diff4th_cpu = DIFF4th(pars['input'], +diff4th_cpu = Diff4th(pars['input'], pars['regularisation_parameter'], pars['edge_parameter'], pars['number_of_iterations'], @@ -604,7 +604,7 @@ plt.title('{}'.format('CPU results')) print ("##############Diff4th GPU##################") start_time = timeit.default_timer() -diff4th_gpu = DIFF4th(pars['input'], +diff4th_gpu = Diff4th(pars['input'], pars['regularisation_parameter'], pars['edge_parameter'], pars['number_of_iterations'], @@ -612,7 +612,7 @@ diff4th_gpu = DIFF4th(pars['input'], rms = rmse(Im, diff4th_gpu) pars['rmse'] = rms -pars['algorithm'] = DIFF4th +pars['algorithm'] = Diff4th txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) diff --git a/Wrappers/Python/demos/demo_gpu_regularisers.py b/Wrappers/Python/demos/demo_gpu_regularisers.py index 2ada559..e1c6575 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, 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 ccpi.filters.regularisers import PatchSelect, NLTV from qualitymetrics import rmse ############################################################################### @@ -223,8 +223,8 @@ pars = {'algorithm' : TGV, \ 'input' : u0,\ 'regularisation_parameter':0.04, \ 'alpha1':1.0,\ - 'alpha0':0.7,\ - 'number_of_iterations' :250 ,\ + 'alpha0':2.0,\ + 'number_of_iterations' :1250 ,\ 'LipshitzConstant' :12 ,\ } @@ -355,13 +355,13 @@ print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot fig = plt.figure() -plt.suptitle('Performance of DIFF4th 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' : DIFF4th, \ +pars = {'algorithm' : Diff4th, \ 'input' : u0,\ 'regularisation_parameter':3.5, \ 'edge_parameter':0.02,\ @@ -371,7 +371,7 @@ pars = {'algorithm' : DIFF4th, \ print ("#############DIFF4th CPU################") start_time = timeit.default_timer() -diff4_gpu = DIFF4th(pars['input'], +diff4_gpu = Diff4th(pars['input'], pars['regularisation_parameter'], pars['edge_parameter'], pars['number_of_iterations'], diff --git a/Wrappers/Python/demos/demo_gpu_regularisers3D.py b/Wrappers/Python/demos/demo_gpu_regularisers3D.py index d5f9a39..b6058d2 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, 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 qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -67,7 +67,7 @@ Im = Im2 del Im2 """ -#%% + slices = 20 filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif") @@ -268,6 +268,53 @@ plt.title('{}'.format('Recovered volume on the GPU using LLT-ROF')) #%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("_______________TGV (3D)_________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure() +plt.suptitle('Performance of TGV 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' : TGV, \ + 'input' : noisyVol,\ + 'regularisation_parameter':0.04, \ + 'alpha1':1.0,\ + 'alpha0':2.0,\ + 'number_of_iterations' :600 ,\ + 'LipshitzConstant' :12 ,\ + } + +print ("#############TGV GPU####################") +start_time = timeit.default_timer() +tgv_gpu3D = TGV(pars['input'], + pars['regularisation_parameter'], + pars['alpha1'], + pars['alpha0'], + pars['number_of_iterations'], + pars['LipshitzConstant'],'gpu') + + +rms = rmse(idealVol, tgv_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(tgv_gpu3D[10,:,:], cmap="gray") +plt.title('{}'.format('Recovered volume on the GPU using TGV')) +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________NDF-TV (3D)_________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") @@ -326,7 +373,7 @@ a.set_title('Noisy Image') imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray") # set parameters -pars = {'algorithm' : DIFF4th, \ +pars = {'algorithm' : Diff4th, \ 'input' : noisyVol,\ 'regularisation_parameter':3.5, \ 'edge_parameter':0.02,\ @@ -336,7 +383,7 @@ pars = {'algorithm' : DIFF4th, \ print ("#############DIFF4th CPU################") start_time = timeit.default_timer() -diff4_gpu3D = DIFF4th(pars['input'], +diff4_gpu3D = Diff4th(pars['input'], pars['regularisation_parameter'], pars['edge_parameter'], pars['number_of_iterations'], diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx index 7d57ed1..11a0617 100644 --- a/Wrappers/Python/src/cpu_regularisers.pyx +++ b/Wrappers/Python/src/cpu_regularisers.pyx @@ -22,7 +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 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 TGV_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, 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); @@ -202,12 +202,8 @@ def TGV_CPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, Lip return TGV_2D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst) elif inputData.ndim == 3: - shape = inputData.shape - out = inputData.copy() - for i in range(shape[0]): - out[i,:,:] = TGV_2D(inputData[i,:,:], regularisation_parameter, - alpha1, alpha0, iterations, LipshitzConst) - return out + return TGV_3D(inputData, regularisation_parameter, alpha1, alpha0, + iterations, LipshitzConst) def TGV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameter, @@ -229,7 +225,30 @@ def TGV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, alpha0, iterationsNumb, LipshitzConst, - dims[1],dims[0]) + dims[1],dims[0],1) + return outputData +def TGV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularisation_parameter, + float alpha1, + float alpha0, + int iterationsNumb, + float LipshitzConst): + + 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 TGV iterations for 3D data */ + TGV_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, + alpha1, + alpha0, + iterationsNumb, + LipshitzConst, + dims[2], dims[1], dims[0]) return outputData #***************************************************************# diff --git a/Wrappers/Python/src/gpu_regularisers.pyx b/Wrappers/Python/src/gpu_regularisers.pyx index 47a6149..b52f669 100644 --- a/Wrappers/Python/src/gpu_regularisers.pyx +++ b/Wrappers/Python/src/gpu_regularisers.pyx @@ -23,7 +23,7 @@ CUDAErrorMessage = 'CUDA error' cdef extern int TV_ROF_GPU_main(float* Input, float* Output, float lambdaPar, int iter, float tau, int N, int M, int Z); cdef extern int 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 int 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 int TGV_GPU_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY); +cdef extern int TGV_GPU_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY, int dimZ); cdef extern int LLT_ROF_GPU_main(float *Input, float *Output, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, int N, int M, int Z); cdef extern int 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 int 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); @@ -102,12 +102,7 @@ def TGV_GPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, Lip if inputData.ndim == 2: return TGV2D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst) elif inputData.ndim == 3: - shape = inputData.shape - out = inputData.copy() - for i in range(shape[0]): - out[i,:,:] = TGV2D(inputData[i,:,:], regularisation_parameter, - alpha1, alpha0, iterations, LipshitzConst) - return out + return TGV3D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst) # Directional Total-variation Fast-Gradient-Projection (FGP) def dTV_FGP_GPU(inputData, refdata, @@ -393,7 +388,6 @@ def LLT_ROF_GPU3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, raise ValueError(CUDAErrorMessage); - #***************************************************************# #***************** Total Generalised Variation *****************# #***************************************************************# @@ -417,11 +411,38 @@ def TGV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, alpha0, iterationsNumb, LipshitzConst, - dims[1],dims[0])==0): + dims[1],dims[0], 1)==0): return outputData else: raise ValueError(CUDAErrorMessage); +def TGV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularisation_parameter, + float alpha1, + float alpha0, + int iterationsNumb, + float LipshitzConst): + + 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 + if (TGV_GPU_main( + &inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, + alpha1, + alpha0, + iterationsNumb, + LipshitzConst, + dims[2], dims[1], dims[0])==0): + return outputData; + else: + raise ValueError(CUDAErrorMessage); + #****************************************************************# #**************Directional Total-variation FGP ******************# |