Nothing
find.clusters<-function(tscores, thresh, chan_hood, thresh_sign, min_nchans){
### NOTE! the code of find cluster is below.
# I embed before the code for follow_clust, an internal function
follow_clust<-function(n_above,current_voxel_id,current_clust_num,above_thresh_ids,above_thresh_times,above_thresh_chans,chan_hood,n_chan){
# NOTE: 25/07/2018 added n_chan, because returned errors
# a = 5
new_members=matrix(0, nrow=1, ncol=n_chan*3)
new_members_ct=0
#
#assign("clust_ids", clust_ids, envir=globalenv())
for (b in current_clust_num : n_above){
if (!clust_ids[above_thresh_ids[b]]){
temp_dist=abs(above_thresh_times[b]-above_thresh_times[current_voxel_id]);
if (above_thresh_chans[current_voxel_id]==above_thresh_chans[b]){
#voxels are at same channel
chan_dist=0;
} else if (chan_hood[above_thresh_chans[current_voxel_id] , above_thresh_chans[b]]) {
#channels are neighbors
chan_dist=1;
} else {
#voxels aren't spatially compatible
chan_dist=2
}
if ((temp_dist+chan_dist)<=1) {
#if voxels are at same time point and neighboring channels OR
#if voxels are at same channel and neighboring time points,
#merge them into the same cluster
clust_ids[above_thresh_ids[b]]<<-current_clust_num;
#keep track of which other voxels are joined to this
#cluster
new_members_ct=new_members_ct+1;
new_members[new_members_ct]=b;
} # if ((temp_dist)
}# if !clust
} # for (b)
return(new_members)
} # end function follow_clust
# https://github.com/dmgroppe/Mass_Univariate_ERP_Toolbox/blob/master/find_clusters.m
#tscores=t(res$t.mat)
#thresh=5
#thresh_sign=-1 # metto un numero positivo
# !!! Data ara arranged nchan*n_tpt
n_chan=dim(tscores)[1] # get number of channels
n_tpt=dim(tscores)[2] # get number of timepoints
# !! I create a clust_ids variable in the Global environment (to be fixed)
clust_ids<<-matrix(0, nrow=n_chan, ncol=n_tpt)
time_matrix=t( replicate(n_chan, 1:n_tpt) ) # matrix 1:nchan * timepoints (chans in riga)
chan_matrix= replicate(n_tpt, 1:n_chan) # crate matrix nchan * 1:timepoints (timepoints in colonna)
if (thresh_sign > 0){
#looking for positive clusters
above_thresh_ids=which(tscores>=thresh);
} else {
#looking for negative clusters
above_thresh_ids=which(tscores<=thresh);
}
above_thresh_times=time_matrix[above_thresh_ids];
above_thresh_chans=chan_matrix[above_thresh_ids];
# Clear chan & time matrix to save memory (just in case)
rm(chan_matrix, time_matrix)
n_above=length(above_thresh_ids);
n_clust=0
# perform the search only if n_above is > 0
# (only if there are actually values above the thresholds)
if (n_above > 0) {
for (a in 1:n_above){
voxel_id=above_thresh_ids[a];
# if there is an NA, it is considered as an already checked voxel
# GIORGIO: not sure this is correct or there is a bug in the code
if (is.na(clust_ids[voxel_id])){
clust_ids[voxel_id]=1
}
## end part added to avoid crashes
if(!clust_ids[voxel_id]){
# this "if" goes if the "voxel" isn't in a cluster yet
n_clust=n_clust+1;
clust_ids[voxel_id]<<-n_clust;
#go through all the remaining voxels and find all the above
#threshold voxels this voxel is neighbors with and give them this
#cluster # G: initially the not_checked voxels are all those above the thresh
voxels_not_checked = rep(1, length(above_thresh_ids));
check_me = rep(0, length(above_thresh_ids));
check_me[a] = 1;
} # close if !clust_ids
while (sum(check_me>0)) {
# find the indices of the non-zero elements and return only the first
first = which(check_me>0)[1]
new = follow_clust(n_above,first,n_clust,above_thresh_ids,above_thresh_times,above_thresh_chans,chan_hood,n_chan)
check_me[new]=1
voxels_not_checked[first]=0 # G: set to 0 the first voxel, now it is checked.
check_me = check_me&voxels_not_checked # G: the voxel to be checked are those "not checked" & with "check me".
} # close while
} # close for( a in 1:n_above)
# NOTE: this works also in the case no cluster above threshold are found
clust_membership=clust_ids;
#clear global clust_ids
# TO BE DONE
# exclude clusters without a minimun number of neighbours
if(!is.null(min_nchans)){
# get cluster values
clust_values=unique(as.numeric(clust_membership))
# exclude 0 (i.e., no cluster)
clust_values=clust_values[clust_values!=0]
# loop and if therea not a least two channels (rows) with
# at least a significant result, exclude that cluster
for (iClust in 1:length(clust_values)){
# initialize
curr_clust=clust_membership
curr_clust[curr_clust!=clust_values[iClust]]=NA
curr_clust=curr_clust!=0 # transform in logical matrix (so 1 means there is a cluster)
curr_check=colSums(curr_clust, na.rm=TRUE)
# if there are not at least min_nchans sign sensors or if for some reasons all check are NA
if (all(curr_check<min_nchans)|all(is.na(curr_check))){
# set 0 to this cluster
clust_membership[clust_membership==clust_values[iClust]]=0
}
}
# transform to NULL if all clusters are excluded
if(all(as.numeric(clust_membership)==0)){
clust_membership=NULL
}
}
# NOTE GIORGIO: now only one out of two results is returned and NOT n_clust
return(clust_membership)
} # end if n_above > 0
} # end function find_cluster
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.