summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rwxr-xr-xWrappers/Python/ccpi/optimisation/operators/BlockOperator.py353
-rwxr-xr-xWrappers/Python/test/test_BlockDataContainer.py178
-rw-r--r--Wrappers/Python/test/test_BlockOperator.py180
3 files changed, 446 insertions, 265 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
index a3ba93e..c9bf794 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
@@ -167,183 +167,184 @@ class BlockLinearOperator(BlockOperator):
if __name__ == '__main__':
+ pass
#from ccpi.optimisation.Algorithms import GradientDescent
- from ccpi.optimisation.algorithms import CGLS
+# 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
+# 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()
+# # 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()
diff --git a/Wrappers/Python/test/test_BlockDataContainer.py b/Wrappers/Python/test/test_BlockDataContainer.py
index 190fb21..ef11a82 100755
--- a/Wrappers/Python/test/test_BlockDataContainer.py
+++ b/Wrappers/Python/test/test_BlockDataContainer.py
@@ -21,7 +21,7 @@ import functools
class TestBlockDataContainer(unittest.TestCase):
def skiptest_BlockDataContainerShape(self):
print ("test block data container")
- ig0 = ImageGeometry(2,3,4)
+ ig0 = ImageGeometry(12,42,55,32)
ig1 = ImageGeometry(12,42,55,32)
data0 = ImageData(geometry=ig0)
@@ -37,7 +37,7 @@ class TestBlockDataContainer(unittest.TestCase):
def skiptest_BlockDataContainerShapeArithmetic(self):
print ("test block data container")
ig0 = ImageGeometry(2,3,4)
- ig1 = ImageGeometry(12,42,55,32)
+ ig1 = ImageGeometry(2,3,4)
data0 = ImageData(geometry=ig0)
data1 = ImageData(geometry=ig1) + 1
@@ -93,7 +93,7 @@ class TestBlockDataContainer(unittest.TestCase):
def test_BlockDataContainer(self):
print ("test block data container")
ig0 = ImageGeometry(2,3,4)
- ig1 = ImageGeometry(12,42,55,32)
+ ig1 = ImageGeometry(2,3,4)
data0 = ImageData(geometry=ig0)
data1 = ImageData(geometry=ig1) + 1
@@ -108,166 +108,166 @@ class TestBlockDataContainer(unittest.TestCase):
print (a[0][0].shape)
#cp2 = BlockDataContainer(*a)
cp2 = cp0.add(cp1)
- assert (cp2.get_item(0,0).as_array()[0][0][0] == 2.)
- assert (cp2.get_item(1,0).as_array()[0][0][0] == 4.)
+ 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,0).as_array()[0][0][0] == 2.)
- assert (cp2.get_item(1,0).as_array()[0][0][0] == 4.)
+ 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,0).as_array()[0][0][0] , 1. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2., decimal = 5)
+ 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,0).as_array()[0][0][0] , 1. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 3., decimal = 5)
+ 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,0).as_array()[0][0][0] , +3. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , +6., decimal = 5)
+ 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,0).as_array()[0][0][0] , +4. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , +7., decimal = 5)
+ 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,0).as_array()[0][0][0] , 2. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 6., decimal = 5)
+ 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,0).as_array()[0][0][0] == -2.)
- assert (cp2.get_item(1,0).as_array()[0][0][0] == -2.)
+ 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,0).as_array()[0][0][0] == -2.)
- assert (cp2.get_item(1,0).as_array()[0][0][0] == -2.)
+ 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,0).as_array()[0][0][0] , -1. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0, decimal = 5)
+ 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,0).as_array()[0][0][0] , -1. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -1., decimal = 5)
+ 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,0).as_array()[0][0][0] , -3. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -4., decimal = 5)
+ 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,0).as_array()[0][0][0] , -4. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -5., decimal = 5)
+ 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,0).as_array()[0][0][0] , -2. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -4., decimal = 5)
+ 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,0).as_array()[0][0][0] == 0.)
- assert (cp2.get_item(1,0).as_array()[0][0][0] == 3.)
+ 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,0).as_array()[0][0][0] == 0.)
- assert (cp2.get_item(1,0).as_array()[0][0][0] == 3.)
+ 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,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2, decimal = 5)
+ 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,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2, decimal = 5)
+ 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,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2., decimal = 5)
+ 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,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2., decimal = 5)
+ 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,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2., decimal = 5)
+ 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,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2., decimal = 5)
+ 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,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2., decimal = 5)
+ 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,0).as_array()[0][0][0] , 0 , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , +6., decimal = 5)
+ 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,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , +6., decimal = 5)
+ 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,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -6., decimal = 5)
+ 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,0).as_array()[0][0][0] == 0.)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1./3., decimal=4)
+ 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,0).as_array()[0][0][0] == 0.)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1./3., decimal=4)
+ 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,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0.5, decimal = 5)
+ 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,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0.5, decimal = 5)
+ 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,0).as_array()[0][0][0] , 0. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0.5, decimal = 5)
+ 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,0).as_array()[0][0][0] , 3. , decimal=5)
- numpy.testing.assert_almost_equal(cp3.get_item(1,0).as_array()[0][0][0] , 1, decimal = 5)
+ 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,0).as_array()[0][0][0] , 1./2 , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 1.5/3., decimal = 5)
+ 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,0).as_array()[0][0][0] , 0.5 , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0.5, decimal = 5)
+ 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,0).as_array()[0][0][0] , -0.5/2. , decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -0.5, decimal = 5)
+ 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,0).as_array()[0][0][0] == 0.)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1., decimal=4)
+ 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,0).as_array()[0][0][0] == 0.)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1., decimal=4)
+ 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,0).as_array()[0][0][0] , 0., decimal=5)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 1., decimal = 5)
+ 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,0).as_array()[0][0][0] == cp1.get_item(0,0).as_array()[0][0][0])
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], cp2.get_item(1,0).as_array()[0][0][0], decimal=4)
+ 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,0).as_array()[0][0][0], 0., decimal=4)
- numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1., decimal=4)
+ 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,0).as_array()[0][0][0], -1., decimal=4)
- numpy.testing.assert_almost_equal(s.get_item(1,0).as_array()[0][0][0], -1., decimal=4)
+ 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,0).as_array()[0][0][0], numpy.sqrt(2), decimal=4)
- numpy.testing.assert_almost_equal(s.get_item(1,0).as_array()[0][0][0], numpy.sqrt(4), decimal=4)
+ 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)
@@ -275,9 +275,9 @@ class TestBlockDataContainer(unittest.TestCase):
numpy.testing.assert_almost_equal(s, 0 + size, decimal=4)
s0 = 1
s1 = 1
- for i in cp0.get_item(0,0).shape:
+ for i in cp0.get_item(0).shape:
s0 *= i
- for i in cp0.get_item(1,0).shape:
+ 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)
diff --git a/Wrappers/Python/test/test_BlockOperator.py b/Wrappers/Python/test/test_BlockOperator.py
index 1896978..8bd673b 100644
--- a/Wrappers/Python/test/test_BlockOperator.py
+++ b/Wrappers/Python/test/test_BlockOperator.py
@@ -108,3 +108,183 @@ class TestBlockOperator(unittest.TestCase):
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()
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 14:33:17 +0000