diff options
Diffstat (limited to 'Wrappers')
| -rw-r--r-- | Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m | 16 | ||||
| -rw-r--r-- | Wrappers/Matlab/demos/demoMatlab_denoise.m | 21 | ||||
| -rw-r--r-- | Wrappers/Matlab/mex_compile/compileCPU_mex.m | 5 | ||||
| -rw-r--r-- | Wrappers/Matlab/mex_compile/compileGPU_mex.m | 12 | ||||
| -rw-r--r-- | Wrappers/Matlab/mex_compile/regularisers_CPU/Diffusion_4thO.c | 76 | ||||
| -rw-r--r-- | Wrappers/Matlab/mex_compile/regularisers_GPU/Diffusion_4thO_GPU.cpp | 76 | ||||
| -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 |
15 files changed, 1229 insertions, 524 deletions
diff --git a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m index c087433..ccd47f1 100644 --- a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m +++ b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m @@ -74,6 +74,22 @@ figure; imshow(u_diff(:,:,15), [0 1]); title('Diffusion denoised volume (CPU)'); % tic; u_diff_g = NonlDiff_GPU(single(vol3D), lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber'); toc; % figure; imshow(u_diff_g(:,:,15), [0 1]); title('Diffusion denoised volume (GPU)'); %% +fprintf('Denoise using Fourth-order anisotropic diffusion model (CPU) \n'); +iter_diff = 300; % number of diffusion iterations +lambda_regDiff = 3.5; % regularisation for the diffusivity +sigmaPar = 0.02; % edge-preserving parameter +tau_param = 0.0025; % time-marching constant +tic; u_diff4 = Diffusion_4thO(single(vol3D), lambda_regDiff, sigmaPar, iter_diff, tau_param); toc; +figure; imshow(u_diff4(:,:,15), [0 1]); title('Diffusion 4thO denoised volume (CPU)'); +%% +% fprintf('Denoise using Fourth-order anisotropic diffusion model (GPU) \n'); +% iter_diff = 300; % number of diffusion iterations +% lambda_regDiff = 3.5; % regularisation for the diffusivity +% sigmaPar = 0.02; % edge-preserving parameter +% tau_param = 0.0025; % time-marching constant +% tic; u_diff4_g = Diffusion_4thO_GPU(single(vol3D), lambda_regDiff, sigmaPar, iter_diff, tau_param); toc; +% figure; imshow(u_diff4_g(:,:,15), [0 1]); title('Diffusion 4thO denoised volume (GPU)'); +%% %>>>>>>>>>>>>>> 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 d93f477..30ad79d 100644 --- a/Wrappers/Matlab/demos/demoMatlab_denoise.m +++ b/Wrappers/Matlab/demos/demoMatlab_denoise.m @@ -53,8 +53,8 @@ figure; imshow(u_sb, [0 1]); title('SB-TV denoised image (CPU)'); %% fprintf('Denoise using Nonlinear-Diffusion model (CPU) \n'); iter_diff = 800; % number of diffusion iterations -lambda_regDiff = 0.06; % regularisation for the diffusivity -sigmaPar = 0.04; % edge-preserving parameter +lambda_regDiff = 0.025; % regularisation for the diffusivity +sigmaPar = 0.015; % edge-preserving parameter tau_param = 0.025; % time-marching constant tic; u_diff = NonlDiff(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber'); toc; figure; imshow(u_diff, [0 1]); title('Diffusion denoised image (CPU)'); @@ -67,6 +67,23 @@ figure; imshow(u_diff, [0 1]); title('Diffusion denoised image (CPU)'); % tic; u_diff_g = NonlDiff_GPU(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber'); toc; % figure; imshow(u_diff_g, [0 1]); title('Diffusion denoised image (GPU)'); %% +fprintf('Denoise using Fourth-order anisotropic diffusion model (CPU) \n'); +iter_diff = 800; % number of diffusion iterations +lambda_regDiff = 3.5; % regularisation for the diffusivity +sigmaPar = 0.02; % edge-preserving parameter +tau_param = 0.005; % time-marching constant +tic; u_diff4 = Diffusion_4thO(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param); toc; +figure; imshow(u_diff4, [0 1]); title('Diffusion 4thO denoised image (CPU)'); +%% +% fprintf('Denoise using Fourth-order anisotropic diffusion model (GPU) \n'); +% iter_diff = 800; % number of diffusion iterations +% lambda_regDiff = 3.5; % regularisation for the diffusivity +% sigmaPar = 0.02; % edge-preserving parameter +% tau_param = 0.005; % time-marching constant +% tic; u_diff4_g = Diffusion_4thO_GPU(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param); toc; +% figure; imshow(u_diff4_g, [0 1]); title('Diffusion 4thO denoised image (GPU)'); +%% + %>>>>>>>>>>>>>> MULTI-CHANNEL priors <<<<<<<<<<<<<<< % fprintf('Denoise using the FGP-dTV model (CPU) \n'); diff --git a/Wrappers/Matlab/mex_compile/compileCPU_mex.m b/Wrappers/Matlab/mex_compile/compileCPU_mex.m index b232f33..19eb301 100644 --- a/Wrappers/Matlab/mex_compile/compileCPU_mex.m +++ b/Wrappers/Matlab/mex_compile/compileCPU_mex.m @@ -31,6 +31,9 @@ movefile('TNV.mex*',Pathmove); mex NonlDiff.c Diffusion_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" movefile('NonlDiff.mex*',Pathmove); +mex Diffusion_4thO.c Diffus4th_order_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" +movefile('Diffusion_4thO.mex*',Pathmove); + mex TV_energy.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" movefile('TV_energy.mex*',Pathmove); @@ -41,7 +44,7 @@ movefile('NonlDiff_Inp.mex*',Pathmove); mex NonlocalMarching_Inpaint.c NonlocalMarching_Inpaint_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" movefile('NonlocalMarching_Inpaint.mex*',Pathmove); -delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* CCPiDefines.h +delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* Diffus4th_order_core* CCPiDefines.h delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core* fprintf('%s \n', 'Regularisers successfully compiled!'); diff --git a/Wrappers/Matlab/mex_compile/compileGPU_mex.m b/Wrappers/Matlab/mex_compile/compileGPU_mex.m index 46d85a6..6b69c34 100644 --- a/Wrappers/Matlab/mex_compile/compileGPU_mex.m +++ b/Wrappers/Matlab/mex_compile/compileGPU_mex.m @@ -42,5 +42,13 @@ movefile('FGP_dTV_GPU.mex*',Pathmove); mex -g -I/usr/local/cuda-7.5/include -L/usr/local/cuda-7.5/lib64 -lcudart -lcufft -lmwgpu NonlDiff_GPU.cpp NonlDiff_GPU_core.o movefile('NonlDiff_GPU.mex*',Pathmove); -delete TV_ROF_GPU_core* TV_FGP_GPU_core* TV_SB_GPU_core* dTV_FGP_GPU_core* NonlDiff_GPU_core* CCPiDefines.h -fprintf('%s \n', 'All successfully compiled!');
\ No newline at end of file +!/usr/local/cuda/bin/nvcc -O0 -c Diffus_4thO_GPU_core.cu -Xcompiler -fPIC -I~/SOFT/MATLAB9/extern/include/ +mex -g -I/usr/local/cuda-7.5/include -L/usr/local/cuda-7.5/lib64 -lcudart -lcufft -lmwgpu Diffusion_4thO_GPU.cpp Diffus_4thO_GPU_core.o +movefile('Diffusion_4thO_GPU.mex*',Pathmove); + +delete TV_ROF_GPU_core* TV_FGP_GPU_core* TV_SB_GPU_core* dTV_FGP_GPU_core* NonlDiff_GPU_core* Diffus_4thO_GPU_core* CCPiDefines.h +fprintf('%s \n', 'All successfully compiled!'); + +pathA2 = sprintf(['..' filesep '..' filesep], 1i); +cd(pathA2); +cd demos
\ No newline at end of file diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/Diffusion_4thO.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/Diffusion_4thO.c new file mode 100644 index 0000000..81c0600 --- /dev/null +++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/Diffusion_4thO.c @@ -0,0 +1,76 @@ +/* + * This work is part of the Core Imaging Library developed by + * Visual Analytics and Imaging System Group of the Science Technology + * Facilities Council, STFC + * + * Copyright 2017 Daniil Kazantsev + * Copyright 2017 Srikanth Nagella, Edoardo Pasca + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * http://www.apache.org/licenses/LICENSE-2.0 + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +#include "matrix.h" +#include "mex.h" +#include "Diffus4th_order_core.h" + +/* C-OMP implementation of fourth-order diffusion scheme [1] for piecewise-smooth recovery (2D/3D case) + * The minimisation is performed using explicit scheme. + * + * Input Parameters: + * 1. Noisy image/volume [REQUIRED] + * 2. lambda - regularization parameter [REQUIRED] + * 3. Edge-preserving parameter (sigma) [REQUIRED] + * 4. Number of iterations, for explicit scheme >= 150 is recommended [OPTIONAL, default 300] + * 5. tau - time-marching step for the explicit scheme [OPTIONAL, default 0.015] + * + * Output: + * [1] Regularized image/volume + * + * This function is based on the paper by + * [1] Hajiaboli, M.R., 2011. An anisotropic fourth-order diffusion filter for image noise removal. International Journal of Computer Vision, 92(2), pp.177-191. + */ + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) + +{ + int number_of_dims, iter_numb, dimX, dimY, dimZ; + const int *dim_array; + float *Input, *Output=NULL, lambda, tau, sigma; + + dim_array = mxGetDimensions(prhs[0]); + number_of_dims = mxGetNumberOfDimensions(prhs[0]); + + /*Handling Matlab input data*/ + Input = (float *) mxGetData(prhs[0]); + lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */ + sigma = (float) mxGetScalar(prhs[2]); /* Edge-preserving parameter */ + iter_numb = 300; /* iterations number */ + tau = 0.01; /* marching step parameter */ + + if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } + if ((nrhs < 3) || (nrhs > 5)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant"); + if ((nrhs == 4) || (nrhs == 5)) iter_numb = (int) mxGetScalar(prhs[3]); /* iterations number */ + if (nrhs == 5) tau = (float) mxGetScalar(prhs[4]); /* marching step parameter */ + + /*Handling Matlab output data*/ + dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; + + /* output arrays*/ + if (number_of_dims == 2) { + dimZ = 1; /*2D case*/ + /* output image/volume */ + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + } + if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + + Diffus4th_CPU_main(Input, Output, lambda, sigma, iter_numb, tau, dimX, dimY, dimZ); +}
\ No newline at end of file diff --git a/Wrappers/Matlab/mex_compile/regularisers_GPU/Diffusion_4thO_GPU.cpp b/Wrappers/Matlab/mex_compile/regularisers_GPU/Diffusion_4thO_GPU.cpp new file mode 100644 index 0000000..0edc067 --- /dev/null +++ b/Wrappers/Matlab/mex_compile/regularisers_GPU/Diffusion_4thO_GPU.cpp @@ -0,0 +1,76 @@ +/* + * This work is part of the Core Imaging Library developed by + * Visual Analytics and Imaging System Group of the Science Technology + * Facilities Council, STFC + * + * Copyright 2017 Daniil Kazantsev + * Copyright 2017 Srikanth Nagella, Edoardo Pasca + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * http://www.apache.org/licenses/LICENSE-2.0 + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +#include "matrix.h" +#include "mex.h" +#include "Diffus_4thO_GPU_core.h" + +/* CUDA implementation of fourth-order diffusion scheme [1] for piecewise-smooth recovery (2D/3D case) + * The minimisation is performed using explicit scheme. + * + * Input Parameters: + * 1. Noisy image/volume [REQUIRED] + * 2. lambda - regularization parameter [REQUIRED] + * 3. Edge-preserving parameter (sigma) [REQUIRED] + * 4. Number of iterations, for explicit scheme >= 150 is recommended [OPTIONAL, default 300] + * 5. tau - time-marching step for the explicit scheme [OPTIONAL, default 0.015] + * + * Output: + * [1] Regularized image/volume + * + * This function is based on the paper by + * [1] Hajiaboli, M.R., 2011. An anisotropic fourth-order diffusion filter for image noise removal. International Journal of Computer Vision, 92(2), pp.177-191. + */ + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) + +{ + int number_of_dims, iter_numb, dimX, dimY, dimZ; + const int *dim_array; + float *Input, *Output=NULL, lambda, tau, sigma; + + dim_array = mxGetDimensions(prhs[0]); + number_of_dims = mxGetNumberOfDimensions(prhs[0]); + + /*Handling Matlab input data*/ + Input = (float *) mxGetData(prhs[0]); + lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */ + sigma = (float) mxGetScalar(prhs[2]); /* Edge-preserving parameter */ + iter_numb = 300; /* iterations number */ + tau = 0.01; /* marching step parameter */ + + if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } + if ((nrhs < 3) || (nrhs > 5)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant"); + if ((nrhs == 4) || (nrhs == 5)) iter_numb = (int) mxGetScalar(prhs[3]); /* iterations number */ + if (nrhs == 5) tau = (float) mxGetScalar(prhs[4]); /* marching step parameter */ + + /*Handling Matlab output data*/ + dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; + + /* output arrays*/ + if (number_of_dims == 2) { + dimZ = 1; /*2D case*/ + /* output image/volume */ + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + } + if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + + Diffus4th_GPU_main(Input, Output, lambda, sigma, iter_numb, tau, dimX, dimY, dimZ); +}
\ No newline at end of file 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 |