デジタル・デザイン・ラボラトリーな日々

アラフィフプログラマーが数学と物理を基礎からやり直す。https://qiita.com/yaju

射影変換(ホモグラフィ)について理解してみる その5

前回までは、代数ライブラリsylvesterを使用して行列計算を求めてきましたが、今回は代数ライブラリを使わずに8次元連立方程式を4x4行列で求める方法を説明します。

ARToolKitを参考に様々な言語に実装し直したARToolKit互換のクラスライブラリNyARToolKitの開発者@nyatlaさんが射影変換の高速化する上で、8x8行列を4x4行列と2x2行列を組み合わせて解く方法を考案されました。その資料(PDF形式)が下記になります。
遠近法の射影変換パラメータ計算の高速化

数学に疎い私には、8次元連立方程式を4x4行列と2x2行列を組み合わせて解くことが出来るなんて、考えもつきませんでした。
3Dライブラリ内で使われる行列ライブラリをみると4x4行列クラスまでは作成されているので、これが出来るんであれば、射影変換がより身近になると思っています。
テクスチャマッピングを理解してみる」で書いたのですが、台形に変形するのに2つの三角形を組み合わせて実現しています。もちろんポリゴンオブジェクトを表示するには、三角形を組み合わせるのが通常なのでいいんですが、ツールとして三角形に分割せずに変形させたいなと思っていました。

射影変換パラメータ計算の高速化をJavascriptで実現しようと目論んでみたわけです。NyARToolKitは、Flashにも対応しているのでActionscriptからJavascriptへのコンバート程度なら簡単に出来るんじゃないかと楽観的な気持ちでした。
しかし、該当するソースリスト(getParam_5)を見ると、もともとAR(拡張現実)用でマーカーのQRコードの歪みを補正するのが目的なので、それ用に最適化されており、PDFの資料と見比べても私の理解力では何が何だか分からないと、もーとほほな状態でした。

それでも資料内容をなんとか理解しながら組んでみました。
実行してみると画像が少しだけ表示されたんです、どこかミスってるかなーと資料とソースを見比べてみたんですが、合ってるんですよね。で、まじまじと資料の(8)の展開したのを見てみると、XとYで番号違いがあることを発見しました。例えばtx21となるところがtx12となっていたりと、B,G,Hのところで複数の記述ミスがありました。(この件はTwitterでnyatlaさんに報告して資料を修正して頂きました。)
それを正しく直したところ、射影変換された画像が正常に表示されたのです。この時はおもわずガッツポーズしてしまいました。

下記ソースリストが、資料内容を参考にした8個のパラメータを求める処理になります。
内部で使用している行列クラスは自作のものを使用しております。

function getParam(src, dest) {

    function Z(val){
        return val == 0 ? 0.5 : val;
    }

    var X1 = Z(src[0][0]);
    var X2 = Z(src[1][0]);
    var X3 = Z(src[2][0]);
    var X4 = Z(src[3][0]);
    var Y1 = Z(src[0][1]);
    var Y2 = Z(src[1][1]);
    var Y3 = Z(src[2][1]);
    var Y4 = Z(src[3][1]);
    var x1 = Z(dest[0][0]);
    var x2 = Z(dest[1][0]);
    var x3 = Z(dest[2][0]);
    var x4 = Z(dest[3][0]);
    var y1 = Z(dest[0][1]);
    var y2 = Z(dest[1][1]);
    var y3 = Z(dest[2][1]);
    var y4 = Z(dest[3][1]);

    //X座標
    var tx = new Matrix44();
    tx.set(X1, Y1, -X1 * x1, -Y1 * x1, X2, Y2, -X2 * x2, -Y2 * x2, X3, Y3, -X3 * x3, -Y3 * x3, X4, Y4, -X4 * x4, -Y4 * x4);
    tx.Invert();
    var kx1 = tx.m11 * x1 + tx.m12 * x2 + tx.m13 * x3 + tx.m14 * x4;
    var kc1 = tx.m11 + tx.m12 + tx.m13 + tx.m14;
    var kx2 = tx.m21 * x1 + tx.m22 * x2 + tx.m23 * x3 + tx.m24 * x4;
    var kc2 = tx.m21 + tx.m22 + tx.m23 + tx.m24;
    var kx3 = tx.m31 * x1 + tx.m32 * x2 + tx.m33 * x3 + tx.m34 * x4;
    var kc3 = tx.m31 + tx.m32 + tx.m33 + tx.m34;
    var kx4 = tx.m41 * x1 + tx.m42 * x2 + tx.m43 * x3 + tx.m44 * x4;
    var kc4 = tx.m41 + tx.m42 + tx.m43 + tx.m44;

    //Y座標
    var ty = new Matrix44();
    ty.set(X1, Y1, -X1 * y1, -Y1 * y1, X2, Y2, -X2 * y2, -Y2 * y2, X3, Y3, -X3 * y3, -Y3 * y3, X4, Y4, -X4 * y4, -Y4 * y4);
    ty.Invert();
    var ky1 = ty.m11 * y1 + ty.m12 * y2 + ty.m13 * y3 + ty.m14 * y4;
    var kf1 = ty.m11 + ty.m12 + ty.m13 + ty.m14;
    var ky2 = ty.m21 * y1 + ty.m22 * y2 + ty.m23 * y3 + ty.m24 * y4;
    var kf2 = ty.m21 + ty.m22 + ty.m23 + ty.m24;
    var ky3 = ty.m31 * y1 + ty.m32 * y2 + ty.m33 * y3 + ty.m34 * y4;
    var kf3 = ty.m31 + ty.m32 + ty.m33 + ty.m34;
    var ky4 = ty.m41 * y1 + ty.m42 * y2 + ty.m43 * y3 + ty.m44 * y4;
    var kf4 = ty.m41 + ty.m42 + ty.m43 + ty.m44;
    var det_1 = kc3 * (-kf4) - (-kf3) * kc4;

    if(det_1 == 0) {
        det_1 = 0.0001;
    }
    det_1 = 1 / det_1;
    var param = new Array(8);
    var C = (-kf4 * det_1) * (kx3 - ky3) + (kf3 * det_1) * (kx4 - ky4);
    var F = (-kc4 * det_1) * (kx3 - ky3) + (kc3 * det_1) * (kx4 - ky4);
    param[2] = C;
    param[5] = F;
    param[6] = kx3 - C * kc3;
    param[7] = kx4 - C * kc4;
    param[0] = kx1 - C * kc1;
    param[1] = kx2 - C * kc2;
    param[3] = ky1 - F * kf1;
    param[4] = ky2 - F * kf2;

    return param;
}

※getParam関数内のクロージャ関数である「Z」は何をしているかというと、逆行列を求める際に対角成分が0の場合、処理の途中でゼロ除算が発生して値が求まらないため、座標に0がある場合は0.5にすることで回避しています。
微妙に違いが出る分は画像補完でなんとかなるかなと思ってます。


射影変換パラメータによる描画処理

function draw(param) {
    var ctx = this.context;
    var imgW = this.image.width;
    var imgH = this.image.height;
    var input = ctx.getImageData(0, 0, imgW, imgH);
    var output = ctx.createImageData(imgW, imgH);
    for(var i = 0; i < imgH; ++i) {
        for(var j = 0; j < imgW; ++j) {
            var tmp = j * param[6] + i * param[7] + param[8];
            var tmpX = (j * param[0] + i * param[1] + param[2]) / tmp;
            var tmpY = (j * param[3] + i * param[4] + param[5]) / tmp;
            var floorX = tmpX | 0;
            var floorY = tmpY | 0;

            if (floorX >= 0 && floorX < imgW && floorY >= 0 && floorY < imgH) {
                var pixelData = this.getPixel(input, floorX, floorY);
                var R = pixelData.R;
                var G = pixelData.G;
                var B = pixelData.B;
                this.setPixel(output, j, i, R, G, B, 255);
            }
        }
    }
    ctx.putImageData(output, 512, 0);
};


射影変換パラメータを取得し、描画用に転送元座標と転送先座標の逆行列を求めて描画する。
逆行列を求めるのにパラメータが8個なので3x3の行列でもいいのですが、大は小を兼ねるので4x4の行列にして求めています。

//射影変換パラメータを取得
var param = this.getParam(origin, markers);

//描画処理用に射影変換パラメータの逆行列を取得
var mx = new Matrix44();
mx.set(param[1], param[2], param[3], 0, param[4], param[5], param[6], 0, param[7], param[8], 1, 0, 0, 0, 0, 1);
mx.Invert();

var inv_param = new Array(9);
inv_param[0] = mx.m11;
inv_param[1] = mx.m12;
inv_param[2] = mx.m13;
inv_param[3] = mx.m21;
inv_param[4] = mx.m22;
inv_param[5] = mx.m23;
inv_param[6] = mx.m31;
inv_param[7] = mx.m32;
inv_param[8] = mx.m33;

//描画処理
draw(inv_param);


代数ライブラリの代わりとして、自作の4x4行列クラスを作成しています。自作するのは数学として理解したいため。
中身としては必要最小限として、単位行列(identity)、掛け算(multi)、コピー(copy)、逆行列(Invert)をもっている。

var Matrix44 = (function () {
    function Matrix44() {
        this.m11 = 1;
        this.m12 = 0;
        this.m13 = 0;
        this.m14 = 0;
        this.m21 = 0;
        this.m22 = 1;
        this.m23 = 0;
        this.m24 = 0;
        this.m31 = 0;
        this.m32 = 0;
        this.m33 = 1;
        this.m34 = 0;
        this.m41 = 0;
        this.m42 = 0;
        this.m43 = 0;
        this.m44 = 1;
        this.identity();
    }

    Matrix44.prototype.identity = function () {
        this.set(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1);
        return this;
    };

    Matrix44.prototype.set = function (m11, m12, m13, m14, m21, m22, m23, m24, m31, m32, m33, m34, m41, m42, m43, m44) {
        this.m11 = m11;
        this.m12 = m12;
        this.m13 = m13;
        this.m14 = m14;
        this.m21 = m21;
        this.m22 = m22;
        this.m23 = m23;
        this.m24 = m24;
        this.m31 = m31;
        this.m32 = m32;
        this.m33 = m33;
        this.m34 = m34;
        this.m41 = m41;
        this.m42 = m42;
        this.m43 = m43;
        this.m44 = m44;
        return this;
    };

    Matrix44.prototype.multi = function (mx) {
        var m11 = this.m11 * mx.m11 + this.m12 * mx.m21 + this.m13 * mx.m31 + this.m14 * mx.m41;
        var m12 = this.m11 * mx.m12 + this.m12 * mx.m22 + this.m13 * mx.m32 + this.m14 * mx.m42;
        var m13 = this.m11 * mx.m13 + this.m12 * mx.m23 + this.m13 * mx.m33 + this.m14 * mx.m43;
        var m14 = this.m11 * mx.m14 + this.m12 * mx.m24 + this.m13 * mx.m34 + this.m14 * mx.m44;
        var m21 = this.m21 * mx.m11 + this.m22 * mx.m21 + this.m23 * mx.m31 + this.m24 * mx.m41;
        var m22 = this.m21 * mx.m12 + this.m22 * mx.m22 + this.m23 * mx.m32 + this.m24 * mx.m42;
        var m23 = this.m21 * mx.m13 + this.m22 * mx.m23 + this.m23 * mx.m33 + this.m24 * mx.m43;
        var m24 = this.m21 * mx.m14 + this.m22 * mx.m24 + this.m23 * mx.m34 + this.m24 * mx.m44;
        var m31 = this.m31 * mx.m11 + this.m32 * mx.m21 + this.m33 * mx.m31 + this.m34 * mx.m41;
        var m32 = this.m31 * mx.m12 + this.m32 * mx.m22 + this.m33 * mx.m32 + this.m34 * mx.m42;
        var m33 = this.m31 * mx.m13 + this.m32 * mx.m23 + this.m33 * mx.m33 + this.m34 * mx.m43;
        var m34 = this.m31 * mx.m14 + this.m32 * mx.m24 + this.m33 * mx.m34 + this.m34 * mx.m44;
        var m41 = this.m41 * mx.m11 + this.m42 * mx.m21 + this.m43 * mx.m31 + this.m44 * mx.m41;
        var m42 = this.m41 * mx.m12 + this.m42 * mx.m22 + this.m43 * mx.m32 + this.m44 * mx.m42;
        var m43 = this.m41 * mx.m13 + this.m42 * mx.m23 + this.m43 * mx.m33 + this.m44 * mx.m43;
        var m44 = this.m41 * mx.m14 + this.m42 * mx.m24 + this.m43 * mx.m34 + this.m44 * mx.m44;
        this.set(m11, m12, m13, m14, m21, m22, m23, m24, m31, m32, m33, m34, m41, m42, m43, m44);
        return this;
    };

    Matrix44.prototype.copy = function (src) {
        this.m11 = src.m11;
        this.m12 = src.m12;
        this.m13 = src.m13;
        this.m14 = src.m14;
        this.m21 = src.m21;
        this.m22 = src.m22;
        this.m23 = src.m23;
        this.m24 = src.m24;
        this.m31 = src.m31;
        this.m32 = src.m32;
        this.m33 = src.m33;
        this.m34 = src.m34;
        this.m41 = src.m41;
        this.m42 = src.m42;
        this.m43 = src.m43;
        this.m44 = src.m44;
        return this;
    };

    Matrix44.prototype.Invert = function () {
        var temp = new Matrix44();
        temp.copy( this);

        var a = temp.m11,  b = temp.m12,  c = temp.m13,  d = temp.m14,
        e = temp.m21,  f = temp.m22,  g = temp.m23,  h = temp.m24,
        i = temp.m31,  j = temp.m32,  k = temp.m33, l = temp.m34,
        m = temp.m41, n = temp.m42, o = temp.m43, p = temp.m44,
        q = a * f - b * e, r = a * g - c * e,
        s = a * h - d * e, t = b * g - c * f,
        u = b * h - d * f, v = c * h - d * g,
        w = i * n - j * m, x = i * o - k * m,
        y = i * p - l * m, z = j * o - k * n,
        A = j * p - l * n, B = k * p - l * o,
        det = 1 / (q * B - r * A + s * z + t * y - u * x + v * w);

        this.m11  = ( f * B - g * A + h * z) * det;
        this.m12  = (-b * B + c * A - d * z) * det;
        this.m13  = ( n * v - o * u + p * t) * det;
        this.m14  = (-j * v + k * u - l * t) * det;
        this.m21  = (-e * B + g * y - h * x) * det;
        this.m22  = ( a * B - c * y + d * x) * det;
        this.m23  = (-m * v + o * s - p * r) * det;
        this.m24  = ( i * v - k * s + l * r) * det;
        this.m31  = ( e * A - f * y + h * w) * det;
        this.m32  = (-a * A + b * y - d * w) * det;
        this.m33 = ( m * u - n * s + p * q) * det;
        this.m34 = (-i * u + j * s - l * q) * det;
        this.m41 = (-e * z + f * x - g * w) * det;
        this.m42 = ( a * z - b * x + c * w) * det;
        this.m43 = (-m * t + n * r - o * q) * det;
        this.m44 = ( i * t - j * r + k * q) * det;

        return this;
    };
})();


次回は、Javascriptの行列ライブラリでは、glMatrix.jsがメジャーなようなので、これに対応してみます。