Hey there Yuxin, 

I agree with Greg's overall approach (reducing 4d to 3d for registration) but I do have a few other thoughts to add for doing this with functional data. 

First I would recommend using an average of a smallish (up to you but i usually use something like 100 frames) number of frames for the 4d-3d average. This might not matter depending on your setup but if there is motion in your imaging setup averaging over many frames can blur the final image. I would also use that same average to do motion correction (gonna make a plug here, but we describe doing this exact procedure with a fly brain back in 2017 in this paper). I recommend using ANTs now for motion correction currently but I've had good results from AFNI's 3dvolreg (it's faster but doesn't do deformations). 

As for the application of the registrations I go by the best practices that human neuroscientists have shared with me. That is that it's best to do as little warping to your functional data as possible. Therefore I always warp the template to the functional data rather than the other way around (This is also what Greg suggested you do as well). Somebody with more technical skills might say it doesn't matter, I honestly don't know. 

Hope that helps. Let me know if you (or Greg) have any other questions. 

-KM