diff options
22 files changed, 1881 insertions, 1295 deletions
diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..d9bd16b --- /dev/null +++ b/.gitattributes @@ -0,0 +1 @@ +*.py text eol=lf diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index d1ad26b..dab2dd9 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -3,7 +3,7 @@ # Visual Analytics and Imaging System Group of the Science Technology # Facilities Council, STFC -# Copyright 2018 Edoardo Pasca +# Copyright 2018-2019 Edoardo Pasca # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. @@ -26,6 +26,7 @@ import numpy import sys from datetime import timedelta, datetime import warnings +from functools import reduce def find_key(dic, val): """return the key of dictionary dic given the value""" @@ -729,6 +730,16 @@ class DataContainer(object): ## reductions def sum(self, out=None, *args, **kwargs): return self.as_array().sum(*args, **kwargs) + def squared_norm(self): + '''return the squared euclidean norm of the DataContainer viewed as a vector''' + shape = self.shape + size = reduce(lambda x,y:x*y, shape, 1) + y = numpy.reshape(self.as_array(), (size, )) + return numpy.dot(y, y.conjugate()) + def norm(self): + '''return the euclidean norm of the DataContainer viewed as a vector''' + return numpy.sqrt(self.squared_norm()) + class ImageData(DataContainer): diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py new file mode 100755 index 0000000..680b268 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -0,0 +1,157 @@ +# -*- 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 2019 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. +import time +from numbers import Integral + +class Algorithm(object): + '''Base class for iterative algorithms + + provides the minimal infrastructure. + Algorithms are iterables so can be easily run in a for loop. They will + stop as soon as the stop cryterion is met. + The user is required to implement the set_up, __init__, update and + and update_objective methods + + A courtesy method run is available to run n iterations. The method accepts + a callback function that receives the current iteration number and the actual objective + value and can be used to trigger print to screens and other user interactions. The run + method will stop when the stopping cryterion is met. + ''' + + def __init__(self): + '''Constructor + + Set the minimal number of parameters: + iteration: current iteration number + max_iteration: maximum number of iterations + memopt: whether to use memory optimisation () + timing: list to hold the times it took to run each iteration + update_objectice_interval: the interval every which we would save the current + objective. 1 means every iteration, 2 every 2 iteration + and so forth. This is by default 1 and should be increased + when evaluating the objective is computationally expensive. + ''' + self.iteration = 0 + self.__max_iteration = 0 + self.__loss = [] + self.memopt = False + self.timing = [] + self.update_objective_interval = 1 + def set_up(self, *args, **kwargs): + '''Set up the algorithm''' + raise NotImplementedError() + def update(self): + '''A single iteration of the algorithm''' + raise NotImplementedError() + + def should_stop(self): + '''default stopping cryterion: number of iterations + + The user can change this in concrete implementatition of iterative algorithms.''' + return self.max_iteration_stop_cryterion() + + def max_iteration_stop_cryterion(self): + '''default stop cryterion for iterative algorithm: max_iteration reached''' + return self.iteration >= self.max_iteration + def __iter__(self): + '''Algorithm is an iterable''' + return self + def next(self): + '''Algorithm is an iterable + + python2 backwards compatibility''' + return self.__next__() + def __next__(self): + '''Algorithm is an iterable + + calling this method triggers update and update_objective + ''' + if self.should_stop(): + raise StopIteration() + else: + time0 = time.time() + self.update() + self.timing.append( time.time() - time0 ) + if self.iteration % self.update_objective_interval == 0: + self.update_objective() + self.iteration += 1 + def get_output(self): + '''Returns the solution found''' + return self.x + def get_last_loss(self): + '''Returns the last stored value of the loss function + + if update_objective_interval is 1 it is the value of the objective at the current + iteration. If update_objective_interval > 1 it is the last stored value. + ''' + return self.__loss[-1] + def get_last_objective(self): + '''alias to get_last_loss''' + return self.get_last_loss() + def update_objective(self): + '''calculates the objective with the current solution''' + raise NotImplementedError() + @property + def loss(self): + '''returns the list of the values of the objective during the iteration + + The length of this list may be shorter than the number of iterations run when + the update_objective_interval > 1 + ''' + return self.__loss + @property + def objective(self): + '''alias of loss''' + return self.loss + @property + def max_iteration(self): + '''gets the maximum number of iterations''' + return self.__max_iteration + @max_iteration.setter + def max_iteration(self, value): + '''sets the maximum number of iterations''' + assert isinstance(value, int) + self.__max_iteration = value + @property + def update_objective_interval(self): + return self.__update_objective_interval + @update_objective_interval.setter + def update_objective_interval(self, value): + if isinstance(value, Integral): + if value >= 1: + self.__update_objective_interval = value + else: + raise ValueError('Update objective interval must be an integer >= 1') + else: + raise ValueError('Update objective interval must be an integer >= 1') + def run(self, iterations, verbose=False, callback=None): + '''run n iterations and update the user with the callback if specified''' + if self.should_stop(): + print ("Stop cryterion has been reached.") + i = 0 + for _ in self: + if verbose: + print ("Iteration {}/{}, objective {}".format(self.iteration, + self.max_iteration, self.get_last_objective()) ) + if callback is not None: + callback(self.iteration, self.get_last_objective()) + i += 1 + if i == iterations: + break + diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py new file mode 100755 index 0000000..7194eb8 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py @@ -0,0 +1,87 @@ +# -*- 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 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. +""" +Created on Thu Feb 21 11:11:23 2019 + +@author: ofn77899 +""" + +from ccpi.optimisation.algorithms import Algorithm +#from collections.abc import Iterable +class CGLS(Algorithm): + + '''Conjugate Gradient Least Squares algorithm + + Parameters: + x_init: initial guess + operator: operator for forward/backward projections + data: data to operate on + ''' + def __init__(self, **kwargs): + super(CGLS, self).__init__() + self.x = kwargs.get('x_init', None) + self.operator = kwargs.get('operator', None) + self.data = kwargs.get('data', None) + if self.x is not None and self.operator is not None and \ + self.data is not None: + print ("Calling from creator") + self.set_up(x_init =kwargs['x_init'], + operator=kwargs['operator'], + data =kwargs['data']) + + def set_up(self, x_init, operator , data ): + + self.r = data.copy() + self.x = x_init.copy() + + self.operator = operator + self.d = operator.adjoint(self.r) + + + self.normr2 = self.d.squared_norm() + #if isinstance(self.normr2, Iterable): + # self.normr2 = sum(self.normr2) + #self.normr2 = numpy.sqrt(self.normr2) + #print ("set_up" , self.normr2) + + def update(self): + + Ad = self.operator.direct(self.d) + #norm = (Ad*Ad).sum() + #if isinstance(norm, Iterable): + # norm = sum(norm) + norm = Ad.squared_norm() + + alpha = self.normr2/norm + self.x += (self.d * alpha) + self.r -= (Ad * alpha) + s = self.operator.adjoint(self.r) + + normr2_new = s.squared_norm() + #if isinstance(normr2_new, Iterable): + # normr2_new = sum(normr2_new) + #normr2_new = numpy.sqrt(normr2_new) + #print (normr2_new) + + beta = normr2_new/self.normr2 + self.normr2 = normr2_new + self.d = s + beta*self.d + + def update_objective(self): + self.loss.append(self.r.squared_norm())
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py b/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py new file mode 100644 index 0000000..798fb61 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py @@ -0,0 +1,86 @@ +# -*- 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 2019 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. +""" +Created on Thu Feb 21 11:09:03 2019 + +@author: ofn77899 +""" + +from ccpi.optimisation.algorithms import Algorithm +from ccpi.optimisation.funcs import ZeroFun + +class FBPD(Algorithm): + '''FBPD Algorithm + + Parameters: + x_init: initial guess + f: constraint + g: data fidelity + h: regularizer + opt: additional algorithm + ''' + constraint = None + data_fidelity = None + regulariser = None + def __init__(self, **kwargs): + pass + def set_up(self, x_init, operator=None, constraint=None, data_fidelity=None,\ + regulariser=None, opt=None): + + # default inputs + if constraint is None: + self.constraint = ZeroFun() + else: + self.constraint = constraint + if data_fidelity is None: + data_fidelity = ZeroFun() + else: + self.data_fidelity = data_fidelity + if regulariser is None: + self.regulariser = ZeroFun() + else: + self.regulariser = regulariser + + # algorithmic parameters + + + # step-sizes + self.tau = 2 / (self.data_fidelity.L + 2) + self.sigma = (1/self.tau - self.data_fidelity.L/2) / self.regulariser.L + + self.inv_sigma = 1/self.sigma + + # initialization + self.x = x_init + self.y = operator.direct(self.x) + + + def update(self): + + # primal forward-backward step + x_old = self.x + self.x = self.x - self.tau * ( self.data_fidelity.grad(self.x) + self.operator.adjoint(self.y) ) + self.x = self.constraint.prox(self.x, self.tau); + + # dual forward-backward step + self.y = self.y + self.sigma * self.operator.direct(2*self.x - x_old); + self.y = self.y - self.sigma * self.regulariser.prox(self.inv_sigma*self.y, self.inv_sigma); + + # time and criterion + self.loss = self.constraint(self.x) + self.data_fidelity(self.x) + self.regulariser(self.operator.direct(self.x)) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py new file mode 100755 index 0000000..bc4489e --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py @@ -0,0 +1,121 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Feb 21 11:07:30 2019 + +@author: ofn77899 +""" + +from ccpi.optimisation.algorithms import Algorithm +from ccpi.optimisation.funcs import ZeroFun +import numpy + +class FISTA(Algorithm): + '''Fast Iterative Shrinkage-Thresholding Algorithm + + Beck, A. and Teboulle, M., 2009. A fast iterative shrinkage-thresholding + algorithm for linear inverse problems. + SIAM journal on imaging sciences,2(1), pp.183-202. + + Parameters: + x_init: initial guess + f: data fidelity + g: regularizer + h: + opt: additional algorithm + ''' + + def __init__(self, **kwargs): + '''initialisation can be done at creation time if all + proper variables are passed or later with set_up''' + super(FISTA, self).__init__() + self.f = None + self.g = None + self.invL = None + self.t_old = 1 + args = ['x_init', 'f', 'g', 'opt'] + for k,v in kwargs.items(): + if k in args: + args.pop(args.index(k)) + if len(args) == 0: + return self.set_up(kwargs['x_init'], + f=kwargs['f'], + g=kwargs['g'], + opt=kwargs['opt']) + + def set_up(self, x_init, f=None, g=None, opt=None): + + # default inputs + if f is None: + self.f = ZeroFun() + else: + self.f = f + if g is None: + g = ZeroFun() + self.g = g + else: + self.g = g + + # algorithmic parameters + if opt is None: + opt = {'tol': 1e-4, 'memopt':False} + + self.tol = opt['tol'] if 'tol' in opt.keys() else 1e-4 + memopt = opt['memopt'] if 'memopt' in opt.keys() else False + self.memopt = memopt + + # initialization + if memopt: + self.y = x_init.clone() + self.x_old = x_init.clone() + self.x = x_init.clone() + self.u = x_init.clone() + else: + self.x_old = x_init.copy() + self.y = x_init.copy() + + #timing = numpy.zeros(max_iter) + #criter = numpy.zeros(max_iter) + + + self.invL = 1/f.L + + self.t_old = 1 + + def update(self): + # algorithm loop + #for it in range(0, max_iter): + + if self.memopt: + # u = y - invL*f.grad(y) + # store the result in x_old + self.f.gradient(self.y, out=self.u) + self.u.__imul__( -self.invL ) + self.u.__iadd__( self.y ) + # x = g.prox(u,invL) + self.g.proximal(self.u, self.invL, out=self.x) + + self.t = 0.5*(1 + numpy.sqrt(1 + 4*(self.t_old**2))) + + # y = x + (t_old-1)/t*(x-x_old) + self.x.subtract(self.x_old, out=self.y) + self.y.__imul__ ((self.t_old-1)/self.t) + self.y.__iadd__( self.x ) + + self.x_old.fill(self.x) + self.t_old = self.t + + + else: + u = self.y - self.invL*self.f.grad(self.y) + + self.x = self.g.prox(u,self.invL) + + self.t = 0.5*(1 + numpy.sqrt(1 + 4*(self.t_old**2))) + + self.y = self.x + (self.t_old-1)/self.t*(self.x-self.x_old) + + self.x_old = self.x.copy() + self.t_old = self.t + + def update_objective(self): + self.loss.append( self.f(self.x) + self.g(self.x) )
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py new file mode 100755 index 0000000..7794b4d --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py @@ -0,0 +1,72 @@ +# -*- 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 2019 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. +""" +Created on Thu Feb 21 11:05:09 2019 + +@author: ofn77899 +""" +from ccpi.optimisation.algorithms import Algorithm + +class GradientDescent(Algorithm): + '''Implementation of Gradient Descent algorithm + ''' + + def __init__(self, **kwargs): + '''initialisation can be done at creation time if all + proper variables are passed or later with set_up''' + super(GradientDescent, self).__init__() + self.x = None + self.rate = 0 + self.objective_function = None + self.regulariser = None + args = ['x_init', 'objective_function', 'rate'] + for k,v in kwargs.items(): + if k in args: + args.pop(args.index(k)) + if len(args) == 0: + return self.set_up(x_init=kwargs['x_init'], + objective_function=kwargs['objective_function'], + rate=kwargs['rate']) + + def should_stop(self): + '''stopping cryterion, currently only based on number of iterations''' + return self.iteration >= self.max_iteration + + def set_up(self, x_init, objective_function, rate): + '''initialisation of the algorithm''' + self.x = x_init.copy() + if self.memopt: + self.x_update = x_init.copy() + self.objective_function = objective_function + self.rate = rate + self.loss.append(objective_function(x_init)) + self.iteration = 0 + + def update(self): + '''Single iteration''' + if self.memopt: + self.objective_function.gradient(self.x, out=self.x_update) + self.x_update *= -self.rate + self.x += self.x_update + else: + self.x += -self.rate * self.objective_function.grad(self.x) + + def update_objective(self): + self.loss.append(self.objective_function(self.x)) +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py new file mode 100644 index 0000000..4a3830f --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py @@ -0,0 +1,29 @@ +# -*- 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 2019 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. +""" +Created on Thu Feb 21 11:03:13 2019 + +@author: ofn77899 +""" + +from .Algorithm import Algorithm +from .CGLS import CGLS +from .GradientDescent import GradientDescent +from .FISTA import FISTA +from .FBPD import FBPD
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/algs.py b/Wrappers/Python/ccpi/optimisation/algs.py index a736277..15638a9 100755 --- a/Wrappers/Python/ccpi/optimisation/algs.py +++ b/Wrappers/Python/ccpi/optimisation/algs.py @@ -158,9 +158,8 @@ def FBPD(x_init, operator=None, constraint=None, data_fidelity=None,\ memopt = opt['memopts'] if 'memopts' in opt.keys() else False # step-sizes - tau = 2 / (data_fidelity.L + 2); - sigma = (1/tau - data_fidelity.L/2) / regulariser.L; - + tau = 2 / (data_fidelity.L + 2) + sigma = (1/tau - data_fidelity.L/2) / regulariser.L inv_sigma = 1/sigma # initialization diff --git a/Wrappers/Python/ccpi/optimisation/funcs.py b/Wrappers/Python/ccpi/optimisation/funcs.py index c7a6143..9b9fc36 100755 --- a/Wrappers/Python/ccpi/optimisation/funcs.py +++ b/Wrappers/Python/ccpi/optimisation/funcs.py @@ -217,10 +217,10 @@ class ZeroFun(Function): class Norm1(Function): def __init__(self,gamma): + super(Norm1, self).__init__() self.gamma = gamma self.L = 1 self.sign_x = None - super(Norm1, self).__init__() def __call__(self,x,out=None): if out is None: diff --git a/Wrappers/Python/ccpi/optimisation/ops.py b/Wrappers/Python/ccpi/optimisation/ops.py index a0e1713..e9e7f44 100755 --- a/Wrappers/Python/ccpi/optimisation/ops.py +++ b/Wrappers/Python/ccpi/optimisation/ops.py @@ -93,9 +93,10 @@ class Identity(Operator): class TomoIdentity(Operator): def __init__(self, geometry, **kwargs): + super(TomoIdentity, self).__init__() self.s1 = 1.0 self.geometry = geometry - super(TomoIdentity, self).__init__() + def direct(self,x,out=None): @@ -236,10 +237,11 @@ def PowerMethodNonsquare(op,numiters , x0=None): # Loop for it in numpy.arange(numiters): x1 = op.adjoint(op.direct(x0)) - x1norm = numpy.sqrt((x1*x1).sum()) + #x1norm = numpy.sqrt((x1*x1).sum()) + x1norm = x1.norm() #print ("x0 **********" ,x0) #print ("x1 **********" ,x1) - s[it] = (x1*x0).sum() / (x0*x0).sum() + s[it] = (x1*x0).sum() / (x0.squared_norm()) x0 = (1.0/x1norm)*x1 return numpy.sqrt(s[-1]), numpy.sqrt(s), x0 diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml index 1b7cae6..8ded429 100644 --- a/Wrappers/Python/conda-recipe/meta.yaml +++ b/Wrappers/Python/conda-recipe/meta.yaml @@ -12,6 +12,15 @@ test: requires: - python-wget - cvxpy # [not win] + + source_files: + - ./test # [win] + - ./ccpi/Wrappers/Python/test # [not win] + + commands: + - python -c "import os; print (os.getcwd())" + - python -m unittest discover # [win] + - python -m unittest discover -s ccpi/Wrappers/Python/test # [not win] requirements: build: diff --git a/Wrappers/Python/conda-recipe/run_test.py b/Wrappers/Python/conda-recipe/run_test.py deleted file mode 100755 index 57afd57..0000000 --- a/Wrappers/Python/conda-recipe/run_test.py +++ /dev/null @@ -1,1145 +0,0 @@ -import unittest
-import numpy
-import numpy as np
-from ccpi.framework import DataContainer
-from ccpi.framework import ImageData
-from ccpi.framework import AcquisitionData
-from ccpi.framework import ImageGeometry
-from ccpi.framework import AcquisitionGeometry
-import sys
-from timeit import default_timer as timer
-from ccpi.optimisation.algs import FISTA
-from ccpi.optimisation.algs import FBPD
-from ccpi.optimisation.funcs import Norm2sq
-from ccpi.optimisation.funcs import ZeroFun
-from ccpi.optimisation.funcs import Norm1
-from ccpi.optimisation.funcs import TV2D
-from ccpi.optimisation.funcs import Norm2
-
-from ccpi.optimisation.ops import LinearOperatorMatrix
-from ccpi.optimisation.ops import TomoIdentity
-from ccpi.optimisation.ops import Identity
-
-
-import numpy.testing
-import wget
-import os
-from ccpi.io.reader import NexusReader
-
-
-try:
- from cvxpy import *
- cvx_not_installable = False
-except ImportError:
- cvx_not_installable = True
-
-
-def aid(x):
- # This function returns the memory
- # block address of an array.
- return x.__array_interface__['data'][0]
-
-
-def dt(steps):
- return steps[-1] - steps[-2]
-
-
-class TestDataContainer(unittest.TestCase):
- def create_simple_ImageData(self):
- N = 64
- ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N)
- Phantom = ImageData(geometry=ig)
-
- x = Phantom.as_array()
-
- x[int(round(N/4)):int(round(3*N/4)),
- int(round(N/4)):int(round(3*N/4))] = 0.5
- x[int(round(N/8)):int(round(7*N/8)),
- int(round(3*N/8)):int(round(5*N/8))] = 1
-
- return (ig, Phantom)
-
- def create_DataContainer(self, X,Y,Z, value=1):
- steps = [timer()]
- a = value * numpy.ones((X, Y, Z), dtype='float32')
- steps.append(timer())
- t0 = dt(steps)
- #print("a refcount " , sys.getrefcount(a))
- ds = DataContainer(a, False, ['X', 'Y', 'Z'])
- return ds
-
- def test_creation_nocopy(self):
- shape = (2, 3, 4, 5)
- size = shape[0]
- for i in range(1, len(shape)):
- size = size * shape[i]
- #print("a refcount " , sys.getrefcount(a))
- a = numpy.asarray([i for i in range(size)])
- #print("a refcount " , sys.getrefcount(a))
- a = numpy.reshape(a, shape)
- #print("a refcount " , sys.getrefcount(a))
- ds = DataContainer(a, False, ['X', 'Y', 'Z', 'W'])
- #print("a refcount " , sys.getrefcount(a))
- self.assertEqual(sys.getrefcount(a), 3)
- self.assertEqual(ds.dimension_labels, {0: 'X', 1: 'Y', 2: 'Z', 3: 'W'})
-
- def testGb_creation_nocopy(self):
- X, Y, Z = 512, 512, 512
- X, Y, Z = 256, 512, 512
- steps = [timer()]
- a = numpy.ones((X, Y, Z), dtype='float32')
- steps.append(timer())
- t0 = dt(steps)
- print("test clone")
- #print("a refcount " , sys.getrefcount(a))
- ds = DataContainer(a, False, ['X', 'Y', 'Z'])
- #print("a refcount " , sys.getrefcount(a))
- self.assertEqual(sys.getrefcount(a), 3)
- ds1 = ds.copy()
- self.assertNotEqual(aid(ds.as_array()), aid(ds1.as_array()))
- ds1 = ds.clone()
- self.assertNotEqual(aid(ds.as_array()), aid(ds1.as_array()))
-
- def testInlineAlgebra(self):
- print("Test Inline Algebra")
- X, Y, Z = 1024, 512, 512
- X, Y, Z = 256, 512, 512
- steps = [timer()]
- a = numpy.ones((X, Y, Z), dtype='float32')
- steps.append(timer())
- t0 = dt(steps)
- print(t0)
- #print("a refcount " , sys.getrefcount(a))
- ds = DataContainer(a, False, ['X', 'Y', 'Z'])
- #ds.__iadd__( 2 )
- ds += 2
- steps.append(timer())
- print(dt(steps))
- self.assertEqual(ds.as_array()[0][0][0], 3.)
- #ds.__isub__( 2 )
- ds -= 2
- steps.append(timer())
- print(dt(steps))
- self.assertEqual(ds.as_array()[0][0][0], 1.)
- #ds.__imul__( 2 )
- ds *= 2
- steps.append(timer())
- print(dt(steps))
- self.assertEqual(ds.as_array()[0][0][0], 2.)
- #ds.__idiv__( 2 )
- ds /= 2
- steps.append(timer())
- print(dt(steps))
- self.assertEqual(ds.as_array()[0][0][0], 1.)
-
- ds1 = ds.copy()
- #ds1.__iadd__( 1 )
- ds1 += 1
- #ds.__iadd__( ds1 )
- ds += ds1
- steps.append(timer())
- print(dt(steps))
- self.assertEqual(ds.as_array()[0][0][0], 3.)
- #ds.__isub__( ds1 )
- ds -= ds1
- steps.append(timer())
- print(dt(steps))
- self.assertEqual(ds.as_array()[0][0][0], 1.)
- #ds.__imul__( ds1 )
- ds *= ds1
- steps.append(timer())
- print(dt(steps))
- self.assertEqual(ds.as_array()[0][0][0], 2.)
- #ds.__idiv__( ds1 )
- ds /= ds1
- steps.append(timer())
- print(dt(steps))
- self.assertEqual(ds.as_array()[0][0][0], 1.)
-
- def test_unary_operations(self):
- print("Test unary operations")
- X, Y, Z = 1024, 512, 512
- X, Y, Z = 256, 512, 512
- steps = [timer()]
- a = -numpy.ones((X, Y, Z), dtype='float32')
- steps.append(timer())
- t0 = dt(steps)
- print(t0)
- #print("a refcount " , sys.getrefcount(a))
- ds = DataContainer(a, False, ['X', 'Y', 'Z'])
-
- ds.sign(out=ds)
- steps.append(timer())
- print(dt(steps))
- self.assertEqual(ds.as_array()[0][0][0], -1.)
-
- ds.abs(out=ds)
- steps.append(timer())
- print(dt(steps))
- self.assertEqual(ds.as_array()[0][0][0], 1.)
-
- ds.__imul__(2)
- ds.sqrt(out=ds)
- steps.append(timer())
- print(dt(steps))
- self.assertEqual(ds.as_array()[0][0][0],
- numpy.sqrt(2., dtype='float32'))
-
- def test_binary_operations(self):
- self.binary_add()
- self.binary_subtract()
- self.binary_multiply()
- self.binary_divide()
-
- def binary_add(self):
- print("Test binary add")
- X, Y, Z = 512, 512, 512
- X, Y, Z = 256, 512, 512
- steps = [timer()]
- a = numpy.ones((X, Y, Z), dtype='float32')
- steps.append(timer())
- t0 = dt(steps)
-
- #print("a refcount " , sys.getrefcount(a))
- ds = DataContainer(a, False, ['X', 'Y', 'Z'])
- ds1 = ds.copy()
-
- steps.append(timer())
- ds.add(ds1, out=ds)
- steps.append(timer())
- t1 = dt(steps)
- print("ds.add(ds1, out=ds)", dt(steps))
- steps.append(timer())
- ds2 = ds.add(ds1)
- steps.append(timer())
- t2 = dt(steps)
- print("ds2 = ds.add(ds1)", dt(steps))
-
- self.assertLess(t1, t2)
- self.assertEqual(ds.as_array()[0][0][0], 2.)
-
- ds0 = ds
- ds0.add(2, out=ds0)
- steps.append(timer())
- print("ds0.add(2,out=ds0)", dt(steps), 3, ds0.as_array()[0][0][0])
- self.assertEqual(4., ds0.as_array()[0][0][0])
-
- dt1 = dt(steps)
- ds3 = ds0.add(2)
- steps.append(timer())
- print("ds3 = ds0.add(2)", dt(steps), 5, ds3.as_array()[0][0][0])
- dt2 = dt(steps)
- self.assertLess(dt1, dt2)
-
- def binary_subtract(self):
- print("Test binary subtract")
- X, Y, Z = 512, 512, 512
- steps = [timer()]
- a = numpy.ones((X, Y, Z), dtype='float32')
- steps.append(timer())
- t0 = dt(steps)
-
- #print("a refcount " , sys.getrefcount(a))
- ds = DataContainer(a, False, ['X', 'Y', 'Z'])
- ds1 = ds.copy()
-
- steps.append(timer())
- ds.subtract(ds1, out=ds)
- steps.append(timer())
- t1 = dt(steps)
- print("ds.subtract(ds1, out=ds)", dt(steps))
- self.assertEqual(0., ds.as_array()[0][0][0])
-
- steps.append(timer())
- ds2 = ds.subtract(ds1)
- self.assertEqual(-1., ds2.as_array()[0][0][0])
-
- steps.append(timer())
- t2 = dt(steps)
- print("ds2 = ds.subtract(ds1)", dt(steps))
-
- self.assertLess(t1, t2)
-
- del ds1
- ds0 = ds.copy()
- steps.append(timer())
- ds0.subtract(2, out=ds0)
- #ds0.__isub__( 2 )
- steps.append(timer())
- print("ds0.subtract(2,out=ds0)", dt(
- steps), -2., ds0.as_array()[0][0][0])
- self.assertEqual(-2., ds0.as_array()[0][0][0])
-
- dt1 = dt(steps)
- ds3 = ds0.subtract(2)
- steps.append(timer())
- print("ds3 = ds0.subtract(2)", dt(steps), 0., ds3.as_array()[0][0][0])
- dt2 = dt(steps)
- self.assertLess(dt1, dt2)
- self.assertEqual(-2., ds0.as_array()[0][0][0])
- self.assertEqual(-4., ds3.as_array()[0][0][0])
-
- def binary_multiply(self):
- print("Test binary multiply")
- X, Y, Z = 1024, 512, 512
- X, Y, Z = 256, 512, 512
- steps = [timer()]
- a = numpy.ones((X, Y, Z), dtype='float32')
- steps.append(timer())
- t0 = dt(steps)
-
- #print("a refcount " , sys.getrefcount(a))
- ds = DataContainer(a, False, ['X', 'Y', 'Z'])
- ds1 = ds.copy()
-
- steps.append(timer())
- ds.multiply(ds1, out=ds)
- steps.append(timer())
- t1 = dt(steps)
- print("ds.multiply(ds1, out=ds)", dt(steps))
- steps.append(timer())
- ds2 = ds.multiply(ds1)
- steps.append(timer())
- t2 = dt(steps)
- print("ds2 = ds.multiply(ds1)", dt(steps))
-
- self.assertLess(t1, t2)
-
- ds0 = ds
- ds0.multiply(2, out=ds0)
- steps.append(timer())
- print("ds0.multiply(2,out=ds0)", dt(
- steps), 2., ds0.as_array()[0][0][0])
- self.assertEqual(2., ds0.as_array()[0][0][0])
-
- dt1 = dt(steps)
- ds3 = ds0.multiply(2)
- steps.append(timer())
- print("ds3 = ds0.multiply(2)", dt(steps), 4., ds3.as_array()[0][0][0])
- dt2 = dt(steps)
- self.assertLess(dt1, dt2)
- self.assertEqual(4., ds3.as_array()[0][0][0])
- self.assertEqual(2., ds.as_array()[0][0][0])
-
- def binary_divide(self):
- print("Test binary divide")
- X, Y, Z = 1024, 512, 512
- X, Y, Z = 256, 512, 512
- steps = [timer()]
- a = numpy.ones((X, Y, Z), dtype='float32')
- steps.append(timer())
- t0 = dt(steps)
-
- #print("a refcount " , sys.getrefcount(a))
- ds = DataContainer(a, False, ['X', 'Y', 'Z'])
- ds1 = ds.copy()
-
- steps.append(timer())
- ds.divide(ds1, out=ds)
- steps.append(timer())
- t1 = dt(steps)
- print("ds.divide(ds1, out=ds)", dt(steps))
- steps.append(timer())
- ds2 = ds.divide(ds1)
- steps.append(timer())
- t2 = dt(steps)
- print("ds2 = ds.divide(ds1)", dt(steps))
-
- self.assertLess(t1, t2)
- self.assertEqual(ds.as_array()[0][0][0], 1.)
-
- ds0 = ds
- ds0.divide(2, out=ds0)
- steps.append(timer())
- print("ds0.divide(2,out=ds0)", dt(steps), 0.5, ds0.as_array()[0][0][0])
- self.assertEqual(0.5, ds0.as_array()[0][0][0])
-
- dt1 = dt(steps)
- ds3 = ds0.divide(2)
- steps.append(timer())
- print("ds3 = ds0.divide(2)", dt(steps), 0.25, ds3.as_array()[0][0][0])
- dt2 = dt(steps)
- self.assertLess(dt1, dt2)
- self.assertEqual(.25, ds3.as_array()[0][0][0])
- self.assertEqual(.5, ds.as_array()[0][0][0])
-
- def test_creation_copy(self):
- shape = (2, 3, 4, 5)
- size = shape[0]
- for i in range(1, len(shape)):
- size = size * shape[i]
- #print("a refcount " , sys.getrefcount(a))
- a = numpy.asarray([i for i in range(size)])
- #print("a refcount " , sys.getrefcount(a))
- a = numpy.reshape(a, shape)
- #print("a refcount " , sys.getrefcount(a))
- ds = DataContainer(a, True, ['X', 'Y', 'Z', 'W'])
- #print("a refcount " , sys.getrefcount(a))
- self.assertEqual(sys.getrefcount(a), 2)
-
- def test_subset(self):
- shape = (4, 3, 2)
- a = [i for i in range(2*3*4)]
- a = numpy.asarray(a)
- a = numpy.reshape(a, shape)
- ds = DataContainer(a, True, ['X', 'Y', 'Z'])
- sub = ds.subset(['X'])
- res = True
- try:
- numpy.testing.assert_array_equal(sub.as_array(),
- numpy.asarray([0, 6, 12, 18]))
- except AssertionError as err:
- res = False
- print(err)
- self.assertTrue(res)
-
- sub = ds.subset(['X'], Y=2, Z=0)
- res = True
- try:
- numpy.testing.assert_array_equal(sub.as_array(),
- numpy.asarray([4, 10, 16, 22]))
- except AssertionError as err:
- res = False
- print(err)
- self.assertTrue(res)
-
- sub = ds.subset(['Y'])
- try:
- numpy.testing.assert_array_equal(
- sub.as_array(), numpy.asarray([0, 2, 4]))
- res = True
- except AssertionError as err:
- res = False
- print(err)
- self.assertTrue(res)
-
- sub = ds.subset(['Z'])
- try:
- numpy.testing.assert_array_equal(
- sub.as_array(), numpy.asarray([0, 1]))
- res = True
- except AssertionError as err:
- res = False
- print(err)
- self.assertTrue(res)
- sub = ds.subset(['Z'], X=1, Y=2)
- try:
- numpy.testing.assert_array_equal(
- sub.as_array(), numpy.asarray([10, 11]))
- res = True
- except AssertionError as err:
- res = False
- print(err)
- self.assertTrue(res)
-
- print(a)
- sub = ds.subset(['X', 'Y'], Z=1)
- res = True
- try:
- numpy.testing.assert_array_equal(sub.as_array(),
- numpy.asarray([[1, 3, 5],
- [7, 9, 11],
- [13, 15, 17],
- [19, 21, 23]]))
- except AssertionError as err:
- res = False
- print(err)
- self.assertTrue(res)
-
- def test_ImageData(self):
- # create ImageData from geometry
- vgeometry = ImageGeometry(voxel_num_x=4, voxel_num_y=3, channels=2)
- vol = ImageData(geometry=vgeometry)
- self.assertEqual(vol.shape, (2, 3, 4))
-
- vol1 = vol + 1
- self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape))
-
- vol1 = vol - 1
- self.assertNumpyArrayEqual(vol1.as_array(), -numpy.ones(vol.shape))
-
- vol1 = 2 * (vol + 1)
- self.assertNumpyArrayEqual(vol1.as_array(), 2 * numpy.ones(vol.shape))
-
- vol1 = (vol + 1) / 2
- self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape) / 2)
-
- vol1 = vol + 1
- self.assertEqual(vol1.sum(), 2*3*4)
- vol1 = (vol + 2) ** 2
- self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape) * 4)
-
- self.assertEqual(vol.number_of_dimensions, 3)
-
- def test_AcquisitionData(self):
- sgeometry = AcquisitionGeometry(dimension=2, angles=numpy.linspace(0, 180, num=10),
- geom_type='parallel', pixel_num_v=3,
- pixel_num_h=5, channels=2)
- sino = AcquisitionData(geometry=sgeometry)
- self.assertEqual(sino.shape, (2, 10, 3, 5))
-
- def assertNumpyArrayEqual(self, first, second):
- res = True
- try:
- numpy.testing.assert_array_equal(first, second)
- except AssertionError as err:
- res = False
- print(err)
- self.assertTrue(res)
-
- def assertNumpyArrayAlmostEqual(self, first, second, decimal=6):
- res = True
- try:
- numpy.testing.assert_array_almost_equal(first, second, decimal)
- except AssertionError as err:
- res = False
- print(err)
- self.assertTrue(res)
- def test_DataContainerChaining(self):
- dc = self.create_DataContainer(256,256,256,1)
-
- dc.add(9,out=dc)\
- .subtract(1,out=dc)
- self.assertEqual(1+9-1,dc.as_array().flatten()[0])
-
-
-
-
-class TestAlgorithms(unittest.TestCase):
- def assertNumpyArrayEqual(self, first, second):
- res = True
- try:
- numpy.testing.assert_array_equal(first, second)
- except AssertionError as err:
- res = False
- print(err)
- self.assertTrue(res)
-
- def assertNumpyArrayAlmostEqual(self, first, second, decimal=6):
- res = True
- try:
- numpy.testing.assert_array_almost_equal(first, second, decimal)
- except AssertionError as err:
- res = False
- print(err)
- self.assertTrue(res)
-
- def test_FISTA_cvx(self):
- if not cvx_not_installable:
- # Problem data.
- m = 30
- n = 20
- np.random.seed(1)
- Amat = np.random.randn(m, n)
- A = LinearOperatorMatrix(Amat)
- bmat = np.random.randn(m)
- bmat.shape = (bmat.shape[0], 1)
-
- # A = Identity()
- # Change n to equal to m.
-
- b = DataContainer(bmat)
-
- # Regularization parameter
- lam = 10
- opt = {'memopt': True}
- # Create object instances with the test data A and b.
- f = Norm2sq(A, b, c=0.5, memopt=True)
- g0 = ZeroFun()
-
- # Initial guess
- x_init = DataContainer(np.zeros((n, 1)))
-
- f.grad(x_init)
-
- # Run FISTA for least squares plus zero function.
- x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0, opt=opt)
-
- # Print solution and final objective/criterion value for comparison
- print("FISTA least squares plus zero function solution and objective value:")
- print(x_fista0.array)
- print(criter0[-1])
-
- # Compare to CVXPY
-
- # Construct the problem.
- x0 = Variable(n)
- objective0 = Minimize(0.5*sum_squares(Amat*x0 - bmat.T[0]))
- prob0 = Problem(objective0)
-
- # The optimal objective is returned by prob.solve().
- result0 = prob0.solve(verbose=False, solver=SCS, eps=1e-9)
-
- # The optimal solution for x is stored in x.value and optimal objective value
- # is in result as well as in objective.value
- print("CVXPY least squares plus zero function solution and objective value:")
- print(x0.value)
- print(objective0.value)
- self.assertNumpyArrayAlmostEqual(
- numpy.squeeze(x_fista0.array), x0.value, 6)
- else:
- self.assertTrue(cvx_not_installable)
-
- def test_FISTA_Norm1_cvx(self):
- if not cvx_not_installable:
- opt = {'memopt': True}
- # Problem data.
- m = 30
- n = 20
- np.random.seed(1)
- Amat = np.random.randn(m, n)
- A = LinearOperatorMatrix(Amat)
- bmat = np.random.randn(m)
- bmat.shape = (bmat.shape[0], 1)
-
- # A = Identity()
- # Change n to equal to m.
-
- b = DataContainer(bmat)
-
- # Regularization parameter
- lam = 10
- opt = {'memopt': True}
- # Create object instances with the test data A and b.
- f = Norm2sq(A, b, c=0.5, memopt=True)
- g0 = ZeroFun()
-
- # Initial guess
- x_init = DataContainer(np.zeros((n, 1)))
-
- # Create 1-norm object instance
- g1 = Norm1(lam)
-
- g1(x_init)
- g1.prox(x_init, 0.02)
-
- # Combine with least squares and solve using generic FISTA implementation
- x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1, opt=opt)
-
- # Print for comparison
- print("FISTA least squares plus 1-norm solution and objective value:")
- print(x_fista1.as_array().squeeze())
- print(criter1[-1])
-
- # Compare to CVXPY
-
- # Construct the problem.
- x1 = Variable(n)
- objective1 = Minimize(
- 0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1, 1))
- prob1 = Problem(objective1)
-
- # The optimal objective is returned by prob.solve().
- result1 = prob1.solve(verbose=False, solver=SCS, eps=1e-9)
-
- # The optimal solution for x is stored in x.value and optimal objective value
- # is in result as well as in objective.value
- print("CVXPY least squares plus 1-norm solution and objective value:")
- print(x1.value)
- print(objective1.value)
-
- self.assertNumpyArrayAlmostEqual(
- numpy.squeeze(x_fista1.array), x1.value, 6)
- else:
- self.assertTrue(cvx_not_installable)
-
- def test_FBPD_Norm1_cvx(self):
- if not cvx_not_installable:
- opt = {'memopt': True}
- # Problem data.
- m = 30
- n = 20
- np.random.seed(1)
- Amat = np.random.randn(m, n)
- A = LinearOperatorMatrix(Amat)
- bmat = np.random.randn(m)
- bmat.shape = (bmat.shape[0], 1)
-
- # A = Identity()
- # Change n to equal to m.
-
- b = DataContainer(bmat)
-
- # Regularization parameter
- lam = 10
- opt = {'memopt': True}
- # Create object instances with the test data A and b.
- f = Norm2sq(A, b, c=0.5, memopt=True)
- g0 = ZeroFun()
-
- # Initial guess
- x_init = DataContainer(np.zeros((n, 1)))
-
- # Create 1-norm object instance
- g1 = Norm1(lam)
-
- # Compare to CVXPY
-
- # Construct the problem.
- x1 = Variable(n)
- objective1 = Minimize(
- 0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1, 1))
- prob1 = Problem(objective1)
-
- # The optimal objective is returned by prob.solve().
- result1 = prob1.solve(verbose=False, solver=SCS, eps=1e-9)
-
- # The optimal solution for x is stored in x.value and optimal objective value
- # is in result as well as in objective.value
- print("CVXPY least squares plus 1-norm solution and objective value:")
- print(x1.value)
- print(objective1.value)
-
- # Now try another algorithm FBPD for same problem:
- x_fbpd1, itfbpd1, timingfbpd1, criterfbpd1 = FBPD(x_init,
- Identity(), None, f, g1)
- print(x_fbpd1)
- print(criterfbpd1[-1])
-
- self.assertNumpyArrayAlmostEqual(
- numpy.squeeze(x_fbpd1.array), x1.value, 6)
- else:
- self.assertTrue(cvx_not_installable)
- # Plot criterion curve to see both FISTA and FBPD converge to same value.
- # Note that FISTA is very efficient for 1-norm minimization so it beats
- # FBPD in this test by a lot. But FBPD can handle a larger class of problems
- # than FISTA can.
-
- # Now try 1-norm and TV denoising with FBPD, first 1-norm.
-
- # Set up phantom size NxN by creating ImageGeometry, initialising the
- # ImageData object with this geometry and empty array and finally put some
- # data into its array, and display as image.
- def test_FISTA_denoise_cvx(self):
- if not cvx_not_installable:
- opt = {'memopt': True}
- N = 64
- ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N)
- Phantom = ImageData(geometry=ig)
-
- x = Phantom.as_array()
-
- x[int(round(N/4)):int(round(3*N/4)),
- int(round(N/4)):int(round(3*N/4))] = 0.5
- x[int(round(N/8)):int(round(7*N/8)),
- int(round(3*N/8)):int(round(5*N/8))] = 1
-
- # Identity operator for denoising
- I = TomoIdentity(ig)
-
- # Data and add noise
- y = I.direct(Phantom)
- y.array = y.array + 0.1*np.random.randn(N, N)
-
- # Data fidelity term
- f_denoise = Norm2sq(I, y, c=0.5, memopt=True)
-
- # 1-norm regulariser
- lam1_denoise = 1.0
- g1_denoise = Norm1(lam1_denoise)
-
- # Initial guess
- x_init_denoise = ImageData(np.zeros((N, N)))
-
- # Combine with least squares and solve using generic FISTA implementation
- x_fista1_denoise, it1_denoise, timing1_denoise, \
- criter1_denoise = \
- FISTA(x_init_denoise, f_denoise, g1_denoise, opt=opt)
-
- print(x_fista1_denoise)
- print(criter1_denoise[-1])
-
- # Now denoise LS + 1-norm with FBPD
- x_fbpd1_denoise, itfbpd1_denoise, timingfbpd1_denoise,\
- criterfbpd1_denoise = \
- FBPD(x_init_denoise, I, None, f_denoise, g1_denoise)
- print(x_fbpd1_denoise)
- print(criterfbpd1_denoise[-1])
-
- # Compare to CVXPY
-
- # Construct the problem.
- x1_denoise = Variable(N**2, 1)
- objective1_denoise = Minimize(
- 0.5*sum_squares(x1_denoise - y.array.flatten()) + lam1_denoise*norm(x1_denoise, 1))
- prob1_denoise = Problem(objective1_denoise)
-
- # The optimal objective is returned by prob.solve().
- result1_denoise = prob1_denoise.solve(
- verbose=False, solver=SCS, eps=1e-12)
-
- # The optimal solution for x is stored in x.value and optimal objective value
- # is in result as well as in objective.value
- print("CVXPY least squares plus 1-norm solution and objective value:")
- print(x1_denoise.value)
- print(objective1_denoise.value)
- self.assertNumpyArrayAlmostEqual(
- x_fista1_denoise.array.flatten(), x1_denoise.value, 5)
-
- self.assertNumpyArrayAlmostEqual(
- x_fbpd1_denoise.array.flatten(), x1_denoise.value, 5)
- x1_cvx = x1_denoise.value
- x1_cvx.shape = (N, N)
-
- # Now TV with FBPD
- lam_tv = 0.1
- gtv = TV2D(lam_tv)
- gtv(gtv.op.direct(x_init_denoise))
-
- opt_tv = {'tol': 1e-4, 'iter': 10000}
-
- x_fbpdtv_denoise, itfbpdtv_denoise, timingfbpdtv_denoise,\
- criterfbpdtv_denoise = \
- FBPD(x_init_denoise, gtv.op, None, f_denoise, gtv, opt=opt_tv)
- print(x_fbpdtv_denoise)
- print(criterfbpdtv_denoise[-1])
-
- # Compare to CVXPY
-
- # Construct the problem.
- xtv_denoise = Variable((N, N))
- objectivetv_denoise = Minimize(
- 0.5*sum_squares(xtv_denoise - y.array) + lam_tv*tv(xtv_denoise))
- probtv_denoise = Problem(objectivetv_denoise)
-
- # The optimal objective is returned by prob.solve().
- resulttv_denoise = probtv_denoise.solve(
- verbose=False, solver=SCS, eps=1e-12)
-
- # The optimal solution for x is stored in x.value and optimal objective value
- # is in result as well as in objective.value
- print("CVXPY least squares plus 1-norm solution and objective value:")
- print(xtv_denoise.value)
- print(objectivetv_denoise.value)
-
- self.assertNumpyArrayAlmostEqual(
- x_fbpdtv_denoise.as_array(), xtv_denoise.value, 1)
-
- else:
- self.assertTrue(cvx_not_installable)
-
-
-class TestFunction(unittest.TestCase):
- def assertNumpyArrayEqual(self, first, second):
- res = True
- try:
- numpy.testing.assert_array_equal(first, second)
- except AssertionError as err:
- res = False
- print(err)
- self.assertTrue(res)
-
- def create_simple_ImageData(self):
- N = 64
- ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N)
- Phantom = ImageData(geometry=ig)
-
- x = Phantom.as_array()
-
- x[int(round(N/4)):int(round(3*N/4)),
- int(round(N/4)):int(round(3*N/4))] = 0.5
- x[int(round(N/8)):int(round(7*N/8)),
- int(round(3*N/8)):int(round(5*N/8))] = 1
-
- return (ig, Phantom)
-
- def _test_Norm2(self):
- print("test Norm2")
- opt = {'memopt': True}
- ig, Phantom = self.create_simple_ImageData()
- x = Phantom.as_array()
- print(Phantom)
- print(Phantom.as_array())
-
- norm2 = Norm2()
- v1 = norm2(x)
- v2 = norm2(Phantom)
- self.assertEqual(v1, v2)
-
- p1 = norm2.prox(Phantom, 1)
- print(p1)
- p2 = norm2.prox(x, 1)
- self.assertNumpyArrayEqual(p1.as_array(), p2)
-
- p3 = norm2.proximal(Phantom, 1)
- p4 = norm2.proximal(x, 1)
- self.assertNumpyArrayEqual(p3.as_array(), p2)
- self.assertNumpyArrayEqual(p3.as_array(), p4)
-
- out = Phantom.copy()
- p5 = norm2.proximal(Phantom, 1, out=out)
- self.assertEqual(id(p5), id(out))
- self.assertNumpyArrayEqual(p5.as_array(), p3.as_array())
-# =============================================================================
-# def test_upper(self):
-# self.assertEqual('foo'.upper(), 'FOO')
-#
-# def test_isupper(self):
-# self.assertTrue('FOO'.isupper())
-# self.assertFalse('Foo'.isupper())
-#
-# def test_split(self):
-# s = 'hello world'
-# self.assertEqual(s.split(), ['hello', 'world'])
-# # check that s.split fails when the separator is not a string
-# with self.assertRaises(TypeError):
-# s.split(2)
-# =============================================================================
-
-class TestNexusReader(unittest.TestCase):
-
- def setUp(self):
- wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs')
- self.filename = '24737_fd.nxs'
-
- def tearDown(self):
- os.remove(self.filename)
-
-
- def testGetDimensions(self):
- nr = NexusReader(self.filename)
- self.assertEqual(nr.get_sinogram_dimensions(), (135, 91, 160), "Sinogram dimensions are not correct")
-
- def testGetProjectionDimensions(self):
- nr = NexusReader(self.filename)
- self.assertEqual(nr.get_projection_dimensions(), (91,135,160), "Projection dimensions are not correct")
-
- def testLoadProjectionWithoutDimensions(self):
- nr = NexusReader(self.filename)
- projections = nr.load_projection()
- self.assertEqual(projections.shape, (91,135,160), "Loaded projection data dimensions are not correct")
-
- def testLoadProjectionWithDimensions(self):
- nr = NexusReader(self.filename)
- projections = nr.load_projection((slice(0,1), slice(0,135), slice(0,160)))
- self.assertEqual(projections.shape, (1,135,160), "Loaded projection data dimensions are not correct")
-
- def testLoadProjectionCompareSingle(self):
- nr = NexusReader(self.filename)
- projections_full = nr.load_projection()
- projections_part = nr.load_projection((slice(0,1), slice(0,135), slice(0,160)))
- numpy.testing.assert_array_equal(projections_part, projections_full[0:1,:,:])
-
- def testLoadProjectionCompareMulti(self):
- nr = NexusReader(self.filename)
- projections_full = nr.load_projection()
- projections_part = nr.load_projection((slice(0,3), slice(0,135), slice(0,160)))
- numpy.testing.assert_array_equal(projections_part, projections_full[0:3,:,:])
-
- def testLoadProjectionCompareRandom(self):
- nr = NexusReader(self.filename)
- projections_full = nr.load_projection()
- projections_part = nr.load_projection((slice(1,8), slice(5,10), slice(8,20)))
- numpy.testing.assert_array_equal(projections_part, projections_full[1:8,5:10,8:20])
-
- def testLoadProjectionCompareFull(self):
- nr = NexusReader(self.filename)
- projections_full = nr.load_projection()
- projections_part = nr.load_projection((slice(None,None), slice(None,None), slice(None,None)))
- numpy.testing.assert_array_equal(projections_part, projections_full[:,:,:])
-
- def testLoadFlatCompareFull(self):
- nr = NexusReader(self.filename)
- flats_full = nr.load_flat()
- flats_part = nr.load_flat((slice(None,None), slice(None,None), slice(None,None)))
- numpy.testing.assert_array_equal(flats_part, flats_full[:,:,:])
-
- def testLoadDarkCompareFull(self):
- nr = NexusReader(self.filename)
- darks_full = nr.load_dark()
- darks_part = nr.load_dark((slice(None,None), slice(None,None), slice(None,None)))
- numpy.testing.assert_array_equal(darks_part, darks_full[:,:,:])
-
- def testProjectionAngles(self):
- nr = NexusReader(self.filename)
- angles = nr.get_projection_angles()
- self.assertEqual(angles.shape, (91,), "Loaded projection number of angles are not correct")
-
- def test_get_acquisition_data_subset(self):
- nr = NexusReader(self.filename)
- key = nr.get_image_keys()
- sl = nr.get_acquisition_data_subset(0,10)
- data = nr.get_acquisition_data().subset(['vertical','horizontal'])
-
- self.assertTrue(sl.shape , (10,data.shape[1]))
-
-class TestCompositeDataContainer(unittest.TestCase):
-
- def test_one(self):
- from ccpi.optimisation.ops import PowerMethodNonsquare
- from ccpi.optimisation.ops import TomoIdentity
- from ccpi.optimisation.funcs import Norm2sq, Norm1
- from ccpi.framework import ImageGeometry, AcquisitionGeometry
- from ccpi.optimisation.operators.CompositeOperator import CompositeDataContainer, CompositeDataContainer
- import matplotlib.pyplot as plt
-
- ig0 = ImageGeometry(2,3,4)
- ig1 = ImageGeometry(12,42,55,32)
-
- data0 = ImageData(geometry=ig0)
- data1 = ImageData(geometry=ig1) + 1
-
- data2 = ImageData(geometry=ig0) + 2
- data3 = ImageData(geometry=ig1) + 3
-
- cp0 = CompositeDataContainer(data0,data1)
- cp1 = CompositeDataContainer(data2,data3)
- #
- a = [ (el, ot) for el,ot in zip(cp0.containers,cp1.containers)]
- print (a[0][0].shape)
- #cp2 = CompositeDataContainer(*a)
- cp2 = cp0.add(cp1)
- assert (cp2.get_item(0,0).as_array()[0][0][0] == 2.)
- assert (cp2.get_item(1,0).as_array()[0][0][0] == 4.)
-
- cp2 = cp0 + cp1
- assert (cp2.get_item(0,0).as_array()[0][0][0] == 2.)
- assert (cp2.get_item(1,0).as_array()[0][0][0] == 4.)
- cp2 = cp0 + 1
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 1. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2., decimal = 5)
- cp2 = cp0 + [1 ,2]
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 1. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 3., decimal = 5)
- cp2 += cp1
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , +3. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , +6., decimal = 5)
-
- cp2 += 1
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , +4. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , +7., decimal = 5)
-
- cp2 += [-2,-1]
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 2. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 6., decimal = 5)
-
-
- cp2 = cp0.subtract(cp1)
- assert (cp2.get_item(0,0).as_array()[0][0][0] == -2.)
- assert (cp2.get_item(1,0).as_array()[0][0][0] == -2.)
- cp2 = cp0 - cp1
- assert (cp2.get_item(0,0).as_array()[0][0][0] == -2.)
- assert (cp2.get_item(1,0).as_array()[0][0][0] == -2.)
-
- cp2 = cp0 - 1
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -1. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0, decimal = 5)
- cp2 = cp0 - [1 ,2]
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -1. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -1., decimal = 5)
-
- cp2 -= cp1
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -3. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -4., decimal = 5)
-
- cp2 -= 1
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -4. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -5., decimal = 5)
-
- cp2 -= [-2,-1]
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -2. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -4., decimal = 5)
-
-
- cp2 = cp0.multiply(cp1)
- assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.)
- assert (cp2.get_item(1,0).as_array()[0][0][0] == 3.)
- cp2 = cp0 * cp1
- assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.)
- assert (cp2.get_item(1,0).as_array()[0][0][0] == 3.)
-
- cp2 = cp0 * 2
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2, decimal = 5)
- cp2 = cp0 * [3 ,2]
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2., decimal = 5)
-
- cp2 *= cp1
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0 , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , +6., decimal = 5)
-
- cp2 *= 1
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , +6., decimal = 5)
-
- cp2 *= [-2,-1]
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -6., decimal = 5)
-
-
- cp2 = cp0.divide(cp1)
- assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1./3., decimal=4)
- cp2 = cp0/cp1
- assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1./3., decimal=4)
-
- cp2 = cp0 / 2
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0.5, decimal = 5)
- cp2 = cp0 / [3 ,2]
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0.5, decimal = 5)
-
- cp2 += 1
- cp2 /= cp1
- # TODO fix inplace division
-
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 1./2 , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 1.5/3., decimal = 5)
-
- cp2 /= 1
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0.5 , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0.5, decimal = 5)
-
- cp2 /= [-2,-1]
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -0.5/2. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -0.5, decimal = 5)
- ####
-
- cp2 = cp0.power(cp1)
- assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1., decimal=4)
- cp2 = cp0**cp1
- assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1., decimal=4)
-
- cp2 = cp0 ** 2
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0., decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 1., decimal = 5)
-
- cp2 = cp0.maximum(cp1)
- assert (cp2.get_item(0,0).as_array()[0][0][0] == cp1.get_item(0,0).as_array()[0][0][0])
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], cp2.get_item(1,0).as_array()[0][0][0], decimal=4)
-
-
- cp2 = cp0.abs()
- numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0], 0., decimal=4)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1., decimal=4)
-
- cp2 = cp0.subtract(cp1)
- s = cp2.sign()
- numpy.testing.assert_almost_equal(s.get_item(0,0).as_array()[0][0][0], -1., decimal=4)
- numpy.testing.assert_almost_equal(s.get_item(1,0).as_array()[0][0][0], -1., decimal=4)
-
- cp2 = cp0.add(cp1)
- s = cp2.sqrt()
- numpy.testing.assert_almost_equal(s.get_item(0,0).as_array()[0][0][0], numpy.sqrt(2), decimal=4)
- numpy.testing.assert_almost_equal(s.get_item(1,0).as_array()[0][0][0], numpy.sqrt(4), decimal=4)
-
- s = cp0.sum()
- numpy.testing.assert_almost_equal(s[0], 0, decimal=4)
- s0 = 1
- s1 = 1
- for i in cp0.get_item(0,0).shape:
- s0 *= i
- for i in cp0.get_item(1,0).shape:
- s1 *= i
-
- numpy.testing.assert_almost_equal(s[1], cp0.get_item(0,0).as_array()[0][0][0]*s0 +cp0.get_item(1,0).as_array()[0][0][0]*s1, decimal=4)
-
-
-if __name__ == '__main__':
- unittest.main()
-
diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py index 00eeed0..85907eb 100644 --- a/Wrappers/Python/setup.py +++ b/Wrappers/Python/setup.py @@ -23,12 +23,16 @@ import os import sys -cil_version='0.11.3' +cil_version=os.environ['CIL_VERSION'] +if cil_version == '': + print("Please set the environmental variable CIL_VERSION") + sys.exit(1) setup( name="ccpi-framework", version=cil_version, - packages=['ccpi' , 'ccpi.io', 'ccpi.optimisation', 'ccpi.optimisation.operators'], + packages=['ccpi' , 'ccpi.io', 'ccpi.optimisation', + 'ccpi.optimisation.operators', 'ccpi.optimisation.algorithms'], # Project uses reStructuredText, so ensure that the docutils get # installed or upgraded on the target machine diff --git a/Wrappers/Python/test/__init__.py b/Wrappers/Python/test/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/Wrappers/Python/test/__init__.py diff --git a/Wrappers/Python/test/TestReader.py b/Wrappers/Python/test/skip_TestReader.py index 51db052..ec7ed58 100644 --- a/Wrappers/Python/test/TestReader.py +++ b/Wrappers/Python/test/skip_TestReader.py @@ -1,134 +1,134 @@ -# -*- 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 Jakob Jorgensen, Daniil Kazantsev, Edoardo Pasca and Srikanth Nagella
-
-# 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.
-
-'''
-Unit tests for Readers
-
-@author: Mr. Srikanth Nagella
-'''
-import unittest
-
-import numpy.testing
-import wget
-import os
-from ccpi.io.reader import NexusReader
-from ccpi.io.reader import XTEKReader
-#@unittest.skip
-class TestNexusReader(unittest.TestCase):
-
- def setUp(self):
- wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs')
- self.filename = '24737_fd.nxs'
-
- def tearDown(self):
- os.remove(self.filename)
-
-
- def testGetDimensions(self):
- nr = NexusReader(self.filename)
- self.assertEqual(nr.getSinogramDimensions(), (135, 91, 160), "Sinogram dimensions are not correct")
-
- def testGetProjectionDimensions(self):
- nr = NexusReader(self.filename)
- self.assertEqual(nr.getProjectionDimensions(), (91,135,160), "Projection dimensions are not correct")
-
- def testLoadProjectionWithoutDimensions(self):
- nr = NexusReader(self.filename)
- projections = nr.loadProjection()
- self.assertEqual(projections.shape, (91,135,160), "Loaded projection data dimensions are not correct")
-
- def testLoadProjectionWithDimensions(self):
- nr = NexusReader(self.filename)
- projections = nr.loadProjection((slice(0,1), slice(0,135), slice(0,160)))
- self.assertEqual(projections.shape, (1,135,160), "Loaded projection data dimensions are not correct")
-
- def testLoadProjectionCompareSingle(self):
- nr = NexusReader(self.filename)
- projections_full = nr.loadProjection()
- projections_part = nr.loadProjection((slice(0,1), slice(0,135), slice(0,160)))
- numpy.testing.assert_array_equal(projections_part, projections_full[0:1,:,:])
-
- def testLoadProjectionCompareMulti(self):
- nr = NexusReader(self.filename)
- projections_full = nr.loadProjection()
- projections_part = nr.loadProjection((slice(0,3), slice(0,135), slice(0,160)))
- numpy.testing.assert_array_equal(projections_part, projections_full[0:3,:,:])
-
- def testLoadProjectionCompareRandom(self):
- nr = NexusReader(self.filename)
- projections_full = nr.loadProjection()
- projections_part = nr.loadProjection((slice(1,8), slice(5,10), slice(8,20)))
- numpy.testing.assert_array_equal(projections_part, projections_full[1:8,5:10,8:20])
-
- def testLoadProjectionCompareFull(self):
- nr = NexusReader(self.filename)
- projections_full = nr.loadProjection()
- projections_part = nr.loadProjection((slice(None,None), slice(None,None), slice(None,None)))
- numpy.testing.assert_array_equal(projections_part, projections_full[:,:,:])
-
- def testLoadFlatCompareFull(self):
- nr = NexusReader(self.filename)
- flats_full = nr.loadFlat()
- flats_part = nr.loadFlat((slice(None,None), slice(None,None), slice(None,None)))
- numpy.testing.assert_array_equal(flats_part, flats_full[:,:,:])
-
- def testLoadDarkCompareFull(self):
- nr = NexusReader(self.filename)
- darks_full = nr.loadDark()
- darks_part = nr.loadDark((slice(None,None), slice(None,None), slice(None,None)))
- numpy.testing.assert_array_equal(darks_part, darks_full[:,:,:])
-
- def testProjectionAngles(self):
- nr = NexusReader(self.filename)
- angles = nr.getProjectionAngles()
- self.assertEqual(angles.shape, (91,), "Loaded projection number of angles are not correct")
-
-class TestXTEKReader(unittest.TestCase):
-
- def setUp(self):
- testpath, filename = os.path.split(os.path.realpath(__file__))
- testpath = os.path.normpath(os.path.join(testpath, '..','..','..'))
- self.filename = os.path.join(testpath,'data','SophiaBeads','SophiaBeads_64_averaged.xtekct')
-
- def tearDown(self):
- pass
-
- def testLoad(self):
- xtekReader = XTEKReader(self.filename)
- self.assertEqual(xtekReader.geometry.pixel_num_h, 500, "Detector pixel X size is not correct")
- self.assertEqual(xtekReader.geometry.pixel_num_v, 500, "Detector pixel Y size is not correct")
- self.assertEqual(xtekReader.geometry.dist_source_center, -80.6392412185669, "Distance from source to center is not correct")
- self.assertEqual(xtekReader.geometry.dist_center_detector, (1007.006 - 80.6392412185669), "Distance from center to detector is not correct")
-
- def testReadAngles(self):
- xtekReader = XTEKReader(self.filename)
- angles = xtekReader.readAngles()
- self.assertEqual(angles.shape, (63,), "Angles doesn't match")
- self.assertAlmostEqual(angles[46], -085.717, 3, "46th Angles doesn't match")
-
- def testLoadProjection(self):
- xtekReader = XTEKReader(self.filename)
- pixels = xtekReader.loadProjection()
- self.assertEqual(pixels.shape, (63, 500, 500), "projections shape doesn't match")
-
-
-
-if __name__ == "__main__":
- #import sys;sys.argv = ['', 'TestNexusReader.testLoad']
+# -*- 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 Jakob Jorgensen, Daniil Kazantsev, Edoardo Pasca and Srikanth Nagella + +# 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. + +''' +Unit tests for Readers + +@author: Mr. Srikanth Nagella +''' +import unittest + +import numpy.testing +import wget +import os +from ccpi.io.reader import NexusReader +from ccpi.io.reader import XTEKReader +#@unittest.skip +class TestNexusReader(unittest.TestCase): + + def setUp(self): + wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs') + self.filename = '24737_fd.nxs' + + def tearDown(self): + os.remove(self.filename) + + + def testGetDimensions(self): + nr = NexusReader(self.filename) + self.assertEqual(nr.getSinogramDimensions(), (135, 91, 160), "Sinogram dimensions are not correct") + + def testGetProjectionDimensions(self): + nr = NexusReader(self.filename) + self.assertEqual(nr.getProjectionDimensions(), (91,135,160), "Projection dimensions are not correct") + + def testLoadProjectionWithoutDimensions(self): + nr = NexusReader(self.filename) + projections = nr.loadProjection() + self.assertEqual(projections.shape, (91,135,160), "Loaded projection data dimensions are not correct") + + def testLoadProjectionWithDimensions(self): + nr = NexusReader(self.filename) + projections = nr.loadProjection((slice(0,1), slice(0,135), slice(0,160))) + self.assertEqual(projections.shape, (1,135,160), "Loaded projection data dimensions are not correct") + + def testLoadProjectionCompareSingle(self): + nr = NexusReader(self.filename) + projections_full = nr.loadProjection() + projections_part = nr.loadProjection((slice(0,1), slice(0,135), slice(0,160))) + numpy.testing.assert_array_equal(projections_part, projections_full[0:1,:,:]) + + def testLoadProjectionCompareMulti(self): + nr = NexusReader(self.filename) + projections_full = nr.loadProjection() + projections_part = nr.loadProjection((slice(0,3), slice(0,135), slice(0,160))) + numpy.testing.assert_array_equal(projections_part, projections_full[0:3,:,:]) + + def testLoadProjectionCompareRandom(self): + nr = NexusReader(self.filename) + projections_full = nr.loadProjection() + projections_part = nr.loadProjection((slice(1,8), slice(5,10), slice(8,20))) + numpy.testing.assert_array_equal(projections_part, projections_full[1:8,5:10,8:20]) + + def testLoadProjectionCompareFull(self): + nr = NexusReader(self.filename) + projections_full = nr.loadProjection() + projections_part = nr.loadProjection((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(projections_part, projections_full[:,:,:]) + + def testLoadFlatCompareFull(self): + nr = NexusReader(self.filename) + flats_full = nr.loadFlat() + flats_part = nr.loadFlat((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(flats_part, flats_full[:,:,:]) + + def testLoadDarkCompareFull(self): + nr = NexusReader(self.filename) + darks_full = nr.loadDark() + darks_part = nr.loadDark((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(darks_part, darks_full[:,:,:]) + + def testProjectionAngles(self): + nr = NexusReader(self.filename) + angles = nr.getProjectionAngles() + self.assertEqual(angles.shape, (91,), "Loaded projection number of angles are not correct") + +class TestXTEKReader(unittest.TestCase): + + def setUp(self): + testpath, filename = os.path.split(os.path.realpath(__file__)) + testpath = os.path.normpath(os.path.join(testpath, '..','..','..')) + self.filename = os.path.join(testpath,'data','SophiaBeads','SophiaBeads_64_averaged.xtekct') + + def tearDown(self): + pass + + def testLoad(self): + xtekReader = XTEKReader(self.filename) + self.assertEqual(xtekReader.geometry.pixel_num_h, 500, "Detector pixel X size is not correct") + self.assertEqual(xtekReader.geometry.pixel_num_v, 500, "Detector pixel Y size is not correct") + self.assertEqual(xtekReader.geometry.dist_source_center, -80.6392412185669, "Distance from source to center is not correct") + self.assertEqual(xtekReader.geometry.dist_center_detector, (1007.006 - 80.6392412185669), "Distance from center to detector is not correct") + + def testReadAngles(self): + xtekReader = XTEKReader(self.filename) + angles = xtekReader.readAngles() + self.assertEqual(angles.shape, (63,), "Angles doesn't match") + self.assertAlmostEqual(angles[46], -085.717, 3, "46th Angles doesn't match") + + def testLoadProjection(self): + xtekReader = XTEKReader(self.filename) + pixels = xtekReader.loadProjection() + self.assertEqual(pixels.shape, (63, 500, 500), "projections shape doesn't match") + + + +if __name__ == "__main__": + #import sys;sys.argv = ['', 'TestNexusReader.testLoad'] unittest.main()
\ No newline at end of file diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py new file mode 100755 index 0000000..05f3fe8 --- /dev/null +++ b/Wrappers/Python/test/test_DataContainer.py @@ -0,0 +1,499 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Feb 27 15:08:21 2019 + +@author: ofn77899 +""" +import sys +import unittest +import numpy +from ccpi.framework import DataContainer +from ccpi.framework import ImageData +from ccpi.framework import AcquisitionData +from ccpi.framework import ImageGeometry +from ccpi.framework import AcquisitionGeometry +from timeit import default_timer as timer + + +def dt(steps): + return steps[-1] - steps[-2] +def aid(x): + # This function returns the memory + # block address of an array. + return x.__array_interface__['data'][0] + + + +class TestDataContainer(unittest.TestCase): + def create_simple_ImageData(self): + N = 64 + ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) + Phantom = ImageData(geometry=ig) + + x = Phantom.as_array() + + x[int(round(N/4)):int(round(3*N/4)), + int(round(N/4)):int(round(3*N/4))] = 0.5 + x[int(round(N/8)):int(round(7*N/8)), + int(round(3*N/8)):int(round(5*N/8))] = 1 + + return (ig, Phantom) + + def create_DataContainer(self, X,Y,Z, value=1): + steps = [timer()] + a = value * numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + return ds + + def test_creation_nocopy(self): + shape = (2, 3, 4, 5) + size = shape[0] + for i in range(1, len(shape)): + size = size * shape[i] + #print("a refcount " , sys.getrefcount(a)) + a = numpy.asarray([i for i in range(size)]) + #print("a refcount " , sys.getrefcount(a)) + a = numpy.reshape(a, shape) + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z', 'W']) + #print("a refcount " , sys.getrefcount(a)) + self.assertEqual(sys.getrefcount(a), 3) + self.assertEqual(ds.dimension_labels, {0: 'X', 1: 'Y', 2: 'Z', 3: 'W'}) + + def testGb_creation_nocopy(self): + X, Y, Z = 512, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + print("test clone") + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + #print("a refcount " , sys.getrefcount(a)) + self.assertEqual(sys.getrefcount(a), 3) + ds1 = ds.copy() + self.assertNotEqual(aid(ds.as_array()), aid(ds1.as_array())) + ds1 = ds.clone() + self.assertNotEqual(aid(ds.as_array()), aid(ds1.as_array())) + + def testInlineAlgebra(self): + print("Test Inline Algebra") + X, Y, Z = 1024, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + print(t0) + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + #ds.__iadd__( 2 ) + ds += 2 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 3.) + #ds.__isub__( 2 ) + ds -= 2 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 1.) + #ds.__imul__( 2 ) + ds *= 2 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 2.) + #ds.__idiv__( 2 ) + ds /= 2 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 1.) + + ds1 = ds.copy() + #ds1.__iadd__( 1 ) + ds1 += 1 + #ds.__iadd__( ds1 ) + ds += ds1 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 3.) + #ds.__isub__( ds1 ) + ds -= ds1 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 1.) + #ds.__imul__( ds1 ) + ds *= ds1 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 2.) + #ds.__idiv__( ds1 ) + ds /= ds1 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 1.) + + def test_unary_operations(self): + print("Test unary operations") + X, Y, Z = 1024, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = -numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + print(t0) + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + + ds.sign(out=ds) + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], -1.) + + ds.abs(out=ds) + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 1.) + + ds.__imul__(2) + ds.sqrt(out=ds) + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], + numpy.sqrt(2., dtype='float32')) + + def test_binary_operations(self): + self.binary_add() + self.binary_subtract() + self.binary_multiply() + self.binary_divide() + + def binary_add(self): + print("Test binary add") + X, Y, Z = 512, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + ds1 = ds.copy() + + steps.append(timer()) + ds.add(ds1, out=ds) + steps.append(timer()) + t1 = dt(steps) + print("ds.add(ds1, out=ds)", dt(steps)) + steps.append(timer()) + ds2 = ds.add(ds1) + steps.append(timer()) + t2 = dt(steps) + print("ds2 = ds.add(ds1)", dt(steps)) + + self.assertLess(t1, t2) + self.assertEqual(ds.as_array()[0][0][0], 2.) + + ds0 = ds + ds0.add(2, out=ds0) + steps.append(timer()) + print("ds0.add(2,out=ds0)", dt(steps), 3, ds0.as_array()[0][0][0]) + self.assertEqual(4., ds0.as_array()[0][0][0]) + + dt1 = dt(steps) + ds3 = ds0.add(2) + steps.append(timer()) + print("ds3 = ds0.add(2)", dt(steps), 5, ds3.as_array()[0][0][0]) + dt2 = dt(steps) + self.assertLess(dt1, dt2) + + def binary_subtract(self): + print("Test binary subtract") + X, Y, Z = 512, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + ds1 = ds.copy() + + steps.append(timer()) + ds.subtract(ds1, out=ds) + steps.append(timer()) + t1 = dt(steps) + print("ds.subtract(ds1, out=ds)", dt(steps)) + self.assertEqual(0., ds.as_array()[0][0][0]) + + steps.append(timer()) + ds2 = ds.subtract(ds1) + self.assertEqual(-1., ds2.as_array()[0][0][0]) + + steps.append(timer()) + t2 = dt(steps) + print("ds2 = ds.subtract(ds1)", dt(steps)) + + self.assertLess(t1, t2) + + del ds1 + ds0 = ds.copy() + steps.append(timer()) + ds0.subtract(2, out=ds0) + #ds0.__isub__( 2 ) + steps.append(timer()) + print("ds0.subtract(2,out=ds0)", dt( + steps), -2., ds0.as_array()[0][0][0]) + self.assertEqual(-2., ds0.as_array()[0][0][0]) + + dt1 = dt(steps) + ds3 = ds0.subtract(2) + steps.append(timer()) + print("ds3 = ds0.subtract(2)", dt(steps), 0., ds3.as_array()[0][0][0]) + dt2 = dt(steps) + self.assertLess(dt1, dt2) + self.assertEqual(-2., ds0.as_array()[0][0][0]) + self.assertEqual(-4., ds3.as_array()[0][0][0]) + + def binary_multiply(self): + print("Test binary multiply") + X, Y, Z = 1024, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + ds1 = ds.copy() + + steps.append(timer()) + ds.multiply(ds1, out=ds) + steps.append(timer()) + t1 = dt(steps) + print("ds.multiply(ds1, out=ds)", dt(steps)) + steps.append(timer()) + ds2 = ds.multiply(ds1) + steps.append(timer()) + t2 = dt(steps) + print("ds2 = ds.multiply(ds1)", dt(steps)) + + self.assertLess(t1, t2) + + ds0 = ds + ds0.multiply(2, out=ds0) + steps.append(timer()) + print("ds0.multiply(2,out=ds0)", dt( + steps), 2., ds0.as_array()[0][0][0]) + self.assertEqual(2., ds0.as_array()[0][0][0]) + + dt1 = dt(steps) + ds3 = ds0.multiply(2) + steps.append(timer()) + print("ds3 = ds0.multiply(2)", dt(steps), 4., ds3.as_array()[0][0][0]) + dt2 = dt(steps) + self.assertLess(dt1, dt2) + self.assertEqual(4., ds3.as_array()[0][0][0]) + self.assertEqual(2., ds.as_array()[0][0][0]) + + def binary_divide(self): + print("Test binary divide") + X, Y, Z = 1024, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + ds1 = ds.copy() + + steps.append(timer()) + ds.divide(ds1, out=ds) + steps.append(timer()) + t1 = dt(steps) + print("ds.divide(ds1, out=ds)", dt(steps)) + steps.append(timer()) + ds2 = ds.divide(ds1) + steps.append(timer()) + t2 = dt(steps) + print("ds2 = ds.divide(ds1)", dt(steps)) + + self.assertLess(t1, t2) + self.assertEqual(ds.as_array()[0][0][0], 1.) + + ds0 = ds + ds0.divide(2, out=ds0) + steps.append(timer()) + print("ds0.divide(2,out=ds0)", dt(steps), 0.5, ds0.as_array()[0][0][0]) + self.assertEqual(0.5, ds0.as_array()[0][0][0]) + + dt1 = dt(steps) + ds3 = ds0.divide(2) + steps.append(timer()) + print("ds3 = ds0.divide(2)", dt(steps), 0.25, ds3.as_array()[0][0][0]) + dt2 = dt(steps) + self.assertLess(dt1, dt2) + self.assertEqual(.25, ds3.as_array()[0][0][0]) + self.assertEqual(.5, ds.as_array()[0][0][0]) + + def test_creation_copy(self): + shape = (2, 3, 4, 5) + size = shape[0] + for i in range(1, len(shape)): + size = size * shape[i] + #print("a refcount " , sys.getrefcount(a)) + a = numpy.asarray([i for i in range(size)]) + #print("a refcount " , sys.getrefcount(a)) + a = numpy.reshape(a, shape) + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, True, ['X', 'Y', 'Z', 'W']) + #print("a refcount " , sys.getrefcount(a)) + self.assertEqual(sys.getrefcount(a), 2) + + def test_subset(self): + shape = (4, 3, 2) + a = [i for i in range(2*3*4)] + a = numpy.asarray(a) + a = numpy.reshape(a, shape) + ds = DataContainer(a, True, ['X', 'Y', 'Z']) + sub = ds.subset(['X']) + res = True + try: + numpy.testing.assert_array_equal(sub.as_array(), + numpy.asarray([0, 6, 12, 18])) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + sub = ds.subset(['X'], Y=2, Z=0) + res = True + try: + numpy.testing.assert_array_equal(sub.as_array(), + numpy.asarray([4, 10, 16, 22])) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + sub = ds.subset(['Y']) + try: + numpy.testing.assert_array_equal( + sub.as_array(), numpy.asarray([0, 2, 4])) + res = True + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + sub = ds.subset(['Z']) + try: + numpy.testing.assert_array_equal( + sub.as_array(), numpy.asarray([0, 1])) + res = True + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + sub = ds.subset(['Z'], X=1, Y=2) + try: + numpy.testing.assert_array_equal( + sub.as_array(), numpy.asarray([10, 11])) + res = True + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + print(a) + sub = ds.subset(['X', 'Y'], Z=1) + res = True + try: + numpy.testing.assert_array_equal(sub.as_array(), + numpy.asarray([[1, 3, 5], + [7, 9, 11], + [13, 15, 17], + [19, 21, 23]])) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def test_ImageData(self): + # create ImageData from geometry + vgeometry = ImageGeometry(voxel_num_x=4, voxel_num_y=3, channels=2) + vol = ImageData(geometry=vgeometry) + self.assertEqual(vol.shape, (2, 3, 4)) + + vol1 = vol + 1 + self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape)) + + vol1 = vol - 1 + self.assertNumpyArrayEqual(vol1.as_array(), -numpy.ones(vol.shape)) + + vol1 = 2 * (vol + 1) + self.assertNumpyArrayEqual(vol1.as_array(), 2 * numpy.ones(vol.shape)) + + vol1 = (vol + 1) / 2 + self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape) / 2) + + vol1 = vol + 1 + self.assertEqual(vol1.sum(), 2*3*4) + vol1 = (vol + 2) ** 2 + self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape) * 4) + + self.assertEqual(vol.number_of_dimensions, 3) + + def test_AcquisitionData(self): + sgeometry = AcquisitionGeometry(dimension=2, angles=numpy.linspace(0, 180, num=10), + geom_type='parallel', pixel_num_v=3, + pixel_num_h=5, channels=2) + sino = AcquisitionData(geometry=sgeometry) + self.assertEqual(sino.shape, (2, 10, 3, 5)) + + def assertNumpyArrayEqual(self, first, second): + res = True + try: + numpy.testing.assert_array_equal(first, second) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def assertNumpyArrayAlmostEqual(self, first, second, decimal=6): + res = True + try: + numpy.testing.assert_array_almost_equal(first, second, decimal) + except AssertionError as err: + res = False + print(err) + print("expected " , second) + print("actual " , first) + + self.assertTrue(res) + def test_DataContainerChaining(self): + dc = self.create_DataContainer(256,256,256,1) + + dc.add(9,out=dc)\ + .subtract(1,out=dc) + self.assertEqual(1+9-1,dc.as_array().flatten()[0]) + def test_reduction(self): + print ("test reductions") + dc = self.create_DataContainer(2,2,2,value=1) + sqnorm = dc.squared_norm() + norm = dc.norm() + self.assertEqual(sqnorm, 8.0) + numpy.testing.assert_almost_equal(norm, numpy.sqrt(8.0), decimal=7) + + + +if __name__ == '__main__': + unittest.main() +
\ No newline at end of file diff --git a/Wrappers/Python/test/test_NexusReader.py b/Wrappers/Python/test/test_NexusReader.py new file mode 100755 index 0000000..55543ba --- /dev/null +++ b/Wrappers/Python/test/test_NexusReader.py @@ -0,0 +1,96 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Feb 27 15:05:00 2019 + +@author: ofn77899 +""" + +import unittest +import wget +import os +from ccpi.io.reader import NexusReader +import numpy + + +class TestNexusReader(unittest.TestCase): + + def setUp(self): + wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs') + self.filename = '24737_fd.nxs' + + def tearDown(self): + os.remove(self.filename) + + + def testGetDimensions(self): + nr = NexusReader(self.filename) + self.assertEqual(nr.get_sinogram_dimensions(), (135, 91, 160), "Sinogram dimensions are not correct") + + def testGetProjectionDimensions(self): + nr = NexusReader(self.filename) + self.assertEqual(nr.get_projection_dimensions(), (91,135,160), "Projection dimensions are not correct") + + def testLoadProjectionWithoutDimensions(self): + nr = NexusReader(self.filename) + projections = nr.load_projection() + self.assertEqual(projections.shape, (91,135,160), "Loaded projection data dimensions are not correct") + + def testLoadProjectionWithDimensions(self): + nr = NexusReader(self.filename) + projections = nr.load_projection((slice(0,1), slice(0,135), slice(0,160))) + self.assertEqual(projections.shape, (1,135,160), "Loaded projection data dimensions are not correct") + + def testLoadProjectionCompareSingle(self): + nr = NexusReader(self.filename) + projections_full = nr.load_projection() + projections_part = nr.load_projection((slice(0,1), slice(0,135), slice(0,160))) + numpy.testing.assert_array_equal(projections_part, projections_full[0:1,:,:]) + + def testLoadProjectionCompareMulti(self): + nr = NexusReader(self.filename) + projections_full = nr.load_projection() + projections_part = nr.load_projection((slice(0,3), slice(0,135), slice(0,160))) + numpy.testing.assert_array_equal(projections_part, projections_full[0:3,:,:]) + + def testLoadProjectionCompareRandom(self): + nr = NexusReader(self.filename) + projections_full = nr.load_projection() + projections_part = nr.load_projection((slice(1,8), slice(5,10), slice(8,20))) + numpy.testing.assert_array_equal(projections_part, projections_full[1:8,5:10,8:20]) + + def testLoadProjectionCompareFull(self): + nr = NexusReader(self.filename) + projections_full = nr.load_projection() + projections_part = nr.load_projection((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(projections_part, projections_full[:,:,:]) + + def testLoadFlatCompareFull(self): + nr = NexusReader(self.filename) + flats_full = nr.load_flat() + flats_part = nr.load_flat((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(flats_part, flats_full[:,:,:]) + + def testLoadDarkCompareFull(self): + nr = NexusReader(self.filename) + darks_full = nr.load_dark() + darks_part = nr.load_dark((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(darks_part, darks_full[:,:,:]) + + def testProjectionAngles(self): + nr = NexusReader(self.filename) + angles = nr.get_projection_angles() + self.assertEqual(angles.shape, (91,), "Loaded projection number of angles are not correct") + + def test_get_acquisition_data_subset(self): + nr = NexusReader(self.filename) + key = nr.get_image_keys() + sl = nr.get_acquisition_data_subset(0,10) + data = nr.get_acquisition_data().subset(['vertical','horizontal']) + + self.assertTrue(sl.shape , (10,data.shape[1])) + + + +if __name__ == '__main__': + unittest.main() + diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py new file mode 100755 index 0000000..b5959b5 --- /dev/null +++ b/Wrappers/Python/test/test_algorithms.py @@ -0,0 +1,123 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Feb 27 15:11:36 2019 + +@author: ofn77899 +""" + +import unittest +import numpy +from ccpi.framework import DataContainer +from ccpi.framework import ImageData +from ccpi.framework import AcquisitionData +from ccpi.framework import ImageGeometry +from ccpi.framework import AcquisitionGeometry +from ccpi.optimisation.ops import TomoIdentity +from ccpi.optimisation.funcs import Norm2sq +from ccpi.optimisation.algorithms import GradientDescent +from ccpi.optimisation.algorithms import CGLS +from ccpi.optimisation.algorithms import FISTA +from ccpi.optimisation.algorithms import FBPD + + + + +class TestAlgorithms(unittest.TestCase): + def setUp(self): + #wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs') + #self.filename = '24737_fd.nxs' + # we use TomoIdentity as the operator and solve the simple least squares + # problem for a random-valued ImageData or AcquisitionData b? + # Then we know the minimiser is b itself + + # || I x -b ||^2 + + # create an ImageGeometry + ig = ImageGeometry(12,13,14) + pass + + def tearDown(self): + #os.remove(self.filename) + pass + + def test_GradientDescent(self): + print ("Test GradientDescent") + ig = ImageGeometry(12,13,14) + x_init = ImageData(geometry=ig) + b = x_init.copy() + # fill with random numbers + b.fill(numpy.random.random(x_init.shape)) + + identity = TomoIdentity(geometry=ig) + + norm2sq = Norm2sq(identity, b) + + alg = GradientDescent(x_init=x_init, + objective_function=norm2sq, + rate=0.3) + alg.max_iteration = 20 + alg.run(20, verbose=True) + self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) + def test_CGLS(self): + print ("Test CGLS") + ig = ImageGeometry(124,153,154) + x_init = ImageData(geometry=ig) + b = x_init.copy() + # fill with random numbers + b.fill(numpy.random.random(x_init.shape)) + + identity = TomoIdentity(geometry=ig) + + alg = CGLS(x_init=x_init, operator=identity, data=b) + alg.max_iteration = 1 + alg.run(20, verbose=True) + self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) + + def test_FISTA(self): + print ("Test FISTA") + ig = ImageGeometry(127,139,149) + x_init = ImageData(geometry=ig) + b = x_init.copy() + # fill with random numbers + b.fill(numpy.random.random(x_init.shape)) + x_init = ImageData(geometry=ig) + x_init.fill(numpy.random.random(x_init.shape)) + + identity = TomoIdentity(geometry=ig) + + norm2sq = Norm2sq(identity, b) + opt = {'tol': 1e-4, 'memopt':False} + alg = FISTA(x_init=x_init, f=norm2sq, g=None, opt=opt) + alg.max_iteration = 2 + alg.run(20, verbose=True) + self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) + alg.run(20, verbose=True) + self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) + + + + def assertNumpyArrayEqual(self, first, second): + res = True + try: + numpy.testing.assert_array_equal(first, second) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def assertNumpyArrayAlmostEqual(self, first, second, decimal=6): + res = True + try: + numpy.testing.assert_array_almost_equal(first, second, decimal) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + + + + +if __name__ == '__main__': + unittest.main() +
\ No newline at end of file diff --git a/Wrappers/Python/test/test_run_test.py b/Wrappers/Python/test/test_run_test.py new file mode 100755 index 0000000..d0b87f5 --- /dev/null +++ b/Wrappers/Python/test/test_run_test.py @@ -0,0 +1,435 @@ +import unittest +import numpy +import numpy as np +from ccpi.framework import DataContainer +from ccpi.framework import ImageData +from ccpi.framework import AcquisitionData +from ccpi.framework import ImageGeometry +from ccpi.framework import AcquisitionGeometry +from ccpi.optimisation.algs import FISTA +from ccpi.optimisation.algs import FBPD +from ccpi.optimisation.funcs import Norm2sq +from ccpi.optimisation.funcs import ZeroFun +from ccpi.optimisation.funcs import Norm1 +from ccpi.optimisation.funcs import TV2D +from ccpi.optimisation.funcs import Norm2 + +from ccpi.optimisation.ops import LinearOperatorMatrix +from ccpi.optimisation.ops import TomoIdentity +from ccpi.optimisation.ops import Identity +from ccpi.optimisation.ops import PowerMethodNonsquare + + +import numpy.testing + +try: + from cvxpy import * + cvx_not_installable = False +except ImportError: + cvx_not_installable = True + + +def aid(x): + # This function returns the memory + # block address of an array. + return x.__array_interface__['data'][0] + + +def dt(steps): + return steps[-1] - steps[-2] + + + + +class TestAlgorithms(unittest.TestCase): + def assertNumpyArrayEqual(self, first, second): + res = True + try: + numpy.testing.assert_array_equal(first, second) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def assertNumpyArrayAlmostEqual(self, first, second, decimal=6): + res = True + try: + numpy.testing.assert_array_almost_equal(first, second, decimal) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def test_FISTA_cvx(self): + if not cvx_not_installable: + # Problem data. + m = 30 + n = 20 + np.random.seed(1) + Amat = np.random.randn(m, n) + A = LinearOperatorMatrix(Amat) + bmat = np.random.randn(m) + bmat.shape = (bmat.shape[0], 1) + + # A = Identity() + # Change n to equal to m. + + b = DataContainer(bmat) + + # Regularization parameter + lam = 10 + opt = {'memopt': True} + # Create object instances with the test data A and b. + f = Norm2sq(A, b, c=0.5, memopt=True) + g0 = ZeroFun() + + # Initial guess + x_init = DataContainer(np.zeros((n, 1))) + + f.grad(x_init) + + # Run FISTA for least squares plus zero function. + x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0, opt=opt) + + # Print solution and final objective/criterion value for comparison + print("FISTA least squares plus zero function solution and objective value:") + print(x_fista0.array) + print(criter0[-1]) + + # Compare to CVXPY + + # Construct the problem. + x0 = Variable(n) + objective0 = Minimize(0.5*sum_squares(Amat*x0 - bmat.T[0])) + prob0 = Problem(objective0) + + # The optimal objective is returned by prob.solve(). + result0 = prob0.solve(verbose=False, solver=SCS, eps=1e-9) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus zero function solution and objective value:") + print(x0.value) + print(objective0.value) + self.assertNumpyArrayAlmostEqual( + numpy.squeeze(x_fista0.array), x0.value, 6) + else: + self.assertTrue(cvx_not_installable) + + def test_FISTA_Norm1_cvx(self): + if not cvx_not_installable: + opt = {'memopt': True} + # Problem data. + m = 30 + n = 20 + np.random.seed(1) + Amat = np.random.randn(m, n) + A = LinearOperatorMatrix(Amat) + bmat = np.random.randn(m) + bmat.shape = (bmat.shape[0], 1) + + # A = Identity() + # Change n to equal to m. + + b = DataContainer(bmat) + + # Regularization parameter + lam = 10 + opt = {'memopt': True} + # Create object instances with the test data A and b. + f = Norm2sq(A, b, c=0.5, memopt=True) + g0 = ZeroFun() + + # Initial guess + x_init = DataContainer(np.zeros((n, 1))) + + # Create 1-norm object instance + g1 = Norm1(lam) + + g1(x_init) + g1.prox(x_init, 0.02) + + # Combine with least squares and solve using generic FISTA implementation + x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1, opt=opt) + + # Print for comparison + print("FISTA least squares plus 1-norm solution and objective value:") + print(x_fista1.as_array().squeeze()) + print(criter1[-1]) + + # Compare to CVXPY + + # Construct the problem. + x1 = Variable(n) + objective1 = Minimize( + 0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1, 1)) + prob1 = Problem(objective1) + + # The optimal objective is returned by prob.solve(). + result1 = prob1.solve(verbose=False, solver=SCS, eps=1e-9) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus 1-norm solution and objective value:") + print(x1.value) + print(objective1.value) + + self.assertNumpyArrayAlmostEqual( + numpy.squeeze(x_fista1.array), x1.value, 6) + else: + self.assertTrue(cvx_not_installable) + + def skip_test_FBPD_Norm1_cvx(self): + print ("test_FBPD_Norm1_cvx") + if not cvx_not_installable: + opt = {'memopt': True} + # Problem data. + m = 30 + n = 20 + np.random.seed(1) + Amat = np.random.randn(m, n) + A = LinearOperatorMatrix(Amat) + bmat = np.random.randn(m) + bmat.shape = (bmat.shape[0], 1) + + # A = Identity() + # Change n to equal to m. + + b = DataContainer(bmat) + + # Regularization parameter + lam = 10 + opt = {'memopt': True} + # Initial guess + x_init = DataContainer(np.random.randn(n, 1)) + + # Create object instances with the test data A and b. + f = Norm2sq(A, b, c=0.5, memopt=True) + f.L = PowerMethodNonsquare(A, 25, x_init)[0] + print ("Lipschitz", f.L) + g0 = ZeroFun() + + + # Create 1-norm object instance + g1 = Norm1(lam) + + # Compare to CVXPY + + # Construct the problem. + x1 = Variable(n) + objective1 = Minimize( + 0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1, 1)) + prob1 = Problem(objective1) + + # The optimal objective is returned by prob.solve(). + result1 = prob1.solve(verbose=False, solver=SCS, eps=1e-9) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus 1-norm solution and objective value:") + print(x1.value) + print(objective1.value) + + # Now try another algorithm FBPD for same problem: + x_fbpd1, itfbpd1, timingfbpd1, criterfbpd1 = FBPD(x_init, + Identity(), None, f, g1) + print(x_fbpd1) + print(criterfbpd1[-1]) + + self.assertNumpyArrayAlmostEqual( + numpy.squeeze(x_fbpd1.array), x1.value, 6) + else: + self.assertTrue(cvx_not_installable) + # Plot criterion curve to see both FISTA and FBPD converge to same value. + # Note that FISTA is very efficient for 1-norm minimization so it beats + # FBPD in this test by a lot. But FBPD can handle a larger class of problems + # than FISTA can. + + # Now try 1-norm and TV denoising with FBPD, first 1-norm. + + # Set up phantom size NxN by creating ImageGeometry, initialising the + # ImageData object with this geometry and empty array and finally put some + # data into its array, and display as image. + def skip_test_FISTA_denoise_cvx(self): + if not cvx_not_installable: + opt = {'memopt': True} + N = 64 + ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) + Phantom = ImageData(geometry=ig) + + x = Phantom.as_array() + + x[int(round(N/4)):int(round(3*N/4)), + int(round(N/4)):int(round(3*N/4))] = 0.5 + x[int(round(N/8)):int(round(7*N/8)), + int(round(3*N/8)):int(round(5*N/8))] = 1 + + # Identity operator for denoising + I = TomoIdentity(ig) + + # Data and add noise + y = I.direct(Phantom) + y.array = y.array + 0.1*np.random.randn(N, N) + + # Data fidelity term + f_denoise = Norm2sq(I, y, c=0.5, memopt=True) + x_init = ImageData(geometry=ig) + f_denoise.L = PowerMethodNonsquare(I, 25, x_init)[0] + + # 1-norm regulariser + lam1_denoise = 1.0 + g1_denoise = Norm1(lam1_denoise) + + # Initial guess + x_init_denoise = ImageData(np.zeros((N, N))) + + # Combine with least squares and solve using generic FISTA implementation + x_fista1_denoise, it1_denoise, timing1_denoise, \ + criter1_denoise = \ + FISTA(x_init_denoise, f_denoise, g1_denoise, opt=opt) + + print(x_fista1_denoise) + print(criter1_denoise[-1]) + + # Now denoise LS + 1-norm with FBPD + x_fbpd1_denoise, itfbpd1_denoise, timingfbpd1_denoise,\ + criterfbpd1_denoise = \ + FBPD(x_init_denoise, I, None, f_denoise, g1_denoise) + print(x_fbpd1_denoise) + print(criterfbpd1_denoise[-1]) + + # Compare to CVXPY + + # Construct the problem. + x1_denoise = Variable(N**2) + objective1_denoise = Minimize( + 0.5*sum_squares(x1_denoise - y.array.flatten()) + lam1_denoise*norm(x1_denoise, 1)) + prob1_denoise = Problem(objective1_denoise) + + # The optimal objective is returned by prob.solve(). + result1_denoise = prob1_denoise.solve( + verbose=False, solver=SCS, eps=1e-12) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus 1-norm solution and objective value:") + print(x1_denoise.value) + print(objective1_denoise.value) + self.assertNumpyArrayAlmostEqual( + x_fista1_denoise.array.flatten(), x1_denoise.value, 5) + + self.assertNumpyArrayAlmostEqual( + x_fbpd1_denoise.array.flatten(), x1_denoise.value, 5) + x1_cvx = x1_denoise.value + x1_cvx.shape = (N, N) + + # Now TV with FBPD + lam_tv = 0.1 + gtv = TV2D(lam_tv) + gtv(gtv.op.direct(x_init_denoise)) + + opt_tv = {'tol': 1e-4, 'iter': 10000} + + x_fbpdtv_denoise, itfbpdtv_denoise, timingfbpdtv_denoise,\ + criterfbpdtv_denoise = \ + FBPD(x_init_denoise, gtv.op, None, f_denoise, gtv, opt=opt_tv) + print(x_fbpdtv_denoise) + print(criterfbpdtv_denoise[-1]) + + # Compare to CVXPY + + # Construct the problem. + xtv_denoise = Variable((N, N)) + objectivetv_denoise = Minimize( + 0.5*sum_squares(xtv_denoise - y.array) + lam_tv*tv(xtv_denoise)) + probtv_denoise = Problem(objectivetv_denoise) + + # The optimal objective is returned by prob.solve(). + resulttv_denoise = probtv_denoise.solve( + verbose=False, solver=SCS, eps=1e-12) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus 1-norm solution and objective value:") + print(xtv_denoise.value) + print(objectivetv_denoise.value) + + self.assertNumpyArrayAlmostEqual( + x_fbpdtv_denoise.as_array(), xtv_denoise.value, 1) + + else: + self.assertTrue(cvx_not_installable) + + +class TestFunction(unittest.TestCase): + def assertNumpyArrayEqual(self, first, second): + res = True + try: + numpy.testing.assert_array_equal(first, second) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def create_simple_ImageData(self): + N = 64 + ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) + Phantom = ImageData(geometry=ig) + + x = Phantom.as_array() + + x[int(round(N/4)):int(round(3*N/4)), + int(round(N/4)):int(round(3*N/4))] = 0.5 + x[int(round(N/8)):int(round(7*N/8)), + int(round(3*N/8)):int(round(5*N/8))] = 1 + + return (ig, Phantom) + + def _test_Norm2(self): + print("test Norm2") + opt = {'memopt': True} + ig, Phantom = self.create_simple_ImageData() + x = Phantom.as_array() + print(Phantom) + print(Phantom.as_array()) + + norm2 = Norm2() + v1 = norm2(x) + v2 = norm2(Phantom) + self.assertEqual(v1, v2) + + p1 = norm2.prox(Phantom, 1) + print(p1) + p2 = norm2.prox(x, 1) + self.assertNumpyArrayEqual(p1.as_array(), p2) + + p3 = norm2.proximal(Phantom, 1) + p4 = norm2.proximal(x, 1) + self.assertNumpyArrayEqual(p3.as_array(), p2) + self.assertNumpyArrayEqual(p3.as_array(), p4) + + out = Phantom.copy() + p5 = norm2.proximal(Phantom, 1, out=out) + self.assertEqual(id(p5), id(out)) + self.assertNumpyArrayEqual(p5.as_array(), p3.as_array()) +# ============================================================================= +# def test_upper(self): +# self.assertEqual('foo'.upper(), 'FOO') +# +# def test_isupper(self): +# self.assertTrue('FOO'.isupper()) +# self.assertFalse('Foo'.isupper()) +# +# def test_split(self): +# s = 'hello world' +# self.assertEqual(s.split(), ['hello', 'world']) +# # check that s.split fails when the separator is not a string +# with self.assertRaises(TypeError): +# s.split(2) +# ============================================================================= + + + +if __name__ == '__main__': + unittest.main() + diff --git a/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py b/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py index 0ec116f..8370c78 100644 --- a/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py +++ b/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py @@ -32,8 +32,8 @@ ProjAngleChannels = np.zeros((totalAngles,totChannels,n,n),dtype='float32') ######################################################################### print ("Loading the data...") -MainPath = '/media/algol/336F96987817D4B4/DATA/IMAT_DATA/' # path to data -pathname0 = '{!s}{!s}'.format(MainPath,'PSI_DATA/DATA/Sample/') +MainPath = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/' # path to data +pathname0 = '{!s}{!s}'.format(MainPath,'PSI_phantom_IMAT/DATA/Sample/') counterFileName = 4675 # A main loop over all available angles for ll in range(0,totalAngles,1): @@ -66,7 +66,7 @@ for ll in range(0,totalAngles,1): counterFileName += 1 ######################################################################### -flat1 = read_fits('{!s}{!s}{!s}'.format(MainPath,'PSI_DATA/DATA/','OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits')) +flat1 = read_fits('{!s}{!s}{!s}'.format(MainPath,'PSI_phantom_IMAT/DATA/','OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits')) nonzero = flat1 > 0 # Apply flat field and take negative log for ll in range(0,totalAngles,1): diff --git a/Wrappers/Python/wip/demo_imat_whitebeam.py b/Wrappers/Python/wip/demo_imat_whitebeam.py index 482c1ae..e0d213e 100644 --- a/Wrappers/Python/wip/demo_imat_whitebeam.py +++ b/Wrappers/Python/wip/demo_imat_whitebeam.py @@ -21,18 +21,18 @@ from ccpi.optimisation.algs import CGLS, FISTA from ccpi.optimisation.funcs import Norm2sq, Norm1 # Load and display a couple of summed projection as examples -pathname0 = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle0/' +pathname0 = '/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle0/' filename0 = 'IMAT00004675_Tomo_test_000_SummedImg.fits' data0 = read_fits(pathname0 + filename0) -pathname10 = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle10/' +pathname10 = '/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle10/' filename10 = 'IMAT00004685_Tomo_test_000_SummedImg.fits' data10 = read_fits(pathname10 + filename10) # Load a flat field (more are available, should we average over them?) -flat1 = read_fits('/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits') +flat1 = read_fits('/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits') # Apply flat field and display after flat-field correction and negative log data0_rel = numpy.zeros(numpy.shape(flat1), dtype = float) @@ -58,7 +58,7 @@ plt.colorbar() plt.show() # Set up for loading all summed images at 250 angles. -pathname = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle{}/' +pathname = '/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle{}/' filename = 'IMAT0000{}_Tomo_test_000_SummedImg.fits' # Dimensions |