surfCoupling/R/kth_neighbors_v3.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('finding 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, 'inflatedneighs.asc', sep='.') )
	locfile = file.path(templatedir, 'surf', paste(hemi, 'inflated.asc', sep='.') )
	system(paste('mris_convert -v', file.path(templatedir, 'surf', paste(hemi, 'inflated', sep='.')), neighfile ) )
	system(paste('mris_convert ', file.path(templatedir, 'surf', paste(hemi, 'inflated', sep='.')), locfile ) )
	nt = read.table(neighfile, col.names=c('vertex', 'nneigh', paste('n', 1:6, sep='')), fill=TRUE)
	loc = read.table(locfile, col.names=c('x', 'y', 'z', 'vertex'), skip=2, nrows=nrow(nt) )
	loc$vertex = nt$vertex

	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='')))
	}

	# compute distances for neighbors
	# Y = unlist(nt[1,grep('neighs[0-9]', names(nt))])
	# x = nt[1, 'vertex']
	dist = function(row, loc){
		x = row['vertex']
		Y = unlist(row[grep('neighs[0-9]', names(nt))])
		Y0 = t(loc[ loc$vertex %in% Y, c('x', 'y', 'z')])
		Y0 = sqrt(colSums((Y0 - unlist(loc[ loc$vertex == x, c('x', 'y', 'z') ]) )^2))
		Y0[ match(Y, as.numeric(names(Y0))-1 )]
	}
	dists = apply(nt, 1, dist, loc=loc)
	maxlen = max(unlist(lapply(dists, length)) )
	dists = lapply(dists, function(x) c(x, rep(NA, maxlen-length(x) ) ) )
	dists = do.call(rbind, dists)
	nt[, paste('n', 1:maxlen, 'dist', sep='') ] = dists
	

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