From 49423f47c7282672fea979cd122c063c0f5ee465 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Wed, 21 Feb 2018 12:06:12 +0000 Subject: ROF_TV demo --- Wrappers/Python/test/test_gpu_regularizers.py | 44 +++++++++++++++++++++++++-- 1 file changed, 41 insertions(+), 3 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/test/test_gpu_regularizers.py b/Wrappers/Python/test/test_gpu_regularizers.py index 735a25d..3346fb2 100644 --- a/Wrappers/Python/test/test_gpu_regularizers.py +++ b/Wrappers/Python/test/test_gpu_regularizers.py @@ -6,14 +6,12 @@ Created on Tue Jan 30 10:24:26 2018 @author: ofn77899 """ - - import matplotlib.pyplot as plt import numpy as np import os from enum import Enum import timeit -from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML +from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, ROF_TV_GPU ############################################################################### def printParametersToString(pars): txt = r'' @@ -146,3 +144,43 @@ imgplot = plt.imshow((nml - u0)**2, cmap="gray") plt.show() + +## Rudin-Osher-Fatemi (ROF) TV regularization +start_time = timeit.default_timer() + +pars = { + 'input' : u0, + 'regularization_parameter': 1,\ + 'time_marching_parameter': 0.003, \ + 'number_of_iterations':300 + } + +rof_tv = ROF_TV_GPU(pars['input'], + pars['number_of_iterations'], + pars['time_marching_parameter'], + pars['number_of_iterations']) + +rms = rmse(Im, rof_tv) +pars['rmse'] = rms +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(2,4,1) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) +# place a text box in upper left in axes coords +a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=12, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(nml, cmap="gray") + +a=fig.add_subplot(2,4,2) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +# place a text box in upper left in axes coords +a.text(0.05, 0.95, 'rof_tv - u0', transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow((rof_tv - u0)**2, cmap="gray") + +plt.show() -- cgit v1.2.3 From 75917255b4f0b8aae2e6ce9492a75e0f749bfb3e Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 21 Feb 2018 12:49:59 +0000 Subject: added TV_ROF --- Wrappers/Python/src/gpu_regularizers.pyx | 8 ++++---- Wrappers/Python/test/test_gpu_regularizers.py | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx index 5a5d274..fcb91cc 100644 --- a/Wrappers/Python/src/gpu_regularizers.pyx +++ b/Wrappers/Python/src/gpu_regularizers.pyx @@ -25,7 +25,7 @@ cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec, int N, int M, int Z, int dimension, int SearchW, int SimilW, int SearchW_real, float denh2, float lambdaf); -cdef extern void TV_ROF_GPU(float* A, float* B, int N, int M, int Z, int iter, float tau, float lambdaf); +cdef extern void TV_ROF_GPU_kernel(float* A, float* B, int N, int M, int Z, int iter, float tau, float lambdaf); cdef extern float pad_crop(float *A, float *Ap, int OldSizeX, int OldSizeY, int OldSizeZ, int NewSizeX, int NewSizeY, int NewSizeZ, @@ -67,7 +67,7 @@ def NML(inputData, h, lambdaf) -def ROF_TV_GPU(inputData, +def GPU_ROF_TV(inputData, iterations, time_marching_parameter, regularization_parameter): @@ -343,7 +343,7 @@ def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, np.zeros([dims[0],dims[1]], dtype='float32') # Running CUDA code here - TV_ROF_GPU( + TV_ROF_GPU_kernel( &inputData[0,0], &B[0,0], dims[0], dims[1], 0, iterations , @@ -366,7 +366,7 @@ def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, np.zeros([dims[0],dims[1],dims[2]], dtype='float32') # Running CUDA code here - TV_ROF_GPU( + TV_ROF_GPU_kernel( &inputData[0,0,0], &B[0,0,0], dims[0], dims[1], dims[2], iterations , diff --git a/Wrappers/Python/test/test_gpu_regularizers.py b/Wrappers/Python/test/test_gpu_regularizers.py index 735a25d..75127a6 100644 --- a/Wrappers/Python/test/test_gpu_regularizers.py +++ b/Wrappers/Python/test/test_gpu_regularizers.py @@ -13,7 +13,7 @@ import numpy as np import os from enum import Enum import timeit -from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML +from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, GPU_ROF_TV ############################################################################### def printParametersToString(pars): txt = r'' -- cgit v1.2.3 From ce313c7dcf45edaf4c824690d8caa7df5df90120 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 21 Feb 2018 14:43:35 +0000 Subject: modified gpu regularizers build --- Wrappers/Python/setup-regularizers.py.in | 6 ++++-- Wrappers/Python/test/test_gpu_regularizers.py | 30 ++++++++++++++------------- 2 files changed, 20 insertions(+), 16 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/setup-regularizers.py.in b/Wrappers/Python/setup-regularizers.py.in index a125261..8655a2e 100644 --- a/Wrappers/Python/setup-regularizers.py.in +++ b/Wrappers/Python/setup-regularizers.py.in @@ -34,7 +34,9 @@ extra_libraries = ['cilreg'] extra_include_dirs += [os.path.join(".." , ".." , "Core"), os.path.join(".." , ".." , "Core", "regularizers_CPU"), - os.path.join(".." , ".." , "Core", "regularizers_GPU") , + os.path.join(".." , ".." , "Core", "regularizers_GPU" , "Diffus_HO" ) , + os.path.join(".." , ".." , "Core", "regularizers_GPU" , "NL_Regul" ) , + os.path.join(".." , ".." , "Core", "regularizers_GPU" , "TV_ROF" ) , "."] if platform.system() == 'Windows': @@ -81,4 +83,4 @@ setup( ) -@SETUP_GPU_WRAPPERS@ \ No newline at end of file +@SETUP_GPU_WRAPPERS@ diff --git a/Wrappers/Python/test/test_gpu_regularizers.py b/Wrappers/Python/test/test_gpu_regularizers.py index ff8af44..9e37627 100644 --- a/Wrappers/Python/test/test_gpu_regularizers.py +++ b/Wrappers/Python/test/test_gpu_regularizers.py @@ -50,14 +50,15 @@ u0 = f(u0).astype('float32') ## plot fig = plt.figure() -a=fig.add_subplot(2,3,1) +a=fig.add_subplot(2,4,1) a.set_title('noise') imgplot = plt.imshow(u0,cmap="gray") ## Diff4thHajiaboli start_time = timeit.default_timer() -pars = {'algorithm' : Diff4thHajiaboli , \ +pars = { +'algorithm' : Diff4thHajiaboli , \ 'input' : u0, 'edge_preserv_parameter':0.1 , \ 'number_of_iterations' :250 ,\ @@ -76,7 +77,7 @@ pars['rmse'] = rms txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) -a=fig.add_subplot(2,3,2) +a=fig.add_subplot(2,4,2) # these are matplotlib.patch.Patch properties props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) @@ -85,14 +86,14 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=12, verticalalignment='top', bbox=props) imgplot = plt.imshow(d4h, cmap="gray") -a=fig.add_subplot(2,3,5) +a=fig.add_subplot(2,4,6) # these are matplotlib.patch.Patch properties props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) # place a text box in upper left in axes coords a.text(0.05, 0.95, 'd4h - u0', transform=a.transAxes, fontsize=12, verticalalignment='top', bbox=props) -imgplot = plt.imshow((d4h - u0)**2, cmap="gray") +imgplot = plt.imshow((d4h - u0), cmap="gray") ## Patch Based Regul NML @@ -107,6 +108,7 @@ pars = {'algorithm' : NML , \ } """ pars = { +'algorithm' : NML , \ 'input' : u0, 'regularization_parameter': 0.01,\ 'searching_window_ratio':3, \ @@ -124,7 +126,7 @@ pars['rmse'] = rms txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) -a=fig.add_subplot(2,3,3) +a=fig.add_subplot(2,4,3) # these are matplotlib.patch.Patch properties props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) @@ -133,22 +135,22 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=12, verticalalignment='top', bbox=props) imgplot = plt.imshow(nml, cmap="gray") -a=fig.add_subplot(2,3,6) +a=fig.add_subplot(2,4,7) # these are matplotlib.patch.Patch properties props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) # place a text box in upper left in axes coords a.text(0.05, 0.95, 'nml - u0', transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) -imgplot = plt.imshow((nml - u0)**2, cmap="gray") +imgplot = plt.imshow((nml - u0), cmap="gray") -plt.show() ## Rudin-Osher-Fatemi (ROF) TV regularization start_time = timeit.default_timer() pars = { +'algorithm' : GPU_ROF_TV , \ 'input' : u0, 'regularization_parameter': 1,\ 'time_marching_parameter': 0.003, \ @@ -158,29 +160,29 @@ pars = { rof_tv = GPU_ROF_TV(pars['input'], pars['number_of_iterations'], pars['time_marching_parameter'], - pars['number_of_iterations']) + pars['regularization_parameter']) rms = rmse(Im, rof_tv) pars['rmse'] = rms txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) -a=fig.add_subplot(2,4,1) +a=fig.add_subplot(2,4,4) # these are matplotlib.patch.Patch properties props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) # place a text box in upper left in axes coords a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=12, verticalalignment='top', bbox=props) -imgplot = plt.imshow(nml, cmap="gray") +imgplot = plt.imshow(rof_tv, cmap="gray") -a=fig.add_subplot(2,4,2) +a=fig.add_subplot(2,4,8) # these are matplotlib.patch.Patch properties props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) # place a text box in upper left in axes coords a.text(0.05, 0.95, 'rof_tv - u0', transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) -imgplot = plt.imshow((rof_tv - u0)**2, cmap="gray") +imgplot = plt.imshow((rof_tv - u0), cmap="gray") plt.show() -- cgit v1.2.3 From d4d12945f5396b10259fdeeb2b09f99b0e2c6afd Mon Sep 17 00:00:00 2001 From: algol Date: Thu, 22 Feb 2018 12:34:39 +0000 Subject: GPU test fixed, cpu vs gpu test added for ROF-TV --- Wrappers/Python/demo/test_cpu_regularizers.py | 10 +- Wrappers/Python/src/gpu_regularizers.pyx | 2 +- .../Python/test/test_cpu_vs_gpu_regularizers.py | 123 +++++++++++++++++++++ Wrappers/Python/test/test_gpu_regularizers.py | 16 +-- 4 files changed, 137 insertions(+), 14 deletions(-) create mode 100644 Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py (limited to 'Wrappers') 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() -- cgit v1.2.3 From 0ebf1f949b7150692e356627c455e905b642af97 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Sun, 25 Feb 2018 17:54:03 +0000 Subject: updates to CPU/GPU core for ROF and demos --- Wrappers/Python/demo/test_cpu_regularizers.py | 6 +++--- Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py | 6 +++--- Wrappers/Python/test/test_gpu_regularizers.py | 6 +++--- 3 files changed, 9 insertions(+), 9 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py index d147b85..b08c418 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':25,\ - 'marching_step': 0.001,\ - 'number_of_iterations': 350 + 'regularization_parameter':0.04,\ + 'marching_step': 0.0025,\ + 'number_of_iterations': 300 } rof = ROF_TV(pars['input'], pars['number_of_iterations'], diff --git a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py index 8c91c73..d742a7f 100644 --- a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py +++ b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py @@ -58,8 +58,8 @@ imgplot = plt.imshow(u0,cmap="gray") # set parameters pars = {'algorithm': ROF_TV , \ 'input' : u0,\ - 'regularization_parameter':12,\ - 'time_marching_parameter': 0.001,\ + 'regularization_parameter':0.04,\ + 'time_marching_parameter': 0.0025,\ 'number_of_iterations': 600 } print ("#################ROF TV CPU#####################") @@ -120,4 +120,4 @@ 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 + print ("Arrays match") diff --git a/Wrappers/Python/test/test_gpu_regularizers.py b/Wrappers/Python/test/test_gpu_regularizers.py index 29f5bad..c473270 100644 --- a/Wrappers/Python/test/test_gpu_regularizers.py +++ b/Wrappers/Python/test/test_gpu_regularizers.py @@ -154,9 +154,9 @@ start_time = timeit.default_timer() pars = { 'algorithm' : GPU_ROF_TV , \ 'input' : u0, - 'regularization_parameter': 25,\ - 'time_marching_parameter': 0.001, \ - 'number_of_iterations':350 + 'regularization_parameter': 0.04,\ + 'time_marching_parameter': 0.0025, \ + 'number_of_iterations':300 } rof_tv = GPU_ROF_TV(pars['input'], -- cgit v1.2.3 From bcdf186ccdca54a3df383512ad5a117004500a60 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Mon, 26 Feb 2018 12:50:58 +0000 Subject: CPU-GPU naming consistency --- Wrappers/Python/src/cpu_regularizers.cpp | 546 ++++++++++----------- Wrappers/Python/src/gpu_regularizers.pyx | 69 ++- Wrappers/Python/test/test_cpu_vs_gpu.py | 10 + .../Python/test/test_cpu_vs_gpu_regularizers.py | 8 +- 4 files changed, 353 insertions(+), 280 deletions(-) create mode 100644 Wrappers/Python/test/test_cpu_vs_gpu.py (limited to 'Wrappers') diff --git a/Wrappers/Python/src/cpu_regularizers.cpp b/Wrappers/Python/src/cpu_regularizers.cpp index e311570..43d5d11 100644 --- a/Wrappers/Python/src/cpu_regularizers.cpp +++ b/Wrappers/Python/src/cpu_regularizers.cpp @@ -27,7 +27,7 @@ limitations under the License. #include "boost/tuple/tuple.hpp" #include "SplitBregman_TV_core.h" -#include "FGP_TV_core.h" +//#include "FGP_TV_core.h" #include "LLT_model_core.h" #include "PatchBased_Regul_core.h" #include "TGV_PD_core.h" @@ -305,289 +305,289 @@ bp::list SplitBregman_TV(np::ndarray input, double d_mu, int iter, double d_epsi -bp::list FGP_TV(np::ndarray input, double d_mu, int iter, double d_epsil, int methTV) { +//bp::list FGP_TV(np::ndarray input, double d_mu, int iter, double d_epsil, int methTV) { - // the result is in the following list - bp::list result; + //// the result is in the following list + //bp::list result; - int number_of_dims, dimX, dimY, dimZ, ll, j, count; - float *A, *D = NULL, *D_old = NULL, *P1 = NULL, *P2 = NULL, *P3 = NULL, *P1_old = NULL, *P2_old = NULL, *P3_old = NULL, *R1 = NULL, *R2 = NULL, *R3 = NULL; - float lambda, tk, tkp1, re, re1, re_old, epsil, funcval; + //int number_of_dims, dimX, dimY, dimZ, ll, j, count; + //float *A, *D = NULL, *D_old = NULL, *P1 = NULL, *P2 = NULL, *P3 = NULL, *P1_old = NULL, *P2_old = NULL, *P3_old = NULL, *R1 = NULL, *R2 = NULL, *R3 = NULL; + //float lambda, tk, tkp1, re, re1, re_old, epsil, funcval; - //number_of_dims = mxGetNumberOfDimensions(prhs[0]); - //dim_array = mxGetDimensions(prhs[0]); - - number_of_dims = input.get_nd(); - int dim_array[3]; + ////number_of_dims = mxGetNumberOfDimensions(prhs[0]); + ////dim_array = mxGetDimensions(prhs[0]); - dim_array[0] = input.shape(0); - dim_array[1] = input.shape(1); - if (number_of_dims == 2) { - dim_array[2] = -1; - } - else { - dim_array[2] = input.shape(2); - } - // Parameter handling is be done in Python - ///*Handling Matlab input data*/ - //if ((nrhs < 2) || (nrhs > 5)) mexErrMsgTxt("At least 2 parameters is required: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1')"); - - ///*Handling Matlab input data*/ - //A = (float *)mxGetData(prhs[0]); /*noisy image (2D/3D) */ - A = reinterpret_cast(input.get_data()); - - //mu = (float)mxGetScalar(prhs[1]); /* regularization parameter */ - lambda = (float)d_mu; - - //iter = 35; /* default iterations number */ + //number_of_dims = input.get_nd(); + //int dim_array[3]; - //epsil = 0.0001; /* default tolerance constant */ - epsil = (float)d_epsil; - //methTV = 0; /* default isotropic TV penalty */ - //if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5)) iter = (int)mxGetScalar(prhs[2]); /* iterations number */ - //if ((nrhs == 4) || (nrhs == 5)) epsil = (float)mxGetScalar(prhs[3]); /* tolerance constant */ - //if (nrhs == 5) { - // char *penalty_type; - // penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */ - // if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',"); - // if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */ - // mxFree(penalty_type); + //dim_array[0] = input.shape(0); + //dim_array[1] = input.shape(1); + //if (number_of_dims == 2) { + //dim_array[2] = -1; //} - //if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input image must be in a single precision"); } - - //plhs[1] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL); - bp::tuple shape1 = bp::make_tuple(dim_array[0], dim_array[1]); - np::dtype dtype = np::dtype::get_builtin(); - np::ndarray out1 = np::zeros(shape1, dtype); + //else { + //dim_array[2] = input.shape(2); + //} + //// Parameter handling is be done in Python + /////*Handling Matlab input data*/ + ////if ((nrhs < 2) || (nrhs > 5)) mexErrMsgTxt("At least 2 parameters is required: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1')"); + + /////*Handling Matlab input data*/ + ////A = (float *)mxGetData(prhs[0]); /*noisy image (2D/3D) */ + //A = reinterpret_cast(input.get_data()); + + ////mu = (float)mxGetScalar(prhs[1]); /* regularization parameter */ + //lambda = (float)d_mu; + + ////iter = 35; /* default iterations number */ + + ////epsil = 0.0001; /* default tolerance constant */ + //epsil = (float)d_epsil; + ////methTV = 0; /* default isotropic TV penalty */ + ////if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5)) iter = (int)mxGetScalar(prhs[2]); /* iterations number */ + ////if ((nrhs == 4) || (nrhs == 5)) epsil = (float)mxGetScalar(prhs[3]); /* tolerance constant */ + ////if (nrhs == 5) { + //// char *penalty_type; + //// penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */ + //// if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',"); + //// if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */ + //// mxFree(penalty_type); + ////} + ////if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input image must be in a single precision"); } + + ////plhs[1] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL); + //bp::tuple shape1 = bp::make_tuple(dim_array[0], dim_array[1]); + //np::dtype dtype = np::dtype::get_builtin(); + //np::ndarray out1 = np::zeros(shape1, dtype); - //float *funcvalA = (float *)mxGetData(plhs[1]); - float * funcvalA = reinterpret_cast(out1.get_data()); - //if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input image must be in a single precision"); } - - /*Handling Matlab output data*/ - dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; - - tk = 1.0f; - tkp1 = 1.0f; - count = 1; - re_old = 0.0f; - - if (number_of_dims == 2) { - dimZ = 1; /*2D case*/ - /*D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - D_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - P1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - P2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - P1_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - P2_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - R1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - R2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));*/ - - bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1]); - np::dtype dtype = np::dtype::get_builtin(); - - - np::ndarray npD = np::zeros(shape, dtype); - np::ndarray npD_old = np::zeros(shape, dtype); - np::ndarray npP1 = np::zeros(shape, dtype); - np::ndarray npP2 = np::zeros(shape, dtype); - np::ndarray npP1_old = np::zeros(shape, dtype); - np::ndarray npP2_old = np::zeros(shape, dtype); - np::ndarray npR1 = np::zeros(shape, dtype); - np::ndarray npR2 = np::zeros(shape, dtype); - - D = reinterpret_cast(npD.get_data()); - D_old = reinterpret_cast(npD_old.get_data()); - P1 = reinterpret_cast(npP1.get_data()); - P2 = reinterpret_cast(npP2.get_data()); - P1_old = reinterpret_cast(npP1_old.get_data()); - P2_old = reinterpret_cast(npP2_old.get_data()); - R1 = reinterpret_cast(npR1.get_data()); - R2 = reinterpret_cast(npR2.get_data()); - - /* begin iterations */ - for (ll = 0; ll(out1.get_data()); + ////if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input image must be in a single precision"); } + + ///*Handling Matlab output data*/ + //dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; + + //tk = 1.0f; + //tkp1 = 1.0f; + //count = 1; + //re_old = 0.0f; + + //if (number_of_dims == 2) { + //dimZ = 1; /*2D case*/ + ///*D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + //D_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + //P1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + //P2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + //P1_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + //P2_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + //R1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + //R2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));*/ + + //bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1]); + //np::dtype dtype = np::dtype::get_builtin(); + + + //np::ndarray npD = np::zeros(shape, dtype); + //np::ndarray npD_old = np::zeros(shape, dtype); + //np::ndarray npP1 = np::zeros(shape, dtype); + //np::ndarray npP2 = np::zeros(shape, dtype); + //np::ndarray npP1_old = np::zeros(shape, dtype); + //np::ndarray npP2_old = np::zeros(shape, dtype); + //np::ndarray npR1 = np::zeros(shape, dtype); + //np::ndarray npR2 = np::zeros(shape, dtype); + + //D = reinterpret_cast(npD.get_data()); + //D_old = reinterpret_cast(npD_old.get_data()); + //P1 = reinterpret_cast(npP1.get_data()); + //P2 = reinterpret_cast(npP2.get_data()); + //P1_old = reinterpret_cast(npP1_old.get_data()); + //P2_old = reinterpret_cast(npP2_old.get_data()); + //R1 = reinterpret_cast(npR1.get_data()); + //R2 = reinterpret_cast(npR2.get_data()); + + ///* begin iterations */ + //for (ll = 0; ll 4) { - Obj_func2D(A, D, P1, P2, lambda, dimX, dimY); - funcval = 0.0f; - for (j = 0; j 2) { - if (re > re_old) { - Obj_func2D(A, D, P1, P2, lambda, dimX, dimY); - funcval = 0.0f; - for (j = 0; j(npD); - result.append(out1); - result.append(ll); - } - if (number_of_dims == 3) { - /*D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - D_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P1_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P2_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P3_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - R1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - R2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - R3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));*/ - bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1], dim_array[2]); - np::dtype dtype = np::dtype::get_builtin(); + ///*updating R and t*/ + //tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f; + //Rupd_func2D(P1, P1_old, P2, P2_old, R1, R2, tkp1, tk, dimX, dimY); + + ///* calculate norm */ + //re = 0.0f; re1 = 0.0f; + //for (j = 0; j 4) { + //Obj_func2D(A, D, P1, P2, lambda, dimX, dimY); + //funcval = 0.0f; + //for (j = 0; j 2) { + //if (re > re_old) { + //Obj_func2D(A, D, P1, P2, lambda, dimX, dimY); + //funcval = 0.0f; + //for (j = 0; j(npD); + //result.append(out1); + //result.append(ll); + //} + //if (number_of_dims == 3) { + ///*D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + //D_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + //P1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + //P2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + //P3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + //P1_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + //P2_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + //P3_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + //R1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + //R2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + //R3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));*/ + //bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1], dim_array[2]); + //np::dtype dtype = np::dtype::get_builtin(); - np::ndarray npD = np::zeros(shape, dtype); - np::ndarray npD_old = np::zeros(shape, dtype); - np::ndarray npP1 = np::zeros(shape, dtype); - np::ndarray npP2 = np::zeros(shape, dtype); - np::ndarray npP3 = np::zeros(shape, dtype); - np::ndarray npP1_old = np::zeros(shape, dtype); - np::ndarray npP2_old = np::zeros(shape, dtype); - np::ndarray npP3_old = np::zeros(shape, dtype); - np::ndarray npR1 = np::zeros(shape, dtype); - np::ndarray npR2 = np::zeros(shape, dtype); - np::ndarray npR3 = np::zeros(shape, dtype); - - D = reinterpret_cast(npD.get_data()); - D_old = reinterpret_cast(npD_old.get_data()); - P1 = reinterpret_cast(npP1.get_data()); - P2 = reinterpret_cast(npP2.get_data()); - P3 = reinterpret_cast(npP3.get_data()); - P1_old = reinterpret_cast(npP1_old.get_data()); - P2_old = reinterpret_cast(npP2_old.get_data()); - P3_old = reinterpret_cast(npP3_old.get_data()); - R1 = reinterpret_cast(npR1.get_data()); - R2 = reinterpret_cast(npR2.get_data()); - R3 = reinterpret_cast(npR3.get_data()); - /* begin iterations */ - for (ll = 0; ll 3) { - Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ); - funcval = 0.0f; - for (j = 0; j 2) { - if (re > re_old) { - Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ); - funcval = 0.0f; - for (j = 0; j(npD); - result.append(out1); - result.append(ll); - } + //np::ndarray npD = np::zeros(shape, dtype); + //np::ndarray npD_old = np::zeros(shape, dtype); + //np::ndarray npP1 = np::zeros(shape, dtype); + //np::ndarray npP2 = np::zeros(shape, dtype); + //np::ndarray npP3 = np::zeros(shape, dtype); + //np::ndarray npP1_old = np::zeros(shape, dtype); + //np::ndarray npP2_old = np::zeros(shape, dtype); + //np::ndarray npP3_old = np::zeros(shape, dtype); + //np::ndarray npR1 = np::zeros(shape, dtype); + //np::ndarray npR2 = np::zeros(shape, dtype); + //np::ndarray npR3 = np::zeros(shape, dtype); + + //D = reinterpret_cast(npD.get_data()); + //D_old = reinterpret_cast(npD_old.get_data()); + //P1 = reinterpret_cast(npP1.get_data()); + //P2 = reinterpret_cast(npP2.get_data()); + //P3 = reinterpret_cast(npP3.get_data()); + //P1_old = reinterpret_cast(npP1_old.get_data()); + //P2_old = reinterpret_cast(npP2_old.get_data()); + //P3_old = reinterpret_cast(npP3_old.get_data()); + //R1 = reinterpret_cast(npR1.get_data()); + //R2 = reinterpret_cast(npR2.get_data()); + //R3 = reinterpret_cast(npR3.get_data()); + ///* begin iterations */ + //for (ll = 0; ll 3) { + //Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ); + //funcval = 0.0f; + //for (j = 0; j 2) { + //if (re > re_old) { + //Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ); + //funcval = 0.0f; + //for (j = 0; j(npD); + //result.append(out1); + //result.append(ll); + //} - return result; -} + //return result; +//} bp::list LLT_model(np::ndarray input, double d_lambda, double d_tau, int iter, double d_epsil, int switcher) { // the result is in the following list diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx index c724471..263fa4a 100644 --- a/Wrappers/Python/src/gpu_regularizers.pyx +++ b/Wrappers/Python/src/gpu_regularizers.pyx @@ -25,7 +25,9 @@ cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec, int N, int M, int Z, int dimension, int SearchW, int SimilW, int SearchW_real, float denh2, float lambdaf); -cdef extern void TV_ROF_GPU_kernel(float* A, float* B, int N, int M, int Z, int iter, float tau, float lambdaf); +cdef extern void TV_ROF_GPU(float* Input, float* Output, int N, int M, int Z, int iter, float tau, float lambdaf); +cdef extern void TV_FGP_GPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z); + cdef extern float pad_crop(float *A, float *Ap, int OldSizeX, int OldSizeY, int OldSizeZ, int NewSizeX, int NewSizeY, int NewSizeZ, @@ -343,7 +345,7 @@ def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, np.zeros([dims[0],dims[1]], dtype='float32') # Running CUDA code here - TV_ROF_GPU_kernel( + TV_ROF_GPU( &inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations , @@ -366,7 +368,7 @@ def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, np.zeros([dims[0],dims[1],dims[2]], dtype='float32') # Running CUDA code here - TV_ROF_GPU_kernel( + TV_ROF_GPU( &inputData[0,0,0], &B[0,0,0], dims[0], dims[1], dims[2], iterations , @@ -374,3 +376,64 @@ def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, regularization_parameter); return B + + +def TVFGP2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularization_parameter, + int iterations, + float tolerance_param, + int methodTV, + int nonneg, + int printM): + + cdef long dims[2] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \ + np.zeros([dims[0],dims[1]], dtype='float32') + + # Running CUDA code here + TV_FGP_GPU( + &inputData[0,0], &B[0,0], + regularization_parameter , + iterations, + tolerance_param, + methodTV, + nonneg, + printM, + dims[0], dims[1], 1); + + return B + +def TVFGP3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularization_parameter, + int iterations, + float tolerance_param, + int methodTV, + int nonneg, + int printM): + + cdef long dims[3] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + dims[2] = inputData.shape[2] + + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \ + np.zeros([dims[0],dims[1],dims[2]], dtype='float32') + + # Running CUDA code here + TV_FGP_GPU( + &inputData[0,0,0], &B[0,0,0], + regularization_parameter , + iterations, + tolerance_param, + methodTV, + nonneg, + printM, + dims[0], dims[1], dims[2]); + + return B + + + diff --git a/Wrappers/Python/test/test_cpu_vs_gpu.py b/Wrappers/Python/test/test_cpu_vs_gpu.py new file mode 100644 index 0000000..74d65dd --- /dev/null +++ b/Wrappers/Python/test/test_cpu_vs_gpu.py @@ -0,0 +1,10 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Feb 21 12:12:22 2018 + +# CPU vs GPU comparison tests + +@author: algol +""" + diff --git a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py index d742a7f..6344021 100644 --- a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py +++ b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py @@ -12,8 +12,8 @@ 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 +from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, TV_ROF_GPU +from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU ############################################################################### def printParametersToString(pars): txt = r'' @@ -64,7 +64,7 @@ pars = {'algorithm': ROF_TV , \ } print ("#################ROF TV CPU#####################") start_time = timeit.default_timer() -rof_cpu = ROF_TV(pars['input'], +rof_cpu = TV_ROF_CPU(pars['input'], pars['number_of_iterations'], pars['regularization_parameter'], pars['time_marching_parameter'] @@ -89,7 +89,7 @@ plt.title('{}'.format('CPU results')) print ("#################ROF TV GPU#####################") start_time = timeit.default_timer() -rof_gpu = GPU_ROF_TV(pars['input'], +rof_gpu = TV_ROF_GPU(pars['input'], pars['number_of_iterations'], pars['time_marching_parameter'], pars['regularization_parameter']) -- cgit v1.2.3 From 8082a76d4dfd9588590bab3fd26eae976b744a94 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Mon, 26 Feb 2018 13:20:33 +0000 Subject: cpu-related cython file updated #33 --- Wrappers/Python/src/cpu_regularizers.pyx | 77 +++++++++++++++++++++++++++----- Wrappers/Python/src/gpu_regularizers.pyx | 2 +- 2 files changed, 68 insertions(+), 11 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx index e7ff78c..b8089a8 100644 --- a/Wrappers/Python/src/cpu_regularizers.pyx +++ b/Wrappers/Python/src/cpu_regularizers.pyx @@ -18,20 +18,20 @@ import cython import numpy as np cimport numpy as np -cdef extern float TV_ROF(float *A, float *B, int dimX, int dimY, int dimZ, - int iterationsNumb, float tau, float flambda); +cdef extern float TV_ROF_CPU(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float flambda); +cdef extern float TV_FGP_CPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ); - -def ROF_TV(inputData, iterations, regularization_parameter, +# Can we use the same name here in "def" as the C function? +def TV_ROF_CPU(inputData, iterations, regularization_parameter, marching_step_parameter): if inputData.ndim == 2: - return ROF_TV_2D(inputData, iterations, regularization_parameter, + return TV_ROF_2D(inputData, iterations, regularization_parameter, marching_step_parameter) elif inputData.ndim == 3: - return ROF_TV_3D(inputData, iterations, regularization_parameter, + return TV_ROF_3D(inputData, iterations, regularization_parameter, marching_step_parameter) -def ROF_TV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, +def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, int iterations, float regularization_parameter, float marching_step_parameter @@ -45,13 +45,13 @@ def ROF_TV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, np.zeros([dims[0],dims[1]], dtype='float32') #/* Run ROF iterations for 2D data */ - TV_ROF(&inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations, + TV_ROF_CPU(&inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations, marching_step_parameter, regularization_parameter) return B -def ROF_TV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, +def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, int iterations, float regularization_parameter, float marching_step_parameter @@ -65,8 +65,65 @@ def ROF_TV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, np.zeros([dims[0],dims[1],dims[2]], dtype='float32') #/* Run ROF iterations for 3D data */ - TV_ROF(&inputData[0,0, 0], &B[0,0, 0], dims[0], dims[1], dims[2], iterations, + TV_FGP_CPU(&inputData[0,0,0], &B[0,0,0], dims[0], dims[1], dims[2], iterations, marching_step_parameter, regularization_parameter) return B +def TV_FGP_CPU(inputData, regularization_parameter, iterations, tolerance_param, methodTV, nonneg, printM): + if inputData.ndim == 2: + return TV_FGP_2D(inputData, regularization_parameter, iterations, tolerance_param, methodTV, nonneg, printM) + elif inputData.ndim == 3: + return TV_FGP_3D(inputData, regularization_parameter, iterations, tolerance_param, methodTV, nonneg, printM) + +def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularization_parameter, + int iterations, + float tolerance_param, + int methodTV, + int nonneg, + int printM): + + cdef long dims[2] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \ + np.zeros([dims[0],dims[1]], dtype='float32') + + #/* Run ROF iterations for 2D data */ + TV_FGP_CPU(&inputData[0,0], &B[0,0], regularization_parameter, + iterations, + tolerance_param, + methodTV, + nonneg, + printM, + dims[0], dims[1], 1) + + return B + + +def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularization_parameter, + int iterations, + float tolerance_param, + int methodTV, + int nonneg, + int printM): + cdef long dims[2] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + dims[2] = inputData.shape[2] + + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \ + np.zeros([dims[0],dims[1],dims[2]], dtype='float32') + + #/* Run ROF iterations for 3D data */ + TV_FGP_CPU(&inputData[0,0, 0], &B[0,0, 0], iterations, + tolerance_param, + methodTV, + nonneg, + printM, + dims[0], dims[1], dims[2]) + + return B diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx index 263fa4a..a14a20d 100644 --- a/Wrappers/Python/src/gpu_regularizers.pyx +++ b/Wrappers/Python/src/gpu_regularizers.pyx @@ -396,7 +396,7 @@ def TVFGP2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, # Running CUDA code here TV_FGP_GPU( &inputData[0,0], &B[0,0], - regularization_parameter , + regularization_parameter, iterations, tolerance_param, methodTV, -- cgit v1.2.3 From 74ff5b5f319b2f7ea3e078c62041ec8a1bb28335 Mon Sep 17 00:00:00 2001 From: algol Date: Mon, 5 Mar 2018 18:12:01 +0000 Subject: Cmake/Cython fixes to compile ROF-FGP --- Wrappers/Python/demo/test_cpu_regularizers.py | 25 ++++++++++++---------- Wrappers/Python/src/cpu_regularizers.cpp | 1 - Wrappers/Python/src/cpu_regularizers.pyx | 25 +++++++++++----------- Wrappers/Python/src/gpu_regularizers.pyx | 2 +- .../Python/test/test_cpu_vs_gpu_regularizers.py | 9 ++++---- 5 files changed, 32 insertions(+), 30 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py index b08c418..53b8538 100644 --- a/Wrappers/Python/demo/test_cpu_regularizers.py +++ b/Wrappers/Python/demo/test_cpu_regularizers.py @@ -11,10 +11,9 @@ import numpy as np import os from enum import Enum import timeit -from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV , FGP_TV ,\ - LLT_model, PatchBased_Regul ,\ +from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV, LLT_model, PatchBased_Regul ,\ TGV_PD -from ccpi.filters.cpu_regularizers_cython import ROF_TV +from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU, TV_FGP_CPU ############################################################################### #https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956 @@ -128,21 +127,25 @@ imgplot = plt.imshow(splitbregman,\ ) ###################### FGP_TV ######################################### -# u = FGP_TV(single(u0), 0.05, 100, 1e-04); + start_time = timeit.default_timer() -pars = {'algorithm' : FGP_TV , \ +pars = {'algorithm' : TV_FGP_CPU , \ 'input' : u0, 'regularization_parameter':0.05, \ 'number_of_iterations' :200 ,\ - 'tolerance_constant':1e-4,\ - 'TV_penalty': 0 + 'tolerance_constant':1e-5,\ + 'methodTV': 0 ,\ + 'nonneg': 0 ,\ + 'printingOut': 0 } -out = FGP_TV (pars['input'], +out = TV_FGP_CPU (pars['input'], pars['regularization_parameter'], pars['number_of_iterations'], pars['tolerance_constant'], - pars['TV_penalty']) + pars['methodTV'], + pars['nonneg'], + pars['printingOut']) fgp = out[0] rms = rmse(Im, fgp) @@ -282,13 +285,13 @@ imgplot = plt.imshow(tgv, cmap="gray") start_time = timeit.default_timer() -pars = {'algorithm': ROF_TV , \ +pars = {'algorithm': TV_ROF_CPU , \ 'input' : u0,\ 'regularization_parameter':0.04,\ 'marching_step': 0.0025,\ 'number_of_iterations': 300 } -rof = ROF_TV(pars['input'], +rof = TV_ROF_CPU(pars['input'], pars['number_of_iterations'], pars['regularization_parameter'], pars['marching_step'] diff --git a/Wrappers/Python/src/cpu_regularizers.cpp b/Wrappers/Python/src/cpu_regularizers.cpp index 43d5d11..b8156d5 100644 --- a/Wrappers/Python/src/cpu_regularizers.cpp +++ b/Wrappers/Python/src/cpu_regularizers.cpp @@ -1040,7 +1040,6 @@ BOOST_PYTHON_MODULE(cpu_regularizers_boost) np::dtype dt2 = np::dtype::get_builtin(); def("SplitBregman_TV", SplitBregman_TV); - def("FGP_TV", FGP_TV); def("LLT_model", LLT_model); def("PatchBased_Regul", PatchBased_Regul); def("TGV_PD", TGV_PD); diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx index b8089a8..448da31 100644 --- a/Wrappers/Python/src/cpu_regularizers.pyx +++ b/Wrappers/Python/src/cpu_regularizers.pyx @@ -18,8 +18,8 @@ import cython import numpy as np cimport numpy as np -cdef extern float TV_ROF_CPU(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float flambda); -cdef extern float TV_FGP_CPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ); +cdef extern float TV_ROF_CPU_main(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float flambda); +cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ); # Can we use the same name here in "def" as the C function? def TV_ROF_CPU(inputData, iterations, regularization_parameter, @@ -45,11 +45,10 @@ def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, np.zeros([dims[0],dims[1]], dtype='float32') #/* Run ROF iterations for 2D data */ - TV_ROF_CPU(&inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations, + TV_ROF_CPU_main(&inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations, marching_step_parameter, regularization_parameter) - return B - + return B def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, int iterations, @@ -65,7 +64,7 @@ def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, np.zeros([dims[0],dims[1],dims[2]], dtype='float32') #/* Run ROF iterations for 3D data */ - TV_FGP_CPU(&inputData[0,0,0], &B[0,0,0], dims[0], dims[1], dims[2], iterations, + TV_ROF_CPU_main(&inputData[0,0,0], &B[0,0,0], dims[0], dims[1], dims[2], iterations, marching_step_parameter, regularization_parameter) return B @@ -88,11 +87,11 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] - cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \ + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \ np.zeros([dims[0],dims[1]], dtype='float32') #/* Run ROF iterations for 2D data */ - TV_FGP_CPU(&inputData[0,0], &B[0,0], regularization_parameter, + TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], regularization_parameter, iterations, tolerance_param, methodTV, @@ -100,8 +99,7 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, printM, dims[0], dims[1], 1) - return B - + return outputData def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularization_parameter, @@ -115,15 +113,16 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, dims[1] = inputData.shape[1] dims[2] = inputData.shape[2] - cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \ + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \ np.zeros([dims[0],dims[1],dims[2]], dtype='float32') #/* Run ROF iterations for 3D data */ - TV_FGP_CPU(&inputData[0,0, 0], &B[0,0, 0], iterations, + TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0, 0], regularization_parameter, + iterations, tolerance_param, methodTV, nonneg, printM, dims[0], dims[1], dims[2]) - return B + return outputData diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx index a14a20d..e99bfa7 100644 --- a/Wrappers/Python/src/gpu_regularizers.pyx +++ b/Wrappers/Python/src/gpu_regularizers.pyx @@ -26,7 +26,7 @@ cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec, int SearchW, int SimilW, int SearchW_real, float denh2, float lambdaf); cdef extern void TV_ROF_GPU(float* Input, float* Output, int N, int M, int Z, int iter, float tau, float lambdaf); -cdef extern void TV_FGP_GPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z); +cdef extern void TV_FGP_GPU(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z); cdef extern float pad_crop(float *A, float *Ap, int OldSizeX, int OldSizeY, int OldSizeZ, diff --git a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py index 6344021..e162afa 100644 --- a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py +++ b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py @@ -12,7 +12,7 @@ import matplotlib.pyplot as plt import numpy as np import os import timeit -from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, TV_ROF_GPU +from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, GPU_ROF_TV from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU ############################################################################### def printParametersToString(pars): @@ -56,11 +56,11 @@ imgplot = plt.imshow(u0,cmap="gray") # set parameters -pars = {'algorithm': ROF_TV , \ +pars = {'algorithm': TV_ROF_CPU , \ 'input' : u0,\ 'regularization_parameter':0.04,\ 'time_marching_parameter': 0.0025,\ - 'number_of_iterations': 600 + 'number_of_iterations': 1200 } print ("#################ROF TV CPU#####################") start_time = timeit.default_timer() @@ -89,13 +89,14 @@ plt.title('{}'.format('CPU results')) print ("#################ROF TV GPU#####################") start_time = timeit.default_timer() -rof_gpu = TV_ROF_GPU(pars['input'], +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 +pars['algorithm'] = GPU_ROF_TV txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) -- cgit v1.2.3 From ccf9b61bba1004af783c6333d58ea9611c0f81f2 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Tue, 6 Mar 2018 11:16:41 +0000 Subject: ROF_FGP_unified syntax --- Wrappers/Python/src/cpu_regularizers.pyx | 59 +++++++++++++++----------------- 1 file changed, 27 insertions(+), 32 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx index 448da31..2654831 100644 --- a/Wrappers/Python/src/cpu_regularizers.pyx +++ b/Wrappers/Python/src/cpu_regularizers.pyx @@ -11,73 +11,68 @@ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. -Author: Edoardo Pasca +Author: Edoardo Pasca, Daniil Kazantsev """ import cython import numpy as np cimport numpy as np -cdef extern float TV_ROF_CPU_main(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float flambda); -cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ); +cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ); +cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ); -# Can we use the same name here in "def" as the C function? -def TV_ROF_CPU(inputData, iterations, regularization_parameter, +def TV_ROF_CPU(inputData, regularization_parameter, iterationsNumb marching_step_parameter): if inputData.ndim == 2: - return TV_ROF_2D(inputData, iterations, regularization_parameter, + return TV_ROF_2D(inputData, regularization_parameter, iterationsNumb marching_step_parameter) elif inputData.ndim == 3: - return TV_ROF_3D(inputData, iterations, regularization_parameter, + return TV_ROF_3D(inputData, regularization_parameter, iterationsNumb marching_step_parameter) def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - int iterations, float regularization_parameter, - float marching_step_parameter - ): - + int iterationsNumb, + float marching_step_parameter): cdef long dims[2] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] - cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \ + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \ np.zeros([dims[0],dims[1]], dtype='float32') - #/* Run ROF iterations for 2D data */ - TV_ROF_CPU_main(&inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations, - marching_step_parameter, regularization_parameter) + # Run ROF iterations for 2D data + TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], regularization_parameter, iterationsNumb, marching_step_parameter, dims[0], dims[1], 1) - return B + return outputData def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, int iterations, float regularization_parameter, float marching_step_parameter ): - cdef long dims[2] + cdef long dims[3] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] dims[2] = inputData.shape[2] - cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \ + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \ np.zeros([dims[0],dims[1],dims[2]], dtype='float32') - #/* Run ROF iterations for 3D data */ - TV_ROF_CPU_main(&inputData[0,0,0], &B[0,0,0], dims[0], dims[1], dims[2], iterations, - marching_step_parameter, regularization_parameter) + # Run ROF iterations for 3D data + TV_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularization_parameter, iterationsNumb, marching_step_parameter, dims[0], dims[1], dims[2]) - return B + return outputData -def TV_FGP_CPU(inputData, regularization_parameter, iterations, tolerance_param, methodTV, nonneg, printM): +def TV_FGP_CPU(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM): if inputData.ndim == 2: - return TV_FGP_2D(inputData, regularization_parameter, iterations, tolerance_param, methodTV, nonneg, printM) + return TV_FGP_2D(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM) elif inputData.ndim == 3: - return TV_FGP_3D(inputData, regularization_parameter, iterations, tolerance_param, methodTV, nonneg, printM) + return TV_FGP_3D(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM) def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularization_parameter, - int iterations, + int iterationsNumb, float tolerance_param, int methodTV, int nonneg, @@ -92,7 +87,7 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, #/* Run ROF iterations for 2D data */ TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], regularization_parameter, - iterations, + iterationsNumb, tolerance_param, methodTV, nonneg, @@ -103,22 +98,22 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularization_parameter, - int iterations, + int iterationsNumb, float tolerance_param, int methodTV, int nonneg, int printM): - cdef long dims[2] + cdef long dims[3] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] dims[2] = inputData.shape[2] cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \ - np.zeros([dims[0],dims[1],dims[2]], dtype='float32') + np.zeros([dims[0], dims[1], dims[2]], dtype='float32') #/* Run ROF iterations for 3D data */ - TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0, 0], regularization_parameter, - iterations, + TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularization_parameter, + iterationsNumb, tolerance_param, methodTV, nonneg, -- cgit v1.2.3 From 8d310478254f3cda63f3663729b416f425ad70b6 Mon Sep 17 00:00:00 2001 From: algol Date: Tue, 6 Mar 2018 11:45:53 +0000 Subject: work on FGP intergration --- Wrappers/Python/demo/test_cpu_regularizers.py | 17 +++++++++-------- Wrappers/Python/src/cpu_regularizers.pyx | 14 +++++--------- 2 files changed, 14 insertions(+), 17 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py index 53b8538..f1eb3c3 100644 --- a/Wrappers/Python/demo/test_cpu_regularizers.py +++ b/Wrappers/Python/demo/test_cpu_regularizers.py @@ -131,9 +131,9 @@ imgplot = plt.imshow(splitbregman,\ start_time = timeit.default_timer() pars = {'algorithm' : TV_FGP_CPU , \ 'input' : u0, - 'regularization_parameter':0.05, \ - 'number_of_iterations' :200 ,\ - 'tolerance_constant':1e-5,\ + 'regularization_parameter':0.07, \ + 'number_of_iterations' :300 ,\ + 'tolerance_constant':0.00001,\ 'methodTV': 0 ,\ 'nonneg': 0 ,\ 'printingOut': 0 @@ -156,7 +156,7 @@ txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) -a=fig.add_subplot(2,4,3) +a=fig.add_subplot(2,4,4) # these are matplotlib.patch.Patch properties props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) @@ -168,8 +168,9 @@ imgplot = plt.imshow(fgp, \ a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) -###################### LLT_model ######################################### +###################### LLT_model ######################################### +""" start_time = timeit.default_timer() pars = {'algorithm': LLT_model , \ @@ -204,7 +205,7 @@ a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(llt,\ cmap="gray" ) - +""" # ###################### PatchBased_Regul ######################################### # # Quick 2D denoising example in Matlab: @@ -292,8 +293,8 @@ pars = {'algorithm': TV_ROF_CPU , \ 'number_of_iterations': 300 } rof = TV_ROF_CPU(pars['input'], - pars['number_of_iterations'], - pars['regularization_parameter'], + pars['regularization_parameter'], + pars['number_of_iterations'], pars['marching_step'] ) #tgv = out diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx index 2654831..d62ca59 100644 --- a/Wrappers/Python/src/cpu_regularizers.pyx +++ b/Wrappers/Python/src/cpu_regularizers.pyx @@ -21,14 +21,11 @@ cimport numpy as np cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ); cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ); -def TV_ROF_CPU(inputData, regularization_parameter, iterationsNumb - marching_step_parameter): +def TV_ROF_CPU(inputData, regularization_parameter, iterationsNumb, marching_step_parameter): if inputData.ndim == 2: - return TV_ROF_2D(inputData, regularization_parameter, iterationsNumb - marching_step_parameter) + return TV_ROF_2D(inputData, regularization_parameter, iterationsNumb, marching_step_parameter) elif inputData.ndim == 3: - return TV_ROF_3D(inputData, regularization_parameter, iterationsNumb - marching_step_parameter) + return TV_ROF_3D(inputData, regularization_parameter, iterationsNumb, marching_step_parameter) def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularization_parameter, @@ -47,10 +44,9 @@ def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, return outputData def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - int iterations, + int iterationsNumb, float regularization_parameter, - float marching_step_parameter - ): + float marching_step_parameter): cdef long dims[3] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] -- cgit v1.2.3 From 309d84445b5ec2980db7c79832958a6481343f17 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Tue, 6 Mar 2018 12:46:09 +0000 Subject: FGP/ROF updates --- .../Matlab/mex_compile/regularizers_CPU/FGP_TV.c | 142 ++------------------- Wrappers/Python/demo/test_cpu_regularizers.py | 10 +- Wrappers/Python/src/cpu_regularizers.pyx | 5 +- 3 files changed, 16 insertions(+), 141 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Matlab/mex_compile/regularizers_CPU/FGP_TV.c b/Wrappers/Matlab/mex_compile/regularizers_CPU/FGP_TV.c index 30cea1a..7cc861a 100644 --- a/Wrappers/Matlab/mex_compile/regularizers_CPU/FGP_TV.c +++ b/Wrappers/Matlab/mex_compile/regularizers_CPU/FGP_TV.c @@ -53,9 +53,9 @@ void mexFunction( int nrhs, const mxArray *prhs[]) { - int number_of_dims, iter, dimX, dimY, dimZ, ll, j, count, methTV; + int number_of_dims, iter, dimX, dimY, dimZ, methTV; const int *dim_array; - float *A, *D=NULL, *D_old=NULL, *P1=NULL, *P2=NULL, *P3=NULL, *P1_old=NULL, *P2_old=NULL, *P3_old=NULL, *R1=NULL, *R2=NULL, *R3=NULL, lambda, tk, tkp1, re, re1, re_old, epsil; + float *Input, *Output, lambda, epsil; number_of_dims = mxGetNumberOfDimensions(prhs[0]); dim_array = mxGetDimensions(prhs[0]); @@ -63,9 +63,9 @@ void mexFunction( /*Handling Matlab input data*/ if ((nrhs < 2) || (nrhs > 5)) mexErrMsgTxt("At least 2 parameters is required: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1')"); - A = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */ + Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */ lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */ - iter = 50; /* default iterations number */ + iter = 300; /* default iterations number */ epsil = 0.0001; /* default tolerance constant */ methTV = 0; /* default isotropic TV penalty */ @@ -78,139 +78,17 @@ void mexFunction( if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */ mxFree(penalty_type); } - /*output function value (last iteration) */ - plhs[1] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL); - float *funcvalA = (float *) mxGetData(plhs[1]); if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } /*Handling Matlab output data*/ - dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; - - tk = 1.0f; - tkp1=1.0f; - count = 0; - re_old = 0.0f; + dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; if (number_of_dims == 2) { dimZ = 1; /*2D case*/ - D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - D_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - P1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - P2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - P1_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - P2_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - R1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - R2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - - /* begin iterations */ - for(ll=0; ll 4) { - Obj_func_CALC2D(A, D, funcvalA, lambda, dimX, dimY); - break; } - - /* check that the residual norm is decreasing */ - if (ll > 2) { - if (re > re_old) { - Obj_func_CALC2D(A, D, funcvalA, lambda, dimX, dimY); - break; }} - re_old = re; - /*printf("%f %i %i \n", re, ll, count); */ - - /*storing old values*/ - copyIm(D, D_old, dimX, dimY, dimZ); - copyIm(P1, P1_old, dimX, dimY, dimZ); - copyIm(P2, P2_old, dimX, dimY, dimZ); - tk = tkp1; - - /* calculating the objective function value */ - if (ll == (iter-1)) Obj_func_CALC2D(A, D, funcvalA, lambda, dimX, dimY); - } - printf("FGP-TV iterations stopped at iteration %i with the function value %f \n", ll, funcvalA[0]); - } - if (number_of_dims == 3) { - D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - D_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P1_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P2_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P3_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - R1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - R2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - R3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - - /* begin iterations */ - for(ll=0; ll 3) { - Obj_func_CALC3D(A, D, funcvalA, lambda, dimX, dimY, dimZ); - break;} - - /* check that the residual norm is decreasing */ - if (ll > 2) { - if (re > re_old) { - Obj_func_CALC3D(A, D, funcvalA, lambda, dimX, dimY, dimZ); - }} - re_old = re; - /*printf("%f %i %i \n", re, ll, count); */ - - /*storing old values*/ - copyIm(D, D_old, dimX, dimY, dimZ); - copyIm(P1, P1_old, dimX, dimY, dimZ); - copyIm(P2, P2_old, dimX, dimY, dimZ); - copyIm(P3, P3_old, dimX, dimY, dimZ); - tk = tkp1; - - if (ll == (iter-1)) Obj_func_CALC3D(A, D, funcvalA, lambda, dimX, dimY, dimZ); - } - printf("FGP-TV iterations stopped at iteration %i with the function value %f \n", ll, funcvalA[0]); - } + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + TV_FGP_CPU_main(Input, Output, lambda, iter, epsil, methTV, 0, 0, dimX, dimY, dimZ); + } + if (number_of_dims == 3) { + } } diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py index f1eb3c3..7f08605 100644 --- a/Wrappers/Python/demo/test_cpu_regularizers.py +++ b/Wrappers/Python/demo/test_cpu_regularizers.py @@ -111,12 +111,10 @@ rms = rmse(Im, splitbregman) pars['rmse'] = rms txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) - +print (txtstr) a=fig.add_subplot(2,4,2) - # these are matplotlib.patch.Patch properties props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) # place a text box in upper left in axes coords @@ -130,14 +128,14 @@ imgplot = plt.imshow(splitbregman,\ start_time = timeit.default_timer() pars = {'algorithm' : TV_FGP_CPU , \ - 'input' : u0, + 'input' : u0,\ 'regularization_parameter':0.07, \ 'number_of_iterations' :300 ,\ 'tolerance_constant':0.00001,\ 'methodTV': 0 ,\ 'nonneg': 0 ,\ - 'printingOut': 0 -} + 'printingOut': 0 + } out = TV_FGP_CPU (pars['input'], pars['regularization_parameter'], diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx index d62ca59..60e8627 100644 --- a/Wrappers/Python/src/cpu_regularizers.pyx +++ b/Wrappers/Python/src/cpu_regularizers.pyx @@ -93,7 +93,7 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, return outputData def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - float regularization_parameter, + float regularization_parameter, int iterationsNumb, float tolerance_param, int methodTV, @@ -109,11 +109,10 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, #/* Run ROF iterations for 3D data */ TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularization_parameter, - iterationsNumb, + iterationsNumb, tolerance_param, methodTV, nonneg, printM, dims[0], dims[1], dims[2]) - return outputData -- cgit v1.2.3 From 69ecdd57434d591eb3fa4afefb72174d3e025fb9 Mon Sep 17 00:00:00 2001 From: algol Date: Tue, 6 Mar 2018 12:48:21 +0000 Subject: FGP_CPU (Cythonized) now works in demo --- Wrappers/Python/demo/test_cpu_regularizers.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py index 7f08605..1d97857 100644 --- a/Wrappers/Python/demo/test_cpu_regularizers.py +++ b/Wrappers/Python/demo/test_cpu_regularizers.py @@ -137,7 +137,7 @@ pars = {'algorithm' : TV_FGP_CPU , \ 'printingOut': 0 } -out = TV_FGP_CPU (pars['input'], +fgp = TV_FGP_CPU(pars['input'], pars['regularization_parameter'], pars['number_of_iterations'], pars['tolerance_constant'], @@ -145,7 +145,6 @@ out = TV_FGP_CPU (pars['input'], pars['nonneg'], pars['printingOut']) -fgp = out[0] rms = rmse(Im, fgp) pars['rmse'] = rms @@ -154,7 +153,7 @@ txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) -a=fig.add_subplot(2,4,4) +a=fig.add_subplot(2,4,3) # these are matplotlib.patch.Patch properties props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) @@ -168,7 +167,6 @@ a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, ###################### LLT_model ######################################### -""" start_time = timeit.default_timer() pars = {'algorithm': LLT_model , \ @@ -203,8 +201,6 @@ a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(llt,\ cmap="gray" ) -""" - # ###################### PatchBased_Regul ######################################### # # Quick 2D denoising example in Matlab: # # Im = double(imread('lena_gray_256.tif'))/255; % loading image @@ -286,7 +282,7 @@ start_time = timeit.default_timer() pars = {'algorithm': TV_ROF_CPU , \ 'input' : u0,\ - 'regularization_parameter':0.04,\ + 'regularization_parameter':0.07,\ 'marching_step': 0.0025,\ 'number_of_iterations': 300 } -- cgit v1.2.3 From 5411ebbd4165c81b398b010d6ad9d11d2e973aad Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Tue, 6 Mar 2018 14:40:11 +0000 Subject: work on cythonization2 --- Wrappers/Python/src/gpu_regularizers.pyx | 97 ++++++++++++++++++-------------- 1 file changed, 56 insertions(+), 41 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx index e99bfa7..cb94e86 100644 --- a/Wrappers/Python/src/gpu_regularizers.pyx +++ b/Wrappers/Python/src/gpu_regularizers.pyx @@ -11,7 +11,7 @@ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. -Author: Edoardo Pasca +Author: Edoardo Pasca, Daniil Kazantsev """ import cython @@ -25,14 +25,16 @@ cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec, int N, int M, int Z, int dimension, int SearchW, int SimilW, int SearchW_real, float denh2, float lambdaf); -cdef extern void TV_ROF_GPU(float* Input, float* Output, int N, int M, int Z, int iter, float tau, float lambdaf); -cdef extern void TV_FGP_GPU(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z); +cdef extern void TV_ROF_GPU_main(float* Input, float* Output, float lambdaPar, int iter, float tau, int N, int M, int Z); +#cdef extern void TV_FGP_GPU(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z); +# correct the function +cdef extern void TV_FGP_GPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z); cdef extern float pad_crop(float *A, float *Ap, int OldSizeX, int OldSizeY, int OldSizeZ, int NewSizeX, int NewSizeY, int NewSizeZ, int padXY, int switchpad_crop); - +#Diffusion 4th order regularizer def Diff4thHajiaboli(inputData, edge_preserv_parameter, iterations, @@ -50,7 +52,7 @@ def Diff4thHajiaboli(inputData, iterations, time_marching_parameter, regularization_parameter) - +# patch-based nonlocal regularization def NML(inputData, SearchW_real, SimilW, @@ -68,23 +70,37 @@ def NML(inputData, SimilW, h, lambdaf) - -def GPU_ROF_TV(inputData, +# Total-variation Rudin-Osher-Fatemi (ROF) +def TV_ROF_GPU(inputData, + regularization_parameter, iterations, - time_marching_parameter, - regularization_parameter): + time_marching_parameter): if inputData.ndim == 2: return ROFTV2D(inputData, - iterations, - time_marching_parameter, - regularization_parameter) + regularization_parameter, + iterations, + time_marching_parameter) elif inputData.ndim == 3: return ROFTV3D(inputData, + regularization_parameter iterations, - time_marching_parameter, - regularization_parameter) - - + time_marching_parameter) +# Total-variation Fast-Gradient-Projection (FGP) +def TV_FGP_GPU(inputData, + regularization_parameter, + iterations, + time_marching_parameter): + if inputData.ndim == 2: + return FGPTV2D(inputData, + regularization_parameter, + iterations, + time_marching_parameter) + elif inputData.ndim == 3: + return FGPTV3D(inputData, + regularization_parameter + iterations, + time_marching_parameter) +#****************************************************************# def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float edge_preserv_parameter, int iterations, @@ -333,52 +349,51 @@ def NML3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, return B def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularization_parameter, int iterations, - float time_marching_parameter, - float regularization_parameter): + float time_marching_parameter): cdef long dims[2] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] - cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \ + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \ np.zeros([dims[0],dims[1]], dtype='float32') # Running CUDA code here - TV_ROF_GPU( - &inputData[0,0], &B[0,0], - dims[0], dims[1], 1, + TV_ROF_GPU_main( + &inputData[0,0], &outputData[0,0], + regularization_parameter, iterations , time_marching_parameter, - regularization_parameter); + dims[0], dims[1], 1); - return B + return outputData def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularization_parameter, int iterations, - float time_marching_parameter, - float regularization_parameter): + float time_marching_parameter): cdef long dims[3] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] dims[2] = inputData.shape[2] - cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \ + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \ np.zeros([dims[0],dims[1],dims[2]], dtype='float32') # Running CUDA code here - TV_ROF_GPU( - &inputData[0,0,0], &B[0,0,0], - dims[0], dims[1], dims[2], + TV_ROF_GPU_main( + &inputData[0,0,0], &outputData[0,0,0], + regularization_parameter, iterations , time_marching_parameter, - regularization_parameter); + dims[0], dims[1], dims[2]); - return B - + return outputData -def TVFGP2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, +def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularization_parameter, int iterations, float tolerance_param, @@ -390,12 +405,12 @@ def TVFGP2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] - cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \ + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \ np.zeros([dims[0],dims[1]], dtype='float32') # Running CUDA code here TV_FGP_GPU( - &inputData[0,0], &B[0,0], + &inputData[0,0], &outputData[0,0], regularization_parameter, iterations, tolerance_param, @@ -404,9 +419,9 @@ def TVFGP2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, printM, dims[0], dims[1], 1); - return B + return outputData -def TVFGP3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, +def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularization_parameter, int iterations, float tolerance_param, @@ -419,12 +434,12 @@ def TVFGP3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, dims[1] = inputData.shape[1] dims[2] = inputData.shape[2] - cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \ + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \ np.zeros([dims[0],dims[1],dims[2]], dtype='float32') # Running CUDA code here TV_FGP_GPU( - &inputData[0,0,0], &B[0,0,0], + &inputData[0,0,0], &outputData[0,0,0], regularization_parameter , iterations, tolerance_param, @@ -433,7 +448,7 @@ def TVFGP3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, printM, dims[0], dims[1], dims[2]); - return B + return outputData -- cgit v1.2.3 From 63ad8306c261a90c572d95084bf1dd8db2b3dce7 Mon Sep 17 00:00:00 2001 From: algol Date: Tue, 6 Mar 2018 16:18:20 +0000 Subject: FGP/ROF-CPU/GPU fully working with new interfaces and regularizers.py script. This closes #39 and #38 --- Wrappers/Python/ccpi/filters/regularizers.py | 44 +++++++ Wrappers/Python/src/gpu_regularizers.pyx | 43 ++++--- Wrappers/Python/test/test_cpu_vs_gpu.py | 10 -- .../Python/test/test_cpu_vs_gpu_regularizers.py | 137 ++++++++++++++++++--- Wrappers/Python/test/test_gpu_regularizers.py | 16 +-- 5 files changed, 196 insertions(+), 54 deletions(-) create mode 100644 Wrappers/Python/ccpi/filters/regularizers.py delete mode 100644 Wrappers/Python/test/test_cpu_vs_gpu.py (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/filters/regularizers.py b/Wrappers/Python/ccpi/filters/regularizers.py new file mode 100644 index 0000000..d6dfa8c --- /dev/null +++ b/Wrappers/Python/ccpi/filters/regularizers.py @@ -0,0 +1,44 @@ +""" +script which assigns a proper device core function based on a flag ('cpu' or 'gpu') +""" + +from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU, TV_FGP_CPU +from ccpi.filters.gpu_regularizers import TV_ROF_GPU, TV_FGP_GPU + +def ROF_TV(inputData, regularization_parameter, iterations, + time_marching_parameter,device='cpu'): + if device == 'cpu': + return TV_ROF_CPU(inputData, + regularization_parameter, + iterations, + time_marching_parameter) + elif device == 'gpu': + return TV_ROF_GPU(inputData, + regularization_parameter, + iterations, + time_marching_parameter) + else: + raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ + .format(device)) + +def FGP_TV(inputData, regularization_parameter,iterations, + tolerance_param, methodTV, nonneg, printM, device='cpu'): + if device == 'cpu': + return TV_FGP_CPU(inputData, + regularization_parameter, + iterations, + tolerance_param, + methodTV, + nonneg, + printM) + elif device == 'gpu': + return TV_FGP_GPU(inputData, + regularization_parameter, + iterations, + tolerance_param, + methodTV, + nonneg, + printM) + else: + raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ + .format(device)) \ No newline at end of file diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx index cb94e86..f96156a 100644 --- a/Wrappers/Python/src/gpu_regularizers.pyx +++ b/Wrappers/Python/src/gpu_regularizers.pyx @@ -26,14 +26,14 @@ cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec, int SearchW, int SimilW, int SearchW_real, float denh2, float lambdaf); cdef extern void TV_ROF_GPU_main(float* Input, float* Output, float lambdaPar, int iter, float tau, int N, int M, int Z); -#cdef extern void TV_FGP_GPU(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z); -# correct the function cdef extern void TV_FGP_GPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z); +# padding function cdef extern float pad_crop(float *A, float *Ap, int OldSizeX, int OldSizeY, int OldSizeZ, int NewSizeX, int NewSizeY, int NewSizeZ, int padXY, int switchpad_crop); + #Diffusion 4th order regularizer def Diff4thHajiaboli(inputData, edge_preserv_parameter, @@ -70,6 +70,7 @@ def NML(inputData, SimilW, h, lambdaf) + # Total-variation Rudin-Osher-Fatemi (ROF) def TV_ROF_GPU(inputData, regularization_parameter, @@ -82,24 +83,34 @@ def TV_ROF_GPU(inputData, time_marching_parameter) elif inputData.ndim == 3: return ROFTV3D(inputData, - regularization_parameter + regularization_parameter, iterations, time_marching_parameter) + # Total-variation Fast-Gradient-Projection (FGP) def TV_FGP_GPU(inputData, regularization_parameter, iterations, - time_marching_parameter): + tolerance_param, + methodTV, + nonneg, + printM): if inputData.ndim == 2: - return FGPTV2D(inputData, + return FGPTV2D(inputData, regularization_parameter, - iterations, - time_marching_parameter) + iterations, + tolerance_param, + methodTV, + nonneg, + printM) elif inputData.ndim == 3: - return FGPTV3D(inputData, - regularization_parameter + return FGPTV3D(inputData, + regularization_parameter, iterations, - time_marching_parameter) + tolerance_param, + methodTV, + nonneg, + printM) #****************************************************************# def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float edge_preserv_parameter, @@ -347,7 +358,8 @@ def NML3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, switchpad_crop) return B - + +# Total-variation ROF def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularization_parameter, int iterations, @@ -393,6 +405,7 @@ def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, return outputData +# Total-variation Fast-Gradient-Projection (FGP) def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularization_parameter, int iterations, @@ -409,8 +422,7 @@ def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, np.zeros([dims[0],dims[1]], dtype='float32') # Running CUDA code here - TV_FGP_GPU( - &inputData[0,0], &outputData[0,0], + TV_FGP_GPU_main(&inputData[0,0], &outputData[0,0], regularization_parameter, iterations, tolerance_param, @@ -438,7 +450,7 @@ def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, np.zeros([dims[0],dims[1],dims[2]], dtype='float32') # Running CUDA code here - TV_FGP_GPU( + TV_FGP_GPU_main( &inputData[0,0,0], &outputData[0,0,0], regularization_parameter , iterations, @@ -449,6 +461,3 @@ def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, dims[0], dims[1], dims[2]); return outputData - - - diff --git a/Wrappers/Python/test/test_cpu_vs_gpu.py b/Wrappers/Python/test/test_cpu_vs_gpu.py deleted file mode 100644 index 74d65dd..0000000 --- a/Wrappers/Python/test/test_cpu_vs_gpu.py +++ /dev/null @@ -1,10 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Wed Feb 21 12:12:22 2018 - -# CPU vs GPU comparison tests - -@author: algol -""" - diff --git a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py index e162afa..a9d0f31 100644 --- a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py +++ b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py @@ -5,15 +5,17 @@ Created on Thu Feb 22 11:39:43 2018 Testing CPU implementation against GPU one -@author: algol +@author: Daniil Kazantsev """ import matplotlib.pyplot as plt import numpy as np -import os +import os import timeit -from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, GPU_ROF_TV -from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU +#from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU, TV_FGP_CPU +#from ccpi.filters.gpu_regularizers import TV_ROF_GPU, TV_FGP_GPU +from ccpi.filters.regularizers import ROF_TV, FGP_TV + ############################################################################### def printParametersToString(pars): txt = r'' @@ -47,6 +49,11 @@ u0 = Im + np.random.normal(loc = Im , f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1) u0 = f(u0).astype('float32') + +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("____________ROF-TV bench___________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + ## plot fig = plt.figure(1) plt.suptitle('Comparison of ROF-TV regularizer using CPU and GPU implementations') @@ -54,22 +61,19 @@ a=fig.add_subplot(1,4,1) a.set_title('Noisy Image') imgplot = plt.imshow(u0,cmap="gray") - # set parameters -pars = {'algorithm': TV_ROF_CPU , \ +pars = {'algorithm': ROF_TV, \ 'input' : u0,\ 'regularization_parameter':0.04,\ - 'time_marching_parameter': 0.0025,\ - 'number_of_iterations': 1200 + 'number_of_iterations': 1200,\ + 'time_marching_parameter': 0.0025 } -print ("#################ROF TV CPU#####################") +print ("#############ROF TV CPU####################") start_time = timeit.default_timer() -rof_cpu = TV_ROF_CPU(pars['input'], - pars['number_of_iterations'], +rof_cpu = ROF_TV(pars['input'], pars['regularization_parameter'], - pars['time_marching_parameter'] - ) -#tgv = out + pars['number_of_iterations'], + pars['time_marching_parameter'],'cpu') rms = rmse(Im, rof_cpu) pars['rmse'] = rms @@ -87,16 +91,16 @@ imgplot = plt.imshow(rof_cpu, cmap="gray") plt.title('{}'.format('CPU results')) -print ("#################ROF TV GPU#####################") +print ("##############ROF TV GPU##################") start_time = timeit.default_timer() -rof_gpu = GPU_ROF_TV(pars['input'], +rof_gpu = ROF_TV(pars['input'], + pars['regularization_parameter'], pars['number_of_iterations'], - pars['time_marching_parameter'], - pars['regularization_parameter']) + pars['time_marching_parameter'],'gpu') rms = rmse(Im, rof_gpu) pars['rmse'] = rms -pars['algorithm'] = GPU_ROF_TV +pars['algorithm'] = ROF_TV txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) @@ -112,7 +116,8 @@ plt.title('{}'.format('GPU results')) print ("--------Compare the results--------") -tolerance = 1e-06 +tolerance = 1e-05 +diff_im = np.zeros(np.shape(rof_cpu)) diff_im = abs(rof_cpu - rof_gpu) diff_im[diff_im > tolerance] = 1 a=fig.add_subplot(1,4,4) @@ -122,3 +127,95 @@ if (diff_im.sum() > 1): print ("Arrays do not match!") else: print ("Arrays match") + +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("____________FGP-TV bench___________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure(2) +plt.suptitle('Comparison of FGP-TV regularizer using CPU and GPU implementations') +a=fig.add_subplot(1,4,1) +a.set_title('Noisy Image') +imgplot = plt.imshow(u0,cmap="gray") + +# set parameters +pars = {'algorithm' : FGP_TV, \ + 'input' : u0,\ + 'regularization_parameter':0.04, \ + 'number_of_iterations' :1200 ,\ + 'tolerance_constant':0.00001,\ + 'methodTV': 0 ,\ + 'nonneg': 0 ,\ + 'printingOut': 0 + } + +print ("#############FGP TV CPU####################") +start_time = timeit.default_timer() +fgp_cpu = FGP_TV(pars['input'], + pars['regularization_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['nonneg'], + pars['printingOut'],'cpu') + + +rms = rmse(Im, fgp_cpu) +pars['rmse'] = rms + +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(1,4,2) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) +# place a text box in upper left in axes coords +a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(fgp_cpu, cmap="gray") +plt.title('{}'.format('CPU results')) + + +print ("##############FGP TV GPU##################") +start_time = timeit.default_timer() +fgp_gpu = FGP_TV(pars['input'], + pars['regularization_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['nonneg'], + pars['printingOut'],'gpu') + +rms = rmse(Im, fgp_gpu) +pars['rmse'] = rms +pars['algorithm'] = FGP_TV +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(1,4,3) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) +# place a text box in upper left in axes coords +a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(fgp_gpu, cmap="gray") +plt.title('{}'.format('GPU results')) + + +print ("--------Compare the results--------") +tolerance = 1e-05 +diff_im = np.zeros(np.shape(rof_cpu)) +diff_im = abs(fgp_cpu - fgp_gpu) +diff_im[diff_im > tolerance] = 1 +a=fig.add_subplot(1,4,4) +imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray") +plt.title('{}'.format('Pixels larger threshold difference')) +if (diff_im.sum() > 1): + print ("Arrays do not match!") +else: + print ("Arrays match") + + diff --git a/Wrappers/Python/test/test_gpu_regularizers.py b/Wrappers/Python/test/test_gpu_regularizers.py index c473270..04aeeb4 100644 --- a/Wrappers/Python/test/test_gpu_regularizers.py +++ b/Wrappers/Python/test/test_gpu_regularizers.py @@ -11,7 +11,8 @@ import numpy as np import os from enum import Enum import timeit -from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, GPU_ROF_TV +from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, TV_ROF_GPU + ############################################################################### def printParametersToString(pars): txt = r'' @@ -152,17 +153,18 @@ plt.colorbar(ticks=[0, 0.03], orientation='vertical') start_time = timeit.default_timer() pars = { -'algorithm' : GPU_ROF_TV , \ +'algorithm' : TV_ROF_GPU , \ 'input' : u0, 'regularization_parameter': 0.04,\ - 'time_marching_parameter': 0.0025, \ - 'number_of_iterations':300 + 'number_of_iterations':300,\ + 'time_marching_parameter': 0.0025 + } -rof_tv = GPU_ROF_TV(pars['input'], +rof_tv = TV_ROF_GPU(pars['input'], + pars['regularization_parameter'], pars['number_of_iterations'], - pars['time_marching_parameter'], - pars['regularization_parameter']) + pars['time_marching_parameter']) rms = rmse(Im, rof_tv) pars['rmse'] = rms -- cgit v1.2.3 From 4593040d228fec5ef1f1a301e511708d47f0d55e Mon Sep 17 00:00:00 2001 From: algol Date: Tue, 6 Mar 2018 16:51:31 +0000 Subject: demos fixed #40 --- Wrappers/Python/demo/test_cpu_regularizers.py | 21 +- Wrappers/Python/src/cpu_regularizers.cpp | 289 --------------------- .../Python/test/test_cpu_vs_gpu_regularizers.py | 4 +- Wrappers/Python/test/test_gpu_regularizers.py | 59 ++++- 4 files changed, 63 insertions(+), 310 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py index 1d97857..76b9de7 100644 --- a/Wrappers/Python/demo/test_cpu_regularizers.py +++ b/Wrappers/Python/demo/test_cpu_regularizers.py @@ -11,10 +11,8 @@ import numpy as np import os from enum import Enum import timeit -from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV, LLT_model, PatchBased_Regul ,\ - TGV_PD -from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU, TV_FGP_CPU - +from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV, LLT_model, PatchBased_Regul, TGV_PD +from ccpi.filters.regularizers import ROF_TV, FGP_TV ############################################################################### #https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956 #NRMSE a normalization of the root of the mean squared error @@ -127,7 +125,7 @@ imgplot = plt.imshow(splitbregman,\ ###################### FGP_TV ######################################### start_time = timeit.default_timer() -pars = {'algorithm' : TV_FGP_CPU , \ +pars = {'algorithm' : FGP_TV , \ 'input' : u0,\ 'regularization_parameter':0.07, \ 'number_of_iterations' :300 ,\ @@ -137,13 +135,13 @@ pars = {'algorithm' : TV_FGP_CPU , \ 'printingOut': 0 } -fgp = TV_FGP_CPU(pars['input'], +fgp = FGP_TV(pars['input'], pars['regularization_parameter'], pars['number_of_iterations'], pars['tolerance_constant'], pars['methodTV'], pars['nonneg'], - pars['printingOut']) + pars['printingOut'], 'cpu') rms = rmse(Im, fgp) pars['rmse'] = rms @@ -280,18 +278,17 @@ imgplot = plt.imshow(tgv, cmap="gray") start_time = timeit.default_timer() -pars = {'algorithm': TV_ROF_CPU , \ +pars = {'algorithm': ROF_TV , \ 'input' : u0,\ 'regularization_parameter':0.07,\ 'marching_step': 0.0025,\ 'number_of_iterations': 300 } -rof = TV_ROF_CPU(pars['input'], +rof = ROF_TV(pars['input'], pars['regularization_parameter'], pars['number_of_iterations'], - pars['marching_step'] - ) -#tgv = out + pars['marching_step'], 'cpu') + rms = rmse(Im, rof) pars['rmse'] = rms diff --git a/Wrappers/Python/src/cpu_regularizers.cpp b/Wrappers/Python/src/cpu_regularizers.cpp index b8156d5..88b7c3c 100644 --- a/Wrappers/Python/src/cpu_regularizers.cpp +++ b/Wrappers/Python/src/cpu_regularizers.cpp @@ -303,292 +303,6 @@ bp::list SplitBregman_TV(np::ndarray input, double d_mu, int iter, double d_epsi } - - -//bp::list FGP_TV(np::ndarray input, double d_mu, int iter, double d_epsil, int methTV) { - - //// the result is in the following list - //bp::list result; - - //int number_of_dims, dimX, dimY, dimZ, ll, j, count; - //float *A, *D = NULL, *D_old = NULL, *P1 = NULL, *P2 = NULL, *P3 = NULL, *P1_old = NULL, *P2_old = NULL, *P3_old = NULL, *R1 = NULL, *R2 = NULL, *R3 = NULL; - //float lambda, tk, tkp1, re, re1, re_old, epsil, funcval; - - ////number_of_dims = mxGetNumberOfDimensions(prhs[0]); - ////dim_array = mxGetDimensions(prhs[0]); - - //number_of_dims = input.get_nd(); - //int dim_array[3]; - - //dim_array[0] = input.shape(0); - //dim_array[1] = input.shape(1); - //if (number_of_dims == 2) { - //dim_array[2] = -1; - //} - //else { - //dim_array[2] = input.shape(2); - //} - //// Parameter handling is be done in Python - /////*Handling Matlab input data*/ - ////if ((nrhs < 2) || (nrhs > 5)) mexErrMsgTxt("At least 2 parameters is required: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1')"); - - /////*Handling Matlab input data*/ - ////A = (float *)mxGetData(prhs[0]); /*noisy image (2D/3D) */ - //A = reinterpret_cast(input.get_data()); - - ////mu = (float)mxGetScalar(prhs[1]); /* regularization parameter */ - //lambda = (float)d_mu; - - ////iter = 35; /* default iterations number */ - - ////epsil = 0.0001; /* default tolerance constant */ - //epsil = (float)d_epsil; - ////methTV = 0; /* default isotropic TV penalty */ - ////if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5)) iter = (int)mxGetScalar(prhs[2]); /* iterations number */ - ////if ((nrhs == 4) || (nrhs == 5)) epsil = (float)mxGetScalar(prhs[3]); /* tolerance constant */ - ////if (nrhs == 5) { - //// char *penalty_type; - //// penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */ - //// if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',"); - //// if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */ - //// mxFree(penalty_type); - ////} - ////if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input image must be in a single precision"); } - - ////plhs[1] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL); - //bp::tuple shape1 = bp::make_tuple(dim_array[0], dim_array[1]); - //np::dtype dtype = np::dtype::get_builtin(); - //np::ndarray out1 = np::zeros(shape1, dtype); - - ////float *funcvalA = (float *)mxGetData(plhs[1]); - //float * funcvalA = reinterpret_cast(out1.get_data()); - ////if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input image must be in a single precision"); } - - ///*Handling Matlab output data*/ - //dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; - - //tk = 1.0f; - //tkp1 = 1.0f; - //count = 1; - //re_old = 0.0f; - - //if (number_of_dims == 2) { - //dimZ = 1; /*2D case*/ - ///*D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - //D_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - //P1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - //P2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - //P1_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - //P2_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - //R1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - //R2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));*/ - - //bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1]); - //np::dtype dtype = np::dtype::get_builtin(); - - - //np::ndarray npD = np::zeros(shape, dtype); - //np::ndarray npD_old = np::zeros(shape, dtype); - //np::ndarray npP1 = np::zeros(shape, dtype); - //np::ndarray npP2 = np::zeros(shape, dtype); - //np::ndarray npP1_old = np::zeros(shape, dtype); - //np::ndarray npP2_old = np::zeros(shape, dtype); - //np::ndarray npR1 = np::zeros(shape, dtype); - //np::ndarray npR2 = np::zeros(shape, dtype); - - //D = reinterpret_cast(npD.get_data()); - //D_old = reinterpret_cast(npD_old.get_data()); - //P1 = reinterpret_cast(npP1.get_data()); - //P2 = reinterpret_cast(npP2.get_data()); - //P1_old = reinterpret_cast(npP1_old.get_data()); - //P2_old = reinterpret_cast(npP2_old.get_data()); - //R1 = reinterpret_cast(npR1.get_data()); - //R2 = reinterpret_cast(npR2.get_data()); - - ///* begin iterations */ - //for (ll = 0; ll 4) { - //Obj_func2D(A, D, P1, P2, lambda, dimX, dimY); - //funcval = 0.0f; - //for (j = 0; j 2) { - //if (re > re_old) { - //Obj_func2D(A, D, P1, P2, lambda, dimX, dimY); - //funcval = 0.0f; - //for (j = 0; j(npD); - //result.append(out1); - //result.append(ll); - //} - //if (number_of_dims == 3) { - ///*D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - //D_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - //P1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - //P2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - //P3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - //P1_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - //P2_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - //P3_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - //R1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - //R2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - //R3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));*/ - //bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1], dim_array[2]); - //np::dtype dtype = np::dtype::get_builtin(); - - //np::ndarray npD = np::zeros(shape, dtype); - //np::ndarray npD_old = np::zeros(shape, dtype); - //np::ndarray npP1 = np::zeros(shape, dtype); - //np::ndarray npP2 = np::zeros(shape, dtype); - //np::ndarray npP3 = np::zeros(shape, dtype); - //np::ndarray npP1_old = np::zeros(shape, dtype); - //np::ndarray npP2_old = np::zeros(shape, dtype); - //np::ndarray npP3_old = np::zeros(shape, dtype); - //np::ndarray npR1 = np::zeros(shape, dtype); - //np::ndarray npR2 = np::zeros(shape, dtype); - //np::ndarray npR3 = np::zeros(shape, dtype); - - //D = reinterpret_cast(npD.get_data()); - //D_old = reinterpret_cast(npD_old.get_data()); - //P1 = reinterpret_cast(npP1.get_data()); - //P2 = reinterpret_cast(npP2.get_data()); - //P3 = reinterpret_cast(npP3.get_data()); - //P1_old = reinterpret_cast(npP1_old.get_data()); - //P2_old = reinterpret_cast(npP2_old.get_data()); - //P3_old = reinterpret_cast(npP3_old.get_data()); - //R1 = reinterpret_cast(npR1.get_data()); - //R2 = reinterpret_cast(npR2.get_data()); - //R3 = reinterpret_cast(npR3.get_data()); - ///* begin iterations */ - //for (ll = 0; ll 3) { - //Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ); - //funcval = 0.0f; - //for (j = 0; j 2) { - //if (re > re_old) { - //Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ); - //funcval = 0.0f; - //for (j = 0; j(npD); - //result.append(out1); - //result.append(ll); - //} - - //return result; -//} - bp::list LLT_model(np::ndarray input, double d_lambda, double d_tau, int iter, double d_epsil, int switcher) { // the result is in the following list bp::list result; @@ -1022,9 +736,6 @@ bp::list TGV_PD(np::ndarray input, double d_lambda, double d_alpha1, double d_al result.append(npU); } - - - return result; } diff --git a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py index a9d0f31..63be1a0 100644 --- a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py +++ b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py @@ -3,7 +3,7 @@ """ Created on Thu Feb 22 11:39:43 2018 -Testing CPU implementation against GPU one +Testing CPU implementation against the GPU one @author: Daniil Kazantsev """ @@ -12,8 +12,6 @@ import matplotlib.pyplot as plt import numpy as np import os import timeit -#from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU, TV_FGP_CPU -#from ccpi.filters.gpu_regularizers import TV_ROF_GPU, TV_FGP_GPU from ccpi.filters.regularizers import ROF_TV, FGP_TV ############################################################################### diff --git a/Wrappers/Python/test/test_gpu_regularizers.py b/Wrappers/Python/test/test_gpu_regularizers.py index 04aeeb4..640b3f9 100644 --- a/Wrappers/Python/test/test_gpu_regularizers.py +++ b/Wrappers/Python/test/test_gpu_regularizers.py @@ -11,8 +11,8 @@ import numpy as np import os from enum import Enum import timeit -from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, TV_ROF_GPU - +from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML +from ccpi.filters.regularizers import ROF_TV, FGP_TV ############################################################################### def printParametersToString(pars): txt = r'' @@ -151,20 +151,19 @@ plt.colorbar(ticks=[0, 0.03], orientation='vertical') ## Rudin-Osher-Fatemi (ROF) TV regularization start_time = timeit.default_timer() - pars = { -'algorithm' : TV_ROF_GPU , \ +'algorithm' : ROF_TV , \ 'input' : u0, 'regularization_parameter': 0.04,\ 'number_of_iterations':300,\ 'time_marching_parameter': 0.0025 } - + rof_tv = TV_ROF_GPU(pars['input'], pars['regularization_parameter'], pars['number_of_iterations'], - pars['time_marching_parameter']) + pars['time_marching_parameter'],'gpu') rms = rmse(Im, rof_tv) pars['rmse'] = rms @@ -190,3 +189,51 @@ a.text(0.05, 0.95, 'rof_tv - u0', transform=a.transAxes, fontsize=14, imgplot = plt.imshow((rof_tv - u0)**2, vmin=0, vmax=0.03, cmap="gray") plt.colorbar(ticks=[0, 0.03], orientation='vertical') plt.show() + +## Fast-Gradient Projection TV regularization +""" +start_time = timeit.default_timer() + +pars = {'algorithm' : FGP_TV, \ + 'input' : u0,\ + 'regularization_parameter':0.04, \ + 'number_of_iterations' :1200 ,\ + 'tolerance_constant':0.00001,\ + 'methodTV': 0 ,\ + 'nonneg': 0 ,\ + 'printingOut': 0 + } + +fgp_gpu = FGP_TV(pars['input'], + pars['regularization_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['nonneg'], + pars['printingOut'],'gpu') + +rms = rmse(Im, fgp_gpu) +pars['rmse'] = rms +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(2,4,4) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) +# place a text box in upper left in axes coords +a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=12, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(fgp_gpu, cmap="gray") + +a=fig.add_subplot(2,4,8) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +# place a text box in upper left in axes coords +a.text(0.05, 0.95, 'fgp_gpu - u0', transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow((fgp_gpu - u0)**2, vmin=0, vmax=0.03, cmap="gray") +plt.colorbar(ticks=[0, 0.03], orientation='vertical') +plt.show() +""" -- cgit v1.2.3 From fbb7a12f7714978c251f02bf84ab9a66c762f428 Mon Sep 17 00:00:00 2001 From: algol Date: Wed, 7 Mar 2018 10:15:13 +0000 Subject: small correct #40 --- Wrappers/Python/src/cpu_regularizers.cpp | 1 - 1 file changed, 1 deletion(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/src/cpu_regularizers.cpp b/Wrappers/Python/src/cpu_regularizers.cpp index 88b7c3c..3529ebd 100644 --- a/Wrappers/Python/src/cpu_regularizers.cpp +++ b/Wrappers/Python/src/cpu_regularizers.cpp @@ -27,7 +27,6 @@ limitations under the License. #include "boost/tuple/tuple.hpp" #include "SplitBregman_TV_core.h" -//#include "FGP_TV_core.h" #include "LLT_model_core.h" #include "PatchBased_Regul_core.h" #include "TGV_PD_core.h" -- cgit v1.2.3