Python tomopy.angles函数代码示例

def main(arg):

    parser = argparse.ArgumentParser()
    parser.add_argument("top", help="top directory where the tiff images are located: /data/")
    parser.add_argument("start", nargs='?', const=0, type=int, default=0, help="index of the first image: 10001 (default 0)")

    args = parser.parse_args()

    top = args.top

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

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

    slider(proj)

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

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

    slider(proj)

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

def rec_test(file_name, sino_start, sino_end):

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

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

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

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

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

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

def main(arg):

    parser = argparse.ArgumentParser()
    parser.add_argument("top", help="top directory where the tiff images are located: /data/")
    parser.add_argument("start", nargs='?', const=1, type=int, default=1, help="index of the first image: 1000 (default 1)")

    args = parser.parse_args()

    top = args.top
    index_start = int(args.start)

    template = os.listdir(top)[0]

    nfile = len(fnmatch.filter(os.listdir(top), '*.tif'))
    index_end = index_start + nfile
    ind_tomo = range(index_start, index_end)
    
    fname = top + template

    print (nfile, index_start, index_end, fname)


    # Select the sinogram range to reconstruct.
    start = 0
    end = 512
    sino=(start, end)

    # Read the tiff raw data.
    ndata = dxchange.read_tiff_stack(fname, ind=ind_tomo, slc=(sino, None))
 
    # Normalize to 1 using the air counts
    ndata = tomopy.normalize_bg(ndata, air=5)

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

    ndata = tomopy.minus_log(ndata)

    # Set binning and number of iterations
    binning = 8
    iters = 21

    print("Original", ndata.shape)
    ndata = tomopy.downsample(ndata, level=binning, axis=1)
#    ndata = tomopy.downsample(ndata, level=binning, axis=2)
    print("Processing:", ndata.shape)

    fdir = 'aligned' + '/noblur_iter_' + str(iters) + '_bin_' + str(binning) 

    print(fdir)
    cprj, sx, sy, conv = alignment.align_seq(ndata, theta, fdir=fdir, iters=iters, pad=(10, 10), blur=False, save=True, debug=True)

    np.save(fdir + '/shift_x', sx)
    np.save(fdir + '/shift_y', sy)

    # Write aligned projections as stack of TIFs.
    dxchange.write_tiff_stack(cprj, fname=fdir + '/radios/image')

def main(arg):

    parser = argparse.ArgumentParser()
    parser.add_argument("top", help="top directory where the tiff images are located: /data/")
    parser.add_argument("start", nargs='?', const=1, type=int, default=1, help="index of the first image: 1000 (default 1)")

    args = parser.parse_args()

    top = args.top
    index_start = int(args.start)

    template = os.listdir(top)[0]

    nfile = len(fnmatch.filter(os.listdir(top), '*.tif'))
    index_end = index_start + nfile
    ind_tomo = range(index_start, index_end)
    
    fname = top + template

    print (nfile, index_start, index_end, fname)


    # Select the sinogram range to reconstruct.
    start = 0
    end = 512
    sino=(start, end)

    # Read the tiff raw data.
    ndata = dxchange.read_tiff_stack(fname, ind=ind_tomo, slc=(sino, None))

    print(ndata.shape)
    binning = 8
    ndata = tomopy.downsample(ndata, level=binning, axis=1)
    print(ndata.shape)
    
    # Normalize to 1 using the air counts
    ndata = tomopy.normalize_bg(ndata, air=5)

    ## slider(ndata)

    # Set data collection angles as equally spaced between 0-180 degrees.
    theta = tomopy.angles(ndata.shape[0])
   
    rot_center = 960
    print("Center of rotation: ", rot_center)

    ndata = tomopy.minus_log(ndata)

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

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

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

def generate(phantom="shepp3d", nsize=512, nangles=360):

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

    return [prj, ang, obj]

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

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

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

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

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

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

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

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

        best_center = 1298
        step_01 = time.time()

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

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

    raise DeprecationWarning

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

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

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

def rec_full(file_name, sino_start, sino_end):

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

    chunks = 10  # number of data chunks for the reconstruction

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

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

        if sino_chunk_end > sino_end:
            break

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

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

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

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

        print output_name

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

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

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

    chunks = 10 # number of data chunks for the reconstruction

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

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

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

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

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

def main(arg):

    fname = '/local/dataraid/elettra/Oak_16bit_slice343_all_repack.h5'
    
    # Read the hdf raw data.
    sino, sflat, sdark, th = dxchange.read_aps_32id(fname)

    slider(sino)
    proj = np.swapaxes(sino,0,1)
    flat = np.swapaxes(sflat,0,1)
    dark = np.swapaxes(sdark,0,1)

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

    print(proj.shape, dark.shape, flat.shape, theta.shape)

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

    # Find rotation center.
    rot_center = 962

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

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

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

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