From 494a857b830fce5e786dfc058f68bf78d9673ba6 Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Tue, 26 Nov 2019 14:51:34 +0000
Subject: PDTV 3D verision

---
 src/Core/CMakeLists.txt                  |   2 +-
 src/Core/regularisers_CPU/FGP_TV_core.c  |  34 --------
 src/Core/regularisers_CPU/FGP_TV_core.h  |   1 -
 src/Core/regularisers_CPU/FGP_dTV_core.c | 141 +++++++++----------------------
 src/Core/regularisers_CPU/FGP_dTV_core.h |   2 -
 src/Core/regularisers_CPU/PD_TV_core.c   |  95 +++++++++++++++++++--
 src/Core/regularisers_CPU/PD_TV_core.h   |   5 +-
 src/Core/regularisers_CPU/utils.c        |  37 +++++++-
 src/Core/regularisers_CPU/utils.h        |   1 +
 src/Python/setup-regularisers.py.in      |   1 -
 src/Python/src/cpu_regularisers.pyx      |  29 ++++++-
 11 files changed, 196 insertions(+), 152 deletions(-)

(limited to 'src')

diff --git a/src/Core/CMakeLists.txt b/src/Core/CMakeLists.txt
index ac7ec91..c1bd7bb 100644
--- a/src/Core/CMakeLists.txt
+++ b/src/Core/CMakeLists.txt
@@ -66,7 +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/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
diff --git a/src/Core/regularisers_CPU/FGP_TV_core.c b/src/Core/regularisers_CPU/FGP_TV_core.c
index ff67af2..12f2254 100644
--- a/src/Core/regularisers_CPU/FGP_TV_core.c
+++ b/src/Core/regularisers_CPU/FGP_TV_core.c
@@ -254,40 +254,6 @@ float Grad_func3D(float *P1, float *P2, float *P3, float *D, float *R1, float *R
             }}}
     return 1;
 }
-float Proj_func3D(float *P1, float *P2, float *P3, int methTV, long DimTotal)
-{
-    float val1, val2, val3, denom, sq_denom;
-    long i;
-    if (methTV == 0) {
-        /* isotropic TV*/
-#pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3,sq_denom)
-        for(i=0; i<DimTotal; i++) {
-            denom = powf(P1[i],2) + powf(P2[i],2) + powf(P3[i],2);
-            if (denom > 1.0f) {
-                sq_denom = 1.0f/sqrtf(denom);
-                P1[i] = P1[i]*sq_denom;
-                P2[i] = P2[i]*sq_denom;
-                P3[i] = P3[i]*sq_denom;
-            }
-        }
-    }
-    else {
-        /* anisotropic TV*/
-#pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3)
-        for(i=0; i<DimTotal; i++) {
-            val1 = fabs(P1[i]);
-            val2 = fabs(P2[i]);
-            val3 = fabs(P3[i]);
-            if (val1 < 1.0f) {val1 = 1.0f;}
-            if (val2 < 1.0f) {val2 = 1.0f;}
-            if (val3 < 1.0f) {val3 = 1.0f;}
-            P1[i] = P1[i]/val1;
-            P2[i] = P2[i]/val2;
-            P3[i] = P3[i]/val3;
-        }
-    }
-    return 1;
-}
 float Rupd_func3D(float *P1, float *P1_old, float *P2, float *P2_old, float *P3, float *P3_old, float *R1, float *R2, float *R3, 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 4466a35..e9c3cde 100644
--- a/src/Core/regularisers_CPU/FGP_TV_core.h
+++ b/src/Core/regularisers_CPU/FGP_TV_core.h
@@ -56,7 +56,6 @@ CCPI_EXPORT float Rupd_func2D(float *P1, float *P1_old, float *P2, float *P2_old
 
 CCPI_EXPORT float Obj_func3D(float *A, float *D, float *R1, float *R2, float *R3, float lambda, long dimX, long dimY, long dimZ);
 CCPI_EXPORT float Grad_func3D(float *P1, float *P2, float *P3, float *D, float *R1, float *R2, float *R3, float lambda, long dimX, long dimY, long dimZ);
-CCPI_EXPORT float Proj_func3D(float *P1, float *P2, float *P3, int methTV, long DimTotal);
 CCPI_EXPORT float Rupd_func3D(float *P1, float *P1_old, float *P2, float *P2_old, float *P3, float *P3_old, float *R1, float *R2, float *R3, float tkp1, float tk, long DimTotal);
 #ifdef __cplusplus
 }
diff --git a/src/Core/regularisers_CPU/FGP_dTV_core.c b/src/Core/regularisers_CPU/FGP_dTV_core.c
index afd7264..e828be6 100644
--- a/src/Core/regularisers_CPU/FGP_dTV_core.c
+++ b/src/Core/regularisers_CPU/FGP_dTV_core.c
@@ -52,11 +52,11 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *info
     float tk = 1.0f;
     float tkp1=1.0f;
     int count = 0;
-    
-    
+
+
     float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL, *InputRef_x=NULL, *InputRef_y=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));
@@ -66,39 +66,39 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *info
     R2 = calloc(DimTotal, sizeof(float));
     InputRef_x = calloc(DimTotal, sizeof(float));
     InputRef_y = calloc(DimTotal, sizeof(float));
-    
+
     if (dimZ <= 1) {
         /*2D case */
         /* calculate gradient field (smoothed) for the reference image */
         GradNorm_func2D(InputRef, InputRef_x, InputRef_y, eta, (long)(dimX), (long)(dimY));
-        
+
         /* begin iterations */
         for(ll=0; ll<iterationsNumb; ll++) {
-            
+
             if ((epsil != 0.0f)  && (ll % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l);
             /*projects a 2D vector field R-1,2 onto the orthogonal complement of another 2D vector field InputRef_xy*/
             ProjectVect_func2D(R1, R2, InputRef_x, InputRef_y, (long)(dimX), (long)(dimY));
-            
+
             /* computing the gradient of the objective function */
             Obj_dfunc2D(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_dfunc2D(P1, P2, Output, R1, R2, InputRef_x, InputRef_y, lambdaPar, (long)(dimX), (long)(dimY));
-            
+
             /* projection step */
-            Proj_dfunc2D(P1, P2, methodTV, DimTotal);
-            
+            Proj_func2D(P1, P2, methodTV, DimTotal);
+
             /*updating R and t*/
             tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f;
             Rupd_dfunc2D(P1, P1_prev, P2, P2_prev, R1, R2, tkp1, tk, DimTotal);
-            
+
             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;
@@ -116,45 +116,45 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *info
     else {
         /*3D case*/
         float *P3=NULL, *P3_prev=NULL, *R3=NULL, *InputRef_z=NULL;
-        
+
         P3 = calloc(DimTotal, sizeof(float));
         P3_prev = calloc(DimTotal, sizeof(float));
         R3 = calloc(DimTotal, sizeof(float));
         InputRef_z = calloc(DimTotal, sizeof(float));
-        
+
         /* calculate gradient field (smoothed) for the reference volume */
         GradNorm_func3D(InputRef, InputRef_x, InputRef_y, InputRef_z, eta, (long)(dimX), (long)(dimY), (long)(dimZ));
-        
+
         /* 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));
-            
+
             /*projects a 3D vector field R-1,2,3 onto the orthogonal complement of another 3D vector field InputRef_xyz*/
             ProjectVect_func3D(R1, R2, R3, InputRef_x, InputRef_y, InputRef_z, (long)(dimX), (long)(dimY), (long)(dimZ));
-            
+
             /* computing the gradient of the objective function */
             Obj_dfunc3D(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_dfunc3D(P1, P2, P3, Output, R1, R2, R3, InputRef_x, InputRef_y, InputRef_z, lambdaPar, (long)(dimX), (long)(dimY), (long)(dimZ));
-            
+
             /* projection step */
-            Proj_dfunc3D(P1, P2, P3, methodTV, DimTotal);
-            
+            Proj_func3D(P1, P2, P3, methodTV, DimTotal);
+
             /*updating R and t*/
             tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f;
             Rupd_dfunc3D(P1, P1_prev, P2, P2_prev, P3, P3_prev, R1, R2, R3, tkp1, tk, DimTotal);
-            
+
             /*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;
-            
+
             /* check early stopping criteria */
             if ((epsil != 0.0f)  && (ll % 5 == 0)) {
                 re = 0.0f; re1 = 0.0f;
@@ -168,16 +168,16 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *info
                 if (count > 3) break;
             }
         }
-        
+
         free(P3); free(P3_prev); free(R3); free(InputRef_z);
     }
     if (epsil != 0.0f) free(Output_prev);
     free(P1); free(P2); free(P1_prev); free(P2_prev); free(R1); free(R2); free(InputRef_x); free(InputRef_y);
-    
+
     /*adding info into info_vector */
     infovector[0] = (float)(ll);  /*iterations number (if stopped earlier based on tolerance)*/
     infovector[1] = re;  /* reached tolerance */
-    
+
     return 0;
 }
 
@@ -185,7 +185,6 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *info
 /********************************************************************/
 /***************************2D Functions*****************************/
 /********************************************************************/
-
 float GradNorm_func2D(float *B, float *B_x, float *B_y, float eta, long dimX, long dimY)
 {
     long i,j,index;
@@ -249,47 +248,17 @@ float Grad_dfunc2D(float *P1, float *P2, float *D, float *R1, float *R2, float *
             /* boundary conditions */
             if (i == dimX-1) val1 = 0.0f; else val1 = D[index] - D[j*dimX + (i+1)];
             if (j == dimY-1) val2 = 0.0f; else val2 = D[index] - D[(j+1)*dimX + i];
-            
+
             in_prod = val1*B_x[index] + val2*B_y[index];   /* calculate inner product */
             val1 = val1 - in_prod*B_x[index];
             val2 = val2 - in_prod*B_y[index];
-            
+
             P1[index] = R1[index] + multip*val1;
             P2[index] = R2[index] + multip*val2;
-            
+
         }}
     return 1;
 }
-float Proj_dfunc2D(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_dfunc2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, long DimTotal)
 {
     long i;
@@ -314,14 +283,14 @@ float GradNorm_func3D(float *B, float *B_x, float *B_y, float *B_z, float eta, l
     for(k=0; k<dimZ; k++) {
         for(j=0; j<dimY; j++) {
             for(i=0; i<dimX; i++) {
-                
+
                 index = (dimX*dimY)*k + j*dimX+i;
-                
+
                 /* zero boundary conditions */
                 if (i == dimX-1) {val1 = 0.0f;} else {val1 = B[(dimX*dimY)*k + j*dimX+(i+1)];}
                 if (j == dimY-1) {val2 = 0.0f;} else {val2 = B[(dimX*dimY)*k + (j+1)*dimX+i];}
                 if (k == dimZ-1) {val3 = 0.0f;} else {val3 = B[(dimX*dimY)*(k+1) + (j)*dimX+i];}
-                
+
                 gradX = val1 - B[index];
                 gradY = val2 - B[index];
                 gradZ = val3 - B[index];
@@ -382,52 +351,18 @@ float Grad_dfunc3D(float *P1, float *P2, float *P3, float *D, float *R1, float *
                 if (i == dimX-1) val1 = 0.0f; else val1 = D[index] - D[(dimX*dimY)*k + j*dimX + (i+1)];
                 if (j == dimY-1) val2 = 0.0f; else val2 = D[index] - D[(dimX*dimY)*k + (j+1)*dimX + i];
                 if (k == dimZ-1) val3 = 0.0f; else val3 = D[index] - D[(dimX*dimY)*(k+1) + j*dimX + i];
-                
+
                 in_prod = val1*B_x[index] + val2*B_y[index] + val3*B_z[index];   /* calculate inner product */
                 val1 = val1 - in_prod*B_x[index];
                 val2 = val2 - in_prod*B_y[index];
                 val3 = val3 - in_prod*B_z[index];
-                
+
                 P1[index] = R1[index] + multip*val1;
                 P2[index] = R2[index] + multip*val2;
                 P3[index] = R3[index] + multip*val3;
             }}}
     return 1;
 }
-float Proj_dfunc3D(float *P1, float *P2, float *P3, int methTV, long DimTotal)
-{
-    float val1, val2, val3, denom, sq_denom;
-    long i;
-    if (methTV == 0) {
-        /* isotropic TV*/
-#pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3,sq_denom)
-        for(i=0; i<DimTotal; i++) {
-            denom = powf(P1[i],2) + powf(P2[i],2) + powf(P3[i],2);
-            if (denom > 1.0f) {
-                sq_denom = 1.0f/sqrtf(denom);
-                P1[i] = P1[i]*sq_denom;
-                P2[i] = P2[i]*sq_denom;
-                P3[i] = P3[i]*sq_denom;
-            }
-        }
-    }
-    else {
-        /* anisotropic TV*/
-#pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3)
-        for(i=0; i<DimTotal; i++) {
-            val1 = fabs(P1[i]);
-            val2 = fabs(P2[i]);
-            val3 = fabs(P3[i]);
-            if (val1 < 1.0f) {val1 = 1.0f;}
-            if (val2 < 1.0f) {val2 = 1.0f;}
-            if (val3 < 1.0f) {val3 = 1.0f;}
-            P1[i] = P1[i]/val1;
-            P2[i] = P2[i]/val2;
-            P3[i] = P3[i]/val3;
-        }
-    }
-    return 1;
-}
 float Rupd_dfunc3D(float *P1, float *P1_old, float *P2, float *P2_old, float *P3, float *P3_old, float *R1, float *R2, float *R3, float tkp1, float tk, long DimTotal)
 {
     long i;
diff --git a/src/Core/regularisers_CPU/FGP_dTV_core.h b/src/Core/regularisers_CPU/FGP_dTV_core.h
index 9ace06d..8582170 100644
--- a/src/Core/regularisers_CPU/FGP_dTV_core.h
+++ b/src/Core/regularisers_CPU/FGP_dTV_core.h
@@ -59,14 +59,12 @@ CCPI_EXPORT float GradNorm_func2D(float *B, float *B_x, float *B_y, float eta, l
 CCPI_EXPORT float ProjectVect_func2D(float *R1, float *R2, float *B_x, float *B_y, long dimX, long dimY);
 CCPI_EXPORT float Obj_dfunc2D(float *A, float *D, float *R1, float *R2, float lambda, long dimX, long dimY);
 CCPI_EXPORT float Grad_dfunc2D(float *P1, float *P2, float *D, float *R1, float *R2, float *B_x, float *B_y, float lambda, long dimX, long dimY);
-CCPI_EXPORT float Proj_dfunc2D(float *P1, float *P2, int methTV, long DimTotal);
 CCPI_EXPORT float Rupd_dfunc2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, long DimTotal);
 
 CCPI_EXPORT float GradNorm_func3D(float *B, float *B_x, float *B_y, float *B_z, float eta, long dimX, long dimY, long dimZ);
 CCPI_EXPORT float ProjectVect_func3D(float *R1, float *R2, float *R3, float *B_x, float *B_y, float *B_z, long dimX, long dimY, long dimZ);
 CCPI_EXPORT float Obj_dfunc3D(float *A, float *D, float *R1, float *R2, float *R3, float lambda, long dimX, long dimY, long dimZ);
 CCPI_EXPORT float Grad_dfunc3D(float *P1, float *P2, float *P3, float *D, float *R1, float *R2, float *R3, float *B_x, float *B_y, float *B_z, float lambda, long dimX, long dimY, long dimZ);
-CCPI_EXPORT float Proj_dfunc3D(float *P1, float *P2, float *P3, int methTV, long DimTotal);
 CCPI_EXPORT float Rupd_dfunc3D(float *P1, float *P1_old, float *P2, float *P2_old, float *P3, float *P3_old, float *R1, float *R2, float *R3, float tkp1, float tk, long DimTotal);
 #ifdef __cplusplus
 }
diff --git a/src/Core/regularisers_CPU/PD_TV_core.c b/src/Core/regularisers_CPU/PD_TV_core.c
index a8c3cfb..f476a8b 100644
--- a/src/Core/regularisers_CPU/PD_TV_core.c
+++ b/src/Core/regularisers_CPU/PD_TV_core.c
@@ -50,13 +50,13 @@ float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar,
     theta = 1.0f;
     lt = tau/lambdaPar;
     ll = 0;
+    DimTotal = (long)(dimX*dimY*dimZ);
 
     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));
@@ -65,8 +65,7 @@ float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar,
         /* 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 */
+            /* computing the the dual P variable */
             DualP2D(U, P1, P2, (long)(dimX), (long)(dimY), sigma);
 
             /* apply nonnegativity */
@@ -94,12 +93,52 @@ float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar,
             }
             /*get updated solution*/
 
-            getX2D(U, U_old, theta, DimTotal);
+            getX(U, U_old, theta, DimTotal);
         }
         free(P1); free(P2); free(U_old);
     }
     else {
-        /*3D case*/
+          /*3D case*/
+        float *U_old=NULL, *P1=NULL, *P2=NULL, *P3=NULL;
+        U_old = calloc(DimTotal, sizeof(float));
+        P1 = calloc(DimTotal, sizeof(float));
+        P2 = calloc(DimTotal, sizeof(float));
+        P3 = calloc(DimTotal, sizeof(float));
+
+        /* begin iterations */
+        for(ll=0; ll<iterationsNumb; ll++) {
+
+         /* computing the the dual P variable */
+            DualP3D(U, P1, P2, P3, (long)(dimX), (long)(dimY),  (long)(dimZ), sigma);
+
+            /* apply nonnegativity */
+            if (nonneg == 1) for(j=0; j<DimTotal; j++) {if (U[j] < 0.0f) U[j] = 0.0f;}
+
+            /* projection step */
+            Proj_func3D(P1, P2, P3, methodTV, DimTotal);
+
+            /* copy U to U_old */
+            copyIm(U, U_old, (long)(dimX), (long)(dimY), (long)(dimZ));
+
+            DivProj3D(U, Input, P1, P2, P3, (long)(dimX), (long)(dimY), (long)(dimZ), 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*/
+
+            getX(U, U_old, theta, DimTotal);
+        }
+        free(P1); free(P2); free(P3); free(U_old);
     }
     /*adding info into info_vector */
     infovector[0] = (float)(ll);  /*iterations number (if stopped earlier based on tolerance)*/
@@ -134,7 +173,7 @@ float DivProj2D(float *U, float *Input, float *P1, float *P2, long dimX, long di
 {
   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)
+  #pragma omp parallel for shared(U,Input,P1,P2) private(index, i, j, P_v1, P_v2, div_var)
   for(j=0; j<dimY; j++) {
     for(i=0; i<dimX; i++) {
             index = j*dimX+i;
@@ -150,7 +189,7 @@ float DivProj2D(float *U, float *Input, float *P1, float *P2, long dimX, long di
 }
 
 /*get the updated solution*/
-float getX2D(float *U, float *U_old, float theta, long DimTotal)
+float getX(float *U, float *U_old, float theta, long DimTotal)
 {
     long i;
     #pragma omp parallel for shared(U,U_old) private(i)
@@ -164,3 +203,45 @@ float getX2D(float *U, float *U_old, float theta, long DimTotal)
 /*****************************************************************/
 /************************3D-case related Functions */
 /*****************************************************************/
+/*Calculating dual variable (using forward differences)*/
+float DualP3D(float *U, float *P1, float *P2, float *P3, long dimX, long dimY, long dimZ, float sigma)
+{
+     long i,j,k,index;
+     #pragma omp parallel for shared(U,P1,P2,P3) private(index,i,j,k)
+     for(k=0; k<dimZ; k++) {
+         for(j=0; j<dimY; j++) {
+           for(i=0; i<dimX; i++) {
+          index = (dimX*dimY)*k + j*dimX+i;
+          /* symmetric boundary conditions (Neuman) */
+          if (i == dimX-1) P1[index] += sigma*(U[(dimX*dimY)*k + j*dimX+(i-1)] - U[index]);
+          else P1[index] += sigma*(U[(dimX*dimY)*k + j*dimX+(i+1)] - U[index]);
+          if (j == dimY-1) P2[index] += sigma*(U[(dimX*dimY)*k + (j-1)*dimX+i] - U[index]);
+          else  P2[index] += sigma*(U[(dimX*dimY)*k + (j+1)*dimX+i] - U[index]);
+          if (k == dimZ-1) P3[index] += sigma*(U[(dimX*dimY)*(k-1) + j*dimX+i] - U[index]);
+          else  P3[index] += sigma*(U[(dimX*dimY)*(k+1) + j*dimX+i] - U[index]);
+        }}}
+     return 1;
+}
+
+/* Divergence for P dual */
+float DivProj3D(float *U, float *Input, float *P1, float *P2, float *P3, long dimX, long dimY, long dimZ, float lt, float tau)
+{
+  long i,j,k,index;
+  float P_v1, P_v2, P_v3, div_var;
+  #pragma omp parallel for shared(U,Input,P1,P2) private(index, i, j, k, P_v1, P_v2, P_v3, div_var)
+  for(k=0; k<dimZ; k++) {
+      for(j=0; j<dimY; j++) {
+        for(i=0; i<dimX; i++) {
+            index = (dimX*dimY)*k + j*dimX+i;
+            /* symmetric boundary conditions (Neuman) */
+            if (i == 0) P_v1 = -P1[index];
+            else P_v1 = -(P1[index] - P1[(dimX*dimY)*k + j*dimX+(i-1)]);
+            if (j == 0) P_v2 = -P2[index];
+            else  P_v2 = -(P2[index] - P2[(dimX*dimY)*k + (j-1)*dimX+i]);
+            if (k == 0) P_v3 = -P3[index];
+            else  P_v3 = -(P3[index] - P3[(dimX*dimY)*(k-1) + j*dimX+i]);
+            div_var = P_v1 + P_v2 + P_v3;
+            U[index] = (U[index] - tau*div_var + lt*Input[index])/(1.0 + lt);
+   }}}
+  return *U;
+}
diff --git a/src/Core/regularisers_CPU/PD_TV_core.h b/src/Core/regularisers_CPU/PD_TV_core.h
index b52689a..b4e8a75 100644
--- a/src/Core/regularisers_CPU/PD_TV_core.h
+++ b/src/Core/regularisers_CPU/PD_TV_core.h
@@ -51,7 +51,10 @@ CCPI_EXPORT float PDTV_CPU_main(float *Input, float *U, float *infovector, float
 
 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);
+CCPI_EXPORT float getX(float *U, float *U_old, float theta, long DimTotal);
+
+CCPI_EXPORT float DualP3D(float *U, float *P1, float *P2, float *P3, long dimX, long dimY, long dimZ, float sigma);
+CCPI_EXPORT float DivProj3D(float *U, float *Input, float *P1, float *P2, float *P3, long dimX, long dimY, long dimZ, float lt, float tau);
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/Core/regularisers_CPU/utils.c b/src/Core/regularisers_CPU/utils.c
index cf4ad72..088ace9 100644
--- a/src/Core/regularisers_CPU/utils.c
+++ b/src/Core/regularisers_CPU/utils.c
@@ -145,7 +145,7 @@ float Im_scale2D(float *Input, float *Scaled, int w, int h, int w2, int h2)
     return *Scaled;
 }
 
-/*2D Projection onto convex set for P*/
+/*2D Projection onto convex set for P (called in PD_TV, FGP_dTV and FGP_TV methods)*/
 float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal)
 {
     float val1, val2, denom, sq_denom;
@@ -176,3 +176,38 @@ float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal)
     }
     return 1;
 }
+/*3D Projection onto convex set for P (called in PD_TV, FGP_TV, FGP_dTV methods)*/
+float Proj_func3D(float *P1, float *P2, float *P3, int methTV, long DimTotal)
+{
+    float val1, val2, val3, denom, sq_denom;
+    long i;
+    if (methTV == 0) {
+        /* isotropic TV*/
+#pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3,sq_denom)
+        for(i=0; i<DimTotal; i++) {
+            denom = powf(P1[i],2) + powf(P2[i],2) + powf(P3[i],2);
+            if (denom > 1.0f) {
+                sq_denom = 1.0f/sqrtf(denom);
+                P1[i] = P1[i]*sq_denom;
+                P2[i] = P2[i]*sq_denom;
+                P3[i] = P3[i]*sq_denom;
+            }
+        }
+    }
+    else {
+        /* anisotropic TV*/
+#pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3)
+        for(i=0; i<DimTotal; i++) {
+            val1 = fabs(P1[i]);
+            val2 = fabs(P2[i]);
+            val3 = fabs(P3[i]);
+            if (val1 < 1.0f) {val1 = 1.0f;}
+            if (val2 < 1.0f) {val2 = 1.0f;}
+            if (val3 < 1.0f) {val3 = 1.0f;}
+            P1[i] = P1[i]/val1;
+            P2[i] = P2[i]/val2;
+            P3[i] = P3[i]/val3;
+        }
+    }
+    return 1;
+}
diff --git a/src/Core/regularisers_CPU/utils.h b/src/Core/regularisers_CPU/utils.h
index e050f86..4b57f71 100644
--- a/src/Core/regularisers_CPU/utils.h
+++ b/src/Core/regularisers_CPU/utils.h
@@ -32,6 +32,7 @@ CCPI_EXPORT float TV_energy3D(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 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);
+CCPI_EXPORT float Proj_func3D(float *P1, float *P2, float *P3, int methTV, long DimTotal);
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/Python/setup-regularisers.py.in b/src/Python/setup-regularisers.py.in
index 9bcd46d..9a5b693 100644
--- a/src/Python/setup-regularisers.py.in
+++ b/src/Python/setup-regularisers.py.in
@@ -38,7 +38,6 @@ 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" ) ,
diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx
index 724634b..08e247c 100644
--- a/src/Python/src/cpu_regularisers.pyx
+++ b/src/Python/src/cpu_regularisers.pyx
@@ -163,7 +163,7 @@ def TV_PD_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_par
     if inputData.ndim == 2:
         return TV_PD_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const)
     elif inputData.ndim == 3:
-        return 0
+        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,
@@ -191,7 +191,34 @@ def TV_PD_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                        methodTV,
                        nonneg,
                        dims[1],dims[0], 1)
+    return (outputData,infovec)
+
+def TV_PD_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+                     float regularisation_parameter,
+                     int iterationsNumb,
+                     float tolerance_param,
+                     int methodTV,
+                     int nonneg,
+                     float lipschitz_const):
+
+    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 FGP-TV iterations for 3D data */
+    PDTV_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter,
+                       iterationsNumb,
+                       tolerance_param,
+                       lipschitz_const,
+                       methodTV,
+                       nonneg,
+                       dims[2], dims[1], dims[0])
     return (outputData,infovec)
 
 #***************************************************************#
-- 
cgit v1.2.3