diff options
5 files changed, 185 insertions, 115 deletions
diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py index a8a0ab4..387b5c1 100755 --- a/Wrappers/Python/ccpi/framework/framework.py +++ b/Wrappers/Python/ccpi/framework/framework.py @@ -760,13 +760,6 @@ class DataContainer(object): return numpy.sqrt(self.squared_norm()) -# def dot(self, other, *args, **kwargs): -# '''return the inner product of 2 DataContainers viewed as vectors''' -# if self.shape == other.shape: -# return numpy.dot(self.as_array().ravel(), other.as_array().ravel()) -# else: -# raise ValueError('Shapes are not aligned: {} != {}'.format(self.shape, other.shape)) - def dot(self, other, *args, **kwargs): '''return the inner product of 2 DataContainers viewed as vectors''' method = kwargs.get('method', 'reduce') diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py index e0b8a32..6ffaf70 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py @@ -111,7 +111,7 @@ class Gradient(LinearOperator): return BlockDataContainer(*mat) - def sum_abs_row(self): + def sum_abs_col(self): tmp = self.gm_range.allocate() res = self.gm_domain.allocate() @@ -120,7 +120,7 @@ class Gradient(LinearOperator): res += spMat.sum_abs_row() return res - def sum_abs_col(self): + def sum_abs_row(self): tmp = self.gm_range.allocate() res = [] diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py index 4903c44..3c295f5 100644 --- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py @@ -17,6 +17,27 @@ # See the License for the specific language governing permissions and # limitations under the License. +""" + +Total Variation Denoising using PDHG algorithm: + + min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) + + +Problem: min_x, x>0 \alpha * ||\nabla x||_{1} + \int x - g * log(x) + + \nabla: Gradient operator + g: Noisy Data with Poisson Noise + \alpha: Regularization parameter + + Method = 0: K = [ \nabla, + Identity] + + Method = 1: K = \nabla + + +""" + from ccpi.framework import ImageData, ImageGeometry import numpy as np diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py index 4189acb..f5d4ce4 100644 --- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py @@ -17,6 +17,27 @@ # See the License for the specific language governing permissions and # limitations under the License. +""" + +Total Variation Denoising using PDHG algorithm: + + min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) + + +Problem: min_x, x>0 \alpha * ||\nabla x||_{1} + ||x-g||_{1} + + \nabla: Gradient operator + g: Noisy Data with Salt & Pepper Noise + \alpha: Regularization parameter + + Method = 0: K = [ \nabla, + Identity] + + Method = 1: K = \nabla + + +""" + from ccpi.framework import ImageData, ImageGeometry import numpy as np diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py b/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py index 0711e91..75286e5 100644 --- a/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py @@ -31,6 +31,25 @@ from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \ from ccpi.astra.ops import AstraProjectorSimple +""" + +Total Variation Denoising using PDHG algorithm: + + min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) + + +Problem: min_x, x>0 \alpha * ||\nabla x||_{1} + int A x -g log(Ax + \eta) + + \nabla: Gradient operator + + A: Projection Matrix + g: Noisy sinogram corrupted with Poisson Noise + + \eta: Background Noise + \alpha: Regularization parameter + +""" + # Create phantom for TV 2D tomography N = 75 x = np.zeros((N,N)) @@ -70,12 +89,26 @@ f = BlockFunction(f1, f2) g = ZeroFunction() -# Compute operator Norm -normK = operator.norm() +diag_precon = True + +if diag_precon: + + def tau_sigma_precond(operator): + + tau = 1/operator.sum_abs_row() + sigma = 1/ operator.sum_abs_col() + + return tau, sigma + + tau, sigma = tau_sigma_precond(operator) + +else: + # Compute operator Norm + normK = operator.norm() + # Primal & dual stepsizes + sigma = 10 + tau = 1/(sigma*normK**2) -# Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) # Setup and run the PDHG algorithm @@ -84,6 +117,8 @@ 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()) @@ -108,104 +143,104 @@ plt.show() #%% Check with CVX solution -from ccpi.optimisation.operators import SparseFiniteDiff -import astra -import numpy - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - - ##Construct problem - u = Variable(N*N) - #q = Variable() - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - - # create matrix representation for Astra operator - - vol_geom = astra.create_vol_geom(N, N) - proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) - - proj_id = astra.create_projector('strip', proj_geom, vol_geom) - - matrix_id = astra.projector.matrix(proj_id) - - ProjMat = astra.matrix.get(matrix_id) - - fidelity = sum( ProjMat * u - noisy_data.as_array().ravel() * log(ProjMat * u)) - #constraints = [q>= fidelity, u>=0] - constraints = [u>=0] - - solver = SCS - obj = Minimize( regulariser + fidelity) - prob = Problem(obj, constraints) - result = prob.solve(verbose = True, solver = solver) - - -##%% Check with CVX solution - -from ccpi.optimisation.operators import SparseFiniteDiff - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - ##Construct problem - u = Variable(ig.shape) - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - # Define Total Variation as a regulariser - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - fidelity = pnorm( u - noisy_data.as_array(),1) - - # choose solver - if 'MOSEK' in installed_solvers(): - solver = MOSEK - else: - solver = SCS - - obj = Minimize( regulariser + fidelity) - prob = Problem(obj) - result = prob.solve(verbose = True, solver = solver) - - - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(pdhg.get_output().as_array()) - plt.title('PDHG solution') - plt.colorbar() - plt.subplot(3,1,2) - plt.imshow(np.reshape(u.value, (N, N))) - plt.title('CVX solution') - plt.colorbar() - plt.subplot(3,1,3) - plt.imshow(diff_cvx) - plt.title('Difference') - plt.colorbar() - plt.show() - - plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') - plt.legend() - plt.title('Middle Line Profiles') - plt.show() - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
\ No newline at end of file +#from ccpi.optimisation.operators import SparseFiniteDiff +#import astra +#import numpy +# +#try: +# from cvxpy import * +# cvx_not_installable = True +#except ImportError: +# cvx_not_installable = False +# +# +#if cvx_not_installable: +# +# +# ##Construct problem +# u = Variable(N*N) +# #q = Variable() +# +# DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') +# DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') +# +# regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) +# +# # create matrix representation for Astra operator +# +# vol_geom = astra.create_vol_geom(N, N) +# proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) +# +# proj_id = astra.create_projector('strip', proj_geom, vol_geom) +# +# matrix_id = astra.projector.matrix(proj_id) +# +# ProjMat = astra.matrix.get(matrix_id) +# +# fidelity = sum( ProjMat * u - noisy_data.as_array().ravel() * log(ProjMat * u)) +# #constraints = [q>= fidelity, u>=0] +# constraints = [u>=0] +# +# solver = SCS +# obj = Minimize( regulariser + fidelity) +# prob = Problem(obj, constraints) +# result = prob.solve(verbose = True, solver = solver) +# +# +###%% Check with CVX solution +# +#from ccpi.optimisation.operators import SparseFiniteDiff +# +#try: +# from cvxpy import * +# cvx_not_installable = True +#except ImportError: +# cvx_not_installable = False +# +# +#if cvx_not_installable: +# +# ##Construct problem +# u = Variable(ig.shape) +# +# DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') +# DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') +# +# # Define Total Variation as a regulariser +# regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) +# fidelity = pnorm( u - noisy_data.as_array(),1) +# +# # choose solver +# if 'MOSEK' in installed_solvers(): +# solver = MOSEK +# else: +# solver = SCS +# +# obj = Minimize( regulariser + fidelity) +# prob = Problem(obj) +# result = prob.solve(verbose = True, solver = solver) +# +# +# plt.figure(figsize=(15,15)) +# plt.subplot(3,1,1) +# plt.imshow(pdhg.get_output().as_array()) +# plt.title('PDHG solution') +# plt.colorbar() +# plt.subplot(3,1,2) +# plt.imshow(np.reshape(u.value, (N, N))) +# plt.title('CVX solution') +# plt.colorbar() +# plt.subplot(3,1,3) +# plt.imshow(diff_cvx) +# plt.title('Difference') +# plt.colorbar() +# plt.show() +# +# plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') +# plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') +# plt.legend() +# plt.title('Middle Line Profiles') +# plt.show() +# +# print('Primal Objective (CVX) {} '.format(obj.value)) +# print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
\ No newline at end of file |