From 0dd1cadcfead9a2a5f225e1500c97cc00a8068d6 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Fri, 11 May 2018 12:44:18 +0100 Subject: some fixes regarding denoising --- Wrappers/Python/ccpi/plugins/regularisers.py | 16 +++++++------ .../Python/wip/demo_compare_RGLTK_TV_denoising.py | 26 +++++++++++++--------- 2 files changed, 25 insertions(+), 17 deletions(-) diff --git a/Wrappers/Python/ccpi/plugins/regularisers.py b/Wrappers/Python/ccpi/plugins/regularisers.py index e9c88a4..6d865cc 100644 --- a/Wrappers/Python/ccpi/plugins/regularisers.py +++ b/Wrappers/Python/ccpi/plugins/regularisers.py @@ -25,7 +25,6 @@ from ccpi.optimisation.ops import Operator import numpy as np - class _ROF_TV_(Operator): def __init__(self,lambdaReg,iterationsTV,tolerance,time_marchstep,device): # set parameters @@ -33,9 +32,10 @@ class _ROF_TV_(Operator): self.iterationsTV = iterationsTV self.time_marchstep = time_marchstep self.device = device # string for 'cpu' or 'gpu' - def __call__(self,x): + def __call__(self,x,x1,typeEnergy): # evaluate objective function of TV gradient - EnergyValTV = TV_ENERGY(np.asarray(x.as_array(), dtype=np.float32), np.asarray(x.as_array(), dtype=np.float32), self.lambdaReg, 2) + # typeEnergy is either 1 (LS + TV for denoising) or 2 (just TV fidelity) + EnergyValTV = TV_ENERGY(np.asarray(x.as_array(), dtype=np.float32), np.asarray(x1.as_array(), dtype=np.float32), self.lambdaReg, typeEnergy) return EnergyValTV def prox(self,x,Lipshitz): pars = {'algorithm' : ROF_TV, \ @@ -60,9 +60,10 @@ class _FGP_TV_(Operator): self.nonnegativity = nonnegativity self.printing = printing self.device = device # string for 'cpu' or 'gpu' - def __call__(self,x): + def __call__(self,x,x1,typeEnergy): # evaluate objective function of TV gradient - EnergyValTV = TV_ENERGY(np.asarray(x.as_array(), dtype=np.float32), np.asarray(x.as_array(), dtype=np.float32), self.lambdaReg, 2) + # typeEnergy is either 1 (LS + TV for denoising) or 2 (just TV fidelity) + EnergyValTV = TV_ENERGY(np.asarray(x.as_array(), dtype=np.float32), np.asarray(x1.as_array(), dtype=np.float32), self.lambdaReg, typeEnergy) return EnergyValTV def prox(self,x,Lipshitz): pars = {'algorithm' : FGP_TV, \ @@ -93,9 +94,10 @@ class _SB_TV_(Operator): self.methodTV = methodTV self.printing = printing self.device = device # string for 'cpu' or 'gpu' - def __call__(self,x): + def __call__(self,x,typeEnergy): # evaluate objective function of TV gradient - EnergyValTV = TV_ENERGY(np.asarray(x.as_array(), dtype=np.float32), np.asarray(x.as_array(), dtype=np.float32), self.lambdaReg, 2) + # typeEnergy is either 1 (LS + TV for denoising) or 2 (just TV fidelity) + EnergyValTV = TV_ENERGY(np.asarray(x.as_array(), dtype=np.float32), np.asarray(x.as_array(), dtype=np.float32), self.lambdaReg, typeEnergy) return EnergyValTV def prox(self,x,Lipshitz): pars = {'algorithm' : SB_TV, \ diff --git a/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py b/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py index 559679e..dd9044e 100644 --- a/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py +++ b/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py @@ -1,24 +1,26 @@ from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, DataContainer from ccpi.optimisation.algs import FISTA, FBPD, CGLS -from ccpi.optimisation.funcs import Norm2sq, ZeroFun, Norm1, TV2D +from ccpi.optimisation.funcs import Norm2sq, ZeroFun, Norm1, TV2D, Identity -from ccpi.optimisation.ops import LinearOperatorMatrix, Identity +from ccpi.optimisation.ops import LinearOperatorMatrix from ccpi.plugins.regularisers import _ROF_TV_, _FGP_TV_ +import numpy as np +import matplotlib.pyplot as plt + + +#%% # Requires CVXPY, see http://www.cvxpy.org/ # CVXPY can be installed in anaconda using # conda install -c cvxgrp cvxpy libgcc - # Whether to use or omit CVXPY use_cvxpy = True if use_cvxpy: from cvxpy import * -import numpy as np -import matplotlib.pyplot as plt - +#%% # Now try 1-norm and TV denoising with FBPD, first 1-norm. @@ -93,6 +95,7 @@ x_fbpdtv_denoise, itfbpdtv_denoise, timingfbpdtv_denoise, criterfbpdtv_denoise = print("CVXPY least squares plus TV solution and objective value:") +plt.figure() plt.imshow(x_fbpdtv_denoise.as_array()) plt.title('FBPD TV') plt.show() @@ -105,27 +108,30 @@ plt.loglog(criterfbpdtv_denoise, label='FBPD TV') plt.show() #%% FISTA with ROF-TV regularisation -g_rof = _ROF_TV_(lambdaReg = lam_tv,iterationsTV=5000,tolerance=1e-5,time_marchstep=0.01,device='cpu') +g_rof = _ROF_TV_(lambdaReg = lam_tv,iterationsTV=5000,tolerance=0,time_marchstep=0.001,device='cpu') xtv_rof = g_rof.prox(y,1.0) print("CCPi-RGL TV ROF:") +plt.figure() plt.imshow(xtv_rof.as_array()) plt.title('ROF TV prox') plt.show() -print(g_rof(xtv_rof)) +print(g_rof(xtv_rof,y,typeEnergy=1)) #%% FISTA with FGP-TV regularisation -g_fgp = _FGP_TV_(lambdaReg = lam_tv,iterationsTV=5000,tolerance=1e-5,methodTV=0,nonnegativity=0,printing=0,device='cpu') +g_fgp = _FGP_TV_(lambdaReg = lam_tv,iterationsTV=5000,tolerance=1e-7,methodTV=0,nonnegativity=0,printing=0,device='cpu') xtv_fgp = g_fgp.prox(y,1.0) print("CCPi-RGL TV FGP:") +plt.figure() plt.imshow(xtv_fgp.as_array()) plt.title('FGP TV prox') plt.show() -print(g_fgp(xtv_fgp)) +print(g_fgp(xtv_fgp,y,typeEnergy=1)) +#%% # Compare all reconstruction clims = (-0.2,1.2) dlims = (-0.2,0.2) -- cgit v1.2.1 From cfd675c4dcc9ca090e6401f70c68bdc34da004db Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Fri, 11 May 2018 15:42:35 +0100 Subject: further corrections, SB-TV added --- Wrappers/Python/ccpi/plugins/regularisers.py | 12 ++++---- .../Python/wip/demo_compare_RGLTK_TV_denoising.py | 35 +++++++++++++++------- 2 files changed, 31 insertions(+), 16 deletions(-) diff --git a/Wrappers/Python/ccpi/plugins/regularisers.py b/Wrappers/Python/ccpi/plugins/regularisers.py index 6d865cc..9f4d3fc 100644 --- a/Wrappers/Python/ccpi/plugins/regularisers.py +++ b/Wrappers/Python/ccpi/plugins/regularisers.py @@ -32,10 +32,10 @@ class _ROF_TV_(Operator): self.iterationsTV = iterationsTV self.time_marchstep = time_marchstep self.device = device # string for 'cpu' or 'gpu' - def __call__(self,x,x1,typeEnergy): + def __call__(self,x): # evaluate objective function of TV gradient # typeEnergy is either 1 (LS + TV for denoising) or 2 (just TV fidelity) - EnergyValTV = TV_ENERGY(np.asarray(x.as_array(), dtype=np.float32), np.asarray(x1.as_array(), dtype=np.float32), self.lambdaReg, typeEnergy) + EnergyValTV = TV_ENERGY(np.asarray(x.as_array(), dtype=np.float32), np.asarray(x.as_array(), dtype=np.float32), self.lambdaReg, 2) return EnergyValTV def prox(self,x,Lipshitz): pars = {'algorithm' : ROF_TV, \ @@ -60,10 +60,10 @@ class _FGP_TV_(Operator): self.nonnegativity = nonnegativity self.printing = printing self.device = device # string for 'cpu' or 'gpu' - def __call__(self,x,x1,typeEnergy): + def __call__(self,x): # evaluate objective function of TV gradient # typeEnergy is either 1 (LS + TV for denoising) or 2 (just TV fidelity) - EnergyValTV = TV_ENERGY(np.asarray(x.as_array(), dtype=np.float32), np.asarray(x1.as_array(), dtype=np.float32), self.lambdaReg, typeEnergy) + EnergyValTV = TV_ENERGY(np.asarray(x.as_array(), dtype=np.float32), np.asarray(x.as_array(), dtype=np.float32), self.lambdaReg, 2) return EnergyValTV def prox(self,x,Lipshitz): pars = {'algorithm' : FGP_TV, \ @@ -94,10 +94,10 @@ class _SB_TV_(Operator): self.methodTV = methodTV self.printing = printing self.device = device # string for 'cpu' or 'gpu' - def __call__(self,x,typeEnergy): + def __call__(self,x): # evaluate objective function of TV gradient # typeEnergy is either 1 (LS + TV for denoising) or 2 (just TV fidelity) - EnergyValTV = TV_ENERGY(np.asarray(x.as_array(), dtype=np.float32), np.asarray(x.as_array(), dtype=np.float32), self.lambdaReg, typeEnergy) + EnergyValTV = TV_ENERGY(np.asarray(x.as_array(), dtype=np.float32), np.asarray(x.as_array(), dtype=np.float32), self.lambdaReg, 2) return EnergyValTV def prox(self,x,Lipshitz): pars = {'algorithm' : SB_TV, \ diff --git a/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py b/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py index dd9044e..07f4090 100644 --- a/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py +++ b/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py @@ -5,7 +5,7 @@ from ccpi.optimisation.funcs import Norm2sq, ZeroFun, Norm1, TV2D, Identity from ccpi.optimisation.ops import LinearOperatorMatrix -from ccpi.plugins.regularisers import _ROF_TV_, _FGP_TV_ +from ccpi.plugins.regularisers import _ROF_TV_, _FGP_TV_, _SB_TV_ import numpy as np import matplotlib.pyplot as plt @@ -72,12 +72,12 @@ if use_cvxpy: # print(xtv_denoise.value) # print(objectivetv_denoise.value) +plt.figure() plt.imshow(xtv_denoise.value) -plt.title('CVX TV') +plt.title('CVX TV with objective equal to {:.2f}'.format(objectivetv_denoise.value)) plt.show() print(objectivetv_denoise.value) - #%% THen FBPD # Data fidelity term @@ -97,7 +97,7 @@ x_fbpdtv_denoise, itfbpdtv_denoise, timingfbpdtv_denoise, criterfbpdtv_denoise = print("CVXPY least squares plus TV solution and objective value:") plt.figure() plt.imshow(x_fbpdtv_denoise.as_array()) -plt.title('FBPD TV') +plt.title('FBPD TV with objective equal to {:.2f}'.format(criterfbpdtv_denoise[-1])) plt.show() print(criterfbpdtv_denoise[-1]) @@ -108,30 +108,45 @@ plt.loglog(criterfbpdtv_denoise, label='FBPD TV') plt.show() #%% FISTA with ROF-TV regularisation -g_rof = _ROF_TV_(lambdaReg = lam_tv,iterationsTV=5000,tolerance=0,time_marchstep=0.001,device='cpu') +g_rof = _ROF_TV_(lambdaReg = lam_tv,iterationsTV=5000,tolerance=0,time_marchstep=0.0009,device='cpu') xtv_rof = g_rof.prox(y,1.0) print("CCPi-RGL TV ROF:") plt.figure() plt.imshow(xtv_rof.as_array()) -plt.title('ROF TV prox') +valObjRof = g_rof(xtv_rof) +plt.title('ROF TV prox with objective equal to {:.2f}'.format(valObjRof[0])) plt.show() -print(g_rof(xtv_rof,y,typeEnergy=1)) +print(valObjRof[0]) #%% FISTA with FGP-TV regularisation -g_fgp = _FGP_TV_(lambdaReg = lam_tv,iterationsTV=5000,tolerance=1e-7,methodTV=0,nonnegativity=0,printing=0,device='cpu') +g_fgp = _FGP_TV_(lambdaReg = lam_tv,iterationsTV=5000,tolerance=0,methodTV=0,nonnegativity=0,printing=0,device='cpu') xtv_fgp = g_fgp.prox(y,1.0) print("CCPi-RGL TV FGP:") plt.figure() plt.imshow(xtv_fgp.as_array()) -plt.title('FGP TV prox') +valObjFGP = g_fgp(xtv_fgp) +plt.title('FGP TV prox with objective equal to {:.2f}'.format(valObjFGP[0])) plt.show() -print(g_fgp(xtv_fgp,y,typeEnergy=1)) +print(valObjFGP[0]) +#%% Split-Bregman-TV regularisation +g_sb = _SB_TV_(lambdaReg = lam_tv,iterationsTV=150,tolerance=0,methodTV=0,printing=0,device='cpu') +xtv_sb = g_sb.prox(y,1.0) + +print("CCPi-RGL TV SB:") +plt.figure() +plt.imshow(xtv_sb.as_array()) +valObjSB = g_sb(xtv_sb) +plt.title('SB TV prox with objective equal to {:.2f}'.format(valObjSB[0])) +plt.show() +print(valObjSB[0]) #%% + + # Compare all reconstruction clims = (-0.2,1.2) dlims = (-0.2,0.2) -- cgit v1.2.1 From 3d939a6139e664c3f8143031d0aaf765298efda5 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Sat, 12 May 2018 17:30:54 +0100 Subject: fixed objective --- .../Python/wip/demo_compare_RGLTK_TV_denoising.py | 29 ++++++++++++---------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py b/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py index 07f4090..0d57e5e 100644 --- a/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py +++ b/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py @@ -1,16 +1,13 @@ from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, DataContainer from ccpi.optimisation.algs import FISTA, FBPD, CGLS -from ccpi.optimisation.funcs import Norm2sq, ZeroFun, Norm1, TV2D, Identity - -from ccpi.optimisation.ops import LinearOperatorMatrix +from ccpi.optimisation.funcs import Norm2sq, ZeroFun, Norm1, TV2D +from ccpi.optimisation.ops import LinearOperatorMatrix, Identity from ccpi.plugins.regularisers import _ROF_TV_, _FGP_TV_, _SB_TV_ import numpy as np import matplotlib.pyplot as plt - - #%% # Requires CVXPY, see http://www.cvxpy.org/ # CVXPY can be installed in anaconda using @@ -108,7 +105,7 @@ plt.loglog(criterfbpdtv_denoise, label='FBPD TV') plt.show() #%% FISTA with ROF-TV regularisation -g_rof = _ROF_TV_(lambdaReg = lam_tv,iterationsTV=5000,tolerance=0,time_marchstep=0.0009,device='cpu') +g_rof = _ROF_TV_(lambdaReg = lam_tv,iterationsTV=2000,tolerance=0,time_marchstep=0.0009,device='cpu') xtv_rof = g_rof.prox(y,1.0) @@ -116,9 +113,11 @@ print("CCPi-RGL TV ROF:") plt.figure() plt.imshow(xtv_rof.as_array()) valObjRof = g_rof(xtv_rof) -plt.title('ROF TV prox with objective equal to {:.2f}'.format(valObjRof[0])) +data_energy = 0.5*np.sum(np.power((xtv_rof.as_array() - y.array),2)) +EnergytotalROF = data_energy + 0.5*valObjRof[0] +plt.title('ROF TV prox with objective equal to {:.2f}'.format(EnergytotalROF)) plt.show() -print(valObjRof[0]) +print(EnergytotalROF) #%% FISTA with FGP-TV regularisation g_fgp = _FGP_TV_(lambdaReg = lam_tv,iterationsTV=5000,tolerance=0,methodTV=0,nonnegativity=0,printing=0,device='cpu') @@ -129,11 +128,13 @@ print("CCPi-RGL TV FGP:") plt.figure() plt.imshow(xtv_fgp.as_array()) valObjFGP = g_fgp(xtv_fgp) -plt.title('FGP TV prox with objective equal to {:.2f}'.format(valObjFGP[0])) +data_energy = 0.5*np.sum(np.power((xtv_fgp.as_array() - y.array),2)) +EnergytotalFGP = data_energy + 0.5*valObjFGP[0] +plt.title('FGP TV prox with objective equal to {:.2f}'.format(EnergytotalFGP)) plt.show() -print(valObjFGP[0]) +print(EnergytotalFGP) #%% Split-Bregman-TV regularisation -g_sb = _SB_TV_(lambdaReg = lam_tv,iterationsTV=150,tolerance=0,methodTV=0,printing=0,device='cpu') +g_sb = _SB_TV_(lambdaReg = lam_tv,iterationsTV=1000,tolerance=0,methodTV=0,printing=0,device='cpu') xtv_sb = g_sb.prox(y,1.0) @@ -141,9 +142,11 @@ print("CCPi-RGL TV SB:") plt.figure() plt.imshow(xtv_sb.as_array()) valObjSB = g_sb(xtv_sb) -plt.title('SB TV prox with objective equal to {:.2f}'.format(valObjSB[0])) +data_energy = 0.5*np.sum(np.power((xtv_sb.as_array() - y.array),2)) +EnergytotalSB = data_energy + 0.5*valObjSB[0] +plt.title('SB TV prox with objective equal to {:.2f}'.format(EnergytotalSB)) plt.show() -print(valObjSB[0]) +print(EnergytotalSB) #%% -- cgit v1.2.1 From d1875172687fc854df35fa9bfc6ac07a148d7f18 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Sat, 12 May 2018 19:03:26 +0100 Subject: fixed objective2 --- Wrappers/Python/ccpi/plugins/regularisers.py | 9 +++------ Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py | 18 +++++++----------- 2 files changed, 10 insertions(+), 17 deletions(-) diff --git a/Wrappers/Python/ccpi/plugins/regularisers.py b/Wrappers/Python/ccpi/plugins/regularisers.py index 9f4d3fc..46464a9 100644 --- a/Wrappers/Python/ccpi/plugins/regularisers.py +++ b/Wrappers/Python/ccpi/plugins/regularisers.py @@ -34,9 +34,8 @@ class _ROF_TV_(Operator): self.device = device # string for 'cpu' or 'gpu' def __call__(self,x): # evaluate objective function of TV gradient - # typeEnergy is either 1 (LS + TV for denoising) or 2 (just TV fidelity) EnergyValTV = TV_ENERGY(np.asarray(x.as_array(), dtype=np.float32), np.asarray(x.as_array(), dtype=np.float32), self.lambdaReg, 2) - return EnergyValTV + return 0.5*EnergyValTV[0] def prox(self,x,Lipshitz): pars = {'algorithm' : ROF_TV, \ 'input' : np.asarray(x.as_array(), dtype=np.float32),\ @@ -62,9 +61,8 @@ class _FGP_TV_(Operator): self.device = device # string for 'cpu' or 'gpu' def __call__(self,x): # evaluate objective function of TV gradient - # typeEnergy is either 1 (LS + TV for denoising) or 2 (just TV fidelity) EnergyValTV = TV_ENERGY(np.asarray(x.as_array(), dtype=np.float32), np.asarray(x.as_array(), dtype=np.float32), self.lambdaReg, 2) - return EnergyValTV + return 0.5*EnergyValTV[0] def prox(self,x,Lipshitz): pars = {'algorithm' : FGP_TV, \ 'input' : np.asarray(x.as_array(), dtype=np.float32),\ @@ -96,9 +94,8 @@ class _SB_TV_(Operator): self.device = device # string for 'cpu' or 'gpu' def __call__(self,x): # evaluate objective function of TV gradient - # typeEnergy is either 1 (LS + TV for denoising) or 2 (just TV fidelity) EnergyValTV = TV_ENERGY(np.asarray(x.as_array(), dtype=np.float32), np.asarray(x.as_array(), dtype=np.float32), self.lambdaReg, 2) - return EnergyValTV + return 0.5*EnergyValTV[0] def prox(self,x,Lipshitz): pars = {'algorithm' : SB_TV, \ 'input' : np.asarray(x.as_array(), dtype=np.float32),\ diff --git a/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py b/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py index 0d57e5e..2bf1286 100644 --- a/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py +++ b/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py @@ -75,11 +75,13 @@ plt.title('CVX TV with objective equal to {:.2f}'.format(objectivetv_denoise.va plt.show() print(objectivetv_denoise.value) -#%% THen FBPD - +#%% # Data fidelity term f_denoise = Norm2sq(I,y,c=0.5) +#%% + +#%% THen FBPD # Initial guess x_init_denoise = ImageData(np.zeros((N,N))) @@ -112,9 +114,7 @@ xtv_rof = g_rof.prox(y,1.0) print("CCPi-RGL TV ROF:") plt.figure() plt.imshow(xtv_rof.as_array()) -valObjRof = g_rof(xtv_rof) -data_energy = 0.5*np.sum(np.power((xtv_rof.as_array() - y.array),2)) -EnergytotalROF = data_energy + 0.5*valObjRof[0] +EnergytotalROF = f_denoise(xtv_rof) + g_rof(xtv_rof) plt.title('ROF TV prox with objective equal to {:.2f}'.format(EnergytotalROF)) plt.show() print(EnergytotalROF) @@ -127,9 +127,7 @@ xtv_fgp = g_fgp.prox(y,1.0) print("CCPi-RGL TV FGP:") plt.figure() plt.imshow(xtv_fgp.as_array()) -valObjFGP = g_fgp(xtv_fgp) -data_energy = 0.5*np.sum(np.power((xtv_fgp.as_array() - y.array),2)) -EnergytotalFGP = data_energy + 0.5*valObjFGP[0] +EnergytotalFGP = f_denoise(xtv_fgp) + g_fgp(xtv_fgp) plt.title('FGP TV prox with objective equal to {:.2f}'.format(EnergytotalFGP)) plt.show() print(EnergytotalFGP) @@ -141,9 +139,7 @@ xtv_sb = g_sb.prox(y,1.0) print("CCPi-RGL TV SB:") plt.figure() plt.imshow(xtv_sb.as_array()) -valObjSB = g_sb(xtv_sb) -data_energy = 0.5*np.sum(np.power((xtv_sb.as_array() - y.array),2)) -EnergytotalSB = data_energy + 0.5*valObjSB[0] +EnergytotalSB = f_denoise(xtv_sb) + g_fgp(xtv_sb) plt.title('SB TV prox with objective equal to {:.2f}'.format(EnergytotalSB)) plt.show() print(EnergytotalSB) -- cgit v1.2.1 From f5aa0dd2976df492374b35d9f6bfcc0933803ebb Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Sat, 12 May 2018 21:07:29 +0100 Subject: Add SB in final display of all. --- Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py b/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py index 2bf1286..fd992de 100644 --- a/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py +++ b/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py @@ -149,7 +149,7 @@ print(EnergytotalSB) # Compare all reconstruction clims = (-0.2,1.2) dlims = (-0.2,0.2) -cols = 3 +cols = 4 rows = 2 current = 1 @@ -171,6 +171,12 @@ a.set_title('FGP') imgplot = plt.imshow(xtv_fgp.as_array(),vmin=clims[0],vmax=clims[1]) plt.axis('off') +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('SB') +imgplot = plt.imshow(xtv_sb.as_array(),vmin=clims[0],vmax=clims[1]) +plt.axis('off') + current = current + 1 a=fig.add_subplot(rows,cols,current) a.set_title('FBPD - CVX') @@ -179,12 +185,18 @@ plt.axis('off') current = current + 1 a=fig.add_subplot(rows,cols,current) -a.set_title('ROF - TV') +a.set_title('ROF - CVX') imgplot = plt.imshow(xtv_rof.as_array()-xtv_denoise.value,vmin=dlims[0],vmax=dlims[1]) plt.axis('off') current = current + 1 a=fig.add_subplot(rows,cols,current) -a.set_title('FGP - TV') +a.set_title('FGP - CVX') imgplot = plt.imshow(xtv_fgp.as_array()-xtv_denoise.value,vmin=dlims[0],vmax=dlims[1]) plt.axis('off') + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('SB - CVX') +imgplot = plt.imshow(xtv_sb.as_array()-xtv_denoise.value,vmin=dlims[0],vmax=dlims[1]) +plt.axis('off') -- cgit v1.2.1 From 653e0adc87c255a392f30590af446fc78043e194 Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Sat, 12 May 2018 21:22:43 +0100 Subject: Tidy up denoising demo. --- .../Python/wip/demo_compare_RGLTK_TV_denoising.py | 133 +++++++++++++-------- 1 file changed, 83 insertions(+), 50 deletions(-) diff --git a/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py b/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py index fd992de..911cff4 100644 --- a/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py +++ b/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py @@ -1,4 +1,10 @@ +# This demo illustrates how the CCPi Regularisation Toolkit can be used +# as TV denoising for use with the FISTA algorithm of the modular +# optimisation framework and compares with the FBPD TV implementation as well +# as CVXPY. + +# All own imports from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, DataContainer from ccpi.optimisation.algs import FISTA, FBPD, CGLS from ccpi.optimisation.funcs import Norm2sq, ZeroFun, Norm1, TV2D @@ -6,12 +12,15 @@ from ccpi.optimisation.ops import LinearOperatorMatrix, Identity from ccpi.plugins.regularisers import _ROF_TV_, _FGP_TV_, _SB_TV_ +# All external imports import numpy as np import matplotlib.pyplot as plt + #%% # Requires CVXPY, see http://www.cvxpy.org/ # CVXPY can be installed in anaconda using # conda install -c cvxgrp cvxpy libgcc + # Whether to use or omit CVXPY use_cvxpy = True if use_cvxpy: @@ -19,8 +28,6 @@ if use_cvxpy: #%% -# Now try 1-norm and TV denoising with FBPD, first 1-norm. - # Set up phantom size NxN by creating ImageGeometry, initialising the # ImageData object with this geometry and empty array and finally put some # data into its array, and display as image. @@ -44,6 +51,7 @@ y = I.direct(Phantom) np.random.seed(0) y.array = y.array + 0.1*np.random.randn(N, N) +# Display noisy image plt.imshow(y.array) plt.title('Noisy image') plt.show() @@ -51,9 +59,8 @@ plt.show() #%% TV parameter lam_tv = 1.0 -#%% Do CVX as high quality ground truth +#%% Do CVX as high quality ground truth for comparison. if use_cvxpy: - # Compare to CVXPY # Construct the problem. xtv_denoise = Variable(N,N) @@ -63,17 +70,15 @@ if use_cvxpy: # The optimal objective is returned by prob.solve(). resulttv_denoise = probtv_denoise.solve(verbose=False,solver=SCS,eps=1e-12) - # The optimal solution for x is stored in x.value and optimal objective value - # is in result as well as in objective.value - print("CVXPY least squares plus TV solution and objective value:") - # print(xtv_denoise.value) - # print(objectivetv_denoise.value) + # The optimal solution for x is stored in x.value and optimal objective + # value is in result as well as in objective.value -plt.figure() -plt.imshow(xtv_denoise.value) -plt.title('CVX TV with objective equal to {:.2f}'.format(objectivetv_denoise.value)) -plt.show() -print(objectivetv_denoise.value) + # Display + plt.figure() + plt.imshow(xtv_denoise.value) + plt.title('CVX TV with objective equal to {:.2f}'.format(objectivetv_denoise.value)) + plt.show() + print(objectivetv_denoise.value) #%% # Data fidelity term @@ -81,19 +86,22 @@ f_denoise = Norm2sq(I,y,c=0.5) #%% -#%% THen FBPD +#%% Then run FBPD algorithm for TV denoising + # Initial guess x_init_denoise = ImageData(np.zeros((N,N))) +# Set up TV function gtv = TV2D(lam_tv) -gtv(gtv.op.direct(x_init_denoise)) -opt_tv = {'tol': 1e-4, 'iter': 10000} +# Evalutate TV of noisy image. +gtv(gtv.op.direct(y)) +# Specify FBPD options and run FBPD. +opt_tv = {'tol': 1e-4, 'iter': 10000} x_fbpdtv_denoise, itfbpdtv_denoise, timingfbpdtv_denoise, criterfbpdtv_denoise = FBPD(x_init_denoise, None, f_denoise, gtv,opt=opt_tv) - -print("CVXPY least squares plus TV solution and objective value:") +print("FBPD least squares plus TV solution and objective value:") plt.figure() plt.imshow(x_fbpdtv_denoise.as_array()) plt.title('FBPD TV with objective equal to {:.2f}'.format(criterfbpdtv_denoise[-1])) @@ -101,16 +109,24 @@ plt.show() print(criterfbpdtv_denoise[-1]) -#%% -plt.loglog([0,opt_tv['iter']], [objectivetv_denoise.value,objectivetv_denoise.value], label='CVX TV') +# Also plot history of criterion vs. CVX +if use_cvxpy: + plt.loglog([0,opt_tv['iter']], [objectivetv_denoise.value,objectivetv_denoise.value], label='CVX TV') plt.loglog(criterfbpdtv_denoise, label='FBPD TV') +plt.legend() plt.show() #%% FISTA with ROF-TV regularisation -g_rof = _ROF_TV_(lambdaReg = lam_tv,iterationsTV=2000,tolerance=0,time_marchstep=0.0009,device='cpu') +g_rof = _ROF_TV_(lambdaReg = lam_tv, + iterationsTV=2000, + tolerance=0, + time_marchstep=0.0009, + device='cpu') +# Evaluating the proximal operator corresponds to denoising. xtv_rof = g_rof.prox(y,1.0) +# Display denoised image and final criterion value. print("CCPi-RGL TV ROF:") plt.figure() plt.imshow(xtv_rof.as_array()) @@ -120,10 +136,18 @@ plt.show() print(EnergytotalROF) #%% FISTA with FGP-TV regularisation -g_fgp = _FGP_TV_(lambdaReg = lam_tv,iterationsTV=5000,tolerance=0,methodTV=0,nonnegativity=0,printing=0,device='cpu') - +g_fgp = _FGP_TV_(lambdaReg = lam_tv, + iterationsTV=5000, + tolerance=0, + methodTV=0, + nonnegativity=0, + printing=0, + device='cpu') + +# Evaluating the proximal operator corresponds to denoising. xtv_fgp = g_fgp.prox(y,1.0) +# Display denoised image and final criterion value. print("CCPi-RGL TV FGP:") plt.figure() plt.imshow(xtv_fgp.as_array()) @@ -131,11 +155,19 @@ EnergytotalFGP = f_denoise(xtv_fgp) + g_fgp(xtv_fgp) plt.title('FGP TV prox with objective equal to {:.2f}'.format(EnergytotalFGP)) plt.show() print(EnergytotalFGP) -#%% Split-Bregman-TV regularisation -g_sb = _SB_TV_(lambdaReg = lam_tv,iterationsTV=1000,tolerance=0,methodTV=0,printing=0,device='cpu') +#%% Split-Bregman-TV regularisation +g_sb = _SB_TV_(lambdaReg = lam_tv, + iterationsTV=1000, + tolerance=0, + methodTV=0, + printing=0, + device='cpu') + +# Evaluating the proximal operator corresponds to denoising. xtv_sb = g_sb.prox(y,1.0) +# Display denoised image and final criterion value. print("CCPi-RGL TV SB:") plt.figure() plt.imshow(xtv_sb.as_array()) @@ -143,8 +175,8 @@ EnergytotalSB = f_denoise(xtv_sb) + g_fgp(xtv_sb) plt.title('SB TV prox with objective equal to {:.2f}'.format(EnergytotalSB)) plt.show() print(EnergytotalSB) -#%% +#%% # Compare all reconstruction clims = (-0.2,1.2) @@ -177,26 +209,27 @@ a.set_title('SB') imgplot = plt.imshow(xtv_sb.as_array(),vmin=clims[0],vmax=clims[1]) plt.axis('off') -current = current + 1 -a=fig.add_subplot(rows,cols,current) -a.set_title('FBPD - CVX') -imgplot = plt.imshow(x_fbpdtv_denoise.as_array()-xtv_denoise.value,vmin=dlims[0],vmax=dlims[1]) -plt.axis('off') - -current = current + 1 -a=fig.add_subplot(rows,cols,current) -a.set_title('ROF - CVX') -imgplot = plt.imshow(xtv_rof.as_array()-xtv_denoise.value,vmin=dlims[0],vmax=dlims[1]) -plt.axis('off') - -current = current + 1 -a=fig.add_subplot(rows,cols,current) -a.set_title('FGP - CVX') -imgplot = plt.imshow(xtv_fgp.as_array()-xtv_denoise.value,vmin=dlims[0],vmax=dlims[1]) -plt.axis('off') - -current = current + 1 -a=fig.add_subplot(rows,cols,current) -a.set_title('SB - CVX') -imgplot = plt.imshow(xtv_sb.as_array()-xtv_denoise.value,vmin=dlims[0],vmax=dlims[1]) -plt.axis('off') +if use_cvxpy: + current = current + 1 + a=fig.add_subplot(rows,cols,current) + a.set_title('FBPD - CVX') + imgplot = plt.imshow(x_fbpdtv_denoise.as_array()-xtv_denoise.value,vmin=dlims[0],vmax=dlims[1]) + plt.axis('off') + + current = current + 1 + a=fig.add_subplot(rows,cols,current) + a.set_title('ROF - CVX') + imgplot = plt.imshow(xtv_rof.as_array()-xtv_denoise.value,vmin=dlims[0],vmax=dlims[1]) + plt.axis('off') + + current = current + 1 + a=fig.add_subplot(rows,cols,current) + a.set_title('FGP - CVX') + imgplot = plt.imshow(xtv_fgp.as_array()-xtv_denoise.value,vmin=dlims[0],vmax=dlims[1]) + plt.axis('off') + + current = current + 1 + a=fig.add_subplot(rows,cols,current) + a.set_title('SB - CVX') + imgplot = plt.imshow(xtv_sb.as_array()-xtv_denoise.value,vmin=dlims[0],vmax=dlims[1]) + plt.axis('off') -- cgit v1.2.1