Точность векторных операций DirectXMath

У меня странные результаты функции XMVector3AngleBetweenVectors. Рассмотрим этот код:

float angle = XMConvertToDegrees(XMVectorGetX(
        XMVector3AngleBetweenVectors(GMathFV(XMFLOAT3(0.0f, 100.0f, 0.0f)), 
        GMathFV(XMFLOAT3(0.0f, 200.0f, 0.0f)))));

Он ищет угол между двумя трехмерными векторами, описанными структурами XMFLOAT3. GMathFV — определяемая пользователем функция, которая преобразует XMFLOAT3 в XMVECTOR следующим образом:

inline XMVECTOR GMathFV(XMFLOAT3& val)
{
    return XMLoadFloat3(&val);    
}

Все остальное — это библиотека directxmath.h. Здесь все в порядке, и результирующий угол равен 0,00000, как и ожидалось.

Но для других векторов с отрицательным значением оси Y, например:

float angle = XMConvertToDegrees(XMVectorGetX(
        XMVector3AngleBetweenVectors(GMathFV(XMFLOAT3(0.0f, -100.0f, 0.0f)), 
        GMathFV(XMFLOAT3(0.0f, -99.0f, 0.0f)))));

Результат 0,0197823402, что я вряд ли могу назвать нулевым углом.

Пожалуйста, помогите мне разобраться в проблеме. Это отрицательная числовая точность, слишком близкие векторные координаты или, может быть, что-то еще?

UPD: Удивительно, но для a(0.0f, 100.0f, 0.0f) x b(0.0f, 99.0f, 0.0f) выдает 0,0197823402, а для a(0.0f, 101.0f, 0.0f) x b(0.0f, 100.0f, 0.0f) 0,000000


person GuardianX    schedule 11.02.2013    source источник
comment
Использует ли DirectXMath радианы или градусы?   -  person Xathereal    schedule 11.02.2013
comment
@Xathereal использует радианы.   -  person GuardianX    schedule 11.02.2013
comment
DirectX по умолчанию устанавливает низкоточный флаг FPU, что может вызвать проблемы, которые вы видите. Чтобы указать ему сохранить флаги FPU внутри вызовов DX, вызовите CreateDevice с D3DCREATE_FPU_PRESERVE в поле поведения или для .Net используйте CreateFlags.FpuPreserve. Это помогает? Следите за тем, чтобы ваши флаги FPU вызывали исключения в некоторых обстоятельствах.   -  person David    schedule 11.02.2013
comment
@DavidM спасибо за ответ. Похоже, что устройство D3D11 уже настроено на двойную точность, так как установить этот флаг невозможно.   -  person GuardianX    schedule 11.02.2013


Ответы (1)


DirectXMath предназначен для 32-битных вычислений с плавающей запятой. Вы видите эскалацию ошибки с плавающей запятой. Вот определение XMVector3AngleBetweenVectors.

inline XMVECTOR XM_CALLCONV XMVector3AngleBetweenVectors(FXMVECTOR V1, FXMVECTOR V2)
{
    XMVECTOR L1 = XMVector3ReciprocalLength(V1);
    XMVECTOR L2 = XMVector3ReciprocalLength(V2);

    XMVECTOR Dot = XMVector3Dot(V1, V2);

    L1 = XMVectorMultiply(L1, L2);

    XMVECTOR CosAngle = XMVectorMultiply(Dot, L1);
    CosAngle = XMVectorClamp(CosAngle, g_XMNegativeOne.v, g_XMOne.v);

    return XMVectorACos(CosAngle);
}

В вашем первом примере CosAngle равен 1.000000000

Во втором примере CosAngle равен 0,999999940.

XMVectorACos (0,999999940) = 0,000345266977

Эта большая ошибка возникает из-за полиномиальной аппроксимации ACos. В общем, вы должны избегать тригонометрических инверсий, когда это возможно. Они медленные и шумные. Вот определение, чтобы вы могли получить представление о его размере.

inline XMVECTOR XM_CALLCONV XMVectorACos (FXMVECTOR V)
{
    __m128 nonnegative = _mm_cmpge_ps(V, g_XMZero);
    __m128 mvalue = _mm_sub_ps(g_XMZero, V);
    __m128 x = _mm_max_ps(V, mvalue);  // |V|

    // Compute (1-|V|), clamp to zero to avoid sqrt of negative number.
    __m128 oneMValue = _mm_sub_ps(g_XMOne, x);
    __m128 clampOneMValue = _mm_max_ps(g_XMZero, oneMValue);
    __m128 root = _mm_sqrt_ps(clampOneMValue);  // sqrt(1-|V|)

    // Compute polynomial approximation
    const XMVECTOR AC1 = g_XMArcCoefficients1;
    XMVECTOR vConstants = XM_PERMUTE_PS( AC1, _MM_SHUFFLE(3, 3, 3, 3) );
    __m128 t0 = _mm_mul_ps(vConstants, x);

    vConstants = XM_PERMUTE_PS( AC1, _MM_SHUFFLE(2, 2, 2, 2) );
    t0 = _mm_add_ps(t0, vConstants);
    t0 = _mm_mul_ps(t0, x);

    vConstants = XM_PERMUTE_PS( AC1, _MM_SHUFFLE(1, 1, 1, 1) );
    t0 = _mm_add_ps(t0, vConstants);
    t0 = _mm_mul_ps(t0, x);

    vConstants = XM_PERMUTE_PS( AC1, _MM_SHUFFLE(0, 0, 0, 0) );
    t0 = _mm_add_ps(t0, vConstants);
    t0 = _mm_mul_ps(t0, x);

    const XMVECTOR AC0 = g_XMArcCoefficients0;
    vConstants = XM_PERMUTE_PS( AC0, _MM_SHUFFLE(3, 3, 3, 3) );
    t0 = _mm_add_ps(t0, vConstants);
    t0 = _mm_mul_ps(t0, x);

    vConstants = XM_PERMUTE_PS( AC0, _MM_SHUFFLE(2, 2, 2, 2) );
    t0 = _mm_add_ps(t0, vConstants);
    t0 = _mm_mul_ps(t0, x);

    vConstants = XM_PERMUTE_PS( AC0, _MM_SHUFFLE(1, 1, 1, 1) );
    t0 = _mm_add_ps(t0, vConstants);
    t0 = _mm_mul_ps(t0, x);

    vConstants = XM_PERMUTE_PS( AC0, _MM_SHUFFLE(0, 0, 0, 0) );
    t0 = _mm_add_ps(t0, vConstants);
    t0 = _mm_mul_ps(t0, root);

    __m128 t1 = _mm_sub_ps(g_XMPi, t0);
    t0 = _mm_and_ps(nonnegative, t0);
    t1 = _mm_andnot_ps(nonnegative, t1);
    t0 = _mm_or_ps(t0, t1);
    return t0;
}
person gradbot    schedule 13.05.2014