Hmm. In calc_shift() I found this:
/*----- enter and sort array of differences -----*/
for (i = 0; i < m; i++)
for (j = 0; j < n; j++)
node_addvalue (&head, yarray[j] - xarray
);
which is done once per voxel. So that's certainly going to take some time.
We sometimes use Ns and Ms of over 100, too, and we need to raise MAX_OBSERVATIONS to do that. Maybe the default could be set to 300?