surfCoupling/R/kth_neighbors_v2.R

In PennBBL/imcoScripts: Functions for Estimating Volumetric Inter-Modal Coupling

# This script will differentiate between neighbors of various distances
# so that weights can be given to neighbors of different degrees
# each set of neighbors is neighbors of only that degree
# hemi = commandArgs(TRUE)[1]
templatedir = commandArgs(TRUE)[1]
# max neighbor number
k = commandArgs(TRUE)[2]
cat('estimating neighborhood structure. k =', k, '\n')
hemis=c('lh', 'rh')
for(hemi in hemis){
	# intermediate file generated by freesurfer with first degree neighbors
	neighfile = file.path(templatedir, 'surf', paste(hemi, 'sphereneighs.asc', sep='.') )
	system(paste('mris_convert -v', file.path(templatedir, 'surf', paste(hemi, 'sphere', sep='.')), neighfile ) )
	nt = read.table(neighfile, col.names=c('vertex', 'nneigh', paste('n', 1:6, sep='')), fill=TRUE)

	nt$neighs1 = apply(nt[,grep('^n[0-9]$', names(nt))], 1, function(x) x[ !is.na(x)] )

	# gets the neighbors of the neighbors
	neighoneigh = function(row, df, neighnum ){
		neighout = unique(c(unlist(df$neighs1[ df$vertex %in% unlist(df[ row, paste('neighs', neighnum-1, sep='') ] ) ]), unlist(df[ row, paste('neighs', 1:(neighnum-1), sep='')] ) ) )
		# remove the neighbors who are in lower degree neighbors
		neighout[ which(! neighout %in% c(df$vertex[row], unlist(df[row, paste('neighs', 1:(neighnum-1), sep='')] ) ) ) ] 
	}

	for( i in 2:k){
		cat(i, '\t')
		nt[ , paste('neighs', i, sep='')] = NA
		temp = lapply(1:nrow(nt), neighoneigh, df=nt, neighnum=i)
		eval(parse(text=paste('nt$neighs', i, ' = temp', sep='')))
	}

	saveRDS(nt, file.path(templatedir, 'surf', paste(hemi, '.', k, 'neighbors_degree.rds', sep='') ) )
}