diff options
16 files changed, 3134 insertions, 0 deletions
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py new file mode 100755 index 0000000..4d6da00 --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py @@ -0,0 +1,221 @@ +#======================================================================== +# 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 Generalised Variation (TGV) Denoising using PDHG algorithm: + + +Problem: min_{x} \alpha * ||\nabla x - w||_{2,1} + + \beta * || E w ||_{2,1} + + fidelity + + where fidelity can be as follows depending on the noise characteristics + of the data: + * Norm2Squared \frac{1}{2} * || x - g ||_{2}^{2} + * KullbackLeibler \int u - g * log(u) + Id_{u>0} + * L1Norm ||u - g||_{1} + + \alpha: Regularization parameter + \beta: Regularization parameter + + \nabla: Gradient operator + E: Symmetrized Gradient operator + + g: Noisy Data with Salt & Pepper Noise + + Method = 0 ( PDHG - split ) : K = [ \nabla, - Identity + ZeroOperator, E + Identity, ZeroOperator] + + + Method = 1 (PDHG - explicit ): K = [ \nabla, - Identity + ZeroOperator, E ] + +""" + +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, KullbackLeibler, L2NormSquared +import sys +if int(numpy.version.version.split('.')[1]) > 12: + from skimage.util import random_noise +else: + from demoutil import random_noise + + + +# user supplied input +if len(sys.argv) > 1: + which_noise = int(sys.argv[1]) +else: + which_noise = 0 +print ("Applying {} noise") + +if len(sys.argv) > 2: + method = sys.argv[2] +else: + method = '1' +print ("method ", method) +# Create phantom for TGV SaltPepper denoising + +N = 100 + +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 + +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +data = ig.allocate() +data.fill(xv/xv.max()) +ag = ig + +# Create noisy data. +# Apply Salt & Pepper noise +# gaussian +# poisson +noises = ['gaussian', 'poisson', 's&p'] +noise = noises[which_noise] +if noise == 's&p': + n1 = random_noise(data.as_array(), mode = noise, salt_vs_pepper = 0.9, amount=0.2, seed=10) +elif noise == 'poisson': + n1 = random_noise(data.as_array(), mode = noise, seed = 10) +elif noise == 'gaussian': + n1 = random_noise(data.as_array(), mode = noise, seed = 10) +else: + raise ValueError('Unsupported Noise ', noise) +noisy_data = ImageData(n1) + +# Show Ground Truth and Noisy Data +plt.figure(figsize=(10,5)) +plt.subplot(1,2,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(1,2,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.show() + +# Regularisation Parameters +if noise == 's&p': + alpha = 0.8 +elif noise == 'poisson': + alpha = .1 +elif noise == 'gaussian': + alpha = .3 + +beta = numpy.sqrt(2)* alpha + +# fidelity +if noise == 's&p': + f3 = L1Norm(b=noisy_data) +elif noise == 'poisson': + f3 = KullbackLeibler(noisy_data) +elif noise == 'gaussian': + f3 = L2NormSquared(b=noisy_data) + +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 depends on the noise characteristics + + 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(f3, 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 = 200 +pdhg.run(2000, verbose = True) + +#%% +plt.figure(figsize=(20,5)) +plt.subplot(1,4,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(1,4,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.subplot(1,4,3) +plt.imshow(pdhg.get_output()[0].as_array()) +plt.title('TGV Reconstruction') +plt.colorbar() +plt.subplot(1,4,4) +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/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py new file mode 100755 index 0000000..0f1effa --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py @@ -0,0 +1,266 @@ +#======================================================================== +# 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,1} + ||x-g||_{1} + + \alpha: Regularization parameter + + \nabla: Gradient operator + + g: Noisy Data with Salt & Pepper Noise + + + Method = 0 ( PDHG - split ) : K = [ \nabla, + Identity] + + + Method = 1 (PDHG - explicit ): K = \nabla + + +""" + +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, L2NormSquared,\ + KullbackLeibler +from ccpi.framework import TestData +import os, sys +if int(numpy.version.version.split('.')[1]) > 12: + from skimage.util import random_noise +else: + from demoutil import random_noise + +# user supplied input +if len(sys.argv) > 1: + which_noise = int(sys.argv[1]) +else: + which_noise = 0 +print ("Applying {} noise") + +if len(sys.argv) > 2: + method = sys.argv[2] +else: + method = '0' +print ("method ", method) +# Create phantom for TV Salt & Pepper denoising +N = 100 + +loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) +data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,N)) +data = loader.load(TestData.PEPPERS, size=(N,N)) +ig = data.geometry +ag = ig + +# Create noisy data. +# Apply Salt & Pepper noise +# gaussian +# poisson +noises = ['gaussian', 'poisson', 's&p'] +noise = noises[which_noise] +if noise == 's&p': + n1 = random_noise(data.as_array(), mode = noise, salt_vs_pepper = 0.9, amount=0.2) +elif noise == 'poisson': + n1 = random_noise(data.as_array(), mode = noise, seed = 10) +elif noise == 'gaussian': + n1 = random_noise(data.as_array(), mode = noise, seed = 10) +else: + raise ValueError('Unsupported Noise ', noise) +noisy_data = ig.allocate() +noisy_data.fill(n1) +#noisy_data = ImageData(n1) + +# Show Ground Truth and Noisy Data +plt.figure(figsize=(10,5)) +plt.subplot(1,2,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(1,2,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.show() + +# Regularisation Parameter +alpha = .2 + +# fidelity +if noise == 's&p': + f2 = L1Norm(b=noisy_data) +elif noise == 'poisson': + f2 = KullbackLeibler(noisy_data) +elif noise == 'gaussian': + f2 = L2NormSquared(b=noisy_data) + +if method == '0': + + # Create operators + op1 = Gradient(ig, correlation=Gradient.CORRELATION_SPACECHANNEL) + 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) + g = f2 + + +# 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) + +if data.geometry.channels > 1: + plt.figure(figsize=(20,15)) + for row in range(data.geometry.channels): + + plt.subplot(3,4,1+row*4) + plt.imshow(data.subset(channel=row).as_array()) + plt.title('Ground Truth') + plt.colorbar() + plt.subplot(3,4,2+row*4) + plt.imshow(noisy_data.subset(channel=row).as_array()) + plt.title('Noisy Data') + plt.colorbar() + plt.subplot(3,4,3+row*4) + plt.imshow(pdhg.get_output().subset(channel=row).as_array()) + plt.title('TV Reconstruction') + plt.colorbar() + plt.subplot(3,4,4+row*4) + plt.plot(np.linspace(0,N,N), data.subset(channel=row).as_array()[int(N/2),:], label = 'GTruth') + plt.plot(np.linspace(0,N,N), pdhg.get_output().subset(channel=row).as_array()[int(N/2),:], label = 'TV reconstruction') + plt.legend() + plt.title('Middle Line Profiles') + plt.show() +else: + plt.figure(figsize=(20,5)) + plt.subplot(1,4,1) + plt.imshow(data.subset(channel=0).as_array()) + plt.title('Ground Truth') + plt.colorbar() + plt.subplot(1,4,2) + plt.imshow(noisy_data.subset(channel=0).as_array()) + plt.title('Noisy Data') + plt.colorbar() + plt.subplot(1,4,3) + plt.imshow(pdhg.get_output().subset(channel=0).as_array()) + plt.title('TV Reconstruction') + plt.colorbar() + plt.subplot(1,4,4) + 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/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py new file mode 100644 index 0000000..5aa334d --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py @@ -0,0 +1,256 @@ +#======================================================================== +# 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. +# +#========================================================================= +""" + +Tikhonov Denoising using PDHG algorithm: + + +Problem: min_{x} \alpha * ||\nabla x||_{2}^{2} + \frac{1}{2} * || x - g ||_{2}^{2} + + \alpha: Regularization parameter + + \nabla: Gradient operator + + g: Noisy Data with Gaussian Noise + + Method = 0 ( PDHG - split ) : K = [ \nabla, + Identity] + + + Method = 1 (PDHG - explicit ): K = \nabla + +""" + +from ccpi.framework import ImageData, TestData + +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, KullbackLeibler, L1Norm + +import sys, os +if int(numpy.version.version.split('.')[1]) > 12: + from skimage.util import random_noise +else: + from demoutil import random_noise + + +# user supplied input +if len(sys.argv) > 1: + which_noise = int(sys.argv[1]) +else: + which_noise = 0 +print ("Applying {} noise") + +if len(sys.argv) > 2: + method = sys.argv[2] +else: + method = '0' +print ("method ", method) + +# Create phantom for TV Salt & Pepper denoising +N = 100 + +loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) +data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,N)) +ig = data.geometry +ag = ig + +# Create noisy data. Apply Salt & Pepper noise +# Create noisy data. +# Apply Salt & Pepper noise +# gaussian +# poisson +noises = ['gaussian', 'poisson', 's&p'] +noise = noises[which_noise] +if noise == 's&p': + n1 = random_noise(data.as_array(), mode = noise, salt_vs_pepper = 0.9, amount=0.2) +elif noise == 'poisson': + n1 = random_noise(data.as_array(), mode = noise, seed = 10) +elif noise == 'gaussian': + n1 = random_noise(data.as_array(), mode = noise, seed = 10) +else: + raise ValueError('Unsupported Noise ', noise) +noisy_data = ImageData(n1) + +# fidelity +if noise == 's&p': + f2 = L1Norm(b=noisy_data) +elif noise == 'poisson': + f2 = KullbackLeibler(noisy_data) +elif noise == 'gaussian': + f2 = 0.5 * L2NormSquared(b=noisy_data) + +# Show Ground Truth and Noisy Data +plt.figure(figsize=(10,5)) +plt.subplot(1,2,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(1,2,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.show() + + +# Regularisation Parameter +# no edge preservation alpha is big +if noise == 's&p': + alpha = 8. +elif noise == 'poisson': + alpha = 8. +elif noise == 'gaussian': + alpha = 8. + +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 must change according to noise + #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 must change according to noise + #g = 0.5 * L2NormSquared(b = noisy_data) + g = f2 + + +# 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(1500) + + +plt.figure(figsize=(20,5)) +plt.subplot(1,4,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(1,4,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.subplot(1,4,3) +plt.imshow(pdhg.get_output().as_array()) +plt.title('Tikhonov Reconstruction') +plt.colorbar() +plt.subplot(1,4,4) +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_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/demos/PDHG_examples/MultiChannel/PDHG_2D_time_denoising.py b/Wrappers/Python/demos/PDHG_examples/MultiChannel/PDHG_2D_time_denoising.py new file mode 100644 index 0000000..045458a --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/MultiChannel/PDHG_2D_time_denoising.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/demos/PDHG_examples/MultiChannel/PDHG_TV_Tomo2D_time.py b/Wrappers/Python/demos/PDHG_examples/MultiChannel/PDHG_TV_Tomo2D_time.py new file mode 100644 index 0000000..045458a --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/MultiChannel/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/demos/PDHG_examples/TGV_Denoising/PDHG_TGV_Denoising_SaltPepper.py b/Wrappers/Python/demos/PDHG_examples/TGV_Denoising/PDHG_TGV_Denoising_SaltPepper.py new file mode 100644 index 0000000..15c0a05 --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/TGV_Denoising/PDHG_TGV_Denoising_SaltPepper.py @@ -0,0 +1,245 @@ +#======================================================================== +# 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 Generalised Variation (TGV) Denoising using PDHG algorithm: + + +Problem: min_{x} \alpha * ||\nabla x - w||_{2,1} + + \beta * || E w ||_{2,1} + + \frac{1}{2} * || x - g ||_{2}^{2} + + \alpha: Regularization parameter + \beta: Regularization parameter + + \nabla: Gradient operator + E: Symmetrized Gradient operator + + g: Noisy Data with Salt & Pepper Noise + + Method = 0 ( PDHG - split ) : K = [ \nabla, - Identity + ZeroOperator, E + Identity, ZeroOperator] + + + Method = 1 (PDHG - explicit ): K = [ \nabla, - Identity + ZeroOperator, E ] + +""" + +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 = 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()) + +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) + +# Show Ground Truth and Noisy Data +plt.figure(figsize=(15,15)) +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 Parameters +alpha = 0.8 +beta = numpy.sqrt(2)* alpha + +method = '1' + +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(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) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 50 +pdhg.run(2000, verbose = False) + +#%% +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 = 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])) + + + + + diff --git a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/.DS_Store b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/.DS_Store Binary files differnew file mode 100644 index 0000000..5008ddf --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/.DS_Store diff --git a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_2D_time.py b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_2D_time.py new file mode 100644 index 0000000..14608db --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_2D_time.py @@ -0,0 +1,192 @@ +#======================================================================== +# 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. +# +#========================================================================= + +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 +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ + 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 = 128 # 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) +ag = ig + +# Create Noisy data. Add Gaussian noise +np.random.seed(10) +noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.25, size=ig.shape) ) + +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 = 0.3 + +# Create operators +#op1 = Gradient(ig) +op1 = Gradient(ig, correlation='Space') +op2 = Gradient(ig, correlation='SpaceChannels') + +op3 = Identity(ig, ag) + +# Create BlockOperator +operator1 = BlockOperator(op1, op3, shape=(2,1) ) +operator2 = BlockOperator(op2, op3, shape=(2,1) ) + +# Create functions + +f1 = alpha * MixedL21Norm() +f2 = 0.5 * L2NormSquared(b = noisy_data) +f = BlockFunction(f1, f2) + +g = ZeroFunction() + +# Compute operator Norm +normK1 = operator1.norm() +normK2 = operator2.norm() + +#%% +# Primal & dual stepsizes +sigma1 = 1 +tau1 = 1/(sigma1*normK1**2) + +sigma2 = 1 +tau2 = 1/(sigma2*normK2**2) + +# Setup and run the PDHG algorithm +pdhg1 = PDHG(f=f,g=g,operator=operator1, tau=tau1, sigma=sigma1) +pdhg1.max_iteration = 2000 +pdhg1.update_objective_interval = 200 +pdhg1.run(2000) + +# Setup and run the PDHG algorithm +pdhg2 = PDHG(f=f,g=g,operator=operator2, tau=tau2, sigma=sigma2) +pdhg2.max_iteration = 2000 +pdhg2.update_objective_interval = 200 +pdhg2.run(2000) + + +#%% + +tindex = [3, 6, 10] +fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8)) + +plt.subplot(3,3,1) +plt.imshow(phantom_2Dt[tindex[0],:,:]) +plt.axis('off') +plt.title('Time {}'.format(tindex[0])) + +plt.subplot(3,3,2) +plt.imshow(phantom_2Dt[tindex[1],:,:]) +plt.axis('off') +plt.title('Time {}'.format(tindex[1])) + +plt.subplot(3,3,3) +plt.imshow(phantom_2Dt[tindex[2],:,:]) +plt.axis('off') +plt.title('Time {}'.format(tindex[2])) + +plt.subplot(3,3,4) +plt.imshow(pdhg1.get_output().as_array()[tindex[0],:,:]) +plt.axis('off') +plt.subplot(3,3,5) +plt.imshow(pdhg1.get_output().as_array()[tindex[1],:,:]) +plt.axis('off') +plt.subplot(3,3,6) +plt.imshow(pdhg1.get_output().as_array()[tindex[2],:,:]) +plt.axis('off') + + +plt.subplot(3,3,7) +plt.imshow(pdhg2.get_output().as_array()[tindex[0],:,:]) +plt.axis('off') +plt.subplot(3,3,8) +plt.imshow(pdhg2.get_output().as_array()[tindex[1],:,:]) +plt.axis('off') +plt.subplot(3,3,9) +plt.imshow(pdhg2.get_output().as_array()[tindex[2],:,:]) +plt.axis('off') + +#%% +im = plt.imshow(pdhg1.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/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Gaussian.py b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Gaussian.py new file mode 100644 index 0000000..9d00ee1 --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Gaussian.py @@ -0,0 +1,225 @@ +#======================================================================== +# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). + +# Copyright 2017 UKRI-STFC +# Copyright 2017 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 + +# 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} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2} + + \alpha: Regularization parameter + + \nabla: Gradient operator + + g: Noisy Data with Gaussian Noise + + Method = 0 ( PDHG - split ) : K = [ \nabla, + Identity] + + + Method = 1 (PDHG - explicit ): K = \nabla + +""" + +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 ccpi.framework import TestData +import os, sys +loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) + +# Load Data +N = 200 +M = 300 + + +# user can change the size of the input data +# you can choose between +# TestData.PEPPERS 2D + Channel +# TestData.BOAT 2D +# TestData.CAMERA 2D +# TestData.RESOLUTION_CHART 2D +# TestData.SIMPLE_PHANTOM_2D 2D +data = loader.load(TestData.BOAT, size=(N,M), scale=(0,1)) + +ig = data.geometry +ag = ig + +# Create Noisy data. Add Gaussian noise +np.random.seed(10) +noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.1, size=data.shape) ) + +print ("min {} max {}".format(data.as_array().min(), data.as_array().max())) + +# Show Ground Truth and Noisy Data +plt.figure() +plt.subplot(1,3,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(1,3,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.subplot(1,3,3) +plt.imshow((data - noisy_data).as_array()) +plt.title('diff') +plt.colorbar() + +plt.show() + +# Regularisation Parameter +alpha = .1 + +method = '0' + +if method == '0': + + # Create operators + op1 = Gradient(ig, correlation=Gradient.CORRELATION_SPACE) + 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) +pdhg.max_iteration = 10000 +pdhg.update_objective_interval = 100 +pdhg.run(1000, verbose=True) + +# Show Results +plt.figure() +plt.subplot(1,3,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.clim(0,1) +plt.subplot(1,3,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.clim(0,1) +plt.subplot(1,3,3) +plt.imshow(pdhg.get_output().as_array()) +plt.title('TV Reconstruction') +plt.clim(0,1) +plt.colorbar() +plt.show() + +plt.plot(np.linspace(0,N,M), noisy_data.as_array()[int(N/2),:], label = 'Noisy data') +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 = '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 = MOSEK) + + 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,M), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') + plt.plot(np.linspace(0,N,M), u.value[int(N/2),:], label = 'CVX') + plt.plot(np.linspace(0,N,M), data.as_array()[int(N/2),:], label = 'Truth') + + 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/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Gaussian_3D.py b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Gaussian_3D.py new file mode 100644 index 0000000..03dc2ef --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Gaussian_3D.py @@ -0,0 +1,181 @@ +#======================================================================== +# 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 (3D) Denoising using PDHG algorithm: + + +Problem: min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2} + + \alpha: Regularization parameter + + \nabla: Gradient operator + + g: Noisy Data with Gaussian Noise + + Method = 0 ( PDHG - split ) : K = [ \nabla, + Identity] + + + Method = 1 (PDHG - explicit ): K = \nabla + +""" + +from ccpi.framework import ImageData, ImageGeometry + +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, verbose = True) + + +#%% +fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8)) +fig.suptitle('TV Reconstruction',fontsize=20) + + +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() + diff --git a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Poisson.py new file mode 100644 index 0000000..1d887c1 --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Poisson.py @@ -0,0 +1,212 @@ +#======================================================================== +# 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,1} + \int x - g * log(x) + + \alpha: Regularization parameter + + \nabla: Gradient operator + + g: Noisy Data with Poisson Noise + + + Method = 0 ( PDHG - split ) : K = [ \nabla, + Identity] + + + Method = 1 (PDHG - explicit ): K = \nabla + +""" + +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, IndicatorBox, \ + MixedL21Norm, BlockFunction + +from skimage.util import random_noise + +# Create phantom for TV Poisson 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 Poisson noise +n1 = random_noise(data.as_array(), mode = 'poisson', seed = 10) +noisy_data = ImageData(n1) + + +# 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 = 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 = IndicatorBox(lower=0) + +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) + +# Setup and Run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) +pdhg.max_iteration = 3000 +pdhg.update_objective_interval = 200 +pdhg.run(3000, verbose=True) + + +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 + u1 = Variable(ig.shape) + q = Variable() + w = 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(kl_div(noisy_data.as_array(), 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/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_SaltPepper.py b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_SaltPepper.py new file mode 100644 index 0000000..c5709c3 --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_SaltPepper.py @@ -0,0 +1,213 @@ +#======================================================================== +# 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,1} + ||x-g||_{1} + + \alpha: Regularization parameter + + \nabla: Gradient operator + + g: Noisy Data with Salt & Pepper Noise + + + Method = 0 ( PDHG - split ) : K = [ \nabla, + Identity] + + + Method = 1 (PDHG - explicit ): K = \nabla + + +""" + +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 = 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 = 's&p', salt_vs_pepper = 0.9, amount=0.2) +noisy_data = ImageData(n1) + +# Show Ground Truth and Noisy Data +plt.figure(figsize=(15,15)) +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 = 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/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py new file mode 100644 index 0000000..bc0081f --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py @@ -0,0 +1,189 @@ +#======================================================================== +# 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 Generalised Variation (TGV) Tomography 2D using PDHG algorithm: + + +Problem: min_{x>0} \alpha * ||\nabla x - w||_{2,1} + + \beta * || E w ||_{2,1} + + int A x - g log(Ax + \eta) + + \alpha: Regularization parameter + \beta: Regularization parameter + + \nabla: Gradient operator + E: Symmetrized Gradient operator + A: Projection Matrix + + g: Noisy Data with Poisson Noise + + K = [ \nabla, - Identity + ZeroOperator, E + A, ZeroOperator] + +""" + +from ccpi.framework import AcquisitionGeometry, AcquisitionData, 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, Gradient, Identity, \ + SymmetrizedGradient, ZeroOperator +from ccpi.optimisation.functions import IndicatorBox, KullbackLeibler, ZeroFunction,\ + MixedL21Norm, BlockFunction + +from ccpi.astra.ops import AstraProjectorSimple +from ccpi.framework import TestData +import os, sys +from skimage.util import random_noise + +loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) + +# Load Data +#N = 50 +#M = 50 +#data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1)) +# +#ig = data.geometry +#ag = ig +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()) + +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +ag = ig + +# Create noisy data. Add Gaussian noise +scale = 0.5 +eta = 0 +n1 = scale * np.random.poisson(eta + sin.as_array()/scale) + +noisy_data = AcquisitionData(n1, ag) + +#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 ') +device = '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 Parameters +alpha = 1 +beta = 5 + +# 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 = BlockFunction(IndicatorBox(lower=0), 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 = 3000 +pdhg.update_objective_interval = 500 +pdhg.run(3000) + +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/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_gaussian.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_gaussian.py new file mode 100644 index 0000000..bd1bb69 --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_gaussian.py @@ -0,0 +1,212 @@ +#======================================================================== +# 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 2D Tomography Reconstruction using PDHG algorithm: + + +Problem: min_x, x>0 \alpha * ||\nabla x||_{2,1} + \frac{1}{2}||Ax - g||^{2} + + \nabla: Gradient operator + + A: Projection Matrix + g: Noisy sinogram corrupted with Gaussian 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, \ + MixedL21Norm, BlockFunction + +from ccpi.astra.ops import AstraProjectorSimple + + + +# Create phantom for TV 2D tomography +N = 100 +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) + +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, dev) +sin = Aop.direct(data) + +# Create noisy data. Apply Poisson noise +n1 = np.random.normal(0, 3, size=ag.shape) +noisy_data = AcquisitionData(n1 + sin.as_array(), 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 = 50 + +# Create operators +op1 = Gradient(ig) +op2 = Aop + +# 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() + +# Compute operator Norm +normK = operator.norm() + +# Primal & dual stepsizes +sigma = 10 +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 = 200 +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) + + 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) + + tmp = noisy_data.as_array().ravel() + + fidelity = 0.5 * sum_squares(ProjMat * u - tmp) + + solver = MOSEK + obj = Minimize( regulariser + fidelity) + prob = Problem(obj) + result = prob.solve(verbose = True, solver = solver) + + diff_cvx = numpy.abs( pdhg.get_output().as_array() - np.reshape(u.value, (N,N) )) + + 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), np.reshape(u.value, (N,N) )[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 diff --git a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_poisson.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_poisson.py new file mode 100644 index 0000000..95fe2c1 --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_poisson.py @@ -0,0 +1,228 @@ +#======================================================================== +# 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,1} + 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 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 KullbackLeibler, \ + MixedL21Norm, BlockFunction, IndicatorBox + +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 = 50 +M = 50 +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 = 2 + +# 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 = IndicatorBox(lower=0) + + +# 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 = 3000 +pdhg.update_objective_interval = 500 +pdhg.run(3000, verbose = True) + +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) + + tmp = noisy_data.as_array().ravel() + + fidelity = sum(kl_div(tmp, ProjMat * u)) + + constraints = [q>=fidelity, u>=0] + solver = SCS + obj = Minimize( regulariser + q) + prob = Problem(obj, constraints) + result = prob.solve(verbose = True, solver = solver) + + diff_cvx = np.abs(pdhg.get_output().as_array() - np.reshape(u.value, (N, N))) + + 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), np.reshape(u.value, (N, N))[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 diff --git a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_Tikhonov_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_Tikhonov_Tomo2D.py new file mode 100644 index 0000000..02cd053 --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/Tomo/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() + + |