# Faster spherical harmonics implemenation #410

Closed
opened this Issue Aug 21, 2014 · 5 comments

Projects
None yet
6 participants
Contributor

### MrBago commented Aug 21, 2014

 @ChantalTax @stefanv @pv @ken-sakaie I've created an issue to track this. I don't want it to hold up merging #404 but I also don't want to forget about this either because all of dipy will benefit greatly. I see no reason why the spherical harmonic implementation and #404 cannot be handled independently.
Contributor

### stefanv commented Aug 21, 2014

 I agree--this can be handled separately.

Closed

Member

### Garyfallidis commented Sep 3, 2014

 Hi @pv and all, @ChantalTax and me tried what you suggested but it seems the results are not the same in my computer too (scipy 0.13.3 and numpy 1.8.1) . Here is a simple script. ``````import numpy as np from scipy.special import sph_harm, lpmv, gammaln from dipy.reconst.shm import real_sph_harm, sph_harm_ind_list from numpy.testing import assert_almost_equal def real_sph_harm2(m, n, theta, phi): val = sph_harm2(np.abs(m), n, theta, phi) real_sh = np.where(m > 0, val.imag, val.real) real_sh *= np.where(m == 0, 1., np.sqrt(2)) return real_sh def sph_harm2(m, n, theta, phi): x = np.cos(phi) val = lpmv(np.abs(m), n, x).astype(complex) val *= np.sqrt((2*n + 1) / 4.0 / np.pi) val *= np.exp(0.5*(gammaln(n-m+1)-gammaln(n+m+1))) val = val * np.exp(1j * m * theta) return val sh_order = 8 m, n = sph_harm_ind_list(sh_order) theta = np.array([1.61491146, 0.76661665, 0.11976141, 1.20198246, 1.74066314, 1.5925956 , 2.13022055, 0.50332859, 1.19868988, 0.78440679, 0.50686938, 0.51739718, 1.80342999, 0.73778957, 2.28559395, 1.29569064, 1.86877091, 0.39239191, 0.54043037, 1.61263047, 0.72695314, 1.90527318, 1.58186125, 0.23130073, 2.51695237, 0.99835604, 1.2883426 , 0.48114057, 1.50079318, 1.07978624, 1.9798903 , 2.36616966, 2.49233299, 2.13116602, 1.36801518, 1.32932608, 0.95926683, 1.070349 , 0.76355762, 2.07148422, 1.50113501, 1.49823314, 0.89248164, 0.22187079, 1.53805373, 1.9765295 , 1.13361568, 1.04908355, 1.68737368, 1.91732452, 1.01937457, 1.45839 , 0.49641525, 0.29087155, 0.52824641, 1.29875871, 1.81023541, 1.17030475, 2.24953206, 1.20280498, 0.76399964, 2.16109722, 0.79780421, 0.87154509]) phi = np.array([-1.5889514 , -3.11092733, -0.61328674, -2.4485381 , 2.88058822, 2.02165946, -1.99783366, 2.71235211, 1.41577992, -2.29413676, -2.24565773, -1.55548635, 2.59318232, -1.84672472, -2.33710739, 2.12111948, 1.87523722, -1.05206575, -2.85381987, -2.22808984, 2.3202034 , -2.19004474, -1.90358372, 2.14818373, 3.1030696 , -2.86620183, -2.19860123, -0.45468447, -3.0034923 , 1.73345011, -2.51716288, 2.49961525, -2.68782986, 2.69699056, 1.78566133, -1.59119705, -2.53378963, -2.02476738, 1.36924987, 2.17600517, 2.38117241, 2.99021511, -1.4218007 , -2.44016802, -2.52868164, 3.01531658, 2.50093627, -1.70745826, -2.7863931 , -2.97359741, 2.17039906, 2.68424643, 1.77896086, 0.45476215, 0.99734418, -2.73107896, 2.28815009, 2.86276506, 3.09450274, -3.09857384, -1.06955885, -2.83826831, 1.81932195, 2.81296654]) print(sph_harm(m[1], n[1], theta[1], phi[1])) print(sph_harm2(m[1], n[1], theta[1], phi[1])) `````` What do you suggest? Can you share your script?

### pv commented Sep 3, 2014

 You have `val = lpmv(np.abs(m), n, x).astype(complex)` instead of `val = lpmv(m, n, x).astype(complex)` in your script. Here's a fixed script: ``````import numpy as np from scipy.special import lpmv, gammaln def sph_harm(m, n, theta, phi): x = np.cos(phi) val = lpmv(m, n, x).astype(complex) val *= np.sqrt((2*n + 1) / 4.0 / np.pi) val *= np.exp(0.5*(gammaln(n-m+1)-gammaln(n+m+1))) val = val * np.exp(1j * m * theta) return val `````` EDIT: oops, sorry, I had the `abs(m)` in my original copypaste in the last PR too...
Contributor

### ChantalTax commented Sep 4, 2014

 @pv @Garyfallidis, yes now it gives the same results, many thanks!
Member

### arokem commented Dec 4, 2014

 Closed through #413