From fb6f44dff00c194e4936f39e4be8d0f22e8fe4da Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 25 Apr 2019 17:43:16 +0100 Subject: TV dynamic tomo --- Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D_time.py | 169 +++++++++++++++++++++++ Wrappers/Python/wip/pdhg_TV_tomography2D_time.py | 4 +- 2 files changed, 171 insertions(+), 2 deletions(-) create mode 100644 Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D_time.py (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D_time.py b/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D_time.py new file mode 100644 index 0000000..045458a --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D_time.py @@ -0,0 +1,169 @@ +# -*- 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, KullbackLeibler, \ + MixedL21Norm, BlockFunction + +from ccpi.astra.ops import AstraProjectorMC + +import os +import tomophantom +from tomophantom import TomoP2D + +# Create phantom for TV 2D dynamic tomography + +model = 102 # note that the selected model is temporal (2D + time) +N = 50 # set dimension of the phantom +# one can specify an exact path to the parameters file +# path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' +path = os.path.dirname(tomophantom.__file__) +path_library2D = os.path.join(path, "Phantom2DLibrary.dat") +#This will generate a N_size x N_size x Time frames phantom (2D + time) +phantom_2Dt = TomoP2D.ModelTemporal(model, N, path_library2D) + +plt.close('all') +plt.figure(1) +plt.rcParams.update({'font.size': 21}) +plt.title('{}''{}'.format('2D+t phantom using model no.',model)) +for sl in range(0,np.shape(phantom_2Dt)[0]): + im = phantom_2Dt[sl,:,:] + plt.imshow(im, vmin=0, vmax=1) + plt.pause(.1) + plt.draw + + +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, channels = np.shape(phantom_2Dt)[0]) +data = ImageData(phantom_2Dt, geometry=ig) + +detectors = N +angles = np.linspace(0,np.pi,N) + +ag = AcquisitionGeometry('parallel','2D', angles, detectors, channels = np.shape(phantom_2Dt)[0]) +Aop = AstraProjectorMC(ig, ag, 'gpu') +sin = Aop.direct(data) + +scale = 2 +n1 = scale * np.random.poisson(sin.as_array()/scale) +noisy_data = AcquisitionData(n1, ag) + +tindex = [3, 6, 10] + +fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 10)) +plt.subplot(1,3,1) +plt.imshow(noisy_data.as_array()[tindex[0],:,:]) +plt.axis('off') +plt.title('Time {}'.format(tindex[0])) +plt.subplot(1,3,2) +plt.imshow(noisy_data.as_array()[tindex[1],:,:]) +plt.axis('off') +plt.title('Time {}'.format(tindex[1])) +plt.subplot(1,3,3) +plt.imshow(noisy_data.as_array()[tindex[2],:,:]) +plt.axis('off') +plt.title('Time {}'.format(tindex[2])) + +fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8, + wspace=0.02, hspace=0.02) + +plt.show() + +#%% +# Regularisation Parameter +alpha = 5 + +# Create operators +#op1 = Gradient(ig) +op1 = Gradient(ig, correlation='SpaceChannels') +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 = 200 +pdhg.run(2000) + + +#%% +fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8)) + +plt.subplot(2,3,1) +plt.imshow(phantom_2Dt[tindex[0],:,:],vmin=0, vmax=1) +plt.axis('off') +plt.title('Time {}'.format(tindex[0])) + +plt.subplot(2,3,2) +plt.imshow(phantom_2Dt[tindex[1],:,:],vmin=0, vmax=1) +plt.axis('off') +plt.title('Time {}'.format(tindex[1])) + +plt.subplot(2,3,3) +plt.imshow(phantom_2Dt[tindex[2],:,:],vmin=0, vmax=1) +plt.axis('off') +plt.title('Time {}'.format(tindex[2])) + + +plt.subplot(2,3,4) +plt.imshow(pdhg.get_output().as_array()[tindex[0],:,:]) +plt.axis('off') +plt.subplot(2,3,5) +plt.imshow(pdhg.get_output().as_array()[tindex[1],:,:]) +plt.axis('off') +plt.subplot(2,3,6) +plt.imshow(pdhg.get_output().as_array()[tindex[2],:,:]) +plt.axis('off') +im = plt.imshow(pdhg.get_output().as_array()[tindex[0],:,:]) + + +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() + diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py b/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py index dea8e5c..5423b22 100644 --- a/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py +++ b/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py @@ -16,7 +16,7 @@ import matplotlib.pyplot as plt from ccpi.optimisation.algorithms import PDHG, PDHG_old from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFun, L2NormSquared, \ +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ MixedL21Norm, BlockFunction, ScaledFunction from ccpi.astra.ops import AstraProjectorSimple, AstraProjectorMC @@ -100,7 +100,7 @@ operator = BlockOperator(op1, op2, shape=(2,1) ) alpha = 50 f = BlockFunction( alpha * MixedL21Norm(), \ 0.5 * L2NormSquared(b = noisy_data) ) -g = ZeroFun() +g = ZeroFunction() # Compute operator Norm normK = operator.norm() -- cgit v1.2.3 From 07e7344926bc8850ef5e2c1580fdff580ca8c9a9 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 25 Apr 2019 17:45:46 +0100 Subject: callback fix --- Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py index 482b3b4..32ab62d 100644 --- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py @@ -87,7 +87,19 @@ opt = {'niter':2000, 'memopt': True} 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) + +def pdgap_print(niter, objective, solution): + + + print( "{:04}/{:04} {:<5} {:.4f} {:<5} {:.4f} {:<5} {:.4f}".\ + format(niter, pdhg.max_iteration,'', \ + objective[0],'',\ + objective[1],'',\ + objective[2])) + +#pdhg.run(2000) + +pdhg.run(2000, callback = pdgap_print) plt.figure(figsize=(15,15)) -- cgit v1.2.3 From 8ef753231e74b4dad339370661b563a57ffe75cf Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 25 Apr 2019 17:48:03 +0100 Subject: fix L2 and ALgo --- .../ccpi/optimisation/algorithms/Algorithm.py | 25 ++++++++++++++-------- .../ccpi/optimisation/functions/L2NormSquared.py | 8 ++++++- 2 files changed, 23 insertions(+), 10 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py index 3c97480..47376a5 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -149,23 +149,30 @@ class Algorithm(object): print("Iteration {:<5} Primal {:<5} Dual {:<5} PDgap".format('','','')) for _ in self: + + if self.iteration % self.update_objective_interval == 0: + if verbose: - if verbose and self.iteration % self.update_objective_interval == 0: +# if verbose and self.iteration % self.update_objective_interval == 0: #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])) + # \t for tab +# print( "{:04}/{:04} {:<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 {}/{}, {}".format(self.iteration, + self.max_iteration, self.get_last_objective()) ) + #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(), self.x) + else: + if callback is not None: + callback(self.iteration, self.get_last_objective(), self.x) i += 1 if i == iterations: break diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py index 20e754e..6cb3d8a 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py @@ -146,7 +146,13 @@ class L2NormSquared(Function): ''' - return ScaledFunction(self, scalar) + return ScaledFunction(self, scalar) + + + def operator_composition(self, operator): + + return FunctionOperatorComposition + if __name__ == '__main__': -- cgit v1.2.3 From 5ae13b1d55da87f4c3f3908fc91ec2424daaf4b3 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Mon, 29 Apr 2019 09:43:18 +0100 Subject: changes to PD aglo --- .../ccpi/optimisation/algorithms/Algorithm.py | 22 ++-- .../Python/wip/Demos/PDHG_TV_Denoising_Poisson.py | 127 ++++++++++----------- 2 files changed, 76 insertions(+), 73 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py index 47376a5..12cbabc 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -146,12 +146,18 @@ class Algorithm(object): print ("Stop cryterion has been reached.") i = 0 - print("Iteration {:<5} Primal {:<5} Dual {:<5} PDgap".format('','','')) +# print("Iteration {:<5} Primal {:<5} Dual {:<5} PDgap".format('','','')) for _ in self: - + if self.iteration % self.update_objective_interval == 0: - if verbose: + + if callback is not None: + callback(self.iteration, self.get_last_objective(), self.x) + + else: + + if verbose: # if verbose and self.iteration % self.update_objective_interval == 0: #pass @@ -163,16 +169,16 @@ class Algorithm(object): # self.get_last_objective()[2])) - print ("Iteration {}/{}, {}".format(self.iteration, - self.max_iteration, self.get_last_objective()) ) + print ("Iteration {}/{}, {}".format(self.iteration, + self.max_iteration, self.get_last_objective()) ) #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(), self.x) +# else: +# if callback is not None: +# callback(self.iteration, self.get_last_objective(), self.x) i += 1 if i == iterations: break diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py index 32ab62d..ccdabb2 100644 --- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py @@ -88,7 +88,7 @@ pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) pdhg.max_iteration = 2000 pdhg.update_objective_interval = 50 -def pdgap_print(niter, objective, solution): +def pdgap_objectives(niter, objective, solution): print( "{:04}/{:04} {:<5} {:.4f} {:<5} {:.4f} {:<5} {:.4f}".\ @@ -97,9 +97,7 @@ def pdgap_print(niter, objective, solution): objective[1],'',\ objective[2])) -#pdhg.run(2000) - -pdhg.run(2000, callback = pdgap_print) +pdhg.run(2000, callback = pdgap_objectives) plt.figure(figsize=(15,15)) @@ -124,66 +122,65 @@ 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])) - - - - - +#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])) +# +# +# +# +# -- cgit v1.2.3