Python dpy.Dpy类代码示例

def tracking_eudx(dir_src, dir_out, verbose=False):

    # Loading FA and evecs data
    fa_name = 'data_' + par_b_tag + '_' + par_dim_tag + '_FA.nii.gz'
    FA, affine = load_nifti(pjoin(dir_src, fa_name), verbose)

    evecs_name = 'data_' + par_b_tag + '_' + par_dim_tag + '_EV.nii.gz'
    evecs, _ = load_nifti(pjoin(dir_src, evecs_name), verbose)

    # Computation of streamlines
    sphere = get_sphere('symmetric724') 
    peak_indices = quantize_evecs(evecs, sphere.vertices)
    streamlines = EuDX(FA.astype('f8'),
                       ind=peak_indices, 
                       seeds=par_eudx_seeds,
                       odf_vertices= sphere.vertices,
                       a_low=par_eudx_threshold)

    # Saving tractography
    voxel_size =  (par_dim_vox,) * 3
    dims = FA.shape[:3]
    hdr = nib.trackvis.empty_header()
    hdr['voxel_size'] = voxel_size
    hdr['voxel_order'] = 'LAS'
    hdr['dim'] = dims
    hdr['vox_to_ras'] = affine
    strm = ((sl, None, None) for sl in streamlines)
    trk_name = 'tractogram_' + par_b_tag + '_' + par_dim_tag + '_' + par_rec_tag + '_' + par_eudx_tag + '.trk'
    trk_out = os.path.join(dir_out, trk_name)
    nib.trackvis.write(trk_out, strm, hdr, points_space='voxel')    

    dpy_out = trk_out.replace('.trk', '.dpy')
    dpy = Dpy(dpy_out, 'w')
    dpy.write_tracks(streamlines)
    dpy.close()

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)

def tracking_prob(dir_src, dir_out, verbose=False):

    wm_name = 'wm_mask_' + par_b_tag + '_' + par_dim_tag + '.nii.gz'
    wm_mask, affine = load_nifti(pjoin(dir_src, wm_name), verbose)

    sh_name = 'sh_' + par_b_tag + '_' + par_dim_tag + '.nii.gz'
    sh, _ = load_nifti(pjoin(dir_src, sh_name), verbose)

    sphere = get_sphere('symmetric724') 

    classifier = ThresholdTissueClassifier(wm_mask.astype('f8'), .5)
    classifier = BinaryTissueClassifier(wm_mask)
    max_dg = ProbabilisticDirectionGetter.from_shcoeff(sh, max_angle=par_trk_max_angle, sphere=sphere)
    seeds = utils.seeds_from_mask(wm_mask, density=2, affine=affine)
    streamlines = LocalTracking(max_dg, classifier, seeds, affine, step_size=par_trk_step_size)
    streamlines = list(streamlines)

    trk_name = 'tractogram_' + par_b_tag + '_' + par_dim_tag + '_' + par_trk_prob_tag + '.trk'
    trk_out = os.path.join(dir_out, trk_name)
 
    save_trk(trk_out, streamlines, affine, wm_mask.shape)

    dpy_out = trk_out.replace('.trk', '.dpy')
    dpy = Dpy(dpy_out, 'w')
    dpy.write_tracks(streamlines)
    dpy.close()

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

def test_picking_trajectories():

    curves=[100*np.random.rand(10,3),100*np.random.rand(5,3),100*np.random.rand(3,3)]
    curves=[100*np.array([[0,0,0],[1,0,0]]), 100*np.array([[0,1,0],[0,1,3]]),100*np.array([[0,2,0],[0,2,3]])]
    
    '''
    from nibabel import trackvis as tv
    #fname='/home/eg309/Data/PROC_MR10032/subj_01/101/1312211075232351192010091419011391228126452ep2dadvdiffDSI25x25x25b4000s003a001_FA_warp.trk'
    fname='/home/eg309/Data/fibers.trk'
    streams,hdr=tv.read(fname)
    T=[s[0] for s in streams]
    curves=T[:200000]
    '''
    fname='/home/eg309/Data/PROC_MR10032/subj_02/101/1312211075232351192010091708112071055601107ep2dadvdiffDSI10125x25x25STs002a001_QA_native.dpy'
    from dipy.io.dpy import Dpy
    dpr=Dpy(fname,'r')
    T=dpr.read_indexed(range(20000))
    
    from dipy.core.track_metrics import length
    curves=[t for t in T if length(t) > 20]
    
    dpr.close()
    #colors=np.random.rand(len(curves),4).astype('f4')
    colors=0.5*np.ones((len(curves),4)).astype('f4')
    for (i,c) in enumerate(curves):        
        orient=c[0]-c[-1]
        orient=np.abs(orient/np.linalg.norm(orient))
        colors[i,:3]=orient
    
    c=InteractiveCurves(curves,colors=colors)
    w=World()
    w.add(c)
    wi=Window()
    wi.attach(w)

def warp_tracks_linearly(flirt_filename,fa_filename, tracks_filename,linear_filename):
    import nibabel as nib
    from dipy.external.fsl import flirt2aff
    
    fsl_ref = '/usr/share/fsl/data/standard/FMRIB58_FA_1mm.nii.gz'
    
    img_fa = nib.load(fa_filename)            

    flirt_affine= np.loadtxt(flirt_filename)    
        
    img_ref =nib.load(fsl_ref)
    
    #create affine matrix from flirt     
    mat=flirt2aff(flirt_affine,img_fa,img_ref)        

    #read tracks    
    tensor_tracks = load_whole_tract(tracks_filename)    
        
    #linear tranform for tractography
    tracks_warped_linear = transform_tracks(tensor_tracks,mat)        

    #save tracks_warped_linear    
    dpr_linear = Dpy(linear_filename, 'w')
    dpr_linear.write_tracks(tracks_warped_linear)
    dpr_linear.close()

def warp_tracks(input_tracks_filename, input_flirt_fmatrix, input_fa_filename,  output_filename = None, input_ref = '/usr/share/fsl/data/standard/FMRIB58_FA_1mm.nii.gz'):

    print 'Loading fa, flirt matrix ...'    
    img_fa = nib.load(input_fa_filename)            

    flirt_affine= np.loadtxt(input_flirt_fmatrix)    
        
    img_ref =nib.load(input_ref)
    
    #create affine matrix from flirt     
    mat=flirt2aff(flirt_affine,img_fa,img_ref)        

    #read tracks    
    print 'Loading tracks ...'
    tensor_tracks = load_tracks(input_tracks_filename)
        
    #linear tranform for tractography
    tracks_warped_linear = transform_tracks(tensor_tracks,mat)        

    if output_filename == None:
        filename_save = input_tracks_filename.split('.')[0]+'_linear.dpy'
    else:
        filename_save = os.path.abspath(output_filename)
        
    #save tracks_warped_linear    
    print 'Saving warped tracks into :', filename_save
    dpr_linear = Dpy(filename_save, 'w')
    dpr_linear.write_tracks(tracks_warped_linear)
    dpr_linear.close()
    
    return filename_save

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)

    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)

def save_dpy(streamlines, filename):
    ''' Save tractography to a .dpy file'''

    print "Save tracks as .dpy"	
    tracks = [track for track in streamlines]
    dpw = Dpy(filename, 'w')
    dpw.write_tracks(tracks)
    dpw.close()

def write_tracks(fdpy,scalar,indices,seed_no=10**6,a_thr=.2,compression=1):
    eudx=EuDX(scalar,indices,seed_no=seed_no,a_low=a_thr)
    #exi=iter(eudx)       
    dpw=Dpy(fdpy,'w',compression=1)    
    #for (i,track) in enumerate(exi):
    for track in eudx:       
        dpw.write_track(track.astype(np.float32))
    dpw.close()

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

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

def create_save_tracks(anisotropy,indices, seeds, low_thresh,filename):
    
    euler = EuDX(anisotropy, 
                        ind=indices, 
                        #odf_vertices=get_sphere('symmetric362'),
                        seeds=seeds, a_low=low_thresh)
                        #odf_vertices=get_sphere('symmetric362').vertices, 
    tracks = [track for track in euler]     
    dpw = Dpy(filename, 'w')
    dpw.write_tracks(tracks)
    dpw.close()

def create_save_tracks(anisotropy, indices, seeds, low_thresh, filename):

    euler = EuDX(
        anisotropy, ind=indices, odf_vertices=get_sphere("symmetric362").vertices, seeds=seeds, a_low=low_thresh
    )

    tensor_tracks_old = [track for track in euler]
    tracks = [track for track in tensor_tracks_old if track.shape[0] > 1]
    dpw = Dpy(filename, "w")
    dpw.write_tracks(tracks)
    dpw.close()