Merge pull request #110 from vais-ral/tol

Tolerance-based stopping criterion, fixes for a new structure, new demos

23 files changed, 364 insertions, 283 deletions

diff --git a/src/Matlab/CMakeLists.txt b/src/Matlab/CMakeLists.txt
index b97f845..0897d7a 100755
--- a/src/Matlab/CMakeLists.txt
+++ b/src/Matlab/CMakeLists.txt
@@ -38,9 +38,9 @@ find_package(Matlab REQUIRED COMPONENTS MAIN_PROGRAM MX_LIBRARY ENG_LIBRARY )
 #set (MEX_TARGETS "CPU_TNV;CPU_ROF")

 #list(APPEND MEX_TARGETS "CPU_TNV")

 #list(APPEND MEX_TARGETS "CPU_ROF")

-

+
 file(GLOB CPU_MEX_FILES

-    "${CMAKE_SOURCE_DIR}/Matlab/mex_compile/regularisers_CPU/*.c"

+    "${CMAKE_SOURCE_DIR}/src/Matlab/mex_compile/regularisers_CPU/*.c"

     #"${CMAKE_SOURCE_DIR}/Matlab/mex_compile/regularisers_GPU/*.c"

 )

 

@@ -85,10 +85,10 @@ foreach(tgt RANGE 0 ${num})
     )

     

 target_include_directories(${current_target}

-   PUBLIC ${CMAKE_SOURCE_DIR}/Core/regularisers_CPU

-   ${CMAKE_SOURCE_DIR}/Core/regularisers_GPU

-   ${CMAKE_SOURCE_DIR}/Core/inpainters_CPU

-   ${CMAKE_SOURCE_DIR}/Core/

+   PUBLIC ${CMAKE_SOURCE_DIR}/src/Core/regularisers_CPU

+   ${CMAKE_SOURCE_DIR}/src/Core/regularisers_GPU

+   ${CMAKE_SOURCE_DIR}/src/Core/inpainters_CPU

+   ${CMAKE_SOURCE_DIR}/src/Core/

    ${MATLAB_INCLUDE_DIR})

    set_property(TARGET ${current_target} PROPERTY C_STANDARD 99)

    list(APPEND CPU_MEX_TARGETS ${current_target})

@@ -101,7 +101,7 @@ if (BUILD_CUDA)
     find_package(CUDA)

     if (CUDA_FOUND)

       file(GLOB GPU_MEX_FILES

-        "${CMAKE_SOURCE_DIR}/Matlab/mex_compile/regularisers_GPU/*.cpp"

+        "${CMAKE_SOURCE_DIR}/src/Matlab/mex_compile/regularisers_GPU/*.cpp"

       )

 

       list(LENGTH GPU_MEX_FILES num)

@@ -131,14 +131,16 @@ message("number of GPU files  " ${num})
           )

         

         target_include_directories(${current_target}

-        PUBLIC ${CMAKE_SOURCE_DIR}/Core/regularisers_CPU

-               ${CMAKE_SOURCE_DIR}/Core/regularisers_GPU

-               ${CMAKE_SOURCE_DIR}/Core/inpainters_CPU

-               ${CMAKE_SOURCE_DIR}/Core/

+        PUBLIC ${CMAKE_SOURCE_DIR}/src/Core/regularisers_CPU

+               ${CMAKE_SOURCE_DIR}/src/Core/regularisers_GPU

+               ${CMAKE_SOURCE_DIR}/src/Core/inpainters_CPU

+               ${CMAKE_SOURCE_DIR}/src/Core/

                ${MATLAB_INCLUDE_DIR})

        

         list(APPEND GPU_MEX_TARGETS ${current_target})

-        INSTALL(TARGETS ${current_target} DESTINATION "${MATLAB_DEST}")

+        INSTALL(TARGETS ${current_target} DESTINATION "${MATLAB_DEST}")
+        
+        

       endforeach()

        

       add_custom_target(MatlabWrapperGPU DEPENDS ${GPU_MEX_TARGETS})

diff --git a/src/Matlab/mex_compile/compileCPU_mex_Linux.m b/src/Matlab/mex_compile/compileCPU_mex_Linux.m
index 72a828e..19fb1a5 100644
--- a/src/Matlab/mex_compile/compileCPU_mex_Linux.m
+++ b/src/Matlab/mex_compile/compileCPU_mex_Linux.m
@@ -1,17 +1,18 @@
-% 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 = '/';
 
-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');
 copyfile(pathcopyFrom2, 'regularisers_CPU');
 
 cd regularisers_CPU
-
+%%
 Pathmove = sprintf(['..' fsep 'installed' fsep], 1i);
 
 fprintf('%s \n', '<<<<<<<<<<<Compiling CPU regularisers>>>>>>>>>>>>>');
@@ -27,15 +28,15 @@ 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);
@@ -47,11 +48,11 @@ 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"
@@ -66,16 +67,15 @@ movefile('TV_energy.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], 1i);
-cd(pathA2);
-cd demos
+pathA2 = sprintf(['..' fsep '..' fsep '..' fsep '..' fsep 'demos'], 1i);
+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..7e15233 100644
--- a/src/Matlab/mex_compile/compileGPU_mex.m
+++ b/src/Matlab/mex_compile/compileGPU_mex.m
@@ -4,17 +4,18 @@
 % In order to compile CUDA modules one needs to have nvcc-compiler
 % installed (see CUDA SDK), check it under MATLAB with !nvcc --version
 
-% In the code bellow we provide a full explicit path to nvcc compiler 
+% In the code bellow we provide a full explicit path to nvcc compiler
 % ! 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');
@@ -67,8 +68,8 @@ 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
+delete PatchSelect_core* Nonlocal_TV_core* shared.h
 fprintf('%s \n', 'All successfully compiled!');
 
-pathA2 = sprintf(['..' fsep '..' fsep], 1i);
+pathA2 = sprintf(['..' fsep '..' fsep '..' fsep '..' fsep 'demos'], 1i);
 cd(pathA2);
-cd demos
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/installed/MEXed_files_location.txt b/src/Matlab/mex_compile/installed/MEXed_files_location.txt
deleted file mode 100644
index e69de29..0000000
--- a/src/Matlab/mex_compile/installed/MEXed_files_location.txt
+++ /dev/null
diff --git a/src/Matlab/mex_compile/regularisers_CPU/Diffusion_4thO.c b/src/Matlab/mex_compile/regularisers_CPU/Diffusion_4thO.c
index 66ea9be..887a76d 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);
-}
\ No newline at end of file
+    mwSize 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);
+}
diff --git a/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c b/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c
index 642362f..251ac52 100644
--- a/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c
+++ b/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c
@@ -29,10 +29,10 @@
  * 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"
@@ -44,54 +44,55 @@ 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;
-    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, 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));
+    }    
+    mwSize vecdim[1];
+    vecdim[0] = 2;
+    infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL));
     
     /* running the function */
-    TV_FGP_CPU_main(Input, Output, lambda, iter, epsil, methTV, nonneg, printswitch, dimX, dimY, dimZ);
-}
\ No newline at end of file
+    TV_FGP_CPU_main(Input, Output, infovec, lambda, iter, epsil, methTV, nonneg, dimX, dimY, dimZ);
+}
diff --git a/src/Matlab/mex_compile/regularisers_CPU/FGP_dTV.c b/src/Matlab/mex_compile/regularisers_CPU/FGP_dTV.c
index 1a0c070..f1b70a8 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));
     
+    mwSize 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);
-}
\ No newline at end of file
+    dTV_FGP_CPU_main(Input, InputRef, Output, infovec, lambda, iter, epsil, eta, methTV, nonneg, dimX, dimY, dimZ);
+}
diff --git a/src/Matlab/mex_compile/regularisers_CPU/LLT_ROF.c b/src/Matlab/mex_compile/regularisers_CPU/LLT_ROF.c
index ab45446..5c6de9d 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));   
+    
+    mwSize 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);    
-}
\ No newline at end of file
+    LLT_ROF_CPU_main(Input, Output, infovec, lambdaROF, lambdaLLT, iterationsNumb, tau, epsil, dimX, dimY, dimZ);    
+}
diff --git a/src/Matlab/mex_compile/regularisers_CPU/NonlDiff.c b/src/Matlab/mex_compile/regularisers_CPU/NonlDiff.c
index ec35b8b..2ca17d2 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);
-}
\ No newline at end of file
+    mwSize 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);
+}
diff --git a/src/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV.c b/src/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV.c
index 014c0a0..34b9915 100644
--- a/src/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV.c
+++ b/src/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV.c
@@ -51,8 +51,8 @@ void mexFunction(
     long number_of_dims,  dimX, dimY, dimZ;
     int IterNumb, NumNeighb = 0;
     unsigned short *H_i, *H_j, *H_k;
-    const int  *dim_array;
-    const int  *dim_array2;
+    const mwSize  *dim_array;
+    const mwSize *dim_array2;
     float *A_orig, *Output=NULL, *Weights, lambda;
     
     dim_array = mxGetDimensions(prhs[0]);
diff --git a/src/Matlab/mex_compile/regularisers_CPU/PatchSelect.c b/src/Matlab/mex_compile/regularisers_CPU/PatchSelect.c
index f942539..1acab29 100644
--- a/src/Matlab/mex_compile/regularisers_CPU/PatchSelect.c
+++ b/src/Matlab/mex_compile/regularisers_CPU/PatchSelect.c
@@ -52,11 +52,11 @@ void mexFunction(
 {
     int number_of_dims,  SearchWindow, SimilarWin, NumNeighb;
     mwSize dimX, dimY, dimZ;
-    unsigned short *H_i=NULL, *H_j=NULL, *H_k=NULL;
-    const int  *dim_array;
+    const mwSize *dim_array;
+    unsigned short *H_i=NULL, *H_j=NULL, *H_k=NULL;    
     float *A, *Weights = NULL, h;
-    int dim_array2[3]; /* for 2D data */
-    int dim_array3[4]; /* for 3D data */
+    mwSize dim_array2[3]; /* for 2D data */
+    mwSize dim_array3[4]; /* for 3D data */
     
     dim_array = mxGetDimensions(prhs[0]);
     number_of_dims = mxGetNumberOfDimensions(prhs[0]);
diff --git a/src/Matlab/mex_compile/regularisers_CPU/ROF_TV.c b/src/Matlab/mex_compile/regularisers_CPU/ROF_TV.c
index 55ef2b1..ffe7b91 100644
--- a/src/Matlab/mex_compile/regularisers_CPU/ROF_TV.c
+++ b/src/Matlab/mex_compile/regularisers_CPU/ROF_TV.c
@@ -29,9 +29,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"
@@ -47,7 +49,8 @@ void mexFunction(
     int number_of_dims, iter_numb;
     mwSize dimX, dimY, dimZ;
     const mwSize *dim_array_i;
-    float *Input, *Output=NULL, lambda, tau;    
+    float *Input, *Output=NULL, lambda, tau, epsil;    
+    float *infovec=NULL;
     
     dim_array_i = mxGetDimensions(prhs[0]);
     number_of_dims = mxGetNumberOfDimensions(prhs[0]);
@@ -57,9 +60,10 @@ 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_i[0]; dimY = dim_array_i[1]; dimZ = dim_array_i[2];        
     
@@ -72,6 +76,10 @@ void mexFunction(
     if (number_of_dims == 3) {
         Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array_i, mxSINGLE_CLASS, mxREAL));
     }
+    
+    mwSize vecdim[1];
+    vecdim[0] = 2;
+    infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL));
      
-    TV_ROF_CPU_main(Input, Output, lambda, iter_numb, tau, dimX, dimY, dimZ);    
-}
\ No newline at end of file
+    TV_ROF_CPU_main(Input, Output, infovec, lambda, iter_numb, tau, epsil, dimX, dimY, dimZ);    
+}
diff --git a/src/Matlab/mex_compile/regularisers_CPU/SB_TV.c b/src/Matlab/mex_compile/regularisers_CPU/SB_TV.c
index 8636322..d1bdb3a 100644
--- a/src/Matlab/mex_compile/regularisers_CPU/SB_TV.c
+++ b/src/Matlab/mex_compile/regularisers_CPU/SB_TV.c
@@ -28,10 +28,10 @@
 * 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
+ * Output:
+ * [1] Regularized image/volume 
+ * [2] Information vector which contains [iteration no., reached tolerance]
 *
 * 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 +42,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 +82,10 @@ void mexFunction(
     }
     if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
     
+    mwSize vecdim[1];
+    vecdim[0] = 2;
+    infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL));    
+    
     /* running the function */
-    SB_TV_CPU_main(Input, Output, lambda, iter, epsil, methTV, printswitch, dimX, dimY, dimZ);
+    SB_TV_CPU_main(Input, Output, infovec, lambda, iter, epsil, methTV, dimX, dimY, dimZ);
 }
diff --git a/src/Matlab/mex_compile/regularisers_CPU/TGV.c b/src/Matlab/mex_compile/regularisers_CPU/TGV.c
index aa4eed4..9e32ae4 100644
--- a/src/Matlab/mex_compile/regularisers_CPU/TGV.c
+++ b/src/Matlab/mex_compile/regularisers_CPU/TGV.c
@@ -20,7 +20,7 @@ limitations under the License.
 #include "mex.h"
 #include "TGV_core.h"
 
-/* C-OMP implementation of Primal-Dual denoising method for 
+/* C-OMP implementation of Primal-Dual denoising method for
  * Total Generilized Variation (TGV)-L2 model [1] (2D/3D)
  *
  * Input Parameters:
@@ -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"
@@ -41,43 +44,51 @@ limitations under the License.
 void mexFunction(
         int nlhs, mxArray *plhs[],
         int nrhs, const mxArray *prhs[])
-        
+
 {
     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]);
-    
+
     /*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");
-    
+    if ((nrhs < 2) || (nrhs > 7)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D), Regularisation parameter, alpha0, alpha1, iterations number, Lipshitz Constant, tolerance");
+
     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 */      
+    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 */
-    
-    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 */
-    
+    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) || (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];
-    
+
     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));
-    }       
+    }
+
+    mwSize 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..42874ef 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);
-}
\ No newline at end of file
+    mwSize 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);
+}
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..d08e50d 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));
+    }    
+    mwSize 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);    
-}
\ No newline at end of file
+    TV_FGP_GPU_main(Input, Output, infovec, lambda, iter, epsil, methTV, nonneg, dimX, dimY, dimZ);    
+}
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..2db4556 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));
     
+    mwSize 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);
-}
\ No newline at end of file
+    dTV_FGP_GPU_main(Input, InputRef, Output, infovec, lambda, iter, epsil, eta, methTV, nonneg, dimX, dimY, dimZ);
+}
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..ff5d577 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));   
+    
+    mwSize 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);    
-}
\ No newline at end of file
+    LLT_ROF_GPU_main(Input, Output, infovec, lambdaROF, lambdaLLT, iterationsNumb, tau, epsil, dimX, dimY, dimZ);    
+}
diff --git a/src/Matlab/mex_compile/regularisers_GPU/NonlDiff_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/NonlDiff_GPU.cpp
index 1cd0cdc..43627c8 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);
-}
\ No newline at end of file
+    mwSize 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);
+}
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..d9b7e83 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));
+    }
+    
+    mwSize 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);    
-}
\ No newline at end of file
+    TV_ROF_GPU_main(Input, Output, infovec, lambda, iter_numb, tau, epsil, dimX, dimY, dimZ);    
+}
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..562dc65 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));
     
+    mwSize 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..eb1f043 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,26 +52,30 @@ 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]);
     
     /*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");
+    if ((nrhs < 2) || (nrhs > 7)) 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));
+    }    
+ 
+    mwSize 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);        
 }