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!
@MrBago yes you are totally right about the axial symmetric response function. Up to now I didn't do this because, if I'm not mistaken, the function sh_to_rh does this before actual CSD. But it would be nicer to output an axially symmetric response function, I totally agree.

I'm not entirely sure what you mean by using a hemisphere, and to which lines of codes you refer to? Could you shed some light on that?
Thanks!