A better way to create a response function for CSD

This implementation includes the recursive calibration of the response function for CSD, Using an FA threshold is not a very robust method.It is dependent on the dataset (non-informed used subjectivity), and still depends on the diffusion tensor (FA and first eigenvector), which has low accuracy at high b-value. This function recursively calibrates the response function, for more information see Tax et al. NeuroImage 2014.

This implementation is slow only on one line (788) which calls real_sph_harm. I would appreciate any help to optimize this function. Thx!
Hey @arokem, thanks a lot for your comments! I'm not sure if this is the case? I want to rotate the main peak direction of every voxel towards the z-axis here. So I actually need to compute the rotation matrix that does this for every voxel independently (this is rotmat). Then I rotate the gradients accordingly, and fit spherical harmonics. Maybe I overlook something, could you perhaps comment on that? Thank you!