Python Dpy.read_tracks方法代码示例

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_tracks [as 别名]
def load_tractogram(filename, lazy_load=False):
    """ Loads tractogram files (*.trk or *.tck or *.dpy)

    Parameters
    ----------
    filename : str
        input trk filename
    lazy_load : {False, True}, optional
        If True, load streamlines in a lazy manner i.e. they will not be kept
        in memory and only be loaded when needed.
        Otherwise, load all streamlines in memory.

    Returns
    -------
    streamlines : list of 2D arrays
        Each 2D array represents a sequence of 3D points (points, 3).
    hdr : dict
        header from a trk file
    """
    if 'dpy' in os.path.splitext(filename)[1].lower():
        dpw = Dpy(filename, 'r')
        streamlines = dpw.read_tracks()
        dpw.close()
        return streamlines, {}

    trk_file = nib.streamlines.load(filename, lazy_load)
    return trk_file.streamlines, trk_file.header

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_tracks [as 别名]
def spherical_rois(fdpy,fsr,sq_radius=4):    
    
    
    R=atlantic_points()    
    dpr=Dpy(fdpy,'r')
    T=dpr.read_tracks()
    dpr.close()
    
    center=R['BCC']
    
    refimg=nib.load(fref)
    aff=refimg.get_affine()
    
    SR={}
    
    for key in R:
        
        center=R[key]
        #back to world space
        centerw=np.dot(aff,np.array(center+(1,)))[:3]        
        centerw.shape=(1,)+centerw.shape   
        centerw=centerw.astype(np.float32)
    
        res= [track_roi_intersection_check(t,centerw,sq_radius) for t in T]
        res= np.array(res,dtype=np.int)
        ind=np.where(res>0)[0]
        
        SR[key]={}
        SR[key]['center']=center
        SR[key]['centerw']=tuple(np.squeeze(centerw))
        SR[key]['radiusw']=np.sqrt(sq_radius)
        SR[key]['indices']=ind
        
    
    save_pickle(fsr,SR)

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_tracks [as 别名]
    def loading_full_tractograpy(self, tracpath=None):
        """
        Loading full tractography and creates StreamlineLabeler to
        show it all.
        """
        # load the tracks registered in MNI space
        self.tracpath=tracpath
        basename = os.path.basename(self.tracpath)
        tracks_basename, tracks_format = os.path.splitext(basename)
        
        if tracks_format == '.dpy': 
            
            dpr = Dpy(self.tracpath, 'r')
            print "Loading", self.tracpath
            self.T = dpr.read_tracks()
            dpr.close()
            self.T = np.array(self.T, dtype=np.object)

            
        elif tracks_format == '.trk':
            print "Loading", self.tracpath
            # Old nibabel API:
            # streams, self.hdr = nib.trackvis.read(self.tracpath, points_space='voxel')
            # self.T = np.array([s[0] for s in streams], dtype=np.object)
            # New nibabel API
            tmp = nib.streamlines.load(self.tracpath)
            streams = tmp.tractogram.apply_affine(np.linalg.inv(tmp.affine)).streamlines
            self.header = tmp.header
            self.T = np.array(streams, dtype=np.object)

        # The following code has been commented out to avoid
        # misalignment between original streamlines IDs and final IDs.
        # print "Removing short streamlines"
        # self.T = np.array([t for t in self.T if length(t)>= 15],  dtype=np.object)
        
        tracks_directoryname = os.path.dirname(self.tracpath) + '/.temp/'
        general_info_filename = tracks_directoryname + tracks_basename + '.spa'
        
        
        # Check if there is the .spa file that contains all the
        # computed information from the tractography anyway and try to
        # load it
        try:
            print "Looking for general information file"
            self.load_info(general_info_filename)
                    
        except (IOError, KeyError):
            print "General information not found, recomputing buffers"
            self.update_info(general_info_filename)
                    
        # create the interaction system for tracks, 
        self.streamlab  = StreamlineLabeler('Bundle Picker',
                                            self.buffers, self.clusters,
                                            vol_shape=self.dims, 
                                            affine=np.copy(self.affine),
                                            clustering_parameter=len(self.clusters),
                                            clustering_parameter_max=len(self.clusters),
                                            full_dissimilarity_matrix=self.full_dissimilarity_matrix)
                
        self.scene.add_actor(self.streamlab)

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_tracks [as 别名]
def skeletonize(fdpy,flsc,points=3):

    dpr=Dpy(fdpy,'r')    
    T=dpr.read_tracks()
    dpr.close()    
    print len(T)
    Td=[downsample(t,points) for t in T]
    C=local_skeleton_clustering(Td,d_thr=10.,points=points)    
    #Tobject=np.array(T,dtype=np.object)
    

    #'''
    #r=fvtk.ren()    
    skeleton=[]    
    for c in C:
        #color=np.random.rand(3)
        if C[c]['N']>0:
            Ttmp=[]
            for i in C[c]['indices']:
                Ttmp.append(T[i])
            si,s=most_similar_track_mam(Ttmp,'avg')
            print si,C[c]['N']    
            C[c]['most']=Ttmp[si]            
            #fvtk.add(r,fvtk.line(Ttmp[si],color))            
    print len(skeleton)
    #r=fos.ren()
    #fos.add(r,fos.line(skeleton,color))    
    #fos.add(r,fos.line(T,fos.red))    
    #fvtk.show(r)
    #'''
    
    save_pickle(flsc,C)

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_tracks [as 别名]
def test_dpy():

    fd,fname = mkstemp()    
    dpw = Dpy(fname,'w')    
    A=np.ones((5,3))
    B=2*A.copy()
    C=3*A.copy()        
    dpw.write_track(A)
    dpw.write_track(B)
    dpw.write_track(C)            
    dpw.write_tracks([C,B,A])    
    dpw.close()    
    dpr = Dpy(fname,'r')    
    assert_equal(dpr.version()=='0.0.1',True)    
    T=dpr.read_tracksi([0,1,2,0,0,2])    
    print(T)    
    T2=dpr.read_tracks()    
    assert_equal(len(T2),6)    
    dpr.close()
    assert_array_equal(A,T[0])
    assert_array_equal(C,T[5])
    os.remove(fname)

    
    
    

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_tracks [as 别名]
def roi_track_counts(fdpy,fref,fatlas,roi_no,dist_transf=True,fres=None):
    
    dpr=Dpy(fdpy,'r')
    T=dpr.read_tracks()
    dpr.close()
    
    img=nib.load(fref)
    affine=img.get_affine()
    zooms = img.get_header().get_zooms()
    iaffine=np.linalg.inv(affine)
    T2=[]
    #go back to volume space
    for t in T:
        T2.append(np.dot(t,iaffine[:3,:3].T)+iaffine[:3,3])
    del T
    
    tcs,tes=track_counts(T2,img.get_shape(),zooms,True)
    
    atlas_img=nib.load(fatlas)
    atlas=atlas_img.get_data()
    roi=atlas.copy()
    roi[atlas!=roi_no]=0
    
    if dist_transf:
        roi2=distance_transform_cdt(roi)
        roi[roi2!=roi2.max()]=0
        I=np.array(np.where(roi==roi_no)).T
    else:
        I=np.array(np.where(roi==roi_no)).T    
    
    """    
    if erosion_level>0:
        roi2=binary_erosion(roi,cross,erosion_level)                 
        I=np.array(np.where(roi2==True)).T
    else:        
        roi2=distance_transform_cdt(roi)        
        I=np.array(np.where(roi==roi_no)).T        
    """
    
    #print I.shape    
    #nib.save(nib.Nifti1Image(roi2,affine),'/tmp/test.nii.gz')
    
    Ttes=[]
    for iroi in I:
        try:
            Ttes.append(tes[tuple(iroi)])
        except KeyError:
            pass
    
    Ttes=list(set(list(chain.from_iterable(Ttes))))    
    T2n=np.array(T2,dtype=np.object)
    res=list(T2n[Ttes])
    
    #back to world space
    res2=[]
    for t in res:
        res2.append(np.dot(t,affine[:3,:3].T)+affine[:3,3])
    np.save(fres,np.array(res2,dtype=np.object))

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_tracks [as 别名]
def load_cst(tracks_filename, cst_index_file, ext):
    from dipy.io.dpy import Dpy
    from dipy.io.pickles import load_pickle
    dpr_tracks = Dpy(tracks_filename, 'r')
    all_tracks=dpr_tracks.read_tracks()
    dpr_tracks.close()
    tracks_id = load_pickle(cst_index_file)
    	
    cst = [all_tracks[i] for i  in tracks_id]    
    
    cst_ext = [all_tracks[i] for i  in tracks_id]
    medoid_cst = []
    #len_dis = 250
    if ext:
        k = np.round(len(cst)*1.2)
        not_cst_fil = []
        min_len = min(len(i) for i in cst)
        #print 'min_len of cst', min_len
        min_len = min_len*2.2/3#2./3.2# - 20
        for i in np.arange(len(all_tracks)):
            if (i not in tracks_id) and (length(all_tracks[i]) > min_len):
                not_cst_fil.append(all_tracks[i])
        
        #for st in all_tracks:
        #    if (length(st)>=min_len) and (st not in cst):
        #        not_cst_fil.append(st)
                
        from dipy.segment.quickbundles import QuickBundles
        
        qb = QuickBundles(cst,200,18)
        
        medoid_cst = qb.centroids[0]
        
        med_notcst_dm = bundles_distances_mam([medoid_cst], not_cst_fil)
        med_cst_dm = bundles_distances_mam([medoid_cst], cst)
        
        cst_rad = med_cst_dm[0][np.argmax(med_cst_dm[0])]
        len_dis = cst_rad * 2.8/2.
        #print med_cst_dm
        #print cst_rad
        #print len_dis
        #k_indices which close to the medoid
        sort = np.argsort(med_notcst_dm,axis = 1)[0]
        #print sort[:k+1]
        while (k>0 and med_notcst_dm[0][sort[k]]>=len_dis):
            k = k - 1
            
        #print med_notcst_dm[0][sort[0:k]]    
        #print k
        #close_indices = np.argsort(cst_dm,axis = 1)[:,0:k][0]
        close_indices = sort[0:k]
        
        for idx in close_indices:
            cst_ext.append(not_cst_fil[idx])            
        
        return cst, cst_ext, medoid_cst

    return cst

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_tracks [as 别名]
def get_luigi_SLFI():
    fseg1 = "/home/eg309/Devel/segmented_bundles/luigi_s1_for_eleftherios/S1_SLFI"
    # fseg2='/home/eg309/Devel/segmented_bundles/luigi_s1_for_eleftherios/S1_SLFII'
    subject = "01"
    fdpyw = "data/subj_" + subject + "/101_32/DTI/tracks_gqi_3M_linear.dpy"
    dpr = Dpy(fdpyw, "r")
    T = dpr.read_tracks()
    dpr.close()
    seg_inds = load_pickle(fseg1)
    T1 = [T[i] for i in seg_inds]
    # seg_inds=load_pickle(fseg2)
    # T2=[T[i] for i in seg_inds]
    return T1

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_tracks [as 别名]
def load_data(figure, data_id):   
    
    if figure=='small_dataset':
        filename = 'ALS_Data/'+ str(data_id) + '/DIFF2DEPI_EKJ_64dirs_14/DTI/tracks_dti_10K.dpy'        
    elif figure=='median_dataset':
        filename = 'ALS_Data/' + str(data_id) + '/DIFF2DEPI_EKJ_64dirs_14/DTI/tracks_dti_1M.dpy'    
    elif figure=='big_dataset':
        filename = 'ALS_Data/' + str(data_id) + '/DIFF2DEPI_EKJ_64dirs_14/DTI/tracks_dti_3M.dpy'
    
    print "Loading tracks."
    dpr = Dpy(filename, 'r')
    tracks = dpr.read_tracks()
    dpr.close()
    tracks = np.array(tracks, dtype=np.object)
    return tracks

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_tracks [as 别名]
def load_data(figure, data_id):

    if figure == "small_dataset":
        filename = "ALS_Data/" + str(data_id) + "/DIFF2DEPI_EKJ_64dirs_14/DTI/tracks_dti_10K.dpy"
    elif figure == "median_dataset":
        filename = "ALS_Data/" + str(data_id) + "/DIFF2DEPI_EKJ_64dirs_14/DTI/tracks_dti_1M.dpy"
    elif figure == "big_dataset":
        filename = "ALS_Data/" + str(data_id) + "/DIFF2DEPI_EKJ_64dirs_14/DTI/tracks_dti_3M.dpy"

    print "Loading tracks."
    dpr = Dpy(filename, "r")
    tracks = dpr.read_tracks()
    dpr.close()
    tracks = np.array(tracks, dtype=np.object)
    return tracks

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_tracks [as 别名]
def load_whole_tract(tracks_filename):
    
    from dipy.io.pickles import load_pickle
    if (tracks_filename[-3:]=='dpy'):
        from dipy.io.dpy import Dpy
        dpr_tracks = Dpy(tracks_filename, 'r')
        all_tracks=dpr_tracks.read_tracks()
        dpr_tracks.close()
    else:
        import nibabel as nib
        streams,hdr=nib.trackvis.read(tracks_filename,points_space='voxel')
        all_tracks = np.array([s[0] for s in streams], dtype=np.object)    

    all_tracks = np.array(all_tracks,dtype=np.object)
    return all_tracks

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_tracks [as 别名]
def create_dataset_from_tractography(size1, size2, same=True):
    if same: assert(size2 >= size1)
    filename = 'data/tracks_dti_10K_linear.dpy'
    print "Loading", filename
    dpr = Dpy(filename, 'r')
    tractography = dpr.read_tracks()
    dpr.close()
    print len(tractography), "streamlines"
    print "Removing streamlines that are too short"
    tractography = filter(lambda x: len(x) > 20, tractography) # remove too short streamlines
    print len(tractography), "streamlines"    
    tractography = np.array(tractography, dtype=np.object)

    print "Creating two simulated tractographies of sizes", size1, "and", size2
    if same:
        ids = fft(tractography, k=max([size1, size2]), distance=bundles_distances_mam)
        tractography1 = tractography[ids[:size1]]
    else:
        # ids1 = np.random.permutation(len(tractography))[:size1]
        # ids1 = sff(tractography, k=size1, distance=bundles_distances_mam)
        ids1 = fft(tractography, k=size1, distance=bundles_distances_mam)
        tractography1 = tractography[ids1[:size1]]

    if same:
        tractography2 = tractography[ids[:size2]]
    else:
        # ids2 = np.random.permutation(len(tractography))[:size2]
        # ids2 = sff(tractography, k=size2, distance=bundles_distances_mam)
        ids2 = fft(tractography, k=size2, distance=bundles_distances_mam)
        tractography2 = tractography[ids2]
        
    print "Done."

    print "Computing the distance matrices for each tractography."
    dm1 = bundles_distances_mam(tractography1, tractography1)
    dm2 = bundles_distances_mam(tractography2, tractography2)

    print("Computing similarity matrices.")
    sigma2 = np.mean([np.median(dm1), np.median(dm2)]) ** 2.0
    print("sigma2 = %f" % sigma2)
    A = np.exp(-dm1 * dm1 / sigma2)
    B = np.exp(-dm2 * dm2 / sigma2)

    # Note: the optimization works even using distance instead of similarity:
    # A = dm1
    # B = dm2

    return A, B

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_tracks [as 别名]
def load_tracks(method="pmt"):
    from nibabel import trackvis as tv

    dname = "/home/eg309/Data/orbital_phantoms/dwi_dir/subject1/"

    if method == "pmt":
        fname = "/home/eg309/Data/orbital_phantoms/dwi_dir/workflow/tractography/_subject_id_subject1/cam2trk_pico_twoten/data_fit_pdfs_tracked.trk"
        streams, hdr = tv.read(fname, points_space="voxel")
        tracks = [s[0] for s in streams]
    if method == "dti":
        fname = dname + "dti_tracks.dpy"
    if method == "dsi":
        fname = dname + "dsi_tracks.dpy"
    if method == "gqs":
        fname = dname + "gqi_tracks.dpy"
    if method == "eit":
        fname = dname + "eit_tracks.dpy"
    if method in ["dti", "dsi", "gqs", "eit"]:
        dpr_linear = Dpy(fname, "r")
        tracks = dpr_linear.read_tracks()
        dpr_linear.close()

    if method != "pmt":
        tracks = [t - np.array([96 / 2.0, 96 / 2.0, 55 / 2.0]) for t in tracks if track_range(t, 100 / 2.5, 150 / 2.5)]
    tracks = [t for t in tracks if track_range(t, 100 / 2.5, 150 / 2.5)]

    print "final no of tracks ", len(tracks)
    qb = QuickBundles(tracks, 25.0 / 2.5, 18)
    # from dipy.viz import fvtk
    # r=fvtk.ren()
    # fvtk.add(r,fvtk.line(qb.virtuals(),fvtk.red))
    # fvtk.show(r)
    # show_tracks(tracks)#qb.exemplars()[0])
    # qb.remove_small_clusters(40)
    del tracks
    # load
    tl = TrackLabeler(qb, qb.downsampled_tracks(), vol_shape=None, tracks_line_width=3.0, tracks_alpha=1)

    # return tracks
    w = World()
    w.add(tl)
    # create window
    wi = Window(caption="Fos", bgcolor=(1.0, 1.0, 1.0, 1.0), width=1600, height=900)
    wi.attach(w)
    # create window manager
    wm = WindowManager()
    wm.add(wi)
    wm.run()

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_tracks [as 别名]
def see_tracks(fdpy,N=2000):
    
    
    dpr=Dpy(fdpy,'r')
    #T=dpr.read_tracksi(range(N))
    T=dpr.read_tracks()
    dpr.close()    
    
    T=[downsample(t,5) for t in T]    

    r=fvtk.ren()
    colors=np.ones((len(T),3)).astype('f4')
    for (i,c) in enumerate(T):        
        orient=c[0]-c[-1]
        orient=np.abs(orient/np.linalg.norm(orient))
        colors[i,:3]=orient    
    fvtk.add(r,fvtk.line(T,colors,opacity=0.5))
    #fos.add(r,fos.sphere((0,0,0),10))
    fvtk.show(r)

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_tracks [as 别名]
def see_spherical_intersections(fdpy,fsr):
    
    dpr=Dpy(fdpy,'r')
    T=dpr.read_tracks()
    dpr.close()
    
    SR=load_pickle(fsr)
    
    r=fvtk.ren()
    
    for key in SR:
        ind=SR[key]['indices']
        intersT=[T[i] for i in ind]
        fvtk.add(r,fvtk.line(intersT,np.random.rand(3)))    
        centerw=SR[key]['centerw']
        radius=SR[key]['radiusw']
        fvtk.add(r,fvtk.sphere(position=centerw,radius=radius))
        
    fvtk.show(r)