diff options
18 files changed, 319 insertions, 248 deletions
diff --git a/demos/demoMatlab_denoise.m b/demos/demoMatlab_denoise.m index a22b40a..7581068 100644 --- a/demos/demoMatlab_denoise.m +++ b/demos/demoMatlab_denoise.m @@ -19,7 +19,7 @@ iter_rof = 2000; % number of ROF iterations tau_rof = 0.001; % time-marching constant epsil_tol = 0.0; % tolerance tic; [u_rof,infovec] = ROF_TV(single(u0), lambda_reg, iter_rof, tau_rof, epsil_tol); toc; -%energyfunc_val_rof = TV_energy(single(u_rof),single(u0),lambda_reg, 1); % get energy function value +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); [ssimval] = ssim(u_rof*255,single(Im)*255); @@ -27,7 +27,7 @@ fprintf('%s %f \n', 'MSSIM error for ROF-TV is:', ssimval); figure; imshow(u_rof, [0 1]); title('ROF-TV denoised image (CPU)'); %% % fprintf('Denoise using the ROF-TV model (GPU) \n'); -% tic; u_rofG = ROF_TV_GPU(single(u0), lambda_reg, iter_rof, tau_rof, epsil_tol); toc; +% tic; [u_rofG,infovec] = ROF_TV_GPU(single(u0), lambda_reg, iter_rof, tau_rof, epsil_tol); toc; % figure; imshow(u_rofG, [0 1]); title('ROF-TV denoised image (GPU)'); %% fprintf('Denoise using the FGP-TV model (CPU) \n'); @@ -59,8 +59,6 @@ fprintf('%s %f \n', 'MSSIM error for SB-TV is:', ssimval); figure; imshow(u_sb, [0 1]); title('SB-TV denoised image (CPU)'); %% % fprintf('Denoise using the SB-TV model (GPU) \n'); -% iter_sb = 80; % number of SB iterations -% epsil_tol = 1.0e-06; % tolerance % 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)'); %% @@ -69,7 +67,8 @@ iter_diff = 450; % number of diffusion iterations lambda_regDiff = 0.025; % regularisation for the diffusivity sigmaPar = 0.015; % edge-preserving parameter tau_param = 0.02; % time-marching constant -tic; u_diff = NonlDiff(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber'); toc; +epsil_tol = 1.0e-06; % tolerance +tic; [u_diff,infovec] = NonlDiff(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber', epsil_tol); toc; rmseDiffus = (RMSE(u_diff(:),Im(:))); fprintf('%s %f \n', 'RMSE error for Nonlinear Diffusion is:', rmseDiffus); [ssimval] = ssim(u_diff*255,single(Im)*255); @@ -85,11 +84,11 @@ figure; imshow(u_diff, [0 1]); title('Diffusion denoised image (CPU)'); % figure; imshow(u_diff_g, [0 1]); title('Diffusion denoised image (GPU)'); %% fprintf('Denoise using the TGV model (CPU) \n'); -lambda_TGV = 0.034; % regularisation parameter +lambda_TGV = 0.045; % regularisation parameter alpha1 = 1.0; % parameter to control the first-order term -alpha0 = 1.0; % 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; +alpha0 = 2.0; % parameter to control the second-order term +iter_TGV = 2500; % number of Primal-Dual iterations for TGV +tic; [u_tgv,infovec] = 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); [ssimval] = ssim(u_tgv*255,single(Im)*255); @@ -109,11 +108,12 @@ figure; imshow(u_tgv, [0 1]); title('TGV denoised image (CPU)'); % figure; imshow(u_tgv_gpu, [0 1]); title('TGV denoised image (GPU)'); %% fprintf('Denoise using the ROF-LLT model (CPU) \n'); -lambda_ROF = 0.016; % ROF regularisation parameter -lambda_LLT = lambda_reg*0.25; % LLT regularisation parameter -iter_LLT = 500; % iterations +lambda_ROF = 0.02; % ROF regularisation parameter +lambda_LLT = 0.01; % LLT regularisation parameter +iter_LLT = 1000; % iterations tau_rof_llt = 0.0025; % time-marching constant -tic; u_rof_llt = LLT_ROF(single(u0), lambda_ROF, lambda_LLT, iter_LLT, tau_rof_llt); toc; +epsil_tol = 1.0e-06; % tolerance +tic; [u_rof_llt,infovec] = LLT_ROF(single(u0), lambda_ROF, lambda_LLT, iter_LLT, tau_rof_llt,epsil_tol); toc; rmseROFLLT = (RMSE(u_rof_llt(:),Im(:))); fprintf('%s %f \n', 'RMSE error for TGV is:', rmseROFLLT); [ssimval] = ssim(u_rof_llt*255,single(Im)*255); @@ -132,10 +132,11 @@ figure; imshow(u_rof_llt, [0 1]); title('ROF-LLT denoised image (CPU)'); %% 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.021; % edge-preserving parameter +lambda_regDiff = 2.5; % regularisation for the diffusivity +sigmaPar = 0.03; % 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; +epsil_tol = 1.0e-06; % tolerance +tic; [u_diff4,infovec] = Diffusion_4thO(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param, epsil_tol); toc; rmseDiffHO = (RMSE(u_diff4(:),Im(:))); fprintf('%s %f \n', 'RMSE error for Fourth-order anisotropic diffusion is:', rmseDiffHO); [ssimval] = ssim(u_diff4*255,single(Im)*255); @@ -174,10 +175,11 @@ fprintf('Denoise using the FGP-dTV model (CPU) \n'); u_ref = Im + .01*randn(size(Im)); u_ref(u_ref < 0) = 0; % u_ref = zeros(size(Im),'single'); % pass zero reference (dTV -> TV) +lambda_reg = 0.04; 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; +tic; [u_fgp_dtv,infovec] = 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)'); diff --git a/src/Matlab/mex_compile/compileCPU_mex_Linux.m b/src/Matlab/mex_compile/compileCPU_mex_Linux.m index d8035f4..19fb1a5 100644 --- a/src/Matlab/mex_compile/compileCPU_mex_Linux.m +++ b/src/Matlab/mex_compile/compileCPU_mex_Linux.m @@ -1,4 +1,5 @@ -% execute this mex file on Linux in Matlab once +% execute this mex file on Linux in Matlab once. See also Cmake driven +% build if this fails fsep = '/'; @@ -27,55 +28,54 @@ movefile('FGP_TV.mex*',Pathmove); fprintf('%s \n', 'Compiling SB-TV...'); mex SB_TV.c SB_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" movefile('SB_TV.mex*',Pathmove); -% -% fprintf('%s \n', 'Compiling dFGP-TV...'); -% mex FGP_dTV.c FGP_dTV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -% movefile('FGP_dTV.mex*',Pathmove); -% -% fprintf('%s \n', 'Compiling TNV...'); -% mex TNV.c TNV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -% movefile('TNV.mex*',Pathmove); -% -% fprintf('%s \n', 'Compiling NonLinear Diffusion...'); -% mex NonlDiff.c Diffusion_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -% movefile('NonlDiff.mex*',Pathmove); -% -% fprintf('%s \n', 'Compiling Anisotropic diffusion of higher order...'); -% mex Diffusion_4thO.c Diffus4th_order_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -% movefile('Diffusion_4thO.mex*',Pathmove); -% -% fprintf('%s \n', 'Compiling TGV...'); -% mex TGV.c TGV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -% movefile('TGV.mex*',Pathmove); -% -% fprintf('%s \n', 'Compiling ROF-LLT...'); -% mex LLT_ROF.c LLT_ROF_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -% movefile('LLT_ROF.mex*',Pathmove); -% -% fprintf('%s \n', 'Compiling NonLocal-TV...'); -% mex PatchSelect.c PatchSelect_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -% mex Nonlocal_TV.c Nonlocal_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -% movefile('Nonlocal_TV.mex*',Pathmove); -% movefile('PatchSelect.mex*',Pathmove); -% -% fprintf('%s \n', 'Compiling additional tools...'); -% mex TV_energy.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -% movefile('TV_energy.mex*',Pathmove); + +fprintf('%s \n', 'Compiling dFGP-TV...'); +mex FGP_dTV.c FGP_dTV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" +movefile('FGP_dTV.mex*',Pathmove); + +fprintf('%s \n', 'Compiling TNV...'); +mex TNV.c TNV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" +movefile('TNV.mex*',Pathmove); + +fprintf('%s \n', 'Compiling NonLinear Diffusion...'); +mex NonlDiff.c Diffusion_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" +movefile('NonlDiff.mex*',Pathmove); + +fprintf('%s \n', 'Compiling Anisotropic diffusion of higher order...'); +mex Diffusion_4thO.c Diffus4th_order_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" +movefile('Diffusion_4thO.mex*',Pathmove); + +fprintf('%s \n', 'Compiling TGV...'); +mex TGV.c TGV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" +movefile('TGV.mex*',Pathmove); + +fprintf('%s \n', 'Compiling ROF-LLT...'); +mex LLT_ROF.c LLT_ROF_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" +movefile('LLT_ROF.mex*',Pathmove); + +fprintf('%s \n', 'Compiling NonLocal-TV...'); +mex PatchSelect.c PatchSelect_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" +mex Nonlocal_TV.c Nonlocal_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" +movefile('Nonlocal_TV.mex*',Pathmove); +movefile('PatchSelect.mex*',Pathmove); + +fprintf('%s \n', 'Compiling additional tools...'); +mex TV_energy.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" +movefile('TV_energy.mex*',Pathmove); %############Inpainters##############% -% fprintf('%s \n', 'Compiling Nonlinear/Linear diffusion inpainting...'); -% mex NonlDiff_Inp.c Diffusion_Inpaint_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -% movefile('NonlDiff_Inp.mex*',Pathmove); -% -% fprintf('%s \n', 'Compiling Nonlocal marching method for inpainting...'); -% mex NonlocalMarching_Inpaint.c NonlocalMarching_Inpaint_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -% movefile('NonlocalMarching_Inpaint.mex*',Pathmove); -% +fprintf('%s \n', 'Compiling Nonlinear/Linear diffusion inpainting...'); +mex NonlDiff_Inp.c Diffusion_Inpaint_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" +movefile('NonlDiff_Inp.mex*',Pathmove); + +fprintf('%s \n', 'Compiling Nonlocal marching method for inpainting...'); +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* TGV_core* LLT_ROF_core* CCPiDefines.h delete PatchSelect_core* Nonlocal_TV_core* delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core* fprintf('%s \n', '<<<<<<< Regularisers successfully compiled! >>>>>>>'); pathA2 = sprintf(['..' fsep '..' fsep '..' fsep '..' fsep 'demos'], 1i); -cd(pathA2); - +cd(pathA2);
\ No newline at end of file diff --git a/src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m b/src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m index 6f7541c..3a9e2af 100644 --- a/src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m +++ b/src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m @@ -10,9 +10,9 @@ fsep = '/'; -pathcopyFrom = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'regularisers_CPU'], 1i); -pathcopyFrom1 = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'CCPiDefines.h'], 1i); -pathcopyFrom2 = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'inpainters_CPU'], 1i); +pathcopyFrom = sprintf(['..' fsep '..' fsep 'Core' fsep 'regularisers_CPU'], 1i); +pathcopyFrom1 = sprintf(['..' fsep '..' fsep 'Core' fsep 'CCPiDefines.h'], 1i); +pathcopyFrom2 = sprintf(['..' fsep '..' fsep 'Core' fsep 'inpainters_CPU'], 1i); copyfile(pathcopyFrom, 'regularisers_CPU'); copyfile(pathcopyFrom1, 'regularisers_CPU'); @@ -79,7 +79,6 @@ fprintf('%s \n', 'Compiling Nonlocal marching method for inpaiting...'); mex NonlocalMarching_Inpaint.c NonlocalMarching_Inpaint_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" movefile('NonlocalMarching_Inpaint.mex*',Pathmove); - %% %%% The second approach to compile using TDM-GCC which follows this %%% discussion: @@ -129,7 +128,5 @@ fprintf('%s \n', 'Regularisers successfully compiled!'); %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%pathA2 = sprintf(['..' fsep '..' fsep], 1i); -%cd(pathA2); -%cd demos +pathA2 = sprintf(['..' fsep '..' fsep '..' fsep '..' fsep 'demos'], 1i); +cd(pathA2); diff --git a/src/Matlab/mex_compile/compileGPU_mex.m b/src/Matlab/mex_compile/compileGPU_mex.m index dd1475c..3a7ac7c 100644 --- a/src/Matlab/mex_compile/compileGPU_mex.m +++ b/src/Matlab/mex_compile/compileGPU_mex.m @@ -8,13 +8,14 @@ % ! paths to matlab and CUDA sdk can be different, modify accordingly ! % Tested on Ubuntu 18.04/MATLAB 2016b/cuda10.0/gcc7.3 - +%>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<< % Installation HAS NOT been tested on Windows, please you Cmake build or % modify the code bellow accordingly fsep = '/'; -pathcopyFrom = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'regularisers_GPU'], 1i); -pathcopyFrom1 = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'CCPiDefines.h'], 1i); +pathcopyFrom = sprintf(['..' fsep '..' fsep 'Core' fsep 'regularisers_GPU'], 1i); +pathcopyFrom1 = sprintf(['..' fsep '..' fsep 'Core' fsep 'CCPiDefines.h'], 1i); + copyfile(pathcopyFrom, 'regularisers_GPU'); copyfile(pathcopyFrom1, 'regularisers_GPU'); @@ -69,6 +70,5 @@ movefile('LLT_ROF_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* TGV_GPU_core* LLT_ROF_GPU_core* CCPiDefines.h fprintf('%s \n', 'All successfully compiled!'); -pathA2 = sprintf(['..' fsep '..' fsep], 1i); -cd(pathA2); -cd demos
\ No newline at end of file +pathA2 = sprintf(['..' fsep '..' fsep '..' fsep '..' fsep 'demos'], 1i); +cd(pathA2);
\ No newline at end of file diff --git a/src/Matlab/mex_compile/regularisers_CPU/Diffusion_4thO.c b/src/Matlab/mex_compile/regularisers_CPU/Diffusion_4thO.c index 66ea9be..a003596 100644 --- a/src/Matlab/mex_compile/regularisers_CPU/Diffusion_4thO.c +++ b/src/Matlab/mex_compile/regularisers_CPU/Diffusion_4thO.c @@ -29,9 +29,11 @@ * 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] + * 6. eplsilon - tolerance constant [OPTIONAL parameter] * * Output: * [1] Regularized image/volume + * [2] Information vector which contains [iteration no., reached tolerance] * * 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. @@ -45,7 +47,8 @@ void mexFunction( int number_of_dims, iter_numb; mwSize dimX, dimY, dimZ; const mwSize *dim_array; - float *Input, *Output=NULL, lambda, tau, sigma; + float *Input, *Output=NULL, lambda, tau, sigma, epsil; + float *infovec=NULL; dim_array = mxGetDimensions(prhs[0]); number_of_dims = mxGetNumberOfDimensions(prhs[0]); @@ -54,13 +57,15 @@ void mexFunction( 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 */ + iter_numb = 500; /* iterations number */ + tau = 0.0025; /* marching step parameter */ + epsil = 1.0e-06; /*tolerance 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 */ + if ((nrhs < 3) || (nrhs > 6)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant, tolerance"); + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6)) iter_numb = (int) mxGetScalar(prhs[3]); /* iterations number */ + if ((nrhs == 5) || (nrhs == 6)) tau = (float) mxGetScalar(prhs[4]); /* marching step parameter */ + if (nrhs == 6) epsil = (float) mxGetScalar(prhs[5]); /* tolerance parameter */ /*Handling Matlab output data*/ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; @@ -73,5 +78,9 @@ void mexFunction( } 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); + int vecdim[1]; + vecdim[0] = 2; + infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL)); + + Diffus4th_CPU_main(Input, Output, infovec, lambda, sigma, iter_numb, tau, epsil, dimX, dimY, dimZ); }
\ No newline at end of file diff --git a/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c b/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c index f160b15..f6db6c8 100644 --- a/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c +++ b/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c @@ -58,7 +58,7 @@ void mexFunction( Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */ lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */ iter = 400; /* default iterations number */ - epsil = 0.0; /* default tolerance constant */ + epsil = 1.0e-06; /* default tolerance constant */ methTV = 0; /* default isotropic TV penalty */ nonneg = 0; /* default nonnegativity switch, off - 0 */ diff --git a/src/Matlab/mex_compile/regularisers_CPU/FGP_dTV.c b/src/Matlab/mex_compile/regularisers_CPU/FGP_dTV.c index 1a0c070..3122610 100644 --- a/src/Matlab/mex_compile/regularisers_CPU/FGP_dTV.c +++ b/src/Matlab/mex_compile/regularisers_CPU/FGP_dTV.c @@ -33,10 +33,10 @@ * 6. eta: smoothing constant to calculate gradient of the reference [OPTIONAL] * * 7. TV-type: methodTV - 'iso' (0) or 'l1' (1) [OPTIONAL] * 8. nonneg: 'nonnegativity (0 is OFF by default) [OPTIONAL] - * 9. print information: 0 (off) or 1 (on) [OPTIONAL] * * Output: - * [1] Filtered/regularized image/volume + * [1] Filtered/regularized image + * [2] Information vector which contains [iteration no., reached tolerance] * * This function is based on the Matlab's codes and papers by * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems" @@ -49,28 +49,27 @@ void mexFunction( int nrhs, const mxArray *prhs[]) { - int number_of_dims, iter, methTV, printswitch, nonneg; + int number_of_dims, iter, methTV, nonneg; mwSize dimX, dimY, dimZ; const mwSize *dim_array; const mwSize *dim_array2; - float *Input, *InputRef, *Output=NULL, lambda, epsil, eta; - + float *Input, *InputRef, *Output=NULL, *infovec=NULL, lambda, epsil, eta; + number_of_dims = mxGetNumberOfDimensions(prhs[0]); dim_array = mxGetDimensions(prhs[0]); dim_array2 = mxGetDimensions(prhs[1]); /*Handling Matlab input data*/ - if ((nrhs < 3) || (nrhs > 9)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Reference(2D/3D), Regularization parameter, iterations number, tolerance, smoothing constant, penalty type ('iso' or 'l1'), nonnegativity switch, print switch"); + if ((nrhs < 3) || (nrhs > 8)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Reference(2D/3D), Regularization parameter, iterations number, tolerance, smoothing constant, penalty type ('iso' or 'l1'), nonnegativity switch"); Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */ InputRef = (float *) mxGetData(prhs[1]); /* reference image (2D/3D) */ lambda = (float) mxGetScalar(prhs[2]); /* regularization parameter */ iter = 300; /* default iterations number */ - epsil = 0.0001; /* default tolerance constant */ + epsil = 1.0e-06; /* default tolerance constant */ eta = 0.01; /* default smoothing constant */ methTV = 0; /* default isotropic TV penalty */ nonneg = 0; /* default nonnegativity switch, off - 0 */ - printswitch = 0; /*default print is switched, off - 0 */ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } @@ -80,35 +79,32 @@ void mexFunction( dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; if (number_of_dims == 2) { if ((dimX != dim_array2[0]) || (dimY != dim_array2[1])) mexErrMsgTxt("The input images have different dimensionalities");} if (number_of_dims == 3) { if ((dimX != dim_array2[0]) || (dimY != dim_array2[1]) || (dimZ != dim_array2[2])) mexErrMsgTxt("The input volumes have different dimensionalities");} - - - if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8) || (nrhs == 9)) iter = (int) mxGetScalar(prhs[3]); /* iterations number */ - if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8) || (nrhs == 9)) epsil = (float) mxGetScalar(prhs[4]); /* tolerance constant */ - if ((nrhs == 6) || (nrhs == 7) || (nrhs == 8) || (nrhs == 9)) { - eta = (float) mxGetScalar(prhs[5]); /* smoothing constant for the gradient of InputRef */ - } - if ((nrhs == 7) || (nrhs == 8) || (nrhs == 9)) { + + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8)) iter = (int) mxGetScalar(prhs[3]); /* iterations number */ + if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8)) epsil = (float) mxGetScalar(prhs[4]); /* tolerance constant */ + if ((nrhs == 6) || (nrhs == 7) || (nrhs == 8)) eta = (float) mxGetScalar(prhs[5]); /* smoothing constant for the gradient of InputRef */ + if ((nrhs == 7) || (nrhs == 8)) { char *penalty_type; penalty_type = mxArrayToString(prhs[6]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */ if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',"); if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */ mxFree(penalty_type); } - if ((nrhs == 8) || (nrhs == 9)) { + if ((nrhs == 8)) { nonneg = (int) mxGetScalar(prhs[7]); if ((nonneg != 0) && (nonneg != 1)) mexErrMsgTxt("Nonnegativity constraint can be enabled by choosing 1 or off - 0"); } - if (nrhs == 9) { - printswitch = (int) mxGetScalar(prhs[8]); - if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0"); - } - + if (number_of_dims == 2) { dimZ = 1; /*2D case*/ 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)); + int vecdim[1]; + vecdim[0] = 2; + infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL)); + /* running the function */ - dTV_FGP_CPU_main(Input, InputRef, Output, lambda, iter, epsil, eta, methTV, nonneg, printswitch, dimX, dimY, dimZ); + dTV_FGP_CPU_main(Input, InputRef, Output, infovec, lambda, iter, epsil, eta, methTV, nonneg, dimX, dimY, dimZ); }
\ No newline at end of file diff --git a/src/Matlab/mex_compile/regularisers_CPU/LLT_ROF.c b/src/Matlab/mex_compile/regularisers_CPU/LLT_ROF.c index ab45446..f630397 100644 --- a/src/Matlab/mex_compile/regularisers_CPU/LLT_ROF.c +++ b/src/Matlab/mex_compile/regularisers_CPU/LLT_ROF.c @@ -32,9 +32,11 @@ * 3. lambdaLLT - LLT-related regularisation parameter * 4. tau - time-marching step * 5. iter - iterations number (for both models) +* 6. eplsilon - tolerance constant [OPTIONAL parameter] * * Output: -* Filtered/regularised image +* [1] Regularized image/volume +* [2] Information vector which contains [iteration no., reached tolerance] * * References: * [1] Lysaker, M., Lundervold, A. and Tai, X.C., 2003. Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time. IEEE Transactions on image processing, 12(12), pp.1579-1590. @@ -49,12 +51,13 @@ void mexFunction( int number_of_dims, iterationsNumb; mwSize dimX, dimY, dimZ; const mwSize *dim_array; - float *Input, *Output=NULL, lambdaROF, lambdaLLT, tau; + float *Input, *Output=NULL, lambdaROF, lambdaLLT, tau, epsil; + float *infovec=NULL; dim_array = mxGetDimensions(prhs[0]); number_of_dims = mxGetNumberOfDimensions(prhs[0]); - if ((nrhs < 3) || (nrhs > 5)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter (ROF), Regularisation parameter (LTT), iterations number, time-marching parameter"); + if ((nrhs < 3) || (nrhs > 6)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter (ROF), Regularisation parameter (LTT), iterations number, time-marching parameter, tolerance"); /*Handling Matlab input data*/ Input = (float *) mxGetData(prhs[0]); @@ -62,10 +65,12 @@ void mexFunction( lambdaLLT = (float) mxGetScalar(prhs[2]); /* ROF regularization parameter */ iterationsNumb = 250; tau = 0.0025; + epsil = 1.0e-06; /*tolerance parameter*/ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } - if ((nrhs == 4) || (nrhs == 5)) iterationsNumb = (int) mxGetScalar(prhs[3]); /* iterations number */ - if (nrhs == 5) tau = (float) mxGetScalar(prhs[4]); /* marching step parameter */ + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6)) iterationsNumb = (int) mxGetScalar(prhs[3]); /* iterations number */ + if ((nrhs == 5) || (nrhs == 6)) tau = (float) mxGetScalar(prhs[4]); /* marching step parameter */ + if (nrhs == 6) epsil = (float) mxGetScalar(prhs[5]); /* epsilon */ /*Handling Matlab output data*/ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; @@ -77,6 +82,10 @@ void mexFunction( 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)); + + int vecdim[1]; + vecdim[0] = 2; + infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL)); - LLT_ROF_CPU_main(Input, Output, lambdaROF, lambdaLLT, iterationsNumb, tau, dimX, dimY, dimZ); + LLT_ROF_CPU_main(Input, Output, infovec, lambdaROF, lambdaLLT, iterationsNumb, tau, epsil, dimX, dimY, dimZ); }
\ No newline at end of file diff --git a/src/Matlab/mex_compile/regularisers_CPU/NonlDiff.c b/src/Matlab/mex_compile/regularisers_CPU/NonlDiff.c index ec35b8b..57c8811 100644 --- a/src/Matlab/mex_compile/regularisers_CPU/NonlDiff.c +++ b/src/Matlab/mex_compile/regularisers_CPU/NonlDiff.c @@ -30,9 +30,11 @@ * 4. Number of iterations, for explicit scheme >= 150 is recommended [OPTIONAL parameter] * 5. tau - time-marching step for explicit scheme [OPTIONAL parameter] * 6. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight [OPTIONAL parameter] + * 7. eplsilon - tolerance constant [OPTIONAL parameter] * * Output: - * [1] Regularized image/volume + * [1] Regularized image/volume + * [2] Information vector which contains [iteration no., reached tolerance] * * This function is based on the paper by * [1] Perona, P. and Malik, J., 1990. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on pattern analysis and machine intelligence, 12(7), pp.629-639. @@ -47,7 +49,8 @@ void mexFunction( mwSize dimX, dimY, dimZ; const mwSize *dim_array; - float *Input, *Output=NULL, lambda, tau, sigma; + float *Input, *Output=NULL, lambda, tau, sigma, epsil; + float *infovec=NULL; dim_array = mxGetDimensions(prhs[0]); number_of_dims = mxGetNumberOfDimensions(prhs[0]); @@ -59,12 +62,13 @@ void mexFunction( iter_numb = 300; /* iterations number */ tau = 0.025; /* marching step parameter */ penaltytype = 1; /* Huber penalty by default */ + epsil = 1.0e-06; /*tolerance parameter*/ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } - if ((nrhs < 3) || (nrhs > 6)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant, penalty type - Huber, PM or Tukey"); - if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6)) iter_numb = (int) mxGetScalar(prhs[3]); /* iterations number */ - if ((nrhs == 5) || (nrhs == 6)) tau = (float) mxGetScalar(prhs[4]); /* marching step parameter */ - if (nrhs == 6) { + if ((nrhs < 3) || (nrhs > 7)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant, penalty type - Huber, PM or Tukey, tolerance"); + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) iter_numb = (int) mxGetScalar(prhs[3]); /* iterations number */ + if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7)) tau = (float) mxGetScalar(prhs[4]); /* marching step parameter */ + if ((nrhs == 6) || (nrhs == 7)) { char *penalty_type; penalty_type = mxArrayToString(prhs[5]); /* Huber, PM or Tukey 'Huber' is the default */ if ((strcmp(penalty_type, "Huber") != 0) && (strcmp(penalty_type, "PM") != 0) && (strcmp(penalty_type, "Tukey") != 0)) mexErrMsgTxt("Choose penalty: 'Huber', 'PM' or 'Tukey',"); @@ -73,6 +77,7 @@ void mexFunction( if (strcmp(penalty_type, "Tukey") == 0) penaltytype = 3; /* enable Tikey Biweight penalty */ mxFree(penalty_type); } + if ((nrhs == 7)) epsil = (float) mxGetScalar(prhs[6]); /* epsilon */ /*Handling Matlab output data*/ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; @@ -85,5 +90,9 @@ void mexFunction( } if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - Diffusion_CPU_main(Input, Output, lambda, sigma, iter_numb, tau, penaltytype, dimX, dimY, dimZ); + int vecdim[1]; + vecdim[0] = 2; + infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL)); + + Diffusion_CPU_main(Input, Output, infovec, lambda, sigma, iter_numb, tau, penaltytype, epsil, dimX, dimY, dimZ); }
\ No newline at end of file diff --git a/src/Matlab/mex_compile/regularisers_CPU/TGV.c b/src/Matlab/mex_compile/regularisers_CPU/TGV.c index aa4eed4..aab01b4 100644 --- a/src/Matlab/mex_compile/regularisers_CPU/TGV.c +++ b/src/Matlab/mex_compile/regularisers_CPU/TGV.c @@ -30,9 +30,12 @@ limitations under the License. * 4. parameter to control the second-order term (alpha0) * 5. Number of Chambolle-Pock (Primal-Dual) iterations * 6. Lipshitz constant (default is 12) + * 7. eplsilon - tolerance constant [OPTIONAL parameter] * * Output: - * Filtered/regulariaed image + * [1] Regularized image/volume + * [2] Information vector which contains [iteration no., reached tolerance] + * * * References: * [1] K. Bredies "Total Generalized Variation" @@ -47,7 +50,8 @@ void mexFunction( mwSize dimX, dimY, dimZ; const mwSize *dim_array; - float *Input, *Output=NULL, lambda, alpha0, alpha1, L2; + float *Input, *Output=NULL, lambda, alpha0, alpha1, L2, epsil; + float *infovec=NULL; number_of_dims = mxGetNumberOfDimensions(prhs[0]); dim_array = mxGetDimensions(prhs[0]); @@ -58,15 +62,17 @@ void mexFunction( Input = (float *) mxGetData(prhs[0]); /*noisy image/volume */ lambda = (float) mxGetScalar(prhs[1]); /* regularisation parameter */ alpha1 = 1.0f; /* parameter to control the first-order term */ - alpha0 = 0.5f; /* parameter to control the second-order term */ - iter = 300; /* Iterations number */ + alpha0 = 2.0f; /* parameter to control the second-order term */ + iter = 500; /* Iterations number */ L2 = 12.0f; /* Lipshitz constant */ + epsil = 1.0e-06; /*tolerance parameter*/ 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 */ + if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) alpha1 = (float) mxGetScalar(prhs[2]); /* parameter to control the first-order term */ + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) alpha0 = (float) mxGetScalar(prhs[3]); /* parameter to control the second-order term */ + if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7)) iter = (int) mxGetScalar(prhs[4]); /* Iterations number */ + if ((nrhs == 6) || (nrhs == 7)) L2 = (float) mxGetScalar(prhs[5]); /* Lipshitz constant */ + if (nrhs == 7) epsil = (float) mxGetScalar(prhs[6]); /* epsilon */ /*Handling Matlab output data*/ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; @@ -77,7 +83,12 @@ void mexFunction( } if (number_of_dims == 3) { Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - } + } + + int vecdim[1]; + vecdim[0] = 2; + infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL)); + /* running the function */ - TGV_main(Input, Output, lambda, alpha1, alpha0, iter, L2, dimX, dimY, dimZ); + TGV_main(Input, Output, infovec, lambda, alpha1, alpha0, iter, L2, epsil, dimX, dimY, dimZ); } diff --git a/src/Matlab/mex_compile/regularisers_GPU/Diffusion_4thO_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/Diffusion_4thO_GPU.cpp index 0cc042b..7b7a220 100644 --- a/src/Matlab/mex_compile/regularisers_GPU/Diffusion_4thO_GPU.cpp +++ b/src/Matlab/mex_compile/regularisers_GPU/Diffusion_4thO_GPU.cpp @@ -23,15 +23,17 @@ /* CUDA implementation of fourth-order diffusion scheme [1] for piecewise-smooth recovery (2D/3D case) * The minimisation is performed using explicit scheme. * - * Input Parameters: + * 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] + * 6. eplsilon - tolerance constant [OPTIONAL parameter] * * Output: * [1] Regularized image/volume + * [2] Information vector which contains [iteration no., reached tolerance] * * 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. @@ -45,7 +47,8 @@ void mexFunction( int number_of_dims, iter_numb; mwSize dimX, dimY, dimZ; const mwSize *dim_array; - float *Input, *Output=NULL, lambda, tau, sigma; + float *Input, *Output=NULL, lambda, tau, sigma, epsil; + float *infovec=NULL; dim_array = mxGetDimensions(prhs[0]); number_of_dims = mxGetNumberOfDimensions(prhs[0]); @@ -54,13 +57,15 @@ void mexFunction( 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 */ + iter_numb = 500; /* iterations number */ + tau = 0.0025; /* marching step parameter */ + epsil = 1.0e-06; /*tolerance 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 */ + if ((nrhs < 3) || (nrhs > 6)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant, tolerance"); + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6)) iter_numb = (int) mxGetScalar(prhs[3]); /* iterations number */ + if ((nrhs == 5) || (nrhs == 6)) tau = (float) mxGetScalar(prhs[4]); /* marching step parameter */ + if (nrhs == 6) epsil = (float) mxGetScalar(prhs[5]); /* tolerance parameter */ /*Handling Matlab output data*/ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; @@ -73,5 +78,9 @@ void mexFunction( } 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); + int vecdim[1]; + vecdim[0] = 2; + infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL)); + + Diffus4th_GPU_main(Input, Output, infovec, lambda, sigma, iter_numb, tau, epsil, dimX, dimY, dimZ); }
\ No newline at end of file diff --git a/src/Matlab/mex_compile/regularisers_GPU/FGP_TV_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/FGP_TV_GPU.cpp index c174e75..5ccc2b2 100644 --- a/src/Matlab/mex_compile/regularisers_GPU/FGP_TV_GPU.cpp +++ b/src/Matlab/mex_compile/regularisers_GPU/FGP_TV_GPU.cpp @@ -22,76 +22,77 @@ /* GPU (CUDA) implementation of FGP-TV [1] denoising/regularization model (2D/3D case) * - * Input Parameters: + * Input Parameters: * 1. Noisy image/volume * 2. lambdaPar - regularization parameter * 3. Number of iterations * 4. eplsilon: tolerance constant * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1) * 6. nonneg: 'nonnegativity (0 is OFF by default) - * 7. print information: 0 (off) or 1 (on) * * Output: * [1] Filtered/regularized image + * [2] Information vector which contains [iteration no., reached tolerance] * * This function is based on the Matlab's code and paper by * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems" */ + void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - int number_of_dims, iter, methTV, printswitch, nonneg; + int number_of_dims, iter, methTV, nonneg; mwSize dimX, dimY, dimZ; const mwSize *dim_array; - - float *Input, *Output=NULL, lambda, epsil; + float *Input, *infovec=NULL, *Output=NULL, lambda, epsil; number_of_dims = mxGetNumberOfDimensions(prhs[0]); dim_array = mxGetDimensions(prhs[0]); /*Handling Matlab input data*/ - if ((nrhs < 2) || (nrhs > 7)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch, print switch"); + if ((nrhs < 2) || (nrhs > 6)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter, Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch"); Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */ lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */ - iter = 300; /* default iterations number */ - epsil = 0.0001; /* default tolerance constant */ + iter = 400; /* default iterations number */ + epsil = 1.0e-06; /* default tolerance constant */ methTV = 0; /* default isotropic TV penalty */ nonneg = 0; /* default nonnegativity switch, off - 0 */ - printswitch = 0; /*default print is switched, off - 0 */ 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) || (nrhs == 7)) iter = (int) mxGetScalar(prhs[2]); /* iterations number */ - if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) epsil = (float) mxGetScalar(prhs[3]); /* tolerance constant */ - if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7)) { + if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6)) iter = (int) mxGetScalar(prhs[2]); /* iterations number */ + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6)) epsil = (float) mxGetScalar(prhs[3]); /* tolerance constant */ + if ((nrhs == 5) || (nrhs == 6)) { char *penalty_type; penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */ if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',"); if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */ mxFree(penalty_type); } - if ((nrhs == 6) || (nrhs == 7)) { + if (nrhs == 6) { nonneg = (int) mxGetScalar(prhs[5]); if ((nonneg != 0) && (nonneg != 1)) mexErrMsgTxt("Nonnegativity constraint can be enabled by choosing 1 or off - 0"); } - if (nrhs == 7) { - printswitch = (int) mxGetScalar(prhs[6]); - if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0"); - } - /*Handling Matlab output data*/ - dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; + /*Handling Matlab output data*/ + dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; + if (number_of_dims == 2) { dimZ = 1; /*2D case*/ - Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + 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)); + if (number_of_dims == 3) { + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + } + int vecdim[1]; + vecdim[0] = 2; + infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL)); /* running the function */ - TV_FGP_GPU_main(Input, Output, lambda, iter, epsil, methTV, nonneg, printswitch, dimX, dimY, dimZ); + TV_FGP_GPU_main(Input, Output, infovec, lambda, iter, epsil, methTV, nonneg, dimX, dimY, dimZ); }
\ No newline at end of file diff --git a/src/Matlab/mex_compile/regularisers_GPU/FGP_dTV_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/FGP_dTV_GPU.cpp index 3f5a4b3..6662e0b 100644 --- a/src/Matlab/mex_compile/regularisers_GPU/FGP_dTV_GPU.cpp +++ b/src/Matlab/mex_compile/regularisers_GPU/FGP_dTV_GPU.cpp @@ -33,43 +33,43 @@ * 6. eta: smoothing constant to calculate gradient of the reference [OPTIONAL] * * 7. TV-type: methodTV - 'iso' (0) or 'l1' (1) [OPTIONAL] * 8. nonneg: 'nonnegativity (0 is OFF by default) [OPTIONAL] - * 9. print information: 0 (off) or 1 (on) [OPTIONAL] * * Output: - * [1] Filtered/regularized image/volume + * [1] Filtered/regularized image + * [2] Information vector which contains [iteration no., reached tolerance] * * This function is based on the Matlab's codes and papers by * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems" * [2] M. J. Ehrhardt and M. M. Betcke, Multi-Contrast MRI Reconstruction with Structure-Guided Total Variation, SIAM Journal on Imaging Sciences 9(3), pp. 1084–1106 */ + + void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - int number_of_dims, iter, methTV, printswitch, nonneg; + int number_of_dims, iter, methTV, nonneg; mwSize dimX, dimY, dimZ; const mwSize *dim_array; - const mwSize *dim_array2; - - float *Input, *InputRef, *Output=NULL, lambda, epsil, eta; - + const mwSize *dim_array2; + float *Input, *InputRef, *Output=NULL, *infovec=NULL, lambda, epsil, eta; + number_of_dims = mxGetNumberOfDimensions(prhs[0]); dim_array = mxGetDimensions(prhs[0]); dim_array2 = mxGetDimensions(prhs[1]); /*Handling Matlab input data*/ - if ((nrhs < 3) || (nrhs > 9)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Reference(2D/3D), Regularization parameter, iterations number, tolerance, smoothing constant, penalty type ('iso' or 'l1'), nonnegativity switch, print switch"); + if ((nrhs < 3) || (nrhs > 8)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Reference(2D/3D), Regularization parameter, iterations number, tolerance, smoothing constant, penalty type ('iso' or 'l1'), nonnegativity switch"); Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */ InputRef = (float *) mxGetData(prhs[1]); /* reference image (2D/3D) */ lambda = (float) mxGetScalar(prhs[2]); /* regularization parameter */ iter = 300; /* default iterations number */ - epsil = 0.0001; /* default tolerance constant */ + epsil = 1.0e-06; /* default tolerance constant */ eta = 0.01; /* default smoothing constant */ methTV = 0; /* default isotropic TV penalty */ nonneg = 0; /* default nonnegativity switch, off - 0 */ - printswitch = 0; /*default print is switched, off - 0 */ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } @@ -79,35 +79,32 @@ void mexFunction( dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; if (number_of_dims == 2) { if ((dimX != dim_array2[0]) || (dimY != dim_array2[1])) mexErrMsgTxt("The input images have different dimensionalities");} if (number_of_dims == 3) { if ((dimX != dim_array2[0]) || (dimY != dim_array2[1]) || (dimZ != dim_array2[2])) mexErrMsgTxt("The input volumes have different dimensionalities");} - - - if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8) || (nrhs == 9)) iter = (int) mxGetScalar(prhs[3]); /* iterations number */ - if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8) || (nrhs == 9)) epsil = (float) mxGetScalar(prhs[4]); /* tolerance constant */ - if ((nrhs == 6) || (nrhs == 7) || (nrhs == 8) || (nrhs == 9)) { - eta = (float) mxGetScalar(prhs[5]); /* smoothing constant for the gradient of InputRef */ - } - if ((nrhs == 7) || (nrhs == 8) || (nrhs == 9)) { + + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8)) iter = (int) mxGetScalar(prhs[3]); /* iterations number */ + if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8)) epsil = (float) mxGetScalar(prhs[4]); /* tolerance constant */ + if ((nrhs == 6) || (nrhs == 7) || (nrhs == 8)) eta = (float) mxGetScalar(prhs[5]); /* smoothing constant for the gradient of InputRef */ + if ((nrhs == 7) || (nrhs == 8)) { char *penalty_type; penalty_type = mxArrayToString(prhs[6]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */ if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',"); if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */ mxFree(penalty_type); } - if ((nrhs == 8) || (nrhs == 9)) { + if ((nrhs == 8)) { nonneg = (int) mxGetScalar(prhs[7]); if ((nonneg != 0) && (nonneg != 1)) mexErrMsgTxt("Nonnegativity constraint can be enabled by choosing 1 or off - 0"); } - if (nrhs == 9) { - printswitch = (int) mxGetScalar(prhs[8]); - if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0"); - } - + if (number_of_dims == 2) { dimZ = 1; /*2D case*/ 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)); + int vecdim[1]; + vecdim[0] = 2; + infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL)); + /* running the function */ - dTV_FGP_GPU_main(Input, InputRef, Output, lambda, iter, epsil, eta, methTV, nonneg, printswitch, dimX, dimY, dimZ); + dTV_FGP_GPU_main(Input, InputRef, Output, infovec, lambda, iter, epsil, eta, methTV, nonneg, dimX, dimY, dimZ); }
\ No newline at end of file diff --git a/src/Matlab/mex_compile/regularisers_GPU/LLT_ROF_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/LLT_ROF_GPU.cpp index e8da4ce..f27767e 100644 --- a/src/Matlab/mex_compile/regularisers_GPU/LLT_ROF_GPU.cpp +++ b/src/Matlab/mex_compile/regularisers_GPU/LLT_ROF_GPU.cpp @@ -32,9 +32,11 @@ * 3. lambdaLLT - LLT-related regularisation parameter * 4. tau - time-marching step * 5. iter - iterations number (for both models) +* 6. eplsilon - tolerance constant [OPTIONAL parameter] * * Output: -* Filtered/regularised image +* [1] Regularized image/volume +* [2] Information vector which contains [iteration no., reached tolerance] * * References: * [1] Lysaker, M., Lundervold, A. and Tai, X.C., 2003. Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time. IEEE Transactions on image processing, 12(12), pp.1579-1590. @@ -48,14 +50,14 @@ void mexFunction( { int number_of_dims, iterationsNumb; mwSize dimX, dimY, dimZ; - const mwSize *dim_array; - - float *Input, *Output=NULL, lambdaROF, lambdaLLT, tau; + const mwSize *dim_array; + float *Input, *Output=NULL, lambdaROF, lambdaLLT, tau, epsil; + float *infovec=NULL; dim_array = mxGetDimensions(prhs[0]); number_of_dims = mxGetNumberOfDimensions(prhs[0]); - if ((nrhs < 3) || (nrhs > 5)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter (ROF), Regularisation parameter (LTT), iterations number, time-marching parameter"); + if ((nrhs < 3) || (nrhs > 6)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter (ROF), Regularisation parameter (LTT), iterations number, time-marching parameter, tolerance"); /*Handling Matlab input data*/ Input = (float *) mxGetData(prhs[0]); @@ -63,10 +65,12 @@ void mexFunction( lambdaLLT = (float) mxGetScalar(prhs[2]); /* ROF regularization parameter */ iterationsNumb = 250; tau = 0.0025; + epsil = 1.0e-06; /*tolerance parameter*/ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } - if ((nrhs == 4) || (nrhs == 5)) iterationsNumb = (int) mxGetScalar(prhs[3]); /* iterations number */ - if (nrhs == 5) tau = (float) mxGetScalar(prhs[4]); /* marching step parameter */ + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6)) iterationsNumb = (int) mxGetScalar(prhs[3]); /* iterations number */ + if ((nrhs == 5) || (nrhs == 6)) tau = (float) mxGetScalar(prhs[4]); /* marching step parameter */ + if (nrhs == 6) epsil = (float) mxGetScalar(prhs[5]); /* epsilon */ /*Handling Matlab output data*/ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; @@ -78,6 +82,10 @@ void mexFunction( 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)); + + int vecdim[1]; + vecdim[0] = 2; + infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL)); - LLT_ROF_GPU_main(Input, Output, lambdaROF, lambdaLLT, iterationsNumb, tau, dimX, dimY, dimZ); + LLT_ROF_GPU_main(Input, Output, infovec, lambdaROF, lambdaLLT, iterationsNumb, tau, epsil, dimX, dimY, dimZ); }
\ No newline at end of file diff --git a/src/Matlab/mex_compile/regularisers_GPU/NonlDiff_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/NonlDiff_GPU.cpp index 1cd0cdc..4ce983f 100644 --- a/src/Matlab/mex_compile/regularisers_GPU/NonlDiff_GPU.cpp +++ b/src/Matlab/mex_compile/regularisers_GPU/NonlDiff_GPU.cpp @@ -26,19 +26,20 @@ * The minimisation is performed using explicit scheme. * * Input Parameters: - * 1. Noisy image/volume + * 1. Noisy image/volume * 2. lambda - regularization parameter * 3. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion - * 4. Number of iterations, for explicit scheme >= 150 is recommended - * 5. tau - time-marching step for explicit scheme - * 6. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight + * 4. Number of iterations, for explicit scheme >= 150 is recommended [OPTIONAL parameter] + * 5. tau - time-marching step for explicit scheme [OPTIONAL parameter] + * 6. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight [OPTIONAL parameter] + * 7. eplsilon - tolerance constant [OPTIONAL parameter] * * Output: * [1] Regularized image/volume + * [2] Information vector which contains [iteration no., reached tolerance] * * This function is based on the paper by * [1] Perona, P. and Malik, J., 1990. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on pattern analysis and machine intelligence, 12(7), pp.629-639. - * [2] Black, M.J., Sapiro, G., Marimont, D.H. and Heeger, D., 1998. Robust anisotropic diffusion. IEEE Transactions on image processing, 7(3), pp.421-432. */ void mexFunction( @@ -48,9 +49,10 @@ void mexFunction( { int number_of_dims, iter_numb, penaltytype; mwSize dimX, dimY, dimZ; - const mwSize *dim_array; + const mwSize *dim_array; - float *Input, *Output=NULL, lambda, tau, sigma; + float *Input, *Output=NULL, lambda, tau, sigma, epsil; + float *infovec=NULL; dim_array = mxGetDimensions(prhs[0]); number_of_dims = mxGetNumberOfDimensions(prhs[0]); @@ -62,12 +64,13 @@ void mexFunction( iter_numb = 300; /* iterations number */ tau = 0.025; /* marching step parameter */ penaltytype = 1; /* Huber penalty by default */ + epsil = 1.0e-06; /*tolerance parameter*/ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } - if ((nrhs < 3) || (nrhs > 6)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant, penalty type - Huber, PM or Tukey"); - if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6)) iter_numb = (int) mxGetScalar(prhs[3]); /* iterations number */ - if ((nrhs == 5) || (nrhs == 6)) tau = (float) mxGetScalar(prhs[4]); /* marching step parameter */ - if (nrhs == 6) { + if ((nrhs < 3) || (nrhs > 7)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant, penalty type - Huber, PM or Tukey, tolerance"); + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) iter_numb = (int) mxGetScalar(prhs[3]); /* iterations number */ + if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7)) tau = (float) mxGetScalar(prhs[4]); /* marching step parameter */ + if ((nrhs == 6) || (nrhs == 7)) { char *penalty_type; penalty_type = mxArrayToString(prhs[5]); /* Huber, PM or Tukey 'Huber' is the default */ if ((strcmp(penalty_type, "Huber") != 0) && (strcmp(penalty_type, "PM") != 0) && (strcmp(penalty_type, "Tukey") != 0)) mexErrMsgTxt("Choose penalty: 'Huber', 'PM' or 'Tukey',"); @@ -76,6 +79,7 @@ void mexFunction( if (strcmp(penalty_type, "Tukey") == 0) penaltytype = 3; /* enable Tikey Biweight penalty */ mxFree(penalty_type); } + if ((nrhs == 7)) epsil = (float) mxGetScalar(prhs[6]); /* epsilon */ /*Handling Matlab output data*/ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; @@ -88,5 +92,9 @@ void mexFunction( } if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - NonlDiff_GPU_main(Input, Output, lambda, sigma, iter_numb, tau, penaltytype, dimX, dimY, dimZ); + int vecdim[1]; + vecdim[0] = 2; + infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL)); + + NonlDiff_GPU_main(Input, Output, infovec, lambda, sigma, iter_numb, tau, penaltytype, epsil, dimX, dimY, dimZ); }
\ No newline at end of file diff --git a/src/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp index bd01d55..4172323 100644 --- a/src/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp +++ b/src/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp @@ -28,9 +28,11 @@ * 2. lambda - regularization parameter [REQUIRED] * 3. Number of iterations, for explicit scheme >= 150 is recommended [REQUIRED] * 4. tau - marching step for explicit scheme, ~1 is recommended [REQUIRED] + * 5. eplsilon: tolerance constant [REQUIRED] * * Output: * [1] Regularized image/volume + * [2] Information vector which contains [iteration no., reached tolerance] * * This function is based on the paper by * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms" @@ -42,13 +44,13 @@ void mexFunction( int nrhs, const mxArray *prhs[]) { - int number_of_dims, iter_numb; + int number_of_dims, iter_numb; mwSize dimX, dimY, dimZ; - const mwSize *dim_array; + const mwSize *dim_array_i; + float *Input, *Output=NULL, lambda, tau, epsil; + float *infovec=NULL; - float *Input, *Output=NULL, lambda, tau; - - dim_array = mxGetDimensions(prhs[0]); + dim_array_i = mxGetDimensions(prhs[0]); number_of_dims = mxGetNumberOfDimensions(prhs[0]); /*Handling Matlab input data*/ @@ -56,19 +58,26 @@ void mexFunction( lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */ iter_numb = (int) mxGetScalar(prhs[2]); /* iterations number */ tau = (float) mxGetScalar(prhs[3]); /* marching step parameter */ + epsil = (float) mxGetScalar(prhs[4]); /* tolerance */ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } - if(nrhs != 4) mexErrMsgTxt("Four inputs reqired: Image(2D,3D), regularization parameter, iterations number, marching step constant"); + if(nrhs != 5) mexErrMsgTxt("Four inputs reqired: Image(2D,3D), regularization parameter, iterations number, marching step constant, tolerance"); /*Handling Matlab output data*/ - dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; + dimX = dim_array_i[0]; dimY = dim_array_i[1]; dimZ = dim_array_i[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)); + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array_i, mxSINGLE_CLASS, mxREAL)); } - if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + if (number_of_dims == 3) { + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array_i, mxSINGLE_CLASS, mxREAL)); + } + + int vecdim[1]; + vecdim[0] = 2; + infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL)); - TV_ROF_GPU_main(Input, Output, lambda, iter_numb, tau, dimX, dimY, dimZ); + TV_ROF_GPU_main(Input, Output, infovec, lambda, iter_numb, tau, epsil, dimX, dimY, dimZ); }
\ No newline at end of file diff --git a/src/Matlab/mex_compile/regularisers_GPU/SB_TV_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/SB_TV_GPU.cpp index 9d1328f..8ec95ab 100644 --- a/src/Matlab/mex_compile/regularisers_GPU/SB_TV_GPU.cpp +++ b/src/Matlab/mex_compile/regularisers_GPU/SB_TV_GPU.cpp @@ -22,17 +22,7 @@ /* CUDA mex-file for implementation of Split Bregman - TV denoising-regularisation model (2D/3D) [1] * -* Input Parameters: -* 1. Noisy image/volume -* 2. lambda - regularisation parameter -* 3. Number of iterations [OPTIONAL parameter] -* 4. eplsilon - tolerance constant [OPTIONAL parameter] -* 5. TV-type: 'iso' or 'l1' [OPTIONAL parameter] -* 6. print information: 0 (off) or 1 (on) [OPTIONAL parameter] -* -* Output: -* 1. Filtered/regularized image -* + * This function is based on the Matlab's code and paper by * [1]. Goldstein, T. and Osher, S., 2009. The split Bregman method for L1-regularized problems. SIAM journal on imaging sciences, 2(2), pp.323-343. */ @@ -42,40 +32,36 @@ void mexFunction( int nrhs, const mxArray *prhs[]) { - int number_of_dims, iter, methTV, printswitch; + int number_of_dims, iter, methTV; mwSize dimX, dimY, dimZ; const mwSize *dim_array; float *Input, *Output=NULL, lambda, epsil; + float *infovec=NULL; 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/3D), Regularization parameter, Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), print switch"); + if ((nrhs < 2) || (nrhs > 5)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter,iterations number, tolerance, penalty type ('iso' or 'l1')"); Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */ lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */ - iter = 100; /* default iterations number */ - epsil = 0.0001; /* default tolerance constant */ + iter = 200; /* default iterations number */ + epsil = 1.0e-06; /* default tolerance constant */ methTV = 0; /* default isotropic TV penalty */ - printswitch = 0; /*default print is switched, off - 0 */ 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)) iter = (int) mxGetScalar(prhs[2]); /* iterations number */ - if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6)) epsil = (float) mxGetScalar(prhs[3]); /* tolerance constant */ - if ((nrhs == 5) || (nrhs == 6)) { + if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5)) iter = (int) mxGetScalar(prhs[2]); /* iterations number */ + if ((nrhs == 4) || (nrhs == 5)) epsil = (float) mxGetScalar(prhs[3]); /* tolerance constant */ + if ((nrhs == 5)) { char *penalty_type; penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */ if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',"); if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */ mxFree(penalty_type); } - if (nrhs == 6) { - printswitch = (int) mxGetScalar(prhs[5]); - if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0"); - } /*Handling Matlab output data*/ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; @@ -86,6 +72,10 @@ void mexFunction( } if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + int vecdim[1]; + vecdim[0] = 2; + infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL)); + /* running the function */ - TV_SB_GPU_main(Input, Output, lambda, iter, epsil, methTV, printswitch, dimX, dimY, dimZ); + TV_SB_GPU_main(Input, Output, infovec, lambda, iter, epsil, methTV, dimX, dimY, dimZ); } diff --git a/src/Matlab/mex_compile/regularisers_GPU/TGV_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/TGV_GPU.cpp index 1173282..bdfd85b 100644 --- a/src/Matlab/mex_compile/regularisers_GPU/TGV_GPU.cpp +++ b/src/Matlab/mex_compile/regularisers_GPU/TGV_GPU.cpp @@ -23,6 +23,9 @@ limitations under the License. /* CUDA implementation of Primal-Dual denoising method for * Total Generilized Variation (TGV)-L2 model [1] (2D/3D) * + /* C-OMP implementation of Primal-Dual denoising method for + * Total Generilized Variation (TGV)-L2 model [1] (2D/3D) + * * Input Parameters: * 1. Noisy image/volume (2D/3D) * 2. lambda - regularisation parameter @@ -30,9 +33,12 @@ limitations under the License. * 4. parameter to control the second-order term (alpha0) * 5. Number of Chambolle-Pock (Primal-Dual) iterations * 6. Lipshitz constant (default is 12) + * 7. eplsilon - tolerance constant [OPTIONAL parameter] * * Output: - * Filtered/regularised image + * [1] Regularized image/volume + * [2] Information vector which contains [iteration no., reached tolerance] + * * * References: * [1] K. Bredies "Total Generalized Variation" @@ -46,7 +52,9 @@ void mexFunction( int number_of_dims, iter; mwSize dimX, dimY, dimZ; const mwSize *dim_array; - float *Input, *Output=NULL, lambda, alpha0, alpha1, L2; + + float *Input, *Output=NULL, lambda, alpha0, alpha1, L2, epsil; + float *infovec=NULL; number_of_dims = mxGetNumberOfDimensions(prhs[0]); dim_array = mxGetDimensions(prhs[0]); @@ -54,18 +62,20 @@ void mexFunction( /*Handling Matlab input data*/ if ((nrhs < 2) || (nrhs > 6)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D), Regularisation parameter, alpha0, alpha1, iterations number, Lipshitz Constant"); - Input = (float *) mxGetData(prhs[0]); /*noisy image (2D) */ + Input = (float *) mxGetData(prhs[0]); /*noisy image/volume */ lambda = (float) mxGetScalar(prhs[1]); /* regularisation parameter */ alpha1 = 1.0f; /* parameter to control the first-order term */ alpha0 = 2.0f; /* parameter to control the second-order term */ iter = 500; /* Iterations number */ L2 = 12.0f; /* Lipshitz constant */ + epsil = 1.0e-06; /*tolerance parameter*/ 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 */ + if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) alpha1 = (float) mxGetScalar(prhs[2]); /* parameter to control the first-order term */ + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) alpha0 = (float) mxGetScalar(prhs[3]); /* parameter to control the second-order term */ + if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7)) iter = (int) mxGetScalar(prhs[4]); /* Iterations number */ + if ((nrhs == 6) || (nrhs == 7)) L2 = (float) mxGetScalar(prhs[5]); /* Lipshitz constant */ + if (nrhs == 7) epsil = (float) mxGetScalar(prhs[6]); /* epsilon */ /*Handling Matlab output data*/ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; @@ -74,8 +84,14 @@ void mexFunction( dimZ = 1; /*2D case*/ 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)); + if (number_of_dims == 3) { + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + } + + int vecdim[1]; + vecdim[0] = 2; + infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL)); /* running the function */ - TGV_GPU_main(Input, Output, lambda, alpha1, alpha0, iter, L2, dimX, dimY, dimZ); + TGV_GPU_main(Input, Output, infovec, lambda, alpha1, alpha0, iter, L2, epsil, dimX, dimY, dimZ); } |