From fb6f44dff00c194e4936f39e4be8d0f22e8fe4da Mon Sep 17 00:00:00 2001
From: epapoutsellis <epapoutsellis@gmail.com>
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')
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 <epapoutsellis@gmail.com>
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')
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 <epapoutsellis@gmail.com>
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')
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 <epapoutsellis@gmail.com>
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')
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