tol work

13 files changed, 217 insertions, 108 deletions

diff --git a/Readme.md b/Readme.md
index 92b4273..cc96b2c 100644
--- a/Readme.md
+++ b/Readme.md
@@ -109,7 +109,7 @@ conda install ccpi-regulariser -c ccpi -c conda-forge
 
 #### Python (conda-build)
 ```
-	export CIL_VERSION=19.02
+	export CIL_VERSION=19.03
 	conda build recipe/ --numpy 1.12 --python 3.5  
 	conda install ccpi-regulariser=${CIL_VERSION} --use-local --force
 	cd demos/
diff --git a/build/run.sh b/build/run.sh
index 332d660..d450299 100644
--- a/build/run.sh
+++ b/build/run.sh
@@ -1,22 +1,24 @@
 #!/bin/bash  
 echo "Building CCPi-regularisation Toolkit using CMake"  
-rm -r build
+rm -r build_proj
 # Requires Cython, install it first: 
 # pip install cython
-mkdir build
-cd build/
+mkdir build_proj
+cd build_proj/
 make clean
+export CIL_VERSION=19.03
 # install Python modules without CUDA
-#cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=OFF -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
+cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=OFF -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
 # install Python modules with CUDA
 # cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
 # install Matlab modules with CUDA
-cmake ../ -DBUILD_PYTHON_WRAPPER=OFF -DMatlab_ROOT_DIR=/dls_sw/apps/matlab/r2014a/ -DBUILD_MATLAB_WRAPPER=ON -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
+# cmake ../ -DBUILD_PYTHON_WRAPPER=OFF -DMatlab_ROOT_DIR=/dls_sw/apps/matlab/r2014a/ -DBUILD_MATLAB_WRAPPER=ON -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
 make install
-# cp install/lib/libcilreg.so install/python/ccpi/filters
-#cd install/python
-#export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:../lib
-# spyder
-# one can also run Matlab in Linux as:
+#### Python
+#cp install/lib/libcilreg.so install/python/ccpi/filters
+cd install/python
+export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:../lib
+spyder
+##### one can also run Matlab in Linux as:
 #PATH="/path/to/mex/:$PATH" LD_LIBRARY_PATH="/path/to/library:$LD_LIBRARY_PATH" matlab
-PATH="/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/build/install/matlab/:$PATH" LD_LIBRARY_PATH="/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/build/install/lib:$LD_LIBRARY_PATH" matlab
+#PATH="/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/build/install/matlab/:$PATH" LD_LIBRARY_PATH="/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/build/install/lib:$LD_LIBRARY_PATH" matlab
diff --git a/demos/SoftwareX_supp/Demo_VolumeDenoise.py b/demos/SoftwareX_supp/Demo_VolumeDenoise.py
index 2387e94..17cdf4d 100644
--- a/demos/SoftwareX_supp/Demo_VolumeDenoise.py
+++ b/demos/SoftwareX_supp/Demo_VolumeDenoise.py
@@ -60,13 +60,14 @@ plt.title('3D Phantom, sagittal view')
 plt.show()
 #%%
 print ("____________________Applying regularisers_______________________")
+print ("________________________________________________________________")
 
 print ("#############ROF TV CPU####################")
 # set parameters
 pars = {'algorithm': ROF_TV, \
         'input' : phantom_noise,\
         'regularisation_parameter':0.04,\
-        'number_of_iterations': 100,\
+        'number_of_iterations': 600,\
         'time_marching_parameter': 0.0025
         }
 
@@ -116,6 +117,65 @@ win2d = win * (win.T)
 ssim_rof = Qtools.ssim(win2d)
 
 print("ROF-TV (gpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE_rof,ssim_rof[0],Run_time_rof))
+#%%
+print ("#############FGP TV CPU####################")
+# set parameters
+pars = {'algorithm':FGP_TV, \
+        'input' : phantom_noise,\
+        'regularisation_parameter':0.04,\
+        'number_of_iterations': 300,\
+        'time_marching_parameter': 0.0025,\
+        'tolerance_constant':1e-05,\
+        }
+
+tic=timeit.default_timer()
+fgp_cpu3D = FGP_TV(pars['input'],
+             pars['regularisation_parameter'],
+             pars['number_of_iterations'],
+             pars['time_marching_parameter'],'cpu')
+toc=timeit.default_timer()
+
+Run_time_fgp = toc - tic
+Qtools = QualityTools(phantom_tm, fgp_cpu3D)
+RMSE_rof = Qtools.rmse()
+
+# SSIM measure
+Qtools = QualityTools(phantom_tm[128,:,:]*255, fgp_cpu3D[128,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim_fgp = Qtools.ssim(win2d)
+
+print("FGP-TV (cpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE_rof,ssim_fgp[0],Run_time_fgp))
+#%%
+print ("#############FGP TV GPU####################")
+# set parameters
+pars = {'algorithm':FGP_TV, \
+        'input' : phantom_noise,\
+        'regularisation_parameter':0.04,\
+        'number_of_iterations': 300,\
+        'time_marching_parameter': 0.0025,\
+        'tolerance_constant':1e-05,\
+        }
 
+tic=timeit.default_timer()
+fgp_gpu3D = FGP_TV(pars['input'],
+             pars['regularisation_parameter'],
+             pars['number_of_iterations'],
+             pars['time_marching_parameter'],'gpu')
+toc=timeit.default_timer()
+
+Run_time_fgp = toc - tic
+Qtools = QualityTools(phantom_tm, fgp_gpu3D)
+RMSE_rof = Qtools.rmse()
+
+# SSIM measure
+Qtools = QualityTools(phantom_tm[128,:,:]*255, fgp_gpu3D[128,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim_fgp = Qtools.ssim(win2d)
+
+print("FGP-TV (gpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE_rof,ssim_fgp[0],Run_time_fgp))
 #%%
 
+
+
diff --git a/demos/demoMatlab_denoise.m b/demos/demoMatlab_denoise.m
index 5e92ee1..fa81f6d 100644
--- a/demos/demoMatlab_denoise.m
+++ b/demos/demoMatlab_denoise.m
@@ -2,9 +2,11 @@
 clear; close all
 fsep = '/';
 
-Path1 = sprintf(['..' fsep 'src' fsep 'Matlab' fsep 'mex_compile' fsep 'installed'], 1i);
-Path2 = sprintf([ data' fsep], 1i);
-Path3 = sprintf(['..' filesep 'src' filesep 'Matlab' filesep 'supp'], 1i);
+%Path1 = sprintf(['..' fsep 'src' fsep 'Matlab' fsep 'mex_compile' fsep 'installed'], 1i);
+Path1 = ('/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/src/Matlab/mex_compile/installed');
+Path2 = sprintf(['data' fsep], 1i);
+%Path3 = sprintf(['..' filesep 'src' filesep 'Matlab' filesep 'supp'], 1i);
+Path3 = '/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/src/Matlab/supp';
 addpath(Path1);
 addpath(Path2);
 addpath(Path3);
@@ -34,8 +36,8 @@ figure; imshow(u_rof, [0 1]); title('ROF-TV denoised image (CPU)');
 %%
 fprintf('Denoise using the FGP-TV model (CPU) \n');
 lambda_reg = 0.033;
-iter_fgp = 300; % number of FGP iterations
-epsil_tol =  1.0e-09; % tolerance
+iter_fgp = 200; % number of FGP iterations
+epsil_tol =  1.0e-05; % tolerance
 tic; u_fgp = FGP_TV(single(u0), lambda_reg, iter_fgp, epsil_tol); toc; 
 energyfunc_val_fgp = TV_energy(single(u_fgp),single(u0),lambda_reg, 1); % get energy function value
 rmseFGP = (RMSE(u_fgp(:),Im(:)));
diff --git a/demos/demo_cpu_regularisers.py b/demos/demo_cpu_regularisers.py
index d34607a..e58ed2f 100644
--- a/demos/demo_cpu_regularisers.py
+++ b/demos/demo_cpu_regularisers.py
@@ -30,8 +30,9 @@ def printParametersToString(pars):
             txt += '\n'
         return txt
 ###############################################################################
-#%%
-filename = os.path.join( "data" ,"lena_gray_512.tif")
+
+# filename = os.path.join( "data" ,"lena_gray_512.tif")
+filename = "/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/test/lena_gray_512.tif"
 
 # read image
 Im = plt.imread(filename)
@@ -127,23 +128,20 @@ imgplot = plt.imshow(u0,cmap="gray")
 pars = {'algorithm' : FGP_TV, \
         'input' : u0,\
         'regularisation_parameter':0.04, \
-        'number_of_iterations' :2000 ,\
+        'number_of_iterations' :200 ,\
         'tolerance_constant':1e-06,\
         'methodTV': 0 ,\
-        'nonneg': 0 ,\
-        'printingOut': 0 
+        'nonneg': 0
         }
         
 print ("#############FGP TV CPU####################")
 start_time = timeit.default_timer()
-fgp_cpu = FGP_TV(pars['input'], 
+fgp_cpu,info_vec = FGP_TV(pars['input'], 
               pars['regularisation_parameter'],
               pars['number_of_iterations'],
               pars['tolerance_constant'], 
               pars['methodTV'],
-              pars['nonneg'],
-              pars['printingOut'],'cpu')  
-             
+              pars['nonneg'],'cpu')
 
 Qtools = QualityTools(Im, fgp_cpu)
 pars['rmse'] = Qtools.rmse()
diff --git a/run.sh b/run.sh
new file mode 100755
index 0000000..c5c396a
--- /dev/null
+++ b/run.sh
@@ -0,0 +1,26 @@
+#!/bin/bash  
+echo "Building CCPi-regularisation Toolkit using CMake"  
+rm -r build_proj
+# Requires Cython, install it first: 
+# pip install cython
+mkdir build_proj
+cd build_proj/
+#make clean
+export CIL_VERSION=19.03
+# install Python modules without CUDA
+cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=OFF -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
+# install Python modules with CUDA
+# cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
+# install Matlab modules without CUDA
+# cmake ../ -DBUILD_PYTHON_WRAPPER=OFF -DMatlab_ROOT_DIR=/dls_sw/apps/matlab/r2014a/ -DBUILD_MATLAB_WRAPPER=ON -DBUILD_CUDA=OFF -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
+# install Matlab modules with CUDA
+# cmake ../ -DBUILD_PYTHON_WRAPPER=OFF -DMatlab_ROOT_DIR=/dls_sw/apps/matlab/r2014a/ -DBUILD_MATLAB_WRAPPER=ON -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
+make install
+#### Python
+#cp install/lib/libcilreg.so install/python/ccpi/filters
+cd install/python
+export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:../lib
+spyder
+##### Matlab (Linux)
+#PATH="/path/to/mex/:$PATH" LD_LIBRARY_PATH="/path/to/library:$LD_LIBRARY_PATH" matlab
+#PATH="/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/build_proj/install/matlab/:$PATH" LD_LIBRARY_PATH="/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/build_proj/install/lib:$LD_LIBRARY_PATH" matlab
diff --git a/src/Core/regularisers_CPU/FGP_TV_core.c b/src/Core/regularisers_CPU/FGP_TV_core.c
index 68d58b7..f613f0e 100644
--- a/src/Core/regularisers_CPU/FGP_TV_core.c
+++ b/src/Core/regularisers_CPU/FGP_TV_core.c
@@ -28,28 +28,33 @@ limitations under the License.
  * 4. eplsilon: tolerance constant 
  * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
  * 6. nonneg: 'nonnegativity (0 is OFF by default) 
- * 7. print information: 0 (off) or 1 (on) 
  *
  * Output:
- * [1] Filtered/regularized image
+ * [1] Filtered/regularized image/volume
+ * [2] Information vector which contains [iteration no., reached tolerance]
  *
  * 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"
  */
  
-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)
+float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int dimX, int dimY, int dimZ)
 {
-	int ll;
+    int ll;
     long j, DimTotal;
-	float re, re1;
-	float tk = 1.0f;
-    float tkp1=1.0f;
+    float re, re1;
+    re = 0.0f; re1 = 0.0f;
+    float tk = 1.0f;
+    float tkp1 =1.0f;
     int count = 0;
+    
+            /*adding info into info_vector */
+//    infovector[0] = 50.1f;  /*iterations number (if stopped earlier based on tolerance)*/
+//    infovector[1] = 0.55f;  /* reached tolerance */
 	
 	if (dimZ <= 1) {
-		/*2D case */
-		float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL;
-		DimTotal = (long)(dimX*dimY);
+	/*2D case */
+  	float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL;
+	DimTotal = (long)(dimX*dimY);
 		
         Output_prev = calloc(DimTotal, sizeof(float));
         P1 = calloc(DimTotal, sizeof(float));
@@ -59,7 +64,7 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio
         R1 = calloc(DimTotal, sizeof(float));
         R2 = calloc(DimTotal, sizeof(float)); 
 		
-		/* begin iterations */
+	/* begin iterations */
         for(ll=0; ll<iterationsNumb; ll++) {
             
             /* computing the gradient of the objective function */
@@ -79,24 +84,32 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio
             Rupd_func2D(P1, P1_prev, P2, P2_prev, R1, R2, tkp1, tk, DimTotal);
             
             /* check early stopping criteria */
-            re = 0.0f; re1 = 0.0f;
-            for(j=0; j<DimTotal; j++)
-            {
-                re += pow(Output[j] - Output_prev[j],2);
-                re1 += pow(Output[j],2);
+            if (epsil != 0.0f) {
+	            for(j=0; j<DimTotal; j++)
+        	    {
+        	        re += pow(Output[j] - Output_prev[j],2);
+        	        re1 += pow(Output[j],2);
+        	    }
+              re = sqrt(re)/sqrt(re1);
+              if (re < epsil)  count++;
+              if (count > 4) break;
+              
+            copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l);
             }
-            re = sqrt(re)/sqrt(re1);
-            if (re < epsil)  count++;
-				if (count > 4) break;
-            
             /*storing old values*/
-            copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l);
             copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), 1l);
             copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), 1l);
             tk = tkp1;
         }
-        if (printM == 1) printf("FGP-TV iterations stopped at iteration %i \n", ll);   
-		free(Output_prev); free(P1); free(P2); free(P1_prev); free(P2_prev); free(R1); free(R2);		
+        /*adding info into info_vector */
+        //info_vector[0] = (float)(ll);  /*iterations number (if stopped earlier based on tolerance)*/
+//        info_vector[1] = 0.55f;  /* reached tolerance */
+	
+	copyIm(Input, infovector, (long)(dimX), (long)(dimY), 1l);
+	
+//	printf("%f\n", infovector[128]);
+	
+	free(Output_prev); free(P1); free(P2); free(P1_prev); free(P2_prev); free(R1); free(R2);		
 	}
 	else {
 		/*3D case*/
@@ -134,26 +147,31 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio
             Rupd_func3D(P1, P1_prev, P2, P2_prev, P3, P3_prev, R1, R2, R3, tkp1, tk, DimTotal);
             
             /* calculate norm - stopping rules*/
-            re = 0.0f; re1 = 0.0f;
-            for(j=0; j<DimTotal; j++)
-            {
-                re += pow(Output[j] - Output_prev[j],2);
-                re1 += pow(Output[j],2);
-            }
+            if (epsil != 0.0f) {
+	            for(j=0; j<DimTotal; j++)
+        	    {
+        	        re += pow(Output[j] - Output_prev[j],2);
+        	        re1 += pow(Output[j],2);
+        	    }
             re = sqrt(re)/sqrt(re1);
             /* stop if the norm residual is less than the tolerance EPS */
             if (re < epsil)  count++;
-            if (count > 4) break;            
-                        
-            /*storing old values*/
-            copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
+            if (count > 4) break; 
+
+            copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));           
+            }
+            
+            /*storing old values*/           
             copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
             copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
             copyIm(P3, P3_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
             tk = tkp1;            
         }	
-		if (printM == 1) printf("FGP-TV iterations stopped at iteration %i \n", ll);   
-		free(Output_prev); free(P1); free(P2); free(P3); free(P1_prev); free(P2_prev); free(P3_prev); free(R1); free(R2); free(R3);
+        /*adding info into info_vector */
+        //infovector[0] = (float)(ll);  /*iterations number (if stopped earlier based on tolerance)*/
+	//infovector[1] = re;  /* reached tolerance */
+        
+	free(Output_prev); free(P1); free(P2); free(P3); free(P1_prev); free(P2_prev); free(P3_prev); free(R1); free(R2); free(R3);
 	}
 	return *Output;
 }
diff --git a/src/Core/regularisers_CPU/FGP_TV_core.h b/src/Core/regularisers_CPU/FGP_TV_core.h
index 3418604..04e6e80 100644
--- a/src/Core/regularisers_CPU/FGP_TV_core.h
+++ b/src/Core/regularisers_CPU/FGP_TV_core.h
@@ -35,10 +35,11 @@ limitations under the License.
  * 4. eplsilon: tolerance constant 
  * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
  * 6. nonneg: 'nonnegativity (0 is OFF by default) 
- * 7. print information: 0 (off) or 1 (on) 
  *
  * Output:
- * [1] Filtered/regularized image
+ * [1] Filtered/regularized image/volume
+ * [2] Information vector which contains [iteration no., reached tolerance]
+
  *
  * 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"
@@ -47,7 +48,7 @@ limitations under the License.
 #ifdef __cplusplus
 extern "C" {
 #endif
-CCPI_EXPORT 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);
+CCPI_EXPORT float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int dimX, int dimY, int dimZ);
 
 CCPI_EXPORT float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, long dimX, long dimY);
 CCPI_EXPORT float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float lambda, long dimX, long dimY);
diff --git a/src/Matlab/mex_compile/compileCPU_mex_Linux.m b/src/Matlab/mex_compile/compileCPU_mex_Linux.m
index 72a828e..f3d9ce1 100644
--- a/src/Matlab/mex_compile/compileCPU_mex_Linux.m
+++ b/src/Matlab/mex_compile/compileCPU_mex_Linux.m
@@ -2,9 +2,9 @@
 
 fsep = '/';
 
-pathcopyFrom = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'regularisers_CPU'], 1i);
-pathcopyFrom1 = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'CCPiDefines.h'], 1i);
-pathcopyFrom2 = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'inpainters_CPU'], 1i);
+pathcopyFrom = sprintf(['..' fsep '..' fsep 'Core' fsep 'regularisers_CPU'], 1i);
+pathcopyFrom1 = sprintf(['..' fsep '..' fsep 'Core' fsep 'CCPiDefines.h'], 1i);
+pathcopyFrom2 = sprintf(['..' fsep '..' fsep 'Core' fsep 'inpainters_CPU'], 1i);
 
 copyfile(pathcopyFrom, 'regularisers_CPU');
 copyfile(pathcopyFrom1, 'regularisers_CPU');
@@ -76,6 +76,7 @@ delete PatchSelect_core* Nonlocal_TV_core*
 delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core*
 fprintf('%s \n', '<<<<<<< Regularisers successfully compiled! >>>>>>>');
 
-pathA2 = sprintf(['..' fsep '..' fsep], 1i);
-cd(pathA2);
-cd demos
+%pathA2 = sprintf(['..' fsep '..' fsep], 1i);
+% cd(pathA2);
+cd('/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/demos')
+%cd demos
diff --git a/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c b/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c
index 642362f..603e0f4 100644
--- a/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c
+++ b/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c
@@ -47,13 +47,13 @@ void mexFunction(
     int number_of_dims, iter, methTV, printswitch, nonneg;
     mwSize dimX, dimY, dimZ;
     const mwSize *dim_array;
-    float *Input, *Output=NULL, lambda, epsil;
+    float *Input, *Output=NULL, *Output2=NULL, lambda, epsil;
     
     number_of_dims = mxGetNumberOfDimensions(prhs[0]);
     dim_array = mxGetDimensions(prhs[0]);
     
     /*Handling Matlab input data*/
-    if ((nrhs < 2) || (nrhs > 7)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter, Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch, print switch");
+    if ((nrhs < 2) || (nrhs > 6)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter, Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch");
     
     Input  = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
     lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
@@ -61,27 +61,24 @@ void mexFunction(
     epsil = 0.0001; /* default tolerance constant */
     methTV = 0;  /* default isotropic TV penalty */
     nonneg = 0; /* default nonnegativity switch, off - 0 */
-    printswitch = 0; /*default print is switched, off - 0 */
+//    printswitch = 0; /*default print is switched, off - 0 */
     
     if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
     
-    if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7))  iter = (int) mxGetScalar(prhs[2]); /* iterations number */
-    if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7))  epsil =  (float) mxGetScalar(prhs[3]); /* tolerance constant */
-    if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7))  {
+    if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6))  iter = (int) mxGetScalar(prhs[2]); /* iterations number */
+    if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6))  epsil =  (float) mxGetScalar(prhs[3]); /* tolerance constant */
+    if ((nrhs == 5) || (nrhs == 6))  {
         char *penalty_type;
         penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */
         if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',");
         if (strcmp(penalty_type, "l1") == 0)  methTV = 1;  /* enable 'l1' penalty */
         mxFree(penalty_type);
     }
-    if ((nrhs == 6) || (nrhs == 7))  {
+    if (nrhs == 6)  {
         nonneg = (int) mxGetScalar(prhs[5]);
         if ((nonneg != 0) && (nonneg != 1)) mexErrMsgTxt("Nonnegativity constraint can be enabled by choosing 1 or off - 0");
     }
-    if (nrhs == 7)  {
-        printswitch = (int) mxGetScalar(prhs[6]);
-        if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0");
-    }
+    
     
     /*Handling Matlab output data*/
     dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
@@ -89,9 +86,12 @@ void mexFunction(
     if (number_of_dims == 2) {
         dimZ = 1; /*2D case*/
         Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+        Output2 = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+    }
+    if (number_of_dims == 3) {    
+    Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
     }
-    if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
     
     /* running the function */
-    TV_FGP_CPU_main(Input, Output, lambda, iter, epsil, methTV, nonneg, printswitch, dimX, dimY, dimZ);
-}
\ No newline at end of file
+    TV_FGP_CPU_main(Input, Output, Output2, lambda, iter, epsil, methTV, nonneg, dimX, dimY, dimZ);
+}
diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py
index 588ea32..fb2c999 100644
--- a/src/Python/ccpi/filters/regularisers.py
+++ b/src/Python/ccpi/filters/regularisers.py
@@ -29,15 +29,14 @@ def ROF_TV(inputData, regularisation_parameter, iterations,
                          .format(device))
 
 def FGP_TV(inputData, regularisation_parameter,iterations,
-                     tolerance_param, methodTV, nonneg, printM, device='cpu'):
+                     tolerance_param, methodTV, nonneg, device='cpu'):
     if device == 'cpu':
         return TV_FGP_CPU(inputData,
                      regularisation_parameter,
                      iterations, 
                      tolerance_param,
                      methodTV,
-                     nonneg,
-                     printM)
+                     nonneg)
     elif device == 'gpu' and gpu_enabled:
         return TV_FGP_GPU(inputData,
                      regularisation_parameter,
@@ -45,7 +44,7 @@ def FGP_TV(inputData, regularisation_parameter,iterations,
                      tolerance_param,
                      methodTV,
                      nonneg,
-                     printM)
+                     1)
     else:
         if not gpu_enabled and device == 'gpu':
             raise ValueError ('GPU is not available')
diff --git a/src/Python/setup-regularisers.py.in b/src/Python/setup-regularisers.py.in
index 82d9f9f..39b820a 100644
--- a/src/Python/setup-regularisers.py.in
+++ b/src/Python/setup-regularisers.py.in
@@ -44,7 +44,7 @@ extra_include_dirs += [os.path.join(".." , "Core"),
                        os.path.join(".." , "Core",  "regularisers_GPU" , "LLTROF" ) ,
                        os.path.join(".." , "Core",  "regularisers_GPU" , "NDF" ) ,
                        os.path.join(".." , "Core",  "regularisers_GPU" , "dTV_FGP" ) ,
-                       os.path.join(".." , "Core",  "regularisers_GPU" , "DIFF4th" ) ,
+                       os.path.join(".." , "Core",  "regularisers_GPU" , "Diff4th" ) ,
                        os.path.join(".." , "Core",  "regularisers_GPU" , "PatchSelect" ) ,
 						   "."]
 
diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx
index 11a0617..b7d029d 100644
--- a/src/Python/src/cpu_regularisers.pyx
+++ b/src/Python/src/cpu_regularisers.pyx
@@ -19,7 +19,7 @@ import numpy as np
 cimport numpy as np
 
 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);
+cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int dimX, int dimY, int dimZ);
 cdef extern float SB_TV_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int printM, int dimX, int dimY, int dimZ);
 cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, int dimX, int dimY, int dimZ);
 cdef extern float TGV_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY, int dimZ);
@@ -45,7 +45,7 @@ def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb, marching_ste
 
 def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
                      float regularisation_parameter,
-                     int iterationsNumb,                     
+                     int iterationsNumb,
                      float marching_step_parameter):
     cdef long dims[2]
     dims[0] = inputData.shape[0]
@@ -80,45 +80,46 @@ def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
 #********************** Total-variation FGP *********************#
 #****************************************************************#
 #******** Total-variation Fast-Gradient-Projection (FGP)*********#
-def TV_FGP_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM):
+def TV_FGP_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg):
     if inputData.ndim == 2:
-        return TV_FGP_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
+        return TV_FGP_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg)
     elif inputData.ndim == 3:
-        return TV_FGP_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
+        return TV_FGP_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg)
 
 def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
                      float regularisation_parameter,
                      int iterationsNumb, 
                      float tolerance_param,
                      int methodTV,
-                     int nonneg,
-                     int printM):
-                         
+                     int nonneg):
+
     cdef long dims[2]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
             np.zeros([dims[0],dims[1]], dtype='float32')
-                   
+        
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] infovec = \
+            np.ones([dims[0],dims[1]], dtype='float32')
+    
     #/* Run FGP-TV iterations for 2D data */
-    TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, 
+    TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0,0], regularisation_parameter, 
                        iterationsNumb, 
                        tolerance_param,
                        methodTV,
                        nonneg,
-                       printM,
                        dims[1],dims[0],1)
     
-    return outputData        
+    return (outputData,infovec)
             
 def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
                      float regularisation_parameter,
                      int iterationsNumb, 
                      float tolerance_param,
                      int methodTV,
-                     int nonneg,
-                     int printM):
+                     int nonneg):
+    
     cdef long dims[3]
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
@@ -126,16 +127,17 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
     
     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 */
-    TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter,
+    TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter,
                        iterationsNumb, 
                        tolerance_param,
                        methodTV,
                        nonneg,
-                       printM,
                        dims[2], dims[1], dims[0])
-    return outputData 
+    return (outputData,infovec)
 
 #***************************************************************#
 #********************** Total-variation SB *********************#