BF(?): cart2sphere and sphere2cart are now invertible.

[arokem]

I am not sure about this one. Are these transformations supposed to be invertible? Have I introduced some other problems somewhere else?

I should say, I am having trouble testing this on my machine. I am getting the following error:

ERROR: test (dipy.core.tests.test_geometry.test_sphere_cart)

Traceback (most recent call last):
File "/Library/Frameworks/Python.framework/Versions/7.1/lib/python2.7/site-packages/nose/case.py", line 133, in run
self.runTest(result)
File "/Library/Frameworks/Python.framework/Versions/7.1/lib/python2.7/site-packages/nose/case.py", line 151, in runTest
test(result)
File "/Library/Frameworks/Python.framework/Versions/7.1/lib/python2.7/unittest/case.py", line 391, in call
return self.run(_args, *_kwds)
File "/Users/arokem/source/dipy/dipy/testing/_paramtestpy2.py", line 76, in run
return self.run_parametric(result, testMethod)
File "/Users/arokem/source/dipy/dipy/testing/_paramtestpy2.py", line 56, in run_parametric
result.addError(self, self._exc_info())
AttributeError: 'test_sphere_cart' object has no attribute '_exc_info'

---

Other than that, the tests pass with this commit included.

[Garyfallidis]

Cool lets monitor this for a while :-) Thank you Ariel.

[MrBago]

Sorry I was out of town so I missed this. Can someone clarify this for me? From what I can tell, the change forces the result of the sphere2cart function to be in [0, 2*pi] instead of [-pi, pi]. Why is this desirable? In terms of being invertible,
x, y, z == sphere2cart(cart2sphere(x, y, z)) should always hold true, but r, t, p == cart2sphere(sphere2cart(r, t, p)) in general cannot be true because the resulting phi from cart2sphere is in a fixed range while the acceptable values for the argument phi are in [-inf inf].

Did I miss something?

[Garyfallidis]

No Bago, that's a good point. Ariel?

[arokem]

If I understand the relevant wikipedia page (http://en.wikipedia.org/wiki/Spherical_coordinate_system#Unique_coordinates) correctly, a choice of conventions needs to be made so that every point on a sphere of given radius has a unique representation. According to the docstring of cart2sphere, the chosen convention for these functions is for theta to be limited to [0,pi] and for phi to be [0,2pi](same as the convention wikipedia mentions to be). I don't think it matters whether we choose that or [-pi,pi], as long as we are consistent about it.

For the second point, the question is whether we should enforce the uniqueness of representation of each point or not. What do you think?

[MrBago]

I'd prefer to be consistent with numpy over Wikipedia but I think it's most important that we're consistent with the documentation. Also the following will map phi in the range [a, a+2*pi) and might be faster than the current implementation:

phi = (phi - a) % (2*np.pi) + a

Bago

[arokem]

I am not sure I understand: what's the inconsistency with numpy? Does numpy already have an implementation of these transformations? Or some convention for spherical coordinates? FWIW, a quick googling revealed that scipy.special.sph_harm follows a similar convention (0 => pi) with the names theta and phi reversed.

[arokem]

And yes - it would seem that the transformation you suggested would be faster.

[MrBago]

You're right there isn't a single convention that numpy/scipy adhere to. I actually use scipy.special.sph_harm in the dipy.reconst.shm module. I don't have a strong preference which convention we use for cart2sphere, would it be more helpful for you if we use the [0, 2*pi) instead of [-pi, pi]?

As to enforcing uniqueness, do you mean that sphere2cart would raise an error if the inputs were outside [0, 2_pi]? I don't think this would be a good idea, and I don't see any reason to do it. Functions such as sin and cos accept values outside [0, 2_pi] and even though scipy.special.sph_harm says in the documentation that theta is in[0, 2*pi], it returns the correct result even when theta is not.

On a slight tangent, I've been using a sphere2cart function of my own from dipy.reconst.sph for some time now and I've wanted to get rid of the code redundancy but I've been putting off do it because I will need to deal with the theta-phi switch you mentioned above. But I think it's about time I fixed it, thanks for bringing this up.

[arokem]

I don't have any preference wrt convention ([-pi,pi] vs [0,2pi]). I chose this convention to maintain consistency with the docs that were there.

As for the other point, I don't think that sphere2cart should throw an error with inputs outside this range, just that cart2sphere should never return an answer outside of this range. That is, the operation:

r2,t2,p2= cart2sphere(sphere2cart(r1, t1, p1))

Will yield (r2,t2,p2)==(r1,t1,p2) only if t1 and p1 are within these ranges.

[MrBago]

If you don't have a preference I suggest we use the same convention as arctan2, [-pi, pi], that way we can avoid the extra operation. We should obviously change the docs to reflect this.

I believe this will meet the condition you described above.

[arokem]

What do you think @Garyfallidis ? Was there any reason to prefer this convention over the other? Is this convention used anywhere else in dipy?

[Garyfallidis]

Ariel you are correct. The operation needs to be inverted. According to our docs theta is 0-pi and phi 0-2pi. This is also the most common polar system used. Lets stick to it. So if you can make both sphere2cart and cart2sphere invertible and constrain phi to be always between 0-2pi in both functions that would in general be very useful. Also, check the following case where theta equals 0 and phi becomes negative.

rtp=np.array([[1,0.,2_pi-pi/6.],[1,0,2_pi-pi/3.]])

print cart2sphere(*sphere2cart(rtp[:,0],rtp[:,1],rtp[:,2]))

The argument about numpy and arctan2 is not correct because if numpy would have a cartesian to polar transform they would definitely make it invertible as it should be. The 1-1 mapping needs to be preserved. As it happens with the FFT for and other invertible transforms for example.

Could I also request a tiny addition from you guys? Can you also make these 2 functions work with Nx3 arrays as well? So if the first component is an Nx3 array then ignore the other two components?

-∞+

[Garyfallidis]

By the way I love seeing coders like you guys who go into the details and try to solve every single small bug. Thank you both! :-)

[MrBago]

Let me go through a few things,

First, apparently the @Parametric mechanism is broken, when there is an error in a test it gives the above nonsensical error instead of the proper error. We either need to update dipy.testing with some patch or stop using it.

I have no idea why still @Parametric is used in these tests. I thought we removed all of them. They should all be removed.

Second, because of the above problem we didn't notice a real bug in the current code, which is:

Sorry, yes I knew about this but forgot to mention it.

In [2]: r, theta, phi = cart2sphere(1, 1, 0)
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
/home/bagrata/<ipython-input-2-9e93cbeeaa08> in <module>()
----> 1 r, theta, phi = cart2sphere(1, 1, 0)

/home/bagrata/python/lib/python2.7/site-packages/dipy/core/geometry.pyc in cart2sphere(x, y, z)
    130     # interval [0,2pi]:

    131     idx = np.where(phi < 0)
--> 132     phi[idx] = 2 * np.pi + phi[idx]
    133     idx = np.where(phi > 2 * np.pi)
    134     phi[idx] = phi[idx] - (2 * np.pi)

IndexError: 0-d arrays can only use a single () or a list of newaxes (and a single ...) as an index

Third, I still don't understand the 'invertible' issue. Who cares if the function is invertible between [0, 2pi] or [-pi, pi]. I vote to keep it simple and not add extra code where it isn't needed. Let's just update the docs to reflect the appropriate range.

OK if we change the docs for both functions it makes sense.

Forth, the -0 issue you mention above would require extra code to fix because -0 < 0 is False and -0 % 2*pi == -0. -0s are common in floating point operations and I don't think there is a good reason to address them explicitly here. Besides -0 holds information, in a limit sense, which might be useful to someone so why should we squash it?

Correct, in general -0 issue is unavoidable so let's not check this.

Fifth, I don't like the idea of accepting Nx3 arguments. Either of the following would work if you an Nx3 array:

>>> x, y, z = A.T
>>> R, theta, phi = cart2sphere(x,y,z)

or in one line:
>>> R, theta, phi = cart2sphere(*A.T)

Good point!

So in sort I feel we should change it back to the way it was and fix the docs.

I agree! This way the code will be simpler and less checks are needed.
Guys please change the docs and remove the @Parametric from the tests.
But also add an example and test for inversion.

Bago

GGT :-) Great!

[Garyfallidis]

Bago and Ariel are you going to make the changes or not? Let me know if you need any help! :-)

[MrBago]

Sorry, I should have said something, I made the changes to the code and the docs and took the @Parametric out of test_geometry.py in the following commit:

95c4b78

Bago

[Garyfallidis]

Brilliant!