2017年10月4日 星期三

atan() 和 atan2() 快速算法

atan() 計算的時候如何得知正確角度;fastAtan2f() 快速計算方法

atan() 在計算的時候是沒有把dx和dy的正負號考慮進去的,也就是說在第3象限的時候會被當作第一象限來處理。
解決方法很容易,有一個已經幫你處好的函式是 atan2() 可以直接輸入 atan2(dy, dx) 計算出正確的結果,不過出來的角度是 -PI~PI 逕度還需要處理一下轉為角度。
這裡提供一個已經處理好上述而且速度更快的方法,是來自於OpenCV內的函式原封不動複製出來的應該足夠安全。

atan2()

static inline float fastAtan2f(float dy, float dx){
    // 快速atan運算
    static const float atan2_p1 = 0.9997878412794807f*(float)(180/M_PI);
    static const float atan2_p3 = -0.3258083974640975f*(float)(180/M_PI);
    static const float atan2_p5 = 0.1555786518463281f*(float)(180/M_PI);
    static const float atan2_p7 = -0.04432655554792128f*(float)(180/M_PI);
    static const float atan2_DBL_EPSILON = 2.2204460492503131e-016;

    float ax = std::abs(dx), ay = std::abs(dy);
    float a, c, c2;
    if (ax >= ay) {
        c = ay/(ax + static_cast<float>(atan2_DBL_EPSILON));
        c2 = c*c;
        a = (((atan2_p7*c2 + atan2_p5)*c2 + atan2_p3)*c2 + atan2_p1)*c;
    } else {
        c = ax/(ay + static_cast<float>(atan2_DBL_EPSILON));
        c2 = c*c;
        a = 90.f - (((atan2_p7*c2 + atan2_p5)*c2 + atan2_p3)*c2 + atan2_p1)*c;
    }
    if (dx < 0)
        a = 180.f - a;
    if (dy < 0)
        a = 360.f - a;
    return a;
}
fastAtan2f() 計算出來就是 0~360 度的角度了,如果需要逕度從代碼內把轉換的移除即可。
M_PI 可能需要自己定義一下,或者引入cmath並啟用宏定義(啟用是一行宏定義的代碼需要再自己查)。
如果只是單單需要拍自己定一個比較快~
#define M_PI 3.14159265358979323846

atan()

這個OpenCV裡面就沒有,不過可以從上面的 atan2() 修改一下即可(簡單來說就是x設定成1)。
這個比較常用到逕度版本的我貼逕度版本的
float fastAtanf_rad(float dy){
    static const float atan2_p1 =  0.9997878412794807f;
    static const float atan2_p3 = -0.3258083974640975f;
    static const float atan2_p5 =  0.1555786518463281f;
    static const float atan2_p7 = -0.04432655554792128f;
    static const float atan2_DBL_EPSILON = 2.2204460492503131e-016;

    float ax = 1.0, ay = std::abs(dy);
    float a, c, c2;
    if (ax >= ay) {
        c = ay/(ax + static_cast<float>(atan2_DBL_EPSILON));
        c2 = c*c;
        a = (((atan2_p7*c2 + atan2_p5)*c2 + atan2_p3)*c2 + atan2_p1)*c;
    } else {
        c = ax/(ay + static_cast<float>(atan2_DBL_EPSILON));
        c2 = c*c;
        a = M_PI*0.5 - (((atan2_p7*c2 + atan2_p5)*c2 + atan2_p3)*c2 + atan2_p1)*c;
    }

    if (dy < 0)
        a = - a;

    return a;
}
這一分我有做簡單的驗證是沒問題的,可以參照底下的代碼有過程,敏感場合建議自己在想方法驗證一次~

其他

全部一共4個版本,atan()和atan2() 以及他們的 逕度版本與角度版本。

沒有留言:

張貼留言