AlphaSim problem - critical bug?

Dear Chao-Gan,

I use AlphaSim in my VBM workflow to take multiple testing into account when I analyze clusters within ROIs.

Everything seemed to work just fine but for the small structure of the amygdala (bilateral only 1109 voxels, voxel size 1.5 mm) the program stopped with an error message indicating that an index was out of bound:

REST Version: 1.8, Release: 20130615
Citation Information:
Xiao-Wei Song, Zhang-Ye Dong, Xiang-Yu Long, Su-Fang Li, Xi-Nian Zuo, Chao-Zhe Zhu, Yong He, Chao-Gan Yan, Yu-Feng Zang. (2011) REST: A Toolkit for Resting-State Functional Magnetic Resonance Imaging Data Processing. PLoS ONE 6(9): e25031. doi:10.1371/journal.pone.0025031
iter=1  pvoxel=0.000000 zthr=0.014561 mc= mean=0.014561
Attempted to access foneimt(2.12295e+06); index out of bounds
because numel(foneimt)=10000.
...

The error only occured with higher estimated smoothness parameters. I tried to enter different values here and for smaller numbers the program worked.

Most of the time the program worked for sereval iterations before the error occured.

I've looked into the code of the function rest_AlphaSim.m and I think that I might have found the reason: The masking step is done right after generating the random numbers and before applying the gauss filter. If the gauss filter is applied to a masked image, it will blur the bounderies of the mask and values in regions that where initially masked are not zero anymore.

Estimation of the mean and standard deviation head into problems because the numbers of voxels (variable "count") is not equivalent to the numbers of voxels that are non zero. For small structures in combination with high smoothness this can lead to negative expressions in the formula of the sd estimation which comes out as imaginary number and is followed by the break down of the program.

I've moved the masking step in the code and also inserted several smaller edits. Some are just alternatives because we done have neither the "Statistical Toolbox" nor the "Imaging processing Toolbox" (we both were in contact about the replacement of bwlabeln back in 2013/01). Further I was not sure why the estimation of mean and sd takes into account the values of the previous iterations. Since the estimation is now before the masking step and is based on the whole volume, I think it will be stable enough to do it seperately for every iteration and each iteration would then be independent from the others.

I'm not sure about the impact of the relocation of the masking step on other structures but first simulation results indicate higher cluster size thresholds which will help to eliminate false positive results in VBM analyses.

I would be very greatful if you would look at this problem since it might be a critical bug. I'll attach my modification of the matlab code. I've tried to documented all my edits in detail in the code. If you have any questions please ask me.

Best wishes from Germany,
Katharina

PS: I attach the edited version of the m-file to this post.