Python Dpy.read_indexed方法代码示例

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_indexed [as 别名]
def skeletonize_both():
    from dipy.viz import fos
    from dipy.core.track_metrics import downsample
    from dipy.core.track_performance import local_skeleton_clustering, most_similar_track_mam
    
    froi='/home/eg309/Data/ICBM_Wmpm/ICBM_WMPM.nii'
    wI=get_roi(froi,9,0) #4 is genu    
    fname='/home/eg309/Data/PROC_MR10032/subj_03/101/1312211075232351192010092217244332311282470ep2dadvdiffDSI10125x25x25STs004a001_QA_warp.dpy'
    #fname='/home/eg309/Data/PROC_MR10032/subj_03/101/1312211075232351192010092217244332311282470ep2dadvdiffDSI10125x25x25STs004a001_QA_native.dpy'
    #fname='/home/eg309/Data/PROC_MR10032/subj_06/101/13122110752323511930000010092916083910900000227ep2dadvdiffDSI10125x25x25STs004a001_QA_native.dpy'
    fname2='/home/eg309/Data/PROC_MR10032/subj_06/101/13122110752323511930000010092916083910900000227ep2dadvdiffDSI10125x25x25STs004a001_QA_warp.dpy'
    r=fos.ren()
    #'''
    dpr=Dpy(fname,'r')    
    T=dpr.read_indexed(range(2*10**4))
    dpr.close()    
    print len(T)    
    Td=[downsample(t,3) for t in T if length(t)>40]
    C=local_skeleton_clustering(Td,d_thr=20.)
    
    for c in C:
        #color=np.random.rand(3)
        color=fos.red
        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']                
            fos.add(r,fos.line(Ttmp[si],color))                

    dpr=Dpy(fname2,'r')    
    T=dpr.read_indexed(range(2*10**4))
    dpr.close()    
    print len(T)    
    Td=[downsample(t,3) for t in T if length(t)>40]
    C=local_skeleton_clustering(Td,d_thr=20.)
    #r=fos.ren()
    for c in C:
        #color=np.random.rand(3)
        color=fos.yellow
        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']                
            fos.add(r,fos.line(Ttmp[si],color))            
    #'''
    fos.add(r,fos.point(wI,fos.green))
    fos.show(r)

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_indexed [as 别名]
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)

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_indexed [as 别名]
def tracks_in_roi():
    
    froi='/home/eg309/Data/ICBM_Wmpm/ICBM_WMPM.nii'
    wI=get_roi(froi,35,1) #4 is genu    
    fname='/home/eg309/Data/PROC_MR10032/subj_03/101/1312211075232351192010092217244332311282470ep2dadvdiffDSI10125x25x25STs004a001_QA_warp.dpy'
    dpr=Dpy(fname,'r')    
    T=dpr.read_indexed(range(2*10**4))
    print len(T)
    Troi=[]
    for t in T:        
        if track_roi_intersection_check(t,wI,.5):
            Troi.append(t)
    print(len(Troi))
    dpr.close()
    
    from dipy.viz import fos
    r=fos.ren()
    fos.add(r,fos.line(Troi,fos.red))
    fos.add(r,fos.point(wI,fos.green))
    fos.show(r)

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_indexed [as 别名]
def skeletonize():
    
    froi='/home/eg309/Data/ICBM_Wmpm/ICBM_WMPM.nii'
    wI=get_roi(froi,35,1) #4 is genu    
    #fname='/home/eg309/Data/PROC_MR10032/subj_03/101/1312211075232351192010092217244332311282470ep2dadvdiffDSI10125x25x25STs004a001_QA_warp.dpy'
    #fname='/home/eg309/Data/PROC_MR10032/subj_03/101/1312211075232351192010092217244332311282470ep2dadvdiffDSI10125x25x25STs004a001_QA_native.dpy'
    #fname='/home/eg309/Data/PROC_MR10032/subj_06/101/13122110752323511930000010092916083910900000227ep2dadvdiffDSI10125x25x25STs004a001_QA_native.dpy'
    fname='/home/eg309/Data/PROC_MR10032/subj_06/101/13122110752323511930000010092916083910900000227ep2dadvdiffDSI10125x25x25STs004a001_QA_warp.dpy'
    
    dpr=Dpy(fname,'r')    
    T=dpr.read_indexed(range(2*10**4))
    dpr.close()
    
    print len(T)    
    from dipy.core.track_metrics import downsample
    from dipy.core.track_performance import local_skeleton_clustering, most_similar_track_mam
    Td=[downsample(t,3) for t in T]
    C=local_skeleton_clustering(Td,d_thr=20.)
    
    #Tobject=np.array(T,dtype=np.object)
    
    from dipy.viz import fos
    r=fos.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']                
            fos.add(r,fos.line(Ttmp[si],color))
            
    #print len(skeleton)
    #fos.add(r,fos.line(skeleton,color))    
    #fos.add(r,fos.line(T,fos.red))    
    fos.show(r)

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_indexed [as 别名]
def warp_tracks():
    dn='/home/eg309/Data/TEST_MR10032/subj_03/101/'
    ffa=dn+'1312211075232351192010092217244332311282470ep2dadvdiffDSI10125x25x25STs004a001_bet_FA.nii.gz'    
    finvw=dn+'1312211075232351192010092217244332311282470ep2dadvdiffDSI10125x25x25STs004a001_warps_in_bet_FA.nii.gz'    
    fqadpy=dn+'1312211075232351192010092217244332311282470ep2dadvdiffDSI10125x25x25STs004a001_QA_native.dpy'
    flaff=dn+'1312211075232351192010092217244332311282470ep2dadvdiffDSI10125x25x25STs004a001_affine_transf.mat'
    fref ='/usr/share/fsl/data/standard/FMRIB58_FA_1mm.nii.gz'    
    fdis =dn+'1312211075232351192010092217244332311282470ep2dadvdiffDSI10125x25x25STs004a001_nonlin_displacements.nii.gz'
    fdis2 =dn+'1312211075232351192010092217244332311282470ep2dadvdiffDSI10125x25x25STs004a001_nonlin_displacements_withaff.nii.gz'
    #read some tracks
    dpr=Dpy(fqadpy,'r')
    T=dpr.read_indexed(range(150))
    dpr.close()
    
    #from fa index to ref index
    res=flirt2aff_files(flaff,ffa,fref)
    
    #load the reference img    
    imgref=ni.load(fref)
    refaff=imgref.get_affine()
    
    #load the invwarp displacements
    imginvw=ni.load(finvw)
    invwdata=imginvw.get_data()
    invwaff = imginvw.get_affine()
    
    #load the forward displacements
    imgdis=ni.load(fdis)
    disdata=imgdis.get_data()
    #load the forward displacements + affine
    imgdis2=ni.load(fdis2)
    disdata2=imgdis2.get_data()
    #from their difference create the affine
    disaff=imgdis2.get_data()-disdata  
    
    shift=np.array([disaff[...,0].mean(),disaff[...,1].mean(),disaff[...,2].mean()])
    
    shape=ni.load(ffa).get_data().shape
    
    disaff0=affine_transform(disaff[...,0],res[:3,:3],res[:3,3],shape,order=1)
    disaff1=affine_transform(disaff[...,1],res[:3,:3],res[:3,3],shape,order=1)
    disaff2=affine_transform(disaff[...,2],res[:3,:3],res[:3,3],shape,order=1)
    
    disdata0=affine_transform(disdata[...,0],res[:3,:3],res[:3,3],shape,order=1)
    disdata1=affine_transform(disdata[...,1],res[:3,:3],res[:3,3],shape,order=1)
    disdata2=affine_transform(disdata[...,2],res[:3,:3],res[:3,3],shape,order=1)
    
    #print disgrad0.shape,disgrad1.shape,disgrad2.shape
    #disdiff=np.empty(invwdata.shape)
    #disdiff[...,0]=disgrad0
    #disdiff[...,1]=disgrad1
    #disdiff[...,2]=disgrad2
    #ni.save(ni.Nifti1Image(disdiff,invwaff),'/tmp/disdiff.nii.gz')
    
    di=disdata0
    dj=disdata1
    dk=disdata2
    
    d2i=invwdata[:,:,:,0] + disaff0
    d2j=invwdata[:,:,:,1] + disaff1
    d2k=invwdata[:,:,:,2] + disaff2
    
    #di=disgrad0
    #dj=disgrad1
    #dk=disgrad2
    
    imgfa=ni.load(ffa)
    fadata=imgfa.get_data()
    faaff =imgfa.get_affine()
    
    Tw=[]
    Tw2=[]
    Tw3=[]
    
    froi='/home/eg309/Data/ICBM_Wmpm/ICBM_WMPM.nii'    
    
    roiI=get_roi(froi,3,1) #3 is GCC     
    roiI2=get_roi(froi,4,1) #4 is BCC
    roiI3=get_roi(froi,5,1) #4 is SCC
    roiI=np.vstack((roiI,roiI2,roiI3))  
  
    for t in T:
        if np.min(t[:,2])>=0:#to be removed
            mci=mc(di,t.T,order=1) #interpolations for i displacement
            mcj=mc(dj,t.T,order=1) #interpolations for j displacement
            mck=mc(dk,t.T,order=1) #interpolations for k displacement            
            D=np.vstack((mci,mcj,mck)).T                        
            WI=np.dot(t,res[:3,:3].T)+res[:3,3]+D#+ shift
            W=np.dot(WI,refaff[:3,:3].T)+refaff[:3,3]
            
            mc2i=mc(d2i,t.T,order=1) #interpolations for i displacement
            mc2j=mc(d2j,t.T,order=1) #interpolations for j displacement
            mc2k=mc(d2k,t.T,order=1) #interpolations for k displacement            
            D2=np.vstack((mc2i,mc2j,mc2k)).T                        
            WI2=np.dot(t,res[:3,:3].T)+res[:3,3]+D2 #+ shift
            W2=np.dot(WI2,refaff[:3,:3].T)+refaff[:3,3]
                        
            WI3=np.dot(t,res[:3,:3].T)+res[:3,3]
            W3=np.dot(WI3,refaff[:3,:3].T)+refaff[:3,3]
            
#.........这里部分代码省略.........

# 需要导入模块: from dipy.io.dpy import Dpy [as 别名]
# 或者: from dipy.io.dpy.Dpy import read_indexed [as 别名]
fdis=d101+'1312211075232351192010092912092080924175865ep2dadvdiffDSI10125x25x25STs005a001_nonlin_displacements.nii.gz'
ffareg=d101+'1312211075232351192010092912092080924175865ep2dadvdiffDSI10125x25x25STs005a001_bet_FA_reg.nii.gz'
flirtaff=d101+'1312211075232351192010092912092080924175865ep2dadvdiffDSI10125x25x25STs005a001_affine_transf.mat'
ftrack=d101+'1312211075232351192010092912092080924175865ep2dadvdiffDSI10125x25x25STs005a001_QA_native.dpy'
froi='/home/eg309/Data/PROC_MR10032/NIFTI_ROIs/AnatomicalROIs/ROI01_GCC.nii'
froi2='/home/eg309/Data/PROC_MR10032/NIFTI_ROIs/AnatomicalROIs/ROI02_BCC.nii'
#froi3='/home/eg309/Data/PROC_MR10032/NIFTI_ROIs/AnatomicalROIs/ROI03_SCC.nii'
froi3='/home/eg309/Downloads/SCC_analyze.nii'

ref_fname = '/usr/share/fsl/data/standard/FMRIB58_FA_1mm.nii.gz'

dpr=Dpy(ftrack,'r')

print dpr.track_no

T=dpr.read_indexed([0,1,2,3,2000,1000000])

for t in T:
    print t.shape
dpr.close()
track=T[4]
im2im = flirt2aff_files(flirtaff, ffa, ref_fname) #ref_name to be replaced by ffareg
print im2im
from dipy.core.track_metrics import length
print len(track)
print length(track)
#ntrack=np.dot(im2im[:3,:3],track.T)+im2im[:3,[3]]
ntrack=np.dot(track,im2im[:3,:3].T)+im2im[:3,3]
print length(ntrack)
#print length(ntrack.T)
print length(ntrack)/length(track)