Python tomopy.find_center函数代码示例

def main():
    normalized = i3.normalize(ct_series, dfs, obs, workdir=os.path.join(workdir, 'normalization'))
    tilt_corrected = i3.correct_tilt(normalized, workdir=os.path.join(workdir, 'tilt-correction'))
    if_corrected = i3.correct_intensity_fluctuation(tilt_corrected, workdir=os.path.join(workdir, 'intensity-fluctuation-correction'))
    angles, sinograms = i3.build_sinograms(if_corrected, workdir=os.path.join(workdir, 'sinogram'))
    # take the middle part to calculate the center of rotation
    sino = [s.data for s in sinograms[900:1100]]
    sino= np.array(sino)
    proj = np.swapaxes(sino, 0, 1)
    rot_center = tomopy.find_center(proj, theta, emission=False, init=1024, tol=0.5)
    rot_center = rot_center[0]
    # reconstruct
    recon = i3.reconstruct(angles, sinograms, workdir=outdir, center=rot_center)
    return

def tomo_reconstruction(sino, omega, algorithm='gridrec',
                        filter_name='shepp', num_iter=1, center=None,
                        refine_center=False, sinogram_order=True):
    '''
    INPUT ->  sino : slice, 2th, x OR 2th, slice, x (with flag sinogram_order=True/False)
    OUTPUT -> tomo : slice, x, y
    '''
    if center is None:
        center = sino.shape[1]/2.
        refine_center = True

    if refine_center:
        center = tomopy.find_center(sino, np.radians(omega), init=center,
                                    ind=0, tol=0.5, sinogram_order=sinogram_order)

    algorithm = algorithm.lower()
    recon_kws = {}
    if algorithm.startswith('gridr'):
        recon_kws['filter_name'] = filter_name
    else:
        recon_kws['num_iter'] = num_iter
    tomo = tomopy.recon(sino, np.radians(omega), algorithm=algorithm,
                        center=center, sinogram_order=sinogram_order, **recon_kws)
    return center, tomo

def recon(io_paras, data_paras, rot_center=None, normalize=True, stripe_removal=10, phase_retrieval=False, 
            opt_center=False, diag_center=False, output="tiff"):
    # 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')

    rot_center_copy = rot_center

    for cycle in xrange(NumCycles):
        # Set start and end of each cycle
        projections_start = 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 cycle #%s (projs %s to %s)" 
            % (cycle, 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:
            print("Normalizing data ...")
            # white = white.mean(axis=0).reshape(-1, *data.shape[1:])
            # dark = dark.mean(axis=0).reshape(-1, *data.shape[1:])
            # data = (data - dark) / (white - dark)
            data = tomopy.normalize(data, white, dark, cutoff=None, ncore=_ncore, nchunk=None)[...]
    
        ## Remove stripes caused by dead pixels in the detector
        if stripe_removal:
            print("Removing stripes ...")
            data = tomopy.remove_stripe_fw(data, level=stripe_removal, wname='db5', sigma=2,
                                    pad=True, ncore=_ncore, 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:
            print("Retrieving phase ...")
            data = tomopy.retrieve_phase(data,
                        pixel_size=1e-4, dist=50, energy=20,
                        alpha=1e-3, pad=True, ncore=_ncore, nchunk=None)
        
        ## Determine and set the center of rotation 
        if opt_center or (rot_center == None):
            ### Using optimization method to automatically find the center
            # d.optimize_center()
            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)
        if diag_center:
            ### Output the reconstruction results using a range of centers,
            ### and then manually find the optimal center.
            # d.diagnose_center()
            if not os.path.exists(diag_cent_dir):
                os.makedirs(diag_cent_dir)
#.........这里部分代码省略.........

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

if __name__ == '__main__':

    # Set path to the micro-CT data to reconstruct.
    fname = 'data_dir/sample_name_prefix'

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

    # Read the APS 1-ID raw data.
    proj, flat, dark = tomopy.read_aps_1id(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)

    # Find rotation center.
    rot_center = tomopy.find_center(proj, theta, emission=False, init=1024, ind=0, tol=0.5)
    print "Center of rotation: ", rot_center

    # Reconstruct object using Gridrec algorithm.
    rec = tomopy.recon(proj, theta, center=rot_center, algorithm='gridrec', emission=False)
        
    # Mask each reconstructed slice with a circle.
    rec = tomopy.circ_mask(rec, axis=0, ratio=0.95)

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

    end = 16

    # APS 26-ID has an x-radia system collecting raw data as xrm.
    proj, flat, metadata = dxchange.read_aps_26id(fname, ind_tomo, ind_flat,
                                                 sino=(start, end))

    # make the darks
    dark = np.zeros((1, proj.shape[1], proj.shape[2]))    

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

    # Find rotation center.
    rot_center = tomopy.find_center(proj, theta, init=1024,
                                    ind=0, tol=0.5)
    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')

#    eng = 31
#    pxl = 0.325e-4
#    rat = 5e-03
#    rat = 1e-03
    #d.phase_retrieval(dist=z, energy=eng, pixel_size=pxl, alpha=rat,padding=True)
    #data = tomopy.retrieve_phase(data, dist=z, energy=eng, pixel_size=pxl, alpha=rat,pad=True)
    
    #if remove_stripe2: d.stripe_removal2()
    if remove_stripe2: data = tomopy.remove_stripe_ti(data)

    #d.downsample2d(level=level) # apply binning on the data
    data = tomopy.downsample(data, level=level) # apply binning on the data
    theta  = tomopy.angles(data.shape[0])
    if 1:
        #if not best_center: d.optimize_center()
        if not best_center: calc_center = tomopy.find_center(data, theta, emission=False, ind=0, tol=0.3)
        else: 
            #d.center=best_center/pow(2,level) # Manage the rotation center
            calc_center = best_center/pow(2,level) # Manage the rotation center
        #d.gridrec(ringWidth=RingW) # Run the reconstruction
        rec = tomopy.recon(data, theta, center=calc_center, algorithm='gridrec', emission=False)
        
        #d.apply_mask(ratio=1)
        rec = tomopy.circ_mask(rec, axis=0)

        # Write data as stack of TIFs.
        #tomopy.xtomo_writer(d.data_recon, output_name, 
        #                    axis=0,
        #                    x_start=slice_first)
        tomopy.io.writer.write_tiff_stack(rec, fname=output_name, axis=0, start=slice_first)

# -*- coding: utf-8 -*-
"""
TomoPy example script to reconstruct the tomography data as
with gridrec.
"""
from __future__ import print_function
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)

#.........这里部分代码省略.........
        sinoused = (0,numslices,1)
    elif sinoused[0]<0:
        sinoused=(int(np.floor(numslices/2.0)-np.ceil(sinoused[1]/2.0)),int(np.floor(numslices/2.0)+np.floor(sinoused[1]/2.0)),1)
    
    num_proj_per_chunk = np.minimum(chunk_proj,projused[1]-projused[0])
    numprojchunks = (projused[1]-projused[0]-1)//num_proj_per_chunk+1
    num_sino_per_chunk = np.minimum(chunk_sino,sinoused[1]-sinoused[0])
    numsinochunks = (sinoused[1]-sinoused[0]-1)//num_sino_per_chunk+1
    numprojused = (projused[1]-projused[0])//projused[2]
    numsinoused = (sinoused[1]-sinoused[0])//sinoused[2]
    
    BeamHardeningCoefficients = (0, 1, 0, 0, 0, .1) if BeamHardeningCoefficients is None else BeamHardeningCoefficients

    if cor is None:
        print("Detecting center of rotation", end="") 
        if angularrange>300:
            lastcor = int(np.floor(numangles/2)-1)
        else:
            lastcor = numangles-1
        #I don't want to see the warnings about the reader using a deprecated variable in dxchange
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            tomo, flat, dark, floc = dxchange.read_als_832h5(inputPath+filename,ind_tomo=(0,lastcor))
        tomo = tomo.astype(np.float32)
        if useNormalize_nf:
            tomopy.normalize_nf(tomo, flat, dark, floc, out=tomo)
        else:
            tomopy.normalize(tomo, flat, dark, out=tomo)

        if corFunction == 'vo':
            # same reason for catching warnings as above
            with warnings.catch_warnings():
                warnings.simplefilter("ignore")
                cor = tomopy.find_center_vo(tomo, ind=voInd, smin=voSMin, smax=voSMax, srad=voSRad, step=voStep,
                                        ratio=voRatio, drop=voDrop)
        elif corFunction == 'nm':
            cor = tomopy.find_center(tomo, tomopy.angles(numangles, angle_offset, angle_offset-angularrange),
                                     ind=nmInd, init=nmInit, tol=nmTol, mask=nmMask, ratio=nmRatio,
                                     sinogram_order=nmSinoOrder)
        elif corFunction == 'pc':
            cor = tomopy.find_center_pc(tomo[0], tomo[1], tol=0.25)
        else:
            raise ValueError("\'corFunction\' must be one of: [ pc, vo, nm ].")
        print(", {}".format(cor))
    else:
        print("using user input center of {}".format(cor))
        
    
    function_list = []

    if doOutliers1D:
        function_list.append('remove_outlier1d')
    if doOutliers2D:
        function_list.append('remove_outlier2d')
    if useNormalize_nf:
        function_list.append('normalize_nf')
    else:
        function_list.append('normalize')
    function_list.append('minus_log')
    if doBeamHardening:
        function_list.append('beam_hardening')
    if doFWringremoval:
        function_list.append('remove_stripe_fw')
    if doTIringremoval:
        function_list.append('remove_stripe_ti')
    if doSFringremoval:

# Read the APS 32-ID or 2-BM raw data
prj, flat, dark = tomopy.io.exchange.read_aps_32id(fname, sino=(start, end))

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

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

# Set the aprox rotation axis location.
# This parameter is the starting angle for auto centering routine
start_center=295 
print "Start Center: ", start_center

# Auto centering
calc_center = tomopy.find_center(prj, theta, emission=False, ind=0, init=start_center, tol=0.3)
print "Calculated Center:", calc_center

# Recon using gridrec
rec = tomopy.recon(prj, theta, center=calc_center, algorithm='gridrec', emission=False)

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

# to save the reconstructed images uncomment and customize the following line:
rec_name = 'rec/tooth'

# Write data as stack of TIFs.
tomopy.io.writer.write_tiff_stack(rec, fname=rec_name)
print "Done!  reconstructions at: ", rec_name