diff options
42 files changed, 4269 insertions, 595 deletions
diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py new file mode 100755 index 0000000..358ba2d --- /dev/null +++ b/Wrappers/Python/ccpi/framework/BlockDataContainer.py @@ -0,0 +1,306 @@ +# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 5 16:04:45 2019
+
+@author: ofn77899
+"""
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import numpy
+from numbers import Number
+import functools
+#from ccpi.framework import AcquisitionData, ImageData
+#from ccpi.optimisation.operators import Operator, LinearOperator
+
+class BlockDataContainer(object):
+ '''Class to hold DataContainers as column vector'''
+ __array_priority__ = 1
+ def __init__(self, *args, **kwargs):
+ ''''''
+ self.containers = args
+ self.index = 0
+ #shape = kwargs.get('shape', None)
+ #if shape is None:
+ # shape = (len(args),1)
+ shape = (len(args),1)
+ self.shape = shape
+ #print (self.shape)
+ n_elements = functools.reduce(lambda x,y: x*y, shape, 1)
+ if len(args) != n_elements:
+ raise ValueError(
+ 'Dimension and size do not match: expected {} got {}'
+ .format(n_elements, len(args)))
+
+
+ def __iter__(self):
+ '''BlockDataContainer is Iterable'''
+ return self
+ def next(self):
+ '''python2 backwards compatibility'''
+ return self.__next__()
+ def __next__(self):
+ try:
+ out = self[self.index]
+ except IndexError as ie:
+ raise StopIteration()
+ self.index+=1
+ return out
+
+ def is_compatible(self, other):
+ '''basic check if the size of the 2 objects fit'''
+ if isinstance(other, Number):
+ return True
+ elif isinstance(other, list):
+ for ot in other:
+ if not isinstance(ot, (Number,\
+ numpy.int, numpy.int8, numpy.int16, numpy.int32, numpy.int64,\
+ numpy.float, numpy.float16, numpy.float32, numpy.float64, \
+ numpy.complex)):
+ raise ValueError('List/ numpy array can only contain numbers {}'\
+ .format(type(ot)))
+ return len(self.containers) == len(other)
+ elif isinstance(other, numpy.ndarray):
+ return self.shape == other.shape
+ return len(self.containers) == len(other.containers)
+
+ def get_item(self, row):
+ if row > self.shape[0]:
+ raise ValueError('Requested row {} > max {}'.format(row, self.shape[0]))
+ return self.containers[row]
+
+ def __getitem__(self, row):
+ return self.get_item(row)
+
+ def add(self, other, *args, **kwargs):
+ assert self.is_compatible(other)
+ out = kwargs.get('out', None)
+ #print ("args" , *args)
+ if isinstance(other, Number):
+ return type(self)(*[ el.add(other, *args, **kwargs) for el in self.containers], shape=self.shape)
+ elif isinstance(other, list) or isinstance(other, numpy.ndarray):
+ return type(self)(*[ el.add(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)
+ return type(self)(
+ *[ el.add(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other.containers)],
+ shape=self.shape)
+
+ def subtract(self, other, *args, **kwargs):
+ assert self.is_compatible(other)
+ out = kwargs.get('out', None)
+ if isinstance(other, Number):
+ return type(self)(*[ el.subtract(other, out, *args, **kwargs) for el in self.containers], shape=self.shape)
+ elif isinstance(other, list) or isinstance(other, numpy.ndarray):
+ return type(self)(*[ el.subtract(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)
+ return type(self)(*[ el.subtract(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other.containers)],
+ shape=self.shape)
+
+ def multiply(self, other, *args, **kwargs):
+ self.is_compatible(other)
+ out = kwargs.get('out', None)
+ if isinstance(other, Number):
+ return type(self)(*[ el.multiply(other, *args, **kwargs) for el in self.containers], shape=self.shape)
+ elif isinstance(other, list):
+ return type(self)(*[ el.multiply(ot, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)
+ elif isinstance(other, numpy.ndarray):
+ return type(self)(*[ el.multiply(ot, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)
+ return type(self)(*[ el.multiply(ot, *args, **kwargs) for el,ot in zip(self.containers,other.containers)],
+ shape=self.shape)
+
+ def divide(self, other, *args, **kwargs):
+ self.is_compatible(other)
+ out = kwargs.get('out', None)
+ if isinstance(other, Number):
+ return type(self)(*[ el.divide(other, *args, **kwargs) for el in self.containers], shape=self.shape)
+ elif isinstance(other, list) or isinstance(other, numpy.ndarray):
+ return type(self)(*[ el.divide(ot, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)
+ return type(self)(*[ el.divide(ot, *args, **kwargs) for el,ot in zip(self.containers,other.containers)],
+ shape=self.shape)
+
+ def power(self, other, *args, **kwargs):
+ assert self.is_compatible(other)
+ out = kwargs.get('out', None)
+ if isinstance(other, Number):
+ return type(self)(*[ el.power(other, *args, **kwargs) for el in self.containers], shape=self.shape)
+ elif isinstance(other, list) or isinstance(other, numpy.ndarray):
+ return type(self)(*[ el.power(ot, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)
+ return type(self)(*[ el.power(ot, *args, **kwargs) for el,ot in zip(self.containers,other.containers)], shape=self.shape)
+
+ def maximum(self,other, *args, **kwargs):
+ assert self.is_compatible(other)
+ out = kwargs.get('out', None)
+ if isinstance(other, Number):
+ return type(self)(*[ el.maximum(other, *args, **kwargs) for el in self.containers], shape=self.shape)
+ elif isinstance(other, list) or isinstance(other, numpy.ndarray):
+ return type(self)(*[ el.maximum(ot, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)
+ return type(self)(*[ el.maximum(ot, *args, **kwargs) for el,ot in zip(self.containers,other.containers)], shape=self.shape)
+
+ ## unary operations
+ def abs(self, *args, **kwargs):
+ out = kwargs.get('out', None)
+ return type(self)(*[ el.abs(*args, **kwargs) for el in self.containers], shape=self.shape)
+ def sign(self, *args, **kwargs):
+ out = kwargs.get('out', None)
+ return type(self)(*[ el.sign(*args, **kwargs) for el in self.containers], shape=self.shape)
+ def sqrt(self, *args, **kwargs):
+ out = kwargs.get('out', None)
+ return type(self)(*[ el.sqrt(*args, **kwargs) for el in self.containers], shape=self.shape)
+ def conjugate(self, out=None):
+ return type(self)(*[el.conjugate() for el in self.containers], shape=self.shape)
+
+ ## reductions
+ def sum(self, *args, **kwargs):
+ return numpy.sum([ el.sum(*args, **kwargs) for el in self.containers])
+ def squared_norm(self):
+ y = numpy.asarray([el.squared_norm() for el in self.containers])
+ return y.sum()
+ def norm(self):
+ return numpy.sqrt(self.squared_norm())
+ def copy(self):
+ '''alias of clone'''
+ return self.clone()
+ def clone(self):
+ return type(self)(*[el.copy() for el in self.containers], shape=self.shape)
+
+ def __add__(self, other):
+ return self.add( other )
+ # __radd__
+
+ def __sub__(self, other):
+ return self.subtract( other )
+ # __rsub__
+
+ def __mul__(self, other):
+ return self.multiply(other)
+ # __rmul__
+
+ def __div__(self, other):
+ return self.divide(other)
+ # __rdiv__
+ def __truediv__(self, other):
+ return self.divide(other)
+
+ def __pow__(self, other):
+ return self.power(other)
+ # reverse operand
+ def __radd__(self, other):
+ '''Reverse addition
+
+ to make sure that this method is called rather than the __mul__ of a numpy array
+ the class constant __array_priority__ must be set > 0
+ https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
+ '''
+ return self + other
+ # __radd__
+
+ def __rsub__(self, other):
+ '''Reverse subtraction
+
+ to make sure that this method is called rather than the __mul__ of a numpy array
+ the class constant __array_priority__ must be set > 0
+ https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
+ '''
+ return (-1 * self) + other
+ # __rsub__
+
+ def __rmul__(self, other):
+ '''Reverse multiplication
+
+ to make sure that this method is called rather than the __mul__ of a numpy array
+ the class constant __array_priority__ must be set > 0
+ https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
+ '''
+ return self * other
+ # __rmul__
+
+ def __rdiv__(self, other):
+ '''Reverse division
+
+ to make sure that this method is called rather than the __mul__ of a numpy array
+ the class constant __array_priority__ must be set > 0
+ https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
+ '''
+ return pow(self / other, -1)
+ # __rdiv__
+ def __rtruediv__(self, other):
+ '''Reverse truedivision
+
+ to make sure that this method is called rather than the __mul__ of a numpy array
+ the class constant __array_priority__ must be set > 0
+ https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
+ '''
+ return self.__rdiv__(other)
+
+ def __rpow__(self, other):
+ '''Reverse power
+
+ to make sure that this method is called rather than the __mul__ of a numpy array
+ the class constant __array_priority__ must be set > 0
+ https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
+ '''
+ return other.power(self)
+
+ def __iadd__(self, other):
+ '''Inline addition'''
+ if isinstance (other, BlockDataContainer):
+ for el,ot in zip(self.containers, other.containers):
+ el += ot
+ elif isinstance(other, Number):
+ for el in self.containers:
+ el += other
+ elif isinstance(other, list) or isinstance(other, numpy.ndarray):
+ self.is_compatible(other)
+ for el,ot in zip(self.containers, other):
+ el += ot
+ return self
+ # __iadd__
+
+ def __isub__(self, other):
+ '''Inline subtraction'''
+ if isinstance (other, BlockDataContainer):
+ for el,ot in zip(self.containers, other.containers):
+ el -= ot
+ elif isinstance(other, Number):
+ for el in self.containers:
+ el -= other
+ elif isinstance(other, list) or isinstance(other, numpy.ndarray):
+ assert self.is_compatible(other)
+ for el,ot in zip(self.containers, other):
+ el -= ot
+ return self
+ # __isub__
+
+ def __imul__(self, other):
+ '''Inline multiplication'''
+ if isinstance (other, BlockDataContainer):
+ for el,ot in zip(self.containers, other.containers):
+ el *= ot
+ elif isinstance(other, Number):
+ for el in self.containers:
+ el *= other
+ elif isinstance(other, list) or isinstance(other, numpy.ndarray):
+ assert self.is_compatible(other)
+ for el,ot in zip(self.containers, other):
+ el *= ot
+ return self
+ # __imul__
+
+ def __idiv__(self, other):
+ '''Inline division'''
+ if isinstance (other, BlockDataContainer):
+ for el,ot in zip(self.containers, other.containers):
+ el /= ot
+ elif isinstance(other, Number):
+ for el in self.containers:
+ el /= other
+ elif isinstance(other, list) or isinstance(other, numpy.ndarray):
+ assert self.is_compatible(other)
+ for el,ot in zip(self.containers, other):
+ el /= ot
+ return self
+ # __rdiv__
+ def __itruediv__(self, other):
+ '''Inline truedivision'''
+ return self.__idiv__(other)
+
diff --git a/Wrappers/Python/ccpi/framework/BlockGeometry.py b/Wrappers/Python/ccpi/framework/BlockGeometry.py new file mode 100755 index 0000000..87dfe92 --- /dev/null +++ b/Wrappers/Python/ccpi/framework/BlockGeometry.py @@ -0,0 +1,34 @@ +from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import numpy
+from numbers import Number
+import functools
+#from ccpi.framework import AcquisitionData, ImageData
+#from ccpi.optimisation.operators import Operator, LinearOperator
+
+class BlockGeometry(object):
+ '''Class to hold Geometry as column vector'''
+ #__array_priority__ = 1
+ def __init__(self, *args, **kwargs):
+ ''''''
+ self.geometries = args
+ self.index = 0
+ #shape = kwargs.get('shape', None)
+ #if shape is None:
+ # shape = (len(args),1)
+ shape = (len(args),1)
+ self.shape = shape
+ #print (self.shape)
+ n_elements = functools.reduce(lambda x,y: x*y, shape, 1)
+ if len(args) != n_elements:
+ raise ValueError(
+ 'Dimension and size do not match: expected {} got {}'
+ .format(n_elements, len(args)))
+
+ def allocate(self):
+ containers = [geom.allocate() for geom in self.geometries]
+ BlockDataContainer(*containers)
+
diff --git a/Wrappers/Python/ccpi/framework/__init__.py b/Wrappers/Python/ccpi/framework/__init__.py new file mode 100755 index 0000000..66e2f56 --- /dev/null +++ b/Wrappers/Python/ccpi/framework/__init__.py @@ -0,0 +1,25 @@ +# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 5 16:00:18 2019
+
+@author: ofn77899
+"""
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import numpy
+import sys
+from datetime import timedelta, datetime
+import warnings
+from functools import reduce
+
+from .framework import DataContainer
+from .framework import ImageData, AcquisitionData
+from .framework import ImageGeometry, AcquisitionGeometry
+from .framework import find_key, message
+from .framework import DataProcessor
+from .framework import AX, PixelByPixelDataProcessor, CastDataContainer
+from .BlockDataContainer import BlockDataContainer
+from .BlockGeometry import BlockGeometry
diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework/framework.py index 24f4ca6..0c43737 100644..100755 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework/framework.py @@ -27,6 +27,7 @@ import sys from datetime import timedelta, datetime import warnings from functools import reduce +from numbers import Number def find_key(dic, val): """return the key of dictionary dic given the value""" @@ -67,6 +68,28 @@ class ImageGeometry(object): self.center_z = center_z self.channels = channels + # this is some code repetition + if self.channels > 1: + if self.voxel_num_z>1: + self.length = 4 + self.shape = (self.channels, self.voxel_num_z, self.voxel_num_y, self.voxel_num_x) + dim_labels = ['channel' ,'vertical' , 'horizontal_y' , 'horizontal_x'] + else: + self.length = 3 + self.shape = (self.channels, self.voxel_num_y, self.voxel_num_x) + dim_labels = ['channel' , 'horizontal_y' , 'horizontal_x'] + else: + if self.voxel_num_z>1: + self.length = 3 + self.shape = (self.voxel_num_z, self.voxel_num_y, self.voxel_num_x) + dim_labels = ['vertical', 'horizontal_y' , 'horizontal_x'] + else: + self.length = 2 + self.shape = (self.voxel_num_y, self.voxel_num_x) + dim_labels = ['horizontal_y' , 'horizontal_x'] + + self.dimension_labels = dim_labels + def get_min_x(self): return self.center_x - 0.5*self.voxel_num_x*self.voxel_size_x @@ -113,9 +136,18 @@ class ImageGeometry(object): return repres def allocate(self, value=0, dimension_labels=None): '''allocates an ImageData according to the size expressed in the instance''' - out = ImageData(geometry=self, dimension_labels=dimension_labels) - if value != 0: - out += value + out = ImageData(geometry=self) + if isinstance(value, Number): + if value != 0: + out += value + else: + if value == 'random': + out.fill(numpy.random.random_sample(self.shape)) + elif value == 'random_int': + out.fill(numpy.random.randint(1, 10 + 1,size=self.shape)) + if dimension_labels is not None: + if dimension_labels != self.dimension_labels: + return out.subset(dimensions=dimension_labels) return out class AcquisitionGeometry(object): @@ -130,7 +162,7 @@ class AcquisitionGeometry(object): dist_source_center=None, dist_center_detector=None, channels=1, - angle_unit='degree' + angle_unit='degree', ): """ General inputs for standard type projection geometries @@ -161,6 +193,7 @@ class AcquisitionGeometry(object): self.geom_type = geom_type # 'parallel' or 'cone' self.dimension = dimension # 2D or 3D self.angles = angles + num_of_angles = len (angles) self.dist_source_center = dist_source_center self.dist_center_detector = dist_center_detector @@ -171,6 +204,24 @@ class AcquisitionGeometry(object): self.pixel_size_v = pixel_size_v self.channels = channels + + if channels > 1: + if pixel_num_v > 1: + shape = (channels, num_of_angles , pixel_num_v, pixel_num_h) + dim_labels = ['channel' , 'angle' , 'vertical' , 'horizontal'] + else: + shape = (channels , num_of_angles, pixel_num_h) + dim_labels = ['channel' , 'angle' , 'horizontal'] + else: + if pixel_num_v > 1: + shape = (num_of_angles, pixel_num_v, pixel_num_h) + dim_labels = ['angle' , 'vertical' , 'horizontal'] + else: + shape = (num_of_angles, pixel_num_h) + dim_labels = ['angle' , 'horizontal'] + self.shape = shape + + self.dimension_labels = dim_labels def clone(self): '''returns a copy of the AcquisitionGeometry''' @@ -198,9 +249,18 @@ class AcquisitionGeometry(object): return repres def allocate(self, value=0, dimension_labels=None): '''allocates an AcquisitionData according to the size expressed in the instance''' - out = AcquisitionData(geometry=self, dimension_labels=dimension_labels) - if value != 0: - out += value + out = AcquisitionData(geometry=self) + if isinstance(value, Number): + if value != 0: + out += value + else: + if value == 'random': + out.fill(numpy.random.random_sample(self.shape)) + elif value == 'random_int': + out.fill(numpy.random.out.fill(numpy.random.randint(1, 10 + 1,size=self.shape))) + if dimension_labels is not None: + if dimension_labels != self.dimension_labels: + return out.subset(dimensions=dimension_labels) return out class DataContainer(object): '''Generic class to hold data @@ -382,7 +442,8 @@ class DataContainer(object): return self.shape == other.shape ## algebra - def __add__(self, other , out=None, *args, **kwargs): + def __add__(self, other, *args, **kwargs): + out = kwargs.get('out', None) if issubclass(type(other), DataContainer): if self.check_dimensions(other): out = self.as_array() + other.as_array() @@ -639,8 +700,8 @@ class DataContainer(object): ## binary operations - def pixel_wise_binary(self,pwop, x2 , out=None, *args, **kwargs): - + def pixel_wise_binary(self, pwop, x2, *args, **kwargs): + out = kwargs.get('out', None) if out is None: if isinstance(x2, (int, float, complex)): out = pwop(self.as_array() , x2 , *args, **kwargs ) @@ -658,7 +719,8 @@ class DataContainer(object): elif issubclass(type(out), DataContainer) and issubclass(type(x2), DataContainer): if self.check_dimensions(out) and self.check_dimensions(x2): - pwop(self.as_array(), x2.as_array(), out=out.as_array(), *args, **kwargs ) + kwargs['out'] = out.as_array() + pwop(self.as_array(), x2.as_array(), *args, **kwargs ) #return type(self)(out.as_array(), # deep_copy=False, # dimension_labels=self.dimension_labels, @@ -668,14 +730,15 @@ class DataContainer(object): raise ValueError(message(type(self),"Wrong size for data memory: ", out.shape,self.shape)) elif issubclass(type(out), DataContainer) and isinstance(x2, (int,float,complex)): if self.check_dimensions(out): - - pwop(self.as_array(), x2, out=out.as_array(), *args, **kwargs ) + kwargs['out']=out.as_array() + pwop(self.as_array(), x2, *args, **kwargs ) return out else: raise ValueError(message(type(self),"Wrong size for data memory: ", out.shape,self.shape)) elif issubclass(type(out), numpy.ndarray): if self.array.shape == out.shape and self.array.dtype == out.dtype: - pwop(self.as_array(), x2 , out=out, *args, **kwargs) + kwargs['out'] = out + pwop(self.as_array(), x2, *args, **kwargs) #return type(self)(out, # deep_copy=False, # dimension_labels=self.dimension_labels, @@ -683,26 +746,27 @@ class DataContainer(object): else: raise ValueError (message(type(self), "incompatible class:" , pwop.__name__, type(out))) - def add(self, other , out=None, *args, **kwargs): - return self.pixel_wise_binary(numpy.add, other, out=out, *args, **kwargs) + def add(self, other, *args, **kwargs): + return self.pixel_wise_binary(numpy.add, other, *args, **kwargs) - def subtract(self, other, out=None , *args, **kwargs): - return self.pixel_wise_binary(numpy.subtract, other, out=out, *args, **kwargs) + def subtract(self, other, *args, **kwargs): + return self.pixel_wise_binary(numpy.subtract, other, *args, **kwargs) - def multiply(self, other , out=None, *args, **kwargs): - return self.pixel_wise_binary(numpy.multiply, other, out=out, *args, **kwargs) + def multiply(self, other, *args, **kwargs): + return self.pixel_wise_binary(numpy.multiply, other, *args, **kwargs) - def divide(self, other , out=None ,*args, **kwargs): - return self.pixel_wise_binary(numpy.divide, other, out=out, *args, **kwargs) + def divide(self, other, *args, **kwargs): + return self.pixel_wise_binary(numpy.divide, other, *args, **kwargs) - def power(self, other , out=None, *args, **kwargs): - return self.pixel_wise_binary(numpy.power, other, out=out, *args, **kwargs) + def power(self, other, *args, **kwargs): + return self.pixel_wise_binary(numpy.power, other, *args, **kwargs) - def maximum(self,x2, out=None, *args, **kwargs): - return self.pixel_wise_binary(numpy.maximum, x2=x2, out=out, *args, **kwargs) + def maximum(self, x2, *args, **kwargs): + return self.pixel_wise_binary(numpy.maximum, x2, *args, **kwargs) ## unary operations - def pixel_wise_unary(self,pwop, out=None, *args, **kwargs): + def pixel_wise_unary(self, pwop, *args, **kwargs): + out = kwargs.get('out', None) if out is None: out = pwop(self.as_array() , *args, **kwargs ) return type(self)(out, @@ -711,38 +775,43 @@ class DataContainer(object): geometry=self.geometry) elif issubclass(type(out), DataContainer): if self.check_dimensions(out): - pwop(self.as_array(), out=out.as_array(), *args, **kwargs ) + kwargs['out'] = out.as_array() + pwop(self.as_array(), *args, **kwargs ) else: raise ValueError(message(type(self),"Wrong size for data memory: ", out.shape,self.shape)) elif issubclass(type(out), numpy.ndarray): if self.array.shape == out.shape and self.array.dtype == out.dtype: - pwop(self.as_array(), out=out, *args, **kwargs) + kwargs['out'] = out + pwop(self.as_array(), *args, **kwargs) else: raise ValueError (message(type(self), "incompatible class:" , pwop.__name__, type(out))) - def abs(self, out=None, *args, **kwargs): - return self.pixel_wise_unary(numpy.abs, out=out, *args, **kwargs) - - def sign(self, out=None, *args, **kwargs): - return self.pixel_wise_unary(numpy.sign , out=out, *args, **kwargs) + def abs(self, *args, **kwargs): + return self.pixel_wise_unary(numpy.abs, *args, **kwargs) - def sqrt(self, out=None, *args, **kwargs): - return self.pixel_wise_unary(numpy.sqrt, out=out, *args, **kwargs) + def sign(self, *args, **kwargs): + return self.pixel_wise_unary(numpy.sign, *args, **kwargs) + def sqrt(self, *args, **kwargs): + return self.pixel_wise_unary(numpy.sqrt, *args, **kwargs) + + def conjugate(self, *args, **kwargs): + return self.pixel_wise_unary(numpy.conjugate, *args, **kwargs) #def __abs__(self): # operation = FM.OPERATION.ABS # return self.callFieldMath(operation, None, self.mask, self.maskOnValue) # __abs__ ## reductions - def sum(self, out=None, *args, **kwargs): + def sum(self, *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()) + #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()) + return self.dot(self) def norm(self): '''return the euclidean norm of the DataContainer viewed as a vector''' return numpy.sqrt(self.squared_norm()) @@ -756,6 +825,7 @@ class DataContainer(object): + class ImageData(DataContainer): '''DataContainer for holding 2D or 3D DataContainer''' def __init__(self, @@ -889,7 +959,7 @@ class AcquisitionData(DataContainer): if channels > 1: if vert > 1: shape = (channels, num_of_angles , vert, horiz) - dim_labels = ['channel' , ' angle' , + dim_labels = ['channel' , 'angle' , 'vertical' , 'horizontal'] else: shape = (channels , num_of_angles, horiz) @@ -918,7 +988,7 @@ class AcquisitionData(DataContainer): elif dim == 'horizontal': shape.append(horiz) if len(shape) != len(dimension_labels): - raise ValueError('Missing {0} axes.\nExpected{1} got {2}}'\ + raise ValueError('Missing {0} axes.\nExpected{1} got {2}'\ .format( len(dimension_labels) - len(shape), dimension_labels, shape) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py index 680b268..cc99344 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -140,7 +140,7 @@ class Algorithm(object): 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): + def run(self, iterations, verbose=True, callback=None): '''run n iterations and update the user with the callback if specified''' if self.should_stop(): print ("Stop cryterion has been reached.") @@ -149,8 +149,9 @@ class Algorithm(object): 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()) + else: + 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/GradientDescent.py b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py index 7794b4d..f1e4132 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py @@ -51,13 +51,17 @@ class GradientDescent(Algorithm): 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 - + try: + self.memopt = self.objective_function.memopt + except AttributeError as ae: + self.memopt = False + if self.memopt: + self.x_update = x_init.copy() + def update(self): '''Single iteration''' if self.memopt: @@ -65,7 +69,7 @@ class GradientDescent(Algorithm): self.x_update *= -self.rate self.x += self.x_update else: - self.x += -self.rate * self.objective_function.grad(self.x) + self.x += -self.rate * self.objective_function.gradient(self.x) def update_objective(self): self.loss.append(self.objective_function(self.x)) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py new file mode 100644 index 0000000..7488310 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Feb 4 16:18:06 2019 + +@author: evangelos +""" + +from ccpi.framework import ImageData +import numpy as np +import matplotlib.pyplot as plt +import time +from Operators.CompositeOperator import CompositeOperator +from Operators.CompositeDataContainer import CompositeDataContainer + +def PDHG(f, g, operator, tau = None, sigma = None, opt = None, **kwargs): + + # algorithmic parameters + if opt is None: + opt = {'tol': 1e-6, 'niter': 500, 'show_iter': 100, \ + 'memopt': False} + + if sigma is None and tau is None: + raise ValueError('Need sigma*tau||K||^2<1') + + niter = opt['niter'] if 'niter' in opt.keys() else 1000 + tol = opt['tol'] if 'tol' in opt.keys() else 1e-4 + memopt = opt['memopt'] if 'memopt' in opt.keys() else False + show_iter = opt['show_iter'] if 'show_iter' in opt.keys() else False + stop_crit = opt['stop_crit'] if 'stop_crit' in opt.keys() else False + + if isinstance(operator, CompositeOperator): +# if isinstance(operator, CompositeOperator_DataContainer): + x_old = operator.alloc_domain_dim() + y_old = operator.alloc_range_dim() + else: + x_old = ImageData(np.zeros(operator.domain_dim())) + y_old = ImageData(np.zeros(operator.range_dim())) + + + xbar = x_old + x_tmp = x_old + x = x_old + + y_tmp = y_old + y = y_tmp + + # relaxation parameter + theta = 1 + + t = time.time() + + objective = [] + + for i in range(niter): + + # Gradient descent, Dual problem solution + y_tmp = y_old + sigma * operator.direct(xbar) + y = f.proximal_conjugate(y_tmp, sigma) + + # Gradient ascent, Primal problem solution + x_tmp = x_old - tau * operator.adjoint(y) + x = g.proximal(x_tmp, tau) + + #Update + xbar = x + theta * (x - x_old) + + x_old = x + y_old = y + +# pdgap + print(f(x) + g(x) + f.convex_conjugate(y) + g.convex_conjugate(-1*operator.adjoint(y)) ) + + + + + +# # TV denoising, pdgap with composite +# +# primal_obj = f.get_item(0).alpha * ImageData(operator.compMat[0][0].direct(x.get_item(0)).power(2).sum(axis=0)).sqrt().sum() +\ +# 0.5*( (operator.compMat[1][0].direct(x.get_item(0)) - f.get_item(1).b).power(2).sum()) +# dual_obj = 0.5 * ((y.get_item(1).power(2)).sum()) + ( y.get_item(1)*f.get_item(1).b ).sum() + + # TV denoising, pdgap with no composite + + + +# primal_obj = f.get_item(0).alpha * ImageData(operator.compMat[0][0].direct(x.get_item(0)).power(2).sum(axis=0)).sqrt().sum() +\ +# 0.5*( (operator.compMat[1][0].direct(x.get_item(0)) - f.get_item(1).b).power(2).sum()) +# dual_obj = 0.5 * ((y.get_item(1).power(2)).sum()) + ( y.get_item(1)*f.get_item(1).b ).sum() + + +# print(primal_obj) +# objective = primal_obj +# + + + + t_end = time.time() + + return x, t_end - t, objective + diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py index 903bc30..7e500e8 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py @@ -27,3 +27,4 @@ from .CGLS import CGLS from .GradientDescent import GradientDescent from .FISTA import FISTA from .FBPD import FBPD + diff --git a/Wrappers/Python/ccpi/optimisation/funcs.py b/Wrappers/Python/ccpi/optimisation/funcs.py index 47ee810..8ce54c7 100755 --- a/Wrappers/Python/ccpi/optimisation/funcs.py +++ b/Wrappers/Python/ccpi/optimisation/funcs.py @@ -20,6 +20,7 @@ from ccpi.optimisation.ops import Identity, FiniteDiff2D import numpy from ccpi.framework import DataContainer +import warnings def isSizeCorrect(data1 ,data2): @@ -40,8 +41,12 @@ class Function(object): def __init__(self): self.L = None def __call__(self,x, out=None): raise NotImplementedError - def grad(self, x): raise NotImplementedError - def prox(self, x, tau): raise NotImplementedError + def grad(self, x): + warnings.warn("grad method is deprecated. use gradient instead", DeprecationWarning) + return self.gradient(x, out=None) + def prox(self, x, tau): + warnings.warn("prox method is deprecated. use proximal instead", DeprecationWarning) + return self.proximal(x,tau,out=None) def gradient(self, x, out=None): raise NotImplementedError def proximal(self, x, tau, out=None): raise NotImplementedError @@ -141,12 +146,20 @@ class Norm2sq(Function): self.A = A # Should be an operator, default identity self.b = b # Default zero DataSet? self.c = c # Default 1. - self.memopt = memopt if memopt: - #self.direct_placehold = A.adjoint(b) - self.direct_placehold = A.allocate_direct() - self.adjoint_placehold = A.allocate_adjoint() - + try: + self.adjoint_placehold = A.range_geometry().allocate() + self.direct_placehold = A.domain_geometry().allocate() + self.memopt = True + except NameError as ne: + warnings.warn(str(ne)) + self.memopt = False + except NotImplementedError as nie: + print (nie) + warnings.warn(str(nie)) + self.memopt = False + else: + self.memopt = False # Compute the Lipschitz parameter from the operator if possible # Leave it initialised to None otherwise @@ -154,11 +167,12 @@ class Norm2sq(Function): self.L = 2.0*self.c*(self.A.get_max_sing_val()**2) except AttributeError as ae: pass + except NotImplementedError as noe: + pass - def grad(self,x): - #return 2*self.c*self.A.adjoint( self.A.direct(x) - self.b ) - return (2.0*self.c)*self.A.adjoint( self.A.direct(x) - self.b ) - + #def grad(self,x): + # return self.gradient(x, out=None) + def __call__(self,x): #return self.c* np.sum(np.square((self.A.direct(x) - self.b).ravel())) #if out is None: @@ -181,12 +195,13 @@ class Norm2sq(Function): self.A.direct(x, out=self.adjoint_placehold) self.adjoint_placehold.__isub__( self.b ) self.A.adjoint(self.adjoint_placehold, out=self.direct_placehold) - self.direct_placehold.__imul__(2.0 * self.c) - # can this be avoided? - out.fill(self.direct_placehold) + #self.direct_placehold.__imul__(2.0 * self.c) + ## can this be avoided? + #out.fill(self.direct_placehold) + self.direct_placehold.multiply(2.0*self.c, out=out) else: - return self.grad(x) - + return (2.0*self.c)*self.A.adjoint( self.A.direct(x) - self.b ) + class ZeroFun(Function): diff --git a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py new file mode 100644 index 0000000..d6c98c4 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py @@ -0,0 +1,70 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 8 10:01:31 2019 + +@author: evangelos +""" + +import numpy as np +#from ccpi.optimisation.funcs import Function +from ccpi.optimisation.functions import Function +from ccpi.framework import BlockDataContainer + +class BlockFunction(Function): + + def __init__(self, operator, *functions): + + self.functions = functions + self.operator = operator + self.length = len(self.functions) + + super(BlockFunction, self).__init__() + + def __call__(self, x): + + tmp = self.operator.direct(x) + + t = 0 + for i in range(tmp.shape[0]): + t += self.functions[i](tmp.get_item(i)) + return t + + def call_adjoint(self, x): + + tmp = operator.adjoint(x) + + t = 0 + for i in range(tmp.shape[0]): + t += self.functions[i](tmp.get_item(i)) + return t + + def convex_conjugate(self, x): + + ''' Convex_conjugate does not take into account the BlockOperator''' + t = 0 + for i in range(x.shape[0]): + t += self.functions[i].convex_conjugate(x.get_item(i)) + return t + + + def proximal_conjugate(self, x, tau, out = None): + + ''' proximal_conjugate does not take into account the BlockOperator''' + out = [None]*self.length + for i in range(self.length): + out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau) + + return BlockDataContainer(*out) + + def proximal(self, x, tau, out = None): + + ''' proximal does not take into account the BlockOperator''' + out = [None]*self.length + for i in range(self.length): + out[i] = self.functions[i].proximal(x.get_item(i), tau) + + return BlockDataContainer(*out) + + def gradient(self,x, out=None): + pass
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/functions/Function.py b/Wrappers/Python/ccpi/optimisation/functions/Function.py new file mode 100644 index 0000000..82f24a6 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/Function.py @@ -0,0 +1,69 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2018-2019 Jakob Jorgensen, Daniil Kazantsev and Edoardo Pasca + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import warnings +from ccpi.optimisation.functions.ScaledFunction import ScaledFunction + +class Function(object): + '''Abstract class representing a function + + Members: + L is the Lipschitz constant of the gradient of the Function + ''' + def __init__(self): + self.L = None + + def __call__(self,x, out=None): + '''Evaluates the function at x ''' + raise NotImplementedError + + def gradient(self, x, out=None): + '''Returns the gradient of the function at x, if the function is differentiable''' + raise NotImplementedError + + def proximal(self, x, tau, out=None): + '''This returns the proximal operator for the function at x, tau''' + raise NotImplementedError + + def convex_conjugate(self, x, out=None): + '''This evaluates the convex conjugate of the function at x''' + raise NotImplementedError + + def proximal_conjugate(self, x, tau, out = None): + '''This returns the proximal operator for the convex conjugate of the function at x, tau''' + raise NotImplementedError + + def grad(self, x): + '''Alias of gradient(x,None)''' + warnings.warn('''This method will disappear in following + versions of the CIL. Use gradient instead''', DeprecationWarning) + return self.gradient(x, out=None) + + def prox(self, x, tau): + '''Alias of proximal(x, tau, None)''' + warnings.warn('''This method will disappear in following + versions of the CIL. Use proximal instead''', DeprecationWarning) + return self.proximal(x, out=None) + + def __rmul__(self, scalar): + '''Defines the multiplication by a scalar on the left + + returns a ScaledFunction''' + return ScaledFunction(self, scalar) + diff --git a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py new file mode 100644 index 0000000..3ac4358 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py @@ -0,0 +1,65 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 8 09:55:36 2019 + +@author: evangelos +""" + +import numpy as np +#from ccpi.optimisation.funcs import Function +from ccpi.optimisation.functions import Function +from ccpi.optimisation.functions import ScaledFunction + + +class FunctionOperatorComposition(Function): + + def __init__(self, operator, function): + super(FunctionOperatorComposition, self).__init__() + self.function = function + self.operator = operator + alpha = 1 + if isinstance (function, ScaledFunction): + alpha = function.scalar + self.L = 2 * alpha * operator.norm()**2 + + + def __call__(self, x): + + return self.function(self.operator.direct(x)) + + def call_adjoint(self, x): + + return self.function(self.operator.adjoint(x)) + + def convex_conjugate(self, x): + + ''' convex_conjugate does not take into account the Operator''' + return self.function.convex_conjugate(x) + + def proximal(self, x, tau): + + ''' proximal does not take into account the Operator''' + + return self.function.proximal(x, tau, out=None) + + def proximal_conjugate(self, x, tau, out=None): + + ''' proximal conjugate does not take into account the Operator''' + + return self.function.proximal_conjugate(x, tau) + + def gradient(self, x, out=None): + + ''' Gradient takes into account the Operator''' + if out is None: + return self.operator.adjoint( + self.function.gradient(self.operator.direct(x)) + ) + else: + self.operator.adjoint( + self.function.gradient(self.operator.direct(x), + out=out) + ) + +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py new file mode 100644 index 0000000..f83de6f --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py @@ -0,0 +1,76 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Mar 6 19:42:34 2019 + +@author: evangelos +""" + +import numpy as np +#from ccpi.optimisation.funcs import Function +from ccpi.optimisation.functions import Function +from ccpi.framework import DataContainer, ImageData, ImageGeometry + + +############################ L1NORM FUNCTIONS ############################# +class SimpleL1Norm(Function): + + def __init__(self, alpha=1): + + super(SimpleL1Norm, self).__init__() + self.alpha = alpha + + def __call__(self, x): + return self.alpha * x.abs().sum() + + def gradient(self,x): + return ValueError('Not Differentiable') + + def convex_conjugate(self,x): + return 0 + + def proximal(self, x, tau): + ''' Soft Threshold''' + return x.sign() * (x.abs() - tau * self.alpha).maximum(0) + + def proximal_conjugate(self, x, tau): + return x.divide((x.abs()/self.alpha).maximum(1.0)) + +class L1Norm(SimpleL1Norm): + + def __init__(self, alpha=1, **kwargs): + + super(L1Norm, self).__init__() + self.alpha = alpha + self.b = kwargs.get('b',None) + + def __call__(self, x): + + if self.b is None: + return SimpleL1Norm.__call__(self, x) + else: + return SimpleL1Norm.__call__(self, x - self.b) + + def gradient(self, x): + return ValueError('Not Differentiable') + + def convex_conjugate(self,x): + if self.b is None: + return SimpleL1Norm.convex_conjugate(self, x) + else: + return SimpleL1Norm.convex_conjugate(self, x) + (self.b * x).sum() + + def proximal(self, x, tau): + + if self.b is None: + return SimpleL1Norm.proximal(self, x, tau) + else: + return self.b + SimpleL1Norm.proximal(self, x - self.b , tau) + + def proximal_conjugate(self, x, tau): + + if self.b is None: + return SimpleL1Norm.proximal_conjugate(self, x, tau) + else: + return SimpleL1Norm.proximal_conjugate(self, x - tau*self.b, tau) +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py new file mode 100644 index 0000000..1baf365 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py @@ -0,0 +1,99 @@ +# -*- coding: utf-8 -*- + +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Thu Feb 7 13:10:56 2019 + +@author: evangelos +""" + +import numpy as np +#from ccpi.optimisation.funcs import Function +from ccpi.optimisation.functions import Function +from ccpi.framework import DataContainer, ImageData, ImageGeometry + + +class SimpleL2NormSq(Function): + + def __init__(self, alpha=1): + + super(SimpleL2NormSq, self).__init__() + # Lispchitz constant of gradient + self.L = 2 + + def __call__(self, x): + return x.power(2).sum() + + def gradient(self,x, out=None): + if out is None: + return 2 * x + else: + out.fill(2*x) + + def convex_conjugate(self,x): + return (1/4) * x.squared_norm() + + def proximal(self, x, tau, out=None): + if out is None: + return x.divide(1+2*tau) + else: + x.divide(1+2*tau, out=out) + + def proximal_conjugate(self, x, tau, out=None): + if out is None: + return x.divide(1 + tau/2) + else: + x.divide(1+tau/2, out=out) + + + +############################ L2NORM FUNCTIONS ############################# +class L2NormSq(SimpleL2NormSq): + + def __init__(self, **kwargs): + super(L2NormSq, self).__init__() + self.b = kwargs.get('b',None) + + def __call__(self, x): + if self.b is None: + return SimpleL2NormSq.__call__(self, x) + else: + return SimpleL2NormSq.__call__(self, x - self.b) + + def gradient(self, x, out=None): + if self.b is None: + return 2 * x + else: + return 2 * (x - self.b) + + def convex_conjugate(self, x): + ''' The convex conjugate corresponds to the simple functional i.e., + f(x) = alpha * ||x - b||_{2}^{2} + ''' + if self.b is None: + return SimpleL2NormSq.convex_conjugate(self, x) + else: + return SimpleL2NormSq.convex_conjugate(self, x) + (self.b * x).sum() + + def proximal(self, x, tau): + + ''' The proximal operator corresponds to the simple functional i.e., + f(x) = alpha * ||x - b||_{2}^{2} + + argmin_x { 0.5||x - u||^{2} + tau f(x) } + ''' + if self.b is None: + return SimpleL2NormSq.proximal(self, x, tau) + else: + return self.b + SimpleL2NormSq.proximal(self, x - self.b , tau) + + def proximal_conjugate(self, x, tau): + ''' The proximal operator corresponds to the simple convex conjugate + functional i.e., f^{*}(x^{) + argmin_x { 0.5||x - u||^{2} + tau f(x) } + ''' + if self.b is None: + return SimpleL2NormSq.proximal_conjugate(self, x, tau) + else: + return SimpleL2NormSq.proximal_conjugate(self, x - tau * self.b, tau) diff --git a/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py b/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py new file mode 100755 index 0000000..8a52566 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py @@ -0,0 +1,66 @@ +from numbers import Number
+import numpy
+
+class ScaledFunction(object):
+ '''ScaledFunction
+
+ A class to represent the scalar multiplication of an Function with a scalar.
+ It holds a function and a scalar. Basically it returns the multiplication
+ of the product of the function __call__, convex_conjugate and gradient with the scalar.
+ For the rest it behaves like the function it holds.
+
+ Args:
+ function (Function): a Function or BlockOperator
+ scalar (Number): a scalar multiplier
+ Example:
+ The scaled operator behaves like the following:
+
+ '''
+ def __init__(self, function, scalar):
+ super(ScaledFunction, self).__init__()
+ self.L = None
+ if not isinstance (scalar, Number):
+ raise TypeError('expected scalar: got {}'.format(type(scalar)))
+ self.scalar = scalar
+ self.function = function
+
+ def __call__(self,x, out=None):
+ '''Evaluates the function at x '''
+ return self.scalar * self.function(x)
+
+ def convex_conjugate(self, x):
+ '''returns the convex_conjugate of the scaled function '''
+ # if out is None:
+ # return self.scalar * self.function.convex_conjugate(x/self.scalar)
+ # else:
+ # out.fill(self.function.convex_conjugate(x/self.scalar))
+ # out *= self.scalar
+ return self.scalar * self.function.convex_conjugate(x/self.scalar)
+
+ def proximal_conjugate(self, x, tau, out = None):
+ '''This returns the proximal operator for the function at x, tau
+
+ TODO check if this is mathematically correct'''
+ return self.function.proximal_conjugate(x, tau, out=out)
+
+ def grad(self, x):
+ '''Alias of gradient(x,None)'''
+ warnings.warn('''This method will disappear in following
+ versions of the CIL. Use gradient instead''', DeprecationWarning)
+ return self.gradient(x, out=None)
+
+ def prox(self, x, tau):
+ '''Alias of proximal(x, tau, None)'''
+ warnings.warn('''This method will disappear in following
+ versions of the CIL. Use proximal instead''', DeprecationWarning)
+ return self.proximal(x, out=None)
+
+ def gradient(self, x, out=None):
+ '''Returns the gradient of the function at x, if the function is differentiable'''
+ return self.scalar * self.function.gradient(x, out=out)
+
+ def proximal(self, x, tau, out=None):
+ '''This returns the proximal operator for the function at x, tau
+
+ TODO check if this is mathematically correct'''
+ return self.function.proximal(x, tau, out=out)
diff --git a/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py b/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py new file mode 100644 index 0000000..9def741 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Mar 6 19:44:10 2019 + +@author: evangelos +""" + +import numpy as np +#from ccpi.optimisation.funcs import Function +from ccpi.optimisation.functions import Function +from ccpi.framework import DataContainer, ImageData +from ccpi.framework import BlockDataContainer + +class ZeroFun(Function): + + def __init__(self): + super(ZeroFun, self).__init__() + + def __call__(self,x): + return 0 + + def convex_conjugate(self, x): + ''' This is the support function sup <x, x^{*}> which in fact is the + indicator function for the set = {0} + So 0 if x=0, or inf if x neq 0 + ''' + + if x.shape[0]==1: + return x.maximum(0).sum() + else: + if isinstance(x, CompositeDataContainer): + return x.get_item(0).maximum(0).sum() + x.get_item(1).maximum(0).sum() + else: + return x.maximum(0).sum() + x.maximum(0).sum() + + def proximal(self,x,tau, out=None): + if out is None: + return x.copy() + else: + out.fill(x) + + def proximal_conjugate(self, x, tau): + return 0
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/functions/__init__.py b/Wrappers/Python/ccpi/optimisation/functions/__init__.py new file mode 100644 index 0000000..9030454 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/__init__.py @@ -0,0 +1,10 @@ +# -*- coding: utf-8 -*- + +from .Function import Function +from .ZeroFun import ZeroFun +from .L1Norm import SimpleL1Norm, L1Norm +from .L2NormSquared import L2NormSq, SimpleL2NormSq +from .mixed_L12Norm import mixed_L12Norm +from .BlockFunction import BlockFunction +from .ScaledFunction import ScaledFunction +from .FunctionOperatorComposition import FunctionOperatorComposition diff --git a/Wrappers/Python/ccpi/optimisation/functions/functions.py b/Wrappers/Python/ccpi/optimisation/functions/functions.py new file mode 100644 index 0000000..8632920 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/functions.py @@ -0,0 +1,312 @@ +# -*- coding: utf-8 -*- + +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Thu Feb 7 13:10:56 2019 + +@author: evangelos +""" + +import numpy as np +#from ccpi.optimisation.funcs import Function +from ccpi.optimisation.functions import Function +from ccpi.framework import DataContainer, ImageData, ImageGeometry +from operators import CompositeDataContainer, Identity, CompositeOperator +from numbers import Number + + +############################ L2NORM FUNCTIONS ############################# +class SimpleL2NormSq(Function): + + def __init__(self, alpha=1): + + super(SimpleL2NormSq, self).__init__() + self.alpha = alpha + + def __call__(self, x): + return self.alpha * x.power(2).sum() + + def gradient(self,x): + return 2 * self.alpha * x + + def convex_conjugate(self,x): + return (1/4*self.alpha) * x.power(2).sum() + + def proximal(self, x, tau): + return x.divide(1+2*tau*self.alpha) + + def proximal_conjugate(self, x, tau): + return x.divide(1 + tau/2*self.alpha ) + + +class L2NormSq(SimpleL2NormSq): + + def __init__(self, A, b = None, alpha=1, **kwargs): + + super(L2NormSq, self).__init__(alpha=alpha) + self.alpha = alpha + self.A = A + self.b = b + + def __call__(self, x): + + if self.b is None: + return SimpleL2NormSq.__call__(self, self.A.direct(x)) + else: + return SimpleL2NormSq.__call__(self, self.A.direct(x) - self.b) + + def convex_conjugate(self, x): + + ''' The convex conjugate corresponds to the simple functional i.e., + f(x) = alpha * ||x - b||_{2}^{2} + ''' + + if self.b is None: + return SimpleL2NormSq.convex_conjugate(self, x) + else: + return SimpleL2NormSq.convex_conjugate(self, x) + (self.b * x).sum() + + def gradient(self, x): + + if self.b is None: + return 2*self.alpha * self.A.adjoint(self.A.direct(x)) + else: + return 2*self.alpha * self.A.adjoint(self.A.direct(x) - self.b) + + def proximal(self, x, tau): + + ''' The proximal operator corresponds to the simple functional i.e., + f(x) = alpha * ||x - b||_{2}^{2} + + argmin_x { 0.5||x - u||^{2} + tau f(x) } + ''' + + if self.b is None: + return SimpleL2NormSq.proximal(self, x, tau) + else: + return self.b + SimpleL2NormSq.proximal(self, x - self.b , tau) + + + def proximal_conjugate(self, x, tau): + + ''' The proximal operator corresponds to the simple convex conjugate + functional i.e., f^{*}(x^{) + argmin_x { 0.5||x - u||^{2} + tau f(x) } + ''' + if self.b is None: + return SimpleL2NormSq.proximal_conjugate(self, x, tau) + else: + return SimpleL2NormSq.proximal_conjugate(self, x - tau * self.b, tau) + + +############################ L1NORM FUNCTIONS ############################# +class SimpleL1Norm(Function): + + def __init__(self, alpha=1): + + super(SimpleL1Norm, self).__init__() + self.alpha = alpha + + def __call__(self, x): + return self.alpha * x.abs().sum() + + def gradient(self,x): + return ValueError('Not Differentiable') + + def convex_conjugate(self,x): + return 0 + + def proximal(self, x, tau): + ''' Soft Threshold''' + return x.sign() * (x.abs() - tau * self.alpha).maximum(1.0) + + def proximal_conjugate(self, x, tau): + return x.divide((x.abs()/self.alpha).maximum(1.0)) + +class L1Norm(SimpleL1Norm): + + def __init__(self, A, b = None, alpha=1, **kwargs): + + super(L1Norm, self).__init__() + self.alpha = alpha + self.A = A + self.b = b + + def __call__(self, x): + + if self.b is None: + return SimpleL1Norm.__call__(self, self.A.direct(x)) + else: + return SimpleL1Norm.__call__(self, self.A.direct(x) - self.b) + + def gradient(self, x): + return ValueError('Not Differentiable') + + def convex_conjugate(self,x): + if self.b is None: + return SimpleL1Norm.convex_conjugate(self, x) + else: + return SimpleL1Norm.convex_conjugate(self, x) + (self.b * x).sum() + + def proximal(self, x, tau): + + if self.b is None: + return SimpleL1Norm.proximal(self, x, tau) + else: + return self.b + SimpleL1Norm.proximal(self, x + self.b , tau) + + def proximal_conjugate(self, x, tau): + + if self.b is None: + return SimpleL1Norm.proximal_conjugate(self, x, tau) + else: + return SimpleL1Norm.proximal_conjugate(self, x - tau*self.b, tau) + + +############################ mixed_L1,2NORM FUNCTIONS ############################# +class mixed_L12Norm(Function): + + def __init__(self, A, b=None, alpha=1, **kwargs): + + super(mixed_L12Norm, self).__init__() + self.alpha = alpha + self.A = A + self.b = b + + self.sym_grad = kwargs.get('sym_grad',False) + + + + def gradient(self,x): + return ValueError('Not Differentiable') + + + def __call__(self,x): + + y = self.A.direct(x) + eucl_norm = ImageData(y.power(2).sum(axis=0)).sqrt() + eucl_norm.__isub__(self.b) + return eucl_norm.sum() * self.alpha + + def convex_conjugate(self,x): + return 0 + + def proximal_conjugate(self, x, tau): + + if self.b is None: + + if self.sym_grad: + tmp2 = np.sqrt(x.as_array()[0]**2 + x.as_array()[1]**2 + 2*x.as_array()[2]**2)/self.alpha + res = x.divide(ImageData(tmp2).maximum(1.0)) + else: + res = x.divide((ImageData(x.power(2).sum(axis=0)).sqrt()/self.alpha).maximum(1.0)) + + else: + res = (x - tau*self.b)/ ((x - tau*self.b)).abs().maximum(1.0) + + return res + + +#%% + +class ZeroFun(Function): + + def __init__(self): + super(ZeroFun, self).__init__() + + def __call__(self,x): + return 0 + + def convex_conjugate(self, x): + ''' This is the support function sup <x, x^{*}> which in fact is the + indicator function for the set = {0} + So 0 if x=0, or inf if x neq 0 + ''' + return x.maximum(0).sum() + + def proximal(self,x,tau): + return x.copy() + + def proximal_conjugate(self, x, tau): + return 0 + + +class CompositeFunction(Function): + + def __init__(self, *args): + self.functions = args + self.length = len(self.functions) + + def get_item(self, ind): + return self.functions[ind] + + def __call__(self,x): + + t = 0 + for i in range(self.length): + for j in range(x.shape[0]): + t +=self.functions[i](x.get_item(j)) + return t + + def convex_conjugate(self, x): + + z = 0 + t = 0 + for i in range(x.shape[0]): + t += self.functions[z].convex_conjugate(x.get_item(i)) + z += 1 + + return t + + def proximal_conjugate(self, x, tau, out = None): + + if isinstance(tau, Number): + tau = CompositeDataContainer(tau) + out = [None]*self.length + for i in range(self.length): + out[i] = self.functions[i].proximal(x.get_item(i), tau.get_item(0)) + return CompositeDataContainer(*out) + + + def proximal_conjugate(self, x, tau, out = None): + + if isinstance(tau, Number): + tau = CompositeDataContainer(tau) + out = [None]*self.length + for i in range(self.length): + out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau.get_item(0)) + return CompositeDataContainer(*out) + + + + +if __name__ == '__main__': + + N = 3 + ig = (N,N) + ag = ig + op1 = Gradient(ig) + op2 = Identity(ig, ag) + + # Form Composite Operator + operator = CompositeOperator((2,1), op1, op2 ) + + # Create functions + alpha = 1 + noisy_data = ImageData(np.random.randint(10, size=ag)) + f = CompositeFunction(L1Norm(op1,alpha), \ + L2NormSq(op2, noisy_data, c = 0.5, memopt = False) ) + + u = ImageData(np.random.randint(10, size=ig)) + uComp = CompositeDataContainer(u) + + print(f(uComp)) # This is f(Kx) = f1(K1*u) + f2(K2*u) + + f1 = L1Norm(op1,alpha) + f2 = L2NormSq(op2, noisy_data, c = 0.5, memopt = False) + + print(f1(u) + f2(u)) + + + diff --git a/Wrappers/Python/ccpi/optimisation/functions/mixed_L12Norm.py b/Wrappers/Python/ccpi/optimisation/functions/mixed_L12Norm.py new file mode 100644 index 0000000..ffeb32e --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/mixed_L12Norm.py @@ -0,0 +1,56 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Mar 6 19:43:12 2019 + +@author: evangelos +""" + +import numpy as np +#from ccpi.optimisation.funcs import Function +from ccpi.optimisation.functions import Function +from ccpi.framework import DataContainer, ImageData, ImageGeometry + +############################ mixed_L1,2NORM FUNCTIONS ############################# +class mixed_L12Norm(Function): + + def __init__(self, alpha, **kwargs): + + super(mixed_L12Norm, self).__init__() + + self.alpha = alpha + self.b = kwargs.get('b',None) + self.sym_grad = kwargs.get('sym_grad',False) + + def __call__(self,x): + + if self.b is None: + tmp1 = x + else: + tmp1 = x - self.b +# + if self.sym_grad: + tmp = np.sqrt(tmp1.as_array()[0]**2 + tmp1.as_array()[1]**2 + 2*tmp1.as_array()[2]**2) + else: + tmp = ImageData(tmp1.power(2).sum(axis=0)).sqrt() + + return self.alpha*tmp.sum() + + def gradient(self,x): + return ValueError('Not Differentiable') + + def convex_conjugate(self,x): + return 0 + + def proximal(self, x, tau): + pass + + def proximal_conjugate(self, x, tau): + + if self.sym_grad: + tmp2 = np.sqrt(x.as_array()[0]**2 + x.as_array()[1]**2 + 2*x.as_array()[2]**2)/self.alpha + res = x.divide(ImageData(tmp2).maximum(1.0)) + else: + res = x.divide((ImageData(x.power(2).sum(axis=0)).sqrt()/self.alpha).maximum(1.0)) + + return res diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py new file mode 100755 index 0000000..4ff38c6 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py @@ -0,0 +1,157 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Feb 14 12:36:40 2019 + +@author: ofn77899 +""" +#from ccpi.optimisation.ops import Operator +import numpy +from numbers import Number +import functools +from ccpi.framework import AcquisitionData, ImageData, BlockDataContainer +from ccpi.optimisation.operators import Operator, LinearOperator +from ccpi.optimisation.operators.BlockScaledOperator import BlockScaledOperator +from ccpi.framework import BlockGeometry + +class BlockOperator(Operator): + '''Class to hold a block operator + + Class to hold a number of Operators in a block. + User may specify the shape of the block, by default is a row vector + + BlockOperators have a generic shape M x N, and when applied on an + Nx1 BlockDataContainer, will yield and Mx1 BlockDataContainer. + Notice: BlockDatacontainer are only allowed to have the shape of N x 1, with + N rows and 1 column. + ''' + __array_priority__ = 1 + def __init__(self, *args, **kwargs): + ''' + Class creator + + Note: + Do not include the `self` parameter in the ``Args`` section. + + Args: + vararg (Operator): Operators in the block. + shape (:obj:`tuple`, optional): If shape is passed the Operators in + vararg are considered input in a row-by-row fashion. + Shape and number of Operators must match. + + Example: + BlockOperator(op0,op1) results in a row block + BlockOperator(op0,op1,shape=(1,2)) results in a column block + ''' + self.operators = args + shape = kwargs.get('shape', None) + if shape is None: + shape = (len(args),1) + self.shape = shape + n_elements = functools.reduce(lambda x,y: x*y, shape, 1) + if len(args) != n_elements: + raise ValueError( + 'Dimension and size do not match: expected {} got {}' + .format(n_elements,len(args))) + def get_item(self, row, col): + if row > self.shape[0]: + raise ValueError('Requested row {} > max {}'.format(row, self.shape[0])) + if col > self.shape[1]: + raise ValueError('Requested col {} > max {}'.format(col, self.shape[1])) + + index = row*self.shape[1]+col + return self.operators[index] + + def norm(self): + norm = [op.norm() for op in self.operators] + b = [] + for i in range(self.shape[0]): + b.append([]) + for j in range(self.shape[1]): + b[-1].append(norm[i*self.shape[1]+j]) + return numpy.asarray(b) + + def direct(self, x, out=None): + '''Direct operation for the BlockOperator + + BlockOperator work on BlockDataContainer, but they will work on DataContainers + and inherited classes by simple wrapping the input in a BlockDataContainer of shape (1,1) + ''' + if not isinstance (x, BlockDataContainer): + x_b = BlockDataContainer(x) + else: + x_b = x + shape = self.get_output_shape(x_b.shape) + res = [] + for row in range(self.shape[0]): + for col in range(self.shape[1]): + if col == 0: + prod = self.get_item(row,col).direct(x_b.get_item(col)) + else: + prod += self.get_item(row,col).direct(x_b.get_item(col)) + res.append(prod) + return BlockDataContainer(*res, shape=shape) + + def adjoint(self, x, out=None): + '''Adjoint operation for the BlockOperator + + BlockOperator may contain both LinearOperator and Operator + This method exists in BlockOperator as it is not known what type of + Operator it will contain. + + BlockOperator work on BlockDataContainer, but they will work on DataContainers + and inherited classes by simple wrapping the input in a BlockDataContainer of shape (1,1) + + Raises: ValueError if the contained Operators are not linear + ''' + if not functools.reduce(lambda x, y: x and y.is_linear(), self.operators, True): + raise ValueError('Not all operators in Block are linear.') + if not isinstance (x, BlockDataContainer): + x_b = BlockDataContainer(x) + else: + x_b = x + shape = self.get_output_shape(x_b.shape, adjoint=True) + res = [] + for row in range(self.shape[1]): + for col in range(self.shape[0]): + if col == 0: + prod = self.get_item(row, col).adjoint(x_b.get_item(col)) + else: + prod += self.get_item(row, col).adjoint(x_b.get_item(col)) + res.append(prod) + return BlockDataContainer(*res, shape=shape) + + def get_output_shape(self, xshape, adjoint=False): + sshape = self.shape[1] + oshape = self.shape[0] + if adjoint: + sshape = self.shape[0] + oshape = self.shape[1] + if sshape != xshape[0]: + raise ValueError('Incompatible shapes {} {}'.format(self.shape, xshape)) + return (oshape, xshape[-1]) + + def __rmul__(self, scalar): + '''Defines the left multiplication with a scalar + + Args: scalar (number or iterable containing numbers): + + Returns: a block operator with Scaled Operators inside''' + if isinstance (scalar, list) or isinstance(scalar, tuple) or \ + isinstance(scalar, numpy.ndarray): + if len(scalar) != len(self.operators): + raise ValueError('dimensions of scalars and operators do not match') + scalars = scalar + else: + scalars = [scalar for _ in self.operators] + # create a list of ScaledOperator-s + ops = [ v * op for v,op in zip(scalars, self.operators)] + #return BlockScaledOperator(self, scalars ,shape=self.shape) + return type(self)(*ops, shape=self.shape) + @property + def T(self): + '''Return the transposed of self''' + shape = (self.shape[1], self.shape[0]) + return type(self)(*self.operators, shape=shape) + +if __name__ == '__main__': + pass diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py new file mode 100644 index 0000000..aeb6c53 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py @@ -0,0 +1,67 @@ +from numbers import Number +import numpy +from ccpi.optimisation.operators import ScaledOperator +import functools + +class BlockScaledOperator(ScaledOperator): + '''ScaledOperator + + A class to represent the scalar multiplication of an Operator with a scalar. + It holds an operator and a scalar. Basically it returns the multiplication + of the result of direct and adjoint of the operator with the scalar. + For the rest it behaves like the operator it holds. + + Args: + operator (Operator): a Operator or LinearOperator + scalar (Number): a scalar multiplier + Example: + The scaled operator behaves like the following: + sop = ScaledOperator(operator, scalar) + sop.direct(x) = scalar * operator.direct(x) + sop.adjoint(x) = scalar * operator.adjoint(x) + sop.norm() = operator.norm() + sop.range_geometry() = operator.range_geometry() + sop.domain_geometry() = operator.domain_geometry() + ''' + def __init__(self, operator, scalar, shape=None): + if shape is None: + shape = operator.shape + + if isinstance(scalar, (list, tuple, numpy.ndarray)): + size = functools.reduce(lambda x,y:x*y, shape, 1) + if len(scalar) != size: + raise ValueError('Scalar and operators size do not match: {}!={}' + .format(len(scalar), len(operator))) + self.scalar = scalar[:] + print ("BlockScaledOperator ", self.scalar) + elif isinstance (scalar, Number): + self.scalar = scalar + else: + raise TypeError('expected scalar to be a number of an iterable: got {}'.format(type(scalar))) + self.operator = operator + self.shape = shape + def direct(self, x, out=None): + print ("BlockScaledOperator self.scalar", self.scalar) + #print ("self.scalar", self.scalar[0]* x.get_item(0).as_array()) + return self.scalar * (self.operator.direct(x, out=out)) + def adjoint(self, x, out=None): + if self.operator.is_linear(): + return self.scalar * self.operator.adjoint(x, out=out) + else: + raise TypeError('Operator is not linear') + def norm(self): + return numpy.abs(self.scalar) * self.operator.norm() + def range_geometry(self): + return self.operator.range_geometry() + def domain_geometry(self): + return self.operator.domain_geometry() + @property + def T(self): + '''Return the transposed of self''' + #print ("transpose before" , self.shape) + #shape = (self.shape[1], self.shape[0]) + ##self.shape = shape + ##self.operator.shape = shape + #print ("transpose" , shape) + #return self + return type(self)(self.operator.T, self.scalar)
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py new file mode 100644 index 0000000..24c4e4b --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py @@ -0,0 +1,322 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 1 22:51:17 2019 + +@author: evangelos +""" + +from ccpi.optimisation.operators import Operator +from ccpi.optimisation.ops import PowerMethodNonsquare +from ccpi.framework import ImageData, BlockDataContainer +import numpy as np + +class FiniteDiff(Operator): + + # Works for Neum/Symmetric & periodic boundary conditions + # TODO add central differences??? + # TODO not very well optimised, too many conditions + # TODO add discretisation step, should get that from imageGeometry + + # Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x'] + # Grad_order = ['channels', 'direction_y', 'direction_x'] + # Grad_order = ['direction_z', 'direction_y', 'direction_x'] + # Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x'] + + def __init__(self, gm_domain, gm_range=None, direction=0, bnd_cond = 'Neumann'): + '''''' + super(FiniteDiff, self).__init__() + '''FIXME: domain and range should be geometries''' + self.gm_domain = gm_domain + self.gm_range = gm_range + self.direction = direction + self.bnd_cond = bnd_cond + + # Domain Geometry = Range Geometry if not stated + if self.gm_range is None: + self.gm_range = self.gm_domain + # check direction and "length" of geometry + if self.direction + 1 > len(self.gm_domain.shape): + raise ValueError('Gradient directions more than geometry domain') + + #self.voxel_size = kwargs.get('voxel_size',1) + # this wrongly assumes a homogeneous voxel size + self.voxel_size = self.gm_domain.voxel_size_x + + + def direct(self, x, out=None): + + x_asarr = x.as_array() + x_sz = len(x.shape) + + if out is None: + out = np.zeros(x.shape) + + fd_arr = out + + ######################## Direct for 2D ############################### + if x_sz == 2: + + if self.direction == 1: + + np.subtract( x_asarr[:,1:], x_asarr[:,0:-1], out = fd_arr[:,0:-1] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0], x_asarr[:,-1], out = fd_arr[:,-1] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 0: + + np.subtract( x_asarr[1:], x_asarr[0:-1], out = fd_arr[0:-1,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:], x_asarr[-1,:], out = fd_arr[-1,:] ) + else: + raise ValueError('No valid boundary conditions') + + ######################## Direct for 3D ############################### + elif x_sz == 3: + + if self.direction == 0: + + np.subtract( x_asarr[1:,:,:], x_asarr[0:-1,:,:], out = fd_arr[0:-1,:,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:,:], x_asarr[-1,:,:], out = fd_arr[-1,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 1: + + np.subtract( x_asarr[:,1:,:], x_asarr[:,0:-1,:], out = fd_arr[:,0:-1,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0,:], x_asarr[:,-1,:], out = fd_arr[:,-1,:] ) + else: + raise ValueError('No valid boundary conditions') + + + if self.direction == 2: + + np.subtract( x_asarr[:,:,1:], x_asarr[:,:,0:-1], out = fd_arr[:,:,0:-1] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,0], x_asarr[:,:,-1], out = fd_arr[:,:,-1] ) + else: + raise ValueError('No valid boundary conditions') + + ######################## Direct for 4D ############################### + elif x_sz == 4: + + if self.direction == 0: + np.subtract( x_asarr[1:,:,:,:], x_asarr[0:-1,:,:,:], out = fd_arr[0:-1,:,:,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:,:,:], x_asarr[-1,:,:,:], out = fd_arr[-1,:,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 1: + np.subtract( x_asarr[:,1:,:,:], x_asarr[:,0:-1,:,:], out = fd_arr[:,0:-1,:,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0,:,:], x_asarr[:,-1,:,:], out = fd_arr[:,-1,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 2: + np.subtract( x_asarr[:,:,1:,:], x_asarr[:,:,0:-1,:], out = fd_arr[:,:,0:-1,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,0,:], x_asarr[:,:,-1,:], out = fd_arr[:,:,-1,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 3: + np.subtract( x_asarr[:,:,:,1:], x_asarr[:,:,:,0:-1], out = fd_arr[:,:,:,0:-1] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,:,0], x_asarr[:,:,:,-1], out = fd_arr[:,:,:,-1] ) + else: + raise ValueError('No valid boundary conditions') + + else: + raise NotImplementedError + + res = out/self.voxel_size + return type(x)(res) + + def adjoint(self, x, out=None): + + x_asarr = x.as_array() + #x_asarr = x + x_sz = len(x.shape) + + if out is None: + out = np.zeros(x.shape) + + fd_arr = out + + ######################## Adjoint for 2D ############################### + if x_sz == 2: + + if self.direction == 1: + + np.subtract( x_asarr[:,1:], x_asarr[:,0:-1], out = fd_arr[:,1:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,0], 0, out = fd_arr[:,0] ) + np.subtract( -x_asarr[:,-2], 0, out = fd_arr[:,-1] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0], x_asarr[:,-1], out = fd_arr[:,0] ) + + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 0: + + np.subtract( x_asarr[1:,:], x_asarr[0:-1,:], out = fd_arr[1:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[0,:], 0, out = fd_arr[0,:] ) + np.subtract( -x_asarr[-2,:], 0, out = fd_arr[-1,:] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:], x_asarr[-1,:], out = fd_arr[0,:] ) + + else: + raise ValueError('No valid boundary conditions') + + ######################## Adjoint for 3D ############################### + elif x_sz == 3: + + if self.direction == 0: + + np.subtract( x_asarr[1:,:,:], x_asarr[0:-1,:,:], out = fd_arr[1:,:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[0,:,:], 0, out = fd_arr[0,:,:] ) + np.subtract( -x_asarr[-2,:,:], 0, out = fd_arr[-1,:,:] ) + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:,:], x_asarr[-1,:,:], out = fd_arr[0,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 1: + np.subtract( x_asarr[:,1:,:], x_asarr[:,0:-1,:], out = fd_arr[:,1:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,0,:], 0, out = fd_arr[:,0,:] ) + np.subtract( -x_asarr[:,-2,:], 0, out = fd_arr[:,-1,:] ) + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0,:], x_asarr[:,-1,:], out = fd_arr[:,0,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 2: + np.subtract( x_asarr[:,:,1:], x_asarr[:,:,0:-1], out = fd_arr[:,:,1:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,:,0], 0, out = fd_arr[:,:,0] ) + np.subtract( -x_asarr[:,:,-2], 0, out = fd_arr[:,:,-1] ) + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,0], x_asarr[:,:,-1], out = fd_arr[:,:,0] ) + else: + raise ValueError('No valid boundary conditions') + + ######################## Adjoint for 4D ############################### + elif x_sz == 4: + + if self.direction == 0: + np.subtract( x_asarr[1:,:,:,:], x_asarr[0:-1,:,:,:], out = fd_arr[1:,:,:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[0,:,:,:], 0, out = fd_arr[0,:,:,:] ) + np.subtract( -x_asarr[-2,:,:,:], 0, out = fd_arr[-1,:,:,:] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:,:,:], x_asarr[-1,:,:,:], out = fd_arr[0,:,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 1: + np.subtract( x_asarr[:,1:,:,:], x_asarr[:,0:-1,:,:], out = fd_arr[:,1:,:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,0,:,:], 0, out = fd_arr[:,0,:,:] ) + np.subtract( -x_asarr[:,-2,:,:], 0, out = fd_arr[:,-1,:,:] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0,:,:], x_asarr[:,-1,:,:], out = fd_arr[:,0,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + + if self.direction == 2: + np.subtract( x_asarr[:,:,1:,:], x_asarr[:,:,0:-1,:], out = fd_arr[:,:,1:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,:,0,:], 0, out = fd_arr[:,:,0,:] ) + np.subtract( -x_asarr[:,:,-2,:], 0, out = fd_arr[:,:,-1,:] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,0,:], x_asarr[:,:,-1,:], out = fd_arr[:,:,0,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 3: + np.subtract( x_asarr[:,:,:,1:], x_asarr[:,:,:,0:-1], out = fd_arr[:,:,:,1:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,:,:,0], 0, out = fd_arr[:,:,:,0] ) + np.subtract( -x_asarr[:,:,:,-2], 0, out = fd_arr[:,:,:,-1] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,:,0], x_asarr[:,:,:,-1], out = fd_arr[:,:,:,0] ) + else: + raise ValueError('No valid boundary conditions') + + else: + raise NotImplementedError + + res = out/self.voxel_size + return type(x)(-res) + + def range_geometry(self): + '''Returns the range geometry''' + return self.gm_range + + def domain_geometry(self): + '''Returns the domain geometry''' + return self.gm_domain + + def norm(self): + x0 = self.gm_domain.allocate() + x0.fill( np.random.random_sample(x0.shape) ) + self.s1, sall, svec = PowerMethodNonsquare(self, 25, x0) + return self.s1 + + + + +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py new file mode 100644 index 0000000..d0d0f43 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py @@ -0,0 +1,136 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 1 22:50:04 2019 + +@author: evangelos +""" + +from ccpi.optimisation.operators import Operator +from ccpi.optimisation.ops import PowerMethodNonsquare +from ccpi.framework import ImageData, BlockDataContainer +import numpy as np +from ccpi.optimisation.operators import FiniteDiff +from ccpi.framework import BlockGeometry + +#%% + +class Gradient(Operator): + + def __init__(self, gm_domain, gm_range=None, bnd_cond = 'Neumann', **kwargs): + + super(Gradient, self).__init__() + + self.gm_domain = gm_domain # Domain of Grad Operator + self.gm_range = gm_range # Range of Grad Operator + self.bnd_cond = bnd_cond # Boundary conditions of Finite Differences + + + if self.gm_range is None: + #FIXME this should be a BlockGeometry + self.gm_range = ((len(self.gm_domain),)+self.gm_domain) + + # Kwargs Default options + self.memopt = kwargs.get('memopt',False) + self.correlation = kwargs.get('correlation','Space') + + #TODO not tested yet, operator norm??? + self.voxel_size = kwargs.get('voxel_size',[1]*len(gm_domain)) + + + def direct(self, x, out=None): + + #tmp = np.zeros(self.gm_range) + tmp = self.gm_range.allocate() + for i in range(len(self.gm_domain)): + #tmp[i] = FiniteDiff(self.gm_domain, direction = i, bnd_cond = self.bnd_cond).direct(x.as_array())/self.voxel_size[i] + if self.correlation == 'Space': + if i == 0 : + i+=1 + tmp[i].fill( + FiniteDiff(self.gm_domain, direction = i, bnd_cond = self.bnd_cond).direct(x.as_array())/self.voxel_size[i] + ) +# return type(x)(tmp) + return type(x)(tmp) + + def adjoint(self, x, out=None): + + tmp = np.zeros(self.gm_domain) + for i in range(len(self.gm_domain)): + tmp+=FiniteDiff(self.gm_domain, direction = i, bnd_cond = self.bnd_cond).adjoint(x.as_array()[i])/self.voxel_size[i] + return type(x)(-tmp) + + def alloc_domain_dim(self): + return ImageData(np.zeros(self.gm_domain)) + + def alloc_range_dim(self): + return ImageData(np.zeros(self.range_dim)) + + def domain_geometry(self): + return self.gm_domain + + def range_geometry(self): + '''fix this''' + return BlockGeometry(self.gm_range, self.gm_range) + + def norm(self): +# return np.sqrt(4*len(self.domainDim())) + #TODO this takes time for big ImageData + # for 2D ||grad|| = sqrt(8), 3D ||grad|| = sqrt(12) + x0 = ImageData(np.random.random_sample(self.domain_dim())) + self.s1, sall, svec = PowerMethodNonsquare(self, 25, x0) + return self.s1 + + +if __name__ == '__main__': + + N, M = (200,300) + ig = (N,M) + G = Gradient(ig) + u = DataContainer(np.random.randint(10, size=G.domain_dim())) + w = DataContainer(np.random.randint(10, size=G.range_dim())) +# w = [DataContainer(np.random.randint(10, size=G.domain_dim())),\ +# DataContainer(np.random.randint(10, size=G.domain_dim()))] + + # domain_dim + print('Domain {}'.format(G.domain_geometry())) + + # range_dim + print('Range {}'.format(G.range_geometry())) + + # Direct + z = G.direct(u) + + # Adjoint + z1 = G.adjoint(w) + + print(z) + print(z1) + + LHS = (G.direct(u)*w).sum() + RHS = (u * G.adjoint(w)).sum() +# + print(LHS,RHS) + print(G.norm()) + +# print(G.adjoint(G.direct(u))) + + + + + + + + + + + + + + + + + + + +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py new file mode 100644 index 0000000..d49cb30 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Mar 6 19:30:51 2019 + +@author: evangelos +""" + +from ccpi.optimisation.operators import Operator + + +class Identity(Operator): + + def __init__(self, gm_domain, gm_range=None): + + self.gm_domain = gm_domain + self.gm_range = gm_range + if self.gm_range is None: + self.gm_range = self.gm_domain + + super(Identity, self).__init__() + + def direct(self,x,out=None): + if out is None: + return x.copy() + else: + out.fill(x) + + def adjoint(self,x, out=None): + if out is None: + return x.copy() + else: + out.fill(x) + + def norm(self): + return 1.0 + + def domain_dim(self): + return self.gm_domain + + def range_dim(self): + return self.gm_range
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py new file mode 100755 index 0000000..e19304f --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py @@ -0,0 +1,22 @@ +# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 5 15:57:52 2019
+
+@author: ofn77899
+"""
+
+from ccpi.optimisation.operators import Operator
+
+
+class LinearOperator(Operator):
+ '''A Linear Operator that maps from a space X <-> Y'''
+ def __init__(self):
+ super(LinearOperator, self).__init__()
+ def is_linear(self):
+ '''Returns if the operator is linear'''
+ return True
+ def adjoint(self,x, out=None):
+ '''returns the adjoint/inverse operation
+
+ only available to linear operators'''
+ raise NotImplementedError
diff --git a/Wrappers/Python/ccpi/optimisation/operators/Operator.py b/Wrappers/Python/ccpi/optimisation/operators/Operator.py new file mode 100755 index 0000000..cdf15a7 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/Operator.py @@ -0,0 +1,30 @@ +# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 5 15:55:56 2019
+
+@author: ofn77899
+"""
+from ccpi.optimisation.operators import ScaledOperator
+
+class Operator(object):
+ '''Operator that maps from a space X -> Y'''
+ def is_linear(self):
+ '''Returns if the operator is linear'''
+ return False
+ def direct(self,x, out=None):
+ '''Returns the application of the Operator on x'''
+ raise NotImplementedError
+ def norm(self):
+ '''Returns the norm of the Operator'''
+ raise NotImplementedError
+ def range_geometry(self):
+ '''Returns the range of the Operator: Y space'''
+ raise NotImplementedError
+ def domain_geometry(self):
+ '''Returns the domain of the Operator: X space'''
+ raise NotImplementedError
+ def __rmul__(self, scalar):
+ '''Defines the multiplication by a scalar on the left
+
+ returns a ScaledOperator'''
+ return ScaledOperator(self, scalar)
diff --git a/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py new file mode 100644 index 0000000..adcc6d9 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py @@ -0,0 +1,42 @@ +from numbers import Number +import numpy + +class ScaledOperator(object): + '''ScaledOperator + + A class to represent the scalar multiplication of an Operator with a scalar. + It holds an operator and a scalar. Basically it returns the multiplication + of the result of direct and adjoint of the operator with the scalar. + For the rest it behaves like the operator it holds. + + Args: + operator (Operator): a Operator or LinearOperator + scalar (Number): a scalar multiplier + Example: + The scaled operator behaves like the following: + sop = ScaledOperator(operator, scalar) + sop.direct(x) = scalar * operator.direct(x) + sop.adjoint(x) = scalar * operator.adjoint(x) + sop.norm() = operator.norm() + sop.range_geometry() = operator.range_geometry() + sop.domain_geometry() = operator.domain_geometry() + ''' + def __init__(self, operator, scalar): + super(ScaledOperator, self).__init__() + if not isinstance (scalar, Number): + raise TypeError('expected scalar: got {}'.format(type(scalar))) + self.scalar = scalar + self.operator = operator + def direct(self, x, out=None): + return self.scalar * self.operator.direct(x, out=out) + def adjoint(self, x, out=None): + if self.operator.is_linear(): + return self.scalar * self.operator.adjoint(x, out=out) + else: + raise TypeError('Operator is not linear') + def norm(self): + return numpy.abs(self.scalar) * self.operator.norm() + def range_geometry(self): + return self.operator.range_geometry() + def domain_geometry(self): + return self.operator.domain_geometry() diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py new file mode 100644 index 0000000..d908e49 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 1 22:53:55 2019 + +@author: evangelos +""" + +from ccpi.optimisation.operators import Operator +from ccpi.optimisation.operators import FiniteDiff +from ccpi.optimisation.ops import PowerMethodNonsquare +from ccpi.framework import ImageData, DataContainer +import numpy as np + + +class SymmetrizedGradient(Operator): + + def __init__(self, gm_domain, gm_range, bnd_cond = 'Neumann', **kwargs): + + super(SymmetrizedGradient, self).__init__() + + self.gm_domain = gm_domain # Domain of Grad Operator + self.gm_range = gm_range # Range of Grad Operator + self.bnd_cond = bnd_cond # Boundary conditions of Finite Differences + + # Kwargs Default options + self.memopt = kwargs.get('memopt',False) + self.correlation = kwargs.get('correlation','Space') + + #TODO not tested yet, operator norm??? + self.voxel_size = kwargs.get('voxel_size',[1]*len(gm_domain)) + + + def direct(self, x, out=None): + + tmp = np.zeros(self.gm_range) + tmp[0] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).adjoint(x.as_array()[0]) + tmp[1] = FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).adjoint(x.as_array()[1]) + tmp[2] = 0.5 * (FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).adjoint(x.as_array()[0]) + + FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).adjoint(x.as_array()[1]) ) + + return type(x)(tmp) + + + def adjoint(self, x, out=None): + + tmp = np.zeros(self.gm_domain) + + tmp[0] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).direct(x.as_array()[0]) + \ + FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).direct(x.as_array()[2]) + + tmp[1] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).direct(x.as_array()[2]) + \ + FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).direct(x.as_array()[1]) + + return type(x)(tmp) + + def alloc_domain_dim(self): + return ImageData(np.zeros(self.gm_domain)) + + def alloc_range_dim(self): + return ImageData(np.zeros(self.range_dim)) + + def domain_dim(self): + return self.gm_domain + + def range_dim(self): + return self.gm_range + + def norm(self): +# return np.sqrt(4*len(self.domainDim())) + #TODO this takes time for big ImageData + # for 2D ||grad|| = sqrt(8), 3D ||grad|| = sqrt(12) + x0 = ImageData(np.random.random_sample(self.domain_dim())) + self.s1, sall, svec = PowerMethodNonsquare(self, 25, x0) + return self.s1 + + + +if __name__ == '__main__': + + ########################################################################### + ## Symmetrized Gradient + + N, M = 2, 3 + ig = (N,M) + ig2 = (2,) + ig + ig3 = (3,) + ig + u1 = DataContainer(np.random.randint(10, size=ig2)) + w1 = DataContainer(np.random.randint(10, size=ig3)) + + E = SymmetrizedGradient(ig2,ig3) + + d1 = E.direct(u1) + d2 = E.adjoint(w1) + + LHS = (d1.as_array()[0]*w1.as_array()[0] + \ + d1.as_array()[1]*w1.as_array()[1] + \ + 2*d1.as_array()[2]*w1.as_array()[2]).sum() + + RHS = (u1.as_array()[0]*d2.as_array()[0] + \ + u1.as_array()[1]*d2.as_array()[1]).sum() + + + print(LHS, RHS, E.norm()) + + +# + + + + + + + + + + +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py new file mode 100644 index 0000000..a7c5f09 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py @@ -0,0 +1,39 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Mar 6 19:25:53 2019 + +@author: evangelos +""" + +import numpy as np +from ccpi.framework import ImageData +from ccpi.optimisation.operators import Operator + +class ZeroOp(Operator): + + def __init__(self, gm_domain, gm_range): + self.gm_domain = gm_domain + self.gm_range = gm_range + super(ZeroOp, self).__init__() + + def direct(self,x,out=None): + if out is None: + return ImageData(np.zeros(self.gm_range)) + else: + return ImageData(np.zeros(self.gm_range)) + + def adjoint(self,x, out=None): + if out is None: + return ImageData(np.zeros(self.gm_domain)) + else: + return ImageData(np.zeros(self.gm_domain)) + + def norm(self): + return 0 + + def domain_dim(self): + return self.gm_domain + + def range_dim(self): + return self.gm_range
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/__init__.py b/Wrappers/Python/ccpi/optimisation/operators/__init__.py new file mode 100755 index 0000000..1e86efc --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/__init__.py @@ -0,0 +1,20 @@ +# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 5 15:56:27 2019
+
+@author: ofn77899
+"""
+
+from .Operator import Operator
+from .LinearOperator import LinearOperator
+from .ScaledOperator import ScaledOperator
+from .BlockOperator import BlockOperator
+from .BlockScaledOperator import BlockScaledOperator
+
+
+from .FiniteDifferenceOperator import FiniteDiff
+from .GradientOperator import Gradient
+from .SymmetrizedGradientOperator import SymmetrizedGradient
+from .IdentityOperator import Identity
+from .ZeroOperator import ZeroOp
+
diff --git a/Wrappers/Python/ccpi/optimisation/ops.py b/Wrappers/Python/ccpi/optimisation/ops.py index e9e7f44..6afb97a 100755 --- a/Wrappers/Python/ccpi/optimisation/ops.py +++ b/Wrappers/Python/ccpi/optimisation/ops.py @@ -49,9 +49,9 @@ class Operator(object): def allocate_adjoint(self): '''Allocates memory on the X space''' raise NotImplementedError - def range_dim(self): + def range_geometry(self): raise NotImplementedError - def domain_dim(self): + def domain_geometry(self): raise NotImplementedError def __rmul__(self, other): '''reverse multiplication of Operator with number sets the variable scalar in the Operator''' @@ -97,7 +97,8 @@ class TomoIdentity(Operator): self.s1 = 1.0 self.geometry = geometry - + def is_linear(self): + return True def direct(self,x,out=None): if out is None: @@ -128,6 +129,10 @@ class TomoIdentity(Operator): raise ValueError("Wrong geometry type: expected ImageGeometry of AcquisitionGeometry, got ", type(self.geometry)) def allocate_adjoint(self): return self.allocate_direct() + def range_geometry(self): + return self.geometry + def domain_geometry(self): + return self.geometry diff --git a/Wrappers/Python/ccpi/processors.py b/Wrappers/Python/ccpi/processors.py index 3a3671a..ccef410 100755 --- a/Wrappers/Python/ccpi/processors.py +++ b/Wrappers/Python/ccpi/processors.py @@ -1,514 +1,514 @@ -# -*- 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
-
-from ccpi.framework import DataProcessor, DataContainer, AcquisitionData,\
- AcquisitionGeometry, ImageGeometry, ImageData
-from ccpi.reconstruction.parallelbeam import alg as pbalg
-import numpy
-from scipy import ndimage
-
-import matplotlib.pyplot as plt
-
-
-class Normalizer(DataProcessor):
- '''Normalization based on flat and dark
-
- This processor read in a AcquisitionData and normalises it based on
- the instrument reading with and without incident photons or neutrons.
-
- Input: AcquisitionData
- Parameter: 2D projection with flat field (or stack)
- 2D projection with dark field (or stack)
- Output: AcquisitionDataSetn
- '''
-
- def __init__(self, flat_field = None, dark_field = None, tolerance = 1e-5):
- kwargs = {
- 'flat_field' : flat_field,
- 'dark_field' : dark_field,
- # very small number. Used when there is a division by zero
- 'tolerance' : tolerance
- }
-
- #DataProcessor.__init__(self, **kwargs)
- super(Normalizer, self).__init__(**kwargs)
- if not flat_field is None:
- self.set_flat_field(flat_field)
- if not dark_field is None:
- self.set_dark_field(dark_field)
-
- def check_input(self, dataset):
- if dataset.number_of_dimensions == 3 or\
- dataset.number_of_dimensions == 2:
- return True
- else:
- raise ValueError("Expected input dimensions is 2 or 3, got {0}"\
- .format(dataset.number_of_dimensions))
-
- def set_dark_field(self, df):
- if type(df) is numpy.ndarray:
- if len(numpy.shape(df)) == 3:
- raise ValueError('Dark Field should be 2D')
- elif len(numpy.shape(df)) == 2:
- self.dark_field = df
- elif issubclass(type(df), DataContainer):
- self.dark_field = self.set_dark_field(df.as_array())
-
- def set_flat_field(self, df):
- if type(df) is numpy.ndarray:
- if len(numpy.shape(df)) == 3:
- raise ValueError('Flat Field should be 2D')
- elif len(numpy.shape(df)) == 2:
- self.flat_field = df
- elif issubclass(type(df), DataContainer):
- self.flat_field = self.set_flat_field(df.as_array())
-
- @staticmethod
- def normalize_projection(projection, flat, dark, tolerance):
- a = (projection - dark)
- b = (flat-dark)
- with numpy.errstate(divide='ignore', invalid='ignore'):
- c = numpy.true_divide( a, b )
- c[ ~ numpy.isfinite( c )] = tolerance # set to not zero if 0/0
- return c
-
- @staticmethod
- def estimate_normalised_error(projection, flat, dark, delta_flat, delta_dark):
- '''returns the estimated relative error of the normalised projection
-
- n = (projection - dark) / (flat - dark)
- Dn/n = (flat-dark + projection-dark)/((flat-dark)*(projection-dark))*(Df/f + Dd/d)
- '''
- a = (projection - dark)
- b = (flat-dark)
- df = delta_flat / flat
- dd = delta_dark / dark
- rel_norm_error = (b + a) / (b * a) * (df + dd)
- return rel_norm_error
-
- def process(self, out=None):
-
- projections = self.get_input()
- dark = self.dark_field
- flat = self.flat_field
-
- if projections.number_of_dimensions == 3:
- if not (projections.shape[1:] == dark.shape and \
- projections.shape[1:] == flat.shape):
- raise ValueError('Flats/Dark and projections size do not match.')
-
-
- a = numpy.asarray(
- [ Normalizer.normalize_projection(
- projection, flat, dark, self.tolerance) \
- for projection in projections.as_array() ]
- )
- elif projections.number_of_dimensions == 2:
- a = Normalizer.normalize_projection(projections.as_array(),
- flat, dark, self.tolerance)
- y = type(projections)( a , True,
- dimension_labels=projections.dimension_labels,
- geometry=projections.geometry)
- return y
-
-
-class CenterOfRotationFinder(DataProcessor):
- '''Processor to find the center of rotation in a parallel beam experiment
-
- This processor read in a AcquisitionDataSet and finds the center of rotation
- based on Nghia Vo's method. https://doi.org/10.1364/OE.22.019078
-
- Input: AcquisitionDataSet
-
- Output: float. center of rotation in pixel coordinate
- '''
-
- def __init__(self):
- kwargs = {
-
- }
-
- #DataProcessor.__init__(self, **kwargs)
- super(CenterOfRotationFinder, self).__init__(**kwargs)
-
- def check_input(self, dataset):
- if dataset.number_of_dimensions == 3:
- if dataset.geometry.geom_type == 'parallel':
- return True
- else:
- raise ValueError('{0} is suitable only for parallel beam geometry'\
- .format(self.__class__.__name__))
- else:
- raise ValueError("Expected input dimensions is 3, got {0}"\
- .format(dataset.number_of_dimensions))
-
-
- # #########################################################################
- # Copyright (c) 2015, UChicago Argonne, LLC. All rights reserved. #
- # #
- # Copyright 2015. UChicago Argonne, LLC. This software was produced #
- # under U.S. Government contract DE-AC02-06CH11357 for Argonne National #
- # Laboratory (ANL), which is operated by UChicago Argonne, LLC for the #
- # U.S. Department of Energy. The U.S. Government has rights to use, #
- # reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR #
- # UChicago Argonne, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR #
- # ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is #
- # modified to produce derivative works, such modified software should #
- # be clearly marked, so as not to confuse it with the version available #
- # from ANL. #
- # #
- # Additionally, redistribution and use in source and binary forms, with #
- # or without modification, are permitted provided that the following #
- # conditions are met: #
- # #
- # * Redistributions of source code must retain the above copyright #
- # notice, this list of conditions and the following disclaimer. #
- # #
- # * Redistributions in binary form must reproduce the above copyright #
- # notice, this list of conditions and the following disclaimer in #
- # the documentation and/or other materials provided with the #
- # distribution. #
- # #
- # * Neither the name of UChicago Argonne, LLC, Argonne National #
- # Laboratory, ANL, the U.S. Government, nor the names of its #
- # contributors may be used to endorse or promote products derived #
- # from this software without specific prior written permission. #
- # #
- # THIS SOFTWARE IS PROVIDED BY UChicago Argonne, LLC AND CONTRIBUTORS #
- # "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT #
- # LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS #
- # FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL UChicago #
- # Argonne, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, #
- # INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, #
- # BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; #
- # LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER #
- # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT #
- # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN #
- # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE #
- # POSSIBILITY OF SUCH DAMAGE. #
- # #########################################################################
-
- @staticmethod
- def as_ndarray(arr, dtype=None, copy=False):
- if not isinstance(arr, numpy.ndarray):
- arr = numpy.array(arr, dtype=dtype, copy=copy)
- return arr
-
- @staticmethod
- def as_dtype(arr, dtype, copy=False):
- if not arr.dtype == dtype:
- arr = numpy.array(arr, dtype=dtype, copy=copy)
- return arr
-
- @staticmethod
- def as_float32(arr):
- arr = CenterOfRotationFinder.as_ndarray(arr, numpy.float32)
- return CenterOfRotationFinder.as_dtype(arr, numpy.float32)
-
-
-
-
- @staticmethod
- def find_center_vo(tomo, ind=None, smin=-40, smax=40, srad=10, step=0.5,
- ratio=2., drop=20):
- """
- Find rotation axis location using Nghia Vo's method. :cite:`Vo:14`.
-
- Parameters
- ----------
- tomo : ndarray
- 3D tomographic data.
- ind : int, optional
- Index of the slice to be used for reconstruction.
- smin, smax : int, optional
- Reference to the horizontal center of the sinogram.
- srad : float, optional
- Fine search radius.
- step : float, optional
- Step of fine searching.
- ratio : float, optional
- The ratio between the FOV of the camera and the size of object.
- It's used to generate the mask.
- drop : int, optional
- Drop lines around vertical center of the mask.
-
- Returns
- -------
- float
- Rotation axis location.
-
- Notes
- -----
- The function may not yield a correct estimate, if:
-
- - the sample size is bigger than the field of view of the camera.
- In this case the ``ratio`` argument need to be set larger
- than the default of 2.0.
-
- - there is distortion in the imaging hardware. If there's
- no correction applied, the center of the projection image may
- yield a better estimate.
-
- - the sample contrast is weak. Paganin's filter need to be applied
- to overcome this.
-
- - the sample was changed during the scan.
- """
- tomo = CenterOfRotationFinder.as_float32(tomo)
-
- if ind is None:
- ind = tomo.shape[1] // 2
- _tomo = tomo[:, ind, :]
-
-
-
- # Reduce noise by smooth filters. Use different filters for coarse and fine search
- _tomo_cs = ndimage.filters.gaussian_filter(_tomo, (3, 1))
- _tomo_fs = ndimage.filters.median_filter(_tomo, (2, 2))
-
- # Coarse and fine searches for finding the rotation center.
- if _tomo.shape[0] * _tomo.shape[1] > 4e6: # If data is large (>2kx2k)
- #_tomo_coarse = downsample(numpy.expand_dims(_tomo_cs,1), level=2)[:, 0, :]
- #init_cen = _search_coarse(_tomo_coarse, smin, smax, ratio, drop)
- #fine_cen = _search_fine(_tomo_fs, srad, step, init_cen*4, ratio, drop)
- init_cen = CenterOfRotationFinder._search_coarse(_tomo_cs, smin,
- smax, ratio, drop)
- fine_cen = CenterOfRotationFinder._search_fine(_tomo_fs, srad,
- step, init_cen,
- ratio, drop)
- else:
- init_cen = CenterOfRotationFinder._search_coarse(_tomo_cs,
- smin, smax,
- ratio, drop)
- fine_cen = CenterOfRotationFinder._search_fine(_tomo_fs, srad,
- step, init_cen,
- ratio, drop)
-
- #logger.debug('Rotation center search finished: %i', fine_cen)
- return fine_cen
-
-
- @staticmethod
- def _search_coarse(sino, smin, smax, ratio, drop):
- """
- Coarse search for finding the rotation center.
- """
- (Nrow, Ncol) = sino.shape
- centerfliplr = (Ncol - 1.0) / 2.0
-
- # Copy the sinogram and flip left right, the purpose is to
- # make a full [0;2Pi] sinogram
- _copy_sino = numpy.fliplr(sino[1:])
-
- # This image is used for compensating the shift of sinogram 2
- temp_img = numpy.zeros((Nrow - 1, Ncol), dtype='float32')
- temp_img[:] = sino[-1]
-
- # Start coarse search in which the shift step is 1
- listshift = numpy.arange(smin, smax + 1)
- listmetric = numpy.zeros(len(listshift), dtype='float32')
- mask = CenterOfRotationFinder._create_mask(2 * Nrow - 1, Ncol,
- 0.5 * ratio * Ncol, drop)
- for i in listshift:
- _sino = numpy.roll(_copy_sino, i, axis=1)
- if i >= 0:
- _sino[:, 0:i] = temp_img[:, 0:i]
- else:
- _sino[:, i:] = temp_img[:, i:]
- listmetric[i - smin] = numpy.sum(numpy.abs(numpy.fft.fftshift(
- #pyfftw.interfaces.numpy_fft.fft2(
- # numpy.vstack((sino, _sino)))
- numpy.fft.fft2(numpy.vstack((sino, _sino)))
- )) * mask)
- minpos = numpy.argmin(listmetric)
- return centerfliplr + listshift[minpos] / 2.0
-
- @staticmethod
- def _search_fine(sino, srad, step, init_cen, ratio, drop):
- """
- Fine search for finding the rotation center.
- """
- Nrow, Ncol = sino.shape
- centerfliplr = (Ncol + 1.0) / 2.0 - 1.0
- # Use to shift the sinogram 2 to the raw CoR.
- shiftsino = numpy.int16(2 * (init_cen - centerfliplr))
- _copy_sino = numpy.roll(numpy.fliplr(sino[1:]), shiftsino, axis=1)
- if init_cen <= centerfliplr:
- lefttake = numpy.int16(numpy.ceil(srad + 1))
- righttake = numpy.int16(numpy.floor(2 * init_cen - srad - 1))
- else:
- lefttake = numpy.int16(numpy.ceil(
- init_cen - (Ncol - 1 - init_cen) + srad + 1))
- righttake = numpy.int16(numpy.floor(Ncol - 1 - srad - 1))
- Ncol1 = righttake - lefttake + 1
- mask = CenterOfRotationFinder._create_mask(2 * Nrow - 1, Ncol1,
- 0.5 * ratio * Ncol, drop)
- numshift = numpy.int16((2 * srad) / step) + 1
- listshift = numpy.linspace(-srad, srad, num=numshift)
- listmetric = numpy.zeros(len(listshift), dtype='float32')
- factor1 = numpy.mean(sino[-1, lefttake:righttake])
- num1 = 0
- for i in listshift:
- _sino = ndimage.interpolation.shift(
- _copy_sino, (0, i), prefilter=False)
- factor2 = numpy.mean(_sino[0,lefttake:righttake])
- _sino = _sino * factor1 / factor2
- sinojoin = numpy.vstack((sino, _sino))
- listmetric[num1] = numpy.sum(numpy.abs(numpy.fft.fftshift(
- #pyfftw.interfaces.numpy_fft.fft2(
- # sinojoin[:, lefttake:righttake + 1])
- numpy.fft.fft2(sinojoin[:, lefttake:righttake + 1])
- )) * mask)
- num1 = num1 + 1
- minpos = numpy.argmin(listmetric)
- return init_cen + listshift[minpos] / 2.0
-
- @staticmethod
- def _create_mask(nrow, ncol, radius, drop):
- du = 1.0 / ncol
- dv = (nrow - 1.0) / (nrow * 2.0 * numpy.pi)
- centerrow = numpy.ceil(nrow / 2) - 1
- centercol = numpy.ceil(ncol / 2) - 1
- # added by Edoardo Pasca
- centerrow = int(centerrow)
- centercol = int(centercol)
- mask = numpy.zeros((nrow, ncol), dtype='float32')
- for i in range(nrow):
- num1 = numpy.round(((i - centerrow) * dv / radius) / du)
- (p1, p2) = numpy.int16(numpy.clip(numpy.sort(
- (-num1 + centercol, num1 + centercol)), 0, ncol - 1))
- mask[i, p1:p2 + 1] = numpy.ones(p2 - p1 + 1, dtype='float32')
- if drop < centerrow:
- mask[centerrow - drop:centerrow + drop + 1,
- :] = numpy.zeros((2 * drop + 1, ncol), dtype='float32')
- mask[:,centercol-1:centercol+2] = numpy.zeros((nrow, 3), dtype='float32')
- return mask
-
- def process(self, out=None):
-
- projections = self.get_input()
-
- cor = CenterOfRotationFinder.find_center_vo(projections.as_array())
-
- return cor
-
-
-class AcquisitionDataPadder(DataProcessor):
- '''Normalization based on flat and dark
-
- This processor read in a AcquisitionData and normalises it based on
- the instrument reading with and without incident photons or neutrons.
-
- Input: AcquisitionData
- Parameter: 2D projection with flat field (or stack)
- 2D projection with dark field (or stack)
- Output: AcquisitionDataSetn
- '''
-
- def __init__(self,
- center_of_rotation = None,
- acquisition_geometry = None,
- pad_value = 1e-5):
- kwargs = {
- 'acquisition_geometry' : acquisition_geometry,
- 'center_of_rotation' : center_of_rotation,
- 'pad_value' : pad_value
- }
-
- super(AcquisitionDataPadder, self).__init__(**kwargs)
-
- def check_input(self, dataset):
- if self.acquisition_geometry is None:
- self.acquisition_geometry = dataset.geometry
- if dataset.number_of_dimensions == 3:
- return True
- else:
- raise ValueError("Expected input dimensions is 2 or 3, got {0}"\
- .format(dataset.number_of_dimensions))
-
- def process(self, out=None):
- projections = self.get_input()
- w = projections.get_dimension_size('horizontal')
- delta = w - 2 * self.center_of_rotation
-
- padded_width = int (
- numpy.ceil(abs(delta)) + w
- )
- delta_pix = padded_width - w
-
- voxel_per_pixel = 1
- geom = pbalg.pb_setup_geometry_from_acquisition(projections.as_array(),
- self.acquisition_geometry.angles,
- self.center_of_rotation,
- voxel_per_pixel )
-
- padded_geometry = self.acquisition_geometry.clone()
-
- padded_geometry.pixel_num_h = geom['n_h']
- padded_geometry.pixel_num_v = geom['n_v']
-
- delta_pix_h = padded_geometry.pixel_num_h - self.acquisition_geometry.pixel_num_h
- delta_pix_v = padded_geometry.pixel_num_v - self.acquisition_geometry.pixel_num_v
-
- if delta_pix_h == 0:
- delta_pix_h = delta_pix
- padded_geometry.pixel_num_h = padded_width
- #initialize a new AcquisitionData with values close to 0
- out = AcquisitionData(geometry=padded_geometry)
- out = out + self.pad_value
-
-
- #pad in the horizontal-vertical plane -> slice on angles
- if delta > 0:
- #pad left of middle
- command = "out.array["
- for i in range(out.number_of_dimensions):
- if out.dimension_labels[i] == 'horizontal':
- value = '{0}:{1}'.format(delta_pix_h, delta_pix_h+w)
- command = command + str(value)
- else:
- if out.dimension_labels[i] == 'vertical' :
- value = '{0}:'.format(delta_pix_v)
- command = command + str(value)
- else:
- command = command + ":"
- if i < out.number_of_dimensions -1:
- command = command + ','
- command = command + '] = projections.array'
- #print (command)
- else:
- #pad right of middle
- command = "out.array["
- for i in range(out.number_of_dimensions):
- if out.dimension_labels[i] == 'horizontal':
- value = '{0}:{1}'.format(0, w)
- command = command + str(value)
- else:
- if out.dimension_labels[i] == 'vertical' :
- value = '{0}:'.format(delta_pix_v)
- command = command + str(value)
- else:
- command = command + ":"
- if i < out.number_of_dimensions -1:
- command = command + ','
- command = command + '] = projections.array'
- #print (command)
- #cleaned = eval(command)
- exec(command)
+# -*- 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 + +from ccpi.framework import DataProcessor, DataContainer, AcquisitionData,\ + AcquisitionGeometry, ImageGeometry, ImageData +from ccpi.reconstruction.parallelbeam import alg as pbalg +import numpy +from scipy import ndimage + +import matplotlib.pyplot as plt + + +class Normalizer(DataProcessor): + '''Normalization based on flat and dark + + This processor read in a AcquisitionData and normalises it based on + the instrument reading with and without incident photons or neutrons. + + Input: AcquisitionData + Parameter: 2D projection with flat field (or stack) + 2D projection with dark field (or stack) + Output: AcquisitionDataSetn + ''' + + def __init__(self, flat_field = None, dark_field = None, tolerance = 1e-5): + kwargs = { + 'flat_field' : flat_field, + 'dark_field' : dark_field, + # very small number. Used when there is a division by zero + 'tolerance' : tolerance + } + + #DataProcessor.__init__(self, **kwargs) + super(Normalizer, self).__init__(**kwargs) + if not flat_field is None: + self.set_flat_field(flat_field) + if not dark_field is None: + self.set_dark_field(dark_field) + + def check_input(self, dataset): + if dataset.number_of_dimensions == 3 or\ + dataset.number_of_dimensions == 2: + return True + else: + raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ + .format(dataset.number_of_dimensions)) + + def set_dark_field(self, df): + if type(df) is numpy.ndarray: + if len(numpy.shape(df)) == 3: + raise ValueError('Dark Field should be 2D') + elif len(numpy.shape(df)) == 2: + self.dark_field = df + elif issubclass(type(df), DataContainer): + self.dark_field = self.set_dark_field(df.as_array()) + + def set_flat_field(self, df): + if type(df) is numpy.ndarray: + if len(numpy.shape(df)) == 3: + raise ValueError('Flat Field should be 2D') + elif len(numpy.shape(df)) == 2: + self.flat_field = df + elif issubclass(type(df), DataContainer): + self.flat_field = self.set_flat_field(df.as_array()) + + @staticmethod + def normalize_projection(projection, flat, dark, tolerance): + a = (projection - dark) + b = (flat-dark) + with numpy.errstate(divide='ignore', invalid='ignore'): + c = numpy.true_divide( a, b ) + c[ ~ numpy.isfinite( c )] = tolerance # set to not zero if 0/0 + return c + + @staticmethod + def estimate_normalised_error(projection, flat, dark, delta_flat, delta_dark): + '''returns the estimated relative error of the normalised projection + + n = (projection - dark) / (flat - dark) + Dn/n = (flat-dark + projection-dark)/((flat-dark)*(projection-dark))*(Df/f + Dd/d) + ''' + a = (projection - dark) + b = (flat-dark) + df = delta_flat / flat + dd = delta_dark / dark + rel_norm_error = (b + a) / (b * a) * (df + dd) + return rel_norm_error + + def process(self, out=None): + + projections = self.get_input() + dark = self.dark_field + flat = self.flat_field + + if projections.number_of_dimensions == 3: + if not (projections.shape[1:] == dark.shape and \ + projections.shape[1:] == flat.shape): + raise ValueError('Flats/Dark and projections size do not match.') + + + a = numpy.asarray( + [ Normalizer.normalize_projection( + projection, flat, dark, self.tolerance) \ + for projection in projections.as_array() ] + ) + elif projections.number_of_dimensions == 2: + a = Normalizer.normalize_projection(projections.as_array(), + flat, dark, self.tolerance) + y = type(projections)( a , True, + dimension_labels=projections.dimension_labels, + geometry=projections.geometry) + return y + + +class CenterOfRotationFinder(DataProcessor): + '''Processor to find the center of rotation in a parallel beam experiment + + This processor read in a AcquisitionDataSet and finds the center of rotation + based on Nghia Vo's method. https://doi.org/10.1364/OE.22.019078 + + Input: AcquisitionDataSet + + Output: float. center of rotation in pixel coordinate + ''' + + def __init__(self): + kwargs = { + + } + + #DataProcessor.__init__(self, **kwargs) + super(CenterOfRotationFinder, self).__init__(**kwargs) + + def check_input(self, dataset): + if dataset.number_of_dimensions == 3: + if dataset.geometry.geom_type == 'parallel': + return True + else: + raise ValueError('{0} is suitable only for parallel beam geometry'\ + .format(self.__class__.__name__)) + else: + raise ValueError("Expected input dimensions is 3, got {0}"\ + .format(dataset.number_of_dimensions)) + + + # ######################################################################### + # Copyright (c) 2015, UChicago Argonne, LLC. All rights reserved. # + # # + # Copyright 2015. UChicago Argonne, LLC. This software was produced # + # under U.S. Government contract DE-AC02-06CH11357 for Argonne National # + # Laboratory (ANL), which is operated by UChicago Argonne, LLC for the # + # U.S. Department of Energy. The U.S. Government has rights to use, # + # reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR # + # UChicago Argonne, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR # + # ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is # + # modified to produce derivative works, such modified software should # + # be clearly marked, so as not to confuse it with the version available # + # from ANL. # + # # + # Additionally, redistribution and use in source and binary forms, with # + # or without modification, are permitted provided that the following # + # conditions are met: # + # # + # * Redistributions of source code must retain the above copyright # + # notice, this list of conditions and the following disclaimer. # + # # + # * Redistributions in binary form must reproduce the above copyright # + # notice, this list of conditions and the following disclaimer in # + # the documentation and/or other materials provided with the # + # distribution. # + # # + # * Neither the name of UChicago Argonne, LLC, Argonne National # + # Laboratory, ANL, the U.S. Government, nor the names of its # + # contributors may be used to endorse or promote products derived # + # from this software without specific prior written permission. # + # # + # THIS SOFTWARE IS PROVIDED BY UChicago Argonne, LLC AND CONTRIBUTORS # + # "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT # + # LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS # + # FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL UChicago # + # Argonne, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, # + # INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, # + # BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; # + # LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER # + # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT # + # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN # + # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # + # POSSIBILITY OF SUCH DAMAGE. # + # ######################################################################### + + @staticmethod + def as_ndarray(arr, dtype=None, copy=False): + if not isinstance(arr, numpy.ndarray): + arr = numpy.array(arr, dtype=dtype, copy=copy) + return arr + + @staticmethod + def as_dtype(arr, dtype, copy=False): + if not arr.dtype == dtype: + arr = numpy.array(arr, dtype=dtype, copy=copy) + return arr + + @staticmethod + def as_float32(arr): + arr = CenterOfRotationFinder.as_ndarray(arr, numpy.float32) + return CenterOfRotationFinder.as_dtype(arr, numpy.float32) + + + + + @staticmethod + def find_center_vo(tomo, ind=None, smin=-40, smax=40, srad=10, step=0.5, + ratio=2., drop=20): + """ + Find rotation axis location using Nghia Vo's method. :cite:`Vo:14`. + + Parameters + ---------- + tomo : ndarray + 3D tomographic data. + ind : int, optional + Index of the slice to be used for reconstruction. + smin, smax : int, optional + Reference to the horizontal center of the sinogram. + srad : float, optional + Fine search radius. + step : float, optional + Step of fine searching. + ratio : float, optional + The ratio between the FOV of the camera and the size of object. + It's used to generate the mask. + drop : int, optional + Drop lines around vertical center of the mask. + + Returns + ------- + float + Rotation axis location. + + Notes + ----- + The function may not yield a correct estimate, if: + + - the sample size is bigger than the field of view of the camera. + In this case the ``ratio`` argument need to be set larger + than the default of 2.0. + + - there is distortion in the imaging hardware. If there's + no correction applied, the center of the projection image may + yield a better estimate. + + - the sample contrast is weak. Paganin's filter need to be applied + to overcome this. + + - the sample was changed during the scan. + """ + tomo = CenterOfRotationFinder.as_float32(tomo) + + if ind is None: + ind = tomo.shape[1] // 2 + _tomo = tomo[:, ind, :] + + + + # Reduce noise by smooth filters. Use different filters for coarse and fine search + _tomo_cs = ndimage.filters.gaussian_filter(_tomo, (3, 1)) + _tomo_fs = ndimage.filters.median_filter(_tomo, (2, 2)) + + # Coarse and fine searches for finding the rotation center. + if _tomo.shape[0] * _tomo.shape[1] > 4e6: # If data is large (>2kx2k) + #_tomo_coarse = downsample(numpy.expand_dims(_tomo_cs,1), level=2)[:, 0, :] + #init_cen = _search_coarse(_tomo_coarse, smin, smax, ratio, drop) + #fine_cen = _search_fine(_tomo_fs, srad, step, init_cen*4, ratio, drop) + init_cen = CenterOfRotationFinder._search_coarse(_tomo_cs, smin, + smax, ratio, drop) + fine_cen = CenterOfRotationFinder._search_fine(_tomo_fs, srad, + step, init_cen, + ratio, drop) + else: + init_cen = CenterOfRotationFinder._search_coarse(_tomo_cs, + smin, smax, + ratio, drop) + fine_cen = CenterOfRotationFinder._search_fine(_tomo_fs, srad, + step, init_cen, + ratio, drop) + + #logger.debug('Rotation center search finished: %i', fine_cen) + return fine_cen + + + @staticmethod + def _search_coarse(sino, smin, smax, ratio, drop): + """ + Coarse search for finding the rotation center. + """ + (Nrow, Ncol) = sino.shape + centerfliplr = (Ncol - 1.0) / 2.0 + + # Copy the sinogram and flip left right, the purpose is to + # make a full [0;2Pi] sinogram + _copy_sino = numpy.fliplr(sino[1:]) + + # This image is used for compensating the shift of sinogram 2 + temp_img = numpy.zeros((Nrow - 1, Ncol), dtype='float32') + temp_img[:] = sino[-1] + + # Start coarse search in which the shift step is 1 + listshift = numpy.arange(smin, smax + 1) + listmetric = numpy.zeros(len(listshift), dtype='float32') + mask = CenterOfRotationFinder._create_mask(2 * Nrow - 1, Ncol, + 0.5 * ratio * Ncol, drop) + for i in listshift: + _sino = numpy.roll(_copy_sino, i, axis=1) + if i >= 0: + _sino[:, 0:i] = temp_img[:, 0:i] + else: + _sino[:, i:] = temp_img[:, i:] + listmetric[i - smin] = numpy.sum(numpy.abs(numpy.fft.fftshift( + #pyfftw.interfaces.numpy_fft.fft2( + # numpy.vstack((sino, _sino))) + numpy.fft.fft2(numpy.vstack((sino, _sino))) + )) * mask) + minpos = numpy.argmin(listmetric) + return centerfliplr + listshift[minpos] / 2.0 + + @staticmethod + def _search_fine(sino, srad, step, init_cen, ratio, drop): + """ + Fine search for finding the rotation center. + """ + Nrow, Ncol = sino.shape + centerfliplr = (Ncol + 1.0) / 2.0 - 1.0 + # Use to shift the sinogram 2 to the raw CoR. + shiftsino = numpy.int16(2 * (init_cen - centerfliplr)) + _copy_sino = numpy.roll(numpy.fliplr(sino[1:]), shiftsino, axis=1) + if init_cen <= centerfliplr: + lefttake = numpy.int16(numpy.ceil(srad + 1)) + righttake = numpy.int16(numpy.floor(2 * init_cen - srad - 1)) + else: + lefttake = numpy.int16(numpy.ceil( + init_cen - (Ncol - 1 - init_cen) + srad + 1)) + righttake = numpy.int16(numpy.floor(Ncol - 1 - srad - 1)) + Ncol1 = righttake - lefttake + 1 + mask = CenterOfRotationFinder._create_mask(2 * Nrow - 1, Ncol1, + 0.5 * ratio * Ncol, drop) + numshift = numpy.int16((2 * srad) / step) + 1 + listshift = numpy.linspace(-srad, srad, num=numshift) + listmetric = numpy.zeros(len(listshift), dtype='float32') + factor1 = numpy.mean(sino[-1, lefttake:righttake]) + num1 = 0 + for i in listshift: + _sino = ndimage.interpolation.shift( + _copy_sino, (0, i), prefilter=False) + factor2 = numpy.mean(_sino[0,lefttake:righttake]) + _sino = _sino * factor1 / factor2 + sinojoin = numpy.vstack((sino, _sino)) + listmetric[num1] = numpy.sum(numpy.abs(numpy.fft.fftshift( + #pyfftw.interfaces.numpy_fft.fft2( + # sinojoin[:, lefttake:righttake + 1]) + numpy.fft.fft2(sinojoin[:, lefttake:righttake + 1]) + )) * mask) + num1 = num1 + 1 + minpos = numpy.argmin(listmetric) + return init_cen + listshift[minpos] / 2.0 + + @staticmethod + def _create_mask(nrow, ncol, radius, drop): + du = 1.0 / ncol + dv = (nrow - 1.0) / (nrow * 2.0 * numpy.pi) + centerrow = numpy.ceil(nrow / 2) - 1 + centercol = numpy.ceil(ncol / 2) - 1 + # added by Edoardo Pasca + centerrow = int(centerrow) + centercol = int(centercol) + mask = numpy.zeros((nrow, ncol), dtype='float32') + for i in range(nrow): + num1 = numpy.round(((i - centerrow) * dv / radius) / du) + (p1, p2) = numpy.int16(numpy.clip(numpy.sort( + (-num1 + centercol, num1 + centercol)), 0, ncol - 1)) + mask[i, p1:p2 + 1] = numpy.ones(p2 - p1 + 1, dtype='float32') + if drop < centerrow: + mask[centerrow - drop:centerrow + drop + 1, + :] = numpy.zeros((2 * drop + 1, ncol), dtype='float32') + mask[:,centercol-1:centercol+2] = numpy.zeros((nrow, 3), dtype='float32') + return mask + + def process(self, out=None): + + projections = self.get_input() + + cor = CenterOfRotationFinder.find_center_vo(projections.as_array()) + + return cor + + +class AcquisitionDataPadder(DataProcessor): + '''Normalization based on flat and dark + + This processor read in a AcquisitionData and normalises it based on + the instrument reading with and without incident photons or neutrons. + + Input: AcquisitionData + Parameter: 2D projection with flat field (or stack) + 2D projection with dark field (or stack) + Output: AcquisitionDataSetn + ''' + + def __init__(self, + center_of_rotation = None, + acquisition_geometry = None, + pad_value = 1e-5): + kwargs = { + 'acquisition_geometry' : acquisition_geometry, + 'center_of_rotation' : center_of_rotation, + 'pad_value' : pad_value + } + + super(AcquisitionDataPadder, self).__init__(**kwargs) + + def check_input(self, dataset): + if self.acquisition_geometry is None: + self.acquisition_geometry = dataset.geometry + if dataset.number_of_dimensions == 3: + return True + else: + raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ + .format(dataset.number_of_dimensions)) + + def process(self, out=None): + projections = self.get_input() + w = projections.get_dimension_size('horizontal') + delta = w - 2 * self.center_of_rotation + + padded_width = int ( + numpy.ceil(abs(delta)) + w + ) + delta_pix = padded_width - w + + voxel_per_pixel = 1 + geom = pbalg.pb_setup_geometry_from_acquisition(projections.as_array(), + self.acquisition_geometry.angles, + self.center_of_rotation, + voxel_per_pixel ) + + padded_geometry = self.acquisition_geometry.clone() + + padded_geometry.pixel_num_h = geom['n_h'] + padded_geometry.pixel_num_v = geom['n_v'] + + delta_pix_h = padded_geometry.pixel_num_h - self.acquisition_geometry.pixel_num_h + delta_pix_v = padded_geometry.pixel_num_v - self.acquisition_geometry.pixel_num_v + + if delta_pix_h == 0: + delta_pix_h = delta_pix + padded_geometry.pixel_num_h = padded_width + #initialize a new AcquisitionData with values close to 0 + out = AcquisitionData(geometry=padded_geometry) + out = out + self.pad_value + + + #pad in the horizontal-vertical plane -> slice on angles + if delta > 0: + #pad left of middle + command = "out.array[" + for i in range(out.number_of_dimensions): + if out.dimension_labels[i] == 'horizontal': + value = '{0}:{1}'.format(delta_pix_h, delta_pix_h+w) + command = command + str(value) + else: + if out.dimension_labels[i] == 'vertical' : + value = '{0}:'.format(delta_pix_v) + command = command + str(value) + else: + command = command + ":" + if i < out.number_of_dimensions -1: + command = command + ',' + command = command + '] = projections.array' + #print (command) + else: + #pad right of middle + command = "out.array[" + for i in range(out.number_of_dimensions): + if out.dimension_labels[i] == 'horizontal': + value = '{0}:{1}'.format(0, w) + command = command + str(value) + else: + if out.dimension_labels[i] == 'vertical' : + value = '{0}:'.format(delta_pix_v) + command = command + str(value) + else: + command = command + ":" + if i < out.number_of_dimensions -1: + command = command + ',' + command = command + '] = projections.array' + #print (command) + #cleaned = eval(command) + exec(command) return out
\ No newline at end of file diff --git a/Wrappers/Python/conda-recipe/conda_build_config.yaml b/Wrappers/Python/conda-recipe/conda_build_config.yaml index 96a211f..30c8e9d 100644 --- a/Wrappers/Python/conda-recipe/conda_build_config.yaml +++ b/Wrappers/Python/conda-recipe/conda_build_config.yaml @@ -4,5 +4,5 @@ python: - 3.6 numpy: # TODO investigage, as it doesn't currently build with cvxp, requires >1.14 - #- 1.12 + - 1.12 - 1.15 diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml index 8ded429..dd3238e 100644 --- a/Wrappers/Python/conda-recipe/meta.yaml +++ b/Wrappers/Python/conda-recipe/meta.yaml @@ -11,7 +11,7 @@ build: test: requires: - python-wget - - cvxpy # [not win] + - cvxpy # [ unix and py36 and np115 ] source_files: - ./test # [win] @@ -24,8 +24,9 @@ test: requirements: build: + - {{ pin_compatible('numpy', max_pin='x.x') }} - python - - numpy {{ numpy }} + - numpy - setuptools run: @@ -33,7 +34,7 @@ requirements: - python - numpy - scipy - - matplotlib + #- matplotlib - h5py about: diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py index eaf124b..87930b5 100644 --- a/Wrappers/Python/setup.py +++ b/Wrappers/Python/setup.py @@ -31,8 +31,11 @@ if cil_version == '': setup( name="ccpi-framework", version=cil_version, - packages=['ccpi' , 'ccpi.io', 'ccpi.optimisation', - 'ccpi.optimisation.algorithms'], + packages=['ccpi' , 'ccpi.io', + 'ccpi.framework', 'ccpi.optimisation', + 'ccpi.optimisation.operators', + 'ccpi.optimisation.algorithms', + 'ccpi.optimisation.functions'], # Project uses reStructuredText, so ensure that the docutils get # installed or upgraded on the target machine diff --git a/Wrappers/Python/test/test_BlockDataContainer.py b/Wrappers/Python/test/test_BlockDataContainer.py new file mode 100755 index 0000000..6c0bede --- /dev/null +++ b/Wrappers/Python/test/test_BlockDataContainer.py @@ -0,0 +1,367 @@ +# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 5 16:08:23 2019
+
+@author: ofn77899
+"""
+
+import unittest
+import numpy
+#from ccpi.plugins.ops import CCPiProjectorSimple
+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.framework import ImageData, AcquisitionData
+#from ccpi.optimisation.algorithms import GradientDescent
+from ccpi.framework import BlockDataContainer
+#from ccpi.optimisation.Algorithms import CGLS
+import functools
+
+from ccpi.optimisation.operators import Gradient, Identity, BlockOperator
+
+class TestBlockDataContainer(unittest.TestCase):
+ def skiptest_BlockDataContainerShape(self):
+ print ("test block data container")
+ ig0 = ImageGeometry(12,42,55,32)
+ 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 = BlockDataContainer(data0,data1)
+ cp1 = BlockDataContainer(data2,data3)
+ transpose_shape = (cp0.shape[1], cp0.shape[0])
+ self.assertTrue(cp0.T.shape == transpose_shape)
+ def skiptest_BlockDataContainerShapeArithmetic(self):
+ print ("test block data container")
+ ig0 = ImageGeometry(2,3,4)
+ ig1 = ImageGeometry(2,3,4)
+
+ data0 = ImageData(geometry=ig0)
+ data1 = ImageData(geometry=ig1) + 1
+
+ data2 = ImageData(geometry=ig0) + 2
+ data3 = ImageData(geometry=ig1) + 3
+
+ cp0 = BlockDataContainer(data0,data1)
+ #cp1 = BlockDataContainer(data2,data3)
+ cp1 = cp0 + 1
+ self.assertTrue(cp1.shape == cp0.shape)
+ cp1 = cp0.T + 1
+
+ transpose_shape = (cp0.shape[1], cp0.shape[0])
+ self.assertTrue(cp1.shape == transpose_shape)
+
+ cp1 = cp0.T - 1
+ transpose_shape = (cp0.shape[1], cp0.shape[0])
+ self.assertTrue(cp1.shape == transpose_shape)
+
+ cp1 = (cp0.T + 1)*2
+ transpose_shape = (cp0.shape[1], cp0.shape[0])
+ self.assertTrue(cp1.shape == transpose_shape)
+
+ cp1 = (cp0.T + 1)/2
+ transpose_shape = (cp0.shape[1], cp0.shape[0])
+ self.assertTrue(cp1.shape == transpose_shape)
+
+ cp1 = cp0.T.power(2.2)
+ transpose_shape = (cp0.shape[1], cp0.shape[0])
+ self.assertTrue(cp1.shape == transpose_shape)
+
+ cp1 = cp0.T.maximum(3)
+ transpose_shape = (cp0.shape[1], cp0.shape[0])
+ self.assertTrue(cp1.shape == transpose_shape)
+
+ cp1 = cp0.T.abs()
+ transpose_shape = (cp0.shape[1], cp0.shape[0])
+ self.assertTrue(cp1.shape == transpose_shape)
+
+ cp1 = cp0.T.sign()
+ transpose_shape = (cp0.shape[1], cp0.shape[0])
+ self.assertTrue(cp1.shape == transpose_shape)
+
+ cp1 = cp0.T.sqrt()
+ transpose_shape = (cp0.shape[1], cp0.shape[0])
+ self.assertTrue(cp1.shape == transpose_shape)
+
+ cp1 = cp0.T.conjugate()
+ transpose_shape = (cp0.shape[1], cp0.shape[0])
+ self.assertTrue(cp1.shape == transpose_shape)
+
+ def test_BlockDataContainer(self):
+ print ("test block data container")
+ ig0 = ImageGeometry(2,3,4)
+ ig1 = ImageGeometry(2,3,4)
+
+ data0 = ImageData(geometry=ig0)
+ data1 = ImageData(geometry=ig1) + 1
+
+ data2 = ImageData(geometry=ig0) + 2
+ data3 = ImageData(geometry=ig1) + 3
+
+ cp0 = BlockDataContainer(data0,data1)
+ cp1 = BlockDataContainer(data2,data3)
+ #
+ a = [ (el, ot) for el,ot in zip(cp0.containers,cp1.containers)]
+ print (a[0][0].shape)
+ #cp2 = BlockDataContainer(*a)
+ cp2 = cp0.add(cp1)
+ assert (cp2.get_item(0).as_array()[0][0][0] == 2.)
+ assert (cp2.get_item(1).as_array()[0][0][0] == 4.)
+
+ cp2 = cp0 + cp1
+ assert (cp2.get_item(0).as_array()[0][0][0] == 2.)
+ assert (cp2.get_item(1).as_array()[0][0][0] == 4.)
+ cp2 = cp0 + 1
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 1. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 2., decimal = 5)
+ cp2 = cp0 + [1 ,2]
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 1. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 3., decimal = 5)
+ cp2 += cp1
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , +3. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , +6., decimal = 5)
+
+ cp2 += 1
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , +4. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , +7., decimal = 5)
+
+ cp2 += [-2,-1]
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 2. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 6., decimal = 5)
+
+
+ cp2 = cp0.subtract(cp1)
+ assert (cp2.get_item(0).as_array()[0][0][0] == -2.)
+ assert (cp2.get_item(1).as_array()[0][0][0] == -2.)
+ cp2 = cp0 - cp1
+ assert (cp2.get_item(0).as_array()[0][0][0] == -2.)
+ assert (cp2.get_item(1).as_array()[0][0][0] == -2.)
+
+ cp2 = cp0 - 1
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , -1. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 0, decimal = 5)
+ cp2 = cp0 - [1 ,2]
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , -1. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , -1., decimal = 5)
+
+ cp2 -= cp1
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , -3. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , -4., decimal = 5)
+
+ cp2 -= 1
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , -4. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , -5., decimal = 5)
+
+ cp2 -= [-2,-1]
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , -2. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , -4., decimal = 5)
+
+
+ cp2 = cp0.multiply(cp1)
+ assert (cp2.get_item(0).as_array()[0][0][0] == 0.)
+ assert (cp2.get_item(1).as_array()[0][0][0] == 3.)
+ cp2 = cp0 * cp1
+ assert (cp2.get_item(0).as_array()[0][0][0] == 0.)
+ assert (cp2.get_item(1).as_array()[0][0][0] == 3.)
+
+ cp2 = cp0 * 2
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 0. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 2, decimal = 5)
+ cp2 = 2 * cp0
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 0. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 2, decimal = 5)
+ cp2 = cp0 * [3 ,2]
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 0. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 2., decimal = 5)
+ cp2 = cp0 * numpy.asarray([3 ,2])
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 0. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 2., decimal = 5)
+
+ cp2 = [3,2] * cp0
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 0. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 2., decimal = 5)
+ cp2 = numpy.asarray([3,2]) * cp0
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 0. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 2., decimal = 5)
+ cp2 = [3,2,3] * cp0
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 0. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 2., decimal = 5)
+
+ cp2 *= cp1
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 0 , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , +6., decimal = 5)
+
+ cp2 *= 1
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 0. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , +6., decimal = 5)
+
+ cp2 *= [-2,-1]
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 0. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , -6., decimal = 5)
+
+
+ cp2 = cp0.divide(cp1)
+ assert (cp2.get_item(0).as_array()[0][0][0] == 0.)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0], 1./3., decimal=4)
+ cp2 = cp0/cp1
+ assert (cp2.get_item(0).as_array()[0][0][0] == 0.)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0], 1./3., decimal=4)
+
+ cp2 = cp0 / 2
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 0. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 0.5, decimal = 5)
+ cp2 = cp0 / [3 ,2]
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 0. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 0.5, decimal = 5)
+ cp2 = cp0 / numpy.asarray([3 ,2])
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 0. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 0.5, decimal = 5)
+ cp3 = numpy.asarray([3 ,2]) / (cp0+1)
+ numpy.testing.assert_almost_equal(cp3.get_item(0).as_array()[0][0][0] , 3. , decimal=5)
+ numpy.testing.assert_almost_equal(cp3.get_item(1).as_array()[0][0][0] , 1, decimal = 5)
+
+ cp2 += 1
+ cp2 /= cp1
+ # TODO fix inplace division
+
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 1./2 , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 1.5/3., decimal = 5)
+
+ cp2 /= 1
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 0.5 , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 0.5, decimal = 5)
+
+ cp2 /= [-2,-1]
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , -0.5/2. , decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , -0.5, decimal = 5)
+ ####
+
+ cp2 = cp0.power(cp1)
+ assert (cp2.get_item(0).as_array()[0][0][0] == 0.)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0], 1., decimal=4)
+ cp2 = cp0**cp1
+ assert (cp2.get_item(0).as_array()[0][0][0] == 0.)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0], 1., decimal=4)
+
+ cp2 = cp0 ** 2
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 0., decimal=5)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 1., decimal = 5)
+
+ cp2 = cp0.maximum(cp1)
+ assert (cp2.get_item(0).as_array()[0][0][0] == cp1.get_item(0).as_array()[0][0][0])
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0], cp2.get_item(1).as_array()[0][0][0], decimal=4)
+
+
+ cp2 = cp0.abs()
+ numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0], 0., decimal=4)
+ numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0], 1., decimal=4)
+
+ cp2 = cp0.subtract(cp1)
+ s = cp2.sign()
+ numpy.testing.assert_almost_equal(s.get_item(0).as_array()[0][0][0], -1., decimal=4)
+ numpy.testing.assert_almost_equal(s.get_item(1).as_array()[0][0][0], -1., decimal=4)
+
+ cp2 = cp0.add(cp1)
+ s = cp2.sqrt()
+ numpy.testing.assert_almost_equal(s.get_item(0).as_array()[0][0][0], numpy.sqrt(2), decimal=4)
+ numpy.testing.assert_almost_equal(s.get_item(1).as_array()[0][0][0], numpy.sqrt(4), decimal=4)
+
+ s = cp0.sum()
+ size = functools.reduce(lambda x,y: x*y, data1.shape, 1)
+ print ("size" , size)
+ numpy.testing.assert_almost_equal(s, 0 + size, decimal=4)
+ s0 = 1
+ s1 = 1
+ for i in cp0.get_item(0).shape:
+ s0 *= i
+ for i in cp0.get_item(1).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)
+ def test_Nested_BlockDataContainer(self):
+ print ("test_Nested_BlockDataContainer")
+ ig0 = ImageGeometry(2,3,4)
+ ig1 = ImageGeometry(2,3,4)
+
+ data0 = ImageData(geometry=ig0)
+ data1 = ImageData(geometry=ig1) + 1
+
+ data2 = ImageData(geometry=ig0) + 2
+ data3 = ImageData(geometry=ig1) + 3
+
+ cp0 = BlockDataContainer(data0,data1)
+ cp1 = BlockDataContainer(data2,data3)
+
+ nbdc = BlockDataContainer(cp0, cp1)
+ nbdc2 = nbdc + 2
+ numpy.testing.assert_almost_equal(nbdc2.get_item(0).get_item(0).as_array()[0][0][0] , 2. , decimal=5)
+ numpy.testing.assert_almost_equal(nbdc2.get_item(0).get_item(1).as_array()[0][0][0] , 3. , decimal=5)
+ numpy.testing.assert_almost_equal(nbdc2.get_item(1).get_item(0).as_array()[0][0][0] , 4. , decimal=5)
+ numpy.testing.assert_almost_equal(nbdc2.get_item(1).get_item(1).as_array()[0][0][0] , 5. , decimal=5)
+
+ nbdc2 = 2 + nbdc
+ numpy.testing.assert_almost_equal(nbdc2.get_item(0).get_item(0).as_array()[0][0][0] , 2. , decimal=5)
+ numpy.testing.assert_almost_equal(nbdc2.get_item(0).get_item(1).as_array()[0][0][0] , 3. , decimal=5)
+ numpy.testing.assert_almost_equal(nbdc2.get_item(1).get_item(0).as_array()[0][0][0] , 4. , decimal=5)
+ numpy.testing.assert_almost_equal(nbdc2.get_item(1).get_item(1).as_array()[0][0][0] , 5. , decimal=5)
+
+
+ nbdc2 = nbdc * 2
+ numpy.testing.assert_almost_equal(nbdc2.get_item(0).get_item(0).as_array()[0][0][0] , 0. , decimal=5)
+ numpy.testing.assert_almost_equal(nbdc2.get_item(0).get_item(1).as_array()[0][0][0] , 2. , decimal=5)
+ numpy.testing.assert_almost_equal(nbdc2.get_item(1).get_item(0).as_array()[0][0][0] , 4. , decimal=5)
+ numpy.testing.assert_almost_equal(nbdc2.get_item(1).get_item(1).as_array()[0][0][0] , 6. , decimal=5)
+
+ nbdc2 = 2 * nbdc
+ numpy.testing.assert_almost_equal(nbdc2.get_item(0).get_item(0).as_array()[0][0][0] , 0. , decimal=5)
+ numpy.testing.assert_almost_equal(nbdc2.get_item(0).get_item(1).as_array()[0][0][0] , 2. , decimal=5)
+ numpy.testing.assert_almost_equal(nbdc2.get_item(1).get_item(0).as_array()[0][0][0] , 4. , decimal=5)
+ numpy.testing.assert_almost_equal(nbdc2.get_item(1).get_item(1).as_array()[0][0][0] , 6. , decimal=5)
+
+ nbdc2 = nbdc / 2
+ numpy.testing.assert_almost_equal(nbdc2.get_item(0).get_item(0).as_array()[0][0][0] , 0. , decimal=5)
+ numpy.testing.assert_almost_equal(nbdc2.get_item(0).get_item(1).as_array()[0][0][0] , .5 , decimal=5)
+ numpy.testing.assert_almost_equal(nbdc2.get_item(1).get_item(0).as_array()[0][0][0] , 1. , decimal=5)
+ numpy.testing.assert_almost_equal(nbdc2.get_item(1).get_item(1).as_array()[0][0][0] , 3./2 , decimal=5)
+
+ c5 = nbdc.get_item(0).power(2).sum()
+ c5a = nbdc.power(2).sum()
+ print ("sum", c5a, c5)
+
+ print ("test_Nested_BlockDataContainer OK")
+ def stest_NestedBlockDataContainer2(self):
+ M, N = 2, 3
+ ig = ImageGeometry(voxel_num_x = M, voxel_num_y = N)
+ ag = ig
+ u = ig.allocate(1)
+ op1 = Gradient(ig)
+ op2 = Identity(ig, ag)
+
+ operator = BlockOperator(op1, op2, shape=(2,1))
+
+ d1 = op1.direct(u)
+ d2 = op2.direct(u)
+
+ d = operator.direct(u)
+
+ dd = operator.domain_geometry()
+ ww = operator.range_geometry()
+
+ print(d.get_item(0).get_item(0).as_array())
+ print(d.get_item(0).get_item(1).as_array())
+ print(d.get_item(1).as_array())
+
+ c1 = d + d
+
+ c2 = 2*d
+
+ c3 = d / (d+0.0001)
+
+
+ c5 = d.get_item(0).power(2).sum()
+
diff --git a/Wrappers/Python/test/test_BlockOperator.py b/Wrappers/Python/test/test_BlockOperator.py new file mode 100644 index 0000000..951aa0a --- /dev/null +++ b/Wrappers/Python/test/test_BlockOperator.py @@ -0,0 +1,313 @@ +import unittest +from ccpi.optimisation.operators import BlockOperator +from ccpi.framework import BlockDataContainer +from ccpi.optimisation.ops import TomoIdentity +from ccpi.framework import ImageGeometry, ImageData +import numpy +from ccpi.optimisation.operators import FiniteDiff + +class TestBlockOperator(unittest.TestCase): + + def test_BlockOperator(self): + print ("test_BlockOperator") + ig = [ ImageGeometry(10,20,30) , \ + ImageGeometry(10,20,30) , \ + ImageGeometry(10,20,30) ] + x = [ g.allocate() for g in ig ] + ops = [ TomoIdentity(g) for g in ig ] + + K = BlockOperator(*ops) + X = BlockDataContainer(x[0]) + Y = K.direct(X) + self.assertTrue(Y.shape == K.shape) + + numpy.testing.assert_array_equal(Y.get_item(0).as_array(),X.get_item(0).as_array()) + numpy.testing.assert_array_equal(Y.get_item(1).as_array(),X.get_item(0).as_array()) + #numpy.testing.assert_array_equal(Y.get_item(2).as_array(),X.get_item(2).as_array()) + + X = BlockDataContainer(*x) + 1 + Y = K.T.direct(X) + # K.T (1,3) X (3,1) => output shape (1,1) + self.assertTrue(Y.shape == (1,1)) + zero = numpy.zeros(X.get_item(0).shape) + numpy.testing.assert_array_equal(Y.get_item(0).as_array(),len(x)+zero) + + + def test_ScaledBlockOperatorSingleScalar(self): + ig = [ ImageGeometry(10,20,30) , \ + ImageGeometry(10,20,30) , \ + ImageGeometry(10,20,30) ] + x = [ g.allocate() for g in ig ] + ops = [ TomoIdentity(g) for g in ig ] + + val = 1 + # test limit as non Scaled + scalar = 1 + k = BlockOperator(*ops) + K = scalar * k + X = BlockDataContainer(*x) + val + + Y = K.T.direct(X) + self.assertTrue(Y.shape == (1,1)) + zero = numpy.zeros(X.get_item(0).shape) + xx = numpy.asarray([val for _ in x]) + numpy.testing.assert_array_equal(Y.get_item(0).as_array(),((scalar*xx).sum()+zero)) + + scalar = 0.5 + k = BlockOperator(*ops) + K = scalar * k + X = BlockDataContainer(*x) + 1 + + Y = K.T.direct(X) + self.assertTrue(Y.shape == (1,1)) + zero = numpy.zeros(X.get_item(0).shape) + numpy.testing.assert_array_equal(Y.get_item(0).as_array(),scalar*(len(x)+zero)) + + + def test_ScaledBlockOperatorScalarList(self): + ig = [ ImageGeometry(2,3) , \ + #ImageGeometry(10,20,30) , \ + ImageGeometry(2,3 ) ] + x = [ g.allocate() for g in ig ] + ops = [ TomoIdentity(g) for g in ig ] + + + # test limit as non Scaled + scalar = numpy.asarray([1 for _ in x]) + k = BlockOperator(*ops) + K = scalar * k + val = 1 + X = BlockDataContainer(*x) + val + + Y = K.T.direct(X) + self.assertTrue(Y.shape == (1,1)) + zero = numpy.zeros(X.get_item(0).shape) + xx = numpy.asarray([val for _ in x]) + numpy.testing.assert_array_equal(Y.get_item(0).as_array(),(scalar*xx).sum()+zero) + + scalar = numpy.asarray([i+1 for i,el in enumerate(x)]) + #scalar = numpy.asarray([6,0]) + k = BlockOperator(*ops) + K = scalar * k + X = BlockDataContainer(*x) + val + Y = K.T.direct(X) + self.assertTrue(Y.shape == (1,1)) + zero = numpy.zeros(X.get_item(0).shape) + xx = numpy.asarray([val for _ in x]) + + + numpy.testing.assert_array_equal(Y.get_item(0).as_array(), + (scalar*xx).sum()+zero) + + + def test_TomoIdentity(self): + ig = ImageGeometry(10,20,30) + img = ig.allocate() + print (img.shape, ig.shape) + self.assertTrue(img.shape == (30,20,10)) + self.assertEqual(img.sum(), 0) + Id = TomoIdentity(ig) + y = Id.direct(img) + numpy.testing.assert_array_equal(y.as_array(), img.as_array()) + + def skiptest_CGLS_tikhonov(self): + from ccpi.optimisation.algorithms import CGLS + + from ccpi.plugins.ops import CCPiProjectorSimple + 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.Algorithms import GradientDescent + #from ccpi.optimisation.Algorithms import CGLS + import matplotlib.pyplot as plt + + + # Set up phantom size N x N x vert by creating ImageGeometry, initialising the + # ImageData object with this geometry and empty array and finally put some + # data into its array, and display one slice as image. + + # Image parameters + N = 128 + vert = 4 + + # Set up image geometry + ig = ImageGeometry(voxel_num_x=N, + voxel_num_y=N, + voxel_num_z=vert) + + # Set up empty image data + Phantom = ImageData(geometry=ig, + dimension_labels=['horizontal_x', + 'horizontal_y', + 'vertical']) + Phantom += 0.05 + # Populate image data by looping over and filling slices + i = 0 + while i < vert: + if vert > 1: + x = Phantom.subset(vertical=i).array + else: + x = Phantom.array + x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 + x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 0.94 + if vert > 1 : + Phantom.fill(x, vertical=i) + i += 1 + + + perc = 0.02 + # Set up empty image data + noise = ImageData(numpy.random.normal(loc = 0.04 , + scale = perc , + size = Phantom.shape), geometry=ig, + dimension_labels=['horizontal_x', + 'horizontal_y', + 'vertical']) + Phantom += noise + + # Set up AcquisitionGeometry object to hold the parameters of the measurement + # setup geometry: # Number of angles, the actual angles from 0 to + # pi for parallel beam, set the width of a detector + # pixel relative to an object pixe and the number of detector pixels. + angles_num = 20 + det_w = 1.0 + det_num = N + + angles = numpy.linspace(0,numpy.pi,angles_num,endpoint=False,dtype=numpy.float32)*\ + 180/numpy.pi + + # Inputs: Geometry, 2D or 3D, angles, horz detector pixel count, + # horz detector pixel size, vert detector pixel count, + # vert detector pixel size. + ag = AcquisitionGeometry('parallel', + '3D', + angles, + N, + det_w, + vert, + det_w) + + # Set up Operator object combining the ImageGeometry and AcquisitionGeometry + # wrapping calls to CCPi projector. + A = CCPiProjectorSimple(ig, ag) + + # Forward and backprojection are available as methods direct and adjoint. Here + # generate test data b and some noise + + b = A.direct(Phantom) + + + #z = A.adjoint(b) + + + # Using the test data b, different reconstruction methods can now be set up as + # demonstrated in the rest of this file. In general all methods need an initial + # guess and some algorithm options to be set. Note that 100 iterations for + # some of the methods is a very low number and 1000 or 10000 iterations may be + # needed if one wants to obtain a converged solution. + x_init = ImageData(geometry=ig, + dimension_labels=['horizontal_x','horizontal_y','vertical']) + X_init = BlockDataContainer(x_init) + B = BlockDataContainer(b, + ImageData(geometry=ig, dimension_labels=['horizontal_x','horizontal_y','vertical'])) + + # setup a tomo identity + Ibig = 1e5 * TomoIdentity(geometry=ig) + Ismall = 1e-5 * TomoIdentity(geometry=ig) + + # composite operator + Kbig = BlockOperator(A, Ibig, shape=(2,1)) + Ksmall = BlockOperator(A, Ismall, shape=(2,1)) + + #out = K.direct(X_init) + + f = Norm2sq(Kbig,B) + f.L = 0.00003 + + fsmall = Norm2sq(Ksmall,B) + f.L = 0.00003 + + simplef = Norm2sq(A, b) + simplef.L = 0.00003 + + gd = GradientDescent( x_init=x_init, objective_function=simplef, + rate=simplef.L) + gd.max_iteration = 10 + + cg = CGLS() + cg.set_up(X_init, Kbig, B ) + cg.max_iteration = 1 + + cgsmall = CGLS() + cgsmall.set_up(X_init, Ksmall, B ) + cgsmall.max_iteration = 1 + + + cgs = CGLS() + cgs.set_up(x_init, A, b ) + cgs.max_iteration = 6 + # + #out.__isub__(B) + #out2 = K.adjoint(out) + + #(2.0*self.c)*self.A.adjoint( self.A.direct(x) - self.b ) + + #for _ in gd: + # print ("iteration {} {}".format(gd.iteration, gd.get_current_loss())) + + #cg.run(10, lambda it,val: print ("iteration {} objective {}".format(it,val)) ) + + #cgs.run(10, lambda it,val: print ("iteration {} objective {}".format(it,val))) + + #cgsmall.run(10, lambda it,val: print ("iteration {} objective {}".format(it,val))) + #cgsmall.run(10, lambda it,val: print ("iteration {} objective {}".format(it,val))) + # for _ in cg: + # print ("iteration {} {}".format(cg.iteration, cg.get_current_loss())) + # + # fig = plt.figure() + # plt.imshow(cg.get_output().get_item(0,0).subset(vertical=0).as_array()) + # plt.title('Composite CGLS') + # plt.show() + # + # for _ in cgs: + # print ("iteration {} {}".format(cgs.iteration, cgs.get_current_loss())) + # + fig = plt.figure() + plt.subplot(1,5,1) + plt.imshow(Phantom.subset(vertical=0).as_array()) + plt.title('Simulated Phantom') + plt.subplot(1,5,2) + plt.imshow(gd.get_output().subset(vertical=0).as_array()) + plt.title('Simple Gradient Descent') + plt.subplot(1,5,3) + plt.imshow(cgs.get_output().subset(vertical=0).as_array()) + plt.title('Simple CGLS') + plt.subplot(1,5,4) + plt.imshow(cg.get_output().get_item(0,0).subset(vertical=0).as_array()) + plt.title('Composite CGLS\nbig lambda') + plt.subplot(1,5,5) + plt.imshow(cgsmall.get_output().get_item(0,0).subset(vertical=0).as_array()) + plt.title('Composite CGLS\nsmall lambda') + plt.show() + + def test_FiniteDiffOperator(self): + N, M = 200, 300 + + + ig = ImageGeometry(voxel_num_x = M, voxel_num_y = N) + u = ig.allocate('random_int') + G = FiniteDiff(ig, direction=0, bnd_cond = 'Neumann') + print(type(u), u.as_array()) + print(G.direct(u).as_array()) + + # Gradient Operator norm, for one direction should be close to 2 + numpy.testing.assert_allclose(G.norm(), numpy.sqrt(4), atol=0.1) + + M1, N1, K1 = 200, 300, 2 + ig1 = ImageGeometry(voxel_num_x = M1, voxel_num_y = N1, channels = K1) + u1 = ig1.allocate('random_int') + G1 = FiniteDiff(ig1, direction=2, bnd_cond = 'Periodic') + print(ig1.shape==u1.shape) + print (G1.norm()) + numpy.testing.assert_allclose(G1.norm(), numpy.sqrt(4), atol=0.1)
\ No newline at end of file diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py index f23179c..7a7e6a0 100755 --- a/Wrappers/Python/test/test_DataContainer.py +++ b/Wrappers/Python/test/test_DataContainer.py @@ -174,7 +174,7 @@ class TestDataContainer(unittest.TestCase): def binary_add(self): print("Test binary add") X, Y, Z = 512, 512, 512 - X, Y, Z = 256, 512, 512 + X, Y, Z = 1024, 512, 512 steps = [timer()] a = numpy.ones((X, Y, Z), dtype='float32') steps.append(timer()) @@ -495,9 +495,10 @@ class TestDataContainer(unittest.TestCase): self.assertEqual(order[1], image.dimension_labels[1]) self.assertEqual(order[2], image.dimension_labels[2]) def test_AcquisitionGeometry_allocate(self): - ageometry = AcquisitionGeometry(dimension=2, angles=numpy.linspace(0, 180, num=10), - geom_type='parallel', pixel_num_v=3, - pixel_num_h=5, channels=2) + ageometry = AcquisitionGeometry(dimension=2, + angles=numpy.linspace(0, 180, num=10), + geom_type='parallel', pixel_num_v=3, + pixel_num_h=5, channels=2) sino = ageometry.allocate() shape = sino.shape print ("shape", shape) @@ -509,8 +510,8 @@ class TestDataContainer(unittest.TestCase): self.assertEqual(1,sino.as_array()[shape[0]-1][shape[1]-1][shape[2]-1][shape[3]-1]) print (sino.dimension_labels, sino.shape, ageometry) - default_order = ['channel' , ' angle' , - 'vertical' , 'horizontal'] + default_order = ['channel' , 'angle' , + 'vertical' , 'horizontal'] self.assertEqual(default_order[0], sino.dimension_labels[0]) self.assertEqual(default_order[1], sino.dimension_labels[1]) self.assertEqual(default_order[2], sino.dimension_labels[2]) @@ -561,4 +562,4 @@ class TestDataContainer(unittest.TestCase): if __name__ == '__main__': unittest.main() -
\ No newline at end of file + diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py new file mode 100644 index 0000000..46e8c7c --- /dev/null +++ b/Wrappers/Python/test/test_Operator.py @@ -0,0 +1,24 @@ +import unittest +#from ccpi.optimisation.operators import Operator +from ccpi.optimisation.ops import TomoIdentity +from ccpi.framework import ImageGeometry, ImageData +import numpy + +class TestOperator(unittest.TestCase): + def test_ScaledOperator(self): + ig = ImageGeometry(10,20,30) + img = ig.allocate() + scalar = 0.5 + sid = scalar * TomoIdentity(ig) + numpy.testing.assert_array_equal(scalar * img.as_array(), sid.direct(img).as_array()) + + + def test_TomoIdentity(self): + ig = ImageGeometry(10,20,30) + img = ig.allocate() + self.assertTrue(img.shape == (30,20,10)) + self.assertEqual(img.sum(), 0) + Id = TomoIdentity(ig) + y = Id.direct(img) + numpy.testing.assert_array_equal(y.as_array(), img.as_array()) + diff --git a/Wrappers/Python/test/test_functions.py b/Wrappers/Python/test/test_functions.py new file mode 100644 index 0000000..6a44641 --- /dev/null +++ b/Wrappers/Python/test/test_functions.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat Mar 2 19:24:37 2019 + +@author: evangelos +""" + + +import numpy as np +#from ccpi.optimisation.funcs import Function +from ccpi.optimisation.functions import Function +from ccpi.framework import DataContainer, ImageData, ImageGeometry +from ccpi.optimisation.operators import Identity +from ccpi.optimisation.operators import BlockOperator +from ccpi.framework import BlockDataContainer +from numbers import Number +from ccpi.optimisation.operators import Gradient + +from ccpi.optimisation.functions import SimpleL2NormSq +from ccpi.optimisation.functions import L2NormSq +from ccpi.optimisation.functions import SimpleL1Norm +from ccpi.optimisation.functions import L1Norm +# from ccpi.optimisation.functions.L2NormSquared import SimpleL2NormSq, L2NormSq +# from ccpi.optimisation.functions.L1Norm import SimpleL1Norm, L1Norm +from ccpi.optimisation.functions import mixed_L12Norm +from ccpi.optimisation.functions import ZeroFun + +from ccpi.optimisation.functions import FunctionOperatorComposition +import unittest + +# + + +class TestFunction(unittest.TestCase): + def test_Function(self): + + + N = 3 + ig = (N,N) + ag = ig + op1 = Gradient(ig) + op2 = Identity(ig, ag) + + # Form Composite Operator + operator = BlockOperator((2,1), op1, op2 ) + + # Create functions + noisy_data = ImageData(np.random.randint(10, size=ag)) + + d = ImageData(np.random.randint(10, size=ag)) + alpha = 0.5 + # scaled function + g = alpha * L2NormSq(b=noisy_data) + + # Compare call of g + a2 = alpha*(d - noisy_data).power(2).sum() + #print(a2, g(d)) + self.assertEqual(a2, g(d)) + + # Compare convex conjugate of g + a3 = 0.5 * d.power(2).sum() + (d*noisy_data).sum() + self.assertEqual(a3, g.convex_conjugate(d)) + #print( a3, g.convex_conjugate(d)) + + + def stest_ScaledFunctin(self): + ig = (N,N) + ag = ig + op1 = Gradient(ig) + op2 = Identity(ig, ag) + + +# +# f1 = L2NormSq(alpha=1, b=noisy_data) +# print(f1(noisy_data)) +# +# f2 = L2NormSq(alpha=5, b=noisy_data).composition_with(op2) +# print(f2(noisy_data)) +# +# print(f1.gradient(noisy_data).as_array()) +# print(f2.gradient(noisy_data).as_array()) +## +# print(f1.proximal(noisy_data,1).as_array()) +# print(f2.proximal(noisy_data,1).as_array()) +# +# +# f3 = mixed_L12Norm(alpha = 1).composition_with(op1) +# print(f3(noisy_data)) +# +# print(ImageData(op1.direct(noisy_data).power(2).sum(axis=0)).sqrt().sum()) +# +# print( 5*(op2.direct(d) - noisy_data).power(2).sum(), f2(d)) +# +# from functions import mixed_L12Norm as mixed_L12Norm_old +# +# print(mixed_L12Norm_old(op1,None,alpha)(noisy_data)) + + + # + + diff --git a/Wrappers/Python/wip/CGLS_tikhonov.py b/Wrappers/Python/wip/CGLS_tikhonov.py new file mode 100644 index 0000000..e9bbcd9 --- /dev/null +++ b/Wrappers/Python/wip/CGLS_tikhonov.py @@ -0,0 +1,196 @@ +from ccpi.optimisation.algorithms import CGLS + +from ccpi.plugins.ops import CCPiProjectorSimple +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, ImageData, AcquisitionData +from ccpi.optimisation.algorithms import GradientDescent +#from ccpi.optimisation.algorithms import CGLS +import matplotlib.pyplot as plt +import numpy +from ccpi.framework import BlockDataContainer +from ccpi.optimisation.operators import BlockOperator + +# Set up phantom size N x N x vert by creating ImageGeometry, initialising the +# ImageData object with this geometry and empty array and finally put some +# data into its array, and display one slice as image. + +# Image parameters +N = 128 +vert = 4 + +# Set up image geometry +ig = ImageGeometry(voxel_num_x=N, + voxel_num_y=N, + voxel_num_z=vert) + +# Set up empty image data +Phantom = ImageData(geometry=ig, + dimension_labels=['horizontal_x', + 'horizontal_y', + 'vertical']) +Phantom += 0.05 +# Populate image data by looping over and filling slices +i = 0 +while i < vert: + if vert > 1: + x = Phantom.subset(vertical=i).array + else: + x = Phantom.array + x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 + x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 0.94 + if vert > 1 : + Phantom.fill(x, vertical=i) + i += 1 + + +perc = 0.02 +# Set up empty image data +noise = ImageData(numpy.random.normal(loc = 0.04 , + scale = perc , + size = Phantom.shape), geometry=ig, + dimension_labels=['horizontal_x', + 'horizontal_y', + 'vertical']) +Phantom += noise + +# Set up AcquisitionGeometry object to hold the parameters of the measurement +# setup geometry: # Number of angles, the actual angles from 0 to +# pi for parallel beam, set the width of a detector +# pixel relative to an object pixe and the number of detector pixels. +angles_num = 20 +det_w = 1.0 +det_num = N + +angles = numpy.linspace(0,numpy.pi,angles_num,endpoint=False,dtype=numpy.float32)*\ + 180/numpy.pi + +# Inputs: Geometry, 2D or 3D, angles, horz detector pixel count, +# horz detector pixel size, vert detector pixel count, +# vert detector pixel size. +ag = AcquisitionGeometry('parallel', + '3D', + angles, + N, + det_w, + vert, + det_w) + +# Set up Operator object combining the ImageGeometry and AcquisitionGeometry +# wrapping calls to CCPi projector. +A = CCPiProjectorSimple(ig, ag) + +# Forward and backprojection are available as methods direct and adjoint. Here +# generate test data b and some noise + +b = A.direct(Phantom) + + +#z = A.adjoint(b) + + +# Using the test data b, different reconstruction methods can now be set up as +# demonstrated in the rest of this file. In general all methods need an initial +# guess and some algorithm options to be set. Note that 100 iterations for +# some of the methods is a very low number and 1000 or 10000 iterations may be +# needed if one wants to obtain a converged solution. +x_init = ImageData(geometry=ig, + dimension_labels=['horizontal_x','horizontal_y','vertical']) +X_init = BlockDataContainer(x_init) +B = BlockDataContainer(b, + ImageData(geometry=ig, dimension_labels=['horizontal_x','horizontal_y','vertical'])) + +# setup a tomo identity +Ibig = 1e5 * TomoIdentity(geometry=ig) +Ismall = 1e-5 * TomoIdentity(geometry=ig) +Iok = 1e1 * TomoIdentity(geometry=ig) + +# composite operator +Kbig = BlockOperator(A, Ibig, shape=(2,1)) +Ksmall = BlockOperator(A, Ismall, shape=(2,1)) +Kok = BlockOperator(A, Iok, shape=(2,1)) + +#out = K.direct(X_init) + +f = Norm2sq(Kbig,B) +f.L = 0.00003 + +fsmall = Norm2sq(Ksmall,B) +fsmall.L = 0.00003 + +fok = Norm2sq(Kok,B) +fok.L = 0.00003 + +simplef = Norm2sq(A, b) +simplef.L = 0.00003 + +gd = GradientDescent( x_init=x_init, objective_function=simplef, + rate=simplef.L) +gd.max_iteration = 50 + +Kbig.direct(X_init) +Kbig.adjoint(B) +cg = CGLS() +cg.set_up(X_init, Kbig, B ) +cg.max_iteration = 10 + +cgsmall = CGLS() +cgsmall.set_up(X_init, Ksmall, B ) +cgsmall.max_iteration = 10 + + +cgs = CGLS() +cgs.set_up(x_init, A, b ) +cgs.max_iteration = 10 + +cgok = CGLS() +cgok.set_up(X_init, Kok, B ) +cgok.max_iteration = 10 +# # +#out.__isub__(B) +#out2 = K.adjoint(out) + +#(2.0*self.c)*self.A.adjoint( self.A.direct(x) - self.b ) + +for _ in gd: + print ("iteration {} {}".format(gd.iteration, gd.get_last_loss())) + +cg.run(10, lambda it,val: print ("iteration {} objective {}".format(it,val))) + +cgs.run(10, lambda it,val: print ("iteration {} objective {}".format(it,val))) + +cgsmall.run(10, lambda it,val: print ("iteration {} objective {}".format(it,val))) +cgsmall.run(10, lambda it,val: print ("iteration {} objective {}".format(it,val))) +cgok.run(10, verbose=True) +# # for _ in cg: +# print ("iteration {} {}".format(cg.iteration, cg.get_current_loss())) +# # +# # fig = plt.figure() +# # plt.imshow(cg.get_output().get_item(0,0).subset(vertical=0).as_array()) +# # plt.title('Composite CGLS') +# # plt.show() +# # +# # for _ in cgs: +# print ("iteration {} {}".format(cgs.iteration, cgs.get_current_loss())) +# # +fig = plt.figure() +plt.subplot(2,3,1) +plt.imshow(Phantom.subset(vertical=0).as_array()) +plt.title('Simulated Phantom') +plt.subplot(2,3,2) +plt.imshow(gd.get_output().subset(vertical=0).as_array()) +plt.title('Simple Gradient Descent') +plt.subplot(2,3,3) +plt.imshow(cgs.get_output().subset(vertical=0).as_array()) +plt.title('Simple CGLS') +plt.subplot(2,3,5) +plt.imshow(cg.get_output().get_item(0).subset(vertical=0).as_array()) +plt.title('Composite CGLS\nbig lambda') +plt.subplot(2,3,6) +plt.imshow(cgsmall.get_output().get_item(0).subset(vertical=0).as_array()) +plt.title('Composite CGLS\nsmall lambda') +plt.subplot(2,3,4) +plt.imshow(cgok.get_output().get_item(0).subset(vertical=0).as_array()) +plt.title('Composite CGLS\nok lambda') +plt.show() diff --git a/Wrappers/Python/wip/CreatePhantom.py b/Wrappers/Python/wip/CreatePhantom.py new file mode 100644 index 0000000..4bf6ea4 --- /dev/null +++ b/Wrappers/Python/wip/CreatePhantom.py @@ -0,0 +1,242 @@ +import numpy +import tomophantom +from tomophantom import TomoP3D +from tomophantom.supp.artifacts import ArtifactsClass as Artifact +import os + +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import CGLS +from ccpi.plugins.ops import CCPiProjectorSimple +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, ImageData, AcquisitionData +from ccpi.optimisation.algorithms import GradientDescent +from ccpi.framework import BlockDataContainer +from ccpi.optimisation.operators import BlockOperator + + +model = 13 # select a model number from tomophantom library +N_size = 64 # Define phantom dimensions using a scalar value (cubic phantom) +path = os.path.dirname(tomophantom.__file__) +path_library3D = os.path.join(path, "Phantom3DLibrary.dat") + +#This will generate a N_size x N_size x N_size phantom (3D) +phantom_tm = TomoP3D.Model(model, N_size, path_library3D) + +# detector column count (horizontal) +detector_horiz = int(numpy.sqrt(2)*N_size) +# detector row count (vertical) (no reason for it to be > N) +detector_vert = N_size +# number of projection angles +angles_num = int(0.5*numpy.pi*N_size) +# angles are expressed in degrees +angles = numpy.linspace(0.0, 179.9, angles_num, dtype='float32') + + +acquisition_data_array = TomoP3D.ModelSino(model, N_size, + detector_horiz, detector_vert, + angles, + path_library3D) + +tomophantom_acquisition_axes_order = ['vertical', 'angle', 'horizontal'] + +artifacts = Artifact(acquisition_data_array) + + +tp_acq_data = AcquisitionData(artifacts.noise(0.2, 'Gaussian'), + dimension_labels=tomophantom_acquisition_axes_order) +#print ("size", acquisition_data.shape) +print ("horiz", detector_horiz) +print ("vert", detector_vert) +print ("angles", angles_num) + +tp_acq_geometry = AcquisitionGeometry(geom_type='parallel', dimension='3D', + angles=angles, + pixel_num_h=detector_horiz, + pixel_num_v=detector_vert, + channels=1, + ) + +acq_data = tp_acq_geometry.allocate() +#print (tp_acq_geometry) +print ("AcquisitionData", acq_data.shape) +print ("TomoPhantom", tp_acq_data.shape, tp_acq_data.dimension_labels) + +default_acquisition_axes_order = ['angle', 'vertical', 'horizontal'] + +acq_data2 = tp_acq_data.subset(dimensions=default_acquisition_axes_order) +print ("AcquisitionData", acq_data2.shape, acq_data2.dimension_labels) +print ("AcquisitionData {} TomoPhantom {}".format(id(acq_data2.as_array()), + id(acquisition_data_array))) + +fig = plt.figure() +plt.subplot(1,2,1) +plt.imshow(acquisition_data_array[20]) +plt.title('Sinogram') +plt.subplot(1,2,2) +plt.imshow(tp_acq_data.as_array()[20]) +plt.title('Sinogram + noise') +plt.show() + +# Set up Operator object combining the ImageGeometry and AcquisitionGeometry +# wrapping calls to CCPi projector. + +ig = ImageGeometry(voxel_num_x=detector_horiz, + voxel_num_y=detector_horiz, + voxel_num_z=detector_vert) +A = CCPiProjectorSimple(ig, tp_acq_geometry) +# Forward and backprojection are available as methods direct and adjoint. Here +# generate test data b and some noise + +#b = A.direct(Phantom) +b = acq_data2 + +#z = A.adjoint(b) + + +# Using the test data b, different reconstruction methods can now be set up as +# demonstrated in the rest of this file. In general all methods need an initial +# guess and some algorithm options to be set. Note that 100 iterations for +# some of the methods is a very low number and 1000 or 10000 iterations may be +# needed if one wants to obtain a converged solution. +x_init = ImageData(geometry=ig, + dimension_labels=['horizontal_x','horizontal_y','vertical']) +X_init = BlockDataContainer(x_init) +B = BlockDataContainer(b, + ImageData(geometry=ig, dimension_labels=['horizontal_x','horizontal_y','vertical'])) + +# setup a tomo identity +Ibig = 4e1 * TomoIdentity(geometry=ig) +Ismall = 1e-3 * TomoIdentity(geometry=ig) +Iok = 7.6e0 * TomoIdentity(geometry=ig) + +# composite operator +Kbig = BlockOperator(A, Ibig, shape=(2,1)) +Ksmall = BlockOperator(A, Ismall, shape=(2,1)) +Kok = BlockOperator(A, Iok, shape=(2,1)) + +#out = K.direct(X_init) +#x0 = x_init.copy() +#x0.fill(numpy.random.randn(*x0.shape)) +#lipschitz = PowerMethodNonsquare(A, 5, x0) +#print("lipschitz", lipschitz) + +#%% + +simplef = Norm2sq(A, b, memopt=False) +#simplef.L = lipschitz[0]/3000. +simplef.L = 0.00003 + +f = Norm2sq(Kbig,B) +f.L = 0.00003 + +fsmall = Norm2sq(Ksmall,B) +fsmall.L = 0.00003 + +fok = Norm2sq(Kok,B) +fok.L = 0.00003 + +print("setup gradient descent") +gd = GradientDescent( x_init=x_init, objective_function=simplef, + rate=simplef.L) +gd.max_iteration = 5 +simplef2 = Norm2sq(A, b, memopt=True) +#simplef.L = lipschitz[0]/3000. +simplef2.L = 0.00003 +print("setup gradient descent") +gd2 = GradientDescent( x_init=x_init, objective_function=simplef2, + rate=simplef2.L) +gd2.max_iteration = 5 + +Kbig.direct(X_init) +Kbig.adjoint(B) +print("setup CGLS") +cg = CGLS() +cg.set_up(X_init, Kbig, B ) +cg.max_iteration = 10 + +print("setup CGLS") +cgsmall = CGLS() +cgsmall.set_up(X_init, Ksmall, B ) +cgsmall.max_iteration = 10 + + +print("setup CGLS") +cgs = CGLS() +cgs.set_up(x_init, A, b ) +cgs.max_iteration = 10 + +print("setup CGLS") +cgok = CGLS() +cgok.set_up(X_init, Kok, B ) +cgok.max_iteration = 10 +# # +#out.__isub__(B) +#out2 = K.adjoint(out) + +#(2.0*self.c)*self.A.adjoint( self.A.direct(x) - self.b ) + + +for _ in gd: + print ("GradientDescent iteration {} {}".format(gd.iteration, gd.get_last_loss())) +#gd2.run(5,verbose=True) +print("CGLS block lambda big") +cg.run(10, lambda it,val: print ("CGLS big iteration {} objective {}".format(it,val))) + +print("CGLS standard") +cgs.run(10, lambda it,val: print ("CGLS standard iteration {} objective {}".format(it,val))) + +print("CGLS block lambda small") +cgsmall.run(10, lambda it,val: print ("CGLS small iteration {} objective {}".format(it,val))) +print("CGLS block lambdaok") +cgok.run(10, verbose=True) +# # for _ in cg: +# print ("iteration {} {}".format(cg.iteration, cg.get_current_loss())) +# # +# # fig = plt.figure() +# # plt.imshow(cg.get_output().get_item(0,0).subset(vertical=0).as_array()) +# # plt.title('Composite CGLS') +# # plt.show() +# # +# # for _ in cgs: +# print ("iteration {} {}".format(cgs.iteration, cgs.get_current_loss())) +# # +Phantom = ImageData(phantom_tm) + +theslice=40 + +fig = plt.figure() +plt.subplot(2,3,1) +plt.imshow(numpy.flip(Phantom.subset(vertical=theslice).as_array(),axis=0), cmap='gray') +plt.clim(0,0.7) +plt.title('Simulated Phantom') +plt.subplot(2,3,2) +plt.imshow(gd.get_output().subset(vertical=theslice).as_array(), cmap='gray') +plt.clim(0,0.7) +plt.title('Simple Gradient Descent') +plt.subplot(2,3,3) +plt.imshow(cgs.get_output().subset(vertical=theslice).as_array(), cmap='gray') +plt.clim(0,0.7) +plt.title('Simple CGLS') +plt.subplot(2,3,5) +plt.imshow(cg.get_output().get_item(0).subset(vertical=theslice).as_array(), cmap='gray') +plt.clim(0,0.7) +plt.title('Composite CGLS\nbig lambda') +plt.subplot(2,3,6) +plt.imshow(cgsmall.get_output().get_item(0).subset(vertical=theslice).as_array(), cmap='gray') +plt.clim(0,0.7) +plt.title('Composite CGLS\nsmall lambda') +plt.subplot(2,3,4) +plt.imshow(cgok.get_output().get_item(0).subset(vertical=theslice).as_array(), cmap='gray') +plt.clim(0,0.7) +plt.title('Composite CGLS\nok lambda') +plt.show() + + +#Ibig = 7e1 * TomoIdentity(geometry=ig) +#Kbig = BlockOperator(A, Ibig, shape=(2,1)) +#cg2 = CGLS(x_init=X_init, operator=Kbig, data=B) +#cg2.max_iteration = 10 +#cg2.run(10, verbose=True) |