diff options
| -rw-r--r-- | Wrappers/Python/ccpi/plugins/regularisers.py | 16 | ||||
| -rw-r--r-- | Wrappers/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) |