プロジェクション行列から内部パラメータを得る

今回の内容

• 内部パラメータを加味したプロジェクション行列の計算方法
• そこから逆算して，プロジェクション行列から内部パラメータを計算する方法

経緯

コンピュータビジョン (CV) 屋さんにとってみれば，カメラは事前に校正しておくものだと考えて内部パラメータは事前計算して，OpenGLで描画するときには内部パラメータからプロジェクション行列を計算する．

一方で，ユーザインタフェース (UI) 屋さんにとっては，プロジェクション行列は手に入るが内部パラメータが分からない，ということがあるらしい．どういう状況だ？と思いますが，HoloLensからプロジェクション行列が手に入るが，内部パラメータは分からないということがある，と言うのである．つまり，以下のような状況だ．

• CV屋さん: 内部パラメータ→プロジェクション行列
• UI屋さん: プロジェクション→内部パラメータ

図 Microsoft社製HoloLens（公式サイトより）

内部パラメータを加味したプロジェクション行列の計算方法

これはどこからか拝借したもののはず．思い出せないです... OpenCVで計算したカメラ位置姿勢は基本的にY-downなので最後のパラメータはtrueをデフォルト引数にしている．

glm::mat4 calcCalibratedProjectionMatrix(
    float nearPlane, float farPlane,
    int w, int h,
    float cx, float cy,
    float fx, float fy,
    float skew,
    bool ydown = true
)
{
    glm::mat4 projMat;

    projMat[0][0] = float(2.0 * fx / w);
    projMat[1][0] = float(-2.0 * skew / w);
    projMat[2][0] = float(w - 2.0 * cx) / w;
    projMat[3][0] = 0.0f;

    projMat[0][1] = 0.0f;
    if (ydown)
    {
        projMat[1][1] = float(2.0 * fy / h);
        projMat[2][1] = float(-h + 2.0 * cy) / h;
    }
    else
    {
        projMat[1][1] = float(-2.0 * fy / h);
        projMat[2][1] = float(h - 2.0 * cy) / h;
    }
    projMat[3][1] = 0.0f;

    projMat[0][2] = 0.0f;
    projMat[1][2] = 0.0f;
    projMat[2][2] = (-farPlane - nearPlane) / (farPlane - nearPlane);
    projMat[3][2] = float(-2.0 * farPlane * nearPlane / (farPlane - nearPlane));

    projMat[0][3] = 0.0f;
    projMat[1][3] = 0.0f;
    projMat[2][3] = -1.0f;
    projMat[3][3] = 0.0f;

    return projMat;
}

プロジェクション行列から内部パラメータを計算する方法

先のプログラムにある計算をもとに逆算する．つまり，以下のようになる．

void calcIntrinsicsFromProjectionMatrix(
    const glm::mat4 &projMat,
    int width, int height,
    float &fx, float &fy, float &cx, float &cy, float &skew,
    bool ydown = true
    )
{
    fx = projMat[0][0] * width / 2.0f;
    fy = projMat[1][1] * height / 2.0f;
    if (!ydown) fy *= -1.0f;
    cx = -(projMat[2][0] * width - width) / 2.0f;
    if (ydown) cy = (projMat[2][1] * height + height) / 2.0f;
    else cy = (height - projMat[2][1] * height) / 2.0f;
    skew = -projMat[1][0] * width / 2.0;
}