[vlc-devel] [PATCH 18/18] viewpoint: use quaternion instead of euler angles

Steve Lhomme robux4 at ycbcr.xyz
Thu Apr 1 06:07:28 UTC 2021


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.

> ---
>   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
> 


More information about the vlc-devel mailing list