glMultMatrixを用いたHomography変換

OpenGLのglMultMatrix関数を使えば、ピンホールカメラモデルであれば任意の投影が扱える。そんな原理は頭では理解しているけれど、いざ実際にプログラムを作ろうとすると色々分からないことがあることに気づく。殆どはOpenGLを正確に理解していないことが原因だった。いったい何を理解していなかったのか？

何がしたかったのか？

プロジェクタを用いて平面に正対した像を投影するプログラムを作ろうとした。平面からプロジェクタへの投影行列であるHomographyは、OpenCV等を用いて正しい物が簡単に求まったにも関わらず、OpenGLのglTranslateとglMultMatrixを用いて平面上の任意の位置にあるマークを表示しようとして時間がかかってしまった。

配列と行列の関係

まずは、glMultMatrixd関数を使う上で気をつけること。この関数は、glMultMatrixd(GLdouble *)となっているので、配列のポインタを渡す必要がある。画像を配列で表す場合は、通常スキャンライン（水平）方向に画素をたどって、順に配列に格納するよう実装するが、4x4の行列を配列で表す場合もこうすることが多いように思う。OpenGLは特別で、下記のように行方向（垂直）に要素をたどり配列に格納する。ただし、これは点の位置をベクトルで表す場合に、縦ベクトル（左から変換行列を掛けるとする）か横ベクトルにする（右側から変換行列を掛ける）かの違いかもしれない。いつもこれを考えるとき混乱します。頼むから統一してくれ・・・

// 行列の定義はこれ。
cv::Mat_<double> m(4, 4);

double m_tmp[16] = {
    m(0, 0),    m(1, 0),    m(2, 0),    m(3, 0),
    m(0, 1),    m(1, 1),    m(2, 1),    m(3, 1),
    m(0, 2),    m(1, 2),    m(2, 2),    m(3, 2),
    m(0, 3),    m(1, 3),    m(2, 3),    m(3, 3)
};
glMultMatrixd(m_tmp);

「行列の積の順序」と「関数の記述順序」の関係

例えば、ユークリッド変換（回転＋並進）は回転行列$R$と並進ベクトル$\mathbf{t}$で次のように書ける。 \[ \mathbf{x}'=R\mathbf{x} + \mathbf{t} \] さらに並進を同次座標表現の行列で表すとこうなる。 \[ \begin{bmatrix}\mathbf{x}'\\1\end{bmatrix}^T= \begin{bmatrix} O& \mathbf{t}\\ \mathbf{0}^T& 1 \end{bmatrix} \begin{bmatrix} R& \mathbf{0}\\ \mathbf{0}^T& 1 \end{bmatrix} \begin{bmatrix} \mathbf{x}\\ 1 \end{bmatrix} \] このとき、式中では左から並進、回転、$\mathbf{x}$の順に現れるが、コードでも下記のようにこの順に書く。計算する順序は、まず回転行列を掛けて、並進行列を掛けるので、イメージ的に逆にも考えることができて、これまた混乱することが多い。

glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glTranslated(t[0], t[1], t[2]);
glRotated(theta, r[0], r[1], r[2]);
glVertex4d(x[0], x[1], x[2], 1.0);

GL_MODELVIEWとGL_PROJECTIONどっち?

OpenGLなどの3D-CG可視化ツールは、カメラの外部パラメータを表すModelView行列とカメラの内部パラメータに相当するProjection行列に分かれている。しかも両者4x4の行列でほぼ同じ。しかし、目的のHomographyはそれら両方を内包している行列なので、どちらにHomography行列をLoadすればいいのか少し迷った。単純な方法は、Projection行列を正射影（単位行列）にし、ModelView行列にHomography行列を代入する方法。mは行列の各要素が格納された配列へのポインタである。

glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glMultMatrixd(m);
glVertex4d(x, y, z, 1.0);

実際のコード

次のコードはC++でOpenCVを使用している。

// homography行列の定義はこれ。値は何からの方法で計算する。
cv::Mat_<double> h;

// ウィンドウ全体をビューポートにする
glViewport(0, 0, w, h);
 
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
 
glOrtho(-0.0f, w, h, 0.0f, -1.0f, 1000.0f);

//3x3のhomography行列を4x4の行列に代入する（結局$z$座標が無効化されいてるだけ）
double h_tmp[16] = {
    h(0, 0),    h(1, 0),    0,    h(2, 0),
    h(0, 1),    h(1, 1),    0,    h(2, 1),
    0,          0,          0,    0,
    h(0, 2),    h(1, 2),    0,    h(2, 2)
};
 
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glMultMatrixd(h_tmp);