Merge pull request #138 from vais-ral/pdtv_gpu_fix

fixes gpu issues with pdtv

18 files changed, 92 insertions, 125 deletions

diff --git a/Readme.md b/Readme.md
index 4d8e7ca..09bf04f 100644
--- a/Readme.md
+++ b/Readme.md
@@ -103,7 +103,8 @@ conda install ccpi-regulariser -c ccpi -c conda-forge
 ```
 	export CIL_VERSION=19.10 (Unix) / set CIL_VERSION=19.10 (Windows)
 	conda build recipe/ --numpy 1.15 --python 3.7  
-	conda install ccpi-regulariser=${CIL_VERSION} --use-local --force-reinstall
+	conda install ccpi-regulariser=${CIL_VERSION} --use-local --force-reinstall # if this don't work (probably you're on Python 3*)
+  conda install -c file://${CONDA_PREFIX}/conda-bld/ ccpi-regulariser=${CIL_VERSION} --force-reinstall # try this one
 	cd demos/
 	python demo_cpu_regularisers.py # to run CPU demo
 	python demo_gpu_regularisers.py # to run GPU demo
diff --git a/demos/demo_cpu_regularisers.py b/demos/demo_cpu_regularisers.py
index 7d66b7f..50a5065 100644
--- a/demos/demo_cpu_regularisers.py
+++ b/demos/demo_cpu_regularisers.py
@@ -176,11 +176,10 @@ pars = {'algorithm' : PD_TV, \
         'input' : u0,\
         'regularisation_parameter':0.02, \
         'number_of_iterations' :1500 ,\
-        'tolerance_constant':1e-08,\
+        'tolerance_constant':1e-06,\
         'methodTV': 0 ,\
         'nonneg': 1,
-        'lipschitz_const' : 8,
-        'tau' : 0.0025}
+        'lipschitz_const' : 8}
         
 print ("#############PD TV CPU####################")
 start_time = timeit.default_timer()
@@ -190,8 +189,7 @@ start_time = timeit.default_timer()
               pars['tolerance_constant'], 
               pars['methodTV'],
               pars['nonneg'],
-              pars['lipschitz_const'],
-              pars['tau'],'cpu')
+              pars['lipschitz_const'],'cpu')
 
 Qtools = QualityTools(Im, pd_cpu)
 pars['rmse'] = Qtools.rmse()
diff --git a/demos/demo_cpu_regularisers3D.py b/demos/demo_cpu_regularisers3D.py
index cfdd2d4..0e7e9be 100644
--- a/demos/demo_cpu_regularisers3D.py
+++ b/demos/demo_cpu_regularisers3D.py
@@ -188,8 +188,7 @@ pars = {'algorithm' : PD_TV, \
         'tolerance_constant':1e-06,\
         'methodTV': 0 ,\
         'nonneg': 0,
-        'lipschitz_const' : 8,
-        'tau' : 0.0025}
+        'lipschitz_const' : 8}
 
 print ("#############FGP TV GPU####################")
 start_time = timeit.default_timer()
@@ -199,8 +198,7 @@ start_time = timeit.default_timer()
               pars['tolerance_constant'], 
               pars['methodTV'],
               pars['nonneg'],
-              pars['lipschitz_const'],
-              pars['tau'],'cpu')
+              pars['lipschitz_const'],'cpu')
              
 Qtools = QualityTools(idealVol, pd_cpu3D)
 pars['rmse'] = Qtools.rmse()
diff --git a/demos/demo_cpu_vs_gpu_regularisers.py b/demos/demo_cpu_vs_gpu_regularisers.py
index 015dfc6..a34bc19 100644
--- a/demos/demo_cpu_vs_gpu_regularisers.py
+++ b/demos/demo_cpu_vs_gpu_regularisers.py
@@ -241,8 +241,7 @@ pars = {'algorithm' : PD_TV, \
         'tolerance_constant':0.0,\
         'methodTV': 0 ,\
         'nonneg': 0,
-        'lipschitz_const' : 8,
-        'tau' : 0.0025}
+        'lipschitz_const' : 8}
         
 print ("#############PD TV CPU####################")
 start_time = timeit.default_timer()
@@ -252,8 +251,7 @@ start_time = timeit.default_timer()
               pars['tolerance_constant'], 
               pars['methodTV'],
               pars['nonneg'],
-              pars['lipschitz_const'],
-              pars['tau'],'cpu')
+              pars['lipschitz_const'],'cpu')
 
 Qtools = QualityTools(Im, pd_cpu)
 pars['rmse'] = Qtools.rmse()
@@ -279,8 +277,7 @@ pars = {'algorithm' : PD_TV, \
         'tolerance_constant':0.0,\
         'methodTV': 0 ,\
         'nonneg': 0,
-        'lipschitz_const' : 8,
-        'tau' : 0.0025}
+        'lipschitz_const' : 8}
         
 print ("#############PD TV CPU####################")
 start_time = timeit.default_timer()
@@ -290,8 +287,7 @@ start_time = timeit.default_timer()
               pars['tolerance_constant'], 
               pars['methodTV'],
               pars['nonneg'],
-              pars['lipschitz_const'],
-              pars['tau'],'gpu')
+              pars['lipschitz_const'],'gpu')
 
 Qtools = QualityTools(Im, pd_gpu)
 pars['rmse'] = Qtools.rmse()
diff --git a/demos/demo_gpu_regularisers.py b/demos/demo_gpu_regularisers.py
index 5131847..c6114db 100644
--- a/demos/demo_gpu_regularisers.py
+++ b/demos/demo_gpu_regularisers.py
@@ -176,8 +176,7 @@ pars = {'algorithm' : PD_TV, \
         'tolerance_constant':1e-06,\
         'methodTV': 0 ,\
         'nonneg': 1,
-        'lipschitz_const' : 8,
-        'tau' : 0.0025}
+        'lipschitz_const' : 8}
         
 print ("#############PD TV CPU####################")
 start_time = timeit.default_timer()
@@ -187,8 +186,7 @@ start_time = timeit.default_timer()
               pars['tolerance_constant'], 
               pars['methodTV'],
               pars['nonneg'],
-              pars['lipschitz_const'],
-              pars['tau'],'gpu')
+              pars['lipschitz_const'],'gpu')
 
 Qtools = QualityTools(Im, pd_gpu)
 pars['rmse'] = Qtools.rmse()
diff --git a/demos/demo_gpu_regularisers3D.py b/demos/demo_gpu_regularisers3D.py
index 2c25d01..18d97e5 100644
--- a/demos/demo_gpu_regularisers3D.py
+++ b/demos/demo_gpu_regularisers3D.py
@@ -192,8 +192,7 @@ pars = {'algorithm' : PD_TV, \
         'tolerance_constant':1e-06,\
         'methodTV': 0 ,\
         'nonneg': 0,
-        'lipschitz_const' : 8,
-        'tau' : 0.0025}
+        'lipschitz_const' : 8}
 
 print ("#############PD TV GPU####################")
 start_time = timeit.default_timer()
@@ -203,8 +202,7 @@ start_time = timeit.default_timer()
               pars['tolerance_constant'], 
               pars['methodTV'],
               pars['nonneg'],
-              pars['lipschitz_const'],
-              pars['tau'],'gpu')
+              pars['lipschitz_const'],'gpu')
 
 Qtools = QualityTools(idealVol, pd_gpu3D)
 pars['rmse'] = Qtools.rmse()
diff --git a/recipe/conda_build_config.yaml b/recipe/conda_build_config.yaml
index 2efefe2..55ceeaf 100644
--- a/recipe/conda_build_config.yaml
+++ b/recipe/conda_build_config.yaml
@@ -8,4 +8,3 @@ numpy:
   - 1.12
   - 1.14
   - 1.15
-  - 1.16
diff --git a/src/Core/regularisers_CPU/PD_TV_core.c b/src/Core/regularisers_CPU/PD_TV_core.c
index 534091b..c1b21e7 100644
--- a/src/Core/regularisers_CPU/PD_TV_core.c
+++ b/src/Core/regularisers_CPU/PD_TV_core.c
@@ -29,7 +29,6 @@
  * 5. lipschitz_const: convergence related parameter
  * 6. TV-type: methodTV - 'iso' (0) or 'l1' (1)
  * 7. nonneg: 'nonnegativity (0 is OFF by default, 1 is ON)
- * 8. tau: time marching parameter
 
  * Output:
  * [1] TV - Filtered/regularized image/volume
@@ -38,16 +37,17 @@
  * [1] Antonin Chambolle, Thomas Pock. "A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging", 2010
  */
 
-float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float lipschitz_const, int methodTV, int nonneg, float tau, int dimX, int dimY, int dimZ)
+float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float lipschitz_const, int methodTV, int nonneg, int dimX, int dimY, int dimZ)
 {
     int ll;
     long j, DimTotal;
-    float re, re1, sigma, theta, lt;
+    float re, re1, sigma, theta, lt, tau;
     re = 0.0f; re1 = 0.0f;
     int count = 0;
 
     //tau = 1.0/powf(lipschitz_const,0.5);
     //sigma = 1.0/powf(lipschitz_const,0.5);
+    tau = lambdaPar*0.1f;
     sigma = 1.0/(lipschitz_const*tau);
     theta = 1.0f;
     lt = tau/lambdaPar;
diff --git a/src/Core/regularisers_CPU/PD_TV_core.h b/src/Core/regularisers_CPU/PD_TV_core.h
index 97edc05..294e75c 100644
--- a/src/Core/regularisers_CPU/PD_TV_core.h
+++ b/src/Core/regularisers_CPU/PD_TV_core.h
@@ -47,7 +47,7 @@ limitations under the License.
 #ifdef __cplusplus
 extern "C" {
 #endif
-float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float lipschitz_const, int methodTV, int nonneg, float tau, int dimX, int dimY, int dimZ);
+float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float lipschitz_const, int methodTV, int nonneg, int dimX, int dimY, int dimZ);
 
 CCPI_EXPORT float DualP2D(float *U, float *P1, float *P2, long dimX, long dimY, float sigma);
 CCPI_EXPORT float DivProj2D(float *U, float *Input, float *P1, float *P2, long dimX, long dimY, float lt, float tau);
diff --git a/src/Core/regularisers_GPU/TV_PD_GPU_core.cu b/src/Core/regularisers_GPU/TV_PD_GPU_core.cu
index e57020a..01e0ab5 100644
--- a/src/Core/regularisers_GPU/TV_PD_GPU_core.cu
+++ b/src/Core/regularisers_GPU/TV_PD_GPU_core.cu
@@ -33,7 +33,6 @@ limitations under the License.
  * 5. lipschitz_const: convergence related parameter
  * 6. TV-type: methodTV - 'iso' (0) or 'l1' (1)
  * 7. nonneg: 'nonnegativity (0 is OFF by default, 1 is ON)
- * 8. tau: time marching parameter
 
  * Output:
  * [1] TV - Filtered/regularized image/volume
@@ -42,8 +41,8 @@ limitations under the License.
  * [1] Antonin Chambolle, Thomas Pock. "A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging", 2010
  */
 
-#define BLKXSIZE2D 8
-#define BLKYSIZE2D 8
+#define BLKXSIZE2D 16
+#define BLKYSIZE2D 16
 
 #define BLKXSIZE 8
 #define BLKYSIZE 8
@@ -322,12 +321,9 @@ __global__ void PDResidCalc3D_kernel(float *Input1, float *Input2, float* Output
         Output[index] = Input1[index] - Input2[index];
     }
 }
-
-
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
-
 ////////////MAIN HOST FUNCTION ///////////////
-extern "C" int TV_PD_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, float lipschitz_const, int methodTV, int nonneg, float tau, int dimX, int dimY, int dimZ)
+extern "C" int TV_PD_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, float lipschitz_const, int methodTV, int nonneg, int dimX, int dimY, int dimZ)
 {
    int deviceCount = -1; // number of devices
    cudaGetDeviceCount(&deviceCount);
@@ -336,9 +332,10 @@ extern "C" int TV_PD_GPU_main(float *Input, float *Output, float *infovector, fl
        return -1;
    }
    int count = 0, i;
-   float re, sigma, theta, lt;
+   float re, sigma, theta, lt, tau;
    re = 0.0f;
 
+   tau = lambdaPar*0.1f;
    sigma = 1.0/(lipschitz_const*tau);
    theta = 1.0f;
    lt = tau/lambdaPar;
diff --git a/src/Core/regularisers_GPU/TV_PD_GPU_core.h b/src/Core/regularisers_GPU/TV_PD_GPU_core.h
index 2b123d9..48e353e 100644
--- a/src/Core/regularisers_GPU/TV_PD_GPU_core.h
+++ b/src/Core/regularisers_GPU/TV_PD_GPU_core.h
@@ -4,6 +4,6 @@
 #include "CCPiDefines.h"
 #include <memory.h>
 
-extern "C" CCPI_EXPORT int TV_PD_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, float lipschitz_const, int methodTV, int nonneg, float tau, int dimX, int dimY, int dimZ);
+extern "C" CCPI_EXPORT int TV_PD_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, float lipschitz_const, int methodTV, int nonneg, int dimX, int dimY, int dimZ);
 
 #endif
diff --git a/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c b/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c
index e5ab1e4..f8f5272 100644
--- a/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c
+++ b/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c
@@ -30,8 +30,7 @@
  * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
  * 6. nonneg: 'nonnegativity (0 is OFF by default, 1 is ON)
  * 7. lipschitz_const: convergence related parameter
- * 8. tau: convergence related parameter
- 
+
  * Output:
  * [1] TV - Filtered/regularized image/volume
  * [2] Information vector which contains [iteration no., reached tolerance]
@@ -41,19 +40,19 @@
 void mexFunction(
         int nlhs, mxArray *plhs[],
         int nrhs, const mxArray *prhs[])
-        
+
 {
     int number_of_dims, iter, methTV, nonneg;
     mwSize dimX, dimY, dimZ;
     const mwSize *dim_array;
-    float *Input, *infovec=NULL, *Output=NULL, lambda, epsil, lipschitz_const, tau;
-    
+    float *Input, *infovec=NULL, *Output=NULL, lambda, epsil, lipschitz_const;
+
     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, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch, lipschitz_const");
-    
+    if ((nrhs < 2) || (nrhs > 6)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch, lipschitz_const");
+
     Input  = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
     lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
     iter = 500; /* default iterations number */
@@ -61,40 +60,39 @@ void mexFunction(
     methTV = 0;  /* default isotropic TV penalty */
     nonneg = 0; /* default nonnegativity switch, off - 0 */
     lipschitz_const = 8.0; /* lipschitz_const */
-    tau = 0.0025; /* tau convergence const */
-    
+
+
     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) || (nrhs == 8))  iter = (int) mxGetScalar(prhs[2]); /* iterations number */
-    if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8))  epsil =  (float) mxGetScalar(prhs[3]); /* tolerance constant */
-    if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8))  {
+
+    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))  {
         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) || (nrhs == 8))  {
+    if ((nrhs == 6) || (nrhs == 7))  {
         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) || (nrhs == 8)) lipschitz_const = (float) mxGetScalar(prhs[6]);
-    if (nrhs == 8) tau = (float) mxGetScalar(prhs[7]);   
-    
+    if (nrhs == 7) lipschitz_const = (float) mxGetScalar(prhs[6]);
+
     /*Handling Matlab output data*/
-    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];    
-       
+    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) {    
+    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 */    
-    PDTV_CPU_main(Input, Output, infovec, lambda, iter,  epsil, lipschitz_const, methTV, nonneg, tau, dimX, dimY, dimZ);
+
+    /* running the function */
+    PDTV_CPU_main(Input, Output, infovec, lambda, iter,  epsil, lipschitz_const, methTV, nonneg, dimX, dimY, dimZ);
 }
diff --git a/src/Matlab/mex_compile/regularisers_GPU/PD_TV_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/PD_TV_GPU.cpp
index e853dd3..2c037a5 100644
--- a/src/Matlab/mex_compile/regularisers_GPU/PD_TV_GPU.cpp
+++ b/src/Matlab/mex_compile/regularisers_GPU/PD_TV_GPU.cpp
@@ -30,8 +30,7 @@
  * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
  * 6. nonneg: 'nonnegativity (0 is OFF by default, 1 is ON)
  * 7. lipschitz_const: convergence related parameter
- * 8. tau: convergence related parameter
- 
+
  * Output:
  * [1] TV - Filtered/regularized image/volume
  * [2] Information vector which contains [iteration no., reached tolerance]
@@ -42,19 +41,19 @@
 void mexFunction(
         int nlhs, mxArray *plhs[],
         int nrhs, const mxArray *prhs[])
-        
+
 {
     int number_of_dims, iter, methTV, nonneg;
     mwSize dimX, dimY, dimZ;
     const mwSize *dim_array;
-    float *Input, *infovec=NULL, *Output=NULL, lambda, epsil, lipschitz_const, tau;
-    
+    float *Input, *infovec=NULL, *Output=NULL, lambda, epsil, lipschitz_const;
+
     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, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch, lipschitz_const");
-    
+    if ((nrhs < 2) || (nrhs > 6)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch, lipschitz_const");
+
     Input  = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
     lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
     iter = 500; /* default iterations number */
@@ -62,40 +61,38 @@ void mexFunction(
     methTV = 0;  /* default isotropic TV penalty */
     nonneg = 0; /* default nonnegativity switch, off - 0 */
     lipschitz_const = 8.0; /* lipschitz_const */
-    tau = 0.0025; /* tau convergence const */
-    
+
     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) || (nrhs == 8))  iter = (int) mxGetScalar(prhs[2]); /* iterations number */
-    if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8))  epsil =  (float) mxGetScalar(prhs[3]); /* tolerance constant */
-    if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8))  {
+
+    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))  {
         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) || (nrhs == 8))  {
+    if ((nrhs == 6) || (nrhs == 7))  {
         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) || (nrhs == 8)) lipschitz_const = (float) mxGetScalar(prhs[6]);
-    if (nrhs == 8) tau = (float) mxGetScalar(prhs[7]);   
-    
+    if (nrhs == 7) lipschitz_const = (float) mxGetScalar(prhs[6]);
+
     /*Handling Matlab output data*/
-    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];    
-       
+    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) {    
+    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_PD_GPU_main(Input, Output, infovec, lambda, iter,  epsil, lipschitz_const, methTV, nonneg, tau, dimX, dimY, dimZ);
+
+    /* running the function */
+    TV_PD_GPU_main(Input, Output, infovec, lambda, iter,  epsil, lipschitz_const, methTV, nonneg, dimX, dimY, dimZ);
 }
diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py
index 5f4001a..e3a984e 100644
--- a/src/Python/ccpi/filters/regularisers.py
+++ b/src/Python/ccpi/filters/regularisers.py
@@ -53,7 +53,7 @@ def FGP_TV(inputData, regularisation_parameter,iterations,
                          .format(device))
 
 def PD_TV(inputData, regularisation_parameter, iterations,
-                     tolerance_param, methodTV, nonneg, lipschitz_const, tau, device='cpu'):
+                     tolerance_param, methodTV, nonneg, lipschitz_const, device='cpu'):
     if device == 'cpu':
         return TV_PD_CPU(inputData,
                      regularisation_parameter,
@@ -61,8 +61,7 @@ def PD_TV(inputData, regularisation_parameter, iterations,
                      tolerance_param,
                      methodTV,
                      nonneg,
-                     lipschitz_const,
-                     tau)
+                     lipschitz_const)
     elif device == 'gpu' and gpu_enabled:
         return TV_PD_GPU(inputData,
                      regularisation_parameter,
@@ -70,8 +69,7 @@ def PD_TV(inputData, regularisation_parameter, iterations,
                      tolerance_param,
                      methodTV,
                      nonneg,
-                     lipschitz_const,
-                     tau)
+                     lipschitz_const)
     else:
         if not gpu_enabled and device == 'gpu':
             raise ValueError ('GPU is not available')
diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx
index 8de6aea..08e247c 100644
--- a/src/Python/src/cpu_regularisers.pyx
+++ b/src/Python/src/cpu_regularisers.pyx
@@ -20,7 +20,7 @@ cimport numpy as np
 
 cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float *infovector, float *lambdaPar, int lambda_is_arr, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ);
 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 PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float lipschitz_const, int methodTV, int nonneg, float tau, int dimX, int dimY, int dimZ);
+cdef extern float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float lipschitz_const, 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 *infovector, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, float epsil, int dimX, int dimY, int dimZ);
@@ -159,11 +159,11 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
 #****************************************************************#
 #****************** Total-variation Primal-dual *****************#
 #****************************************************************#
-def TV_PD_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const, tau):
+def TV_PD_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const):
     if inputData.ndim == 2:
-        return TV_PD_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const, tau)
+        return TV_PD_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const)
     elif inputData.ndim == 3:
-        return TV_PD_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const, tau)
+        return TV_PD_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const)
 
 def TV_PD_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float regularisation_parameter,
@@ -171,8 +171,7 @@ def TV_PD_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float tolerance_param,
                      int methodTV,
                      int nonneg,
-                     float lipschitz_const,
-                     float tau):
+                     float lipschitz_const):
 
     cdef long dims[2]
     dims[0] = inputData.shape[0]
@@ -191,7 +190,6 @@ def TV_PD_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                        lipschitz_const,
                        methodTV,
                        nonneg,
-                       tau,
                        dims[1],dims[0], 1)
     return (outputData,infovec)
 
@@ -200,9 +198,8 @@ def TV_PD_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      int iterationsNumb,
                      float tolerance_param,
                      int methodTV,
-                     int nonneg,                     
-                     float lipschitz_const,
-                     float tau):
+                     int nonneg,
+                     float lipschitz_const):
 
     cdef long dims[3]
     dims[0] = inputData.shape[0]
@@ -221,7 +218,6 @@ def TV_PD_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                        lipschitz_const,
                        methodTV,
                        nonneg,
-                       tau,
                        dims[2], dims[1], dims[0])
     return (outputData,infovec)
 
diff --git a/src/Python/src/gpu_regularisers.pyx b/src/Python/src/gpu_regularisers.pyx
index b22d15e..8a4568e 100644
--- a/src/Python/src/gpu_regularisers.pyx
+++ b/src/Python/src/gpu_regularisers.pyx
@@ -22,7 +22,7 @@ CUDAErrorMessage = 'CUDA error'
 
 cdef extern int TV_ROF_GPU_main(float* Input, float* Output, float *infovector, float *lambdaPar, int lambda_is_arr, int iter, float tau, float epsil, int N, int M, int Z);
 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_PD_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, float lipschitz_const, int methodTV, int nonneg, float tau, int dimX, int dimY, int dimZ);
+cdef extern int TV_PD_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, float lipschitz_const, int methodTV, int nonneg, int dimX, int dimY, int dimZ);
 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 *infovector, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, float epsil, int dimX, int dimY, int dimZ);
@@ -72,11 +72,11 @@ def TV_FGP_GPU(inputData,
                      methodTV,
                      nonneg)
 # Total-variation Primal-Dual (PD)
-def TV_PD_GPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const, tau):
+def TV_PD_GPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const):
     if inputData.ndim == 2:
-        return TVPD2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const, tau)
+        return TVPD2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const)
     elif inputData.ndim == 3:
-        return TVPD3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const, tau)
+        return TVPD3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const)
 
 def TVPD2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float regularisation_parameter,
@@ -84,8 +84,7 @@ def TVPD2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                      float tolerance_param,
                      int methodTV,
                      int nonneg,
-                     float lipschitz_const,
-                     float tau):
+                     float lipschitz_const):
 
     cdef long dims[2]
     dims[0] = inputData.shape[0]
@@ -103,7 +102,6 @@ def TVPD2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                        lipschitz_const,
                        methodTV,
                        nonneg,
-                       tau,
                        dims[1],dims[0], 1) ==0):
         return (outputData,infovec)
     else:
@@ -115,8 +113,7 @@ def TVPD3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                      float tolerance_param,
                      int methodTV,
                      int nonneg,
-                     float lipschitz_const,
-                     float tau):
+                     float lipschitz_const):
 
     cdef long dims[3]
     dims[0] = inputData.shape[0]
@@ -134,7 +131,6 @@ def TVPD3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                        lipschitz_const,
                        methodTV,
                        nonneg,
-                       tau,
                        dims[2], dims[1], dims[0]) ==0):
         return (outputData,infovec)
     else:
diff --git a/test/test_CPU_regularisers.py b/test/test_CPU_regularisers.py
index 266ca8a..1eba479 100644
--- a/test/test_CPU_regularisers.py
+++ b/test/test_CPU_regularisers.py
@@ -42,7 +42,7 @@ class TestRegularisers(unittest.TestCase):
     def test_PD_TV_CPU(self):
         Im,input,ref = self.getPars()
 
-        pd_cpu,info = PD_TV(input, 0.02, 300, 0.0, 0, 0, 8, 0.0025, 'cpu');
+        pd_cpu,info = PD_TV(input, 0.02, 300, 0.0, 0, 0, 8, 'cpu');
 
         rms = rmse(Im, pd_cpu)
         
diff --git a/test/test_run_test.py b/test/test_run_test.py
index 1707aec..f562593 100755
--- a/test/test_run_test.py
+++ b/test/test_run_test.py
@@ -200,8 +200,7 @@ class TestRegularisers(unittest.TestCase):
                 'tolerance_constant':0.0,\

                 'methodTV': 0 ,\

                 'nonneg': 0,

-                'lipschitz_const' : 8,

-                'tau' : 0.0025}

+                'lipschitz_const' : 8}

                 

         print ("#############PD TV CPU####################")

         start_time = timeit.default_timer()

@@ -211,8 +210,7 @@ class TestRegularisers(unittest.TestCase):
                       pars['tolerance_constant'], 

                       pars['methodTV'],

                       pars['nonneg'],

-                      pars['lipschitz_const'],

-                      pars['tau'],'cpu')

+                      pars['lipschitz_const'],'cpu')

 

         rms = rmse(Im, pd_cpu)

         pars['rmse'] = rms

@@ -230,8 +228,7 @@ class TestRegularisers(unittest.TestCase):
               pars['tolerance_constant'], 

               pars['methodTV'],

               pars['nonneg'],

-              pars['lipschitz_const'],

-              pars['tau'],'gpu')

+              pars['lipschitz_const'],'gpu')

 

         except ValueError as ve:

             self.skipTest("Results not comparable. GPU computing error.")