diff options
Diffstat (limited to 'Wrappers/Matlab')
| -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 |
6 files changed, 210 insertions, 8 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");} +} |