From f12e27ba088462858f35315193cf6bf624325d1e Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Wed, 24 Apr 2019 10:49:06 +0100 Subject: fix memopt proximal --- .../ccpi/optimisation/functions/BlockFunction.py | 28 +++++++++++++++------- Wrappers/Python/wip/pdhg_TV_denoising.py | 2 +- 2 files changed, 21 insertions(+), 9 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py index bf627a5..0836324 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py +++ b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py @@ -93,16 +93,28 @@ class BlockFunction(Function): ''' + + if out is None: - out = [None]*self.length - if isinstance(tau, Number): - for i in range(self.length): - out[i] = self.functions[i].proximal(x.get_item(i), tau) + out = [None]*self.length + if isinstance(tau, Number): + for i in range(self.length): + out[i] = self.functions[i].proximal(x.get_item(i), tau) + else: + for i in range(self.length): + out[i] = self.functions[i].proximal(x.get_item(i), tau.get_item(i)) + + return BlockDataContainer(*out) + else: - for i in range(self.length): - out[i] = self.functions[i].proximal(x.get_item(i), tau.get_item(i)) - - return BlockDataContainer(*out) + if isinstance(tau, Number): + for i in range(self.length): + self.functions[i].proximal(x.get_item(i), tau, out[i]) + else: + for i in range(self.length): + self.functions[i].proximal(x.get_item(i), tau.get_item(i), out[i]) + + def gradient(self,x, out=None): diff --git a/Wrappers/Python/wip/pdhg_TV_denoising.py b/Wrappers/Python/wip/pdhg_TV_denoising.py index 872d832..a8e25d2 100755 --- a/Wrappers/Python/wip/pdhg_TV_denoising.py +++ b/Wrappers/Python/wip/pdhg_TV_denoising.py @@ -6,7 +6,7 @@ Created on Fri Feb 22 14:53:03 2019 @author: evangelos """ -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer +from ccpi.framework import ImageData, ImageGeometry import numpy as np import numpy -- cgit v1.2.3 From b0a2d779602df9f8abfcb68083f9de0df43ed5b0 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Wed, 24 Apr 2019 10:49:41 +0100 Subject: fix TGV example --- .../Python/ccpi/optimisation/algorithms/PDHG.py | 7 +- .../operators/SymmetrizedGradientOperator.py | 29 ++-- Wrappers/Python/wip/pdhg_TGV_denoising.py | 184 ++++++++++++--------- Wrappers/Python/wip/pdhg_TV_tomography2D.py | 6 +- 4 files changed, 130 insertions(+), 96 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py index 7631e29..110f998 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -104,6 +104,7 @@ class PDHG(Algorithm): self.y_old = self.y def update_objective(self): + p1 = self.f(self.operator.direct(self.x)) + self.g(self.x) d1 = -(self.f.convex_conjugate(self.y) + self.g.convex_conjugate(-1*self.operator.adjoint(self.y))) @@ -166,7 +167,7 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs): x_old.fill(x) y_old.fill(y) - if i%10==0: + if i%50==0: p1 = f(operator.direct(x)) + g(x) d1 = - ( f.convex_conjugate(y) + g.convex_conjugate(-1*operator.adjoint(y)) ) @@ -194,8 +195,8 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs): x_old.fill(x) y_old.fill(y) - - if i%10==0: + + if i%50==0: p1 = f(operator.direct(x)) + g(x) d1 = - ( f.convex_conjugate(y) + g.convex_conjugate(-1*operator.adjoint(y)) ) diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py index 3c15fa1..243809b 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py @@ -69,9 +69,8 @@ class SymmetrizedGradient(Gradient): self.FD.direction = i self.FD.adjoint(x.get_item(j), out=out[ind]) ind+=1 - out1 = [out[i] for i in self.order_ind] - res = [0.5 * sum(x) for x in zip(out, out1)] - out.fill(BlockDataContainer(*res)) + out1 = BlockDataContainer(*[out[i] for i in self.order_ind]) + out.fill( 0.5 * (out + out1) ) def adjoint(self, x, out=None): @@ -102,10 +101,11 @@ class SymmetrizedGradient(Gradient): self.FD.direct(x[i], out=tmp[j]) i+=1 tmp1+=tmp[j] - out[k].fill(tmp1) + out[k].fill(tmp1) +# tmp = self.adjoint(x) +# out.fill(tmp) - - + def domain_geometry(self): return self.gm_domain @@ -114,12 +114,12 @@ class SymmetrizedGradient(Gradient): def norm(self): - #TODO need dot method for BlockDataContainer - return numpy.sqrt(4*self.gm_domain.shape[0]) +# TODO need dot method for BlockDataContainer +# return numpy.sqrt(4*self.gm_domain.shape[0]) -# x0 = self.gm_domain.allocate('random_int') -# self.s1, sall, svec = LinearOperator.PowerMethod(self, 10, x0) -# return self.s1 +# x0 = self.gm_domain.allocate('random') + self.s1, sall, svec = LinearOperator.PowerMethod(self, 50) + return self.s1 @@ -226,6 +226,13 @@ if __name__ == '__main__': print(LHS_out, RHS_out) + out = E1.range_geometry().allocate() + E1.direct(u1, out=out) + E1.adjoint(out, out=out1) + + print(E1.norm()) + print(E2.norm()) + diff --git a/Wrappers/Python/wip/pdhg_TGV_denoising.py b/Wrappers/Python/wip/pdhg_TGV_denoising.py index 0b6e594..1c570cb 100644 --- a/Wrappers/Python/wip/pdhg_TGV_denoising.py +++ b/Wrappers/Python/wip/pdhg_TGV_denoising.py @@ -105,100 +105,126 @@ normK = operator.norm() sigma = 1 tau = 1/(sigma*normK**2) ## -opt = {'niter':2000} -opt1 = {'niter':2000, 'memopt': True} +opt = {'niter':500} +opt1 = {'niter':500, 'memopt': True} # t1 = timer() res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) t2 = timer() # +t3 = timer() +res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) +t4 = timer() +# +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) plt.imshow(res[0].as_array()) +plt.title('no memopt') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(res1[0].as_array()) +plt.title('memopt') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow((res1[0] - res[0]).abs().as_array()) +plt.title('diff') +plt.colorbar() plt.show() +print("NoMemopt/Memopt is {}/{}".format(t2-t1, t4-t3)) + + +###### + +#%% -#t3 = timer() -#res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -#t4 = timer() +#def update_plot(it_update, objective, x): +# +## sol = pdhg.get_output() +# plt.figure() +# plt.imshow(x[0].as_array()) +# plt.show() +# +# +##def stop_criterion(x,) # -#plt.figure(figsize=(15,15)) -#plt.subplot(3,1,1) -#plt.imshow(res[0].as_array()) -#plt.title('no memopt') -#plt.colorbar() -#plt.subplot(3,1,2) -#plt.imshow(res1[0].as_array()) -#plt.title('memopt') -#plt.colorbar() -#plt.subplot(3,1,3) -#plt.imshow((res1[0] - res[0]).abs().as_array()) -#plt.title('diff') -#plt.colorbar() -#plt.show() +#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +#pdhg.max_iteration = 2000 +#pdhg.update_objective_interval = 100 # -#print("NoMemopt/Memopt is {}/{}".format(t2-t1, t4-t3)) - +#pdhg.run(4000, verbose=False, callback=update_plot) -#%% Check with CVX solution -from ccpi.optimisation.operators import SparseFiniteDiff +#%% -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False -if cvx_not_installable: - - u = Variable(ig.shape) - w1 = Variable((N, N)) - w2 = Variable((N, N)) - - # create TGV regulariser - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u) - vec(w1), \ - DY.matrix() * vec(u) - vec(w2)]), 2, axis = 0)) + \ - beta * sum(norm(vstack([ DX.matrix().transpose() * vec(w1), DY.matrix().transpose() * vec(w2), \ - 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ), \ - 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ) ]), 2, axis = 0 ) ) - - constraints = [] - fidelity = 0.5 * sum_squares(u - noisy_data.as_array()) - solver = MOSEK - obj = Minimize( regulariser + fidelity) - prob = Problem(obj) - result = prob.solve(verbose = True, solver = solver) - - diff_cvx = numpy.abs( res[0].as_array() - u.value ) - - # Show result - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(res[0].as_array()) - plt.title('PDHG solution') - plt.colorbar() - - plt.subplot(3,1,2) - plt.imshow(u.value) - plt.title('CVX solution') - plt.colorbar() - - plt.subplot(3,1,3) - plt.imshow(diff_cvx) - plt.title('Difference') - plt.colorbar() - plt.show() - - plt.plot(np.linspace(0,N,N), res[0].as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') - plt.legend() - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(primal[-1])) - print('Min/Max of absolute difference {}/{}'.format(diff_cvx.min(), diff_cvx.max())) - + + + + + +#%% Check with CVX solution + +#from ccpi.optimisation.operators import SparseFiniteDiff +# +#try: +# from cvxpy import * +# cvx_not_installable = True +#except ImportError: +# cvx_not_installable = False +# +#if cvx_not_installable: +# +# u = Variable(ig.shape) +# w1 = Variable((N, N)) +# w2 = Variable((N, N)) +# +# # create TGV regulariser +# DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') +# DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') +# +# regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u) - vec(w1), \ +# DY.matrix() * vec(u) - vec(w2)]), 2, axis = 0)) + \ +# beta * sum(norm(vstack([ DX.matrix().transpose() * vec(w1), DY.matrix().transpose() * vec(w2), \ +# 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ), \ +# 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ) ]), 2, axis = 0 ) ) +# +# constraints = [] +# fidelity = 0.5 * sum_squares(u - noisy_data.as_array()) +# solver = MOSEK +# +# obj = Minimize( regulariser + fidelity) +# prob = Problem(obj) +# result = prob.solve(verbose = True, solver = solver) +# +# diff_cvx = numpy.abs( res[0].as_array() - u.value ) +# +# # Show result +# plt.figure(figsize=(15,15)) +# plt.subplot(3,1,1) +# plt.imshow(res[0].as_array()) +# plt.title('PDHG solution') +# plt.colorbar() +# +# plt.subplot(3,1,2) +# plt.imshow(u.value) +# plt.title('CVX solution') +# plt.colorbar() +# +# plt.subplot(3,1,3) +# plt.imshow(diff_cvx) +# plt.title('Difference') +# plt.colorbar() +# plt.show() +# +# plt.plot(np.linspace(0,N,N), res[0].as_array()[int(N/2),:], label = 'PDHG') +# plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') +# plt.legend() +# +# print('Primal Objective (CVX) {} '.format(obj.value)) +# print('Primal Objective (PDHG) {} '.format(primal[-1])) +# print('Min/Max of absolute difference {}/{}'.format(diff_cvx.min(), diff_cvx.max())) +# diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D.py b/Wrappers/Python/wip/pdhg_TV_tomography2D.py index 2713cfd..91c48c7 100644 --- a/Wrappers/Python/wip/pdhg_TV_tomography2D.py +++ b/Wrappers/Python/wip/pdhg_TV_tomography2D.py @@ -8,16 +8,16 @@ Created on Fri Feb 22 14:53:03 2019 @author: evangelos """ -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, AcquisitionGeometry, AcquisitionData +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData import numpy as np import matplotlib.pyplot as plt from ccpi.optimisation.algorithms import PDHG, PDHG_old -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient +from ccpi.optimisation.operators import BlockOperator, Gradient from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction, ScaledFunction + MixedL21Norm, BlockFunction from ccpi.astra.ops import AstraProjectorSimple from skimage.util import random_noise -- cgit v1.2.3 From 0604ebff2fccfe1fa48defcecfba6dd75ef3249b Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Wed, 24 Apr 2019 10:49:59 +0100 Subject: fix TGV tomo2D example --- Wrappers/Python/wip/pdhg_TGV_tomography2D.py | 200 +++++++++++++++++++++++++++ 1 file changed, 200 insertions(+) create mode 100644 Wrappers/Python/wip/pdhg_TGV_tomography2D.py (limited to 'Wrappers') diff --git a/Wrappers/Python/wip/pdhg_TGV_tomography2D.py b/Wrappers/Python/wip/pdhg_TGV_tomography2D.py new file mode 100644 index 0000000..ee3b089 --- /dev/null +++ b/Wrappers/Python/wip/pdhg_TGV_tomography2D.py @@ -0,0 +1,200 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Feb 22 14:53:03 2019 + +@author: evangelos +""" + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG, PDHG_old + +from ccpi.optimisation.operators import BlockOperator, Identity, \ + Gradient, SymmetrizedGradient, ZeroOperator +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ + MixedL21Norm, BlockFunction + +from skimage.util import random_noise + +from timeit import default_timer as timer +from ccpi.astra.ops import AstraProjectorSimple + +#def dt(steps): +# return steps[-1] - steps[-2] + +# Create phantom for TGV Gaussian denoising + +N = 100 + +data = np.zeros((N,N)) + +x1 = np.linspace(0, int(N/2), N) +x2 = np.linspace(int(N/2), 0., N) +xv, yv = np.meshgrid(x1, x2) + +xv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1] = yv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1].T + +data = xv +data = ImageData(data/data.max()) + +plt.imshow(data.as_array()) +plt.show() + +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) + +detectors = N +angles = np.linspace(0,np.pi,N) + +ag = AcquisitionGeometry('parallel','2D',angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +plt.imshow(sin.as_array()) +plt.title('Sinogram') +plt.colorbar() +plt.show() + +# Add Gaussian noise to the sinogram data +np.random.seed(10) +n1 = np.random.random(sin.shape) + +noisy_data = sin + ImageData(5*n1) + +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Sinogram') +plt.colorbar() +plt.show() + +#%% + +alpha, beta = 20, 50 + + +# Create operators +op11 = Gradient(ig) +op12 = Identity(op11.range_geometry()) + +op22 = SymmetrizedGradient(op11.domain_geometry()) +op21 = ZeroOperator(ig, op22.range_geometry()) + +op31 = Aop +op32 = ZeroOperator(op22.domain_geometry(), ag) + +operator = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) + + +f1 = alpha * MixedL21Norm() +f2 = beta * MixedL21Norm() +f3 = 0.5 * L2NormSquared(b = noisy_data) +f = BlockFunction(f1, f2, f3) +g = ZeroFunction() + +## Compute operator Norm +normK = operator.norm() +# +## Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) +## +opt = {'niter':5000} +opt1 = {'niter':5000, 'memopt': True} +# +t1 = timer() +res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) +t2 = timer() +# +plt.imshow(res[0].as_array()) +plt.show() + + +#t3 = timer() +#res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) +#t4 = timer() +# +#plt.figure(figsize=(15,15)) +#plt.subplot(3,1,1) +#plt.imshow(res[0].as_array()) +#plt.title('no memopt') +#plt.colorbar() +#plt.subplot(3,1,2) +#plt.imshow(res1[0].as_array()) +#plt.title('memopt') +#plt.colorbar() +#plt.subplot(3,1,3) +#plt.imshow((res1[0] - res[0]).abs().as_array()) +#plt.title('diff') +#plt.colorbar() +#plt.show() +# +#print("NoMemopt/Memopt is {}/{}".format(t2-t1, t4-t3)) + + +#%% Check with CVX solution + +#from ccpi.optimisation.operators import SparseFiniteDiff +# +#try: +# from cvxpy import * +# cvx_not_installable = True +#except ImportError: +# cvx_not_installable = False +# +#if cvx_not_installable: +# +# u = Variable(ig.shape) +# w1 = Variable((N, N)) +# w2 = Variable((N, N)) +# +# # create TGV regulariser +# DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') +# DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') +# +# regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u) - vec(w1), \ +# DY.matrix() * vec(u) - vec(w2)]), 2, axis = 0)) + \ +# beta * sum(norm(vstack([ DX.matrix().transpose() * vec(w1), DY.matrix().transpose() * vec(w2), \ +# 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ), \ +# 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ) ]), 2, axis = 0 ) ) +# +# constraints = [] +# fidelity = 0.5 * sum_squares(u - noisy_data.as_array()) +# solver = MOSEK +# +# obj = Minimize( regulariser + fidelity) +# prob = Problem(obj) +# result = prob.solve(verbose = True, solver = solver) +# +# diff_cvx = numpy.abs( res[0].as_array() - u.value ) +# +# # Show result +# plt.figure(figsize=(15,15)) +# plt.subplot(3,1,1) +# plt.imshow(res[0].as_array()) +# plt.title('PDHG solution') +# plt.colorbar() +# +# plt.subplot(3,1,2) +# plt.imshow(u.value) +# plt.title('CVX solution') +# plt.colorbar() +# +# plt.subplot(3,1,3) +# plt.imshow(diff_cvx) +# plt.title('Difference') +# plt.colorbar() +# plt.show() +# +# plt.plot(np.linspace(0,N,N), res[0].as_array()[int(N/2),:], label = 'PDHG') +# plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') +# plt.legend() +# +# print('Primal Objective (CVX) {} '.format(obj.value)) +# print('Primal Objective (PDHG) {} '.format(primal[-1])) +# print('Min/Max of absolute difference {}/{}'.format(diff_cvx.min(), diff_cvx.max())) + + + -- cgit v1.2.3 From 2e8b2def18c2d57920ad063056540873cb30d292 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Wed, 24 Apr 2019 13:49:21 +0100 Subject: test pdhg tomo TV 2D --- .../Python/ccpi/optimisation/algorithms/PDHG.py | 78 ++++++++++++++-------- Wrappers/Python/wip/pdhg_TV_tomography2D.py | 2 +- 2 files changed, 51 insertions(+), 29 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py index 110f998..0f5e8ef 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -47,34 +47,43 @@ class PDHG(Algorithm): self.y_old = self.operator.range_geometry().allocate() self.xbar = self.x_old.copy() - + self.x_tmp = self.x_old.copy() self.x = self.x_old.copy() - self.y = self.y_old.copy() - if self.memopt: - self.y_tmp = self.y_old.copy() - self.x_tmp = self.x_old.copy() - #y = y_tmp + + self.y_tmp = self.y_old.copy() + self.y = self.y_tmp.copy() + + + + #self.y = self.y_old.copy() + + + #if self.memopt: + # self.y_tmp = self.y_old.copy() + # self.x_tmp = self.x_old.copy() + # relaxation parameter self.theta = 1 def update(self): + if self.memopt: # Gradient descent, Dual problem solution # self.y_old += self.sigma * self.operator.direct(self.xbar) self.operator.direct(self.xbar, out=self.y_tmp) self.y_tmp *= self.sigma - self.y_old += self.y_tmp + self.y_tmp += self.y_old #self.y = self.f.proximal_conjugate(self.y_old, self.sigma) - self.f.proximal_conjugate(self.y_old, self.sigma, out=self.y) + self.f.proximal_conjugate(self.y_tmp, self.sigma, out=self.y) # Gradient ascent, Primal problem solution self.operator.adjoint(self.y, out=self.x_tmp) - self.x_tmp *= self.tau - self.x_old -= self.x_tmp + self.x_tmp *= -1*self.tau + self.x_tmp += self.x_old - self.g.proximal(self.x_old, self.tau, out=self.x) + self.g.proximal(self.x_tmp, self.tau, out=self.x) #Update self.x.subtract(self.x_old, out=self.xbar) @@ -83,7 +92,8 @@ class PDHG(Algorithm): self.x_old.fill(self.x) self.y_old.fill(self.y) - + + else: # Gradient descent, Dual problem solution self.y_old += self.sigma * self.operator.direct(self.xbar) @@ -94,14 +104,23 @@ class PDHG(Algorithm): self.x = self.g.proximal(self.x_old, self.tau) #Update - #xbar = x + theta * (x - x_old) - self.xbar.fill(self.x) - self.xbar -= self.x_old + self.x.subtract(self.x_old, out=self.xbar) self.xbar *= self.theta self.xbar += self.x - self.x_old = self.x - self.y_old = self.y + self.x_old.fill(self.x) + self.y_old.fill(self.y) + + #xbar = x + theta * (x - x_old) +# self.xbar.fill(self.x) +# self.xbar -= self.x_old +# self.xbar *= self.theta +# self.xbar += self.x + +# self.x_old.fill(self.x) +# self.y_old.fill(self.y) + + def update_objective(self): @@ -153,7 +172,7 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs): if not memopt: - + y_tmp = y_old + sigma * operator.direct(xbar) y = f.proximal_conjugate(y_tmp, sigma) @@ -163,10 +182,7 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs): x.subtract(x_old, out=xbar) xbar *= theta xbar += x - - x_old.fill(x) - y_old.fill(y) - + if i%50==0: p1 = f(operator.direct(x)) + g(x) @@ -174,7 +190,11 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs): primal.append(p1) dual.append(d1) pdgap.append(p1-d1) - print(p1, d1, p1-d1) + print(p1, d1, p1-d1) + + x_old.fill(x) + y_old.fill(y) + else: @@ -192,10 +212,7 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs): x.subtract(x_old, out=xbar) xbar *= theta xbar += x - - x_old.fill(x) - y_old.fill(y) - + if i%50==0: p1 = f(operator.direct(x)) + g(x) @@ -203,7 +220,12 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs): primal.append(p1) dual.append(d1) pdgap.append(p1-d1) - print(p1, d1, p1-d1) + print(p1, d1, p1-d1) + + x_old.fill(x) + y_old.fill(y) + + t_end = time.time() diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D.py b/Wrappers/Python/wip/pdhg_TV_tomography2D.py index 91c48c7..0e167e3 100644 --- a/Wrappers/Python/wip/pdhg_TV_tomography2D.py +++ b/Wrappers/Python/wip/pdhg_TV_tomography2D.py @@ -126,7 +126,7 @@ t3 = timer() res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) t4 = timer() # - +# plt.figure(figsize=(15,15)) plt.subplot(3,1,1) plt.imshow(res.as_array()) -- cgit v1.2.3 From 90431756d8c385e51ae4c4de0c7a280e466cf915 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Wed, 24 Apr 2019 15:44:52 +0100 Subject: add complete Demos PDHG --- .../wip/Demos/PDHG_TGV_Denoising_SaltPepper.py | 196 +++++++++++++++++++++ .../Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py | 177 +++++++++++++++++++ .../Python/wip/Demos/PDHG_TV_Denoising_Poisson.py | 177 +++++++++++++++++++ .../wip/Demos/PDHG_TV_Denoising_SaltPepper.py | 177 +++++++++++++++++++ .../wip/Demos/sinogram_demo_tomography2D.npy | Bin 0 -> 60128 bytes Wrappers/Python/wip/Demos/untitled4.py | 87 +++++++++ 6 files changed, 814 insertions(+) create mode 100644 Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py create mode 100644 Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py create mode 100644 Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py create mode 100644 Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py create mode 100644 Wrappers/Python/wip/Demos/sinogram_demo_tomography2D.npy create mode 100644 Wrappers/Python/wip/Demos/untitled4.py (limited to 'Wrappers') diff --git a/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py b/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py new file mode 100644 index 0000000..dffb6d8 --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py @@ -0,0 +1,196 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Feb 22 14:53:03 2019 + +@author: evangelos +""" + +from ccpi.framework import ImageData, ImageGeometry + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG + +from ccpi.optimisation.operators import BlockOperator, Identity, \ + Gradient, SymmetrizedGradient, ZeroOperator +from ccpi.optimisation.functions import ZeroFunction, L1Norm, \ + MixedL21Norm, BlockFunction + +from skimage.util import random_noise + +# Create phantom for TGV SaltPepper denoising + +N = 300 + +data = np.zeros((N,N)) + +x1 = np.linspace(0, int(N/2), N) +x2 = np.linspace(int(N/2), 0., N) +xv, yv = np.meshgrid(x1, x2) + +xv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1] = yv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1].T + +data = xv +data = ImageData(data/data.max()) + +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +ag = ig + +# Create noisy data. Add Gaussian noise +n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2) +noisy_data = ImageData(n1) + +# Regularisation Parameters +alpha = 0.8 +beta = numpy.sqrt(2)* alpha + +method = '0' + +if method == '0': + + # Create operators + op11 = Gradient(ig) + op12 = Identity(op11.range_geometry()) + + op22 = SymmetrizedGradient(op11.domain_geometry()) + op21 = ZeroOperator(ig, op22.range_geometry()) + + op31 = Identity(ig, ag) + op32 = ZeroOperator(op22.domain_geometry(), ag) + + operator = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) + + f1 = alpha * MixedL21Norm() + f2 = beta * MixedL21Norm() + f3 = L1Norm(b=noisy_data) + f = BlockFunction(f1, f2, f3) + g = ZeroFunction() + +else: + + # Create operators + op11 = Gradient(ig) + op12 = Identity(op11.range_geometry()) + op22 = SymmetrizedGradient(op11.domain_geometry()) + op21 = ZeroOperator(ig, op22.range_geometry()) + + operator = BlockOperator(op11, -1*op12, op21, op22, shape=(2,2) ) + + f1 = alpha * MixedL21Norm() + f2 = beta * MixedL21Norm() + + f = BlockFunction(f1, f2) + g = BlockFunction(0.5 * L1Norm(b=noisy_data), ZeroFunction()) + +## Compute operator Norm +normK = operator.norm() +# +# Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) + + +# Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 50 +pdhg.run(2000) + +#%% +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(pdhg.get_output()[0].as_array()) +plt.title('TGV Reconstruction') +plt.colorbar() +plt.show() +## +plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') +plt.plot(np.linspace(0,N,N), pdhg.get_output()[0].as_array()[int(N/2),:], label = 'TV reconstruction') +plt.legend() +plt.title('Middle Line Profiles') +plt.show() + + + + + + + + + + +#%% Check with CVX solution + +#from ccpi.optimisation.operators import SparseFiniteDiff +# +#try: +# from cvxpy import * +# cvx_not_installable = True +#except ImportError: +# cvx_not_installable = False +# +#if cvx_not_installable: +# +# u = Variable(ig.shape) +# w1 = Variable((N, N)) +# w2 = Variable((N, N)) +# +# # create TGV regulariser +# DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') +# DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') +# +# regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u) - vec(w1), \ +# DY.matrix() * vec(u) - vec(w2)]), 2, axis = 0)) + \ +# beta * sum(norm(vstack([ DX.matrix().transpose() * vec(w1), DY.matrix().transpose() * vec(w2), \ +# 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ), \ +# 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ) ]), 2, axis = 0 ) ) +# +# constraints = [] +# fidelity = 0.5 * sum_squares(u - noisy_data.as_array()) +# solver = MOSEK +# +# obj = Minimize( regulariser + fidelity) +# prob = Problem(obj) +# result = prob.solve(verbose = True, solver = solver) +# +# diff_cvx = numpy.abs( res[0].as_array() - u.value ) +# +# # Show result +# plt.figure(figsize=(15,15)) +# plt.subplot(3,1,1) +# plt.imshow(res[0].as_array()) +# plt.title('PDHG solution') +# plt.colorbar() +# +# plt.subplot(3,1,2) +# plt.imshow(u.value) +# plt.title('CVX solution') +# plt.colorbar() +# +# plt.subplot(3,1,3) +# plt.imshow(diff_cvx) +# plt.title('Difference') +# plt.colorbar() +# plt.show() +# +# plt.plot(np.linspace(0,N,N), res[0].as_array()[int(N/2),:], label = 'PDHG') +# plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') +# plt.legend() +# +# print('Primal Objective (CVX) {} '.format(obj.value)) +# print('Primal Objective (PDHG) {} '.format(primal[-1])) +# print('Min/Max of absolute difference {}/{}'.format(diff_cvx.min(), diff_cvx.max())) +# + + diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py new file mode 100644 index 0000000..10e5621 --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py @@ -0,0 +1,177 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.framework import ImageData, ImageGeometry + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG + +from ccpi.optimisation.operators import BlockOperator, Identity, Gradient +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ + MixedL21Norm, BlockFunction + +from skimage.util import random_noise + +# Create phantom for TV Gaussian denoising +N = 300 + +data = np.zeros((N,N)) +data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 +data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 +data = ImageData(data) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +ag = ig + +# Create noisy data. Add Gaussian noise +n1 = random_noise(data.as_array(), mode = 'gaussian', mean=0, var = 0.05, seed=10) +noisy_data = ImageData(n1) + +# Regularisation Parameter +alpha = 0.5 + +method = '1' + +if method == '0': + + # Create operators + op1 = Gradient(ig) + op2 = Identity(ig, ag) + + # Create BlockOperator + operator = BlockOperator(op1, op2, shape=(2,1) ) + + # Create functions + + f1 = alpha * MixedL21Norm() + f2 = 0.5 * L2NormSquared(b = noisy_data) + f = BlockFunction(f1, f2) + + g = ZeroFunction() + +else: + + # Without the "Block Framework" + operator = Gradient(ig) + f = alpha * MixedL21Norm() + g = 0.5 * L2NormSquared(b = noisy_data) + + +# Compute operator Norm +normK = operator.norm() + +# Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) + + +# Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 50 +pdhg.run(2000) + + +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(pdhg.get_output().as_array()) +plt.title('TV Reconstruction') +plt.colorbar() +plt.show() +## +plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') +plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction') +plt.legend() +plt.title('Middle Line Profiles') +plt.show() + + +##%% Check with CVX solution + +from ccpi.optimisation.operators import SparseFiniteDiff + +try: + from cvxpy import * + cvx_not_installable = True +except ImportError: + cvx_not_installable = False + + +if cvx_not_installable: + + ##Construct problem + u = Variable(ig.shape) + + DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') + DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') + + # Define Total Variation as a regulariser + regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) + fidelity = 0.5 * sum_squares(u - noisy_data.as_array()) + + # choose solver + if 'MOSEK' in installed_solvers(): + solver = MOSEK + else: + solver = SCS + + obj = Minimize( regulariser + fidelity) + prob = Problem(obj) + result = prob.solve(verbose = True, solver = solver) + + diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value ) + + plt.figure(figsize=(15,15)) + plt.subplot(3,1,1) + plt.imshow(pdhg.get_output().as_array()) + plt.title('PDHG solution') + plt.colorbar() + plt.subplot(3,1,2) + plt.imshow(u.value) + plt.title('CVX solution') + plt.colorbar() + plt.subplot(3,1,3) + plt.imshow(diff_cvx) + plt.title('Difference') + plt.colorbar() + plt.show() + + plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') + plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') + plt.legend() + plt.title('Middle Line Profiles') + plt.show() + + print('Primal Objective (CVX) {} '.format(obj.value)) + print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) + + + + + diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py new file mode 100644 index 0000000..f2c2023 --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py @@ -0,0 +1,177 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.framework import ImageData, ImageGeometry + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG + +from ccpi.optimisation.operators import BlockOperator, Identity, Gradient +from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \ + MixedL21Norm, BlockFunction + +from skimage.util import random_noise + +# Create phantom for TV Poisson denoising +N = 300 + +data = np.zeros((N,N)) +data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 +data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 +data = ImageData(data) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +ag = ig + +# Create noisy data. Apply Poisson noise +n1 = random_noise(data.as_array(), mode = 'poisson', seed = 10) +noisy_data = ImageData(n1) + +# Regularisation Parameter +alpha = 2 + +method = '0' + +if method == '0': + + # Create operators + op1 = Gradient(ig) + op2 = Identity(ig, ag) + + # Create BlockOperator + operator = BlockOperator(op1, op2, shape=(2,1) ) + + # Create functions + + f1 = alpha * MixedL21Norm() + f2 = KullbackLeibler(noisy_data) + f = BlockFunction(f1, f2) + + g = ZeroFunction() + +else: + + # Without the "Block Framework" + operator = Gradient(ig) + f = alpha * MixedL21Norm() + g = KullbackLeibler(noisy_data) + + +# Compute operator Norm +normK = operator.norm() + +# Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) +opt = {'niter':2000, 'memopt': True} + +# Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 50 +pdhg.run(2000) + + +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(pdhg.get_output().as_array()) +plt.title('TV Reconstruction') +plt.colorbar() +plt.show() +## +plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') +plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction') +plt.legend() +plt.title('Middle Line Profiles') +plt.show() + + +##%% Check with CVX solution +#%% Check with CVX solution + +from ccpi.optimisation.operators import SparseFiniteDiff + +try: + from cvxpy import * + cvx_not_installable = True +except ImportError: + cvx_not_installable = False + + +if cvx_not_installable: + + ##Construct problem + u1 = Variable(ig.shape) + q = Variable() + + DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') + DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') + + # Define Total Variation as a regulariser + regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u1), DY.matrix() * vec(u1)]), 2, axis = 0)) + + fidelity = sum( u1 - multiply(noisy_data.as_array(), log(u1)) ) + constraints = [q>= fidelity, u1>=0] + + solver = ECOS + obj = Minimize( regulariser + q) + prob = Problem(obj, constraints) + result = prob.solve(verbose = True, solver = solver) + + + diff_cvx = numpy.abs( pdhg.get_output().as_array() - u1.value ) + + plt.figure(figsize=(15,15)) + plt.subplot(3,1,1) + plt.imshow(pdhg.get_output().as_array()) + plt.title('PDHG solution') + plt.colorbar() + plt.subplot(3,1,2) + plt.imshow(u1.value) + plt.title('CVX solution') + plt.colorbar() + plt.subplot(3,1,3) + plt.imshow(diff_cvx) + plt.title('Difference') + plt.colorbar() + plt.show() + + plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') + plt.plot(np.linspace(0,N,N), u1.value[int(N/2),:], label = 'CVX') + plt.legend() + plt.title('Middle Line Profiles') + plt.show() + + print('Primal Objective (CVX) {} '.format(obj.value)) + print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) + + + + + diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py new file mode 100644 index 0000000..f96acd5 --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py @@ -0,0 +1,177 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.framework import ImageData, ImageGeometry + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG + +from ccpi.optimisation.operators import BlockOperator, Identity, Gradient +from ccpi.optimisation.functions import ZeroFunction, L1Norm, \ + MixedL21Norm, BlockFunction + +from skimage.util import random_noise + +# Create phantom for TV Salt & Pepper denoising +N = 300 + +data = np.zeros((N,N)) +data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 +data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 +data = ImageData(data) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +ag = ig + +# Create noisy data. Apply Salt & Pepper noise +n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2) +noisy_data = ImageData(n1) + +# Regularisation Parameter +alpha = 2 + +method = '0' + +if method == '0': + + # Create operators + op1 = Gradient(ig) + op2 = Identity(ig, ag) + + # Create BlockOperator + operator = BlockOperator(op1, op2, shape=(2,1) ) + + # Create functions + + f1 = alpha * MixedL21Norm() + f2 = L1Norm(b = noisy_data) + f = BlockFunction(f1, f2) + + g = ZeroFunction() + +else: + + # Without the "Block Framework" + operator = Gradient(ig) + f = alpha * MixedL21Norm() + g = L1Norm(b = noisy_data) + + +# Compute operator Norm +normK = operator.norm() + +# Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) +opt = {'niter':2000, 'memopt': True} + +# Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 50 +pdhg.run(2000) + + +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(pdhg.get_output().as_array()) +plt.title('TV Reconstruction') +plt.colorbar() +plt.show() +## +plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') +plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction') +plt.legend() +plt.title('Middle Line Profiles') +plt.show() + + +##%% Check with CVX solution + +from ccpi.optimisation.operators import SparseFiniteDiff + +try: + from cvxpy import * + cvx_not_installable = True +except ImportError: + cvx_not_installable = False + + +if cvx_not_installable: + + ##Construct problem + u = Variable(ig.shape) + + DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') + DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') + + # Define Total Variation as a regulariser + regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) + fidelity = pnorm( u - noisy_data.as_array(),1) + + # choose solver + if 'MOSEK' in installed_solvers(): + solver = MOSEK + else: + solver = SCS + + obj = Minimize( regulariser + fidelity) + prob = Problem(obj) + result = prob.solve(verbose = True, solver = solver) + + diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value ) + + plt.figure(figsize=(15,15)) + plt.subplot(3,1,1) + plt.imshow(pdhg.get_output().as_array()) + plt.title('PDHG solution') + plt.colorbar() + plt.subplot(3,1,2) + plt.imshow(u.value) + plt.title('CVX solution') + plt.colorbar() + plt.subplot(3,1,3) + plt.imshow(diff_cvx) + plt.title('Difference') + plt.colorbar() + plt.show() + + plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') + plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') + plt.legend() + plt.title('Middle Line Profiles') + plt.show() + + print('Primal Objective (CVX) {} '.format(obj.value)) + print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) + + + + + diff --git a/Wrappers/Python/wip/Demos/sinogram_demo_tomography2D.npy b/Wrappers/Python/wip/Demos/sinogram_demo_tomography2D.npy new file mode 100644 index 0000000..f37fd4b Binary files /dev/null and b/Wrappers/Python/wip/Demos/sinogram_demo_tomography2D.npy differ diff --git a/Wrappers/Python/wip/Demos/untitled4.py b/Wrappers/Python/wip/Demos/untitled4.py new file mode 100644 index 0000000..0cacbd7 --- /dev/null +++ b/Wrappers/Python/wip/Demos/untitled4.py @@ -0,0 +1,87 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Apr 24 14:21:08 2019 + +@author: evangelos +""" + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData +import numpy +import numpy as np +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG, PDHG_old + +from ccpi.optimisation.operators import BlockOperator, Gradient +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ + MixedL21Norm, BlockFunction + +from ccpi.astra.ops import AstraProjectorSimple, AstraProjector3DSimple +from skimage.util import random_noise +from timeit import default_timer as timer + +#N = 75 +#x = np.zeros((N,N)) +#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 + +#data = ImageData(x) + +N = 75 +#x = np.zeros((N,N)) + +vert = 4 +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, voxel_num_z=vert) +data = ig.allocate() +Phantom = data +# Populate image data by looping over and filling slices +i = 0 +while i < vert: + if vert > 1: + x = Phantom.subset(vertical=i).array + else: + x = Phantom.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)] = 0.98 + if vert > 1 : + Phantom.fill(x, vertical=i) + i += 1 + +angles_num = 100 +det_w = 1.0 +det_num = N + +angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*\ + 180/np.pi + +# Inputs: Geometry, 2D or 3D, angles, horz detector pixel count, +# horz detector pixel size, vert detector pixel count, +# vert detector pixel size. +ag = AcquisitionGeometry('parallel', + '3D', + angles, + N, + det_w, + vert, + det_w) + +sino = numpy.load("sinogram_demo_tomography2D.npy", mmap_mode='r') +plt.imshow(sino) +plt.title('Sinogram CCPi') +plt.colorbar() +plt.show() + +#%% +Aop = AstraProjector3DSimple(ig, ag) +sin = Aop.direct(data) + +plt.imshow(sin.as_array()) + +plt.title('Sinogram Astra') +plt.colorbar() +plt.show() + + + +#%% \ No newline at end of file -- cgit v1.2.3 From e05bed848d826d219efba6f78354b4c3f76161cc Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 25 Apr 2019 10:06:59 +0100 Subject: add denoising demos --- .../wip/Demos/PDHG_TGV_Denoising_SaltPepper.py | 132 ++++++++++----------- .../Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py | 2 +- .../Python/wip/Demos/PDHG_TV_Denoising_Poisson.py | 4 +- .../wip/Demos/PDHG_TV_Denoising_SaltPepper.py | 4 +- 4 files changed, 70 insertions(+), 72 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py b/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py index dffb6d8..7b65c31 100644 --- a/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py +++ b/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py @@ -23,7 +23,7 @@ from skimage.util import random_noise # Create phantom for TGV SaltPepper denoising -N = 300 +N = 100 data = np.zeros((N,N)) @@ -47,7 +47,7 @@ noisy_data = ImageData(n1) alpha = 0.8 beta = numpy.sqrt(2)* alpha -method = '0' +method = '1' if method == '0': @@ -83,7 +83,7 @@ else: f2 = beta * MixedL21Norm() f = BlockFunction(f1, f2) - g = BlockFunction(0.5 * L1Norm(b=noisy_data), ZeroFunction()) + g = BlockFunction(L1Norm(b=noisy_data), ZeroFunction()) ## Compute operator Norm normK = operator.norm() @@ -122,75 +122,73 @@ plt.title('Middle Line Profiles') plt.show() +#%% Check with CVX solution +from ccpi.optimisation.operators import SparseFiniteDiff +try: + from cvxpy import * + cvx_not_installable = True +except ImportError: + cvx_not_installable = False +if cvx_not_installable: + + u = Variable(ig.shape) + w1 = Variable((N, N)) + w2 = Variable((N, N)) + + # create TGV regulariser + DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') + DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') + + regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u) - vec(w1), \ + DY.matrix() * vec(u) - vec(w2)]), 2, axis = 0)) + \ + beta * sum(norm(vstack([ DX.matrix().transpose() * vec(w1), DY.matrix().transpose() * vec(w2), \ + 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ), \ + 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ) ]), 2, axis = 0 ) ) + + constraints = [] + fidelity = pnorm(u - noisy_data.as_array(),1) + solver = MOSEK + + # choose solver + if 'MOSEK' in installed_solvers(): + solver = MOSEK + else: + solver = SCS + + obj = Minimize( regulariser + fidelity) + prob = Problem(obj) + result = prob.solve(verbose = True, solver = solver) + + diff_cvx = numpy.abs( pdhg.get_output()[0].as_array() - u.value ) + + plt.figure(figsize=(15,15)) + plt.subplot(3,1,1) + plt.imshow(pdhg.get_output()[0].as_array()) + plt.title('PDHG solution') + plt.colorbar() + plt.subplot(3,1,2) + plt.imshow(u.value) + plt.title('CVX solution') + plt.colorbar() + plt.subplot(3,1,3) + plt.imshow(diff_cvx) + plt.title('Difference') + plt.colorbar() + plt.show() + + plt.plot(np.linspace(0,N,N), pdhg.get_output()[0].as_array()[int(N/2),:], label = 'PDHG') + plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') + plt.legend() + plt.title('Middle Line Profiles') + plt.show() + + print('Primal Objective (CVX) {} '.format(obj.value)) + print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) -#%% Check with CVX solution - -#from ccpi.optimisation.operators import SparseFiniteDiff -# -#try: -# from cvxpy import * -# cvx_not_installable = True -#except ImportError: -# cvx_not_installable = False -# -#if cvx_not_installable: -# -# u = Variable(ig.shape) -# w1 = Variable((N, N)) -# w2 = Variable((N, N)) -# -# # create TGV regulariser -# DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') -# DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') -# -# regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u) - vec(w1), \ -# DY.matrix() * vec(u) - vec(w2)]), 2, axis = 0)) + \ -# beta * sum(norm(vstack([ DX.matrix().transpose() * vec(w1), DY.matrix().transpose() * vec(w2), \ -# 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ), \ -# 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ) ]), 2, axis = 0 ) ) -# -# constraints = [] -# fidelity = 0.5 * sum_squares(u - noisy_data.as_array()) -# solver = MOSEK -# -# obj = Minimize( regulariser + fidelity) -# prob = Problem(obj) -# result = prob.solve(verbose = True, solver = solver) -# -# diff_cvx = numpy.abs( res[0].as_array() - u.value ) -# -# # Show result -# plt.figure(figsize=(15,15)) -# plt.subplot(3,1,1) -# plt.imshow(res[0].as_array()) -# plt.title('PDHG solution') -# plt.colorbar() -# -# plt.subplot(3,1,2) -# plt.imshow(u.value) -# plt.title('CVX solution') -# plt.colorbar() -# -# plt.subplot(3,1,3) -# plt.imshow(diff_cvx) -# plt.title('Difference') -# plt.colorbar() -# plt.show() -# -# plt.plot(np.linspace(0,N,N), res[0].as_array()[int(N/2),:], label = 'PDHG') -# plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') -# plt.legend() -# -# print('Primal Objective (CVX) {} '.format(obj.value)) -# print('Primal Objective (PDHG) {} '.format(primal[-1])) -# print('Min/Max of absolute difference {}/{}'.format(diff_cvx.min(), diff_cvx.max())) -# - - diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py index 10e5621..1a3e0df 100644 --- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py @@ -48,7 +48,7 @@ noisy_data = ImageData(n1) # Regularisation Parameter alpha = 0.5 -method = '1' +method = '0' if method == '0': diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py index f2c2023..482b3b4 100644 --- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py @@ -32,7 +32,7 @@ from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \ from skimage.util import random_noise # Create phantom for TV Poisson denoising -N = 300 +N = 100 data = np.zeros((N,N)) data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 @@ -48,7 +48,7 @@ noisy_data = ImageData(n1) # Regularisation Parameter alpha = 2 -method = '0' +method = '1' if method == '0': diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py index f96acd5..484fe42 100644 --- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py @@ -32,7 +32,7 @@ from ccpi.optimisation.functions import ZeroFunction, L1Norm, \ from skimage.util import random_noise # Create phantom for TV Salt & Pepper denoising -N = 300 +N = 100 data = np.zeros((N,N)) data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 @@ -48,7 +48,7 @@ noisy_data = ImageData(n1) # Regularisation Parameter alpha = 2 -method = '0' +method = '1' if method == '0': -- cgit v1.2.3 From f1da1e2bda23a978c92b91a5d533c7d4b331acc3 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 25 Apr 2019 10:08:08 +0100 Subject: change verbose format --- .../Python/ccpi/optimisation/algorithms/Algorithm.py | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py index ed95c3f..3c97480 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -145,13 +145,27 @@ class Algorithm(object): if self.should_stop(): print ("Stop cryterion has been reached.") i = 0 + + print("Iteration {:<5} Primal {:<5} Dual {:<5} PDgap".format('','','')) for _ in self: + + if verbose and self.iteration % self.update_objective_interval == 0: - print ("Iteration {}/{}, objective {}".format(self.iteration, - self.max_iteration, self.get_last_objective()) ) + #pass + print( "{}/{} {:<5} {:.4f} {:<5} {:.4f} {:<5} {:.4f}".\ + format(self.iteration, self.max_iteration,'', \ + self.get_last_objective()[0],'',\ + self.get_last_objective()[1],'',\ + self.get_last_objective()[2])) + + + #print ("Iteration {}/{}, Primal, Dual, PDgap = {}".format(self.iteration, + # self.max_iteration, self.get_last_objective()) ) + + else: if callback is not None: - callback(self.iteration, self.get_last_objective()) + callback(self.iteration, self.get_last_objective(), self.x) i += 1 if i == iterations: break -- cgit v1.2.3 From e021cf74ab461bc72ecc42543852d119bfcd4549 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 25 Apr 2019 10:09:23 +0100 Subject: add dot method --- Wrappers/Python/ccpi/framework/BlockDataContainer.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py index 823a1bd..166014b 100755 --- a/Wrappers/Python/ccpi/framework/BlockDataContainer.py +++ b/Wrappers/Python/ccpi/framework/BlockDataContainer.py @@ -443,10 +443,10 @@ class BlockDataContainer(object): '''Inline truedivision''' return self.__idiv__(other) -# def dot(self, other): + def dot(self, other): # -# tmp = [ self.containers[i].dot(other.containers[i]) for i in range(self.shape[0])] -# return sum(tmp) + tmp = [ self.containers[i].dot(other.containers[i]) for i in range(self.shape[0])] + return sum(tmp) -- cgit v1.2.3 From f6d46771664363b4448c53c8399c3fe1da425f59 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 25 Apr 2019 10:09:52 +0100 Subject: Symmetric BlockGeometry --- Wrappers/Python/ccpi/framework/BlockGeometry.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/framework/BlockGeometry.py b/Wrappers/Python/ccpi/framework/BlockGeometry.py index 7fc5cb8..ed44d99 100755 --- a/Wrappers/Python/ccpi/framework/BlockGeometry.py +++ b/Wrappers/Python/ccpi/framework/BlockGeometry.py @@ -37,16 +37,15 @@ class BlockGeometry(object): containers = [geom.allocate(value) for geom in self.geometries] if symmetry == True: - - # TODO works but needs better coding - + # for 2x2 # [ ig11, ig12\ # ig21, ig22] + # Row-wise Order if len(containers)==4: - containers[1] = containers[2] + containers[1]=containers[2] # for 3x3 # [ ig11, ig12, ig13\ -- cgit v1.2.3 From e7bfeab8ef7c723022d0f36a3be945f6b8a056c8 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 25 Apr 2019 11:17:50 +0100 Subject: Tikhonov demos --- .../Python/wip/Demos/PDHG_Tikhonov_Denoising.py | 176 +++++++++++++++++++++ Wrappers/Python/wip/Demos/PDHG_Tikhonov_Tomo2D.py | 108 +++++++++++++ 2 files changed, 284 insertions(+) create mode 100644 Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py create mode 100644 Wrappers/Python/wip/Demos/PDHG_Tikhonov_Tomo2D.py (limited to 'Wrappers') diff --git a/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py b/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py new file mode 100644 index 0000000..3f275e2 --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py @@ -0,0 +1,176 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.framework import ImageData, ImageGeometry + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG + +from ccpi.optimisation.operators import BlockOperator, Identity, Gradient +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, BlockFunction + +from skimage.util import random_noise + +# Create phantom for TV Salt & Pepper denoising +N = 100 + +data = np.zeros((N,N)) +data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 +data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 +data = ImageData(data) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +ag = ig + +# Create noisy data. Apply Salt & Pepper noise +n1 = random_noise(data.as_array(), mode = 'gaussian', mean=0, var = 0.05, seed=10) +noisy_data = ImageData(n1) + +# Regularisation Parameter +alpha = 4 + +method = '1' + +if method == '0': + + # Create operators + op1 = Gradient(ig) + op2 = Identity(ig, ag) + + # Create BlockOperator + operator = BlockOperator(op1, op2, shape=(2,1) ) + + # Create functions + + f1 = alpha * L2NormSquared() + f2 = 0.5 * L2NormSquared(b = noisy_data) + f = BlockFunction(f1, f2) + g = ZeroFunction() + +else: + + # Without the "Block Framework" + operator = Gradient(ig) + f = alpha * L2NormSquared() + g = 0.5 * L2NormSquared(b = noisy_data) + + +# Compute operator Norm +normK = operator.norm() + +# Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) +opt = {'niter':2000, 'memopt': True} + +# Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 50 +pdhg.run(2000) + + +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(pdhg.get_output().as_array()) +plt.title('Tikhonov Reconstruction') +plt.colorbar() +plt.show() +## +plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') +plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'Tikhonov reconstruction') +plt.legend() +plt.title('Middle Line Profiles') +plt.show() + + +##%% Check with CVX solution + +from ccpi.optimisation.operators import SparseFiniteDiff + +try: + from cvxpy import * + cvx_not_installable = True +except ImportError: + cvx_not_installable = False + + +if cvx_not_installable: + + ##Construct problem + u = Variable(ig.shape) + + DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') + DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') + + # Define Total Variation as a regulariser + + regulariser = alpha * sum_squares(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) + fidelity = 0.5 * sum_squares(u - noisy_data.as_array()) + + # choose solver + if 'MOSEK' in installed_solvers(): + solver = MOSEK + else: + solver = SCS + + obj = Minimize( regulariser + fidelity) + prob = Problem(obj) + result = prob.solve(verbose = True, solver = solver) + + diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value ) + + plt.figure(figsize=(15,15)) + plt.subplot(3,1,1) + plt.imshow(pdhg.get_output().as_array()) + plt.title('PDHG solution') + plt.colorbar() + plt.subplot(3,1,2) + plt.imshow(u.value) + plt.title('CVX solution') + plt.colorbar() + plt.subplot(3,1,3) + plt.imshow(diff_cvx) + plt.title('Difference') + plt.colorbar() + plt.show() + + plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') + plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') + plt.legend() + plt.title('Middle Line Profiles') + plt.show() + + print('Primal Objective (CVX) {} '.format(obj.value)) + print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) + + + + + diff --git a/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Tomo2D.py b/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Tomo2D.py new file mode 100644 index 0000000..5c03362 --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Tomo2D.py @@ -0,0 +1,108 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG + +from ccpi.optimisation.operators import BlockOperator, Gradient +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, BlockFunction +from skimage.util import random_noise +from ccpi.astra.ops import AstraProjectorSimple + +# Create phantom for TV 2D tomography +N = 75 +x = np.zeros((N,N)) +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 + +data = ImageData(x) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) + +detectors = N +angles = np.linspace(0, np.pi, N, dtype=np.float32) + +ag = AcquisitionGeometry('parallel','2D',angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +# Create noisy data. Apply Gaussian noise + +np.random.seed(10) +noisy_data = sin + AcquisitionData(np.random.normal(0, 3, sin.shape)) + +# Regularisation Parameter +alpha = 500 + +# Create operators +op1 = Gradient(ig) +op2 = Aop + +# Create BlockOperator +operator = BlockOperator(op1, op2, shape=(2,1) ) + +# Create functions + +f1 = alpha * L2NormSquared() +f2 = 0.5 * L2NormSquared(b=noisy_data) +f = BlockFunction(f1, f2) + +g = ZeroFunction() + +# Compute operator Norm +normK = operator.norm() + +# Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) + + +# Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 5000 +pdhg.update_objective_interval = 50 +pdhg.run(2000) + +#%% +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(pdhg.get_output().as_array()) +plt.title('Tikhonov Reconstruction') +plt.colorbar() +plt.show() +## +plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') +plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'Tikhonov reconstruction') +plt.legend() +plt.title('Middle Line Profiles') +plt.show() + + -- cgit v1.2.3 From b79143152345bc0e2c49320bc1d7607cedda51cc Mon Sep 17 00:00:00 2001 From: Vaggelis Papoutsellis <22398586+epapoutsellis@users.noreply.github.com> Date: Thu, 25 Apr 2019 11:21:08 +0100 Subject: Delete sinogram_demo_tomography2D.npy --- .../Python/wip/Demos/sinogram_demo_tomography2D.npy | Bin 60128 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 Wrappers/Python/wip/Demos/sinogram_demo_tomography2D.npy (limited to 'Wrappers') diff --git a/Wrappers/Python/wip/Demos/sinogram_demo_tomography2D.npy b/Wrappers/Python/wip/Demos/sinogram_demo_tomography2D.npy deleted file mode 100644 index f37fd4b..0000000 Binary files a/Wrappers/Python/wip/Demos/sinogram_demo_tomography2D.npy and /dev/null differ -- cgit v1.2.3 From d8cbae5c0ed33ea0c30f1e0f3518e2e97c0a86ba Mon Sep 17 00:00:00 2001 From: Vaggelis Papoutsellis <22398586+epapoutsellis@users.noreply.github.com> Date: Thu, 25 Apr 2019 11:21:24 +0100 Subject: Delete untitled4.py --- Wrappers/Python/wip/Demos/untitled4.py | 87 ---------------------------------- 1 file changed, 87 deletions(-) delete mode 100644 Wrappers/Python/wip/Demos/untitled4.py (limited to 'Wrappers') diff --git a/Wrappers/Python/wip/Demos/untitled4.py b/Wrappers/Python/wip/Demos/untitled4.py deleted file mode 100644 index 0cacbd7..0000000 --- a/Wrappers/Python/wip/Demos/untitled4.py +++ /dev/null @@ -1,87 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Wed Apr 24 14:21:08 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData -import numpy -import numpy as np -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction - -from ccpi.astra.ops import AstraProjectorSimple, AstraProjector3DSimple -from skimage.util import random_noise -from timeit import default_timer as timer - -#N = 75 -#x = np.zeros((N,N)) -#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 - -#data = ImageData(x) - -N = 75 -#x = np.zeros((N,N)) - -vert = 4 -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, voxel_num_z=vert) -data = ig.allocate() -Phantom = data -# Populate image data by looping over and filling slices -i = 0 -while i < vert: - if vert > 1: - x = Phantom.subset(vertical=i).array - else: - x = Phantom.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)] = 0.98 - if vert > 1 : - Phantom.fill(x, vertical=i) - i += 1 - -angles_num = 100 -det_w = 1.0 -det_num = N - -angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*\ - 180/np.pi - -# Inputs: Geometry, 2D or 3D, angles, horz detector pixel count, -# horz detector pixel size, vert detector pixel count, -# vert detector pixel size. -ag = AcquisitionGeometry('parallel', - '3D', - angles, - N, - det_w, - vert, - det_w) - -sino = numpy.load("sinogram_demo_tomography2D.npy", mmap_mode='r') -plt.imshow(sino) -plt.title('Sinogram CCPi') -plt.colorbar() -plt.show() - -#%% -Aop = AstraProjector3DSimple(ig, ag) -sin = Aop.direct(data) - -plt.imshow(sin.as_array()) - -plt.title('Sinogram Astra') -plt.colorbar() -plt.show() - - - -#%% \ No newline at end of file -- cgit v1.2.3 From 3be687f3d78b2edcbfec19bb24c3cd0493e7259a Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 25 Apr 2019 11:22:14 +0100 Subject: TV, TGV Tomo --- Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py | 124 ++++++++++++++++ Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py | 211 +++++++++++++++++++++++++++ 2 files changed, 335 insertions(+) create mode 100644 Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py create mode 100644 Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py (limited to 'Wrappers') diff --git a/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py b/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py new file mode 100644 index 0000000..26578bb --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py @@ -0,0 +1,124 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG + +from ccpi.optimisation.operators import BlockOperator, Gradient, Identity, \ + SymmetrizedGradient, ZeroOperator +from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \ + MixedL21Norm, BlockFunction + +from ccpi.astra.ops import AstraProjectorSimple + +# Create phantom for TV 2D tomography +N = 75 + +data = np.zeros((N,N)) + +x1 = np.linspace(0, int(N/2), N) +x2 = np.linspace(int(N/2), 0., N) +xv, yv = np.meshgrid(x1, x2) + +xv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1] = yv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1].T +data = xv +data = ImageData(data/data.max()) + +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) + +detectors = N +angles = np.linspace(0, np.pi, N, dtype=np.float32) + +ag = AcquisitionGeometry('parallel','2D',angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +# Create noisy data. Apply Poisson noise +scale = 0.1 +np.random.seed(5) +n1 = scale * np.random.poisson(sin.as_array()/scale) +noisy_data = AcquisitionData(n1, ag) + + +plt.imshow(noisy_data.as_array()) +plt.show() +#%% +# Regularisation Parameters +alpha = 0.7 +beta = 2 + +# Create Operators +op11 = Gradient(ig) +op12 = Identity(op11.range_geometry()) + +op22 = SymmetrizedGradient(op11.domain_geometry()) +op21 = ZeroOperator(ig, op22.range_geometry()) + +op31 = Aop +op32 = ZeroOperator(op22.domain_geometry(), ag) + +operator = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) + +f1 = alpha * MixedL21Norm() +f2 = beta * MixedL21Norm() +f3 = KullbackLeibler(noisy_data) +f = BlockFunction(f1, f2, f3) +g = ZeroFunction() + +# Compute operator Norm +normK = operator.norm() + +# Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) + + +# Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 50 +pdhg.run(2000) + +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(pdhg.get_output()[0].as_array()) +plt.title('TGV Reconstruction') +plt.colorbar() +plt.show() +## +plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') +plt.plot(np.linspace(0,N,N), pdhg.get_output()[0].as_array()[int(N/2),:], label = 'TGV reconstruction') +plt.legend() +plt.title('Middle Line Profiles') +plt.show() + + diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py b/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py new file mode 100644 index 0000000..0711e91 --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py @@ -0,0 +1,211 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG + +from ccpi.optimisation.operators import BlockOperator, Identity, Gradient +from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \ + MixedL21Norm, BlockFunction + +from ccpi.astra.ops import AstraProjectorSimple + +# Create phantom for TV 2D tomography +N = 75 +x = np.zeros((N,N)) +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 + +data = ImageData(x) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) + +detectors = N +angles = np.linspace(0, np.pi, N, dtype=np.float32) + +ag = AcquisitionGeometry('parallel','2D',angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +# Create noisy data. Apply Poisson noise +scale = 2 +n1 = scale * np.random.poisson(sin.as_array()/scale) +noisy_data = AcquisitionData(n1, ag) + +# Regularisation Parameter +alpha = 5 + +# Create operators +op1 = Gradient(ig) +op2 = Aop + +# Create BlockOperator +operator = BlockOperator(op1, op2, shape=(2,1) ) + +# Create functions + +f1 = alpha * MixedL21Norm() +f2 = KullbackLeibler(noisy_data) +f = BlockFunction(f1, f2) + +g = ZeroFunction() + +# Compute operator Norm +normK = operator.norm() + +# Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) + + +# Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 50 +pdhg.run(2000) + +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(pdhg.get_output().as_array()) +plt.title('TV Reconstruction') +plt.colorbar() +plt.show() +## +plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') +plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction') +plt.legend() +plt.title('Middle Line Profiles') +plt.show() + + +#%% Check with CVX solution + +from ccpi.optimisation.operators import SparseFiniteDiff +import astra +import numpy + +try: + from cvxpy import * + cvx_not_installable = True +except ImportError: + cvx_not_installable = False + + +if cvx_not_installable: + + + ##Construct problem + u = Variable(N*N) + #q = Variable() + + DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') + DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') + + regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) + + # create matrix representation for Astra operator + + vol_geom = astra.create_vol_geom(N, N) + proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) + + proj_id = astra.create_projector('strip', proj_geom, vol_geom) + + matrix_id = astra.projector.matrix(proj_id) + + ProjMat = astra.matrix.get(matrix_id) + + fidelity = sum( ProjMat * u - noisy_data.as_array().ravel() * log(ProjMat * u)) + #constraints = [q>= fidelity, u>=0] + constraints = [u>=0] + + solver = SCS + obj = Minimize( regulariser + fidelity) + prob = Problem(obj, constraints) + result = prob.solve(verbose = True, solver = solver) + + +##%% Check with CVX solution + +from ccpi.optimisation.operators import SparseFiniteDiff + +try: + from cvxpy import * + cvx_not_installable = True +except ImportError: + cvx_not_installable = False + + +if cvx_not_installable: + + ##Construct problem + u = Variable(ig.shape) + + DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') + DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') + + # Define Total Variation as a regulariser + regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) + fidelity = pnorm( u - noisy_data.as_array(),1) + + # choose solver + if 'MOSEK' in installed_solvers(): + solver = MOSEK + else: + solver = SCS + + obj = Minimize( regulariser + fidelity) + prob = Problem(obj) + result = prob.solve(verbose = True, solver = solver) + + + plt.figure(figsize=(15,15)) + plt.subplot(3,1,1) + plt.imshow(pdhg.get_output().as_array()) + plt.title('PDHG solution') + plt.colorbar() + plt.subplot(3,1,2) + plt.imshow(np.reshape(u.value, (N, N))) + plt.title('CVX solution') + plt.colorbar() + plt.subplot(3,1,3) + plt.imshow(diff_cvx) + plt.title('Difference') + plt.colorbar() + plt.show() + + plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') + plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') + plt.legend() + plt.title('Middle Line Profiles') + plt.show() + + print('Primal Objective (CVX) {} '.format(obj.value)) + print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) \ No newline at end of file -- cgit v1.2.3 From 36c36aa4395eb7625b28180bdd6bd376ae2017a7 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 25 Apr 2019 13:10:03 +0100 Subject: add 3D TV denoising --- .../wip/Demos/PDHG_TV_Denoising_Gaussian_3D.py | 155 +++++++++++++++++++++ 1 file changed, 155 insertions(+) create mode 100644 Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_3D.py (limited to 'Wrappers') diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_3D.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_3D.py new file mode 100644 index 0000000..c86ddc9 --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_3D.py @@ -0,0 +1,155 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.framework import ImageData, ImageGeometry + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG + +from ccpi.optimisation.operators import BlockOperator, Identity, Gradient +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ + MixedL21Norm, BlockFunction + +from skimage.util import random_noise + +# Create phantom for TV Gaussian denoising +import timeit +import os +from tomophantom import TomoP3D +import tomophantom + +print ("Building 3D phantom using TomoPhantom software") +tic=timeit.default_timer() +model = 13 # select a model number from the library +N = 64 # Define phantom dimensions using a scalar value (cubic phantom) +path = os.path.dirname(tomophantom.__file__) +path_library3D = os.path.join(path, "Phantom3DLibrary.dat") + +#This will generate a N x N x N phantom (3D) +phantom_tm = TomoP3D.Model(model, N, path_library3D) + +# Create noisy data. Add Gaussian noise +ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N, voxel_num_z=N) +ag = ig +n1 = random_noise(phantom_tm, mode = 'gaussian', mean=0, var = 0.001, seed=10) +noisy_data = ImageData(n1) + +sliceSel = int(0.5*N) +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(noisy_data.as_array()[sliceSel,:,:],vmin=0, vmax=1) +plt.title('Axial View') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(noisy_data.as_array()[:,sliceSel,:],vmin=0, vmax=1) +plt.title('Coronal View') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(noisy_data.as_array()[:,:,sliceSel],vmin=0, vmax=1) +plt.title('Sagittal View') +plt.colorbar() +plt.show() + + +# Regularisation Parameter +alpha = 0.05 + +method = '0' + +if method == '0': + + # Create operators + op1 = Gradient(ig) + op2 = Identity(ig, ag) + + # Create BlockOperator + operator = BlockOperator(op1, op2, shape=(2,1) ) + + # Create functions + + f1 = alpha * MixedL21Norm() + f2 = 0.5 * L2NormSquared(b = noisy_data) + f = BlockFunction(f1, f2) + + g = ZeroFunction() + +else: + + # Without the "Block Framework" + operator = Gradient(ig) + f = alpha * MixedL21Norm() + g = 0.5 * L2NormSquared(b = noisy_data) + + +# Compute operator Norm +normK = operator.norm() + +# Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) + + +# Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 200 +pdhg.run(2000) + +fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8)) + +plt.subplot(2,3,1) +plt.imshow(noisy_data.as_array()[sliceSel,:,:],vmin=0, vmax=1) +plt.axis('off') +plt.title('Axial View') + +plt.subplot(2,3,2) +plt.imshow(noisy_data.as_array()[:,sliceSel,:],vmin=0, vmax=1) +plt.axis('off') +plt.title('Coronal View') + +plt.subplot(2,3,3) +plt.imshow(noisy_data.as_array()[:,:,sliceSel],vmin=0, vmax=1) +plt.axis('off') +plt.title('Sagittal View') + + +plt.subplot(2,3,4) +plt.imshow(pdhg.get_output().as_array()[sliceSel,:,:],vmin=0, vmax=1) +plt.axis('off') +plt.subplot(2,3,5) +plt.imshow(pdhg.get_output().as_array()[:,sliceSel,:],vmin=0, vmax=1) +plt.axis('off') +plt.subplot(2,3,6) +plt.imshow(pdhg.get_output().as_array()[:,:,sliceSel],vmin=0, vmax=1) +plt.axis('off') +im = plt.imshow(pdhg.get_output().as_array()[:,:,sliceSel],vmin=0, vmax=1) + + +fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8, + wspace=0.02, hspace=0.02) + +cb_ax = fig.add_axes([0.83, 0.1, 0.02, 0.8]) +cbar = fig.colorbar(im, cax=cb_ax) + + +plt.show() + -- cgit v1.2.3 From d34028a672b22f7c4a02464736d74c70fb354362 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 25 Apr 2019 16:23:57 +0100 Subject: memopt fix prox conjugate --- .../ccpi/optimisation/functions/KullbackLeibler.py | 39 ++++++---------------- 1 file changed, 11 insertions(+), 28 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py index 22d21fd..14b5ea0 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py @@ -62,19 +62,6 @@ class KullbackLeibler(Function): if out is None: return 1 - self.b/(x + self.bnoise) else: -#<<<<<<< HEAD -# self.b.divide(x+self.bnoise, out=out) -# out.subtract(1, out=out) -# -# def convex_conjugate(self, x): -# -# tmp = self.b/( 1 - x ) -# ind = tmp.as_array()>0 -# -# sel -# -# return (self.b * ( ImageData( numpy.log(tmp) ) - 1 ) - self.bnoise * (x - 1)).sum() -#======= x.add(self.bnoise, out=out) self.b.divide(out, out=out) out.subtract(1, out=out) @@ -116,23 +103,19 @@ class KullbackLeibler(Function): if out is None: z = x + tau * self.bnoise - return 0.5*((z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt()) else: - z = x + tau * self.bnoise - res = 0.5*((z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt()) - out.fill(res) -# else: -# z_m = x + tau * self.bnoise -1 -# self.b.multiply(4*tau, out=out) -# z_m.multiply(z_m, out=z_m) -# out += z_m -# out.sqrt(out=out) -# z = z_m + 2 -# z_m.sqrt(out=z_m) -# z_m += 2 -# out *= -1 -# out += z_m + z_m = x + tau * self.bnoise -1 + self.b.multiply(4*tau, out=out) + z_m.multiply(z_m, out=z_m) + out += z_m + out.sqrt(out=out) + z_m.sqrt(out=z_m) + z_m += 2 + out *= -1 + out += z_m + out *= 0.5 + def __rmul__(self, scalar): -- cgit v1.2.3 From b726c09f16c256daf92544668d2ff51bed3cf06f Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 25 Apr 2019 16:30:33 +0100 Subject: L2Norm fix memopt --- Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py index 5490782..20e754e 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py @@ -94,12 +94,12 @@ class L2NormSquared(Function): return x/(1+2*tau) else: -# tmp = x.subtract(self.b) + tmp = x.subtract(self.b) # tmp -= self.b -# tmp /= (1+2*tau) -# tmp += self.b -# return tmp - return (x-self.b)/(1+2*tau) + self.b + tmp /= (1+2*tau) + tmp += self.b + return tmp +# return (x-self.b)/(1+2*tau) + self.b # if self.b is not None: # out=x -- cgit v1.2.3