def reconstruct(h5fname, sino, rot_center, binning, algorithm='gridrec'):

    sample_detector_distance = 30       # Propagation distance of the wavefront in cm
    detector_pixel_size_x = 1.17e-4     # Detector pixel size in cm (5x: 1.17e-4, 2X: 2.93e-4)
    monochromator_energy = 25.74        # Energy of incident wave in keV
    alpha = 1e-02                       # Phase retrieval coeff.
    zinger_level = 1000                 # Zinger level for projections
    zinger_level_w = 1000               # Zinger level for white

    miss_angles = [141,226]

    # Read APS 32-BM raw data.
    proj, flat, dark, theta = dxchange.read_aps_32id(h5fname, sino=sino)
        

    print (theta)
    # Manage the missing angles:
    #proj_size = np.shape(proj)
    #theta = np.linspace(0,180,proj_size[0])
    proj = np.concatenate((proj[0:miss_angles[0],:,:], proj[miss_angles[1]+1:-1,:,:]), axis=0)
    theta = np.concatenate((theta[0:miss_angles[0]], theta[miss_angles[1]+1:-1]))

    # zinger_removal
    #proj = tomopy.misc.corr.remove_outlier(proj, zinger_level, size=15, axis=0)
    #flat = tomopy.misc.corr.remove_outlier(flat, zinger_level_w, size=15, axis=0)

    # Flat-field correction of raw data.
    data = tomopy.normalize(proj, flat, dark, cutoff=0.8)

    # remove stripes
    data = tomopy.remove_stripe_fw(data,level=7,wname='sym16',sigma=1,pad=True)

    # phase retrieval
    # data = tomopy.prep.phase.retrieve_phase(data,pixel_size=detector_pixel_size_x,dist=sample_detector_distance,energy=monochromator_energy,alpha=alpha,pad=True)

    print("Raw data: ", h5fname)
    print("Center: ", rot_center)

    data = tomopy.minus_log(data)

    data = tomopy.remove_nan(data, val=0.0)
    data = tomopy.remove_neg(data, val=0.00)
    data[np.where(data == np.inf)] = 0.00

    rot_center = rot_center/np.power(2, float(binning))
    data = tomopy.downsample(data, level=binning) 
    data = tomopy.downsample(data, level=binning, axis=1)

    # Reconstruct object.
    if algorithm == 'sirtfbp':
        rec = rec_sirtfbp(data, theta, rot_center)
    else:
        rec = tomopy.recon(data, theta, center=rot_center, algorithm=algorithm, filter_name='parzen')
        
    print("Algorithm: ", algorithm)

    # Mask each reconstructed slice with a circle.
    ##rec = tomopy.circ_mask(rec, axis=0, ratio=0.95)
    
    return rec

def rec_try(h5fname, nsino, rot_center, center_search_width, algorithm, binning):
    zinger_level = 800                  # Zinger level for projections
    zinger_level_w = 1000               # Zinger level for white
    
    data_shape = get_dx_dims(h5fname, 'data')
    print(data_shape)
    ssino = int(data_shape[1] * nsino)

    center_range = (rot_center-center_search_width, rot_center+center_search_width, 0.5)
    #print(sino,ssino, center_range)
    #print(center_range[0], center_range[1], center_range[2])

    # Select sinogram range to reconstruct
    sino = None
        
    start = ssino
    end = start + 1
    sino = (start, end)

    # Read APS 32-BM raw data.
    proj, flat, dark, theta = dxchange.read_aps_32id(h5fname, sino=sino)

    # zinger_removal
    proj = tomopy.misc.corr.remove_outlier(proj, zinger_level, size=15, axis=0)
    flat = tomopy.misc.corr.remove_outlier(flat, zinger_level_w, size=15, axis=0)
        
    # Flat-field correction of raw data.
    data = tomopy.normalize(proj, flat, dark, cutoff=1.4)

    # remove stripes
    data = tomopy.remove_stripe_fw(data,level=7,wname='sym16',sigma=1,pad=True)


    print("Raw data: ", h5fname)
    print("Center: ", rot_center)

    data = tomopy.minus_log(data)

    stack = np.empty((len(np.arange(*center_range)), data_shape[0], data_shape[2]))

    index = 0
    for axis in np.arange(*center_range):
        stack[index] = data[:, 0, :]
        index = index + 1

    # Reconstruct the same slice with a range of centers.
    rec = tomopy.recon(stack, theta, center=np.arange(*center_range), sinogram_order=True, algorithm='gridrec', filter_name='parzen', nchunk=1)

    # Mask each reconstructed slice with a circle.
    rec = tomopy.circ_mask(rec, axis=0, ratio=0.95)

    index = 0
    # Save images to a temporary folder.
    fname = os.path.dirname(h5fname) + '/' + 'try_rec/' + 'recon_' + os.path.splitext(os.path.basename(h5fname))[0]    
    for axis in np.arange(*center_range):
        rfname = fname + '_' + str('{0:.2f}'.format(axis) + '.tiff')
        dxchange.write_tiff(rec[index], fname=rfname, overwrite=True)
        index = index + 1

    print("Reconstructions: ", fname)

def reconstruct(h5fname, sino, rot_center, binning, algorithm='gridrec'):

    sample_detector_distance = 8        # Propagation distance of the wavefront in cm
    detector_pixel_size_x = 2.247e-4    # Detector pixel size in cm (5x: 1.17e-4, 2X: 2.93e-4)
    monochromator_energy = 24.9         # Energy of incident wave in keV
    alpha = 1e-02                       # Phase retrieval coeff.
    zinger_level = 800                  # Zinger level for projections
    zinger_level_w = 1000               # Zinger level for white

    # Read APS 32-BM raw data.
    proj, flat, dark, theta = dxchange.read_aps_32id(h5fname, sino=sino)
        
    # zinger_removal
    # proj = tomopy.misc.corr.remove_outlier(proj, zinger_level, size=15, axis=0)
    # flat = tomopy.misc.corr.remove_outlier(flat, zinger_level_w, size=15, axis=0)

    # Flat-field correction of raw data.
    ##data = tomopy.normalize(proj, flat, dark, cutoff=0.8)
    data = tomopy.normalize(proj, flat, dark)

    # remove stripes
    data = tomopy.remove_stripe_fw(data,level=7,wname='sym16',sigma=1,pad=True)

    # data = tomopy.remove_stripe_ti(data, alpha=1.5)
    # data = tomopy.remove_stripe_sf(data, size=150)

    # phase retrieval
    #data = tomopy.prep.phase.retrieve_phase(data,pixel_size=detector_pixel_size_x,dist=sample_detector_distance,energy=monochromator_energy,alpha=alpha,pad=True)

    print("Raw data: ", h5fname)
    print("Center: ", rot_center)

    data = tomopy.minus_log(data)

    data = tomopy.remove_nan(data, val=0.0)
    data = tomopy.remove_neg(data, val=0.00)
    data[np.where(data == np.inf)] = 0.00

    rot_center = rot_center/np.power(2, float(binning))
    data = tomopy.downsample(data, level=binning) 
    data = tomopy.downsample(data, level=binning, axis=1)

    # Reconstruct object.
    if algorithm == 'sirtfbp':
        rec = rec_sirtfbp(data, theta, rot_center)
    elif algorithm == 'astrasirt':
        extra_options ={'MinConstraint':0}
        options = {'proj_type':'cuda', 'method':'SIRT_CUDA', 'num_iter':200, 'extra_options':extra_options}
        rec = tomopy.recon(data, theta, center=rot_center, algorithm=tomopy.astra, options=options)
    else:
        rec = tomopy.recon(data, theta, center=rot_center, algorithm=algorithm, filter_name='parzen')
        
    print("Algorithm: ", algorithm)

    # Mask each reconstructed slice with a circle.
    rec = tomopy.circ_mask(rec, axis=0, ratio=0.95)
    
    return rec

def main(arg):

    parser = argparse.ArgumentParser()
    parser.add_argument("top", help="top directory where the tiff images are located: /data/")
    parser.add_argument("start", nargs='?', const=1, type=int, default=1, help="index of the first image: 1000 (default 1)")

    args = parser.parse_args()

    top = args.top
    index_start = int(args.start)

    template = os.listdir(top)[0]

    nfile = len(fnmatch.filter(os.listdir(top), '*.tif'))
    index_end = index_start + nfile
    ind_tomo = range(index_start, index_end)
    
    fname = top + template

    print (nfile, index_start, index_end, fname)


    # Select the sinogram range to reconstruct.
    start = 0
    end = 512
    sino=(start, end)

    # Read the tiff raw data.
    ndata = dxchange.read_tiff_stack(fname, ind=ind_tomo, slc=(sino, None))
 
    # Normalize to 1 using the air counts
    ndata = tomopy.normalize_bg(ndata, air=5)

    # Set data collection angles as equally spaced between 0-180 degrees.
    theta = tomopy.angles(ndata.shape[0])

    ndata = tomopy.minus_log(ndata)

    # Set binning and number of iterations
    binning = 8
    iters = 21

    print("Original", ndata.shape)
    ndata = tomopy.downsample(ndata, level=binning, axis=1)
#    ndata = tomopy.downsample(ndata, level=binning, axis=2)
    print("Processing:", ndata.shape)

    fdir = 'aligned' + '/noblur_iter_' + str(iters) + '_bin_' + str(binning) 

    print(fdir)
    cprj, sx, sy, conv = alignment.align_seq(ndata, theta, fdir=fdir, iters=iters, pad=(10, 10), blur=False, save=True, debug=True)

    np.save(fdir + '/shift_x', sx)
    np.save(fdir + '/shift_y', sy)

    # Write aligned projections as stack of TIFs.
    dxchange.write_tiff_stack(cprj, fname=fdir + '/radios/image')

def main(arg):

    parser = argparse.ArgumentParser()
    parser.add_argument("top", help="top directory where the tiff images are located: /data/")
    parser.add_argument("start", nargs='?', const=1, type=int, default=1, help="index of the first image: 1000 (default 1)")

    args = parser.parse_args()

    top = args.top
    index_start = int(args.start)

    template = os.listdir(top)[0]

    nfile = len(fnmatch.filter(os.listdir(top), '*.tif'))
    index_end = index_start + nfile
    ind_tomo = range(index_start, index_end)
    
    fname = top + template

    print (nfile, index_start, index_end, fname)


    # Select the sinogram range to reconstruct.
    start = 0
    end = 512
    sino=(start, end)

    # Read the tiff raw data.
    ndata = dxchange.read_tiff_stack(fname, ind=ind_tomo, slc=(sino, None))

    print(ndata.shape)
    binning = 8
    ndata = tomopy.downsample(ndata, level=binning, axis=1)
    print(ndata.shape)
    
    # Normalize to 1 using the air counts
    ndata = tomopy.normalize_bg(ndata, air=5)

    ## slider(ndata)

    # Set data collection angles as equally spaced between 0-180 degrees.
    theta = tomopy.angles(ndata.shape[0])
   
    rot_center = 960
    print("Center of rotation: ", rot_center)

    ndata = tomopy.minus_log(ndata)

    # Reconstruct object using Gridrec algorithm.
    rec = tomopy.recon(ndata, theta, center=rot_center, algorithm='gridrec')

    # Mask each reconstructed slice with a circle.
    rec = tomopy.circ_mask(rec, axis=0, ratio=0.95)

    # Write data as stack of TIFs.
    dxchange.write_tiff_stack(rec, fname='/local/dataraid/mark/rec/recon')

def reconstruct(h5fname, sino, rot_center, binning, algorithm='gridrec'):

    sample_detector_distance = 8        # Propagation distance of the wavefront in cm
    detector_pixel_size_x = 2.247e-4    # Detector pixel size in cm (5x: 1.17e-4, 2X: 2.93e-4)
    monochromator_energy = 24.9         # Energy of incident wave in keV
    alpha = 1e-02                       # Phase retrieval coeff.
    zinger_level = 800                  # Zinger level for projections
    zinger_level_w = 1000               # Zinger level for white

    # h5fname_norm = '/local/data/2019-02/Burke/C47M_0015.h5'
    h5fname_norm = '/local/data/2019-02/Burke/kc78_Menardii_0003.h5'
    proj1, flat, dark, theta1 = dxchange.read_aps_32id(h5fname_norm, sino=sino)
    proj, dummy, dummy1, theta = dxchange.read_aps_32id(h5fname, sino=sino)
        
    # zinger_removal
    proj = tomopy.misc.corr.remove_outlier(proj, zinger_level, size=15, axis=0)
    flat = tomopy.misc.corr.remove_outlier(flat, zinger_level_w, size=15, axis=0)

    # Flat-field correction of raw data.
    ##data = tomopy.normalize(proj, flat, dark, cutoff=0.8)
    data = tomopy.normalize(proj, flat, dark)

    # remove stripes
    data = tomopy.remove_stripe_fw(data,level=7,wname='sym16',sigma=1,pad=True)

    #data = tomopy.remove_stripe_ti(data, alpha=1.5)
    data = tomopy.remove_stripe_sf(data, size=20)

    # phase retrieval
    #data = tomopy.prep.phase.retrieve_phase(data,pixel_size=detector_pixel_size_x,dist=sample_detector_distance,energy=monochromator_energy,alpha=alpha,pad=True)

    print("Raw data: ", h5fname)
    print("Center: ", rot_center)

    data = tomopy.minus_log(data)

    data = tomopy.remove_nan(data, val=0.0)
    data = tomopy.remove_neg(data, val=0.00)
    data[np.where(data == np.inf)] = 0.00

    rot_center = rot_center/np.power(2, float(binning))
    data = tomopy.downsample(data, level=binning) 
    data = tomopy.downsample(data, level=binning, axis=1)

    # Reconstruct object.
    if algorithm == 'sirtfbp':
        rec = rec_sirtfbp(data, theta, rot_center)
    else:
        rec = tomopy.recon(data, theta, center=rot_center, algorithm=algorithm, filter_name='parzen')
        
    print("Algorithm: ", algorithm)

    # Mask each reconstructed slice with a circle.
    rec = tomopy.circ_mask(rec, axis=0, ratio=0.95)
    
    return rec

def reconstruct(h5fname, sino, rot_center, args, blocked_views=None):

    # Read APS 32-BM raw data.
    proj, flat, dark, theta = dxchange.read_aps_32id(h5fname, sino=sino)

    # Manage the missing angles:
    if blocked_views is not None:
        print("Blocked Views: ", blocked_views)
        proj = np.concatenate((proj[0:blocked_views[0], :, :],
                               proj[blocked_views[1]+1:-1, :, :]), axis=0)
        theta = np.concatenate((theta[0:blocked_views[0]],
                                theta[blocked_views[1]+1: -1]))

    # Flat-field correction of raw data.
    data = tomopy.normalize(proj, flat, dark, cutoff=1.4)

    # remove stripes
    data = tomopy.remove_stripe_fw(data, level=7, wname='sym16', sigma=1,
                                   pad=True)

    print("Raw data: ", h5fname)
    print("Center: ", rot_center)

    data = tomopy.minus_log(data)

    data = tomopy.remove_nan(data, val=0.0)
    data = tomopy.remove_neg(data, val=0.00)
    data[np.where(data == np.inf)] = 0.00

    algorithm = args.algorithm
    ncores = args.ncores
    nitr = args.num_iter

    # always add algorithm
    _kwargs = {"algorithm": algorithm}

    # assign number of cores
    _kwargs["ncore"] = ncores

    # don't assign "num_iter" if gridrec or fbp
    if algorithm not in ["fbp", "gridrec"]:
        _kwargs["num_iter"] = nitr

    # Reconstruct object.
    with timemory.util.auto_timer(
        "[tomopy.recon(algorithm='{}')]".format(algorithm)):
        rec = tomopy.recon(proj, theta, **_kwargs)

    # Mask each reconstructed slice with a circle.
    rec = tomopy.circ_mask(rec, axis=0, ratio=0.95)

    return rec

def main(arg):

    parser = argparse.ArgumentParser()
    parser.add_argument("top", help="top directory where the tiff images are located: /data/")
    parser.add_argument("start", nargs='?', const=1, type=int, default=1, help="index of the first image: 10001 (default 1)")

    args = parser.parse_args()

    top = args.top
    index_start = int(args.start)

    template = os.listdir(top)[1]

    nfile = len(fnmatch.filter(os.listdir(top), '*.tif'))
    index_end = index_start + nfile
    ind_tomo = range(index_start, index_end)

    fname = top + template

    # Read the tiff raw data.
    rdata = dxchange.read_tiff_stack(fname, ind=ind_tomo)
    particle_bed_reference = particle_bed_location(rdata[0], plot=False)
    print("Particle bed location: ", particle_bed_reference)
    
    # Cut the images to remove the particle bed
    cdata = rdata[:, 0:particle_bed_reference, :]

    # Find the image when the shutter starts to close
    dark_index = shutter_off(rdata)
    print("shutter closes on image: ", dark_index)
    # Set the [start, end] index of the blocked images, flat and dark.
    flat_range = [0, 1]
    data_range = [48, dark_index]
    dark_range = [dark_index, nfile]

    # # for fast testing
    # data_range = [48, dark_index]

    flat = cdata[flat_range[0]:flat_range[1], :, :]
    proj = cdata[data_range[0]:data_range[1], :, :]
    dark = np.zeros((dark_range[1]-dark_range[0], proj.shape[1], proj.shape[2]))  

    # if you want to use the shutter closed images as dark uncomment this:
    #dark = cdata[dark_range[0]:dark_range[1], :, :]  

    ndata = tomopy.normalize(proj, flat, dark)
    ndata = tomopy.normalize_bg(ndata, air=ndata.shape[2]/2.5)
    ndata = tomopy.minus_log(ndata)
    sharpening(ndata)
    slider(ndata)

def main(arg):

    fname = '/local/dataraid/elettra/Oak_16bit_slice343_all_repack.h5'
    
    # Read the hdf raw data.
    sino, sflat, sdark, th = dxchange.read_aps_32id(fname)

    slider(sino)
    proj = np.swapaxes(sino,0,1)
    flat = np.swapaxes(sflat,0,1)
    dark = np.swapaxes(sdark,0,1)

    # Set data collection angles as equally spaced between 0-180 degrees.
    theta = tomopy.angles(proj.shape[0], ang1=0.0, ang2=180.0)

    print(proj.shape, dark.shape, flat.shape, theta.shape)

    # Flat-field correction of raw data.
    ndata = tomopy.normalize(proj, flat, dark)
    #slider(ndata)

    # Find rotation center.
    rot_center = 962

    binning = 1
    ndata = tomopy.downsample(ndata, level=int(binning))
    rot_center = rot_center/np.power(2, float(binning))    

    ndata = tomopy.minus_log(ndata)
    
    # Reconstruct object using Gridrec algorithm.
    rec = tomopy.recon(ndata, theta, center=rot_center, algorithm='gridrec')

    # Mask each reconstructed slice with a circle.
    rec = tomopy.circ_mask(rec, axis=0, ratio=0.95)

    # Write data as stack of TIFs.
    dxchange.write_tiff_stack(rec, fname='recon_dir/recon')

def reconstruct(sname, rot_center, ovlpfind, s_start, s_end):
    fname = dfolder + sname + '.h5'
    print (fname)
    start = s_start  
    end =   s_end
    chunks = 24 
    num_sino = (end - start) // chunks
    for m in range(chunks):
        sino_start = start + num_sino * m
        sino_end = start + num_sino * (m + 1)
        start_read_time = time.time()
        proj, flat, dark, thetat = dxchange.read_aps_2bm(fname, sino=(sino_start, sino_end))
        print('   done read in %0.1f min' % ((time.time() - start_read_time)/60))
        dark = proj[9001:9002]
        flat = proj[0:1]
        proj = proj[1:9000]
        theta = tomopy.angles(proj.shape[0], 0., 360.)
        proj = tomopy.sino_360_to_180(proj, overlap=ovlpfind, rotation='right')
        proj = tomopy.remove_outlier(proj, dif=0.4)
        proj = tomopy.normalize_bg(proj, air=10)
        proj = tomopy.minus_log(proj)
        center = rot_center
        start_ring_time = time.time()
        proj = tomopy.remove_stripe_fw(proj, wname='sym5', sigma=4, pad=False)
        proj = tomopy.remove_stripe_sf(proj, size=3)
        print('   done pre-process in %0.1f min' % ((time.time() - start_ring_time)/60))
        start_phase_time = time.time()
        proj = tomopy.retrieve_phase(proj, pixel_size=detector_pixel_size_x, dist=sample_detector_distance, energy=energy, alpha=alpha, pad=True, ncore=None, nchunk=None)
        print('   done phase retrieval in %0.1f min' % ((time.time() - start_phase_time)/60))
        start_recon_time = time.time()
        rec = tomopy.recon(proj, theta, center=center, algorithm='gridrec', filter_name='ramalk')
        tomopy.circ_mask(rec, axis=0, ratio=0.95)
        print ("Reconstructed", rec.shape)
        dxchange.write_tiff_stack(rec, fname = dfolder + '/' + sname + '/' + sname, overwrite=True, start=sino_start)
        print('   Chunk reconstruction done in %0.1f min' % ((time.time() - start_recon_time)/60))
    print ("Done!")

def reconstruct(h5fname, sino, rot_center, args, blocked_views=None):

    # Read APS 32-BM raw data.
    proj, flat, dark, theta = dxchange.read_aps_32id(h5fname, sino=sino)

    # Manage the missing angles:
    if blocked_views is not None:
        print("Blocked Views: ", blocked_views)
        proj = np.concatenate((proj[0:blocked_views[0], :, :],
                               proj[blocked_views[1]+1:-1, :, :]), axis=0)
        theta = np.concatenate((theta[0:blocked_views[0]],
                                theta[blocked_views[1]+1: -1]))

    # Flat-field correction of raw data.
    data = tomopy.normalize(proj, flat, dark, cutoff=1.4)

    # remove stripes
    data = tomopy.remove_stripe_fw(data, level=7, wname='sym16', sigma=1,
                                   pad=True)

    print("Raw data: ", h5fname)
    print("Center: ", rot_center)

    data = tomopy.minus_log(data)

    data = tomopy.remove_nan(data, val=0.0)
    data = tomopy.remove_neg(data, val=0.00)
    data[np.where(data == np.inf)] = 0.00

    algorithm = args.algorithm
    ncores = args.ncores
    nitr = args.num_iter

    # always add algorithm
    _kwargs = {"algorithm": algorithm}

    # assign number of cores
    _kwargs["ncore"] = ncores

    # use the accelerated version
    if algorithm in ["mlem", "sirt"]:
        _kwargs["accelerated"] = True

    # don't assign "num_iter" if gridrec or fbp
    if algorithm not in ["fbp", "gridrec"]:
        _kwargs["num_iter"] = nitr

    sname = os.path.join(args.output_dir, 'proj_{}'.format(args.algorithm))
    print(proj.shape)
    tmp = np.zeros((proj.shape[0], proj.shape[2]))
    tmp[:,:] = proj[:,0,:]
    output_image(tmp, sname + "." + args.format)

    # Reconstruct object.
    with timemory.util.auto_timer(
        "[tomopy.recon(algorithm='{}')]".format(algorithm)):
        print("Starting reconstruction with kwargs={}...".format(_kwargs))
        rec = tomopy.recon(data, theta, **_kwargs)
    print("Completed reconstruction...")

    # Mask each reconstructed slice with a circle.
    rec = tomopy.circ_mask(rec, axis=0, ratio=0.95)

    obj = np.zeros(rec.shape, dtype=rec.dtype)
    label = "{} @ {}".format(algorithm.upper(), h5fname)
    quantify_difference(label, obj, rec)

    return rec

def reconstruct(h5fname, sino, rot_center, binning, algorithm='gridrec'):

    sample_detector_distance = 8        # Propagation distance of the wavefront in cm
    detector_pixel_size_x = 2.247e-4    # Detector pixel size in cm (5x: 1.17e-4, 2X: 2.93e-4)
    monochromator_energy = 24.9         # Energy of incident wave in keV
    alpha = 1e-02                       # Phase retrieval coeff.
    zinger_level = 800                  # Zinger level for projections
    zinger_level_w = 1000               # Zinger level for white

    # Read APS 32-BM raw data.
    proj, flat, dark, theta = dxchange.read_aps_32id(h5fname, sino=sino)
        
    # zinger_removal
    # proj = tomopy.misc.corr.remove_outlier(proj, zinger_level, size=15, axis=0)
    # flat = tomopy.misc.corr.remove_outlier(flat, zinger_level_w, size=15, axis=0)

    # Flat-field correction of raw data.
    ##data = tomopy.normalize(proj, flat, dark, cutoff=0.8)
    data = tomopy.normalize(proj, flat, dark)

    # remove stripes
    #data = tomopy.remove_stripe_fw(data,level=7,wname='sym16',sigma=1,pad=True)

    #data = tomopy.remove_stripe_ti(data, alpha=1.5)
    #data = tomopy.remove_stripe_sf(data, size=150)

    # phase retrieval
    #data = tomopy.prep.phase.retrieve_phase(data,pixel_size=detector_pixel_size_x,dist=sample_detector_distance,energy=monochromator_energy,alpha=alpha,pad=True)

    print("Raw data: ", h5fname)
    print("Center: ", rot_center)

    data = tomopy.minus_log(data)

    data = tomopy.remove_nan(data, val=0.0)
    data = tomopy.remove_neg(data, val=0.00)
    data[np.where(data == np.inf)] = 0.00

    rot_center = rot_center/np.power(2, float(binning))
    data = tomopy.downsample(data, level=binning) 
    data = tomopy.downsample(data, level=binning, axis=1)
    # padding 
    N = data.shape[2]
    data_pad = np.zeros([data.shape[0],data.shape[1],3*N//2],dtype = "float32")
    data_pad[:,:,N//4:5*N//4] = data
    data_pad[:,:,0:N//4] = np.tile(np.reshape(data[:,:,0],[data.shape[0],data.shape[1],1]),(1,1,N//4))
    data_pad[:,:,5*N//4:] = np.tile(np.reshape(data[:,:,-1],[data.shape[0],data.shape[1],1]),(1,1,N//4))

    data = data_pad
    rot_center = rot_center+N//4
 
    # Reconstruct object.
    if algorithm == 'sirtfbp':
        rec = rec_sirtfbp(data, theta, rot_center)
    else:
        rec = tomopy.recon(data, theta, center=rot_center, algorithm=algorithm, filter_name='parzen')
    rec = rec[:,N//4:5*N//4,N//4:5*N//4]
        
    print("Algorithm: ", algorithm)

    # Mask each reconstructed slice with a circle.
#   rec = tomopy.circ_mask(rec, axis=0, ratio=0.95)
    
    return rec

	inputPath = '{}_{:d}{}'.format(fn,y,fileextension)
	tomo[y] = dxchange.reader.read_tiff(inputPath,slc = (sinoused, raysused))

print('loading flat images')
for y in range(0,len(floc)):
	inputPath = '{}{}_{:d}{}'.format(fn,flatextension,floc[y],fileextension)
	flat[y] = dxchange.reader.read_tiff(inputPath,slc = (sinoused, raysused))

print('loading dark images')
for y in range(0,numdrk):
	inputPath = '{}{}_{:d}{}'.format(fn,darkextension,y,fileextension)
	dark[y] = dxchange.reader.read_tiff(inputPath,slc = (sinoused, raysused))	
	
print('normalizing')
tomo = tomo.astype(np.float32)
tomopy.normalize_nf(tomo, flat, dark, floc, out=tomo)

tomopy.minus_log(tomo, out=tomo)

tomo = tomopy.pad(tomo, 2, npad=npad, mode='edge')
rec = tomopy.recon(tomo, tomopy.angles(numangles, angle_offset, angle_offset-angularrange), center=cor+npad, algorithm='gridrec', filter_name='butterworth', filter_par=[.25, 2])
rec = rec[:, npad:-npad, npad:-npad]
rec /= pxsize  # convert reconstructed voxel values from 1/pixel to 1/cm
rec = tomopy.circ_mask(rec, 0)



print('writing recon')
dxchange.write_tiff_stack(rec, fname='rec/'+fn, start=sinoused[0])

# Flat field correct data
logger.info("Flat field correcting data")
proj.scatter(0)
tomopy.normalize(proj.local_arr, flat, dark, ncore=1, out=proj.local_arr)
np.clip(proj.local_arr, 1e-6, 1.0, proj.local_arr)
del flat, dark

# Remove Stripe
# NOTE: we need to change remove_strip_fw to take sinogram order data, since it internally rotates the data
#proj.scatter(1)
#proj.local_arr = tomopy.remove_stripe_fw(proj.local_arr, ncore=1)

# Take the minus log to prepare for reconstruction
#NOTE: no scatter required since minus_log doesn't care about order
tomopy.minus_log(proj.local_arr, ncore=1, out=proj.local_arr)

# Find rotation center per set of sinograms
logger.info("Finding center of rotation")
proj.scatter(1)
# NOTE: center finding doesn't work for my datasets :-(
#center = tomopy.find_center(proj.local_arr, theta, sinogram_order=True)
center = proj.shape[2] // 2
logger.info("Center for sinograms [%d:%d] is %f" % (proj.offset, proj.offset+proj.size, center))

alg = 'gridrec'
logger.info("Reconstructing using: %s" % alg)
# Reconstruct object using algorithm
proj.scatter(1)
rec = tomopy.recon(proj.local_arr,
                   theta,

			
		if useNormalize_nf:
			logging.info('Doing normalize (nearest flats)')
			tomo = tomopy.normalize_nf(projs, flat, dark, floc)
		else:
			logging.info('Doing normalize')
			tomo = tomopy.normalize(projs, flat, dark)
		
		
		#sinofilenametowrite = odirectory+'/rec'+iname[x]+'/'+iname[x]+'sino'
		#dxchange.write_tiff_stack(tomo, fname=sinofilenametowrite, start=sinorange[0]+y*num_sino_per_chunk,axis=1)
		projs = None
		flat = None
					
		logging.info('Doing -log')
		tomo = tomopy.minus_log(np.maximum(tomo,0.000000000001), out=tomo) # in place logarithm 
		
		angularrange = float(gdata['arange'])
		logging.info('angular range: %f', angularrange)


	
		# Use padding to remove halo in reconstruction if present
		if pad_sino:
			npad = int(np.ceil(tomo.shape[2] * np.sqrt(2)) - tomo.shape[2])//2
			tomo = tomopy.pad(tomo, 2, npad=npad, mode='edge')
			cor_rec = cor + npad # account for padding
		else:
			cor_rec = cor
	

    ## Set path (without file suffix) to the micro-CT data to reconstruct.
    fname = 'data_dir/sample'

    ## Import Data.
    proj, flat, dark, theta = dx.exchange.read_aps_13bm(fname, format = 'netcdf4')

    ## Flat-field correction of raw data.
    proj = tp.normalize(proj, flat = flat, dark = dark)

    ## Additional flat-field correction of raw data to negate need to mask.
    proj = tp.normalize_bg(proj, air = 10)

    ## Set rotation center.
    rot_center = tp.find_center_vo(proj)
    print('Center of rotation: ', rot_center)

    tp.minus_log(proj, out = proj)

    # Reconstruct object using Gridrec algorith.
    rec = tp.recon(proj, theta, center = rot_center, sinogram_order = False, algorithm = 'gridrec', filter_name = 'hann')
    rec = tp.remove_nan(rec)

    ## Writing data in netCDF3 .volume.
    ncfile = Dataset('filename.volume', 'w', format = 'NETCDF3_64BIT', clobber = True)
    NX = ncfile.createDimension('NX', rec.shape[2])
    NY = ncfile.createDimension('NY', rec.shape[1])
    NZ = ncfile.createDimension('NZ', rec.shape[0])
    volume = ncfile.createVariable('VOLUME', 'f4', ('NZ','NY','NX'))
    volume[:] = rec
    ncfile.close()

    # Select the sinogram range to reconstruct.
    start = 290
    end = 294

    # Read the APS 5-BM raw data
    proj, flat, dark = dxchange.read_aps_5bm(fname, sino=(start, end))

    # Set data collection angles as equally spaced between 0-180 degrees.
    theta = tomopy.angles(proj.shape[0])

    # Flat-field correction of raw data.
    proj = tomopy.normalize(proj, flat, dark)

    # remove stripes
    proj = tomopy.remove_stripe_fw(proj,level=7,wname='sym16',sigma=1,pad=True)

    # Set rotation center.
    rot_center = proj.shape[2] / 2.0
    print("Center of rotation: ", rot_center)

    proj = tomopy.minus_log(proj)

    # Reconstruct object using Gridrec algorithm.
    rec = tomopy.recon(proj, theta, center=rot_center, algorithm='gridrec')

    # Mask each reconstructed slice with a circle.
    rec = tomopy.circ_mask(rec, axis=0, ratio=0.95)

    # Write data as stack of TIFs.
    dxchange.write_tiff_stack(rec, fname='recon_dir/recon')

def reconstruct(h5fname, sino, rot_center, binning, algorithm='gridrec'):

    sample_detector_distance = 31      # Propagation distance of the wavefront in cm
    detector_pixel_size_x = 1.17e-4    # Detector pixel size in cm (5x: 1.17e-4, 2X: 2.93e-4)
    monochromator_energy = 65    # Energy of incident wave in keV
    # used pink beam

    alpha = 4*1e-4                       # Phase retrieval coeff.
    zinger_level = 800                  # Zinger level for projections
    zinger_level_w = 1000               # Zinger level for white

    # Read APS 2-BM raw data.
    # DIMAX saves 3 files: proj, flat, dark
    # when loading the data set select the proj file (larger size)

    fname = os.path.splitext(h5fname)[0]    
 
    fbase = fname.rsplit('_', 1)[0]
    fnum = fname.rsplit('_', 1)[1]
    fext = os.path.splitext(h5fname)[1]  

    fnum_flat = str("%4.4d" % (int(fnum)+1))
    fnum_dark = str("%4.4d" % (int(fnum)+2))

    fnproj = fbase + '_' + fnum + fext
    fnflat = fbase + '_' + fnum_flat + fext
    fndark = fbase + '_' + fnum_dark + fext
    
    print('proj', fnproj)
    print('flat', fnflat)
    print('dark', fndark)
    # Read APS 2-BM DIMAX raw data.
    proj, dum, dum2, theta = dxchange.read_aps_32id(fnproj, sino=sino)
    dum3, flat, dum4, dum5 = dxchange.read_aps_32id(fnflat, sino=sino)
    #flat, dum3, dum4, dum5 = dxchange.read_aps_32id(fnflat, sino=sino)          
    dum6, dum7, dark, dum8 = dxchange.read_aps_32id(fndark, sino=sino)

    # Flat-field correction of raw data.
    data = tomopy.normalize(proj, flat, dark, cutoff=1.4)

    # remove stripes
    data = tomopy.remove_stripe_fw(data,level=7,wname='sym16',sigma=1,pad=True)        
    
    # zinger_removal
    proj = tomopy.misc.corr.remove_outlier(proj, zinger_level, size=15, axis=0)
    flat = tomopy.misc.corr.remove_outlier(flat, zinger_level_w, size=15, axis=0)

    # Flat-field correction of raw data.
    ##data = tomopy.normalize(proj, flat, dark, cutoff=0.8)
    data = tomopy.normalize(proj, flat, dark)

    # remove stripes
    #data = tomopy.remove_stripe_fw(data,level=7,wname='sym16',sigma=1,pad=True)

    #data = tomopy.remove_stripe_ti(data, alpha=1.5)
    data = tomopy.remove_stripe_sf(data, size=150)

    # phase retrieval
    #data = tomopy.prep.phase.retrieve_phase(data,pixel_size=detector_pixel_size_x,dist=sample_detector_distance,energy=monochromator_energy,alpha=alpha,pad=True)

    print("Raw data: ", h5fname)
    print("Center: ", rot_center)

    data = tomopy.minus_log(data)

    data = tomopy.remove_nan(data, val=0.0)
    data = tomopy.remove_neg(data, val=0.00)
    data[np.where(data == np.inf)] = 0.00

    rot_center = rot_center/np.power(2, float(binning))
    data = tomopy.downsample(data, level=binning) 
    data = tomopy.downsample(data, level=binning, axis=1)

    # padding 
    N = data.shape[2]
    data_pad = np.zeros([data.shape[0],data.shape[1],3*N//2],dtype = "float32")
    data_pad[:,:,N//4:5*N//4] = data
    data_pad[:,:,0:N//4] = np.tile(np.reshape(data[:,:,0],[data.shape[0],data.shape[1],1]),(1,1,N//4))
    data_pad[:,:,5*N//4:] = np.tile(np.reshape(data[:,:,-1],[data.shape[0],data.shape[1],1]),(1,1,N//4))

    data = data_pad
    rot_center = rot_center+N//4

    nframes = 8 
    nproj = 1500
    theta = np.linspace(0, np.pi*nframes, nproj*nframes, endpoint=False)
    rec = np.zeros(
            (nframes, data.shape[1], data.shape[2], data.shape[2]), dtype='float32')
    for time_frame in range(0, nframes):
        rec0 = tomopy.recon(data[time_frame*nproj:(time_frame+1)*nproj], theta[time_frame*nproj:(
               time_frame+1)*nproj], center=rot_center, algorithm='gridrec')
        # Mask each reconstructed slice with a circle.
        rec[time_frame] = tomopy.circ_mask(rec0, axis=0, ratio=0.95)
    rec = rec[:,:,N//4:5*N//4,N//4:5*N//4]

        
    print("Algorithm: ", algorithm)
    
    return rec

def rec_try(h5fname, nsino, rot_center, center_search_width, algorithm, binning, dark_file):
    
    data_shape = get_dx_dims(h5fname, 'data')
    print(data_shape)
    ssino = int(data_shape[1] * nsino)

    center_range = (rot_center-center_search_width, rot_center+center_search_width, 0.5)
    #print(sino,ssino, center_range)
    #print(center_range[0], center_range[1], center_range[2])

    # Select sinogram range to reconstruct
    sino = None
        
    start = ssino
    end = start + 1
    sino = (start, end)

    # Read APS 32-BM raw data.
    proj, flat, dark, theta = dxchange.read_aps_32id(h5fname, sino=sino)
    if dark_file is not None:
        print('Reading white/dark from {}'.format(dark_file))
        proj_, flat, dark, theta_ = dxchange.read_aps_32id(dark_file, sino=sino)
        del proj_, theta_

    print(proj.shape, flat.shape, dark.shape)
        
    # Flat-field correction of raw data.
    data = tomopy.normalize(proj, flat, dark, cutoff=1.4)    

    # remove stripes
    # data = tomopy.remove_stripe_fw(data,level=7,wname='sym16',sigma=1,pad=True)


    print("Raw data: ", h5fname)
    print("Center: ", rot_center)

    data = tomopy.minus_log(data)

    data = tomopy.remove_nan(data, val=0.0)
    data = tomopy.remove_neg(data, val=0.00)
    data[np.where(data == np.inf)] = 0.00

    stack = np.empty((len(np.arange(*center_range)), data_shape[0], data_shape[2]))

    index = 0
    for axis in np.arange(*center_range):
        stack[index] = data[:, 0, :]
        index = index + 1

    # Reconstruct the same slice with a range of centers.
    rec = tomopy.recon(stack, theta, center=np.arange(*center_range), sinogram_order=True, algorithm='gridrec', filter_name='parzen', nchunk=1)

    # Mask each reconstructed slice with a circle.
    rec = tomopy.circ_mask(rec, axis=0, ratio=0.95)

    index = 0
    # Save images to a temporary folder.
    #fname = os.path.dirname(h5fname) + os.sep + 'try_rec/' + path_base_name(h5fname) + os.sep + 'recon_' + os.path.splitext(os.path.basename(h5fname))[0]    
    fname = os.path.dirname(h5fname) + os.sep + 'centers/' + path_base_name(h5fname) + os.sep + 'recon_' + os.path.splitext(os.path.basename(h5fname))[0]    
    for axis in np.arange(*center_range):
        rfname = fname + '_' + str('{0:.2f}'.format(axis) + '.tiff')
        dxchange.write_tiff(rec[index], fname=rfname, overwrite=True)
        index = index + 1

    print("Reconstructions: ", fname)

def preprocess_data(prj, flat, dark, FF_norm=flat_field_norm, remove_rings = remove_rings, medfilt_size=medfilt_size, FF_drift_corr=flat_field_drift_corr, downspling=binning):

    if FF_norm:
        # normalize the prj
        print('\n*** Applying flat field correction:') 
        start_norm_time = time.time()
        prj = tomo