# 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='') ) )
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.