diff options
| -rw-r--r-- | Wrappers/Python/ccpi/plugins/regularisers.py | 115 | ||||
| -rw-r--r-- | Wrappers/Python/wip/demo_compare_RGLTK_TV_denoising.py | 170 | ||||
| -rw-r--r-- | Wrappers/Python/wip/demo_simple_RGLTK.md | 214 |
3 files changed, 499 insertions, 0 deletions
diff --git a/Wrappers/Python/ccpi/plugins/regularisers.py b/Wrappers/Python/ccpi/plugins/regularisers.py new file mode 100644 index 0000000..e9c88a4 --- /dev/null +++ b/Wrappers/Python/ccpi/plugins/regularisers.py @@ -0,0 +1,115 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2018 Jakob Jorgensen, Daniil Kazantsev and Edoardo Pasca + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# This requires CCPi-Regularisation toolbox to be installed +from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV +from ccpi.filters.cpu_regularisers import TV_ENERGY +from ccpi.framework import DataContainer +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 + self.lambdaReg = lambdaReg + self.iterationsTV = iterationsTV + self.time_marchstep = time_marchstep + self.device = device # string for 'cpu' or 'gpu' + def __call__(self,x): + # 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) + return EnergyValTV + def prox(self,x,Lipshitz): + pars = {'algorithm' : ROF_TV, \ + 'input' : np.asarray(x.as_array(), dtype=np.float32),\ + 'regularization_parameter':self.lambdaReg*Lipshitz, \ + 'number_of_iterations' :self.iterationsTV ,\ + 'time_marching_parameter':self.time_marchstep} + + out = ROF_TV(pars['input'], + pars['regularization_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'], self.device) + return DataContainer(out) + +class _FGP_TV_(Operator): + def __init__(self,lambdaReg,iterationsTV,tolerance,methodTV,nonnegativity,printing,device): + # set parameters + self.lambdaReg = lambdaReg + self.iterationsTV = iterationsTV + self.tolerance = tolerance + self.methodTV = methodTV + self.nonnegativity = nonnegativity + self.printing = printing + self.device = device # string for 'cpu' or 'gpu' + def __call__(self,x): + # 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) + return EnergyValTV + def prox(self,x,Lipshitz): + pars = {'algorithm' : FGP_TV, \ + 'input' : np.asarray(x.as_array(), dtype=np.float32),\ + 'regularization_parameter':self.lambdaReg*Lipshitz, \ + 'number_of_iterations' :self.iterationsTV ,\ + 'tolerance_constant':self.tolerance,\ + 'methodTV': self.methodTV ,\ + 'nonneg': self.nonnegativity ,\ + 'printingOut': self.printing} + + out = FGP_TV(pars['input'], + pars['regularization_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['nonneg'], + pars['printingOut'], self.device) + return DataContainer(out) + + +class _SB_TV_(Operator): + def __init__(self,lambdaReg,iterationsTV,tolerance,methodTV,printing,device): + # set parameters + self.lambdaReg = lambdaReg + self.iterationsTV = iterationsTV + self.tolerance = tolerance + self.methodTV = methodTV + self.printing = printing + self.device = device # string for 'cpu' or 'gpu' + def __call__(self,x): + # 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) + return EnergyValTV + def prox(self,x,Lipshitz): + pars = {'algorithm' : SB_TV, \ + 'input' : np.asarray(x.as_array(), dtype=np.float32),\ + 'regularization_parameter':self.lambdaReg*Lipshitz, \ + 'number_of_iterations' :self.iterationsTV ,\ + 'tolerance_constant':self.tolerance,\ + 'methodTV': self.methodTV ,\ + 'printingOut': self.printing} + + out = SB_TV(pars['input'], + pars['regularization_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['printingOut'], self.device) + return DataContainer(out) 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..559679e --- /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('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 - TV') +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') +imgplot = plt.imshow(xtv_fgp.as_array()-xtv_denoise.value,vmin=dlims[0],vmax=dlims[1]) +plt.axis('off') diff --git a/Wrappers/Python/wip/demo_simple_RGLTK.md b/Wrappers/Python/wip/demo_simple_RGLTK.md new file mode 100644 index 0000000..9f0a4c3 --- /dev/null +++ b/Wrappers/Python/wip/demo_simple_RGLTK.md @@ -0,0 +1,214 @@ + +from ccpi.framework import ImageData , ImageGeometry, AcquisitionGeometry +from ccpi.optimisation.algs import FISTA, FBPD, CGLS +from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D +from ccpi.astra.ops import AstraProjectorSimple +from ccpi.plugins.regularisers import _ROF_TV_, _FGP_TV_ + +import numpy as np +import matplotlib.pyplot as plt + +test_case = 1 # 1=parallel2D, 2=cone2D + +# Set up phantom +N = 128 + + +vg = ImageGeometry(voxel_num_x=N,voxel_num_y=N) +Phantom = ImageData(geometry=vg) + +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.show() + +# Set up measurement geometry +angles_num = 20; # angles number + +if test_case==1: + angles = np.linspace(0,np.pi,angles_num,endpoint=False) +elif test_case==2: + angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) +else: + NotImplemented + +det_w = 1.0 +det_num = N +SourceOrig = 200 +OrigDetec = 0 + +# Parallelbeam geometry test +if test_case==1: + pg = AcquisitionGeometry('parallel', + '2D', + angles, + det_num,det_w) +elif test_case==2: + pg = AcquisitionGeometry('cone', + '2D', + angles, + det_num, + det_w, + dist_source_center=SourceOrig, + dist_center_detector=OrigDetec) + +# ASTRA operator using volume and sinogram geometries +Aop = AstraProjectorSimple(vg, pg, 'cpu') + +# Unused old astra projector without geometry +# Aop_old = AstraProjector(det_w, det_num, SourceOrig, +# OrigDetec, angles, +# N,'fanbeam','gpu') + +# Try forward and backprojection +b = Aop.direct(Phantom) +out2 = Aop.adjoint(b) + +#plt.imshow(b.array) +#plt.show() + +#plt.imshow(out2.array) +#plt.show() + +# Create least squares object instance with projector and data. +f = Norm2sq(Aop,b,c=0.5) + +# Initial guess +x_init = ImageData(np.zeros(x.shape),geometry=vg) +#%% +# FISTA with ROF-TV regularisation +g_rof = _ROF_TV_(lambdaReg = 10.0,iterationsTV=50,tolerance=1e-5,time_marchstep=0.01,device='cpu') + +opt = {'tol': 1e-4, 'iter': 100} + +x_fista_rof, it1, timing1, criter_rof = FISTA(x_init, f, g_rof,opt) + +plt.figure() +plt.subplot(121) +plt.imshow(x_fista_rof.array,cmap="BuPu") +plt.title('FISTA-ROF-TV') +plt.subplot(122) +plt.semilogy(criter_rof) +plt.show() +#%% +# FISTA with FGP-TV regularisation +g_fgp = _FGP_TV_(lambdaReg = 10.0,iterationsTV=50,tolerance=1e-5,methodTV=0,nonnegativity=0,printing=0,device='cpu') + +x_fista_fgp, it1, timing1, criter_fgp = FISTA(x_init, f, g_fgp,opt) + +plt.figure() +plt.subplot(121) +plt.imshow(x_fista_fgp.array,cmap="BuPu") +plt.title('FISTA-FGP-TV') +plt.subplot(122) +plt.semilogy(criter_fgp) +plt.show() +#%% +# Run FISTA for least squares without regularization +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt) + +plt.imshow(x_fista0.array) +plt.title('FISTA0') +plt.show() +#%% +# Now least squares plus 1-norm regularization +lam = 0.1 +g0 = Norm1(lam) + +# Run FISTA for least squares plus 1-norm function. +x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0) + +plt.imshow(x_fista1.array) +plt.title('FISTA1') +plt.show() + +plt.semilogy(criter1) +plt.show() +#%% +# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm +opt = {'tol': 1e-4, 'iter': 100} +x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt) + +plt.imshow(x_fbpd1.array) +plt.title('FBPD1') +plt.show() + +plt.semilogy(criter_fbpd1) +plt.show() +#%% +opt_FBPD = {'tol': 1e-4, 'iter': 10000} +# Now FBPD for least squares plus TV +lamtv = 10.0 +gtv = TV2D(lamtv) + +x_fbpdtv, it_fbpdtv, timing_fbpdtv, criter_fbpdtv = FBPD(x_init,None,f,gtv,opt=opt_FBPD) + +plt.imshow(x_fbpdtv.array) +plt.show() + +plt.semilogy(criter_fbpdtv) +plt.show() + + +# Run CGLS, which should agree with the FISTA0 +x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt ) + +plt.imshow(x_CGLS.array) +plt.title('CGLS') +#plt.title('CGLS recon, compare FISTA0') +plt.show() + +plt.semilogy(criter_CGLS) +plt.title('CGLS criterion') +plt.show() +#%% + +clims = (0,1) +cols = 3 +rows = 2 +current = 1 +fig = plt.figure() +# projections row +a=fig.add_subplot(rows,cols,current) +a.set_title('phantom {0}'.format(np.shape(Phantom.as_array()))) + +imgplot = plt.imshow(Phantom.as_array(),vmin=clims[0],vmax=clims[1]) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FISTA0') +imgplot = plt.imshow(x_fista0.as_array(),vmin=clims[0],vmax=clims[1]) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FISTA1') +imgplot = plt.imshow(x_fista1.as_array(),vmin=clims[0],vmax=clims[1]) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FBPD1') +imgplot = plt.imshow(x_fbpd1.as_array(),vmin=clims[0],vmax=clims[1]) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('CGLS') +imgplot = plt.imshow(x_CGLS.as_array(),vmin=clims[0],vmax=clims[1]) + +#current = current + 1 +#a=fig.add_subplot(rows,cols,current) +#a.set_title('FBPD TV') +#imgplot = plt.imshow(x_fbpdtv.as_array(),vmin=clims[0],vmax=clims[1]) + +fig = plt.figure() +# projections row +b=fig.add_subplot(1,1,1) +b.set_title('criteria') +imgplot = plt.loglog(criter0 , label='FISTA0') +imgplot = plt.loglog(criter1 , label='FISTA1') +imgplot = plt.loglog(criter_fbpd1, label='FBPD1') +imgplot = plt.loglog(criter_CGLS, label='CGLS') +b.legend(loc='right') +plt.show() +#%% |