handled dimension change from nifti to nibabel packages

group_binfile_parcellation.py

@@ -122,13 +122,17 @@ def group_binfile_parcellate( infiles, outfile, K, n_voxels ):
     print "finished reading in data and calculating connectivity\n"
 
     # perform clustering
+
+
+
     eigenval,eigenvec = ncut(W,K)
     eigenvec_discrete = discretisation(eigenvec)
 
     ## write out the results
     group_img=eigenvec_discrete[:,0]
 
     for i in range(1,K):
+    	if not i%10: print i
         group_img=group_img+(i+1)*eigenvec_discrete[:,i]
 
     save(outfile,group_img.todense())

make_image_from_bin.py

@@ -88,9 +88,11 @@ def make_image_from_bin( image, binfile, mask ):
     imdat[imdat>0]=short(a[0:sum(imdat)].flatten())
 
     # write out the image as nifti
-    thdr=nim.header
+    thdr=nim.get_header()
     thdr['scl_slope']=1
+
+    nim_aff = nim.get_affine()
 
-    nim_out = nb.Nifti1Image(imdat, nim.get_affine())
+    nim_out = nb.Nifti1Image(imdat, nim_aff, thdr)
     nim_out.set_data_dtype('int16')
     nim_out.to_filename(image)

make_local_connectivity_ones.py

@@ -102,7 +102,7 @@ def make_local_connectivity_ones( maskfile, outfile ):
 
     # read in the mask
     msk=nb.load(maskfile)
-    msz=shape(msk.get_data())
+    msz=msk.shape
 
     # convert the 3D mask array into a 1D vector
     mskdat=reshape(msk.get_data(),prod(msz))
@@ -111,7 +111,7 @@ def make_local_connectivity_ones( maskfile, outfile ):
     # elements of the mask
     iv=nonzero(mskdat)[0]
     m=len(iv)
-
+    print m, ' # of non-zero voxels in the mask'
     # construct a sparse matrix from the mask
     msk=csc_matrix((range(1,m+1),(iv,zeros(m))),shape=(prod(msz),1))
 
@@ -121,6 +121,7 @@ def make_local_connectivity_ones( maskfile, outfile ):
 
     # loop over all of the voxels in the mask 	
     for i in range(0,m):
+    	if i % 100 == 0: print 'voxel# ', i
 
         # calculate the voxels that are in the 3D neighborhood
         # of the center voxel
@@ -155,8 +156,8 @@ def make_local_connectivity_ones( maskfile, outfile ):
 
     # concatenate the i, j and w_ij into a single vector		
     outlist=sparse_i
-    outlist=append(outlist,sparse_j)
-    outlist=append(outlist,sparse_w)
+    outlist=append(outlist, sparse_j)
+    outlist=append(outlist, sparse_w)
 
     # save the output file to a .NPY file
     save(outfile,outlist)

make_local_connectivity_scorr.py

@@ -108,54 +108,63 @@ def make_local_connectivity_scorr( infile, maskfile, outfile, thresh ):
 
     # read in the mask
     msk=nb.load(maskfile)
-    msz=shape(msk.get_data())
+    msz=msk.shape
 
     # convert the 3D mask array into a 1D vector
     mskdat=reshape(msk.get_data(),prod(msz))
 
     # determine the 1D coordinates of the non-zero 
     # elements of the mask
     iv=nonzero(mskdat)[0]
-
+    
     # read in the fmri data
+    # NOTE the format of x,y,z axes and time dimension after reading
+    # nb.load('x.nii.gz').shape -> (x,y,z,t)
     nim=nb.load(infile)
-    sz=shape(nim.get_data())
+    sz=nim.shape
+    print sz, ' dimensions of the 4D fMRI data'
+
 
     # reshape fmri data to a num_voxels x num_timepoints array
     imdat=reshape(nim.get_data(),(prod(sz[:3]),sz[3]))
-
+    
     # mask the datset to only then in-mask voxels
     imdat=imdat[iv,:]
-
+    imdat_sz = imdat.shape
+
     #zscore fmri time courses, this makes calculation of the
     # correlation coefficient a simple matrix product
-    imdat_s=tile(std(imdat,1),(148,1)).T
+    imdat_s=tile(std(imdat,1),(imdat_sz[1],1)).T
     # replace 0 with really large number to avoid div by zero
     imdat_s[imdat_s==0]=1000000
-    imdat_m=tile(mean(imdat,1),(148,1)).T
+    imdat_m=tile(mean(imdat,1),(imdat_sz[1],1)).T
     imdat=(imdat-imdat_m)/imdat_s
+
     # set values with no variance to zero
     imdat[imdat_s==0]=0
     imdat[isnan(imdat)]=0
-
+    
     # remove voxels with zero variance, do this here
     # so that the mapping will be consistent across
     # subjects
     vndx=nonzero(var(imdat,1)!=0)[0]
     iv=iv[vndx]
+
+    m = len(iv)
+    print m , ' # of non-zero valued or non-zero variance voxels in the mask'
 
     # construct a sparse matrix from the mask
-    msk=csc_matrix((vndx+1,(iv,zeros(len(iv)))),shape=(prod(msz),1))
+    msk=csc_matrix((vndx+1,(iv,zeros(m))),shape=(prod(msz),1))
 
     sparse_i=[]
     sparse_j=[]
     sparse_w=[[]]
 
-    for i in range(0,len(iv)):
+    for i in range(0,m):
         # convert index into 3D and calculate neighbors
-        ndx3d=indx_1dto3d(iv[i],sz[1:])+neighbors
+        ndx3d=indx_1dto3d(iv[i],sz[:-1])+neighbors
         # convert resulting 3D indices into 1D
-        ndx1d=indx_3dto1d(ndx3d,sz[1:])
+        ndx1d=indx_3dto1d(ndx3d,sz[:-1])
         # convert 1D indices into masked versions
         ondx1d=msk[ndx1d].todense()
         # exclude indices not in the mask

make_local_connectivity_tcorr.py

@@ -119,17 +119,17 @@ def make_local_connectivity_tcorr( infile, maskfile, outfile, thresh ):
     # elements of the mask
     iv=nonzero(mskdat)[0]
     m=len(iv)
-
+    print m, '# of non-zero voxels in the mask'
     # read in the fmri data
+    # NOTE the format of x,y,z axes and time dimension after reading
+    # nb.load('x.nii.gz').shape -> (x,y,z,t)
     nim=nb.load(infile)
-    sz=shape(nim.get_data())
-
+    sz=nim.shape
     # reshape fmri data to a num_voxels x num_timepoints array	
     imdat=reshape(nim.get_data(),(prod(sz[:3]),sz[3]))
 
     # construct a sparse matrix from the mask
-    msk=csc_matrix((range(1,m+1),(iv,zeros(m))),shape=(prod(sz[1:]),1))
-
+    msk=csc_matrix((range(1,m+1),(iv,zeros(m))),shape=(prod(sz[:-1]),1))
     sparse_i=[]
     sparse_j=[]
     sparse_w=[]
@@ -138,12 +138,11 @@ def make_local_connectivity_tcorr( infile, maskfile, outfile, thresh ):
 
     # loop over all of the voxels in the mask 	
     for i in range(0,m):
-
         # calculate the voxels that are in the 3D neighborhood
         # of the center voxel
-        ndx3d=indx_1dto3d(iv[i],sz[1:])+neighbors
-        ndx1d=indx_3dto1d(ndx3d,sz[1:])
-
+        ndx3d=indx_1dto3d(iv[i],sz[:-1])+neighbors
+        ndx1d=indx_3dto1d(ndx3d,sz[:-1])
+        #acd
         # restrict the neigborhood using the mask
         ondx1d=msk[ndx1d].todense()
         ndx1d=ndx1d[nonzero(ondx1d)[0]]
@@ -153,11 +152,14 @@ def make_local_connectivity_tcorr( infile, maskfile, outfile, thresh ):
 
         # determine the index of the seed voxel in the neighborhood
         nndx=nonzero(ndx1d==iv[i])[0]
-
+	#print nndx,
+	#print nndx.shape
+	#print ndx1d.shape
+	#print ndx1d
         # exctract the timecourses for all of the voxels in the 
         # neighborhood
         tc=matrix(imdat[ndx1d,:])
-
+	 
         # make sure that the "seed" has variance, if not just
         # skip it
         if var(tc[nndx,:]) == 0:

parcel_naming.py

@@ -138,7 +138,7 @@ def read_and_conform_atlas(atlas_file,atlas_label_file,\
 def main():
 
     try:
-        fsl_path=os.environ['FSL_DIR']
+        fsl_path=os.environ['FSLDIR']
     except KeyError:
         print "FSL_DIR is not set in the environment, is FSL installed?"
         sys.exit()
@@ -171,15 +171,16 @@ def main():
     parcel_filename=sys.argv[1]
     parcel_outname=sys.argv[2]
     parcel_vals=[int(n) for n in sys.argv[3].split(',')]
-
+    print parcel_vals
     print "%s called with %s, %s, %s"%(sys.argv[0],parcel_filename,\
         parcel_outname,",".join([str(i) for i in parcel_vals]))
 
     print "Read in the parcellation results %s"%(parcel_filename)
     # lets read in the parcellation results that we want to label
     parcels_nii=nb.load(parcel_filename)
     parcels_img=parcels_nii.get_data()
-
+    print np.shape(parcels_img)
+    print len(parcel_vals)
     if len(parcel_vals) != np.shape(parcels_img)[3]:
         print "Length of parcel values (%d) != number of parcel images (%d)"%( \
             len(parcel_vals),np.shape(parcels_img)[3])