diff options
Diffstat (limited to 'Wrappers')
| -rw-r--r-- | Wrappers/Python/demos/PDHG_TV_Denoising_Gaussian_3D.py | 155 | ||||
| -rw-r--r-- | Wrappers/Python/demos/PDHG_TV_Tomo2D.py | 245 | ||||
| -rw-r--r-- | Wrappers/Python/demos/PDHG_TV_Tomo2D_gaussian.py | 204 |
3 files changed, 0 insertions, 604 deletions
diff --git a/Wrappers/Python/demos/PDHG_TV_Denoising_Gaussian_3D.py b/Wrappers/Python/demos/PDHG_TV_Denoising_Gaussian_3D.py deleted file mode 100644 index c86ddc9..0000000 --- a/Wrappers/Python/demos/PDHG_TV_Denoising_Gaussian_3D.py +++ /dev/null @@ -1,155 +0,0 @@ -# -*- coding: utf-8 -*- -# This work is part of the Core Imaging Library developed by -# Visual Analytics and Imaging System Group of the Science Technology -# Facilities Council, STFC - -# Copyright 2018-2019 Evangelos Papoutsellis and Edoardo Pasca - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -from ccpi.framework import ImageData, ImageGeometry - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - 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) - -fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8)) - -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_TV_Tomo2D.py b/Wrappers/Python/demos/PDHG_TV_Tomo2D.py deleted file mode 100644 index 87d5328..0000000 --- a/Wrappers/Python/demos/PDHG_TV_Tomo2D.py +++ /dev/null @@ -1,245 +0,0 @@ -# -*- 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, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \ - MixedL21Norm, BlockFunction - -from ccpi.astra.ops import AstraProjectorSimple - -""" - -Total Variation Denoising using PDHG algorithm: - - min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) - - -Problem: min_x, x>0 \alpha * ||\nabla x||_{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 - -""" - -# Create phantom for TV 2D tomography -N = 75 -x = np.zeros((N,N)) -x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 -x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 - -data = ImageData(x) -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) - -detectors = N -angles = np.linspace(0, np.pi, N, dtype=np.float32) - -ag = AcquisitionGeometry('parallel','2D',angles, detectors) -Aop = AstraProjectorSimple(ig, ag, 'cpu') -sin = Aop.direct(data) - -# Create noisy data. Apply Poisson noise -scale = 2 -n1 = scale * np.random.poisson(sin.as_array()/scale) -noisy_data = AcquisitionData(n1, ag) - -# Regularisation Parameter -alpha = 5 - -# 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 = ZeroFunction() - -diag_precon = True - -if diag_precon: - - def tau_sigma_precond(operator): - - tau = 1/operator.sum_abs_row() - sigma = 1/ operator.sum_abs_col() - - return tau, sigma - - tau, sigma = tau_sigma_precond(operator) - -else: - # 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, 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 -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) - - fidelity = sum( ProjMat * u - noisy_data.as_array().ravel() * log(ProjMat * u)) - #constraints = [q>= fidelity, u>=0] - constraints = [u>=0] - - solver = SCS - obj = Minimize( regulariser + fidelity) - prob = Problem(obj, constraints) - result = prob.solve(verbose = True, solver = solver) - - -##%% 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) - - - 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), 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]))
\ No newline at end of file diff --git a/Wrappers/Python/demos/PDHG_TV_Tomo2D_gaussian.py b/Wrappers/Python/demos/PDHG_TV_Tomo2D_gaussian.py deleted file mode 100644 index dc473a8..0000000 --- a/Wrappers/Python/demos/PDHG_TV_Tomo2D_gaussian.py +++ /dev/null @@ -1,204 +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} + \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 = 200 -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) -Aop = AstraProjectorSimple(ig, ag, 'cpu') -sin = Aop.direct(data) - -# Create noisy data. Apply Poisson noise -n1 = np.random.normal(0, 3, size=ig.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 |