diff options
Diffstat (limited to 'Wrappers/Python')
18 files changed, 0 insertions, 3484 deletions
diff --git a/Wrappers/Python/demos/PDHG_TGV_Tomo2D.py b/Wrappers/Python/demos/PDHG_TGV_Tomo2D.py deleted file mode 100644 index 19cf684..0000000 --- a/Wrappers/Python/demos/PDHG_TGV_Tomo2D.py +++ /dev/null @@ -1,127 +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, Identity, \ - SymmetrizedGradient, ZeroOperator -from ccpi.optimisation.functions import ZeroFunction, IndicatorBox, KullbackLeibler, \ - MixedL21Norm, BlockFunction - -from ccpi.astra.ops import AstraProjectorSimple - -# Create phantom for TV 2D tomography -N = 75 - -data = np.zeros((N,N)) - -x1 = np.linspace(0, int(N/2), N) -x2 = np.linspace(int(N/2), 0., N) -xv, yv = np.meshgrid(x1, x2) - -xv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1] = yv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1].T -data = xv -data = ImageData(data/data.max()) - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) - -detectors = N -angles = np.linspace(0, np.pi, N, dtype=np.float32) - -ag = AcquisitionGeometry('parallel','2D',angles, detectors) -Aop = AstraProjectorSimple(ig, ag, 'gpu') -sin = Aop.direct(data) - -# Create noisy data. Apply Poisson noise -scale = 0.1 -np.random.seed(5) -n1 = scale * np.random.poisson(sin.as_array()/scale) -noisy_data = AcquisitionData(n1, ag) - - -plt.imshow(noisy_data.as_array()) -plt.show() -#%% -# Regularisation Parameters -alpha = 0.7 -beta = 2 - -# Create Operators -op11 = Gradient(ig) -op12 = Identity(op11.range_geometry()) - -op22 = SymmetrizedGradient(op11.domain_geometry()) -op21 = ZeroOperator(ig, op22.range_geometry()) - -op31 = Aop -op32 = ZeroOperator(op22.domain_geometry(), ag) - -operator = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) - -f1 = alpha * MixedL21Norm() -f2 = beta * MixedL21Norm() -f3 = KullbackLeibler(noisy_data) -f = BlockFunction(f1, f2, f3) - -g = BlockFunction(-1 * IndicatorBox(lower=0), ZeroFunction()) -#g = IndicatorBox(lower=0) -#g = ZeroFunction() - -# Compute operator Norm -normK = operator.norm() - -# Primal & dual stepsizes -sigma = 10 -tau = 1/(sigma*normK**2) - - -# Setup and run the PDHG algorithm -pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) -pdhg.max_iteration = 2000 -pdhg.update_objective_interval = 50 -pdhg.run(2000) - -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow(pdhg.get_output()[0].as_array()) -plt.title('TGV Reconstruction') -plt.colorbar() -plt.show() -## -plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') -plt.plot(np.linspace(0,N,N), pdhg.get_output()[0].as_array()[int(N/2),:], label = 'TGV reconstruction') -plt.legend() -plt.title('Middle Line Profiles') -plt.show() - - diff --git a/Wrappers/Python/demos/PDHG_TV_Tomo2D_poisson.py b/Wrappers/Python/demos/PDHG_TV_Tomo2D_poisson.py deleted file mode 100644 index 72d0670..0000000 --- a/Wrappers/Python/demos/PDHG_TV_Tomo2D_poisson.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. -# -#========================================================================= - -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, IndicatorBox - -from ccpi.astra.ops import AstraProjectorSimple - -""" - -Total Variation Denoising using PDHG algorithm: - - -Problem: min_x, x>0 \alpha * ||\nabla x||_{2,1} + int A x -g log(Ax + \eta) - - \nabla: Gradient operator - - A: Projection Matrix - g: Noisy sinogram corrupted with Poisson Noise - - \eta: Background Noise - \alpha: Regularization parameter - -""" - -# Create phantom for TV 2D tomography -N = 50 -x = np.zeros((N,N)) -x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 -x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 - -data = ImageData(x) -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) - -detectors = N -angles = np.linspace(0, np.pi, N) - -ag = AcquisitionGeometry('parallel','2D',angles, detectors) -Aop = AstraProjectorSimple(ig, ag, 'cpu') -sin = Aop.direct(data) - -# Create noisy data. Apply Poisson noise -scale = 0.25 -eta = 0 #np.random.randint(0, sin.as_array().max()/2, ag.shape) -n1 = scale * np.random.poisson(eta + sin.as_array()/scale) - -noisy_data = AcquisitionData(n1, ag) - - -# Show Ground Truth and Noisy Data -plt.figure(figsize=(10,10)) -plt.subplot(2,1,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(2,1,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.show() - -#%% - -# Regularisation Parameter -alpha = 2 - -# Create operators -op1 = Gradient(ig) -op2 = Aop - -# Create BlockOperator -operator = BlockOperator(op1, op2, shape=(2,1) ) - -# Create functions - -f1 = alpha * MixedL21Norm() -f2 = KullbackLeibler(noisy_data) -f = BlockFunction(f1, f2) - -#g = ZeroFunction() -g = IndicatorBox(lower=0) - -# Compute operator Norm -normK = operator.norm() - -# Primal & dual stepsizes -sigma = 2 -tau = 1/(sigma*normK**2) -#sigma = 1/normK -#tau = 1/normK - -# 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, verbose = True) - -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow(pdhg.get_output().as_array()) -plt.title('TV Reconstruction') -plt.colorbar() -plt.show() -## -plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') -plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction') -plt.legend() -plt.title('Middle Line Profiles') -plt.show() - - -##%% Check with CVX solution -# -#from ccpi.optimisation.operators import SparseFiniteDiff -#import astra -#import numpy -# -#try: -# from cvxpy import * -# cvx_not_installable = True -#except ImportError: -# cvx_not_installable = False -# -# -#if cvx_not_installable: -# -# -# ##Construct problem -# u = Variable(N*N) -# #q = Variable() -# -# DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') -# DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') -# -# regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) -# -# # create matrix representation for Astra operator -# -# vol_geom = astra.create_vol_geom(N, N) -# proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) -# -# proj_id = astra.create_projector('strip', proj_geom, vol_geom) -# -# matrix_id = astra.projector.matrix(proj_id) -# -# ProjMat = astra.matrix.get(matrix_id) -# -# tmp = noisy_data.as_array().ravel('F') -# -## fidelity = sum( ProjMat * u - tmp * log(ProjMat * u + 1e-6)) -# #constraints = [q>= fidelity, u>=0] -# constraints = [] -# -# fidelity = sum(kl_div(tmp, ProjMat * u + 1e-6)) -## fidelity = kl_div(cp.multiply(alpha, W), -## cp.multiply(alpha, W + cp.multiply(beta, P))) - \ -## cp.multiply(alpha, cp.multiply(beta, P)) -# -# -# -# solver = SCS -# obj = Minimize( regulariser + fidelity) -# prob = Problem(obj, constraints) -# result = prob.solve(verbose = True, solver = solver) -# -# -####%% Check with CVX solution -## -##from ccpi.optimisation.operators import SparseFiniteDiff -## -##try: -## from cvxpy import * -## cvx_not_installable = True -##except ImportError: -## cvx_not_installable = False -## -## -##if cvx_not_installable: -## -## ##Construct problem -## u = Variable(ig.shape) -## -## DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') -## DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') -## -## # Define Total Variation as a regulariser -## regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) -## fidelity = pnorm( u - noisy_data.as_array(),1) -## -## # choose solver -## if 'MOSEK' in installed_solvers(): -## solver = MOSEK -## else: -## solver = SCS -## -## obj = Minimize( regulariser + fidelity) -## prob = Problem(obj) -## result = prob.solve(verbose = True, solver = solver) -## -## -# plt.figure(figsize=(15,15)) -# plt.subplot(3,1,1) -# plt.imshow(pdhg.get_output().as_array()) -# plt.title('PDHG solution') -# plt.colorbar() -# plt.subplot(3,1,2) -# plt.imshow(np.reshape(u.value, (N, N))) -# plt.title('CVX solution') -# plt.colorbar() -# plt.subplot(3,1,3) -# plt.imshow(diff_cvx) -# plt.title('Difference') -# plt.colorbar() -# plt.show() -# -# plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') -# plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') -# plt.legend() -# plt.title('Middle Line Profiles') -# plt.show() -# -# print('Primal Objective (CVX) {} '.format(obj.value)) -# print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
\ No newline at end of file diff --git a/Wrappers/Python/demos/PDHG_TV_Tomo2D_time.py b/Wrappers/Python/demos/PDHG_TV_Tomo2D_time.py deleted file mode 100644 index 045458a..0000000 --- a/Wrappers/Python/demos/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_Tikhonov_Tomo2D.py b/Wrappers/Python/demos/PDHG_Tikhonov_Tomo2D.py deleted file mode 100644 index 02cd053..0000000 --- a/Wrappers/Python/demos/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_examples/.DS_Store b/Wrappers/Python/demos/PDHG_examples/.DS_Store Binary files differnew file mode 100644 index 0000000..141508d --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/.DS_Store diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_2D_time_denoising.py b/Wrappers/Python/demos/PDHG_examples/PDHG_2D_time_denoising.py deleted file mode 100644 index 045458a..0000000 --- a/Wrappers/Python/demos/PDHG_examples/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/PDHG_TGV_Denoising.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising.py deleted file mode 100755 index 4d6da00..0000000 --- a/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising.py +++ /dev/null @@ -1,221 +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_{x} \alpha * ||\nabla x - w||_{2,1} + - \beta * || E w ||_{2,1} + - fidelity - - where fidelity can be as follows depending on the noise characteristics - of the data: - * Norm2Squared \frac{1}{2} * || x - g ||_{2}^{2} - * KullbackLeibler \int u - g * log(u) + Id_{u>0} - * L1Norm ||u - g||_{1} - - \alpha: Regularization parameter - \beta: Regularization parameter - - \nabla: Gradient operator - E: Symmetrized Gradient operator - - g: Noisy Data with Salt & Pepper Noise - - Method = 0 ( PDHG - split ) : K = [ \nabla, - Identity - ZeroOperator, E - Identity, ZeroOperator] - - - Method = 1 (PDHG - explicit ): K = [ \nabla, - Identity - ZeroOperator, E ] - -""" - -from ccpi.framework import ImageData, ImageGeometry - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG - -from ccpi.optimisation.operators import BlockOperator, Identity, \ - Gradient, SymmetrizedGradient, ZeroOperator -from ccpi.optimisation.functions import ZeroFunction, L1Norm, \ - MixedL21Norm, BlockFunction, KullbackLeibler, L2NormSquared -import sys -if int(numpy.version.version.split('.')[1]) > 12: - from skimage.util import random_noise -else: - from demoutil import random_noise - - - -# user supplied input -if len(sys.argv) > 1: - which_noise = int(sys.argv[1]) -else: - which_noise = 0 -print ("Applying {} noise") - -if len(sys.argv) > 2: - method = sys.argv[2] -else: - method = '1' -print ("method ", method) -# Create phantom for TGV SaltPepper denoising - -N = 100 - -x1 = np.linspace(0, int(N/2), N) -x2 = np.linspace(int(N/2), 0., N) -xv, yv = np.meshgrid(x1, x2) - -xv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1] = yv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1].T - -#data = xv - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -data = ig.allocate() -data.fill(xv/xv.max()) -ag = ig - -# Create noisy data. -# Apply Salt & Pepper noise -# gaussian -# poisson -noises = ['gaussian', 'poisson', 's&p'] -noise = noises[which_noise] -if noise == 's&p': - n1 = random_noise(data.as_array(), mode = noise, salt_vs_pepper = 0.9, amount=0.2, seed=10) -elif noise == 'poisson': - n1 = random_noise(data.as_array(), mode = noise, seed = 10) -elif noise == 'gaussian': - n1 = random_noise(data.as_array(), mode = noise, seed = 10) -else: - raise ValueError('Unsupported Noise ', noise) -noisy_data = ImageData(n1) - -# Show Ground Truth and Noisy Data -plt.figure(figsize=(10,5)) -plt.subplot(1,2,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(1,2,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.show() - -# Regularisation Parameters -if noise == 's&p': - alpha = 0.8 -elif noise == 'poisson': - alpha = .1 -elif noise == 'gaussian': - alpha = .3 - -beta = numpy.sqrt(2)* alpha - -# fidelity -if noise == 's&p': - f3 = L1Norm(b=noisy_data) -elif noise == 'poisson': - f3 = KullbackLeibler(noisy_data) -elif noise == 'gaussian': - f3 = L2NormSquared(b=noisy_data) - -if method == '0': - - # Create operators - op11 = Gradient(ig) - op12 = Identity(op11.range_geometry()) - - op22 = SymmetrizedGradient(op11.domain_geometry()) - op21 = ZeroOperator(ig, op22.range_geometry()) - - op31 = Identity(ig, ag) - op32 = ZeroOperator(op22.domain_geometry(), ag) - - operator = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) - - f1 = alpha * MixedL21Norm() - f2 = beta * MixedL21Norm() - # f3 depends on the noise characteristics - - f = BlockFunction(f1, f2, f3) - g = ZeroFunction() - -else: - - # Create operators - op11 = Gradient(ig) - op12 = Identity(op11.range_geometry()) - op22 = SymmetrizedGradient(op11.domain_geometry()) - op21 = ZeroOperator(ig, op22.range_geometry()) - - operator = BlockOperator(op11, -1*op12, op21, op22, shape=(2,2) ) - - f1 = alpha * MixedL21Norm() - f2 = beta * MixedL21Norm() - - f = BlockFunction(f1, f2) - g = BlockFunction(f3, ZeroFunction()) - -## Compute operator Norm -normK = operator.norm() -# -# Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) - - -# Setup and run the PDHG algorithm -pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) -pdhg.max_iteration = 2000 -pdhg.update_objective_interval = 200 -pdhg.run(2000, verbose = True) - -#%% -plt.figure(figsize=(20,5)) -plt.subplot(1,4,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(1,4,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.subplot(1,4,3) -plt.imshow(pdhg.get_output()[0].as_array()) -plt.title('TGV Reconstruction') -plt.colorbar() -plt.subplot(1,4,4) -plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') -plt.plot(np.linspace(0,N,N), pdhg.get_output()[0].as_array()[int(N/2),:], label = 'TGV reconstruction') -plt.legend() -plt.title('Middle Line Profiles') -plt.show() - - diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising_SaltPepper.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising_SaltPepper.py deleted file mode 100644 index 15c0a05..0000000 --- a/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising_SaltPepper.py +++ /dev/null @@ -1,245 +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_{x} \alpha * ||\nabla x - w||_{2,1} + - \beta * || E w ||_{2,1} + - \frac{1}{2} * || x - g ||_{2}^{2} - - \alpha: Regularization parameter - \beta: Regularization parameter - - \nabla: Gradient operator - E: Symmetrized Gradient operator - - g: Noisy Data with Salt & Pepper Noise - - Method = 0 ( PDHG - split ) : K = [ \nabla, - Identity - ZeroOperator, E - Identity, ZeroOperator] - - - Method = 1 (PDHG - explicit ): K = [ \nabla, - Identity - ZeroOperator, E ] - -""" - -from ccpi.framework import ImageData, ImageGeometry - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG - -from ccpi.optimisation.operators import BlockOperator, Identity, \ - Gradient, SymmetrizedGradient, ZeroOperator -from ccpi.optimisation.functions import ZeroFunction, L1Norm, \ - MixedL21Norm, BlockFunction - -from skimage.util import random_noise - -# Create phantom for TGV SaltPepper denoising - -N = 100 - -data = np.zeros((N,N)) - -x1 = np.linspace(0, int(N/2), N) -x2 = np.linspace(int(N/2), 0., N) -xv, yv = np.meshgrid(x1, x2) - -xv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1] = yv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1].T - -data = xv -data = ImageData(data/data.max()) - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -ag = ig - -# Create noisy data. Add Gaussian noise -n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2) -noisy_data = ImageData(n1) - -# Show Ground Truth and Noisy Data -plt.figure(figsize=(15,15)) -plt.subplot(2,1,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(2,1,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.show() - -# Regularisation Parameters -alpha = 0.8 -beta = numpy.sqrt(2)* alpha - -method = '1' - -if method == '0': - - # Create operators - op11 = Gradient(ig) - op12 = Identity(op11.range_geometry()) - - op22 = SymmetrizedGradient(op11.domain_geometry()) - op21 = ZeroOperator(ig, op22.range_geometry()) - - op31 = Identity(ig, ag) - op32 = ZeroOperator(op22.domain_geometry(), ag) - - operator = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) - - f1 = alpha * MixedL21Norm() - f2 = beta * MixedL21Norm() - f3 = L1Norm(b=noisy_data) - f = BlockFunction(f1, f2, f3) - g = ZeroFunction() - -else: - - # Create operators - op11 = Gradient(ig) - op12 = Identity(op11.range_geometry()) - op22 = SymmetrizedGradient(op11.domain_geometry()) - op21 = ZeroOperator(ig, op22.range_geometry()) - - operator = BlockOperator(op11, -1*op12, op21, op22, shape=(2,2) ) - - f1 = alpha * MixedL21Norm() - f2 = beta * MixedL21Norm() - - f = BlockFunction(f1, f2) - g = BlockFunction(L1Norm(b=noisy_data), ZeroFunction()) - -## Compute operator Norm -normK = operator.norm() -# -# Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) - - -# Setup and run the PDHG algorithm -pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) -pdhg.max_iteration = 2000 -pdhg.update_objective_interval = 50 -pdhg.run(2000, verbose = False) - -#%% -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow(pdhg.get_output()[0].as_array()) -plt.title('TGV Reconstruction') -plt.colorbar() -plt.show() -## -plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') -plt.plot(np.linspace(0,N,N), pdhg.get_output()[0].as_array()[int(N/2),:], label = 'TV reconstruction') -plt.legend() -plt.title('Middle Line Profiles') -plt.show() - - -#%% Check with CVX solution - -from ccpi.optimisation.operators import SparseFiniteDiff - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - -if cvx_not_installable: - - u = Variable(ig.shape) - w1 = Variable((N, N)) - w2 = Variable((N, N)) - - # create TGV regulariser - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u) - vec(w1), \ - DY.matrix() * vec(u) - vec(w2)]), 2, axis = 0)) + \ - beta * sum(norm(vstack([ DX.matrix().transpose() * vec(w1), DY.matrix().transpose() * vec(w2), \ - 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ), \ - 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ) ]), 2, axis = 0 ) ) - - constraints = [] - fidelity = pnorm(u - noisy_data.as_array(),1) - solver = MOSEK - - # choose solver - if 'MOSEK' in installed_solvers(): - solver = MOSEK - else: - solver = SCS - - obj = Minimize( regulariser + fidelity) - prob = Problem(obj) - result = prob.solve(verbose = True, solver = solver) - - diff_cvx = numpy.abs( pdhg.get_output()[0].as_array() - u.value ) - - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(pdhg.get_output()[0].as_array()) - plt.title('PDHG solution') - plt.colorbar() - plt.subplot(3,1,2) - plt.imshow(u.value) - plt.title('CVX solution') - plt.colorbar() - plt.subplot(3,1,3) - plt.imshow(diff_cvx) - plt.title('Difference') - plt.colorbar() - plt.show() - - plt.plot(np.linspace(0,N,N), pdhg.get_output()[0].as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') - plt.legend() - plt.title('Middle Line Profiles') - plt.show() - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) - - - - - diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising.py deleted file mode 100755 index 0f1effa..0000000 --- a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising.py +++ /dev/null @@ -1,266 +0,0 @@ -#======================================================================== -# Copyright 2019 Science Technology Facilities Council -# Copyright 2019 University of Manchester -# -# This work is part of the Core Imaging Library developed by Science Technology -# Facilities Council and University of Manchester -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0.txt -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -#========================================================================= - -""" - -Total Variation Denoising using PDHG algorithm: - - -Problem: min_x, x>0 \alpha * ||\nabla x||_{2,1} + ||x-g||_{1} - - \alpha: Regularization parameter - - \nabla: Gradient operator - - g: Noisy Data with Salt & Pepper Noise - - - Method = 0 ( PDHG - split ) : K = [ \nabla, - Identity] - - - Method = 1 (PDHG - explicit ): K = \nabla - - -""" - -from ccpi.framework import ImageData, ImageGeometry - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, L1Norm, \ - MixedL21Norm, BlockFunction, L2NormSquared,\ - KullbackLeibler -from ccpi.framework import TestData -import os, sys -if int(numpy.version.version.split('.')[1]) > 12: - from skimage.util import random_noise -else: - from demoutil import random_noise - -# user supplied input -if len(sys.argv) > 1: - which_noise = int(sys.argv[1]) -else: - which_noise = 0 -print ("Applying {} noise") - -if len(sys.argv) > 2: - method = sys.argv[2] -else: - method = '0' -print ("method ", method) -# Create phantom for TV Salt & Pepper denoising -N = 100 - -loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) -data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,N)) -data = loader.load(TestData.PEPPERS, size=(N,N)) -ig = data.geometry -ag = ig - -# Create noisy data. -# Apply Salt & Pepper noise -# gaussian -# poisson -noises = ['gaussian', 'poisson', 's&p'] -noise = noises[which_noise] -if noise == 's&p': - n1 = random_noise(data.as_array(), mode = noise, salt_vs_pepper = 0.9, amount=0.2) -elif noise == 'poisson': - n1 = random_noise(data.as_array(), mode = noise, seed = 10) -elif noise == 'gaussian': - n1 = random_noise(data.as_array(), mode = noise, seed = 10) -else: - raise ValueError('Unsupported Noise ', noise) -noisy_data = ig.allocate() -noisy_data.fill(n1) -#noisy_data = ImageData(n1) - -# Show Ground Truth and Noisy Data -plt.figure(figsize=(10,5)) -plt.subplot(1,2,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(1,2,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.show() - -# Regularisation Parameter -alpha = .2 - -# fidelity -if noise == 's&p': - f2 = L1Norm(b=noisy_data) -elif noise == 'poisson': - f2 = KullbackLeibler(noisy_data) -elif noise == 'gaussian': - f2 = L2NormSquared(b=noisy_data) - -if method == '0': - - # Create operators - op1 = Gradient(ig, correlation=Gradient.CORRELATION_SPACECHANNEL) - op2 = Identity(ig, ag) - - # Create BlockOperator - operator = BlockOperator(op1, op2, shape=(2,1) ) - - # Create functions - f1 = alpha * MixedL21Norm() - #f2 = L1Norm(b = noisy_data) - f = BlockFunction(f1, f2) - - g = ZeroFunction() - -else: - - # Without the "Block Framework" - operator = Gradient(ig) - f = alpha * MixedL21Norm() - #g = L1Norm(b = noisy_data) - g = f2 - - -# Compute operator Norm -normK = operator.norm() - -# Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) -opt = {'niter':2000, 'memopt': True} - -# Setup and run the PDHG algorithm -pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) -pdhg.max_iteration = 2000 -pdhg.update_objective_interval = 50 -pdhg.run(2000) - -if data.geometry.channels > 1: - plt.figure(figsize=(20,15)) - for row in range(data.geometry.channels): - - plt.subplot(3,4,1+row*4) - plt.imshow(data.subset(channel=row).as_array()) - plt.title('Ground Truth') - plt.colorbar() - plt.subplot(3,4,2+row*4) - plt.imshow(noisy_data.subset(channel=row).as_array()) - plt.title('Noisy Data') - plt.colorbar() - plt.subplot(3,4,3+row*4) - plt.imshow(pdhg.get_output().subset(channel=row).as_array()) - plt.title('TV Reconstruction') - plt.colorbar() - plt.subplot(3,4,4+row*4) - plt.plot(np.linspace(0,N,N), data.subset(channel=row).as_array()[int(N/2),:], label = 'GTruth') - plt.plot(np.linspace(0,N,N), pdhg.get_output().subset(channel=row).as_array()[int(N/2),:], label = 'TV reconstruction') - plt.legend() - plt.title('Middle Line Profiles') - plt.show() -else: - plt.figure(figsize=(20,5)) - plt.subplot(1,4,1) - plt.imshow(data.subset(channel=0).as_array()) - plt.title('Ground Truth') - plt.colorbar() - plt.subplot(1,4,2) - plt.imshow(noisy_data.subset(channel=0).as_array()) - plt.title('Noisy Data') - plt.colorbar() - plt.subplot(1,4,3) - plt.imshow(pdhg.get_output().subset(channel=0).as_array()) - plt.title('TV Reconstruction') - plt.colorbar() - plt.subplot(1,4,4) - plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') - plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction') - plt.legend() - plt.title('Middle Line Profiles') - plt.show() - - -##%% Check with CVX solution - -from ccpi.optimisation.operators import SparseFiniteDiff - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - ##Construct problem - u = Variable(ig.shape) - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - # Define Total Variation as a regulariser - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - fidelity = pnorm( u - noisy_data.as_array(),1) - - # choose solver - if 'MOSEK' in installed_solvers(): - solver = MOSEK - else: - solver = SCS - - obj = Minimize( regulariser + fidelity) - prob = Problem(obj) - result = prob.solve(verbose = True, solver = solver) - - diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value ) - - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(pdhg.get_output().as_array()) - plt.title('PDHG solution') - plt.colorbar() - plt.subplot(3,1,2) - plt.imshow(u.value) - plt.title('CVX solution') - plt.colorbar() - plt.subplot(3,1,3) - plt.imshow(diff_cvx) - plt.title('Difference') - plt.colorbar() - plt.show() - - plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') - plt.legend() - plt.title('Middle Line Profiles') - plt.show() - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_2D_time.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_2D_time.py deleted file mode 100644 index 14608db..0000000 --- a/Wrappers/Python/demos/PDHG_examples/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/PDHG_TV_Denoising_Gaussian.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Gaussian.py deleted file mode 100644 index 9d00ee1..0000000 --- a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Gaussian.py +++ /dev/null @@ -1,225 +0,0 @@ -#======================================================================== -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). - -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -#========================================================================= -""" - -Total Variation Denoising using PDHG algorithm: - - -Problem: min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2} - - \alpha: Regularization parameter - - \nabla: Gradient operator - - g: Noisy Data with Gaussian Noise - - Method = 0 ( PDHG - split ) : K = [ \nabla, - Identity] - - - Method = 1 (PDHG - explicit ): K = \nabla - -""" - -from ccpi.framework import ImageData, ImageGeometry - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction - -from ccpi.framework import TestData -import os, sys -loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) - -# Load Data -N = 200 -M = 300 - - -# user can change the size of the input data -# you can choose between -# TestData.PEPPERS 2D + Channel -# TestData.BOAT 2D -# TestData.CAMERA 2D -# TestData.RESOLUTION_CHART 2D -# TestData.SIMPLE_PHANTOM_2D 2D -data = loader.load(TestData.BOAT, size=(N,M), scale=(0,1)) - -ig = data.geometry -ag = ig - -# Create Noisy data. Add Gaussian noise -np.random.seed(10) -noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.1, size=data.shape) ) - -print ("min {} max {}".format(data.as_array().min(), data.as_array().max())) - -# Show Ground Truth and Noisy Data -plt.figure() -plt.subplot(1,3,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(1,3,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.subplot(1,3,3) -plt.imshow((data - noisy_data).as_array()) -plt.title('diff') -plt.colorbar() - -plt.show() - -# Regularisation Parameter -alpha = .1 - -method = '0' - -if method == '0': - - # Create operators - op1 = Gradient(ig, correlation=Gradient.CORRELATION_SPACE) - op2 = Identity(ig, ag) - - # Create BlockOperator - operator = BlockOperator(op1, op2, shape=(2,1) ) - - # Create functions - f1 = alpha * MixedL21Norm() - f2 = 0.5 * L2NormSquared(b = noisy_data) - f = BlockFunction(f1, f2) - - g = ZeroFunction() - -else: - - # Without the "Block Framework" - operator = Gradient(ig) - f = alpha * MixedL21Norm() - g = 0.5 * L2NormSquared(b = noisy_data) - -# Compute Operator Norm -normK = operator.norm() - -# Primal & Dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) - -# Setup and Run the PDHG algorithm -pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) -pdhg.max_iteration = 10000 -pdhg.update_objective_interval = 100 -pdhg.run(1000, verbose=True) - -# Show Results -plt.figure() -plt.subplot(1,3,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.clim(0,1) -plt.subplot(1,3,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.clim(0,1) -plt.subplot(1,3,3) -plt.imshow(pdhg.get_output().as_array()) -plt.title('TV Reconstruction') -plt.clim(0,1) -plt.colorbar() -plt.show() - -plt.plot(np.linspace(0,N,M), noisy_data.as_array()[int(N/2),:], label = 'Noisy data') -plt.plot(np.linspace(0,N,M), data.as_array()[int(N/2),:], label = 'GTruth') -plt.plot(np.linspace(0,N,M), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction') - -plt.legend() -plt.title('Middle Line Profiles') -plt.show() - - -#%% Check with CVX solution - -from ccpi.optimisation.operators import SparseFiniteDiff - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - ##Construct problem - u = Variable(ig.shape) - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - # Define Total Variation as a regulariser - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - fidelity = 0.5 * sum_squares(u - noisy_data.as_array()) - - # choose solver - if 'MOSEK' in installed_solvers(): - solver = MOSEK - else: - solver = SCS - - obj = Minimize( regulariser + fidelity) - prob = Problem(obj) - result = prob.solve(verbose = True, solver = MOSEK) - - diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value ) - - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(pdhg.get_output().as_array()) - plt.title('PDHG solution') - plt.colorbar() - plt.subplot(3,1,2) - plt.imshow(u.value) - plt.title('CVX solution') - plt.colorbar() - plt.subplot(3,1,3) - plt.imshow(diff_cvx) - plt.title('Difference') - plt.colorbar() - plt.show() - - plt.plot(np.linspace(0,N,M), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,M), u.value[int(N/2),:], label = 'CVX') - plt.plot(np.linspace(0,N,M), data.as_array()[int(N/2),:], label = 'Truth') - - plt.legend() - plt.title('Middle Line Profiles') - plt.show() - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) - diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Gaussian_3D.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Gaussian_3D.py deleted file mode 100644 index 03dc2ef..0000000 --- a/Wrappers/Python/demos/PDHG_examples/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/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Poisson.py deleted file mode 100644 index 1d887c1..0000000 --- a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Poisson.py +++ /dev/null @@ -1,212 +0,0 @@ -#======================================================================== -# Copyright 2019 Science Technology Facilities Council -# Copyright 2019 University of Manchester -# -# This work is part of the Core Imaging Library developed by Science Technology -# Facilities Council and University of Manchester -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0.txt -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -#========================================================================= - -""" - -Total Variation Denoising using PDHG algorithm: - -Problem: min_x, x>0 \alpha * ||\nabla x||_{2,1} + \int x - g * log(x) - - \alpha: Regularization parameter - - \nabla: Gradient operator - - g: Noisy Data with Poisson Noise - - - Method = 0 ( PDHG - split ) : K = [ \nabla, - Identity] - - - Method = 1 (PDHG - explicit ): K = \nabla - -""" - -from ccpi.framework import ImageData, ImageGeometry - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, IndicatorBox, \ - MixedL21Norm, BlockFunction - -from skimage.util import random_noise - -# Create phantom for TV Poisson denoising -N = 100 - -data = np.zeros((N,N)) -data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 -data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 -data = ImageData(data) -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -ag = ig - -# Create noisy data. Apply Poisson noise -n1 = random_noise(data.as_array(), mode = 'poisson', seed = 10) -noisy_data = ImageData(n1) - - -# Show Ground Truth and Noisy Data -plt.figure(figsize=(10,10)) -plt.subplot(2,1,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(2,1,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.show() - - -# Regularisation Parameter -alpha = 2 - -method = '0' - -if method == '0': - - # Create operators - op1 = Gradient(ig) - op2 = Identity(ig, ag) - - # Create BlockOperator - operator = BlockOperator(op1, op2, shape=(2,1) ) - - # Create functions - - f1 = alpha * MixedL21Norm() - f2 = KullbackLeibler(noisy_data) - f = BlockFunction(f1, f2) - g = IndicatorBox(lower=0) - -else: - - # Without the "Block Framework" - operator = Gradient(ig) - f = alpha * MixedL21Norm() - g = KullbackLeibler(noisy_data) - - -# Compute operator Norm -normK = operator.norm() - -# Primal & Dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) - -# Setup and Run the PDHG algorithm -pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) -pdhg.max_iteration = 3000 -pdhg.update_objective_interval = 200 -pdhg.run(3000, verbose=True) - - -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow(pdhg.get_output().as_array()) -plt.title('TV Reconstruction') -plt.colorbar() -plt.show() -## -plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') -plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction') -plt.legend() -plt.title('Middle Line Profiles') -plt.show() - - -#%% Check with CVX solution - -from ccpi.optimisation.operators import SparseFiniteDiff - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - ##Construct problem - u1 = Variable(ig.shape) - q = Variable() - w = Variable() - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - # Define Total Variation as a regulariser - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u1), DY.matrix() * vec(u1)]), 2, axis = 0)) - fidelity = sum(kl_div(noisy_data.as_array(), u1)) - - constraints = [q>=fidelity, u1>=0] - - solver = ECOS - obj = Minimize( regulariser + q) - prob = Problem(obj, constraints) - result = prob.solve(verbose = True, solver = solver) - - - diff_cvx = numpy.abs( pdhg.get_output().as_array() - u1.value ) - - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(pdhg.get_output().as_array()) - plt.title('PDHG solution') - plt.colorbar() - plt.subplot(3,1,2) - plt.imshow(u1.value) - plt.title('CVX solution') - plt.colorbar() - plt.subplot(3,1,3) - plt.imshow(diff_cvx) - plt.title('Difference') - plt.colorbar() - plt.show() - - plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), u1.value[int(N/2),:], label = 'CVX') - plt.legend() - plt.title('Middle Line Profiles') - plt.show() - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) - - - - - diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_SaltPepper.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_SaltPepper.py deleted file mode 100644 index c5709c3..0000000 --- a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_SaltPepper.py +++ /dev/null @@ -1,213 +0,0 @@ -#======================================================================== -# Copyright 2019 Science Technology Facilities Council -# Copyright 2019 University of Manchester -# -# This work is part of the Core Imaging Library developed by Science Technology -# Facilities Council and University of Manchester -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0.txt -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -#========================================================================= - -""" - -Total Variation Denoising using PDHG algorithm: - - -Problem: min_x, x>0 \alpha * ||\nabla x||_{2,1} + ||x-g||_{1} - - \alpha: Regularization parameter - - \nabla: Gradient operator - - g: Noisy Data with Salt & Pepper Noise - - - Method = 0 ( PDHG - split ) : K = [ \nabla, - Identity] - - - Method = 1 (PDHG - explicit ): K = \nabla - - -""" - -from ccpi.framework import ImageData, ImageGeometry - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, L1Norm, \ - MixedL21Norm, BlockFunction - -from skimage.util import random_noise - -# Create phantom for TV Salt & Pepper denoising -N = 100 - -data = np.zeros((N,N)) -data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 -data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 -data = ImageData(data) -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -ag = ig - -# Create noisy data. Apply Salt & Pepper noise -n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2) -noisy_data = ImageData(n1) - -# Show Ground Truth and Noisy Data -plt.figure(figsize=(15,15)) -plt.subplot(2,1,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(2,1,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.show() - -# Regularisation Parameter -alpha = 2 - -method = '0' - -if method == '0': - - # Create operators - op1 = Gradient(ig) - op2 = Identity(ig, ag) - - # Create BlockOperator - operator = BlockOperator(op1, op2, shape=(2,1) ) - - # Create functions - f1 = alpha * MixedL21Norm() - f2 = L1Norm(b = noisy_data) - f = BlockFunction(f1, f2) - - g = ZeroFunction() - -else: - - # Without the "Block Framework" - operator = Gradient(ig) - f = alpha * MixedL21Norm() - g = L1Norm(b = noisy_data) - - -# Compute operator Norm -normK = operator.norm() - -# Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) -opt = {'niter':2000, 'memopt': True} - -# Setup and run the PDHG algorithm -pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) -pdhg.max_iteration = 2000 -pdhg.update_objective_interval = 50 -pdhg.run(2000) - - -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow(pdhg.get_output().as_array()) -plt.title('TV Reconstruction') -plt.colorbar() -plt.show() -## -plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') -plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction') -plt.legend() -plt.title('Middle Line Profiles') -plt.show() - - -##%% Check with CVX solution - -from ccpi.optimisation.operators import SparseFiniteDiff - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - ##Construct problem - u = Variable(ig.shape) - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - # Define Total Variation as a regulariser - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - fidelity = pnorm( u - noisy_data.as_array(),1) - - # choose solver - if 'MOSEK' in installed_solvers(): - solver = MOSEK - else: - solver = SCS - - obj = Minimize( regulariser + fidelity) - prob = Problem(obj) - result = prob.solve(verbose = True, solver = solver) - - diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value ) - - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(pdhg.get_output().as_array()) - plt.title('PDHG solution') - plt.colorbar() - plt.subplot(3,1,2) - plt.imshow(u.value) - plt.title('CVX solution') - plt.colorbar() - plt.subplot(3,1,3) - plt.imshow(diff_cvx) - plt.title('Difference') - plt.colorbar() - plt.show() - - plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') - plt.legend() - plt.title('Middle Line Profiles') - plt.show() - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) - - - - - diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_gaussian.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_gaussian.py deleted file mode 100644 index 6acbfcc..0000000 --- a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_gaussian.py +++ /dev/null @@ -1,212 +0,0 @@ -#======================================================================== -# Copyright 2019 Science Technology Facilities Council -# Copyright 2019 University of Manchester -# -# This work is part of the Core Imaging Library developed by Science Technology -# Facilities Council and University of Manchester -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0.txt -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -#========================================================================= - -""" - -Total Variation 2D Tomography Reconstruction using PDHG algorithm: - - -Problem: min_x, x>0 \alpha * ||\nabla x||_{2,1} + \frac{1}{2}||Ax - g||^{2} - - \nabla: Gradient operator - - A: Projection Matrix - g: Noisy sinogram corrupted with Gaussian Noise - - \alpha: Regularization parameter - -""" - -from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG - -from ccpi.optimisation.operators import BlockOperator, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction - -from ccpi.astra.ops import AstraProjectorSimple - - - -# Create phantom for TV 2D tomography -N = 100 -x = np.zeros((N,N)) -x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 -x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 - -data = ImageData(x) -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) - -detectors = N -angles = np.linspace(0, np.pi, N) - -ag = AcquisitionGeometry('parallel','2D',angles, detectors) - -device = input('Available device: GPU==1 / CPU==0 ') - -if device=='1': - dev = 'gpu' -else: - dev = 'cpu' - -Aop = AstraProjectorSimple(ig, ag, dev) -sin = Aop.direct(data) - -# Create noisy data. Apply Poisson noise -n1 = np.random.normal(0, 3, size=ig.shape) -noisy_data = AcquisitionData(n1 + sin.as_array(), ag) - -# Show Ground Truth and Noisy Data -plt.figure(figsize=(10,10)) -plt.subplot(2,1,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(2,1,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.show() - -# Regularisation Parameter -alpha = 50 - -# Create operators -op1 = Gradient(ig) -op2 = Aop - -# Create BlockOperator -operator = BlockOperator(op1, op2, shape=(2,1) ) - -# Create functions - -f1 = alpha * MixedL21Norm() -f2 = 0.5 * L2NormSquared(b=noisy_data) -f = BlockFunction(f1, f2) - -g = ZeroFunction() - -# Compute operator Norm -normK = operator.norm() - -# Primal & dual stepsizes -sigma = 10 -tau = 1/(sigma*normK**2) - -# Setup and run the PDHG algorithm -pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) -pdhg.max_iteration = 2000 -pdhg.update_objective_interval = 200 -pdhg.run(2000) - -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow(pdhg.get_output().as_array()) -plt.title('TV Reconstruction') -plt.colorbar() -plt.show() -## -plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') -plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction') -plt.legend() -plt.title('Middle Line Profiles') -plt.show() - - -#%% Check with CVX solution - -from ccpi.optimisation.operators import SparseFiniteDiff -import astra -import numpy - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - -if cvx_not_installable: - - ##Construct problem - u = Variable(N*N) - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - - # create matrix representation for Astra operator - vol_geom = astra.create_vol_geom(N, N) - proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) - - proj_id = astra.create_projector('strip', proj_geom, vol_geom) - - matrix_id = astra.projector.matrix(proj_id) - - ProjMat = astra.matrix.get(matrix_id) - - tmp = noisy_data.as_array().ravel() - - fidelity = 0.5 * sum_squares(ProjMat * u - tmp) - - solver = MOSEK - obj = Minimize( regulariser + fidelity) - prob = Problem(obj) - result = prob.solve(verbose = True, solver = solver) - - diff_cvx = numpy.abs( pdhg.get_output().as_array() - np.reshape(u.value, (N,N) )) - - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(pdhg.get_output().as_array()) - plt.title('PDHG solution') - plt.colorbar() - plt.subplot(3,1,2) - plt.imshow(np.reshape(u.value, (N, N))) - plt.title('CVX solution') - plt.colorbar() - plt.subplot(3,1,3) - plt.imshow(diff_cvx) - plt.title('Difference') - plt.colorbar() - plt.show() - - plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), np.reshape(u.value, (N,N) )[int(N/2),:], label = 'CVX') - plt.legend() - plt.title('Middle Line Profiles') - plt.show() - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
\ No newline at end of file diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py deleted file mode 100644 index c44f393..0000000 --- a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py +++ /dev/null @@ -1,226 +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, AcquisitionData - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG - -from ccpi.optimisation.operators import BlockOperator, Gradient -from ccpi.optimisation.functions import KullbackLeibler, \ - MixedL21Norm, BlockFunction, IndicatorBox - -from ccpi.astra.ops import AstraProjectorSimple -from ccpi.framework import TestData -import os, sys - -""" - -Total Variation Denoising using PDHG algorithm: - - -Problem: min_x, x>0 \alpha * ||\nabla x||_{2,1} + int A x -g log(Ax + \eta) - - \nabla: Gradient operator - - A: Projection Matrix - g: Noisy sinogram corrupted with Poisson Noise - - \eta: Background Noise - \alpha: Regularization parameter - -""" - -loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) - -# Load Data -N = 50 -M = 50 -data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1)) - -ig = data.geometry -ag = ig - -#Create Acquisition Data and apply poisson noise - -detectors = N -angles = np.linspace(0, np.pi, N) - -ag = AcquisitionGeometry('parallel','2D',angles, detectors) - -device = input('Available device: GPU==1 / CPU==0 ') - -if device=='1': - dev = 'gpu' -else: - dev = 'cpu' - -Aop = AstraProjectorSimple(ig, ag, 'cpu') -sin = Aop.direct(data) - -# Create noisy data. Apply Poisson noise -scale = 0.5 -eta = 0 -n1 = scale * np.random.poisson(eta + sin.as_array()/scale) - -noisy_data = AcquisitionData(n1, ag) - -# Show Ground Truth and Noisy Data -plt.figure(figsize=(10,10)) -plt.subplot(2,1,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(2,1,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.show() - -# Regularisation Parameter -alpha = 2 - -# Create operators -op1 = Gradient(ig) -op2 = Aop - -# Create BlockOperator -operator = BlockOperator(op1,op2, shape=(2,1) ) - -# Create functions - -f1 = alpha * MixedL21Norm() -f2 = KullbackLeibler(noisy_data) -f = BlockFunction(f1, f2) - -g = IndicatorBox(lower=0) - - -# Compute operator Norm -normK = operator.norm() - -# Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) - - -# Setup and run the PDHG algorithm -pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) -pdhg.max_iteration = 3000 -pdhg.update_objective_interval = 500 -pdhg.run(3000, verbose = True) - -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow(pdhg.get_output().as_array()) -plt.title('TV Reconstruction') -plt.colorbar() -plt.show() -## -plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') -plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction') -plt.legend() -plt.title('Middle Line Profiles') -plt.show() - - -#%% Check with CVX solution -# -from ccpi.optimisation.operators import SparseFiniteDiff -import astra -import numpy - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - - ##Construct problem - u = Variable(N*N) - q = Variable() - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - - # create matrix representation for Astra operator - - vol_geom = astra.create_vol_geom(N, N) - proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) - - proj_id = astra.create_projector('strip', proj_geom, vol_geom) - - matrix_id = astra.projector.matrix(proj_id) - - ProjMat = astra.matrix.get(matrix_id) - - tmp = noisy_data.as_array().ravel() - - fidelity = sum(kl_div(tmp, ProjMat * u)) - - constraints = [q>=fidelity, u>=0] - solver = SCS - obj = Minimize( regulariser + q) - prob = Problem(obj, constraints) - result = prob.solve(verbose = True, solver = solver) - - diff_cvx = np.abs(pdhg.get_output().as_array() - np.reshape(u.value, (N, N))) - - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(pdhg.get_output().as_array()) - plt.title('PDHG solution') - plt.colorbar() - plt.subplot(3,1,2) - plt.imshow(np.reshape(u.value, (N, N))) - plt.title('CVX solution') - plt.colorbar() - plt.subplot(3,1,3) - plt.imshow(diff_cvx) - plt.title('Difference') - plt.colorbar() - plt.show() - - plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), np.reshape(u.value, (N, N))[int(N/2),:], label = 'CVX') - plt.legend() - plt.title('Middle Line Profiles') - plt.show() - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
\ No newline at end of file diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_Tikhonov_Denoising.py b/Wrappers/Python/demos/PDHG_examples/PDHG_Tikhonov_Denoising.py deleted file mode 100644 index f00f1cc..0000000 --- a/Wrappers/Python/demos/PDHG_examples/PDHG_Tikhonov_Denoising.py +++ /dev/null @@ -1,256 +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 Denoising using PDHG algorithm: - - -Problem: min_{x} \alpha * ||\nabla x||_{2}^{2} + \frac{1}{2} * || x - g ||_{2}^{2} - - \alpha: Regularization parameter - - \nabla: Gradient operator - - g: Noisy Data with Gaussian Noise - - Method = 0 ( PDHG - split ) : K = [ \nabla, - Identity] - - - Method = 1 (PDHG - explicit ): K = \nabla - -""" - -from ccpi.framework import ImageData, ImageGeometry, TestData - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared,\ - BlockFunction, KullbackLeibler, L1Norm - -import sys, os -if int(numpy.version.version.split('.')[1]) > 12: - from skimage.util import random_noise -else: - from demoutil import random_noise - - -# user supplied input -if len(sys.argv) > 1: - which_noise = int(sys.argv[1]) -else: - which_noise = 0 -print ("Applying {} noise") - -if len(sys.argv) > 2: - method = sys.argv[2] -else: - method = '0' -print ("method ", method) - -# Create phantom for TV Salt & Pepper denoising -N = 100 - -loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) -data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,N)) -ig = data.geometry -ag = ig - -# Create noisy data. Apply Salt & Pepper noise -# Create noisy data. -# Apply Salt & Pepper noise -# gaussian -# poisson -noises = ['gaussian', 'poisson', 's&p'] -noise = noises[which_noise] -if noise == 's&p': - n1 = random_noise(data.as_array(), mode = noise, salt_vs_pepper = 0.9, amount=0.2) -elif noise == 'poisson': - n1 = random_noise(data.as_array(), mode = noise, seed = 10) -elif noise == 'gaussian': - n1 = random_noise(data.as_array(), mode = noise, seed = 10) -else: - raise ValueError('Unsupported Noise ', noise) -noisy_data = ImageData(n1) - -# fidelity -if noise == 's&p': - f2 = L1Norm(b=noisy_data) -elif noise == 'poisson': - f2 = KullbackLeibler(noisy_data) -elif noise == 'gaussian': - f2 = 0.5 * L2NormSquared(b=noisy_data) - -# Show Ground Truth and Noisy Data -plt.figure(figsize=(10,5)) -plt.subplot(1,2,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(1,2,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.show() - - -# Regularisation Parameter -# no edge preservation alpha is big -if noise == 's&p': - alpha = 8. -elif noise == 'poisson': - alpha = 8. -elif noise == 'gaussian': - alpha = 8. - -if method == '0': - - # Create operators - op1 = Gradient(ig) - op2 = Identity(ig, ag) - - # Create BlockOperator - operator = BlockOperator(op1, op2, shape=(2,1) ) - - # Create functions - - f1 = alpha * L2NormSquared() - # f2 must change according to noise - #f2 = 0.5 * L2NormSquared(b = noisy_data) - f = BlockFunction(f1, f2) - g = ZeroFunction() - -else: - - # Without the "Block Framework" - operator = Gradient(ig) - f = alpha * L2NormSquared() - # g must change according to noise - #g = 0.5 * L2NormSquared(b = noisy_data) - g = f2 - - -# Compute operator Norm -normK = operator.norm() - -# Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) -opt = {'niter':2000, 'memopt': True} - -# Setup and run the PDHG algorithm -pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) -pdhg.max_iteration = 2000 -pdhg.update_objective_interval = 50 -pdhg.run(1500) - - -plt.figure(figsize=(20,5)) -plt.subplot(1,4,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(1,4,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.subplot(1,4,3) -plt.imshow(pdhg.get_output().as_array()) -plt.title('Tikhonov Reconstruction') -plt.colorbar() -plt.subplot(1,4,4) -plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') -plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction') -plt.legend() -plt.title('Middle Line Profiles') -plt.show() - - - -##%% Check with CVX solution - -from ccpi.optimisation.operators import SparseFiniteDiff - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - ##Construct problem - u = Variable(ig.shape) - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - # Define Total Variation as a regulariser - - regulariser = alpha * sum_squares(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - fidelity = 0.5 * sum_squares(u - noisy_data.as_array()) - - # choose solver - if 'MOSEK' in installed_solvers(): - solver = MOSEK - else: - solver = SCS - - obj = Minimize( regulariser + fidelity) - prob = Problem(obj) - result = prob.solve(verbose = True, solver = solver) - - diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value ) - - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(pdhg.get_output().as_array()) - plt.title('PDHG solution') - plt.colorbar() - plt.subplot(3,1,2) - plt.imshow(u.value) - plt.title('CVX solution') - plt.colorbar() - plt.subplot(3,1,3) - plt.imshow(diff_cvx) - plt.title('Difference') - plt.colorbar() - plt.show() - - plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') - plt.legend() - plt.title('Middle Line Profiles') - plt.show() - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) - - - - - diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_Tikhonov_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/PDHG_Tikhonov_Tomo2D.py deleted file mode 100644 index 02cd053..0000000 --- a/Wrappers/Python/demos/PDHG_examples/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() - - |