Torsten may know if there is a cleverer way to do this, but one simple option is to generate 3D stacks for each timepoint and then apply the registration to each of them. This is how munger handles 4D (3D + colour) multichannel stacks. If you had provided the following stacks in your images/ folder

stack01.nrrd - time averaged sample stack
stack02.nrrd - timepoint 1
stack03.nrrd - timepoint 2
…

Then munger would handle this automatically if you pass in options

-r "010203"

I'm sure you could also script something this for your own data fairly easily.

As an alternative, if your analysis requires some kind of mask in reference space, have you considered moving that into the sample space? You would do a lot less reformatting.

Best,

Greg.