diff options
| -rw-r--r-- | Wrappers/Python/demos/PDHG_Tikhonov_Tomo2D.py | 12 | ||||
| -rw-r--r-- | Wrappers/Python/demos/PDHG_examples/PDHG_Tikhonov_Tomo2D.py | 156 |
2 files changed, 162 insertions, 6 deletions
diff --git a/Wrappers/Python/demos/PDHG_Tikhonov_Tomo2D.py b/Wrappers/Python/demos/PDHG_Tikhonov_Tomo2D.py index 22972da..02cd053 100644 --- a/Wrappers/Python/demos/PDHG_Tikhonov_Tomo2D.py +++ b/Wrappers/Python/demos/PDHG_Tikhonov_Tomo2D.py @@ -47,8 +47,8 @@ import matplotlib.pyplot as plt from ccpi.optimisation.algorithms import PDHG from ccpi.optimisation.operators import BlockOperator, Gradient -from ccpi.optimisation.functions import IndicatorBox, L2NormSquared, BlockFunction -from skimage.util import random_noise +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, BlockFunction + from ccpi.astra.ops import AstraProjectorSimple from ccpi.framework import TestData import os, sys @@ -101,7 +101,7 @@ plt.show() # Regularisation Parameter -alpha = 500 +alpha = 1000 # Create operators op1 = Gradient(ig) @@ -115,8 +115,8 @@ operator = BlockOperator(op1, op2, shape=(2,1) ) f1 = alpha * L2NormSquared() f2 = 0.5 * L2NormSquared(b=noisy_data) f = BlockFunction(f1, f2) - -g = IndicatorBox(lower=0) + +g = ZeroFunction() # Compute operator Norm normK = operator.norm() @@ -127,7 +127,7 @@ 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 = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) pdhg.max_iteration = 2000 pdhg.update_objective_interval = 500 pdhg.run(2000) diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_Tikhonov_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/PDHG_Tikhonov_Tomo2D.py new file mode 100644 index 0000000..02cd053 --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/PDHG_Tikhonov_Tomo2D.py @@ -0,0 +1,156 @@ +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# 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.txt +# +# 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. +# +#========================================================================= + +""" + +Total Variation Denoising using PDHG algorithm: + +Problem: min_x, x>0 \alpha * ||\nabla x||_{2}^{2} + int A x -g log(Ax + \eta) + + \nabla: Gradient operator + + A: Projection Matrix + g: Noisy sinogram corrupted with Poisson Noise + + \eta: Background Noise + \alpha: Regularization parameter + + + +""" + + +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 ccpi.astra.ops import AstraProjectorSimple +from ccpi.framework import TestData +import os, sys + +loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) + +# Load Data +N = 100 +M = 100 +data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1)) + +ig = data.geometry +ag = ig + +#Create Acquisition Data and apply poisson noise + +detectors = N +angles = np.linspace(0, np.pi, N) + +ag = AcquisitionGeometry('parallel','2D',angles, detectors) + +device = input('Available device: GPU==1 / CPU==0 ') + +if device=='1': + dev = 'gpu' +else: + dev = 'cpu' + +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +# Create noisy data. Apply Poisson noise +scale = 0.5 +eta = 0 +n1 = scale * np.random.poisson(eta + sin.as_array()/scale) + +noisy_data = AcquisitionData(n1, ag) + +# Show Ground Truth and Noisy Data +plt.figure(figsize=(10,10)) +plt.subplot(2,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(2,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.show() + + +# Regularisation Parameter +alpha = 1000 + +# 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) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 500 +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,M), data.as_array()[int(N/2),:], label = 'GTruth') +plt.plot(np.linspace(0,N,M), pdhg.get_output().as_array()[int(N/2),:], label = 'Tikhonov reconstruction') +plt.legend() +plt.title('Middle Line Profiles') +plt.show() + + |