summaryrefslogtreecommitdiffstats
path: root/Wrappers
diff options
context:
space:
mode:
authorepapoutsellis <epapoutsellis@gmail.com>2019-06-12 10:44:43 +0100
committerepapoutsellis <epapoutsellis@gmail.com>2019-06-12 10:44:43 +0100
commit5c3621aeb63e6465f1884be81ebd12d8a39dcdbd (patch)
tree35eae75af08f24d054e9073cd2ea049a3350dc47 /Wrappers
parentec479c78c8956a1645604a780c2eb8c0d601165b (diff)
downloadframework-5c3621aeb63e6465f1884be81ebd12d8a39dcdbd.tar.gz
framework-5c3621aeb63e6465f1884be81ebd12d8a39dcdbd.tar.bz2
framework-5c3621aeb63e6465f1884be81ebd12d8a39dcdbd.tar.xz
framework-5c3621aeb63e6465f1884be81ebd12d8a39dcdbd.zip
delete old demos from wip
Diffstat (limited to 'Wrappers')
-rw-r--r--Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py124
-rw-r--r--Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py225
-rw-r--r--Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_3D.py155
-rw-r--r--Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py207
-rw-r--r--Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py198
-rw-r--r--Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py176
-rw-r--r--Wrappers/Python/wip/Demos/fista_test.py127
7 files changed, 0 insertions, 1212 deletions
diff --git a/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py b/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py
deleted file mode 100644
index 49d4db6..0000000
--- a/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py
+++ /dev/null
@@ -1,124 +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, 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 = 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 = 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/wip/Demos/PDHG_TV_Denoising_Gaussian.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py
deleted file mode 100644
index 5df02b1..0000000
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py
+++ /dev/null
@@ -1,225 +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:
-
- min_{x} max_{y} < K x, y > + g(x) - f^{*}(y)
-
-
-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: K = [ \nabla,
- Identity]
-
- Method = 1: 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 Data import *
-
-#%%
-
-
-data = ImageData(plt.imread('camera.png'))
-
-#
-## Create phantom for TV Gaussian 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
-#data = ImageData(data)
-#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-#ag = ig
-#
-#
-#
-## Replace with http://sipi.usc.edu/database/database.php?volume=misc&image=36#top
-
-
-
-# Create noisy data. Add Gaussian noise
-np.random.seed(10)
-noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.05, size=ig.shape) )
-
-# 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 = 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 = 3000
-pdhg.update_objective_interval = 200
-pdhg.run(3000, verbose=False)
-
-# Show Results
-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 = 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,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.plot(np.linspace(0,N,N), 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/wip/Demos/PDHG_TV_Denoising_Gaussian_3D.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_3D.py
deleted file mode 100644
index c86ddc9..0000000
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_3D.py
+++ /dev/null
@@ -1,155 +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
-
-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 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)
-
-fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8))
-
-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/wip/Demos/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
deleted file mode 100644
index 70f6b9b..0000000
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
+++ /dev/null
@@ -1,207 +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 STFC, 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:
-
- 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
-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, \
- 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)
-
-# Regularisation Parameter
-alpha = 2
-
-method = '1'
-
-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 = ZeroFunction()
-
-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)
-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
-
-def pdgap_objectives(niter, objective, solution):
-
-
- print( "{:04}/{:04} {:<5} {:.4f} {:<5} {:.4f} {:<5} {:.4f}".\
- format(niter, pdhg.max_iteration,'', \
- objective[0],'',\
- objective[1],'',\
- objective[2]))
-
-pdhg.run(2000, callback = pdgap_objectives)
-
-
-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()
-
- 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( 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/wip/Demos/PDHG_TV_Denoising_SaltPepper.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py
deleted file mode 100644
index f5d4ce4..0000000
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py
+++ /dev/null
@@ -1,198 +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.
-
-"""
-
-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
-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)
-
-# 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/wip/Demos/PDHG_Tikhonov_Denoising.py b/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py
deleted file mode 100644
index 041d4ee..0000000
--- a/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py
+++ /dev/null
@@ -1,176 +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
-
-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
-
-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 = 'gaussian', mean=0, var = 0.05, seed=10)
-noisy_data = ImageData(n1)
-
-# Regularisation Parameter
-alpha = 4
-
-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 * L2NormSquared()
- 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 = 0.5 * L2NormSquared(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('Tikhonov 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 = '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]))
-
-
-
-
-
diff --git a/Wrappers/Python/wip/Demos/fista_test.py b/Wrappers/Python/wip/Demos/fista_test.py
deleted file mode 100644
index dd1f6fa..0000000
--- a/Wrappers/Python/wip/Demos/fista_test.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
-
-import numpy as np
-import numpy
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import FISTA, PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Gradient, Identity
-from ccpi.optimisation.functions import L2NormSquared, L1Norm, \
- MixedL21Norm, FunctionOperatorComposition, BlockFunction, ZeroFunction
-
-from skimage.util import random_noise
-
-# 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
-data = ImageData(data)
-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)
-
-# Regularisation Parameter
-alpha = 5
-
-operator = Gradient(ig)
-
-#fidelity = L1Norm(b=noisy_data)
-#regulariser = FunctionOperatorComposition(alpha * L2NormSquared(), operator)
-
-fidelity = FunctionOperatorComposition(alpha * MixedL21Norm(), operator)
-regulariser = 0.5 * L2NormSquared(b = noisy_data)
-
-x_init = ig.allocate()
-
-## Setup and run the PDHG algorithm
-opt = {'tol': 1e-4, 'memopt':True}
-fista = FISTA(x_init=x_init , f=regulariser, g=fidelity, opt=opt)
-fista.max_iteration = 2000
-fista.update_objective_interval = 50
-fista.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(fista.get_output().as_array())
-plt.title('TV Reconstruction')
-plt.colorbar()
-plt.show()
-
-# Compare with PDHG
-method = '0'
-#
-if method == '0':
-#
-# # Create operators
- op1 = Gradient(ig)
- op2 = Identity(ig, ag)
-#
-# # Create BlockOperator
- operator = BlockOperator(op1, op2, shape=(2,1) )
- f = BlockFunction(alpha * L2NormSquared(), fidelity)
- 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 = 50
-pdhg.run(2000)
-#
-#%%
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(fista.get_output().as_array())
-plt.title('FISTA')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(pdhg.get_output().as_array())
-plt.title('PDHG')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow(np.abs(pdhg.get_output().as_array()-fista.get_output().as_array()))
-plt.title('Diff FISTA-PDHG')
-plt.colorbar()
-plt.show()
-
-