diff options
| -rw-r--r-- | Wrappers/Python/wip/fista_TV_denoising.py | 72 | ||||
| -rw-r--r-- | Wrappers/Python/wip/pdhg_TGV_denoising.py | 230 | ||||
| -rw-r--r-- | Wrappers/Python/wip/pdhg_TGV_tomography2D.py | 200 | ||||
| -rwxr-xr-x | Wrappers/Python/wip/pdhg_TV_denoising.py | 204 | ||||
| -rw-r--r-- | Wrappers/Python/wip/pdhg_TV_denoising3D.py | 360 | ||||
| -rw-r--r-- | Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py | 268 | ||||
| -rw-r--r-- | Wrappers/Python/wip/pdhg_TV_tomography2D.py | 231 | ||||
| -rw-r--r-- | Wrappers/Python/wip/pdhg_TV_tomography2D_time.py | 152 | ||||
| -rw-r--r-- | Wrappers/Python/wip/pdhg_tv_denoising_poisson.py | 187 |
9 files changed, 0 insertions, 1904 deletions
diff --git a/Wrappers/Python/wip/fista_TV_denoising.py b/Wrappers/Python/wip/fista_TV_denoising.py deleted file mode 100644 index a9712fc..0000000 --- a/Wrappers/Python/wip/fista_TV_denoising.py +++ /dev/null @@ -1,72 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer - -import numpy as np -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old, FISTA - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, FunctionOperatorComposition, BlockFunction, ScaledFunction - -from ccpi.optimisation.algs import FISTA - -from skimage.util import random_noise - -from timeit import default_timer as timer -def dt(steps): - return steps[-1] - steps[-2] - -# Create phantom for TV 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 - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -ag = ig - -# Create noisy data. Add Gaussian noise -n1 = random_noise(data, mode = 'gaussian', mean=0, var = 0.05, seed=10) -noisy_data = ImageData(n1) - - -plt.imshow(noisy_data.as_array()) -plt.title('Noisy data') -plt.show() - -# Regularisation Parameter -alpha = 2 - -operator = Gradient(ig) -g = alpha * MixedL21Norm() -f = 0.5 * L2NormSquared(b = noisy_data) - -x_init = ig.allocate() -opt = {'niter':2000} - - -x = FISTA(x_init, f, g, opt) - -#fista = FISTA() -#fista.set_up(x_init, f, g, opt ) -#fista.max_iteration = 10 -# -#fista.run(2000) -#plt.figure(figsize=(15,15)) -#plt.subplot(3,1,1) -#plt.imshow(fista.get_output().as_array()) -#plt.title('no memopt class') - - - diff --git a/Wrappers/Python/wip/pdhg_TGV_denoising.py b/Wrappers/Python/wip/pdhg_TGV_denoising.py deleted file mode 100644 index 1c570cb..0000000 --- a/Wrappers/Python/wip/pdhg_TGV_denoising.py +++ /dev/null @@ -1,230 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, \ - Gradient, SymmetrizedGradient, ZeroOperator -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction - -from skimage.util import random_noise - -from timeit import default_timer as timer -#def dt(steps): -# return steps[-1] - steps[-2] - -# Create phantom for TGV Gaussian 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 = data/data.max() - -plt.imshow(data) -plt.show() - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -ag = ig - -# Create noisy data. Add Gaussian noise -n1 = random_noise(data, mode = 'gaussian', mean=0, var = 0.005, seed=10) -noisy_data = ImageData(n1) - - -plt.imshow(noisy_data.as_array()) -plt.title('Noisy data') -plt.show() - -alpha, beta = 0.1, 0.5 - -#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") -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 = 0.5 * L2NormSquared(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(0.5 * L2NormSquared(b = noisy_data), ZeroFunction()) - -## Compute operator Norm -normK = operator.norm() -# -## Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) -## -opt = {'niter':500} -opt1 = {'niter':500, 'memopt': True} -# -t1 = timer() -res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -t2 = timer() -# -t3 = timer() -res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -t4 = timer() -# -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(res[0].as_array()) -plt.title('no memopt') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(res1[0].as_array()) -plt.title('memopt') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow((res1[0] - res[0]).abs().as_array()) -plt.title('diff') -plt.colorbar() -plt.show() - -print("NoMemopt/Memopt is {}/{}".format(t2-t1, t4-t3)) - - -###### - -#%% - -#def update_plot(it_update, objective, x): -# -## sol = pdhg.get_output() -# plt.figure() -# plt.imshow(x[0].as_array()) -# plt.show() -# -# -##def stop_criterion(x,) -# -#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) -#pdhg.max_iteration = 2000 -#pdhg.update_objective_interval = 100 -# -#pdhg.run(4000, verbose=False, callback=update_plot) - - -#%% - - - - - - - - -#%% 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 = 0.5 * sum_squares(u - noisy_data.as_array()) -# solver = MOSEK -# -# obj = Minimize( regulariser + fidelity) -# prob = Problem(obj) -# result = prob.solve(verbose = True, solver = solver) -# -# diff_cvx = numpy.abs( res[0].as_array() - u.value ) -# -# # Show result -# plt.figure(figsize=(15,15)) -# plt.subplot(3,1,1) -# plt.imshow(res[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), res[0].as_array()[int(N/2),:], label = 'PDHG') -# plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') -# plt.legend() -# -# print('Primal Objective (CVX) {} '.format(obj.value)) -# print('Primal Objective (PDHG) {} '.format(primal[-1])) -# print('Min/Max of absolute difference {}/{}'.format(diff_cvx.min(), diff_cvx.max())) -# - - diff --git a/Wrappers/Python/wip/pdhg_TGV_tomography2D.py b/Wrappers/Python/wip/pdhg_TGV_tomography2D.py deleted file mode 100644 index ee3b089..0000000 --- a/Wrappers/Python/wip/pdhg_TGV_tomography2D.py +++ /dev/null @@ -1,200 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, \ - Gradient, SymmetrizedGradient, ZeroOperator -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction - -from skimage.util import random_noise - -from timeit import default_timer as timer -from ccpi.astra.ops import AstraProjectorSimple - -#def dt(steps): -# return steps[-1] - steps[-2] - -# Create phantom for TGV Gaussian 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()) - -plt.imshow(data.as_array()) -plt.show() - -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) - -plt.imshow(sin.as_array()) -plt.title('Sinogram') -plt.colorbar() -plt.show() - -# Add Gaussian noise to the sinogram data -np.random.seed(10) -n1 = np.random.random(sin.shape) - -noisy_data = sin + ImageData(5*n1) - -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Sinogram') -plt.colorbar() -plt.show() - -#%% - -alpha, beta = 20, 50 - - -# 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 = 0.5 * L2NormSquared(b = noisy_data) -f = BlockFunction(f1, f2, f3) -g = ZeroFunction() - -## Compute operator Norm -normK = operator.norm() -# -## Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) -## -opt = {'niter':5000} -opt1 = {'niter':5000, 'memopt': True} -# -t1 = timer() -res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -t2 = timer() -# -plt.imshow(res[0].as_array()) -plt.show() - - -#t3 = timer() -#res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -#t4 = timer() -# -#plt.figure(figsize=(15,15)) -#plt.subplot(3,1,1) -#plt.imshow(res[0].as_array()) -#plt.title('no memopt') -#plt.colorbar() -#plt.subplot(3,1,2) -#plt.imshow(res1[0].as_array()) -#plt.title('memopt') -#plt.colorbar() -#plt.subplot(3,1,3) -#plt.imshow((res1[0] - res[0]).abs().as_array()) -#plt.title('diff') -#plt.colorbar() -#plt.show() -# -#print("NoMemopt/Memopt is {}/{}".format(t2-t1, t4-t3)) - - -#%% 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 = 0.5 * sum_squares(u - noisy_data.as_array()) -# solver = MOSEK -# -# obj = Minimize( regulariser + fidelity) -# prob = Problem(obj) -# result = prob.solve(verbose = True, solver = solver) -# -# diff_cvx = numpy.abs( res[0].as_array() - u.value ) -# -# # Show result -# plt.figure(figsize=(15,15)) -# plt.subplot(3,1,1) -# plt.imshow(res[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), res[0].as_array()[int(N/2),:], label = 'PDHG') -# plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') -# plt.legend() -# -# print('Primal Objective (CVX) {} '.format(obj.value)) -# print('Primal Objective (PDHG) {} '.format(primal[-1])) -# print('Min/Max of absolute difference {}/{}'.format(diff_cvx.min(), diff_cvx.max())) - - - diff --git a/Wrappers/Python/wip/pdhg_TV_denoising.py b/Wrappers/Python/wip/pdhg_TV_denoising.py deleted file mode 100755 index b16e8f9..0000000 --- a/Wrappers/Python/wip/pdhg_TV_denoising.py +++ /dev/null @@ -1,204 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction - -from skimage.util import random_noise - -from timeit import default_timer as timer -#def dt(steps): -# return steps[-1] - steps[-2] - -# Create phantom for TV Gaussian 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 - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -ag = ig - -# Create noisy data. Add Gaussian noise -n1 = random_noise(data, mode = 'gaussian', mean=0, var = 0.05, seed=10) -noisy_data = ImageData(n1) - - -plt.imshow(noisy_data.as_array()) -plt.title('Noisy data') -plt.show() - -# Regularisation Parameter -alpha = 0.5 - -#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") -method = '1' - -if method == '0': - - # Create operators - op1 = Gradient(ig) - op2 = Identity(ig, ag) - - # Form Composite Operator - 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: - - ########################################################################### - # No Composite # - ########################################################################### - operator = Gradient(ig) - f = alpha * MixedL21Norm() - g = 0.5 * L2NormSquared(b = noisy_data) - - -diag_precon = False - -if diag_precon: - - def tau_sigma_precond(operator): - - tau = 1/operator.sum_abs_row() - sigma = 1/ operator.sum_abs_col() - - return tau, sigma - - tau, sigma = tau_sigma_precond(operator) - -else: - # Compute operator Norm - normK = operator.norm() - - # Primal & dual stepsizes - sigma = 1 - tau = 1/(sigma*normK**2) - - -opt = {'niter':5000} -opt1 = {'niter':5000, 'memopt': True} - -t1 = timer() -res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -t2 = timer() - -t3 = timer() -res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -t4 = timer() - -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(res.as_array()) -plt.title('no memopt') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(res1.as_array()) -plt.title('memopt') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow((res1 - res).abs().as_array()) -plt.title('diff') -plt.colorbar() -plt.show() -# -plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt') -plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt') -plt.legend() -plt.show() -# -print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(t2-t1, t4 -t3)) -diff = (res1 - res).abs().as_array().max() -# -print(" Max of abs difference is {}".format(diff)) - - -#%% 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 = solver) - - diff_cvx = numpy.abs( res.as_array() - u.value ) - - # Show result - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(res.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), res1.as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') - plt.legend() - - - - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(primal[-1])) - - - - diff --git a/Wrappers/Python/wip/pdhg_TV_denoising3D.py b/Wrappers/Python/wip/pdhg_TV_denoising3D.py deleted file mode 100644 index 06ecfa2..0000000 --- a/Wrappers/Python/wip/pdhg_TV_denoising3D.py +++ /dev/null @@ -1,360 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer - -import numpy as np -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, FunctionOperatorComposition, BlockFunction - -from skimage.util import random_noise - -from timeit import default_timer as timer -def dt(steps): - return steps[-1] - steps[-2] - -#%% - -# Create phantom for TV 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_size = 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_size x N_size x N_size phantom (3D) -phantom_tm = TomoP3D.Model(model, N_size, path_library3D) -#toc=timeit.default_timer() -#Run_time = toc - tic -#print("Phantom has been built in {} seconds".format(Run_time)) -# -#sliceSel = int(0.5*N_size) -##plt.gray() -#plt.figure() -#plt.subplot(131) -#plt.imshow(phantom_tm[sliceSel,:,:],vmin=0, vmax=1) -#plt.title('3D Phantom, axial view') -# -#plt.subplot(132) -#plt.imshow(phantom_tm[:,sliceSel,:],vmin=0, vmax=1) -#plt.title('3D Phantom, coronal view') -# -#plt.subplot(133) -#plt.imshow(phantom_tm[:,:,sliceSel],vmin=0, vmax=1) -#plt.title('3D Phantom, sagittal view') -#plt.show() - -#%% - -N = N_size -ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N, voxel_num_z=N) - -n1 = random_noise(phantom_tm, mode = 'gaussian', mean=0, var = 0.001, seed=10) -noisy_data = ImageData(n1) -#plt.imshow(noisy_data.as_array()[:,:,32]) - -#%% - -# Regularisation Parameter -alpha = 0.02 - -#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") -method = '0' - -if method == '0': - - # Create operators - op1 = Gradient(ig) - op2 = Identity(ig) - - # Form Composite Operator - 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: - - ########################################################################### - # No Composite # - ########################################################################### - operator = Gradient(ig) - f = alpha * FunctionOperatorComposition(operator, MixedL21Norm()) - g = L2NormSquared(b = noisy_data) - - ########################################################################### -#%% - -# Compute operator Norm -normK = operator.norm() - -# Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) - -opt = {'niter':2000} -opt1 = {'niter':2000, 'memopt': True} - -#t1 = timer() -#res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -#t2 = timer() - - -t3 = timer() -res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -t4 = timer() - -#import cProfile -#cProfile.run('res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) ') -### -print ("No memopt in {}s, memopt in {}s ".format(t2-t1, t4 -t3)) -# -## -##%% -# -#plt.figure(figsize=(10,10)) -#plt.subplot(311) -#plt.imshow(res1.as_array()[sliceSel,:,:]) -#plt.colorbar() -#plt.title('3D Phantom, axial view') -# -#plt.subplot(312) -#plt.imshow(res1.as_array()[:,sliceSel,:]) -#plt.colorbar() -#plt.title('3D Phantom, coronal view') -# -#plt.subplot(313) -#plt.imshow(res1.as_array()[:,:,sliceSel]) -#plt.colorbar() -#plt.title('3D Phantom, sagittal view') -#plt.show() -# -#plt.figure(figsize=(10,10)) -#plt.subplot(311) -#plt.imshow(res.as_array()[sliceSel,:,:]) -#plt.colorbar() -#plt.title('3D Phantom, axial view') -# -#plt.subplot(312) -#plt.imshow(res.as_array()[:,sliceSel,:]) -#plt.colorbar() -#plt.title('3D Phantom, coronal view') -# -#plt.subplot(313) -#plt.imshow(res.as_array()[:,:,sliceSel]) -#plt.colorbar() -#plt.title('3D Phantom, sagittal view') -#plt.show() -# -#diff = (res1 - res).abs() -# -#plt.figure(figsize=(10,10)) -#plt.subplot(311) -#plt.imshow(diff.as_array()[sliceSel,:,:]) -#plt.colorbar() -#plt.title('3D Phantom, axial view') -# -#plt.subplot(312) -#plt.imshow(diff.as_array()[:,sliceSel,:]) -#plt.colorbar() -#plt.title('3D Phantom, coronal view') -# -#plt.subplot(313) -#plt.imshow(diff.as_array()[:,:,sliceSel]) -#plt.colorbar() -#plt.title('3D Phantom, sagittal view') -#plt.show() -# -# -# -# -##%% -#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) -#pdhg.max_iteration = 2000 -#pdhg.update_objective_interval = 100 -#### -#pdhgo = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) -#pdhgo.max_iteration = 2000 -#pdhgo.update_objective_interval = 100 -#### -#steps = [timer()] -#pdhgo.run(2000) -#steps.append(timer()) -#t1 = dt(steps) -## -#pdhg.run(2000) -#steps.append(timer()) -#t2 = dt(steps) -# -#print ("Time difference {}s {}s {}s Speedup {:.2f}".format(t1,t2,t2-t1, t2/t1)) -#res = pdhg.get_output() -#res1 = pdhgo.get_output() - -#%% -#plt.figure(figsize=(15,15)) -#plt.subplot(3,1,1) -#plt.imshow(res.as_array()) -#plt.title('no memopt') -#plt.colorbar() -#plt.subplot(3,1,2) -#plt.imshow(res1.as_array()) -#plt.title('memopt') -#plt.colorbar() -#plt.subplot(3,1,3) -#plt.imshow((res1 - res).abs().as_array()) -#plt.title('diff') -#plt.colorbar() -#plt.show() - - -#plt.figure(figsize=(15,15)) -#plt.subplot(3,1,1) -#plt.imshow(pdhg.get_output().as_array()) -#plt.title('no memopt class') -#plt.colorbar() -#plt.subplot(3,1,2) -#plt.imshow(res.as_array()) -#plt.title('no memopt') -#plt.colorbar() -#plt.subplot(3,1,3) -#plt.imshow((pdhg.get_output() - res).abs().as_array()) -#plt.title('diff') -#plt.colorbar() -#plt.show() -# -# -# -#plt.figure(figsize=(15,15)) -#plt.subplot(3,1,1) -#plt.imshow(pdhgo.get_output().as_array()) -#plt.title('memopt class') -#plt.colorbar() -#plt.subplot(3,1,2) -#plt.imshow(res1.as_array()) -#plt.title('no memopt') -#plt.colorbar() -#plt.subplot(3,1,3) -#plt.imshow((pdhgo.get_output() - res1).abs().as_array()) -#plt.title('diff') -#plt.colorbar() -#plt.show() - - - - - -# print ("Time difference {}s {}s {}s Speedup {:.2f}".format(t1,t2,t2-t1, t2/t1)) -# res = pdhg.get_output() -# res1 = pdhgo.get_output() -# -# diff = (res-res1) -# print ("diff norm {} max {}".format(diff.norm(), diff.abs().as_array().max())) -# print ("Sum ( abs(diff) ) {}".format(diff.abs().sum())) -# -# -# plt.figure(figsize=(5,5)) -# plt.subplot(1,3,1) -# plt.imshow(res.as_array()) -# plt.colorbar() -# #plt.show() -# -# #plt.figure(figsize=(5,5)) -# plt.subplot(1,3,2) -# plt.imshow(res1.as_array()) -# plt.colorbar() - -#plt.show() - - - -#======= -## opt = {'niter':2000, 'memopt': True} -# -## res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -# -#>>>>>>> origin/pdhg_fix -# -# -## opt = {'niter':2000, 'memopt': False} -## res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -# -## plt.figure(figsize=(5,5)) -## plt.subplot(1,3,1) -## plt.imshow(res.as_array()) -## plt.title('memopt') -## plt.colorbar() -## plt.subplot(1,3,2) -## plt.imshow(res1.as_array()) -## plt.title('no memopt') -## plt.colorbar() -## plt.subplot(1,3,3) -## plt.imshow((res1 - res).abs().as_array()) -## plt.title('diff') -## plt.colorbar() -## plt.show() -#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) -#pdhg.max_iteration = 2000 -#pdhg.update_objective_interval = 100 -# -# -#pdhgo = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) -#pdhgo.max_iteration = 2000 -#pdhgo.update_objective_interval = 100 -# -#steps = [timer()] -#pdhgo.run(200) -#steps.append(timer()) -#t1 = dt(steps) -# -#pdhg.run(200) -#steps.append(timer()) -#t2 = dt(steps) -# -#print ("Time difference {} {} {}".format(t1,t2,t2-t1)) -#sol = pdhg.get_output().as_array() -##sol = result.as_array() -## -#fig = plt.figure() -#plt.subplot(1,3,1) -#plt.imshow(noisy_data.as_array()) -#plt.colorbar() -#plt.subplot(1,3,2) -#plt.imshow(sol) -#plt.colorbar() -#plt.subplot(1,3,3) -#plt.imshow(pdhgo.get_output().as_array()) -#plt.colorbar() -# -#plt.show() -### -## -#### -##plt.plot(np.linspace(0,N,N), data[int(N/2),:], label = 'GTruth') -##plt.plot(np.linspace(0,N,N), sol[int(N/2),:], label = 'Recon') -##plt.legend() -##plt.show() -# -# -##%% -## diff --git a/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py b/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py deleted file mode 100644 index bd330fc..0000000 --- a/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py +++ /dev/null @@ -1,268 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, L1Norm, \ - MixedL21Norm, FunctionOperatorComposition, BlockFunction, ScaledFunction - -from skimage.util import random_noise - -from timeit import default_timer as timer - - - -# ############################################################################ -# Create phantom for TV 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 - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -ag = ig - -# Create noisy data. Add Gaussian noise -n1 = random_noise(data, mode = 's&p', salt_vs_pepper = 0.9, amount=0.2) -noisy_data = ImageData(n1) - -plt.imshow(noisy_data.as_array()) -plt.colorbar() -plt.show() - -#%% - -# Regularisation Parameter -alpha = 2 - -#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") -method = '0' -if method == '0': - - # Create operators - op1 = Gradient(ig) - op2 = Identity(ig, ag) - - operator = BlockOperator(op1, op2, shape=(2,1) ) - - f1 = alpha * MixedL21Norm() - f2 = L1Norm(b = noisy_data) - - f = BlockFunction(f1, f2 ) - g = ZeroFunction() - -else: - - ########################################################################### - # No Composite # - ########################################################################### - operator = Gradient(ig) - f = alpha * MixedL21Norm() - g = L1Norm(b = noisy_data) - ########################################################################### -#%% - -diag_precon = False - -if diag_precon: - - def tau_sigma_precond(operator): - - tau = 1/operator.sum_abs_row() - sigma = 1/ operator.sum_abs_col() - - return tau, sigma - - tau, sigma = tau_sigma_precond(operator) - -else: - # Compute operator Norm - normK = operator.norm() - - # Primal & dual stepsizes - sigma = 1 - tau = 1/(sigma*normK**2) - - -opt = {'niter':5000} -opt1 = {'niter':5000, 'memopt': True} - -t1 = timer() -res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -t2 = timer() - - -t3 = timer() -res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -t4 = timer() - -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(res.as_array()) -plt.title('no memopt') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(res1.as_array()) -plt.title('memopt') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow((res1 - res).abs().as_array()) -plt.title('diff') -plt.colorbar() -plt.show() -# -plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt') -plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt') -plt.legend() -plt.show() -# -print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(t2-t1, t4 -t3)) -diff = (res1 - res).abs().as_array().max() -# -print(" Max of abs difference is {}".format(diff)) - -#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) -#pdhg.max_iteration = 2000 -#pdhg.update_objective_interval = 10 -# -#pdhg.run(2000) - - - -#sol = pdhg.get_output().as_array() -##sol = result.as_array() -## -#fig = plt.figure() -#plt.subplot(1,2,1) -#plt.imshow(noisy_data.as_array()) -##plt.colorbar() -#plt.subplot(1,2,2) -#plt.imshow(sol) -##plt.colorbar() -#plt.show() -## - -## -#plt.plot(np.linspace(0,N,N), data[int(N/2),:], label = 'GTruth') -#plt.plot(np.linspace(0,N,N), sol[int(N/2),:], label = 'Recon') -#plt.legend() -#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) - - 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( res.as_array() - u.value ) - -# Show result - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(res.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), res1.as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') - plt.legend() - - - - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(primal[-1])) - - - -#try_cvx = input("Do you want CVX comparison (0/1)") -# -#if try_cvx=='0': -# -# from cvxpy import * -# import sys -# sys.path.insert(0,'/Users/evangelos/Desktop/Projects/CCPi/CCPi-Framework/Wrappers/Python/ccpi/optimisation/cvx_scripts') -# from cvx_functions import TV_cvx -# -# u = Variable((N, N)) -# fidelity = pnorm( u - noisy_data.as_array(),1) -# regulariser = alpha * TV_cvx(u) -# solver = MOSEK -# obj = Minimize( regulariser + fidelity) -# constraints = [] -# prob = Problem(obj, constraints) -# -# # Choose solver (SCS is fast but less accurate than MOSEK) -# result = prob.solve(verbose = True, solver = solver) -# -# print('Objective value is {} '.format(obj.value)) -# -# diff_pdhg_cvx = np.abs(u.value - res.as_array()) -# plt.imshow(diff_pdhg_cvx) -# plt.colorbar() -# plt.title('|CVX-PDHG|') -# plt.show() -# -# plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') -# plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'PDHG') -# plt.legend() -# plt.show() -# -#else: -# print('No CVX solution available') - - - - diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D.py b/Wrappers/Python/wip/pdhg_TV_tomography2D.py deleted file mode 100644 index e123739..0000000 --- a/Wrappers/Python/wip/pdhg_TV_tomography2D.py +++ /dev/null @@ -1,231 +0,0 @@ -# -*- coding: utf-8 -*- - -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData - -import numpy as np -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction - -from ccpi.astra.operators import AstraProjectorSimple -from timeit import default_timer as timer - - -#%%############################################################################### -# Create phantom for TV tomography - -#import os -#import tomophantom -#from tomophantom import TomoP2D -#from tomophantom.supp.qualitymetrics import QualityTools - -#model = 1 # select a model number from the library -#N = 150 # set dimension of the phantom -## one can specify an exact path to the parameters file -## path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' -#path = os.path.dirname(tomophantom.__file__) -#path_library2D = os.path.join(path, "Phantom2DLibrary.dat") -##This will generate a N_size x N_size phantom (2D) -#phantom_2D = TomoP2D.Model(model, N, path_library2D) -#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -#data = ImageData(phantom_2D, geometry=ig) - -N = 75 -x = np.zeros((N,N)) -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) - -plt.imshow(sin.as_array()) -plt.title('Sinogram') -plt.colorbar() -plt.show() - -# Add Gaussian noise to the sinogram data -np.random.seed(10) -n1 = np.random.random(sin.shape) - -noisy_data = sin + ImageData(5*n1) - -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Sinogram') -plt.colorbar() -plt.show() - -# Create operators -op1 = Gradient(ig) -op2 = Aop - -# Form Composite Operator -operator = BlockOperator(op1, op2, shape=(2,1) ) - -alpha = 10 -f = BlockFunction( alpha * MixedL21Norm(), \ - 0.5 * L2NormSquared(b = noisy_data) ) -g = ZeroFunction() - -# Compute operator Norm -normK = operator.norm() - -## Primal & dual stepsizes -diag_precon = True - -if diag_precon: - - def tau_sigma_precond(operator): - - tau = 1/operator.sum_abs_row() - sigma = 1/ operator.sum_abs_col() - - return tau, sigma - - tau, sigma = tau_sigma_precond(operator) - -else: - normK = operator.norm() - sigma = 1 - tau = 1/(sigma*normK**2) - -# Compute operator Norm - - -# Primal & dual stepsizes - -opt = {'niter':200} -opt1 = {'niter':200, 'memopt': True} - -t1 = timer() -res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -t2 = timer() - - -t3 = timer() -res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -t4 = timer() -# -# -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(res.as_array()) -plt.title('no memopt') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(res1.as_array()) -plt.title('memopt') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow((res1 - res).abs().as_array()) -plt.title('diff') -plt.colorbar() -plt.show() -# -plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt') -plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt') -plt.legend() -plt.show() -# -print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(t2-t1, t4 -t3)) -diff = (res1 - res).abs().as_array().max() -# -print(" Max of abs difference is {}".format(diff)) - - -#%% 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: -# -# -# 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) -# -# fidelity = 0.5 * sum_squares(ProjMat * u - noisy_data.as_array().ravel()) -# -# # 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) -# -##%% -# -# u_rs = np.reshape(u.value, (N,N)) -# -# diff_cvx = numpy.abs( res.as_array() - u_rs ) -# -# # Show result -# plt.figure(figsize=(15,15)) -# plt.subplot(3,1,1) -# plt.imshow(res.as_array()) -# plt.title('PDHG solution') -# plt.colorbar() -# -# plt.subplot(3,1,2) -# plt.imshow(u_rs) -# 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), res1.as_array()[int(N/2),:], label = 'PDHG') -# plt.plot(np.linspace(0,N,N), u_rs[int(N/2),:], label = 'CVX') -# plt.legend() -# -# -# print('Primal Objective (CVX) {} '.format(obj.value)) -# print('Primal Objective (PDHG) {} '.format(primal[-1]))
\ No newline at end of file diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py b/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py deleted file mode 100644 index 5423b22..0000000 --- a/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py +++ /dev/null @@ -1,152 +0,0 @@ -# -*- coding: utf-8 -*- - -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, AcquisitionGeometry, AcquisitionData - -import numpy as np -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction, ScaledFunction - -from ccpi.astra.ops import AstraProjectorSimple, AstraProjectorMC -from skimage.util import random_noise - - -#%%############################################################################### -# Create phantom for TV tomography - -import numpy as np -import matplotlib.pyplot as plt -import os -import tomophantom -from tomophantom import TomoP2D - -model = 102 # note that the selected model is temporal (2D + time) -N = 150 # set dimension of the phantom -# one can specify an exact path to the parameters file -# path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' -path = os.path.dirname(tomophantom.__file__) -path_library2D = os.path.join(path, "Phantom2DLibrary.dat") -#This will generate a N_size x N_size 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 - -#N = 150 -#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 - -#%% -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, channels = np.shape(phantom_2Dt)[0]) -data = ImageData(phantom_2Dt, geometry=ig) - - - -detectors = 150 -angles = np.linspace(0,np.pi,100) - -ag = AcquisitionGeometry('parallel','2D',angles, detectors, channels = np.shape(phantom_2Dt)[0]) -Aop = AstraProjectorMC(ig, ag, 'gpu') -sin = Aop.direct(data) - -plt.imshow(sin.as_array()[10]) -plt.title('Sinogram') -plt.colorbar() -plt.show() - -# Add Gaussian noise to the sinogram data -np.random.seed(10) -n1 = np.random.random(sin.shape) - -noisy_data = sin + ImageData(5*n1) - -plt.imshow(noisy_data.as_array()[10]) -plt.title('Noisy Sinogram') -plt.colorbar() -plt.show() - - -#%% Works only with Composite Operator Structure of PDHG - -#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) - -# Create operators -op1 = Gradient(ig) -op2 = Aop - -# Form Composite Operator -operator = BlockOperator(op1, op2, shape=(2,1) ) - -alpha = 50 -f = BlockFunction( alpha * MixedL21Norm(), \ - 0.5 * L2NormSquared(b = noisy_data) ) -g = ZeroFunction() - -# Compute operator Norm -normK = operator.norm() - -## Primal & dual stepsizes - -sigma = 1 -tau = 1/(sigma*normK**2) - -#sigma = 1/normK -#tau = 1/normK - -opt = {'niter':2000} - -res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) - -plt.figure(figsize=(5,5)) -plt.imshow(res.as_array()) -plt.colorbar() -plt.show() - -#sigma = 10 -#tau = 1/(sigma*normK**2) -# -#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) -#pdhg.max_iteration = 5000 -#pdhg.update_objective_interval = 20 -# -#pdhg.run(5000) -# -##%% -#sol = pdhg.get_output().as_array() -#fig = plt.figure() -#plt.subplot(1,2,1) -#plt.imshow(noisy_data.as_array()) -##plt.colorbar() -#plt.subplot(1,2,2) -#plt.imshow(sol) -##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), sol[int(N/2),:], label = 'Recon') -plt.legend() -plt.show() - - diff --git a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py b/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py deleted file mode 100644 index cb37598..0000000 --- a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py +++ /dev/null @@ -1,187 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.framework import ImageData, ImageGeometry - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \ - MixedL21Norm, BlockFunction - - -from skimage.util import random_noise -from timeit import default_timer as timer - - - -# ############################################################################ -# Create phantom for TV Poisson denoising - -N = 200 -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 - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -ag = ig - -# Create noisy data. Add Gaussian noise -n1 = random_noise(data, mode = 'poisson', seed = 10) -noisy_data = ImageData(n1) - -plt.imshow(noisy_data.as_array()) -plt.colorbar() -plt.show() - -# Regularisation Parameter -alpha = 2 - -#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") - -method = '0' -if method == '0': - - # Create operators - op1 = Gradient(ig) - op2 = Identity(ig, ag) - - # Form Composite Operator - operator = BlockOperator(op1, op2, shape=(2,1) ) - - f1 = alpha * MixedL21Norm() - f2 = KullbackLeibler(noisy_data) - - f = BlockFunction(f1, f2 ) - g = ZeroFunction() - -else: - - ########################################################################### - # No Composite # - ########################################################################### - 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) - -opt = {'niter':2000} -opt1 = {'niter':2000, 'memopt': True} - -t1 = timer() -res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -t2 = timer() - -t3 = timer() -res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -t4 = timer() - -print(pdgap[-1]) - - -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(res.as_array()) -plt.title('no memopt') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(res1.as_array()) -plt.title('memopt') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow((res1 - res).abs().as_array()) -plt.title('diff') -plt.colorbar() -plt.show() -# -plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt') -plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt') -plt.legend() -plt.show() - -print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(t2-t1, t4 -t3)) -diff = (res1 - res).abs().as_array().max() - -print(" Max of abs difference is {}".format(diff)) - - -#%% Check with CVX solution - -from ccpi.optimisation.operators import SparseFiniteDiff - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - ##Construct problem - u1 = Variable(ig.shape) - q = Variable() - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - # Define Total Variation as a regulariser - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u1), DY.matrix() * vec(u1)]), 2, axis = 0)) - - fidelity = sum( u1 - multiply(noisy_data.as_array(), log(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( res.as_array() - u1.value ) - - # Show result - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(res.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), res1.as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), u1.value[int(N/2),:], label = 'CVX') - plt.legend() - - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(primal[-1])) - - - - - |