Add regularization step

Add regularization step
OS seems to work

3 files changed, 27 insertions, 7 deletions

diff --git a/src/Python/ccpi/reconstruction/FISTAReconstructor.py b/src/Python/ccpi/reconstruction/FISTAReconstructor.py
index f43966c..c903712 100644
--- a/src/Python/ccpi/reconstruction/FISTAReconstructor.py
+++ b/src/Python/ccpi/reconstruction/FISTAReconstructor.py
@@ -363,6 +363,9 @@ class FISTAReconstructor():
             except Exception():
                 subsets = 0
             #return subsets
+        else:
+            self.setParameter(subsets=subsets)
+            
 
         angles = self.getParameter('projector_geometry')['ProjectionAngles'] 
         
@@ -371,7 +374,7 @@ class FISTAReconstructor():
         #                          subsets + 1)
         binsDiscr, binEdges = numpy.histogram(angles, bins=subsets)
         # get rearranged subset indices
-        IndicesReorg = numpy.zeros((numpy.shape(angles)))
+        IndicesReorg = numpy.zeros((numpy.shape(angles)), dtype=numpy.int32)
         counterM = 0
         for ii in range(binsDiscr.max()):
             counter = 0
diff --git a/src/Python/test/test_reconstructor-os.py b/src/Python/test/test_reconstructor-os.py
index a36feda..6c82ae0 100644
--- a/src/Python/test/test_reconstructor-os.py
+++ b/src/Python/test/test_reconstructor-os.py
@@ -12,6 +12,7 @@ import numpy
 from ccpi.reconstruction.FISTAReconstructor import FISTAReconstructor
 import astra
 import matplotlib.pyplot as plt
+from ccpi.imaging.Regularizer import Regularizer
 
 def RMSE(signal1, signal2):
     '''RMSE Root Mean Squared Error'''
@@ -77,6 +78,12 @@ fistaRecon.setParameter(ring_alpha = 21)
 fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
 
 
+reg = Regularizer(Regularizer.Algorithm.LLT_model)
+reg.setParameter(regularization_parameter=25,
+                          time_step=0.0003,
+                          tolerance_constant=0.0001,
+                          number_of_iterations=300)
+
 ## Ordered subset
 if True:
     subsets = 16
@@ -210,11 +217,12 @@ if True:
             # the number of projections per subset
             numProjSub = fistaRecon.getParameter('os_bins')[ss]
             CurrSubIndices = fistaRecon.getParameter('os_indices')\
-                             [counterInd:counterInd+numProjSub-1]
+                             [counterInd:counterInd+numProjSub]
+            #print ("Len CurrSubIndices {0}".format(numProjSub))
             mask = numpy.zeros(numpy.shape(angles), dtype=bool)
             cc = 0
-            for i in range(len(CurrSubIndices)):
-                mask[int(CurrSubIndices[i])] = True
+            for j in range(len(CurrSubIndices)):
+                mask[int(CurrSubIndices[j])] = True
             proj_geomSUB['ProjectionAngles'] = angles[mask]
 
             shape = list(numpy.shape(fistaRecon.getParameter('input_sinogram')))
@@ -254,7 +262,7 @@ if True:
     ##                    sino_updt_FULL(:,indC,:) = squeeze(sino_updt_Sub(:,kkk,:)); % filling the full sinogram
     ##                end
                  for kkk in range(numProjSub):
-                     print ("ring removal indC ... {0}".format(kkk))
+                     #print ("ring removal indC ... {0}".format(kkk))
                      indC = int(CurrSubIndices[kkk])
                      residualSub[:,kkk,:] = weights[:,indC,:].squeeze() * \
                             (sino_updt_Sub[:,kkk,:].squeeze() - \
@@ -297,7 +305,8 @@ if True:
         # regularizer = fistaRecon.getParameter('regularizer')
         # for slices:
         # out = regularizer(input=X)
-        print ("skipping regularizer")
+        print ("regularizer")
+        #X = reg(input=X)
 
 
         ## FINAL
@@ -321,7 +330,8 @@ if True:
             Resid_error[i] = RMSE(X*ROI, X_ideal*ROI)
             string = 'Iteration Number {0} | RMS Error {1} | Objective {2} \n'
             print (string.format(i,Resid_error[i], objective[i]))
-            
+
+    numpy.save("X_out_os.npy", X)
 
 else:
     fistaRecon = FISTAReconstructor(proj_geom,
diff --git a/src/Python/test/test_reconstructor.py b/src/Python/test/test_reconstructor.py
index be28624..3342301 100644
--- a/src/Python/test/test_reconstructor.py
+++ b/src/Python/test/test_reconstructor.py
@@ -76,6 +76,13 @@ fistaRecon.setParameter(Lipschitz_constant = 767893952.0)
 fistaRecon.setParameter(ring_alpha = 21)
 fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
 
+reg = Regularizer(Regularizer.Algorithm.LLT_model)
+reg.setParameter(regularization_parameter=25,
+                          time_step=0.0003,
+                          tolerance_constant=0.0001,
+                          number_of_iterations=300)
+fistaRecon.setParameter(regularizer = reg)
+
 ## Ordered subset
 if False:
     subsets = 16