fix for matlab wrappers

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);        
 }