cleaning stage, leaving only FGP/ROF and related compilers, demos

7 files changed, 20 insertions, 1040 deletions

diff --git a/Wrappers/Python/test/astra_test.py b/Wrappers/Python/test/astra_test.py
deleted file mode 100644
index 42c375a..0000000
--- a/Wrappers/Python/test/astra_test.py
+++ /dev/null
@@ -1,85 +0,0 @@
-import astra
-import numpy
-import filefun
-
-
-# read in the same data as the DemoRD2
-angles = filefun.dlmread("DemoRD2/angles.csv")
-darks_ar = filefun.dlmread("DemoRD2/darks_ar.csv", separator=",")
-flats_ar = filefun.dlmread("DemoRD2/flats_ar.csv", separator=",")
-
-if True:
-    Sino3D = numpy.load("DemoRD2/Sino3D.npy")
-else:
-    sino = filefun.dlmread("DemoRD2/sino_01.csv", separator=",")
-    a = map (lambda x:x, numpy.shape(sino))
-    a.append(20)
-
-    Sino3D = numpy.zeros(tuple(a), dtype="float")
-
-    for i in range(1,numpy.shape(Sino3D)[2]+1):
-        print("Read file DemoRD2/sino_%02d.csv" % i)
-        sino = filefun.dlmread("DemoRD2/sino_%02d.csv" % i, separator=",")
-        Sino3D.T[i-1] = sino.T
-    
-Weights3D = numpy.asarray(Sino3D, dtype="float")
-
-##angles_rad = angles*(pi/180); % conversion to radians
-##size_det = size(data_raw3D,1); % detectors dim
-##angSize = size(data_raw3D, 2); % angles dim
-##slices_tot = size(data_raw3D, 3); % no of slices
-##recon_size = 950; % reconstruction size
-
-
-angles_rad = angles * numpy.pi /180.
-size_det, angSize, slices_tot = numpy.shape(Sino3D)
-size_det, angSize, slices_tot = [int(i) for i in numpy.shape(Sino3D)]
-recon_size = 950
-Z_slices = 3;
-det_row_count = Z_slices;
-
-#proj_geom = astra_create_proj_geom('parallel3d', 1, 1,
-# det_row_count, size_det, angles_rad);
-
-detectorSpacingX = 1.0
-detectorSpacingY = detectorSpacingX
-proj_geom = astra.create_proj_geom('parallel3d',
-                                            detectorSpacingX,
-                                            detectorSpacingY,
-                                            det_row_count,
-                                            size_det,
-                                            angles_rad)
-
-#vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices);
-vol_geom = astra.create_vol_geom(recon_size,recon_size,Z_slices);
-
-sino = numpy.zeros((size_det, angSize, slices_tot), dtype="float")
-
-#weights = ones(size(sino));
-weights = numpy.ones(numpy.shape(sino))
-
-#####################################################################
-## PowerMethod for Lipschitz constant
-
-N = vol_geom['GridColCount']
-x1 = numpy.random.rand(1,N,N)
-#sqweight = sqrt(weights(:,:,1));
-sqweight = numpy.sqrt(weights.T[0]).T
-##proj_geomT = proj_geom;
-proj_geomT = proj_geom.copy()
-##proj_geomT.DetectorRowCount = 1;
-proj_geomT['DetectorRowCount'] = 1
-##vol_geomT = vol_geom;
-vol_geomT = vol_geom.copy()
-##vol_geomT.GridSliceCount = 1;
-vol_geomT['GridSliceCount'] = 1
-
-##[sino_id, y] = astra_create_sino3d_cuda(x1, proj_geomT, vol_geomT);
-
-#sino_id, y = astra.create_sino3d_gpu(x1, proj_geomT, vol_geomT);
-sino_id, y = astra.create_sino(x1, proj_geomT, vol_geomT);
-
-##y = sqweight.*y;
-##astra_mex_data3d('delete', sino_id);
-        
-
diff --git a/Wrappers/Python/test/create_phantom_projections.py b/Wrappers/Python/test/create_phantom_projections.py
deleted file mode 100644
index 20a9278..0000000
--- a/Wrappers/Python/test/create_phantom_projections.py
+++ /dev/null
@@ -1,49 +0,0 @@
-from ccpi.reconstruction.AstraDevice import AstraDevice
-from ccpi.reconstruction.DeviceModel import DeviceModel
-import h5py
-import numpy
-import matplotlib.pyplot as plt
-
-nx = h5py.File('phant3D_256.h5', "r")
-phantom = numpy.asarray(nx.get('/dataset1'))
-pX,pY,pZ = numpy.shape(phantom)
-
-filename = r'/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/demos/DendrData.h5'
-nxa = h5py.File(filename, "r")
-#getEntry(nx, '/')
-# I have exported the entries as children of /
-entries = [entry for entry in nxa['/'].keys()]
-print (entries)
-
-angles_rad = numpy.asarray(nxa.get('/angles_rad'), dtype="float32")
-
-
-device = AstraDevice(
-    DeviceModel.DeviceType.PARALLEL3D.value,
-    [ pX , pY , 1., 1., angles_rad],
-    [ pX, pY, pZ ] )
-
-
-proj = device.doForwardProject(phantom)
-stack = [proj[:,i,:] for i in range(len(angles_rad))]
-stack = numpy.asarray(stack)
-
-
-fig = plt.figure()
-a=fig.add_subplot(1,2,1)
-a.set_title('proj')
-imgplot = plt.imshow(proj[:,100,:])
-a=fig.add_subplot(1,2,2)
-a.set_title('stack')
-imgplot = plt.imshow(stack[100])
-plt.show()
-
-pf = h5py.File("phantom3D256_projections.h5" , "w")
-pf.create_dataset("/projections", data=stack)
-pf.create_dataset("/sinogram", data=proj)
-pf.create_dataset("/angles", data=angles_rad)
-pf.create_dataset("/reconstruction_volume" , data=numpy.asarray([pX, pY, pZ]))
-pf.create_dataset("/camera/size" , data=numpy.asarray([pX , pY ]))
-pf.create_dataset("/camera/spacing" , data=numpy.asarray([1.,1.]))
-pf.flush()
-pf.close()
diff --git a/Wrappers/Python/test/metrics.py b/Wrappers/Python/test/metrics.py
new file mode 100644
index 0000000..53f68fb
--- /dev/null
+++ b/Wrappers/Python/test/metrics.py
@@ -0,0 +1,20 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Feb 21 13:34:32 2018
+# quality metrics
+@author: algol
+"""
+import numpy as np
+
+def nrmse(im1, im2):
+    a, b = im1.shape
+    rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b))
+    max_val = max(np.max(im1), np.max(im2))
+    min_val = min(np.min(im1), np.min(im2))
+    return 1 - (rmse / (max_val - min_val))
+    
+def rmse(im1, im2):
+    a, b = im1.shape
+    rmse = np.sqrt(np.sum((im1 - im2) ** 2) / float(a * b))    
+    return rmse
\ No newline at end of file
diff --git a/Wrappers/Python/test/readhd5.py b/Wrappers/Python/test/readhd5.py
deleted file mode 100644
index eff6c43..0000000
--- a/Wrappers/Python/test/readhd5.py
+++ /dev/null
@@ -1,42 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Wed Aug 23 16:34:49 2017
-
-@author: ofn77899
-"""
-
-import h5py
-import numpy
-
-def getEntry(nx, location):
-    for item in nx[location].keys():
-        print (item)
-        
-filename = r'/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/demos/DendrData.h5'
-nx = h5py.File(filename, "r")
-#getEntry(nx, '/')
-# I have exported the entries as children of /
-entries = [entry for entry in nx['/'].keys()]
-print (entries)
-
-Sino3D = numpy.asarray(nx.get('/Sino3D'))
-Weights3D = numpy.asarray(nx.get('/Weights3D'))
-angSize = numpy.asarray(nx.get('/angSize'), dtype=int)[0]
-angles_rad = numpy.asarray(nx.get('/angles_rad'))
-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]
-
-#from ccpi.viewer.CILViewer2D import CILViewer2D
-#v = CILViewer2D()
-#v.setInputAsNumpy(Weights3D)
-#v.startRenderLoop()
-
-import matplotlib.pyplot as plt
-fig = plt.figure()
-
-a=fig.add_subplot(1,1,1)
-a.set_title('noise')
-imgplot = plt.imshow(Weights3D[0].T)
-plt.show()
diff --git a/Wrappers/Python/test/simple_astra_test.py b/Wrappers/Python/test/simple_astra_test.py
deleted file mode 100644
index 905eeea..0000000
--- a/Wrappers/Python/test/simple_astra_test.py
+++ /dev/null
@@ -1,25 +0,0 @@
-import astra
-import numpy
-
-detectorSpacingX = 1.0
-detectorSpacingY = 1.0
-det_row_count = 128
-det_col_count = 128
-
-angles_rad = numpy.asarray([i for i in range(360)], dtype=float) / 180. * numpy.pi
-
-proj_geom = astra.creators.create_proj_geom('parallel3d',
-                                            detectorSpacingX,
-                                            detectorSpacingY,
-                                            det_row_count,
-                                            det_col_count,
-                                            angles_rad)
-
-image_size_x = 64
-image_size_y = 64
-image_size_z = 32
-
-vol_geom = astra.creators.create_vol_geom(image_size_x,image_size_y,image_size_z)
-
-x1 = numpy.random.rand(image_size_z,image_size_y,image_size_x)
-sino_id, y = astra.creators.create_sino3d_gpu(x1, proj_geom, vol_geom)
diff --git a/Wrappers/Python/test/test_reconstructor-os_phantom.py b/Wrappers/Python/test/test_reconstructor-os_phantom.py
deleted file mode 100644
index 01f1354..0000000
--- a/Wrappers/Python/test/test_reconstructor-os_phantom.py
+++ /dev/null
@@ -1,480 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Wed Aug 23 16:34:49 2017
-
-@author: ofn77899
-Based on DemoRD2.m
-"""
-
-import h5py
-import numpy
-
-from ccpi.reconstruction.FISTAReconstructor import FISTAReconstructor
-import astra
-import matplotlib.pyplot as plt
-from ccpi.imaging.Regularizer import Regularizer
-from ccpi.reconstruction.AstraDevice import AstraDevice
-from ccpi.reconstruction.DeviceModel import DeviceModel
-
-#from ccpi.viewer.CILViewer2D import *
-
-
-def RMSE(signal1, signal2):
-    '''RMSE Root Mean Squared Error'''
-    if numpy.shape(signal1) == numpy.shape(signal2):
-        err = (signal1 - signal2)
-        err = numpy.sum( err * err )/numpy.size(signal1);  # MSE
-        err = sqrt(err);                                   # RMSE
-        return err
-    else:
-        raise Exception('Input signals must have the same shape')
-
-filename = r'/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/src/Python/test/phantom3D256_projections.h5'
-nx = h5py.File(filename, "r")
-#getEntry(nx, '/')
-# I have exported the entries as children of /
-entries = [entry for entry in nx['/'].keys()]
-print (entries)
-
-projections = numpy.asarray(nx.get('/projections'), dtype="float32")
-#Weights3D = numpy.asarray(nx.get('/Weights3D'), dtype="float32")
-#angSize = numpy.asarray(nx.get('/angSize'), dtype=int)[0]
-angles_rad = numpy.asarray(nx.get('/angles'), dtype="float32")
-angSize = numpy.size(angles_rad)
-image_size_x, image_size_y, image_size_z = \
-              numpy.asarray(nx.get('/reconstruction_volume'), dtype=int)
-det_col_count, det_row_count = \
-              numpy.asarray(nx.get('/camera/size'))
-#slices_tot = numpy.asarray(nx.get('/slices_tot'), dtype=int)[0]
-detectorSpacingX, detectorSpacingY = numpy.asarray(nx.get('/camera/spacing'), dtype=int)
-
-Z_slices = 20
-#det_row_count = image_size_y
-# next definition is just for consistency of naming
-#det_col_count = image_size_x
-
-detectorSpacingX = 1.0
-detectorSpacingY = detectorSpacingX
-
-
-proj_geom = astra.creators.create_proj_geom('parallel3d',
-                                            detectorSpacingX,
-                                            detectorSpacingY,
-                                            det_row_count,
-                                            det_col_count,
-                                            angles_rad)
-
-#vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices);
-##image_size_x = recon_size
-##image_size_y = recon_size
-##image_size_z = Z_slices
-vol_geom = astra.creators.create_vol_geom( image_size_x,
-                                           image_size_y,
-                                           image_size_z)
-
-## First pass the arguments to the FISTAReconstructor and test the
-## Lipschitz constant
-astradevice = AstraDevice(DeviceModel.DeviceType.PARALLEL3D.value,
-                [proj_geom['DetectorRowCount'] ,
-                 proj_geom['DetectorColCount'] ,
-                 proj_geom['DetectorSpacingX'] ,
-                 proj_geom['DetectorSpacingY'] ,
-                 proj_geom['ProjectionAngles']
-                 ],
-                [
-                    vol_geom['GridColCount'],
-                    vol_geom['GridRowCount'], 
-                    vol_geom['GridSliceCount'] ] )
-## create the sinogram 
-Sino3D = numpy.transpose(projections, axes=[1,0,2])
-
-fistaRecon = FISTAReconstructor(proj_geom,
-                                vol_geom,
-                                Sino3D ,
-                                #weights=Weights3D,
-                                device=astradevice)
-
-print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
-fistaRecon.setParameter(number_of_iterations = 4)
-#fistaRecon.setParameter(Lipschitz_constant = 767893952.0)
-fistaRecon.setParameter(ring_alpha = 21)
-fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
-#fistaRecon.setParameter(ring_lambda_R_L1 = 0)
-subsets = 8
-fistaRecon.setParameter(subsets=subsets)
-
-
-#reg = Regularizer(Regularizer.Algorithm.FGP_TV)
-#reg.setParameter(regularization_parameter=0.005,
-#                          number_of_iterations=50)
-reg = Regularizer(Regularizer.Algorithm.FGP_TV)
-reg.setParameter(regularization_parameter=5e6,
-                          tolerance_constant=0.0001,
-                          number_of_iterations=50)
-
-#fistaRecon.setParameter(regularizer=reg)
-#lc = fistaRecon.getParameter('Lipschitz_constant')
-#reg.setParameter(regularization_parameter=5e6/lc)
-
-## Ordered subset
-if True:
-    #subsets = 8
-    fistaRecon.setParameter(subsets=subsets)
-    fistaRecon.createOrderedSubsets()
-else:
-    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:
-    print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
-    print ("prepare for iteration")
-    fistaRecon.prepareForIteration()
-    
-    
-
-    print("initializing  ...")
-    if True:
-        # if X doesn't exist
-        #N = params.vol_geom.GridColCount
-        N = vol_geom['GridColCount']
-        print ("N " + str(N))
-        X = numpy.asarray(numpy.ones((image_size_x,image_size_y,image_size_z)),
-                        dtype=numpy.float) * 0.001
-        X = numpy.asarray(numpy.zeros((image_size_x,image_size_y,image_size_z)),
-                        dtype=numpy.float) 
-    else:
-        #X = fistaRecon.initialize()
-        X = numpy.load("X.npy")
-    
-    print (numpy.shape(X))
-    X_t = X.copy()
-    print ("initialized")
-    proj_geom , vol_geom, sino , \
-        SlicesZ, weights , alpha_ring = fistaRecon.getParameter(
-            ['projector_geometry' , 'output_geometry',
-             'input_sinogram', 'SlicesZ' ,  'weights', 'ring_alpha'])
-    lambdaR_L1 , alpha_ring , weights , L_const= \
-                       fistaRecon.getParameter(['ring_lambda_R_L1',
-                                               'ring_alpha' , 'weights',
-                                               'Lipschitz_constant'])
-
-    #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)); 
-
-      
-    t = 1
-    
-
-    ## additional for 
-    proj_geomSUB = proj_geom.copy()
-    fistaRecon.residual2 = numpy.zeros(numpy.shape(fistaRecon.pars['input_sinogram']))
-    residual2 = fistaRecon.residual2
-    sino_updt_FULL = fistaRecon.residual.copy()
-    r_x = fistaRecon.r.copy()
-
-    results = []
-    print ("starting iterations")
-##    % Outer FISTA iterations loop
-    for i in range(fistaRecon.getParameter('number_of_iterations')):
-##        % With OS approach it becomes trickier to correlate independent subsets, hence additional work is required
-##        % one solution is to work with a full sinogram at times
-##        if ((i >= 3) && (lambdaR_L1 > 0))
-##            [sino_id2, sino_updt2] = astra_create_sino3d_cuda(X, proj_geom, vol_geom);
-##            astra_mex_data3d('delete', sino_id2);
-##        end
-        # With OS approach it becomes trickier to correlate independent subsets,
-        # hence additional work is required one solution is to work with a full
-        # sinogram at times
-
-        
-        #t_old = t
-        SlicesZ, anglesNumb, Detectors = \
-                    numpy.shape(fistaRecon.getParameter('input_sinogram'))
-        ## https://github.com/vais-ral/CCPi-FISTA_Reconstruction/issues/4
-        r_old = fistaRecon.r.copy()
-            
-        if (i > 1 and lambdaR_L1 > 0) :
-            for kkk in range(anglesNumb):
-                 
-                 residual2[:,kkk,:] = (weights[:,kkk,:]).squeeze() * \
-                                       ((sino_updt_FULL[:,kkk,:]).squeeze() - \
-                                        (sino[:,kkk,:]).squeeze() -\
-                                        (alpha_ring * r_x)
-                                        )
-            #r_old = fistaRecon.r.copy()
-            vec = fistaRecon.residual.sum(axis = 1)
-            #if SlicesZ > 1:
-            #    vec = vec[:,1,:] # 1 or 0?
-            r_x = fistaRecon.r_x
-            # update ring variable
-            fistaRecon.r = (r_x - (1./L_const) * vec)
-
-        # subset loop
-        counterInd = 1
-        geometry_type = fistaRecon.getParameter('projector_geometry')['type']
-        angles = fistaRecon.getParameter('projector_geometry')['ProjectionAngles']
-    
-##        if geometry_type == 'parallel' or \
-##           geometry_type == 'fanflat' or \
-##           geometry_type == 'fanflat_vec' :
-##            
-##            for kkk in range(SlicesZ):
-##                sino_id, sinoT[kkk] = \
-##                         astra.creators.create_sino3d_gpu(
-##                             X_t[kkk:kkk+1], proj_geomSUB, vol_geom)
-##                sino_updt_Sub[kkk] = sinoT.T.copy()
-##                
-##        else:
-##            sino_id, sino_updt_Sub = \
-##                astra.creators.create_sino3d_gpu(X_t, proj_geomSUB, vol_geom)
-##        
-##        astra.matlab.data3d('delete', sino_id)
-  
-        for ss in range(fistaRecon.getParameter('subsets')):
-            print ("Subset {0}".format(ss))
-            X_old = X.copy()
-            t_old = t
-            print ("X[0][0][0] {0} t {1}".format(X[0][0][0], t))
-            
-            # the number of projections per subset
-            numProjSub = fistaRecon.getParameter('os_bins')[ss]
-            CurrSubIndices = fistaRecon.getParameter('os_indices')\
-                             [counterInd:counterInd+numProjSub]
-            shape = list(numpy.shape(fistaRecon.getParameter('input_sinogram')))
-            shape[1] = numProjSub
-            sino_updt_Sub = numpy.zeros(shape)
-            
-            #print ("Len CurrSubIndices {0}".format(numProjSub))
-            mask = numpy.zeros(numpy.shape(angles), dtype=bool)
-            cc = 0
-            for j in range(len(CurrSubIndices)):
-                mask[int(CurrSubIndices[j])] = True
-
-            ## this is a reduced device
-            rdev = fistaRecon.getParameter('device_model')\
-                   .createReducedDevice(proj_par={'angles' : angles[mask]},
-                                        vol_par={})
-            proj_geomSUB['ProjectionAngles'] = angles[mask]
-
-            
-
-            if geometry_type == 'parallel' or \
-               geometry_type == 'fanflat' or \
-               geometry_type == 'fanflat_vec' :
-
-                for kkk in range(SlicesZ):
-                    sino_id, sinoT = astra.creators.create_sino3d_gpu (
-                        X_t[kkk:kkk+1] , proj_geomSUB, vol_geom)
-                    sino_updt_Sub[kkk] = sinoT.T.copy()
-                    astra.matlab.data3d('delete', sino_id)
-            else:
-                # for 3D geometry (watch the GPU memory overflow in ASTRA < 1.8)
-                sino_id, sino_updt_Sub = \
-                     astra.creators.create_sino3d_gpu (X_t,
-                                                       proj_geomSUB,
-                                                       vol_geom)
-                
-                astra.matlab.data3d('delete', sino_id)
-                
-            
-
-
-            ## RING REMOVAL
-            residual = fistaRecon.residual
-            
-            
-            if lambdaR_L1 > 0 :
-                 print ("ring removal")
-                 residualSub = numpy.zeros(shape)
-    ##             for a chosen subset
-    ##                for kkk = 1:numProjSub
-    ##                    indC = CurrSubIndeces(kkk);
-    ##                    residualSub(:,kkk,:) =  squeeze(weights(:,indC,:)).*(squeeze(sino_updt_Sub(:,kkk,:)) - (squeeze(sino(:,indC,:)) - alpha_ring.*r_x));
-    ##                    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))
-                     indC = int(CurrSubIndices[kkk])
-                     residualSub[:,kkk,:] = weights[:,indC,:].squeeze() * \
-                            (sino_updt_Sub[:,kkk,:].squeeze() - \
-                             sino[:,indC,:].squeeze() - alpha_ring * r_x)
-                     # filling the full sinogram
-                     sino_updt_FULL[:,indC,:] = sino_updt_Sub[:,kkk,:].squeeze()
-
-            else:
-                #PWLS model
-                # I guess we need to use mask here instead
-                residualSub = weights[:,CurrSubIndices,:] * \
-                              ( sino_updt_Sub - \
-                                sino[:,CurrSubIndices,:].squeeze() )
-            # it seems that in the original code the following like is not
-            # calculated in the case of ring removal
-            objective[i] = 0.5 * numpy.linalg.norm(residualSub)
-
-            #backprojection
-            if geometry_type == 'parallel' or \
-               geometry_type == 'fanflat' or \
-               geometry_type == 'fanflat_vec' :
-                # if geometry is 2D use slice-by-slice projection-backprojection
-                # routine
-                x_temp = numpy.zeros(numpy.shape(X), dtype=numpy.float32)
-                for kkk in range(SlicesZ):
-                    
-                    x_id, x_temp[kkk] = \
-                             astra.creators.create_backprojection3d_gpu(
-                                 residualSub[kkk:kkk+1],
-                                 proj_geomSUB, vol_geom)
-                    astra.matlab.data3d('delete', x_id)
-                    
-            else:
-                x_id, x_temp = \
-                      astra.creators.create_backprojection3d_gpu(
-                          residualSub, proj_geomSUB, vol_geom)
-
-                astra.matlab.data3d('delete', x_id)
-                
-            X = X_t - (1/L_const) * x_temp
-
-        
-
-            ## REGULARIZATION
-            ## SKIPPING FOR NOW
-            ## Should be simpli
-            # regularizer = fistaRecon.getParameter('regularizer')
-            # for slices:
-            # out = regularizer(input=X)
-            print ("regularizer")
-            reg = fistaRecon.getParameter('regularizer')
-
-            if reg is not None:
-                X = reg(input=X,
-                        output_all=False)
-
-            t = (1 + numpy.sqrt(1 + 4 * t **2))/2
-            X_t = X + (((t_old -1)/t) * (X-X_old))
-            counterInd = counterInd + numProjSub - 1
-        if i == 1:
-            r_old = fistaRecon.r.copy()
-            
-        ## FINAL
-        print ("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)
-            # updating r
-            r_x = fistaRecon.r + ((t_old-1)/t) * (fistaRecon.r - r_old)
-        
-
-        if fistaRecon.getParameter('region_of_interest') is None:
-            string = 'Iteration Number {0} | Objective {1} \n'
-            print (string.format( i, objective[i]))
-        else:
-            ROI , X_ideal = fistaRecon.getParameter('region_of_interest',
-                                                    'ideal_image')
-            
-            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]))
-
-        results.append(X[10])
-    numpy.save("X_out_os.npy", X)
-
-else:
-    
-    
-    
-    astradevice = AstraDevice(DeviceModel.DeviceType.PARALLEL3D.value,
-                [proj_geom['DetectorRowCount'] ,
-                 proj_geom['DetectorColCount'] ,
-                 proj_geom['DetectorSpacingX'] ,
-                 proj_geom['DetectorSpacingY'] ,
-                 proj_geom['ProjectionAngles']
-                 ],
-                [
-                    vol_geom['GridColCount'],
-                    vol_geom['GridRowCount'], 
-                    vol_geom['GridSliceCount'] ] )
-    regul = Regularizer(Regularizer.Algorithm.FGP_TV)
-    regul.setParameter(regularization_parameter=5e6,
-                       number_of_iterations=50,
-                       tolerance_constant=1e-4,
-                       TV_penalty=Regularizer.TotalVariationPenalty.isotropic)
-
-    fistaRecon = FISTAReconstructor(proj_geom,
-                                vol_geom,
-                                Sino3D ,
-                                weights=Weights3D,
-                                device=astradevice,
-                                #regularizer = regul,
-                                subsets=8)
-
-    print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
-    fistaRecon.setParameter(number_of_iterations = 1)
-    fistaRecon.setParameter(Lipschitz_constant = 767893952.0)
-    fistaRecon.setParameter(ring_alpha = 21)
-    fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
-    #fistaRecon.setParameter(subsets=8)
-    
-    #lc = fistaRecon.getParameter('Lipschitz_constant')
-    #fistaRecon.getParameter('regularizer').setParameter(regularization_parameter=5e6/lc)
-    
-    fistaRecon.prepareForIteration()
-    X = fistaRecon.iterate(numpy.load("X.npy"))
-    
-
-# plot
-fig = plt.figure()
-cols = 3
-
-## add the difference
-rd = []
-for i in range(1,len(results)):
-    rd.append(results[i-1])
-    rd.append(results[i])
-    rd.append(results[i] - results[i-1])
-
-rows = (lambda x: int(numpy.floor(x/cols) + 1) if x%cols != 0  else int(x/cols)) \
-       (len (rd))
-for i in range(len (results)):
-    a=fig.add_subplot(rows,cols,i+1)
-    imgplot = plt.imshow(results[i], vmin=0, vmax=1)
-    a.text(0.05, 0.95, "iteration {0}".format(i),
-               verticalalignment='top')
-##    i = i + 1
-##    a=fig.add_subplot(rows,cols,i+1)
-##    imgplot = plt.imshow(results[i], vmin=0, vmax=10)
-##    a.text(0.05, 0.95, "iteration {0}".format(i),
-##               verticalalignment='top')
-    
-##    a=fig.add_subplot(rows,cols,i+2)
-##    imgplot = plt.imshow(results[i]-results[i-1], vmin=0, vmax=10)
-##    a.text(0.05, 0.95, "difference {0}-{1}".format(i, i-1),
-##               verticalalignment='top')
-        
-        
-
-plt.show()
diff --git a/Wrappers/Python/test/test_reconstructor.py b/Wrappers/Python/test/test_reconstructor.py
deleted file mode 100644
index 40065e7..0000000
--- a/Wrappers/Python/test/test_reconstructor.py
+++ /dev/null
@@ -1,359 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Wed Aug 23 16:34:49 2017
-
-@author: ofn77899
-Based on DemoRD2.m
-"""
-
-import h5py
-import numpy
-
-from ccpi.reconstruction.FISTAReconstructor import FISTAReconstructor
-import astra
-import matplotlib.pyplot as plt
-from ccpi.imaging.Regularizer import Regularizer
-from ccpi.reconstruction.AstraDevice import AstraDevice
-from ccpi.reconstruction.DeviceModel import DeviceModel
-
-def RMSE(signal1, signal2):
-    '''RMSE Root Mean Squared Error'''
-    if numpy.shape(signal1) == numpy.shape(signal2):
-        err = (signal1 - signal2)
-        err = numpy.sum( err * err )/numpy.size(signal1);  # MSE
-        err = sqrt(err);                                   # RMSE
-        return err
-    else:
-        raise Exception('Input signals must have the same shape')
-
-def createAstraDevice(projector_geometry, output_geometry):
-    '''TODO remove'''
-    
-    device = AstraDevice(DeviceModel.DeviceType.PARALLEL3D.value,
-                [projector_geometry['DetectorRowCount'] ,
-                 projector_geometry['DetectorColCount'] ,
-                 projector_geometry['DetectorSpacingX'] ,
-                 projector_geometry['DetectorSpacingY'] ,
-                 projector_geometry['ProjectionAngles']
-                 ],
-                [
-                    output_geometry['GridColCount'],
-                    output_geometry['GridRowCount'], 
-                    output_geometry['GridSliceCount'] ] )
-    return device
-
-filename = r'/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/demos/DendrData.h5'
-nx = h5py.File(filename, "r")
-#getEntry(nx, '/')
-# I have exported the entries as children of /
-entries = [entry for entry in nx['/'].keys()]
-print (entries)
-
-Sino3D = numpy.asarray(nx.get('/Sino3D'), dtype="float32")
-Weights3D = numpy.asarray(nx.get('/Weights3D'), dtype="float32")
-angSize = numpy.asarray(nx.get('/angSize'), dtype=int)[0]
-angles_rad = numpy.asarray(nx.get('/angles_rad'), dtype="float32")
-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 = 20
-det_row_count = Z_slices
-# next definition is just for consistency of naming
-det_col_count = size_det
-
-detectorSpacingX = 1.0
-detectorSpacingY = detectorSpacingX
-
-
-proj_geom = astra.creators.create_proj_geom('parallel3d',
-                                            detectorSpacingX,
-                                            detectorSpacingY,
-                                            det_row_count,
-                                            det_col_count,
-                                            angles_rad)
-
-#vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices);
-image_size_x = recon_size
-image_size_y = recon_size
-image_size_z = Z_slices
-vol_geom = astra.creators.create_vol_geom( image_size_x,
-                                           image_size_y,
-                                           image_size_z)
-
-## First pass the arguments to the FISTAReconstructor and test the
-## Lipschitz constant
-
-##fistaRecon = FISTAReconstructor(proj_geom,
-##                                vol_geom,
-##                                Sino3D ,
-##                                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)
-##
-##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
-    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 False:
-    print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
-    print ("prepare for iteration")
-    fistaRecon.prepareForIteration()
-    
-    
-
-    print("initializing  ...")
-    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 ("initialized")
-    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)); 
-
-      
-    t = 1
-
-
-    print ("starting iterations")
-##    % 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'] == 'fanflat' or \
-           fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat_vec' :
-            # 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):
-                sino_id, sino_updt[kkk] = \
-                         astra.creators.create_sino3d_gpu(
-                             X_t[kkk:kkk+1], proj_geom, vol_geom)
-                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.creators.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 :
-             print ("ring removal")
-             for kkk in range(anglesNumb):
-                 
-                 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
-
-        
-
-        # Projection/Backprojection Routine
-        if fistaRecon.getParameter('projector_geometry')['type'] == 'parallel' or \
-           fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat' or\
-           fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat_vec':
-            x_temp = numpy.zeros(numpy.shape(X),dtype=numpy.float32)
-            print ("Projection/Backprojection Routine")
-            for kkk in range(SlicesZ):
-                
-                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)
-        print ("skipping regularizer")
-
-
-        ## FINAL
-        print ("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('region_of_interest') is None:
-            string = 'Iteration Number {0} | Objective {1} \n'
-            print (string.format( i, objective[i]))
-        else:
-            ROI , X_ideal = fistaRecon.getParameter('region_of_interest',
-                                                    'ideal_image')
-            
-            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]))
-            
-##        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
-else:
-    
-    # create a device for forward/backprojection
-    #astradevice = createAstraDevice(proj_geom, vol_geom)
-    
-    astradevice = AstraDevice(DeviceModel.DeviceType.PARALLEL3D.value,
-                [proj_geom['DetectorRowCount'] ,
-                 proj_geom['DetectorColCount'] ,
-                 proj_geom['DetectorSpacingX'] ,
-                 proj_geom['DetectorSpacingY'] ,
-                 proj_geom['ProjectionAngles']
-                 ],
-                [
-                    vol_geom['GridColCount'],
-                    vol_geom['GridRowCount'], 
-                    vol_geom['GridSliceCount'] ] )
-
-    regul = Regularizer(Regularizer.Algorithm.FGP_TV)
-    regul.setParameter(regularization_parameter=5e6,
-                       number_of_iterations=50,
-                       tolerance_constant=1e-4,
-                       TV_penalty=Regularizer.TotalVariationPenalty.isotropic)
-    
-    fistaRecon = FISTAReconstructor(proj_geom,
-                                vol_geom,
-                                Sino3D ,
-                                device = astradevice,
-                                weights=Weights3D,
-                                regularizer = regul
-                                )
-
-    print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
-    fistaRecon.setParameter(number_of_iterations = 18)
-    fistaRecon.setParameter(Lipschitz_constant = 767893952.0)
-    fistaRecon.setParameter(ring_alpha = 21)
-    fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
-
-    
-    
-    fistaRecon.prepareForIteration()
-    X = numpy.load("X.npy")
-
-    
-    X = fistaRecon.iterate(X)
-    #numpy.save("X_out.npy", X)