Merge pull request #40 from vais-ral/GPU_test

ROF/FGP updated via Cython

8 files changed, 588 insertions, 397 deletions

diff --git a/Wrappers/Python/ccpi/filters/regularizers.py b/Wrappers/Python/ccpi/filters/regularizers.py
new file mode 100644
index 0000000..d6dfa8c
--- /dev/null
+++ b/Wrappers/Python/ccpi/filters/regularizers.py
@@ -0,0 +1,44 @@
+"""
+script which assigns a proper device core function based on a flag ('cpu' or 'gpu')
+"""
+
+from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU, TV_FGP_CPU
+from ccpi.filters.gpu_regularizers import TV_ROF_GPU, TV_FGP_GPU
+
+def ROF_TV(inputData, regularization_parameter, iterations,
+                     time_marching_parameter,device='cpu'):
+    if device == 'cpu':
+        return TV_ROF_CPU(inputData,
+                     regularization_parameter,
+                     iterations, 
+                     time_marching_parameter)
+    elif device == 'gpu':
+        return TV_ROF_GPU(inputData,
+                     regularization_parameter,
+                     iterations, 
+                     time_marching_parameter)
+    else:
+        raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
+                         .format(device))
+
+def FGP_TV(inputData, regularization_parameter,iterations,
+                     tolerance_param, methodTV, nonneg, printM, device='cpu'):
+    if device == 'cpu':
+        return TV_FGP_CPU(inputData,
+                     regularization_parameter,
+                     iterations, 
+                     tolerance_param,
+                     methodTV,
+                     nonneg,
+                     printM)
+    elif device == 'gpu':
+        return TV_FGP_GPU(inputData,
+                     regularization_parameter,
+                     iterations, 
+                     tolerance_param,
+                     methodTV,
+                     nonneg,
+                     printM)
+    else:
+        raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
+                         .format(device))
\ No newline at end of file
diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py
index 5908c3c..76b9de7 100644
--- a/Wrappers/Python/demo/test_cpu_regularizers.py
+++ b/Wrappers/Python/demo/test_cpu_regularizers.py
@@ -11,11 +11,8 @@ import numpy as np
 import os    
 from enum import Enum
 import timeit
-from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV , FGP_TV ,\
-                                                 LLT_model, PatchBased_Regul ,\
-                                                 TGV_PD
-from ccpi.filters.cpu_regularizers_cython import ROF_TV
-
+from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV, LLT_model, PatchBased_Regul, TGV_PD
+from ccpi.filters.regularizers import ROF_TV, FGP_TV
 ###############################################################################
 #https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956
 #NRMSE a normalization of the root of the mean squared error
@@ -112,12 +109,10 @@ rms = rmse(Im, splitbregman)
 pars['rmse'] = rms
 txtstr = printParametersToString(pars) 
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-    
+print (txtstr)   
 
 a=fig.add_subplot(2,4,2)
 
-
 # these are matplotlib.patch.Patch properties
 props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
 # place a text box in upper left in axes coords
@@ -128,23 +123,26 @@ imgplot = plt.imshow(splitbregman,\
                      )
 
 ###################### FGP_TV #########################################
-# u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+
 start_time = timeit.default_timer()
 pars = {'algorithm' : FGP_TV , \
-        'input' : u0,
-        'regularization_parameter':0.05, \
-        'number_of_iterations' :200 ,\
-        'tolerance_constant':1e-4,\
-        'TV_penalty': 0
-}
+        'input' : u0,\
+        'regularization_parameter':0.07, \
+        'number_of_iterations' :300 ,\
+        'tolerance_constant':0.00001,\
+        'methodTV': 0 ,\
+        'nonneg': 0 ,\
+        'printingOut': 0 
+        }
 
-out = FGP_TV (pars['input'], 
+fgp = FGP_TV(pars['input'], 
               pars['regularization_parameter'],
               pars['number_of_iterations'],
               pars['tolerance_constant'], 
-              pars['TV_penalty'])  
+              pars['methodTV'],
+              pars['nonneg'],
+              pars['printingOut'], 'cpu')  
 
-fgp = out[0]
 rms = rmse(Im, fgp)
 pars['rmse'] = rms
 
@@ -165,8 +163,8 @@ imgplot = plt.imshow(fgp, \
 a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
         verticalalignment='top', bbox=props)
 
-###################### LLT_model #########################################
 
+###################### LLT_model #########################################
 start_time = timeit.default_timer()
 
 pars = {'algorithm': LLT_model , \
@@ -201,8 +199,6 @@ a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
 imgplot = plt.imshow(llt,\
                      cmap="gray"
                      )
-
-
 # ###################### PatchBased_Regul #########################################
 # # Quick 2D denoising example in Matlab:   
 # #   Im = double(imread('lena_gray_256.tif'))/255;  % loading image
@@ -284,16 +280,15 @@ start_time = timeit.default_timer()
 
 pars = {'algorithm': ROF_TV , \
         'input' : u0,\
-        'regularization_parameter':1,\
-        'marching_step': 0.003,\
+        'regularization_parameter':0.07,\
+        'marching_step': 0.0025,\
         'number_of_iterations': 300
         }
 rof = ROF_TV(pars['input'],
-             pars['number_of_iterations'],
-             pars['regularization_parameter'],
-             pars['marching_step'] 
-             )
-#tgv = out
+			 pars['regularization_parameter'],
+             pars['number_of_iterations'],             
+             pars['marching_step'], 'cpu')
+
 rms = rmse(Im, rof)
 pars['rmse'] = rms
 
@@ -307,9 +302,7 @@ props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
 # place a text box in upper left in axes coords
 a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
          verticalalignment='top', bbox=props)
-imgplot = plt.imshow(tgv, cmap="gray")
-
-
+imgplot = plt.imshow(rof, cmap="gray")
 
 plt.show()
 
diff --git a/Wrappers/Python/setup-regularizers.py.in b/Wrappers/Python/setup-regularizers.py.in
index a125261..8655a2e 100644
--- a/Wrappers/Python/setup-regularizers.py.in
+++ b/Wrappers/Python/setup-regularizers.py.in
@@ -34,7 +34,9 @@ extra_libraries = ['cilreg']
 
 extra_include_dirs += [os.path.join(".." , ".." , "Core"),
                        os.path.join(".." , ".." , "Core",  "regularizers_CPU"),
-                       os.path.join(".." , ".." , "Core",  "regularizers_GPU") , 
+                       os.path.join(".." , ".." , "Core",  "regularizers_GPU" , "Diffus_HO" ) , 
+                       os.path.join(".." , ".." , "Core",  "regularizers_GPU" , "NL_Regul" ) , 
+                       os.path.join(".." , ".." , "Core",  "regularizers_GPU" , "TV_ROF" ) , 
 						   "."]
 
 if platform.system() == 'Windows':				   
@@ -81,4 +83,4 @@ setup(
 )
 
 
-@SETUP_GPU_WRAPPERS@
\ No newline at end of file
+@SETUP_GPU_WRAPPERS@
diff --git a/Wrappers/Python/src/cpu_regularizers.cpp b/Wrappers/Python/src/cpu_regularizers.cpp
index e311570..3529ebd 100644
--- a/Wrappers/Python/src/cpu_regularizers.cpp
+++ b/Wrappers/Python/src/cpu_regularizers.cpp
@@ -27,7 +27,6 @@ limitations under the License.
 #include "boost/tuple/tuple.hpp"
 
 #include "SplitBregman_TV_core.h"
-#include "FGP_TV_core.h"
 #include "LLT_model_core.h"
 #include "PatchBased_Regul_core.h"
 #include "TGV_PD_core.h"
@@ -303,292 +302,6 @@ bp::list SplitBregman_TV(np::ndarray input, double d_mu, int iter, double d_epsi
 
 	}
 
-
-
-bp::list FGP_TV(np::ndarray input, double d_mu, int iter, double d_epsil, int methTV) {
-
-	// the result is in the following list
-	bp::list result;
-
-	int number_of_dims, dimX, dimY, dimZ, ll, j, count;
-	float *A, *D = NULL, *D_old = NULL, *P1 = NULL, *P2 = NULL, *P3 = NULL, *P1_old = NULL, *P2_old = NULL, *P3_old = NULL, *R1 = NULL, *R2 = NULL, *R3 = NULL;
-	float lambda, tk, tkp1, re, re1, re_old, epsil, funcval;
-
-	//number_of_dims = mxGetNumberOfDimensions(prhs[0]);
-	//dim_array = mxGetDimensions(prhs[0]);
-
-	number_of_dims = input.get_nd();
-	int dim_array[3];
-
-	dim_array[0] = input.shape(0);
-	dim_array[1] = input.shape(1);
-	if (number_of_dims == 2) {
-		dim_array[2] = -1;
-	}
-	else {
-		dim_array[2] = input.shape(2);
-	}
-	// Parameter handling is be done in Python
-	///*Handling Matlab input data*/
-	//if ((nrhs < 2) || (nrhs > 5)) mexErrMsgTxt("At least 2 parameters is required: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1')");
-
-	///*Handling Matlab input data*/
-	//A = (float *)mxGetData(prhs[0]); /*noisy image (2D/3D) */
-	A = reinterpret_cast<float *>(input.get_data());
-
-	//mu = (float)mxGetScalar(prhs[1]); /* regularization parameter */
-	lambda = (float)d_mu;
-
-	//iter = 35; /* default iterations number */
-
-	//epsil = 0.0001; /* default tolerance constant */
-	epsil = (float)d_epsil;
-	//methTV = 0;  /* default isotropic TV penalty */
-	//if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5))  iter = (int)mxGetScalar(prhs[2]); /* iterations number */
-	//if ((nrhs == 4) || (nrhs == 5))  epsil = (float)mxGetScalar(prhs[3]); /* tolerance constant */
-	//if (nrhs == 5) {
-	//	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 (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input image must be in a single precision"); }
-
-	//plhs[1] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL);
-	bp::tuple shape1 = bp::make_tuple(dim_array[0], dim_array[1]);
-	np::dtype dtype = np::dtype::get_builtin<float>();
-	np::ndarray out1 = np::zeros(shape1, dtype);
-	
-	//float *funcvalA = (float *)mxGetData(plhs[1]);
-	float * funcvalA = reinterpret_cast<float *>(out1.get_data());
-	//if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input image must be in a single precision"); }
-
-	/*Handling Matlab output data*/
-	dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
-
-	tk = 1.0f;
-	tkp1 = 1.0f;
-	count = 1;
-	re_old = 0.0f;
-
-	if (number_of_dims == 2) {
-		dimZ = 1; /*2D case*/
-		/*D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
-		D_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
-		P1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
-		P2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
-		P1_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
-		P2_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
-		R1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
-		R2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));*/
-
-		bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1]);
-		np::dtype dtype = np::dtype::get_builtin<float>();
-
-
-		np::ndarray npD      = np::zeros(shape, dtype);
-		np::ndarray npD_old  = np::zeros(shape, dtype);
-		np::ndarray npP1     = np::zeros(shape, dtype);
-		np::ndarray npP2     = np::zeros(shape, dtype);
-		np::ndarray npP1_old = np::zeros(shape, dtype);
-		np::ndarray npP2_old = np::zeros(shape, dtype);
-		np::ndarray npR1     = np::zeros(shape, dtype);
-		np::ndarray npR2     = np::zeros(shape, dtype);
-
-		D      = reinterpret_cast<float *>(npD.get_data());
-		D_old  = reinterpret_cast<float *>(npD_old.get_data());
-		P1     = reinterpret_cast<float *>(npP1.get_data());
-		P2     = reinterpret_cast<float *>(npP2.get_data());
-		P1_old = reinterpret_cast<float *>(npP1_old.get_data());
-		P2_old = reinterpret_cast<float *>(npP2_old.get_data());
-		R1     = reinterpret_cast<float *>(npR1.get_data());
-		R2     = reinterpret_cast<float *>(npR2.get_data());
-
-		/* begin iterations */
-		for (ll = 0; ll<iter; ll++) {
-			/* computing the gradient of the objective function */
-			Obj_func2D(A, D, R1, R2, lambda, dimX, dimY);
-
-			/*Taking a step towards minus of the gradient*/
-			Grad_func2D(P1, P2, D, R1, R2, lambda, dimX, dimY);
-  
-            /* projection step */
-            Proj_func2D(P1, P2, methTV, dimX, dimY);
-            
-            /*updating R and t*/
-            tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f;
-            Rupd_func2D(P1, P1_old, P2, P2_old, R1, R2, tkp1, tk, dimX, dimY);
-
-			/* calculate norm */
-			re = 0.0f; re1 = 0.0f;
-			for (j = 0; j<dimX*dimY*dimZ; j++)
-			{
-				re += pow(D[j] - D_old[j], 2);
-				re1 += pow(D[j], 2);
-			}
-			re = sqrt(re) / sqrt(re1);
-			if (re < epsil)  count++;
-			if (count > 4) {
-				Obj_func2D(A, D, P1, P2, lambda, dimX, dimY);
-				funcval = 0.0f;
-				for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
-				//funcvalA[0] = sqrt(funcval);
-				float fv = sqrt(funcval);
-				std::memcpy(funcvalA, &fv, sizeof(float));
-				break;
-			}
-
-			/* check that the residual norm is decreasing */
-			if (ll > 2) {
-				if (re > re_old) {
-					Obj_func2D(A, D, P1, P2, lambda, dimX, dimY);
-					funcval = 0.0f;
-					for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
-					//funcvalA[0] = sqrt(funcval);
-					float fv = sqrt(funcval);
-					std::memcpy(funcvalA, &fv, sizeof(float));
-					break;
-				}
-			}
-			re_old = re;
-			/*printf("%f %i %i \n", re, ll, count); */
-
-			/*storing old values*/
-			copyIm(D, D_old, dimX, dimY, dimZ);
-			copyIm(P1, P1_old, dimX, dimY, dimZ);
-			copyIm(P2, P2_old, dimX, dimY, dimZ);
-			tk = tkp1;
-
-			/* calculating the objective function value */
-			if (ll == (iter - 1)) {
-				Obj_func2D(A, D, P1, P2, lambda, dimX, dimY);
-				funcval = 0.0f;
-				for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
-				//funcvalA[0] = sqrt(funcval);
-				float fv = sqrt(funcval);
-				std::memcpy(funcvalA, &fv, sizeof(float));
-			}
-		}
-		//printf("FGP-TV iterations stopped at iteration %i with the function value %f \n", ll, funcvalA[0]);
-		result.append<np::ndarray>(npD);
-		result.append<np::ndarray>(out1);
-		result.append<int>(ll);
-	}
-	if (number_of_dims == 3) {
-		/*D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-		D_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-		P1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-		P2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-		P3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-		P1_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-		P2_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-		P3_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-		R1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-		R2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-		R3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));*/
-		bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1], dim_array[2]);
-		np::dtype dtype = np::dtype::get_builtin<float>();
-		
-		np::ndarray npD      = np::zeros(shape, dtype);
-		np::ndarray npD_old  = np::zeros(shape, dtype);
-		np::ndarray npP1     = np::zeros(shape, dtype);
-		np::ndarray npP2     = np::zeros(shape, dtype);
-		np::ndarray npP3     = np::zeros(shape, dtype);
-		np::ndarray npP1_old = np::zeros(shape, dtype);
-		np::ndarray npP2_old = np::zeros(shape, dtype);
-		np::ndarray npP3_old = np::zeros(shape, dtype);
-		np::ndarray npR1     = np::zeros(shape, dtype);
-		np::ndarray npR2     = np::zeros(shape, dtype);
-		np::ndarray npR3     = np::zeros(shape, dtype);
-
-		D      = reinterpret_cast<float *>(npD.get_data());
-		D_old  = reinterpret_cast<float *>(npD_old.get_data());
-		P1     = reinterpret_cast<float *>(npP1.get_data());
-		P2     = reinterpret_cast<float *>(npP2.get_data());
-		P3     = reinterpret_cast<float *>(npP3.get_data());
-		P1_old = reinterpret_cast<float *>(npP1_old.get_data());
-		P2_old = reinterpret_cast<float *>(npP2_old.get_data());
-		P3_old = reinterpret_cast<float *>(npP3_old.get_data());
-		R1     = reinterpret_cast<float *>(npR1.get_data());
-		R2     = reinterpret_cast<float *>(npR2.get_data());
-		R3     = reinterpret_cast<float *>(npR3.get_data());
-		/* begin iterations */
-		for (ll = 0; ll<iter; ll++) {
-			/* computing the gradient of the objective function */
-			Obj_func3D(A, D, R1, R2, R3, lambda, dimX, dimY, dimZ);
-			/*Taking a step towards minus of the gradient*/
-			Grad_func3D(P1, P2, P3, D, R1, R2, R3, lambda, dimX, dimY, dimZ);
-
-			/* projection step */
-			Proj_func3D(P1, P2, P3, dimX, dimY, dimZ);
-
-			/*updating R and t*/
-			tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f;
-			Rupd_func3D(P1, P1_old, P2, P2_old, P3, P3_old, R1, R2, R3, tkp1, tk, dimX, dimY, dimZ);
-
-			/* calculate norm - stopping rules*/
-			re = 0.0f; re1 = 0.0f;
-			for (j = 0; j<dimX*dimY*dimZ; j++)
-			{
-				re += pow(D[j] - D_old[j], 2);
-				re1 += pow(D[j], 2);
-			}
-			re = sqrt(re) / sqrt(re1);
-			/* stop if the norm residual is less than the tolerance EPS */
-			if (re < epsil)  count++;
-			if (count > 3) {
-				Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ);
-				funcval = 0.0f;
-				for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
-				//funcvalA[0] = sqrt(funcval);
-				float fv = sqrt(funcval);
-				std::memcpy(funcvalA, &fv, sizeof(float));
-				break;
-			}
-
-			/* check that the residual norm is decreasing */
-			if (ll > 2) {
-				if (re > re_old) {
-					Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ);
-					funcval = 0.0f;
-					for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
-					//funcvalA[0] = sqrt(funcval);
-					float fv = sqrt(funcval);
-					std::memcpy(funcvalA, &fv, sizeof(float));
-					break;
-				}
-			}
-
-			re_old = re;
-			/*printf("%f %i %i \n", re, ll, count); */
-
-			/*storing old values*/
-			copyIm(D, D_old, dimX, dimY, dimZ);
-			copyIm(P1, P1_old, dimX, dimY, dimZ);
-			copyIm(P2, P2_old, dimX, dimY, dimZ);
-			copyIm(P3, P3_old, dimX, dimY, dimZ);
-			tk = tkp1;
-
-			if (ll == (iter - 1)) {
-				Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ);
-				funcval = 0.0f;
-				for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
-				//funcvalA[0] = sqrt(funcval);
-				float fv = sqrt(funcval);
-				std::memcpy(funcvalA, &fv, sizeof(float));
-			}
-
-		}
-		//printf("FGP-TV iterations stopped at iteration %i with the function value %f \n", ll, funcvalA[0]);
-		result.append<np::ndarray>(npD);
-		result.append<np::ndarray>(out1);
-		result.append<int>(ll);
-	}
-
-	return result;
-}
-
 bp::list LLT_model(np::ndarray input, double d_lambda, double d_tau, int iter, double d_epsil, int switcher) {
 	// the result is in the following list
 	bp::list result;
@@ -1022,9 +735,6 @@ bp::list TGV_PD(np::ndarray input, double d_lambda, double d_alpha1, double d_al
 		result.append<np::ndarray>(npU);
 	}
 	
-
-	
-	
 	return result;
 }
 
@@ -1040,7 +750,6 @@ BOOST_PYTHON_MODULE(cpu_regularizers_boost)
 	np::dtype dt2 = np::dtype::get_builtin<uint16_t>();
 
 	def("SplitBregman_TV", SplitBregman_TV);
-	def("FGP_TV", FGP_TV);
 	def("LLT_model", LLT_model);
 	def("PatchBased_Regul", PatchBased_Regul);
 	def("TGV_PD", TGV_PD);
diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx
index e7ff78c..60e8627 100644
--- a/Wrappers/Python/src/cpu_regularizers.pyx
+++ b/Wrappers/Python/src/cpu_regularizers.pyx
@@ -11,62 +11,108 @@ 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.
 
-Author: Edoardo Pasca
+Author: Edoardo Pasca, Daniil Kazantsev
 """
 
 import cython
 import numpy as np
 cimport numpy as np
 
-cdef extern float TV_ROF(float *A, float *B, int dimX, int dimY, int dimZ, 
-                         int iterationsNumb, float tau, float flambda);
+cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ);
+cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
 
+def TV_ROF_CPU(inputData, regularization_parameter, iterationsNumb, marching_step_parameter):
+    if inputData.ndim == 2:
+        return TV_ROF_2D(inputData, regularization_parameter, iterationsNumb, marching_step_parameter)
+    elif inputData.ndim == 3:
+        return TV_ROF_3D(inputData, regularization_parameter, iterationsNumb, marching_step_parameter)
+
+def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularization_parameter,
+                     int iterationsNumb,                     
+                     float marching_step_parameter):
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')
+                   
+    # Run ROF iterations for 2D data 
+    TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], regularization_parameter, iterationsNumb, marching_step_parameter, dims[0], dims[1], 1)
+    
+    return outputData
+            
+def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     int iterationsNumb,
+                     float regularization_parameter,
+                     float marching_step_parameter):
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+    
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+           
+    # Run ROF iterations for 3D data 
+    TV_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularization_parameter, iterationsNumb, marching_step_parameter, dims[0], dims[1], dims[2])
 
-def ROF_TV(inputData, iterations, regularization_parameter, 
-           marching_step_parameter):
+    return outputData    
+                     
+def TV_FGP_CPU(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM):
     if inputData.ndim == 2:
-        return ROF_TV_2D(inputData, iterations, regularization_parameter, 
-                         marching_step_parameter)
+        return TV_FGP_2D(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
     elif inputData.ndim == 3:
-        return ROF_TV_3D(inputData, iterations, regularization_parameter, 
-                         marching_step_parameter)
+        return TV_FGP_3D(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
 
-def ROF_TV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                     int iterations,
+def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
                      float regularization_parameter,
-                     float marching_step_parameter
-                     ):
+                     int iterationsNumb, 
+                     float tolerance_param,
+                     int methodTV,
+                     int nonneg,
+                     int printM):
                          
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
             np.zeros([dims[0],dims[1]], dtype='float32')
                    
     #/* Run ROF iterations for 2D data */
-    TV_ROF(&inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations,
-               marching_step_parameter, regularization_parameter)
+    TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], regularization_parameter, 
+                       iterationsNumb, 
+                       tolerance_param,
+                       methodTV,
+                       nonneg,
+                       printM,
+                       dims[0], dims[1], 1)
     
-    return B
-        
+    return outputData        
             
-def ROF_TV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     int iterations,
+def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
                      float regularization_parameter,
-                     float marching_step_parameter
-                     ):
-    cdef long dims[2]
+                     int iterationsNumb, 
+                     float tolerance_param,
+                     int methodTV,
+                     int nonneg,
+                     int printM):
+    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"] B = \
-            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+            np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
            
     #/* Run ROF iterations for 3D data */
-    TV_ROF(&inputData[0,0, 0], &B[0,0, 0], dims[0], dims[1], dims[2], iterations,
-             marching_step_parameter, regularization_parameter)
-
-    return B    
-                     
+    TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularization_parameter,
+                       iterationsNumb, 
+                       tolerance_param,
+                       methodTV,
+                       nonneg,
+                       printM,
+                       dims[0], dims[1], dims[2])
+    return outputData 
diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx
index 5a5d274..f96156a 100644
--- a/Wrappers/Python/src/gpu_regularizers.pyx
+++ b/Wrappers/Python/src/gpu_regularizers.pyx
@@ -11,7 +11,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.
 
-Author: Edoardo Pasca
+Author: Edoardo Pasca, Daniil Kazantsev
 """
 
 import cython
@@ -25,12 +25,16 @@ cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec,
                                 int N, int M,  int Z, int dimension, 
                                 int SearchW, int SimilW, 
                                 int SearchW_real, float denh2, float lambdaf);
-cdef extern void TV_ROF_GPU(float* A, float* B, int N, int M, int Z, int iter, float tau, float lambdaf);
+cdef extern void TV_ROF_GPU_main(float* Input, float* Output, float lambdaPar, int iter, float tau, int N, int M, int Z);
+cdef extern void TV_FGP_GPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z);
+
+# padding function
 cdef extern float pad_crop(float *A, float *Ap, 
                            int OldSizeX, int OldSizeY, int OldSizeZ, 
                            int NewSizeX, int NewSizeY, int NewSizeZ, 
                            int padXY, int switchpad_crop);
-
+                           
+#Diffusion 4th order regularizer
 def Diff4thHajiaboli(inputData, 
                      edge_preserv_parameter, 
                      iterations, 
@@ -48,7 +52,7 @@ def Diff4thHajiaboli(inputData,
                      iterations, 
                      time_marching_parameter,
                      regularization_parameter)
-        
+# patch-based nonlocal regularization
 def NML(inputData, 
                      SearchW_real, 
                      SimilW, 
@@ -66,23 +70,48 @@ def NML(inputData,
                      SimilW, 
                      h,
                      lambdaf)
-
-def ROF_TV_GPU(inputData,
+                     
+# Total-variation Rudin-Osher-Fatemi (ROF)
+def TV_ROF_GPU(inputData,
+                     regularization_parameter,
                      iterations, 
-                     time_marching_parameter,
-                     regularization_parameter):
+                     time_marching_parameter):
     if inputData.ndim == 2:
         return ROFTV2D(inputData, 
-                     iterations, 
-                     time_marching_parameter,
-                     regularization_parameter)
+                     regularization_parameter,
+                     iterations,
+                     time_marching_parameter)
     elif inputData.ndim == 3:
         return ROFTV3D(inputData, 
+                     regularization_parameter,
                      iterations, 
-                     time_marching_parameter,
-                     regularization_parameter)
-
-                    
+                     time_marching_parameter)
+                     
+# Total-variation Fast-Gradient-Projection (FGP)
+def TV_FGP_GPU(inputData,
+                     regularization_parameter,
+                     iterations, 
+                     tolerance_param,
+                     methodTV,
+                     nonneg,
+                     printM):
+    if inputData.ndim == 2:
+        return FGPTV2D(inputData,
+                     regularization_parameter,
+                     iterations, 
+                     tolerance_param,
+                     methodTV,
+                     nonneg,
+                     printM)
+    elif inputData.ndim == 3:
+        return FGPTV3D(inputData,
+                     regularization_parameter,
+                     iterations, 
+                     tolerance_param,
+                     methodTV,
+                     nonneg,
+                     printM)
+#****************************************************************#
 def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
                      float edge_preserv_parameter, 
                      int iterations, 
@@ -329,48 +358,106 @@ def NML3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
              switchpad_crop)
     
     return B                   
-                     
+  
+# Total-variation ROF 
 def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularization_parameter,
                      int iterations, 
-                     float time_marching_parameter,
-                     float regularization_parameter):
+                     float time_marching_parameter):
     
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
 
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
 		    np.zeros([dims[0],dims[1]], dtype='float32')
           
     # Running CUDA code here    
-    TV_ROF_GPU(            
-            &inputData[0,0], &B[0,0], 
-                       dims[0], dims[1], 0, 
+    TV_ROF_GPU_main(            
+            &inputData[0,0], &outputData[0,0], 
+                       regularization_parameter,
                        iterations , 
                        time_marching_parameter, 
-                       regularization_parameter);   
+                       dims[0], dims[1], 1);   
      
-    return B
+    return outputData
     
 def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularization_parameter,
                      int iterations, 
-                     float time_marching_parameter,
-                     float regularization_parameter):
+                     float time_marching_parameter):
     
     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"] B = \
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
 		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
           
     # Running CUDA code here    
-    TV_ROF_GPU(            
-            &inputData[0,0,0], &B[0,0,0], 
-                       dims[0], dims[1], dims[2], 
+    TV_ROF_GPU_main(            
+            &inputData[0,0,0], &outputData[0,0,0], 
+                       regularization_parameter,
                        iterations , 
                        time_marching_parameter, 
-                       regularization_parameter);   
+                       dims[0], dims[1], dims[2]);   
      
-    return B
+    return outputData
+
+# Total-variation Fast-Gradient-Projection (FGP)
+def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularization_parameter,
+                     int iterations, 
+                     float tolerance_param,
+                     int methodTV,
+                     int nonneg,
+                     int printM):
+    
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+		    np.zeros([dims[0],dims[1]], dtype='float32')
+          
+    # Running CUDA code here    
+    TV_FGP_GPU_main(&inputData[0,0], &outputData[0,0],                        
+                       regularization_parameter, 
+                       iterations, 
+                       tolerance_param,
+                       methodTV,
+                       nonneg,
+                       printM,
+                       dims[0], dims[1], 1);   
+     
+    return outputData
+    
+def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularization_parameter,
+                     int iterations, 
+                     float tolerance_param,
+                     int methodTV,
+                     int nonneg,
+                     int printM):
+    
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+          
+    # Running CUDA code here    
+    TV_FGP_GPU_main(            
+            &inputData[0,0,0], &outputData[0,0,0], 
+                       regularization_parameter , 
+                       iterations, 
+                       tolerance_param,
+                       methodTV,
+                       nonneg,
+                       printM,
+                       dims[0], dims[1], dims[2]);   
+     
+    return outputData    
diff --git a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py
new file mode 100644
index 0000000..63be1a0
--- /dev/null
+++ b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py
@@ -0,0 +1,219 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Feb 22 11:39:43 2018
+
+Testing CPU implementation against the GPU one
+
+@author: Daniil Kazantsev
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+import timeit
+from ccpi.filters.regularizers import ROF_TV, FGP_TV
+
+###############################################################################
+def printParametersToString(pars):
+        txt = r''
+        for key, value in pars.items():
+            if key== 'algorithm' :
+                txt += "{0} = {1}".format(key, value.__name__)
+            elif key == 'input':
+                txt += "{0} = {1}".format(key, np.shape(value))
+            else:
+                txt += "{0} = {1}".format(key, value)
+            txt += '\n'
+        return txt
+###############################################################################
+def rmse(im1, im2):
+    a, b = im1.shape
+    rmse = np.sqrt(np.sum((im1 - im2) ** 2) / float(a * b))    
+    return rmse
+        
+filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
+
+# read image
+Im = plt.imread(filename)                     
+Im = np.asarray(Im, dtype='float32')
+
+Im = Im/255
+perc = 0.075
+u0 = Im + np.random.normal(loc = Im ,
+                                  scale = perc * Im , 
+                                  size = np.shape(Im))
+# map the u0 u0->u0>0
+f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+u0 = f(u0).astype('float32')
+
+
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print ("____________ROF-TV bench___________________")
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+
+## plot 
+fig = plt.figure(1)
+plt.suptitle('Comparison of ROF-TV regularizer using CPU and GPU implementations')
+a=fig.add_subplot(1,4,1)
+a.set_title('Noisy Image')
+imgplot = plt.imshow(u0,cmap="gray")
+
+# set parameters
+pars = {'algorithm': ROF_TV, \
+        'input' : u0,\
+        'regularization_parameter':0.04,\
+        'number_of_iterations': 1200,\
+        'time_marching_parameter': 0.0025        
+        }
+print ("#############ROF TV CPU####################")
+start_time = timeit.default_timer()
+rof_cpu = ROF_TV(pars['input'],
+             pars['regularization_parameter'],
+             pars['number_of_iterations'],
+             pars['time_marching_parameter'],'cpu')
+rms = rmse(Im, rof_cpu)
+pars['rmse'] = rms
+
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,4,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(rof_cpu, cmap="gray")
+plt.title('{}'.format('CPU results'))
+
+
+print ("##############ROF TV GPU##################")
+start_time = timeit.default_timer()
+rof_gpu = ROF_TV(pars['input'], 
+                     pars['regularization_parameter'],
+                     pars['number_of_iterations'], 
+                     pars['time_marching_parameter'],'gpu')
+                     
+rms = rmse(Im, rof_gpu)
+pars['rmse'] = rms
+pars['algorithm'] = ROF_TV
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,4,3)
+
+# 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(rof_gpu, cmap="gray")
+plt.title('{}'.format('GPU results'))
+
+
+print ("--------Compare the results--------")
+tolerance = 1e-05
+diff_im = np.zeros(np.shape(rof_cpu))
+diff_im = abs(rof_cpu - rof_gpu)
+diff_im[diff_im > tolerance] = 1
+a=fig.add_subplot(1,4,4)
+imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray")
+plt.title('{}'.format('Pixels larger threshold difference'))
+if (diff_im.sum() > 1):
+    print ("Arrays do not match!")
+else:
+    print ("Arrays match")
+
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print ("____________FGP-TV bench___________________")
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+
+## plot 
+fig = plt.figure(2)
+plt.suptitle('Comparison of FGP-TV regularizer using CPU and GPU implementations')
+a=fig.add_subplot(1,4,1)
+a.set_title('Noisy Image')
+imgplot = plt.imshow(u0,cmap="gray")
+
+# set parameters
+pars = {'algorithm' : FGP_TV, \
+        'input' : u0,\
+        'regularization_parameter':0.04, \
+        'number_of_iterations' :1200 ,\
+        'tolerance_constant':0.00001,\
+        'methodTV': 0 ,\
+        'nonneg': 0 ,\
+        'printingOut': 0 
+        }
+        
+print ("#############FGP TV CPU####################")
+start_time = timeit.default_timer()
+fgp_cpu = FGP_TV(pars['input'], 
+              pars['regularization_parameter'],
+              pars['number_of_iterations'],
+              pars['tolerance_constant'], 
+              pars['methodTV'],
+              pars['nonneg'],
+              pars['printingOut'],'cpu')  
+             
+             
+rms = rmse(Im, fgp_cpu)
+pars['rmse'] = rms
+
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,4,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(fgp_cpu, cmap="gray")
+plt.title('{}'.format('CPU results'))
+
+
+print ("##############FGP TV GPU##################")
+start_time = timeit.default_timer()
+fgp_gpu = FGP_TV(pars['input'], 
+              pars['regularization_parameter'],
+              pars['number_of_iterations'],
+              pars['tolerance_constant'], 
+              pars['methodTV'],
+              pars['nonneg'],
+              pars['printingOut'],'gpu')
+                                   
+rms = rmse(Im, fgp_gpu)
+pars['rmse'] = rms
+pars['algorithm'] = FGP_TV
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,4,3)
+
+# 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(fgp_gpu, cmap="gray")
+plt.title('{}'.format('GPU results'))
+
+
+print ("--------Compare the results--------")
+tolerance = 1e-05
+diff_im = np.zeros(np.shape(rof_cpu))
+diff_im = abs(fgp_cpu - fgp_gpu)
+diff_im[diff_im > tolerance] = 1
+a=fig.add_subplot(1,4,4)
+imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray")
+plt.title('{}'.format('Pixels larger threshold difference'))
+if (diff_im.sum() > 1):
+    print ("Arrays do not match!")
+else:
+    print ("Arrays match")
+
+
diff --git a/Wrappers/Python/test/test_gpu_regularizers.py b/Wrappers/Python/test/test_gpu_regularizers.py
index 735a25d..640b3f9 100644
--- a/Wrappers/Python/test/test_gpu_regularizers.py
+++ b/Wrappers/Python/test/test_gpu_regularizers.py
@@ -6,14 +6,13 @@ Created on Tue Jan 30 10:24:26 2018
 @author: ofn77899
 """
 
-
-
 import matplotlib.pyplot as plt
 import numpy as np
 import os    
 from enum import Enum
 import timeit
 from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML
+from ccpi.filters.regularizers import ROF_TV, FGP_TV
 ###############################################################################
 def printParametersToString(pars):
         txt = r''
@@ -52,14 +51,15 @@ u0 = f(u0).astype('float32')
 ## plot 
 fig = plt.figure()
 
-a=fig.add_subplot(2,3,1)
+a=fig.add_subplot(2,4,1)
 a.set_title('noise')
 imgplot = plt.imshow(u0,cmap="gray")
 
 
 ## Diff4thHajiaboli
 start_time = timeit.default_timer()
-pars = {'algorithm' : Diff4thHajiaboli , \
+pars = {
+'algorithm' : Diff4thHajiaboli , \
         'input' : u0,
         'edge_preserv_parameter':0.1 , \
 'number_of_iterations' :250 ,\
@@ -78,7 +78,7 @@ pars['rmse'] = rms
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
 print (txtstr)
-a=fig.add_subplot(2,3,2)
+a=fig.add_subplot(2,4,2)
 
 # these are matplotlib.patch.Patch properties
 props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
@@ -87,14 +87,15 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=12,
          verticalalignment='top', bbox=props)
 imgplot = plt.imshow(d4h, cmap="gray")
 
-a=fig.add_subplot(2,3,5)
+a=fig.add_subplot(2,4,6)
 
 # these are matplotlib.patch.Patch properties
 props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
 # place a text box in upper left in axes coords
 a.text(0.05, 0.95, 'd4h - u0', transform=a.transAxes, fontsize=12,
          verticalalignment='top', bbox=props)
-imgplot = plt.imshow((d4h - u0)**2, cmap="gray")
+imgplot = plt.imshow((d4h - u0)**2,  vmin=0, vmax=0.03, cmap="gray")
+plt.colorbar(ticks=[0, 0.03], orientation='vertical')
 
 
 ## Patch Based Regul NML
@@ -109,6 +110,7 @@ pars = {'algorithm' : NML , \
 }
 """
 pars = {
+'algorithm' : NML , \
         'input' : u0,
         'regularization_parameter': 0.01,\
         'searching_window_ratio':3, \
@@ -126,7 +128,7 @@ pars['rmse'] = rms
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
 print (txtstr)
-a=fig.add_subplot(2,3,3)
+a=fig.add_subplot(2,4,3)
 
 # these are matplotlib.patch.Patch properties
 props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
@@ -135,14 +137,103 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=12,
          verticalalignment='top', bbox=props)
 imgplot = plt.imshow(nml, cmap="gray")
 
-a=fig.add_subplot(2,3,6)
+a=fig.add_subplot(2,4,7)
 
 # these are matplotlib.patch.Patch properties
 props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
 # place a text box in upper left in axes coords
 a.text(0.05, 0.95, 'nml - u0', transform=a.transAxes, fontsize=14,
          verticalalignment='top', bbox=props)
-imgplot = plt.imshow((nml - u0)**2, cmap="gray")
+imgplot = plt.imshow((nml - u0)**2,  vmin=0, vmax=0.03, cmap="gray")
+plt.colorbar(ticks=[0, 0.03], orientation='vertical')
 
-plt.show()
         
+        
+## Rudin-Osher-Fatemi (ROF) TV regularization
+start_time = timeit.default_timer()
+pars = {
+'algorithm' : ROF_TV , \
+        'input' : u0,
+        'regularization_parameter': 0.04,\
+        'number_of_iterations':300,\
+        'time_marching_parameter': 0.0025
+        
+	}
+	
+rof_tv = TV_ROF_GPU(pars['input'], 
+                     pars['regularization_parameter'],
+                     pars['number_of_iterations'], 
+                     pars['time_marching_parameter'],'gpu')
+                     
+rms = rmse(Im, rof_tv)
+pars['rmse'] = rms
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(2,4,4)
+
+# 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=12,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow(rof_tv, cmap="gray")
+
+a=fig.add_subplot(2,4,8)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, 'rof_tv - u0', transform=a.transAxes, fontsize=14,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow((rof_tv - u0)**2, vmin=0, vmax=0.03, cmap="gray")
+plt.colorbar(ticks=[0, 0.03], orientation='vertical')
+plt.show()
+
+## Fast-Gradient Projection TV regularization
+"""
+start_time = timeit.default_timer()
+
+pars = {'algorithm' : FGP_TV, \
+        'input' : u0,\
+        'regularization_parameter':0.04, \
+        'number_of_iterations' :1200 ,\
+        'tolerance_constant':0.00001,\
+        'methodTV': 0 ,\
+        'nonneg': 0 ,\
+        'printingOut': 0 
+        }
+
+fgp_gpu = FGP_TV(pars['input'], 
+              pars['regularization_parameter'],
+              pars['number_of_iterations'],
+              pars['tolerance_constant'], 
+              pars['methodTV'],
+              pars['nonneg'],
+              pars['printingOut'],'gpu')
+                                   
+rms = rmse(Im, fgp_gpu)
+pars['rmse'] = rms
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(2,4,4)
+
+# 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=12,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow(fgp_gpu, cmap="gray")
+
+a=fig.add_subplot(2,4,8)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, 'fgp_gpu - u0', transform=a.transAxes, fontsize=14,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow((fgp_gpu - u0)**2, vmin=0, vmax=0.03, cmap="gray")
+plt.colorbar(ticks=[0, 0.03], orientation='vertical')
+plt.show()
+"""