added PD_TV_2D_CPU version

14 files changed, 541 insertions, 104 deletions

diff --git a/build/run.sh b/build/run.sh
index 285476b..bbbbc6a 100755
--- a/build/run.sh
+++ b/build/run.sh
@@ -6,11 +6,11 @@ rm -r ../build_proj
 mkdir ../build_proj
 cd ../build_proj/
 #make clean
-export CIL_VERSION=19.09
+export CIL_VERSION=19.10
 # install Python modules without CUDA
-#cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=OFF -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
+cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=OFF -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
 # install Python modules with CUDA
-cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
+# cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
 # install Matlab modules without CUDA
 # cmake ../ -DBUILD_PYTHON_WRAPPER=OFF -DMatlab_ROOT_DIR=/dls_sw/apps/matlab/r2014a/ -DBUILD_MATLAB_WRAPPER=ON -DBUILD_CUDA=OFF -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
 # install Matlab modules with CUDA
diff --git a/demos/demo_cpu_regularisers.py b/demos/demo_cpu_regularisers.py
index 8655623..d0bbe63 100644
--- a/demos/demo_cpu_regularisers.py
+++ b/demos/demo_cpu_regularisers.py
@@ -12,7 +12,7 @@ import matplotlib.pyplot as plt
 import numpy as np
 import os
 import timeit
-from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, TNV, NDF, Diff4th
+from ccpi.filters.regularisers import ROF_TV, FGP_TV, PD_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, TNV, NDF, Diff4th
 from ccpi.filters.regularisers import PatchSelect, NLTV
 from ccpi.supp.qualitymetrics import QualityTools
 ###############################################################################
@@ -129,7 +129,7 @@ imgplot = plt.imshow(u0,cmap="gray")
 pars = {'algorithm' : FGP_TV, \
         'input' : u0,\
         'regularisation_parameter':0.02, \
-        'number_of_iterations' :400 ,\
+        'number_of_iterations' :1500 ,\
         'tolerance_constant':1e-06,\
         'methodTV': 0 ,\
         'nonneg': 0}
@@ -161,6 +161,53 @@ plt.title('{}'.format('CPU results'))
 
 #%%
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print ("_______________PD-TV (2D)__________________")
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+
+## plot 
+fig = plt.figure()
+plt.suptitle('Performance of PD-TV regulariser using the CPU')
+a=fig.add_subplot(1,2,1)
+a.set_title('Noisy Image')
+imgplot = plt.imshow(u0,cmap="gray")
+
+# set parameters
+pars = {'algorithm' : PD_TV, \
+        'input' : u0,\
+        'regularisation_parameter':0.02, \
+        'number_of_iterations' :1500 ,\
+        'tolerance_constant':1e-06,\
+        'methodTV': 0 ,\
+        'nonneg': 1,
+        'lipschitz_const' : 12}
+        
+print ("#############PD TV CPU####################")
+start_time = timeit.default_timer()
+(pd_cpu,info_vec_cpu) = PD_TV(pars['input'], 
+              pars['regularisation_parameter'],
+              pars['number_of_iterations'],
+              pars['tolerance_constant'], 
+              pars['methodTV'],
+              pars['nonneg'],
+              pars['lipschitz_const'], 'cpu')
+
+Qtools = QualityTools(Im, pd_cpu)
+pars['rmse'] = Qtools.rmse()
+
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,2,2)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
+# place a text box in upper left in axes coords
+a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow(pd_cpu, cmap="gray")
+plt.title('{}'.format('CPU results'))
+#%%
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("_______________SB-TV (2D)__________________")
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 
diff --git a/src/Core/CMakeLists.txt b/src/Core/CMakeLists.txt
index eea0d63..ac7ec91 100644
--- a/src/Core/CMakeLists.txt
+++ b/src/Core/CMakeLists.txt
@@ -20,7 +20,7 @@ if (OPENMP_FOUND)
     set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}")
    set (CMAKE_SHARED_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_SHARED_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}")
    set (CMAKE_STATIC_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_STATIC_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}")
-   
+
 endif()
 
 ## Build the regularisers package as a library
@@ -39,21 +39,21 @@ if(WIN32)
   set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}")
   set (CMAKE_C_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}")
   set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} /NODEFAULTLIB:MSVCRT.lib")
-  
+
   set (EXTRA_LIBRARIES)
-		
+
   message("library lib: ${LIBRARY_LIB}")
-  
+
 elseif(UNIX)
-   set (FLAGS "-O2 -funsigned-char -Wall  -Wl,--no-undefined  -DCCPiReconstructionIterative_EXPORTS ")  
+   set (FLAGS "-O2 -funsigned-char -Wall  -Wl,--no-undefined  -DCCPiReconstructionIterative_EXPORTS ")
    set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}")
    set (CMAKE_C_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}")
-  
-   set (EXTRA_LIBRARIES 
+
+   set (EXTRA_LIBRARIES
 		"gomp"
 		"m"
 		)
-   
+
 endif()
 message("CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS}")
 
@@ -66,6 +66,7 @@ message("Adding regularisers as a shared library")
 add_library(cilreg SHARED
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/FGP_TV_core.c
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/SB_TV_core.c
+      ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/PD_TV_core.c
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/TGV_core.c
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffusion_core.c
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffus4th_order_core.c
@@ -80,8 +81,8 @@ add_library(cilreg SHARED
 	    ${CMAKE_CURRENT_SOURCE_DIR}/inpainters_CPU/NonlocalMarching_Inpaint_core.c
 	    )
 target_link_libraries(cilreg ${EXTRA_LIBRARIES} )
-include_directories(cilreg PUBLIC 
-                      ${LIBRARY_INC}/include 
+include_directories(cilreg PUBLIC
+                      ${LIBRARY_INC}/include
 					  ${CMAKE_CURRENT_SOURCE_DIR}
 		              ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/
 		              ${CMAKE_CURRENT_SOURCE_DIR}/inpainters_CPU/  )
@@ -92,14 +93,14 @@ if (UNIX)
 message ("I'd install into ${CMAKE_INSTALL_PREFIX}/lib")
 install(TARGETS cilreg
 	LIBRARY DESTINATION lib
-	CONFIGURATIONS ${CMAKE_BUILD_TYPE} 
+	CONFIGURATIONS ${CMAKE_BUILD_TYPE}
 	)
 elseif(WIN32)
 message ("I'd install into ${CMAKE_INSTALL_PREFIX} lib bin")
-  install(TARGETS cilreg 
+  install(TARGETS cilreg
 	RUNTIME DESTINATION bin
 	ARCHIVE DESTINATION lib
-	CONFIGURATIONS ${CMAKE_BUILD_TYPE} 
+	CONFIGURATIONS ${CMAKE_BUILD_TYPE}
 	)
 endif()
 
@@ -126,14 +127,14 @@ if (BUILD_CUDA)
         message ("I'd install into ${CMAKE_INSTALL_PREFIX}/lib")
         install(TARGETS cilregcuda
         LIBRARY DESTINATION lib
-        CONFIGURATIONS ${CMAKE_BUILD_TYPE} 
+        CONFIGURATIONS ${CMAKE_BUILD_TYPE}
         )
       elseif(WIN32)
         message ("I'd install into ${CMAKE_INSTALL_PREFIX} lib bin")
         install(TARGETS cilregcuda
         RUNTIME DESTINATION bin
         ARCHIVE DESTINATION lib
-        CONFIGURATIONS ${CMAKE_BUILD_TYPE} 
+        CONFIGURATIONS ${CMAKE_BUILD_TYPE}
         )
       endif()
     else()
diff --git a/src/Core/regularisers_CPU/FGP_TV_core.c b/src/Core/regularisers_CPU/FGP_TV_core.c
index a16a2e5..ff67af2 100644
--- a/src/Core/regularisers_CPU/FGP_TV_core.c
+++ b/src/Core/regularisers_CPU/FGP_TV_core.c
@@ -46,12 +46,12 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb
     float tk = 1.0f;
     float tkp1 =1.0f;
     int count = 0;
-    
+
     if (dimZ <= 1) {
         /*2D case */
         float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL;
         DimTotal = (long)(dimX*dimY);
-        
+
         if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
         P1 = calloc(DimTotal, sizeof(float));
         P2 = calloc(DimTotal, sizeof(float));
@@ -59,32 +59,32 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb
         P2_prev = calloc(DimTotal, sizeof(float));
         R1 = calloc(DimTotal, sizeof(float));
         R2 = calloc(DimTotal, sizeof(float));
-        
+
         /* begin iterations */
         for(ll=0; ll<iterationsNumb; ll++) {
-            
+
             if ((epsil != 0.0f)  && (ll % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l);
             /* computing the gradient of the objective function */
             Obj_func2D(Input, Output, R1, R2, lambdaPar, (long)(dimX), (long)(dimY));
-            
+
             /* apply nonnegativity */
             if (nonneg == 1) for(j=0; j<DimTotal; j++) {if (Output[j] < 0.0f) Output[j] = 0.0f;}
-            
+
             /*Taking a step towards minus of the gradient*/
             Grad_func2D(P1, P2, Output, R1, R2, lambdaPar, (long)(dimX), (long)(dimY));
-            
+
             /* projection step */
             Proj_func2D(P1, P2, methodTV, DimTotal);
-            
+
             /*updating R and t*/
             tkp1 = (1.0f + sqrtf(1.0f + 4.0f*tk*tk))*0.5f;
             Rupd_func2D(P1, P1_prev, P2, P2_prev, R1, R2, tkp1, tk, DimTotal);
-            
+
             /*storing old values*/
             copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), 1l);
             copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), 1l);
             tk = tkp1;
-            
+
             /* check early stopping criteria */
             if ((epsil != 0.0f)  && (ll % 5 == 0)) {
                 re = 0.0f; re1 = 0.0f;
@@ -105,7 +105,7 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb
         /*3D case*/
         float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P3=NULL, *P1_prev=NULL, *P2_prev=NULL, *P3_prev=NULL, *R1=NULL, *R2=NULL, *R3=NULL;
         DimTotal = (long)(dimX*dimY*dimZ);
-        
+
         if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
         P1 = calloc(DimTotal, sizeof(float));
         P2 = calloc(DimTotal, sizeof(float));
@@ -116,28 +116,28 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb
         R1 = calloc(DimTotal, sizeof(float));
         R2 = calloc(DimTotal, sizeof(float));
         R3 = calloc(DimTotal, sizeof(float));
-        
+
         /* begin iterations */
         for(ll=0; ll<iterationsNumb; ll++) {
-            
+
             if ((epsil != 0.0f)  && (ll % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
-            
+
             /* computing the gradient of the objective function */
             Obj_func3D(Input, Output, R1, R2, R3, lambdaPar, (long)(dimX), (long)(dimY), (long)(dimZ));
-            
+
             /* apply nonnegativity */
             if (nonneg == 1) for(j=0; j<DimTotal; j++) {if (Output[j] < 0.0f) Output[j] = 0.0f;}
-            
+
             /*Taking a step towards minus of the gradient*/
             Grad_func3D(P1, P2, P3, Output, R1, R2, R3, lambdaPar, (long)(dimX), (long)(dimY), (long)(dimZ));
-            
+
             /* projection step */
             Proj_func3D(P1, P2, P3, methodTV, DimTotal);
-            
+
             /*updating R and t*/
             tkp1 = (1.0f + sqrtf(1.0f + 4.0f*tk*tk))*0.5f;
             Rupd_func3D(P1, P1_prev, P2, P2_prev, P3, P3_prev, R1, R2, R3, tkp1, tk, DimTotal);
-            
+
             /* calculate norm - stopping rules*/
             if ((epsil != 0.0f)  && (ll % 5 == 0)) {
                 re = 0.0f; re1 = 0.0f;
@@ -151,22 +151,22 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb
                 if (re < epsil)  count++;
                 if (count > 3) break;
             }
-            
+
             /*storing old values*/
             copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
             copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
             copyIm(P3, P3_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
             tk = tkp1;
         }
-        
+
         if (epsil != 0.0f) free(Output_prev);
         free(P1); free(P2); free(P3); free(P1_prev); free(P2_prev); free(P3_prev); free(R1); free(R2); free(R3);
     }
-    
+
     /*adding info into info_vector */
     infovector[0] = (float)(ll);  /*iterations number (if stopped earlier based on tolerance)*/
     infovector[1] = re;  /* reached tolerance */
-    
+
     return 0;
 }
 
@@ -202,36 +202,6 @@ float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float la
         }}
     return 1;
 }
-float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal)
-{
-    float val1, val2, denom, sq_denom;
-    long i;
-    if (methTV == 0) {
-        /* isotropic TV*/
-#pragma omp parallel for shared(P1,P2) private(i,denom,sq_denom)
-        for(i=0; i<DimTotal; i++) {
-            denom = powf(P1[i],2) +  powf(P2[i],2);
-            if (denom > 1.0f) {
-                sq_denom = 1.0f/sqrtf(denom);
-                P1[i] = P1[i]*sq_denom;
-                P2[i] = P2[i]*sq_denom;
-            }
-        }
-    }
-    else {
-        /* anisotropic TV*/
-#pragma omp parallel for shared(P1,P2) private(i,val1,val2)
-        for(i=0; i<DimTotal; i++) {
-            val1 = fabs(P1[i]);
-            val2 = fabs(P2[i]);
-            if (val1 < 1.0f) {val1 = 1.0f;}
-            if (val2 < 1.0f) {val2 = 1.0f;}
-            P1[i] = P1[i]/val1;
-            P2[i] = P2[i]/val2;
-        }
-    }
-    return 1;
-}
 float Rupd_func2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, long DimTotal)
 {
     long i;
diff --git a/src/Core/regularisers_CPU/FGP_TV_core.h b/src/Core/regularisers_CPU/FGP_TV_core.h
index 04e6e80..4466a35 100644
--- a/src/Core/regularisers_CPU/FGP_TV_core.h
+++ b/src/Core/regularisers_CPU/FGP_TV_core.h
@@ -29,12 +29,12 @@ limitations under the License.
 /* C-OMP implementation of FGP-TV [1] denoising/regularization model (2D/3D case)
  *
  * Input Parameters:
- * 1. Noisy image/volume 
- * 2. lambda - regularization parameter 
+ * 1. Noisy image/volume
+ * 2. lambda - regularization parameter
  * 3. Number of iterations
- * 4. eplsilon: tolerance constant 
+ * 4. eplsilon: tolerance constant
  * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
- * 6. nonneg: 'nonnegativity (0 is OFF by default) 
+ * 6. nonneg: 'nonnegativity (0 is OFF by default)
  *
  * Output:
  * [1] Filtered/regularized image/volume
@@ -44,7 +44,7 @@ limitations under the License.
  * This function is based on the Matlab's code and paper by
  * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems"
  */
- 
+
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -52,7 +52,6 @@ CCPI_EXPORT float TV_FGP_CPU_main(float *Input, float *Output, float *infovector
 
 CCPI_EXPORT float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, long dimX, long dimY);
 CCPI_EXPORT float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float lambda, long dimX, long dimY);
-CCPI_EXPORT float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal);
 CCPI_EXPORT float Rupd_func2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, long DimTotal);
 
 CCPI_EXPORT float Obj_func3D(float *A, float *D, float *R1, float *R2, float *R3, float lambda, long dimX, long dimY, long dimZ);
diff --git a/src/Core/regularisers_CPU/PD_TV_core.c b/src/Core/regularisers_CPU/PD_TV_core.c
new file mode 100644
index 0000000..a8c3cfb
--- /dev/null
+++ b/src/Core/regularisers_CPU/PD_TV_core.c
@@ -0,0 +1,166 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2019 Daniil Kazantsev
+ * Copyright 2019 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "PD_TV_core.h"
+
+/* C-OMP implementation of Primal-Dual TV [1] by Chambolle Pock denoising/regularization model (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume
+ * 2. lambdaPar - regularization parameter
+ * 3. Number of iterations
+ * 4. eplsilon: tolerance constant
+ * 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)
+
+ * Output:
+ * [1] TV - Filtered/regularized image/volume
+ * [2] Information vector which contains [iteration no., reached tolerance]
+ *
+ * [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, int dimX, int dimY, int dimZ)
+{
+    int ll;
+    long j, DimTotal;
+    float re, re1, tau, sigma, theta, lt;
+    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);
+    theta = 1.0f;
+    lt = tau/lambdaPar;
+    ll = 0;
+
+    copyIm(Input, U, (long)(dimX), (long)(dimY), (long)(dimZ));
+
+    if (dimZ <= 1) {
+        /*2D case */
+        float *U_old=NULL, *P1=NULL, *P2=NULL;
+        DimTotal = (long)(dimX*dimY);
+
+        U_old = calloc(DimTotal, sizeof(float));
+        P1 = calloc(DimTotal, sizeof(float));
+        P2 = calloc(DimTotal, sizeof(float));
+
+        /* begin iterations */
+        for(ll=0; ll<iterationsNumb; ll++) {
+
+            //if ((epsil != 0.0f)  && (ll % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l);
+            /* computing the gradient of the objective function */
+            DualP2D(U, P1, P2, (long)(dimX), (long)(dimY), sigma);
+
+            /* apply nonnegativity */
+            if (nonneg == 1) for(j=0; j<DimTotal; j++) {if (U[j] < 0.0f) U[j] = 0.0f;}
+
+            /* projection step */
+            Proj_func2D(P1, P2, methodTV, DimTotal);
+
+            /* copy U to U_old */
+            copyIm(U, U_old, (long)(dimX), (long)(dimY), 1l);
+
+            DivProj2D(U, Input, P1, P2,(long)(dimX), (long)(dimY), lt, tau);
+
+            /* check early stopping criteria */
+            if ((epsil != 0.0f)  && (ll % 5 == 0)) {
+                re = 0.0f; re1 = 0.0f;
+                for(j=0; j<DimTotal; j++)
+                {
+                    re += powf(U[j] - U_old[j],2);
+                    re1 += powf(U[j],2);
+                }
+                re = sqrtf(re)/sqrtf(re1);
+                if (re < epsil)  count++;
+                if (count > 3) break;
+            }
+            /*get updated solution*/
+
+            getX2D(U, U_old, theta, DimTotal);
+        }
+        free(P1); free(P2); free(U_old);
+    }
+    else {
+        /*3D case*/
+    }
+    /*adding info into info_vector */
+    infovector[0] = (float)(ll);  /*iterations number (if stopped earlier based on tolerance)*/
+    infovector[1] = re;  /* reached tolerance */
+
+    return 0;
+}
+
+/*****************************************************************/
+/************************2D-case related Functions */
+/*****************************************************************/
+
+/*Calculating dual variable (using forward differences)*/
+float DualP2D(float *U, float *P1, float *P2, long dimX, long dimY, float sigma)
+{
+     long i,j,index;
+     #pragma omp parallel for shared(U,P1,P2) private(index,i,j)
+     for(j=0; j<dimY; j++) {
+       for(i=0; i<dimX; i++) {
+          index = j*dimX+i;
+          /* symmetric boundary conditions (Neuman) */
+          if (i == dimX-1) P1[index] += sigma*(U[j*dimX+(i-1)] - U[index]);
+          else P1[index] += sigma*(U[j*dimX+(i+1)] - U[index]);
+          if (j == dimY-1) P2[index] += sigma*(U[(j-1)*dimX+i] - U[index]);
+          else  P2[index] += sigma*(U[(j+1)*dimX+i] - U[index]);
+        }}
+     return 1;
+}
+
+/* Divergence for P dual */
+float DivProj2D(float *U, float *Input, float *P1, float *P2, long dimX, long dimY, float lt, float tau)
+{
+  long i,j,index;
+  float P_v1, P_v2, div_var;
+  #pragma omp parallel for shared(U,Input,P1,P2) private(i, j, P_v1, P_v2, div_var)
+  for(j=0; j<dimY; j++) {
+    for(i=0; i<dimX; i++) {
+            index = j*dimX+i;
+            /* symmetric boundary conditions (Neuman) */
+            if (i == 0) P_v1 = -P1[index];
+            else P_v1 = -(P1[index] - P1[j*dimX+(i-1)]);
+            if (j == 0) P_v2 = -P2[index];
+            else  P_v2 = -(P2[index] - P2[(j-1)*dimX+i]);
+            div_var = P_v1 + P_v2;
+            U[index] = (U[index] - tau*div_var + lt*Input[index])/(1.0 + lt);
+          }}
+  return *U;
+}
+
+/*get the updated solution*/
+float getX2D(float *U, float *U_old, float theta, long DimTotal)
+{
+    long i;
+    #pragma omp parallel for shared(U,U_old) private(i)
+    for(i=0; i<DimTotal; i++) {
+          U[i] +=  theta*(U[i] - U_old[i]);
+          }
+    return *U;
+}
+
+
+/*****************************************************************/
+/************************3D-case related Functions */
+/*****************************************************************/
diff --git a/src/Core/regularisers_CPU/PD_TV_core.h b/src/Core/regularisers_CPU/PD_TV_core.h
new file mode 100644
index 0000000..b52689a
--- /dev/null
+++ b/src/Core/regularisers_CPU/PD_TV_core.h
@@ -0,0 +1,57 @@
+/*
+This work is part of the Core Imaging Library developed by
+Visual Analytics and Imaging System Group of the Science Technology
+Facilities Council, STFC
+
+Copyright 2019 Daniil Kazantsev
+Copyright 2019 Srikanth Nagella, Edoardo Pasca
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+http://www.apache.org/licenses/LICENSE-2.0
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+//#include <matrix.h>
+#include <math.h>
+#include <stdlib.h>
+#include <memory.h>
+#include <stdio.h>
+#include "omp.h"
+#include "utils.h"
+#include "CCPiDefines.h"
+
+/* C-OMP implementation of Primal-Dual TV [1] by Chambolle Pock denoising/regularization model (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume
+ * 2. lambdaPar - regularization parameter
+ * 3. Number of iterations
+ * 4. eplsilon: tolerance constant
+ * 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)
+
+ * Output:
+ * [1] TV - Filtered/regularized image/volume
+ * [2] Information vector which contains [iteration no., reached tolerance]
+ *
+ * [1] Antonin Chambolle, Thomas Pock. "A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging", 2010
+ */
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+CCPI_EXPORT 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);
+CCPI_EXPORT float getX2D(float *U, float *U_old, float theta, long DimTotal);
+#ifdef __cplusplus
+}
+#endif
diff --git a/src/Core/regularisers_CPU/utils.c b/src/Core/regularisers_CPU/utils.c
index 5bb3a5c..cf4ad72 100644
--- a/src/Core/regularisers_CPU/utils.c
+++ b/src/Core/regularisers_CPU/utils.c
@@ -65,7 +65,7 @@ float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int type, int
 {
     int i, j, i1, j1, index;
     float NOMx_2, NOMy_2, E_Grad=0.0f, E_Data=0.0f;
-    
+
     /* first calculate \grad U_xy*/
     for(j=0; j<dimY; j++) {
         for(i=0; i<dimX; i++) {
@@ -73,7 +73,7 @@ float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int type, int
             /* boundary conditions */
             i1 = i + 1; if (i == dimX-1) i1 = i;
             j1 = j + 1; if (j == dimY-1) j1 = j;
-            
+
             /* Forward differences */
             NOMx_2 = powf((float)(U[j1*dimX + i] - U[index]),2); /* x+ */
             NOMy_2 = powf((float)(U[j*dimX + i1] - U[index]),2); /* y+ */
@@ -90,7 +90,7 @@ float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, int
 {
     long i, j, k, i1, j1, k1, index;
     float NOMx_2, NOMy_2, NOMz_2, E_Grad=0.0f, E_Data=0.0f;
-    
+
     /* first calculate \grad U_xy*/
     for(j=0; j<(long)(dimY); j++) {
         for(i=0; i<(long)(dimX); i++) {
@@ -100,12 +100,12 @@ float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, int
                 i1 = i + 1; if (i == (long)(dimX-1)) i1 = i;
                 j1 = j + 1; if (j == (long)(dimY-1)) j1 = j;
                 k1 = k + 1; if (k == (long)(dimZ-1)) k1 = k;
-                
+
                 /* Forward differences */
                 NOMx_2 = powf((float)(U[(dimX*dimY)*k + j1*dimX+i] - U[index]),2); /* x+ */
                 NOMy_2 = powf((float)(U[(dimX*dimY)*k + j*dimX+i1] - U[index]),2); /* y+ */
                 NOMz_2 = powf((float)(U[(dimX*dimY)*k1 + j*dimX+i] - U[index]),2); /* z+ */
-                
+
                 E_Grad += 2.0f*lambda*sqrtf((float)(NOMx_2) + (float)(NOMy_2) + (float)(NOMz_2)); /* gradient term energy */
                 E_Data += (powf((float)(U[index]-U0[index]),2)); /* fidelity term energy */
             }
@@ -131,12 +131,12 @@ float Im_scale2D(float *Input, float *Scaled, int w, int h, int w2, int h2)
             x_diff = (x_ratio * j) - x;
             y_diff = (y_ratio * i) - y;
             index = y*w+x ;
-            
+
             A = Input[index];
             B = Input[index+1];
             C = Input[index+w];
             D = Input[index+w+1];
-            
+
             gray = (float)(A*(1.0 - x_diff)*(1.0 - y_diff) +  B*(x_diff)*(1.0 - y_diff) +
                     C*(y_diff)*(1.0 - x_diff) +  D*(x_diff*y_diff));
 
@@ -144,3 +144,35 @@ float Im_scale2D(float *Input, float *Scaled, int w, int h, int w2, int h2)
         }}
     return *Scaled;
 }
+
+/*2D Projection onto convex set for P*/
+float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal)
+{
+    float val1, val2, denom, sq_denom;
+    long i;
+    if (methTV == 0) {
+        /* isotropic TV*/
+#pragma omp parallel for shared(P1,P2) private(i,denom,sq_denom)
+        for(i=0; i<DimTotal; i++) {
+            denom = powf(P1[i],2) +  powf(P2[i],2);
+            if (denom > 1.0f) {
+                sq_denom = 1.0f/sqrtf(denom);
+                P1[i] = P1[i]*sq_denom;
+                P2[i] = P2[i]*sq_denom;
+            }
+        }
+    }
+    else {
+        /* anisotropic TV*/
+#pragma omp parallel for shared(P1,P2) private(i,val1,val2)
+        for(i=0; i<DimTotal; i++) {
+            val1 = fabs(P1[i]);
+            val2 = fabs(P2[i]);
+            if (val1 < 1.0f) {val1 = 1.0f;}
+            if (val2 < 1.0f) {val2 = 1.0f;}
+            P1[i] = P1[i]/val1;
+            P2[i] = P2[i]/val2;
+        }
+    }
+    return 1;
+}
diff --git a/src/Core/regularisers_CPU/utils.h b/src/Core/regularisers_CPU/utils.h
index 8f0ba59..e050f86 100644
--- a/src/Core/regularisers_CPU/utils.h
+++ b/src/Core/regularisers_CPU/utils.h
@@ -31,6 +31,7 @@ CCPI_EXPORT float TV_energy2D(float *U, float *U0, float *E_val, float lambda, i
 CCPI_EXPORT float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY, int dimZ);
 CCPI_EXPORT float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY, int dimZ);
 CCPI_EXPORT float Im_scale2D(float *Input, float *Scaled, int w, int h, int w2, int h2);
+CCPI_EXPORT float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal);
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/Matlab/mex_compile/compileCPU_mex_Linux.m b/src/Matlab/mex_compile/compileCPU_mex_Linux.m
index 2ed7ea6..5330d7f 100644
--- a/src/Matlab/mex_compile/compileCPU_mex_Linux.m
+++ b/src/Matlab/mex_compile/compileCPU_mex_Linux.m
@@ -28,6 +28,10 @@ movefile('FGP_TV.mex*',Pathmove);
 fprintf('%s \n', 'Compiling SB-TV...');
 mex SB_TV.c SB_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
 movefile('SB_TV.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling PD-TV...');
+mex PD_TV.c PD_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('PD_TV.mex*',Pathmove);
  
 fprintf('%s \n', 'Compiling dFGP-TV...');
 mex FGP_dTV.c FGP_dTV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
@@ -75,7 +79,8 @@ movefile('NonlocalMarching_Inpaint.mex*',Pathmove);
 delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* Diffus4th_order_core* TGV_core* LLT_ROF_core* CCPiDefines.h
 delete PatchSelect_core* Nonlocal_TV_core*
 delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core*
-fprintf('%s \n', '<<<<<<< Regularisers successfully compiled! >>>>>>>');
+delete PD_TV_core*
+fprintf('%s \n', '<<<<<<< CPU regularisers were successfully compiled! >>>>>>>');
 
 pathA2 = sprintf(['..' fsep '..' fsep '..' fsep '..' fsep 'demos' fsep 'Matlab_demos'], 1i);
 cd(pathA2);
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c b/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c
new file mode 100644
index 0000000..eac2d18
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c
@@ -0,0 +1,98 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2019 Daniil Kazantsev
+ * Copyright 2019 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "PD_TV_core.h"
+
+/* C-OMP implementation of Primal-Dual TV [1] by Chambolle Pock denoising/regularization model (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume
+ * 2. lambdaPar - regularization parameter
+ * 3. Number of iterations
+ * 4. eplsilon: tolerance constant
+ * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
+ * 6. nonneg: 'nonnegativity (0 is OFF by default, 1 is ON)
+ * 7. lipschitz_const: convergence related parameter
+ 
+ * Output:
+ * [1] TV - Filtered/regularized image/volume
+ * [2] Information vector which contains [iteration no., reached tolerance]
+ *
+ * [1] Antonin Chambolle, Thomas Pock. "A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging", 2010
+ */
+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;
+    
+    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");
+    
+    Input  = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
+    lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
+    iter = 400; /* default iterations number */
+    epsil = 1.0e-06; /* default tolerance constant */
+    methTV = 0;  /* default isotropic TV penalty */
+    nonneg = 0; /* default nonnegativity switch, off - 0 */
+    lipschitz_const = 12.0; /* lipschitz_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))  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))  {
+        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) lipschitz_const = (float) mxGetScalar(prhs[6]);
+    
+    
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];    
+       
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));        
+    }
+    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, dimX, dimY, dimZ);
+}
diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py
index 0b5b2ee..d65c0b9 100644
--- a/src/Python/ccpi/filters/regularisers.py
+++ b/src/Python/ccpi/filters/regularisers.py
@@ -2,7 +2,7 @@
 script which assigns a proper device core function based on a flag ('cpu' or 'gpu')
 """
 
-from ccpi.filters.cpu_regularisers import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, Diff4th_CPU, TGV_CPU, LLT_ROF_CPU, PATCHSEL_CPU, NLTV_CPU
+from ccpi.filters.cpu_regularisers import TV_ROF_CPU, TV_FGP_CPU, TV_PD_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, Diff4th_CPU, TGV_CPU, LLT_ROF_CPU, PATCHSEL_CPU, NLTV_CPU
 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
@@ -51,6 +51,31 @@ def FGP_TV(inputData, regularisation_parameter,iterations,
             raise ValueError ('GPU is not available')
         raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
                          .format(device))
+
+def PD_TV(inputData, regularisation_parameter, iterations,
+                     tolerance_param, methodTV, nonneg, lipschitz_const, device='cpu'):
+    if device == 'cpu':
+        return TV_PD_CPU(inputData,
+                     regularisation_parameter,
+                     iterations,
+                     tolerance_param,
+                     methodTV,
+                     nonneg,
+                     lipschitz_const)
+    elif device == 'gpu' and gpu_enabled:
+        return TV_PD_CPU(inputData,
+                     regularisation_parameter,
+                     iterations,
+                     tolerance_param,
+                     methodTV,
+                     nonneg,
+                     lipschitz_const)
+    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 SB_TV(inputData, regularisation_parameter, iterations,
                      tolerance_param, methodTV, device='cpu'):
     if device == 'cpu':
@@ -212,4 +237,3 @@ def NDF_INP(inputData, maskData, regularisation_parameter, edge_parameter, itera
 
 def NVM_INP(inputData, maskData, SW_increment, iterations):
         return NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterations)
-
diff --git a/src/Python/setup-regularisers.py.in b/src/Python/setup-regularisers.py.in
index 4c578e3..9bcd46d 100644
--- a/src/Python/setup-regularisers.py.in
+++ b/src/Python/setup-regularisers.py.in
@@ -8,13 +8,13 @@ from Cython.Distutils import build_ext
 import os
 import sys
 import numpy
-import platform	
+import platform
 
 cil_version=os.environ['CIL_VERSION']
 if  cil_version == '':
     print("Please set the environmental variable CIL_VERSION")
     sys.exit(1)
-	
+
 library_include_path = ""
 library_lib_path = ""
 try:
@@ -23,7 +23,7 @@ try:
 except:
     library_include_path = os.environ['PREFIX']+'/include'
     pass
-    
+
 extra_include_dirs = [numpy.get_include(), library_include_path]
 #extra_library_dirs = [os.path.join(library_include_path, "..", "lib")]
 extra_compile_args = []
@@ -38,6 +38,7 @@ extra_include_dirs += [os.path.join(".." , "Core"),
                        os.path.join(".." , "Core",  "regularisers_CPU"),
                        os.path.join(".." , "Core",  "inpainters_CPU"),
                        os.path.join(".." , "Core",  "regularisers_GPU" , "TV_FGP" ) ,
+                       os.path.join(".." , "Core",  "regularisers_GPU" , "TV_PD" ) ,
                        os.path.join(".." , "Core",  "regularisers_GPU" , "TV_ROF" ) ,
                        os.path.join(".." , "Core",  "regularisers_GPU" , "TV_SB" ) ,
                        os.path.join(".." , "Core",  "regularisers_GPU" , "TGV" ) ,
@@ -48,12 +49,12 @@ extra_include_dirs += [os.path.join(".." , "Core"),
                        os.path.join(".." , "Core",  "regularisers_GPU" , "PatchSelect" ) ,
 						   "."]
 
-if platform.system() == 'Windows':				   
-    extra_compile_args[0:] = ['/DWIN32','/EHsc','/DBOOST_ALL_NO_LIB' , '/openmp' ]   
+if platform.system() == 'Windows':
+    extra_compile_args[0:] = ['/DWIN32','/EHsc','/DBOOST_ALL_NO_LIB' , '/openmp' ]
 else:
     extra_compile_args = ['-fopenmp','-O2', '-funsigned-char', '-Wall', '-std=c++0x']
     extra_libraries += [@EXTRA_OMP_LIB@]
-    
+
 setup(
     name='ccpi',
 	description='CCPi Core Imaging Library - Image regularisers',
@@ -61,13 +62,13 @@ setup(
     cmdclass = {'build_ext': build_ext},
     ext_modules = [Extension("ccpi.filters.cpu_regularisers",
                              sources=[os.path.join("." , "src", "cpu_regularisers.pyx" ) ],
-                             include_dirs=extra_include_dirs, 
-							 library_dirs=extra_library_dirs, 
-							 extra_compile_args=extra_compile_args, 
-							 libraries=extra_libraries ), 
-    
+                             include_dirs=extra_include_dirs,
+							 library_dirs=extra_library_dirs,
+							 extra_compile_args=extra_compile_args,
+							 libraries=extra_libraries ),
+
     ],
-	zip_safe = False,	
+	zip_safe = False,
 	packages = {'ccpi', 'ccpi.filters', 'ccpi.supp'},
 )
 
diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx
index 4917d06..724634b 100644
--- a/src/Python/src/cpu_regularisers.pyx
+++ b/src/Python/src/cpu_regularisers.pyx
@@ -20,6 +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, 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);
@@ -58,9 +59,6 @@ def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
             np.ones([2], dtype='float32')
 
-    # 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)
-     # Run ROF iterations for 2D data
     if isinstance (regularisation_parameter, np.ndarray):
         reg = regularisation_parameter.copy()
         TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], &reg[0,0],  1, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1)
@@ -158,6 +156,44 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
                        dims[2], dims[1], dims[0])
     return (outputData,infovec)
 
+#****************************************************************#
+#****************** Total-variation Primal-dual *****************#
+#****************************************************************#
+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)
+    elif inputData.ndim == 3:
+        return 0
+
+def TV_PD_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+                     float regularisation_parameter,
+                     int iterationsNumb,
+                     float tolerance_param,
+                     int methodTV,
+                     int nonneg,
+                     float lipschitz_const):
+
+    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 */
+    PDTV_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter,
+                       iterationsNumb,
+                       tolerance_param,
+                       lipschitz_const,
+                       methodTV,
+                       nonneg,
+                       dims[1],dims[0], 1)
+
+    return (outputData,infovec)
+
 #***************************************************************#
 #********************** Total-variation SB *********************#
 #***************************************************************#