Python tomopy.remove_stripe_fw函数代码示例

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 = 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 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 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 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 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 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 center(io_paras, data_paras, center_start, center_end, center_step, diag_cycle=0, 
            mode='diag', normalize=True, stripe_removal=10, phase_retrieval=False):
    
    # Input and output
    datafile = io_paras.get('datafile')
    path2white = io_paras.get('path2white', datafile)
    path2dark = io_paras.get('path2dark', path2white)
    out_dir = io_paras.get('out_dir')
    diag_cent_dir = io_paras.get('diag_cent_dir', out_dir+"/center_diagnose/")
    recon_dir = io_paras.get('recon_dir', out_dir+"/recon/")
    out_prefix = io_paras.get('out_prefix', "recon_")

    # Parameters of dataset
    NumCycles = data_paras.get('NumCycles', 1) # Number of cycles used for recon
    ProjPerCycle = data_paras.get('ProjPerCycle') # Number of projections per cycle, N_theta
    cycle_offset = data_paras.get('cycle_offset', 0) # Offset in output cycle number
    proj_start = data_paras.get('proj_start', 0) # Starting projection of reconstruction 
    proj_step = data_paras.get('proj_step')
    z_start = data_paras.get('z_start', 0)
    z_end = data_paras.get('z_end', z_start+1)
    z_step = data_paras.get('z_step')
    x_start = data_paras.get('x_start')
    x_end = data_paras.get('x_end', x_start+1)
    x_step = data_paras.get('x_step')
    white_start = data_paras.get('white_start')
    white_end = data_paras.get('white_end')
    dark_start = data_paras.get('dark_start')
    dark_end = data_paras.get('dark_end')

    # Set start and end of each subcycle
    projections_start = diag_cycle * ProjPerCycle + proj_start
    projections_end = projections_start + ProjPerCycle
    slice1 = slice(projections_start, projections_end, proj_step)
    slice2 = slice(z_start, z_end, z_step)
    slice3 = slice(x_start, x_end, x_step)
    slices = (slice1, slice2, slice3)
    white_slices = (slice(white_start, white_end), slice2, slice3)
    dark_slices = (slice(dark_start, dark_end), slice2, slice3)
    print("Running center diagnosis (projs %s to %s)" 
        % (projections_start, projections_end))
    
    # Read HDF5 file.
    print("Reading datafile %s..." % datafile, end="")
    sys.stdout.flush()
    data, white, dark = reader.read_aps_2bm(datafile, slices, white_slices, dark_slices, 
                                    path2white=path2white, path2dark=path2dark)
    theta = gen_theta(data.shape[0])
    print("Done!")
    print("Data shape = %s;\nwhite shape = %s;\ndark shape = %s." 
        % (data.shape, white.shape, dark.shape))
    
    ## Normalize dataset using data_white and data_dark
    if normalize:
        data = tomopy.normalize(data, white, dark, cutoff=None, ncore=_ncore, nchunk=None)

    ## Remove stripes caused by dead pixels in the detector
    if stripe_removal:
        data = tomopy.remove_stripe_fw(data, level=stripe_removal, wname='db5', 
                                        sigma=2, pad=True, ncore=None, nchunk=None)
        # data = tomopy.remove_stripe_ti(data, nblock=0, alpha=1.5, 
        #                                 ncore=None, nchunk=None)
    
#        # Show preprocessed projection
#        plt.figure("%s-prep" % projections_start)
#        plt.imshow(d.data[0,:,:], cmap=cm.Greys_r)
#        plt.savefig(out_dir+"/preprocess/%s-prep.jpg" 
#                    % projections_start)
#        # plt.show()
#        continue

    ## Phase retrieval
    if phase_retrieval:
        data = tomopy.retrieve_phase(data,
                    pixel_size=6.5e-5, dist=33, energy=30,
                    alpha=1e-3, pad=True, ncore=_ncore, nchunk=None)
    
    ## Determine and set the center of rotation
    ### Using optimization method to automatically find the center
    # d.optimize_center()
    if 'opti' in mode:
        print("Optimizing center ...", end="")
        sys.stdout.flush()
        rot_center = tomopy.find_center(data, theta, ind=None, emission=True, init=None,
                                        tol=0.5, mask=True, ratio=1.)
        print("Done!")
        print("center = %s" % rot_center)
    ### Output the reconstruction results using a range of centers,
    ### and then manually find the optimal center.
    if 'diag' in mode:
        if not os.path.exists(diag_cent_dir):
            os.makedirs(diag_cent_dir)
        print("Testing centers ...", end="")
        sys.stdout.flush()
        tomopy.write_center(data, theta, dpath=diag_cent_dir, 
                            cen_range=[center_start, center_end, center_step], 
                            ind=None, emission=False, mask=False, ratio=1.)
        print("Done!")

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

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

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

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

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

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

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

    proj = tomopy.minus_log(proj)

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

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

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

def fast_tomo_recon(argv):
    """
    Reconstruct subset slices (sinograms) equally spaced within tomographic
    dataset
    """

    logger = logging.getLogger('fast_tomopy.fast_tomo_recon')

    # Parse arguments passed to function
    parser = argparse.ArgumentParser()
    parser.add_argument('-i', '--input', type=str, help='path to input raw '
                        'dataset', required=True)
    parser.add_argument('-o', '--output-file', type=str, help='full path to h5 output '
                        'file', default=os.path.join(os.getcwd(), "fast-tomopy.h5"))
    parser.add_argument('-sn', '--sino-num', type=int, help='Number of slices '
                        'to reconstruct', default=5)
    parser.add_argument('-a', '--algorithm', type=str, help='Reconstruction'
                        ' algorithm', default='gridrec',
                        choices=['art', 'bart', 'fbp', 'gridrec', 'mlem',
                                 'ospml_hybrid', 'ospml_quad', 'pml_hybrid',
                                 'pml_quad', 'sirt'])
    parser.add_argument('-c', '--center', type=float, help='Center of rotation',
                        default=None)
    parser.add_argument('-fn', '--filter-name', type=str, help='Name of filter'
                        ' used for reconstruction',
                        choices=['none', 'shepp', 'cosine', 'hann', 'hamming',
                                 'ramlak', 'parzen', 'butterworth'],
                        default='butterworth')
    parser.add_argument('-rr', '--ring-remove', type=str, help='Ring removal '
                        'method', choices=['Octopus', 'Tomopy-FW', 'Tomopy-T'],
                        default='Tomopy-FW')
    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 = tomopy.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_nf(first_last, flats, darks, floc)
        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(tomo.shape[0], 270, 90)

    logger.info('Reconstructing normalized data')
    # Reconstruct sinograms
    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
#.........这里部分代码省略.........

    start = 500
    end = 501
    sino = (start, end)

    # Read APS 32-BM raw data.
    proj, flat, dark, theta = dxchange.read_aps_32id(fname, 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=5,wname='sym16',sigma=1,pad=True)
    #data = tomopy.prep.stripe.remove_stripe_ti(data,alpha=7)
    #data = tomopy.prep.stripe.remove_stripe_sf(data,size=51)

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

    # Find rotation center
    # rot_center = 955
    # rot_center = 953.25
    # rot_center = tomopy.find_center(data, theta, init=rot_center, ind=0, tol=0.5)
    rot_center = tomopy.find_center_vo(data)   
    print(h5name, rot_center)

    data = tomopy.minus_log(data)

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

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

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

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

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

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

    fnproj = fbase + '_' + fnum + fext
    fnflat = fbase + '_' + fnum_flat + fext
    fndark = fbase + '_' + fnum_dark + fext

    fnflat = '/local/data/2018-11/Chawla/1G_A/1G_A_0002.hdf'
    fndark = '/local/data/2018-11/Chawla/1G_A/1G_A_0003.hdf'

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

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

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

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

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

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

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

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

    data = tomopy.minus_log(data)

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

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

    # 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)
    #rec = np.swapaxes(rec,0,2)
    
    return rec