summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--Wrappers/Python/ccpi/plugins/regularisers.py16
-rw-r--r--Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py26
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)