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/Python') 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/Python') 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/Python') 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/Python') 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/Python') 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 cc9e5d510378e1810d0a1a6b0fe70279389255be Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 17 Apr 2019 16:33:13 +0100 Subject: fixed conflict --- Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py | 4 ---- 1 file changed, 4 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py index cc3356a..c277482 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py @@ -20,12 +20,8 @@ import numpy from ccpi.optimisation.functions import Function from ccpi.optimisation.functions.ScaledFunction import ScaledFunction -<<<<<<< HEAD -from ccpi.framework import ImageData -======= from ccpi.framework import ImageData, ImageGeometry import functools ->>>>>>> composite_operator_datacontainer class KullbackLeibler(Function): -- cgit v1.2.3 From 9734788e8bc5088040ca8958db3ce9811a845758 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 23 Apr 2019 12:36:30 +0100 Subject: pass current solution to callback, fix PDHG setup --- Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py | 2 +- Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py | 7 ++++++- 2 files changed, 7 insertions(+), 2 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py index ed95c3f..7604c8a 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -151,7 +151,7 @@ class Algorithm(object): self.max_iteration, self.get_last_objective()) ) else: if callback is not None: - callback(self.iteration, self.get_last_objective()) + callback(self.iteration, self.get_last_objective(), self.x) i += 1 if i == iterations: break diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py index bc080f8..bc3059b 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -37,7 +37,12 @@ class PDHG(Algorithm): def set_up(self, f, g, operator, tau = None, sigma = None, opt = None, **kwargs): # algorithmic parameters - + self.operator = operator + self.f = f + self.g = g + self.tau = tau + self.sigma = sigma + self.opt = opt if sigma is None and tau is None: raise ValueError('Need sigma*tau||K||^2<1') -- cgit v1.2.3 From 541bbba333d0fd36e0a7a050879d40e6bd279306 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 23 Apr 2019 12:38:45 +0100 Subject: add tomography type of example with ccpi projector --- Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py | 273 ++++++++++++++++++++++ 1 file changed, 273 insertions(+) create mode 100644 Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py b/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py new file mode 100644 index 0000000..ea1ad08 --- /dev/null +++ b/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py @@ -0,0 +1,273 @@ +# -*- coding: utf-8 -*- + +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Feb 22 14:53:03 2019 + +@author: evangelos +""" + +from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, \ + AcquisitionGeometry, AcquisitionData + +import numpy as np +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG, PDHG_old + +from ccpi.optimisation.operators import BlockOperator, Identity, Gradient +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ + MixedL21Norm, BlockFunction, ScaledFunction + +#from ccpi.astra.ops import AstraProjectorSimple +from ccpi.plugins.ops import CCPiProjectorSimple +#from skimage.util import random_noise +from timeit import default_timer as timer + +#%% + +#%%############################################################################### +# Create phantom for TV tomography + +#import os +#import tomophantom +#from tomophantom import TomoP2D +#from tomophantom.supp.qualitymetrics import QualityTools + +#model = 1 # select a model number from the library +#N = 150 # set dimension of the phantom +## one can specify an exact path to the parameters file +## path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' +#path = os.path.dirname(tomophantom.__file__) +#path_library2D = os.path.join(path, "Phantom2DLibrary.dat") +##This will generate a N_size x N_size phantom (2D) +#phantom_2D = TomoP2D.Model(model, N, path_library2D) +#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +#data = ImageData(phantom_2D, geometry=ig) + +N = 150 +#x = np.zeros((N,N)) + +vert = 4 +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, voxel_num_z=vert) +data = ig.allocate() +Phantom = data +# Populate image data by looping over and filling slices +i = 0 +while i < vert: + if vert > 1: + x = Phantom.subset(vertical=i).array + else: + x = Phantom.array + x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 + x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 0.98 + if vert > 1 : + Phantom.fill(x, vertical=i) + i += 1 + + +#%% +#detectors = N +#angles = np.linspace(0,np.pi,100) +angles_num = 100 +det_w = 1.0 +det_num = N + +angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*\ + 180/np.pi + +# Inputs: Geometry, 2D or 3D, angles, horz detector pixel count, +# horz detector pixel size, vert detector pixel count, +# vert detector pixel size. +ag = AcquisitionGeometry('parallel', + '3D', + angles, + N, + det_w, + vert, + det_w) + +from ccpi.reconstruction.parallelbeam import alg as pbalg +#from ccpi.plugins.processors import setupCCPiGeometries +def setupCCPiGeometries(ig, ag, counter): + + #vg = ImageGeometry(voxel_num_x=voxel_num_x,voxel_num_y=voxel_num_y, voxel_num_z=voxel_num_z) + #Phantom_ccpi = ImageData(geometry=vg, + # dimension_labels=['horizontal_x','horizontal_y','vertical']) + ##.subset(['horizontal_x','horizontal_y','vertical']) + ## ask the ccpi code what dimensions it would like + Phantom_ccpi = ig.allocate(dimension_labels=[ImageGeometry.HORIZONTAL_X, + ImageGeometry.HORIZONTAL_Y, + ImageGeometry.VERTICAL]) + + voxel_per_pixel = 1 + angles = ag.angles + geoms = pbalg.pb_setup_geometry_from_image(Phantom_ccpi.as_array(), + angles, + voxel_per_pixel ) + + pg = AcquisitionGeometry('parallel', + '3D', + angles, + geoms['n_h'], 1.0, + geoms['n_v'], 1.0 #2D in 3D is a slice 1 pixel thick + ) + + center_of_rotation = Phantom_ccpi.get_dimension_size('horizontal_x') / 2 + #ad = AcquisitionData(geometry=pg,dimension_labels=['angle','vertical','horizontal']) + ad = pg.allocate(dimension_labels=[AcquisitionGeometry.ANGLE, + AcquisitionGeometry.VERTICAL, + AcquisitionGeometry.HORIZONTAL]) + geoms_i = pbalg.pb_setup_geometry_from_acquisition(ad.as_array(), + angles, + center_of_rotation, + voxel_per_pixel ) + + counter+=1 + + if counter < 4: + print (geoms, geoms_i) + if (not ( geoms_i == geoms )): + print ("not equal and {} {} {}".format(counter, geoms['output_volume_z'], geoms_i['output_volume_z'])) + X = max(geoms['output_volume_x'], geoms_i['output_volume_x']) + Y = max(geoms['output_volume_y'], geoms_i['output_volume_y']) + Z = max(geoms['output_volume_z'], geoms_i['output_volume_z']) + return setupCCPiGeometries(X,Y,Z,angles, counter) + else: + print ("happy now {} {} {}".format(counter, geoms['output_volume_z'], geoms_i['output_volume_z'])) + + return geoms + else: + return geoms_i + + + +#voxel_num_x, voxel_num_y, voxel_num_z, angles, counter +print ("###############################################") +print (ig) +print (ag) +g = setupCCPiGeometries(ig, ag, 0) +print (g) +print ("###############################################") +print ("###############################################") +#ag = AcquisitionGeometry('parallel','2D',angles, detectors) +#Aop = AstraProjectorSimple(ig, ag, 'cpu') +Aop = CCPiProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +plt.imshow(sin.subset(vertical=0).as_array()) +plt.title('Sinogram') +plt.colorbar() +plt.show() + + +#%% +# Add Gaussian noise to the sinogram data +np.random.seed(10) +n1 = np.random.random(sin.shape) + +noisy_data = sin + ImageData(5*n1) + +plt.imshow(noisy_data.subset(vertical=0).as_array()) +plt.title('Noisy Sinogram') +plt.colorbar() +plt.show() + + +#%% Works only with Composite Operator Structure of PDHG + +#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) + +# Create operators +op1 = Gradient(ig) +op2 = Aop + +# Form Composite Operator +operator = BlockOperator(op1, op2, shape=(2,1) ) + +alpha = 50 +f = BlockFunction( alpha * MixedL21Norm(), \ + 0.5 * L2NormSquared(b = noisy_data) ) +g = ZeroFunction() + +normK = Aop.norm() + +# Compute operator Norm +normK = operator.norm() + +## Primal & dual stepsizes +diag_precon = False + +if diag_precon: + + def tau_sigma_precond(operator): + + tau = 1/operator.sum_abs_row() + sigma = 1/ operator.sum_abs_col() + + return tau, sigma + + tau, sigma = tau_sigma_precond(operator) + +else: + sigma = 1 + tau = 1/(sigma*normK**2) + +# Compute operator Norm +normK = operator.norm() + +# Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) + +opt = {'niter':2000} +opt1 = {'niter':2000, 'memopt': True} + + +pdhg1 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma) +pdhg1.max_iteration = 2000 +pdhg1.update_objective_interval = 200 +pdhg2 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg2.max_iteration = 2000 +pdhg2.update_objective_interval = 200 +t1 = timer() +#res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) +pdhg1.run(200) +res = pdhg1.get_output().subset(vertical=0) +t2 = timer() + +t3 = timer() +#res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) +pdhg2.run(200) +res1 = pdhg2.get_output().subset(vertical=0) +t4 = timer() +# +print ("No memopt in {}s, memopt in {}s ".format(t2-t1, t4 -t3)) + +#%% +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(res.as_array()) +plt.title('no memopt') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(res1.as_array()) +plt.title('memopt') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow((res1 - res).abs().as_array()) +plt.title('diff') +plt.colorbar() +plt.show() +# +plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt') +plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt') +plt.legend() +plt.show() +# +print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(t2-t1, t4 -t3)) +diff = (res1 - res).abs().as_array().max() +# +print(" Max of abs difference is {}".format(diff)) + -- 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/Python') 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 cf4f909599f945c1af34daf00a9928dfeff4d041 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 26 Apr 2019 11:16:49 +0100 Subject: fixed prints --- .../ccpi/optimisation/algorithms/Algorithm.py | 23 +++++++--------------- 1 file changed, 7 insertions(+), 16 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py index 3c97480..bd48e13 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -34,7 +34,7 @@ class Algorithm(object): method will stop when the stopping cryterion is met. ''' - def __init__(self): + def __init__(self, **kwargs): '''Constructor Set the minimal number of parameters: @@ -48,7 +48,7 @@ class Algorithm(object): when evaluating the objective is computationally expensive. ''' self.iteration = 0 - self.__max_iteration = 0 + self.__max_iteration = kwargs.get('max_iteration', 0) self.__loss = [] self.memopt = False self.timing = [] @@ -91,9 +91,11 @@ class Algorithm(object): if self.iteration % self.update_objective_interval == 0: self.update_objective() self.iteration += 1 + def get_output(self): '''Returns the solution found''' return self.x + def get_last_loss(self): '''Returns the last stored value of the loss function @@ -146,21 +148,10 @@ class Algorithm(object): print ("Stop cryterion has been reached.") i = 0 - print("Iteration {:<5} Primal {:<5} Dual {:<5} PDgap".format('','','')) for _ in self: - - - if verbose and self.iteration % self.update_objective_interval == 0: - #pass - print( "{}/{} {:<5} {:.4f} {:<5} {:.4f} {:<5} {:.4f}".\ - format(self.iteration, self.max_iteration,'', \ - self.get_last_objective()[0],'',\ - self.get_last_objective()[1],'',\ - self.get_last_objective()[2])) - - - #print ("Iteration {}/{}, Primal, Dual, PDgap = {}".format(self.iteration, - # self.max_iteration, self.get_last_objective()) ) + if verbose and (self.iteration -1) % self.update_objective_interval == 0: + print ("Iteration {}/{}, = {}".format(self.iteration-1, + self.max_iteration, self.get_last_objective()) ) else: -- cgit v1.2.3 From a0c532bc1e55574aa5262901614545bdfb006dc7 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 26 Apr 2019 11:17:09 +0100 Subject: removed memopt as not needed --- Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py index 2cbd146..c0b774d 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -19,14 +19,13 @@ class PDHG(Algorithm): '''Primal Dual Hybrid Gradient''' def __init__(self, **kwargs): - super(PDHG, self).__init__() + super(PDHG, self).__init__(max_iteration=kwargs.get('max_iteration',0)) self.f = kwargs.get('f', None) self.operator = kwargs.get('operator', None) self.g = kwargs.get('g', None) self.tau = kwargs.get('tau', None) self.sigma = kwargs.get('sigma', None) - self.memopt = kwargs.get('memopt', False) - + if self.f is not None and self.operator is not None and \ self.g is not None: print ("Calling from creator") @@ -46,25 +45,22 @@ class PDHG(Algorithm): self.opt = opt if sigma is None and tau is None: raise ValueError('Need sigma*tau||K||^2<1') - - + self.x_old = self.operator.domain_geometry().allocate() self.x_tmp = self.x_old.copy() self.x = self.x_old.copy() - + self.y_old = self.operator.range_geometry().allocate() self.y_tmp = self.y_old.copy() self.y = self.y_old.copy() - + self.xbar = self.x_old.copy() - # relaxation parameter self.theta = 1 def update(self): - # Gradient descent, Dual problem solution self.operator.direct(self.xbar, out=self.y_tmp) self.y_tmp *= self.sigma @@ -90,7 +86,7 @@ class PDHG(Algorithm): self.y_old.fill(self.y) def update_objective(self): - + p1 = self.f(self.operator.direct(self.x)) + self.g(self.x) d1 = -(self.f.convex_conjugate(self.y) + self.g.convex_conjugate(-1*self.operator.adjoint(self.y))) -- cgit v1.2.3 From c39df4e997039c61f8a3bb883bf135d88db2498e Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 26 Apr 2019 11:17:32 +0100 Subject: demo updated --- Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py b/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py index 70ca57a..9de48a5 100644 --- a/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py +++ b/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py @@ -21,7 +21,8 @@ from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ MixedL21Norm, BlockFunction, ScaledFunction #from ccpi.astra.ops import AstraProjectorSimple -from ccpi.plugins.ops import CCPiProjectorSimple +#from ccpi.plugins.ops import CCPiProjectorSimple +from ccpi.plugins.operators import CCPiProjectorSimple #from skimage.util import random_noise from timeit import default_timer as timer @@ -90,8 +91,8 @@ ag = AcquisitionGeometry('parallel', det_w) from ccpi.reconstruction.parallelbeam import alg as pbalg -#from ccpi.plugins.processors import setupCCPiGeometries -def setupCCPiGeometries(ig, ag, counter): +from ccpi.plugins.processors import setupCCPiGeometries +def ssetupCCPiGeometries(ig, ag, counter): #vg = ImageGeometry(voxel_num_x=voxel_num_x,voxel_num_y=voxel_num_y, voxel_num_z=voxel_num_z) #Phantom_ccpi = ImageData(geometry=vg, @@ -221,16 +222,16 @@ normK = operator.norm() # Primal & dual stepsizes sigma = 1 tau = 1/(sigma*normK**2) -niter = 3 +niter = 50 opt = {'niter':niter} opt1 = {'niter':niter, 'memopt': True} -pdhg1 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma) -pdhg1.max_iteration = 2000 -pdhg1.update_objective_interval = 10 -pdhg2 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg1 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma, memopt=True, max_iteration=niter) +#pdhg1.max_iteration = 2000 +pdhg1.update_objective_interval = 100 +pdhg2 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma, memopt=False) pdhg2.max_iteration = 2000 pdhg2.update_objective_interval = 100 @@ -238,12 +239,14 @@ t1_old = timer() resold, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) t2_old = timer() +print ("memopt = False, shouldn't matter") pdhg1.run(niter) print (sum(pdhg1.timing)) res = pdhg1.get_output().subset(vertical=0) - +print (pdhg1.objective) t3 = timer() #res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) +print ("memopt = True, shouldn't matter") pdhg2.run(niter) print (sum(pdhg2.timing)) res1 = pdhg2.get_output().subset(vertical=0) -- 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/Python') 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 3a2fcc7e379f4761bad029dd99edf7d806526a10 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 26 Apr 2019 13:48:00 +0100 Subject: removed matplotlib dependence --- Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py | 1 - 1 file changed, 1 deletion(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py index c0b774d..39b092b 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -13,7 +13,6 @@ import time from ccpi.optimisation.operators import BlockOperator from ccpi.framework import BlockDataContainer from ccpi.optimisation.functions import FunctionOperatorComposition -import matplotlib.pyplot as plt class PDHG(Algorithm): '''Primal Dual Hybrid Gradient''' -- cgit v1.2.3 From 43d0a144e4d8766b8275d9c13e06f2d0423d198e Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 30 Apr 2019 13:55:28 +0100 Subject: removed numpy from build requirements --- Wrappers/Python/conda-recipe/meta.yaml | 2 -- 1 file changed, 2 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml index dd3238e..ac5f763 100644 --- a/Wrappers/Python/conda-recipe/meta.yaml +++ b/Wrappers/Python/conda-recipe/meta.yaml @@ -26,7 +26,6 @@ requirements: build: - {{ pin_compatible('numpy', max_pin='x.x') }} - python - - numpy - setuptools run: @@ -34,7 +33,6 @@ requirements: - python - numpy - scipy - #- matplotlib - h5py about: -- cgit v1.2.3 From 64ad60acb5c2a993e9b61682c18105723d9ae912 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 30 Apr 2019 13:56:46 +0100 Subject: restructured processors --- .../ccpi/processors/CenterOfRotationFinder.py | 408 +++++++++++++++++++++ Wrappers/Python/ccpi/processors/Normalizer.py | 124 +++++++ Wrappers/Python/ccpi/processors/__init__.py | 9 + Wrappers/Python/setup.py | 1 + 4 files changed, 542 insertions(+) create mode 100755 Wrappers/Python/ccpi/processors/CenterOfRotationFinder.py create mode 100755 Wrappers/Python/ccpi/processors/Normalizer.py create mode 100755 Wrappers/Python/ccpi/processors/__init__.py (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/ccpi/processors/CenterOfRotationFinder.py b/Wrappers/Python/ccpi/processors/CenterOfRotationFinder.py new file mode 100755 index 0000000..936dc05 --- /dev/null +++ b/Wrappers/Python/ccpi/processors/CenterOfRotationFinder.py @@ -0,0 +1,408 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2018 Edoardo Pasca + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License + +from ccpi.framework import DataProcessor, DataContainer, AcquisitionData,\ + AcquisitionGeometry, ImageGeometry, ImageData +import numpy +from scipy import ndimage + +class CenterOfRotationFinder(DataProcessor): + '''Processor to find the center of rotation in a parallel beam experiment + + This processor read in a AcquisitionDataSet and finds the center of rotation + based on Nghia Vo's method. https://doi.org/10.1364/OE.22.019078 + + Input: AcquisitionDataSet + + Output: float. center of rotation in pixel coordinate + ''' + + def __init__(self): + kwargs = { + + } + + #DataProcessor.__init__(self, **kwargs) + super(CenterOfRotationFinder, self).__init__(**kwargs) + + def check_input(self, dataset): + if dataset.number_of_dimensions == 3: + if dataset.geometry.geom_type == 'parallel': + return True + else: + raise ValueError('{0} is suitable only for parallel beam geometry'\ + .format(self.__class__.__name__)) + else: + raise ValueError("Expected input dimensions is 3, got {0}"\ + .format(dataset.number_of_dimensions)) + + + # ######################################################################### + # Copyright (c) 2015, UChicago Argonne, LLC. All rights reserved. # + # # + # Copyright 2015. UChicago Argonne, LLC. This software was produced # + # under U.S. Government contract DE-AC02-06CH11357 for Argonne National # + # Laboratory (ANL), which is operated by UChicago Argonne, LLC for the # + # U.S. Department of Energy. The U.S. Government has rights to use, # + # reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR # + # UChicago Argonne, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR # + # ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is # + # modified to produce derivative works, such modified software should # + # be clearly marked, so as not to confuse it with the version available # + # from ANL. # + # # + # Additionally, redistribution and use in source and binary forms, with # + # or without modification, are permitted provided that the following # + # conditions are met: # + # # + # * Redistributions of source code must retain the above copyright # + # notice, this list of conditions and the following disclaimer. # + # # + # * Redistributions in binary form must reproduce the above copyright # + # notice, this list of conditions and the following disclaimer in # + # the documentation and/or other materials provided with the # + # distribution. # + # # + # * Neither the name of UChicago Argonne, LLC, Argonne National # + # Laboratory, ANL, the U.S. Government, nor the names of its # + # contributors may be used to endorse or promote products derived # + # from this software without specific prior written permission. # + # # + # THIS SOFTWARE IS PROVIDED BY UChicago Argonne, LLC AND CONTRIBUTORS # + # "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT # + # LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS # + # FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL UChicago # + # Argonne, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, # + # INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, # + # BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; # + # LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER # + # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT # + # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN # + # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # + # POSSIBILITY OF SUCH DAMAGE. # + # ######################################################################### + + @staticmethod + def as_ndarray(arr, dtype=None, copy=False): + if not isinstance(arr, numpy.ndarray): + arr = numpy.array(arr, dtype=dtype, copy=copy) + return arr + + @staticmethod + def as_dtype(arr, dtype, copy=False): + if not arr.dtype == dtype: + arr = numpy.array(arr, dtype=dtype, copy=copy) + return arr + + @staticmethod + def as_float32(arr): + arr = CenterOfRotationFinder.as_ndarray(arr, numpy.float32) + return CenterOfRotationFinder.as_dtype(arr, numpy.float32) + + + + + @staticmethod + def find_center_vo(tomo, ind=None, smin=-40, smax=40, srad=10, step=0.5, + ratio=2., drop=20): + """ + Find rotation axis location using Nghia Vo's method. :cite:`Vo:14`. + + Parameters + ---------- + tomo : ndarray + 3D tomographic data. + ind : int, optional + Index of the slice to be used for reconstruction. + smin, smax : int, optional + Reference to the horizontal center of the sinogram. + srad : float, optional + Fine search radius. + step : float, optional + Step of fine searching. + ratio : float, optional + The ratio between the FOV of the camera and the size of object. + It's used to generate the mask. + drop : int, optional + Drop lines around vertical center of the mask. + + Returns + ------- + float + Rotation axis location. + + Notes + ----- + The function may not yield a correct estimate, if: + + - the sample size is bigger than the field of view of the camera. + In this case the ``ratio`` argument need to be set larger + than the default of 2.0. + + - there is distortion in the imaging hardware. If there's + no correction applied, the center of the projection image may + yield a better estimate. + + - the sample contrast is weak. Paganin's filter need to be applied + to overcome this. + + - the sample was changed during the scan. + """ + tomo = CenterOfRotationFinder.as_float32(tomo) + + if ind is None: + ind = tomo.shape[1] // 2 + _tomo = tomo[:, ind, :] + + + + # Reduce noise by smooth filters. Use different filters for coarse and fine search + _tomo_cs = ndimage.filters.gaussian_filter(_tomo, (3, 1)) + _tomo_fs = ndimage.filters.median_filter(_tomo, (2, 2)) + + # Coarse and fine searches for finding the rotation center. + if _tomo.shape[0] * _tomo.shape[1] > 4e6: # If data is large (>2kx2k) + #_tomo_coarse = downsample(numpy.expand_dims(_tomo_cs,1), level=2)[:, 0, :] + #init_cen = _search_coarse(_tomo_coarse, smin, smax, ratio, drop) + #fine_cen = _search_fine(_tomo_fs, srad, step, init_cen*4, ratio, drop) + init_cen = CenterOfRotationFinder._search_coarse(_tomo_cs, smin, + smax, ratio, drop) + fine_cen = CenterOfRotationFinder._search_fine(_tomo_fs, srad, + step, init_cen, + ratio, drop) + else: + init_cen = CenterOfRotationFinder._search_coarse(_tomo_cs, + smin, smax, + ratio, drop) + fine_cen = CenterOfRotationFinder._search_fine(_tomo_fs, srad, + step, init_cen, + ratio, drop) + + #logger.debug('Rotation center search finished: %i', fine_cen) + return fine_cen + + + @staticmethod + def _search_coarse(sino, smin, smax, ratio, drop): + """ + Coarse search for finding the rotation center. + """ + (Nrow, Ncol) = sino.shape + centerfliplr = (Ncol - 1.0) / 2.0 + + # Copy the sinogram and flip left right, the purpose is to + # make a full [0;2Pi] sinogram + _copy_sino = numpy.fliplr(sino[1:]) + + # This image is used for compensating the shift of sinogram 2 + temp_img = numpy.zeros((Nrow - 1, Ncol), dtype='float32') + temp_img[:] = sino[-1] + + # Start coarse search in which the shift step is 1 + listshift = numpy.arange(smin, smax + 1) + listmetric = numpy.zeros(len(listshift), dtype='float32') + mask = CenterOfRotationFinder._create_mask(2 * Nrow - 1, Ncol, + 0.5 * ratio * Ncol, drop) + for i in listshift: + _sino = numpy.roll(_copy_sino, i, axis=1) + if i >= 0: + _sino[:, 0:i] = temp_img[:, 0:i] + else: + _sino[:, i:] = temp_img[:, i:] + listmetric[i - smin] = numpy.sum(numpy.abs(numpy.fft.fftshift( + #pyfftw.interfaces.numpy_fft.fft2( + # numpy.vstack((sino, _sino))) + numpy.fft.fft2(numpy.vstack((sino, _sino))) + )) * mask) + minpos = numpy.argmin(listmetric) + return centerfliplr + listshift[minpos] / 2.0 + + @staticmethod + def _search_fine(sino, srad, step, init_cen, ratio, drop): + """ + Fine search for finding the rotation center. + """ + Nrow, Ncol = sino.shape + centerfliplr = (Ncol + 1.0) / 2.0 - 1.0 + # Use to shift the sinogram 2 to the raw CoR. + shiftsino = numpy.int16(2 * (init_cen - centerfliplr)) + _copy_sino = numpy.roll(numpy.fliplr(sino[1:]), shiftsino, axis=1) + if init_cen <= centerfliplr: + lefttake = numpy.int16(numpy.ceil(srad + 1)) + righttake = numpy.int16(numpy.floor(2 * init_cen - srad - 1)) + else: + lefttake = numpy.int16(numpy.ceil( + init_cen - (Ncol - 1 - init_cen) + srad + 1)) + righttake = numpy.int16(numpy.floor(Ncol - 1 - srad - 1)) + Ncol1 = righttake - lefttake + 1 + mask = CenterOfRotationFinder._create_mask(2 * Nrow - 1, Ncol1, + 0.5 * ratio * Ncol, drop) + numshift = numpy.int16((2 * srad) / step) + 1 + listshift = numpy.linspace(-srad, srad, num=numshift) + listmetric = numpy.zeros(len(listshift), dtype='float32') + factor1 = numpy.mean(sino[-1, lefttake:righttake]) + num1 = 0 + for i in listshift: + _sino = ndimage.interpolation.shift( + _copy_sino, (0, i), prefilter=False) + factor2 = numpy.mean(_sino[0,lefttake:righttake]) + _sino = _sino * factor1 / factor2 + sinojoin = numpy.vstack((sino, _sino)) + listmetric[num1] = numpy.sum(numpy.abs(numpy.fft.fftshift( + #pyfftw.interfaces.numpy_fft.fft2( + # sinojoin[:, lefttake:righttake + 1]) + numpy.fft.fft2(sinojoin[:, lefttake:righttake + 1]) + )) * mask) + num1 = num1 + 1 + minpos = numpy.argmin(listmetric) + return init_cen + listshift[minpos] / 2.0 + + @staticmethod + def _create_mask(nrow, ncol, radius, drop): + du = 1.0 / ncol + dv = (nrow - 1.0) / (nrow * 2.0 * numpy.pi) + centerrow = numpy.ceil(nrow / 2) - 1 + centercol = numpy.ceil(ncol / 2) - 1 + # added by Edoardo Pasca + centerrow = int(centerrow) + centercol = int(centercol) + mask = numpy.zeros((nrow, ncol), dtype='float32') + for i in range(nrow): + num1 = numpy.round(((i - centerrow) * dv / radius) / du) + (p1, p2) = numpy.int16(numpy.clip(numpy.sort( + (-num1 + centercol, num1 + centercol)), 0, ncol - 1)) + mask[i, p1:p2 + 1] = numpy.ones(p2 - p1 + 1, dtype='float32') + if drop < centerrow: + mask[centerrow - drop:centerrow + drop + 1, + :] = numpy.zeros((2 * drop + 1, ncol), dtype='float32') + mask[:,centercol-1:centercol+2] = numpy.zeros((nrow, 3), dtype='float32') + return mask + + def process(self, out=None): + + projections = self.get_input() + + cor = CenterOfRotationFinder.find_center_vo(projections.as_array()) + + return cor + + +class AcquisitionDataPadder(DataProcessor): + '''Normalization based on flat and dark + + This processor read in a AcquisitionData and normalises it based on + the instrument reading with and without incident photons or neutrons. + + Input: AcquisitionData + Parameter: 2D projection with flat field (or stack) + 2D projection with dark field (or stack) + Output: AcquisitionDataSetn + ''' + + def __init__(self, + center_of_rotation = None, + acquisition_geometry = None, + pad_value = 1e-5): + kwargs = { + 'acquisition_geometry' : acquisition_geometry, + 'center_of_rotation' : center_of_rotation, + 'pad_value' : pad_value + } + + super(AcquisitionDataPadder, self).__init__(**kwargs) + + def check_input(self, dataset): + if self.acquisition_geometry is None: + self.acquisition_geometry = dataset.geometry + if dataset.number_of_dimensions == 3: + return True + else: + raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ + .format(dataset.number_of_dimensions)) + + def process(self, out=None): + projections = self.get_input() + w = projections.get_dimension_size('horizontal') + delta = w - 2 * self.center_of_rotation + + padded_width = int ( + numpy.ceil(abs(delta)) + w + ) + delta_pix = padded_width - w + + voxel_per_pixel = 1 + geom = pbalg.pb_setup_geometry_from_acquisition(projections.as_array(), + self.acquisition_geometry.angles, + self.center_of_rotation, + voxel_per_pixel ) + + padded_geometry = self.acquisition_geometry.clone() + + padded_geometry.pixel_num_h = geom['n_h'] + padded_geometry.pixel_num_v = geom['n_v'] + + delta_pix_h = padded_geometry.pixel_num_h - self.acquisition_geometry.pixel_num_h + delta_pix_v = padded_geometry.pixel_num_v - self.acquisition_geometry.pixel_num_v + + if delta_pix_h == 0: + delta_pix_h = delta_pix + padded_geometry.pixel_num_h = padded_width + #initialize a new AcquisitionData with values close to 0 + out = AcquisitionData(geometry=padded_geometry) + out = out + self.pad_value + + + #pad in the horizontal-vertical plane -> slice on angles + if delta > 0: + #pad left of middle + command = "out.array[" + for i in range(out.number_of_dimensions): + if out.dimension_labels[i] == 'horizontal': + value = '{0}:{1}'.format(delta_pix_h, delta_pix_h+w) + command = command + str(value) + else: + if out.dimension_labels[i] == 'vertical' : + value = '{0}:'.format(delta_pix_v) + command = command + str(value) + else: + command = command + ":" + if i < out.number_of_dimensions -1: + command = command + ',' + command = command + '] = projections.array' + #print (command) + else: + #pad right of middle + command = "out.array[" + for i in range(out.number_of_dimensions): + if out.dimension_labels[i] == 'horizontal': + value = '{0}:{1}'.format(0, w) + command = command + str(value) + else: + if out.dimension_labels[i] == 'vertical' : + value = '{0}:'.format(delta_pix_v) + command = command + str(value) + else: + command = command + ":" + if i < out.number_of_dimensions -1: + command = command + ',' + command = command + '] = projections.array' + #print (command) + #cleaned = eval(command) + exec(command) + return out \ No newline at end of file diff --git a/Wrappers/Python/ccpi/processors/Normalizer.py b/Wrappers/Python/ccpi/processors/Normalizer.py new file mode 100755 index 0000000..da65e5c --- /dev/null +++ b/Wrappers/Python/ccpi/processors/Normalizer.py @@ -0,0 +1,124 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2018 Edoardo Pasca + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License + +from ccpi.framework import DataProcessor, DataContainer, AcquisitionData,\ + AcquisitionGeometry, ImageGeometry, ImageData +import numpy + +class Normalizer(DataProcessor): + '''Normalization based on flat and dark + + This processor read in a AcquisitionData and normalises it based on + the instrument reading with and without incident photons or neutrons. + + Input: AcquisitionData + Parameter: 2D projection with flat field (or stack) + 2D projection with dark field (or stack) + Output: AcquisitionDataSetn + ''' + + def __init__(self, flat_field = None, dark_field = None, tolerance = 1e-5): + kwargs = { + 'flat_field' : flat_field, + 'dark_field' : dark_field, + # very small number. Used when there is a division by zero + 'tolerance' : tolerance + } + + #DataProcessor.__init__(self, **kwargs) + super(Normalizer, self).__init__(**kwargs) + if not flat_field is None: + self.set_flat_field(flat_field) + if not dark_field is None: + self.set_dark_field(dark_field) + + def check_input(self, dataset): + if dataset.number_of_dimensions == 3 or\ + dataset.number_of_dimensions == 2: + return True + else: + raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ + .format(dataset.number_of_dimensions)) + + def set_dark_field(self, df): + if type(df) is numpy.ndarray: + if len(numpy.shape(df)) == 3: + raise ValueError('Dark Field should be 2D') + elif len(numpy.shape(df)) == 2: + self.dark_field = df + elif issubclass(type(df), DataContainer): + self.dark_field = self.set_dark_field(df.as_array()) + + def set_flat_field(self, df): + if type(df) is numpy.ndarray: + if len(numpy.shape(df)) == 3: + raise ValueError('Flat Field should be 2D') + elif len(numpy.shape(df)) == 2: + self.flat_field = df + elif issubclass(type(df), DataContainer): + self.flat_field = self.set_flat_field(df.as_array()) + + @staticmethod + def normalize_projection(projection, flat, dark, tolerance): + a = (projection - dark) + b = (flat-dark) + with numpy.errstate(divide='ignore', invalid='ignore'): + c = numpy.true_divide( a, b ) + c[ ~ numpy.isfinite( c )] = tolerance # set to not zero if 0/0 + return c + + @staticmethod + def estimate_normalised_error(projection, flat, dark, delta_flat, delta_dark): + '''returns the estimated relative error of the normalised projection + + n = (projection - dark) / (flat - dark) + Dn/n = (flat-dark + projection-dark)/((flat-dark)*(projection-dark))*(Df/f + Dd/d) + ''' + a = (projection - dark) + b = (flat-dark) + df = delta_flat / flat + dd = delta_dark / dark + rel_norm_error = (b + a) / (b * a) * (df + dd) + return rel_norm_error + + def process(self, out=None): + + projections = self.get_input() + dark = self.dark_field + flat = self.flat_field + + if projections.number_of_dimensions == 3: + if not (projections.shape[1:] == dark.shape and \ + projections.shape[1:] == flat.shape): + raise ValueError('Flats/Dark and projections size do not match.') + + + a = numpy.asarray( + [ Normalizer.normalize_projection( + projection, flat, dark, self.tolerance) \ + for projection in projections.as_array() ] + ) + elif projections.number_of_dimensions == 2: + a = Normalizer.normalize_projection(projections.as_array(), + flat, dark, self.tolerance) + y = type(projections)( a , True, + dimension_labels=projections.dimension_labels, + geometry=projections.geometry) + return y + \ No newline at end of file diff --git a/Wrappers/Python/ccpi/processors/__init__.py b/Wrappers/Python/ccpi/processors/__init__.py new file mode 100755 index 0000000..f8d050e --- /dev/null +++ b/Wrappers/Python/ccpi/processors/__init__.py @@ -0,0 +1,9 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Apr 30 13:51:09 2019 + +@author: ofn77899 +""" + +from .CenterOfRotationFinder import CenterOfRotationFinder +from .Normalizer import Normalizer diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py index a3fde59..78f88dd 100644 --- a/Wrappers/Python/setup.py +++ b/Wrappers/Python/setup.py @@ -36,6 +36,7 @@ setup( 'ccpi.optimisation.operators', 'ccpi.optimisation.algorithms', 'ccpi.optimisation.functions', + 'ccpi.optimisation.processors', 'ccpi.contrib','ccpi.contrib.optimisation', 'ccpi.contrib.optimisation.algorithms'], -- cgit v1.2.3 From 136d908afbac2b1010b56f8b96081df956f2b8bf Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 30 Apr 2019 13:57:42 +0100 Subject: added test for center of rotation finder --- Wrappers/Python/test/test_DataProcessor.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index 1c1de3a..a5bd9c6 100755 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -11,8 +11,29 @@ from timeit import default_timer as timer from ccpi.framework import AX, CastDataContainer, PixelByPixelDataProcessor +from ccpi.io.reader import NexusReader +from ccpi.processors import CenterOfRotationFinder +import wget +import os + class TestDataProcessor(unittest.TestCase): + def setUp(self): + wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs') + self.filename = '24737_fd.nxs' + + def tearDown(self): + os.remove(self.filename) + def test_CenterOfRotation(self): + reader = NexusReader(self.filename) + ad = reader.get_acquisition_data_whole() + print (ad.geometry) + cf = CenterOfRotationFinder() + cf.set_input(ad) + print ("Center of rotation", cf.get_output()) + + + def test_DataProcessorChaining(self): shape = (2,3,4,5) size = shape[0] -- cgit v1.2.3 From d4abd5cec7097caeda5d30a7c891fd6a47162a8a Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 30 Apr 2019 14:21:52 +0100 Subject: added test and removed processors.py --- Wrappers/Python/ccpi/processors.py | 514 ----------------------------- Wrappers/Python/setup.py | 2 +- Wrappers/Python/test/test_DataProcessor.py | 3 + 3 files changed, 4 insertions(+), 515 deletions(-) delete mode 100755 Wrappers/Python/ccpi/processors.py (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/ccpi/processors.py b/Wrappers/Python/ccpi/processors.py deleted file mode 100755 index ccef410..0000000 --- a/Wrappers/Python/ccpi/processors.py +++ /dev/null @@ -1,514 +0,0 @@ -# -*- coding: utf-8 -*- -# This work is part of the Core Imaging Library developed by -# Visual Analytics and Imaging System Group of the Science Technology -# Facilities Council, STFC - -# Copyright 2018 Edoardo Pasca - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License - -from ccpi.framework import DataProcessor, DataContainer, AcquisitionData,\ - AcquisitionGeometry, ImageGeometry, ImageData -from ccpi.reconstruction.parallelbeam import alg as pbalg -import numpy -from scipy import ndimage - -import matplotlib.pyplot as plt - - -class Normalizer(DataProcessor): - '''Normalization based on flat and dark - - This processor read in a AcquisitionData and normalises it based on - the instrument reading with and without incident photons or neutrons. - - Input: AcquisitionData - Parameter: 2D projection with flat field (or stack) - 2D projection with dark field (or stack) - Output: AcquisitionDataSetn - ''' - - def __init__(self, flat_field = None, dark_field = None, tolerance = 1e-5): - kwargs = { - 'flat_field' : flat_field, - 'dark_field' : dark_field, - # very small number. Used when there is a division by zero - 'tolerance' : tolerance - } - - #DataProcessor.__init__(self, **kwargs) - super(Normalizer, self).__init__(**kwargs) - if not flat_field is None: - self.set_flat_field(flat_field) - if not dark_field is None: - self.set_dark_field(dark_field) - - def check_input(self, dataset): - if dataset.number_of_dimensions == 3 or\ - dataset.number_of_dimensions == 2: - return True - else: - raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ - .format(dataset.number_of_dimensions)) - - def set_dark_field(self, df): - if type(df) is numpy.ndarray: - if len(numpy.shape(df)) == 3: - raise ValueError('Dark Field should be 2D') - elif len(numpy.shape(df)) == 2: - self.dark_field = df - elif issubclass(type(df), DataContainer): - self.dark_field = self.set_dark_field(df.as_array()) - - def set_flat_field(self, df): - if type(df) is numpy.ndarray: - if len(numpy.shape(df)) == 3: - raise ValueError('Flat Field should be 2D') - elif len(numpy.shape(df)) == 2: - self.flat_field = df - elif issubclass(type(df), DataContainer): - self.flat_field = self.set_flat_field(df.as_array()) - - @staticmethod - def normalize_projection(projection, flat, dark, tolerance): - a = (projection - dark) - b = (flat-dark) - with numpy.errstate(divide='ignore', invalid='ignore'): - c = numpy.true_divide( a, b ) - c[ ~ numpy.isfinite( c )] = tolerance # set to not zero if 0/0 - return c - - @staticmethod - def estimate_normalised_error(projection, flat, dark, delta_flat, delta_dark): - '''returns the estimated relative error of the normalised projection - - n = (projection - dark) / (flat - dark) - Dn/n = (flat-dark + projection-dark)/((flat-dark)*(projection-dark))*(Df/f + Dd/d) - ''' - a = (projection - dark) - b = (flat-dark) - df = delta_flat / flat - dd = delta_dark / dark - rel_norm_error = (b + a) / (b * a) * (df + dd) - return rel_norm_error - - def process(self, out=None): - - projections = self.get_input() - dark = self.dark_field - flat = self.flat_field - - if projections.number_of_dimensions == 3: - if not (projections.shape[1:] == dark.shape and \ - projections.shape[1:] == flat.shape): - raise ValueError('Flats/Dark and projections size do not match.') - - - a = numpy.asarray( - [ Normalizer.normalize_projection( - projection, flat, dark, self.tolerance) \ - for projection in projections.as_array() ] - ) - elif projections.number_of_dimensions == 2: - a = Normalizer.normalize_projection(projections.as_array(), - flat, dark, self.tolerance) - y = type(projections)( a , True, - dimension_labels=projections.dimension_labels, - geometry=projections.geometry) - return y - - -class CenterOfRotationFinder(DataProcessor): - '''Processor to find the center of rotation in a parallel beam experiment - - This processor read in a AcquisitionDataSet and finds the center of rotation - based on Nghia Vo's method. https://doi.org/10.1364/OE.22.019078 - - Input: AcquisitionDataSet - - Output: float. center of rotation in pixel coordinate - ''' - - def __init__(self): - kwargs = { - - } - - #DataProcessor.__init__(self, **kwargs) - super(CenterOfRotationFinder, self).__init__(**kwargs) - - def check_input(self, dataset): - if dataset.number_of_dimensions == 3: - if dataset.geometry.geom_type == 'parallel': - return True - else: - raise ValueError('{0} is suitable only for parallel beam geometry'\ - .format(self.__class__.__name__)) - else: - raise ValueError("Expected input dimensions is 3, got {0}"\ - .format(dataset.number_of_dimensions)) - - - # ######################################################################### - # Copyright (c) 2015, UChicago Argonne, LLC. All rights reserved. # - # # - # Copyright 2015. UChicago Argonne, LLC. This software was produced # - # under U.S. Government contract DE-AC02-06CH11357 for Argonne National # - # Laboratory (ANL), which is operated by UChicago Argonne, LLC for the # - # U.S. Department of Energy. The U.S. Government has rights to use, # - # reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR # - # UChicago Argonne, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR # - # ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is # - # modified to produce derivative works, such modified software should # - # be clearly marked, so as not to confuse it with the version available # - # from ANL. # - # # - # Additionally, redistribution and use in source and binary forms, with # - # or without modification, are permitted provided that the following # - # conditions are met: # - # # - # * Redistributions of source code must retain the above copyright # - # notice, this list of conditions and the following disclaimer. # - # # - # * Redistributions in binary form must reproduce the above copyright # - # notice, this list of conditions and the following disclaimer in # - # the documentation and/or other materials provided with the # - # distribution. # - # # - # * Neither the name of UChicago Argonne, LLC, Argonne National # - # Laboratory, ANL, the U.S. Government, nor the names of its # - # contributors may be used to endorse or promote products derived # - # from this software without specific prior written permission. # - # # - # THIS SOFTWARE IS PROVIDED BY UChicago Argonne, LLC AND CONTRIBUTORS # - # "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT # - # LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS # - # FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL UChicago # - # Argonne, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, # - # INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, # - # BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; # - # LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER # - # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT # - # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN # - # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # - # POSSIBILITY OF SUCH DAMAGE. # - # ######################################################################### - - @staticmethod - def as_ndarray(arr, dtype=None, copy=False): - if not isinstance(arr, numpy.ndarray): - arr = numpy.array(arr, dtype=dtype, copy=copy) - return arr - - @staticmethod - def as_dtype(arr, dtype, copy=False): - if not arr.dtype == dtype: - arr = numpy.array(arr, dtype=dtype, copy=copy) - return arr - - @staticmethod - def as_float32(arr): - arr = CenterOfRotationFinder.as_ndarray(arr, numpy.float32) - return CenterOfRotationFinder.as_dtype(arr, numpy.float32) - - - - - @staticmethod - def find_center_vo(tomo, ind=None, smin=-40, smax=40, srad=10, step=0.5, - ratio=2., drop=20): - """ - Find rotation axis location using Nghia Vo's method. :cite:`Vo:14`. - - Parameters - ---------- - tomo : ndarray - 3D tomographic data. - ind : int, optional - Index of the slice to be used for reconstruction. - smin, smax : int, optional - Reference to the horizontal center of the sinogram. - srad : float, optional - Fine search radius. - step : float, optional - Step of fine searching. - ratio : float, optional - The ratio between the FOV of the camera and the size of object. - It's used to generate the mask. - drop : int, optional - Drop lines around vertical center of the mask. - - Returns - ------- - float - Rotation axis location. - - Notes - ----- - The function may not yield a correct estimate, if: - - - the sample size is bigger than the field of view of the camera. - In this case the ``ratio`` argument need to be set larger - than the default of 2.0. - - - there is distortion in the imaging hardware. If there's - no correction applied, the center of the projection image may - yield a better estimate. - - - the sample contrast is weak. Paganin's filter need to be applied - to overcome this. - - - the sample was changed during the scan. - """ - tomo = CenterOfRotationFinder.as_float32(tomo) - - if ind is None: - ind = tomo.shape[1] // 2 - _tomo = tomo[:, ind, :] - - - - # Reduce noise by smooth filters. Use different filters for coarse and fine search - _tomo_cs = ndimage.filters.gaussian_filter(_tomo, (3, 1)) - _tomo_fs = ndimage.filters.median_filter(_tomo, (2, 2)) - - # Coarse and fine searches for finding the rotation center. - if _tomo.shape[0] * _tomo.shape[1] > 4e6: # If data is large (>2kx2k) - #_tomo_coarse = downsample(numpy.expand_dims(_tomo_cs,1), level=2)[:, 0, :] - #init_cen = _search_coarse(_tomo_coarse, smin, smax, ratio, drop) - #fine_cen = _search_fine(_tomo_fs, srad, step, init_cen*4, ratio, drop) - init_cen = CenterOfRotationFinder._search_coarse(_tomo_cs, smin, - smax, ratio, drop) - fine_cen = CenterOfRotationFinder._search_fine(_tomo_fs, srad, - step, init_cen, - ratio, drop) - else: - init_cen = CenterOfRotationFinder._search_coarse(_tomo_cs, - smin, smax, - ratio, drop) - fine_cen = CenterOfRotationFinder._search_fine(_tomo_fs, srad, - step, init_cen, - ratio, drop) - - #logger.debug('Rotation center search finished: %i', fine_cen) - return fine_cen - - - @staticmethod - def _search_coarse(sino, smin, smax, ratio, drop): - """ - Coarse search for finding the rotation center. - """ - (Nrow, Ncol) = sino.shape - centerfliplr = (Ncol - 1.0) / 2.0 - - # Copy the sinogram and flip left right, the purpose is to - # make a full [0;2Pi] sinogram - _copy_sino = numpy.fliplr(sino[1:]) - - # This image is used for compensating the shift of sinogram 2 - temp_img = numpy.zeros((Nrow - 1, Ncol), dtype='float32') - temp_img[:] = sino[-1] - - # Start coarse search in which the shift step is 1 - listshift = numpy.arange(smin, smax + 1) - listmetric = numpy.zeros(len(listshift), dtype='float32') - mask = CenterOfRotationFinder._create_mask(2 * Nrow - 1, Ncol, - 0.5 * ratio * Ncol, drop) - for i in listshift: - _sino = numpy.roll(_copy_sino, i, axis=1) - if i >= 0: - _sino[:, 0:i] = temp_img[:, 0:i] - else: - _sino[:, i:] = temp_img[:, i:] - listmetric[i - smin] = numpy.sum(numpy.abs(numpy.fft.fftshift( - #pyfftw.interfaces.numpy_fft.fft2( - # numpy.vstack((sino, _sino))) - numpy.fft.fft2(numpy.vstack((sino, _sino))) - )) * mask) - minpos = numpy.argmin(listmetric) - return centerfliplr + listshift[minpos] / 2.0 - - @staticmethod - def _search_fine(sino, srad, step, init_cen, ratio, drop): - """ - Fine search for finding the rotation center. - """ - Nrow, Ncol = sino.shape - centerfliplr = (Ncol + 1.0) / 2.0 - 1.0 - # Use to shift the sinogram 2 to the raw CoR. - shiftsino = numpy.int16(2 * (init_cen - centerfliplr)) - _copy_sino = numpy.roll(numpy.fliplr(sino[1:]), shiftsino, axis=1) - if init_cen <= centerfliplr: - lefttake = numpy.int16(numpy.ceil(srad + 1)) - righttake = numpy.int16(numpy.floor(2 * init_cen - srad - 1)) - else: - lefttake = numpy.int16(numpy.ceil( - init_cen - (Ncol - 1 - init_cen) + srad + 1)) - righttake = numpy.int16(numpy.floor(Ncol - 1 - srad - 1)) - Ncol1 = righttake - lefttake + 1 - mask = CenterOfRotationFinder._create_mask(2 * Nrow - 1, Ncol1, - 0.5 * ratio * Ncol, drop) - numshift = numpy.int16((2 * srad) / step) + 1 - listshift = numpy.linspace(-srad, srad, num=numshift) - listmetric = numpy.zeros(len(listshift), dtype='float32') - factor1 = numpy.mean(sino[-1, lefttake:righttake]) - num1 = 0 - for i in listshift: - _sino = ndimage.interpolation.shift( - _copy_sino, (0, i), prefilter=False) - factor2 = numpy.mean(_sino[0,lefttake:righttake]) - _sino = _sino * factor1 / factor2 - sinojoin = numpy.vstack((sino, _sino)) - listmetric[num1] = numpy.sum(numpy.abs(numpy.fft.fftshift( - #pyfftw.interfaces.numpy_fft.fft2( - # sinojoin[:, lefttake:righttake + 1]) - numpy.fft.fft2(sinojoin[:, lefttake:righttake + 1]) - )) * mask) - num1 = num1 + 1 - minpos = numpy.argmin(listmetric) - return init_cen + listshift[minpos] / 2.0 - - @staticmethod - def _create_mask(nrow, ncol, radius, drop): - du = 1.0 / ncol - dv = (nrow - 1.0) / (nrow * 2.0 * numpy.pi) - centerrow = numpy.ceil(nrow / 2) - 1 - centercol = numpy.ceil(ncol / 2) - 1 - # added by Edoardo Pasca - centerrow = int(centerrow) - centercol = int(centercol) - mask = numpy.zeros((nrow, ncol), dtype='float32') - for i in range(nrow): - num1 = numpy.round(((i - centerrow) * dv / radius) / du) - (p1, p2) = numpy.int16(numpy.clip(numpy.sort( - (-num1 + centercol, num1 + centercol)), 0, ncol - 1)) - mask[i, p1:p2 + 1] = numpy.ones(p2 - p1 + 1, dtype='float32') - if drop < centerrow: - mask[centerrow - drop:centerrow + drop + 1, - :] = numpy.zeros((2 * drop + 1, ncol), dtype='float32') - mask[:,centercol-1:centercol+2] = numpy.zeros((nrow, 3), dtype='float32') - return mask - - def process(self, out=None): - - projections = self.get_input() - - cor = CenterOfRotationFinder.find_center_vo(projections.as_array()) - - return cor - - -class AcquisitionDataPadder(DataProcessor): - '''Normalization based on flat and dark - - This processor read in a AcquisitionData and normalises it based on - the instrument reading with and without incident photons or neutrons. - - Input: AcquisitionData - Parameter: 2D projection with flat field (or stack) - 2D projection with dark field (or stack) - Output: AcquisitionDataSetn - ''' - - def __init__(self, - center_of_rotation = None, - acquisition_geometry = None, - pad_value = 1e-5): - kwargs = { - 'acquisition_geometry' : acquisition_geometry, - 'center_of_rotation' : center_of_rotation, - 'pad_value' : pad_value - } - - super(AcquisitionDataPadder, self).__init__(**kwargs) - - def check_input(self, dataset): - if self.acquisition_geometry is None: - self.acquisition_geometry = dataset.geometry - if dataset.number_of_dimensions == 3: - return True - else: - raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ - .format(dataset.number_of_dimensions)) - - def process(self, out=None): - projections = self.get_input() - w = projections.get_dimension_size('horizontal') - delta = w - 2 * self.center_of_rotation - - padded_width = int ( - numpy.ceil(abs(delta)) + w - ) - delta_pix = padded_width - w - - voxel_per_pixel = 1 - geom = pbalg.pb_setup_geometry_from_acquisition(projections.as_array(), - self.acquisition_geometry.angles, - self.center_of_rotation, - voxel_per_pixel ) - - padded_geometry = self.acquisition_geometry.clone() - - padded_geometry.pixel_num_h = geom['n_h'] - padded_geometry.pixel_num_v = geom['n_v'] - - delta_pix_h = padded_geometry.pixel_num_h - self.acquisition_geometry.pixel_num_h - delta_pix_v = padded_geometry.pixel_num_v - self.acquisition_geometry.pixel_num_v - - if delta_pix_h == 0: - delta_pix_h = delta_pix - padded_geometry.pixel_num_h = padded_width - #initialize a new AcquisitionData with values close to 0 - out = AcquisitionData(geometry=padded_geometry) - out = out + self.pad_value - - - #pad in the horizontal-vertical plane -> slice on angles - if delta > 0: - #pad left of middle - command = "out.array[" - for i in range(out.number_of_dimensions): - if out.dimension_labels[i] == 'horizontal': - value = '{0}:{1}'.format(delta_pix_h, delta_pix_h+w) - command = command + str(value) - else: - if out.dimension_labels[i] == 'vertical' : - value = '{0}:'.format(delta_pix_v) - command = command + str(value) - else: - command = command + ":" - if i < out.number_of_dimensions -1: - command = command + ',' - command = command + '] = projections.array' - #print (command) - else: - #pad right of middle - command = "out.array[" - for i in range(out.number_of_dimensions): - if out.dimension_labels[i] == 'horizontal': - value = '{0}:{1}'.format(0, w) - command = command + str(value) - else: - if out.dimension_labels[i] == 'vertical' : - value = '{0}:'.format(delta_pix_v) - command = command + str(value) - else: - command = command + ":" - if i < out.number_of_dimensions -1: - command = command + ',' - command = command + '] = projections.array' - #print (command) - #cleaned = eval(command) - exec(command) - return out \ No newline at end of file diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py index 78f88dd..95c0dea 100644 --- a/Wrappers/Python/setup.py +++ b/Wrappers/Python/setup.py @@ -36,7 +36,7 @@ setup( 'ccpi.optimisation.operators', 'ccpi.optimisation.algorithms', 'ccpi.optimisation.functions', - 'ccpi.optimisation.processors', + 'ccpi.processors', 'ccpi.contrib','ccpi.contrib.optimisation', 'ccpi.contrib.optimisation.algorithms'], diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index a5bd9c6..3e6a83e 100755 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -31,6 +31,9 @@ class TestDataProcessor(unittest.TestCase): cf = CenterOfRotationFinder() cf.set_input(ad) print ("Center of rotation", cf.get_output()) + self.assertAlmostEqual(86.25, cf.get_output()) + def test_Normalizer(self): + pass -- cgit v1.2.3 From 590a593b692b07af76e3e87a95e4cc01a2843aea Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 30 Apr 2019 14:23:50 +0100 Subject: change numpy versions takes into consideration #271 --- Wrappers/Python/conda-recipe/conda_build_config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/conda-recipe/conda_build_config.yaml b/Wrappers/Python/conda-recipe/conda_build_config.yaml index 30c8e9d..393ae18 100644 --- a/Wrappers/Python/conda-recipe/conda_build_config.yaml +++ b/Wrappers/Python/conda-recipe/conda_build_config.yaml @@ -4,5 +4,5 @@ python: - 3.6 numpy: # TODO investigage, as it doesn't currently build with cvxp, requires >1.14 + - 1.11 - 1.12 - - 1.15 -- 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/Python') 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 f0204684c5d9c731fcc9a49ffab2a49638a38ef1 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 1 May 2019 10:58:34 +0100 Subject: update demo --- Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py | 201 ++++++++-------------- 1 file changed, 71 insertions(+), 130 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py b/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py index 9de48a5..854f645 100644 --- a/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py +++ b/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py @@ -20,11 +20,17 @@ from ccpi.optimisation.operators import BlockOperator, Identity, Gradient from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ MixedL21Norm, BlockFunction, ScaledFunction -#from ccpi.astra.ops import AstraProjectorSimple -#from ccpi.plugins.ops import CCPiProjectorSimple from ccpi.plugins.operators import CCPiProjectorSimple -#from skimage.util import random_noise from timeit import default_timer as timer +from ccpi.reconstruction.parallelbeam import alg as pbalg +import os + +try: + import tomophantom + from tomophantom import TomoP3D + no_tomophantom = False +except ImportError as ie: + no_tomophantom = True #%% @@ -52,32 +58,11 @@ N = 75 vert = 4 ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, voxel_num_z=vert) -data = ig.allocate() -Phantom = data -# Populate image data by looping over and filling slices -i = 0 -while i < vert: - if vert > 1: - x = Phantom.subset(vertical=i).array - else: - x = Phantom.array - x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 - x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 0.98 - if vert > 1 : - Phantom.fill(x, vertical=i) - i += 1 - -#%% -#detectors = N -#angles = np.linspace(0,np.pi,100) -#angles_num = 100 -angles_num = N +angles_num = 100 det_w = 1.0 det_num = N -angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*\ - 180/np.pi angles = np.linspace(-90.,90.,N, dtype=np.float32) # Inputs: Geometry, 2D or 3D, angles, horz detector pixel count, # horz detector pixel size, vert detector pixel count, @@ -90,73 +75,59 @@ ag = AcquisitionGeometry('parallel', vert, det_w) -from ccpi.reconstruction.parallelbeam import alg as pbalg -from ccpi.plugins.processors import setupCCPiGeometries -def ssetupCCPiGeometries(ig, ag, counter): +#no_tomophantom = True +if no_tomophantom: + data = ig.allocate() + Phantom = data + # Populate image data by looping over and filling slices + i = 0 + while i < vert: + if vert > 1: + x = Phantom.subset(vertical=i).array + else: + x = Phantom.array + x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 + x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 0.98 + if vert > 1 : + Phantom.fill(x, vertical=i) + i += 1 - #vg = ImageGeometry(voxel_num_x=voxel_num_x,voxel_num_y=voxel_num_y, voxel_num_z=voxel_num_z) - #Phantom_ccpi = ImageData(geometry=vg, - # dimension_labels=['horizontal_x','horizontal_y','vertical']) - ##.subset(['horizontal_x','horizontal_y','vertical']) - ## ask the ccpi code what dimensions it would like - Phantom_ccpi = ig.allocate(dimension_labels=[ImageGeometry.HORIZONTAL_X, - ImageGeometry.HORIZONTAL_Y, - ImageGeometry.VERTICAL]) + Aop = CCPiProjectorSimple(ig, ag, 'cpu') + sin = Aop.direct(data) +else: + + model = 13 # select a model number from the library + N_size = N # Define phantom dimensions using a scalar value (cubic phantom) + path = os.path.dirname(tomophantom.__file__) + path_library3D = os.path.join(path, "Phantom3DLibrary.dat") + #This will generate a N_size x N_size x N_size phantom (3D) + phantom_tm = TomoP3D.Model(model, N_size, path_library3D) - voxel_per_pixel = 1 - angles = ag.angles - geoms = pbalg.pb_setup_geometry_from_image(Phantom_ccpi.as_array(), - angles, - voxel_per_pixel ) + #%% + Horiz_det = int(np.sqrt(2)*N_size) # detector column count (horizontal) + Vert_det = N_size # detector row count (vertical) (no reason for it to be > N) + #angles_num = int(0.5*np.pi*N_size); # angles number + #angles = np.linspace(0.0,179.9,angles_num,dtype='float32') # in degrees - pg = AcquisitionGeometry('parallel', - '3D', - angles, - geoms['n_h'], 1.0, - geoms['n_v'], 1.0 #2D in 3D is a slice 1 pixel thick - ) + print ("Building 3D analytical projection data with TomoPhantom") + projData3D_analyt = TomoP3D.ModelSino(model, + N_size, + Horiz_det, + Vert_det, + angles, + path_library3D) - center_of_rotation = Phantom_ccpi.get_dimension_size('horizontal_x') / 2 - #ad = AcquisitionData(geometry=pg,dimension_labels=['angle','vertical','horizontal']) - ad = pg.allocate(dimension_labels=[AcquisitionGeometry.ANGLE, - AcquisitionGeometry.VERTICAL, - AcquisitionGeometry.HORIZONTAL]) - geoms_i = pbalg.pb_setup_geometry_from_acquisition(ad.as_array(), - angles, - center_of_rotation, - voxel_per_pixel ) + # tomophantom outputs in [vert,angles,horiz] + # we want [angle,vert,horiz] + data = np.transpose(projData3D_analyt, [1,0,2]) + ag.pixel_num_h = Horiz_det + ag.pixel_num_v = Vert_det + sin = ag.allocate() + sin.fill(data) + ig.voxel_num_y = Vert_det - counter+=1 + Aop = CCPiProjectorSimple(ig, ag, 'cpu') - if counter < 4: - print (geoms, geoms_i) - if (not ( geoms_i == geoms )): - print ("not equal and {} {} {}".format(counter, geoms['output_volume_z'], geoms_i['output_volume_z'])) - X = max(geoms['output_volume_x'], geoms_i['output_volume_x']) - Y = max(geoms['output_volume_y'], geoms_i['output_volume_y']) - Z = max(geoms['output_volume_z'], geoms_i['output_volume_z']) - return setupCCPiGeometries(X,Y,Z,angles, counter) - else: - print ("happy now {} {} {}".format(counter, geoms['output_volume_z'], geoms_i['output_volume_z'])) - - return geoms - else: - return geoms_i - - - -#voxel_num_x, voxel_num_y, voxel_num_z, angles, counter -print ("###############################################") -print (ig) -print (ag) -g = setupCCPiGeometries(ig, ag, 0) -print (g) -print ("###############################################") -print ("###############################################") -#ag = AcquisitionGeometry('parallel','2D',angles, detectors) -#Aop = AstraProjectorSimple(ig, ag, 'cpu') -Aop = CCPiProjectorSimple(ig, ag, 'cpu') -sin = Aop.direct(data) plt.imshow(sin.subset(vertical=0).as_array()) plt.title('Sinogram') @@ -228,70 +199,40 @@ opt1 = {'niter':niter, 'memopt': True} -pdhg1 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma, memopt=True, max_iteration=niter) +pdhg1 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma, max_iteration=niter) #pdhg1.max_iteration = 2000 pdhg1.update_objective_interval = 100 -pdhg2 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma, memopt=False) -pdhg2.max_iteration = 2000 -pdhg2.update_objective_interval = 100 t1_old = timer() resold, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) t2_old = timer() -print ("memopt = False, shouldn't matter") pdhg1.run(niter) print (sum(pdhg1.timing)) res = pdhg1.get_output().subset(vertical=0) -print (pdhg1.objective) -t3 = timer() -#res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -print ("memopt = True, shouldn't matter") -pdhg2.run(niter) -print (sum(pdhg2.timing)) -res1 = pdhg2.get_output().subset(vertical=0) -t4 = timer() -# -print ("No memopt in {}s, memopt in {}/{}s old {}s".format(sum(pdhg1.timing), - sum(pdhg2.timing),t4-t3, t2_old-t1_old)) - -t1_old = timer() -resold1, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -t2_old = timer() #%% plt.figure() -plt.subplot(2,3,1) +plt.subplot(1,4,1) plt.imshow(res.as_array()) -plt.title('no memopt') -plt.colorbar() -plt.subplot(2,3,2) -plt.imshow(res1.as_array()) -plt.title('memopt') +plt.title('Algorithm') plt.colorbar() -plt.subplot(2,3,3) -plt.imshow((res1 - resold1.subset(vertical=0)).abs().as_array()) -plt.title('diff') -plt.colorbar() -plt.subplot(2,3,4) +plt.subplot(1,4,2) plt.imshow(resold.subset(vertical=0).as_array()) -plt.title('old nomemopt') +plt.title('function') plt.colorbar() -plt.subplot(2,3,5) -plt.imshow(resold1.subset(vertical=0).as_array()) -plt.title('old memopt') -plt.colorbar() -plt.subplot(2,3,6) -plt.imshow((resold1 - resold).subset(vertical=0).as_array()) -plt.title('diff old') +plt.subplot(1,4,3) +plt.imshow((res - resold.subset(vertical=0)).abs().as_array()) +plt.title('diff') plt.colorbar() -#plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt') -#plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt') -#plt.legend() +plt.subplot(1,4,4) +plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'Algorithm') +plt.plot(np.linspace(0,N,N), resold.subset(vertical=0).as_array()[int(N/2),:], label = 'function') +plt.legend() plt.show() # -print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(sum(pdhg1.timing), t4 -t3)) -diff = (res1 - res).abs().as_array().max() +print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(sum(pdhg1.timing), t2_old -t1_old)) +diff = (res - resold.subset(vertical=0)).abs().as_array().max() # print(" Max of abs difference is {}".format(diff)) -- 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/Python') 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/Python') 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/Python') 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/Python') 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 From 640feb256a172c22ef4ab78f792e214f62ef5beb Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 7 May 2019 10:10:10 +0100 Subject: allows passing update_objective_interval in creator --- Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py index c923a30..a14378c 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -52,7 +52,7 @@ class Algorithm(object): self.__loss = [] self.memopt = False self.timing = [] - self.update_objective_interval = 1 + self.update_objective_interval = kwargs.get('update_objective_interval', 1) def set_up(self, *args, **kwargs): '''Set up the algorithm''' raise NotImplementedError() -- cgit v1.2.3 From 6a8e5a212679dea33e767a2f30bac4e88bce21fe Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 7 May 2019 11:17:00 +0100 Subject: update astra calls --- Wrappers/Python/wip/pdhg_TV_tomography2D.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D.py b/Wrappers/Python/wip/pdhg_TV_tomography2D.py index cd91409..1be3dfa 100644 --- a/Wrappers/Python/wip/pdhg_TV_tomography2D.py +++ b/Wrappers/Python/wip/pdhg_TV_tomography2D.py @@ -19,8 +19,7 @@ from ccpi.optimisation.operators import BlockOperator, Gradient from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ MixedL21Norm, BlockFunction -from ccpi.astra.ops import AstraProjectorSimple -from skimage.util import random_noise +from ccpi.astra.operators import AstraProjectorSimple from timeit import default_timer as timer -- cgit v1.2.3 From d1fbb8b98862eeaddcda29ebf76e590212103ad8 Mon Sep 17 00:00:00 2001 From: Vaggelis Papoutsellis <22398586+epapoutsellis@users.noreply.github.com> Date: Tue, 7 May 2019 11:48:47 +0100 Subject: Update yaml Add matplotlib --- Wrappers/Python/conda-recipe/meta.yaml | 1 + 1 file changed, 1 insertion(+) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml index ac5f763..6564014 100644 --- a/Wrappers/Python/conda-recipe/meta.yaml +++ b/Wrappers/Python/conda-recipe/meta.yaml @@ -33,6 +33,7 @@ requirements: - python - numpy - scipy + - matplotlib - h5py about: -- cgit v1.2.3