diff options
| -rwxr-xr-x | Wrappers/Python/ccpi/framework/framework.py | 4 | ||||
| -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 | 19 | ||||
| -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 |
6 files changed, 331 insertions, 72 deletions
diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py index 7516447..ffc91ae 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) 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..89b5519 100755 --- a/Wrappers/Python/ccpi/optimisation/algs.py +++ b/Wrappers/Python/ccpi/optimisation/algs.py @@ -280,10 +280,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 +289,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/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 |