delete/move demos

3 files changed, 0 insertions, 604 deletions

diff --git a/Wrappers/Python/demos/PDHG_TV_Denoising_Gaussian_3D.py b/Wrappers/Python/demos/PDHG_TV_Denoising_Gaussian_3D.py
deleted file mode 100644
index c86ddc9..0000000
--- a/Wrappers/Python/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/demos/PDHG_TV_Tomo2D.py b/Wrappers/Python/demos/PDHG_TV_Tomo2D.py
deleted file mode 100644
index 87d5328..0000000
--- a/Wrappers/Python/demos/PDHG_TV_Tomo2D.py
+++ /dev/null
@@ -1,245 +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, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \
-                      MixedL21Norm, BlockFunction
-
-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))
-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, dtype=np.float32)
-
-ag = AcquisitionGeometry('parallel','2D',angles, detectors)
-Aop = AstraProjectorSimple(ig, ag, 'cpu')
-sin = Aop.direct(data)
-
-# Create noisy data. Apply Poisson noise
-scale = 2
-n1 = scale * np.random.poisson(sin.as_array()/scale)
-noisy_data = AcquisitionData(n1, ag)
-
-# Regularisation Parameter
-alpha = 5
-
-# Create operators
-op1 = Gradient(ig)
-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()
-    
-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)
-
-
-# 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
-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
diff --git a/Wrappers/Python/demos/PDHG_TV_Tomo2D_gaussian.py b/Wrappers/Python/demos/PDHG_TV_Tomo2D_gaussian.py
deleted file mode 100644
index dc473a8..0000000
--- a/Wrappers/Python/demos/PDHG_TV_Tomo2D_gaussian.py
+++ /dev/null
@@ -1,204 +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:
-
-
-Problem:     min_x, x>0  \alpha * ||\nabla x||_{2,1} + \frac{1}{2}||Ax - g||^{2}
-
-             \nabla: Gradient operator 
-             
-             A: Projection Matrix
-             g: Noisy sinogram corrupted with Gaussian Noise
-             
-             \alpha: Regularization parameter
- 
-"""
-
-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, L2NormSquared, \
-                      MixedL21Norm, BlockFunction
-
-from ccpi.astra.ops import AstraProjectorSimple
-
-
-
-# Create phantom for TV 2D tomography 
-N = 200
-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)
-
-# Create noisy data. Apply Poisson noise
-n1 = np.random.normal(0, 3, size=ig.shape)
-noisy_data = AcquisitionData(n1 + sin.as_array(), ag)
-
-# Show Ground Truth and Noisy Data
-plt.figure(figsize=(10,10))
-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 = 50
-
-# Create operators
-op1 = Gradient(ig)
-op2 = Aop
-
-# 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()
-
-# Compute operator Norm
-normK = operator.norm()
-    
-# Primal & dual stepsizes
-sigma = 10
-tau = 1/(sigma*normK**2)
-
-# Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
-pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 200
-pdhg.run(2000)
-
-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
-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)
-    
-    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
-    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
-
-    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
-    
-    # create matrix representation for Astra operator
-    vol_geom = astra.create_vol_geom(N, N)
-    proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles)
-
-    proj_id = astra.create_projector('strip', proj_geom, vol_geom)
-
-    matrix_id = astra.projector.matrix(proj_id)
-
-    ProjMat = astra.matrix.get(matrix_id)
-    
-    tmp = noisy_data.as_array().ravel()
-    
-    fidelity = 0.5 * sum_squares(ProjMat * u - tmp)
-
-    solver = MOSEK
-    obj =  Minimize( regulariser +  fidelity)
-    prob = Problem(obj)
-    result = prob.solve(verbose = True, solver = solver)    
-         
-    diff_cvx = numpy.abs( pdhg.get_output().as_array() - np.reshape(u.value, (N,N) ))
-           
-    plt.figure(figsize=(15,15))
-    plt.subplot(3,1,1)
-    plt.imshow(pdhg.get_output().as_array())
-    plt.title('PDHG solution')
-    plt.colorbar()
-    plt.subplot(3,1,2)
-    plt.imshow(np.reshape(u.value, (N, N)))
-    plt.title('CVX solution')
-    plt.colorbar()
-    plt.subplot(3,1,3)
-    plt.imshow(diff_cvx)
-    plt.title('Difference')
-    plt.colorbar()
-    plt.show()    
-    
-    plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,N,N), np.reshape(u.value, (N,N) )[int(N/2),:], label = 'CVX')
-    plt.legend()
-    plt.title('Middle Line Profiles')
-    plt.show()
-            
-    print('Primal Objective (CVX) {} '.format(obj.value))
-    print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
\ No newline at end of file