前回までは、代数ライブラリsylvesterを使用して行列計算を求めてきましたが、今回は代数ライブラリを使わずに8次元連立方程式を4x4行列で求める方法を説明します。
ARToolKitを参考に様々な言語に実装し直したARToolKit互換のクラスライブラリNyARToolKitの開発者@nyatlaさんが射影変換の高速化する上で、8x8行列を4x4行列と2x2行列を組み合わせて解く方法を考案されました。その資料(PDF形式)が下記になります。
遠近法の射影変換パラメータ計算の高速化
数学に疎い私には、8次元連立方程式を4x4行列と2x2行列を組み合わせて解くことが出来るなんて、考えもつきませんでした。
3Dライブラリ内で使われる行列ライブラリをみると4x4行列クラスまでは作成されているので、これが出来るんであれば、射影変換がより身近になると思っています。
「テクスチャマッピングを理解してみる」で書いたのですが、台形に変形するのに２つの三角形を組み合わせて実現しています。もちろんポリゴンオブジェクトを表示するには、三角形を組み合わせるのが通常なのでいいんですが、ツールとして三角形に分割せずに変形させたいなと思っていました。
射影変換パラメータ計算の高速化を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がメジャーなようなので、これに対応してみます。
