diff options
Diffstat (limited to 'Wrappers')
| -rw-r--r-- | Wrappers/Python/demos/PDHG_TV_Tomo2D_poisson.py | 205 | ||||
| -rw-r--r-- | Wrappers/Python/demos/PDHG_examples/PDHG_2D_time_denoising.py | 169 | ||||
| -rw-r--r-- | Wrappers/Python/wip/demo_box_constraints_FISTA.py | 2 |
3 files changed, 276 insertions, 100 deletions
diff --git a/Wrappers/Python/demos/PDHG_TV_Tomo2D_poisson.py b/Wrappers/Python/demos/PDHG_TV_Tomo2D_poisson.py index b6d7725..72d0670 100644 --- a/Wrappers/Python/demos/PDHG_TV_Tomo2D_poisson.py +++ b/Wrappers/Python/demos/PDHG_TV_Tomo2D_poisson.py @@ -29,7 +29,7 @@ from ccpi.optimisation.algorithms import PDHG from ccpi.optimisation.operators import BlockOperator, Gradient from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \ - MixedL21Norm, BlockFunction + MixedL21Norm, BlockFunction, IndicatorBox from ccpi.astra.ops import AstraProjectorSimple @@ -67,10 +67,13 @@ Aop = AstraProjectorSimple(ig, ag, 'cpu') sin = Aop.direct(data) # Create noisy data. Apply Poisson noise -scale = 0.5 -n1 = scale * np.random.poisson(sin.as_array()/scale) +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) @@ -83,9 +86,10 @@ plt.title('Noisy Data') plt.colorbar() plt.show() +#%% # Regularisation Parameter -alpha = 0.5 +alpha = 2 # Create operators op1 = Gradient(ig) @@ -100,20 +104,23 @@ f1 = alpha * MixedL21Norm() f2 = KullbackLeibler(noisy_data) f = BlockFunction(f1, f2) -g = ZeroFunction() +#g = ZeroFunction() +g = IndicatorBox(lower=0) # Compute operator Norm normK = operator.norm() # Primal & dual stepsizes -sigma = 10 +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 = 200 -pdhg.run(2000) +pdhg.update_objective_interval = 500 +pdhg.run(2000, verbose = True) plt.figure(figsize=(15,15)) plt.subplot(3,1,1) @@ -137,64 +144,11 @@ 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 +##%% Check with CVX solution # #from ccpi.optimisation.operators import SparseFiniteDiff +#import astra +#import numpy # #try: # from cvxpy import * @@ -204,48 +158,101 @@ if cvx_not_installable: # # #if cvx_not_installable: +# # # ##Construct problem -# u = Variable(ig.shape) +# u = Variable(N*N) +# #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(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 +# # 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) -# result = prob.solve(verbose = True, solver = solver) +# 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.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 +# 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_examples/PDHG_2D_time_denoising.py b/Wrappers/Python/demos/PDHG_examples/PDHG_2D_time_denoising.py new file mode 100644 index 0000000..045458a --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/PDHG_2D_time_denoising.py @@ -0,0 +1,169 @@ +# -*- 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/wip/demo_box_constraints_FISTA.py b/Wrappers/Python/wip/demo_box_constraints_FISTA.py index 2f9e0c6..b15dd45 100644 --- a/Wrappers/Python/wip/demo_box_constraints_FISTA.py +++ b/Wrappers/Python/wip/demo_box_constraints_FISTA.py @@ -72,7 +72,7 @@ else: # Set up Operator object combining the ImageGeometry and AcquisitionGeometry # wrapping calls to ASTRA as well as specifying whether to use CPU or GPU. -Aop = AstraProjectorSimple(ig, ag, 'gpu') +Aop = AstraProjectorSimple(ig, ag, 'cpu') Aop = Identity(ig,ig) |