summaryrefslogtreecommitdiffstats
path: root/Wrappers
diff options
context:
space:
mode:
authorepapoutsellis <epapoutsellis@gmail.com>2019-04-25 11:17:50 +0100
committerepapoutsellis <epapoutsellis@gmail.com>2019-04-25 11:17:50 +0100
commite7bfeab8ef7c723022d0f36a3be945f6b8a056c8 (patch)
treef482a711cc0e2481f11b2903e742344c955f6f63 /Wrappers
parentf6d46771664363b4448c53c8399c3fe1da425f59 (diff)
downloadframework-e7bfeab8ef7c723022d0f36a3be945f6b8a056c8.tar.gz
framework-e7bfeab8ef7c723022d0f36a3be945f6b8a056c8.tar.bz2
framework-e7bfeab8ef7c723022d0f36a3be945f6b8a056c8.tar.xz
framework-e7bfeab8ef7c723022d0f36a3be945f6b8a056c8.zip
Tikhonov demos
Diffstat (limited to 'Wrappers')
-rw-r--r--Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py176
-rw-r--r--Wrappers/Python/wip/Demos/PDHG_Tikhonov_Tomo2D.py108
2 files changed, 284 insertions, 0 deletions
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()
+
+