summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorDaniil Kazantsev <dkazanc@hotmail.com>2019-03-11 22:27:57 +0000
committerDaniil Kazantsev <dkazanc@hotmail.com>2019-03-11 22:27:57 +0000
commitb901525766f8e3473ef58a19bf3fadc178d3778c (patch)
treeb2eaa40fea9cb9155784db2d13ae30405d990056
parentcd2844b6384cb74fae5ee3f7f0fdb91be752653e (diff)
downloadregularization-b901525766f8e3473ef58a19bf3fadc178d3778c.tar.gz
regularization-b901525766f8e3473ef58a19bf3fadc178d3778c.tar.bz2
regularization-b901525766f8e3473ef58a19bf3fadc178d3778c.tar.xz
regularization-b901525766f8e3473ef58a19bf3fadc178d3778c.zip
fix for matlab wrappers
-rw-r--r--demos/demoMatlab_denoise.m36
-rw-r--r--src/Matlab/mex_compile/compileCPU_mex_Linux.m90
-rw-r--r--src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m13
-rw-r--r--src/Matlab/mex_compile/compileGPU_mex.m12
-rw-r--r--src/Matlab/mex_compile/regularisers_CPU/Diffusion_4thO.c23
-rw-r--r--src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c2
-rw-r--r--src/Matlab/mex_compile/regularisers_CPU/FGP_dTV.c42
-rw-r--r--src/Matlab/mex_compile/regularisers_CPU/LLT_ROF.c21
-rw-r--r--src/Matlab/mex_compile/regularisers_CPU/NonlDiff.c23
-rw-r--r--src/Matlab/mex_compile/regularisers_CPU/TGV.c31
-rw-r--r--src/Matlab/mex_compile/regularisers_GPU/Diffusion_4thO_GPU.cpp25
-rw-r--r--src/Matlab/mex_compile/regularisers_GPU/FGP_TV_GPU.cpp45
-rw-r--r--src/Matlab/mex_compile/regularisers_GPU/FGP_dTV_GPU.cpp47
-rw-r--r--src/Matlab/mex_compile/regularisers_GPU/LLT_ROF_GPU.cpp24
-rw-r--r--src/Matlab/mex_compile/regularisers_GPU/NonlDiff_GPU.cpp32
-rw-r--r--src/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp29
-rw-r--r--src/Matlab/mex_compile/regularisers_GPU/SB_TV_GPU.cpp38
-rw-r--r--src/Matlab/mex_compile/regularisers_GPU/TGV_GPU.cpp34
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);
}