Merge remote-tracking branch 'remotes/origin/master' into newdirstructure

Conflicts:
	demos/demoMatlab_denoise.m
	demos/qualitymetrics.py

20 files changed, 827 insertions, 856 deletions

diff --git a/Readme.md b/Readme.md
index 187f8ac..3a39066 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=0.10.4
+	export CIL_VERSION=19.02
 	conda build Wrappers/Python/conda-recipe --numpy 1.12 --python 3.5 
 	conda install ccpi-regulariser=${CIL_VERSION} --use-local --force
 	cd demos/
diff --git a/Wrappers/Python/ccpi/supp/__init__.py b/Wrappers/Python/ccpi/supp/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/Wrappers/Python/ccpi/supp/__init__.py
diff --git a/Wrappers/Python/ccpi/supp/qualitymetrics.py b/Wrappers/Python/ccpi/supp/qualitymetrics.py
new file mode 100644
index 0000000..f44d832
--- /dev/null
+++ b/Wrappers/Python/ccpi/supp/qualitymetrics.py
@@ -0,0 +1,65 @@
+#!/usr/bin/env python2
+# -*- coding: utf-8 -*-
+"""
+A class for some standard image quality metrics
+"""
+import numpy as np
+
+class QualityTools:
+    def __init__(self, im1, im2):
+        if im1.size != im2.size:
+            print ('Error: Sizes of images/volumes are different')
+            raise SystemExit
+        self.im1 = im1 # image or volume - 1
+        self.im2 = im2 # image or volume - 2
+    def nrmse(self):
+        """ Normalised Root Mean Square Error """
+        rmse = np.sqrt(np.sum((self.im2 - self.im1) ** 2) / float(self.im1.size))
+        max_val = max(np.max(self.im1), np.max(self.im2))
+        min_val = min(np.min(self.im1), np.min(self.im2))
+        return 1 - (rmse / (max_val - min_val))
+    def rmse(self):
+        """ Root Mean Square Error """
+        rmse = np.sqrt(np.sum((self.im1 - self.im2) ** 2) / float(self.im1.size))
+        return rmse
+    def ssim(self, window, k=(0.01, 0.03), l=255):
+        from scipy.signal import fftconvolve
+        """See https://ece.uwaterloo.ca/~z70wang/research/ssim/"""
+        # Check if the window is smaller than the images.
+        for a, b in zip(window.shape, self.im1.shape):
+            if a > b:
+                return None, None
+        # Values in k must be positive according to the base implementation.
+        for ki in k:
+            if ki < 0:
+                return None, None
+    
+        c1 = (k[0] * l) ** 2
+        c2 = (k[1] * l) ** 2
+        window = window/np.sum(window)
+    
+        mu1 = fftconvolve(self.im1, window, mode='valid')
+        mu2 = fftconvolve(self.im2, window, mode='valid')
+        mu1_sq = mu1 * mu1
+        mu2_sq = mu2 * mu2
+        mu1_mu2 = mu1 * mu2
+        sigma1_sq = fftconvolve(self.im1 * self.im1, window, mode='valid') - mu1_sq
+        sigma2_sq = fftconvolve(self.im2 * self.im2, window, mode='valid') - mu2_sq
+        sigma12 = fftconvolve(self.im1 * self.im2, window, mode='valid') - mu1_mu2
+    
+        if c1 > 0 and c2 > 0:
+            num = (2 * mu1_mu2 + c1) * (2 * sigma12 + c2)
+            den = (mu1_sq + mu2_sq + c1) * (sigma1_sq + sigma2_sq + c2)
+            ssim_map = num / den
+        else:
+            num1 = 2 * mu1_mu2 + c1
+            num2 = 2 * sigma12 + c2
+            den1 = mu1_sq + mu2_sq + c1
+            den2 = sigma1_sq + sigma2_sq + c2
+            ssim_map = np.ones(np.shape(mu1))
+            index = (den1 * den2) > 0
+            ssim_map[index] = (num1[index] * num2[index]) / (den1[index] * den2[index])
+            index = (den1 != 0) & (den2 == 0)
+            ssim_map[index] = num1[index] / den1[index]    
+        mssim = ssim_map.mean()
+        return mssim, ssim_map
diff --git a/build/run.sh b/build/run.sh
index a8e5555..332d660 100644
--- a/build/run.sh
+++ b/build/run.sh
@@ -1,19 +1,22 @@
 #!/bin/bash  
 echo "Building CCPi-regularisation Toolkit using CMake"  
-# rm -r build
+rm -r build
 # Requires Cython, install it first: 
 # pip install cython
-# mkdir build
+mkdir build
 cd build/
 make clean
-# install Python modules only 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 only with CUDA
+# 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 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
 # cp install/lib/libcilreg.so install/python/ccpi/filters
-cd install/python
-export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:../lib
+#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="/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
diff --git a/demos/demoMatlab_3Ddenoise.m b/demos/demoMatlab_3Ddenoise.m
index cdd3117..cf2c88a 100644
--- a/demos/demoMatlab_3Ddenoise.m
+++ b/demos/demoMatlab_3Ddenoise.m
@@ -8,7 +8,7 @@ addpath(Path2);
 addpath(Path3);
 
 N = 512; 
-slices = 7;
+slices = 15;
 vol3D = zeros(N,N,slices, 'single');
 Ideal3D = zeros(N,N,slices, 'single');
 Im = double(imread('lena_gray_512.tif'))/255;  % loading image
@@ -17,9 +17,7 @@ vol3D(:,:,i) = Im + .05*randn(size(Im));
 Ideal3D(:,:,i) = Im;
 end
 vol3D(vol3D < 0) = 0;
-figure; imshow(vol3D(:,:,15), [0 1]); title('Noisy image');
-
-
+figure; imshow(vol3D(:,:,7), [0 1]); title('Noisy image');
 lambda_reg = 0.03; % regularsation parameter for all methods
 %%
 fprintf('Denoise a volume using the ROF-TV model (CPU) \n');
@@ -143,6 +141,16 @@ rmseTGV = RMSE(Ideal3D(:),u_tgv(:));
 fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV);
 figure; imshow(u_tgv(:,:,3), [0 1]); title('TGV denoised volume (CPU)');
 %%
+% fprintf('Denoise using the TGV model (GPU) \n');
+% lambda_TGV = 0.03; % regularisation parameter
+% alpha1 = 1.0; % parameter to control the first-order term
+% alpha0 = 2.0; % parameter to control the second-order term
+% iter_TGV = 500; % number of Primal-Dual iterations for TGV
+% tic; u_tgv_gpu = TGV_GPU(single(vol3D), lambda_TGV, alpha1, alpha0, iter_TGV); toc; 
+% rmseTGV = RMSE(Ideal3D(:),u_tgv_gpu(:));
+% fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV);
+% figure; imshow(u_tgv_gpu(:,:,3), [0 1]); title('TGV denoised volume (GPU)');
+%%
 %>>>>>>>>>>>>>> MULTI-CHANNEL priors <<<<<<<<<<<<<<< %
 fprintf('Denoise a volume using the FGP-dTV model (CPU) \n');
 
diff --git a/demos/demoMatlab_denoise.m b/demos/demoMatlab_denoise.m
index 2031853..5135129 100644
--- a/demos/demoMatlab_denoise.m
+++ b/demos/demoMatlab_denoise.m
@@ -5,7 +5,9 @@ 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);
-addpath(Path1); addpath(Path2); addpath(Path3);
+addpath(Path1);
+addpath(Path2);
+addpath(Path3);
 
 Im = double(imread('lena_gray_512.tif'))/255;  % loading image
 u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
@@ -29,7 +31,7 @@ figure; imshow(u_rof, [0 1]); title('ROF-TV denoised image (CPU)');
 % figure; imshow(u_rofG, [0 1]); title('ROF-TV denoised image (GPU)');
 %%
 fprintf('Denoise using the FGP-TV model (CPU) \n');
-iter_fgp = 1000; % number of FGP iterations
+iter_fgp = 1300; % number of FGP iterations
 epsil_tol =  1.0e-06; % 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
@@ -39,8 +41,8 @@ figure; imshow(u_fgp, [0 1]); title('FGP-TV denoised image (CPU)');
 
 %%
 % fprintf('Denoise using the FGP-TV model (GPU) \n');
-% iter_fgp = 1000; % number of FGP iterations
-% epsil_tol =  1.0e-05; % tolerance
+% iter_fgp = 1300; % number of FGP iterations
+% epsil_tol =  1.0e-06; % tolerance
 % tic; u_fgpG = FGP_TV_GPU(single(u0), lambda_reg, iter_fgp, epsil_tol); toc; 
 % figure; imshow(u_fgpG, [0 1]); title('FGP-TV denoised image (GPU)');
 %%
@@ -63,17 +65,17 @@ fprintf('Denoise using the TGV model (CPU) \n');
 lambda_TGV = 0.045; % regularisation parameter
 alpha1 = 1.0; % parameter to control the first-order term
 alpha0 = 2.0; % parameter to control the second-order term
-iter_TGV = 2000; % number of Primal-Dual iterations for TGV
+iter_TGV = 1500; % number of Primal-Dual iterations for TGV
 tic; u_tgv = TGV(single(u0), lambda_TGV, alpha1, alpha0, iter_TGV); toc; 
 rmseTGV = (RMSE(u_tgv(:),Im(:)));
 fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV);
 figure; imshow(u_tgv, [0 1]); title('TGV denoised image (CPU)');
-%%
+
 % fprintf('Denoise using the TGV model (GPU) \n');
 % lambda_TGV = 0.045; % regularisation parameter
 % alpha1 = 1.0; % parameter to control the first-order term
 % alpha0 = 2.0; % parameter to control the second-order term
-% iter_TGV = 2000; % number of Primal-Dual iterations for TGV
+% iter_TGV = 1500; % number of Primal-Dual iterations for TGV
 % tic; u_tgv_gpu = TGV_GPU(single(u0), lambda_TGV, alpha1, alpha0, iter_TGV); toc; 
 % rmseTGV_gpu = (RMSE(u_tgv_gpu(:),Im(:)));
 % fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV_gpu);
diff --git a/demos/demo_cpu_inpainters.py b/demos/demo_cpu_inpainters.py
index d07e74a..2e6ccf2 100644
--- a/demos/demo_cpu_inpainters.py
+++ b/demos/demo_cpu_inpainters.py
@@ -11,7 +11,7 @@ import os
 import timeit
 from scipy import io
 from ccpi.filters.regularisers import NDF_INP, NVM_INP
-from qualitymetrics import rmse
+from ccpi.supp.qualitymetrics import QualityTools
 ###############################################################################
 def printParametersToString(pars):
         txt = r''
@@ -85,9 +85,9 @@ ndf_inp_linear = NDF_INP(pars['input'],
               pars['number_of_iterations'],
               pars['time_marching_parameter'], 
               pars['penalty_type'])
-             
-rms = rmse(sino_full, ndf_inp_linear)
-pars['rmse'] = rms
+
+Qtools = QualityTools(sino_full, ndf_inp_linear)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -133,8 +133,9 @@ ndf_inp_nonlinear = NDF_INP(pars['input'],
               pars['time_marching_parameter'], 
               pars['penalty_type'])
              
-rms = rmse(sino_full, ndf_inp_nonlinear)
-pars['rmse'] = rms
+
+Qtools = QualityTools(sino_full, ndf_inp_nonlinear)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -174,8 +175,9 @@ start_time = timeit.default_timer()
               pars['SW_increment'],
               pars['number_of_iterations'])
              
-rms = rmse(sino_full, nvm_inp)
-pars['rmse'] = rms
+
+Qtools = QualityTools(sino_full, nvm_inp)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
diff --git a/demos/demo_cpu_regularisers.py b/demos/demo_cpu_regularisers.py
index 373502b..d34607a 100644
--- a/demos/demo_cpu_regularisers.py
+++ b/demos/demo_cpu_regularisers.py
@@ -14,7 +14,7 @@ import os
 import timeit
 from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, TNV, NDF, Diff4th
 from ccpi.filters.regularisers import PatchSelect, NLTV
-from qualitymetrics import rmse
+from ccpi.supp.qualitymetrics import QualityTools
 ###############################################################################
 def printParametersToString(pars):
         txt = r''
@@ -84,9 +84,9 @@ imgplot = plt.imshow(u0,cmap="gray")
 # set parameters
 pars = {'algorithm': ROF_TV, \
         'input' : u0,\
-        'regularisation_parameter':0.04,\
-        'number_of_iterations': 1200,\
-        'time_marching_parameter': 0.0025        
+        'regularisation_parameter':0.02,\
+        'number_of_iterations': 2000,\
+        'time_marching_parameter': 0.0025
         }
 print ("#############ROF TV CPU####################")
 start_time = timeit.default_timer()
@@ -94,8 +94,9 @@ rof_cpu = ROF_TV(pars['input'],
              pars['regularisation_parameter'],
              pars['number_of_iterations'],
              pars['time_marching_parameter'],'cpu')
-rms = rmse(Im, rof_cpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, rof_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -143,9 +144,9 @@ fgp_cpu = FGP_TV(pars['input'],
               pars['nonneg'],
               pars['printingOut'],'cpu')  
              
-             
-rms = rmse(Im, fgp_cpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, fgp_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -191,9 +192,8 @@ sb_cpu = SB_TV(pars['input'],
               pars['methodTV'],
               pars['printingOut'],'cpu')  
              
-             
-rms = rmse(Im, sb_cpu)
-pars['rmse'] = rms
+Qtools = QualityTools(Im, sb_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -240,8 +240,8 @@ tgv_cpu = TGV(pars['input'],
               pars['LipshitzConstant'],'cpu')
              
              
-rms = rmse(Im, tgv_cpu)
-pars['rmse'] = rms
+Qtools = QualityTools(Im, tgv_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -286,8 +286,8 @@ lltrof_cpu = LLT_ROF(pars['input'],
               pars['number_of_iterations'],
               pars['time_marching_parameter'],'cpu')
 
-rms = rmse(Im, lltrof_cpu)
-pars['rmse'] = rms
+Qtools = QualityTools(Im, lltrof_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -335,8 +335,8 @@ ndf_cpu = NDF(pars['input'],
               pars['time_marching_parameter'], 
               pars['penalty_type'],'cpu')  
              
-rms = rmse(Im, ndf_cpu)
-pars['rmse'] = rms
+Qtools = QualityTools(Im, ndf_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -380,8 +380,8 @@ diff4_cpu = Diff4th(pars['input'],
               pars['number_of_iterations'],
               pars['time_marching_parameter'],'cpu')
              
-rms = rmse(Im, diff4_cpu)
-pars['rmse'] = rms
+Qtools = QualityTools(Im, diff4_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -452,8 +452,8 @@ nltv_cpu = NLTV(pars2['input'],
               pars2['regularisation_parameter'],
               pars2['iterations'])
 
-rms = rmse(Im, nltv_cpu)
-pars['rmse'] = rms
+Qtools = QualityTools(Im, nltv_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -505,8 +505,8 @@ fgp_dtv_cpu = FGP_dTV(pars['input'],
               pars['nonneg'],
               pars['printingOut'],'cpu')
              
-rms = rmse(Im, fgp_dtv_cpu)
-pars['rmse'] = rms
+Qtools = QualityTools(Im, fgp_dtv_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -554,9 +554,9 @@ tnv_cpu = TNV(pars['input'],
               pars['regularisation_parameter'],
               pars['number_of_iterations'],
               pars['tolerance_constant'])
-             
-rms = rmse(idealVol, tnv_cpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(idealVol, tnv_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
diff --git a/demos/demo_cpu_regularisers3D.py b/demos/demo_cpu_regularisers3D.py
index 56baf13..fd6c545 100644
--- a/demos/demo_cpu_regularisers3D.py
+++ b/demos/demo_cpu_regularisers3D.py
@@ -13,7 +13,7 @@ import numpy as np
 import os
 import timeit
 from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th
-from qualitymetrics import rmse
+from ccpi.supp.qualitymetrics import QualityTools
 ###############################################################################
 def printParametersToString(pars):
         txt = r''
@@ -104,8 +104,9 @@ rof_cpu3D = ROF_TV(pars['input'],
              pars['regularisation_parameter'],
              pars['number_of_iterations'],
              pars['time_marching_parameter'],'cpu')
-rms = rmse(idealVol, rof_cpu3D)
-pars['rmse'] = rms
+
+Qtools = QualityTools(idealVol, rof_cpu3D)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -153,9 +154,8 @@ fgp_cpu3D = FGP_TV(pars['input'],
               pars['nonneg'],
               pars['printingOut'],'cpu')  
              
-             
-rms = rmse(idealVol, fgp_cpu3D)
-pars['rmse'] = rms
+Qtools = QualityTools(idealVol, fgp_cpu3D)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -201,8 +201,10 @@ sb_cpu3D = SB_TV(pars['input'],
               pars['methodTV'],
               pars['printingOut'],'cpu')
              
-rms = rmse(idealVol, sb_cpu3D)
-pars['rmse'] = rms
+
+Qtools = QualityTools(idealVol, sb_cpu3D)
+pars['rmse'] = Qtools.rmse()
+
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -246,8 +248,9 @@ lltrof_cpu3D = LLT_ROF(pars['input'],
               pars['number_of_iterations'],
               pars['time_marching_parameter'],'cpu')
 
-rms = rmse(idealVol, lltrof_cpu3D)
-pars['rmse'] = rms
+
+Qtools = QualityTools(idealVol, lltrof_cpu3D)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -294,8 +297,8 @@ tgv_cpu3D = TGV(pars['input'],
               pars['LipshitzConstant'],'cpu')
              
 
-rms = rmse(idealVol, tgv_cpu3D)
-pars['rmse'] = rms
+Qtools = QualityTools(idealVol, tgv_cpu3D)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -341,8 +344,9 @@ ndf_cpu3D = NDF(pars['input'],
               pars['time_marching_parameter'], 
               pars['penalty_type'])  
              
-rms = rmse(idealVol, ndf_cpu3D)
-pars['rmse'] = rms
+
+Qtools = QualityTools(idealVol, ndf_cpu3D)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -386,8 +390,9 @@ diff4th_cpu3D = Diff4th(pars['input'],
               pars['number_of_iterations'],
               pars['time_marching_parameter'])  
              
-rms = rmse(idealVol, diff4th_cpu3D)
-pars['rmse'] = rms
+
+Qtools = QualityTools(idealVol, diff4th_cpu3D)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -439,9 +444,9 @@ fgp_dTV_cpu3D = FGP_dTV(pars['input'],
               pars['nonneg'],
               pars['printingOut'],'cpu')
              
-             
-rms = rmse(idealVol, fgp_dTV_cpu3D)
-pars['rmse'] = rms
+
+Qtools = QualityTools(idealVol, fgp_dTV_cpu3D)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
diff --git a/demos/demo_cpu_vs_gpu_regularisers.py b/demos/demo_cpu_vs_gpu_regularisers.py
index 5ce8da4..e1eb91f 100644
--- a/demos/demo_cpu_vs_gpu_regularisers.py
+++ b/demos/demo_cpu_vs_gpu_regularisers.py
@@ -14,7 +14,7 @@ import os
 import timeit
 from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th
 from ccpi.filters.regularisers import PatchSelect
-from qualitymetrics import rmse
+from ccpi.supp.qualitymetrics import QualityTools
 ###############################################################################
 def printParametersToString(pars):
         txt = r''
@@ -76,8 +76,9 @@ rof_cpu = ROF_TV(pars['input'],
              pars['regularisation_parameter'],
              pars['number_of_iterations'],
              pars['time_marching_parameter'],'cpu')
-rms = rmse(Im, rof_cpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, rof_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -98,9 +99,10 @@ rof_gpu = ROF_TV(pars['input'],
                      pars['regularisation_parameter'],
                      pars['number_of_iterations'], 
                      pars['time_marching_parameter'],'gpu')
-                     
-rms = rmse(Im, rof_gpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, rof_gpu)
+pars['rmse'] = Qtools.rmse()
+
 pars['algorithm'] = ROF_TV
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -162,9 +164,9 @@ fgp_cpu = FGP_TV(pars['input'],
               pars['nonneg'],
               pars['printingOut'],'cpu')  
              
-             
-rms = rmse(Im, fgp_cpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, fgp_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -189,9 +191,10 @@ fgp_gpu = FGP_TV(pars['input'],
               pars['methodTV'],
               pars['nonneg'],
               pars['printingOut'],'gpu')
-                                   
-rms = rmse(Im, fgp_gpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, fgp_gpu)
+pars['rmse'] = Qtools.rmse()
+
 pars['algorithm'] = FGP_TV
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -251,9 +254,9 @@ sb_cpu = SB_TV(pars['input'],
               pars['methodTV'],
               pars['printingOut'],'cpu')  
              
-             
-rms = rmse(Im, sb_cpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, sb_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -277,9 +280,9 @@ sb_gpu = SB_TV(pars['input'],
               pars['tolerance_constant'], 
               pars['methodTV'],
               pars['printingOut'],'gpu')
-                                   
-rms = rmse(Im, sb_gpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, sb_gpu)
+pars['rmse'] = Qtools.rmse()
 pars['algorithm'] = SB_TV
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -337,8 +340,8 @@ tgv_cpu = TGV(pars['input'],
               pars['number_of_iterations'],
               pars['LipshitzConstant'],'cpu')
              
-rms = rmse(Im, tgv_cpu)
-pars['rmse'] = rms
+Qtools = QualityTools(Im, tgv_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -362,8 +365,8 @@ tgv_gpu = TGV(pars['input'],
               pars['number_of_iterations'],
               pars['LipshitzConstant'],'gpu')
                                    
-rms = rmse(Im, tgv_gpu)
-pars['rmse'] = rms
+Qtools = QualityTools(Im, tgv_gpu)
+pars['rmse'] = Qtools.rmse()
 pars['algorithm'] = TGV
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -419,8 +422,8 @@ lltrof_cpu = LLT_ROF(pars['input'],
               pars['number_of_iterations'],
               pars['time_marching_parameter'],'cpu')
 
-rms = rmse(Im, lltrof_cpu)
-pars['rmse'] = rms
+Qtools = QualityTools(Im, lltrof_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -443,8 +446,9 @@ lltrof_gpu = LLT_ROF(pars['input'],
               pars['number_of_iterations'],
               pars['time_marching_parameter'],'gpu')
 
-rms = rmse(Im, lltrof_gpu)
-pars['rmse'] = rms
+Qtools = QualityTools(Im, lltrof_gpu)
+pars['rmse'] = Qtools.rmse()
+
 pars['algorithm'] = LLT_ROF
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -501,9 +505,9 @@ ndf_cpu = NDF(pars['input'],
               pars['number_of_iterations'],
               pars['time_marching_parameter'], 
               pars['penalty_type'],'cpu')
-             
-rms = rmse(Im, ndf_cpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, ndf_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -527,9 +531,9 @@ ndf_gpu = NDF(pars['input'],
               pars['number_of_iterations'],
               pars['time_marching_parameter'], 
               pars['penalty_type'],'gpu')
-             
-rms = rmse(Im, ndf_gpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, ndf_gpu)
+pars['rmse'] = Qtools.rmse()
 pars['algorithm'] = NDF
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -585,9 +589,9 @@ diff4th_cpu = Diff4th(pars['input'],
               pars['edge_parameter'], 
               pars['number_of_iterations'],
               pars['time_marching_parameter'],'cpu')
-             
-rms = rmse(Im, diff4th_cpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, diff4th_cpu)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -609,9 +613,9 @@ diff4th_gpu = Diff4th(pars['input'],
               pars['edge_parameter'], 
               pars['number_of_iterations'],
               pars['time_marching_parameter'], 'gpu')
-             
-rms = rmse(Im, diff4th_gpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, diff4th_gpu)
+pars['rmse'] = Qtools.rmse()
 pars['algorithm'] = Diff4th
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -675,10 +679,10 @@ fgp_dtv_cpu = FGP_dTV(pars['input'],
               pars['methodTV'],
               pars['nonneg'],
               pars['printingOut'],'cpu')
-             
-             
-rms = rmse(Im, fgp_dtv_cpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, fgp_dtv_cpu)
+pars['rmse'] = Qtools.rmse()
+
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -704,8 +708,8 @@ fgp_dtv_gpu = FGP_dTV(pars['input'],
               pars['methodTV'],
               pars['nonneg'],
               pars['printingOut'],'gpu')
-rms = rmse(Im, fgp_dtv_gpu)
-pars['rmse'] = rms
+Qtools = QualityTools(Im, fgp_dtv_gpu)
+pars['rmse'] = Qtools.rmse()
 pars['algorithm'] = FGP_dTV
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
diff --git a/demos/demo_gpu_regularisers.py b/demos/demo_gpu_regularisers.py
index bc9baf2..89bb948 100644
--- a/demos/demo_gpu_regularisers.py
+++ b/demos/demo_gpu_regularisers.py
@@ -14,7 +14,7 @@ import os
 import timeit
 from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th
 from ccpi.filters.regularisers import PatchSelect, NLTV
-from qualitymetrics import rmse
+from ccpi.supp.qualitymetrics import QualityTools
 ###############################################################################
 def printParametersToString(pars):
         txt = r''
@@ -93,9 +93,9 @@ rof_gpu = ROF_TV(pars['input'],
                      pars['regularisation_parameter'],
                      pars['number_of_iterations'], 
                      pars['time_marching_parameter'],'gpu')
-                     
-rms = rmse(Im, rof_gpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, rof_gpu)
+pars['rmse'] = Qtools.rmse()
 pars['algorithm'] = ROF_TV
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -142,9 +142,8 @@ fgp_gpu = FGP_TV(pars['input'],
               pars['methodTV'],
               pars['nonneg'],
               pars['printingOut'],'gpu')
-                                   
-rms = rmse(Im, fgp_gpu)
-pars['rmse'] = rms
+Qtools = QualityTools(Im, fgp_gpu)
+pars['rmse'] = Qtools.rmse()
 pars['algorithm'] = FGP_TV
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -189,9 +188,9 @@ sb_gpu = SB_TV(pars['input'],
               pars['tolerance_constant'], 
               pars['methodTV'],
               pars['printingOut'],'gpu')
-                                   
-rms = rmse(Im, sb_gpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, sb_gpu)
+pars['rmse'] = Qtools.rmse()
 pars['algorithm'] = SB_TV
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -236,11 +235,9 @@ tgv_gpu = TGV(pars['input'],
               pars['alpha0'],
               pars['number_of_iterations'],
               pars['LipshitzConstant'],'gpu')  
-             
-             
-rms = rmse(Im, tgv_gpu)
-pars['rmse'] = rms
 
+Qtools = QualityTools(Im, tgv_gpu)
+pars['rmse'] = Qtools.rmse()
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
 print (txtstr)
@@ -284,10 +281,8 @@ lltrof_gpu = LLT_ROF(pars['input'],
               pars['number_of_iterations'],
               pars['time_marching_parameter'],'gpu')
              
-             
-rms = rmse(Im, lltrof_gpu)
-pars['rmse'] = rms
-
+Qtools = QualityTools(Im, lltrof_gpu)
+pars['rmse'] = Qtools.rmse()
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
 print (txtstr)
@@ -331,9 +326,9 @@ ndf_gpu = NDF(pars['input'],
               pars['number_of_iterations'],
               pars['time_marching_parameter'], 
               pars['penalty_type'],'gpu')  
-             
-rms = rmse(Im, ndf_gpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, ndf_gpu)
+pars['rmse'] = Qtools.rmse()
 pars['algorithm'] = NDF
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -376,10 +371,10 @@ diff4_gpu = Diff4th(pars['input'],
               pars['edge_parameter'], 
               pars['number_of_iterations'],
               pars['time_marching_parameter'],'gpu')
-             
-rms = rmse(Im, diff4_gpu)
-pars['rmse'] = rms
 
+Qtools = QualityTools(Im, diff4_gpu)
+pars['algorithm'] = Diff4th
+pars['rmse'] = Qtools.rmse()
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
 print (txtstr)
@@ -449,9 +444,8 @@ nltv_cpu = NLTV(pars2['input'],
               pars2['regularisation_parameter'],
               pars2['iterations'])
 
-rms = rmse(Im, nltv_cpu)
-pars['rmse'] = rms
-
+Qtools = QualityTools(Im, nltv_cpu)
+pars['rmse'] = Qtools.rmse()
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
 print (txtstr)
@@ -500,9 +494,9 @@ fgp_dtv_gpu = FGP_dTV(pars['input'],
               pars['methodTV'],
               pars['nonneg'],
               pars['printingOut'],'gpu')
-                                   
-rms = rmse(Im, fgp_dtv_gpu)
-pars['rmse'] = rms
+
+Qtools = QualityTools(Im, fgp_dtv_gpu)
+pars['rmse'] = Qtools.rmse()
 pars['algorithm'] = FGP_dTV
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
diff --git a/demos/demo_gpu_regularisers3D.py b/demos/demo_gpu_regularisers3D.py
index 2f49cb9..be16921 100644
--- a/demos/demo_gpu_regularisers3D.py
+++ b/demos/demo_gpu_regularisers3D.py
@@ -13,7 +13,7 @@ import numpy as np
 import os
 import timeit
 from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th
-from qualitymetrics import rmse
+from ccpi.supp.qualitymetrics import QualityTools
 ###############################################################################
 def printParametersToString(pars):
         txt = r''
@@ -111,9 +111,9 @@ rof_gpu3D = ROF_TV(pars['input'],
              pars['regularisation_parameter'],
              pars['number_of_iterations'],
              pars['time_marching_parameter'],'gpu')
-rms = rmse(idealVol, rof_gpu3D)
-pars['rmse'] = rms
 
+Qtools = QualityTools(idealVol, rof_gpu3D)
+pars['rmse'] = Qtools.rmse()
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
 print (txtstr)
@@ -159,9 +159,8 @@ fgp_gpu3D = FGP_TV(pars['input'],
               pars['nonneg'],
               pars['printingOut'],'gpu')
 
-rms = rmse(idealVol, fgp_gpu3D)
-pars['rmse'] = rms
-
+Qtools = QualityTools(idealVol, fgp_gpu3D)
+pars['rmse'] = Qtools.rmse()
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
 print (txtstr)
@@ -206,8 +205,8 @@ sb_gpu3D = SB_TV(pars['input'],
               pars['methodTV'],
               pars['printingOut'],'gpu')
 
-rms = rmse(idealVol, sb_gpu3D)
-pars['rmse'] = rms
+Qtools = QualityTools(idealVol, sb_gpu3D)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -250,8 +249,8 @@ lltrof_gpu3D = LLT_ROF(pars['input'],
               pars['number_of_iterations'],
               pars['time_marching_parameter'],'gpu')
 
-rms = rmse(idealVol, lltrof_gpu3D)
-pars['rmse'] = rms
+Qtools = QualityTools(idealVol, lltrof_gpu3D)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
@@ -296,11 +295,9 @@ tgv_gpu3D = TGV(pars['input'],
               pars['alpha0'],
               pars['number_of_iterations'],
               pars['LipshitzConstant'],'gpu')
-             
-
-rms = rmse(idealVol, tgv_gpu3D)
-pars['rmse'] = rms
 
+Qtools = QualityTools(idealVol, tgv_gpu3D)
+pars['rmse'] = Qtools.rmse()
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
 print (txtstr)
@@ -344,9 +341,8 @@ ndf_gpu3D = NDF(pars['input'],
               pars['time_marching_parameter'], 
               pars['penalty_type'],'gpu')
 
-rms = rmse(idealVol, ndf_gpu3D)
-pars['rmse'] = rms
-
+Qtools = QualityTools(idealVol, ndf_gpu3D)
+pars['rmse'] = Qtools.rmse()
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
 print (txtstr)
@@ -388,10 +384,9 @@ diff4_gpu3D = Diff4th(pars['input'],
               pars['edge_parameter'], 
               pars['number_of_iterations'],
               pars['time_marching_parameter'],'gpu')
-             
-rms = rmse(idealVol, diff4_gpu3D)
-pars['rmse'] = rms
 
+Qtools = QualityTools(idealVol, diff4_gpu3D)
+pars['rmse'] = Qtools.rmse()
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
 print (txtstr)
@@ -442,8 +437,8 @@ fgp_dTV_gpu3D = FGP_dTV(pars['input'],
               pars['nonneg'],
               pars['printingOut'],'gpu')
 
-rms = rmse(idealVol, fgp_dTV_gpu3D)
-pars['rmse'] = rms
+Qtools = QualityTools(idealVol, fgp_dTV_gpu3D)
+pars['rmse'] = Qtools.rmse()
 
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
diff --git a/demos/qualitymetrics.py b/demos/qualitymetrics.py
index 850829e..e69de29 100644
--- a/demos/qualitymetrics.py
+++ b/demos/qualitymetrics.py
@@ -1,18 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Wed Feb 21 13:34:32 2018
-# quality metrics
-@authors: Daniil Kazantsev, Edoardo Pasca
-"""
-import numpy as np
-
-def nrmse(im1, im2):
-    rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(im1.size))
-    max_val = max(np.max(im1), np.max(im2))
-    min_val = min(np.min(im1), np.min(im2))
-    return 1 - (rmse / (max_val - min_val))
-    
-def rmse(im1, im2):
-    rmse = np.sqrt(np.sum((im1 - im2) ** 2) / float(im1.size))
-    return rmse
diff --git a/src/Core/regularisers_CPU/TGV_core.c b/src/Core/regularisers_CPU/TGV_core.c
index 805c3d4..136e0bd 100644
--- a/src/Core/regularisers_CPU/TGV_core.c
+++ b/src/Core/regularisers_CPU/TGV_core.c
@@ -1,25 +1,25 @@
 /*
-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.
-*/
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2019 Daniil Kazantsev
+ * Copyright 2019 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
 
 #include "TGV_core.h"
 
-/* C-OMP implementation of Primal-Dual denoising method for 
+/* C-OMP implementation of Primal-Dual denoising method for
  * Total Generilized Variation (TGV)-L2 model [1] (2D/3D case)
  *
  * Input Parameters:
@@ -29,44 +29,44 @@ limitations under the License.
  * 4. parameter to control the second-order term (alpha0)
  * 5. Number of Chambolle-Pock (Primal-Dual) iterations
  * 6. Lipshitz constant (default is 12)
- * 
+ *
  * Output:
  * Filtered/regularised image/volume
  *
  * References:
  * [1] K. Bredies "Total Generalized Variation"
- * 
+ *
  */
- 
+
 float TGV_main(float *U0, float *U, float lambda, float alpha1, float alpha0, int iter, float L2, int dimX, int dimY, int dimZ)
 {
-	long DimTotal;
-	int ll;
-	float *U_old, *P1, *P2, *Q1, *Q2, *Q3, *V1, *V1_old, *V2, *V2_old, tau, sigma;
-
-	DimTotal = (long)(dimX*dimY*dimZ);
-	copyIm(U0, U, (long)(dimX), (long)(dimY), (long)(dimZ)); /* initialize */
-        tau = pow(L2,-0.5);
-        sigma = pow(L2,-0.5);
-
-        /* dual variables */
-        P1 = calloc(DimTotal, sizeof(float));
-        P2 = calloc(DimTotal, sizeof(float));
-        
-        Q1 = calloc(DimTotal, sizeof(float));
-        Q2 = calloc(DimTotal, sizeof(float));
-        Q3 = calloc(DimTotal, sizeof(float));
-        
-        U_old = calloc(DimTotal, sizeof(float));
+    long DimTotal;
+    int ll;
+    float *U_old, *P1, *P2, *Q1, *Q2, *Q3, *V1, *V1_old, *V2, *V2_old, tau, sigma;
+    
+    DimTotal = (long)(dimX*dimY*dimZ);
+    copyIm(U0, U, (long)(dimX), (long)(dimY), (long)(dimZ)); /* initialize */
+    tau = pow(L2,-0.5);
+    sigma = pow(L2,-0.5);
+    
+    /* dual variables */
+    P1 = calloc(DimTotal, sizeof(float));
+    P2 = calloc(DimTotal, sizeof(float));
+    
+    Q1 = calloc(DimTotal, sizeof(float));
+    Q2 = calloc(DimTotal, sizeof(float));
+    Q3 = calloc(DimTotal, sizeof(float));
+    
+    U_old = calloc(DimTotal, sizeof(float));
+    
+    V1 = calloc(DimTotal, sizeof(float));
+    V1_old = calloc(DimTotal, sizeof(float));
+    V2 = calloc(DimTotal, sizeof(float));
+    V2_old = calloc(DimTotal, sizeof(float));
+    
+    if (dimZ == 1) {
+        /*2D case*/
         
-        V1 = calloc(DimTotal, sizeof(float));
-        V1_old = calloc(DimTotal, sizeof(float));
-        V2 = calloc(DimTotal, sizeof(float));
-        V2_old = calloc(DimTotal, sizeof(float));
-	
-	if (dimZ == 1) {
-	/*2D case*/
-	
         /* Primal-dual iterations begin here */
         for(ll = 0; ll < iter; ll++) {
             
@@ -102,8 +102,8 @@ float TGV_main(float *U0, float *U, float lambda, float alpha1, float alpha0, in
             newU(V1, V1_old, (long)(dimX), (long)(dimY));
             newU(V2, V2_old, (long)(dimX), (long)(dimY));
         } /*end of iterations*/
-        	}
-        else {
+    }
+    else {
         /*3D case*/
         float *P3, *Q4, *Q5, *Q6, *V3, *V3_old;
         
@@ -114,7 +114,7 @@ float TGV_main(float *U0, float *U, float lambda, float alpha1, float alpha0, in
         V3 = calloc(DimTotal, sizeof(float));
         V3_old = calloc(DimTotal, sizeof(float));
         
-         /* Primal-dual iterations begin here */
+        /* Primal-dual iterations begin here */
         for(ll = 0; ll < iter; ll++) {
             
             /* Calculate Dual Variable P */
@@ -145,21 +145,20 @@ float TGV_main(float *U0, float *U, float lambda, float alpha1, float alpha0, in
             UpdV_3D(V1, V2, V3, P1, P2, P3, Q1, Q2, Q3, Q4, Q5, Q6, (long)(dimX), (long)(dimY), (long)(dimZ), tau);
             
             /*get new V*/
-            newU3D_3Ar(V1, V2, V3, V1_old, V2_old, V3_old, (long)(dimX), (long)(dimY), (long)(dimZ));           
-	        } /*end of iterations*/
+            newU3D_3Ar(V1, V2, V3, V1_old, V2_old, V3_old, (long)(dimX), (long)(dimY), (long)(dimZ));
+        } /*end of iterations*/
         free(P3);free(Q4);free(Q5);free(Q6);free(V3);free(V3_old);
-        }     
-
+    }
+    
     /*freeing*/
     free(P1);free(P2);free(Q1);free(Q2);free(Q3);free(U_old);
     free(V1);free(V2);free(V1_old);free(V2_old);
-	return *U;
+    return *U;
 }
 
 /********************************************************************/
 /***************************2D Functions*****************************/
 /********************************************************************/
-
 /*Calculating dual variable P (using forward differences)*/
 float DualP_2D(float *U, float *V1, float *V2, float *P1, float *P2, long dimX, long dimY, float sigma)
 {
@@ -167,12 +166,13 @@ float DualP_2D(float *U, float *V1, float *V2, float *P1, float *P2, long dimX,
 #pragma omp parallel for shared(U,V1,V2,P1,P2) private(i,j,index)
     for(i=0; i<dimX; i++) {
         for(j=0; j<dimY; j++) {
-			 index = j*dimX+i;
+            index = j*dimX+i;
             /* symmetric boundary conditions (Neuman) */
-            if (i == dimX-1) P1[index] += sigma*((U[j*dimX+(i-1)] - U[index]) - V1[index]); 
-            else P1[index] += sigma*((U[j*dimX+(i+1)] - U[index])  - V1[index]); 
-            if (j == dimY-1) P2[index] += sigma*((U[(j-1)*dimX+i] - U[index])  - V2[index]);
+            if (i == dimX-1) P1[index] += sigma*(-V1[index]);
+            else P1[index] += sigma*((U[j*dimX+(i+1)] - U[index])  - V1[index]);
+            if (j == dimY-1) P2[index] += sigma*(-V2[index]);
             else  P2[index] += sigma*((U[(j+1)*dimX+i] - U[index])  - V2[index]);
+            
         }}
     return 1;
 }
@@ -184,7 +184,7 @@ float ProjP_2D(float *P1, float *P2, long dimX, long dimY, float alpha1)
 #pragma omp parallel for shared(P1,P2) private(i,j,index,grad_magn)
     for(i=0; i<dimX; i++) {
         for(j=0; j<dimY; j++) {
-	    index = j*dimX+i;
+            index = j*dimX+i;
             grad_magn = (sqrtf(pow(P1[index],2) + pow(P2[index],2)))/alpha1;
             if (grad_magn > 1.0f) {
                 P1[index] /= grad_magn;
@@ -201,8 +201,8 @@ float DualQ_2D(float *V1, float *V2, float *Q1, float *Q2, float *Q3, long dimX,
 #pragma omp parallel for shared(Q1,Q2,Q3,V1,V2) private(i,j,index,q1,q2,q11,q22)
     for(i=0; i<dimX; i++) {
         for(j=0; j<dimY; j++) {
-    	    index = j*dimX+i;
-    	    q1 = 0.0f; q11 = 0.0f; q2 = 0.0f; q22 = 0.0f;
+            index = j*dimX+i;
+            q1 = 0.0f; q11 = 0.0f; q2 = 0.0f; q22 = 0.0f;
             /* boundary conditions (Neuman) */
             if (i != dimX-1){
                 q1 = V1[j*dimX+(i+1)] - V1[index];
@@ -225,7 +225,7 @@ float ProjQ_2D(float *Q1, float *Q2, float *Q3, long dimX, long dimY, float alph
 #pragma omp parallel for shared(Q1,Q2,Q3) private(i,j,index,grad_magn)
     for(i=0; i<dimX; i++) {
         for(j=0; j<dimY; j++) {
-	   index = j*dimX+i;
+            index = j*dimX+i;
             grad_magn = sqrtf(pow(Q1[index],2) + pow(Q2[index],2) + 2*pow(Q3[index],2));
             grad_magn = grad_magn/alpha0;
             if (grad_magn > 1.0f) {
@@ -236,7 +236,7 @@ float ProjQ_2D(float *Q1, float *Q2, float *Q3, long dimX, long dimY, float alph
         }}
     return 1;
 }
-/* Divergence and projection for P*/
+/* Divergence and projection for P (backward differences)*/
 float DivProjP_2D(float *U, float *U0, float *P1, float *P2, long dimX, long dimY, float lambda, float tau)
 {
     long i,j,index;
@@ -244,11 +244,16 @@ float DivProjP_2D(float *U, float *U0, float *P1, float *P2, long dimX, long dim
 #pragma omp parallel for shared(U,U0,P1,P2) private(i,j,index,P_v1,P_v2,div)
     for(i=0; i<dimX; i++) {
         for(j=0; j<dimY; j++) {
-	    index = j*dimX+i;
+            index = j*dimX+i;
+            
             if (i == 0) P_v1 = P1[index];
+            else if (i == dimX-1) P_v1 = -P1[j*dimX+(i-1)];
             else P_v1 = P1[index] - P1[j*dimX+(i-1)];
+            
             if (j == 0) P_v2 = P2[index];
-            else  P_v2 = P2[index] - P2[(j-1)*dimX+i];
+            else if (j == dimY-1) P_v2 = -P2[(j-1)*dimX+i];
+            else P_v2 = P2[index] - P2[(j-1)*dimX+i];
+            
             div = P_v1 + P_v2;
             U[index] = (lambda*(U[index] + tau*div) + tau*U0[index])/(lambda + tau);
         }}
@@ -262,7 +267,7 @@ float newU(float *U, float *U_old, long dimX, long dimY)
     for(i=0; i<dimX*dimY; i++) U[i] = 2*U[i] - U_old[i];
     return *U;
 }
-/*get update for V*/
+/*get update for V (backward differences)*/
 float UpdV_2D(float *V1, float *V2, float *P1, float *P2, float *Q1, float *Q2, float *Q3, long dimX, long dimY, float tau)
 {
     long i, j, index;
@@ -270,17 +275,30 @@ float UpdV_2D(float *V1, float *V2, float *P1, float *P2, float *Q1, float *Q2,
 #pragma omp parallel for shared(V1,V2,P1,P2,Q1,Q2,Q3) private(i, j, index, q1, q3_x, q3_y, q2, div1, div2)
     for(i=0; i<dimX; i++) {
         for(j=0; j<dimY; j++) {
-	    index = j*dimX+i;
-              q2 = 0.0f;  q3_y = 0.0f; q1 = 0.0f; q3_x = 0.0;
+            index = j*dimX+i;
+            
             /* boundary conditions (Neuman) */
-            if (i != 0) {
+            if (i == 0) {
+                q1 = Q1[index];
+                q3_x = Q3[index]; }
+            else if (i == dimX-1) {
+                q1 = -Q1[j*dimX+(i-1)];
+                q3_x = -Q3[j*dimX+(i-1)];  }
+            else {
                 q1 = Q1[index] - Q1[j*dimX+(i-1)];
-                q3_x = Q3[index] - Q3[j*dimX+(i-1)];
-            }
-            if (j != 0) {
+                q3_x = Q3[index] - Q3[j*dimX+(i-1)];  }
+            
+            if (j == 0) {
+                q2 = Q2[index];
+                q3_y = Q3[index]; }
+            else if (j == dimY-1) {
+                q2 = -Q2[(j-1)*dimX+i];
+                q3_y = -Q3[(j-1)*dimX+i]; }
+            else {
                 q2 = Q2[index] - Q2[(j-1)*dimX+i];
-                q3_y = Q3[index] - Q3[(j-1)*dimX+i];
-            }
+                q3_y = Q3[index] - Q3[(j-1)*dimX+i]; }
+            
+            
             div1 = q1 + q3_y;
             div2 = q3_x + q2;
             V1[index] += tau*(P1[index] + div1);
@@ -299,16 +317,16 @@ float DualP_3D(float *U, float *V1, float *V2, float *V3, float *P1, float *P2,
 #pragma omp parallel for shared(U,V1,V2,V3,P1,P2,P3) private(i,j,k,index)
     for(i=0; i<dimX; i++) {
         for(j=0; j<dimY; j++) {
-          for(k=0; k<dimZ; k++) {             	   
-    	   index = (dimX*dimY)*k + j*dimX+i;    	   
-            /* symmetric boundary conditions (Neuman) */
-            if (i == dimX-1) P1[index] += sigma*((U[(dimX*dimY)*k + j*dimX+(i-1)] - U[index]) - V1[index]); 
-            else P1[index] += sigma*((U[(dimX*dimY)*k + j*dimX+(i+1)] - U[index])  - V1[index]); 
-            if (j == dimY-1) P2[index] += sigma*((U[(dimX*dimY)*k + (j-1)*dimX+i] - U[index])  - V2[index]);
-            else  P2[index] += sigma*((U[(dimX*dimY)*k + (j+1)*dimX+i] - U[index])  - V2[index]);
-            if (k == dimZ-1) P3[index] += sigma*((U[(dimX*dimY)*(k-1) + j*dimX+i] - U[index])  - V3[index]);
-            else  P3[index] += sigma*((U[(dimX*dimY)*(k+1) + j*dimX+i] - U[index])  - V3[index]);
-        }}}
+            for(k=0; k<dimZ; k++) {
+                index = (dimX*dimY)*k + j*dimX+i;
+                /* symmetric boundary conditions (Neuman) */
+                if (i == dimX-1) P1[index] += sigma*(-V1[index]);
+                else P1[index] += sigma*((U[(dimX*dimY)*k + j*dimX+(i+1)] - U[index])  - V1[index]);
+                if (j == dimY-1) P2[index] += sigma*(-V2[index]);
+                else  P2[index] += sigma*((U[(dimX*dimY)*k + (j+1)*dimX+i] - U[index])  - V2[index]);
+                if (k == dimZ-1) P3[index] += sigma*(-V3[index]);
+                else  P3[index] += sigma*((U[(dimX*dimY)*(k+1) + j*dimX+i] - U[index])  - V3[index]);
+            }}}
     return 1;
 }
 /*Projection onto convex set for P*/
@@ -319,15 +337,15 @@ float ProjP_3D(float *P1, float *P2, float *P3, long dimX, long dimY, long dimZ,
 #pragma omp parallel for shared(P1,P2,P3) private(i,j,k,index,grad_magn)
     for(i=0; i<dimX; i++) {
         for(j=0; j<dimY; j++) {
-	  for(k=0; k<dimZ; k++) {   	
-   	    index = (dimX*dimY)*k + j*dimX+i;
-            grad_magn = (sqrtf(pow(P1[index],2) + pow(P2[index],2) + pow(P3[index],2)))/alpha1;
-            if (grad_magn > 1.0f) {
-                P1[index] /= grad_magn;
-                P2[index] /= grad_magn;
-                P3[index] /= grad_magn;
-            }
-        }}}
+            for(k=0; k<dimZ; k++) {
+                index = (dimX*dimY)*k + j*dimX+i;
+                grad_magn = (sqrtf(pow(P1[index],2) + pow(P2[index],2) + pow(P3[index],2)))/alpha1;
+                if (grad_magn > 1.0f) {
+                    P1[index] /= grad_magn;
+                    P2[index] /= grad_magn;
+                    P3[index] /= grad_magn;
+                }
+            }}}
     return 1;
 }
 /*Calculating dual variable Q (using forward differences)*/
@@ -338,33 +356,33 @@ float DualQ_3D(float *V1, float *V2, float *V3, float *Q1, float *Q2, float *Q3,
 #pragma omp parallel for shared(Q1,Q2,Q3,Q4,Q5,Q6,V1,V2,V3) private(i,j,k,index,q1,q2,q3,q11,q22,q33,q44,q55,q66)
     for(i=0; i<dimX; i++) {
         for(j=0; j<dimY; j++) {
-       	  for(k=0; k<dimZ; k++) {   	
-	    index = (dimX*dimY)*k + j*dimX+i;
-	    q1 = 0.0f; q11 = 0.0f; q33 = 0.0f; q2 = 0.0f; q22 = 0.0f; q55 = 0.0f; q3 = 0.0f; q44 = 0.0f; q66 = 0.0f;
-            /* symmetric boundary conditions (Neuman) */
-            if (i != dimX-1){ 
-                q1 = V1[(dimX*dimY)*k + j*dimX+(i+1)] - V1[index];              
-                q11 = V2[(dimX*dimY)*k + j*dimX+(i+1)] - V2[index];
-                q33 = V3[(dimX*dimY)*k + j*dimX+(i+1)] - V3[index];
-            }
-            if (j != dimY-1) {
-                q2 = V2[(dimX*dimY)*k + (j+1)*dimX+i] - V2[index];                
-                q22 = V1[(dimX*dimY)*k + (j+1)*dimX+i] - V1[index];
-                q55 = V3[(dimX*dimY)*k + (j+1)*dimX+i] - V3[index];
-            }
-            if (k != dimZ-1) {
-                q3 = V3[(dimX*dimY)*(k+1) + j*dimX+i] - V3[index];
-                q44 = V1[(dimX*dimY)*(k+1) + j*dimX+i] - V1[index];
-                q66 = V2[(dimX*dimY)*(k+1) + j*dimX+i] - V2[index];
-            }
-            
-            Q1[index] += sigma*(q1); /*Q11*/
-            Q2[index] += sigma*(q2); /*Q22*/            
-            Q3[index] += sigma*(q3); /*Q33*/
-            Q4[index] += sigma*(0.5f*(q11 + q22)); /* Q21 / Q12 */
-            Q5[index] += sigma*(0.5f*(q33 + q44)); /* Q31 / Q13 */
-            Q6[index] += sigma*(0.5f*(q55 + q66)); /* Q32 / Q23 */
-        }}}
+            for(k=0; k<dimZ; k++) {
+                index = (dimX*dimY)*k + j*dimX+i;
+                q1 = 0.0f; q11 = 0.0f; q33 = 0.0f; q2 = 0.0f; q22 = 0.0f; q55 = 0.0f; q3 = 0.0f; q44 = 0.0f; q66 = 0.0f;
+                /* symmetric boundary conditions (Neuman) */
+                if (i != dimX-1){
+                    q1 = V1[(dimX*dimY)*k + j*dimX+(i+1)] - V1[index];
+                    q11 = V2[(dimX*dimY)*k + j*dimX+(i+1)] - V2[index];
+                    q33 = V3[(dimX*dimY)*k + j*dimX+(i+1)] - V3[index];
+                }
+                if (j != dimY-1) {
+                    q2 = V2[(dimX*dimY)*k + (j+1)*dimX+i] - V2[index];
+                    q22 = V1[(dimX*dimY)*k + (j+1)*dimX+i] - V1[index];
+                    q55 = V3[(dimX*dimY)*k + (j+1)*dimX+i] - V3[index];
+                }
+                if (k != dimZ-1) {
+                    q3 = V3[(dimX*dimY)*(k+1) + j*dimX+i] - V3[index];
+                    q44 = V1[(dimX*dimY)*(k+1) + j*dimX+i] - V1[index];
+                    q66 = V2[(dimX*dimY)*(k+1) + j*dimX+i] - V2[index];
+                }
+                
+                Q1[index] += sigma*(q1); /*Q11*/
+                Q2[index] += sigma*(q2); /*Q22*/
+                Q3[index] += sigma*(q3); /*Q33*/
+                Q4[index] += sigma*(0.5f*(q11 + q22)); /* Q21 / Q12 */
+                Q5[index] += sigma*(0.5f*(q33 + q44)); /* Q31 / Q13 */
+                Q6[index] += sigma*(0.5f*(q55 + q66)); /* Q32 / Q23 */
+            }}}
     return 1;
 }
 float ProjQ_3D(float *Q1, float *Q2, float *Q3, float *Q4, float *Q5, float *Q6, long dimX, long dimY, long dimZ, float alpha0)
@@ -374,19 +392,19 @@ float ProjQ_3D(float *Q1, float *Q2, float *Q3, float *Q4, float *Q5, float *Q6,
 #pragma omp parallel for shared(Q1,Q2,Q3,Q4,Q5,Q6) private(i,j,k,index,grad_magn)
     for(i=0; i<dimX; i++) {
         for(j=0; j<dimY; j++) {
-       	  for(k=0; k<dimZ; k++) {   	
-	    index = (dimX*dimY)*k + j*dimX+i;           
-            grad_magn = sqrtf(pow(Q1[index],2) + pow(Q2[index],2) + pow(Q3[index],2) + 2.0f*pow(Q4[index],2) + 2.0f*pow(Q5[index],2) + 2.0f*pow(Q6[index],2));
-            grad_magn = grad_magn/alpha0;
-            if (grad_magn > 1.0f) {
-                Q1[index] /= grad_magn;
-                Q2[index] /= grad_magn;
-                Q3[index] /= grad_magn;
-                Q4[index] /= grad_magn;
-                Q5[index] /= grad_magn;
-                Q6[index] /= grad_magn;
-            }
-        }}}
+            for(k=0; k<dimZ; k++) {
+                index = (dimX*dimY)*k + j*dimX+i;
+                grad_magn = sqrtf(pow(Q1[index],2) + pow(Q2[index],2) + pow(Q3[index],2) + 2.0f*pow(Q4[index],2) + 2.0f*pow(Q5[index],2) + 2.0f*pow(Q6[index],2));
+                grad_magn = grad_magn/alpha0;
+                if (grad_magn > 1.0f) {
+                    Q1[index] /= grad_magn;
+                    Q2[index] /= grad_magn;
+                    Q3[index] /= grad_magn;
+                    Q4[index] /= grad_magn;
+                    Q5[index] /= grad_magn;
+                    Q6[index] /= grad_magn;
+                }
+            }}}
     return 1;
 }
 /* Divergence and projection for P*/
@@ -397,18 +415,22 @@ float DivProjP_3D(float *U, float *U0, float *P1, float *P2, float *P3, long dim
 #pragma omp parallel for shared(U,U0,P1,P2,P3) private(i,j,k,index,P_v1,P_v2,P_v3,div)
     for(i=0; i<dimX; i++) {
         for(j=0; j<dimY; j++) {
-       	  for(k=0; k<dimZ; k++) {   	
-	    index = (dimX*dimY)*k + j*dimX+i; 	    
-            if (i == 0) P_v1 = P1[index];
-            else P_v1 = P1[index] - P1[(dimX*dimY)*k + j*dimX+(i-1)];
-            if (j == 0) P_v2 = P2[index];
-            else P_v2 = P2[index] - P2[(dimX*dimY)*k + (j-1)*dimX+i];
-            if (k == 0) P_v3 = P3[index];
-            else P_v3 = P3[index] - P3[(dimX*dimY)*(k-1) + (j)*dimX+i];              
-                      
-            div = P_v1 + P_v2 + P_v3;
-            U[index] = (lambda*(U[index] + tau*div) + tau*U0[index])/(lambda + tau); 
-        }}}
+            for(k=0; k<dimZ; k++) {
+                index = (dimX*dimY)*k + j*dimX+i;
+                
+                if (i == 0) P_v1 = P1[index];
+                else if (i == dimX-1)  P_v1 = -P1[(dimX*dimY)*k + j*dimX+(i-1)];
+                else P_v1 = P1[index] - P1[(dimX*dimY)*k + j*dimX+(i-1)];
+                if (j == 0) P_v2 = P2[index];
+                else if (j == dimY-1) P_v2 = -P2[(dimX*dimY)*k + (j-1)*dimX+i];
+                else P_v2 = P2[index] - P2[(dimX*dimY)*k + (j-1)*dimX+i];
+                if (k == 0) P_v3 = P3[index];
+                else if (k == dimZ-1) P_v3 = -P3[(dimX*dimY)*(k-1) + (j)*dimX+i];
+                else P_v3 = P3[index] - P3[(dimX*dimY)*(k-1) + (j)*dimX+i];
+                
+                div = P_v1 + P_v2 + P_v3;
+                U[index] = (lambda*(U[index] + tau*div) + tau*U0[index])/(lambda + tau);
+            }}}
     return *U;
 }
 /*get update for V*/
@@ -419,47 +441,70 @@ float UpdV_3D(float *V1, float *V2, float *V3, float *P1, float *P2, float *P3,
 #pragma omp parallel for shared(V1,V2,V3,P1,P2,P3,Q1,Q2,Q3,Q4,Q5,Q6) private(i,j,k,index,q1,q4x,q5x,q2,q4y,q6y,q6z,q5z,q3,div1,div2,div3)
     for(i=0; i<dimX; i++) {
         for(j=0; j<dimY; j++) {
-       	  for(k=0; k<dimZ; k++) {   	
-	    index = (dimX*dimY)*k + j*dimX+i; 	
-	    q1 = 0.0f; q4x= 0.0f; q5x= 0.0f; q2= 0.0f; q4y= 0.0f; q6y= 0.0f; q6z= 0.0f; q5z= 0.0f; q3= 0.0f;
-            /* Q1 - Q11, Q2 - Q22, Q3 -  Q33, Q4 - Q21/Q12, Q5 - Q31/Q13, Q6 - Q32/Q23*/            
-            /* symmetric boundary conditions (Neuman) */
-            if (i != 0) {
-                q1 = Q1[index] - Q1[(dimX*dimY)*k + j*dimX+(i-1)];
-                q4x = Q4[index] - Q4[(dimX*dimY)*k + j*dimX+(i-1)];                
-                q5x = Q5[index] - Q5[(dimX*dimY)*k + j*dimX+(i-1)];
-            }
-            if (j != 0) {
-                q2 = Q2[index] - Q2[(dimX*dimY)*k + (j-1)*dimX+i];
-                q4y = Q4[index] - Q4[(dimX*dimY)*k + (j-1)*dimX+i];
-                q6y = Q6[index] - Q6[(dimX*dimY)*k + (j-1)*dimX+i];
-            }
-             if (k != 0) {
-                q6z = Q6[index] - Q6[(dimX*dimY)*(k-1) + (j)*dimX+i];
-                q5z = Q5[index] - Q5[(dimX*dimY)*(k-1) + (j)*dimX+i];
-                q3 = Q3[index] - Q3[(dimX*dimY)*(k-1) + (j)*dimX+i];
-            }
-            div1 = q1 + q4y + q5z;
-            div2 = q4x + q2 + q6z;            
-            div3 = q5x + q6y + q3;
-            
-            V1[index] += tau*(P1[index] + div1);
-            V2[index] += tau*(P2[index] + div2);
-            V3[index] += tau*(P3[index] + div3);
-        }}}
+            for(k=0; k<dimZ; k++) {
+                index = (dimX*dimY)*k + j*dimX+i;
+                q1 = 0.0f; q4x= 0.0f; q5x= 0.0f; q2= 0.0f; q4y= 0.0f; q6y= 0.0f; q6z= 0.0f; q5z= 0.0f; q3= 0.0f;
+                /* Q1 - Q11, Q2 - Q22, Q3 -  Q33, Q4 - Q21/Q12, Q5 - Q31/Q13, Q6 - Q32/Q23*/
+                /* symmetric boundary conditions (Neuman) */
+                
+                if (i == 0) {
+                    q1 = Q1[index];
+                    q4x = Q4[index];
+                    q5x = Q5[index]; }
+                else if (i == dimX-1) {
+                    q1 = -Q1[(dimX*dimY)*k + j*dimX+(i-1)];
+                    q4x = -Q4[(dimX*dimY)*k + j*dimX+(i-1)];
+                    q5x = -Q5[(dimX*dimY)*k + j*dimX+(i-1)]; }
+                else {
+                    q1 = Q1[index] - Q1[(dimX*dimY)*k + j*dimX+(i-1)];
+                    q4x = Q4[index] - Q4[(dimX*dimY)*k + j*dimX+(i-1)];
+                    q5x = Q5[index] - Q5[(dimX*dimY)*k + j*dimX+(i-1)]; }
+                if (j == 0) {
+                    q2 = Q2[index];
+                    q4y = Q4[index];
+                    q6y = Q6[index]; }
+                else if (j == dimY-1) {
+                    q2 = -Q2[(dimX*dimY)*k + (j-1)*dimX+i];
+                    q4y = -Q4[(dimX*dimY)*k + (j-1)*dimX+i];
+                    q6y = -Q6[(dimX*dimY)*k + (j-1)*dimX+i]; }
+                else {
+                    q2 = Q2[index] - Q2[(dimX*dimY)*k + (j-1)*dimX+i];
+                    q4y = Q4[index] - Q4[(dimX*dimY)*k + (j-1)*dimX+i];
+                    q6y = Q6[index] - Q6[(dimX*dimY)*k + (j-1)*dimX+i]; }
+                if (k == 0) {
+                    q6z = Q6[index];
+                    q5z = Q5[index];
+                    q3 = Q3[index]; }
+                else if (k == dimZ-1) {
+                    q6z = -Q6[(dimX*dimY)*(k-1) + (j)*dimX+i];
+                    q5z =  -Q5[(dimX*dimY)*(k-1) + (j)*dimX+i];
+                    q3 =  -Q3[(dimX*dimY)*(k-1) + (j)*dimX+i]; }
+                else {
+                    q6z = Q6[index] - Q6[(dimX*dimY)*(k-1) + (j)*dimX+i];
+                    q5z = Q5[index] - Q5[(dimX*dimY)*(k-1) + (j)*dimX+i];
+                    q3 = Q3[index] - Q3[(dimX*dimY)*(k-1) + (j)*dimX+i]; }
+                
+                div1 = q1 + q4y + q5z;
+                div2 = q4x + q2 + q6z;
+                div3 = q5x + q6y + q3;
+                
+                V1[index] += tau*(P1[index] + div1);
+                V2[index] += tau*(P2[index] + div2);
+                V3[index] += tau*(P3[index] + div3);
+            }}}
     return 1;
 }
 
 float copyIm_3Ar(float *V1, float *V2, float *V3, float *V1_old, float *V2_old, float *V3_old, long dimX, long dimY, long dimZ)
 {
-	long j;
+    long j;
 #pragma omp parallel for shared(V1, V2, V3, V1_old, V2_old, V3_old) private(j)
-	for (j = 0; j<dimX*dimY*dimZ; j++)  {	
-	V1_old[j] = V1[j];
-	V2_old[j] = V2[j];
-	V3_old[j] = V3[j];	
-	}
-	return 1;
+    for (j = 0; j<dimX*dimY*dimZ; j++)  {
+        V1_old[j] = V1[j];
+        V2_old[j] = V2[j];
+        V3_old[j] = V3[j];
+    }
+    return 1;
 }
 
 /*get updated solution U*/
@@ -478,9 +523,9 @@ float newU3D_3Ar(float *V1, float *V2, float *V3, float *V1_old, float *V2_old,
     long i;
 #pragma omp parallel for shared(V1, V2, V3, V1_old, V2_old, V3_old) private(i)
     for(i=0; i<dimX*dimY*dimZ; i++) {
-    V1[i] = 2.0f*V1[i] - V1_old[i];
-    V2[i] = 2.0f*V2[i] - V2_old[i];
-    V3[i] = 2.0f*V3[i] - V3_old[i];
+        V1[i] = 2.0f*V1[i] - V1_old[i];
+        V2[i] = 2.0f*V2[i] - V2_old[i];
+        V3[i] = 2.0f*V3[i] - V3_old[i];
     }
     return 1;
 }
diff --git a/src/Core/regularisers_GPU/TGV_GPU_core.cu b/src/Core/regularisers_GPU/TGV_GPU_core.cu
index 58b2c41..e4abf72 100644
--- a/src/Core/regularisers_GPU/TGV_GPU_core.cu
+++ b/src/Core/regularisers_GPU/TGV_GPU_core.cu
@@ -3,8 +3,8 @@ 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
+Copyright 2019 Daniil Kazantsev
+Copyright 2019 Srikanth Nagella, Edoardo Pasca
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
@@ -32,7 +32,7 @@ limitations under the License.
  * 6. Lipshitz constant (default is 12)
  *
  * Output:
- * Filtered/regulariaed image 
+ * Filtered/regularised image 
  *
  * References:
  * [1] K. Bredies "Total Generalized Variation"
@@ -42,8 +42,8 @@ limitations under the License.
 #define BLKYSIZE 8
 #define BLKZSIZE 8    
     
-#define BLKXSIZE2D 16
-#define BLKYSIZE2D 16
+#define BLKXSIZE2D 8
+#define BLKYSIZE2D 8
 #define EPS 1.0e-7
 #define idivup(a, b) ( ((a)%(b) != 0) ? (a)/(b)+1 : (a)/(b) )
 
@@ -53,80 +53,84 @@ limitations under the License.
 /********************************************************************/
 __global__ void DualP_2D_kernel(float *U, float *V1, float *V2, float *P1, float *P2, int dimX, int dimY, float sigma)
 {    
-	int i = blockDim.x * blockIdx.x + threadIdx.x;
-        int j = blockDim.y * blockIdx.y + threadIdx.y;
+	int num_total = dimX*dimY;
+	const int i = blockDim.x * blockIdx.x + threadIdx.x;
+        const int j = blockDim.y * blockIdx.y + threadIdx.y;
         
-        int index = i + dimX*j;
+        int index = i + dimX*j;       
         
-        if ((i >= 0) && (i < dimX) && (j >= 0) && (j < dimY)) {
-            /* symmetric boundary conditions (Neuman) */
-            if (i == dimX-1) P1[index] += sigma*((U[j*dimX+(i-1)] - U[index]) - V1[index]); 
-            else P1[index] += sigma*((U[j*dimX+(i+1)] - U[index])  - V1[index]); 
-            if (j == dimY-1) P2[index] += sigma*((U[(j-1)*dimX+i] - U[index])  - V2[index]);
-            else  P2[index] += sigma*((U[(j+1)*dimX+i] - U[index])  - V2[index]);
-		}
+        if (index < num_total) {  
+        /* symmetric boundary conditions (Neuman) */            
+        if ((i >= 0) && (i < dimX-1))  P1[index] += sigma*((U[(i+1) + dimX*j] - U[index])  - V1[index]); 
+        else P1[index] += sigma*(-V1[index]); 
+        if ((j >= 0) && (j < dimY-1))  P2[index] += sigma*((U[i + dimX*(j+1)] - U[index])  - V2[index]);
+        else P2[index] += sigma*(-V2[index]);                    
+	}
 	return;
 } 
 
 __global__ void ProjP_2D_kernel(float *P1, float *P2, int dimX, int dimY, float alpha1)
 {
    	float grad_magn;
-
-	int i = blockDim.x * blockIdx.x + threadIdx.x;
-        int j = blockDim.y * blockIdx.y + threadIdx.y;
+	int num_total = dimX*dimY;
+	
+	const int i = blockDim.x * blockIdx.x + threadIdx.x;
+        const int j = blockDim.y * blockIdx.y + threadIdx.y;
         
         int index = i + dimX*j;
-        
-        if ((i >= 0) && (i < dimX) && (j >= 0) && (j < dimY)) {
-            
-            grad_magn = sqrt(pow(P1[index],2) + pow(P2[index],2));
+                
+        if (index < num_total) {            
+            grad_magn = sqrtf(pow(P1[index],2) + pow(P2[index],2));
             grad_magn = grad_magn/alpha1;
             if (grad_magn > 1.0f) {
                 P1[index] /= grad_magn;
                 P2[index] /= grad_magn;
             }
-		}
+	}
 	return;
 } 
 
 __global__ void DualQ_2D_kernel(float *V1, float *V2, float *Q1, float *Q2, float *Q3, int dimX, int dimY, float sigma)
 {
         float q1, q2, q11, q22;
-
-	int i = blockDim.x * blockIdx.x + threadIdx.x;
-        int j = blockDim.y * blockIdx.y + threadIdx.y;
+	int num_total = dimX*dimY;
+	
+	const int i = blockDim.x * blockIdx.x + threadIdx.x;
+        const int j = blockDim.y * blockIdx.y + threadIdx.y;
         
-        int index = i + dimX*j;
+        int index = i + dimX*j;      
         
-        if ((i >= 0) && (i < dimX) && (j >= 0) && (j < dimY)) {            
-            /* symmetric boundary conditions (Neuman) */
-    	    q1 = 0.0f; q11 = 0.0f; q2 = 0.0f; q22 = 0.0f;
-            /* boundary conditions (Neuman) */
-            if (i != dimX-1){
-                q1 = V1[j*dimX+(i+1)] - V1[index];
-                q11 = V2[j*dimX+(i+1)] - V2[index];
-            }
-            if (j != dimY-1) {
-                q2 = V2[(j+1)*dimX+i] - V2[index];
-                q22 = V1[(j+1)*dimX+i] - V1[index];
-            }
+        if (index < num_total) {
+         q1 = 0.0f; q2  = 0.0f; q11  = 0.0f; q22  = 0.0f;
+         
+	        if ((i >= 0) && (i < dimX-1))  {            
+        	    /* boundary conditions (Neuman) */
+        	    q1 = V1[(i+1) + dimX*j] - V1[index];
+        	    q11 = V2[(i+1) + dimX*j] - V2[index];
+	        }
+        	if ((j >= 0) && (j < dimY-1)) {
+        	    q2 = V2[i + dimX*(j+1)] - V2[index];
+        	    q22 = V1[i + dimX*(j+1)] - V1[index];
+        	}
+        	
             Q1[index] += sigma*(q1);
             Q2[index] += sigma*(q2);
             Q3[index] += sigma*(0.5f*(q11 + q22));
-	}
+	 }            
 	return;
 } 
 
 __global__ void ProjQ_2D_kernel(float *Q1, float *Q2, float *Q3, int dimX, int dimY, float alpha0)
 {
 	float grad_magn;
+        int num_total = dimX*dimY;
 
-	int i = blockDim.x * blockIdx.x + threadIdx.x;
-        int j = blockDim.y * blockIdx.y + threadIdx.y;
-        
-        int index = i + dimX*j;
+	const int i = blockDim.x * blockIdx.x + threadIdx.x;
+        const int j = blockDim.y * blockIdx.y + threadIdx.y;
         
-        if ((i >= 0) && (i < dimX) && (j >= 0) && (j < dimY)) {            
+        int index = i + dimX*j;        
+       
+        if (index < num_total) {
             grad_magn = sqrt(pow(Q1[index],2) + pow(Q2[index],2) + 2*pow(Q3[index],2));
             grad_magn = grad_magn/alpha0;
             if (grad_magn > 1.0f) {
@@ -141,44 +145,75 @@ __global__ void ProjQ_2D_kernel(float *Q1, float *Q2, float *Q3, int dimX, int d
 __global__ void DivProjP_2D_kernel(float *U, float *U0, float *P1, float *P2, int dimX, int dimY, float lambda, float tau)
 {
 	float P_v1, P_v2, div;
-
-	int i = blockDim.x * blockIdx.x + threadIdx.x;
-        int j = blockDim.y * blockIdx.y + threadIdx.y;
+	int num_total = dimX*dimY;
+	
+	const int i = blockDim.x * blockIdx.x + threadIdx.x;
+        const int j = blockDim.y * blockIdx.y + threadIdx.y;       
         
         int index = i + dimX*j;
+
+        if (index < num_total) { 
+	P_v1 = 0.0f; P_v2 = 0.0f;
         
-        if ((i >= 0) && (i < dimX) && (j >= 0) && (j < dimY)) {
-			
-            if (i == 0) P_v1 = P1[index];
-            else P_v1 = P1[index] - P1[j*dimX+(i-1)];
-            if (j == 0) P_v2 = P2[index];
-            else  P_v2 = P2[index] - P2[(j-1)*dimX+i];
-            div = P_v1 + P_v2;
-            U[index] = (lambda*(U[index] + tau*div) + tau*U0[index])/(lambda + tau);
-		}
+        if (i == 0) P_v1 = P1[index];
+        if (i == dimX-1) P_v1 = -P1[(i-1) + dimX*j];
+        if ((i > 0) && (i < dimX-1)) P_v1 = P1[index] - P1[(i-1) + dimX*j];
+
+        if (j == 0) P_v2 = P2[index];
+        if (j == dimY-1) P_v2 = -P2[i + dimX*(j-1)];
+      	if ((j > 0) && (j < dimY-1))  P_v2 = P2[index] - P2[i + dimX*(j-1)];
+         
+        div = P_v1 + P_v2;
+        U[index] = (lambda*(U[index] + tau*div) + tau*U0[index])/(lambda + tau);
+	}
 	return;
 } 
 
 __global__ void UpdV_2D_kernel(float *V1, float *V2, float *P1, float *P2, float *Q1, float *Q2, float *Q3, int dimX, int dimY, float tau)
 {
 	float q1, q3_x, q2, q3_y, div1, div2;
-
-	int i = blockDim.x * blockIdx.x + threadIdx.x;
-        int j = blockDim.y * blockIdx.y + threadIdx.y;
-        
-        int index = i + dimX*j;
+	int num_total = dimX*dimY;
+	int i1, j1;
+	
+	const int i = blockDim.x * blockIdx.x + threadIdx.x;
+        const int j = blockDim.y * blockIdx.y + threadIdx.y;
         
-        if ((i >= 0) && (i < dimX) && (j >= 0) && (j < dimY)) {			
-   	    q2 = 0.0f;  q3_y = 0.0f; q1 = 0.0f; q3_x = 0.0;
-            /* boundary conditions (Neuman) */
-            if (i != 0) {
-                q1 = Q1[index] - Q1[j*dimX+(i-1)];
-                q3_x = Q3[index] - Q3[j*dimX+(i-1)];
-            }
-            if (j != 0) {
-                q2 = Q2[index] - Q2[(j-1)*dimX+i];
-                q3_y = Q3[index] - Q3[(j-1)*dimX+i];
-            }
+        int index = i + dimX*j;      
+      
+	if (index < num_total) {   
+	
+	    i1 = (i-1) + dimX*j;
+            j1 = (i) + dimX*(j-1);
+
+            /* boundary conditions (Neuman) */        
+            if ((i > 0) && (i < dimX-1)) {
+            q1 = Q1[index] - Q1[i1];
+            q3_x = Q3[index] - Q3[i1];  }            
+            else if (i == 0) {
+            q1 = Q1[index];
+            q3_x = Q3[index]; } 
+            else if (i == dimX-1) {
+            q1 = -Q1[i1];
+            q3_x = -Q3[i1];  }
+            else {
+            q1 = 0.0f;
+            q3_x = 0.0f;
+            }    
+            
+            if ((j > 0) && (j < dimY-1)) {
+            q2 = Q2[index] - Q2[j1];
+            q3_y = Q3[index] - Q3[j1]; } 
+            else if (j == dimY-1) {
+            q2 = -Q2[j1];
+            q3_y = -Q3[j1]; }
+            else if (j == 0) {
+            q2 = Q2[index];
+            q3_y = Q3[index]; }
+            else {
+            q2 = 0.0f;
+            q3_y = 0.0f;
+            }       
+            
             div1 = q1 + q3_y;
             div2 = q3_x + q2;
             V1[index] += tau*(P1[index] + div1);
@@ -243,21 +278,22 @@ __global__ void newU_kernel_ar2(float *V1, float *V2, float *V1_old, float *V2_o
 __global__ void DualP_3D_kernel(float *U, float *V1, float *V2, float *V3, float *P1, float *P2, float *P3, int dimX, int dimY, int dimZ, float sigma)
 {    
 	int index;
-	int i = blockDim.x * blockIdx.x + threadIdx.x;
-        int j = blockDim.y * blockIdx.y + threadIdx.y;
-        int k = blockDim.z * blockIdx.z + threadIdx.z;
+	const int i = blockDim.x * blockIdx.x + threadIdx.x;
+        const int j = blockDim.y * blockIdx.y + threadIdx.y;
+        const int k = blockDim.z * blockIdx.z + threadIdx.z;
+
+   	int num_total = dimX*dimY*dimZ;
         
-        if ((i >= 0) && (i < dimX) && (j >= 0) && (j < dimY) && (k >= 0) && (k < dimZ)) {
-	
-	    index = (dimX*dimY)*k + j*dimX+i;
-            /* symmetric boundary conditions (Neuman) */
-            if (i == dimX-1) P1[index] += sigma*((U[(dimX*dimY)*k + j*dimX+(i-1)] - U[index]) - V1[index]); 
-            else P1[index] += sigma*((U[(dimX*dimY)*k + j*dimX+(i+1)] - U[index])  - V1[index]); 
-            if (j == dimY-1) P2[index] += sigma*((U[(dimX*dimY)*k + (j-1)*dimX+i] - U[index])  - V2[index]);
-            else  P2[index] += sigma*((U[(dimX*dimY)*k + (j+1)*dimX+i] - U[index])  - V2[index]);
-            if (k == dimZ-1) P3[index] += sigma*((U[(dimX*dimY)*(k-1) + j*dimX+i] - U[index])  - V3[index]);
-            else  P3[index] += sigma*((U[(dimX*dimY)*(k+1) + j*dimX+i] - U[index])  - V3[index]);
-	}
+        index = (dimX*dimY)*k + i*dimX+j;    
+        if (index < num_total) {                
+            /* symmetric boundary conditions (Neuman) */            
+            if ((i >= 0) && (i < dimX-1)) P1[index] += sigma*((U[(dimX*dimY)*k + (i+1)*dimX+j] - U[index])  - V1[index]);  
+	    else P1[index] += sigma*(-V1[index]); 
+	    if ((j >= 0) && (j < dimY-1)) P2[index] += sigma*((U[(dimX*dimY)*k + i*dimX+(j+1)] - U[index])  - V2[index]);        
+	    else P2[index] += sigma*(-V2[index]);                
+      	    if ((k >= 0) && (k < dimZ-1)) P3[index] += sigma*((U[(dimX*dimY)*(k+1) + i*dimX+(j)] - U[index])  - V3[index]);        
+	    else P3[index] += sigma*(-V3[index]);                	    
+	 }	
 	return;
 } 
 
@@ -265,14 +301,14 @@ __global__ void ProjP_3D_kernel(float *P1, float *P2, float *P3, int dimX, int d
 {
    	float grad_magn;
    	int index;
+   	int num_total = dimX*dimY*dimZ;
    	
 	int i = blockDim.x * blockIdx.x + threadIdx.x;
         int j = blockDim.y * blockIdx.y + threadIdx.y;
         int k = blockDim.z * blockIdx.z + threadIdx.z;
         
-        if ((i >= 0) && (i < dimX) && (j >= 0) && (j < dimY) && (k >= 0) && (k < dimZ)) {	
-	    index = (dimX*dimY)*k + j*dimX+i;
-            
+        index = (dimX*dimY)*k + i*dimX+j;    
+        if (index < num_total) {            
             grad_magn = (sqrtf(pow(P1[index],2) + pow(P2[index],2) + pow(P3[index],2)))/alpha1;
             if (grad_magn > 1.0f) {
                 P1[index] /= grad_magn;
@@ -288,55 +324,57 @@ __global__ void DualQ_3D_kernel(float *V1, float *V2, float *V3, float *Q1, floa
 	int index; 
         float q1, q2, q3, q11, q22, q33, q44, q55, q66;
 
+   	int num_total = dimX*dimY*dimZ;
+   	
 	int i = blockDim.x * blockIdx.x + threadIdx.x;
         int j = blockDim.y * blockIdx.y + threadIdx.y;
         int k = blockDim.z * blockIdx.z + threadIdx.z;
         
-        if ((i >= 0) && (i < dimX) && (j >= 0) && (j < dimY) && (k >= 0) && (k < dimZ)) {	
-	    
-	    index = (dimX*dimY)*k + j*dimX+i;	    
-	    q1 = 0.0f; q11 = 0.0f; q33 = 0.0f; q2 = 0.0f; q22 = 0.0f; q55 = 0.0f; q3 = 0.0f; q44 = 0.0f; q66 = 0.0f;
-            /* symmetric boundary conditions (Neuman) */
-            if (i != dimX-1){ 
-                q1 = V1[(dimX*dimY)*k + j*dimX+(i+1)] - V1[index];              
-                q11 = V2[(dimX*dimY)*k + j*dimX+(i+1)] - V2[index];
-                q33 = V3[(dimX*dimY)*k + j*dimX+(i+1)] - V3[index];
-            }
-            if (j != dimY-1) {
-                q2 = V2[(dimX*dimY)*k + (j+1)*dimX+i] - V2[index];                
-                q22 = V1[(dimX*dimY)*k + (j+1)*dimX+i] - V1[index];
-                q55 = V3[(dimX*dimY)*k + (j+1)*dimX+i] - V3[index];
-            }
-            if (k != dimZ-1) {
-                q3 = V3[(dimX*dimY)*(k+1) + j*dimX+i] - V3[index];
-                q44 = V1[(dimX*dimY)*(k+1) + j*dimX+i] - V1[index];
-                q66 = V2[(dimX*dimY)*(k+1) + j*dimX+i] - V2[index];
-            }
-            
+        index = (dimX*dimY)*k + i*dimX+j; 
+        int i1 = (dimX*dimY)*k + (i+1)*dimX+j;
+        int j1 = (dimX*dimY)*k + (i)*dimX+(j+1);
+        int k1 = (dimX*dimY)*(k+1) + (i)*dimX+(j);
+                
+        if (index < num_total) {
+ 	q1 = 0.0f; q11 = 0.0f; q33 = 0.0f; q2 = 0.0f; q22 = 0.0f; q55 = 0.0f; q3 = 0.0f; q44 = 0.0f; q66 = 0.0f;
+         
+	        /* boundary conditions (Neuman) */
+	        if ((i >= 0) && (i < dimX-1))  {                    	
+                q1 = V1[i1] - V1[index];              
+                q11 = V2[i1] - V2[index];
+                q33 = V3[i1] - V3[index];  }
+        	if ((j >= 0) && (j < dimY-1)) {
+                q2 = V2[j1] - V2[index];                
+                q22 = V1[j1] - V1[index];
+                q55 = V3[j1] - V3[index];  }
+        	if ((k >= 0) && (k < dimZ-1)) {
+                q3 = V3[k1] - V3[index];
+                q44 = V1[k1] - V1[index];
+                q66 = V2[k1] - V2[index]; }
+        	
             Q1[index] += sigma*(q1); /*Q11*/
             Q2[index] += sigma*(q2); /*Q22*/            
             Q3[index] += sigma*(q3); /*Q33*/
             Q4[index] += sigma*(0.5f*(q11 + q22)); /* Q21 / Q12 */
             Q5[index] += sigma*(0.5f*(q33 + q44)); /* Q31 / Q13 */
             Q6[index] += sigma*(0.5f*(q55 + q66)); /* Q32 / Q23 */
-	}
+	 }
 	return;
 }
 
-
 __global__ void ProjQ_3D_kernel(float *Q1, float *Q2, float *Q3, float *Q4, float *Q5, float *Q6, int dimX, int dimY, int dimZ, float alpha0)
 {
 	float grad_magn;
 	int index;
-
+   	int num_total = dimX*dimY*dimZ;
+   	
 	int i = blockDim.x * blockIdx.x + threadIdx.x;
         int j = blockDim.y * blockIdx.y + threadIdx.y;
-        int k = blockDim.z * blockIdx.z + threadIdx.z;
+        int k = blockDim.z * blockIdx.z + threadIdx.z;        
+
+        index = (dimX*dimY)*k + i*dimX+j; 
         
-        if ((i >= 0) && (i < dimX) && (j >= 0) && (j < dimY) && (k >= 0) && (k < dimZ)) {	
-	    
-        index = (dimX*dimY)*k + j*dimX+i;	
-	
+        if (index < num_total) {	
 	grad_magn = sqrtf(pow(Q1[index],2) + pow(Q2[index],2) + pow(Q3[index],2) + 2.0f*pow(Q4[index],2) + 2.0f*pow(Q5[index],2) + 2.0f*pow(Q6[index],2));
             grad_magn = grad_magn/alpha0;
             if (grad_magn > 1.0f) {
@@ -354,21 +392,33 @@ __global__ void DivProjP_3D_kernel(float *U, float *U0, float *P1, float *P2, fl
 {
 	float P_v1, P_v2, P_v3, div;
 	int index;
-
+   	int num_total = dimX*dimY*dimZ;
+   	
+   	
 	int i = blockDim.x * blockIdx.x + threadIdx.x;
         int j = blockDim.y * blockIdx.y + threadIdx.y;
-        int k = blockDim.z * blockIdx.z + threadIdx.z;
+        int k = blockDim.z * blockIdx.z + threadIdx.z;       
+     
+        index = (dimX*dimY)*k + i*dimX+j; 
+        int i1 = (dimX*dimY)*k + (i-1)*dimX+j;
+        int j1 = (dimX*dimY)*k + (i)*dimX+(j-1);
+        int k1 = (dimX*dimY)*(k-1) + (i)*dimX+(j);        
+        
+        if (index < num_total) { 
+	P_v1 = 0.0f; P_v2 = 0.0f; P_v3 = 0.0f;
         
-        if ((i >= 0) && (i < dimX) && (j >= 0) && (j < dimY) && (k >= 0) && (k < dimZ)) {	
-
-        index = (dimX*dimY)*k + j*dimX+i;	
-			
         if (i == 0) P_v1 = P1[index];
-        else P_v1 = P1[index] - P1[(dimX*dimY)*k + j*dimX+(i-1)];
+        if (i == dimX-1) P_v1 = -P1[i1];
+        if ((i > 0) && (i < dimX-1)) P_v1 = P1[index] - P1[i1];
+
         if (j == 0) P_v2 = P2[index];
-        else P_v2 = P2[index] - P2[(dimX*dimY)*k + (j-1)*dimX+i];
+        if (j == dimY-1) P_v2 = -P2[j1];
+      	if ((j > 0) && (j < dimY-1))  P_v2 = P2[index] - P2[j1];
+        
         if (k == 0) P_v3 = P3[index];
-        else P_v3 = P3[index] - P3[(dimX*dimY)*(k-1) + (j)*dimX+i];              
+        if (k == dimZ-1) P_v3 = -P3[k1];
+      	if ((k > 0) && (k < dimZ-1))  P_v3 = P3[index] - P3[k1];
+     
                       
         div = P_v1 + P_v2 + P_v3;
         U[index] = (lambda*(U[index] + tau*div) + tau*U0[index])/(lambda + tau);             
@@ -379,33 +429,73 @@ __global__ void UpdV_3D_kernel(float *V1, float *V2, float *V3, float *P1, float
 {
 	float q1, q4x, q5x, q2, q4y, q6y, q6z, q5z, q3, div1, div2, div3;
 	int index;
+	int num_total = dimX*dimY*dimZ;
 	
 	int i = blockDim.x * blockIdx.x + threadIdx.x;
         int j = blockDim.y * blockIdx.y + threadIdx.y;
         int k = blockDim.z * blockIdx.z + threadIdx.z;
         
-        if ((i >= 0) && (i < dimX) && (j >= 0) && (j < dimY) && (k >= 0) && (k < dimZ)) {	
-
-        index = (dimX*dimY)*k + j*dimX+i;	
+        index = (dimX*dimY)*k + i*dimX+j; 
+        int i1 = (dimX*dimY)*k + (i-1)*dimX+j;
+        int j1 = (dimX*dimY)*k + (i)*dimX+(j-1);
+        int k1 = (dimX*dimY)*(k-1) + (i)*dimX+(j);    
         
-	q1 = 0.0f; q4x= 0.0f; q5x= 0.0f; q2= 0.0f; q4y= 0.0f; q6y= 0.0f; q6z= 0.0f; q5z= 0.0f; q3= 0.0f;
-        /* Q1 - Q11, Q2 - Q22, Q3 -  Q33, Q4 - Q21/Q12, Q5 - Q31/Q13, Q6 - Q32/Q23*/            
-        /* symmetric boundary conditions (Neuman) */
-        if (i != 0) {
-                q1 = Q1[index] - Q1[(dimX*dimY)*k + j*dimX+(i-1)];
-                q4x = Q4[index] - Q4[(dimX*dimY)*k + j*dimX+(i-1)];                
-                q5x = Q5[index] - Q5[(dimX*dimY)*k + j*dimX+(i-1)];
-        }
-       if (j != 0) {
-                q2 = Q2[index] - Q2[(dimX*dimY)*k + (j-1)*dimX+i];
-                q4y = Q4[index] - Q4[(dimX*dimY)*k + (j-1)*dimX+i];
-                q6y = Q6[index] - Q6[(dimX*dimY)*k + (j-1)*dimX+i];
-       }
-       if (k != 0) {
-                q6z = Q6[index] - Q6[(dimX*dimY)*(k-1) + (j)*dimX+i];
-                q5z = Q5[index] - Q5[(dimX*dimY)*(k-1) + (j)*dimX+i];
-                q3 = Q3[index] - Q3[(dimX*dimY)*(k-1) + (j)*dimX+i];
-       }
+        /* Q1 - Q11, Q2 - Q22, Q3 -  Q33, Q4 - Q21/Q12, Q5 - Q31/Q13, Q6 - Q32/Q23*/       
+        if (index < num_total) {  
+
+  	 /* boundary conditions (Neuman) */        
+            if ((i > 0) && (i < dimX-1)) {
+                q1 = Q1[index] - Q1[i1];
+                q4x = Q4[index] - Q4[i1];                
+                q5x = Q5[index] - Q5[i1]; }            
+            else if (i == 0) {
+                q1 = Q1[index];
+                q4x = Q4[index];                
+                q5x = Q5[index]; } 
+            else if (i == dimX-1) {
+                q1 = -Q1[i1];
+                q4x = -Q4[i1];                
+                q5x = -Q5[i1]; }
+            else {
+                q1 = 0.0f;
+                q4x = 0.0f;
+                q5x = 0.0f;  }    
+            
+            if ((j > 0) && (j < dimY-1)) {
+                q2 = Q2[index] - Q2[j1];
+                q4y = Q4[index] - Q4[j1];
+                q6y = Q6[index] - Q6[j1]; } 
+            else if (j == dimY-1) {
+                q2 = -Q2[j1];
+                q4y = -Q4[j1];
+                q6y = -Q6[j1]; }
+            else if (j == 0) {
+                q2 = Q2[index];
+                q4y = Q4[index];
+                q6y = Q6[index]; }
+            else {
+                q2 =  0.0f;
+                q4y = 0.0f;
+                q6y = 0.0f;
+               }       
+
+            if ((k > 0) && (k < dimZ-1)) {
+                q6z = Q6[index] - Q6[k1];
+                q5z = Q5[index] - Q5[k1];
+                q3 = Q3[index] - Q3[k1]; } 
+            else if (k == dimZ-1) {
+                q6z = -Q6[k1];
+                q5z = -Q5[k1];
+                q3 = -Q3[k1]; }
+            else if (k == 0) {
+                q6z = Q6[index];
+                q5z = Q5[index];
+                q3 = Q3[index]; }
+            else {
+                q6z = 0.0f;
+                q5z = 0.0f;
+                q3 = 0.0f; }
+
        div1 = q1 + q4y + q5z;
        div2 = q4x + q2 + q6z;            
        div3 = q5x + q6y + q3;
@@ -510,7 +600,14 @@ extern "C" int TGV_GPU_main(float *U0, float *U, float lambda, float alpha1, flo
         CHECK(cudaMalloc((void**)&V2_old,dimTotal*sizeof(float)));
         
         CHECK(cudaMemcpy(d_U0,U0,dimTotal*sizeof(float),cudaMemcpyHostToDevice));
-        CHECK(cudaMemcpy(d_U,U0,dimTotal*sizeof(float),cudaMemcpyHostToDevice));      
+        CHECK(cudaMemcpy(d_U,U0,dimTotal*sizeof(float),cudaMemcpyHostToDevice));   
+        cudaMemset(P1, 0, dimTotal*sizeof(float));
+        cudaMemset(P2, 0, dimTotal*sizeof(float));
+        cudaMemset(Q1, 0, dimTotal*sizeof(float));
+        cudaMemset(Q2, 0, dimTotal*sizeof(float));
+        cudaMemset(Q3, 0, dimTotal*sizeof(float));
+        cudaMemset(V1, 0, dimTotal*sizeof(float));
+        cudaMemset(V2, 0, dimTotal*sizeof(float));           
         
         if (dimZ == 1) {
 	/*2D case */
@@ -521,14 +618,14 @@ extern "C" int TGV_GPU_main(float *U0, float *U, float lambda, float alpha1, flo
 			
 	    /* Calculate Dual Variable P */
             DualP_2D_kernel<<<dimGrid,dimBlock>>>(d_U, V1, V2, P1, P2, dimX, dimY, sigma);
-	    CHECK(cudaDeviceSynchronize());
+      	    CHECK(cudaDeviceSynchronize());
             /*Projection onto convex set for P*/
             ProjP_2D_kernel<<<dimGrid,dimBlock>>>(P1, P2, dimX, dimY, alpha1);
             CHECK(cudaDeviceSynchronize());
             /* Calculate Dual Variable Q */
             DualQ_2D_kernel<<<dimGrid,dimBlock>>>(V1, V2, Q1, Q2, Q3, dimX, dimY, sigma);
             CHECK(cudaDeviceSynchronize());
-             /*Projection onto convex set for Q*/
+            /*Projection onto convex set for Q*/
             ProjQ_2D_kernel<<<dimGrid,dimBlock>>>(Q1, Q2, Q3, dimX, dimY, alpha0);
             CHECK(cudaDeviceSynchronize());
             /*saving U into U_old*/
@@ -549,7 +646,7 @@ extern "C" int TGV_GPU_main(float *U0, float *U, float lambda, float alpha1, flo
             /*get new V*/
             newU_kernel_ar2<<<dimGrid,dimBlock>>>(V1, V2, V1_old, V2_old, dimX, dimY, dimTotal);
             CHECK(cudaDeviceSynchronize());            
-	        }
+	    }
         }
         else {
         /*3D case */
@@ -565,6 +662,12 @@ extern "C" int TGV_GPU_main(float *U0, float *U, float lambda, float alpha1, flo
         CHECK(cudaMalloc((void**)&V3,dimTotal*sizeof(float)));
         CHECK(cudaMalloc((void**)&V3_old,dimTotal*sizeof(float)));
         
+        cudaMemset(Q4, 0.0f, dimTotal*sizeof(float));
+        cudaMemset(Q5, 0.0f, dimTotal*sizeof(float));
+        cudaMemset(Q6, 0.0f, dimTotal*sizeof(float));
+        cudaMemset(P3, 0.0f, dimTotal*sizeof(float));
+        cudaMemset(V3, 0.0f, dimTotal*sizeof(float));        
+        
         for(int n=0; n < iterationsNumb; n++) {
 			
 	    /* Calculate Dual Variable P */
diff --git a/src/Matlab/mex_compile/regularisers_GPU/TGV_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/TGV_GPU.cpp
index edb551d..1173282 100644
--- a/src/Matlab/mex_compile/regularisers_GPU/TGV_GPU.cpp
+++ b/src/Matlab/mex_compile/regularisers_GPU/TGV_GPU.cpp
@@ -21,18 +21,18 @@ limitations under the License.
 #include "TGV_GPU_core.h"
 
 /* CUDA implementation of Primal-Dual denoising method for 
- * Total Generilized Variation (TGV)-L2 model [1] (2D case only)
+ * Total Generilized Variation (TGV)-L2 model [1] (2D/3D)
  *
  * Input Parameters:
- * 1. Noisy image (2D) (required)
- * 2. lambda - regularisation parameter (required)
- * 3. parameter to control the first-order term (alpha1) (default - 1)
- * 4. parameter to control the second-order term (alpha0) (default - 0.5)
- * 5. Number of Chambolle-Pock (Primal-Dual) iterations (default is 300)
+ * 1. Noisy image/volume (2D/3D)
+ * 2. lambda - regularisation parameter
+ * 3. parameter to control the first-order term (alpha1)
+ * 4. parameter to control the second-order term (alpha0)
+ * 5. Number of Chambolle-Pock (Primal-Dual) iterations
  * 6. Lipshitz constant (default is 12)
  *
  * Output:
- * Filtered/regulariaed image 
+ * Filtered/regularised image 
  *
  * References:
  * [1] K. Bredies "Total Generalized Variation"
@@ -44,7 +44,7 @@ void mexFunction(
         
 {
     int number_of_dims, iter;
-    mwSize dimX, dimY;
+    mwSize dimX, dimY, dimZ;
     const mwSize *dim_array;
     float *Input, *Output=NULL, lambda, alpha0, alpha1, L2;
     
@@ -57,8 +57,8 @@ void mexFunction(
     Input  = (float *) mxGetData(prhs[0]); /*noisy image (2D) */
     lambda =  (float) mxGetScalar(prhs[1]); /* regularisation parameter */
     alpha1 =  1.0f; /* parameter to control the first-order term */ 
-    alpha0 =  0.5f; /* parameter to control the second-order term */
-    iter =  300; /* Iterations number */      
+    alpha0 =  2.0f; /* parameter to control the second-order term */
+    iter =  500; /* Iterations number */      
     L2 =  12.0f; /* Lipshitz constant */
     
     if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }   
@@ -68,12 +68,14 @@ void mexFunction(
     if (nrhs == 6)  L2 =  (float) mxGetScalar(prhs[5]); /* Lipshitz constant */
     
     /*Handling Matlab output data*/
-    dimX = dim_array[0]; dimY = dim_array[1];
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
     
     if (number_of_dims == 2) {
-        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
-        /* running the function */
-        TGV_GPU_main(Input, Output, lambda, alpha1, alpha0, iter, L2, dimX, dimY);        
+        dimZ = 1; /*2D case*/
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));        
     }
-    if (number_of_dims == 3) {mexErrMsgTxt("Only 2D images accepted");}       
+    if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+    
+    /* running the function */
+    TGV_GPU_main(Input, Output, lambda, alpha1, alpha0, iter, L2, dimX, dimY, dimZ);        
 }
diff --git a/src/Python/setup-regularisers.py.in b/src/Python/setup-regularisers.py.in
index 59be768..82d9f9f 100644
--- a/src/Python/setup-regularisers.py.in
+++ b/src/Python/setup-regularisers.py.in
@@ -68,7 +68,7 @@ setup(
     
     ],
 	zip_safe = False,	
-	packages = {'ccpi','ccpi.filters'},
+	packages = {'ccpi','ccpi.filters', 'ccpi.supp'},
 )
 
 
diff --git a/test/test_CPU_regularisers.py b/test/test_CPU_regularisers.py
index 42e4735..8940926 100644
--- a/test/test_CPU_regularisers.py
+++ b/test/test_CPU_regularisers.py
@@ -9,7 +9,7 @@ from testroutines import *
 
 class TestRegularisers(unittest.TestCase):
 
-    def getPars(self,alg,noi=1200):
+    def getPars(self):
         filename = os.path.join("lena_gray_512.tif")
         plt = TiffReader()
         # read image
@@ -28,64 +28,101 @@ class TestRegularisers(unittest.TestCase):
         u0 = u0.astype('float32')
         u_ref = u_ref.astype('float32')
         # set parameters
-        pars = {'algorithm': alg, \
-                'input': u0, \
-                'regularisation_parameter': 0.04, \
-                'number_of_iterations': noi, \
-                'tolerance_constant': 0.00001, \
-                'methodTV': 0, \
-                'nonneg': 0, \
-                'printingOut': 0, \
-                'time_marching_parameter': 0.00002
-                }
-        return Im, pars
+        #pars = {'algorithm': alg, \
+        #        'input': u0, \
+        #        'regularisation_parameter': 0.04, \
+        #        'number_of_iterations': noi, \
+        #        'tolerance_constant': 0.00001, \
+        #        'methodTV': 0, \
+        #        'nonneg': 0, \
+        #        'printingOut': 0, \
+        #        'time_marching_parameter': 0.00002
+        #        }
+        return Im,u0,u_ref
 
 
     def test_FGP_TV_CPU(self):
-        Im, pars = self.getPars(FGP_TV)
+        Im,input,ref = self.getPars()
 
-        fgp_cpu = FGP_TV(pars['input'],
-                         pars['regularisation_parameter'],
-                         pars['number_of_iterations'],
-                         pars['tolerance_constant'],
-                         pars['methodTV'],
-                         pars['nonneg'],
-                         pars['printingOut'], 'cpu')
+        fgp_cpu = FGP_TV(input,0.04,1200,1e-5,0,0,0,'cpu');
 
         rms = rmse(Im, fgp_cpu)
-        pars['rmse'] = rms
+
         self.assertAlmostEqual(rms,0.02,delta=0.01)
 
     def test_TV_ROF_CPU(self):
         # set parameters
-        Im, pars = self.getPars(ROF_TV)
+        Im, input,ref = self.getPars()
         # call routine
-        fgp_cpu = ROF_TV(pars['input'],
-                         pars['regularisation_parameter'],
-                         pars['number_of_iterations'],
-                         pars['time_marching_parameter'], 'cpu')
+        fgp_cpu = ROF_TV(input,0.04,1200,2e-5, 'cpu')
 
         rms = rmse(Im, fgp_cpu)
-        pars['rmse'] = rms
 
-        #txtstr = printParametersToString(pars)
-        #print(txtstr)
         # now test that it generates some expected output
         self.assertAlmostEqual(rms,0.02,delta=0.01)
 
     def test_SB_TV_CPU(self):
         # set parameters
-        Im, pars = self.getPars(SB_TV)
+        Im, input,ref = self.getPars()
         # call routine
-        fgp_cpu = SB_TV(pars['input'],
-                         pars['regularisation_parameter'],
-                         pars['number_of_iterations'],
-                         pars['time_marching_parameter'], 'cpu')
+        sb_cpu = SB_TV(input,0.04,150,1e-5,0,0,'cpu')
 
-        rms = rmse(Im, fgp_cpu)
-        pars['rmse'] = rms
+        rms = rmse(Im, sb_cpu)
+
+        # now test that it generates some expected output
+        self.assertAlmostEqual(rms,0.02,delta=0.01)
+
+    def test_TGV_CPU(self):
+        # set parameters
+        Im, input,ref = self.getPars()
+        # call routine
+        sb_cpu = TGV(input,0.04,1.0,2.0,250,12,'cpu')
+
+        rms = rmse(Im, sb_cpu)
 
-        #txtstr = printParametersToString(pars)
-        #print(txtstr)
         # now test that it generates some expected output
         self.assertAlmostEqual(rms,0.02,delta=0.01)
+
+    def test_LLT_ROF_CPU(self):
+        # set parameters
+        Im, input,ref = self.getPars()
+        # call routine
+        sb_cpu = LLT_ROF(input,0.04,0.01,1000,1e-4,'cpu')
+
+        rms = rmse(Im, sb_cpu)
+
+        # now test that it generates some expected output
+        self.assertAlmostEqual(rms,0.02,delta=0.01)
+
+    def test_NDF_CPU(self):
+        # set parameters
+        Im, input,ref = self.getPars()
+        # call routine
+        sb_cpu = NDF(input, 0.06, 0.04,1000,0.025,1, 'cpu')
+
+        rms = rmse(Im, sb_cpu)
+
+        # now test that it generates some expected output
+        self.assertAlmostEqual(rms, 0.02, delta=0.01)
+
+    def test_Diff4th_CPU(self):
+        # set parameters
+        Im, input,ref = self.getPars()
+        # call routine
+        sb_cpu = Diff4th(input, 3.5,0.02,500,0.001, 'cpu')
+
+        rms = rmse(Im, sb_cpu)
+
+        # now test that it generates some expected output
+        self.assertAlmostEqual(rms, 0.02, delta=0.01)
+
+    def test_FGP_dTV_CPU(self):
+        # set parameters
+        Im, input,ref = self.getPars()
+        # call routine
+        sb_cpu = FGP_dTV(input,ref,0.04,1000,1e-7,0.2,0,0,0, 'cpu')
+
+        rms = rmse(Im, sb_cpu)
+
+        # now test that it generates some expected output
+        self.assertAlmostEqual(rms, 0.02, delta=0.01)
diff --git a/test/test_FGP_TV.py b/test/test_FGP_TV.py
deleted file mode 100644
index f0dc540..0000000
--- a/test/test_FGP_TV.py
+++ /dev/null
@@ -1,152 +0,0 @@
-import unittest
-import math
-import os
-import timeit
-from ccpi.filters.regularisers import FGP_TV
-#, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th
-from testroutines import *
-
-###############################################################################
-
-class TestRegularisers(unittest.TestCase):
-
-    def test_FGP_TV_CPU(self):
-        print(__name__)
-        filename = os.path.join("lena_gray_512.tif")
-        plt = TiffReader()
-        # read image
-        Im = plt.imread(filename)
-        Im = np.asarray(Im, dtype='float32')
-
-        Im = Im / 255
-        perc = 0.05
-        u0 = Im + np.random.normal(loc=0,
-                                   scale=perc * Im,
-                                   size=np.shape(Im))
-        u_ref = Im + np.random.normal(loc=0,
-                                      scale=0.01 * 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 = u0.astype('float32')
-        u_ref = u_ref.astype('float32')
-
-        print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-        print("____________FGP-TV bench___________________")
-        print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-        # set parameters
-        pars = {'algorithm': FGP_TV, \
-                'input': u0, \
-                'regularisation_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['regularisation_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)
-        self.assertTrue(math.isclose(rms,0.02,rel_tol=1e-1))
-
-    def test_FGP_TV_CPU_vs_GPU(self):
-        print(__name__)
-        filename = os.path.join("lena_gray_512.tif")
-        plt = TiffReader()
-        # read image
-        Im = plt.imread(filename)
-        Im = np.asarray(Im, dtype='float32')
-
-        Im = Im / 255
-        perc = 0.05
-        u0 = Im + np.random.normal(loc=0,
-                                   scale=perc * Im,
-                                   size=np.shape(Im))
-        u_ref = Im + np.random.normal(loc=0,
-                                      scale=0.01 * 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 = u0.astype('float32')
-        u_ref = u_ref.astype('float32')
-
-        print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-        print("____________FGP-TV bench___________________")
-        print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-        # set parameters
-        pars = {'algorithm': FGP_TV, \
-                'input': u0, \
-                'regularisation_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['regularisation_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)
-
-        print("##############FGP TV GPU##################")
-        start_time = timeit.default_timer()
-        try:
-            fgp_gpu = FGP_TV(pars['input'],
-                             pars['regularisation_parameter'],
-                             pars['number_of_iterations'],
-                             pars['tolerance_constant'],
-                             pars['methodTV'],
-                             pars['nonneg'],
-                             pars['printingOut'], 'gpu')
-
-        except ValueError as ve:
-            self.skipTest("Results not comparable. GPU computing error.")
-
-        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)
-
-        print("--------Compare the results--------")
-        tolerance = 1e-05
-        diff_im = np.zeros(np.shape(fgp_cpu))
-        diff_im = abs(fgp_cpu - fgp_gpu)
-        diff_im[diff_im > tolerance] = 1
-
-        self.assertLessEqual(diff_im.sum(), 1)
-
-if __name__ == '__main__':
-    unittest.main()
diff --git a/test/test_ROF_TV.py b/test/test_ROF_TV.py
deleted file mode 100644
index fa35680..0000000
--- a/test/test_ROF_TV.py
+++ /dev/null
@@ -1,124 +0,0 @@
-import unittest
-import math
-import os
-import timeit
-from ccpi.filters.regularisers import ROF_TV
-#, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th
-from testroutines import *
-
-class TestRegularisers(unittest.TestCase):
-
-    def test_ROF_TV_CPU(self):
-        filename = os.path.join("lena_gray_512.tif")
-        plt = TiffReader()
-        # read image
-        Im = plt.imread(filename)
-        Im = np.asarray(Im, dtype='float32')
-
-        Im = Im / 255
-        perc = 0.05
-        u0 = Im + np.random.normal(loc=0,
-                                   scale=perc * Im,
-                                   size=np.shape(Im))
-        u_ref = Im + np.random.normal(loc=0,
-                                      scale=0.01 * 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 = u0.astype('float32')
-        u_ref = u_ref.astype('float32')
-
-
-        # set parameters
-        pars = {'algorithm': ROF_TV, \
-                'input': u0, \
-                'regularisation_parameter': 0.04, \
-                'number_of_iterations': 2500, \
-                'time_marching_parameter': 0.00002
-                }
-        print("#############ROF TV CPU####################")
-        start_time = timeit.default_timer()
-        rof_cpu = ROF_TV(pars['input'],
-                         pars['regularisation_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)
-
-        self.assertTrue(math.isclose(rms,0.02067839,rel_tol=1e-2))
-
-
-    def test_ROF_TV_CPU_vs_GPU(self):
-        filename = os.path.join("lena_gray_512.tif")
-        plt = TiffReader()
-        # read image
-        Im = plt.imread(filename)
-        Im = np.asarray(Im, dtype='float32')
-
-        Im = Im / 255
-        perc = 0.05
-        u0 = Im + np.random.normal(loc=0,
-                                   scale=perc * Im,
-                                   size=np.shape(Im))
-        u_ref = Im + np.random.normal(loc=0,
-                                      scale=0.01 * 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 = u0.astype('float32')
-        u_ref = u_ref.astype('float32')
-
-        print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-        print("____________ROF-TV bench___________________")
-        print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-        # set parameters
-        pars = {'algorithm': ROF_TV, \
-                'input': u0, \
-                'regularisation_parameter': 0.04, \
-                'number_of_iterations': 2500, \
-                'time_marching_parameter': 0.00002
-                }
-        print("##############ROF TV GPU##################")
-        start_time = timeit.default_timer()
-        try:
-            rof_gpu = ROF_TV(pars['input'],
-                             pars['regularisation_parameter'],
-                             pars['number_of_iterations'],
-                             pars['time_marching_parameter'], 'gpu')
-        except ValueError as ve:
-            self.skipTest("Results not comparable. GPU computing error.")
-
-        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)
-
-        print("#############ROF TV CPU####################")
-        start_time = timeit.default_timer()
-        rof_cpu = ROF_TV(pars['input'],
-                         pars['regularisation_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)
-        print("--------Compare the results--------")
-        tolerance = 1e-04
-        diff_im = np.zeros(np.shape(rof_cpu))
-        diff_im = abs(rof_cpu - rof_gpu)
-        diff_im[diff_im > tolerance] = 1
-        self.assertLessEqual(diff_im.sum(), 1)
-
-if __name__ == '__main__':
-    unittest.main()