diff options
Diffstat (limited to 'Wrappers')
| -rw-r--r-- | Wrappers/Matlab/demos/demoMatlab_denoise.m | 42 | ||||
| -rw-r--r-- | Wrappers/Matlab/mex_compile/compileCPU_mex_Linux.m | 5 | ||||
| -rw-r--r-- | Wrappers/Matlab/mex_compile/compileCPU_mex_WINDOWS.m | 9 | ||||
| -rw-r--r-- | Wrappers/Matlab/mex_compile/compileGPU_mex.m | 6 | ||||
| -rw-r--r-- | Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c | 78 | ||||
| -rw-r--r-- | Wrappers/Matlab/mex_compile/regularisers_GPU/TGV_GPU.cpp | 78 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/filters/regularisers.py | 24 | ||||
| -rw-r--r-- | Wrappers/Python/demos/demo_cpu_regularisers.py | 67 | ||||
| -rw-r--r-- | Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py | 106 | ||||
| -rw-r--r-- | Wrappers/Python/demos/demo_gpu_regularisers.py | 70 | ||||
| -rw-r--r-- | Wrappers/Python/setup-regularisers.py.in | 1 | ||||
| -rw-r--r-- | Wrappers/Python/src/cpu_regularisers.pyx | 34 | ||||
| -rw-r--r-- | Wrappers/Python/src/gpu_regularisers.pyx | 34 |
13 files changed, 516 insertions, 38 deletions
diff --git a/Wrappers/Matlab/demos/demoMatlab_denoise.m b/Wrappers/Matlab/demos/demoMatlab_denoise.m index 8289f41..3f0ca54 100644 --- a/Wrappers/Matlab/demos/demoMatlab_denoise.m +++ b/Wrappers/Matlab/demos/demoMatlab_denoise.m @@ -1,9 +1,11 @@ % Image (2D) denoising demo using CCPi-RGL clear; close all -Path1 = sprintf(['..' filesep 'mex_compile' filesep 'installed'], 1i); -Path2 = sprintf(['..' filesep '..' filesep '..' filesep 'data' filesep], 1i); -addpath(Path1); -addpath(Path2); +fsep = '/'; + +Path1 = sprintf(['..' fsep 'mex_compile' fsep 'installed'], 1i); +Path2 = sprintf(['..' fsep '..' fsep '..' fsep 'data' fsep], 1i); +Path3 = sprintf(['..' fsep 'supp'], 1i); +addpath(Path1); addpath(Path2); addpath(Path3); Im = double(imread('lena_gray_512.tif'))/255; % loading image u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; @@ -16,6 +18,8 @@ tau_rof = 0.0025; % time-marching constant iter_rof = 750; % number of ROF iterations tic; u_rof = ROF_TV(single(u0), lambda_reg, iter_rof, tau_rof); toc; energyfunc_val_rof = TV_energy(single(u_rof),single(u0),lambda_reg, 1); % get energy function value +rmseROF = (RMSE(u_rof(:),Im(:))); +fprintf('%s %f \n', 'RMSE error for ROF-TV is:', rmseROF); figure; imshow(u_rof, [0 1]); title('ROF-TV denoised image (CPU)'); %% % fprintf('Denoise using the ROF-TV model (GPU) \n'); @@ -29,6 +33,8 @@ iter_fgp = 1000; % number of FGP iterations epsil_tol = 1.0e-06; % tolerance tic; u_fgp = FGP_TV(single(u0), lambda_reg, iter_fgp, epsil_tol); toc; energyfunc_val_fgp = TV_energy(single(u_fgp),single(u0),lambda_reg, 1); % get energy function value +rmseFGP = (RMSE(u_fgp(:),Im(:))); +fprintf('%s %f \n', 'RMSE error for FGP-TV is:', rmseFGP); figure; imshow(u_fgp, [0 1]); title('FGP-TV denoised image (CPU)'); %% @@ -43,6 +49,8 @@ iter_sb = 150; % number of SB iterations epsil_tol = 1.0e-06; % tolerance tic; u_sb = SB_TV(single(u0), lambda_reg, iter_sb, epsil_tol); toc; energyfunc_val_sb = TV_energy(single(u_sb),single(u0),lambda_reg, 1); % get energy function value +rmseSB = (RMSE(u_sb(:),Im(:))); +fprintf('%s %f \n', 'RMSE error for SB-TV is:', rmseSB); figure; imshow(u_sb, [0 1]); title('SB-TV denoised image (CPU)'); %% % fprintf('Denoise using the SB-TV model (GPU) \n'); @@ -51,12 +59,34 @@ figure; imshow(u_sb, [0 1]); title('SB-TV denoised image (CPU)'); % tic; u_sbG = SB_TV_GPU(single(u0), lambda_reg, iter_sb, epsil_tol); toc; % 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 +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 +% 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); +% figure; imshow(u_tgv_gpu, [0 1]); title('TGV denoised image (GPU)'); +%% fprintf('Denoise using Nonlinear-Diffusion model (CPU) \n'); iter_diff = 800; % number of diffusion iterations 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; +rmseDiffus = (RMSE(u_diff(:),Im(:))); +fprintf('%s %f \n', 'RMSE error for Nonlinear Diffusion is:', rmseDiffus); figure; imshow(u_diff, [0 1]); title('Diffusion denoised image (CPU)'); %% % fprintf('Denoise using Nonlinear-Diffusion model (GPU) \n'); @@ -73,6 +103,8 @@ lambda_regDiff = 3.5; % regularisation for the diffusivity sigmaPar = 0.02; % edge-preserving parameter tau_param = 0.0015; % time-marching constant tic; u_diff4 = Diffusion_4thO(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param); toc; +rmseDiffHO = (RMSE(u_diff4(:),Im(:))); +fprintf('%s %f \n', 'RMSE error for Fourth-order anisotropic diffusion is:', rmseDiffHO); figure; imshow(u_diff4, [0 1]); title('Diffusion 4thO denoised image (CPU)'); %% % fprintf('Denoise using Fourth-order anisotropic diffusion model (GPU) \n'); @@ -95,6 +127,8 @@ iter_fgp = 1000; % number of FGP iterations epsil_tol = 1.0e-06; % tolerance eta = 0.2; % Reference image gradient smoothing constant tic; u_fgp_dtv = FGP_dTV(single(u0), single(u_ref), lambda_reg, iter_fgp, epsil_tol, eta); toc; +rmse_dTV= (RMSE(u_fgp_dtv(:),Im(:))); +fprintf('%s %f \n', 'RMSE error for Directional Total Variation (dTV) is:', rmse_dTV); figure; imshow(u_fgp_dtv, [0 1]); title('FGP-dTV denoised image (CPU)'); %% % fprintf('Denoise using the FGP-dTV model (GPU) \n'); diff --git a/Wrappers/Matlab/mex_compile/compileCPU_mex_Linux.m b/Wrappers/Matlab/mex_compile/compileCPU_mex_Linux.m index 82b681a..8acc1b7 100644 --- a/Wrappers/Matlab/mex_compile/compileCPU_mex_Linux.m +++ b/Wrappers/Matlab/mex_compile/compileCPU_mex_Linux.m @@ -36,6 +36,9 @@ 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 TGV.c TGV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" +movefile('TGV.mex*',Pathmove); + mex TV_energy.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" movefile('TV_energy.mex*',Pathmove); @@ -46,7 +49,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* Diffus4th_order_core* CCPiDefines.h +delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* Diffus4th_order_core* TGV_core* CCPiDefines.h delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core* fprintf('%s \n', 'Regularisers successfully compiled!'); diff --git a/Wrappers/Matlab/mex_compile/compileCPU_mex_WINDOWS.m b/Wrappers/Matlab/mex_compile/compileCPU_mex_WINDOWS.m index 629a346..ea1ad7d 100644 --- a/Wrappers/Matlab/mex_compile/compileCPU_mex_WINDOWS.m +++ b/Wrappers/Matlab/mex_compile/compileCPU_mex_WINDOWS.m @@ -44,6 +44,9 @@ movefile('NonlDiff.mex*',Pathmove); mex Diffusion_4thO.c Diffus4th_order_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" movefile('Diffusion_4thO.mex*',Pathmove); +mex TGV.c TGV_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" +movefile('TGV.mex*',Pathmove); + mex TV_energy.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" movefile('TV_energy.mex*',Pathmove); @@ -54,7 +57,7 @@ movefile('NonlDiff_Inp.mex*',Pathmove); mex NonlocalMarching_Inpaint.c NonlocalMarching_Inpaint_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" movefile('NonlocalMarching_Inpaint.mex*',Pathmove); -delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* Diffus4th_order_core* CCPiDefines.h +delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* Diffus4th_order_core* TGV_core* CCPiDefines.h delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core* fprintf('%s \n', 'Regularisers successfully compiled!'); %% @@ -82,13 +85,15 @@ fprintf('%s \n', 'Regularisers successfully compiled!'); % movefile('NonlDiff.mex*',Pathmove); % mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" Diffusion_4thO.c Diffus4th_order_core.c utils.c % movefile('Diffusion_4thO.mex*',Pathmove); +% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TGV.c TGV_core.c utils.c +% movefile('TGV.mex*',Pathmove); % mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TV_energy.c utils.c % movefile('TV_energy.mex*',Pathmove); % mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" NonlDiff_Inp.c Diffusion_Inpaint_core.c utils.c % movefile('NonlDiff_Inp.mex*',Pathmove); % mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" NonlocalMarching_Inpaint.c NonlocalMarching_Inpaint_core.c utils.c % movefile('NonlocalMarching_Inpaint.mex*',Pathmove); -% delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* Diffus4th_order_core* CCPiDefines.h +% delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* Diffus4th_order_core* TGV_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 a4a636e..003c6ec 100644 --- a/Wrappers/Matlab/mex_compile/compileGPU_mex.m +++ b/Wrappers/Matlab/mex_compile/compileGPU_mex.m @@ -37,6 +37,10 @@ movefile('FGP_TV_GPU.mex*',Pathmove); mex -g -I/usr/local/cuda-7.5/include -L/usr/local/cuda-7.5/lib64 -lcudart -lcufft -lmwgpu SB_TV_GPU.cpp TV_SB_GPU_core.o movefile('SB_TV_GPU.mex*',Pathmove); +!/usr/local/cuda/bin/nvcc -O0 -c TGV_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 TGV_GPU.cpp TGV_GPU_core.o +movefile('TGV_GPU.mex*',Pathmove); + !/usr/local/cuda/bin/nvcc -O0 -c dTV_FGP_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 FGP_dTV_GPU.cpp dTV_FGP_GPU_core.o movefile('FGP_dTV_GPU.mex*',Pathmove); @@ -49,7 +53,7 @@ movefile('NonlDiff_GPU.mex*',Pathmove); 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 +delete TV_ROF_GPU_core* TV_FGP_GPU_core* TV_SB_GPU_core* dTV_FGP_GPU_core* NonlDiff_GPU_core* Diffus_4thO_GPU_core* TGV_GPU_core* CCPiDefines.h fprintf('%s \n', 'All successfully compiled!'); pathA2 = sprintf(['..' fsep '..' fsep], 1i); diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c new file mode 100644 index 0000000..9516869 --- /dev/null +++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c @@ -0,0 +1,78 @@ +/* +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 "mex.h" +#include "TGV_core.h" + +/* C-OMP implementation of Primal-Dual denoising method for + * Total Generilized Variation (TGV)-L2 model [1] (2D case only) + * + * 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) + * 6. Lipshitz constant (default is 12) + * + * Output: + * Filtered/regulariaed image + * + * References: + * [1] K. Bredies "Total Generalized Variation" + */ + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) + +{ + int number_of_dims, iter, dimX, dimY; + const int *dim_array; + float *Input, *Output=NULL, lambda, alpha0, alpha1, L2; + + number_of_dims = mxGetNumberOfDimensions(prhs[0]); + dim_array = mxGetDimensions(prhs[0]); + + /*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) */ + 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 */ + iter = 300; /* Iterations number */ + L2 = 12.0f; /* Lipshitz constant */ + + if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } + if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6)) alpha1 = (float) mxGetScalar(prhs[2]); /* parameter to control the first-order term */ + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6)) alpha0 = (float) mxGetScalar(prhs[3]); /* parameter to control the second-order term */ + if ((nrhs == 5) || (nrhs == 6)) iter = (int) mxGetScalar(prhs[4]); /* Iterations number */ + if (nrhs == 6) L2 = (float) mxGetScalar(prhs[5]); /* Lipshitz constant */ + + /*Handling Matlab output data*/ + dimX = dim_array[0]; dimY = dim_array[1]; + + 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); + } + if (number_of_dims == 3) {mexErrMsgTxt("Only 2D images accepted");} +} diff --git a/Wrappers/Matlab/mex_compile/regularisers_GPU/TGV_GPU.cpp b/Wrappers/Matlab/mex_compile/regularisers_GPU/TGV_GPU.cpp new file mode 100644 index 0000000..5a0df5b --- /dev/null +++ b/Wrappers/Matlab/mex_compile/regularisers_GPU/TGV_GPU.cpp @@ -0,0 +1,78 @@ +/* +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 "mex.h" +#include "TGV_GPU_core.h" + +/* CUDA implementation of Primal-Dual denoising method for + * Total Generilized Variation (TGV)-L2 model [1] (2D case only) + * + * 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) + * 6. Lipshitz constant (default is 12) + * + * Output: + * Filtered/regulariaed image + * + * References: + * [1] K. Bredies "Total Generalized Variation" + */ + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) + +{ + int number_of_dims, iter, dimX, dimY; + const int *dim_array; + float *Input, *Output=NULL, lambda, alpha0, alpha1, L2; + + number_of_dims = mxGetNumberOfDimensions(prhs[0]); + dim_array = mxGetDimensions(prhs[0]); + + /*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) */ + 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 */ + iter = 300; /* Iterations number */ + L2 = 12.0f; /* Lipshitz constant */ + + if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } + if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6)) alpha1 = (float) mxGetScalar(prhs[2]); /* parameter to control the first-order term */ + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6)) alpha0 = (float) mxGetScalar(prhs[3]); /* parameter to control the second-order term */ + if ((nrhs == 5) || (nrhs == 6)) iter = (int) mxGetScalar(prhs[4]); /* Iterations number */ + if (nrhs == 6) L2 = (float) mxGetScalar(prhs[5]); /* Lipshitz constant */ + + /*Handling Matlab output data*/ + dimX = dim_array[0]; dimY = dim_array[1]; + + if (number_of_dims == 2) { + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + /* running the function */ + TGV_GPU_main(Input, Output, lambda, alpha1, alpha0, iter, L2, dimX, dimY); + } + if (number_of_dims == 3) {mexErrMsgTxt("Only 2D images accepted");} +} diff --git a/Wrappers/Python/ccpi/filters/regularisers.py b/Wrappers/Python/ccpi/filters/regularisers.py index 0b79dac..0e435a6 100644 --- a/Wrappers/Python/ccpi/filters/regularisers.py +++ b/Wrappers/Python/ccpi/filters/regularisers.py @@ -2,8 +2,8 @@ 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, 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 TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, Diff4th_CPU, TGV_CPU +from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU, Diff4th_GPU, TGV_GPU from ccpi.filters.cpu_regularisers import NDF_INPAINT_CPU, NVM_INPAINT_CPU def ROF_TV(inputData, regularisation_parameter, iterations, @@ -128,6 +128,26 @@ def DIFF4th(inputData, regularisation_parameter, edge_parameter, iterations, else: raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ .format(device)) +def TGV(inputData, regularisation_parameter, alpha1, alpha0, iterations, + LipshitzConst, device='cpu'): + if device == 'cpu': + return TGV_CPU(inputData, + regularisation_parameter, + alpha1, + alpha0, + iterations, + LipshitzConst) + elif device == 'gpu': + return TGV_GPU(inputData, + regularisation_parameter, + alpha1, + alpha0, + iterations, + LipshitzConst) + 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 ff500ae..5c20244 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, DIFF4th +from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, FGP_dTV, TNV, NDF, DIFF4th from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -74,7 +74,7 @@ print ("_______________ROF-TV (2D)_________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(1) +fig = plt.figure() plt.suptitle('Performance of ROF-TV regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -109,13 +109,13 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(rof_cpu, cmap="gray") plt.title('{}'.format('CPU results')) - +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________FGP-TV (2D)__________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(2) +fig = plt.figure() plt.suptitle('Performance of FGP-TV regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -159,12 +159,13 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(fgp_cpu, cmap="gray") plt.title('{}'.format('CPU results')) +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________SB-TV (2D)__________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(3) +fig = plt.figure() plt.suptitle('Performance of SB-TV regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -205,14 +206,62 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) imgplot = plt.imshow(sb_cpu, cmap="gray") plt.title('{}'.format('CPU results')) +#%% + +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("_____Total Generalised Variation (2D)______") +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(u0,cmap="gray") + +# set parameters +pars = {'algorithm' : TGV, \ + 'input' : u0,\ + 'regularisation_parameter':0.04, \ + 'alpha1':1.0,\ + 'alpha0':0.7,\ + 'number_of_iterations' :250 ,\ + 'LipshitzConstant' :12 ,\ + } + +print ("#############TGV CPU####################") +start_time = timeit.default_timer() +tgv_cpu = TGV(pars['input'], + pars['regularisation_parameter'], + pars['alpha1'], + pars['alpha0'], + pars['number_of_iterations'], + pars['LipshitzConstant'],'cpu') + + +rms = rmse(Im, tgv_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(tgv_cpu, cmap="gray") +plt.title('{}'.format('CPU results')) +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("________________NDF (2D)___________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(4) +fig = plt.figure() plt.suptitle('Performance of NDF regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -259,7 +308,7 @@ print ("___Anisotropic Diffusion 4th Order (2D)____") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(5) +fig = plt.figure() plt.suptitle('Performance of DIFF4th regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -304,7 +353,7 @@ print ("_____________FGP-dTV (2D)__________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(6) +fig = plt.figure() plt.suptitle('Performance of FGP-dTV regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -356,7 +405,7 @@ print ("__________Total nuclear Variation__________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(7) +fig = plt.figure() plt.suptitle('Performance of TNV regulariser using the CPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') diff --git a/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py index 4611522..46b8ffc 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, DIFF4th +from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, FGP_dTV, NDF, DIFF4th from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -50,12 +50,13 @@ u_ref = Im + np.random.normal(loc = 0 , u0 = u0.astype('float32') u_ref = u_ref.astype('float32') +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("____________ROF-TV bench___________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(1) +fig = plt.figure() plt.suptitle('Comparison of ROF-TV regulariser using CPU and GPU implementations') a=fig.add_subplot(1,4,1) a.set_title('Noisy Image') @@ -90,7 +91,6 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(rof_cpu, cmap="gray") plt.title('{}'.format('CPU results')) - print ("##############ROF TV GPU##################") start_time = timeit.default_timer() rof_gpu = ROF_TV(pars['input'], @@ -128,12 +128,13 @@ if (diff_im.sum() > 1): else: print ("Arrays match") +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("____________FGP-TV bench___________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(2) +fig = plt.figure() plt.suptitle('Comparison of FGP-TV regulariser using CPU and GPU implementations') a=fig.add_subplot(1,4,1) a.set_title('Noisy Image') @@ -218,12 +219,13 @@ if (diff_im.sum() > 1): else: print ("Arrays match") +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("____________SB-TV bench___________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(3) +fig = plt.figure() plt.suptitle('Comparison of SB-TV regulariser using CPU and GPU implementations') a=fig.add_subplot(1,4,1) a.set_title('Noisy Image') @@ -303,14 +305,98 @@ if (diff_im.sum() > 1): print ("Arrays do not match!") else: print ("Arrays match") +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("____________TGV bench___________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure() +plt.suptitle('Comparison of TGV 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' : TGV, \ + 'input' : u0,\ + 'regularisation_parameter':0.04, \ + 'alpha1':1.0,\ + 'alpha0':0.7,\ + 'number_of_iterations' :250 ,\ + 'LipshitzConstant' :12 ,\ + } + +print ("#############TGV CPU####################") +start_time = timeit.default_timer() +tgv_cpu = TGV(pars['input'], + pars['regularisation_parameter'], + pars['alpha1'], + pars['alpha0'], + pars['number_of_iterations'], + pars['LipshitzConstant'],'cpu') + +rms = rmse(Im, tgv_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(tgv_cpu, cmap="gray") +plt.title('{}'.format('CPU results')) + +print ("##############SB TV GPU##################") +start_time = timeit.default_timer() +tgv_gpu = TGV(pars['input'], + pars['regularisation_parameter'], + pars['alpha1'], + pars['alpha0'], + pars['number_of_iterations'], + pars['LipshitzConstant'],'gpu') + +rms = rmse(Im, tgv_gpu) +pars['rmse'] = rms +pars['algorithm'] = TGV +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(tgv_gpu, cmap="gray") +plt.title('{}'.format('GPU results')) +print ("--------Compare the results--------") +tolerance = 1e-05 +diff_im = np.zeros(np.shape(tgv_gpu)) +diff_im = abs(tgv_cpu - tgv_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 ("_______________NDF bench___________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(4) +fig = plt.figure() plt.suptitle('Comparison of NDF regulariser using CPU and GPU implementations') a=fig.add_subplot(1,4,1) a.set_title('Noisy Image') @@ -390,13 +476,13 @@ if (diff_im.sum() > 1): else: print ("Arrays match") - +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("___Anisotropic Diffusion 4th Order (2D)____") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(5) +fig = plt.figure() plt.suptitle('Comparison of Diff4th regulariser using CPU and GPU implementations') a=fig.add_subplot(1,4,1) a.set_title('Noisy Image') @@ -472,12 +558,13 @@ if (diff_im.sum() > 1): else: print ("Arrays match") +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("____________FGP-dTV bench___________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(6) +fig = plt.figure() plt.suptitle('Comparison of FGP-dTV regulariser using CPU and GPU implementations') a=fig.add_subplot(1,4,1) a.set_title('Noisy Image') @@ -565,3 +652,4 @@ if (diff_im.sum() > 1): print ("Arrays do not match!") else: print ("Arrays match") +#%%
\ No newline at end of file diff --git a/Wrappers/Python/demos/demo_gpu_regularisers.py b/Wrappers/Python/demos/demo_gpu_regularisers.py index 3179428..792a019 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, DIFF4th +from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, FGP_dTV, NDF, DIFF4th from qualitymetrics import rmse ############################################################################### def printParametersToString(pars): @@ -66,13 +66,14 @@ Im2[:,0:M] = Im[:,0:M] Im = Im2 del Im2 """ +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("____________ROF-TV regulariser_____________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(1) +fig = plt.figure() plt.suptitle('Performance of the ROF-TV regulariser using the GPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -108,13 +109,13 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(rof_gpu, cmap="gray") plt.title('{}'.format('GPU results')) - +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("____________FGP-TV regulariser_____________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(2) +fig = plt.figure() plt.suptitle('Performance of the FGP-TV regulariser using the GPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -157,13 +158,13 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(fgp_gpu, cmap="gray") plt.title('{}'.format('GPU results')) - +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("____________SB-TV regulariser______________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(3) +fig = plt.figure() plt.suptitle('Performance of the SB-TV regulariser using the GPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -203,14 +204,62 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) imgplot = plt.imshow(sb_gpu, cmap="gray") plt.title('{}'.format('GPU results')) +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("_____Total Generalised Variation (2D)______") +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(u0,cmap="gray") + +# set parameters +pars = {'algorithm' : TGV, \ + 'input' : u0,\ + 'regularisation_parameter':0.04, \ + 'alpha1':1.0,\ + 'alpha0':0.7,\ + 'number_of_iterations' :250 ,\ + 'LipshitzConstant' :12 ,\ + } + +print ("#############TGV CPU####################") +start_time = timeit.default_timer() +tgv_gpu = TGV(pars['input'], + pars['regularisation_parameter'], + pars['alpha1'], + pars['alpha0'], + pars['number_of_iterations'], + pars['LipshitzConstant'],'gpu') + + +rms = rmse(Im, tgv_gpu) +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_gpu, cmap="gray") +plt.title('{}'.format('GPU results')) + +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________NDF regulariser_____________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(4) +fig = plt.figure() plt.suptitle('Performance of the NDF regulariser using the GPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -251,13 +300,13 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(ndf_gpu, cmap="gray") plt.title('{}'.format('GPU results')) - +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("___Anisotropic Diffusion 4th Order (2D)____") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(5) +fig = plt.figure() plt.suptitle('Performance of DIFF4th regulariser using the GPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') @@ -296,12 +345,13 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(diff4_gpu, cmap="gray") plt.title('{}'.format('GPU results')) +#%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("____________FGP-dTV bench___________________") print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") ## plot -fig = plt.figure(6) +fig = plt.figure() plt.suptitle('Performance of the FGP-dTV regulariser using the GPU') a=fig.add_subplot(1,2,1) a.set_title('Noisy Image') diff --git a/Wrappers/Python/setup-regularisers.py.in b/Wrappers/Python/setup-regularisers.py.in index 76dfecf..89ebaf9 100644 --- a/Wrappers/Python/setup-regularisers.py.in +++ b/Wrappers/Python/setup-regularisers.py.in @@ -38,6 +38,7 @@ extra_include_dirs += [os.path.join(".." , ".." , "Core"), os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_FGP" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_ROF" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TV_SB" ) , + os.path.join(".." , ".." , "Core", "regularisers_GPU" , "TGV" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "NDF" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "dTV_FGP" ) , os.path.join(".." , ".." , "Core", "regularisers_GPU" , "DIFF4th" ) , diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx index bdb1eff..cf81bec 100644 --- a/Wrappers/Python/src/cpu_regularisers.pyx +++ b/Wrappers/Python/src/cpu_regularisers.pyx @@ -21,6 +21,7 @@ cimport numpy as np cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ); 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 TGV_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY); 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); @@ -189,6 +190,39 @@ def TV_SB_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, printM, dims[2], dims[1], dims[0]) return outputData + +#***************************************************************# +#***************** Total Generalised Variation *****************# +#***************************************************************# +def TGV_CPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst): + if inputData.ndim == 2: + return TGV_2D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst) + elif inputData.ndim == 3: + return 0 + +def TGV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularisation_parameter, + float alpha1, + float alpha0, + int iterationsNumb, + float LipshitzConst): + + 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 TGV iterations for 2D data */ + TGV_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, + alpha1, + alpha0, + iterationsNumb, + LipshitzConst, + dims[1],dims[0]) + return outputData + #****************************************************************# #**************Directional Total-variation FGP ******************# #****************************************************************# diff --git a/Wrappers/Python/src/gpu_regularisers.pyx b/Wrappers/Python/src/gpu_regularisers.pyx index b67e62b..4a202d7 100644 --- a/Wrappers/Python/src/gpu_regularisers.pyx +++ b/Wrappers/Python/src/gpu_regularisers.pyx @@ -21,6 +21,7 @@ cimport numpy as np cdef extern void TV_ROF_GPU_main(float* Input, float* Output, float lambdaPar, int iter, float tau, int N, int M, int Z); cdef extern void 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 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 TGV_GPU_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY); 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); @@ -86,6 +87,12 @@ def TV_SB_GPU(inputData, tolerance_param, methodTV, printM) +# Total Generilised Variation (TGV) +def TGV_GPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst): + if inputData.ndim == 2: + return TGV2D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst) + elif inputData.ndim == 3: + return 0 # Directional Total-variation Fast-Gradient-Projection (FGP) def dTV_FGP_GPU(inputData, refdata, @@ -315,6 +322,33 @@ def SBTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, dims[2], dims[1], dims[0]); return outputData + +#***************************************************************# +#***************** Total Generalised Variation *****************# +#***************************************************************# +def TGV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularisation_parameter, + float alpha1, + float alpha0, + int iterationsNumb, + float LipshitzConst): + + 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 TGV iterations for 2D data */ + TGV_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, + alpha1, + alpha0, + iterationsNumb, + LipshitzConst, + dims[1],dims[0]) + return outputData + #****************************************************************# #**************Directional Total-variation FGP ******************# #****************************************************************# |