GPU test fixed, cpu vs gpu test added for ROF-TV

5 files changed, 138 insertions, 15 deletions

diff --git a/Core/regularizers_GPU/TV_ROF/TV_ROF_GPU.cu b/Core/regularizers_GPU/TV_ROF/TV_ROF_GPU.cu
index 73a52e1..897b5d0 100755
--- a/Core/regularizers_GPU/TV_ROF/TV_ROF_GPU.cu
+++ b/Core/regularizers_GPU/TV_ROF/TV_ROF_GPU.cu
@@ -312,7 +312,7 @@ extern "C" void TV_ROF_GPU_kernel(float* Input, float* Output, int N, int M, int
 		int dev = 0;
 		CHECK(cudaSetDevice(dev));
 		
-        float *d_input, *d_update, *d_D1, *d_D2;
+        float *d_input, *d_update, *d_D1, *d_D2;        
         
         CHECK(cudaMalloc((void**)&d_input,N*M*Z*sizeof(float)));
         CHECK(cudaMalloc((void**)&d_update,N*M*Z*sizeof(float)));
diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py
index 5908c3c..d147b85 100644
--- a/Wrappers/Python/demo/test_cpu_regularizers.py
+++ b/Wrappers/Python/demo/test_cpu_regularizers.py
@@ -284,9 +284,9 @@ start_time = timeit.default_timer()
 
 pars = {'algorithm': ROF_TV , \
         'input' : u0,\
-        'regularization_parameter':1,\
-        'marching_step': 0.003,\
-        'number_of_iterations': 300
+        'regularization_parameter':25,\
+        'marching_step': 0.001,\
+        'number_of_iterations': 350
         }
 rof = ROF_TV(pars['input'],
              pars['number_of_iterations'],
@@ -307,9 +307,7 @@ props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
 # place a text box in upper left in axes coords
 a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
          verticalalignment='top', bbox=props)
-imgplot = plt.imshow(tgv, cmap="gray")
-
-
+imgplot = plt.imshow(rof, cmap="gray")
 
 plt.show()
 
diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx
index fcb91cc..c724471 100644
--- a/Wrappers/Python/src/gpu_regularizers.pyx
+++ b/Wrappers/Python/src/gpu_regularizers.pyx
@@ -345,7 +345,7 @@ def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
     # Running CUDA code here    
     TV_ROF_GPU_kernel(            
             &inputData[0,0], &B[0,0], 
-                       dims[0], dims[1], 0, 
+                       dims[0], dims[1], 1, 
                        iterations , 
                        time_marching_parameter, 
                        regularization_parameter);   
diff --git a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py
new file mode 100644
index 0000000..8c91c73
--- /dev/null
+++ b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py
@@ -0,0 +1,123 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Feb 22 11:39:43 2018
+
+Testing CPU implementation against GPU one
+
+@author: algol
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import os    
+import timeit
+from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, GPU_ROF_TV
+from ccpi.filters.cpu_regularizers_cython import ROF_TV
+###############################################################################
+def printParametersToString(pars):
+        txt = r''
+        for key, value in pars.items():
+            if key== 'algorithm' :
+                txt += "{0} = {1}".format(key, value.__name__)
+            elif key == 'input':
+                txt += "{0} = {1}".format(key, np.shape(value))
+            else:
+                txt += "{0} = {1}".format(key, value)
+            txt += '\n'
+        return txt
+###############################################################################
+def rmse(im1, im2):
+    a, b = im1.shape
+    rmse = np.sqrt(np.sum((im1 - im2) ** 2) / float(a * b))    
+    return rmse
+        
+filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
+
+# read image
+Im = plt.imread(filename)                     
+Im = np.asarray(Im, dtype='float32')
+
+Im = Im/255
+perc = 0.075
+u0 = Im + np.random.normal(loc = Im ,
+                                  scale = perc * Im , 
+                                  size = np.shape(Im))
+# map the u0 u0->u0>0
+f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+u0 = f(u0).astype('float32')
+
+## plot 
+fig = plt.figure(1)
+plt.suptitle('Comparison of ROF-TV regularizer using CPU and GPU implementations')
+a=fig.add_subplot(1,4,1)
+a.set_title('Noisy Image')
+imgplot = plt.imshow(u0,cmap="gray")
+
+
+# set parameters
+pars = {'algorithm': ROF_TV , \
+        'input' : u0,\
+        'regularization_parameter':12,\
+        'time_marching_parameter': 0.001,\
+        'number_of_iterations': 600
+        }
+print ("#################ROF TV CPU#####################")
+start_time = timeit.default_timer()
+rof_cpu = ROF_TV(pars['input'],
+             pars['number_of_iterations'],
+             pars['regularization_parameter'],
+             pars['time_marching_parameter'] 
+             )
+#tgv = out
+rms = rmse(Im, rof_cpu)
+pars['rmse'] = rms
+
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,4,2)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
+# place a text box in upper left in axes coords
+a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow(rof_cpu, cmap="gray")
+plt.title('{}'.format('CPU results'))
+
+
+print ("#################ROF TV GPU#####################")
+start_time = timeit.default_timer()
+rof_gpu = GPU_ROF_TV(pars['input'], 
+                     pars['number_of_iterations'], 
+                     pars['time_marching_parameter'], 
+                     pars['regularization_parameter'])
+                     
+rms = rmse(Im, rof_gpu)
+pars['rmse'] = rms
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,4,3)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
+# place a text box in upper left in axes coords
+a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow(rof_gpu, cmap="gray")
+plt.title('{}'.format('GPU results'))
+
+
+print ("--------Compare the results--------")
+tolerance = 1e-06
+diff_im = abs(rof_cpu - rof_gpu)
+diff_im[diff_im > tolerance] = 1
+a=fig.add_subplot(1,4,4)
+imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray")
+plt.title('{}'.format('Pixels larger threshold difference'))
+if (diff_im.sum() > 1):
+    print ("Arrays do not match!")
+else:
+    print ("Arrays match")
\ No newline at end of file
diff --git a/Wrappers/Python/test/test_gpu_regularizers.py b/Wrappers/Python/test/test_gpu_regularizers.py
index 9e37627..29f5bad 100644
--- a/Wrappers/Python/test/test_gpu_regularizers.py
+++ b/Wrappers/Python/test/test_gpu_regularizers.py
@@ -93,7 +93,8 @@ props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
 # place a text box in upper left in axes coords
 a.text(0.05, 0.95, 'd4h - u0', transform=a.transAxes, fontsize=12,
          verticalalignment='top', bbox=props)
-imgplot = plt.imshow((d4h - u0), cmap="gray")
+imgplot = plt.imshow((d4h - u0)**2,  vmin=0, vmax=0.03, cmap="gray")
+plt.colorbar(ticks=[0, 0.03], orientation='vertical')
 
 
 ## Patch Based Regul NML
@@ -142,7 +143,8 @@ props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
 # place a text box in upper left in axes coords
 a.text(0.05, 0.95, 'nml - u0', transform=a.transAxes, fontsize=14,
          verticalalignment='top', bbox=props)
-imgplot = plt.imshow((nml - u0), cmap="gray")
+imgplot = plt.imshow((nml - u0)**2,  vmin=0, vmax=0.03, cmap="gray")
+plt.colorbar(ticks=[0, 0.03], orientation='vertical')
 
         
         
@@ -152,9 +154,9 @@ start_time = timeit.default_timer()
 pars = {
 'algorithm' : GPU_ROF_TV , \
         'input' : u0,
-        'regularization_parameter': 1,\
-        'time_marching_parameter': 0.003, \
-        'number_of_iterations':300        
+        'regularization_parameter': 25,\
+        'time_marching_parameter': 0.001, \
+        'number_of_iterations':350        
 	}
 
 rof_tv = GPU_ROF_TV(pars['input'], 
@@ -183,6 +185,6 @@ props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
 # place a text box in upper left in axes coords
 a.text(0.05, 0.95, 'rof_tv - u0', transform=a.transAxes, fontsize=14,
          verticalalignment='top', bbox=props)
-imgplot = plt.imshow((rof_tv - u0), cmap="gray")
-
+imgplot = plt.imshow((rof_tv - u0)**2, vmin=0, vmax=0.03, cmap="gray")
+plt.colorbar(ticks=[0, 0.03], orientation='vertical')
 plt.show()