def astrarecon(tilt_data,tilt_angles,iterations=1,geometry='parallel3d',SO_dist=1.0,OD_dist=1.0):
    proj_shape = np.shape(tilt_data)
    recon_shape = (proj_shape[2],proj_shape[2],proj_shape[1])
    
    vol_geom = astra.create_vol_geom(recon_shape)

    angles = np.pi*tilt_angles/180
    
    if geometry == 'parallel3d':
        proj_geom = astra.create_proj_geom(geometry, 1.0, 1.0, proj_shape[1], proj_shape[2], angles)
        cfg = astra.astra_dict('SIRT3D_CUDA')
    elif geometry == 'cone':
        proj_geom = astra.create_proj_geom(geometry, 1.0, 1.0, proj_shape[1], proj_shape[2], angles, SO_dist, OD_dist)
        cfg = astra.astra_dict('FDK_CUDA')
        
    proj_id = astra.data3d.create('-proj3d', proj_geom, np.swapaxes(tilt_data,0,1))
    
    rec_id = astra.data3d.create('-vol', vol_geom)
    
    cfg['ReconstructionDataId'] = rec_id
    cfg['ProjectionDataId'] = proj_id
    
    alg_id = astra.algorithm.create(cfg)
    
    astra.algorithm.run(alg_id, iterations)
    
    rec = astra.data3d.get(rec_id)
    
    astra.algorithm.delete(alg_id)
    astra.data3d.delete(rec_id)
    astra.data3d.delete(proj_id)
    
    return(rec)

def astra_rec_cpu(tomo, center, recon, theta, vol_geom, niter, proj_type, opts):
    # Lazy import ASTRA
    import astra as astra_mod
    nslices, nang, ndet = tomo.shape
    cfg = astra_mod.astra_dict(opts['method'])
    if 'extra_options' in opts:
        cfg['option'] = opts['extra_options']
    proj_geom = astra_mod.create_proj_geom('parallel', 1.0, ndet, theta.astype(np.float64))
    pid = astra_mod.create_projector(proj_type, proj_geom, vol_geom)
    sino = np.zeros((nang, ndet), dtype=np.float32)
    sid = astra_mod.data2d.link('-sino', proj_geom, sino)
    cfg['ProjectorId'] = pid
    cfg['ProjectionDataId'] = sid
    for i in range(nslices):
        shft = int(np.round(ndet / 2. - center[i]))
        if not shft == 0:
            sino[:] = np.roll(tomo[i], shft)
            l = shft
            r = ndet + shft
            if l < 0:
                l = 0
            if r > ndet:
                r = ndet
            sino[:, :l] = 0
            sino[:,  r:] = 0
        else:
            sino[:] = tomo[i]
        vid = astra_mod.data2d.link('-vol', vol_geom, recon[i])
        cfg['ReconstructionDataId'] = vid
        alg_id = astra_mod.algorithm.create(cfg)
        astra_mod.algorithm.run(alg_id, niter)
        astra_mod.algorithm.delete(alg_id)
        astra_mod.data2d.delete(vid)
    astra_mod.data2d.delete(sid)
    astra_mod.projector.delete(pid)

def astra_rec_cuda(tomo, center, recon, theta, vol_geom, niter, proj_type, gpu_index, opts):
    # Lazy import ASTRA
    import astra as astra_mod
    nslices, nang, ndet = tomo.shape
    cfg = astra_mod.astra_dict(opts['method'])
    if 'extra_options' in opts:
        # NOTE: we are modifying 'extra_options' and so need to make a copy
        cfg['option'] = copy.deepcopy(opts['extra_options'])
    else:
        cfg['option'] = {}
    if gpu_index is not None:
        cfg['option']['GPUindex'] = gpu_index
    oc = None
    const_theta = np.ones(nang)
    proj_geom = astra_mod.create_proj_geom(
        'parallel', 1.0, ndet, theta.astype(np.float64))
    for i in range(nslices):
        if center[i] != oc:
            oc = center[i]
            proj_geom['option'] = {
                'ExtraDetectorOffset':
                (center[i] - ndet / 2.) * const_theta}
        pid = astra_mod.create_projector(proj_type, proj_geom, vol_geom)
        cfg['ProjectorId'] = pid
        sid = astra_mod.data2d.link('-sino', proj_geom, tomo[i])
        cfg['ProjectionDataId'] = sid
        vid = astra_mod.data2d.link('-vol', vol_geom, recon[i])
        cfg['ReconstructionDataId'] = vid
        alg_id = astra_mod.algorithm.create(cfg)
        astra_mod.algorithm.run(alg_id, niter)
        astra_mod.algorithm.delete(alg_id)
        astra_mod.data2d.delete(vid)
        astra_mod.data2d.delete(sid)
        astra_mod.projector.delete(pid)

 def init(self, volume_data=1, projection_data=1):
     # Create volume geometry
     self.volume_geom = astra.create_vol_geom(self.num_voxel)
 
     # Create projection geometry
     self.projection_geom = astra.create_proj_geom(self.geometry_type,
                                        self.detector_spacing_x, self.detector_spacing_y,
                                        self.det_row_count, self.det_col_count,
                                        self.angles,
                                        self.source_origin, self.origin_detector)
 
     # Allocate and store volume data in ASTRA memory
     self.volume_id = astra.data3d.create('-vol', self.volume_geom, volume_data)
 
     # Allocate and store projection data in ASTRA memeory
     self.projection_id = astra.data3d.create('-sino', self.projection_geom, projection_data)
     
     # Create algorithm object: forward projector
     cfg = astra.astra_dict('FP3D_CUDA')
     cfg['option'] = {'GPUindex': self.gpu_index}
     cfg['ProjectionDataId'] = self.projection_id
     cfg['VolumeDataId'] = self.volume_id
     self.forward_alg_id = astra.algorithm.create(cfg)
     
     # Create algorithm object: backward projector
     cfg = astra.astra_dict('BP3D_CUDA')
     cfg['option'] = {'GPUindex': self.gpu_index}
     cfg['ProjectionDataId'] = self.projection_id
     cfg['ReconstructionDataId'] = self.volume_id
     self.backward_alg_id = astra.algorithm.create(cfg)

 def run(self, iterations):
     self.fn = getFilterFile(self.fd, self.pg, iterations, self.rg, self.osz)
     if os.path.exists(self.fn):
         flt = np.load(self.fn)
         self.v[:] = self.customFBP(flt)
         return
     nd = self.nd
     if self.osz:
         nds = self.osz
     else:
         nds = nd
     na = len(self.ang)
     pgc = astra.create_proj_geom('parallel',1.0,nd,self.ang)
     vgc = astra.create_vol_geom((nds,nds))
     pidc = astra.create_projector('strip',pgc,vgc)
     x = np.zeros((nds,nds))
     xs = np.zeros((nds,nds))
     sf = np.zeros((na,nd))
     vid = astra.data2d.create('-vol',vgc)
     sid = astra.data2d.create('-sino',pgc)
     cfg = astra.astra_dict('FP')
     cfg['ProjectorId']=pidc
     cfg['ProjectionDataId']=sid
     cfg['VolumeDataId']=vid
     fpid = astra.algorithm.create(cfg)
     cfg = astra.astra_dict('BP')
     cfg['ProjectorId']=pidc
     cfg['ProjectionDataId']=sid
     cfg['ReconstructionDataId']=vid
     bpid = astra.algorithm.create(cfg)
     vc = astra.data2d.get_shared(vid)
     sc = astra.data2d.get_shared(sid)
     x[nds//2,nds//2]=1
     alp = 1./(na*nds)
     if self.rg:
         if self.rg*alp >=0.1:
             alp = 0.1/self.rg
     astra.log.info('Computing filter...')
     for i in range(iterations):
         if i%10==0: astra.log.info('{:.2f} % done'.format(100*float(i)/iterations))
         xs+=x
         vc[:] = x
         astra.algorithm.run(fpid)
         astra.algorithm.run(bpid)
         if self.rg:
             dx = x[:-1,:] - x[1:,:]
             dy = x[:,:-1] - x[:,1:]
             x[:-1,:] -= self.rg*dx*alp
             x[1:,:] += self.rg*dx*alp
             x[:,:-1] -= self.rg*dy*alp
             x[:,1:] += self.rg*dy*alp
         x -= vc*alp
     vc[:] = xs
     astra.algorithm.run(fpid)
     flt = sc.copy()*alp
     astra.algorithm.delete([fpid,bpid])
     astra.algorithm.delete([vid,sid])
     np.save(self.fn,flt)
     self.v[:] = self.customFBP(flt)
     astra.projector.delete(pidc)

 def geom_setup_2D(self, sino, angles, shape, cors):
     p_low, p_high = self.array_pad(cors, self.nCols)
     sino = np.pad(sino, ((0, 0), (p_low, p_high)), mode='reflect')
     vol_geom = astra.create_vol_geom(shape[0], shape[1])
     proj_geom = astra.create_proj_geom('parallel', 1.0, sino.shape[1],
                                        np.deg2rad(angles))
     return sino, vol_geom, proj_geom

    def __init__(self,
                 geometry_obj=Geometry(1),
                 vol_vector=None,
                 proj_vector=None, gpu_index=0):
        self.geom = geometry_obj
        if vol_vector is None:
            self.vol = Rn(self.geom.vol_size)
        else:
            self.vol = vol_vector
        if proj_vector is None:
            self.proj = Rn(self.geom.proj_size)
        else:
            self.proj = proj_vector
        self.gpu_index = gpu_index
        self.bp_id = None
        self.fp_id = None

        # Create volume geometry
        self.vol_geom = astra.create_vol_geom(self.geom.vol_shape)

        # Create projection geometry
        self.proj_geom = astra.create_proj_geom(
            self.geom.geom_type,
            self.geom.detector_spacing_x, self.geom.detector_spacing_y,
            self.geom.det_row_count, self.geom.det_col_count,
            self.geom.angles,
            self.geom.source_origin, self.geom.origin_detector)

        # Allocate ASTRA memory for volume data
        self.volume_id = astra.data3d.create('-vol', self.vol_geom)

        # Allocate ASTRA memory for projection data
        self.proj_id = astra.data3d.create('-sino', self.proj_geom)

def gen_random_geometry_parallel():
  if not NONUNITDET:
    w = 1.0
  else:
    w = 0.8 + 0.4 * np.random.random()
  pg = astra.create_proj_geom('parallel', w, np.random.randint(*range2d), np.linspace(0, 2*np.pi, np.random.randint(*range2d), endpoint=False))
  return pg

def gen_random_geometry_fanflat():
  if not NONUNITDET:
    w = 1.0
  else:
    w = 0.6 + 0.8 * np.random.random()
  pg = astra.create_proj_geom('fanflat', w, np.random.randint(*range2d), np.linspace(0, 2*np.pi, np.random.randint(*range2d), endpoint=False), 256 * (0.5 + np.random.random()), 256 * np.random.random())
  return pg

def bp0(sino, det_col=111, num_angles=222, voxel_size_mm=1):
    """Wrapper for astra forward projector

    :param sino:
    :param projector_id:
    :return: backprojected sinogram
    """

    vol_geom = astra.create_vol_geom(sino.shape)
    proj_id = astra.create_projector(
        'cuda',
        astra.create_proj_geom('parallel', 1.0, det_col,
                               np.linspace(0, np.pi, num_angles, False)),
        vol_geom)


    rec_id, backprojection = astra.create_backprojection(sino * voxel_size_mm,
                                                   proj_id)

    # rec_id, backprojection = astra.create_backprojection(sino, projector_id)

    # backprojection /= sino.shape[0]
    # backprojection *= np.pi
    astra.data2d.delete(rec_id)
    astra.projector.delete(proj_id)

    return backprojection

def pid(number_of_angles):
    """Return projector ID.
    :param number_of_angles:
    :return:
    """
    return astra.create_projector('cuda', astra.create_proj_geom('parallel',
    1.0, vshape[0], np.linspace(0, np.pi, number_of_angles, False)), vol_geom)

def recon_sirt_fbp(im, angles, iter, temppath ):
	"""Reconstruct a sinogram with the SIRT-FBP algorithm (Pelt, 2015).

	Parameters
	----------
	im : array_like
		Sinogram image data as numpy array.

    iter : int
		Number of iterations to be used for the computation of SIRT filter.

	angles : double
		Value in radians representing the number of angles of the input sinogram.
	
	"""
	# Create ASTRA geometries:
	vol_geom = astra.create_vol_geom(im.shape[1] , im.shape[1])
	proj_geom = astra.create_proj_geom('parallel', 1.0, im.shape[1], linspace(0,angles,im.shape[0],False))
	proj = astra.create_projector('cuda', proj_geom, vol_geom)
	p = astra.OpTomo(proj)

	# Register plugin with ASTRA
	astra.plugin.register(sirtfbp.plugin)

	# Create the ASTRA projector
	im_rec = p.reconstruct('SIRT-FBP',im, iter, extraOptions={'filter_dir':temppath})
	
	return im_rec.astype(float32)

def recon_mr_fbp(im, angles):
	"""Reconstruct a sinogram with the Minimum Residual FBP algorithm (Pelt, 2013).

	Parameters
	----------
	im : array_like
		Sinogram image data as numpy array.

	angles : double
		Value in radians representing the number of angles of the input sinogram.
	
	"""
	# Create ASTRA geometries:
	vol_geom = astra.create_vol_geom(im.shape[1] , im.shape[1])
	proj_geom = astra.create_proj_geom('parallel', 1.0, im.shape[1], linspace(0,angles,im.shape[0],False))

	# Create the ASTRA projector:
	p = mrfbp.ASTRAProjector.ASTRAProjector2D(proj_geom,vol_geom)
	
	# Create the MR-FBP Reconstructor:
	rec = mrfbp.Reconstructor(p)

	# Reconstruct the image using MR-FBP:
	im_rec = rec.reconstruct(im)
	
	return im_rec.astype(float32)

def astrafp(im, ang, prj="cuda"):
    proj_geom = astra.create_proj_geom("parallel", 1.0, im.shape[0], np.array([ang, 0]))
    vol_geom = astra.create_vol_geom(im.shape)
    pid = astra.create_projector(prj, proj_geom, vol_geom)
    w = astra.OpTomo(pid)
    fpim = w * im
    astra.projector.delete(pid)
    return fpim[0 : im.shape[0]]

 def __init__(self, n_pixels, n_angles, rayperdetec=None):
     '''
     Initialize the ASTRA toolbox with a simple parallel configuration.
     The image is assumed to be square, and the detector count is equal to the number of rows/columns.
     '''
     self.vol_geom = astra.create_vol_geom(n_pixels, n_pixels)
     self.proj_geom = astra.create_proj_geom('parallel', 1.0, n_pixels, np.linspace(0,np.pi,n_angles,False))
     self.proj_id = astra.create_projector('cuda', self.proj_geom, self.vol_geom)

    def reconstruct(self, sinogram, centre_of_rotation, angles, shape, center):

        ctr = centre_of_rotation
        width = sinogram.shape[1]
        pad = 50

        sino = np.nan_to_num(sinogram)

        # pad the array so that the centre of rotation is in the middle
        alen = ctr
        blen = width - ctr
        mid = width / 2.0

        if (ctr > mid):
            plow = pad
            phigh = (alen - blen) + pad
        else:
            plow = (blen - alen) + pad
            phigh = pad

        logdata = np.log(sino+1)
        sinogram = np.pad(logdata, ((0, 0), (int(plow), int(phigh))),
                          mode='reflect')

        width = sinogram.shape[1]

        vol_geom = astra.create_vol_geom(shape[0], shape[1])
        proj_geom = astra.create_proj_geom('parallel', 1.0, width,
                                           np.deg2rad(angles))

        sinogram_id = astra.data2d.create("-sino", proj_geom, sinogram)

        # Create a data object for the reconstruction
        rec_id = astra.data2d.create('-vol', vol_geom)

        proj_id = astra.create_projector('strip', proj_geom, vol_geom)

        cfg = astra.astra_dict(self.parameters['reconstruction_type'])
        cfg['ReconstructionDataId'] = rec_id
        cfg['ProjectionDataId'] = sinogram_id
        cfg['ProjectorId'] = proj_id

        # Create the algorithm object from the configuration structure
        alg_id = astra.algorithm.create(cfg)
        # Run 20 iterations of the algorithm
        itterations = int(self.parameters['number_of_iterations'])
        # This will have a runtime in the order of 10 seconds.
        astra.algorithm.run(alg_id, itterations)
        # Get the result
        rec = astra.data2d.get(rec_id)

        # Clean up.
        astra.algorithm.delete(alg_id)
        astra.data2d.delete(rec_id)
        astra.data2d.delete(sinogram_id)

        return rec

def _basic_par3d_fp():
  import astra
  import numpy as np
  vg = astra.create_vol_geom(2, 32, 32)
  pg = astra.create_proj_geom('parallel3d', 1, 1, 32, 32, [0])
  vol = np.random.rand(32, 2, 32)
  (sino_id, sino) = astra.create_sino3d_gpu(vol, pg, vg)
  astra.data3d.delete(sino_id)
  err = np.max(np.abs(sino[:,0,:] - np.sum(vol,axis=1)))
  return err < 1e-6

def _basic_par2d_fp(type):
  import astra
  import numpy as np
  vg = astra.create_vol_geom(2, 32)
  pg = astra.create_proj_geom('parallel', 1, 32, [0])
  proj_id = astra.create_projector(type, pg, vg)
  vol = np.random.rand(2, 32)
  (sino_id, sino) = astra.create_sino(vol, proj_id)
  astra.data2d.delete(sino_id)
  astra.projector.delete(proj_id)
  err = np.max(np.abs(sino[0,:] - np.sum(vol,axis=0)))
  return err < 1e-6

    def geom_setup_3D(self, sino, angles, shape, cors):
        nSinos = sino.shape[self.slice_dim]
        length = len(angles)
        angles = np.deg2rad(angles)

        vectors = np.zeros((length, 12))
        for i in range(len(angles)):
            # ray direction
            vectors[i, 0] = np.sin(angles[i])
            vectors[i, 1] = -np.cos(angles[i])
            vectors[i, 2] = 0

            # center of detector
            vectors[i, 3:6] = 0

            # vector from detector pixel (0,0) to (0,1)
            vectors[i, 6] = np.cos(angles[i])
            vectors[i, 7] = np.sin(angles[i])
            vectors[i, 8] = 0

            # vector from detector pixel (0,0) to (1,0)
            vectors[i, 9] = 0
            vectors[i, 10] = 0
            vectors[i, 11] = 1

#        i = range(length)
#        # ray direction
#        vectors[i, 0] = np.sin(theta[i])
#        vectors[i, 1] = -np.cos(theta[i])
#        vectors[i, 2] = 0
#        # detector centre (use this for translation)
#        # assuming all sinograms are translated by the same amount for now
#        #det_vec = [cors[0], 0, 0]
#        det_vec = [0, 0, 0]
#        vectors[i, 3:6] = det_vec
#        # (use the following vectors for rotation)
#        # vector from detector pixel (0,0) to (0,1)
#        vectors[i, 6] = np.cos(theta[i])
#        vectors[i, 7] = np.sin(theta[i])
#        vectors[i, 8] = 0
#        # vector from detector pixel (0,0) to (1,0)
#        vectors[i, 9:12] = [0, 0, 1]

        # Parameters: #rows, #columns, vectors
        vol_geom = astra.create_vol_geom(nSinos, shape[0], shape[2])
        proj_geom = astra.create_proj_geom('parallel3d_vec', sino.shape[1],
                                           sino.shape[2], vectors)
        return vol_geom, proj_geom

    def astra_2D_recon(self, data):
        sino = data[0]
        cor, angles, vol_shape, init = self.get_frame_params()
        angles = np.deg2rad(angles)
        if self.res:
            res = np.zeros(self.len_res)
        # create volume geom
        vol_geom = astra.create_vol_geom(vol_shape)
        # create projection geom
        det_width = sino.shape[self.dim_detX]
        proj_geom = astra.create_proj_geom('parallel', 1.0, det_width, angles)
        sino = np.transpose(sino, (self.dim_rot, self.dim_detX))

        # create sinogram id
        sino_id = astra.data2d.create("-sino", proj_geom, sino)
        # create reconstruction id
        if init is not None:
            rec_id = astra.data2d.create('-vol', vol_geom, init)
        else:
            rec_id = astra.data2d.create('-vol', vol_geom)

#        if self.mask_id:
#            self.mask_id = astra.data2d.create('-vol', vol_geom, self.mask)
        # setup configuration options
        cfg = self.set_config(rec_id, sino_id, proj_geom, vol_geom)
        # create algorithm id
        alg_id = astra.algorithm.create(cfg)
        # run algorithm
        if self.res:
            for j in range(self.iters):
                # Run a single iteration
                astra.algorithm.run(alg_id, 1)
                res[j] = astra.algorithm.get_res_norm(alg_id)
        else:
            astra.algorithm.run(alg_id, self.iters)
        # get reconstruction matrix

        if self.manual_mask is not False:
            recon = self.manual_mask*astra.data2d.get(rec_id)
        else:
            recon = astra.data2d.get(rec_id)

        # delete geometry
        self.delete(alg_id, sino_id, rec_id, False)
        return [recon, res] if self.res else recon