From b9fafd363d1d181a4a8b42ea4038924097207913 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Mon, 9 Apr 2018 13:41:06 +0100 Subject: major renaming and new 3D demos for Matlab --- Wrappers/Python/test/test_cpu_regularisers.py | 442 +++++++++++++++++++++ Wrappers/Python/test/test_cpu_regularizers.py | 442 --------------------- .../Python/test/test_cpu_vs_gpu_regularisers.py | 219 ++++++++++ .../Python/test/test_cpu_vs_gpu_regularizers.py | 219 ---------- Wrappers/Python/test/test_gpu_regularisers.py | 143 +++++++ Wrappers/Python/test/test_gpu_regularizers.py | 239 ----------- Wrappers/Python/test/test_regularisers_3d.py | 425 ++++++++++++++++++++ Wrappers/Python/test/test_regularizers.py | 395 ------------------ Wrappers/Python/test/test_regularizers_3d.py | 425 -------------------- 9 files changed, 1229 insertions(+), 1720 deletions(-) create mode 100644 Wrappers/Python/test/test_cpu_regularisers.py delete mode 100644 Wrappers/Python/test/test_cpu_regularizers.py create mode 100644 Wrappers/Python/test/test_cpu_vs_gpu_regularisers.py delete mode 100644 Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py create mode 100644 Wrappers/Python/test/test_gpu_regularisers.py delete mode 100644 Wrappers/Python/test/test_gpu_regularizers.py create mode 100644 Wrappers/Python/test/test_regularisers_3d.py delete mode 100644 Wrappers/Python/test/test_regularizers.py delete mode 100644 Wrappers/Python/test/test_regularizers_3d.py (limited to 'Wrappers/Python/test') diff --git a/Wrappers/Python/test/test_cpu_regularisers.py b/Wrappers/Python/test/test_cpu_regularisers.py new file mode 100644 index 0000000..9713baa --- /dev/null +++ b/Wrappers/Python/test/test_cpu_regularisers.py @@ -0,0 +1,442 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Aug 4 11:10:05 2017 + +@author: ofn77899 +""" + + +import matplotlib.pyplot as plt +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 ,\ + TGV_PD + +############################################################################### +#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 +#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum +# intensity from the two images being compared, and respectively the same for +# minval. RMSE is given by the square root of MSE: +# sqrt[(sum(A - B) ** 2) / |A|], +# where |A| means the number of elements in A. By doing this, the maximum value +# given by RMSE is maxval. + +def nrmse(im1, im2): + a, b = im1.shape + rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b)) + max_val = max(np.max(im1), np.max(im2)) + min_val = min(np.min(im1), np.min(im2)) + return 1 - (rmse / (max_val - min_val)) +############################################################################### +def 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 +############################################################################### +# +# 2D Regularizers +# +############################################################################### +#Example: +# figure; +# Im = double(imread('lena_gray_256.tif'))/255; % loading image +# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; +# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); + +# assumes the script is launched from the test directory +filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif") +#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif" +#filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif" +#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif' + +Im = plt.imread(filename) +Im = np.asarray(Im, dtype='float32') + +perc = 0.15 +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() + +a=fig.add_subplot(2,3,1) +a.set_title('noise') +imgplot = plt.imshow(u0#,cmap="gray" + ) + +reg_output = [] +############################################################################## +# Call regularizer + +####################### SplitBregman_TV ##################################### +# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); + +start_time = timeit.default_timer() +pars = {'algorithm' : SplitBregman_TV , \ + 'input' : u0, + 'regularization_parameter':10. , \ +'number_of_iterations' :35 ,\ +'tolerance_constant':0.0001 , \ +'TV_penalty': 0 +} + +out = SplitBregman_TV (pars['input'], pars['regularization_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['TV_penalty']) +splitbregman = out[0] +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) + + +a=fig.add_subplot(2,3,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, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(splitbregman,\ + #cmap="gray" + ) + +###################### FGP_TV ######################################### +# u = FGP_TV(single(u0), 0.05, 100, 1e-04); +start_time = timeit.default_timer() +pars = {'algorithm' : FGP_TV , \ + 'input' : u0, + 'regularization_parameter':5e-4, \ + 'number_of_iterations' :10 ,\ + 'tolerance_constant':0.001,\ + 'TV_penalty': 0 +} + +out = FGP_TV (pars['input'], + pars['regularization_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['TV_penalty']) + +fgp = out[0] +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) + + +a=fig.add_subplot(2,3,3) + +# 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 +imgplot = plt.imshow(fgp, \ + #cmap="gray" + ) +# 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) + +###################### LLT_model ######################################### + +start_time = timeit.default_timer() + +pars = {'algorithm': LLT_model , \ + 'input' : u0, + 'regularization_parameter': 25,\ + 'time_step':0.0003, \ + 'number_of_iterations' :300,\ + 'tolerance_constant':0.001,\ + 'restrictive_Z_smoothing': 0 +} +out = LLT_model(pars['input'], + pars['regularization_parameter'], + pars['time_step'] , + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['restrictive_Z_smoothing'] ) + +llt = out[0] +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(2,3,4) + +# 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, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(llt,\ + #cmap="gray" + ) + + +# ###################### PatchBased_Regul ######################################### +# # Quick 2D denoising example in Matlab: +# # Im = double(imread('lena_gray_256.tif'))/255; % loading image +# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise +# # ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); +start_time = timeit.default_timer() + +pars = {'algorithm': PatchBased_Regul , \ + 'input' : u0, + 'regularization_parameter': 0.05,\ + 'searching_window_ratio':3, \ + 'similarity_window_ratio':1,\ + 'PB_filtering_parameter': 0.08 +} +out = PatchBased_Regul(pars['input'], + pars['regularization_parameter'], + pars['searching_window_ratio'] , + pars['similarity_window_ratio'] , + pars['PB_filtering_parameter']) +pbr = out[0] +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) + +a=fig.add_subplot(2,3,5) + + +# 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, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(pbr #,cmap="gray" + ) + + +# ###################### TGV_PD ######################################### +# # Quick 2D denoising example in Matlab: +# # Im = double(imread('lena_gray_256.tif'))/255; % loading image +# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise +# # u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); + +start_time = timeit.default_timer() + +pars = {'algorithm': TGV_PD , \ + 'input' : u0,\ + 'regularization_parameter':0.05,\ + 'first_order_term': 1.3,\ + 'second_order_term': 1, \ + 'number_of_iterations': 550 + } +out = TGV_PD(pars['input'], + pars['regularization_parameter'], + pars['first_order_term'] , + pars['second_order_term'] , + pars['number_of_iterations']) +tgv = out[0] +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(2,3,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, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(tgv #, cmap="gray") + ) + + +plt.show() + +################################################################################ +## +## 3D Regularizers +## +################################################################################ +##Example: +## figure; +## Im = double(imread('lena_gray_256.tif'))/255; % loading image +## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; +## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); +# +##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha" +#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha" +# +#reader = vtk.vtkMetaImageReader() +#reader.SetFileName(os.path.normpath(filename)) +#reader.Update() +##vtk returns 3D images, let's take just the one slice there is as 2D +#Im = Converter.vtk2numpy(reader.GetOutput()) +#Im = Im.astype('float32') +##imgplot = plt.imshow(Im) +#perc = 0.05 +#u0 = Im + (perc* np.random.normal(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') +#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(), +# reader.GetOutput().GetOrigin()) +#converter.Update() +#writer = vtk.vtkMetaImageWriter() +#writer.SetInputData(converter.GetOutput()) +#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha") +##writer.Write() +# +# +### plot +#fig3D = plt.figure() +##a=fig.add_subplot(3,3,1) +##a.set_title('Original') +##imgplot = plt.imshow(Im) +#sliceNo = 32 +# +#a=fig3D.add_subplot(2,3,1) +#a.set_title('noise') +#imgplot = plt.imshow(u0.T[sliceNo]) +# +#reg_output3d = [] +# +############################################################################### +## Call regularizer +# +######################## SplitBregman_TV ##################################### +## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); +# +##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV) +# +##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30, +## #tolerance_constant=1e-4, +## TV_Penalty=Regularizer.TotalVariationPenalty.l1) +# +#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30, +# tolerance_constant=1e-4, +# TV_Penalty=Regularizer.TotalVariationPenalty.l1) +# +# +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## 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, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) +# +####################### FGP_TV ######################################### +## u = FGP_TV(single(u0), 0.05, 100, 1e-04); +#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005, +# number_of_iterations=200) +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## 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, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) +# +####################### LLT_model ######################################### +## * u0 = Im + .03*randn(size(Im)); % adding noise +## [Den] = LLT_model(single(u0), 10, 0.1, 1); +##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); +##input, regularization_parameter , time_step, number_of_iterations, +## tolerance_constant, restrictive_Z_smoothing=0 +#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25, +# time_step=0.0003, +# tolerance_constant=0.0001, +# number_of_iterations=300) +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## 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, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) +# +####################### PatchBased_Regul ######################################### +## Quick 2D denoising example in Matlab: +## Im = double(imread('lena_gray_256.tif'))/255; % loading image +## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise +## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); +# +#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05, +# searching_window_ratio=3, +# similarity_window_ratio=1, +# PB_filtering_parameter=0.08) +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## 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, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) +# + +###################### TGV_PD ######################################### +# Quick 2D denoising example in Matlab: +# Im = double(imread('lena_gray_256.tif'))/255; % loading image +# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise +# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); + + +#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05, +# first_order_term=1.3, +# second_order_term=1, +# number_of_iterations=550) +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## 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, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) diff --git a/Wrappers/Python/test/test_cpu_regularizers.py b/Wrappers/Python/test/test_cpu_regularizers.py deleted file mode 100644 index 9713baa..0000000 --- a/Wrappers/Python/test/test_cpu_regularizers.py +++ /dev/null @@ -1,442 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Fri Aug 4 11:10:05 2017 - -@author: ofn77899 -""" - - -import matplotlib.pyplot as plt -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 ,\ - TGV_PD - -############################################################################### -#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 -#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum -# intensity from the two images being compared, and respectively the same for -# minval. RMSE is given by the square root of MSE: -# sqrt[(sum(A - B) ** 2) / |A|], -# where |A| means the number of elements in A. By doing this, the maximum value -# given by RMSE is maxval. - -def nrmse(im1, im2): - a, b = im1.shape - rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b)) - max_val = max(np.max(im1), np.max(im2)) - min_val = min(np.min(im1), np.min(im2)) - return 1 - (rmse / (max_val - min_val)) -############################################################################### -def 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 -############################################################################### -# -# 2D Regularizers -# -############################################################################### -#Example: -# figure; -# Im = double(imread('lena_gray_256.tif'))/255; % loading image -# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; -# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); - -# assumes the script is launched from the test directory -filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif") -#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif" -#filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif" -#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif' - -Im = plt.imread(filename) -Im = np.asarray(Im, dtype='float32') - -perc = 0.15 -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() - -a=fig.add_subplot(2,3,1) -a.set_title('noise') -imgplot = plt.imshow(u0#,cmap="gray" - ) - -reg_output = [] -############################################################################## -# Call regularizer - -####################### SplitBregman_TV ##################################### -# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); - -start_time = timeit.default_timer() -pars = {'algorithm' : SplitBregman_TV , \ - 'input' : u0, - 'regularization_parameter':10. , \ -'number_of_iterations' :35 ,\ -'tolerance_constant':0.0001 , \ -'TV_penalty': 0 -} - -out = SplitBregman_TV (pars['input'], pars['regularization_parameter'], - pars['number_of_iterations'], - pars['tolerance_constant'], - pars['TV_penalty']) -splitbregman = out[0] -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) - - -a=fig.add_subplot(2,3,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, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(splitbregman,\ - #cmap="gray" - ) - -###################### FGP_TV ######################################### -# u = FGP_TV(single(u0), 0.05, 100, 1e-04); -start_time = timeit.default_timer() -pars = {'algorithm' : FGP_TV , \ - 'input' : u0, - 'regularization_parameter':5e-4, \ - 'number_of_iterations' :10 ,\ - 'tolerance_constant':0.001,\ - 'TV_penalty': 0 -} - -out = FGP_TV (pars['input'], - pars['regularization_parameter'], - pars['number_of_iterations'], - pars['tolerance_constant'], - pars['TV_penalty']) - -fgp = out[0] -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) - - -a=fig.add_subplot(2,3,3) - -# 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 -imgplot = plt.imshow(fgp, \ - #cmap="gray" - ) -# 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) - -###################### LLT_model ######################################### - -start_time = timeit.default_timer() - -pars = {'algorithm': LLT_model , \ - 'input' : u0, - 'regularization_parameter': 25,\ - 'time_step':0.0003, \ - 'number_of_iterations' :300,\ - 'tolerance_constant':0.001,\ - 'restrictive_Z_smoothing': 0 -} -out = LLT_model(pars['input'], - pars['regularization_parameter'], - pars['time_step'] , - pars['number_of_iterations'], - pars['tolerance_constant'], - pars['restrictive_Z_smoothing'] ) - -llt = out[0] -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(2,3,4) - -# 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, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(llt,\ - #cmap="gray" - ) - - -# ###################### PatchBased_Regul ######################################### -# # Quick 2D denoising example in Matlab: -# # Im = double(imread('lena_gray_256.tif'))/255; % loading image -# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -# # ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); -start_time = timeit.default_timer() - -pars = {'algorithm': PatchBased_Regul , \ - 'input' : u0, - 'regularization_parameter': 0.05,\ - 'searching_window_ratio':3, \ - 'similarity_window_ratio':1,\ - 'PB_filtering_parameter': 0.08 -} -out = PatchBased_Regul(pars['input'], - pars['regularization_parameter'], - pars['searching_window_ratio'] , - pars['similarity_window_ratio'] , - pars['PB_filtering_parameter']) -pbr = out[0] -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) - -a=fig.add_subplot(2,3,5) - - -# 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, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(pbr #,cmap="gray" - ) - - -# ###################### TGV_PD ######################################### -# # Quick 2D denoising example in Matlab: -# # Im = double(imread('lena_gray_256.tif'))/255; % loading image -# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -# # u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); - -start_time = timeit.default_timer() - -pars = {'algorithm': TGV_PD , \ - 'input' : u0,\ - 'regularization_parameter':0.05,\ - 'first_order_term': 1.3,\ - 'second_order_term': 1, \ - 'number_of_iterations': 550 - } -out = TGV_PD(pars['input'], - pars['regularization_parameter'], - pars['first_order_term'] , - pars['second_order_term'] , - pars['number_of_iterations']) -tgv = out[0] -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(2,3,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, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(tgv #, cmap="gray") - ) - - -plt.show() - -################################################################################ -## -## 3D Regularizers -## -################################################################################ -##Example: -## figure; -## Im = double(imread('lena_gray_256.tif'))/255; % loading image -## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; -## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); -# -##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha" -#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha" -# -#reader = vtk.vtkMetaImageReader() -#reader.SetFileName(os.path.normpath(filename)) -#reader.Update() -##vtk returns 3D images, let's take just the one slice there is as 2D -#Im = Converter.vtk2numpy(reader.GetOutput()) -#Im = Im.astype('float32') -##imgplot = plt.imshow(Im) -#perc = 0.05 -#u0 = Im + (perc* np.random.normal(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') -#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(), -# reader.GetOutput().GetOrigin()) -#converter.Update() -#writer = vtk.vtkMetaImageWriter() -#writer.SetInputData(converter.GetOutput()) -#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha") -##writer.Write() -# -# -### plot -#fig3D = plt.figure() -##a=fig.add_subplot(3,3,1) -##a.set_title('Original') -##imgplot = plt.imshow(Im) -#sliceNo = 32 -# -#a=fig3D.add_subplot(2,3,1) -#a.set_title('noise') -#imgplot = plt.imshow(u0.T[sliceNo]) -# -#reg_output3d = [] -# -############################################################################### -## Call regularizer -# -######################## SplitBregman_TV ##################################### -## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); -# -##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV) -# -##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30, -## #tolerance_constant=1e-4, -## TV_Penalty=Regularizer.TotalVariationPenalty.l1) -# -#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30, -# tolerance_constant=1e-4, -# TV_Penalty=Regularizer.TotalVariationPenalty.l1) -# -# -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## 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, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# -####################### FGP_TV ######################################### -## u = FGP_TV(single(u0), 0.05, 100, 1e-04); -#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005, -# number_of_iterations=200) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## 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, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# -####################### LLT_model ######################################### -## * u0 = Im + .03*randn(size(Im)); % adding noise -## [Den] = LLT_model(single(u0), 10, 0.1, 1); -##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); -##input, regularization_parameter , time_step, number_of_iterations, -## tolerance_constant, restrictive_Z_smoothing=0 -#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25, -# time_step=0.0003, -# tolerance_constant=0.0001, -# number_of_iterations=300) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## 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, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# -####################### PatchBased_Regul ######################################### -## Quick 2D denoising example in Matlab: -## Im = double(imread('lena_gray_256.tif'))/255; % loading image -## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); -# -#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05, -# searching_window_ratio=3, -# similarity_window_ratio=1, -# PB_filtering_parameter=0.08) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## 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, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# - -###################### TGV_PD ######################################### -# Quick 2D denoising example in Matlab: -# Im = double(imread('lena_gray_256.tif'))/255; % loading image -# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); - - -#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05, -# first_order_term=1.3, -# second_order_term=1, -# number_of_iterations=550) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## 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, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) diff --git a/Wrappers/Python/test/test_cpu_vs_gpu_regularisers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularisers.py new file mode 100644 index 0000000..15e9042 --- /dev/null +++ b/Wrappers/Python/test/test_cpu_vs_gpu_regularisers.py @@ -0,0 +1,219 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Thu Feb 22 11:39:43 2018 + +Testing CPU implementation against the GPU one + +@author: Daniil Kazantsev +""" + +import matplotlib.pyplot as plt +import numpy as np +import os +import timeit +from ccpi.filters.regularisers import ROF_TV, FGP_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') + + +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("____________ROF-TV bench___________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure(1) +plt.suptitle('Comparison of ROF-TV regulariser 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,\ + 'regularisation_parameter':0.04,\ + 'number_of_iterations': 1200,\ + 'time_marching_parameter': 0.0025 + } +print ("#############ROF TV CPU####################") +start_time = timeit.default_timer() +rof_cpu = ROF_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'cpu') +rms = rmse(Im, rof_cpu) +pars['rmse'] = rms + +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +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 = ROF_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'gpu') + +rms = rmse(Im, rof_gpu) +pars['rmse'] = rms +pars['algorithm'] = ROF_TV +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +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-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) +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") + +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("____________FGP-TV bench___________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure(2) +plt.suptitle('Comparison of FGP-TV regulariser 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,\ + 'regularisation_parameter':0.04, \ + 'number_of_iterations' :1200 ,\ + 'tolerance_constant':0.00001,\ + 'methodTV': 0 ,\ + 'nonneg': 0 ,\ + 'printingOut': 0 + } + +print ("#############FGP TV CPU####################") +start_time = timeit.default_timer() +fgp_cpu = FGP_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['nonneg'], + pars['printingOut'],'cpu') + + +rms = rmse(Im, fgp_cpu) +pars['rmse'] = rms + +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +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['regularisation_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_cpu_vs_gpu_regularizers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py deleted file mode 100644 index 63be1a0..0000000 --- a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py +++ /dev/null @@ -1,219 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Thu Feb 22 11:39:43 2018 - -Testing CPU implementation against the GPU one - -@author: Daniil Kazantsev -""" - -import matplotlib.pyplot as plt -import numpy as np -import os -import timeit -from ccpi.filters.regularizers import ROF_TV, FGP_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') - - -print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("____________ROF-TV bench___________________") -print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") - -## 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':0.04,\ - 'number_of_iterations': 1200,\ - 'time_marching_parameter': 0.0025 - } -print ("#############ROF TV CPU####################") -start_time = timeit.default_timer() -rof_cpu = ROF_TV(pars['input'], - pars['regularization_parameter'], - pars['number_of_iterations'], - pars['time_marching_parameter'],'cpu') -rms = rmse(Im, rof_cpu) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -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 = ROF_TV(pars['input'], - pars['regularization_parameter'], - pars['number_of_iterations'], - pars['time_marching_parameter'],'gpu') - -rms = rmse(Im, rof_gpu) -pars['rmse'] = rms -pars['algorithm'] = ROF_TV -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -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-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) -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") - -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_regularisers.py b/Wrappers/Python/test/test_gpu_regularisers.py new file mode 100644 index 0000000..2103c0e --- /dev/null +++ b/Wrappers/Python/test/test_gpu_regularisers.py @@ -0,0 +1,143 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +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.regularisers import ROF_TV, FGP_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") + +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() + +a=fig.add_subplot(2,4,1) +a.set_title('noise') +imgplot = plt.imshow(u0,cmap="gray") + + + +## Rudin-Osher-Fatemi (ROF) TV regularisation +start_time = timeit.default_timer() +pars = { +'algorithm' : ROF_TV , \ + 'input' : u0, + 'regularisation_parameter': 0.04,\ + 'number_of_iterations':300,\ + 'time_marching_parameter': 0.0025 + + } + +rof_tv = TV_ROF_GPU(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'gpu') + +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,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(rof_tv, 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, 'rof_tv - u0', transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +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 regularisation +""" +start_time = timeit.default_timer() + +pars = {'algorithm' : FGP_TV, \ + 'input' : u0,\ + 'regularisation_parameter':0.04, \ + 'number_of_iterations' :1200 ,\ + 'tolerance_constant':0.00001,\ + 'methodTV': 0 ,\ + 'nonneg': 0 ,\ + 'printingOut': 0 + } + +fgp_gpu = FGP_TV(pars['input'], + pars['regularisation_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() +""" diff --git a/Wrappers/Python/test/test_gpu_regularizers.py b/Wrappers/Python/test/test_gpu_regularizers.py deleted file mode 100644 index 640b3f9..0000000 --- a/Wrappers/Python/test/test_gpu_regularizers.py +++ /dev/null @@ -1,239 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -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.regularizers import ROF_TV, FGP_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") -#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif" -#filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif" -#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif' - -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() - -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 , \ - 'input' : u0, - 'edge_preserv_parameter':0.1 , \ -'number_of_iterations' :250 ,\ -'time_marching_parameter':0.03 ,\ -'regularization_parameter':0.7 -} - - -d4h = Diff4thHajiaboli(pars['input'], - pars['edge_preserv_parameter'], - pars['number_of_iterations'], - pars['time_marching_parameter'], - pars['regularization_parameter']) -rms = rmse(Im, d4h) -pars['rmse'] = rms -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(2,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=12, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(d4h, cmap="gray") - -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, vmin=0, vmax=0.03, cmap="gray") -plt.colorbar(ticks=[0, 0.03], orientation='vertical') - - -## Patch Based Regul NML -start_time = timeit.default_timer() -""" -pars = {'algorithm' : NML , \ - 'input' : u0, - 'SearchW_real':3 , \ -'SimilW' :1,\ -'h':0.05 ,# -'lambda' : 0.08 -} -""" -pars = { -'algorithm' : NML , \ - 'input' : u0, - 'regularization_parameter': 0.01,\ - 'searching_window_ratio':3, \ - 'similarity_window_ratio':1,\ - 'PB_filtering_parameter': 0.2 -} - -nml = NML(pars['input'], - pars['searching_window_ratio'], - pars['similarity_window_ratio'], - pars['PB_filtering_parameter'], - pars['regularization_parameter']) -rms = rmse(Im, nml) -pars['rmse'] = rms -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(2,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=12, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(nml, cmap="gray") - -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, vmin=0, vmax=0.03, cmap="gray") -plt.colorbar(ticks=[0, 0.03], orientation='vertical') - - - -## Rudin-Osher-Fatemi (ROF) TV regularization -start_time = timeit.default_timer() -pars = { -'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'],'gpu') - -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,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(rof_tv, 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, 'rof_tv - u0', transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -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() -""" diff --git a/Wrappers/Python/test/test_regularisers_3d.py b/Wrappers/Python/test/test_regularisers_3d.py new file mode 100644 index 0000000..2d11a7e --- /dev/null +++ b/Wrappers/Python/test/test_regularisers_3d.py @@ -0,0 +1,425 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Aug 4 11:10:05 2017 + +@author: ofn77899 +""" + +#from ccpi.viewer.CILViewer2D import Converter +#import vtk + +import matplotlib.pyplot as plt +import numpy as np +import os +from enum import Enum +import timeit +#from PIL import Image +#from Regularizer import Regularizer +from ccpi.imaging.Regularizer import Regularizer + +############################################################################### +#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 +#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum +# intensity from the two images being compared, and respectively the same for +# minval. RMSE is given by the square root of MSE: +# sqrt[(sum(A - B) ** 2) / |A|], +# where |A| means the number of elements in A. By doing this, the maximum value +# given by RMSE is maxval. + +def nrmse(im1, im2): + a, b = im1.shape + rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b)) + max_val = max(np.max(im1), np.max(im2)) + min_val = min(np.min(im1), np.min(im2)) + return 1 - (rmse / (max_val - min_val)) +############################################################################### + +############################################################################### +# +# 2D Regularizers +# +############################################################################### +#Example: +# figure; +# Im = double(imread('lena_gray_256.tif'))/255; % loading image +# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; +# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); + + +#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif" +filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif" +#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif' + +#reader = vtk.vtkTIFFReader() +#reader.SetFileName(os.path.normpath(filename)) +#reader.Update() +Im = plt.imread(filename) +#Im = Image.open('/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif')/255 +#img.show() +Im = np.asarray(Im, dtype='float32') + +# create a 3D image by stacking N of this images + + +#imgplot = plt.imshow(Im) +perc = 0.05 +u_n = Im + (perc* np.random.normal(size=np.shape(Im))) +y,z = np.shape(u_n) +u_n = np.reshape(u_n , (1,y,z)) + +u0 = u_n.copy() +for i in range (19): + u_n = Im + (perc* np.random.normal(size=np.shape(Im))) + u_n = np.reshape(u_n , (1,y,z)) + + u0 = np.vstack ( (u0, u_n) ) + +# map the u0 u0->u0>0 +f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1) +u0 = f(u0).astype('float32') + +print ("Passed image shape {0}".format(np.shape(u0))) + +## plot +fig = plt.figure() +#a=fig.add_subplot(3,3,1) +#a.set_title('Original') +#imgplot = plt.imshow(Im) +sliceno = 10 + +a=fig.add_subplot(2,3,1) +a.set_title('noise') +imgplot = plt.imshow(u0[sliceno],cmap="gray") + +reg_output = [] +############################################################################## +# Call regularizer + +####################### SplitBregman_TV ##################################### +# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); + +use_object = True +if use_object: + reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV) + print (reg.pars) + reg.setParameter(input=u0) + reg.setParameter(regularization_parameter=10.) + # or + # reg.setParameter(input=u0, regularization_parameter=10., #number_of_iterations=30, + #tolerance_constant=1e-4, + #TV_Penalty=Regularizer.TotalVariationPenalty.l1) + plotme = reg() [0] + pars = reg.pars + textstr = reg.printParametersToString() + + #out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30, + #tolerance_constant=1e-4, + # TV_Penalty=Regularizer.TotalVariationPenalty.l1) + +#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30, +# tolerance_constant=1e-4, +# TV_Penalty=Regularizer.TotalVariationPenalty.l1) + +else: + out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10. ) + pars = out2[2] + reg_output.append(out2) + plotme = reg_output[-1][0] + textstr = out2[-1] + +a=fig.add_subplot(2,3,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, textstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(plotme[sliceno],cmap="gray") + +###################### FGP_TV ######################################### +# u = FGP_TV(single(u0), 0.05, 100, 1e-04); +out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.0005, + number_of_iterations=50) +pars = out2[-2] + +reg_output.append(out2) + +a=fig.add_subplot(2,3,3) + +textstr = out2[-1] + +# 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, textstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(reg_output[-1][0][sliceno]) +# place a text box in upper left in axes coords +a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray") + +###################### LLT_model ######################################### +# * u0 = Im + .03*randn(size(Im)); % adding noise +# [Den] = LLT_model(single(u0), 10, 0.1, 1); +#Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); +#input, regularization_parameter , time_step, number_of_iterations, +# tolerance_constant, restrictive_Z_smoothing=0 +out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25, + time_step=0.0003, + tolerance_constant=0.0001, + number_of_iterations=300) +pars = out2[-2] + +reg_output.append(out2) + +a=fig.add_subplot(2,3,4) + +textstr = out2[-1] + +# 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, textstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray") + + +# ###################### PatchBased_Regul ######################################### +# # Quick 2D denoising example in Matlab: +# # Im = double(imread('lena_gray_256.tif'))/255; % loading image +# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise +# # ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); + +out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05, + searching_window_ratio=3, + similarity_window_ratio=1, + PB_filtering_parameter=0.08) +pars = out2[-2] +reg_output.append(out2) + +a=fig.add_subplot(2,3,5) + + +textstr = out2[-1] + +# 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, textstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray") + + +# ###################### TGV_PD ######################################### +# # Quick 2D denoising example in Matlab: +# # Im = double(imread('lena_gray_256.tif'))/255; % loading image +# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise +# # u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); + + +out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05, + first_order_term=1.3, + second_order_term=1, + number_of_iterations=550) +pars = out2[-2] +reg_output.append(out2) + +a=fig.add_subplot(2,3,6) + + +textstr = out2[-1] + + +# 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, textstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray") + + +plt.show() + +################################################################################ +## +## 3D Regularizers +## +################################################################################ +##Example: +## figure; +## Im = double(imread('lena_gray_256.tif'))/255; % loading image +## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; +## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); +# +##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha" +#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha" +# +#reader = vtk.vtkMetaImageReader() +#reader.SetFileName(os.path.normpath(filename)) +#reader.Update() +##vtk returns 3D images, let's take just the one slice there is as 2D +#Im = Converter.vtk2numpy(reader.GetOutput()) +#Im = Im.astype('float32') +##imgplot = plt.imshow(Im) +#perc = 0.05 +#u0 = Im + (perc* np.random.normal(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') +#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(), +# reader.GetOutput().GetOrigin()) +#converter.Update() +#writer = vtk.vtkMetaImageWriter() +#writer.SetInputData(converter.GetOutput()) +#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha") +##writer.Write() +# +# +### plot +#fig3D = plt.figure() +##a=fig.add_subplot(3,3,1) +##a.set_title('Original') +##imgplot = plt.imshow(Im) +#sliceNo = 32 +# +#a=fig3D.add_subplot(2,3,1) +#a.set_title('noise') +#imgplot = plt.imshow(u0.T[sliceNo]) +# +#reg_output3d = [] +# +############################################################################### +## Call regularizer +# +######################## SplitBregman_TV ##################################### +## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); +# +##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV) +# +##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30, +## #tolerance_constant=1e-4, +## TV_Penalty=Regularizer.TotalVariationPenalty.l1) +# +#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30, +# tolerance_constant=1e-4, +# TV_Penalty=Regularizer.TotalVariationPenalty.l1) +# +# +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## 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, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) +# +####################### FGP_TV ######################################### +## u = FGP_TV(single(u0), 0.05, 100, 1e-04); +#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005, +# number_of_iterations=200) +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## 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, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) +# +####################### LLT_model ######################################### +## * u0 = Im + .03*randn(size(Im)); % adding noise +## [Den] = LLT_model(single(u0), 10, 0.1, 1); +##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); +##input, regularization_parameter , time_step, number_of_iterations, +## tolerance_constant, restrictive_Z_smoothing=0 +#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25, +# time_step=0.0003, +# tolerance_constant=0.0001, +# number_of_iterations=300) +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## 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, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) +# +####################### PatchBased_Regul ######################################### +## Quick 2D denoising example in Matlab: +## Im = double(imread('lena_gray_256.tif'))/255; % loading image +## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise +## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); +# +#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05, +# searching_window_ratio=3, +# similarity_window_ratio=1, +# PB_filtering_parameter=0.08) +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## 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, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) +# + +###################### TGV_PD ######################################### +# Quick 2D denoising example in Matlab: +# Im = double(imread('lena_gray_256.tif'))/255; % loading image +# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise +# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); + + +#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05, +# first_order_term=1.3, +# second_order_term=1, +# number_of_iterations=550) +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## 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, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) diff --git a/Wrappers/Python/test/test_regularizers.py b/Wrappers/Python/test/test_regularizers.py deleted file mode 100644 index cf5da2b..0000000 --- a/Wrappers/Python/test/test_regularizers.py +++ /dev/null @@ -1,395 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Fri Aug 4 11:10:05 2017 - -@author: ofn77899 -""" - -#from ccpi.viewer.CILViewer2D import Converter -#import vtk - -import matplotlib.pyplot as plt -import numpy as np -import os -from enum import Enum -import timeit -#from PIL import Image -#from Regularizer import Regularizer -from ccpi.filters.Regularizer import Regularizer - -############################################################################### -#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 -#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum -# intensity from the two images being compared, and respectively the same for -# minval. RMSE is given by the square root of MSE: -# sqrt[(sum(A - B) ** 2) / |A|], -# where |A| means the number of elements in A. By doing this, the maximum value -# given by RMSE is maxval. - -def nrmse(im1, im2): - a, b = im1.shape - rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b)) - max_val = max(np.max(im1), np.max(im2)) - min_val = min(np.min(im1), np.min(im2)) - return 1 - (rmse / (max_val - min_val)) -############################################################################### - -############################################################################### -# -# 2D Regularizers -# -############################################################################### -#Example: -# figure; -# Im = double(imread('lena_gray_256.tif'))/255; % loading image -# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; -# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); - - -filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif" -#filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif" -#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif' - -#reader = vtk.vtkTIFFReader() -#reader.SetFileName(os.path.normpath(filename)) -#reader.Update() -Im = plt.imread(filename) -#Im = Image.open('/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif')/255 -#img.show() -Im = np.asarray(Im, dtype='float32') - - - - -#imgplot = plt.imshow(Im) -perc = 0.05 -u0 = Im + (perc* np.random.normal(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() -#a=fig.add_subplot(3,3,1) -#a.set_title('Original') -#imgplot = plt.imshow(Im) - -a=fig.add_subplot(2,3,1) -a.set_title('noise') -imgplot = plt.imshow(u0,cmap="gray") - -reg_output = [] -############################################################################## -# Call regularizer - -####################### SplitBregman_TV ##################################### -# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); - -start_time = timeit.default_timer() -reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV) -print (reg.pars) -reg.setParameter(input=u0) -reg.setParameter(regularization_parameter=10.) -# or -# reg.setParameter(input=u0, regularization_parameter=10., #number_of_iterations=30, - #tolerance_constant=1e-4, - #TV_Penalty=Regularizer.TotalVariationPenalty.l1) -plotme = reg(output_all=True) [0] -pars = reg.pars -txtstr = reg.printParametersToString() -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) - - -a=fig.add_subplot(2,3,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, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(plotme,cmap="gray") - -###################### FGP_TV ######################################### -# u = FGP_TV(single(u0), 0.05, 100, 1e-04); -start_time = timeit.default_timer() -reg = Regularizer(Regularizer.Algorithm.FGP_TV) -out2 = reg(input=u0, regularization_parameter=5e-4, - number_of_iterations=10) -txtstr = reg.printParametersToString() -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) - - -a=fig.add_subplot(2,3,3) - -# 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 -imgplot = plt.imshow(out2,cmap="gray") -# 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) - -###################### LLT_model ######################################### -# * u0 = Im + .03*randn(size(Im)); % adding noise -# [Den] = LLT_model(single(u0), 10, 0.1, 1); -#Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); -#input, regularization_parameter , time_step, number_of_iterations, -# tolerance_constant, restrictive_Z_smoothing=0 - -del out2 -start_time = timeit.default_timer() -reg = Regularizer(Regularizer.Algorithm.LLT_model) -out2 = reg(input=u0, regularization_parameter=25, - time_step=0.0003, - tolerance_constant=0.001, - number_of_iterations=300) -txtstr = reg.printParametersToString() -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(2,3,4) - -# 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, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(out2,cmap="gray") - - -# ###################### PatchBased_Regul ######################################### -# # Quick 2D denoising example in Matlab: -# # Im = double(imread('lena_gray_256.tif'))/255; % loading image -# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -# # ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); -start_time = timeit.default_timer() -reg = Regularizer(Regularizer.Algorithm.PatchBased_Regul) -out2 = reg(input=u0, regularization_parameter=0.05, - searching_window_ratio=3, - similarity_window_ratio=1, - PB_filtering_parameter=0.08) -txtstr = reg.printParametersToString() -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) - -a=fig.add_subplot(2,3,5) - - -# 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, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(out2,cmap="gray") - - -# ###################### TGV_PD ######################################### -# # Quick 2D denoising example in Matlab: -# # Im = double(imread('lena_gray_256.tif'))/255; % loading image -# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -# # u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); - -start_time = timeit.default_timer() -reg = Regularizer(Regularizer.Algorithm.TGV_PD) -out2 = reg(input=u0, regularization_parameter=0.05, - first_order_term=1.3, - second_order_term=1, - number_of_iterations=550) -txtstr = reg.printParametersToString() -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(2,3,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, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(out2,cmap="gray") - - -plt.show() - -################################################################################ -## -## 3D Regularizers -## -################################################################################ -##Example: -## figure; -## Im = double(imread('lena_gray_256.tif'))/255; % loading image -## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; -## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); -# -##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha" -#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha" -# -#reader = vtk.vtkMetaImageReader() -#reader.SetFileName(os.path.normpath(filename)) -#reader.Update() -##vtk returns 3D images, let's take just the one slice there is as 2D -#Im = Converter.vtk2numpy(reader.GetOutput()) -#Im = Im.astype('float32') -##imgplot = plt.imshow(Im) -#perc = 0.05 -#u0 = Im + (perc* np.random.normal(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') -#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(), -# reader.GetOutput().GetOrigin()) -#converter.Update() -#writer = vtk.vtkMetaImageWriter() -#writer.SetInputData(converter.GetOutput()) -#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha") -##writer.Write() -# -# -### plot -#fig3D = plt.figure() -##a=fig.add_subplot(3,3,1) -##a.set_title('Original') -##imgplot = plt.imshow(Im) -#sliceNo = 32 -# -#a=fig3D.add_subplot(2,3,1) -#a.set_title('noise') -#imgplot = plt.imshow(u0.T[sliceNo]) -# -#reg_output3d = [] -# -############################################################################### -## Call regularizer -# -######################## SplitBregman_TV ##################################### -## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); -# -##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV) -# -##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30, -## #tolerance_constant=1e-4, -## TV_Penalty=Regularizer.TotalVariationPenalty.l1) -# -#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30, -# tolerance_constant=1e-4, -# TV_Penalty=Regularizer.TotalVariationPenalty.l1) -# -# -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## 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, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# -####################### FGP_TV ######################################### -## u = FGP_TV(single(u0), 0.05, 100, 1e-04); -#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005, -# number_of_iterations=200) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## 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, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# -####################### LLT_model ######################################### -## * u0 = Im + .03*randn(size(Im)); % adding noise -## [Den] = LLT_model(single(u0), 10, 0.1, 1); -##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); -##input, regularization_parameter , time_step, number_of_iterations, -## tolerance_constant, restrictive_Z_smoothing=0 -#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25, -# time_step=0.0003, -# tolerance_constant=0.0001, -# number_of_iterations=300) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## 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, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# -####################### PatchBased_Regul ######################################### -## Quick 2D denoising example in Matlab: -## Im = double(imread('lena_gray_256.tif'))/255; % loading image -## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); -# -#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05, -# searching_window_ratio=3, -# similarity_window_ratio=1, -# PB_filtering_parameter=0.08) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## 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, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# - -###################### TGV_PD ######################################### -# Quick 2D denoising example in Matlab: -# Im = double(imread('lena_gray_256.tif'))/255; % loading image -# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); - - -#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05, -# first_order_term=1.3, -# second_order_term=1, -# number_of_iterations=550) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## 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, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) diff --git a/Wrappers/Python/test/test_regularizers_3d.py b/Wrappers/Python/test/test_regularizers_3d.py deleted file mode 100644 index 2d11a7e..0000000 --- a/Wrappers/Python/test/test_regularizers_3d.py +++ /dev/null @@ -1,425 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Fri Aug 4 11:10:05 2017 - -@author: ofn77899 -""" - -#from ccpi.viewer.CILViewer2D import Converter -#import vtk - -import matplotlib.pyplot as plt -import numpy as np -import os -from enum import Enum -import timeit -#from PIL import Image -#from Regularizer import Regularizer -from ccpi.imaging.Regularizer import Regularizer - -############################################################################### -#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 -#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum -# intensity from the two images being compared, and respectively the same for -# minval. RMSE is given by the square root of MSE: -# sqrt[(sum(A - B) ** 2) / |A|], -# where |A| means the number of elements in A. By doing this, the maximum value -# given by RMSE is maxval. - -def nrmse(im1, im2): - a, b = im1.shape - rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b)) - max_val = max(np.max(im1), np.max(im2)) - min_val = min(np.min(im1), np.min(im2)) - return 1 - (rmse / (max_val - min_val)) -############################################################################### - -############################################################################### -# -# 2D Regularizers -# -############################################################################### -#Example: -# figure; -# Im = double(imread('lena_gray_256.tif'))/255; % loading image -# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; -# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); - - -#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif" -filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif" -#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif' - -#reader = vtk.vtkTIFFReader() -#reader.SetFileName(os.path.normpath(filename)) -#reader.Update() -Im = plt.imread(filename) -#Im = Image.open('/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif')/255 -#img.show() -Im = np.asarray(Im, dtype='float32') - -# create a 3D image by stacking N of this images - - -#imgplot = plt.imshow(Im) -perc = 0.05 -u_n = Im + (perc* np.random.normal(size=np.shape(Im))) -y,z = np.shape(u_n) -u_n = np.reshape(u_n , (1,y,z)) - -u0 = u_n.copy() -for i in range (19): - u_n = Im + (perc* np.random.normal(size=np.shape(Im))) - u_n = np.reshape(u_n , (1,y,z)) - - u0 = np.vstack ( (u0, u_n) ) - -# map the u0 u0->u0>0 -f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1) -u0 = f(u0).astype('float32') - -print ("Passed image shape {0}".format(np.shape(u0))) - -## plot -fig = plt.figure() -#a=fig.add_subplot(3,3,1) -#a.set_title('Original') -#imgplot = plt.imshow(Im) -sliceno = 10 - -a=fig.add_subplot(2,3,1) -a.set_title('noise') -imgplot = plt.imshow(u0[sliceno],cmap="gray") - -reg_output = [] -############################################################################## -# Call regularizer - -####################### SplitBregman_TV ##################################### -# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); - -use_object = True -if use_object: - reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV) - print (reg.pars) - reg.setParameter(input=u0) - reg.setParameter(regularization_parameter=10.) - # or - # reg.setParameter(input=u0, regularization_parameter=10., #number_of_iterations=30, - #tolerance_constant=1e-4, - #TV_Penalty=Regularizer.TotalVariationPenalty.l1) - plotme = reg() [0] - pars = reg.pars - textstr = reg.printParametersToString() - - #out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30, - #tolerance_constant=1e-4, - # TV_Penalty=Regularizer.TotalVariationPenalty.l1) - -#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30, -# tolerance_constant=1e-4, -# TV_Penalty=Regularizer.TotalVariationPenalty.l1) - -else: - out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10. ) - pars = out2[2] - reg_output.append(out2) - plotme = reg_output[-1][0] - textstr = out2[-1] - -a=fig.add_subplot(2,3,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, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(plotme[sliceno],cmap="gray") - -###################### FGP_TV ######################################### -# u = FGP_TV(single(u0), 0.05, 100, 1e-04); -out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.0005, - number_of_iterations=50) -pars = out2[-2] - -reg_output.append(out2) - -a=fig.add_subplot(2,3,3) - -textstr = out2[-1] - -# 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, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(reg_output[-1][0][sliceno]) -# place a text box in upper left in axes coords -a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray") - -###################### LLT_model ######################################### -# * u0 = Im + .03*randn(size(Im)); % adding noise -# [Den] = LLT_model(single(u0), 10, 0.1, 1); -#Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); -#input, regularization_parameter , time_step, number_of_iterations, -# tolerance_constant, restrictive_Z_smoothing=0 -out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25, - time_step=0.0003, - tolerance_constant=0.0001, - number_of_iterations=300) -pars = out2[-2] - -reg_output.append(out2) - -a=fig.add_subplot(2,3,4) - -textstr = out2[-1] - -# 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, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray") - - -# ###################### PatchBased_Regul ######################################### -# # Quick 2D denoising example in Matlab: -# # Im = double(imread('lena_gray_256.tif'))/255; % loading image -# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -# # ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); - -out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05, - searching_window_ratio=3, - similarity_window_ratio=1, - PB_filtering_parameter=0.08) -pars = out2[-2] -reg_output.append(out2) - -a=fig.add_subplot(2,3,5) - - -textstr = out2[-1] - -# 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, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray") - - -# ###################### TGV_PD ######################################### -# # Quick 2D denoising example in Matlab: -# # Im = double(imread('lena_gray_256.tif'))/255; % loading image -# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -# # u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); - - -out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05, - first_order_term=1.3, - second_order_term=1, - number_of_iterations=550) -pars = out2[-2] -reg_output.append(out2) - -a=fig.add_subplot(2,3,6) - - -textstr = out2[-1] - - -# 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, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray") - - -plt.show() - -################################################################################ -## -## 3D Regularizers -## -################################################################################ -##Example: -## figure; -## Im = double(imread('lena_gray_256.tif'))/255; % loading image -## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; -## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); -# -##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha" -#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha" -# -#reader = vtk.vtkMetaImageReader() -#reader.SetFileName(os.path.normpath(filename)) -#reader.Update() -##vtk returns 3D images, let's take just the one slice there is as 2D -#Im = Converter.vtk2numpy(reader.GetOutput()) -#Im = Im.astype('float32') -##imgplot = plt.imshow(Im) -#perc = 0.05 -#u0 = Im + (perc* np.random.normal(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') -#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(), -# reader.GetOutput().GetOrigin()) -#converter.Update() -#writer = vtk.vtkMetaImageWriter() -#writer.SetInputData(converter.GetOutput()) -#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha") -##writer.Write() -# -# -### plot -#fig3D = plt.figure() -##a=fig.add_subplot(3,3,1) -##a.set_title('Original') -##imgplot = plt.imshow(Im) -#sliceNo = 32 -# -#a=fig3D.add_subplot(2,3,1) -#a.set_title('noise') -#imgplot = plt.imshow(u0.T[sliceNo]) -# -#reg_output3d = [] -# -############################################################################### -## Call regularizer -# -######################## SplitBregman_TV ##################################### -## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); -# -##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV) -# -##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30, -## #tolerance_constant=1e-4, -## TV_Penalty=Regularizer.TotalVariationPenalty.l1) -# -#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30, -# tolerance_constant=1e-4, -# TV_Penalty=Regularizer.TotalVariationPenalty.l1) -# -# -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## 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, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# -####################### FGP_TV ######################################### -## u = FGP_TV(single(u0), 0.05, 100, 1e-04); -#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005, -# number_of_iterations=200) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## 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, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# -####################### LLT_model ######################################### -## * u0 = Im + .03*randn(size(Im)); % adding noise -## [Den] = LLT_model(single(u0), 10, 0.1, 1); -##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); -##input, regularization_parameter , time_step, number_of_iterations, -## tolerance_constant, restrictive_Z_smoothing=0 -#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25, -# time_step=0.0003, -# tolerance_constant=0.0001, -# number_of_iterations=300) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## 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, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# -####################### PatchBased_Regul ######################################### -## Quick 2D denoising example in Matlab: -## Im = double(imread('lena_gray_256.tif'))/255; % loading image -## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); -# -#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05, -# searching_window_ratio=3, -# similarity_window_ratio=1, -# PB_filtering_parameter=0.08) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## 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, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# - -###################### TGV_PD ######################################### -# Quick 2D denoising example in Matlab: -# Im = double(imread('lena_gray_256.tif'))/255; % loading image -# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); - - -#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05, -# first_order_term=1.3, -# second_order_term=1, -# number_of_iterations=550) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## 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, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -- cgit v1.2.3