From 162129a93b0a1d42cff1c1c44f1ce01e67e3b1fe Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Tue, 16 Apr 2019 16:12:05 +0100 Subject: Transferred from move_sirt_algorithm branch --- Wrappers/Python/ccpi/framework/framework.py | 4 + .../Python/ccpi/optimisation/algorithms/SIRT.py | 89 +++++++++++++++++++++ .../ccpi/optimisation/algorithms/__init__.py | 1 + Wrappers/Python/ccpi/optimisation/algs.py | 9 +-- Wrappers/Python/wip/demo_test_sirt.py | 90 +++++++++++++++++++--- 5 files changed, 176 insertions(+), 17 deletions(-) create mode 100644 Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py 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..389ec99 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py @@ -0,0 +1,89 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Apr 10 13:39:35 2019 + @author: jakob +""" + + # -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + + # Copyright 2018 Edoardo Pasca + + # Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + + # http://www.apache.org/licenses/LICENSE-2.0 + + # Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +""" +Created on Thu Feb 21 11:11:23 2019 + @author: ofn77899 +""" + + from ccpi.optimisation.algorithms import Algorithm +from ccpi.framework import ImageData, AcquisitionData + + #from collections.abc import Iterable +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. + ''' + 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('data', None) + if self.x is not None and self.operator is not None and \ + self.data is not None: + print ("Calling from creator") + self.set_up(x_init =kwargs['x_init'], + operator=kwargs['operator'], + data =kwargs['data']) + + def set_up(self, x_init, operator , data, 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. + im1 = ImageData(geometry=self.x.geometry) + im1.array[:] = 1.0 + self.M = 1/operator.direct(im1) + aq1 = AcquisitionData(geometry=self.M.geometry) + aq1.array[:] = 1.0 + self.D = 1/operator.adjoint(aq1) + + + 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..5a95c5c 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) @@ -307,7 +303,10 @@ def SIRT(x_init, operator , data , opt=None, constraint=None): 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_test_sirt.py index 6f5a44d..c7a3c76 100644 --- a/Wrappers/Python/wip/demo_test_sirt.py +++ b/Wrappers/Python/wip/demo_test_sirt.py @@ -8,6 +8,9 @@ from ccpi.optimisation.algs import FISTA, FBPD, CGLS, SIRT from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D, IndicatorBox from ccpi.astra.ops import AstraProjectorSimple +from ccpi.optimisation.algorithms import CGLS as CGLSALG +from ccpi.optimisation.algorithms import SIRT as SIRTALG + import numpy as np import matplotlib.pyplot as plt @@ -25,6 +28,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,10 +73,12 @@ Aop = AstraProjectorSimple(ig, ag, 'gpu') 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() @@ -81,96 +87,156 @@ 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.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() + +my_CGLS_alg = CGLSALG() +my_CGLS_alg.set_up(x_init, Aop, b ) +my_CGLS_alg.max_iteration = 2000 +my_CGLS_alg.run(opt['iter']) +x_CGLS_alg = my_CGLS_alg.get_output() + +plt.figure() +plt.imshow(x_CGLS_alg.array) +plt.title('CGLS ALG') +plt.colorbar() +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.figure() plt.imshow(x_SIRT.array) plt.title('SIRT unconstrained') plt.colorbar() plt.show() +plt.figure() plt.semilogy(criter_SIRT) plt.title('SIRT unconstrained criterion') plt.show() + + +my_SIRT_alg = SIRTALG() +my_SIRT_alg.set_up(x_init, Aop, b ) +my_SIRT_alg.max_iteration = 2000 +my_SIRT_alg.run(opt['iter']) +x_SIRT_alg = my_SIRT_alg.get_output() + +plt.figure() +plt.imshow(x_SIRT_alg.array) +plt.title('SIRT ALG') +plt.colorbar() +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)) - +plt.figure() plt.imshow(x_SIRT0.array) plt.title('SIRT nonneg') plt.colorbar() plt.show() +plt.figure() plt.semilogy(criter_SIRT0) plt.title('SIRT nonneg criterion') plt.show() + +my_SIRT_alg0 = SIRTALG() +my_SIRT_alg0.set_up(x_init, Aop, b, constraint=IndicatorBox(lower=0) ) +my_SIRT_alg0.max_iteration = 2000 +my_SIRT_alg0.run(opt['iter']) +x_SIRT_alg0 = my_SIRT_alg0.get_output() + +plt.figure() +plt.imshow(x_SIRT_alg0.array) +plt.title('SIRT ALG0') +plt.colorbar() +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.figure() plt.imshow(x_SIRT01.array) plt.title('SIRT box(0,1)') plt.colorbar() plt.show() +plt.figure() plt.semilogy(criter_SIRT01) plt.title('SIRT box(0,1) criterion') plt.show() +my_SIRT_alg01 = SIRTALG() +my_SIRT_alg01.set_up(x_init, Aop, b, constraint=IndicatorBox(lower=0,upper=1) ) +my_SIRT_alg01.max_iteration = 2000 +my_SIRT_alg01.run(opt['iter']) +x_SIRT_alg01 = my_SIRT_alg01.get_output() + +plt.figure() +plt.imshow(x_SIRT_alg01.array) +plt.title('SIRT ALG01') +plt.colorbar() +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.figure() plt.imshow(x_fista.array) plt.title('FISTA Least squares') plt.show() - +plt.figure() 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) - +plt.figure() plt.imshow(x_fista0.array) plt.title('FISTA Least squares nonneg') plt.show() - +plt.figure() plt.semilogy(criter0) plt.title('FISTA Least squares nonneg 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) - +plt.figure() plt.imshow(x_fista01.array) plt.title('FISTA Least squares box(0,1)') plt.show() - +plt.figure() plt.semilogy(criter01) plt.title('FISTA Least squares box(0,1) criterion') -plt.show() \ No newline at end of file +plt.show() +''' \ No newline at end of file -- cgit v1.2.3 From 490aedb50630a162b7914e27917c074c10bada8e Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Tue, 16 Apr 2019 16:17:12 +0100 Subject: bugfixes, whitespace --- Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py | 8 ++++---- Wrappers/Python/wip/demo_test_sirt.py | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py index 389ec99..9fc3b8e 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py @@ -28,15 +28,15 @@ Created on Thu Feb 21 11:11:23 2019 @author: ofn77899 """ - from ccpi.optimisation.algorithms import Algorithm +from ccpi.optimisation.algorithms import Algorithm from ccpi.framework import ImageData, AcquisitionData #from collections.abc import Iterable class SIRT(Algorithm): - '''Simultaneous Iterative Reconstruction Technique - Parameters: - x_init: initial guess + '''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 diff --git a/Wrappers/Python/wip/demo_test_sirt.py b/Wrappers/Python/wip/demo_test_sirt.py index c7a3c76..d3f27cf 100644 --- a/Wrappers/Python/wip/demo_test_sirt.py +++ b/Wrappers/Python/wip/demo_test_sirt.py @@ -5,7 +5,7 @@ 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.funcs import Norm2sq, Norm1, IndicatorBox from ccpi.astra.ops import AstraProjectorSimple from ccpi.optimisation.algorithms import CGLS as CGLSALG -- cgit v1.2.3 From fdc6c564722d6b93f62310d5ba57f03173028b7e Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Tue, 16 Apr 2019 16:19:39 +0100 Subject: Copied SIRT manually --- .../Python/ccpi/optimisation/algorithms/SIRT.py | 58 +++++++++++----------- 1 file changed, 30 insertions(+), 28 deletions(-) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py index 9fc3b8e..cb8b731 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py @@ -2,41 +2,44 @@ # -*- coding: utf-8 -*- """ Created on Wed Apr 10 13:39:35 2019 - @author: jakob + +@author: jakob """ - # -*- coding: utf-8 -*- +# -*- 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 +# Copyright 2018 Edoardo Pasca - # Licensed under the Apache License, Version 2.0 (the "License"); +# 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 +# http://www.apache.org/licenses/LICENSE-2.0 - # Unless required by applicable law or agreed to in writing, software +# Unless required by applicable law or agreed to in writing, software # distributed under the License is distributed on an "AS IS" BASIS, # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. """ Created on Thu Feb 21 11:11:23 2019 - @author: ofn77899 + +@author: ofn77899 """ from ccpi.optimisation.algorithms import Algorithm from ccpi.framework import ImageData, AcquisitionData - #from collections.abc import Iterable +#from collections.abc import Iterable class SIRT(Algorithm): '''Simultaneous Iterative Reconstruction Technique + Parameters: - x_init: initial guess + x_init: initial guess operator: operator for forward/backward projections data: data to operate on constraint: Function with prox-method, for example IndicatorBox to @@ -55,18 +58,18 @@ class SIRT(Algorithm): operator=kwargs['operator'], data =kwargs['data']) - def set_up(self, x_init, operator , data, constraint=None ): - - self.x = x_init.copy() + 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.r = data.copy() + + self.relax_par = 1.0 + + # Set up scaling matrices D and M. im1 = ImageData(geometry=self.x.geometry) im1.array[:] = 1.0 self.M = 1/operator.direct(im1) @@ -75,15 +78,14 @@ class SIRT(Algorithm): self.D = 1/operator.adjoint(aq1) - 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: + 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 + def update_objective(self): + self.loss.append(self.r.squared_norm()) \ No newline at end of file -- cgit v1.2.3 From 124bed2643557ba90b1df6cf4b09f4a6a5f30b47 Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Wed, 17 Apr 2019 14:53:38 +0100 Subject: Implemented using geometries in SIRT --- .../Python/ccpi/optimisation/algorithms/SIRT.py | 24 ++++------- Wrappers/Python/ccpi/optimisation/algs.py | 10 +---- Wrappers/Python/wip/demo_test_sirt.py | 46 +++------------------- 3 files changed, 16 insertions(+), 64 deletions(-) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py index cb8b731..beba913 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py @@ -24,11 +24,6 @@ Created on Wed Apr 10 13:39:35 2019 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. -""" -Created on Thu Feb 21 11:11:23 2019 - -@author: ofn77899 -""" from ccpi.optimisation.algorithms import Algorithm from ccpi.framework import ImageData, AcquisitionData @@ -43,20 +38,21 @@ class SIRT(Algorithm): operator: operator for forward/backward projections data: data to operate on constraint: Function with prox-method, for example IndicatorBox to - enforce box constraints. + 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('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']) + 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 ): @@ -70,12 +66,8 @@ class SIRT(Algorithm): self.relax_par = 1.0 # Set up scaling matrices D and M. - im1 = ImageData(geometry=self.x.geometry) - im1.array[:] = 1.0 - self.M = 1/operator.direct(im1) - aq1 = AcquisitionData(geometry=self.M.geometry) - aq1.array[:] = 1.0 - self.D = 1/operator.adjoint(aq1) + 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): diff --git a/Wrappers/Python/ccpi/optimisation/algs.py b/Wrappers/Python/ccpi/optimisation/algs.py index 5a95c5c..89b5519 100755 --- a/Wrappers/Python/ccpi/optimisation/algs.py +++ b/Wrappers/Python/ccpi/optimisation/algs.py @@ -289,14 +289,8 @@ 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): diff --git a/Wrappers/Python/wip/demo_test_sirt.py b/Wrappers/Python/wip/demo_test_sirt.py index d3f27cf..8f65f39 100644 --- a/Wrappers/Python/wip/demo_test_sirt.py +++ b/Wrappers/Python/wip/demo_test_sirt.py @@ -11,6 +11,8 @@ from ccpi.astra.ops import AstraProjectorSimple from ccpi.optimisation.algorithms import CGLS as CGLSALG from ccpi.optimisation.algorithms import SIRT as SIRTALG +from ccpi.optimisation.operators import Identity + import numpy as np import matplotlib.pyplot as plt @@ -68,6 +70,8 @@ else: # 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) @@ -89,6 +93,7 @@ plt.show() x_init = ImageData(np.zeros(x.shape),geometry=ig) 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) @@ -132,7 +137,6 @@ plt.title('SIRT unconstrained criterion') plt.show() - my_SIRT_alg = SIRTALG() my_SIRT_alg.set_up(x_init, Aop, b ) my_SIRT_alg.max_iteration = 2000 @@ -145,6 +149,7 @@ plt.title('SIRT ALG') plt.colorbar() 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: @@ -201,42 +206,3 @@ plt.imshow(x_SIRT_alg01.array) plt.title('SIRT ALG01') plt.colorbar() 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.figure() -plt.imshow(x_fista.array) -plt.title('FISTA Least squares') -plt.show() -plt.figure() -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) -plt.figure() -plt.imshow(x_fista0.array) -plt.title('FISTA Least squares nonneg') -plt.show() -plt.figure() -plt.semilogy(criter0) -plt.title('FISTA Least squares nonneg 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) -plt.figure() -plt.imshow(x_fista01.array) -plt.title('FISTA Least squares box(0,1)') -plt.show() -plt.figure() -plt.semilogy(criter01) -plt.title('FISTA Least squares box(0,1) criterion') -plt.show() -''' \ No newline at end of file -- cgit v1.2.3 From ae686e0f5b8baaab8e9dae00debf889c1fbd5921 Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Wed, 17 Apr 2019 15:32:22 +0100 Subject: Add demo FISTA with constraints previously in SIRT demo --- Wrappers/Python/wip/demo_box_constraints_FISTA.py | 158 ++++++++++++++++++++++ 1 file changed, 158 insertions(+) create mode 100644 Wrappers/Python/wip/demo_box_constraints_FISTA.py 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 -- cgit v1.2.3 From 90b2626c06ba69f002dabb1aae0238344646312e Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Thu, 25 Apr 2019 21:38:52 +0100 Subject: Tidy SIRT and demo --- .../Python/ccpi/optimisation/algorithms/SIRT.py | 9 - Wrappers/Python/wip/demo_SIRT.py | 205 ++++++++++++++++++++ Wrappers/Python/wip/demo_test_sirt.py | 208 --------------------- 3 files changed, 205 insertions(+), 217 deletions(-) create mode 100644 Wrappers/Python/wip/demo_SIRT.py delete mode 100644 Wrappers/Python/wip/demo_test_sirt.py diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py index beba913..30584d4 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py @@ -1,11 +1,4 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Wed Apr 10 13:39:35 2019 - -@author: jakob -""" - # -*- coding: utf-8 -*- # This work is part of the Core Imaging Library developed by # Visual Analytics and Imaging System Group of the Science Technology @@ -26,9 +19,7 @@ Created on Wed Apr 10 13:39:35 2019 # limitations under the License. from ccpi.optimisation.algorithms import Algorithm -from ccpi.framework import ImageData, AcquisitionData -#from collections.abc import Iterable class SIRT(Algorithm): '''Simultaneous Iterative Reconstruction Technique diff --git a/Wrappers/Python/wip/demo_SIRT.py b/Wrappers/Python/wip/demo_SIRT.py new file mode 100644 index 0000000..66b82a2 --- /dev/null +++ b/Wrappers/Python/wip/demo_SIRT.py @@ -0,0 +1,205 @@ +# 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.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 + +# 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, 'gpu') + +# 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} + + +# 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.array) +plt.title('CGLS ALG') +plt.colorbar() +plt.show() + +plt.figure() +plt.semilogy(CGLS_alg.loss) +plt.title('CGLS criterion') +plt.show() + + +# 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.array) +plt.title('SIRT unconstrained') +plt.colorbar() +plt.show() + +plt.figure() +plt.semilogy(SIRT_alg.loss) +plt.title('SIRT unconstrained criterion') +plt.show() + +# 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.figure() +plt.imshow(x_SIRT_alg2.array) +plt.title('SIRT unconstrained, extra iterations') +plt.colorbar() +plt.show() + +plt.figure() +plt.semilogy(SIRT_alg.loss) +plt.title('SIRT unconstrained criterion, extra iterations') +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. 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.array) +plt.title('SIRT nonnegativity constrained') +plt.colorbar() +plt.show() + +plt.figure() +plt.semilogy(SIRT_alg0.loss) +plt.title('SIRT nonnegativity criterion') +plt.show() + + +# 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.figure() +plt.imshow(x_SIRT_alg01.array) +plt.title('SIRT boc(0,1)') +plt.colorbar() +plt.show() + +plt.figure() +plt.semilogy(SIRT_alg01.loss) +plt.title('SIRT box(0,1) criterion') +plt.show() + +# 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.figure() +plt.imshow(x_SIRT_alg0208.array) +plt.title('SIRT boc(0.2,0.8)') +plt.colorbar() +plt.show() + +plt.figure() +plt.semilogy(SIRT_alg0208.loss) +plt.title('SIRT box(0.2,0.8) criterion') +plt.show() \ No newline at end of file diff --git a/Wrappers/Python/wip/demo_test_sirt.py b/Wrappers/Python/wip/demo_test_sirt.py deleted file mode 100644 index 8f65f39..0000000 --- a/Wrappers/Python/wip/demo_test_sirt.py +++ /dev/null @@ -1,208 +0,0 @@ -# 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.optimisation.funcs import Norm2sq, Norm1, IndicatorBox -from ccpi.astra.ops import AstraProjectorSimple - -from ccpi.optimisation.algorithms import CGLS as CGLSALG -from ccpi.optimisation.algorithms import SIRT as SIRTALG - -from ccpi.optimisation.operators import Identity - -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, '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} - - -# First a CGLS reconstruction 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() - - -my_CGLS_alg = CGLSALG() -my_CGLS_alg.set_up(x_init, Aop, b ) -my_CGLS_alg.max_iteration = 2000 -my_CGLS_alg.run(opt['iter']) -x_CGLS_alg = my_CGLS_alg.get_output() - -plt.figure() -plt.imshow(x_CGLS_alg.array) -plt.title('CGLS ALG') -plt.colorbar() -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.figure() -plt.imshow(x_SIRT.array) -plt.title('SIRT unconstrained') -plt.colorbar() -plt.show() - -plt.figure() -plt.semilogy(criter_SIRT) -plt.title('SIRT unconstrained criterion') -plt.show() - - -my_SIRT_alg = SIRTALG() -my_SIRT_alg.set_up(x_init, Aop, b ) -my_SIRT_alg.max_iteration = 2000 -my_SIRT_alg.run(opt['iter']) -x_SIRT_alg = my_SIRT_alg.get_output() - -plt.figure() -plt.imshow(x_SIRT_alg.array) -plt.title('SIRT ALG') -plt.colorbar() -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)) -plt.figure() -plt.imshow(x_SIRT0.array) -plt.title('SIRT nonneg') -plt.colorbar() -plt.show() - -plt.figure() -plt.semilogy(criter_SIRT0) -plt.title('SIRT nonneg criterion') -plt.show() - - -my_SIRT_alg0 = SIRTALG() -my_SIRT_alg0.set_up(x_init, Aop, b, constraint=IndicatorBox(lower=0) ) -my_SIRT_alg0.max_iteration = 2000 -my_SIRT_alg0.run(opt['iter']) -x_SIRT_alg0 = my_SIRT_alg0.get_output() - -plt.figure() -plt.imshow(x_SIRT_alg0.array) -plt.title('SIRT ALG0') -plt.colorbar() -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.figure() -plt.imshow(x_SIRT01.array) -plt.title('SIRT box(0,1)') -plt.colorbar() -plt.show() - -plt.figure() -plt.semilogy(criter_SIRT01) -plt.title('SIRT box(0,1) criterion') -plt.show() - -my_SIRT_alg01 = SIRTALG() -my_SIRT_alg01.set_up(x_init, Aop, b, constraint=IndicatorBox(lower=0,upper=1) ) -my_SIRT_alg01.max_iteration = 2000 -my_SIRT_alg01.run(opt['iter']) -x_SIRT_alg01 = my_SIRT_alg01.get_output() - -plt.figure() -plt.imshow(x_SIRT_alg01.array) -plt.title('SIRT ALG01') -plt.colorbar() -plt.show() -- cgit v1.2.3 From 347376fc17ddd6709e160f0c0c74577328afbbf1 Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Fri, 26 Apr 2019 12:57:48 +0100 Subject: Change array->as_array and loss->objective --- Wrappers/Python/wip/demo_SIRT.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/Wrappers/Python/wip/demo_SIRT.py b/Wrappers/Python/wip/demo_SIRT.py index 66b82a2..5a85d41 100644 --- a/Wrappers/Python/wip/demo_SIRT.py +++ b/Wrappers/Python/wip/demo_SIRT.py @@ -71,12 +71,12 @@ b = Aop.direct(Phantom) z = Aop.adjoint(b) plt.figure() -plt.imshow(b.array) +plt.imshow(b.as_array()) plt.title('Simulated data') plt.show() plt.figure() -plt.imshow(z.array) +plt.imshow(z.as_array()) plt.title('Backprojected data') plt.show() @@ -95,13 +95,13 @@ CGLS_alg.run(opt['iter']) x_CGLS_alg = CGLS_alg.get_output() plt.figure() -plt.imshow(x_CGLS_alg.array) +plt.imshow(x_CGLS_alg.as_array()) plt.title('CGLS ALG') plt.colorbar() plt.show() plt.figure() -plt.semilogy(CGLS_alg.loss) +plt.semilogy(CGLS_alg.objective) plt.title('CGLS criterion') plt.show() @@ -115,13 +115,13 @@ SIRT_alg.run(opt['iter']) x_SIRT_alg = SIRT_alg.get_output() plt.figure() -plt.imshow(x_SIRT_alg.array) +plt.imshow(x_SIRT_alg.as_array()) plt.title('SIRT unconstrained') plt.colorbar() plt.show() plt.figure() -plt.semilogy(SIRT_alg.loss) +plt.semilogy(SIRT_alg.objective) plt.title('SIRT unconstrained criterion') plt.show() @@ -132,13 +132,13 @@ SIRT_alg.run(opt['iter']) x_SIRT_alg2 = SIRT_alg.get_output() plt.figure() -plt.imshow(x_SIRT_alg2.array) +plt.imshow(x_SIRT_alg2.as_array()) plt.title('SIRT unconstrained, extra iterations') plt.colorbar() plt.show() plt.figure() -plt.semilogy(SIRT_alg.loss) +plt.semilogy(SIRT_alg.objective) plt.title('SIRT unconstrained criterion, extra iterations') plt.show() @@ -154,13 +154,13 @@ SIRT_alg0.run(opt['iter']) x_SIRT_alg0 = SIRT_alg0.get_output() plt.figure() -plt.imshow(x_SIRT_alg0.array) +plt.imshow(x_SIRT_alg0.as_array()) plt.title('SIRT nonnegativity constrained') plt.colorbar() plt.show() plt.figure() -plt.semilogy(SIRT_alg0.loss) +plt.semilogy(SIRT_alg0.objective) plt.title('SIRT nonnegativity criterion') plt.show() @@ -173,13 +173,13 @@ SIRT_alg01.run(opt['iter']) x_SIRT_alg01 = SIRT_alg01.get_output() plt.figure() -plt.imshow(x_SIRT_alg01.array) +plt.imshow(x_SIRT_alg01.as_array()) plt.title('SIRT boc(0,1)') plt.colorbar() plt.show() plt.figure() -plt.semilogy(SIRT_alg01.loss) +plt.semilogy(SIRT_alg01.objective) plt.title('SIRT box(0,1) criterion') plt.show() @@ -194,12 +194,12 @@ SIRT_alg0208.run(opt['iter']) x_SIRT_alg0208 = SIRT_alg0208.get_output() plt.figure() -plt.imshow(x_SIRT_alg0208.array) +plt.imshow(x_SIRT_alg0208.as_array()) plt.title('SIRT boc(0.2,0.8)') plt.colorbar() plt.show() plt.figure() -plt.semilogy(SIRT_alg0208.loss) +plt.semilogy(SIRT_alg0208.objective) plt.title('SIRT box(0.2,0.8) criterion') plt.show() \ No newline at end of file -- cgit v1.2.3