def build_panorama(src_folder, file_grid, shift_grid, frame=0, method='max', method2=None, blend_options={}, blend_options2={},
                   blur=None, color_correction=False, margin=100):

    t00 = time.time()
    root = os.getcwd()
    os.chdir(src_folder)
    cam_size = g_shapes(file_grid[0, 0])
    cam_size = cam_size[1:3]
    img_size = shift_grid[-1, -1] + cam_size
    buff = np.zeros([1, 1])
    last_none = False
    if method2 is None:
        for (y, x), value in np.ndenumerate(file_grid):
            if (value != None and frame < g_shapes(value)[0]):
                prj, flt, drk = read_aps_32id_adaptive(value, proj=(frame, frame + 1))
                prj = tomopy.normalize(prj, flt, drk)
                prj = preprocess(prj, blur=blur)
                t0 = time.time()
                buff = blend(buff, np.squeeze(prj), shift_grid[y, x, :], method=method, color_correction=color_correction, **blend_options)
                print('Rank: {:d}; Frame: {:d}; Pos: ({:d}, {:d}); Method: {:s}; Color Corr.:{:b}; Tile stitched in '
                      '{:.2f} s.'.format(rank, frame, y, x, method, color_correction, time.time()-t0))
                if last_none:
                    buff[margin:, margin:-margin][np.isnan(buff[margin:, margin:-margin])] = 0
                    last_none = False
            else:
                last_none = True
    else:
        for y in range(file_grid.shape[0]):
            temp_grid = file_grid[y:y+1, :]
            temp_shift = np.copy(shift_grid[y:y+1, :, :])
            offset = np.min(temp_shift[:, :, 0])
            temp_shift[:, :, 0] = temp_shift[:, :, 0] - offset
            row_buff = np.zeros([1, 1])
            prj, flt, drk = read_aps_32id_adaptive(temp_grid[0, 0], proj=(frame, frame + 1))
            prj = tomopy.normalize(prj, flt, drk)
            prj = preprocess(prj, blur=blur)
            row_buff, _ = arrange_image(row_buff, np.squeeze(prj), temp_shift[0, 0, :], order=1)
            for x in range(1, temp_grid.shape[1]):
                value = temp_grid[0, x]
                if (value != None and frame < g_shapes(value)[0]):
                    prj, flt, drk = read_aps_32id_adaptive(value, proj=(frame, frame + 1))
                    prj = tomopy.normalize(prj, flt, drk)
                    prj = preprocess(prj, blur=blur)
                    t0 = time.time()
                    row_buff = blend(row_buff, np.squeeze(prj), temp_shift[0, x, :], method=method, color_correction=color_correction, **blend_options)
                    print('Rank: {:d}; Frame: {:d}; Pos: ({:d}, {:d}); Method: {:s}; Color Corr.:{:b}; Tile stitched in '
                          '{:.2f} s.'.format(rank, frame, y, x, method, color_correction, time.time() - t0))
                    if last_none:
                        row_buff[margin:, margin:-margin][np.isnan(row_buff[margin:, margin:-margin])] = 0
                        last_none = False
                else:
                    last_none = True
            t0 = time.time()
            buff = blend(buff, row_buff, [offset, 0], method=method2, color_correction=False, **blend_options2)
            print('Rank: {:d}; Frame: {:d}; Row: {:d}; Row stitched in {:.2f} s.'.format(rank, frame, y, time.time()-t0))
    print('Rank: {:d}; Frame: {:d}; Panorama built in {:.2f} s.'.format(rank, frame, time.time()-t00))
    os.chdir(root)
    return buff

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 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_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 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 find_rotation_axis(h5fname, nsino):
    
    data_size = get_dx_dims(h5fname, 'data')
    ssino = int(data_size[1] * nsino)

    # 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)
        
    # Flat-field correction of raw data
    data = tomopy.normalize(proj, flat, dark, cutoff=1.4)

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

    # find rotation center
    rot_center = tomopy.find_center_vo(data)   

    return rot_center

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 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 save_partial_frames(file_grid, save_folder, prefix, frame=0):
    for (y, x), value in np.ndenumerate(file_grid):
        print(value)
        if (value != None):
            prj, flt, drk = read_aps_32id_adaptive(value, proj=(frame, frame + 1))
            prj = tomopy.normalize(prj, flt, drk)
            prj = preprocess(prj)
            fname = prefix + 'Y' + str(y).zfill(2) + '_X' + str(x).zfill(2)
            dxchange.write_tiff(np.squeeze(prj), fname=os.path.join(save_folder, 'partial_frames', 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

    # 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 load_sino(filename, sino_n, normalize=True):
    print('Loading {:s}, slice {:d}'.format(filename, sino_n))
    sino_n = int(sino_n)
    sino, flt, drk = read_aps_32id_adaptive(filename, sino=(sino_n, sino_n + 1))
    if not normalize:
        flt[:, :, :] = flt.max()
        drk[:, :, :] = 0
    sino = tomopy.normalize(sino, flt, drk)
    # 1st slice of each tile of some samples contains mostly abnormally large values which should be removed.
    if sino.max() > 1e2:
        sino[np.abs(sino) > 1] = 1
    return np.squeeze(sino)

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 write_first_frames(ui):

    root = os.getcwd()
    os.chdir(ui.raw_folder)
    try:
        os.mkdir('first_frames')
    except:
        pass
    for i in ui.filelist:
        ui.boxMetaOut.insert(END, i + '\n')
        prj, flt, drk = dxchange.read_aps_32id(i, proj=(0, 1))
        prj = tomopy.normalize(prj, flt, drk)
        dxchange.write_tiff(prj, os.path.join('first_frames', os.path.splitext(i)[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 = 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 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):

    parser = argparse.ArgumentParser()
    parser.add_argument("fname", help="file name of a single dataset to normalize: /data/sample.h5")

    args = parser.parse_args()

    fname = args.fname

    if os.path.isfile(fname):
        data_shape = get_dx_dims(fname, 'data')

        # Select projgram range to reconstruct.
        proj_start = 0
        proj_end = data_shape[0]

        chunks = 6          # number of projgram chunks to reconstruct
                            # only one chunk at the time is converted
                            # allowing for limited RAM machines to complete a full reconstruction

        nProj_per_chunk = (proj_end - proj_start)/chunks
        print("Normalizing [%d] slices from slice [%d] to [%d] in [%d] chunks of [%d] slices each" % ((proj_end - proj_start), proj_start, proj_end, chunks, nProj_per_chunk))            

        strt = 0
        for iChunk in range(0,chunks):
            print('\n  -- chunk # %i' % (iChunk+1))
            proj_chunk_start = np.int(proj_start + nProj_per_chunk*iChunk)
            proj_chunk_end = np.int(proj_start + nProj_per_chunk*(iChunk+1))
            print('\n  --------> [%i, %i]' % (proj_chunk_start, proj_chunk_end))
                    
            if proj_chunk_end > proj_end: 
                break

            nproj = (int(proj_chunk_start), int(proj_chunk_end))
            # Reconstruct.
            proj, flat, dark, dummy = dxchange.read_aps_32id(fname, proj=nproj)

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

            # Write data as stack of TIFs.
            tifffname = os.path.dirname(fname) + os.sep + os.path.splitext(os.path.basename(fname))[0]+ '_tiff' + os.sep + os.path.splitext(os.path.basename(fname))[0]
            print("Converted files: ", tifffname)
            dxchange.write_tiff_stack(data, fname=tifffname, start=strt)
            strt += nproj[1] - nproj[0]
    else:
        print("File Name does not exist: ", fname)

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