search

cryoem / eman2

A scientific image processing software suite with a focus on CryoEM and CryoET
143 stars 50 forks source link

corner cases of Transform class #543

Open wjiang opened 2 years ago

wjiang commented 2 years ago

The following code demonstrates that the Transform class does not properly handle the conversion to spin/quaternion for some Euler angles. The two distinct Euler angles with phi and 360-phi will produce identical spin/quaternion. The scipy.spatial.transform.Rotation class is fine for these angles. 

import EMAN2
t = EMAN2.Transform()
for alt in [179.9, 180]:    # 179 works but 180 fails
    for phi in [37, 360-37]:    # 37 or any angle here will reproduce the issue
        t.set_rotation(dict(type='eman', az=0, alt=alt, phi=phi))
        print(t.get_rotation('eman'))
        #print("\t%s" % (t.get_rotation('matrix'))) # ok
        #print("\t%s" % (t.get_rotation('spin')))   # problematic
        print("\t%s" % (t.get_rotation('quaternion')))  # problematic
ghost commented 2 years ago

Hi,

I worked on these issues quite a bit and wanted to help, but couldn’t figure from your letter what the problems were. What dies it mean “it fails”. Did the program crash without giving the answer? What “issue” are you referring to? What “problematic” is supposed to mean?

I would gladly help but I am not sure what are the issues you encounter.

Regards, Paul

On Feb 17, 2022, at 6:56 PM, Wen Jiang @.***> wrote:

 The following code demonstrates that the Transform class does not properly handle the conversion to spin/quaternion for some Euler angles. The two distinct Euler angles with phi and 360-phi will produce identical spin/quaternion. The scipy.spatial.transform.Rotation class is fine for these angles.

import EMAN2 t = EMAN2.Transform() for alt in [179.9, 180]: # 179 works but 180 fails for phi in [37, 360-37]: # 37 or any angle here will reproduce the issue t.set_rotation(dict(type='eman', az=0, alt=alt, phi=phi)) print(t.get_rotation('eman'))

print("\t%s" % (t.get_rotation('matrix'))) # ok

  #print("\t%s" % (t.get_rotation('spin')))   # problematic
  print("\t%s" % (t.get_rotation('quaternion')))  # problematic

— Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you are subscribed to this thread.

wjiang commented 2 years ago

I should have included the outputs:

{'alt': 179.0000000534512, 'az': 0.0, 'phi': 37.00000140230337, 'type': 'eman'}
    {'e0': 0.00827471061039176, 'e1': 0.9483871831915285, 'e2': -0.31732592276992155, 'e3': 0.0027696419845010937, 'type': 'quaternion'}
{'alt': 179.0000000534512, 'az': 0.0, 'phi': 322.9999985976966, 'type': 'eman'}
    {'e0': 0.00827471061039176, 'e1': 0.9483871831915285, 'e2': 0.31732592276992155, 'e3': -0.0027696419845010937, 'type': 'quaternion'}
{'alt': 180.0, 'az': 0.0, 'phi': 37.00000193417059, 'type': 'eman'}
    {'e0': 0.0, 'e1': 0.9483236480200433, 'e2': 0.3173046778822478, 'e3': 0.0, 'type': 'quaternion'}
{'alt': 180.0, 'az': 0.0, 'phi': 322.9999980658294, 'type': 'eman'}
    {'e0': 0.0, 'e1': 0.9483236480200433, 'e2': 0.3173046778822478, 'e3': 0.0, 'type': 'quaternion'}

The quaternions are different for the two Eulers (alt=179, az=0, phi=37) and (alt=179, az=0, phi=323) as expected, but are identical for Eulers (alt=180, az=0, phi=37) and Eulers (alt=180, az=0, phi=37) that should not.

ghost commented 2 years ago

Hi,

I am not that familiar with EMAN convention and I do not have my laptop handy to check it. However, by setting az=0 you likely create a degenerate case and despite have two different values for alt and phi there is actually only one effective angle. Meaning (Alt, 0, phi) Is equivalent to (0,0, phi + alt) Or something like that.

Based on that, I suspect the second pair of results is actually correct.

To break the degeneracy please try to set az in your tests to a very small number. Say 1.0e-5 and check what results you are getting.

Another test is set transform to (Alt, 0, phi) And then recover the angle. You should get (0, 0, something).

Finally, transform class used to break for degenerate cases, one reason is that internally it uses single precision, which makes recovery of angles unstable for az->0.

Regards, Paul

On Feb 18, 2022, at 4:11 PM, Wen Jiang @.***> wrote:

 I should have included the outputs:

{'alt': 179.0000000534512, 'az': 0.0, 'phi': 37.00000140230337, 'type': 'eman'} {'e0': 0.00827471061039176, 'e1': 0.9483871831915285, 'e2': -0.31732592276992155, 'e3': 0.0027696419845010937, 'type': 'quaternion'} {'alt': 179.0000000534512, 'az': 0.0, 'phi': 322.9999985976966, 'type': 'eman'} {'e0': 0.00827471061039176, 'e1': 0.9483871831915285, 'e2': 0.31732592276992155, 'e3': -0.0027696419845010937, 'type': 'quaternion'} {'alt': 180.0, 'az': 0.0, 'phi': 37.00000193417059, 'type': 'eman'} {'e0': 0.0, 'e1': 0.9483236480200433, 'e2': 0.3173046778822478, 'e3': 0.0, 'type': 'quaternion'} {'alt': 180.0, 'az': 0.0, 'phi': 322.9999980658294, 'type': 'eman'} {'e0': 0.0, 'e1': 0.9483236480200433, 'e2': 0.3173046778822478, 'e3': 0.0, 'type': 'quaternion'} The quaternions are different for the two Eulers (alt=179, az=0, phi=37) and (alt=179, az=0, phi=323) as expected, but are identical for Eulers (alt=180, az=0, phi=37) and Eulers (alt=180, az=0, phi=37) that should not.

— Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you commented.

ghost commented 2 years ago

Hi,

I finally checked your problem using my laptop. My example is attached below.

I added there printout using spider convention angles, which I understand much better.

As you can see, for theta=0 (this is angle with respect to Z-axis) the result is correct, As in this case phi and psi angles are equivalent.

For theta=180 the readout angles are also correct. As you certainly know, for each 3D rotation There are two equivalent sets of Euler angles that results in the same effective rotation (The “path” is different though). So the program recovers angles not as phi, theta, but as theta, psi+180. This is still correct. However, as you pointed out, quaternions for these two sets are the same.
To answer the question why, I would have to study quaternion representation more. Generally, if you think about degenerate rotations, i.e., those for which theta is 0 or 180, representations that involve normal vector will have a problem. In naive terms the reason is that rotation of a vector about its axis results in the same vector. Therefore, it matters how one think about the 3D rotation, i.e., which rotation is applied first.

In short, there is no easy way to fix the problem.

Regards, Pawel.

In [10]: ...: t = EMAN2.Transform() ...: for alt in [0,179.9, 180]:^I# 179 works but 180 fails ...: ^Ifor phi in [37, 360-37]:^I# 37 or any angle here will reproduce the issue ...: ^I^Iprint(" ",alt,phi) ...: ^I^It.set_rotation(dict(type='eman', az=0, alt=alt, phi=phi)) ...: ^I^I#t.set_rotation(dict(type='spider', theta=0, psi=alt, phi=phi)) ...: ^I^Iprint(t.get_rotation('spider')) ...: ^I^I#print("\t%s" % (t.get_rotation('matrix')))^I# ok ...: ^I^I#print("\t%s" % (t.get_rotation('spin')))^I# problematic ...: ^I^Iprint("\t%s" % (t.get_rotation('quaternion')))^I# problematic ...: ...: ...: 0 37 {'phi': 0.0, 'psi': 37.00000193417059, 'theta': 0.0, 'type': 'spider'} {'e0': 0.9483236480200433, 'e1': 0.0, 'e2': 0.0, 'e3': 0.3173046702654797, 'type': 'quaternion'} 0 323 {'phi': 0.0, 'psi': 322.9999980658294, 'theta': 0.0, 'type': 'spider'} {'e0': 0.9483236480200433, 'e1': 0.0, 'e2': 0.0, 'e3': -0.3173046702654797, 'type': 'quaternion'} 179.9 37 {'phi': 270.0, 'psi': 127.00000037025436, 'theta': 179.8999999940324, 'type': 'spider'} {'e0': 0.0008279211405182212, 'e1': 0.9479192263265416, 'e2': -0.3171693404042028, 'e3': 0.00028797257061503347, 'type': 'quaternion'} 179.9 323 {'phi': 270.0, 'psi': 52.999999629745616, 'theta': 179.8999999940324, 'type': 'spider'} {'e0': 0.0008279211405182212, 'e1': 0.9479192263265416, 'e2': 0.3171693404042028, 'e3': -0.00028797257061503347, 'type': 'quaternion'} 180 37 {'phi': 0.0, 'psi': 217.0000019341706, 'theta': 180.0, 'type': 'spider'} {'e0': 0.0, 'e1': 0.9483236480200433, 'e2': 0.3173046778822478, 'e3': 0.0, 'type': 'quaternion'} 180 323 {'phi': 0.0, 'psi': 142.9999980658294, 'theta': 180.0, 'type': 'spider'} {'e0': 0.0, 'e1': 0.9483236480200433, 'e2': 0.3173046778822478, 'e3': 0.0, 'type': 'quaternion'}

On Feb 18, 2022, at 5:11 PM, Wen Jiang @.***> wrote:

I should have included the outputs:

{'alt': 179.0000000534512, 'az': 0.0, 'phi': 37.00000140230337, 'type': 'eman'} {'e0': 0.00827471061039176, 'e1': 0.9483871831915285, 'e2': -0.31732592276992155, 'e3': 0.0027696419845010937, 'type': 'quaternion'} {'alt': 179.0000000534512, 'az': 0.0, 'phi': 322.9999985976966, 'type': 'eman'} {'e0': 0.00827471061039176, 'e1': 0.9483871831915285, 'e2': 0.31732592276992155, 'e3': -0.0027696419845010937, 'type': 'quaternion'} {'alt': 180.0, 'az': 0.0, 'phi': 37.00000193417059, 'type': 'eman'} {'e0': 0.0, 'e1': 0.9483236480200433, 'e2': 0.3173046778822478, 'e3': 0.0, 'type': 'quaternion'} {'alt': 180.0, 'az': 0.0, 'phi': 322.9999980658294, 'type': 'eman'} {'e0': 0.0, 'e1': 0.9483236480200433, 'e2': 0.3173046778822478, 'e3': 0.0, 'type': 'quaternion'} The quaternions are different for the two Eulers (alt=179, az=0, phi=37) and (alt=179, az=0, phi=323) as expected, but are identical for Eulers (alt=180, az=0, phi=37) and Eulers (alt=180, az=0, phi=37) that should not.

— Reply to this email directly, view it on GitHub https://github.com/cryoem/eman2/issues/543#issuecomment-1045259722, or unsubscribe https://github.com/notifications/unsubscribe-auth/AD4OFUUJ5HEWVHCX43C6SXDU327ZJANCNFSM5OWKUARA. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub. You are receiving this because you commented.

wjiang commented 2 years ago

a few more tests with the follow code suggest that the quaternions are not self consistent as setting rotation with the returned quaternions would result in a different Transform as measured by the rotation matrix (v2.91). Interestingly, the same tests are fine with an older version (2.23).

import EMAN2
import numpy as np

phi = np.random.random()*360
theta = np.random.random()*180 
psi = np.random.random()*360 

def get_matrix_numpy(transform):
    m = transform.get_rotation('matrix')        
    m = np.array([
        [ m['m11'], m['m12'], m['m13'] ],
        [ m['m21'], m['m22'], m['m23'] ],
        [ m['m31'], m['m32'], m['m33'] ]
        ])
    return m

t = EMAN2.Transform(dict(type='spider', phi=phi, theta=theta, psi=psi))
print(t.get_rotation('spider'))

m = get_matrix_numpy(t)

t.set_rotation(t.get_rotation('spider'))
m2 = get_matrix_numpy(t)
print("m spider\n", m)
print("m2 spider\n", m2)
assert(np.allclose(m, m2, atol=1e-6))   # always pass

t.set_rotation(t.get_rotation('matrix'))
m2 = get_matrix_numpy(t)
print("m2 matrix\n", m2)
assert(np.allclose(m, m2, atol=1e-6))   # always pass

t.set_rotation(t.get_rotation('quaternion'))
m2 = get_matrix_numpy(t)
print("m2 quaternion\n", m2)
assert(np.allclose(m, m2, atol=1e-6))   # pass (v2.23) but fail (v2.91)

Outputs for v2.23:

{'theta': 20.821667151647016, 'phi': 74.98372483565657, 'psi': 336.0864770394823, 'type': 'spider'}
('m spider\n', array([[ 0.61289924,  0.72025269, -0.32494712],
       [-0.7847755 ,  0.6027984 , -0.14408851],
       [ 0.09209748,  0.34332228,  0.93469131]]))
('m2 spider\n', array([[ 0.61289918,  0.72025269, -0.32494712],
       [-0.78477556,  0.60279834, -0.14408851],
       [ 0.09209746,  0.34332228,  0.93469131]]))
('m2 matrix\n', array([[ 0.61289918,  0.72025269, -0.32494712],
       [-0.78477556,  0.60279834, -0.14408851],
       [ 0.09209746,  0.34332228,  0.93469131]]))
('m2 quaternion\n', array([[ 0.61289912,  0.72025269, -0.32494712],
       [-0.78477556,  0.60279828, -0.14408849],
       [ 0.09209745,  0.34332228,  0.93469131]]))

Outputs for v.2.91:

{'phi': 262.6421402967897, 'psi': 264.7865592745186, 'theta': 92.20078133895466, 'type': 'spider'}
m spider
 [[-0.98810965  0.12407576  0.09079918]
 [-0.08522039  0.04956456 -0.99512857]
 [-0.12797174 -0.99103409 -0.03840144]]
m2 spider
 [[-0.98810965  0.12407575  0.09079918]
 [-0.08522039  0.04956456 -0.99512857]
 [-0.12797174 -0.99103409 -0.03840144]]
m2 matrix
 [[-0.98810965  0.12407575  0.09079918]
 [-0.08522039  0.04956456 -0.99512857]
 [-0.12797174 -0.99103409 -0.03840144]]
m2 quaternion
 [[-0.99863003 -0.04771598  0.02147503]
 [ 0.05219268 -0.87907547  0.47381665]
 [-0.00373045  0.4742884   0.88036157]]
ghost commented 2 years ago

Hi,

Good point. Since the results are correct in the older version the problem clearly is not with the code itself but with something that happened between versions.

If so, somebody from Steve’s group should be able to fix the problem.

Regards, Paul

On Feb 22, 2022, at 11:18 AM, Wen Jiang @.***> wrote:

 a few more tests with the follow code suggest that the quaternions are not self consistent as setting rotation with the returned quaternions would result in a different Transform as measured by the rotation matrix (v2.91). Interestingly, the same tests are fine with an older version (2.23).

import EMAN2 import numpy as np

phi = np.random.random()360 theta = np.random.random()180 psi = np.random.random()*360

def get_matrix_numpy(transform): m = transform.get_rotation('matrix')
m = np.array([ [ m['m11'], m['m12'], m['m13'] ], [ m['m21'], m['m22'], m['m23'] ], [ m['m31'], m['m32'], m['m33'] ] ]) return m

t = EMAN2.Transform(dict(type='spider', phi=phi, theta=theta, psi=psi)) print(t.get_rotation('spider'))

m = get_matrix_numpy(t)

t.set_rotation(t.get_rotation('spider')) m2 = get_matrix_numpy(t) print("m spider\n", m) print("m2 spider\n", m2) assert(np.allclose(m, m2, atol=1e-6)) # always pass

t.set_rotation(t.get_rotation('matrix')) m2 = get_matrix_numpy(t) print("m2 matrix\n", m2) assert(np.allclose(m, m2, atol=1e-6)) # always pass

t.set_rotation(t.get_rotation('quaternion')) m2 = get_matrix_numpy(t) print("m2 quaternion\n", m2) assert(np.allclose(m, m2, atol=1e-6)) # pass (v2.23) but fail (v2.91)

Outputs for v2.23:

{'theta': 20.821667151647016, 'phi': 74.98372483565657, 'psi': 336.0864770394823, 'type': 'spider'} ('m spider\n', array([[ 0.61289924, 0.72025269, -0.32494712], [-0.7847755 , 0.6027984 , -0.14408851], [ 0.09209748, 0.34332228, 0.93469131]])) ('m2 spider\n', array([[ 0.61289918, 0.72025269, -0.32494712], [-0.78477556, 0.60279834, -0.14408851], [ 0.09209746, 0.34332228, 0.93469131]])) ('m2 matrix\n', array([[ 0.61289918, 0.72025269, -0.32494712], [-0.78477556, 0.60279834, -0.14408851], [ 0.09209746, 0.34332228, 0.93469131]])) ('m2 quaternion\n', array([[ 0.61289912, 0.72025269, -0.32494712], [-0.78477556, 0.60279828, -0.14408849], [ 0.09209745, 0.34332228, 0.93469131]]))

Outputs for v.2.91:

{'phi': 262.6421402967897, 'psi': 264.7865592745186, 'theta': 92.20078133895466, 'type': 'spider'} m spider [[-0.98810965 0.12407576 0.09079918] [-0.08522039 0.04956456 -0.99512857] [-0.12797174 -0.99103409 -0.03840144]] m2 spider [[-0.98810965 0.12407575 0.09079918] [-0.08522039 0.04956456 -0.99512857] [-0.12797174 -0.99103409 -0.03840144]] m2 matrix [[-0.98810965 0.12407575 0.09079918] [-0.08522039 0.04956456 -0.99512857] [-0.12797174 -0.99103409 -0.03840144]] m2 quaternion [[-0.99863003 -0.04771598 0.02147503] [ 0.05219268 -0.87907547 0.47381665] [-0.00373045 0.4742884 0.88036157]]

— Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you commented.

prbprb2 commented 2 years ago

Hi Everyone,

I recently changed a couple of lines in the Transform class which made the quaternion stuff better for small twists. I wrote some notes on it. There were some issues for pi rotations, which I fixed now, I believe.

The issue is sort of interesting, so I wrote up some notes.

Best, Phil

On Tue, Feb 22, 2022 at 10:31 AM Pawel @.***> wrote:

Hi,

Good point. Since the results are correct in the older version the problem clearly is not with the code itself but with something that happened between versions.

If so, somebody from Steve’s group should be able to fix the problem.

Regards, Paul

On Feb 22, 2022, at 11:18 AM, Wen Jiang @.***> wrote:

 a few more tests with the follow code suggest that the quaternions are not self consistent as setting rotation with the returned quaternions would result in a different Transform as measured by the rotation matrix (v2.91). Interestingly, the same tests are fine with an older version (2.23).

import EMAN2 import numpy as np

phi = np.random.random()360 theta = np.random.random()180 psi = np.random.random()*360

def get_matrix_numpy(transform): m = transform.get_rotation('matrix') m = np.array([ [ m['m11'], m['m12'], m['m13'] ], [ m['m21'], m['m22'], m['m23'] ], [ m['m31'], m['m32'], m['m33'] ] ]) return m

t = EMAN2.Transform(dict(type='spider', phi=phi, theta=theta, psi=psi)) print(t.get_rotation('spider'))

m = get_matrix_numpy(t)

t.set_rotation(t.get_rotation('spider')) m2 = get_matrix_numpy(t) print("m spider\n", m) print("m2 spider\n", m2) assert(np.allclose(m, m2, atol=1e-6)) # always pass

t.set_rotation(t.get_rotation('matrix')) m2 = get_matrix_numpy(t) print("m2 matrix\n", m2) assert(np.allclose(m, m2, atol=1e-6)) # always pass

t.set_rotation(t.get_rotation('quaternion')) m2 = get_matrix_numpy(t) print("m2 quaternion\n", m2) assert(np.allclose(m, m2, atol=1e-6)) # pass (v2.23) but fail (v2.91)

Outputs for v2.23:

{'theta': 20.821667151647016, 'phi': 74.98372483565657, 'psi': 336.0864770394823, 'type': 'spider'} ('m spider\n', array([[ 0.61289924, 0.72025269, -0.32494712], [-0.7847755 , 0.6027984 , -0.14408851], [ 0.09209748, 0.34332228, 0.93469131]])) ('m2 spider\n', array([[ 0.61289918, 0.72025269, -0.32494712], [-0.78477556, 0.60279834, -0.14408851], [ 0.09209746, 0.34332228, 0.93469131]])) ('m2 matrix\n', array([[ 0.61289918, 0.72025269, -0.32494712], [-0.78477556, 0.60279834, -0.14408851], [ 0.09209746, 0.34332228, 0.93469131]])) ('m2 quaternion\n', array([[ 0.61289912, 0.72025269, -0.32494712], [-0.78477556, 0.60279828, -0.14408849], [ 0.09209745, 0.34332228, 0.93469131]]))

Outputs for v.2.91:

{'phi': 262.6421402967897, 'psi': 264.7865592745186, 'theta': 92.20078133895466, 'type': 'spider'} m spider [[-0.98810965 0.12407576 0.09079918] [-0.08522039 0.04956456 -0.99512857] [-0.12797174 -0.99103409 -0.03840144]] m2 spider [[-0.98810965 0.12407575 0.09079918] [-0.08522039 0.04956456 -0.99512857] [-0.12797174 -0.99103409 -0.03840144]] m2 matrix [[-0.98810965 0.12407575 0.09079918] [-0.08522039 0.04956456 -0.99512857] [-0.12797174 -0.99103409 -0.03840144]] m2 quaternion [[-0.99863003 -0.04771598 0.02147503] [ 0.05219268 -0.87907547 0.47381665] [-0.00373045 0.4742884 0.88036157]]

— Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you commented.

— Reply to this email directly, view it on GitHub https://github.com/cryoem/eman2/issues/543#issuecomment-1047982541, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABWHXQRIGTRPPINOP4LBH7TU4O26JANCNFSM5OWKUARA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you are subscribed to this thread.Message ID: @.***>