diff options
7 files changed, 0 insertions, 1265 deletions
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 deleted file mode 100644 index 9d00ee1..0000000 --- a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Gaussian.py +++ /dev/null @@ -1,225 +0,0 @@ -#======================================================================== -# 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_Poisson.py b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Poisson.py deleted file mode 100644 index 1d887c1..0000000 --- a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Poisson.py +++ /dev/null @@ -1,212 +0,0 @@ -#======================================================================== -# 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 deleted file mode 100644 index c5709c3..0000000 --- a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_SaltPepper.py +++ /dev/null @@ -1,213 +0,0 @@ -#======================================================================== -# 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_TV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D.py deleted file mode 100644 index f179e95..0000000 --- a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D.py +++ /dev/null @@ -1,173 +0,0 @@ -#======================================================================== -# 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_u \alpha * ||\nabla u||_{2,1} + \frac{1}{2}||Au - g||^{2} - min_u, u>0 \alpha * ||\nabla u||_{2,1} + \int A u - g log (Au + \eta) - - \nabla: Gradient operator - A: System Matrix - g: Noisy sinogram - \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, \ - MixedL21Norm, BlockFunction, KullbackLeibler, IndicatorBox - -from ccpi.astra.ops import AstraProjectorSimple -from ccpi.framework import TestData -from PIL import Image -import os, sys -if int(numpy.version.version.split('.')[1]) > 12: - from skimage.util import random_noise -else: - from demoutil import random_noise - -import scipy.io - -# user supplied input -if len(sys.argv) > 1: - which_noise = int(sys.argv[1]) -else: - which_noise = 1 - -# Load 256 shepp-logan -data256 = scipy.io.loadmat('phantom.mat')['phantom256'] -data = ImageData(numpy.array(Image.fromarray(data256).resize((256,256)))) -N, M = data.shape -ig = ImageGeometry(voxel_num_x=N, voxel_num_y=M) - -# Add it to testdata or use tomophantom -#loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) -#data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(50, 50)) -#ig = data.geometry - -# Create acquisition data and geometry -detectors = N -angles = np.linspace(0, np.pi, 180) -ag = AcquisitionGeometry('parallel','2D',angles, detectors) - -# Select device -device = '0' -#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 Gaussian noise -noises = ['gaussian', 'poisson'] -noise = noises[which_noise] - -if noise == 'poisson': - scale = 5 - eta = 0 - noisy_data = AcquisitionData(np.random.poisson( scale * (eta + sin.as_array()))/scale, ag) -elif noise == 'gaussian': - n1 = np.random.normal(0, 1, size = ag.shape) - noisy_data = AcquisitionData(n1 + sin.as_array(), ag) -else: - raise ValueError('Unsupported Noise ', noise) - -# Show Ground Truth and Noisy Data -plt.figure(figsize=(10,10)) -plt.subplot(1,2,2) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(1,2,1) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.show() - -# Create operators -op1 = Gradient(ig) -op2 = Aop - -# Create BlockOperator -operator = BlockOperator(op1, op2, shape=(2,1) ) - -# Compute operator Norm -normK = operator.norm() - -# Create functions -if noise == 'poisson': - alpha = 3 - f2 = KullbackLeibler(noisy_data) - g = IndicatorBox(lower=0) - sigma = 1 - tau = 1/(sigma*normK**2) - -elif noise == 'gaussian': - alpha = 20 - f2 = 0.5 * L2NormSquared(b=noisy_data) - g = ZeroFunction() - sigma = 10 - tau = 1/(sigma*normK**2) - -f1 = alpha * MixedL21Norm() -f = BlockFunction(f1, f2) - -# 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,ig.shape[1],ig.shape[1]), data.as_array()[int(N/2),:], label = 'GTruth') -plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(N/2),:], label = 'TV 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 deleted file mode 100644 index bd1bb69..0000000 --- a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_gaussian.py +++ /dev/null @@ -1,212 +0,0 @@ -#======================================================================== -# 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 deleted file mode 100644 index e7e33cb..0000000 --- a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_poisson.py +++ /dev/null @@ -1,230 +0,0 @@ -#======================================================================== -# 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 = np.random.poisson(eta + sin.as_array()) - -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/phantom.mat b/Wrappers/Python/demos/PDHG_examples/Tomo/phantom.mat Binary files differdeleted file mode 100755 index c465bbe..0000000 --- a/Wrappers/Python/demos/PDHG_examples/Tomo/phantom.mat +++ /dev/null |