Paul,
It seems that what you want isn't exactly the full residuals, but the data with the covariates projected out. That is, using Gang's notation, after the voxel-wise fit is made (giving the b
.. coefficients)
y
i1 = b
01 + b
11 * cov
i1 + b
21 * cov
i2 + e
i1
for the data in the i
th subject in group 1, then you want the modified data
z
i1 = b
01 + e
i1
to be written out to a dataset for further messing around. Is that right?