Python tomopy.circ_mask函数代码示例

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

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

#.........这里部分代码省略.........
    parser.add_argument("-lf", "--log-file", type=str, help="log file name", default="fast-tomopy.log")

    args = parser.parse_args()

    fh = logging.FileHandler(args.log_file)
    fh.setLevel(logging.INFO)
    fh.setFormatter(formatter)
    logger.addHandler(fh)

    if os.path.isdir(os.path.dirname(args.output_file)) is False:
        raise IOError(2, "Directory of output file does not exist", args.output_file)

    # Read file metadata
    logger.info("Reading input file metadata")
    fdata, gdata = read_als_832h5_metadata(args.input)
    proj_total = int(gdata["nangles"])
    last = proj_total - 1
    sino_total = int(gdata["nslices"])
    ray_total = int(gdata["nrays"])
    px_size = float(gdata["pxsize"]) / 10  # cm

    # Set parameters for sinograms to read
    step = sino_total // (args.sino_num + 2)
    start = step
    end = step * (args.sino_num + 1)
    sino = (start, end, step)

    # Read full first and last projection to determine center of rotation
    if args.center is None:
        logger.info("Reading full first and last projection for COR")
        first_last, flats, darks, floc = tomopy.read_als_832h5(args.input, ind_tomo=(0, last))
        first_last = tomopy.normalize(first_last, flats, darks)
        args.center = tomopy.find_center_pc(first_last[0, :, :], first_last[1, :, :], tol=0.1)
        logger.info("Detected center: %f", args.center)

    # Read and normalize raw sinograms
    logger.info("Reading raw data")
    tomo, flats, darks, floc = tomopy.read_als_832h5(args.input, sino=sino)
    logger.info("Normalizing raw data")
    tomo = tomopy.normalize_nf(tomo, flats, darks, floc)

    # Remove stripes from sinograms (remove rings)
    logger.info("Preprocessing normalized data")
    if args.ring_remove == "Tomopy-FW":
        logger.info("Removing stripes from sinograms with %s", args.ring_remove)
        tomo = tomopy.remove_stripe_fw(tomo)
    elif args.ring_remove == "Tomopy-T":
        logger.info("Removing stripes from sinograms with %s", args.ring_remove)
        tomo = tomopy.remove_stripe_ti(tomo)

    # Pad sinograms with edge values
    npad = int(np.ceil(ray_total * np.sqrt(2)) - ray_total) // 2
    tomo = tomopy.pad(tomo, 2, npad=npad, mode="edge")
    args.center += npad  # account for padding

    filter_name = np.array(args.filter_name, dtype=(str, 16))
    theta = tomopy.angles(proj_total, 270, 90)

    logger.info("Reconstructing normalized data")
    # Reconstruct sinograms
    # rec = tomopy.minus_log(tomo, out=tomo)
    rec = tomopy.recon(
        tomo, theta, center=args.center, emission=False, algorithm=args.algorithm, filter_name=filter_name
    )
    rec = tomopy.circ_mask(rec[:, npad:-npad, npad:-npad], 0)
    rec = rec / px_size

    # Remove rings from reconstruction
    if args.ring_remove == "Octopus":
        logger.info("Removing rings from reconstructions with %s", args.ring_remove)
        thresh = float(gdata["ring_threshold"])
        thresh_max = float(gdata["upp_ring_value"])
        thresh_min = float(gdata["low_ring_value"])
        theta_min = int(gdata["max_arc_length"])
        rwidth = int(gdata["max_ring_size"])
        rec = tomopy.remove_rings(
            rec,
            center_x=args.center,
            thresh=thresh,
            thresh_max=thresh_max,
            thresh_min=thresh_min,
            theta_min=theta_min,
            rwidth=rwidth,
        )

    # Write reconstruction data to new hdf5 file
    fdata["stage"] = "fast-tomopy"
    fdata["stage_flow"] = "/raw/" + fdata["stage"]
    fdata["stage_version"] = "fast-tomopy-0.1"
    # Generate a new uuid based on host ID and current time
    fdata["uuid"] = str(uuid.uuid1())

    gdata["Reconstruction_Type"] = "tomopy-gridrec"
    gdata["ring_removal_method"] = args.ring_remove
    gdata["rfilter"] = args.filter_name

    logger.info("Writing reconstructed data to h5 file")
    write_als_832h5(rec, args.input, fdata, gdata, args.output_file, step)

    return

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 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))
#.........这里部分代码省略.........

import tomopy
import dxchange

if __name__ == '__main__':

    # Set path to the micro-CT data to reconstruct.
    fname = '../../../tomopy/data/tooth.h5'

    # Select the sinogram range to reconstruct.
    start = 0
    end = 2

    # Read the APS 2-BM 0r 32-ID raw data.
    proj, flat, dark = dxchange.read_aps_32id(fname, sino=(start, end))

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

    # Set data collection angles as equally spaced between 0-180 degrees.
    proj = tomopy.normalize(proj, flat, dark)

    # Set data collection angles as equally spaced between 0-180 degrees.
    rot_center = tomopy.find_center(proj, theta, init=290, ind=0, tol=0.5)

    tomopy.minus_log(proj)

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

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

def main(arg):

    parser = argparse.ArgumentParser()
    parser.add_argument("fname", help="Full file name: /data/fname.raw")
    parser.add_argument("--start", nargs='?', type=int, default=0, help="First image to read")
    parser.add_argument("--nimg", nargs='?', type=int, default=1, help="Number of images to read")
    parser.add_argument("--ndark", nargs='?', type=int, default=10, help="Number of dark images")
    parser.add_argument("--nflat", nargs='?', type=int, default=10, help="Number of white images")

    args = parser.parse_args()

    fname = args.fname
    start = args.start
    end = args.start + args.nimg

    nflat, ndark, nimg, height, width = read_adimec_header(fname)
    print("Image Size:", width, height)
    print("Dataset metadata (nflat, ndark, nimg:", nflat, ndark, nimg)

    # override nflat and ndark from header with the passed parameter
    # comment the two lines below if the meta data in the binary 
    # file for nflat and ndark is correct
    nflat = args.nflat
    ndark = args.ndark

    proj = read_adimec_stack(fname, img=(start, end))
    print("Projection:", proj.shape)
    # slider(proj)

    flat = read_adimec_stack(fname, img=(nimg-ndark-nflat, nimg-ndark))
    print("Flat:", flat.shape)
    # slider(flat)

    dark = read_adimec_stack(fname, img=(nimg-ndark, nimg))
    print("Dark:", dark.shape)
    # slider(dark)

    nproj = tomopy.normalize(proj, flat, dark)
    print("Normalized projection:", nproj.shape)
    # slider(proj)

    
    proj = nproj[:,100:110, :]
    print("Sino chunk:", proj.shape)
    slider(proj)
    
    theta = tomopy.angles(proj.shape[0])
    print(theta.shape)

    proj = tomopy.minus_log(proj)

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

    rot_center = 1280
    # 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')

#.........这里部分代码省略.........
                elif func_name == 'normalize':
                    tomo = tomo.astype(np.float32,copy=False)
                    tomopy.normalize(tomo, flat, dark, out=tomo)
                elif func_name == 'minus_log':
                    mx = np.float32(0.00000000000000000001)
                    ne.evaluate('where(tomo>mx, tomo, mx)', out=tomo)
                    tomopy.minus_log(tomo, out=tomo)
                elif func_name == 'beam_hardening':
                    loc_dict = {'a{}'.format(i):np.float32(val) for i,val in enumerate(BeamHardeningCoefficients)}
                    tomo = ne.evaluate('a0 + a1*tomo + a2*tomo**2 + a3*tomo**3 + a4*tomo**4 + a5*tomo**5', local_dict=loc_dict, out=tomo)
                elif func_name == 'remove_stripe_fw':
                    tomo = tomopy.remove_stripe_fw(tomo, sigma=ringSigma, level=ringLevel, pad=True, wname=ringWavelet)
                elif func_name == 'remove_stripe_ti':
                    tomo = tomopy.remove_stripe_ti(tomo, nblock=ringNBlock, alpha=ringAlpha)
                elif func_name == 'remove_stripe_sf':
                    tomo = tomopy.remove_stripe_sf(tomo, size=ringSize)
                elif func_name == 'correcttilt':
                    if tiltcenter_slice is None:
                        tiltcenter_slice = numslices/2.
                    if tiltcenter_det is None:
                        tiltcenter_det = tomo.shape[2]/2
                    new_center = tiltcenter_slice - 0.5 - sinoused[0]
                    center_det = tiltcenter_det - 0.5
                    
                    #add padding of 10 pixels, to be unpadded right after tilt correction. This makes the tilted image not have zeros at certain edges, which matters in cases where sample is bigger than the field of view. For the small amounts we are generally tilting the images, 10 pixels is sufficient.
#                    tomo = tomopy.pad(tomo, 2, npad=10, mode='edge')
#                    center_det = center_det + 10
                    
                    cntr = (center_det, new_center)
                    for b in range(tomo.shape[0]):
                        tomo[b] = st.rotate(tomo[b], correcttilt, center=cntr, preserve_range=True, order=1, mode='edge', clip=True) #center=None means image is rotated around its center; order=1 is default, order of spline interpolation
#                    tomo = tomo[:, :, 10:-10]    
                        
                elif func_name == 'do_360_to_180':
                    
                    # Keep values around for processing the next chunk in the list
                    keepvalues = [angularrange, numangles, projused, num_proj_per_chunk, numprojchunks, numprojused, numrays, anglelist]
                    
                    #why -.5 on one and not on the other?
                    if tomo.shape[0]%2>0:
                        tomo = sino_360_to_180(tomo[0:-1,:,:], overlap=int(np.round((tomo.shape[2]-cor-.5))*2), rotation='right')
                        angularrange = angularrange/2 - angularrange/(tomo.shape[0]-1)
                    else:
                        tomo = sino_360_to_180(tomo[:,:,:], overlap=int(np.round((tomo.shape[2]-cor))*2), rotation='right')
                        angularrange = angularrange/2
                    numangles = int(numangles/2)
                    projused = (0,numangles-1,1)
                    num_proj_per_chunk = np.minimum(chunk_proj,projused[1]-projused[0])
                    numprojchunks = (projused[1]-projused[0]-1)//num_proj_per_chunk+1
                    numprojused = (projused[1]-projused[0])//projused[2]
                    numrays = tomo.shape[2]
                    
                    anglelist = anglelist[:numangles]
                
                elif func_name == 'phase_retrieval':
                    tomo = tomopy.retrieve_phase(tomo, pixel_size=pxsize, dist=propagation_dist, energy=kev, alpha=alphaReg, pad=True)
                
                elif func_name == 'translation_correction':
                    tomo = linear_translation_correction(tomo,dx=xshift,dy=yshift,interpolation=False):
                    
                elif func_name == 'recon_mask':
                    tomo = tomopy.pad(tomo, 2, npad=npad, mode='edge')

                    if projIgnoreList is not None:
                        for badproj in projIgnoreList:
                            tomo[badproj] = 0

                    rec = tomopy.recon(tomo, anglelist, center=cor+npad, algorithm='gridrec', filter_name='butterworth', filter_par=[butterworth_cutoff, butterworth_order])
                    rec = rec[:, npad:-npad, npad:-npad]
                    rec /= pxsize  # convert reconstructed voxel values from 1/pixel to 1/cm
                    rec = tomopy.circ_mask(rec, 0)
                elif func_name == 'polar_ring':
                    rec = np.ascontiguousarray(rec, dtype=np.float32)
                    rec = tomopy.remove_ring(rec, theta_min=Rarc, rwidth=Rmaxwidth, thresh_max=Rtmax, thresh=Rthr, thresh_min=Rtmin,out=rec)
                elif func_name == 'castTo8bit':
                    rec = convert8bit(rec, cast8bit_min, cast8bit_max)
                elif func_name == 'bilateral_filter':
                    rec = pyF3D.run_BilateralFilter(rec, spatialRadius=bilateral_srad, rangeRadius=bilateral_rrad)
                elif func_name == 'write_output':
                    dxchange.write_tiff_stack(rec, fname=filenametowrite, start=y*num_sino_per_chunk + sinoused[0])
                print('(took {:.2f} seconds)'.format(time.time()-curtime))
                dofunc+=1
                if dofunc==len(function_list):
                    break
            if y<niter-1 and keepvalues: # Reset original values for next chunk
                angularrange, numangles, projused, num_proj_per_chunk, numprojchunks, numprojused, numrays, anglelist = keepvalues
                
        curtemp = 1 - curtemp
        curfunc = dofunc
        if curfunc==len(function_list):
            break
        axis = slice_dir[function_list[curfunc]]
    print("cleaning up temp files")
    for tmpfile in tempfilenames:
        try:
            os.remove(tmpfile)
        except OSError:
            pass
    print("End Time: "+time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime()))
    print('It took {:.3f} s to process {}'.format(time.time()-start_time,inputPath+filename))

			logging.info('...with center of rotation shifted %f',k)
			rec = tomopy.recon(tomo, tomopy.angles(tomo.shape[0], 270, 270-angularrange), center=cor_rec+k, algorithm='gridrec', filter_name='butterworth', filter_par=butterworthpars)
			
			rec /= pxsize  # intensity values in cm^-1
			

			
			CORtoWrite =cor_rec+k
			
			if doPolarRing:
				logging.info('Doing ring removal (polar mean filter)')
				rec = tomopy.remove_ring(rec, theta_min=Rarc, rwidth=Rmaxwidth, thresh_max=Rtmax, thresh=Rthr, thresh_min=Rtmin)
			
			if pad_sino:
				logging.info('Unpadding...')
				rec = tomopy.circ_mask(rec[:, npad:-npad, npad:-npad], 0)
				CORtoWrite = CORtoWrite - npad
			else:
				rec = tomopy.circ_mask(rec, 0, ratio=1.0, val=0.0)
		


			logging.info('Writing reconstruction slices to %s', iname[x])
			
			if testCOR_insteps:
				filenametowrite = odirectory+'/rec'+iname[x]+'/'+'cor'+str(CORtoWrite)+'_'+iname[x]
			else:
				filenametowrite = odirectory+'/rec'+iname[x]+'/'+iname[x]

			if castTo8bit:
				rec = convert8bit(rec,data_min,data_max)