Steps 1 and 2 make sense, or you could use FreeSurfer in Step 1 (it's slow but accurate in segmentation).
I'm not sure what you want from Steps 3 and 4. You could use the mask from Step 2 to extract all the time series from CSF and/or WM via 3dmaskdump, and then compute a low dimensional approximation to these using the first few principal components from 1dsvd. Then you could use these as 'baseline' (or 'null hypothesis') regressors in 3dDeconvolve/3dREMLfit, or any other linear regression program of your choice.
