Progress in pythonization

2 files changed, 298 insertions, 35 deletions

diff --git a/src/Python/ccpi/fista/FISTAReconstructor.py b/src/Python/ccpi/fista/FISTAReconstructor.py
index 8318ea6..87dd2c0 100644
--- a/src/Python/ccpi/fista/FISTAReconstructor.py
+++ b/src/Python/ccpi/fista/FISTAReconstructor.py
@@ -81,7 +81,7 @@ class FISTAReconstructor():
         self.pars['projector_geometry'] = projector_geometry # proj_geom
         self.pars['output_geometry'] = output_geometry       # vol_geom
         self.pars['input_sinogram'] = input_sinogram         # sino
-        detectors, nangles, sliceZ = numpy.shape(input_sinogram)
+        sliceZ, nangles, detectors = numpy.shape(input_sinogram)
         self.pars['detectors'] = detectors
         self.pars['number_of_angles'] = nangles
         self.pars['SlicesZ'] = sliceZ
@@ -108,7 +108,9 @@ class FISTAReconstructor():
               'regularizer' , 
               'ring_lambda_R_L1',
               'ring_alpha',
-              'subsets')
+              'subsets',
+              'use_studentt_fidelity',
+              'studentt')
         self.acceptedInputKeywords = list(kw)
         
         # handle keyworded parameters
@@ -143,16 +145,18 @@ class FISTAReconstructor():
             else:
                 self.pars['region_of_interest'] = numpy.nonzero(
                     self.pars['ideal_image']>0.0)
-            
+                
+        # the regularizer must be a correctly instantiated object    
         if not 'regularizer' in kwargs.keys() :
             self.pars['regularizer'] = None
-        else:
-            # the regularizer must be a correctly instantiated object
-            if not 'ring_lambda_R_L1' in kwargs.keys():
-                self.pars['ring_lambda_R_L1'] = 0
-            if not 'ring_alpha' in kwargs.keys():
-                self.pars['ring_alpha'] = 1
 
+        #RING REMOVAL
+        if not 'ring_lambda_R_L1' in kwargs.keys():
+            self.pars['ring_lambda_R_L1'] = 0
+        if not 'ring_alpha' in kwargs.keys():
+            self.pars['ring_alpha'] = 1
+
+        # ORDERED SUBSET
         if not 'subsets' in kwargs.keys():
             self.pars['subsets'] = 0
         else:
@@ -160,6 +164,15 @@ class FISTAReconstructor():
 
         if not 'initialize' in kwargs.keys():
             self.pars['initialize'] = False
+
+        if not 'use_studentt_fidelity' in kwargs.keys():
+            self.setParameter(studentt=False)
+        else:
+            print ("studentt {0}".format(kwargs['use_studentt_fidelity']))
+            if kwargs['use_studentt_fidelity']:
+                raise Exception('Not implemented')
+            
+            self.setParameter(studentt=kwargs['use_studentt_fidelity'])
             
             
     def setParameter(self, **kwargs):
@@ -170,6 +183,8 @@ class FISTAReconstructor():
         '''
         for key , value in kwargs.items():
             if key in self.acceptedInputKeywords:
+                if key == 'use_studentt_fidelity':
+                    raise Exception('use_studentt_fidelity Not implemented')
                 self.pars[key] = value
             else:
                 raise Exception('Wrong parameter {0} for '.format(key) +
@@ -354,10 +369,28 @@ class FISTAReconstructor():
             #return subsets
 
         angles = self.getParameter('projector_geometry')['ProjectionAngles'] 
-            
-            
-            
         
+        #binEdges = numpy.linspace(angles.min(),
+        #                          angles.max(),
+        #                          subsets + 1)
+        binsDiscr, binEdges = numpy.histogram(angles, bins=subsets)
+        # get rearranged subset indices
+        IndicesReorg = numpy.zeros((numpy.shape(angles)))
+        counterM = 0
+        for ii in range(binsDiscr.max()):
+            counter = 0
+            for jj in range(subsets):
+                curr_index = ii + jj  + counter
+                #print ("{0} {1} {2}".format(binsDiscr[jj] , ii, counterM))
+                if binsDiscr[jj] > ii:
+                    if (counterM < numpy.size(IndicesReorg)):
+                        IndicesReorg[counterM] = curr_index
+                    counterM = counterM + 1
+                    
+                counter = counter + binsDiscr[jj] - 1    
+                
+            
+        return IndicesReorg
             
 
     def prepareForIteration(self):
@@ -368,23 +401,24 @@ class FISTAReconstructor():
         detectors, nangles, sliceZ  = numpy.shape(self.pars['input_sinogram'])
         self.r = numpy.zeros((detectors, sliceZ), dtype=numpy.float)
         # another ring variable
-        self.rx = self.r.copy()
+        self.r_x = self.r.copy()
 
         self.residual = numpy.zeros(numpy.shape(self.pars['input_sinogram']))
         
         if self.getParameter('Lipschitz_constant') is None:
             self.pars['Lipschitz_constant'] = \
                             self.calculateLipschitzConstantWithPowerMethod()
+        
 
     # prepareForIteration
 
     def iterate(self, Xin=None):
         # convenience variable storage
         proj_geom , vol_geom, sino , \
-                  SlicesZ = self.getParameter(['projector_geometry' ,
-                                                        'output_geometry',
-                                                        'input_sinogram',
-                                                        'SlicesZ'])
+                  SlicesZ = self.getParameter([ 'projector_geometry' ,
+                                                'output_geometry',
+                                                'input_sinogram',
+                                                'SlicesZ'])
                         
         t = 1
         if Xin is None:    
@@ -394,7 +428,8 @@ class FISTAReconstructor():
                 N = vol_geom['GridColCount']
                 X = numpy.zeros((N,N,SlicesZ), dtype=numpy.float)
         else:
-            X = Xin.copy()
+            # copy by reference
+            X = Xin
 
         X_t = X.copy()
         
@@ -402,28 +437,31 @@ class FISTAReconstructor():
             X_old = X.copy()
             t_old = t
             r_old = self.r.copy()
-            if self.pars['projector_geometry']['type'] == 'parallel' or \
-               self.pars['projector_geometry']['type'] == 'parallel3d':
+            if self.getParameter('projector_geometry')['type'] == 'parallel' or \
+               self.getParameter('projector_geometry')['type'] == 'parallel3d':
                 # if the geometry is parallel use slice-by-slice
                 # projection-backprojection routine
                 #sino_updt = zeros(size(sino),'single');
+                proj_geomT = proj_geom.copy()
+                proj_geomT['DetectorRowCount'] = 1
+                vol_geomT = vol_geom.copy()
+                vol_geomT['GridSliceCount'] = 1;
                 sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float)
-                
-                #for kkk = 1:SlicesZ
-                #    [sino_id, sino_updt(:,:,kkk)] =
-                # astra_create_sino3d_cuda(X_t(:,:,kkk), proj_geomT, vol_geomT);
-                #    astra_mex_data3d('delete', sino_id);
                 for kkk in range(SlicesZ):
+                    print (kkk)
                     sino_id, sino_updt[kkk] = \
                              astra.creators.create_sino3d_gpu(
-                                 X_t[kkk], proj_geomT, vol_geomT)
-                    
+                                 X_t[kkk:kkk+1], proj_geomT, vol_geomT)
+                    astra.matlab.data3d('delete', sino_id)
             else:
-                # for divergent 3D geometry (watch GPU memory overflow in
-                # Astra < 1.8
-                sino_id, y = astra.creators.create_sino3d_gpu(X_t, 
-                                                              proj_geom, 
-                                                              vol_geom)
-    
+                # for divergent 3D geometry (watch the GPU memory overflow in
+                # ASTRA versions < 1.8)
+                #[sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom);
+                sino_id, sino_updt = astra.matlab.create_sino3d_gpu(
+                    X_t, proj_geom, vol_geom)
+
+
+            ## RING REMOVAL
 
+            ## REGULARIZATION
             
diff --git a/src/Python/test_reconstructor.py b/src/Python/test_reconstructor.py
index a338d34..f8f6b3c 100644
--- a/src/Python/test_reconstructor.py
+++ b/src/Python/test_reconstructor.py
@@ -31,7 +31,7 @@ recon_size = numpy.asarray(nx.get('/recon_size'), dtype=int)[0]
 size_det = numpy.asarray(nx.get('/size_det'), dtype=int)[0]
 slices_tot = numpy.asarray(nx.get('/slices_tot'), dtype=int)[0]
 
-Z_slices = 3
+Z_slices = 20
 det_row_count = Z_slices
 # next definition is just for consistency of naming
 det_col_count = size_det
@@ -64,5 +64,230 @@ fistaRecon = FISTAReconstructor(proj_geom,
                                 weights=Weights3D)
 
 print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
+fistaRecon.setParameter(number_of_iterations = 12)
+fistaRecon.setParameter(Lipschitz_constant = 767893952.0)
+fistaRecon.setParameter(ring_alpha = 21)
+fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
+#fistaRecon.setParameter(use_studentt_fidelity= True)
 
-## the calculation of the lipschitz constant should not start by itself
+## Ordered subset
+if False:
+    subsets = 16
+    angles = fistaRecon.getParameter('projector_geometry')['ProjectionAngles']
+    #binEdges = numpy.linspace(angles.min(),
+    #                          angles.max(),
+    #                          subsets + 1)
+    binsDiscr, binEdges = numpy.histogram(angles, bins=subsets)
+    # get rearranged subset indices
+    IndicesReorg = numpy.zeros((numpy.shape(angles)))
+    counterM = 0
+    for ii in range(binsDiscr.max()):
+        counter = 0
+        for jj in range(subsets):
+            curr_index = ii + jj  + counter
+            #print ("{0} {1} {2}".format(binsDiscr[jj] , ii, counterM))
+            if binsDiscr[jj] > ii:
+                if (counterM < numpy.size(IndicesReorg)):
+                    IndicesReorg[counterM] = curr_index
+                counterM = counterM + 1
+                
+            counter = counter + binsDiscr[jj] - 1
+
+
+if True:
+    fistaRecon.prepareForIteration()
+    print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
+
+    
+
+    proj_geom , vol_geom, sino , \
+                      SlicesZ = fistaRecon.getParameter(['projector_geometry' ,
+                                                            'output_geometry',
+                                                            'input_sinogram',
+                                                         'SlicesZ'])
+
+    fistaRecon.setParameter(number_of_iterations = 3)
+    iterFISTA = fistaRecon.getParameter('number_of_iterations')
+    # errors vector (if the ground truth is given)
+    Resid_error = numpy.zeros((iterFISTA));
+    # objective function values vector
+    objective = numpy.zeros((iterFISTA)); 
+
+      
+    print ("line")                        
+    t = 1
+    print ("line")
+
+    if False:
+        # if X doesn't exist
+        #N = params.vol_geom.GridColCount
+        N = vol_geom['GridColCount']
+        print ("N " + str(N))
+        X = numpy.zeros((N,N,SlicesZ), dtype=numpy.float)
+    else:
+        #X = fistaRecon.initialize()
+        X = numpy.load("X.npy")
+    
+    print (numpy.shape(X))
+    X_t = X.copy()
+    print ("X_t copy")
+##    % Outer FISTA iterations loop
+    for i in range(fistaRecon.getParameter('number_of_iterations')):
+        X_old = X.copy()
+        t_old = t
+        r_old = fistaRecon.r.copy()
+        if fistaRecon.getParameter('projector_geometry')['type'] == 'parallel' or \
+           fistaRecon.getParameter('projector_geometry')['type'] == 'parallel3d':
+            # if the geometry is parallel use slice-by-slice
+            # projection-backprojection routine
+            #sino_updt = zeros(size(sino),'single');
+            proj_geomT = proj_geom.copy()
+            proj_geomT['DetectorRowCount'] = 1
+            vol_geomT = vol_geom.copy()
+            vol_geomT['GridSliceCount'] = 1;
+            sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float)
+            for kkk in range(SlicesZ):
+                print (kkk)
+                sino_id, sino_updt[kkk] = \
+                         astra.creators.create_sino3d_gpu(
+                             X_t[kkk:kkk+1], proj_geomT, vol_geomT)
+                astra.matlab.data3d('delete', sino_id)
+        else:
+            # for divergent 3D geometry (watch the GPU memory overflow in
+            # ASTRA versions < 1.8)
+            #[sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom);
+            sino_id, sino_updt = astra.matlab.create_sino3d_gpu(
+                X_t, proj_geom, vol_geom)
+
+        ## RING REMOVAL
+        residual = fistaRecon.residual
+        lambdaR_L1 , alpha_ring , weights , L_const= \
+                   fistaRecon.getParameter(['ring_lambda_R_L1',
+                                           'ring_alpha' , 'weights',
+                                           'Lipschitz_constant'])
+        r_x = fistaRecon.r_x
+        SlicesZ, anglesNumb, Detectors = \
+                    numpy.shape(fistaRecon.getParameter('input_sinogram'))
+        if lambdaR_L1 > 0 :
+             for kkk in range(anglesNumb):
+                 print ("angles {0}".format(kkk))
+                 residual[:,kkk,:] = (weights[:,kkk,:]).squeeze() * \
+                                       ((sino_updt[:,kkk,:]).squeeze() - \
+                                        (sino[:,kkk,:]).squeeze() -\
+                                        (alpha_ring * r_x)
+                                        )
+             vec = residual.sum(axis = 1)
+             #if SlicesZ > 1:
+             #    vec = vec[:,1,:].squeeze()
+             fistaRecon.r = (r_x - (1./L_const) * vec).copy()
+             objective[i] = (0.5 * (residual ** 2).sum())
+##            % the ring removal part (Group-Huber fidelity)
+##            for kkk = 1:anglesNumb
+##                residual(:,kkk,:) =  squeeze(weights(:,kkk,:)).*
+##                 (squeeze(sino_updt(:,kkk,:)) -
+##                 (squeeze(sino(:,kkk,:)) - alpha_ring.*r_x));
+##            end
+##            vec = sum(residual,2);
+##            if (SlicesZ > 1)
+##                vec = squeeze(vec(:,1,:));
+##            end
+##            r = r_x - (1./L_const).*vec;
+##            objective(i) = (0.5*sum(residual(:).^2)); % for the objective function output
+
+        else:
+            if fistaRecon.getParameter('use_studentt_fidelity'):
+                residual = weights * (sino_updt - sino)
+                for kkk in range(SlicesZ):
+                    # reshape(residual(:,:,kkk), Detectors*anglesNumb, 1)
+                    # 1D
+                    res_vec = numpy.reshape(residual[kkk], (Detectors * anglesNumb,1))
+                
+##            else            
+##			if (studentt == 1)
+##				% artifacts removal with Students t penalty
+##				residual = weights.*(sino_updt - sino);
+##				for kkk = 1:SlicesZ
+##				   res_vec = reshape(residual(:,:,kkk), Detectors*anglesNumb, 1); % 1D vectorized sinogram
+##                                 %s = 100;
+##                                 %gr = (2)*res_vec./(s*2 + conj(res_vec).*res_vec);
+##                                 [ff, gr] = studentst(res_vec, 1);
+##                                 residual(:,:,kkk) = reshape(gr, Detectors, anglesNumb);
+##				end            
+##				objective(i) = ff; % for the objective function output
+##            else
+##            % no ring removal (LS model)
+##            residual = weights.*(sino_updt - sino);
+##            objective(i) = (0.5*sum(residual(:).^2)); % for the objective function output
+##            end	
+##        end       
+
+        # Projection/Backprojection Routine
+        if fistaRecon.getParameter('projector_geometry')['type'] == 'parallel' or \
+           fistaRecon.getParameter('projector_geometry')['type'] == 'parallel3d':
+            x_temp = numpy.zeros(numpy.shape(X),dtype=numpy.float32)
+            for kkk in range(SlicesZ):
+                print ("Projection/Backprojection Routine {0}".format( kkk ))
+                x_id, x_temp[kkk] = \
+                         astra.creators.create_backprojection3d_gpu(
+                             residual[kkk:kkk+1],
+                             proj_geomT, vol_geomT)
+                astra.matlab.data3d('delete', x_id)
+        else:
+            x_id, x_temp = \
+                  astra.creators.create_backprojection3d_gpu(
+                      residual, proj_geom, vol_geom)            
+
+        X = X_t - (1/L_const) * x_temp
+        astra.matlab.data3d('delete', sino_id)
+        astra.matlab.data3d('delete', x_id)
+        
+
+        ## REGULARIZATION
+        ## SKIPPING FOR NOW
+        ## Should be simpli
+        # regularizer = fistaRecon.getParameter('regularizer')
+        # for slices:
+        # out = regularizer(input=X)
+
+
+        ## FINAL
+        lambdaR_L1 = fistaRecon.getParameter('ring_lambda_R_L1')
+        if lambdaR_L1 > 0:
+            fistaRecon.r = numpy.max(
+                numpy.abs(fistaRecon.r) - lambdaR_L1 , 0) * \
+                numpy.sign(fistaRecon.r)
+        t = (1 + numpy.sqrt(1 + 4 * t**2))/2
+        X_t = X + (((t_old -1)/t) * (X - X_old))
+
+        if lambdaR_L1 > 0:
+            fistaRecon.r_x = fistaRecon.r + \
+                             (((t_old-1)/t) * (fistaRecon.r - r_old))
+
+        if fistaRecon.getParameter('ideal_image') is None:
+            string = 'Iteration Number {0} | Objective {1} \n'
+            print (string.format( i, objective[i]))
+            
+##        if (lambdaR_L1 > 0)
+##            r =  max(abs(r)-lambdaR_L1, 0).*sign(r); % soft-thresholding operator for ring vector
+##        end
+##        
+##        t = (1 + sqrt(1 + 4*t^2))/2; % updating t
+##        X_t = X + ((t_old-1)/t).*(X - X_old); % updating X
+##        
+##        if (lambdaR_L1 > 0)
+##            r_x = r + ((t_old-1)/t).*(r - r_old); % updating r
+##        end
+##        
+##        if (show == 1)
+##            figure(10); imshow(X(:,:,slice), [0 maxvalplot]);
+##            if (lambdaR_L1 > 0)
+##                figure(11); plot(r); title('Rings offset vector')
+##            end
+##            pause(0.01);
+##        end
+##        if (strcmp(X_ideal, 'none' ) == 0)
+##            Resid_error(i) = RMSE(X(ROI), X_ideal(ROI));
+##            fprintf('%s %i %s %s %.4f  %s %s %f \n', 'Iteration Number:', i, '|', 'Error RMSE:', Resid_error(i), '|', 'Objective:', objective(i));
+##        else
+##            fprintf('%s %i  %s %s %f \n', 'Iteration Number:', i, '|', 'Objective:', objective(i));
+##        end