all python code updated

3 files changed, 594 insertions, 525 deletions

diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py
index 09b465a..398e11c 100644
--- a/src/Python/ccpi/filters/regularisers.py
+++ b/src/Python/ccpi/filters/regularisers.py
@@ -7,7 +7,7 @@ try:
     from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU, Diff4th_GPU, TGV_GPU, LLT_ROF_GPU, PATCHSEL_GPU
     gpu_enabled = True
 except ImportError:
-    gpu_enabled = False    
+    gpu_enabled = False
 from ccpi.filters.cpu_regularisers import NDF_INPAINT_CPU, NVM_INPAINT_CPU
 
 def ROF_TV(inputData, regularisation_parameter, iterations,
@@ -15,13 +15,13 @@ def ROF_TV(inputData, regularisation_parameter, iterations,
     if device == 'cpu':
         return TV_ROF_CPU(inputData,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      time_marching_parameter,
                      tolerance_param)
     elif device == 'gpu' and gpu_enabled:
         return TV_ROF_GPU(inputData,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      time_marching_parameter,
                      tolerance_param)
     else:
@@ -35,14 +35,14 @@ def FGP_TV(inputData, regularisation_parameter,iterations,
     if device == 'cpu':
         return TV_FGP_CPU(inputData,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      tolerance_param,
                      methodTV,
                      nonneg)
     elif device == 'gpu' and gpu_enabled:
         return TV_FGP_GPU(inputData,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      tolerance_param,
                      methodTV,
                      nonneg)
@@ -56,13 +56,13 @@ def SB_TV(inputData, regularisation_parameter, iterations,
     if device == 'cpu':
         return TV_SB_CPU(inputData,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      tolerance_param,
                      methodTV)
     elif device == 'gpu' and gpu_enabled:
         return TV_SB_GPU(inputData,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      tolerance_param,
                      methodTV)
     else:
@@ -81,91 +81,115 @@ def LLT_ROF(inputData, regularisation_parameterROF, regularisation_parameterLLT,
             raise ValueError ('GPU is not available')
         raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
                          .format(device))
-def FGP_dTV(inputData, refdata, regularisation_parameter, iterations,
-                     tolerance_param, eta_const, methodTV, nonneg, printM, device='cpu'):
+def TGV(inputData, regularisation_parameter, alpha1, alpha0, iterations,
+                     LipshitzConst, tolerance_param, device='cpu'):
     if device == 'cpu':
-        return dTV_FGP_CPU(inputData,
-                     refdata,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     eta_const,
-                     methodTV,
-                     nonneg,
-                     printM)
+        return TGV_CPU(inputData,
+					regularisation_parameter,
+					alpha1,
+					alpha0,
+					iterations,
+                    LipshitzConst,
+                    tolerance_param)
     elif device == 'gpu' and gpu_enabled:
-        return dTV_FGP_GPU(inputData,
-                     refdata,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     eta_const,
-                     methodTV,
-                     nonneg,
-                     printM)
+        return TGV_GPU(inputData,
+					regularisation_parameter,
+					alpha1,
+					alpha0,
+					iterations,
+                    LipshitzConst,
+                    tolerance_param)
     else:
         if not gpu_enabled and device == 'gpu':
             raise ValueError ('GPU is not available')
         raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
                          .format(device))
-def TNV(inputData, regularisation_parameter, iterations, tolerance_param):
-        return TNV_CPU(inputData,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param)
 def NDF(inputData, regularisation_parameter, edge_parameter, iterations,
-                     time_marching_parameter, penalty_type, device='cpu'):
+                     time_marching_parameter, penalty_type, tolerance_param, device='cpu'):
     if device == 'cpu':
         return NDF_CPU(inputData,
                      regularisation_parameter,
                      edge_parameter,
-                     iterations, 
+                     iterations,
                      time_marching_parameter,
-                     penalty_type)
+                     penalty_type,
+                     tolerance_param)
     elif device == 'gpu' and gpu_enabled:
         return NDF_GPU(inputData,
                      regularisation_parameter,
                      edge_parameter,
-                     iterations, 
+                     iterations,
                      time_marching_parameter,
-                     penalty_type)
+                     penalty_type,
+                     tolerance_param)
     else:
         if not gpu_enabled and device == 'gpu':
     	    raise ValueError ('GPU is not available')
         raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
                          .format(device))
 def Diff4th(inputData, regularisation_parameter, edge_parameter, iterations,
-                     time_marching_parameter, device='cpu'):
+                     time_marching_parameter, tolerance_param, device='cpu'):
     if device == 'cpu':
         return Diff4th_CPU(inputData,
                      regularisation_parameter,
                      edge_parameter,
-                     iterations, 
-                     time_marching_parameter)
+                     iterations,
+                     time_marching_parameter,
+                     tolerance_param)
     elif device == 'gpu' and gpu_enabled:
         return Diff4th_GPU(inputData,
                      regularisation_parameter,
                      edge_parameter,
-                     iterations, 
-                     time_marching_parameter)
+                     iterations,
+                     time_marching_parameter,
+                     tolerance_param)
     else:
         if not gpu_enabled and device == 'gpu':
             raise ValueError ('GPU is not available')
         raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
                          .format(device))
-        
+def FGP_dTV(inputData, refdata, regularisation_parameter, iterations,
+                     tolerance_param, eta_const, methodTV, nonneg, device='cpu'):
+    if device == 'cpu':
+        return dTV_FGP_CPU(inputData,
+                     refdata,
+                     regularisation_parameter,
+                     iterations,
+                     tolerance_param,
+                     eta_const,
+                     methodTV,
+                     nonneg)
+    elif device == 'gpu' and gpu_enabled:
+        return dTV_FGP_GPU(inputData,
+                     refdata,
+                     regularisation_parameter,
+                     iterations,
+                     tolerance_param,
+                     eta_const,
+                     methodTV,
+                     nonneg)
+    else:
+        if not gpu_enabled and device == 'gpu':
+            raise ValueError ('GPU is not available')
+        raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
+                         .format(device))
+def TNV(inputData, regularisation_parameter, iterations, tolerance_param):
+        return TNV_CPU(inputData,
+                     regularisation_parameter,
+                     iterations,
+                     tolerance_param)
 def PatchSelect(inputData, searchwindow, patchwindow, neighbours, edge_parameter, device='cpu'):
     if device == 'cpu':
         return PATCHSEL_CPU(inputData,
                      searchwindow,
                      patchwindow,
-                     neighbours, 
+                     neighbours,
                      edge_parameter)
     elif device == 'gpu' and gpu_enabled:
         return PATCHSEL_GPU(inputData,
                      searchwindow,
                      patchwindow,
-                     neighbours, 
+                     neighbours,
                      edge_parameter)
     else:
         if not gpu_enabled and device == 'gpu':
@@ -177,36 +201,14 @@ def NLTV(inputData, H_i, H_j, H_k, Weights, regularisation_parameter, iterations
     return NLTV_CPU(inputData,
                      H_i,
                      H_j,
-                     H_k, 
+                     H_k,
                      Weights,
                      regularisation_parameter,
                      iterations)
-
-def TGV(inputData, regularisation_parameter, alpha1, alpha0, iterations,
-                     LipshitzConst, device='cpu'):
-    if device == 'cpu':
-        return TGV_CPU(inputData, 
-					regularisation_parameter, 
-					alpha1, 
-					alpha0, 
-					iterations,
-                    LipshitzConst)
-    elif device == 'gpu' and gpu_enabled:
-        return TGV_GPU(inputData, 
-					regularisation_parameter, 
-					alpha1, 
-					alpha0, 
-					iterations,
-                    LipshitzConst)
-    else:
-        if not gpu_enabled and device == 'gpu':
-            raise ValueError ('GPU is not available')
-        raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
-                         .format(device))
 def NDF_INP(inputData, maskData, regularisation_parameter, edge_parameter, iterations,
                      time_marching_parameter, penalty_type):
-        return NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter, 
+        return NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter,
         edge_parameter, iterations, time_marching_parameter, penalty_type)
-        
+
 def NVM_INP(inputData, maskData, SW_increment, iterations):
         return NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterations)
diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx
index f2276bb..add641b 100644
--- a/src/Python/src/cpu_regularisers.pyx
+++ b/src/Python/src/cpu_regularisers.pyx
@@ -22,11 +22,11 @@ cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float *infovector
 cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int dimX, int dimY, int dimZ);
 cdef extern float SB_TV_CPU_main(float *Input, float *Output, float *infovector, float mu, int iter, float epsil, int methodTV, int dimX, int dimY, int dimZ);
 cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float *infovector, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ);
-cdef extern float TGV_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY, int dimZ);
-cdef extern float Diffusion_CPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ);
-cdef extern float Diffus4th_CPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ);
+cdef extern float TGV_main(float *Input, float *Output, float *infovector, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, float epsil, int dimX, int dimY, int dimZ);
+cdef extern float Diffusion_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
+cdef extern float Diffus4th_CPU_main(float *Input, float *Output,  float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ);
+cdef extern float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int dimX, int dimY, int dimZ);
 cdef extern float TNV_CPU_main(float *Input, float *u, float lambdaPar, int maxIter, float tol, int dimX, int dimY, int dimZ);
-cdef extern float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
 cdef extern float PatchSelect_CPU_main(float *Input, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int SearchWindow, int SimilarWin, int NumNeighb, float h, int switchM);
 cdef extern float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int NumNeighb, float lambdaReg, int IterNumb);
 
@@ -43,7 +43,7 @@ def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb, marching_ste
     elif inputData.ndim == 3:
         return TV_ROF_3D(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter,tolerance_param)
 
-def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float regularisation_parameter,
                      int iterationsNumb,
                      float marching_step_parameter,
@@ -51,18 +51,18 @@ def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
-    
+
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
             np.zeros([dims[0],dims[1]], dtype='float32')
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
             np.ones([2], dtype='float32')
-            
-    # Run ROF iterations for 2D data 
+
+    # Run ROF iterations for 2D data
     TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1)
-    
+
     return (outputData,infovec)
-            
-def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+
+def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      float regularisation_parameter,
                      int iterationsNumb,
                      float marching_step_parameter,
@@ -71,13 +71,13 @@ def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     dims[2] = inputData.shape[2]
-    
+
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
             np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
             np.ones([2], dtype='float32')
-           
-    # Run ROF iterations for 3D data 
+
+    # Run ROF iterations for 3D data
     TV_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter, iterationsNumb, marching_step_parameter, tolerance_param, dims[2], dims[1], dims[0])
 
     return (outputData,infovec)
@@ -92,9 +92,9 @@ def TV_FGP_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_pa
     elif inputData.ndim == 3:
         return TV_FGP_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg)
 
-def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float regularisation_parameter,
-                     int iterationsNumb, 
+                     int iterationsNumb,
                      float tolerance_param,
                      int methodTV,
                      int nonneg):
@@ -102,35 +102,35 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
-    
+
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
             np.zeros([dims[0],dims[1]], dtype='float32')
-        
+
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
             np.ones([2], dtype='float32')
-    
+
     #/* Run FGP-TV iterations for 2D data */
-    TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter, 
-                       iterationsNumb, 
+    TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter,
+                       iterationsNumb,
                        tolerance_param,
                        methodTV,
                        nonneg,
                        dims[1],dims[0],1)
-    
+
     return (outputData,infovec)
-            
-def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+
+def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      float regularisation_parameter,
-                     int iterationsNumb, 
+                     int iterationsNumb,
                      float tolerance_param,
                      int methodTV,
                      int nonneg):
-    
+
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     dims[2] = inputData.shape[2]
-    
+
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
             np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
@@ -138,7 +138,7 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
 
     #/* Run FGP-TV iterations for 3D data */
     TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter,
-                       iterationsNumb, 
+                       iterationsNumb,
                        tolerance_param,
                        methodTV,
                        nonneg,
@@ -155,54 +155,54 @@ def TV_SB_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_par
     elif inputData.ndim == 3:
         return TV_SB_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV)
 
-def TV_SB_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+def TV_SB_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float regularisation_parameter,
-                     int iterationsNumb, 
+                     int iterationsNumb,
                      float tolerance_param,
                      int methodTV):
-                         
+
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
-    
+
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
             np.zeros([dims[0],dims[1]], dtype='float32')
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
             np.zeros([2], dtype='float32')
-                   
+
     #/* Run SB-TV iterations for 2D data */
     SB_TV_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0],
-                       regularisation_parameter, 
-                       iterationsNumb, 
+                       regularisation_parameter,
+                       iterationsNumb,
                        tolerance_param,
                        methodTV,
                        dims[1],dims[0], 1)
-    
+
     return (outputData,infovec)
-            
-def TV_SB_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+
+def TV_SB_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      float regularisation_parameter,
-                     int iterationsNumb, 
+                     int iterationsNumb,
                      float tolerance_param,
                      int methodTV):
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     dims[2] = inputData.shape[2]
-    
+
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
             np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
             np.zeros([2], dtype='float32')
-            
+
     #/* Run SB-TV iterations for 3D data */
     SB_TV_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0],
                        regularisation_parameter,
-                       iterationsNumb, 
+                       iterationsNumb,
                        tolerance_param,
                        methodTV,
                        dims[2], dims[1], dims[0])
-    return (outputData,infovec) 
+    return (outputData,infovec)
 #***************************************************************#
 #******************* ROF - LLT regularisation ******************#
 #***************************************************************#
@@ -212,288 +212,321 @@ def LLT_ROF_CPU(inputData, regularisation_parameterROF, regularisation_parameter
     elif inputData.ndim == 3:
         return LLT_ROF_3D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, tolerance_param)
 
-def LLT_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+def LLT_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float regularisation_parameterROF,
                      float regularisation_parameterLLT,
-                     int iterations, 
+                     int iterations,
                      float time_marching_parameter,
                      float tolerance_param):
-                         
+
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
-    
+
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
             np.zeros([dims[0],dims[1]], dtype='float32')
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
             np.zeros([2], dtype='float32')
-            
+
     #/* Run ROF-LLT iterations for 2D data */
-    LLT_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, 
+    LLT_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter,
                      tolerance_param,
                      dims[1],dims[0],1)
-    return (outputData,infovec) 
+    return (outputData,infovec)
 
-def LLT_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+def LLT_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      float regularisation_parameterROF,
                      float regularisation_parameterLLT,
-                     int iterations, 
+                     int iterations,
                      float time_marching_parameter,
                      float tolerance_param):
-						 
+
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     dims[2] = inputData.shape[2]
-    
+
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
             np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
             np.zeros([2], dtype='float32')
-            
+
     #/* Run ROF-LLT iterations for 3D data */
-    LLT_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameterROF, regularisation_parameterLLT, iterations, 
+    LLT_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameterROF, regularisation_parameterLLT, iterations,
                      time_marching_parameter,
-                     tolerance_param, 
+                     tolerance_param,
                      dims[2], dims[1], dims[0])
     return (outputData,infovec)
 #***************************************************************#
 #***************** Total Generalised Variation *****************#
 #***************************************************************#
-def TGV_CPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst):
+def TGV_CPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst, tolerance_param):
     if inputData.ndim == 2:
-        return TGV_2D(inputData, regularisation_parameter, alpha1, alpha0, 
-                      iterations, LipshitzConst)
+        return TGV_2D(inputData, regularisation_parameter, alpha1, alpha0,
+                      iterations, LipshitzConst, tolerance_param)
     elif inputData.ndim == 3:
-        return TGV_3D(inputData, regularisation_parameter, alpha1, alpha0, 
-                      iterations, LipshitzConst)
+        return TGV_3D(inputData, regularisation_parameter, alpha1, alpha0,
+                      iterations, LipshitzConst, tolerance_param)
 
-def TGV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+def TGV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float regularisation_parameter,
                      float alpha1,
                      float alpha0,
-                     int iterationsNumb, 
-                     float LipshitzConst):
-                         
+                     int iterationsNumb,
+                     float LipshitzConst,
+                     float tolerance_param):
+
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
-    
+
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
             np.zeros([dims[0],dims[1]], dtype='float32')
-                   
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+                np.zeros([2], dtype='float32')
+
     #/* Run TGV iterations for 2D data */
-    TGV_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, 
+    TGV_main(&inputData[0,0], &outputData[0,0],  &infovec[0],  regularisation_parameter,
                        alpha1,
                        alpha0,
-                       iterationsNumb, 
+                       iterationsNumb,
                        LipshitzConst,
+                       tolerance_param,
                        dims[1],dims[0],1)
-    return outputData
-def TGV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+    return (outputData,infovec)
+def TGV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      float regularisation_parameter,
                      float alpha1,
                      float alpha0,
-                     int iterationsNumb, 
-                     float LipshitzConst):
-                         
+                     int iterationsNumb,
+                     float LipshitzConst,
+                     float tolerance_param):
+
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     dims[2] = inputData.shape[2]
-    
+
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
             np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
-                   
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+                np.zeros([2], dtype='float32')
+
     #/* Run TGV iterations for 3D data */
-    TGV_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, 
+    TGV_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter,
                        alpha1,
                        alpha0,
-                       iterationsNumb, 
+                       iterationsNumb,
                        LipshitzConst,
+                       tolerance_param,
                        dims[2], dims[1], dims[0])
-    return outputData
-
+    return (outputData,infovec)
 
 #****************************************************************#
-#**************Directional Total-variation FGP ******************#
+#***************Nonlinear (Isotropic) Diffusion******************#
 #****************************************************************#
-#******** Directional TV Fast-Gradient-Projection (FGP)*********#
-def dTV_FGP_CPU(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg, printM):
+def NDF_CPU(inputData, regularisation_parameter, edge_parameter, iterationsNumb,time_marching_parameter, penalty_type,tolerance_param):
     if inputData.ndim == 2:
-        return dTV_FGP_2D(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg, printM)
+        return NDF_2D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, tolerance_param)
     elif inputData.ndim == 3:
-        return dTV_FGP_3D(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg, printM)
+        return NDF_3D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, tolerance_param)
 
-def dTV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-               np.ndarray[np.float32_t, ndim=2, mode="c"] refdata,
+def NDF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float regularisation_parameter,
-                     int iterationsNumb, 
-                     float tolerance_param,
-                     float eta_const,
-                     int methodTV,
-                     int nonneg,
-                     int printM):
-                         
+                     float edge_parameter,
+                     int iterationsNumb,
+                     float time_marching_parameter,
+                     int penalty_type,
+                     float tolerance_param):
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
-    
+
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
             np.zeros([dims[0],dims[1]], dtype='float32')
-                   
-    #/* Run FGP-dTV iterations for 2D data */
-    dTV_FGP_CPU_main(&inputData[0,0], &refdata[0,0], &outputData[0,0], regularisation_parameter, 
-                       iterationsNumb, 
-                       tolerance_param,
-                       eta_const,
-                       methodTV,                       
-                       nonneg,
-                       printM,
-                       dims[1], dims[0], 1)
-    
-    return outputData        
-            
-def dTV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
-               np.ndarray[np.float32_t, ndim=3, mode="c"] refdata,
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+                np.zeros([2], dtype='float32')
+
+    # Run Nonlinear Diffusion iterations for 2D data
+    Diffusion_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0],
+    regularisation_parameter, edge_parameter, iterationsNumb,
+    time_marching_parameter, penalty_type,
+    tolerance_param,
+    dims[1], dims[0], 1)
+    return (outputData,infovec)
+
+def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      float regularisation_parameter,
-                     int iterationsNumb, 
-                     float tolerance_param,
-                     float eta_const,
-                     int methodTV,
-                     int nonneg,
-                     int printM):
+                     float edge_parameter,
+                     int iterationsNumb,
+                     float time_marching_parameter,
+                     int penalty_type,
+                     float tolerance_param):
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     dims[2] = inputData.shape[2]
-    
+
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-            np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
-           
-    #/* Run FGP-dTV iterations for 3D data */
-    dTV_FGP_CPU_main(&inputData[0,0,0], &refdata[0,0,0], &outputData[0,0,0], regularisation_parameter,
-                       iterationsNumb, 
-                       tolerance_param,
-                       eta_const,
-                       methodTV,
-                       nonneg,
-                       printM,
-                       dims[2], dims[1], dims[0])
-    return outputData
-    
+            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+                np.zeros([2], dtype='float32')
+
+    # Run Nonlinear Diffusion iterations for  3D data
+    Diffusion_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0],
+    regularisation_parameter, edge_parameter, iterationsNumb,
+    time_marching_parameter, penalty_type,
+    tolerance_param,
+    dims[2], dims[1], dims[0])
+    return (outputData,infovec)
 #****************************************************************#
-#*********************Total Nuclear Variation********************#
+#*************Anisotropic Fourth-Order diffusion*****************#
 #****************************************************************#
-def TNV_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param):
+def Diff4th_CPU(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter,tolerance_param):
     if inputData.ndim == 2:
-        return 
+        return Diff4th_2D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter,tolerance_param)
     elif inputData.ndim == 3:
-        return TNV_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param)
+        return Diff4th_3D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter,tolerance_param)
+
+def Diff4th_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+                     float regularisation_parameter,
+                     float edge_parameter,
+                     int iterationsNumb,
+                     float time_marching_parameter,
+                     float tolerance_param):
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
 
-def TNV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+                np.zeros([2], dtype='float32')
+
+    # Run Anisotropic Fourth-Order diffusion for 2D data
+    Diffus4th_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0],
+    regularisation_parameter,
+    edge_parameter, iterationsNumb,
+    time_marching_parameter,
+    tolerance_param,
+    dims[1], dims[0], 1)
+    return (outputData,infovec)
+
+def Diff4th_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      float regularisation_parameter,
+                     float edge_parameter,
                      int iterationsNumb,
+                     float time_marching_parameter,
                      float tolerance_param):
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     dims[2] = inputData.shape[2]
-    
+
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
             np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-           
-    # Run TNV iterations for 3D (X,Y,Channels) data 
-    TNV_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, iterationsNumb, tolerance_param, dims[2], dims[1], dims[0])
-    return outputData
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+                    np.zeros([2], dtype='float32')
+
+    # Run Anisotropic Fourth-Order diffusion for  3D data
+    Diffus4th_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0],
+    regularisation_parameter, edge_parameter,
+    iterationsNumb, time_marching_parameter,
+    tolerance_param,
+    dims[2], dims[1], dims[0])
+    return (outputData,infovec)
 #****************************************************************#
-#***************Nonlinear (Isotropic) Diffusion******************#
+#**************Directional Total-variation FGP ******************#
 #****************************************************************#
-def NDF_CPU(inputData, regularisation_parameter, edge_parameter, iterationsNumb,time_marching_parameter, penalty_type):
+#******** Directional TV Fast-Gradient-Projection (FGP)*********#
+def dTV_FGP_CPU(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg):
     if inputData.ndim == 2:
-        return NDF_2D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
+        return dTV_FGP_2D(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg)
     elif inputData.ndim == 3:
-        return NDF_3D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
+        return dTV_FGP_3D(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg)
 
-def NDF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+def dTV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+               np.ndarray[np.float32_t, ndim=2, mode="c"] refdata,
                      float regularisation_parameter,
-                     float edge_parameter,
-                     int iterationsNumb,                     
-                     float time_marching_parameter,
-                     int penalty_type):
+                     int iterationsNumb,
+                     float tolerance_param,
+                     float eta_const,
+                     int methodTV,
+                     int nonneg):
+
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
-    
+
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')   
-    
-    # Run Nonlinear Diffusion iterations for 2D data 
-    Diffusion_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1)
-    return outputData
-            
-def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+            np.zeros([dims[0],dims[1]], dtype='float32')
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+                    np.zeros([2], dtype='float32')
+
+    #/* Run FGP-dTV iterations for 2D data */
+    dTV_FGP_CPU_main(&inputData[0,0], &refdata[0,0], &outputData[0,0], &infovec[0],
+                       regularisation_parameter,
+                       iterationsNumb,
+                       tolerance_param,
+                       eta_const,
+                       methodTV,
+                       nonneg,
+                       dims[1], dims[0], 1)
+    return (outputData,infovec)
+
+def dTV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+               np.ndarray[np.float32_t, ndim=3, mode="c"] refdata,
                      float regularisation_parameter,
-                     float edge_parameter,
-                     int iterationsNumb,                     
-                     float time_marching_parameter,
-                     int penalty_type):
+                     int iterationsNumb,
+                     float tolerance_param,
+                     float eta_const,
+                     int methodTV,
+                     int nonneg):
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     dims[2] = inputData.shape[2]
-    
+
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-    
-    # Run Nonlinear Diffusion iterations for  3D data 
-    Diffusion_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0])
+            np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+                    np.zeros([2], dtype='float32')
 
-    return outputData
+    #/* Run FGP-dTV iterations for 3D data */
+    dTV_FGP_CPU_main(&inputData[0,0,0], &refdata[0,0,0], &outputData[0,0,0], &infovec[0],
+                       regularisation_parameter,
+                       iterationsNumb,
+                       tolerance_param,
+                       eta_const,
+                       methodTV,
+                       nonneg,
+                       dims[2], dims[1], dims[0])
+    return (outputData,infovec)
 
 #****************************************************************#
-#*************Anisotropic Fourth-Order diffusion*****************#
+#*********************Total Nuclear Variation********************#
 #****************************************************************#
-def Diff4th_CPU(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter):
+def TNV_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param):
     if inputData.ndim == 2:
-        return Diff4th_2D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter)
+        return
     elif inputData.ndim == 3:
-        return Diff4th_3D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter)
+        return TNV_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param)
 
-def Diff4th_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+def TNV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      float regularisation_parameter,
-                     float edge_parameter,
-                     int iterationsNumb,                     
-                     float time_marching_parameter):
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')   
-    
-    # Run Anisotropic Fourth-Order diffusion for 2D data 
-    Diffus4th_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[1], dims[0], 1)
-    return outputData
-          
-def Diff4th_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     float edge_parameter,
                      int iterationsNumb,
-                     float time_marching_parameter):
+                     float tolerance_param):
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     dims[2] = inputData.shape[2]
-    
+
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
             np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-    
-    # Run Anisotropic Fourth-Order diffusion for  3D data 
-    Diffus4th_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[2], dims[1], dims[0])
 
+    # Run TNV iterations for 3D (X,Y,Channels) data
+    TNV_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, iterationsNumb, tolerance_param, dims[2], dims[1], dims[0])
     return outputData
-
 #****************************************************************#
 #***************Patch-based weights calculation******************#
 #****************************************************************#
@@ -511,14 +544,14 @@ def PatchSel_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
     dims[0] = neighbours
     dims[1] = inputData.shape[0]
     dims[2] = inputData.shape[1]
-    
-    
+
+
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Weights = \
             np.zeros([dims[0], dims[1],dims[2]], dtype='float32')
-    
+
     cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_i = \
             np.zeros([dims[0], dims[1],dims[2]], dtype='uint16')
-            
+
     cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_j = \
             np.zeros([dims[0], dims[1],dims[2]], dtype='uint16')
 
@@ -536,16 +569,16 @@ def PatchSel_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
     dims[1] = inputData.shape[1]
     dims[2] = inputData.shape[2]
     dims[3] = neighbours
-    
+
     cdef np.ndarray[np.float32_t, ndim=4, mode="c"] Weights = \
             np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='float32')
-    
+
     cdef np.ndarray[np.uint16_t, ndim=4, mode="c"] H_i = \
             np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='uint16')
-            
+
     cdef np.ndarray[np.uint16_t, ndim=4, mode="c"] H_j = \
             np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='uint16')
-            
+
     cdef np.ndarray[np.uint16_t, ndim=4, mode="c"] H_k = \
             np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='uint16')
 
@@ -573,10 +606,10 @@ def NLTV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     neighbours = H_i.shape[0]
-    
+
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
             np.zeros([dims[0],dims[1]], dtype='float32')
-    
+
     # Run nonlocal TV regularisation
     Nonlocal_TV_CPU_main(&inputData[0,0], &outputData[0,0], &H_i[0,0,0], &H_j[0,0,0], &H_i[0,0,0], &Weights[0,0,0], dims[1], dims[0], 0, neighbours, regularisation_parameter, iterations)
     return outputData
@@ -590,7 +623,7 @@ def NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter, edge_paramete
     elif inputData.ndim == 3:
         return NDF_INP_3D(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
 
-def NDF_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+def NDF_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData,
                      float regularisation_parameter,
                      float edge_parameter,
@@ -605,12 +638,12 @@ def NDF_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
 
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
             np.zeros([dims[0],dims[1]], dtype='float32')
-    
-    # Run Inpaiting by Diffusion iterations for 2D data 
+
+    # Run Inpaiting by Diffusion iterations for 2D data
     Diffusion_Inpaint_CPU_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1)
     return outputData
-            
-def NDF_INP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+
+def NDF_INP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      np.ndarray[np.uint8_t, ndim=3, mode="c"] maskData,
                      float regularisation_parameter,
                      float edge_parameter,
@@ -621,11 +654,11 @@ def NDF_INP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     dims[2] = inputData.shape[2]
-    
+
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
             np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-    
-    # Run Inpaiting by Diffusion iterations for 3D data 
+
+    # Run Inpaiting by Diffusion iterations for 3D data
     Diffusion_Inpaint_CPU_main(&inputData[0,0,0], &maskData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0])
 
     return outputData
@@ -636,27 +669,27 @@ def NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterationsNumb):
     if inputData.ndim == 2:
         return NVM_INP_2D(inputData, maskData, SW_increment, iterationsNumb)
     elif inputData.ndim == 3:
-        return 
+        return
 
-def NVM_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+def NVM_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData,
                      int SW_increment,
                      int iterationsNumb):
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
-    
+
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')   
-    
+            np.zeros([dims[0],dims[1]], dtype='float32')
+
     cdef np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData_upd = \
             np.zeros([dims[0],dims[1]], dtype='uint8')
-    
-    # Run Inpaiting by Nonlocal vertical marching method for 2D data 
-    NonlocalMarching_Inpaint_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], 
+
+    # Run Inpaiting by Nonlocal vertical marching method for 2D data
+    NonlocalMarching_Inpaint_main(&inputData[0,0], &maskData[0,0], &outputData[0,0],
                                   &maskData_upd[0,0],
                                   SW_increment, iterationsNumb, 1, dims[1], dims[0], 1)
-    
+
     return (outputData, maskData_upd)
 
 
@@ -669,36 +702,36 @@ def TV_ENERGY(inputData, inputData0, regularisation_parameter, typeFunctional):
     elif inputData.ndim == 3:
         return TV_ENERGY_3D(inputData, inputData0, regularisation_parameter, typeFunctional)
 
-def TV_ENERGY_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                 np.ndarray[np.float32_t, ndim=2, mode="c"] inputData0, 
+def TV_ENERGY_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+                 np.ndarray[np.float32_t, ndim=2, mode="c"] inputData0,
                      float regularisation_parameter,
                      int typeFunctional):
-    
+
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
-    
+
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] outputData = \
             np.zeros([1], dtype='float32')
-                   
-    # run function    
+
+    # run function
     TV_energy2D(&inputData[0,0], &inputData0[0,0], &outputData[0], regularisation_parameter, typeFunctional, dims[1], dims[0])
-    
+
     return outputData
-            
+
 def TV_ENERGY_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
-                 np.ndarray[np.float32_t, ndim=3, mode="c"] inputData0, 
+                 np.ndarray[np.float32_t, ndim=3, mode="c"] inputData0,
                      float regularisation_parameter,
                      int typeFunctional):
-						 
+
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     dims[2] = inputData.shape[2]
-    
+
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] outputData = \
             np.zeros([1], dtype='float32')
-           
+
     # Run function
     TV_energy3D(&inputData[0,0,0], &inputData0[0,0,0], &outputData[0], regularisation_parameter, typeFunctional, dims[2], dims[1], dims[0])
 
diff --git a/src/Python/src/gpu_regularisers.pyx b/src/Python/src/gpu_regularisers.pyx
index 6fc4135..84ee981 100644
--- a/src/Python/src/gpu_regularisers.pyx
+++ b/src/Python/src/gpu_regularisers.pyx
@@ -24,68 +24,68 @@ cdef extern int TV_ROF_GPU_main(float* Input, float* Output, float *infovector,
 cdef extern int TV_FGP_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int N, int M, int Z);
 cdef extern int TV_SB_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, int methodTV, int N, int M, int Z);
 cdef extern int LLT_ROF_GPU_main(float *Input, float *Output, float *infovector, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau,  float epsil, int N, int M, int Z);
-cdef extern int TGV_GPU_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY, int dimZ);
-cdef extern int NonlDiff_GPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int N, int M, int Z);
-cdef extern int dTV_FGP_GPU_main(float *Input, float *InputRef, float *Output, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int printM, int N, int M, int Z);
-cdef extern int Diffus4th_GPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int N, int M, int Z);
+cdef extern int TGV_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, float epsil, int dimX, int dimY, int dimZ);
+cdef extern int NonlDiff_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int N, int M, int Z);
+cdef extern int Diffus4th_GPU_main(float *Input, float *Output,  float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, float epsil, int N, int M, int Z);
+cdef extern int dTV_FGP_GPU_main(float *Input, float *InputRef, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int N, int M, int Z);
 cdef extern int PatchSelect_GPU_main(float *Input, unsigned short *H_i, unsigned short *H_j, float *Weights, int N, int M, int SearchWindow, int SimilarWin, int NumNeighb, float h);
 
 # Total-variation Rudin-Osher-Fatemi (ROF)
 def TV_ROF_GPU(inputData,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      time_marching_parameter,
                      tolerance_param):
     if inputData.ndim == 2:
-        return ROFTV2D(inputData, 
+        return ROFTV2D(inputData,
                      regularisation_parameter,
                      iterations,
                      time_marching_parameter,
                      tolerance_param)
     elif inputData.ndim == 3:
-        return ROFTV3D(inputData, 
+        return ROFTV3D(inputData,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      time_marching_parameter,
                      tolerance_param)
-                     
+
 # Total-variation Fast-Gradient-Projection (FGP)
 def TV_FGP_GPU(inputData,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      tolerance_param,
                      methodTV,
                      nonneg):
     if inputData.ndim == 2:
         return FGPTV2D(inputData,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      tolerance_param,
                      methodTV,
                      nonneg)
     elif inputData.ndim == 3:
         return FGPTV3D(inputData,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      tolerance_param,
                      methodTV,
                      nonneg)
 # Total-variation Split Bregman (SB)
 def TV_SB_GPU(inputData,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      tolerance_param,
                      methodTV):
     if inputData.ndim == 2:
         return SBTV2D(inputData,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      tolerance_param,
                      methodTV)
     elif inputData.ndim == 3:
         return SBTV3D(inputData,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      tolerance_param,
                      methodTV)
 # LLT-ROF model
@@ -95,90 +95,93 @@ def LLT_ROF_GPU(inputData, regularisation_parameterROF, regularisation_parameter
     elif inputData.ndim == 3:
         return LLT_ROF_GPU3D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, tolerance_param)
 # Total Generilised Variation (TGV)
-def TGV_GPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst):
+def TGV_GPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst, tolerance_param):
     if inputData.ndim == 2:
-        return TGV2D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst)
+        return TGV2D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst, tolerance_param)
     elif inputData.ndim == 3:
-        return TGV3D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst)
+        return TGV3D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst, tolerance_param)
 # Directional Total-variation Fast-Gradient-Projection (FGP)
 def dTV_FGP_GPU(inputData,
                      refdata,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      tolerance_param,
                      eta_const,
                      methodTV,
-                     nonneg,
-                     printM):
+                     nonneg):
     if inputData.ndim == 2:
         return FGPdTV2D(inputData,
                      refdata,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      tolerance_param,
                      eta_const,
                      methodTV,
-                     nonneg,
-                     printM)
+                     nonneg)
     elif inputData.ndim == 3:
         return FGPdTV3D(inputData,
                      refdata,
                      regularisation_parameter,
-                     iterations, 
+                     iterations,
                      tolerance_param,
                      eta_const,
                      methodTV,
-                     nonneg,
-                     printM)
+                     nonneg)
 # Nonlocal Isotropic Diffusion (NDF)
 def NDF_GPU(inputData,
                      regularisation_parameter,
                      edge_parameter,
-                     iterations, 
+                     iterations,
                      time_marching_parameter,
-                     penalty_type):
+                     penalty_type,
+                     tolerance_param):
     if inputData.ndim == 2:
         return NDF_GPU_2D(inputData,
                      regularisation_parameter,
                      edge_parameter,
-                     iterations, 
+                     iterations,
                      time_marching_parameter,
-                     penalty_type)
+                     penalty_type,
+                     tolerance_param)
     elif inputData.ndim == 3:
         return NDF_GPU_3D(inputData,
                      regularisation_parameter,
                      edge_parameter,
-                     iterations, 
+                     iterations,
                      time_marching_parameter,
-                     penalty_type)
+                     penalty_type,
+                     tolerance_param)
 # Anisotropic Fourth-Order diffusion
 def Diff4th_GPU(inputData,
                      regularisation_parameter,
                      edge_parameter,
-                     iterations, 
-                     time_marching_parameter):
+                     iterations,
+                     time_marching_parameter,
+                     tolerance_param):
     if inputData.ndim == 2:
         return Diff4th_2D(inputData,
                      regularisation_parameter,
                      edge_parameter,
-                     iterations, 
-                     time_marching_parameter)
+                     iterations,
+                     time_marching_parameter,
+                     tolerance_param)
     elif inputData.ndim == 3:
         return Diff4th_3D(inputData,
                      regularisation_parameter,
                      edge_parameter,
-                     iterations, 
-                     time_marching_parameter)
-                     
+                     iterations,
+                     time_marching_parameter,
+                     tolerance_param)
+
 #****************************************************************#
 #********************** Total-variation ROF *********************#
 #****************************************************************#
-def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float regularisation_parameter,
-                     int iterations, 
+                     int iterations,
                      float time_marching_parameter,
                      float tolerance_param):
-    
+
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
@@ -187,25 +190,25 @@ def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
 		    np.zeros([dims[0],dims[1]], dtype='float32')
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
             np.ones([2], dtype='float32')
-          
+
     # Running CUDA code here
     if (TV_ROF_GPU_main(
             &inputData[0,0], &outputData[0,0], &infovec[0],
                        regularisation_parameter,
-                       iterations, 
-                       time_marching_parameter, 
+                       iterations,
+                       time_marching_parameter,
                        tolerance_param,
                        dims[1], dims[0], 1)==0):
         return (outputData,infovec)
     else:
         raise ValueError(CUDAErrorMessage);
-    
-def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+
+def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      float regularisation_parameter,
-                     int iterations, 
+                     int iterations,
                      float time_marching_parameter,
                      float tolerance_param):
-    
+
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
@@ -215,13 +218,13 @@ def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
 		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
             np.ones([2], dtype='float32')
-          
-    # Running CUDA code here    
+
+    # Running CUDA code here
     if (TV_ROF_GPU_main(
             &inputData[0,0,0], &outputData[0,0,0], &infovec[0],
                        regularisation_parameter,
-                       iterations, 
-                       time_marching_parameter, 
+                       iterations,
+                       time_marching_parameter,
                        tolerance_param,
                        dims[2], dims[1], dims[0])==0):
         return (outputData,infovec)
@@ -231,13 +234,13 @@ def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
 #********************** Total-variation FGP *********************#
 #****************************************************************#
 #******** Total-variation Fast-Gradient-Projection (FGP)*********#
-def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float regularisation_parameter,
-                     int iterations, 
+                     int iterations,
                      float tolerance_param,
                      int methodTV,
                      int nonneg):
-    
+
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
@@ -246,10 +249,10 @@ def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
 		    np.zeros([dims[0],dims[1]], dtype='float32')
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
             np.ones([2], dtype='float32')
-          
-    # Running CUDA code here    
+
+    # Running CUDA code here
     if (TV_FGP_GPU_main(&inputData[0,0], &outputData[0,0], &infovec[0],
-                       regularisation_parameter, 
+                       regularisation_parameter,
                        iterations,
                        tolerance_param,
                        methodTV,
@@ -258,14 +261,14 @@ def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
         return (outputData,infovec)
     else:
         raise ValueError(CUDAErrorMessage);
-    
-def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+
+def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      float regularisation_parameter,
-                     int iterations, 
+                     int iterations,
                      float tolerance_param,
                      int methodTV,
                      int nonneg):
-    
+
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
@@ -276,11 +279,11 @@ def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
 		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
             np.ones([2], dtype='float32')
-          
-    # Running CUDA code here    
+
+    # Running CUDA code here
     if (TV_FGP_GPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0],
-                       regularisation_parameter, 
-                       iterations, 
+                       regularisation_parameter,
+                       iterations,
                        tolerance_param,
                        methodTV,
                        nonneg,
@@ -293,12 +296,12 @@ def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
 #********************** Total-variation SB *********************#
 #***************************************************************#
 #*************** Total-variation Split Bregman (SB)*************#
-def SBTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+def SBTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float regularisation_parameter,
-                     int iterations, 
+                     int iterations,
                      float tolerance_param,
                      int methodTV):
-    
+
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
@@ -307,11 +310,11 @@ def SBTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
 		    np.zeros([dims[0],dims[1]], dtype='float32')
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
             np.ones([2], dtype='float32')
-          
+
     # Running CUDA code here
     if (TV_SB_GPU_main(&inputData[0,0], &outputData[0,0],&infovec[0],
-                       regularisation_parameter, 
-                       iterations, 
+                       regularisation_parameter,
+                       iterations,
                        tolerance_param,
                        methodTV,
                        dims[1], dims[0], 1)==0):
@@ -319,13 +322,13 @@ def SBTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
     else:
         raise ValueError(CUDAErrorMessage);
 
-    
-def SBTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+
+def SBTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      float regularisation_parameter,
-                     int iterations, 
+                     int iterations,
                      float tolerance_param,
                      int methodTV):
-    
+
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
@@ -335,11 +338,11 @@ def SBTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
 		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
             np.ones([2], dtype='float32')
-            
-    # Running CUDA code here    
+
+    # Running CUDA code here
     if (TV_SB_GPU_main(&inputData[0,0,0], &outputData[0,0,0],&infovec[0],
-                       regularisation_parameter , 
-                       iterations, 
+                       regularisation_parameter ,
+                       iterations,
                        tolerance_param,
                        methodTV,
                        dims[2], dims[1], dims[0])==0):
@@ -352,13 +355,13 @@ def SBTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
 #************************ LLT-ROF model ************************#
 #***************************************************************#
 #************Joint LLT-ROF model for higher order **************#
-def LLT_ROF_GPU2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+def LLT_ROF_GPU2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float regularisation_parameterROF,
                      float regularisation_parameterLLT,
-                     int iterations, 
+                     int iterations,
                      float time_marching_parameter,
                      float tolerance_param):
-    
+
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
@@ -367,24 +370,24 @@ def LLT_ROF_GPU2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
 		    np.zeros([dims[0],dims[1]], dtype='float32')
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
             np.ones([2], dtype='float32')
-            
-    # Running CUDA code here    
-    if (LLT_ROF_GPU_main(&inputData[0,0], &outputData[0,0],&infovec[0],regularisation_parameterROF, regularisation_parameterLLT, iterations, 
-                         time_marching_parameter, 
+
+    # Running CUDA code here
+    if (LLT_ROF_GPU_main(&inputData[0,0], &outputData[0,0],&infovec[0],regularisation_parameterROF, regularisation_parameterLLT, iterations,
+                         time_marching_parameter,
                          tolerance_param,
                          dims[1],dims[0],1)==0):
         return (outputData,infovec)
     else:
         raise ValueError(CUDAErrorMessage);
 
-    
-def LLT_ROF_GPU3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+
+def LLT_ROF_GPU3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      float regularisation_parameterROF,
                      float regularisation_parameterLLT,
-                     int iterations, 
+                     int iterations,
                      float time_marching_parameter,
                      float tolerance_param):
-    
+
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
@@ -394,11 +397,11 @@ def LLT_ROF_GPU3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
 		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
             np.ones([2], dtype='float32')
-          
+
     # Running CUDA code here
-    if (LLT_ROF_GPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameterROF, regularisation_parameterLLT, 
-                         iterations, 
-                         time_marching_parameter, 
+    if (LLT_ROF_GPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameterROF, regularisation_parameterLLT,
+                         iterations,
+                         time_marching_parameter,
                          tolerance_param,
                          dims[2], dims[1], dims[0])==0):
         return (outputData,infovec)
@@ -409,38 +412,43 @@ def LLT_ROF_GPU3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
 #***************************************************************#
 #***************** Total Generalised Variation *****************#
 #***************************************************************#
-def TGV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+def TGV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float regularisation_parameter,
                      float alpha1,
                      float alpha0,
-                     int iterationsNumb, 
-                     float LipshitzConst):
-                         
+                     int iterationsNumb,
+                     float LipshitzConst,
+                     float tolerance_param):
+
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
-    
+
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
             np.zeros([dims[0],dims[1]], dtype='float32')
-                   
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+            np.ones([2], dtype='float32')
+
     #/* Run TGV iterations for 2D data */
-    if (TGV_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter,
+    if (TGV_GPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter,
                        alpha1,
                        alpha0,
-                       iterationsNumb, 
+                       iterationsNumb,
                        LipshitzConst,
+                       tolerance_param,
                        dims[1],dims[0], 1)==0):
-        return outputData
+        return (outputData,infovec)
     else:
         raise ValueError(CUDAErrorMessage);
 
-def TGV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+def TGV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      float regularisation_parameter,
                      float alpha1,
                      float alpha0,
-                     int iterationsNumb, 
-                     float LipshitzConst):
-    
+                     int iterationsNumb,
+                     float LipshitzConst,
+                     float tolerance_param):
+
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
@@ -448,178 +456,205 @@ def TGV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
 
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
 		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-          
-    # Running CUDA code here    
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+            np.ones([2], dtype='float32')
+
+    # Running CUDA code here
     if (TGV_GPU_main(
-            &inputData[0,0,0], &outputData[0,0,0], regularisation_parameter,
+            &inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter,
                        alpha1,
                        alpha0,
-                       iterationsNumb, 
+                       iterationsNumb,
                        LipshitzConst,
+                       tolerance_param,
                        dims[2], dims[1], dims[0])==0):
-        return outputData;
+        return (outputData,infovec)
     else:
         raise ValueError(CUDAErrorMessage);
 
-
 #****************************************************************#
-#**************Directional Total-variation FGP ******************#
+#***************Nonlinear (Isotropic) Diffusion******************#
 #****************************************************************#
-#******** Directional TV Fast-Gradient-Projection (FGP)*********#
-def FGPdTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-             np.ndarray[np.float32_t, ndim=2, mode="c"] refdata,
+def NDF_GPU_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float regularisation_parameter,
-                     int iterations, 
-                     float tolerance_param,
-                     float eta_const,
-                     int methodTV,
-                     int nonneg,
-                     int printM):
-    
+                     float edge_parameter,
+                     int iterationsNumb,
+                     float time_marching_parameter,
+                     int penalty_type,
+                     float tolerance_param):
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
 
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-		    np.zeros([dims[0],dims[1]], dtype='float32')
-          
-    # Running CUDA code here    
-    if (dTV_FGP_GPU_main(&inputData[0,0], &refdata[0,0], &outputData[0,0],
-                       regularisation_parameter, 
-                       iterations, 
-                       tolerance_param,
-                       eta_const,
-                       methodTV,
-                       nonneg,
-                       printM,
-                       dims[1], dims[0], 1)==0):
-        return outputData
+            np.zeros([dims[0],dims[1]], dtype='float32')
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+            np.ones([2], dtype='float32')
+
+    #rangecheck = penalty_type < 1 and penalty_type > 3
+    #if not rangecheck:
+#        raise ValueError('Choose penalty type as 1 for Huber, 2 - Perona-Malik, 3 - Tukey Biweight')
+
+    # Run Nonlinear Diffusion iterations for 2D data
+    # Running CUDA code here
+    if (NonlDiff_GPU_main(&inputData[0,0], &outputData[0,0],&infovec[0],
+    regularisation_parameter,
+    edge_parameter, iterationsNumb,
+    time_marching_parameter, penalty_type,
+    tolerance_param,
+    dims[1], dims[0], 1)==0):
+        return (outputData,infovec)
     else:
         raise ValueError(CUDAErrorMessage);
-
-    
-def FGPdTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-             np.ndarray[np.float32_t, ndim=3, mode="c"] refdata, 
+def NDF_GPU_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      float regularisation_parameter,
-                     int iterations, 
-                     float tolerance_param,
-                     float eta_const,
-                     int methodTV,
-                     int nonneg,
-                     int printM):
-    
+                     float edge_parameter,
+                     int iterationsNumb,
+                     float time_marching_parameter,
+                     int penalty_type,
+                     float tolerance_param):
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     dims[2] = inputData.shape[2]
 
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-          
-    # Running CUDA code here    
-    if (dTV_FGP_GPU_main(&inputData[0,0,0], &refdata[0,0,0], &outputData[0,0,0],
-                       regularisation_parameter , 
-                       iterations, 
-                       tolerance_param,
-                       eta_const,
-                       methodTV,
-                       nonneg,
-                       printM,
-                       dims[2], dims[1], dims[0])==0):
-        return outputData;
+            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+            np.ones([2], dtype='float32')
+
+    # Run Nonlinear Diffusion iterations for  3D data
+    # Running CUDA code here
+    if (NonlDiff_GPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0],
+    regularisation_parameter, edge_parameter,
+    iterationsNumb, time_marching_parameter,
+    penalty_type,
+    tolerance_param,
+    dims[2], dims[1], dims[0])==0):
+        return (outputData,infovec)
     else:
         raise ValueError(CUDAErrorMessage);
 
-
 #****************************************************************#
-#***************Nonlinear (Isotropic) Diffusion******************#
+#************Anisotropic Fourth-Order diffusion******************#
 #****************************************************************#
-def NDF_GPU_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+def Diff4th_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float regularisation_parameter,
                      float edge_parameter,
-                     int iterationsNumb,                     
+                     int iterationsNumb,
                      float time_marching_parameter,
-                     int penalty_type):
+                     float tolerance_param):
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
-    
+
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
             np.zeros([dims[0],dims[1]], dtype='float32')
-    
-    #rangecheck = penalty_type < 1 and penalty_type > 3
-    #if not rangecheck:
-#        raise ValueError('Choose penalty type as 1 for Huber, 2 - Perona-Malik, 3 - Tukey Biweight')
-    
-    # Run Nonlinear Diffusion iterations for 2D data 
-    # Running CUDA code here  
-    if (NonlDiff_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1)==0):
-        return outputData;
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+            np.ones([2], dtype='float32')
+
+    # Run Anisotropic Fourth-Order diffusion for 2D data
+    # Running CUDA code here
+    if (Diffus4th_GPU_main(&inputData[0,0], &outputData[0,0], &infovec[0],
+    regularisation_parameter, edge_parameter, iterationsNumb,
+    time_marching_parameter,
+    tolerance_param,
+    dims[1], dims[0], 1)==0):
+        return (outputData,infovec)
     else:
         raise ValueError(CUDAErrorMessage);
 
-            
-def NDF_GPU_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+def Diff4th_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      float regularisation_parameter,
                      float edge_parameter,
-                     int iterationsNumb,                     
+                     int iterationsNumb,
                      float time_marching_parameter,
-                     int penalty_type):
+                     float tolerance_param):
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     dims[2] = inputData.shape[2]
-    
+
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')    
-       
-    # Run Nonlinear Diffusion iterations for  3D data 
-    # Running CUDA code here  
-    if (NonlDiff_GPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0])==0):
-        return outputData;
+            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+            np.ones([2], dtype='float32')
+
+    # Run Anisotropic Fourth-Order diffusion for  3D data
+    # Running CUDA code here
+    if (Diffus4th_GPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0],
+    regularisation_parameter, edge_parameter,
+    iterationsNumb, time_marching_parameter,
+    tolerance_param,
+    dims[2], dims[1], dims[0])==0):
+        return (outputData,infovec)
     else:
         raise ValueError(CUDAErrorMessage);
-
 #****************************************************************#
-#************Anisotropic Fourth-Order diffusion******************#
+#**************Directional Total-variation FGP ******************#
 #****************************************************************#
-def Diff4th_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+#******** Directional TV Fast-Gradient-Projection (FGP)*********#
+def FGPdTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+             np.ndarray[np.float32_t, ndim=2, mode="c"] refdata,
                      float regularisation_parameter,
-                     float edge_parameter,
-                     int iterationsNumb,
-                     float time_marching_parameter):
+                     int iterations,
+                     float tolerance_param,
+                     float eta_const,
+                     int methodTV,
+                     int nonneg):
+
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
-    
+
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')
-    
-    # Run Anisotropic Fourth-Order diffusion for 2D data 
-    # Running CUDA code here  
-    if (Diffus4th_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[1], dims[0], 1)==0):
-        return outputData
+		    np.zeros([dims[0],dims[1]], dtype='float32')
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+            np.ones([2], dtype='float32')
+
+    # Running CUDA code here
+    if (dTV_FGP_GPU_main(&inputData[0,0], &refdata[0,0], &outputData[0,0], &infovec[0],
+                       regularisation_parameter,
+                       iterations,
+                       tolerance_param,
+                       eta_const,
+                       methodTV,
+                       nonneg,
+                       dims[1], dims[0], 1)==0):
+        return (outputData,infovec)
     else:
         raise ValueError(CUDAErrorMessage);
 
-            
-def Diff4th_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+
+def FGPdTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+             np.ndarray[np.float32_t, ndim=3, mode="c"] refdata,
                      float regularisation_parameter,
-                     float edge_parameter,
-                     int iterationsNumb,
-                     float time_marching_parameter):
+                     int iterations,
+                     float tolerance_param,
+                     float eta_const,
+                     int methodTV,
+                     int nonneg):
+
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     dims[2] = inputData.shape[2]
-    
+
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')    
-       
-    # Run Anisotropic Fourth-Order diffusion for  3D data 
-    # Running CUDA code here  
-    if (Diffus4th_GPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[2], dims[1], dims[0])==0):
-        return outputData;
+		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+            np.ones([2], dtype='float32')
+
+    # Running CUDA code here
+    if (dTV_FGP_GPU_main(&inputData[0,0,0], &refdata[0,0,0], &outputData[0,0,0], &infovec[0],
+                       regularisation_parameter ,
+                       iterations,
+                       tolerance_param,
+                       eta_const,
+                       methodTV,
+                       nonneg,
+                       dims[2], dims[1], dims[0])==0):
+        return (outputData,infovec)
     else:
         raise ValueError(CUDAErrorMessage);
 
@@ -639,14 +674,14 @@ def PatchSel_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
     cdef long dims[3]
     dims[0] = neighbours
     dims[1] = inputData.shape[0]
-    dims[2] = inputData.shape[1]    
-    
+    dims[2] = inputData.shape[1]
+
     cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Weights = \
             np.zeros([dims[0], dims[1],dims[2]], dtype='float32')
-    
+
     cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_i = \
             np.zeros([dims[0], dims[1],dims[2]], dtype='uint16')
-            
+
     cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_j = \
             np.zeros([dims[0], dims[1],dims[2]], dtype='uint16')
 
@@ -655,4 +690,3 @@ def PatchSel_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
         return H_i, H_j, Weights;
     else:
         raise ValueError(CUDAErrorMessage);
-