diff options
author | jakobsj <jakobsj@users.noreply.github.com> | 2018-05-10 15:08:33 +0100 |
---|---|---|
committer | GitHub <noreply@github.com> | 2018-05-10 15:08:33 +0100 |
commit | 01ade10ddedd4da1c13fc4dcae4a0dea400550f0 (patch) | |
tree | df951b61d296b770412023ae1b98cd21383b0c7f /Wrappers | |
parent | 3c36f3b10e2b747eb5b1bf32c9320885a5494b56 (diff) | |
parent | 6bdd0ebce34b10ed299b986641830b664c140007 (diff) | |
download | framework-plugins-01ade10ddedd4da1c13fc4dcae4a0dea400550f0.tar.gz framework-plugins-01ade10ddedd4da1c13fc4dcae4a0dea400550f0.tar.bz2 framework-plugins-01ade10ddedd4da1c13fc4dcae4a0dea400550f0.tar.xz framework-plugins-01ade10ddedd4da1c13fc4dcae4a0dea400550f0.zip |
Merge pull request #6 from vais-ral/RGLTK_TV_denoising_demo
Add RGLTK TV denoising demo/test
Diffstat (limited to 'Wrappers')
-rw-r--r-- | Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py | 170 |
1 files changed, 170 insertions, 0 deletions
diff --git a/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py b/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py new file mode 100644 index 0000000..b8ccea9 --- /dev/null +++ b/Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py @@ -0,0 +1,170 @@ + +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.ops import LinearOperatorMatrix, Identity + +from ccpi.plugins.regularisers import _ROF_TV_, _FGP_TV_ + +# 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. + +# 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. +N = 64 +ig = ImageGeometry(voxel_num_x=N,voxel_num_y=N) +Phantom = ImageData(geometry=ig) + +x = Phantom.as_array() +x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 +x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 + +plt.imshow(x) +plt.title('Phantom image') +plt.show() + +# Identity operator for denoising +I = Identity() + +# Data and add noise +y = I.direct(Phantom) +np.random.seed(0) +y.array = y.array + 0.1*np.random.randn(N, N) + +plt.imshow(y.array) +plt.title('Noisy image') +plt.show() + +#%% TV parameter +lam_tv = 1.0 + +#%% Do CVX as high quality ground truth +if use_cvxpy: + # Compare to CVXPY + + # Construct the problem. + xtv_denoise = Variable(N,N) + objectivetv_denoise = Minimize(0.5*sum_squares(xtv_denoise - y.array) + lam_tv*tv(xtv_denoise) ) + probtv_denoise = Problem(objectivetv_denoise) + + # 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) + +plt.imshow(xtv_denoise.value) +plt.title('CVX TV') +plt.show() +print(objectivetv_denoise.value) + + +#%% THen FBPD + +# Data fidelity term +f_denoise = Norm2sq(I,y,c=0.5) + +# Initial guess +x_init_denoise = ImageData(np.zeros((N,N))) + +gtv = TV2D(lam_tv) +gtv(gtv.op.direct(x_init_denoise)) + +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:") +plt.imshow(x_fbpdtv_denoise.as_array()) +plt.title('FBPD TV') +plt.show() + +print(criterfbpdtv_denoise[-1]) + +#%% +plt.loglog([0,opt_tv['iter']], [objectivetv_denoise.value,objectivetv_denoise.value], label='CVX TV') +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') + +xtv_rof = g_rof.prox(y,1.0) + +print("CCPi-RGL TV ROF:") +plt.imshow(xtv_rof.as_array()) +plt.title('ROF TV prox') +plt.show() +print(g_rof(xtv_rof)) + +#%% 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') + +xtv_fgp = g_fgp.prox(y,1.0) + +print("CCPi-RGL TV FGP:") +plt.imshow(xtv_fgp.as_array()) +plt.title('FGP TV prox') +plt.show() +print(g_fgp(xtv_fgp)) + +# Compare all reconstruction +clims = (-0.2,1.2) +dlims = (-0.2,0.2) +cols = 3 +rows = 2 +current = 1 + +fig = plt.figure() +a=fig.add_subplot(rows,cols,current) +a.set_title('FBPD') +imgplot = plt.imshow(x_fbpdtv_denoise.as_array(),vmin=clims[0],vmax=clims[1]) +plt.axis('off') + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('ROF') +imgplot = plt.imshow(xtv_rof.as_array(),vmin=clims[0],vmax=clims[1]) +plt.axis('off') + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +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('FISTA LS+1') +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('FBPD LS+1') +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('FBPD TV') +imgplot = plt.imshow(xtv_fgp.as_array()-xtv_denoise.value,vmin=dlims[0],vmax=dlims[1]) +plt.axis('off') |