QuadWarp用の射影変換(ホモグラフィ)行列の導出

# 前置き
プロジェクションマッピングなどで矩形を歪めて変形するときに、
三角形ポリゴンで行うとuv値の補間が各三角形の中だけで行われるのでテクスチャマッピングが意図した通りになりません

<img src="https://qiita-image-store.s3.amazonaws.com/0/15038/5ab7eade-08a5-4abc-8c1a-597e196b1660.png" width=50%>
これが

<img src="https://qiita-image-store.s3.amazonaws.com/0/15038/31110f63-fbe4-cb61-5d46-bcb9227a6a88.png" width=50%>!
こうなるけど

<img src="https://qiita-image-store.s3.amazonaws.com/0/15038/d4b520e6-d8f0-2713-9c7e-28c797ad3661.png" width=50%>!
ほんとはこうなって欲しい

こんなときに意図通りマッピングするuv値を求めるHomography行列を求めます
答えだけ書いてるところは見たけど導出は見なかったので書いてみます


# 導出
$(0,0)(0,1)(1,1)(1,0)$の正方形を任意形状の四角形$(x_0',y_0')(x_1',y_1')(x_2',y_2')(x_3',y_3')$にする変換行列を求める


```math
\begin{pmatrix}
x_i'' \\
y_i'' \\
z_i'' \\
\end{pmatrix}
=
\begin{pmatrix}
a & b & c \\
d & e & f \\
g & h & 1
\end{pmatrix}
\begin{pmatrix}
x_i \\
y_i \\
1 \\
\end{pmatrix}
\\
x_i'=\frac{x_i''}{z_i''},\quad y_i'=\frac{y_i''}{z_i''}
```


方程式の形に直すと

```math
x_i' = \frac{ax_i + by_i + c}{gx_i + hy_i + 1}\quad
y_i' = \frac{dx_i + ey_i + f}{gx_i + hy_i + 1}
```
$(x_0,y_0)=(0,0), (x_1,y_1)=(0,1), (x_2,y_2)=(1,1), (x_3,y_3)=(1,0)$として、
$(x_0',y_0')～(x_3',y_3')$を入力したときそれを満たす$a～h$が求まればよい
言い換えると$a～h$を$(x_0',y_0')～(x_3',y_3')$で表せられればよい
未知変数$a～h$は８つ、立てられる連立方程式は$(x_0',y_0')～(x_3',y_3')$で８つ、求まりそう

$iが0～3$のときそれぞれ考えると、

$i=0$のとき
$(x_0,y_0)=(0,0)$なので、

```math
x_0' = \frac{a*0 + b*0 + c}{g*0 + h*0 + 1}\quad
y_0' = \frac{d*0 + e * 0 + f}{g*0 + h*0 + 1}\\
```

```math
c = x_0' \quad
f = y_0' \qquad\cdots①
```


$i=1$のとき
$(x_1,y_1)=(0,1)$なので

```math
x_1' = \frac{a*0 + b*1 + c}{g*0 + h*1 + 1}\quad
y_1' = \frac{d*0 + e*1 + f}{g*0 + h*1 + 1}\\
```

```math
x_1' = \frac{b + c}{h + 1}\quad
y_1' = \frac{e + f}{h + 1}\qquad\cdots②
```

以降同様に、
$i=2$のとき

```math
x_2' = \frac{a + b + c}{g + h + 1}\quad
y_2' = \frac{d + e + f}{g + h + 1}\qquad\cdots③
```
$i=3$のとき

```math
x_3' = \frac{a + c}{g + 1}\quad
y_3' = \frac{d + f}{g + 1}\qquad\cdots④
```

①～④で連立方程式８つが出た（いきなり$c,f$もわかった！）ので、あとはゴリゴリ解く

```math
①②より、
b=x_1'h-x_0'+x_1'\quad
e=y_1'h-y_0'+y_1'\qquad\cdots⑤\\
③より
a+b-x_2'g-x_2'h+x_0'-x_2'=0\quad
d+e-y_2'g-y_2'h+y_0'-y_2'=0\qquad\cdots⑥\\
①④より
a=x_3'g-x_0'+x_3'\quad
d=y_3'g-y_0'+y_3\qquad\cdots⑦\\
```
⑤、⑦を⑥に代入して、$g,h$の式にすると

```math
(x_3'g-x_0'+x_3')+(x_1'h-x_0'+x_1')-x_2'g-x_2'h+x_0'-x_2'=0\\
(y_3'g-y_0'+y_3)+(y_1'h-y_0'+y_1')-y_2'g-y_2'h+y_0'-y_2'=0\\
```

整理して

```math
(x_1'-x_2')h+(x_3'-x_2')g-x_0'+x_1'-x_2'+x_3'=0\qquad(y_1'-y_2')h+(y_3'-y_2')g-y_0'+y_1'-y_2'+y_3'=0\qquad\cdots⑧
```

ここでおもむろに以下のように定義して

```math
sx = x_0'-x_1'+x_2'-x_3'\\
sy = y_0'-y_1'+y_2'-y_3'\\
dx1 = x_1'-x_2'\\
dx2 = x_3'-x_2'\\
dy1 = y_1'-y_2'\\
dy2 = y_3'-y_2'
```

⑧に代入すると

```math
dx1*h+dx2*g-sx=0\qquad dy1*h+dy2*g-sy=0
```

$g,h$について連立方程式を解いて、

```math
g=\frac{sx*dy1-sy*dx1}{dy1*dx2-dx1*dy2}\quad
h=\frac{sy*dx2-sx*dy2}{dy1*dx2-dx1*dy2}\qquad\cdots⑨
```

$g,h$が求まったので、⑤,⑦に代入すれば、$a,b,d,e$が求まり、$a～h$のすべてが求まる


めでたしめでたし

# 実装
上記の変換行列は$(0,0)(0,1)(1,1)(1,0)$を任意の位置に移動させる行列ですが、
uv値はこの逆で、任意の位置を$(0,0)(0,1)(1,1)(1,0)$へ変換すればいいので逆行列をポジションに乗算する実装になります


#ソースコード（Unity）
https://github.com/fuqunaga/QuadWarp


# 参考
http://yaju3d.hatenablog.jp/entry/2013/08/04/152524

連立方程式を行列の式の形にして逆行列で求めるやり方
http://mf-atelier.sakura.ne.jp/mf-atelier/modules/tips/index.php/program/algorithm/a6.html

OpenCVは４点に限らない高度なのやってるっぽい
http://opencv.jp/opencv-2svn/cpp/camera_calibration_and_3d_reconstruction.html#cv-findhomography