Hi Ivana,
I’d suggest you use a single (set of) response function(s) for your analysis. Think of it as a common basis that allows you to directly compare the FOD coefficients. If you have reason to believe that the stroke adversely affects single fibre voxels (see output of dwi2response -voxels) in one of the time-points you might want to use the other time-point only, otherwise I’d use all time-points and subjects to calculate a (set of) average response functions (see docs). If there were genuine difference in the response function between the time-points, these would translate to differences in the FOD coefficients if you used one set of response functions. This way you can focus your analysis on the FODs only. If you used two sets of response functions you’d need to related FODs and response functions to each other making the interpretation of the results very hard.
Registration is typically performed and the resulting transformation applied to the FOD or FOD-derived images, not the DWI images as non-rigid reorientation requires a tissue-specific model. So the order is, 1. preprocess your data (as you have done), 2. estimate a group-specific (set of) response function(s), 3. use these to get the ODFs, 4. register the ODFs.
If you want to perform your analysis in the space after registration or in subject-native space and propagate the results to the joint space or a mixture of the two depends on the analysis you want to perform. You’ll need to think about spatial correspondence but also about how (nonlinear) transformations affect the properties of your data. Here is some food for thought.