Merge pull request #102 from vais-ral/qualmet

Adds qualitymetrics class and updates demos

10 files changed, 216 insertions, 169 deletions

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/Wrappers/Python/demos/demo_cpu_inpainters.py b/Wrappers/Python/demos/demo_cpu_inpainters.py
index 3b4191b..c61ea50 100644
--- a/Wrappers/Python/demos/demo_cpu_inpainters.py
+++ b/Wrappers/Python/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/Wrappers/Python/demos/demo_cpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_regularisers.py
index e6befa9..b8dadf5 100644
--- a/Wrappers/Python/demos/demo_cpu_regularisers.py
+++ b/Wrappers/Python/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/Wrappers/Python/demos/demo_cpu_regularisers3D.py b/Wrappers/Python/demos/demo_cpu_regularisers3D.py
index 2d2fc22..df8af27 100644
--- a/Wrappers/Python/demos/demo_cpu_regularisers3D.py
+++ b/Wrappers/Python/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/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py
index 230a761..6c4ab5e 100644
--- a/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py
+++ b/Wrappers/Python/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/Wrappers/Python/demos/demo_gpu_regularisers.py b/Wrappers/Python/demos/demo_gpu_regularisers.py
index e1c6575..54a1c14 100644
--- a/Wrappers/Python/demos/demo_gpu_regularisers.py
+++ b/Wrappers/Python/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/Wrappers/Python/demos/demo_gpu_regularisers3D.py b/Wrappers/Python/demos/demo_gpu_regularisers3D.py
index b6058d2..d50c08e 100644
--- a/Wrappers/Python/demos/demo_gpu_regularisers3D.py
+++ b/Wrappers/Python/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/Wrappers/Python/demos/qualitymetrics.py b/Wrappers/Python/demos/qualitymetrics.py
deleted file mode 100644
index 850829e..0000000
--- a/Wrappers/Python/demos/qualitymetrics.py
+++ /dev/null
@@ -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/Wrappers/Python/setup-regularisers.py.in b/Wrappers/Python/setup-regularisers.py.in
index 462edda..16fb380 100644
--- a/Wrappers/Python/setup-regularisers.py.in
+++ b/Wrappers/Python/setup-regularisers.py.in
@@ -68,7 +68,7 @@ setup(
     
     ],
 	zip_safe = False,	
-	packages = {'ccpi','ccpi.filters'},
+	packages = {'ccpi','ccpi.filters', 'ccpi.supp'},
 )