diff options
| -rw-r--r-- | Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py | 24 | ||||
| -rwxr-xr-x | Wrappers/Python/ccpi/optimisation/algs.py | 10 | ||||
| -rw-r--r-- | Wrappers/Python/wip/demo_test_sirt.py | 46 |
3 files changed, 16 insertions, 64 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py index cb8b731..beba913 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py @@ -24,11 +24,6 @@ Created on Wed Apr 10 13:39:35 2019 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. -""" -Created on Thu Feb 21 11:11:23 2019 - -@author: ofn77899 -""" from ccpi.optimisation.algorithms import Algorithm from ccpi.framework import ImageData, AcquisitionData @@ -43,20 +38,21 @@ class SIRT(Algorithm): operator: operator for forward/backward projections data: data to operate on constraint: Function with prox-method, for example IndicatorBox to - enforce box constraints. + enforce box constraints, default is None). ''' def __init__(self, **kwargs): super(SIRT, self).__init__() self.x = kwargs.get('x_init', None) self.operator = kwargs.get('operator', None) self.data = kwargs.get('data', None) - self.constraint = kwargs.get('data', None) + self.constraint = kwargs.get('constraint', None) if self.x is not None and self.operator is not None and \ self.data is not None: print ("Calling from creator") - self.set_up(x_init =kwargs['x_init'], - operator=kwargs['operator'], - data =kwargs['data']) + self.set_up(x_init=kwargs['x_init'], + operator=kwargs['operator'], + data=kwargs['data'], + constraint=kwargs['constraint']) def set_up(self, x_init, operator , data, constraint=None ): @@ -70,12 +66,8 @@ class SIRT(Algorithm): self.relax_par = 1.0 # Set up scaling matrices D and M. - im1 = ImageData(geometry=self.x.geometry) - im1.array[:] = 1.0 - self.M = 1/operator.direct(im1) - aq1 = AcquisitionData(geometry=self.M.geometry) - aq1.array[:] = 1.0 - self.D = 1/operator.adjoint(aq1) + self.M = 1/self.operator.direct(self.operator.domain_geometry().allocate(value=1.0)) + self.D = 1/self.operator.adjoint(self.operator.range_geometry().allocate(value=1.0)) def update(self): diff --git a/Wrappers/Python/ccpi/optimisation/algs.py b/Wrappers/Python/ccpi/optimisation/algs.py index 5a95c5c..89b5519 100755 --- a/Wrappers/Python/ccpi/optimisation/algs.py +++ b/Wrappers/Python/ccpi/optimisation/algs.py @@ -289,14 +289,8 @@ def SIRT(x_init, operator , data , opt=None, constraint=None): relax_par = 1.0 # Set up scaling matrices D and M. - im1 = ImageData(geometry=x_init.geometry) - im1.array[:] = 1.0 - M = 1/operator.direct(im1) - del im1 - aq1 = AcquisitionData(geometry=M.geometry) - aq1.array[:] = 1.0 - D = 1/operator.adjoint(aq1) - del aq1 + M = 1/operator.direct(operator.domain_geometry().allocate(value=1.0)) + D = 1/operator.adjoint(operator.range_geometry().allocate(value=1.0)) # algorithm loop for it in range(0, max_iter): diff --git a/Wrappers/Python/wip/demo_test_sirt.py b/Wrappers/Python/wip/demo_test_sirt.py index d3f27cf..8f65f39 100644 --- a/Wrappers/Python/wip/demo_test_sirt.py +++ b/Wrappers/Python/wip/demo_test_sirt.py @@ -11,6 +11,8 @@ from ccpi.astra.ops import AstraProjectorSimple from ccpi.optimisation.algorithms import CGLS as CGLSALG from ccpi.optimisation.algorithms import SIRT as SIRTALG +from ccpi.optimisation.operators import Identity + import numpy as np import matplotlib.pyplot as plt @@ -68,6 +70,8 @@ else: # wrapping calls to ASTRA as well as specifying whether to use CPU or GPU. Aop = AstraProjectorSimple(ig, ag, 'gpu') +Aop = Identity(ig,ig) + # Forward and backprojection are available as methods direct and adjoint. Here # generate test data b and do simple backprojection to obtain z. b = Aop.direct(Phantom) @@ -89,6 +93,7 @@ plt.show() x_init = ImageData(np.zeros(x.shape),geometry=ig) opt = {'tol': 1e-4, 'iter': 100} + # First a CGLS reconstruction can be done: x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt) @@ -132,7 +137,6 @@ plt.title('SIRT unconstrained criterion') plt.show() - my_SIRT_alg = SIRTALG() my_SIRT_alg.set_up(x_init, Aop, b ) my_SIRT_alg.max_iteration = 2000 @@ -145,6 +149,7 @@ plt.title('SIRT ALG') plt.colorbar() plt.show() + # A SIRT nonnegativity constrained reconstruction can be done using the # additional input "constraint" set to a box indicator function with 0 as the # lower bound and the default upper bound of infinity: @@ -201,42 +206,3 @@ plt.imshow(x_SIRT_alg01.array) plt.title('SIRT ALG01') plt.colorbar() plt.show() - -# The indicator function can also be used with the FISTA algorithm to do -# least squares with nonnegativity constraint. - -''' -# Create least squares object instance with projector, test data and a constant -# coefficient of 0.5: -f = Norm2sq(Aop,b,c=0.5) -# Run FISTA for least squares without constraints -x_fista, it, timing, criter = FISTA(x_init, f, None,opt) -plt.figure() -plt.imshow(x_fista.array) -plt.title('FISTA Least squares') -plt.show() -plt.figure() -plt.semilogy(criter) -plt.title('FISTA Least squares criterion') -plt.show() -# Run FISTA for least squares with nonnegativity constraint -x_fista0, it0, timing0, criter0 = FISTA(x_init, f, IndicatorBox(lower=0),opt) -plt.figure() -plt.imshow(x_fista0.array) -plt.title('FISTA Least squares nonneg') -plt.show() -plt.figure() -plt.semilogy(criter0) -plt.title('FISTA Least squares nonneg criterion') -plt.show() -# Run FISTA for least squares with box constraint [0,1] -x_fista01, it01, timing01, criter01 = FISTA(x_init, f, IndicatorBox(lower=0,upper=1),opt) -plt.figure() -plt.imshow(x_fista01.array) -plt.title('FISTA Least squares box(0,1)') -plt.show() -plt.figure() -plt.semilogy(criter01) -plt.title('FISTA Least squares box(0,1) criterion') -plt.show() -'''
\ No newline at end of file |