[vlc-devel] [PATCH 18/18] viewpoint: use quaternion instead of euler angles
Alexandre Janniaux
ajanni at videolabs.io
Thu Apr 1 07:58:56 UTC 2021
Hi,
On Thu, Apr 01, 2021 at 08:07:28AM +0200, Steve Lhomme wrote:
> On 2021-03-31 11:25, Alexandre Janniaux wrote:
> > It enables viewpoint producers to provide an orientation without
> > singularities at north and south poles.
> >
> > Refs #18089, #18760
>
> These 2 issues are about HMD integration and say nothing about a singularity
> pole issue.
HMD headsets are using quaternions, and are sensible to the
gimbal lock issue, hence the referencing of those issues.
This indeed does not solve the issues, hence why I didn´t
write that it fixes them. ;)
> > + // TODO: normalize quaternion
> What happens with non normalized input if we push this patch ?
Sorry, this is a stray comment. Normalization is done at the
conversion sites. You can search for the comment on unit or for
square_norm. I´ll remove the comment locally.
> > ---
> > include/vlc_viewpoint.h | 18 +--
> > src/misc/viewpoint.c | 240 ++++++++++++++++++++++++++++++++--------
> > 2 files changed, 205 insertions(+), 53 deletions(-)
> >
> > diff --git a/include/vlc_viewpoint.h b/include/vlc_viewpoint.h
> > index b3fbe812c7..b3336b5ae0 100644
> > --- a/include/vlc_viewpoint.h
> > +++ b/include/vlc_viewpoint.h
> > @@ -37,24 +37,28 @@
> > /**
> > * Viewpoints
> > */
> > +
> > struct vlc_viewpoint_t {
> > - float yaw; /* yaw in degrees */
> > - float pitch; /* pitch in degrees */
> > - float roll; /* roll in degrees */
> > + /**
> > + * orientation quaternion with the following properties:
> > + * 1/ use ijk = -1 for the operations
> > + * 2/ memory layout is [x, y, z, w] (like GLSL)
> > + * 3/ system is right-handed
> > + */
> > + float quat[4];
> > float fov; /* field of view in degrees */
> > };
> >
> > static inline void vlc_viewpoint_init( vlc_viewpoint_t *p_vp )
> > {
> > - p_vp->yaw = p_vp->pitch = p_vp->roll = 0.0f;
> > + p_vp->quat[3] = 1;
> > + p_vp->quat[0] = p_vp->quat[1] = p_vp->quat[2] = 0;
> > p_vp->fov = FIELD_OF_VIEW_DEGREES_DEFAULT;
> > }
> >
> > static inline void vlc_viewpoint_clip( vlc_viewpoint_t *p_vp )
> > {
> > - p_vp->yaw = fmodf( p_vp->yaw, 360.f );
> > - p_vp->pitch = fmodf( p_vp->pitch, 360.f );
> > - p_vp->roll = fmodf( p_vp->roll, 360.f );
> > + // TODO: normalize quaternion
>
> What happens with non normalized input if we push this patch ?
>
> > p_vp->fov = VLC_CLIP( p_vp->fov, FIELD_OF_VIEW_DEGREES_MIN,
> > FIELD_OF_VIEW_DEGREES_MAX );
> > }
> > diff --git a/src/misc/viewpoint.c b/src/misc/viewpoint.c
> > index 72f68bb40f..ea03c6ced1 100644
> > --- a/src/misc/viewpoint.c
> > +++ b/src/misc/viewpoint.c
> > @@ -25,63 +25,211 @@
> > #endif
> >
> > #include <vlc_viewpoint.h>
> > +#include <assert.h>
> >
> > -void vlc_viewpoint_to_4x4( const vlc_viewpoint_t *vp, float *m )
> > +/* Quaternion to/from Euler conversion. */
> > +
> > +static void QuaternionToEuler(float *yaw, float *pitch, float *roll, const float *q)
> > {
> > - float yaw = vp->yaw * (float)M_PI / 180.f;
> > - float pitch = vp->pitch * (float)M_PI / 180.f;
> > - float roll = vp->roll * (float)M_PI / 180.f;
> > -
> > - float s, c;
> > -
> > - s = sinf(pitch);
> > - c = cosf(pitch);
> > - const float x_rot[4][4] = {
> > - { 1.f, 0.f, 0.f, 0.f },
> > - { 0.f, c, -s, 0.f },
> > - { 0.f, s, c, 0.f },
> > - { 0.f, 0.f, 0.f, 1.f } };
> > -
> > - s = sinf(yaw);
> > - c = cosf(yaw);
> > - const float y_rot[4][4] = {
> > - { c, 0.f, s, 0.f },
> > - { 0.f, 1.f, 0.f, 0.f },
> > - { -s, 0.f, c, 0.f },
> > - { 0.f, 0.f, 0.f, 1.f } };
> > -
> > - s = sinf(roll);
> > - c = cosf(roll);
> > - const float z_rot[4][4] = {
> > - { c, s, 0.f, 0.f },
> > - { -s, c, 0.f, 0.f },
> > - { 0.f, 0.f, 1.f, 0.f },
> > - { 0.f, 0.f, 0.f, 1.f } };
> > -
> > - /**
> > - * Column-major matrix multiplication mathematically equal to
> > - * z_rot * x_rot * y_rot
> > + /* The matrix built from the angles is made from the multiplication of the
> > + * following matrices:
> > + * ⎡ cos(yaw) 0 -sin(yaw) ⎤
> > + * m_yaw (y_rot) = ⎢ 0 1 0 ⎥
> > + * ⎣ sin(yaw) 0 cos(yaw) ⎦
> > + *
> > + * ⎡ 1 0 0 ⎤
> > + * m_pitch (x_rot) = ⎢ 0 cos(pitch) sin(pitch) ⎥
> > + * ⎣ 0 -sin(pitch) cos(pitch) ⎦
> > + *
> > + * ⎡ cos(roll) sin(roll) 0 ⎤
> > + * m_roll (z_rot) = ⎢ -sin(roll) cos(roll) 0 ⎥
> > + * ⎣ 0 0 1 ⎦
> > + *
> > + * Which, multiplied in the correct order will bring, with the symbols
> > + * rewritten: sin = s , cos = c, yaw = y, pitch = p, roll = r
> > + *
> > + * ⎡s(p)⋅s(r)⋅s(y) + c(r)⋅c(y) s(r)⋅c(p) s(p)⋅s(r)⋅c(y) - s(y)⋅c(r)⎤
> > + * V = ⎢s(p)⋅s(y)⋅c(r) - s(r)⋅c(y) c(p)⋅c(r) s(p)⋅c(r)⋅c(y) + s(r)⋅s(y)⎥
> > + * ⎣ s(y)⋅c(p) -s(p) c(p)⋅c(y) ⎦
> > + *
> > + * We can first extract pitch = atan2( -V_32, sqrt(V_31^2 + V_33^2) )
> > + *
> > + * By taking the case |pitch| = 90 degree, it simplify c(y) and s(y) and:
> > + * roll = atan2( V_11, -V_21 )
> > + * yaw = atan2( V_11, -V_13 )
> > + *
> > + * Otherwise, |pitch| != 90 degree and we can get:
> > + * roll = atan2( V_12, V_22 )
> > + * yaw = atan2( V_31, V_33 )
> > + *
> > + * The equivalent matrix obtained by converting the equivalent quaternion
> > + * Q = ((x, y, z), w) into a tranform matrix (\ref vlc_viewpoint_to_4x4)
> > + * is equal to:
> > + *
> > + * ⎡ xx + ww -yy - zz 2*(xy - zw) 2*(xz + yw) ⎤
> > + * U = ⎢ 2*(xy + zw) 1 - 2*(xx + zz) 2*(yz - xw) ⎥
> > + * ⎣ 2*(xz - yw) 2*(yz + xw) 1 - 2*(xx + yy) ⎦
> > + *
> > + * By identifying the coefficient in the matrix V computed above and the
> > + * matrix U obtained from converting the quaternion to 3x3 matrix, we get
> > + * the following results.
> > */
> > - memset(m, 0, 16 * sizeof(float));
> > - for (int i=0; i<4; ++i)
> > - for (int j=0; j<4; ++j)
> > - for (int k=0; k<4; ++k)
> > - for (int l=0; l<4; ++l)
> > - m[4*i+l] += y_rot[i][j] * x_rot[j][k] * z_rot[k][l];
> > +
> > + /* Rename variables and precompute square values to improve readability, as
> > + * they are used multiple times. */
> > + float qw = q[3];
> > + float qx = q[0];
> > + float qy = q[1];
> > + float qz = q[2];
> > +
> > + float sqx = qx * qx;
> > + float sqy = qy * qy;
> > + float sqz = qz * qz;
> > + float sqw = qw * qw;
> > +
> > + /* The test value is extracted from V_23 = -sin(pitch), which is also
> > + * computed as V_23 = 2 * (qw*qx + qy*qz). When abs(V_23 / 2) > = 0.4999,
> > + * the cos(pitch) will be 0 and we need a fallback method to get the other
> > + * angles in the singularity. If the quaternion is scaled, the scale is
> > + * squared because we multiply components by pair, so there's no need for
> > + * square root here when computing unit. */
> > + float test = - (qw*qx + qy*qz);
> > + float unit = sqx + sqy + sqz + sqw;
> > +
> > + /* Diagonal values. */
> > + float V_11 = 1 - 2 * (sqy + sqz);
> > + float V_22 = 1 - 2 * (sqx + sqz);
> > + float V_33 = 1 - 2 * (sqx + sqy);
> > +
> > + /* Values with a minus sign. */
> > + float V_13 = 2 * (qx * qz - qw * qy);
> > + float V_21 = 2 * (qx * qy - qw * qz);
> > +
> > + /* Values with only plus sign. */
> > + float V_23 = 2 * (qw * qx + qy * qz);
> > +
> > + /**
> > + * Original code for this part was from:
> > + * http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToEuler/
> > + **/
> > + if (test > 0.499f * unit)
> > + {
> > + /* singularity at north pole */
> > + *yaw = -asin(V_11);
> > + *pitch = - (float)M_PI_2;
> > + *roll = 0;
> > + }
> > + else if (test < -0.499f * unit)
> > + {
> > + /* singularity at south pole */
> > + *yaw = asin(V_11);
> > + *pitch = (float)M_PI_2;
> > + *roll = 0;
> > + }
> > + else
> > + {
> > + *yaw = atan2( V_13, V_33 );
> > + *pitch = atan2( -V_23, sqrt( V_21*V_21 + V_22*V_22 ) );
> > + *roll = -atan2( V_21, V_22 );
> > + }
> > +}
> > +
> > +static void EulerToQuaternion(float *q, float yaw, float pitch, float roll)
> > +{
> > + const float c_yaw = cos(yaw / 2.f);
> > + const float s_yaw = sin(yaw / 2.f);
> > + const float c_pitch = cos(pitch / 2.f);
> > + const float s_pitch = sin(pitch / 2.f);
> > + const float c_roll = cos(roll / 2.f);
> > + const float s_roll = sin(roll / 2.f);
> > +
> > + /* The values below are the result of converting each Euler rotation
> > + * into it's equivalent quaternion, and then multiplying the quaternions
> > + * together. It also means that any additional Euler rotation order
> > + * conventions can easly be added with the same method. */
> > +
> > + q[3] = c_yaw * c_pitch * c_roll - s_yaw * s_pitch * s_roll;
> > + q[0] = s_yaw * c_pitch * s_roll - c_yaw * s_pitch * c_roll;
> > + q[1] = -s_yaw * c_pitch * c_roll - c_yaw * s_pitch * s_roll;
> > + q[2] = c_yaw * c_pitch * s_roll + s_yaw * s_pitch * c_roll;
> > +}
> > +
> > +void vlc_viewpoint_to_4x4( const vlc_viewpoint_t *vp, float *m )
> > +{
> > + const float *q = vp->quat;
> > +
> > + /* Scaling factor is squared when multiplying the components
> > + * pair-wise, so we don't need to use square roots here. */
> > + const float square_norm = q[0]*q[0] + q[1]*q[1]
> > + + q[2]*q[2] + q[3]*q[3];
> > +
> > + const float xx = q[0] * q[0] / square_norm;
> > + const float yy = q[1] * q[1] / square_norm;
> > + const float zz = q[2] * q[2] / square_norm;
> > + const float ww = q[3] * q[3] / square_norm;
> > +
> > + const float xy = q[0] * q[1] / square_norm;
> > + const float zw = q[2] * q[3] / square_norm;
> > + const float yz = q[1] * q[2] / square_norm;
> > + const float xw = q[0] * q[3] / square_norm;
> > +
> > + const float xz = q[0] * q[2] / square_norm;
> > + const float yw = q[1] * q[3] / square_norm;
> > +
> > + /* The quaternion is the opposite rotation of the view.
> > + * We need to inverse the matrix at the same time. */
> > +
> > + /* ⎡ xx + ww -yy - zz 2*(xy - zw) 2*(xz + yw) ⎤
> > + * U = ⎢ 2*(xy + zw) 1 - 2*(xx + zz) 2*(yz - xw) ⎥
> > + * ⎣ 2*(xz - yw) 2*(yz + xw) 1 - 2*(xx + yy) ⎦
> > + *
> > + * Written column by column, as a usual transformation matrix in the
> > + * projective space,
> > + *
> > + * M = [ R ][ T ]
> > + * [ 0 ][ 1 ]
> > + *
> > + * with R being the rotation matrix from the quaternion and T being a
> > + * translation vector, { 0 } here.
> > + **/
> > +
> > + m[0] = xx + ww - yy - zz;
> > + m[1] = 2 * (xy + zw);
> > + m[2] = 2 * (xz - yw);
> > + m[3] = 0;
> > +
> > + m[4] = 2 * (xy - zw);
> > + m[5] = 1 - 2 * (xx + zz);
> > + m[6] = 2 * (yz + xw);
> > + m[7] = 0;
> > +
> > + m[8] = 2 * (xz + yw);
> > + m[9] = 2 * (yz - xw);
> > + m[10] = 1 - 2 * (xx + yy);
> > + m[11] = 0;
> > +
> > + m[12] = m[13] = m[14] = 0;
> > + m[15] = 1;
> > }
> >
> > void vlc_viewpoint_from_euler(vlc_viewpoint_t *vp,
> > float yaw, float pitch, float roll)
> > {
> > - vp->yaw = -yaw;
> > - vp->pitch = -pitch;
> > - vp->roll = -roll;
> > + /* convert angles from degrees into radians */
>
> You could do the normalization here.
>
> > + yaw = -yaw * (float)M_PI / 180.f;
> > + pitch = -pitch * (float)M_PI / 180.f;
> > + roll = -roll * (float)M_PI / 180.f;
> > +
> > + EulerToQuaternion(vp->quat, yaw, pitch, roll);
> > }
> >
> > void vlc_viewpoint_to_euler(const vlc_viewpoint_t *vp,
> > float *yaw, float *pitch, float *roll)
> > {
> > - *yaw = -vp->yaw;
> > - *pitch = -vp->pitch;
> > - *roll = -vp->roll;
> > + QuaternionToEuler(yaw, pitch, roll, vp->quat);
> > +
> > + /* convert angles from radian into degrees */
> > + *yaw = -180.f / (float)M_PI * (*yaw);
> > + *pitch = -180.f / (float)M_PI * (*pitch);
> > + *roll = -180.f / (float)M_PI * (*roll);
> > }
> > --
> > 2.31.0
> > _______________________________________________
> > vlc-devel mailing list
> > To unsubscribe or modify your subscription options:
> > https://mailman.videolan.org/listinfo/vlc-devel
> >
> _______________________________________________
> vlc-devel mailing list
> To unsubscribe or modify your subscription options:
> https://mailman.videolan.org/listinfo/vlc-devel
More information about the vlc-devel
mailing list