diff options
23 files changed, 0 insertions, 4526 deletions
diff --git a/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py b/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py deleted file mode 100644 index 2ac050f..0000000 --- a/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py +++ /dev/null @@ -1,211 +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. -# -#========================================================================= - -""" -Conjugate Gradient for (Regularized) Least Squares for Tomography - - -Problem: min_u alpha * || L x ||^{2}_{2} + || A u - g ||_{2}^{2} - - A: Projection operator - g: Sinogram - L: Identity or Gradient Operator - -""" - - -from ccpi.framework import ImageGeometry, ImageData, \ - AcquisitionGeometry, BlockDataContainer, AcquisitionData - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import CGLS -from ccpi.optimisation.operators import BlockOperator, Gradient, Identity - -import tomophantom -from tomophantom import TomoP2D -from ccpi.astra.operators import AstraProjectorSimple -import os - - -# Load Shepp-Logan phantom -model = 1 # select a model number from the library -N = 64 # set dimension of the phantom -path = os.path.dirname(tomophantom.__file__) -path_library2D = os.path.join(path, "Phantom2DLibrary.dat") -phantom_2D = TomoP2D.Model(model, N, path_library2D) - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -data = ImageData(phantom_2D) - -detectors = N -angles = np.linspace(0, np.pi, 180, dtype=np.float32) - -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) - -np.random.seed(10) -noisy_data = AcquisitionData( sin.as_array() + np.random.normal(0,1,ag.shape)) - -# 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() - -# Setup and run the simple CGLS algorithm -x_init = ig.allocate() - -cgls1 = CGLS(x_init = x_init, operator = Aop, data = noisy_data) -cgls1.max_iteration = 20 -cgls1.update_objective_interval = 5 -cgls1.run(20, verbose = True) - -# Setup and run the regularised CGLS algorithm (Tikhonov with Identity) - -x_init = ig.allocate() -alpha1 = 50 -op1 = Identity(ig) - -block_op1 = BlockOperator( Aop, alpha1 * op1, shape=(2,1)) -block_data1 = BlockDataContainer(noisy_data, op1.range_geometry().allocate()) - -cgls2 = CGLS(x_init = x_init, operator = block_op1, data = block_data1) -cgls2.max_iteration = 200 -cgls2.update_objective_interval = 10 -cgls2.run(200, verbose=True) - -# Setup and run the regularised CGLS algorithm (Tikhonov with Gradient) - -x_init = ig.allocate() -alpha2 = 25 -op2 = Gradient(ig) - -block_op2 = BlockOperator( Aop, alpha2 * op2, shape=(2,1)) -block_data2 = BlockDataContainer(noisy_data, op2.range_geometry().allocate()) - -cgls3 = CGLS(x_init=x_init, operator = block_op2, data = block_data2) -cgls3.max_iteration = 200 -cgls3.update_objective_interval = 5 -cgls3.run(200, verbose = True) - -# Show results -plt.figure(figsize=(8,8)) - -plt.subplot(2,2,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') - -plt.subplot(2,2,2) -plt.imshow(cgls1.get_output().as_array()) -plt.title('GGLS') - -plt.subplot(2,2,3) -plt.imshow(cgls2.get_output().as_array()) -plt.title('Regularised GGLS L = {} * Identity'.format(alpha1)) - -plt.subplot(2,2,4) -plt.imshow(cgls3.get_output().as_array()) -plt.title('Regularised GGLS L = {} * Gradient'.format(alpha2)) - -plt.show() - - - -print('Compare CVX vs Regularised CG with L = Gradient') - -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 = alpha2**2 * sum_squares(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('line', proj_geom, vol_geom) - - matrix_id = astra.projector.matrix(proj_id) - - ProjMat = astra.matrix.get(matrix_id) - - fidelity = sum_squares( ProjMat * u - noisy_data.as_array().ravel()) - - solver = MOSEK - obj = Minimize( regulariser + fidelity) - prob = Problem(obj) - result = prob.solve(verbose = True, solver = solver) - - -plt.figure(figsize=(10,20)) - -plt.subplot(1,3,1) -plt.imshow(np.reshape(u.value, (N, N))) -plt.colorbar() - -plt.subplot(1,3,2) -plt.imshow(cgls3.get_output().as_array()) -plt.colorbar() - -plt.subplot(1,3,3) -plt.imshow(np.abs(cgls3.get_output().as_array() - np.reshape(u.value, (N, N)) )) -plt.colorbar() - -plt.show() - -print('Primal Objective (CVX) {} '.format(obj.value)) -print('Primal Objective (CGLS) {} '.format(cgls3.objective[-1])) - diff --git a/Wrappers/Python/demos/CGLS_examples/LinearSystem.py b/Wrappers/Python/demos/CGLS_examples/LinearSystem.py deleted file mode 100644 index cc398cb..0000000 --- a/Wrappers/Python/demos/CGLS_examples/LinearSystem.py +++ /dev/null @@ -1,82 +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. -# -#========================================================================= - -import numpy -from ccpi.optimisation.operators import * -from ccpi.optimisation.algorithms import * -from ccpi.optimisation.functions import * -from ccpi.framework import * - -# Problem dimension. -m = 200 -n = 200 - -numpy.random.seed(10) - -# Create matrix A and data b -Amat = numpy.asarray( numpy.random.randn(m, n), dtype = numpy.float32) -bmat = numpy.asarray( numpy.random.randn(m), dtype=numpy.float32) - -# Geometries for data and solution -vgb = VectorGeometry(m) -vgx = VectorGeometry(n) - -b = vgb.allocate(0, dtype=numpy.float32) -b.fill(bmat) -A = LinearOperatorMatrix(Amat) - -# Setup and run CGLS -x_init = vgx.allocate() - -cgls = CGLS(x_init = x_init, operator = A, data = b) -cgls.max_iteration = 2000 -cgls.update_objective_interval = 200 -cgls.run(2000, verbose = True) - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - -# Compare with CVX -x = Variable(n) -obj = Minimize(sum_squares(A.A*x - b.as_array())) -prob = Problem(obj) - -# choose solver -if 'MOSEK' in installed_solvers(): - solver = MOSEK -else: - solver = SCS - -result = prob.solve(solver = MOSEK) - - -diff_sol = x.value - cgls.get_output().as_array() - -print('Error |CVX - CGLS| = {}'.format(numpy.sum(numpy.abs(diff_sol)))) -print('CVX objective = {}'.format(obj.value)) -print('CGLS objective = {}'.format(cgls.objective[-1])) - - - - diff --git a/Wrappers/Python/demos/CGLS_examples/Regularised_CGLS_Denoising.py b/Wrappers/Python/demos/CGLS_examples/Regularised_CGLS_Denoising.py deleted file mode 100644 index c124fa1..0000000 --- a/Wrappers/Python/demos/CGLS_examples/Regularised_CGLS_Denoising.py +++ /dev/null @@ -1,146 +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. -# -#========================================================================= - -""" -Conjugate Gradient for (Regularized) Least Squares for Tomography - - -Problem: min_u alpha * || L x ||^{2}_{2} + || A u - g ||_{2}^{2} - - A: Identity operator - g: Sinogram - L: Identity or Gradient Operator - -""" - - -from ccpi.framework import ImageGeometry, ImageData, \ - AcquisitionGeometry, BlockDataContainer, AcquisitionData - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import CGLS -from ccpi.optimisation.operators import BlockOperator, Gradient, Identity -from ccpi.framework import TestData - -import os, sys - -loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) -data = loader.load(TestData.SHAPES) -ig = data.geometry -ag = ig - -noisy_data = ImageData(TestData.random_noise(data.as_array(), mode = 'gaussian', seed = 1)) -#noisy_data = ImageData(data.as_array()) - -# 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() - -# Setup and run the regularised CGLS algorithm (Tikhonov with Gradient) -x_init = ig.allocate() -alpha = 2 -op = Gradient(ig) - -block_op = BlockOperator( Identity(ig), alpha * op, shape=(2,1)) -block_data = BlockDataContainer(noisy_data, op.range_geometry().allocate()) - -cgls = CGLS(x_init=x_init, operator = block_op, data = block_data) -cgls.max_iteration = 200 -cgls.update_objective_interval = 5 -cgls.run(200, verbose = True) - -# Show results -plt.figure(figsize=(20,10)) -plt.subplot(3,1,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.subplot(3,1,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy') -plt.subplot(3,1,3) -plt.imshow(cgls.get_output().as_array()) -plt.title('Regularised GGLS with Gradient') -plt.show() - -#%% - -print('Compare CVX vs Regularised CG with L = Gradient') - -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) - #q = Variable() - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - regulariser = alpha**2 * sum_squares(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - - fidelity = 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) - - -plt.figure(figsize=(20,20)) -plt.subplot(3,1,1) -plt.imshow(np.reshape(u.value, ig.shape)) -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(cgls.get_output().as_array()) -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow(np.abs(cgls.get_output().as_array() - np.reshape(u.value, ig.shape) )) -plt.colorbar() -plt.show() - -print('Primal Objective (CVX) {} '.format(obj.value)) -print('Primal Objective (CGLS) {} '.format(cgls.objective[-1])) - diff --git a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py deleted file mode 100644 index 09e350b..0000000 --- a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py +++ /dev/null @@ -1,199 +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. -# -#========================================================================= - -""" -Compare solutions of FISTA & PDHG & CGLS - & Astra Built-in algorithms for Least Squares - - -Problem: min_x || A x - g ||_{2}^{2} - - A: Projection operator - g: Sinogram - -""" - - -from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, CGLS, FISTA - -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, FunctionOperatorComposition -from ccpi.astra.ops import AstraProjectorSimple -import astra - -import tomophantom -from tomophantom import TomoP2D -import os - -# Load Shepp-Logan phantom -model = 1 # select a model number from the library -N = 128 # set dimension of the phantom -path = os.path.dirname(tomophantom.__file__) -path_library2D = os.path.join(path, "Phantom2DLibrary.dat") -phantom_2D = TomoP2D.Model(model, N, path_library2D) - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -data = ImageData(phantom_2D) - -detectors = N -angles = np.linspace(0, np.pi, 180, dtype=np.float32) - -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) -sinogram = Aop.direct(data) - - - -############################################################################### -# Setup and run Astra CGLS algorithm -vol_geom = astra.create_vol_geom(N, N) -proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) -proj_id = astra.create_projector('line', proj_geom, vol_geom) - -# Create a sinogram id -sinogram_id = astra.data2d.create('-sino', proj_geom, sinogram.as_array()) - -# Create a data id -rec_id = astra.data2d.create('-vol', vol_geom) - -cgls_astra = astra.astra_dict('CGLS') -cgls_astra['ReconstructionDataId'] = rec_id -cgls_astra['ProjectionDataId'] = sinogram_id -cgls_astra['ProjectorId'] = proj_id - -# Create the algorithm object from the configuration structure -alg_id = astra.algorithm.create(cgls_astra) - -astra.algorithm.run(alg_id, 500) -recon_cgls_astra = ImageData(astra.data2d.get(rec_id)) - - -#%% - -# Setup and run the simple CGLS algorithm -x_init = ig.allocate() - -cgls = CGLS(x_init = x_init, operator = Aop, data = sinogram) -cgls.max_iteration = 500 -cgls.update_objective_interval = 100 -cgls.run(500, verbose = True) - -#%% - -############################################################################### -# Setup and run the PDHG algorithm -operator = Aop -f = L2NormSquared(b = sinogram) -g = ZeroFunction() - -## Compute operator Norm -normK = operator.norm() - -## Primal & dual stepsizes -sigma = 0.02 -tau = 1/(sigma*normK**2) - - -pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) -pdhg.max_iteration = 1000 -pdhg.update_objective_interval = 100 -pdhg.run(1000, verbose=True) - -#%% -############################################################################### -# Setup and run the FISTA algorithm -fidelity = FunctionOperatorComposition(L2NormSquared(b=sinogram), Aop) -regularizer = ZeroFunction() - -fista = FISTA(x_init=x_init , f=fidelity, g=regularizer) -fista.max_iteration = 500 -fista.update_objective_interval = 100 -fista.run(500, verbose = True) - -#%% Show results - -plt.figure(figsize=(10,10)) -plt.suptitle('Reconstructions ', fontsize=16) - -plt.subplot(2,2,1) -plt.imshow(cgls.get_output().as_array()) -plt.colorbar() -plt.title('CGLS reconstruction') - -plt.subplot(2,2,2) -plt.imshow(fista.get_output().as_array()) -plt.colorbar() -plt.title('FISTA reconstruction') - -plt.subplot(2,2,3) -plt.imshow(pdhg.get_output().as_array()) -plt.colorbar() -plt.title('PDHG reconstruction') - -plt.subplot(2,2,4) -plt.imshow(recon_cgls_astra.as_array()) -plt.colorbar() -plt.title('CGLS astra') - -diff1 = pdhg.get_output() - recon_cgls_astra -diff2 = fista.get_output() - recon_cgls_astra -diff3 = cgls.get_output() - recon_cgls_astra - -plt.figure(figsize=(15,15)) - -plt.subplot(3,1,1) -plt.imshow(diff1.abs().as_array()) -plt.title('Diff PDHG vs CGLS astra') -plt.colorbar() - -plt.subplot(3,1,2) -plt.imshow(diff2.abs().as_array()) -plt.title('Diff FISTA vs CGLS astra') -plt.colorbar() - -plt.subplot(3,1,3) -plt.imshow(diff3.abs().as_array()) -plt.title('Diff CLGS vs CGLS astra') -plt.colorbar() - - -cgls_astra_obj = fidelity(ImageData(recon_cgls_astra)) - -print('Primal Objective (FISTA) {} '.format(fista.objective[-1])) -print('Primal Objective (CGLS) {} '.format(cgls.objective[-1])) -print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) -print('Primal Objective (CGLS_astra) {} '.format(cgls_astra_obj)) - - - diff --git a/Wrappers/Python/demos/CompareAlgorithms/CGLS_PDHG_Tikhonov.py b/Wrappers/Python/demos/CompareAlgorithms/CGLS_PDHG_Tikhonov.py deleted file mode 100644 index 9b6d10f..0000000 --- a/Wrappers/Python/demos/CompareAlgorithms/CGLS_PDHG_Tikhonov.py +++ /dev/null @@ -1,133 +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. -# -#========================================================================= - -""" -Compare solutions of PDHG & "Block CGLS" algorithms for - - -Problem: min_x alpha * ||\grad x ||^{2}_{2} + || A x - g ||_{2}^{2} - - - A: Projection operator - g: Sinogram - -""" - - -from ccpi.framework import AcquisitionGeometry, BlockDataContainer, AcquisitionData - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, CGLS -from ccpi.optimisation.operators import BlockOperator, Gradient - -from ccpi.optimisation.functions import ZeroFunction, BlockFunction, L2NormSquared -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')) - -# Create Ground truth phantom and Sinogram -N = 150 -M = 150 -data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1)) -ig = data.geometry - -detectors = N -angles = np.linspace(0, np.pi, N, dtype=np.float32) -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) - -noisy_data = AcquisitionData( sin.as_array() + np.random.normal(0,3,ig.shape)) - -# Setup and run the CGLS algorithm -alpha = 50 -Grad = Gradient(ig) - -# Form Tikhonov as a Block CGLS structure -op_CGLS = BlockOperator( Aop, alpha * Grad, shape=(2,1)) -block_data = BlockDataContainer(noisy_data, Grad.range_geometry().allocate()) - -x_init = ig.allocate() -cgls = CGLS(x_init=x_init, operator=op_CGLS, data=block_data) -cgls.max_iteration = 1000 -cgls.update_objective_interval = 200 -cgls.run(1000,verbose=False) - - -#Setup and run the PDHG algorithm - -# Create BlockOperator -op_PDHG = BlockOperator(Grad, Aop, shape=(2,1) ) -# Create functions -f1 = 0.5 * alpha**2 * L2NormSquared() -f2 = 0.5 * L2NormSquared(b = noisy_data) -f = BlockFunction(f1, f2) -g = ZeroFunction() - -## Compute operator Norm -normK = op_PDHG.norm() - -## Primal & dual stepsizes -sigma = 10 -tau = 1/(sigma*normK**2) - -pdhg = PDHG(f=f,g=g,operator=op_PDHG, tau=tau, sigma=sigma) -pdhg.max_iteration = 1000 -pdhg.update_objective_interval = 200 -pdhg.run(1000, verbose=False) - -# Show results -plt.figure(figsize=(10,10)) - -plt.subplot(2,1,1) -plt.imshow(cgls.get_output().as_array()) -plt.title('CGLS reconstruction') - -plt.subplot(2,1,2) -plt.imshow(pdhg.get_output().as_array()) -plt.title('PDHG reconstruction') - -plt.show() - -diff1 = pdhg.get_output() - cgls.get_output() - -plt.imshow(diff1.abs().as_array()) -plt.title('Diff PDHG vs CGLS') -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), cgls.get_output().as_array()[int(N/2),:], label = 'CGLS') -plt.legend() -plt.title('Middle Line Profiles') -plt.show() diff --git a/Wrappers/Python/demos/CompareAlgorithms/FISTA_vs_PDHG.py b/Wrappers/Python/demos/CompareAlgorithms/FISTA_vs_PDHG.py deleted file mode 100644 index c24ebac..0000000 --- a/Wrappers/Python/demos/CompareAlgorithms/FISTA_vs_PDHG.py +++ /dev/null @@ -1,124 +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. -# -#========================================================================= - - -""" -Compare solutions of FISTA & PDHG algorithms - - -Problem: min_x alpha * ||\grad x ||^{2}_{2} + || x - g ||_{1} - - \alpha: Regularization parameter - - \nabla: Gradient operator - - g: Noisy Data with Salt & Pepper Noise - -""" - -from ccpi.framework import ImageData, ImageGeometry - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import FISTA, PDHG - -from ccpi.optimisation.operators import BlockOperator, Gradient, Identity -from ccpi.optimisation.functions import L2NormSquared, L1Norm, \ - FunctionOperatorComposition, BlockFunction, ZeroFunction - -from skimage.util import random_noise - - -# Create Ground truth and Noisy data -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 - -n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2) -noisy_data = ImageData(n1) - -# Regularisation Parameter -alpha = 5 - -############################################################################### -# Setup and run the FISTA algorithm -operator = Gradient(ig) -fidelity = L1Norm(b=noisy_data) -regulariser = FunctionOperatorComposition(alpha * L2NormSquared(), operator) - -x_init = ig.allocate() -opt = {'memopt':True} -fista = FISTA(x_init=x_init , f=regulariser, g=fidelity, opt=opt) -fista.max_iteration = 2000 -fista.update_objective_interval = 50 -fista.run(2000, verbose=False) -############################################################################### - - -############################################################################### -# Setup and run the PDHG algorithm -op1 = Gradient(ig) -op2 = Identity(ig, ag) - -operator = BlockOperator(op1, op2, shape=(2,1) ) -f = BlockFunction(alpha * L2NormSquared(), fidelity) -g = ZeroFunction() - -normK = operator.norm() - -sigma = 1 -tau = 1/(sigma*normK**2) - -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=False) -############################################################################### - -# Show results - -plt.figure(figsize=(10,10)) - -plt.subplot(2,1,1) -plt.imshow(pdhg.get_output().as_array()) -plt.title('PDHG reconstruction') - -plt.subplot(2,1,2) -plt.imshow(fista.get_output().as_array()) -plt.title('FISTA reconstruction') - -plt.show() - -diff1 = pdhg.get_output() - fista.get_output() - -plt.imshow(diff1.abs().as_array()) -plt.title('Diff PDHG vs FISTA') -plt.colorbar() -plt.show() - - diff --git a/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py b/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py deleted file mode 100644 index 1a96a2d..0000000 --- a/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py +++ /dev/null @@ -1,215 +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. -# -#========================================================================= - -from ccpi.framework import AcquisitionGeometry -from ccpi.optimisation.algorithms import FISTA -from ccpi.optimisation.functions import IndicatorBox, ZeroFunction, \ - L2NormSquared, FunctionOperatorComposition -from ccpi.astra.operators import AstraProjectorSimple - -import numpy as np -import matplotlib.pyplot as plt -from ccpi.framework import TestData -import os, sys -from ccpi.optimisation.funcs import Norm2sq - -# Load Data -loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) - -N = 50 -M = 50 -data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1)) - -ig = data.geometry - -# Show Ground Truth -plt.figure(figsize=(5,5)) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.show() - -#%% -# Set up AcquisitionGeometry object to hold the parameters of the measurement -# setup geometry: # Number of angles, the actual angles from 0 to -# pi for parallel beam and 0 to 2pi for fanbeam, set the width of a detector -# pixel relative to an object pixel, the number of detector pixels, and the -# source-origin and origin-detector distance (here the origin-detector distance -# set to 0 to simulate a "virtual detector" with same detector pixel size as -# object pixel size). - -#device = input('Available device: GPU==1 / CPU==0 ') - -#if device=='1': -# dev = 'gpu' -#else: -# dev = 'cpu' - -test_case = 1 -dev = 'cpu' - -if test_case==1: - - detectors = N - angles_num = N - det_w = 1.0 - - angles = np.linspace(0, np.pi, angles_num, endpoint=False) - ag = AcquisitionGeometry('parallel', - '2D', - angles, - detectors,det_w) - -elif test_case==2: - - SourceOrig = 200 - OrigDetec = 0 - angles = np.linspace(0,2*np.pi,angles_num) - ag = AcquisitionGeometry('cone', - '2D', - angles, - detectors, - det_w, - dist_source_center=SourceOrig, - dist_center_detector=OrigDetec) -else: - NotImplemented - -# Set up Operator object combining the ImageGeometry and AcquisitionGeometry -# wrapping calls to ASTRA as well as specifying whether to use CPU or GPU. -Aop = AstraProjectorSimple(ig, ag, dev) - -# Forward and backprojection are available as methods direct and adjoint. Here -# generate test data b and do simple backprojection to obtain z. -sin = Aop.direct(data) -back_proj = Aop.adjoint(sin) - -plt.figure(figsize=(5,5)) -plt.imshow(sin.array) -plt.title('Simulated data') -plt.show() - -plt.figure(figsize=(5,5)) -plt.imshow(back_proj.array) -plt.title('Backprojected data') -plt.show() - -#%% - -# Using the test data b, different reconstruction methods can now be set up as -# demonstrated in the rest of this file. In general all methods need an initial -# guess and some algorithm options to be set: - -f = FunctionOperatorComposition(L2NormSquared(b=sin), Aop) -#f = Norm2sq(Aop, sin, c=0.5, memopt=True) -g = ZeroFunction() - -x_init = ig.allocate() -fista = FISTA(x_init=x_init, f=f, g=g) -fista.max_iteration = 1000 -fista.update_objective_interval = 100 -fista.run(1000, verbose=True) - -plt.figure() -plt.imshow(fista.get_output().as_array()) -plt.title('FISTA unconstrained') -plt.colorbar() -plt.show() - - -# Run FISTA for least squares with lower/upper bound -fista0 = FISTA(x_init=x_init, f=f, g=IndicatorBox(lower=0,upper=1)) -fista0.max_iteration = 1000 -fista0.update_objective_interval = 100 -fista0.run(1000, verbose=True) - -plt.figure() -plt.imshow(fista0.get_output().as_array()) -plt.title('FISTA constrained in [0,1]') -plt.colorbar() -plt.show() - -#plt.figure() -#plt.semilogy(fista0.objective) -#plt.title('FISTA constrained in [0,1]') -#plt.show() - -#%% Check with CVX solution - -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) - - # 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('linear', proj_geom, vol_geom) - - matrix_id = astra.projector.matrix(proj_id) - - ProjMat = astra.matrix.get(matrix_id) - - tmp = sin.as_array().ravel() - - fidelity = 0.5 * sum_squares(ProjMat * u - tmp) - - solver = MOSEK - obj = Minimize(fidelity) - constraints = [u>=0, u<=1] - prob = Problem(obj, constraints=constraints) - result = prob.solve(verbose = True, solver = solver) - - diff_cvx = numpy.abs( fista0.get_output().as_array() - np.reshape(u.value, (N,M) )) - - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(fista0.get_output().as_array()) - plt.title('FISTA solution') - plt.colorbar() - plt.subplot(3,1,2) - plt.imshow(np.reshape(u.value, (N, M))) - 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), fista0.get_output().as_array()[int(N/2),:], label = 'FISTA') - plt.plot(np.linspace(0,N,M), np.reshape(u.value, (N,M) )[int(N/2),:], label = 'CVX') - plt.legend() - plt.title('Middle Line Profiles') - plt.show() - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (FISTA) {} '.format(fista0.loss[1]))
\ No newline at end of file diff --git a/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py b/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py deleted file mode 100644 index 6007990..0000000 --- a/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py +++ /dev/null @@ -1,201 +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. -# -#========================================================================= - -""" - -"Tikhonov regularization" for Poisson denoising using FISTA algorithm: - -Problem: min_x, x>0 \alpha * ||\nabla x||_{2}^{2} + \int x - g * log(x) - - \alpha: Regularization parameter - - \nabla: Gradient operator - - g: Noisy Data with Poisson Noise - - -""" - -from ccpi.framework import ImageData - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import FISTA - -from ccpi.optimisation.operators import Gradient -from ccpi.optimisation.functions import KullbackLeibler, L2NormSquared, FunctionOperatorComposition - -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 = 100 -M = 100 -data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1)) - -ig = data.geometry -ag = ig - -# Create Noisy data with 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 = 10 - -# Setup and run the FISTA algorithm -operator = Gradient(ig) -fid = KullbackLeibler(noisy_data) - -def KL_Prox_PosCone(x, tau, out=None): - - if out is None: - tmp = 0.5 *( (x - fid.bnoise - tau) + ( (x + fid.bnoise - tau)**2 + 4*tau*fid.b ) .sqrt() ) - return tmp.maximum(0) - else: - tmp = 0.5 *( (x - fid.bnoise - tau) + - ( (x + fid.bnoise - tau)**2 + 4*tau*fid.b ) .sqrt() - ) - x.add(fid.bnoise, out=out) - out -= tau - out *= out - tmp = fid.b * (4 * tau) - out.add(tmp, out=out) - out.sqrt(out=out) - - x.subtract(fid.bnoise, out=tmp) - tmp -= tau - - out += tmp - - out *= 0.5 - - # ADD the constraint here - out.maximum(0, out=out) - -fid.proximal = KL_Prox_PosCone - -reg = FunctionOperatorComposition(alpha * L2NormSquared(), operator) - -x_init = ig.allocate() -fista = FISTA(x_init=x_init , f=reg, g=fid) -fista.max_iteration = 2000 -fista.update_objective_interval = 500 -fista.run(2000, verbose=True) - -# Show results -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(fista.get_output().as_array()) -plt.title('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), fista.get_output().as_array()[int(N/2),:], label = '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() - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - regulariser = alpha * sum_squares(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 = SCS - obj = Minimize( regulariser + q) - prob = Problem(obj, constraints) - result = prob.solve(verbose = True, solver = solver) - - diff_cvx = numpy.abs( fista.get_output().as_array() - u1.value ) - - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(fista.get_output().as_array()) - plt.title('FISTA 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,M), fista.get_output().as_array()[int(N/2),:], label = 'FISTA') - plt.plot(np.linspace(0,N,M), 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 (FISTA) {} '.format(fista.loss[1])) - - - diff --git a/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py b/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py deleted file mode 100644 index e69060f..0000000 --- a/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py +++ /dev/null @@ -1,273 +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. -# -#========================================================================= - - -from ccpi.framework import ImageGeometry, ImageData, AcquisitionGeometry, AcquisitionData, BlockDataContainer - -import numpy as numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, CGLS -from ccpi.optimisation.algs import CGLS as CGLS_old - -from ccpi.optimisation.operators import BlockOperator, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction - -from ccpi.astra.operators import AstraProjectorMC -from scipy.io import loadmat -import h5py - -#%% - -phantom = 'powder' - -if phantom == 'carbon': - pathname = '/media/newhd/shared/Data/ColourBay/spectral_data_sets/CarbonPd/' - filename = 'carbonPd_full_sinogram_stripes_removed.mat' - X = loadmat(pathname + filename) - X = numpy.transpose(X['SS'],(3,1,2,0)) - X = X[80:100] # delete this to take all channels - -elif phantom == 'powder': - pathname = '/media/newhd/shared/DataProcessed/' - filename = 'S_180.mat' - path = pathname + filename - arrays = {} - f = h5py.File(path) - for k, v in f.items(): - arrays[k] = numpy.array(v) - XX = arrays['S'] - X = numpy.transpose(XX,(0,2,1,3)) - X = X[100:120] - - -#%% Setup Geometry of Colorbay - -num_channels = X.shape[0] -num_pixels_h = X.shape[3] -num_pixels_v = X.shape[2] -num_angles = X.shape[1] - -# Display a single projection in a single channel -plt.imshow(X[5,5,:,:]) -plt.title('Example of a projection image in one channel' ) -plt.show() - -# Set angles to use -angles = numpy.linspace(-numpy.pi,numpy.pi,num_angles,endpoint=False) - -# Define full 3D acquisition geometry and data container. -# Geometric info is taken from the txt-file in the same dir as the mat-file -ag = AcquisitionGeometry('cone', - '3D', - angles, - pixel_num_h=num_pixels_h, - pixel_size_h=0.25, - pixel_num_v=num_pixels_v, - pixel_size_v=0.25, - dist_source_center=233.0, - dist_center_detector=245.0, - channels=num_channels) -data = AcquisitionData(X, geometry=ag) - -# Reduce to central slice by extracting relevant parameters from data and its -# geometry. Perhaps create function to extract central slice automatically? -data2d = data.subset(vertical=40) -ag2d = AcquisitionGeometry('cone', - '2D', - ag.angles, - pixel_num_h=ag.pixel_num_h, - pixel_size_h=ag.pixel_size_h, - pixel_num_v=1, - pixel_size_v=ag.pixel_size_h, - dist_source_center=ag.dist_source_center, - dist_center_detector=ag.dist_center_detector, - channels=ag.channels) -data2d.geometry = ag2d - -# Set up 2D Image Geometry. -# First need the geometric magnification to scale the voxel size relative -# to the detector pixel size. -mag = (ag.dist_source_center + ag.dist_center_detector)/ag.dist_source_center -ig2d = ImageGeometry(voxel_num_x=ag2d.pixel_num_h, - voxel_num_y=ag2d.pixel_num_h, - voxel_size_x=ag2d.pixel_size_h/mag, - voxel_size_y=ag2d.pixel_size_h/mag, - channels=X.shape[0]) - -# Create GPU multichannel projector/backprojector operator with ASTRA. -Aall = AstraProjectorMC(ig2d,ag2d,'gpu') - -# Compute and simple backprojction and display one channel as image. -Xbp = Aall.adjoint(data2d) -plt.imshow(Xbp.subset(channel=5).array) -plt.show() - -#%% CGLS - -def callback(iteration, objective, x): - plt.imshow(x.as_array()[5]) - plt.colorbar() - plt.show() - -x_init = ig2d.allocate() -cgls1 = CGLS(x_init=x_init, operator=Aall, data=data2d) -cgls1.max_iteration = 100 -cgls1.update_objective_interval = 1 -cgls1.run(5,verbose=True, callback = callback) - -plt.imshow(cgls1.get_output().subset(channel=5).array) -plt.title('CGLS') -plt.show() - -#%% Tikhonov - -alpha = 2.5 -Grad = Gradient(ig2d, correlation=Gradient.CORRELATION_SPACE) # use also CORRELATION_SPACECHANNEL - -# Form Tikhonov as a Block CGLS structure -op_CGLS = BlockOperator( Aall, alpha * Grad, shape=(2,1)) -block_data = BlockDataContainer(data2d, Grad.range_geometry().allocate()) - -cgls2 = CGLS(x_init=x_init, operator=op_CGLS, data=block_data) -cgls2.max_iteration = 100 -cgls2.update_objective_interval = 1 - -cgls2.run(10,verbose=True, callback=callback) - -plt.imshow(cgls2.get_output().subset(channel=5).array) -plt.title('Tikhonov') -plt.show() - -#%% Total Variation - -# Regularisation Parameter -#alpha_TV = 0.08 # for carbon -alpha_TV = 0.08 # for powder - -# Create operators -op1 = Gradient(ig2d, correlation=Gradient.CORRELATION_SPACE) -op2 = Aall - -# Create BlockOperator -operator = BlockOperator(op1, op2, shape=(2,1) ) - -# Create functions -f1 = alpha_TV * MixedL21Norm() -f2 = 0.5 * L2NormSquared(b=data2d) -f = BlockFunction(f1, f2) -g = ZeroFunction() - -# Compute operator Norm -#normK = 8.70320267279591 # For powder Run one time no need to compute again takes time -normK = 14.60320657253632 # for carbon - -# 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 = 100 -pdhg.run(1000, verbose =True, callback=callback) - - -#%% Show sinograms -channel_ind = [10,15,19] - -plt.figure(figsize=(15,15)) - -plt.subplot(4,3,1) -plt.imshow(data2d.subset(channel = channel_ind[0]).as_array()) -plt.title('Channel {}'.format(channel_ind[0])) -plt.colorbar() - -plt.subplot(4,3,2) -plt.imshow(data2d.subset(channel = channel_ind[1]).as_array()) -plt.title('Channel {}'.format(channel_ind[1])) -plt.colorbar() - -plt.subplot(4,3,3) -plt.imshow(data2d.subset(channel = channel_ind[2]).as_array()) -plt.title('Channel {}'.format(channel_ind[2])) -plt.colorbar() - -############################################################################### -# Show CGLS -plt.subplot(4,3,4) -plt.imshow(cgls1.get_output().subset(channel = channel_ind[0]).as_array()) -plt.colorbar() - -plt.subplot(4,3,5) -plt.imshow(cgls1.get_output().subset(channel = channel_ind[1]).as_array()) -plt.colorbar() - -plt.subplot(4,3,6) -plt.imshow(cgls1.get_output().subset(channel = channel_ind[2]).as_array()) -plt.colorbar() - -############################################################################### -# Show Tikhonov - -plt.subplot(4,3,7) -plt.imshow(cgls2.get_output().subset(channel = channel_ind[0]).as_array()) -plt.colorbar() - -plt.subplot(4,3,8) -plt.imshow(cgls2.get_output().subset(channel = channel_ind[1]).as_array()) -plt.colorbar() - -plt.subplot(4,3,9) -plt.imshow(cgls2.get_output().subset(channel = channel_ind[2]).as_array()) -plt.colorbar() - - -############################################################################### -# Show Total variation - -plt.subplot(4,3,10) -plt.imshow(pdhg.get_output().subset(channel = channel_ind[0]).as_array()) -plt.colorbar() - -plt.subplot(4,3,11) -plt.imshow(pdhg.get_output().subset(channel = channel_ind[1]).as_array()) -plt.colorbar() - -plt.subplot(4,3,12) -plt.imshow(pdhg.get_output().subset(channel = channel_ind[2]).as_array()) -plt.colorbar() - - -############################################################################### - - - - - - - - - - - - diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py deleted file mode 100755 index 9dbcf3e..0000000 --- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py +++ /dev/null @@ -1,282 +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 Generalised Variation (TGV) Denoising using PDHG algorithm: - - -Problem: min_{u} \alpha * ||\nabla u - w||_{2,1} + - \beta * || E u ||_{2,1} + - Fidelity(u, g) - - \nabla: Gradient operator - E: Symmetrized Gradient operator - \alpha: Regularization parameter - \beta: Regularization parameter - - g: Noisy Data - - Fidelity = 1) L2NormSquarred ( \frac{1}{2} * || u - g ||_{2}^{2} ) if Noise is Gaussian - 2) L1Norm ( ||u - g||_{1} )if Noise is Salt & Pepper - 3) Kullback Leibler (\int u - g * log(u) + Id_{u>0}) if Noise is Poisson - - - Method = 0 ( PDHG - split ) : K = [ \nabla, - Identity - ZeroOperator, E - Identity, ZeroOperator] - - - Method = 1 (PDHG - explicit ): K = [ \nabla, - Identity - ZeroOperator, E ] - - Default: TGV denoising - noise = Gaussian - Fidelity = L2NormSquarred - method = 0 - -""" - -from ccpi.framework import ImageData - -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 - -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) - - -loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) -data = loader.load(TestData.SHAPES) -ig = data.geometry -ag = ig - -# Create noisy data. -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': - scale = 5 - n1 = random_noise( data.as_array()/scale, mode = noise, seed = 10)*scale -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 Parameter depending on the noise distribution -if noise == 's&p': - alpha = 0.8 -elif noise == 'poisson': - alpha = .3 -elif noise == 'gaussian': - alpha = .2 - -# TODO add ref why this choice -beta = 2 * alpha - -# Fidelity -if noise == 's&p': - f3 = L1Norm(b=noisy_data) -elif noise == 'poisson': - f3 = KullbackLeibler(noisy_data) -elif noise == 'gaussian': - f3 = 0.5 * 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() - - 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 = 100 -pdhg.run(2000) - -# Show results -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()[0].as_array()) -plt.title('TGV Reconstruction') -plt.colorbar() -plt.subplot(1,4,4) -plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[int(ig.shape[0]/2),:], label = 'GTruth') -plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output()[0].as_array()[int(ig.shape[0]/2),:], label = 'TGV 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(ig.shape) - w2 = Variable(ig.shape) - - # 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 = [] - - # choose solver - if 'MOSEK' in installed_solvers(): - solver = MOSEK - else: - solver = SCS - - # fidelity - if noise == 's&p': - fidelity = pnorm( u - noisy_data.as_array(),1) - elif noise == 'poisson': - fidelity = sum(kl_div(noisy_data.as_array(), u)) - solver = SCS - elif noise == 'gaussian': - fidelity = 0.5 * sum_squares(noisy_data.as_array() - u) - - 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,ig.shape[1],ig.shape[1]), pdhg.get_output()[0].as_array()[int(ig.shape[0]/2),:], label = 'PDHG') - plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), u.value[int(ig.shape[0]/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/GatherAll/PDHG_TV_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py deleted file mode 100755 index 6937fa0..0000000 --- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py +++ /dev/null @@ -1,280 +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_{u}, \alpha * ||\nabla u||_{2,1} + Fidelity(u, g) - - \alpha: Regularization parameter - - \nabla: Gradient operator - - g: Noisy Data - - Fidelity = 1) L2NormSquarred ( \frac{1}{2} * || u - g ||_{2}^{2} ) if Noise is Gaussian - 2) L1Norm ( ||u - g||_{1} )if Noise is Salt & Pepper - 3) Kullback Leibler (\int u - g * log(u) + Id_{u>0}) if Noise is Poisson - - Method = 0 ( PDHG - split ) : K = [ \nabla, - Identity] - - - Method = 1 (PDHG - explicit ): K = \nabla - - - Default: ROF denoising - noise = Gaussian - Fidelity = L2NormSquarred - method = 0 - - -""" - -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 -#import scipy.io -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) - - -loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) -data = loader.load(TestData.SHAPES, size=(50,50)) -ig = data.geometry -ag = ig - -# Create noisy data. -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': - scale = 5 - n1 = random_noise( data.as_array()/scale, mode = noise, seed = 10)*scale -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) - -# 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 depending on the noise distribution -if noise == 's&p': - alpha = 0.8 -elif noise == 'poisson': - alpha = 1 -elif noise == 'gaussian': - alpha = .3 - -# 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) - -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 - f = BlockFunction(alpha * MixedL21Norm(), f2) - g = ZeroFunction() - -else: - - operator = Gradient(ig) - f = alpha * MixedL21Norm() - g = f2 - - -# 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 = 100 -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,ig.shape[1],ig.shape[1]), data.subset(channel=row).as_array()[int(N/2),:], label = 'GTruth') - plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), 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,ig.shape[1],ig.shape[1]), data.as_array()[int(ig.shape[0]/2),:], label = 'GTruth') - plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(ig.shape[0]/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([Constant(DX.matrix()) * vec(u), Constant(DY.matrix()) * vec(u)]), 2, axis = 0)) - - # choose solver - if 'MOSEK' in installed_solvers(): - solver = MOSEK - else: - solver = SCS - - # fidelity - if noise == 's&p': - fidelity = pnorm( u - noisy_data.as_array(),1) - elif noise == 'poisson': - fidelity = sum(kl_div(noisy_data.as_array(), u)) - solver = SCS - elif noise == 'gaussian': - fidelity = 0.5 * sum_squares(noisy_data.as_array() - u) - - 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,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(ig.shape[0]/2),:], label = 'PDHG') - plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), u.value[int(ig.shape[0]/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_TV_Denoising_2D_time.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_2D_time.py deleted file mode 100644 index febe76d..0000000 --- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_2D_time.py +++ /dev/null @@ -1,243 +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 (Dynamic) Denoising using PDHG algorithm and Tomophantom: - - -Problem: min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2} - - \alpha: Regularization parameter - - \nabla: Gradient operator - - g: 2D Dynamic noisy data with Gaussian Noise - - 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, 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 - -path = os.path.dirname(tomophantom.__file__) -path_library2D = os.path.join(path, "Phantom2DLibrary.dat") -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 - -# Setup geometries -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. Apply Gaussian noise -np.random.seed(10) -noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.25, size=ig.shape) ) - -# time-frames index -tindex = [8, 16, 24] - -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, 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 = [8, 16, 24] -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() - -#%% -import matplotlib.animation as animation -fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 30)) -ims1 = [] -ims2 = [] -ims3 = [] -for sl in range(0,np.shape(phantom_2Dt)[0]): - - plt.subplot(1,3,1) - im1 = plt.imshow(phantom_2Dt[sl,:,:], animated=True) - - plt.subplot(1,3,2) - im2 = plt.imshow(pdhg1.get_output().as_array()[sl,:,:]) - - plt.subplot(1,3,3) - im3 = plt.imshow(pdhg2.get_output().as_array()[sl,:,:]) - - ims1.append([im1]) - ims2.append([im2]) - ims3.append([im3]) - - -ani1 = animation.ArtistAnimation(fig, ims1, interval=500, - repeat_delay=10) - -ani2 = animation.ArtistAnimation(fig, ims2, interval=500, - repeat_delay=10) - -ani3 = animation.ArtistAnimation(fig, ims3, interval=500, - repeat_delay=10) -plt.show() -# plt.pause(0.25) -# plt.show() - - - - - - - - diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_Gaussian_3D.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_Gaussian_3D.py deleted file mode 100644 index 15709cd..0000000 --- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_Gaussian_3D.py +++ /dev/null @@ -1,147 +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 (3D) Denoising using PDHG algorithm and Tomophantom: - - -Problem: min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2} - - \alpha: Regularization parameter - - \nabla: Gradient operator - - g: 3D 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 Gradient -from ccpi.optimisation.functions import L2NormSquared, MixedL21Norm - -from skimage.util import random_noise - -import timeit -import os -from tomophantom import TomoP3D -import tomophantom - -# Create a phantom from 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. Apply 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) - -# Show results -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 - -# Setup and run the PDHG algorithm -operator = Gradient(ig) -f = alpha * MixedL21Norm() -g = 0.5 * L2NormSquared(b = noisy_data) - -normK = operator.norm() - -sigma = 1 -tau = 1/(sigma*normK**2) - -pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) -pdhg.max_iteration = 1000 -pdhg.update_objective_interval = 200 -pdhg.run(1000, verbose = True) - -# Show results -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/GatherAll/PDHG_TV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Tomo2D.py deleted file mode 100644 index 4f7639e..0000000 --- a/Wrappers/Python/demos/PDHG_examples/GatherAll/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/GatherAll/PDHG_Tikhonov_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py deleted file mode 100644 index 78d4980..0000000 --- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py +++ /dev/null @@ -1,258 +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. -# -#========================================================================= - -""" - -``Tikhonov`` Regularization Denoising using PDHG algorithm: - - -Problem: min_{u}, \alpha * ||\nabla u||_{2,2} + Fidelity(u, g) - - \alpha: Regularization parameter - - \nabla: Gradient operator - - g: Noisy Data - - Fidelity = 1) L2NormSquarred ( \frac{1}{2} * || u - g ||_{2}^{2} ) if Noise is Gaussian - 2) L1Norm ( ||u - g||_{1} )if Noise is Salt & Pepper - 3) Kullback Leibler (\int u - g * log(u) + Id_{u>0}) if Noise is Poisson - - Method = 0 ( PDHG - split ) : K = [ \nabla, - Identity] - - - Method = 1 (PDHG - explicit ): K = \nabla - - - Default: Tikhonov denoising - noise = Gaussian - Fidelity = L2NormSquarred - method = 0 - -""" - -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 = 2 -print ("Applying {} noise") - -if len(sys.argv) > 2: - method = sys.argv[2] -else: - method = '0' -print ("method ", method) - -loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) -data = loader.load(TestData.SHAPES) -ig = data.geometry -ag = ig - -# Create noisy data. -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': - scale = 5 - n1 = random_noise( data.as_array()/scale, mode = noise, seed = 10)*scale -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 Parameter depending on the noise distribution -if noise == 's&p': - alpha = 20 -elif noise == 'poisson': - alpha = 10 -elif noise == 'gaussian': - alpha = 5 - -# 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) - -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() - f = BlockFunction(f1, f2) - g = ZeroFunction() - -else: - - operator = Gradient(ig) - g = f2 - - -# 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 = 100 -pdhg.run(2000) - - -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('TV Reconstruction') -plt.colorbar() -plt.subplot(1,4,4) -plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[int(ig.shape[0]/2),:], label = 'GTruth') -plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(ig.shape[0]/2),:], label = 'Tikhonov reconstruction') -plt.legend() -plt.title('Middle Line Profiles') -plt.show() - - - -##%% Check with CVX solution - -from ccpi.optimisation.operators import SparseFiniteDiff - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - -if cvx_not_installable: - - ##Construct problem - u = Variable(ig.shape) - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - # Define Total Variation as a regulariser - - regulariser = alpha * sum_squares(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - - # choose solver - if 'MOSEK' in installed_solvers(): - solver = MOSEK - else: - solver = SCS - - # fidelity - if noise == 's&p': - fidelity = pnorm( u - noisy_data.as_array(),1) - elif noise == 'poisson': - fidelity = sum(kl_div(noisy_data.as_array(), u)) - solver = SCS - elif noise == 'gaussian': - fidelity = 0.5 * sum_squares(noisy_data.as_array() - u) - - 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,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(ig.shape[0]/2),:], label = 'PDHG') - plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), u.value[int(ig.shape[0]/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/IMATDemo.py b/Wrappers/Python/demos/PDHG_examples/IMATDemo.py deleted file mode 100644 index 2051860..0000000 --- a/Wrappers/Python/demos/PDHG_examples/IMATDemo.py +++ /dev/null @@ -1,339 +0,0 @@ - - -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Mon Mar 25 12:50:27 2019 - -@author: vaggelis -""" - -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, AcquisitionGeometry, AcquisitionData -from astropy.io import fits -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import CGLS, PDHG -from ccpi.optimisation.functions import MixedL21Norm, L2NormSquared, BlockFunction, ZeroFunction, KullbackLeibler, IndicatorBox -from ccpi.optimisation.operators import Gradient, BlockOperator - -from ccpi.astra.operators import AstraProjectorMC, AstraProjectorSimple - -import pickle - - -# load file - -#filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_282.fits' -#filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_564.fits' -#filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_141.fits' -filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_80_channels.fits' - -sino_handler = fits.open(filename_sino) -sino = numpy.array(sino_handler[0].data, dtype=float) - -# change axis order: channels, angles, detectors -sino_new = numpy.rollaxis(sino, 2) -sino_handler.close() - - -sino_shape = sino_new.shape - -num_channels = sino_shape[0] # channelss -num_pixels_h = sino_shape[2] # detectors -num_pixels_v = sino_shape[2] # detectors -num_angles = sino_shape[1] # angles - - -ig = ImageGeometry(voxel_num_x = num_pixels_h, voxel_num_y = num_pixels_v, channels = num_channels) - -with open("/media/newhd/vaggelis/CCPi/IMAT_reconstruction/CCPi-Framework/Wrappers/Python/ccpi/optimisation/IMAT_data/golden_angles_new.txt") as f: - angles_string = [line.rstrip() for line in f] - angles = numpy.array(angles_string).astype(float) - - -ag = AcquisitionGeometry('parallel', '2D', angles * numpy.pi / 180, pixel_num_h = num_pixels_h, channels = num_channels) -op_MC = AstraProjectorMC(ig, ag, 'gpu') - -sino_aqdata = AcquisitionData(sino_new, ag) -result_bp = op_MC.adjoint(sino_aqdata) - -#%% - -channel = [40, 60] -for j in range(2): - z4 = sino_aqdata.as_array()[channel[j]] - plt.figure(figsize=(10,6)) - plt.imshow(z4, cmap='viridis') - plt.axis('off') - plt.savefig('Sino_141/Sinogram_ch_{}_.png'.format(channel[j]), bbox_inches='tight', transparent=True) - plt.show() - -#%% - -def callback(iteration, objective, x): - plt.imshow(x.as_array()[40]) - plt.colorbar() - plt.show() - -#%% -# CGLS - -x_init = ig.allocate() -cgls1 = CGLS(x_init=x_init, operator=op_MC, data=sino_aqdata) -cgls1.max_iteration = 100 -cgls1.update_objective_interval = 2 -cgls1.run(20,verbose=True, callback=callback) - -plt.imshow(cgls1.get_output().subset(channel=20).array) -plt.title('CGLS') -plt.colorbar() -plt.show() - -#%% -with open('Sino_141/CGLS/CGLS_{}_iter.pkl'.format(20), 'wb') as f: - z = cgls1.get_output() - pickle.dump(z, f) - -#%% -#% Tikhonov Space - -x_init = ig.allocate() -alpha = [1,3,5,10,20,50] - -for a in alpha: - - Grad = Gradient(ig, correlation = Gradient.CORRELATION_SPACE) - operator = BlockOperator(op_MC, a * Grad, shape=(2,1)) - blockData = BlockDataContainer(sino_aqdata, \ - Grad.range_geometry().allocate()) - cgls2 = CGLS() - cgls2.max_iteration = 500 - cgls2.set_up(x_init, operator, blockData) - cgls2.update_objective_interval = 50 - cgls2.run(100,verbose=True) - - with open('Sino_141/CGLS_Space/CGLS_a_{}.pkl'.format(a), 'wb') as f: - z = cgls2.get_output() - pickle.dump(z, f) - -#% Tikhonov SpaceChannels - -for a1 in alpha: - - Grad1 = Gradient(ig, correlation = Gradient.CORRELATION_SPACECHANNEL) - operator1 = BlockOperator(op_MC, a1 * Grad1, shape=(2,1)) - blockData1 = BlockDataContainer(sino_aqdata, \ - Grad1.range_geometry().allocate()) - cgls3 = CGLS() - cgls3.max_iteration = 500 - cgls3.set_up(x_init, operator1, blockData1) - cgls3.update_objective_interval = 10 - cgls3.run(100, verbose=True) - - with open('Sino_141/CGLS_SpaceChannels/CGLS_a_{}.pkl'.format(a1), 'wb') as f1: - z1 = cgls3.get_output() - pickle.dump(z1, f1) - - - -#%% -# -ig_tmp = ImageGeometry(voxel_num_x = num_pixels_h, voxel_num_y = num_pixels_v) -ag_tmp = AcquisitionGeometry('parallel', '2D', angles * numpy.pi / 180, pixel_num_h = num_pixels_h) -op_tmp = AstraProjectorSimple(ig_tmp, ag_tmp, 'gpu') -normK1 = op_tmp.norm() - -alpha_TV = [2, 5, 10] # for powder - -# Create operators -op1 = Gradient(ig, correlation=Gradient.CORRELATION_SPACECHANNEL) -op2 = op_MC - -# Create BlockOperator -operator = BlockOperator(op1, op2, shape=(2,1) ) - - -for alpha in alpha_TV: -# Create functions - f1 = alpha * MixedL21Norm() - - f2 = KullbackLeibler(sino_aqdata) - f = BlockFunction(f1, f2) - g = IndicatorBox(lower=0) - - # Compute operator Norm - normK = numpy.sqrt(8 + normK1**2) - - # 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 = 5000 - pdhg.update_objective_interval = 500 -# pdhg.run(2000, verbose=True, callback=callback) - pdhg.run(5000, verbose=True, callback=callback) -# - with open('Sino_141/TV_SpaceChannels/TV_a = {}.pkl'.format(alpha), 'wb') as f3: - z3 = pdhg.get_output() - pickle.dump(z3, f3) -# -# -# -# -##%% -# -#ig_tmp = ImageGeometry(voxel_num_x = num_pixels_h, voxel_num_y = num_pixels_v) -#ag_tmp = AcquisitionGeometry('parallel', '2D', angles * numpy.pi / 180, pixel_num_h = num_pixels_h) -#op_tmp = AstraProjectorSimple(ig_tmp, ag_tmp, 'gpu') -#normK1 = op_tmp.norm() -# -#alpha_TV = 10 # for powder -# -## Create operators -#op1 = Gradient(ig, correlation=Gradient.CORRELATION_SPACECHANNEL) -#op2 = op_MC -# -## Create BlockOperator -#operator = BlockOperator(op1, op2, shape=(2,1) ) -# -# -## Create functions -#f1 = alpha_TV * MixedL21Norm() -#f2 = 0.5 * L2NormSquared(b=sino_aqdata) -#f = BlockFunction(f1, f2) -#g = ZeroFunction() -# -## Compute operator Norm -##normK = 8.70320267279591 # For powder Run one time no need to compute again takes time -#normK = numpy.sqrt(8 + normK1**2) # for carbon -# -## Primal & dual stepsizes -#sigma = 0.1 -#tau = 1/(sigma*normK**2) -# -#def callback(iteration, objective, x): -# plt.imshow(x.as_array()[100]) -# plt.colorbar() -# plt.show() -# -## 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 = 100 -#pdhg.run(2000, verbose=True) -# -# -# -# -# -# - - - - - - - - - - -#%% - -#with open('/media/newhd/vaggelis/CCPi/IMAT_reconstruction/CCPi-Framework/Wrappers/Python/ccpi/optimisation/CGLS_Tikhonov/CGLS_Space/CGLS_Space_a = 50.pkl', 'wb') as f: -# z = cgls2.get_output() -# pickle.dump(z, f) -# - - #%% -with open('Sino_141/CGLS_Space/CGLS_Space_a_20.pkl', 'rb') as f1: - x = pickle.load(f1) - -with open('Sino_141/CGLS_SpaceChannels/CGLS_SpaceChannels_a_20.pkl', 'rb') as f1: - x1 = pickle.load(f1) - - - -# -plt.imshow(x.as_array()[40]*mask) -plt.colorbar() -plt.show() - -plt.imshow(x1.as_array()[40]*mask) -plt.colorbar() -plt.show() - -plt.plot(x.as_array()[40,100,:]) -plt.plot(x1.as_array()[40,100,:]) -plt.show() - -#%% - -# Show results - -def circ_mask(h, w, center=None, radius=None): - - if center is None: # use the middle of the image - center = [int(w/2), int(h/2)] - if radius is None: # use the smallest distance between the center and image walls - radius = min(center[0], center[1], w-center[0], h-center[1]) - - Y, X = numpy.ogrid[:h, :w] - dist_from_center = numpy.sqrt((X - center[0])**2 + (Y-center[1])**2) - - mask = dist_from_center <= radius - return mask - -mask = circ_mask(141, 141, center=None, radius = 55) -plt.imshow(numpy.multiply(x.as_array()[40],mask)) -plt.show() -#%% -#channel = [100, 200, 300] -# -#for i in range(3): -# tmp = cgls1.get_output().as_array()[channel[i]] -# -# z = tmp * mask -# plt.figure(figsize=(10,6)) -# plt.imshow(z, vmin=0, cmap='viridis') -# plt.axis('off') -## plt.clim(0, 0.02) -## plt.colorbar() -## del z -# plt.savefig('CGLS_282/CGLS_Chan_{}.png'.format(channel[i]), bbox_inches='tight', transparent=True) -# plt.show() -# -# -##%% Line Profiles -# -#n1, n2, n3 = cgs.get_output().as_array().shape -#mask = circ_mask(564, 564, center=None, radius = 220) -#material = ['Cu', 'Fe', 'Ni'] -#ycoords = [200, 300, 380] -# -#for i in range(3): -# z = cgs.get_output().as_array()[channel[i]] * mask -# -# for k1 in range(len(ycoords)): -# plt.plot(numpy.arange(0,n2), z[ycoords[k1],:]) -# plt.title('Channel {}: {}'.format(channel[i], material[k1])) -# plt.savefig('CGLS/line_profile_chan_{}_material_{}.png'.\ -# format(channel[i], material[k1]), bbox_inches='tight') -# plt.show() -# -# -# -# -# -##%% -# -#%% - - - -#%% - -#plt.imshow(pdhg.get_output().subset(channel=100).as_array()) -#plt.show() 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 deleted file mode 100644 index 045458a..0000000 --- a/Wrappers/Python/demos/PDHG_examples/MultiChannel/PDHG_2D_time_denoising.py +++ /dev/null @@ -1,169 +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, 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 deleted file mode 100644 index 045458a..0000000 --- a/Wrappers/Python/demos/PDHG_examples/MultiChannel/PDHG_TV_Tomo2D_time.py +++ /dev/null @@ -1,169 +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, 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/PDHG_TV_Color_Denoising.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Color_Denoising.py deleted file mode 100644 index ddf5ace..0000000 --- a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Color_Denoising.py +++ /dev/null @@ -1,115 +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 - - -""" - -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG - -from ccpi.optimisation.operators import Gradient, BlockOperator, FiniteDiff -from ccpi.optimisation.functions import MixedL21Norm, MixedL11Norm, L2NormSquared, BlockFunction, L1Norm -from ccpi.framework import TestData, ImageGeometry -import os, sys -if int(numpy.version.version.split('.')[1]) > 12: - from skimage.util import random_noise -else: - from demoutil import random_noise - -loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) -data = loader.load(TestData.PEPPERS, size=(256,256)) -ig = data.geometry -ag = ig - -# Create noisy data. -n1 = random_noise(data.as_array(), mode = 'gaussian', var = 0.15, seed = 50) -noisy_data = ig.allocate() -noisy_data.fill(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 -operator = Gradient(ig, correlation=Gradient.CORRELATION_SPACE) -f1 = 5 * 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 -pdhg1 = PDHG(f=f1,g=g,operator=operator, tau=tau, sigma=sigma) -pdhg1.max_iteration = 2000 -pdhg1.update_objective_interval = 200 -pdhg1.run(1000) - - -# Show results -plt.figure(figsize=(10,10)) -plt.subplot(2,2,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(2,2,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.subplot(2,2,3) -plt.imshow(pdhg1.get_output().as_array()) -plt.title('TV Reconstruction') -plt.colorbar() -plt.subplot(2,2,4) -plt.imshow(pdhg2.get_output().as_array()) -plt.title('TV Reconstruction') -plt.colorbar() -plt.show() - 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 deleted file mode 100644 index 14608db..0000000 --- a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_2D_time.py +++ /dev/null @@ -1,192 +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. -# -#========================================================================= - -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_3D.py b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Gaussian_3D.py deleted file mode 100644 index 03dc2ef..0000000 --- a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Gaussian_3D.py +++ /dev/null @@ -1,181 +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 (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/Tomo/PDHG_Tikhonov_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_Tikhonov_Tomo2D.py deleted file mode 100644 index 02cd053..0000000 --- a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_Tikhonov_Tomo2D.py +++ /dev/null @@ -1,156 +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}^{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() - - diff --git a/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py b/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py deleted file mode 100644 index 854f645..0000000 --- a/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py +++ /dev/null @@ -1,238 +0,0 @@ -# -*- coding: utf-8 -*- - -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, \ - AcquisitionGeometry, AcquisitionData - -import numpy as np -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction, ScaledFunction - -from ccpi.plugins.operators import CCPiProjectorSimple -from timeit import default_timer as timer -from ccpi.reconstruction.parallelbeam import alg as pbalg -import os - -try: - import tomophantom - from tomophantom import TomoP3D - no_tomophantom = False -except ImportError as ie: - no_tomophantom = True - -#%% - -#%%############################################################################### -# Create phantom for TV tomography - -#import os -#import tomophantom -#from tomophantom import TomoP2D -#from tomophantom.supp.qualitymetrics import QualityTools - -#model = 1 # select a model number from the library -#N = 150 # 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 phantom (2D) -#phantom_2D = TomoP2D.Model(model, N, path_library2D) -#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -#data = ImageData(phantom_2D, geometry=ig) - -N = 75 -#x = np.zeros((N,N)) - -vert = 4 -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, voxel_num_z=vert) - -angles_num = 100 -det_w = 1.0 -det_num = N - -angles = np.linspace(-90.,90.,N, dtype=np.float32) -# Inputs: Geometry, 2D or 3D, angles, horz detector pixel count, -# horz detector pixel size, vert detector pixel count, -# vert detector pixel size. -ag = AcquisitionGeometry('parallel', - '3D', - angles, - N, - det_w, - vert, - det_w) - -#no_tomophantom = True -if no_tomophantom: - data = ig.allocate() - Phantom = data - # Populate image data by looping over and filling slices - i = 0 - while i < vert: - if vert > 1: - x = Phantom.subset(vertical=i).array - else: - x = Phantom.array - 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)] = 0.98 - if vert > 1 : - Phantom.fill(x, vertical=i) - i += 1 - - Aop = CCPiProjectorSimple(ig, ag, 'cpu') - sin = Aop.direct(data) -else: - - model = 13 # select a model number from the library - N_size = N # 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_size x N_size x N_size phantom (3D) - phantom_tm = TomoP3D.Model(model, N_size, path_library3D) - - #%% - Horiz_det = int(np.sqrt(2)*N_size) # detector column count (horizontal) - Vert_det = N_size # detector row count (vertical) (no reason for it to be > N) - #angles_num = int(0.5*np.pi*N_size); # angles number - #angles = np.linspace(0.0,179.9,angles_num,dtype='float32') # in degrees - - print ("Building 3D analytical projection data with TomoPhantom") - projData3D_analyt = TomoP3D.ModelSino(model, - N_size, - Horiz_det, - Vert_det, - angles, - path_library3D) - - # tomophantom outputs in [vert,angles,horiz] - # we want [angle,vert,horiz] - data = np.transpose(projData3D_analyt, [1,0,2]) - ag.pixel_num_h = Horiz_det - ag.pixel_num_v = Vert_det - sin = ag.allocate() - sin.fill(data) - ig.voxel_num_y = Vert_det - - Aop = CCPiProjectorSimple(ig, ag, 'cpu') - - -plt.imshow(sin.subset(vertical=0).as_array()) -plt.title('Sinogram') -plt.colorbar() -plt.show() - - -#%% -# Add Gaussian noise to the sinogram data -np.random.seed(10) -n1 = np.random.random(sin.shape) - -noisy_data = sin + ImageData(5*n1) - -plt.imshow(noisy_data.subset(vertical=0).as_array()) -plt.title('Noisy Sinogram') -plt.colorbar() -plt.show() - - -#%% Works only with Composite Operator Structure of PDHG - -#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) - -# Create operators -op1 = Gradient(ig) -op2 = Aop - -# Form Composite Operator -operator = BlockOperator(op1, op2, shape=(2,1) ) - -alpha = 50 -f = BlockFunction( alpha * MixedL21Norm(), \ - 0.5 * L2NormSquared(b = noisy_data) ) -g = ZeroFunction() - -normK = Aop.norm() - -# Compute operator Norm -normK = operator.norm() - -## Primal & dual stepsizes -diag_precon = False - -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: - sigma = 1 - tau = 1/(sigma*normK**2) - -# Compute operator Norm -normK = operator.norm() - -# Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) -niter = 50 -opt = {'niter':niter} -opt1 = {'niter':niter, 'memopt': True} - - - -pdhg1 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma, max_iteration=niter) -#pdhg1.max_iteration = 2000 -pdhg1.update_objective_interval = 100 - -t1_old = timer() -resold, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -t2_old = timer() - -pdhg1.run(niter) -print (sum(pdhg1.timing)) -res = pdhg1.get_output().subset(vertical=0) - -#%% -plt.figure() -plt.subplot(1,4,1) -plt.imshow(res.as_array()) -plt.title('Algorithm') -plt.colorbar() -plt.subplot(1,4,2) -plt.imshow(resold.subset(vertical=0).as_array()) -plt.title('function') -plt.colorbar() -plt.subplot(1,4,3) -plt.imshow((res - resold.subset(vertical=0)).abs().as_array()) -plt.title('diff') -plt.colorbar() -plt.subplot(1,4,4) -plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'Algorithm') -plt.plot(np.linspace(0,N,N), resold.subset(vertical=0).as_array()[int(N/2),:], label = 'function') -plt.legend() -plt.show() -# -print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(sum(pdhg1.timing), t2_old -t1_old)) -diff = (res - resold.subset(vertical=0)).abs().as_array().max() -# -print(" Max of abs difference is {}".format(diff)) - |