summaryrefslogtreecommitdiffstats
path: root/Wrappers/Python
diff options
context:
space:
mode:
authorepapoutsellis <epapoutsellis@gmail.com>2019-06-12 15:22:05 +0100
committerepapoutsellis <epapoutsellis@gmail.com>2019-06-12 15:22:05 +0100
commitc2ec8d85841b059437d9e97a46540ee4e712b593 (patch)
tree69503b36d2288cfc3125d822605ed8dd1c752bc6 /Wrappers/Python
parent144c23b09281a4bdd767ea89db70d028cda05b40 (diff)
downloadframework-c2ec8d85841b059437d9e97a46540ee4e712b593.tar.gz
framework-c2ec8d85841b059437d9e97a46540ee4e712b593.tar.bz2
framework-c2ec8d85841b059437d9e97a46540ee4e712b593.tar.xz
framework-c2ec8d85841b059437d9e97a46540ee4e712b593.zip
delete old demo
Diffstat (limited to 'Wrappers/Python')
-rwxr-xr-xWrappers/Python/ccpi/optimisation/operators/BlockOperator.py2
-rw-r--r--Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py194
2 files changed, 1 insertions, 195 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
index 5f04363..cbdc420 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
@@ -104,7 +104,7 @@ class BlockOperator(Operator):
index = row*self.shape[1]+col
return self.operators[index]
- def calculate_norm(self, **kwargs):
+ def norm(self, **kwargs):
norm = [op.norm(**kwargs)**2 for op in self.operators]
return numpy.sqrt(sum(norm))
diff --git a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py
deleted file mode 100644
index e74e1c6..0000000
--- a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py
+++ /dev/null
@@ -1,194 +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 Generalised Variation (TGV) Tomography 2D using PDHG algorithm:
-
-
-Problem: min_{x>0} \alpha * ||\nabla x - w||_{2,1} + \beta * || E w ||_{2,1} +
- \frac{1}{2}||Au - g||^{2}
-
- min_{u>0} \alpha * ||\nabla u - w||_{2,1} + \beta * || E w ||_{2,1} +
- int A u - g log(Au + \eta)
-
- \alpha: Regularization parameter
- \beta: Regularization parameter
-
- \nabla: Gradient operator
- E: Symmetrized Gradient operator
- A: System Matrix
-
- g: Noisy Sinogram
-
- K = [ \nabla, - Identity
- ZeroOperator, E
- A, ZeroOperator]
-
-"""
-
-from ccpi.framework import AcquisitionGeometry, AcquisitionData, 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, Gradient, Identity, \
- SymmetrizedGradient, ZeroOperator
-from ccpi.optimisation.functions import IndicatorBox, KullbackLeibler, ZeroFunction,\
- MixedL21Norm, BlockFunction, L2NormSquared
-
-from ccpi.astra.ops import AstraProjectorSimple
-import os, sys
-
-
-import tomophantom
-from tomophantom import TomoP2D
-
-# user supplied input
-if len(sys.argv) > 1:
- which_noise = int(sys.argv[1])
-else:
- which_noise = 1
-
-# Load Piecewise smooth Shepp-Logan phantom
-model = 2 # select a model number from the library
-N = 128 # 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)
-
-#Create Acquisition Data
-detectors = N
-angles = np.linspace(0, np.pi, N)
-ag = AcquisitionGeometry('parallel','2D',angles, detectors)
-
-#device = input('Available device: GPU==1 / CPU==0 ')
-device = '1'
-if device=='1':
- dev = 'gpu'
-else:
- dev = 'cpu'
-
-Aop = AstraProjectorSimple(ig, ag, 'cpu')
-sin = Aop.direct(data)
-
-# Create noisy sinogram.
-noises = ['gaussian', 'poisson']
-noise = noises[which_noise]
-
-if noise == 'poisson':
- scale = 5
- eta = 0
- noisy_data = AcquisitionData(np.random.poisson( scale * (eta + sin.as_array()))/scale, ag)
-elif noise == 'gaussian':
- n1 = np.random.normal(0, 1, size = ag.shape)
- noisy_data = AcquisitionData(n1 + sin.as_array(), ag)
-else:
- raise ValueError('Unsupported Noise ', noise)
-
-# Show Ground Truth and Noisy Data
-plt.figure(figsize=(10,10))
-plt.subplot(1,2,2)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(1,2,1)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.show()
-
-# 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) )
-normK = operator.norm()
-
-# Create functions
-if noise == 'poisson':
- alpha = 3
- beta = 6
- f3 = KullbackLeibler(noisy_data)
- g = BlockFunction(IndicatorBox(lower=0), ZeroFunction())
-
- # Primal & dual stepsizes
- sigma = 1
- tau = 1/(sigma*normK**2)
-
-elif noise == 'gaussian':
- alpha = 20
- beta = 50
- f3 = 0.5 * L2NormSquared(b=noisy_data)
- g = BlockFunction(ZeroFunction(), ZeroFunction())
-
- # Primal & dual stepsizes
- sigma = 10
- tau = 1/(sigma*normK**2)
-
-f1 = alpha * MixedL21Norm()
-f2 = beta * MixedL21Norm()
-f = BlockFunction(f1, f2, f3)
-
-# Compute operator Norm
-normK = operator.norm()
-
-# 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 = 500
-pdhg.run(3000)
-#%%
-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,ig.shape[1],ig.shape[1]), data.as_array()[:, int(N/2)], label = 'GTruth')
-plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output()[0].as_array()[:, int(N/2)], label = 'TGV reconstruction')
-plt.legend()
-plt.title('Middle Line Profiles')
-plt.show()
-
-