Merge pull request #8 from vais-ral/pythonize

Pythonize

17 files changed, 249 insertions, 59 deletions

diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m
index fa98360..3d22b97 100644
--- a/main_func/FISTA_REC.m
+++ b/main_func/FISTA_REC.m
@@ -517,6 +517,18 @@ else
         
         % subsets loop
         counterInd = 1;
+		if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec'))
+			% if geometry is 2D use slice-by-slice projection-backprojection routine                    
+			for kkk = 1:SlicesZ
+				[sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom);
+				sino_updt_Sub(:,:,kkk) = sinoT';
+				astra_mex_data2d('delete', sino_id);
+			end
+		else
+			% for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8)
+			[sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom);
+			astra_mex_data3d('delete', sino_id);
+		end
         for ss = 1:subsets
             X_old = X;
             t_old = t;
@@ -541,6 +553,7 @@ else
             
             if (lambdaR_L1 > 0)
                 % Group-Huber fidelity (ring removal)
+
                 
                 residualSub = zeros(Detectors, numProjSub, SlicesZ,'single'); % residual for a chosen subset
                 for kkk = 1:numProjSub
@@ -551,6 +564,7 @@ else
                 
             elseif (studentt > 0)
                 % student t data fidelity
+
                 % artifacts removal with Students t penalty
                 residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:)));
                 
@@ -564,9 +578,11 @@ else
                 objective(i) = ff; % for the objective function output
             else
                 % PWLS model
+
                 residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:)));
                 objective(i) = 0.5*norm(residualSub(:)); % for the objective function output
             end
+
             
             % perform backprojection of a subset
             if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec'))
diff --git a/main_func/regularizers_CPU/FGP_TV.c b/main_func/regularizers_CPU/FGP_TV.c
index e0dae57..30cea1a 100644
--- a/main_func/regularizers_CPU/FGP_TV.c
+++ b/main_func/regularizers_CPU/FGP_TV.c
@@ -3,7 +3,7 @@ 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 2017 Daniil Kazanteev
+Copyright 2017 Daniil Kazantsev
 Copyright 2017 Srikanth Nagella, Edoardo Pasca
 
 Licensed under the Apache License, Version 2.0 (the "License");
@@ -16,6 +16,7 @@ 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 "FGP_TV_core.h"
 
diff --git a/main_func/regularizers_CPU/FGP_TV_core.c b/main_func/regularizers_CPU/FGP_TV_core.c
index 9cde327..03cd445 100644
--- a/main_func/regularizers_CPU/FGP_TV_core.c
+++ b/main_func/regularizers_CPU/FGP_TV_core.c
@@ -3,7 +3,7 @@ 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 2017 Daniil Kazanteev
+Copyright 2017 Daniil Kazantsev
 Copyright 2017 Srikanth Nagella, Edoardo Pasca
 
 Licensed under the Apache License, Version 2.0 (the "License");
@@ -263,13 +263,4 @@ float Rupd_func3D(float *P1, float *P1_old, float *P2, float *P2_old, float *P3,
 	return 1;
 }
 
-/* General Functions */
-/*****************************************************************/
-/* Copy Image */
-float copyIm(float *A, float *B, int dimX, int dimY, int dimZ)
-{
-	int j;
-#pragma omp parallel for shared(A, B) private(j)
-	for (j = 0; j<dimX*dimY*dimZ; j++)  B[j] = A[j];
-	return *B;
-}
+
diff --git a/main_func/regularizers_CPU/FGP_TV_core.h b/main_func/regularizers_CPU/FGP_TV_core.h
index 697fd84..6430bf2 100644
--- a/main_func/regularizers_CPU/FGP_TV_core.h
+++ b/main_func/regularizers_CPU/FGP_TV_core.h
@@ -3,7 +3,7 @@ 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 2017 Daniil Kazanteev
+Copyright 2017 Daniil Kazantsev
 Copyright 2017 Srikanth Nagella, Edoardo Pasca
 
 Licensed under the Apache License, Version 2.0 (the "License");
@@ -17,14 +17,44 @@ See the License for the specific language governing permissions and
 limitations under the License.
 */
 
-#include <matrix.h>
+//#include <matrix.h>
 #include <math.h>
 #include <stdlib.h>
 #include <memory.h>
 #include <stdio.h>
 #include "omp.h"
+#include "utils.h"
 
-float copyIm(float *A, float *B, int dimX, int dimY, int dimZ);
+/* C-OMP implementation of FGP-TV [1] denoising/regularization model (2D/3D case)
+*
+* Input Parameters:
+* 1. Noisy image/volume [REQUIRED]
+* 2. lambda - regularization parameter [REQUIRED]
+* 3. Number of iterations [OPTIONAL parameter]
+* 4. eplsilon: tolerance constant [OPTIONAL parameter]
+* 5. TV-type: 'iso' or 'l1' [OPTIONAL parameter]
+*
+* Output:
+* [1] Filtered/regularized image
+* [2] last function value
+*
+* Example of image denoising:
+* figure;
+* Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+* u0 = Im + .05*randn(size(Im)); % adding noise
+* u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+*
+* to compile with OMP support: mex FGP_TV.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+* 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"
+*
+* D. Kazantsev, 2016-17
+*
+*/
+#ifdef __cplusplus
+extern "C" {
+#endif
+//float copyIm(float *A, float *B, int dimX, int dimY, int dimZ);
 float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, int dimX, int dimY);
 float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float lambda, int dimX, int dimY);
 float Proj_func2D(float *P1, float *P2, int methTV, int dimX, int dimY);
@@ -36,3 +66,6 @@ float Grad_func3D(float *P1, float *P2, float *P3, float *D, float *R1, float *R
 float Proj_func3D(float *P1, float *P2, float *P3, int dimX, int dimY, int dimZ);
 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, int dimX, int dimY, int dimZ);
 float Obj_func_CALC3D(float *A, float *D, float *funcvalA, float lambda, int dimX, int dimY, int dimZ);
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file
diff --git a/main_func/regularizers_CPU/LLT_model.c b/main_func/regularizers_CPU/LLT_model.c
index 47146a5..0b07b47 100644
--- a/main_func/regularizers_CPU/LLT_model.c
+++ b/main_func/regularizers_CPU/LLT_model.c
@@ -3,7 +3,7 @@ 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 2017 Daniil Kazanteev
+Copyright 2017 Daniil Kazantsev
 Copyright 2017 Srikanth Nagella, Edoardo Pasca
 
 Licensed under the Apache License, Version 2.0 (the "License");
@@ -18,12 +18,13 @@ limitations under the License.
 */
 
 #include "mex.h"
+#include "matrix.h"
 #include "LLT_model_core.h"
 
 /* C-OMP implementation of Lysaker, Lundervold and Tai (LLT) model of higher order regularization penalty
 *
 * Input Parameters:
-* 1. U0 - origanal noise image/volume
+* 1. U0 - original noise image/volume
 * 2. lambda - regularization parameter
 * 3. tau - time-step  for explicit scheme
 * 4. iter - iterations number
@@ -46,7 +47,6 @@ limitations under the License.
 * 28.11.16/Harwell
 */
 
-
 void mexFunction(
         int nlhs, mxArray *plhs[],
         int nrhs, const mxArray *prhs[])
diff --git a/main_func/regularizers_CPU/LLT_model_core.c b/main_func/regularizers_CPU/LLT_model_core.c
index 7a1cdbe..3a853d2 100644
--- a/main_func/regularizers_CPU/LLT_model_core.c
+++ b/main_func/regularizers_CPU/LLT_model_core.c
@@ -3,7 +3,7 @@ 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 2017 Daniil Kazanteev
+Copyright 2017 Daniil Kazantsev
 Copyright 2017 Srikanth Nagella, Edoardo Pasca
 
 Licensed under the Apache License, Version 2.0 (the "License");
@@ -314,12 +314,5 @@ float cleanMap(unsigned short *Map, int dimX, int dimY, int dimZ)
 	return *Map;
 }
 
-/* Copy Image */
-float copyIm(float *A, float *U, int dimX, int dimY, int dimZ)
-{
-	int j;
-#pragma omp parallel for shared(A, U) private(j)
-	for (j = 0; j<dimX*dimY*dimZ; j++)  U[j] = A[j];
-	return *U;
-}
+
 /*********************3D *********************/
\ No newline at end of file
diff --git a/main_func/regularizers_CPU/LLT_model_core.h b/main_func/regularizers_CPU/LLT_model_core.h
index 10f52fe..13fce5a 100644
--- a/main_func/regularizers_CPU/LLT_model_core.h
+++ b/main_func/regularizers_CPU/LLT_model_core.h
@@ -3,7 +3,7 @@ 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 2017 Daniil Kazanteev
+Copyright 2017 Daniil Kazantsev
 Copyright 2017 Srikanth Nagella, Edoardo Pasca
 
 Licensed under the Apache License, Version 2.0 (the "License");
@@ -17,16 +17,20 @@ See the License for the specific language governing permissions and
 limitations under the License.
 */
 
-#include <matrix.h>
+//#include <matrix.h>
 #include <math.h>
 #include <stdlib.h>
 #include <memory.h>
 #include <stdio.h>
 #include "omp.h"
+#include "utils.h"
 
 #define EPS 0.01
 
 /* 2D functions */
+#ifdef __cplusplus
+extern "C" {
+#endif
 float der2D(float *U, float *D1, float *D2, int dimX, int dimY, int dimZ);
 float div_upd2D(float *U0, float *U, float *D1, float *D2, int dimX, int dimY, int dimZ, float lambda, float tau);
 
@@ -36,4 +40,7 @@ float div_upd3D(float *U0, float *U, float *D1, float *D2, float *D3, unsigned s
 float calcMap(float *U, unsigned short *Map, int dimX, int dimY, int dimZ);
 float cleanMap(unsigned short *Map, int dimX, int dimY, int dimZ);
 
-float copyIm(float *A, float *U, int dimX, int dimY, int dimZ);
+//float copyIm(float *A, float *U, int dimX, int dimY, int dimZ);
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file
diff --git a/main_func/regularizers_CPU/PatchBased_Regul.c b/main_func/regularizers_CPU/PatchBased_Regul.c
index 59eb3b3..9c925df 100644
--- a/main_func/regularizers_CPU/PatchBased_Regul.c
+++ b/main_func/regularizers_CPU/PatchBased_Regul.c
@@ -3,7 +3,7 @@ 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 2017 Daniil Kazanteev

+Copyright 2017 Daniil Kazantsev

 Copyright 2017 Srikanth Nagella, Edoardo Pasca

 

 Licensed under the Apache License, Version 2.0 (the "License");

@@ -18,6 +18,7 @@ limitations under the License.
 */

 

 #include "mex.h"

+#include "matrix.h"

 #include "PatchBased_Regul_core.h"

 

 

diff --git a/main_func/regularizers_CPU/PatchBased_Regul_core.h b/main_func/regularizers_CPU/PatchBased_Regul_core.h
index 5aa6415..d4a8a46 100644
--- a/main_func/regularizers_CPU/PatchBased_Regul_core.h
+++ b/main_func/regularizers_CPU/PatchBased_Regul_core.h
@@ -19,13 +19,51 @@ limitations under the License.
 
 #define _USE_MATH_DEFINES
 
-#include <matrix.h>
+//#include <matrix.h>
 #include <math.h>
 #include <stdlib.h>
 #include <memory.h>
 #include <stdio.h>
 #include "omp.h"
 
+/* C-OMP implementation of  patch-based (PB) regularization (2D and 3D cases).
+* This method finds self-similar patches in data and performs one fixed point iteration to mimimize the PB penalty function
+*
+* References: 1. Yang Z. & Jacob M. "Nonlocal Regularization of Inverse Problems"
+*             2. Kazantsev D. et al. "4D-CT reconstruction with unified spatial-temporal patch-based regularization"
+*
+* Input Parameters (mandatory):
+* 1. Image (2D or 3D)
+* 2. ratio of the searching window (e.g. 3 = (2*3+1) = 7 pixels window)
+* 3. ratio of the similarity window (e.g. 1 = (2*1+1) = 3 pixels window)
+* 4. h - parameter for the PB penalty function
+* 5. lambda - regularization parameter
+
+* Output:
+* 1. regularized (denoised) Image (N x N)/volume (N x N x N)
+*
+* Quick 2D denoising example in Matlab:
+Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
+*
+* Please see more tests in a file:
+TestTemporalSmoothing.m
+
+*
+* Matlab + C/mex compilers needed
+* to compile with OMP support: mex PB_Regul_CPU.c CFLAGS="\$CFLAGS -fopenmp -Wall" LDFLAGS="\$LDFLAGS -fopenmp"
+*
+* D. Kazantsev *
+* 02/07/2014
+* Harwell, UK
+*/
+#ifdef __cplusplus
+extern "C" {
+#endif
 float pad_crop(float *A, float *Ap, int OldSizeX, int OldSizeY, int OldSizeZ, int NewSizeX, int NewSizeY, int NewSizeZ, int padXY, int switchpad_crop);
 float PB_FUNC2D(float *A, float *B, int dimX, int dimY, int padXY, int SearchW, int SimilW, float h, float lambda);
 float PB_FUNC3D(float *A, float *B, int dimX, int dimY, int dimZ, int padXY, int SearchW, int SimilW, float h, float lambda);
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file
diff --git a/main_func/regularizers_CPU/SplitBregman_TV.c b/main_func/regularizers_CPU/SplitBregman_TV.c
index 0dc638d..38f6a9d 100644
--- a/main_func/regularizers_CPU/SplitBregman_TV.c
+++ b/main_func/regularizers_CPU/SplitBregman_TV.c
@@ -3,7 +3,7 @@ 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 2017 Daniil Kazanteev
+Copyright 2017 Daniil Kazantsev
 Copyright 2017 Srikanth Nagella, Edoardo Pasca
 
 Licensed under the Apache License, Version 2.0 (the "License");
@@ -18,6 +18,7 @@ limitations under the License.
 */
 
 #include "mex.h"
+#include <matrix.h>
 #include "SplitBregman_TV_core.h"
 
 /* C-OMP implementation of Split Bregman - TV denoising-regularization model (2D/3D)
diff --git a/main_func/regularizers_CPU/SplitBregman_TV_core.c b/main_func/regularizers_CPU/SplitBregman_TV_core.c
index 283dd43..4109a4b 100644
--- a/main_func/regularizers_CPU/SplitBregman_TV_core.c
+++ b/main_func/regularizers_CPU/SplitBregman_TV_core.c
@@ -3,7 +3,7 @@ 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 2017 Daniil Kazanteev
+Copyright 2017 Daniil Kazantsev
 Copyright 2017 Srikanth Nagella, Edoardo Pasca
 
 Licensed under the Apache License, Version 2.0 (the "License");
@@ -257,13 +257,3 @@ float updBxByBz3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *B
             }}}
     return 1;
 }
-/* General Functions */
-/*****************************************************************/
-/* Copy Image */
-float copyIm(float *A, float *B, int dimX, int dimY, int dimZ)
-{
-    int j;
-#pragma omp parallel for shared(A, B) private(j)
-    for(j=0; j<dimX*dimY*dimZ; j++)  B[j] = A[j];
-    return *B;
-}
\ No newline at end of file
diff --git a/main_func/regularizers_CPU/SplitBregman_TV_core.h b/main_func/regularizers_CPU/SplitBregman_TV_core.h
index a7aaabb..6ed3ff9 100644
--- a/main_func/regularizers_CPU/SplitBregman_TV_core.h
+++ b/main_func/regularizers_CPU/SplitBregman_TV_core.h
@@ -3,7 +3,7 @@ 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 2017 Daniil Kazanteev
+Copyright 2017 Daniil Kazantsev
 Copyright 2017 Srikanth Nagella, Edoardo Pasca
 
 Licensed under the Apache License, Version 2.0 (the "License");
@@ -16,14 +16,44 @@ 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 <matrix.h>
 #include <math.h>
 #include <stdlib.h>
 #include <memory.h>
 #include <stdio.h>
 #include "omp.h"
 
-float copyIm(float *A, float *B, int dimX, int dimY, int dimZ);
+#include "utils.h"
+
+/* C-OMP implementation of Split Bregman - TV denoising-regularization model (2D/3D)
+*
+* Input Parameters:
+* 1. Noisy image/volume
+* 2. lambda - regularization parameter
+* 3. Number of iterations [OPTIONAL parameter]
+* 4. eplsilon - tolerance constant [OPTIONAL parameter]
+* 5. TV-type: 'iso' or 'l1' [OPTIONAL parameter]
+*
+* Output:
+* Filtered/regularized image
+*
+* Example:
+* figure;
+* Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+* u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
+* u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+*
+* to compile with OMP support: mex SplitBregman_TV.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+* References:
+* The Split Bregman Method for L1 Regularized Problems, by Tom Goldstein and Stanley Osher.
+* D. Kazantsev, 2016*
+*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+//float copyIm(float *A, float *B, int dimX, int dimY, int dimZ);
 float gauss_seidel2D(float *U, float *A, float *Dx, float *Dy, float *Bx, float *By, int dimX, int dimY, float lambda, float mu);
 float updDxDy_shrinkAniso2D(float *U, float *Dx, float *Dy, float *Bx, float *By, int dimX, int dimY, float lambda);
 float updDxDy_shrinkIso2D(float *U, float *Dx, float *Dy, float *Bx, float *By, int dimX, int dimY, float lambda);
@@ -33,3 +63,7 @@ float gauss_seidel3D(float *U, float *A, float *Dx, float *Dy, float *Dz, float
 float updDxDyDz_shrinkAniso3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, int dimX, int dimY, int dimZ, float lambda);
 float updDxDyDz_shrinkIso3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, int dimX, int dimY, int dimZ, float lambda);
 float updBxByBz3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, int dimX, int dimY, int dimZ);
+
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file
diff --git a/main_func/regularizers_CPU/TGV_PD.c b/main_func/regularizers_CPU/TGV_PD.c
index 6593d8e..c9cb440 100644
--- a/main_func/regularizers_CPU/TGV_PD.c
+++ b/main_func/regularizers_CPU/TGV_PD.c
@@ -3,7 +3,7 @@ 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 2017 Daniil Kazanteev
+Copyright 2017 Daniil Kazantsev
 Copyright 2017 Srikanth Nagella, Edoardo Pasca
 
 Licensed under the Apache License, Version 2.0 (the "License");
diff --git a/main_func/regularizers_CPU/TGV_PD_core.c b/main_func/regularizers_CPU/TGV_PD_core.c
index 1164b73..4139d10 100644
--- a/main_func/regularizers_CPU/TGV_PD_core.c
+++ b/main_func/regularizers_CPU/TGV_PD_core.c
@@ -186,14 +186,6 @@ float UpdV_2D(float *V1, float *V2, float *P1, float *P2, float *Q1, float *Q2,
         }}
     return 1;
 }
-/* Copy Image */
-float copyIm(float *A, float *U, int dimX, int dimY, int dimZ)
-{
-    int j;
-#pragma omp parallel for shared(A, U) private(j)
-    for(j=0; j<dimX*dimY*dimZ; j++)  U[j] = A[j];
-    return *U;
-}
 /*********************3D *********************/
 
 /*Calculating dual variable P (using forward differences)*/
diff --git a/main_func/regularizers_CPU/TGV_PD_core.h b/main_func/regularizers_CPU/TGV_PD_core.h
index 04ba95c..d5378df 100644
--- a/main_func/regularizers_CPU/TGV_PD_core.h
+++ b/main_func/regularizers_CPU/TGV_PD_core.h
@@ -3,7 +3,7 @@ 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 2017 Daniil Kazanteev
+Copyright 2017 Daniil Kazantsev
 Copyright 2017 Srikanth Nagella, Edoardo Pasca
 
 Licensed under the Apache License, Version 2.0 (the "License");
@@ -17,13 +17,42 @@ See the License for the specific language governing permissions and
 limitations under the License.
 */
 
-#include <matrix.h>
+//#include <matrix.h>
 #include <math.h>
 #include <stdlib.h>
 #include <memory.h>
 #include <stdio.h>
 #include "omp.h"
+#include "utils.h"
 
+/* C-OMP implementation of Primal-Dual denoising method for
+* Total Generilized Variation (TGV)-L2 model (2D case only)
+*
+* Input Parameters:
+* 1. Noisy image/volume (2D)
+* 2. lambda - regularization parameter
+* 3. parameter to control first-order term (alpha1)
+* 4. parameter to control the second-order term (alpha0)
+* 5. Number of CP iterations
+*
+* Output:
+* Filtered/regularized image
+*
+* Example:
+* figure;
+* Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+* u0 = Im + .03*randn(size(Im)); % adding noise
+* tic; u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); toc;
+*
+* to compile with OMP support: mex TGV_PD.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+* References:
+* K. Bredies "Total Generalized Variation"
+*
+* 28.11.16/Harwell
+*/
+#ifdef __cplusplus
+extern "C" {
+#endif
 /* 2D functions */
 float DualP_2D(float *U, float *V1, float *V2, float *P1, float *P2, int dimX, int dimY, int dimZ, float sigma);
 float ProjP_2D(float *P1, float *P2, int dimX, int dimY, int dimZ, float alpha1);
@@ -32,4 +61,7 @@ float ProjQ_2D(float *Q1, float *Q2, float *Q3, int dimX, int dimY, int dimZ, fl
 float DivProjP_2D(float *U, float *A, float *P1, float *P2, int dimX, int dimY, int dimZ, float lambda, float tau);
 float UpdV_2D(float *V1, float *V2, float *P1, float *P2, float *Q1, float *Q2, float *Q3, int dimX, int dimY, int dimZ, float tau);
 float newU(float *U, float *U_old, int dimX, int dimY, int dimZ);
-float copyIm(float *A, float *U, int dimX, int dimY, int dimZ);
+//float copyIm(float *A, float *U, int dimX, int dimY, int dimZ);
+#ifdef __cplusplus
+}
+#endif
diff --git a/main_func/regularizers_CPU/utils.c b/main_func/regularizers_CPU/utils.c
new file mode 100644
index 0000000..0e83d2c
--- /dev/null
+++ b/main_func/regularizers_CPU/utils.c
@@ -0,0 +1,29 @@
+/*
+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 2017 Daniil Kazanteev
+Copyright 2017 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 "utils.h"
+
+/* Copy Image */
+float copyIm(float *A, float *U, int dimX, int dimY, int dimZ)
+{
+	int j;
+#pragma omp parallel for shared(A, U) private(j)
+	for (j = 0; j<dimX*dimY*dimZ; j++)  U[j] = A[j];
+	return *U;
+}
\ No newline at end of file
diff --git a/main_func/regularizers_CPU/utils.h b/main_func/regularizers_CPU/utils.h
new file mode 100644
index 0000000..53463a3
--- /dev/null
+++ b/main_func/regularizers_CPU/utils.h
@@ -0,0 +1,32 @@
+/*
+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 2017 Daniil Kazantsev
+Copyright 2017 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"
+#ifdef __cplusplus
+extern "C" {
+#endif
+float copyIm(float *A, float *U, int dimX, int dimY, int dimZ);
+#ifdef __cplusplus
+}
+#endif