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=0, type=int, default=0, help="index of the first image: 10001 (default 0)")

    args = parser.parse_args()

    top = args.top

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

    print(top)
    # Read the Australian Synchrotron Facility data
    proj, flat, dark = dxchange.read_aps_5bm(top)
#    proj, flat, dark = dxchange.read_aps_5bm(fname, sino=(start, end))

    slider(proj)

    # 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)

    slider(proj)

def main():
    #****************************************************************************
    file_name = '/local/dataraid/databank/dataExchange/tmp/Australian_rank3.h5'
    output_name = '/local/dataraid/databank/dataExchange/tmp/rec/Australian_rank3'    
    sino_start = 290    
    sino_end = 294    

    # Read HDF5 file.
    exchange_rank = 3;
    prj, flat, dark = tomopy.io.exchange.read_aps_32id(file_name, exchange_rank, sino=(sino_start, sino_end))
    theta  = tomopy.angles(prj.shape[0])

    # normalize the data
    prj = tomopy.normalize(prj, flat, dark)

    best_center=1184
    print "Best Center: ", best_center
    calc_center = best_center
    #calc_center = tomopy.find_center(prj, theta, emission=False, ind=0, init=best_center, tol=0.8)
    print "Calculated Center:", calc_center
    
    # reconstruct 
    rec = tomopy.recon(prj, theta, center=calc_center, algorithm='gridrec', emission=False)
    #rec = tomopy.circ_mask(rec, axis=0)
    
    # Write data as stack of TIFs.
    tomopy.io.writer.write_tiff_stack(rec, fname=output_name)
    plt.gray()
    plt.axis('off')
    plt.imshow(rec[0])

def main(argv):
    try:
        opts, args = getopt.getopt(argv,"hc:s:",["core=","sino="])
    except getopt.GetoptError:
        print 'test.py -c <ncore> -s <nsino>'
        sys.exit(2)
    for opt, arg in opts:
        if opt == '-h':
            print 'test.py -c <ncore> -s <nsino>'
            sys.exit()
        elif opt in ("-c", "--core"):
            ncore = int(arg)
        elif opt in ("-s", "--sino"):
            nsino = int(arg)
    file_name = '/local/decarlo/data/proj_10.hdf'
    output_name = './recon/proj10_rec'
    sino_start = 200


    # Read HDF5 file.
    prj, flat, dark = tomopy.io.exchange.read_aps_32id(file_name, sino=(sino_start, sino_start+nsino))

    # Fix flats because sample did not move
    flat = np.full((flat.shape[0], flat.shape[1], flat.shape[2]), 1000)

    # Set angles
    theta  = tomopy.angles(prj.shape[0])

def generate(phantom, args):
    """Return the simulated data for the given phantom."""
    with timemory.util.auto_timer("[tomopy.misc.phantom.{}]".format(phantom)):
        obj = getattr(tomopy.misc.phantom, phantom)(size=args.size)
        obj = tomopy.misc.morph.pad(obj, axis=1, mode='constant')
        obj = tomopy.misc.morph.pad(obj, axis=2, mode='constant')

        if args.partial:
            data_size = obj.shape[0]
            subset = list(args.subset)
            subset.sort()
            nbeg, nend = subset[0], subset[1]
            if nbeg == nend:
                nend += 1
            if not args.no_center:
                ndiv = (nend - nbeg) // 2
                offset = data_size // 2
                nbeg = (offset - ndiv)
                nend = (offset + ndiv)
            print("[partial]> slices = {} ({}, {}) of {}".format(
                nend - nbeg, nbeg, nend, data_size))
            obj = obj[nbeg:nend,:,:]

    with timemory.util.auto_timer("[tomopy.angles]"):
        ang = tomopy.angles(args.angles)

    with timemory.util.auto_timer("[tomopy.project]"):
        prj = tomopy.project(obj, ang)

    print("[dims]> projection = {}, angles = {}, object = {}".format(
        prj.shape, ang.shape, obj.shape))
    return [prj, ang, obj]

def rec_test(file_name, sino_start, sino_end, astra_method, extra_options, num_iter=1):

    print '\n#### Processing '+ file_name
    sino_start = sino_start + 200
    sino_end = sino_start + 2
    print "Test reconstruction of slice [%d]" % sino_start
    # Read HDF5 file.
    prj, flat, dark = tomopy.io.exchange.read_aps_32id(file_name, sino=(sino_start, sino_end))

    # Manage the missing angles:
    theta  = tomopy.angles(prj.shape[0])
    prj = np.concatenate((prj[0:miss_angles[0],:,:], prj[miss_angles[1]+1:-1,:,:]), axis=0)
    theta = np.concatenate((theta[0:miss_angles[0]], theta[miss_angles[1]+1:-1]))

    # normalize the prj
    prj = tomopy.normalize(prj, flat, dark)
    
    # remove ring artefacts
    prjn = tomopy.remove_stripe_fw(prj)

    # reconstruct 
    rec = tomopy.recon(prj[:,::reduce_amount,::reduce_amount], theta, center=float(best_center)/reduce_amount, algorithm=tomopy.astra, options={'proj_type':proj_type,'method':astra_method,'extra_options':extra_options,'num_iter':num_iter}, emission=False)
        
    # Write data as stack of TIFs.
    tomopy.io.writer.write_tiff_stack(rec, fname=output_name)

    print "Slice saved as [%s_00000.tiff]" % output_name

def rec_test(file_name, sino_start, sino_end):

    print "\n#### Processing " + file_name
    sino_start = sino_start + 200
    sino_end = sino_start + 2
    print "Test reconstruction of slice [%d]" % sino_start
    # Read HDF5 file.
    prj, flat, dark = tomopy.io.exchange.read_aps_32id(file_name, sino=(sino_start, sino_end))

    # Manage the missing angles:
    theta = tomopy.angles(prj.shape[0])
    prj = np.concatenate((prj[0 : miss_angles[0], :, :], prj[miss_angles[1] + 1 : -1, :, :]), axis=0)
    theta = np.concatenate((theta[0 : miss_angles[0]], theta[miss_angles[1] + 1 : -1]))

    # normalize the prj
    prj = tomopy.normalize(prj, flat, dark)

    # reconstruct
    rec = tomopy.recon(prj, theta, center=best_center, algorithm="gridrec", emission=False)

    # Write data as stack of TIFs.
    tomopy.io.writer.write_tiff_stack(rec, fname=output_name)

    print "Slice saved as [%s_00000.tiff]" % output_name
    # show the reconstructed slice
    pl.gray()
    pl.axis("off")
    pl.imshow(rec[0])

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 generate(phantom="shepp3d", nsize=512, nangles=360):

    with timemory.util.auto_timer("[tomopy.misc.phantom.{}]".format(phantom)):
        obj = getattr(tomopy.misc.phantom, phantom)(size=nsize)
    with timemory.util.auto_timer("[tomopy.angles]"):
        ang = tomopy.angles(nangles)
    with timemory.util.auto_timer("[tomopy.project]"):
        prj = tomopy.project(obj, ang)

    return [prj, ang, obj]

def recon_slice(row_sino, center_pos, sinogram_order=False, algorithm=None,
        init_recon=None, ncore=None, nchunk=None, **kwargs):
    t = time.time()
    ang = tomopy.angles(row_sino.shape[0])
    print(row_sino.shape)
    row_sino = row_sino.astype('float32')
    # row_sino = tomopy.normalize_bg(row_sino) # WARNING: normalize_bg can unpredicatably give bad results for some slices
    row_sino = tomopy.remove_stripe_ti(row_sino, alpha=4)
    rec = tomopy.recon(row_sino, ang, center=center_pos, sinogram_order=sinogram_order, algorithm=algorithm,
        init_recon=init_recon, ncore=ncore, nchunk=nchunk, **kwargs)

    print('recon:           ' + str(time.time() - t))
    return rec

def main(argv):
    try:
        opts, args = getopt.getopt(argv,"hc:s:",["core=","sino="])
    except getopt.GetoptError:
        print 'test.py -c <ncore> -s <nsino>'
        sys.exit(2)
    for opt, arg in opts:
        if opt == '-h':
            print 'test.py -c <ncore> -s <nsino>'
            sys.exit()
        elif opt in ("-c", "--core"):
            ncore = int(arg)
        elif opt in ("-s", "--sino"):
            nsino = int(arg)

    # **********************************************
    #file_name = '/local/decarlo/data/proj_10.hdf'
    #output_name = './recon/proj10_rec'
    #sino_start = 0
    #sino_end = 2048
    # **********************************************
    file_name = '/local/decarlo/data/Hornby_APS_2011.h5'
    output_name = './recon/Hornby_APS_2011_'
    best_center=1024
    sino_start = 0
    sino_end = 1792
    # **********************************************

    step_00 = time.time()
    step_02_delta_total = 0
    
    count = 0
    while (sino_start <= (sino_end - nsino)):
        # Read HDF5 file.
        prj, flat, dark = tomopy.io.exchange.read_aps_32id(file_name, sino=(sino_start, sino_start+nsino))

        # Fix flats because sample did not move
        flat = np.full((flat.shape[0], flat.shape[1], flat.shape[2]), 1000)

        # Set angles
        theta  = tomopy.angles(prj.shape[0])

        # normalize the prj
        prj = tomopy.normalize(prj, flat, dark)

        best_center = 1298
        step_01 = time.time()

        # reconstruct 
        rec = tomopy.recon(prj, theta, center=best_center, algorithm='gridrec', emission=False, ncore = ncore)

def recon_hdf5_mpi(src_fanme, dest_folder, sino_range, sino_step, center_vec, shift_grid, dtype='float32',
               algorithm='gridrec', tolerance=1, save_sino=False, sino_blur=None, **kwargs):
    """
    Reconstruct a single tile, or fused HDF5 created using util/total_fusion. MPI supported.
    """

    raise DeprecationWarning

    if rank == 0:
        if not os.path.exists(dest_folder):
            os.mkdir(dest_folder)
    sino_ini = int(sino_range[0])
    sino_end = int(sino_range[1])
    f = h5py.File(src_fanme)
    dset = f['exchange/data']
    full_shape = dset.shape
    theta = tomopy.angles(full_shape[0])
    center_vec = np.asarray(center_vec)
    sino_ls = np.arange(sino_ini, sino_end, sino_step, dtype='int')
    grid_bins = np.ceil(shift_grid[:, 0, 0])

    t0 = time.time()
    alloc_set = allocate_mpi_subsets(sino_ls.size, size, task_list=sino_ls)
    for slice in alloc_set[rank]:
        print('    Rank {:d}: reconstructing {:d}'.format(rank, slice))
        grid_line = np.digitize(slice, grid_bins)
        grid_line = grid_line - 1
        center = center_vec[grid_line]
        data = dset[:, slice, :]
        if sino_blur is not None:
            data = gaussian_filter(data, sino_blur)
        data = data.reshape([full_shape[0], 1, full_shape[2]])
        data[np.isnan(data)] = 0
        data = data.astype('float32')
        if save_sino:
            dxchange.write_tiff(data[:, slice, :], fname=os.path.join(dest_folder, 'sino/recon_{:05d}_{:d}.tiff').format(slice, center))
        # data = tomopy.remove_stripe_ti(data)
        rec = tomopy.recon(data, theta, center=center, algorithm=algorithm, **kwargs)
        # rec = tomopy.remove_ring(rec)
        rec = tomopy.remove_outlier(rec, tolerance)
        rec = tomopy.circ_mask(rec, axis=0, ratio=0.95)
        dxchange.write_tiff(rec, fname='{:s}/recon/recon_{:05d}_{:d}'.format(dest_folder, slice, center), dtype=dtype)

    print('Rank {:d} finished in {:.2f} s.'.format(rank, time.time()-t0))
    return

def rec_full(file_name, sino_start, sino_end):

    print "\n#### Processing " + file_name

    chunks = 10  # number of data chunks for the reconstruction

    nSino_per_chunk = (sino_end - sino_start) / chunks
    print "Reconstructing [%d] slices from slice [%d] to [%d] in [%d] chunks of [%d] slices each" % (
        (sino_end - sino_start),
        sino_start,
        sino_end,
        chunks,
        nSino_per_chunk,
    )

    for iChunk in range(0, chunks):
        print "\n  -- chunk # %i" % (iChunk + 1)
        sino_chunk_start = sino_start + nSino_per_chunk * iChunk
        sino_chunk_end = sino_start + nSino_per_chunk * (iChunk + 1)
        print "\n  --------> [%i, %i]" % (sino_chunk_start, sino_chunk_end)

        if sino_chunk_end > sino_end:
            break

        # Read HDF5 file.
        prj, flat, dark = tomopy.io.exchange.read_aps_32id(file_name, sino=(sino_chunk_start, sino_chunk_end))

        # Manage the missing angles:
        theta = tomopy.angles(prj.shape[0])
        prj = np.concatenate((prj[0 : miss_angles[0], :, :], prj[miss_angles[1] + 1 : -1, :, :]), axis=0)
        theta = np.concatenate((theta[0 : miss_angles[0]], theta[miss_angles[1] + 1 : -1]))

        # normalize the prj
        prj = tomopy.normalize(prj, flat, dark)

        # reconstruct
        rec = tomopy.recon(prj, theta, center=best_center, algorithm="gridrec", emission=False)

        print output_name

        # Write data as stack of TIFs.
        tomopy.io.writer.write_tiff_stack(rec, fname=output_name, start=sino_chunk_start)

def rec_full(file_name, sino_start, sino_end, astra_method, extra_options, num_iter=1):

    print '\n#### Processing '+ file_name

    chunks = 10 # number of data chunks for the reconstruction

    nSino_per_chunk = (sino_end - sino_start)/chunks
    print "Reconstructing [%d] slices from slice [%d] to [%d] in [%d] chunks of [%d] slices each" % ((sino_end - sino_start), sino_start, sino_end, chunks, nSino_per_chunk)
    strt = 0
    for iChunk in range(0,chunks):
        print '\n  -- chunk # %i' % (iChunk+1)
        sino_chunk_start = sino_start + nSino_per_chunk*iChunk 
        sino_chunk_end = sino_start + nSino_per_chunk*(iChunk+1)
        print '\n  --------> [%i, %i]' % (sino_chunk_start, sino_chunk_end)
        
        if sino_chunk_end > sino_end: 
            break
                
        # Read HDF5 file.
        prj, flat, dark = tomopy.io.exchange.read_aps_32id(file_name, sino=(sino_chunk_start, sino_chunk_end))

        # Manage the missing angles:
        theta  = tomopy.angles(prj.shape[0])
        prj = np.concatenate((prj[0:miss_angles[0],:,:], prj[miss_angles[1]+1:-1,:,:]), axis=0)
        theta = np.concatenate((theta[0:miss_angles[0]], theta[miss_angles[1]+1:-1]))

        # normalize the prj
        prj = tomopy.normalize(prj, flat, dark)
        
        # remove ring artefacts
        prj = tomopy.remove_stripe_fw(prj)

        # reconstruct 
        rec = tomopy.recon(prj[:,::reduce_amount,::reduce_amount], theta, center=float(best_center)/reduce_amount, algorithm=tomopy.astra, options={'proj_type':proj_type,'method':astra_method,'extra_options':extra_options,'num_iter':num_iter}, emission=False)
        
        print output_name

        # Write data as stack of TIFs.
        tomopy.io.writer.write_tiff_stack(rec, fname=output_name, start=strt)
        strt += prj[:,::reduce_amount,:].shape[1]

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!")

# -*- coding: utf-8 -*-
"""
Created on Fri Nov  6 14:41:29 2015

@author: lbluque
"""

import tomopy
import numpy as np

proj_num = 45
input_path = '/home/lbluque/TestDataSets/PhantomSets/projections832small' + str(proj_num)
basename = 'projections832small' + str(proj_num)
filename = basename + '_0000.tif'
input_name = input_path + '/' + filename
ind = range(proj_num)
digits = 4

tomo = tomopy.read_tiff_stack(input_name, ind, digits)
tomo_padded = np.pad(tomo, ((0, 0), (0,0), (20,20)), 'constant', constant_values=0)
center = (tomo.shape[2] - 1)/2.0
theta = tomopy.angles(tomo.shape[0], 90, 270 - 180/proj_num)

rec = tomopy.recon(tomo_padded, theta, center=center, algorithm='gridrec', emission=False)

output_path = '/home/lbluque/TestRecons/PhantomRecons/' + basename + '_reverse/'
outname = output_path + basename
tomopy.write_tiff_stack(rec, fname=outname, digit=digits)

    step = 100;

    index = dataProj.shape[1] / step+1

    for i in range(0, index):
        sliceStart = i*step;
        sliceEnd = (i+1)*step;
        if(sliceStart >= dataProj.shape[1]):
            sys.exit(0);
        if(sliceEnd > dataProj.shape[1]):
            sliceEnd = dataProj.shape[1];

        proj = dataProj[:,sliceStart:sliceEnd,:]

        start1 = timeit.default_timer()
        theta = tomopy.angles(proj.shape[0], 11, 168)

        stop1 = timeit.default_timer()
        print("end angles", (stop1 - start1))

        start2 = timeit.default_timer()
        # # Flat-field correction of raw data.
        proj = tomopy.normalize(proj, flat, dark)
        stop2 = timeit.default_timer()
        print("end normalize", (stop2 - start2))

        start3 = timeit.default_timer()
        # # Find rotation center.
        rot_center = tomopy.find_center(proj, theta, emission=False, ind=0, init=1024, tol=0.5)
        stop3 = timeit.default_timer()
        print("end find_center", (stop3 - start3))

def recon_hdf5(src_fanme, dest_folder, sino_range, sino_step, shift_grid, center_vec=None, center_eq=None, dtype='float32',
               algorithm='gridrec', tolerance=1, chunk_size=20, save_sino=False, sino_blur=None, flattened_radius=120,
               mode='180', test_mode=False, phase_retrieval=None, ring_removal=True, **kwargs):
    """
    center_eq: a and b parameters in fitted center position equation center = a*slice + b.
    """

    if not os.path.exists(dest_folder):
        try:
            os.mkdir(dest_folder)
        except:
            pass
    sino_ini = int(sino_range[0])
    sino_end = int(sino_range[1])
    sino_ls_all = np.arange(sino_ini, sino_end, sino_step, dtype='int')
    alloc_set = allocate_mpi_subsets(sino_ls_all.size, size, task_list=sino_ls_all)
    sino_ls = alloc_set[rank]

    # prepare metadata
    f = h5py.File(src_fanme)
    dset = f['exchange/data']
    full_shape = dset.shape
    theta = tomopy.angles(full_shape[0])
    if center_eq is not None:
        a, b = center_eq
        center_ls = sino_ls * a + b
        center_ls = np.round(center_ls)
        for iblock in range(int(sino_ls.size/chunk_size)+1):
            print('Beginning block {:d}.'.format(iblock))
            t0 = time.time()
            istart = iblock*chunk_size
            iend = np.min([(iblock+1)*chunk_size, sino_ls.size])
            fstart = sino_ls[istart]
            fend = sino_ls[iend]
            center = center_ls[istart:iend]
            data = dset[:, fstart:fend:sino_step, :]
            data[np.isnan(data)] = 0
            data = data.astype('float32')
            data = tomopy.remove_stripe_ti(data, alpha=4)
            if sino_blur is not None:
                for i in range(data.shape[1]):
                    data[:, i, :] = gaussian_filter(data[:, i, :], sino_blur)
            rec = tomopy.recon(data, theta, center=center, algorithm=algorithm, **kwargs)
            rec = tomopy.remove_ring(rec)
            rec = tomopy.remove_outlier(rec, tolerance)
            rec = tomopy.circ_mask(rec, axis=0, ratio=0.95)
            for i in range(rec.shape[0]):
                slice = fstart + i*sino_step
                dxchange.write_tiff(rec[i, :, :], fname=os.path.join(dest_folder, 'recon/recon_{:05d}_{:05d}.tiff').format(slice, sino_ini))
                if save_sino:
                    dxchange.write_tiff(data[:, i, :], fname=os.path.join(dest_folder, 'sino/recon_{:05d}_{:d}.tiff').format(slice, int(center[i])))
            iblock += 1
            print('Block {:d} finished in {:.2f} s.'.format(iblock, time.time()-t0))
    else:
        # divide chunks
        grid_bins = np.append(np.ceil(shift_grid[:, 0, 0]), full_shape[1])
        chunks = []
        center_ls = []
        istart = 0
        counter = 0
        # irow should be 0 for slice 0
        irow = np.searchsorted(grid_bins, sino_ls[0], side='right')-1

        for i in range(sino_ls.size):
            counter += 1
            sino_next = i+1 if i != sino_ls.size-1 else i
            if counter >= chunk_size or sino_ls[sino_next] >= grid_bins[irow+1] or sino_next == i:
                iend = i+1
                chunks.append((istart, iend))
                istart = iend
                center_ls.append(center_vec[irow])
                if sino_ls[sino_next] >= grid_bins[irow+1]:
                    irow += 1
                counter = 0

        # reconstruct chunks
        iblock = 1
        for (istart, iend), center in izip(chunks, center_ls):
            print('Beginning block {:d}.'.format(iblock))
            t0 = time.time()
            fstart = sino_ls[istart]
            fend = sino_ls[iend-1]
            print('Reading data...')
            data = dset[:, fstart:fend+1:sino_step, :]
            if mode == '360':
                overlap = 2 * (dset.shape[2] - center)
                data = tomosaic.morph.sino_360_to_180(data, overlap=overlap, rotation='right')
                theta = tomopy.angles(data.shape[0])
            data[np.isnan(data)] = 0
            data = data.astype('float32')
            if sino_blur is not None:
                for i in range(data.shape[1]):
                    data[:, i, :] = gaussian_filter(data[:, i, :], sino_blur)
            if ring_removal:
                data = tomopy.remove_stripe_ti(data, alpha=4)
                if phase_retrieval:
                    data = tomopy.retrieve_phase(data, kwargs['pixel_size'], kwargs['dist'], kwargs['energy'],
                                                 kwargs['alpha'])
                rec0 = tomopy.recon(data, theta, center=center, algorithm=algorithm, **kwargs)
                rec = tomopy.remove_ring(np.copy(rec0))
#.........这里部分代码省略.........

def reconstruct(filename,inputPath="", outputPath="", COR=COR, doOutliers=doOutliers, outlier_diff=outlier_diff, outlier_size=outlier_size, doFWringremoval=doFWringremoval, ringSigma=ringSigma,ringLevel=ringLevel, ringWavelet=ringWavelet,pad_sino=pad_sino,  doPhaseRetrieval=doPhaseRetrieval, propagation_dist=propagation_dist, kev=kev,alphaReg=alphaReg, butterworthpars=butterworthpars, doPolarRing=doPolarRing,Rarc=Rarc, Rmaxwidth=Rmaxwidth, Rtmax=Rtmax, Rthr=Rthr, Rtmin=Rtmin, useAutoCOR=useAutoCOR, use360to180=use360to180, num_substacks=num_substacks,recon_slice=recon_slice):

	# Convert filename to list type if only one file name is given
	if type(filename) != list:
		filename=[filename]

	# If useAutoCor == true, a list of COR will be automatically calculated for all files
	# If a list of COR is given, only entries with boolean False will use automatic COR calculation
	if useAutoCOR==True or (len(COR) != len(filename)):
		logging.info('using auto COR for all input files')
		COR = [False]*len(filename)

	for x in range(len(filename)):
		logging.info('opening data set, checking metadata')

		fdata, gdata = read_als_832h5_metadata(inputPath[x]+filename[x]+'.h5')
		pxsize = float(gdata['pxsize'])/10.0 # convert from metadata (mm) to this script (cm)
		numslices = int(gdata['nslices'])

		# recon_slice == True, only center slice will be reconstructed
		# if integer is given, a specific 		
		if recon_slice != False:
			if (type(recon_slice) == int) and (recon_slice <= numslices):
				sinorange [recon_slice-1, recon_slice]
			else:
				sinorange = [numslices//2-1, numslices//2]
		else:
			sinorange = [0, numslices]

		# Calculate number of substacks (chunks)
		substacks = num_substacks #(sinorange[1]-sinorange[0]-1)//num_sino_per_substack+1

		if (sinorange[1]-sinorange[0]) >= substacks:
			num_sino_per_substack = (sinorange[1]-sinorange[0])//num_substacks
		else:
			num_sino_per_substack = 1

	
		firstcor, lastcor = 0, int(gdata['nangles'])-1
		projs, flat, dark, floc = dxchange.read_als_832h5(inputPath[x]+filename[x]+'.h5', ind_tomo=(firstcor, lastcor))
		projs = tomopy.normalize_nf(projs, flat, dark, floc)
		autocor = tomopy.find_center_pc(projs[0], projs[1], tol=0.25)


		if (type(COR[x]) == bool) or (COR[x]<0) or (COR[x]=='auto'):
			firstcor, lastcor = 0, int(gdata['nangles'])-1
			projs, flat, dark, floc = dxchange.read_als_832h5(inputPath[x]+filename[x]+'.h5', ind_tomo=(firstcor, lastcor))
			projs = tomopy.normalize_nf(projs, flat, dark, floc)
			cor = tomopy.find_center_pc(projs[0], projs[1], tol=0.25)
		else:
			cor = COR[x]

		logging.info('Dataset %s, has %d total slices, reconstructing slices %d through %d in %d substack(s), using COR: %f',filename[x], int(gdata['nslices']), sinorange[0], sinorange[1]-1, substacks, cor)
		
		for y in range(0, substacks):
			logging.info('Starting dataset %s (%d of %d), substack %d of %d',filename[x], x+1, len(filename), y+1, substacks)

			logging.info('Reading sinograms...')
			projs, flat, dark, floc = dxchange.read_als_832h5(inputPath[x]+filename[x]+'.h5', sino=(sinorange[0]+y*num_sino_per_substack, sinorange[0]+(y+1)*num_sino_per_substack, 1)) 

			logging.info('Doing remove outliers, norm (nearest flats), and -log...')
			if doOutliers:
				projs = tomopy.remove_outlier(projs, outlier_diff, size=outlier_size, axis=0)
				flat = tomopy.remove_outlier(flat, outlier_diff, size=outlier_size, axis=0)
			tomo = tomopy.normalize_nf(projs, flat, dark, floc)
			tomo = tomopy.minus_log(tomo, out=tomo) # in place logarithm 
		
			# 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
		
			if doFWringremoval:
				logging.info('Doing ring (Fourier-wavelet) function...')
				tomo = tomopy.remove_stripe_fw(tomo, sigma=ringSigma, level=ringLevel, pad=True, wname=ringWavelet)		

			if doPhaseRetrieval:
				logging.info('Doing Phase retrieval...')
				#tomo = tomopy.retrieve_phase(tomo, pixel_size=pxsize, dist=propagation_dist, energy=kev, alpha=alphaReg, pad=True)	
				tomo = tomopy.retrieve_phase(tomo, pixel_size=pxsize, dist=propagation_dist, energy=kev, alpha=alphaReg, pad=True)		

			logging.info('Doing recon (gridrec) function and scaling/masking, with cor %f...',cor_rec)
			rec = tomopy.recon(tomo, tomopy.angles(tomo.shape[0], 270, 90), center=cor_rec, algorithm='gridrec', filter_name='butterworth', filter_par=butterworthpars)
			#rec = tomopy.recon(tomo, tomopy.angles(tomo.shape[0], 180+angularrange/2, 180-angularrange/2), center=cor_rec, algorithm='gridrec', filter_name='butterworth', filter_par=butterworthpars)		
			rec /= pxsize  # intensity values in cm^-1
			if pad_sino:
				rec = tomopy.circ_mask(rec[:, npad:-npad, npad:-npad], 0)
			else:
				rec = tomopy.circ_mask(rec, 0, ratio=1.0, val=0.0)
			
			if doPolarRing:
				logging.info('Doing ring (polar mean filter) function...')
				rec = tomopy.remove_ring(rec, theta_min=Rarc, rwidth=Rmaxwidth, thresh_max=Rtmax, thresh=Rthr, thresh_min=Rtmin)

			logging.info('Writing reconstruction slices to %s', filename[x])
			#dxchange.write_tiff_stack(rec, fname=outputPath+'alpha'+str(alphaReg)+'/rec'+filename[x]+'/rec'+filename[x], start=sinorange[0]+y*num_sino_per_substack)
			dxchange.write_tiff_stack(rec, fname=outputPath + 'recon_'+filename[x]+'/recon_'+filename[x], start=sinorange[0]+y*num_sino_per_substack)
#.........这里部分代码省略.........