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 (limited to 'Wrappers') 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(-) (limited to 'Wrappers') 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(-) (limited to 'Wrappers') 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(-) (limited to 'Wrappers') 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 (limited to 'Wrappers') 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 (limited to 'Wrappers') 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(-) (limited to 'Wrappers') 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 From a9ec78c1669f460ebaf5227600f8b0c082cfcf56 Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Tue, 30 Apr 2019 22:38:20 +0100 Subject: Fix dot product bug --- Wrappers/Python/ccpi/framework/framework.py | 3 +- Wrappers/Python/wip/compare_CGLS_algos.py | 127 ++++++++++++++++++++++++++++ 2 files changed, 129 insertions(+), 1 deletion(-) create mode 100644 Wrappers/Python/wip/compare_CGLS_algos.py (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py index ffc91ae..e278795 100755 --- a/Wrappers/Python/ccpi/framework/framework.py +++ b/Wrappers/Python/ccpi/framework/framework.py @@ -765,7 +765,8 @@ class DataContainer(object): def dot(self, other, *args, **kwargs): '''return the inner product of 2 DataContainers viewed as vectors''' if self.shape == other.shape: - return numpy.dot(self.as_array().ravel(), other.as_array().ravel()) + return (self*other).sum() + #return numpy.dot(self.as_array().ravel(), other.as_array().ravel()) else: raise ValueError('Shapes are not aligned: {} != {}'.format(self.shape, other.shape)) diff --git a/Wrappers/Python/wip/compare_CGLS_algos.py b/Wrappers/Python/wip/compare_CGLS_algos.py new file mode 100644 index 0000000..333805d --- /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.ops 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 -- cgit v1.2.3 From 6d1c24c8fc389365f2dd83ee480c63969b08ce9f Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 1 May 2019 16:13:40 +0100 Subject: removed useless imports --- Wrappers/Python/ccpi/optimisation/algs.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/algs.py b/Wrappers/Python/ccpi/optimisation/algs.py index 89b5519..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 -- cgit v1.2.3 From da581e6061ebe206e007fe4378cc9d449b67d791 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 1 May 2019 16:14:02 +0100 Subject: reimplements dot product following discussion in #273 an implementation of the dot product is made where we rely on Python to choose an appropriate type for the result of the reduction (e.g. float64 if the data is float32) --- Wrappers/Python/ccpi/framework/framework.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py index e278795..66420b9 100755 --- a/Wrappers/Python/ccpi/framework/framework.py +++ b/Wrappers/Python/ccpi/framework/framework.py @@ -764,13 +764,23 @@ 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 self.shape == other.shape: - return (self*other).sum() - #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)) - - + -- cgit v1.2.3 From 67fa9b5a985912232d8b5aeeec12451ed55a1a7a Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 1 May 2019 16:16:17 +0100 Subject: fix import --- Wrappers/Python/test/test_DataContainer.py | 5 +++++ Wrappers/Python/wip/compare_CGLS_algos.py | 2 +- 2 files changed, 6 insertions(+), 1 deletion(-) (limited to 'Wrappers') 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 index 333805d..119752c 100644 --- a/Wrappers/Python/wip/compare_CGLS_algos.py +++ b/Wrappers/Python/wip/compare_CGLS_algos.py @@ -5,7 +5,7 @@ from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, \ AcquisitionData from ccpi.optimisation.algs import FISTA, FBPD, CGLS, SIRT -from ccpi.astra.ops import AstraProjectorSimple +from ccpi.astra.operators import AstraProjectorSimple from ccpi.optimisation.algorithms import CGLS as CGLSalg -- cgit v1.2.3 From 82d94d608ea639c0aa8aefb80cc97c5d8b1ba2cb Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 1 May 2019 16:29:46 +0100 Subject: raise error if method is not recognised --- Wrappers/Python/ccpi/framework/framework.py | 3 +++ 1 file changed, 3 insertions(+) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py index 66420b9..7236e0e 100755 --- a/Wrappers/Python/ccpi/framework/framework.py +++ b/Wrappers/Python/ccpi/framework/framework.py @@ -765,6 +765,9 @@ class DataContainer(object): 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 (self*other).sum() if method == 'numpy': -- cgit v1.2.3