diff options
3 files changed, 221 insertions, 236 deletions
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py index 4d6da00..e9bad7e 100755 --- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py +++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py @@ -23,23 +23,21 @@ 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 - +Problem: min_{u} \alpha * ||\nabla u - w||_{2,1} + + \beta * || E u ||_{2,1} + + Fidelity(u, g) + \nabla: Gradient operator - E: Symmetrized Gradient operator + E: Symmetrized Gradient operator + \alpha: Regularization parameter + \beta: Regularization parameter - g: Noisy Data with Salt & Pepper Noise + g: Noisy Data + + Fidelity = 1) L2NormSquarred ( \frac{1}{2} * || u - g ||_{2}^{2} ) if Noise is Gaussian + 2) L1Norm ( ||u - g||_{1} )if Noise is Salt & Pepper + 3) Kullback Leibler (\int u - g * log(u) + Id_{u>0}) if Noise is Poisson + Method = 0 ( PDHG - split ) : K = [ \nabla, - Identity ZeroOperator, E @@ -48,10 +46,15 @@ Problem: min_{x} \alpha * ||\nabla x - w||_{2,1} + Method = 1 (PDHG - explicit ): K = [ \nabla, - Identity ZeroOperator, E ] + + Default: TGV denoising + noise = Gaussian + Fidelity = L2NormSquarred + method = 0 """ -from ccpi.framework import ImageData, ImageGeometry +from ccpi.framework import ImageData import numpy as np import numpy @@ -63,14 +66,14 @@ from ccpi.optimisation.operators import BlockOperator, Identity, \ Gradient, SymmetrizedGradient, ZeroOperator from ccpi.optimisation.functions import ZeroFunction, L1Norm, \ MixedL21Norm, BlockFunction, KullbackLeibler, L2NormSquared -import sys + +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]) @@ -81,35 +84,23 @@ print ("Applying {} noise") if len(sys.argv) > 2: method = sys.argv[2] else: - method = '1' + method = '0' 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()) +loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) +data = loader.load(TestData.SHAPES) +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, seed=10) elif noise == 'poisson': - n1 = random_noise(data.as_array(), mode = noise, seed = 10) + scale = 5 + n1 = random_noise( data.as_array()/scale, mode = noise, seed = 10)*scale elif noise == 'gaussian': n1 = random_noise(data.as_array(), mode = noise, seed = 10) else: @@ -128,23 +119,24 @@ plt.title('Noisy Data') plt.colorbar() plt.show() -# Regularisation Parameters +# Regularisation Parameter depending on the noise distribution if noise == 's&p': alpha = 0.8 elif noise == 'poisson': - alpha = .1 -elif noise == 'gaussian': alpha = .3 +elif noise == 'gaussian': + alpha = .2 -beta = numpy.sqrt(2)* alpha +# TODO add ref why this choice +beta = 2 * alpha -# fidelity +# 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) + f3 = 0.5 * L2NormSquared(b=noisy_data) if method == '0': @@ -162,7 +154,6 @@ if method == '0': f1 = alpha * MixedL21Norm() f2 = beta * MixedL21Norm() - # f3 depends on the noise characteristics f = BlockFunction(f1, f2, f3) g = ZeroFunction() @@ -183,28 +174,27 @@ else: f = BlockFunction(f1, f2) g = BlockFunction(f3, ZeroFunction()) -## Compute operator Norm +# 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) +pdhg.update_objective_interval = 100 +pdhg.run(2000) -#%% +# Show results plt.figure(figsize=(20,5)) plt.subplot(1,4,1) -plt.imshow(data.as_array()) +plt.imshow(data.subset(channel=0).as_array()) plt.title('Ground Truth') plt.colorbar() plt.subplot(1,4,2) -plt.imshow(noisy_data.as_array()) +plt.imshow(noisy_data.subset(channel=0).as_array()) plt.title('Noisy Data') plt.colorbar() plt.subplot(1,4,3) @@ -212,10 +202,8 @@ 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.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[int(ig.shape[0]/2),:], label = 'GTruth') +plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output()[0].as_array()[int(ig.shape[0]/2),:], label = 'TGV reconstruction') plt.legend() plt.title('Middle Line Profiles') -plt.show() - - +plt.show()
\ No newline at end of file diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py index 5aa334d..8a9920c 100644 --- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py +++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py @@ -18,26 +18,37 @@ # limitations under the License. # #========================================================================= + """ -Tikhonov Denoising using PDHG algorithm: +``Tikhonov`` Regularization Denoising using PDHG algorithm: -Problem: min_{x} \alpha * ||\nabla x||_{2}^{2} + \frac{1}{2} * || x - g ||_{2}^{2} +Problem: min_{u}, \alpha * ||\nabla u||_{2,2} + Fidelity(u, g) \alpha: Regularization parameter \nabla: Gradient operator - g: Noisy Data with Gaussian Noise + g: Noisy Data + Fidelity = 1) L2NormSquarred ( \frac{1}{2} * || u - g ||_{2}^{2} ) if Noise is Gaussian + 2) L1Norm ( ||u - g||_{1} )if Noise is Salt & Pepper + 3) Kullback Leibler (\int u - g * log(u) + Id_{u>0}) if Noise is Poisson + Method = 0 ( PDHG - split ) : K = [ \nabla, Identity] - Method = 1 (PDHG - explicit ): K = \nabla - -""" + Method = 1 (PDHG - explicit ): K = \nabla + + + Default: Tikhonov denoising + noise = Gaussian + Fidelity = L2NormSquarred + method = 0 + +""" from ccpi.framework import ImageData, TestData @@ -62,7 +73,7 @@ else: if len(sys.argv) > 1: which_noise = int(sys.argv[1]) else: - which_noise = 0 + which_noise = 2 print ("Applying {} noise") if len(sys.argv) > 2: @@ -71,39 +82,25 @@ 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.SHAPES) 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) + scale = 5 + n1 = random_noise( data.as_array()/scale, mode = noise, seed = 10)*scale elif noise == 'gaussian': n1 = random_noise(data.as_array(), mode = noise, seed = 10) else: raise ValueError('Unsupported Noise ', noise) noisy_data = ImageData(n1) -# 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) @@ -116,15 +113,21 @@ plt.title('Noisy Data') plt.colorbar() plt.show() - -# Regularisation Parameter -# no edge preservation alpha is big +# Regularisation Parameter depending on the noise distribution +if noise == 's&p': + alpha = 20 +elif noise == 'poisson': + alpha = 10 +elif noise == 'gaussian': + alpha = 5 + +# fidelity if noise == 's&p': - alpha = 8. + f2 = L1Norm(b=noisy_data) elif noise == 'poisson': - alpha = 8. + f2 = KullbackLeibler(noisy_data) elif noise == 'gaussian': - alpha = 8. + f2 = 0.5 * L2NormSquared(b=noisy_data) if method == '0': @@ -135,21 +138,14 @@ if method == '0': # Create BlockOperator operator = BlockOperator(op1, op2, shape=(2,1) ) - # Create functions - + # 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 @@ -159,13 +155,12 @@ 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 = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) pdhg.max_iteration = 2000 -pdhg.update_objective_interval = 50 -pdhg.run(1500) +pdhg.update_objective_interval = 100 +pdhg.run(2000) plt.figure(figsize=(20,5)) @@ -179,78 +174,78 @@ plt.title('Noisy Data') plt.colorbar() plt.subplot(1,4,3) plt.imshow(pdhg.get_output().as_array()) -plt.title('Tikhonov Reconstruction') +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.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[int(ig.shape[0]/2),:], label = 'GTruth') +plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(ig.shape[0]/2),:], label = 'Tikhonov reconstruction') plt.legend() plt.title('Middle Line Profiles') plt.show() -##%% Check with CVX solution - -from ccpi.optimisation.operators import SparseFiniteDiff - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - ##Construct problem - u = Variable(ig.shape) - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - # Define Total Variation as a regulariser - - regulariser = alpha * sum_squares(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - 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])) - - - - - +###%% 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/Tomo/PDHG_TV_Tomo2D_poisson.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_poisson.py index 95fe2c1..e7e33cb 100644 --- a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_poisson.py +++ b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_poisson.py @@ -84,7 +84,7 @@ 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) +n1 = np.random.poisson(eta + sin.as_array()) noisy_data = AcquisitionData(n1, ag) @@ -100,6 +100,8 @@ plt.title('Noisy Data') plt.colorbar() plt.show() + +#%% # Regularisation Parameter alpha = 2 @@ -157,72 +159,72 @@ 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 +#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 |