math - 複数の四元数の「平均」?

Question

OpenGL プログラムでスケルトン アニメーションのマトリックスからクォータニオンに切り替えようとしていますが、問題が発生しました。

単位クォータニオンの数が与えられた場合、ベクトルを変換するために使用すると、各クォータニオンによって個別に変換されたベクトルの平均であるベクトルが得られるクォータニオンを取得する必要があります。(行列の場合、単純に行列を足し合わせて、行列の数で割ります)

Answer 1

コンピュータ グラフィックス業界で一般に信じられていることとは反対に、この問題を解決するための堅牢で正確、かつ単純なアルゴリズムが航空宇宙業界から提供されています。これは、平均化されるクォータニオンの数に (大きな) 定数係数を加えた時間線形で実行されます。

Q = [ a ₁ q ₁ , a ₂ q ₂ , ..., a n _q n _] ,

ここで、a _iはi番目の四元数の重みで、q _iは列ベクトルとして平均化される i番目の四元数です。したがって、 Qは 4 × N行列です。

QQ ^Tの最大固有値に対応する正規化された固有ベクトルは加重平均です。QQ ^Tは自己共役であり、少なくとも正の半正定であるため、その固有問題を解くための高速で堅牢な方法が利用可能です。行列と行列の積を計算することは、平均化される要素の数とともに増加する唯一のステップです。

Journal of Guidance, Control, and Dynamics from 2007 のこのテクニカル ノートを参照してください。これは、この方法と他の方法の要約論文です。現代では、上で引用した方法は、実装の信頼性と堅牢性の良いトレードオフとなり、1978 年には既に教科書に掲載されていました。

Answer 2

残念ながら、それを行うのはそれほど簡単ではありませんが、可能です。その背後にある数学を説明するホワイトペーパーは次のとおりです。

Unity3D Wiki ページをチェックしてください (以下のコード サンプルは同じ記事からのものです): http://wiki.unity3d.com/index.php/Averaging_Quaternions_and_Vectors

//Get an average (mean) from more then two quaternions (with two, slerp would be used).
//Note: this only works if all the quaternions are relatively close together.
//Usage: 
//-Cumulative is an external Vector4 which holds all the added x y z and w components.
//-newRotation is the next rotation to be added to the average pool
//-firstRotation is the first quaternion of the array to be averaged
//-addAmount holds the total amount of quaternions which are currently added
//This function returns the current average quaternion
public static Quaternion AverageQuaternion(ref Vector4 cumulative, Quaternion newRotation, Quaternion firstRotation, int addAmount){

    float w = 0.0f;
    float x = 0.0f;
    float y = 0.0f;
    float z = 0.0f;

    //Before we add the new rotation to the average (mean), we have to check whether the quaternion has to be inverted. Because
    //q and -q are the same rotation, but cannot be averaged, we have to make sure they are all the same.
    if(!Math3d.AreQuaternionsClose(newRotation, firstRotation)){

        newRotation = Math3d.InverseSignQuaternion(newRotation);    
    }

    //Average the values
    float addDet = 1f/(float)addAmount;
    cumulative.w += newRotation.w;
    w = cumulative.w * addDet;
    cumulative.x += newRotation.x;
    x = cumulative.x * addDet;
    cumulative.y += newRotation.y;
    y = cumulative.y * addDet;
    cumulative.z += newRotation.z;
    z = cumulative.z * addDet;      

    //note: if speed is an issue, you can skip the normalization step
    return NormalizeQuaternion(x, y, z, w);
}

public static Quaternion NormalizeQuaternion(float x, float y, float z, float w){

    float lengthD = 1.0f / (w*w + x*x + y*y + z*z);
    w *= lengthD;
    x *= lengthD;
    y *= lengthD;
    z *= lengthD;

    return new Quaternion(x, y, z, w);
}

//Changes the sign of the quaternion components. This is not the same as the inverse.
public static Quaternion InverseSignQuaternion(Quaternion q){

    return new Quaternion(-q.x, -q.y, -q.z, -q.w);
}

//Returns true if the two input quaternions are close to each other. This can
//be used to check whether or not one of two quaternions which are supposed to
//be very similar but has its component signs reversed (q has the same rotation as
//-q)
public static bool AreQuaternionsClose(Quaternion q1, Quaternion q2){

    float dot = Quaternion.Dot(q1, q2);

    if(dot < 0.0f){

        return false;                   
    }

    else{

        return true;
    }
}

この投稿も: http://forum.unity3d.com/threads/86898-Average-quaternions

Answer 3

ここで提案されているようにクォータニオンを Slerping しようとしましたが、それは私がやろうとしていることではうまくいきませんでした (モデルは歪んでいました)。より良い解決策）。

Answer 4

四元数を追加することはできません。できることは、途中を含む 2 つの角度の間で連続的に回転するクォータニオンを見つけることです。クォータニオン補間は「slerp」として知られており、ウィキペディアのページがあります。これは、アニメーションに非常に役立つトリックです。いくつかの点で、コンピュータ グラフィックスで四元数を使用する主な理由は slerp です。

Answer 5

クォータニオンでも同じことができますが、小さな修正があります。2. ライブラリが単位四元数で動作する場合は、平均化の終わりである平均四元数を正規化します。

平均クォータニオンは、おおよその平均回転を表します (最大誤差は約 5 度)。

警告: 回転があまりにも異なる場合、異なる方向の平均行列が壊れる可能性があります。

Answer 6

クォータニオンは、制約のない平均を計算するときに回転に使用する DOF の理想的なセットではありません。

私がよく使うのはこちら(

[MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static Vector3 ToAngularVelocity( this Quaternion q )
    {
        if ( abs(q.w) > 1023.5f / 1024.0f)
            return new Vector3();
            var angle = acos( abs(q.w) );
            var gain = Sign(q.w)*2.0f * angle / Sin(angle);

        return new Vector3(q.x * gain, q.y * gain, q.z * gain);
    }


    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static Quaternion FromAngularVelocity( this Vector3 w )
    {
        var mag = w.magnitude;
        if (mag <= 0)
            return Quaternion.identity;
        var cs = cos(mag * 0.5f);
        var siGain = sin(mag * 0.5f) / mag;
        return new Quaternion(w.x * siGain, w.y * siGain, w.z * siGain, cs);

    }

    internal static Quaternion Average(this Quaternion refence, Quaternion[] source)
    {

        var refernceInverse = refence.Inverse();
        Assert.IsFalse(source.IsNullOrEmpty());
        Vector3 result = new Vector3();
        foreach (var q in source)
        {
            result += (refernceInverse*q).ToAngularVelocity();
        }

        return reference*((result / source.Length).FromAngularVelocity());
    }
     internal static Quaternion Average(Quaternion[] source)
    {
        Assert.IsFalse(source.IsNullOrEmpty());
        Vector3 result = new Vector3();
        foreach (var q in source)
        {
            result += q.ToAngularVelocity();
        }

        return (result / source.Length).FromAngularVelocity();
    }
     internal static Quaternion Average(Quaternion[] source, int iterations)
    {
        Assert.IsFalse(source.IsNullOrEmpty());
        var reference = Quaternion.identity;
        for(int i = 0;i < iterations;i++)
        {
            reference = Average(reference,source);

        }
        return reference;

    }`

Answer 7

クォータニオンを平均化するときにおそらく実際に必要な別のスラープベースの方法があります。

まず、固有分析ベースの平均化と比較してみましょう。

重み w_A = 2 および w_B = 1 を使用して、X 軸の周りの 0° および 90° の回転に対応する 2 つのクォータニオン A および B の平均化を検討してください。予想される加重平均は、X 軸の周りの 30° の回転に対応する必要があります。 . Slerp ベースの加重平均は期待値を返します。固有分析に基づく加重平均は、26.56° の回転を返します。

固有解析ベースの方法は、対応する回転行列のフロベニウス ノルムを最小化するクォータニオンを返します。代わりに、slerp ベースの方法は、クォータニオンに関して 3D 回転空間の平均を計算します。

import math
import numpy as np
import quaternion # (pip install numpy-quaternion)
d2r = math.pi/180
r2d = 180/math.pi

def recover_angle(mat):
    return np.arctan2(mat[1,0], mat[0,0])

# ground truth
angles = np.array([0,90])
weights = np.array([2,1])
mean_angle = np.sum(angles*(weights/weights.sum()))
quaternion_mean = quaternion.from_euler_angles(mean_angle*d2r,0,0)

# eigen analysis
Q = quaternion.as_float_array(
    [
        (weight**0.5) * quaternion.from_euler_angles(0,0,angle*d2r) 
        for angle,weight 
        in zip(angles,weights)
    ]
).T
QQT = Q @ Q.T
eigval,eigvec = np.linalg.eig(QQT)
quaternion_mean_eig = quaternion.from_float_array( eigvec[:,np.argmax(eigval)] )

# slerp based
def slerp_avg(quaternions, weights):
    # welford's mean in terms of linear mix operations:
    toqt = quaternion.from_float_array
    mix = lambda a,b,k: quaternion.slerp(toqt(a),toqt(b),0,1,k)
    wmean, wsum, num = quaternions[0], weights[0], len(weights)
    for i in range(1,num):
        wsum += weights[i]
        k = (weights[i]/wsum)
        # wmean := wmean*(1-k) + quaternions[i]*k
        wmean = mix(wmean,quaternions[i],k) 
    return wmean

quaternion_mean_slerp = slerp_avg(quaternion.as_float_array(
    [quaternion.from_euler_angles(0,0,angle*d2r) for angle in angles]
), weights)

mat_mean = quaternion.as_rotation_matrix(quaternion_mean)
mat_mean_eig = quaternion.as_rotation_matrix(quaternion_mean_eig)
mat_mean_slerp = quaternion.as_rotation_matrix(quaternion_mean_slerp)

print("expected", recover_angle(mat_mean)*r2d)
print("eigen", recover_angle(mat_mean_eig)*r2d)
print("slerp", recover_angle(mat_mean_slerp)*r2d)

出力:

expected 29.999999999999996
eigen 26.565051177077994
slerp 30.00000000000001

このためのゼロ依存 C++ ライブラリは、 https://github.com/xaedes/average_affine_transform_mat/にあります。

Answer 8

これは、この提案されたアルゴリズムの実装を含む GitHub リポジトリです:) https://github.com/christophhagen/averaging-quaternions

christophhagen ofc に感謝とクレジット ;)