Let T be the number of points in each time series, N the number of voxels in the volume analysed, and M the size of the empirical distribution. The regression procedures used in Phase 1 are linear in T, and insertion of each of the resulting correlations into an ordered binary tree is log(N). Each of these steps is performed once for each voxel analysed. Thus Phase 1 is O(NT + N log(N)).
In Phase 2, the step of randomisation of the time series is again linear in T [Durstenfeld 1964]. Computing correlations for each voxel takes time proportional to N.T, as above. Saving the maximal correlation in an ordered binary tree takes log(M) time. Each of these steps is performed once for each of the M entries in the empirical distribution. Thus Phase 2 is O(MNT + M log(M)).
In Phase 3, extracting the maximal correlation from the binary tree that was constructed in Phase 1 takes log(N) time. Finding the rank of this correlation within the empirical distribution that was constructed in Phase 2 takes log(M) time. Deleting from the empirical distribution the correlations produced by the activated voxel again takes log(M) time. (The tree structure that stores the distribution is keyed both on correlation values and on coordinates, so that both the rank-ordering and coordinate-deletion operations can be performed in logarithmic time.) Each of these steps is performed once for each voxel whose activation is significant. This number of significantly activated voxels generally will be some fraction of N, the total number of voxels analysed. Thus Phase 3 is O(N log(N) + N log(M)).
The entire algorithm is thus O(MNT + N log(N) + M log(M)). For any practical values of these parameters the first term dominates, and so the algorithm behaves linearly in M, N, and T. It is thus desirable to hold M, N, and T to the minimal values necessary to produce adequate statistics. T cannot be usefully reduced, else we would lose a great deal of information in the correlations. M should not be taken much lower than 104; otherwise the empirical distribution, and the resulting adjusted probabilities, would be too grainy. N, though, can be decreased without losing any of the information that we care about, by computing the set of voxels that represent brain tissue and applying the procedure only to those voxels. (We implement this selection of brain voxels using a combination of intensity thresholding and region-growing.) This restriction has the beneficial side effect of making the computed distribution more representative of the tissue under consideration.