Implemented using geometries in SIRT

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