diff options
| -rwxr-xr-x | Wrappers/Python/ccpi/framework/framework.py | 24 | ||||
| -rwxr-xr-x | Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py | 2 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py | 74 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/optimisation/algorithms/__init__.py | 1 | ||||
| -rwxr-xr-x | Wrappers/Python/ccpi/optimisation/algs.py | 28 | ||||
| -rw-r--r-- | Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py | 201 | ||||
| -rwxr-xr-x | Wrappers/Python/test/test_DataContainer.py | 5 | ||||
| -rw-r--r-- | Wrappers/Python/wip/compare_CGLS_algos.py | 127 | ||||
| -rw-r--r-- | Wrappers/Python/wip/demo_SIRT.py (renamed from Wrappers/Python/wip/demo_test_sirt.py) | 147 | ||||
| -rw-r--r-- | Wrappers/Python/wip/demo_box_constraints_FISTA.py | 158 | ||||
| -rw-r--r-- | Wrappers/Python/wip/pdhg_TV_tomography2D.py | 3 |
11 files changed, 555 insertions, 215 deletions
diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py index 7516447..7236e0e 100755 --- a/Wrappers/Python/ccpi/framework/framework.py +++ b/Wrappers/Python/ccpi/framework/framework.py @@ -707,6 +707,10 @@ class DataContainer(object): def maximum(self, x2, *args, **kwargs): return self.pixel_wise_binary(numpy.maximum, x2, *args, **kwargs) + def minimum(self,x2, out=None, *args, **kwargs): + return self.pixel_wise_binary(numpy.minimum, x2=x2, out=out, *args, **kwargs) + + ## unary operations def pixel_wise_unary(self, pwop, *args, **kwargs): out = kwargs.get('out', None) @@ -760,12 +764,26 @@ class DataContainer(object): return numpy.sqrt(self.squared_norm()) def dot(self, other, *args, **kwargs): '''return the inner product of 2 DataContainers viewed as vectors''' + method = kwargs.get('method', 'reduce') + if method not in ['numpy','reduce']: + raise ValueError('dot: specified method not valid. Expecting numpy or reduce got {} '.format( + method)) if self.shape == other.shape: - return numpy.dot(self.as_array().ravel(), other.as_array().ravel()) + # return (self*other).sum() + if method == 'numpy': + return numpy.dot(self.as_array().ravel(), other.as_array()) + elif method == 'reduce': + # see https://github.com/vais-ral/CCPi-Framework/pull/273 + # notice that Python seems to be smart enough to use + # the appropriate type to hold the result of the reduction + sf = reduce(lambda x,y: x + y[0]*y[1], + zip(self.as_array().ravel(), + other.as_array().ravel()), + 0) + return sf else: raise ValueError('Shapes are not aligned: {} != {}'.format(self.shape, other.shape)) - - + diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py index c923a30..a14378c 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -52,7 +52,7 @@ class Algorithm(object): self.__loss = [] self.memopt = False self.timing = [] - self.update_objective_interval = 1 + self.update_objective_interval = kwargs.get('update_objective_interval', 1) def set_up(self, *args, **kwargs): '''Set up the algorithm''' raise NotImplementedError() diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py new file mode 100644 index 0000000..30584d4 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py @@ -0,0 +1,74 @@ +#!/usr/bin/env python3 +# -*- 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.optimisation.algorithms import Algorithm + +class SIRT(Algorithm): + + '''Simultaneous Iterative Reconstruction Technique + + Parameters: + x_init: initial guess + operator: operator for forward/backward projections + data: data to operate on + constraint: Function with prox-method, for example IndicatorBox to + enforce box constraints, default is None). + ''' + def __init__(self, **kwargs): + super(SIRT, self).__init__() + self.x = kwargs.get('x_init', None) + self.operator = kwargs.get('operator', None) + self.data = kwargs.get('data', None) + self.constraint = kwargs.get('constraint', None) + if self.x is not None and self.operator is not None and \ + self.data is not None: + print ("Calling from creator") + self.set_up(x_init=kwargs['x_init'], + operator=kwargs['operator'], + data=kwargs['data'], + constraint=kwargs['constraint']) + + def set_up(self, x_init, operator , data, constraint=None ): + + self.x = x_init.copy() + self.operator = operator + self.data = data + self.constraint = constraint + + self.r = data.copy() + + self.relax_par = 1.0 + + # Set up scaling matrices D and M. + self.M = 1/self.operator.direct(self.operator.domain_geometry().allocate(value=1.0)) + self.D = 1/self.operator.adjoint(self.operator.range_geometry().allocate(value=1.0)) + + + def update(self): + + self.r = self.data - self.operator.direct(self.x) + + self.x += self.relax_par * (self.D*self.operator.adjoint(self.M*self.r)) + + if self.constraint != None: + self.x = self.constraint.prox(self.x,None) + + def update_objective(self): + self.loss.append(self.r.squared_norm())
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py index f562973..2dbacfc 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py @@ -24,6 +24,7 @@ Created on Thu Feb 21 11:03:13 2019 from .Algorithm import Algorithm from .CGLS import CGLS +from .SIRT import SIRT from .GradientDescent import GradientDescent from .FISTA import FISTA from .FBPD import FBPD diff --git a/Wrappers/Python/ccpi/optimisation/algs.py b/Wrappers/Python/ccpi/optimisation/algs.py index 2f819d3..f5ba85e 100755 --- a/Wrappers/Python/ccpi/optimisation/algs.py +++ b/Wrappers/Python/ccpi/optimisation/algs.py @@ -20,13 +20,8 @@ import numpy import time -from ccpi.optimisation.functions import Function -from ccpi.optimisation.functions import ZeroFunction -from ccpi.framework import ImageData -from ccpi.framework import AcquisitionData -from ccpi.optimisation.spdhg import spdhg -from ccpi.optimisation.spdhg import KullbackLeibler -from ccpi.optimisation.spdhg import KullbackLeiblerConvexConjugate + + def FISTA(x_init, f=None, g=None, opt=None): '''Fast Iterative Shrinkage-Thresholding Algorithm @@ -280,10 +275,6 @@ def SIRT(x_init, operator , data , opt=None, constraint=None): tol = opt['tol'] max_iter = opt['iter'] - # Set default constraint to unconstrained - if constraint==None: - constraint = Function() - x = x_init.clone() timing = numpy.zeros(max_iter) @@ -293,21 +284,18 @@ def SIRT(x_init, operator , data , opt=None, constraint=None): relax_par = 1.0 # Set up scaling matrices D and M. - im1 = ImageData(geometry=x_init.geometry) - im1.array[:] = 1.0 - M = 1/operator.direct(im1) - del im1 - aq1 = AcquisitionData(geometry=M.geometry) - aq1.array[:] = 1.0 - D = 1/operator.adjoint(aq1) - del aq1 + M = 1/operator.direct(operator.domain_geometry().allocate(value=1.0)) + D = 1/operator.adjoint(operator.range_geometry().allocate(value=1.0)) # algorithm loop for it in range(0, max_iter): t = time.time() r = data - operator.direct(x) - x = constraint.prox(x + relax_par * (D*operator.adjoint(M*r)),None) + x = x + relax_par * (D*operator.adjoint(M*r)) + + if constraint != None: + x = constraint.prox(x,None) timing[it] = time.time() - t if it > 0: diff --git a/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py b/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py index 9de48a5..854f645 100644 --- a/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py +++ b/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py @@ -20,11 +20,17 @@ from ccpi.optimisation.operators import BlockOperator, Identity, Gradient from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ MixedL21Norm, BlockFunction, ScaledFunction -#from ccpi.astra.ops import AstraProjectorSimple -#from ccpi.plugins.ops import CCPiProjectorSimple from ccpi.plugins.operators import CCPiProjectorSimple -#from skimage.util import random_noise from timeit import default_timer as timer +from ccpi.reconstruction.parallelbeam import alg as pbalg +import os + +try: + import tomophantom + from tomophantom import TomoP3D + no_tomophantom = False +except ImportError as ie: + no_tomophantom = True #%% @@ -52,32 +58,11 @@ N = 75 vert = 4 ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, voxel_num_z=vert) -data = ig.allocate() -Phantom = data -# 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.98 - if vert > 1 : - Phantom.fill(x, vertical=i) - i += 1 - -#%% -#detectors = N -#angles = np.linspace(0,np.pi,100) -#angles_num = 100 -angles_num = N +angles_num = 100 det_w = 1.0 det_num = N -angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*\ - 180/np.pi angles = np.linspace(-90.,90.,N, dtype=np.float32) # Inputs: Geometry, 2D or 3D, angles, horz detector pixel count, # horz detector pixel size, vert detector pixel count, @@ -90,73 +75,59 @@ ag = AcquisitionGeometry('parallel', vert, det_w) -from ccpi.reconstruction.parallelbeam import alg as pbalg -from ccpi.plugins.processors import setupCCPiGeometries -def ssetupCCPiGeometries(ig, ag, counter): +#no_tomophantom = True +if no_tomophantom: + data = ig.allocate() + Phantom = data + # 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.98 + if vert > 1 : + Phantom.fill(x, vertical=i) + i += 1 - #vg = ImageGeometry(voxel_num_x=voxel_num_x,voxel_num_y=voxel_num_y, voxel_num_z=voxel_num_z) - #Phantom_ccpi = ImageData(geometry=vg, - # dimension_labels=['horizontal_x','horizontal_y','vertical']) - ##.subset(['horizontal_x','horizontal_y','vertical']) - ## ask the ccpi code what dimensions it would like - Phantom_ccpi = ig.allocate(dimension_labels=[ImageGeometry.HORIZONTAL_X, - ImageGeometry.HORIZONTAL_Y, - ImageGeometry.VERTICAL]) + Aop = CCPiProjectorSimple(ig, ag, 'cpu') + sin = Aop.direct(data) +else: + + model = 13 # select a model number from the library + N_size = N # 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) - voxel_per_pixel = 1 - angles = ag.angles - geoms = pbalg.pb_setup_geometry_from_image(Phantom_ccpi.as_array(), - angles, - voxel_per_pixel ) + #%% + Horiz_det = int(np.sqrt(2)*N_size) # detector column count (horizontal) + Vert_det = N_size # detector row count (vertical) (no reason for it to be > N) + #angles_num = int(0.5*np.pi*N_size); # angles number + #angles = np.linspace(0.0,179.9,angles_num,dtype='float32') # in degrees - pg = AcquisitionGeometry('parallel', - '3D', - angles, - geoms['n_h'], 1.0, - geoms['n_v'], 1.0 #2D in 3D is a slice 1 pixel thick - ) + print ("Building 3D analytical projection data with TomoPhantom") + projData3D_analyt = TomoP3D.ModelSino(model, + N_size, + Horiz_det, + Vert_det, + angles, + path_library3D) - center_of_rotation = Phantom_ccpi.get_dimension_size('horizontal_x') / 2 - #ad = AcquisitionData(geometry=pg,dimension_labels=['angle','vertical','horizontal']) - ad = pg.allocate(dimension_labels=[AcquisitionGeometry.ANGLE, - AcquisitionGeometry.VERTICAL, - AcquisitionGeometry.HORIZONTAL]) - geoms_i = pbalg.pb_setup_geometry_from_acquisition(ad.as_array(), - angles, - center_of_rotation, - voxel_per_pixel ) + # tomophantom outputs in [vert,angles,horiz] + # we want [angle,vert,horiz] + data = np.transpose(projData3D_analyt, [1,0,2]) + ag.pixel_num_h = Horiz_det + ag.pixel_num_v = Vert_det + sin = ag.allocate() + sin.fill(data) + ig.voxel_num_y = Vert_det - counter+=1 + Aop = CCPiProjectorSimple(ig, ag, 'cpu') - if counter < 4: - print (geoms, geoms_i) - if (not ( geoms_i == geoms )): - print ("not equal and {} {} {}".format(counter, geoms['output_volume_z'], geoms_i['output_volume_z'])) - X = max(geoms['output_volume_x'], geoms_i['output_volume_x']) - Y = max(geoms['output_volume_y'], geoms_i['output_volume_y']) - Z = max(geoms['output_volume_z'], geoms_i['output_volume_z']) - return setupCCPiGeometries(X,Y,Z,angles, counter) - else: - print ("happy now {} {} {}".format(counter, geoms['output_volume_z'], geoms_i['output_volume_z'])) - - return geoms - else: - return geoms_i - - - -#voxel_num_x, voxel_num_y, voxel_num_z, angles, counter -print ("###############################################") -print (ig) -print (ag) -g = setupCCPiGeometries(ig, ag, 0) -print (g) -print ("###############################################") -print ("###############################################") -#ag = AcquisitionGeometry('parallel','2D',angles, detectors) -#Aop = AstraProjectorSimple(ig, ag, 'cpu') -Aop = CCPiProjectorSimple(ig, ag, 'cpu') -sin = Aop.direct(data) plt.imshow(sin.subset(vertical=0).as_array()) plt.title('Sinogram') @@ -228,70 +199,40 @@ opt1 = {'niter':niter, 'memopt': True} -pdhg1 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma, memopt=True, max_iteration=niter) +pdhg1 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma, max_iteration=niter) #pdhg1.max_iteration = 2000 pdhg1.update_objective_interval = 100 -pdhg2 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma, memopt=False) -pdhg2.max_iteration = 2000 -pdhg2.update_objective_interval = 100 t1_old = timer() resold, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) t2_old = timer() -print ("memopt = False, shouldn't matter") pdhg1.run(niter) print (sum(pdhg1.timing)) res = pdhg1.get_output().subset(vertical=0) -print (pdhg1.objective) -t3 = timer() -#res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -print ("memopt = True, shouldn't matter") -pdhg2.run(niter) -print (sum(pdhg2.timing)) -res1 = pdhg2.get_output().subset(vertical=0) -t4 = timer() -# -print ("No memopt in {}s, memopt in {}/{}s old {}s".format(sum(pdhg1.timing), - sum(pdhg2.timing),t4-t3, t2_old-t1_old)) - -t1_old = timer() -resold1, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -t2_old = timer() #%% plt.figure() -plt.subplot(2,3,1) +plt.subplot(1,4,1) plt.imshow(res.as_array()) -plt.title('no memopt') -plt.colorbar() -plt.subplot(2,3,2) -plt.imshow(res1.as_array()) -plt.title('memopt') +plt.title('Algorithm') plt.colorbar() -plt.subplot(2,3,3) -plt.imshow((res1 - resold1.subset(vertical=0)).abs().as_array()) -plt.title('diff') -plt.colorbar() -plt.subplot(2,3,4) +plt.subplot(1,4,2) plt.imshow(resold.subset(vertical=0).as_array()) -plt.title('old nomemopt') +plt.title('function') plt.colorbar() -plt.subplot(2,3,5) -plt.imshow(resold1.subset(vertical=0).as_array()) -plt.title('old memopt') -plt.colorbar() -plt.subplot(2,3,6) -plt.imshow((resold1 - resold).subset(vertical=0).as_array()) -plt.title('diff old') +plt.subplot(1,4,3) +plt.imshow((res - resold.subset(vertical=0)).abs().as_array()) +plt.title('diff') plt.colorbar() -#plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt') -#plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt') -#plt.legend() +plt.subplot(1,4,4) +plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'Algorithm') +plt.plot(np.linspace(0,N,N), resold.subset(vertical=0).as_array()[int(N/2),:], label = 'function') +plt.legend() plt.show() # -print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(sum(pdhg1.timing), t4 -t3)) -diff = (res1 - res).abs().as_array().max() +print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(sum(pdhg1.timing), t2_old -t1_old)) +diff = (res - resold.subset(vertical=0)).abs().as_array().max() # print(" Max of abs difference is {}".format(diff)) diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py index 8e8ab87..e92d4c6 100755 --- a/Wrappers/Python/test/test_DataContainer.py +++ b/Wrappers/Python/test/test_DataContainer.py @@ -455,6 +455,11 @@ class TestDataContainer(unittest.TestCase): self.assertTrue(False) except ValueError as ve: self.assertTrue(True) + + print ("test dot numpy") + n0 = (ds0 * ds1).sum() + n1 = ds0.as_array().ravel().dot(ds1.as_array().ravel()) + self.assertEqual(n0, n1) diff --git a/Wrappers/Python/wip/compare_CGLS_algos.py b/Wrappers/Python/wip/compare_CGLS_algos.py new file mode 100644 index 0000000..119752c --- /dev/null +++ b/Wrappers/Python/wip/compare_CGLS_algos.py @@ -0,0 +1,127 @@ +# This demo illustrates how to use the SIRT algorithm without and with +# nonnegativity and box constraints. The ASTRA 2D projectors are used. + +# First make all imports +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, \ + AcquisitionData +from ccpi.optimisation.algs import FISTA, FBPD, CGLS, SIRT +from ccpi.astra.operators import AstraProjectorSimple + +from ccpi.optimisation.algorithms import CGLS as CGLSalg + +import numpy as np +import matplotlib.pyplot as plt + +# Choose either a parallel-beam (1=parallel2D) or fan-beam (2=cone2D) test case +test_case = 1 + +# Set up phantom size NxN by creating ImageGeometry, initialising the +# ImageData object with this geometry and empty array and finally put some +# data into its array, and display as image. +N = 128 +ig = ImageGeometry(voxel_num_x=N,voxel_num_y=N) +Phantom = ImageData(geometry=ig) + +x = Phantom.as_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)] = 1 + +#plt.figure() +#plt.imshow(x) +#plt.title('Phantom image') +#plt.show() + +# 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 and 0 to 2pi for fanbeam, set the width of a detector +# pixel relative to an object pixel, the number of detector pixels, and the +# source-origin and origin-detector distance (here the origin-detector distance +# set to 0 to simulate a "virtual detector" with same detector pixel size as +# object pixel size). +angles_num = 20 +det_w = 1.0 +det_num = N +SourceOrig = 200 +OrigDetec = 0 + +if test_case==1: + angles = np.linspace(0,np.pi,angles_num,endpoint=False) + ag = AcquisitionGeometry('parallel', + '2D', + angles, + det_num,det_w) +elif test_case==2: + angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) + ag = AcquisitionGeometry('cone', + '2D', + angles, + det_num, + det_w, + dist_source_center=SourceOrig, + dist_center_detector=OrigDetec) +else: + NotImplemented + +# Set up Operator object combining the ImageGeometry and AcquisitionGeometry +# wrapping calls to ASTRA as well as specifying whether to use CPU or GPU. +Aop = AstraProjectorSimple(ig, ag, 'cpu') + +# Forward and backprojection are available as methods direct and adjoint. Here +# generate test data b and do simple backprojection to obtain z. +b = Aop.direct(Phantom) +z = Aop.adjoint(b) + +#plt.figure() +#plt.imshow(b.array) +#plt.title('Simulated data') +#plt.show() + +#plt.figure() +#plt.imshow(z.array) +#plt.title('Backprojected data') +#plt.show() + +# 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: +x_init = ImageData(np.zeros(x.shape),geometry=ig) +opt = {'tol': 1e-4, 'iter': 7} + +# First a CGLS reconstruction using the function version of CGLS can be done: +x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt) + +#plt.figure() +#plt.imshow(x_CGLS.array) +#plt.title('CGLS') +#plt.colorbar() +#plt.show() + +#plt.figure() +#plt.semilogy(criter_CGLS) +#plt.title('CGLS criterion') +#plt.show() + + + +# Now CLGS using the algorithm class +CGLS_alg = CGLSalg() +CGLS_alg.set_up(x_init, Aop, b ) +CGLS_alg.max_iteration = 2000 +CGLS_alg.run(opt['iter']) +x_CGLS_alg = CGLS_alg.get_output() + +#plt.figure() +#plt.imshow(x_CGLS_alg.as_array()) +#plt.title('CGLS ALG') +#plt.colorbar() +#plt.show() + +#plt.figure() +#plt.semilogy(CGLS_alg.objective) +#plt.title('CGLS criterion') +#plt.show() + +print(criter_CGLS) +print(CGLS_alg.objective) + +print((x_CGLS - x_CGLS_alg).norm())
\ No newline at end of file diff --git a/Wrappers/Python/wip/demo_test_sirt.py b/Wrappers/Python/wip/demo_SIRT.py index 6f5a44d..5a85d41 100644 --- a/Wrappers/Python/wip/demo_test_sirt.py +++ b/Wrappers/Python/wip/demo_SIRT.py @@ -4,9 +4,9 @@ # First make all imports from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, \ AcquisitionData -from ccpi.optimisation.algs import FISTA, FBPD, CGLS, SIRT -from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D, IndicatorBox +from ccpi.optimisation.functions import IndicatorBox from ccpi.astra.ops import AstraProjectorSimple +from ccpi.optimisation.algorithms import SIRT, CGLS import numpy as np import matplotlib.pyplot as plt @@ -25,6 +25,7 @@ x = Phantom.as_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)] = 1 +plt.figure() plt.imshow(x) plt.title('Phantom image') plt.show() @@ -69,11 +70,13 @@ Aop = AstraProjectorSimple(ig, ag, 'gpu') b = Aop.direct(Phantom) z = Aop.adjoint(b) -plt.imshow(b.array) +plt.figure() +plt.imshow(b.as_array()) plt.title('Simulated data') plt.show() -plt.imshow(z.array) +plt.figure() +plt.imshow(z.as_array()) plt.title('Backprojected data') plt.show() @@ -81,96 +84,122 @@ plt.show() # demonstrated in the rest of this file. In general all methods need an initial # guess and some algorithm options to be set: x_init = ImageData(np.zeros(x.shape),geometry=ig) -opt = {'tol': 1e-4, 'iter': 1000} +opt = {'tol': 1e-4, 'iter': 100} -# First a CGLS reconstruction can be done: -x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt) -plt.imshow(x_CGLS.array) -plt.title('CGLS') +# First run a simple CGLS reconstruction: +CGLS_alg = CGLS() +CGLS_alg.set_up(x_init, Aop, b ) +CGLS_alg.max_iteration = 2000 +CGLS_alg.run(opt['iter']) +x_CGLS_alg = CGLS_alg.get_output() + +plt.figure() +plt.imshow(x_CGLS_alg.as_array()) +plt.title('CGLS ALG') plt.colorbar() plt.show() -plt.semilogy(criter_CGLS) +plt.figure() +plt.semilogy(CGLS_alg.objective) plt.title('CGLS criterion') plt.show() -# A SIRT unconstrained reconstruction can be done: similarly: -x_SIRT, it_SIRT, timing_SIRT, criter_SIRT = SIRT(x_init, Aop, b, opt) -plt.imshow(x_SIRT.array) +# A SIRT reconstruction can be done simply by replacing CGLS by SIRT. +# In the first instance, no constraints are enforced. +SIRT_alg = SIRT() +SIRT_alg.set_up(x_init, Aop, b ) +SIRT_alg.max_iteration = 2000 +SIRT_alg.run(opt['iter']) +x_SIRT_alg = SIRT_alg.get_output() + +plt.figure() +plt.imshow(x_SIRT_alg.as_array()) plt.title('SIRT unconstrained') plt.colorbar() plt.show() -plt.semilogy(criter_SIRT) +plt.figure() +plt.semilogy(SIRT_alg.objective) plt.title('SIRT unconstrained criterion') plt.show() -# A SIRT nonnegativity constrained reconstruction can be done using the -# additional input "constraint" set to a box indicator function with 0 as the -# lower bound and the default upper bound of infinity: -x_SIRT0, it_SIRT0, timing_SIRT0, criter_SIRT0 = SIRT(x_init, Aop, b, opt, - constraint=IndicatorBox(lower=0)) +# The SIRT algorithm is stopped after the specified number of iterations has +# been run. It can be resumed by calling the run command again, which will run +# it for the specificed number of iterations +SIRT_alg.run(opt['iter']) +x_SIRT_alg2 = SIRT_alg.get_output() -plt.imshow(x_SIRT0.array) -plt.title('SIRT nonneg') +plt.figure() +plt.imshow(x_SIRT_alg2.as_array()) +plt.title('SIRT unconstrained, extra iterations') plt.colorbar() plt.show() -plt.semilogy(criter_SIRT0) -plt.title('SIRT nonneg criterion') +plt.figure() +plt.semilogy(SIRT_alg.objective) +plt.title('SIRT unconstrained criterion, extra iterations') plt.show() -# A SIRT reconstruction with box constraints on [0,1] can also be done: -x_SIRT01, it_SIRT01, timing_SIRT01, criter_SIRT01 = SIRT(x_init, Aop, b, opt, - constraint=IndicatorBox(lower=0,upper=1)) -plt.imshow(x_SIRT01.array) -plt.title('SIRT box(0,1)') +# A SIRT nonnegativity constrained reconstruction can be done using the +# additional input "constraint" set to a box indicator function with 0 as the +# lower bound and the default upper bound of infinity. First setup a new +# instance of the SIRT algorithm. +SIRT_alg0 = SIRT() +SIRT_alg0.set_up(x_init, Aop, b, constraint=IndicatorBox(lower=0) ) +SIRT_alg0.max_iteration = 2000 +SIRT_alg0.run(opt['iter']) +x_SIRT_alg0 = SIRT_alg0.get_output() + +plt.figure() +plt.imshow(x_SIRT_alg0.as_array()) +plt.title('SIRT nonnegativity constrained') plt.colorbar() plt.show() -plt.semilogy(criter_SIRT01) -plt.title('SIRT box(0,1) criterion') -plt.show() - -# The indicator function can also be used with the FISTA algorithm to do -# least squares with nonnegativity constraint. - -# Create least squares object instance with projector, test data and a constant -# coefficient of 0.5: -f = Norm2sq(Aop,b,c=0.5) - -# Run FISTA for least squares without constraints -x_fista, it, timing, criter = FISTA(x_init, f, None,opt) - -plt.imshow(x_fista.array) -plt.title('FISTA Least squares') +plt.figure() +plt.semilogy(SIRT_alg0.objective) +plt.title('SIRT nonnegativity criterion') plt.show() -plt.semilogy(criter) -plt.title('FISTA Least squares criterion') -plt.show() -# Run FISTA for least squares with nonnegativity constraint -x_fista0, it0, timing0, criter0 = FISTA(x_init, f, IndicatorBox(lower=0),opt) +# A SIRT reconstruction with box constraints on [0,1] can also be done. +SIRT_alg01 = SIRT() +SIRT_alg01.set_up(x_init, Aop, b, constraint=IndicatorBox(lower=0,upper=1) ) +SIRT_alg01.max_iteration = 2000 +SIRT_alg01.run(opt['iter']) +x_SIRT_alg01 = SIRT_alg01.get_output() -plt.imshow(x_fista0.array) -plt.title('FISTA Least squares nonneg') +plt.figure() +plt.imshow(x_SIRT_alg01.as_array()) +plt.title('SIRT boc(0,1)') +plt.colorbar() plt.show() -plt.semilogy(criter0) -plt.title('FISTA Least squares nonneg criterion') +plt.figure() +plt.semilogy(SIRT_alg01.objective) +plt.title('SIRT box(0,1) criterion') plt.show() -# Run FISTA for least squares with box constraint [0,1] -x_fista01, it01, timing01, criter01 = FISTA(x_init, f, IndicatorBox(lower=0,upper=1),opt) +# The test image has values in the range [0,1], so enforcing values in the +# reconstruction to be within this interval improves a lot. Just for fun +# we can also easily see what happens if we choose a narrower interval as +# constrint in the reconstruction, lower bound 0.2, upper bound 0.8. +SIRT_alg0208 = SIRT() +SIRT_alg0208.set_up(x_init,Aop,b,constraint=IndicatorBox(lower=0.2,upper=0.8)) +SIRT_alg0208.max_iteration = 2000 +SIRT_alg0208.run(opt['iter']) +x_SIRT_alg0208 = SIRT_alg0208.get_output() -plt.imshow(x_fista01.array) -plt.title('FISTA Least squares box(0,1)') +plt.figure() +plt.imshow(x_SIRT_alg0208.as_array()) +plt.title('SIRT boc(0.2,0.8)') +plt.colorbar() plt.show() -plt.semilogy(criter01) -plt.title('FISTA Least squares box(0,1) criterion') +plt.figure() +plt.semilogy(SIRT_alg0208.objective) +plt.title('SIRT box(0.2,0.8) criterion') plt.show()
\ No newline at end of file diff --git a/Wrappers/Python/wip/demo_box_constraints_FISTA.py b/Wrappers/Python/wip/demo_box_constraints_FISTA.py new file mode 100644 index 0000000..2f9e0c6 --- /dev/null +++ b/Wrappers/Python/wip/demo_box_constraints_FISTA.py @@ -0,0 +1,158 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Apr 17 14:46:21 2019 + +@author: jakob + +Demonstrate the use of box constraints in FISTA +""" + +# First make all imports +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, \ + AcquisitionData +from ccpi.optimisation.algorithms import FISTA +from ccpi.optimisation.functions import Norm2sq, IndicatorBox +from ccpi.astra.ops import AstraProjectorSimple + +from ccpi.optimisation.operators import Identity + +import numpy as np +import matplotlib.pyplot as plt + + +# Set up phantom size NxN by creating ImageGeometry, initialising the +# ImageData object with this geometry and empty array and finally put some +# data into its array, and display as image. +N = 128 +ig = ImageGeometry(voxel_num_x=N,voxel_num_y=N) +Phantom = ImageData(geometry=ig) + +x = Phantom.as_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)] = 1 + +plt.figure() +plt.imshow(x) +plt.title('Phantom image') +plt.show() + +# 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 and 0 to 2pi for fanbeam, set the width of a detector +# pixel relative to an object pixel, the number of detector pixels, and the +# source-origin and origin-detector distance (here the origin-detector distance +# set to 0 to simulate a "virtual detector" with same detector pixel size as +# object pixel size). +angles_num = 20 +det_w = 1.0 +det_num = N +SourceOrig = 200 +OrigDetec = 0 + +test_case = 1 + +if test_case==1: + angles = np.linspace(0,np.pi,angles_num,endpoint=False) + ag = AcquisitionGeometry('parallel', + '2D', + angles, + det_num,det_w) +elif test_case==2: + angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) + ag = AcquisitionGeometry('cone', + '2D', + angles, + det_num, + det_w, + dist_source_center=SourceOrig, + dist_center_detector=OrigDetec) +else: + NotImplemented + +# Set up Operator object combining the ImageGeometry and AcquisitionGeometry +# wrapping calls to ASTRA as well as specifying whether to use CPU or GPU. +Aop = AstraProjectorSimple(ig, ag, 'gpu') + +Aop = Identity(ig,ig) + +# Forward and backprojection are available as methods direct and adjoint. Here +# generate test data b and do simple backprojection to obtain z. +b = Aop.direct(Phantom) +z = Aop.adjoint(b) + +plt.figure() +plt.imshow(b.array) +plt.title('Simulated data') +plt.show() + +plt.figure() +plt.imshow(z.array) +plt.title('Backprojected data') +plt.show() + +# 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: +x_init = ImageData(np.zeros(x.shape),geometry=ig) +opt = {'tol': 1e-4, 'iter': 100} + + + +# Create least squares object instance with projector, test data and a constant +# coefficient of 0.5: +f = Norm2sq(Aop,b,c=0.5) + +# Run FISTA for least squares without constraints +FISTA_alg = FISTA() +FISTA_alg.set_up(x_init=x_init, f=f, opt=opt) +FISTA_alg.max_iteration = 2000 +FISTA_alg.run(opt['iter']) +x_FISTA = FISTA_alg.get_output() + +plt.figure() +plt.imshow(x_FISTA.array) +plt.title('FISTA unconstrained') +plt.colorbar() +plt.show() + +plt.figure() +plt.semilogy(FISTA_alg.objective) +plt.title('FISTA unconstrained criterion') +plt.show() + +# Run FISTA for least squares with lower bound 0.1 +FISTA_alg0 = FISTA() +FISTA_alg0.set_up(x_init=x_init, f=f, g=IndicatorBox(lower=0.1), opt=opt) +FISTA_alg0.max_iteration = 2000 +FISTA_alg0.run(opt['iter']) +x_FISTA0 = FISTA_alg0.get_output() + +plt.figure() +plt.imshow(x_FISTA0.array) +plt.title('FISTA lower bound 0.1') +plt.colorbar() +plt.show() + +plt.figure() +plt.semilogy(FISTA_alg0.objective) +plt.title('FISTA criterion, lower bound 0.1') +plt.show() + +# Run FISTA for least squares with box constraint [0.1,0.8] +FISTA_alg0 = FISTA() +FISTA_alg0.set_up(x_init=x_init, f=f, g=IndicatorBox(lower=0.1,upper=0.8), opt=opt) +FISTA_alg0.max_iteration = 2000 +FISTA_alg0.run(opt['iter']) +x_FISTA0 = FISTA_alg0.get_output() + +plt.figure() +plt.imshow(x_FISTA0.array) +plt.title('FISTA box(0.1,0.8) constrained') +plt.colorbar() +plt.show() + +plt.figure() +plt.semilogy(FISTA_alg0.objective) +plt.title('FISTA criterion, box(0.1,0.8) constrained criterion') +plt.show()
\ No newline at end of file diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D.py b/Wrappers/Python/wip/pdhg_TV_tomography2D.py index cd91409..1be3dfa 100644 --- a/Wrappers/Python/wip/pdhg_TV_tomography2D.py +++ b/Wrappers/Python/wip/pdhg_TV_tomography2D.py @@ -19,8 +19,7 @@ from ccpi.optimisation.operators import BlockOperator, Gradient from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ MixedL21Norm, BlockFunction -from ccpi.astra.ops import AstraProjectorSimple -from skimage.util import random_noise +from ccpi.astra.operators import AstraProjectorSimple from timeit import default_timer as timer |