結婚式で余興を頼まれた時の画像処理をPythonでやってみる

#モチベーション

画像処理は色々なところで使える技術です。
そうたとえ結婚式の余興であってもです。
20代後半になると周りが結婚しているので、いやでも意識してしまう（汗

<img width="400" alt="couple_afin.jpg" src="https://qiita-image-store.s3.amazonaws.com/0/10496/5d18ba25-321c-e375-90cc-0097bb5382b5.jpeg">


[結婚式・ウェディングのテーブルセット](
 https://www.pakutaso.com/20130406108post-2650.html"
結婚式・ウェディングのテーブルセット") 

そして結婚式といえば余興

時間がない中でもクオリティを上げたい・・・
そんなニーズに応えてくれます。
そう画像処理ならね。

本記事で扱う範囲としては、「画像間の写像」を対象とします([実践 コンピュータビジョンの3章](https://www.oreilly.co.jp/books/9784873116075/)に相当)。なお、上級画像処理者の方は誤りなどあれば指摘していただければ幸いです。

#画像間の写像に必要な処理

余興用の写真を用意するときに必要になるのが写真の合成です。

写真の合成も単純に行うと違和感のある写真になってしまいます。

そこで必要になってくる処理は下記になります。

・画像の変形
・画像中に上手く画像を埋め込む
・位置合わせ
・パノラマ用の位置合わせ

#各処理に必要な画像処理技術

画像の変換処理に必要な要素は大きく分けると下記になります。

<img width="850" alt="Screen Shot 2016-02-04 at 18.14.31.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/e8ca9783-f535-24a1-e5f6-513a33a5a043.png">

[EECS 442 – Computer vision Conversation hour](
http://wwweb.eecs.umich.edu/vision/teaching/EECS442_2011/lectures/discussion2.pdf "EECS 442 – Computer vision Conversation hour") 

・画像のスケーリング
・画像のスケーリング（非対称）
・回転
・言葉では説明できないので一番右の図をご覧下さい。

これらの処理を組み合わせて画像を変換していきます。

ここで重要になってくる数式があります。

```math

\begin{pmatrix}
x'  \\
y'  \\
w'
\end{pmatrix}

=

\begin{pmatrix}
h_{11} & h_{21} & h_{31}\\
h_{12} & h_{22} & h_{32}\\
h_{13} & h_{23} & h_{33}
\end{pmatrix}


\begin{pmatrix}
x  \\
y  \\
w
\end{pmatrix}
```

```math
\begin{matrix}
\boldsymbol{x'} = \boldsymbol{H}\cdot{\boldsymbol{x}}
\end{matrix}
```


上記の式はある点がベクトルで与えられたときにその点を変換する行列をかけて別の空間に写像していることを表しています。
言葉でいうと意味不明なので、下記に図を載せます。

<img width="799" alt="Screen Shot 2016-02-04 at 18.45.46.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/b7e0d449-00e9-f012-fd6d-64b3a26b0b58.png">

[Lecture 16: Planar Homographies](
http://www.cse.psu.edu/~rtc12/CSE486/lecture16.pdf "Lecture 16:
Planar Homographies") 

Pinholeカメラから物体を覗いたときにPinhole camera上でどう見えるかを表したり、Pinhole cameraから平面のイメージで表すときに行列Hをかけて変換していることを表しています。

最初に述べた基本的な画像変換の処理を行列Hに網羅しています。
網羅された行列Hを求めれば、Pinholeカメラから物体を覗いたときにPinhole camera上でどう見えるかを表したり、Pinhole cameraから平面上に画像を表すことができます。
この行列Hをホモグラフィー行列と呼びます。

ここで一旦まとめます。

合成写真を作成するには画像変換処理が必要。
画像変換処理はホモグラフィー行列を使用すれば行える。
ホモグラフィー行列を使用すると平面から斜めに移った画像または斜めから平面の画像などを再現することが可能

#ホモグラフィー行列

ここでホモグラフィー行列を求める方法について述べます。

Direct Linear Transfer(直接線形変換)アルゴリズムを使用します。

では先ほどの数式をもう一度確認します。


```math

\begin{pmatrix}
x'  \\
y'  \\
w'
\end{pmatrix}

=

\begin{pmatrix}
h_{11} & h_{21} & h_{31}\\
h_{12} & h_{22} & h_{32}\\
h_{13} & h_{23} & h_{33}
\end{pmatrix}


\begin{pmatrix}
x  \\
y  \\
w
\end{pmatrix}
```

今回は画像なので２次元の座標を表していますが、同時座標系と呼ばれる座標を採用しており、x,yの値に加えてwの値があります。

wの値は通常は正規化して１にしますが、この値があることの効果は平行移動も表せることです。

画像の変換において、上記４つだけなく位置が変わる平行移動が表現できる。
これは重要なポイントです。

分かりやすい図解は下記をご覧ください。

http://d.hatena.ne.jp/Zellij/20120523/p1

次は正規化したパターンを確認してみましょう。
この数式から求めたいのは平面の画像を斜めの画像にしたり、斜めの画像を平面にしたりしてくれるホモグラフィー行列です。

```math

\begin{pmatrix}
x'  \\
y'  \\
1'
\end{pmatrix}

=

\begin{pmatrix}
h_{11} & h_{21} & h_{31}\\
h_{12} & h_{22} & h_{32}\\
h_{13} & h_{23} & h_{33}
\end{pmatrix}


\begin{pmatrix}
x  \\
y  \\
1
\end{pmatrix}
```

元の画像と写像した画像の点における関係性を数式で表すと

```math
  x' = \frac{h_{11}x + h_{21}y + h_{31}}{h_{13}x + h_{23}y + h_{33}}\\

  y' = \frac{h_{12}x + h_{22}y + h_{32}}{h_{13}x + h_{23}y + h_{33}}

```

行列Hを求めるには9次元の値を求める必要がありますがここで制約を加えます。
h33の値を1である制約をかけます。
なぜここで制約をかけても大丈夫なのかですが、

```math
  x' = \frac{kh_{11}x + kh_{21}y + kh_{31}}{kh_{13}x + kh_{23}y + kh_{33}}\\

  y' = \frac{kh_{12}x + kh_{22}y + kh_{32}}{kh_{13}x + kh_{23}y + kh_{33}}\\
\\
```

h33を1にした場合  

```math
  x' = \frac{h_{11}x + h_{21}y + h_{31}}{h_{13}x + h_{23}y + 1}\\

  y' = \frac{h_{12}x + h_{22}y + h_{32}}{h_{13}x + h_{23}y + 1}

```

k倍した値を求めても分母、分子で割るので実質的に無意味なので、さらにh33はx,yの両者で使用されています。ここでkはh33を1にするような値として想定しても問題ない。
そのため、h33の値を1である制約をかけて９次元の値を求める問題から8次元の値を求める問題にできました。

ここから直接的線形変換法の直接的な解法

```math

  (h_{13}x + h_{23}y + 1)x' = h_{11}x + h_{21}y + h_{31}\\

  ({h_{13}x + h_{23}y + 1})y' = h_{12}x + h_{22}y + h_{32}
```

求められたx',y'

```math

 　x' = h_{11}x + h_{21}y + h_{31} - h_{13}xx' - h_{23}yx' 　\\

  y' = h_{12}x + h_{22}y + h_{32} - h_{13}xy' - h_{23}yy' 

```

写像された１点のデータから２つの数式が得られることがわかりました。
つまり数式的には４つの点が求められれば８つの値を求めるのに必要な連立方程式が得られることがわかりました。
ここで４点が得られたときの式を表してみます。

```math

\begin{pmatrix}
x'_{1}  \\
y'_{1}  \\
x'_{2}  \\
y'_{2}  \\
x'_{3}  \\
y'_{3}  \\
x'_{4}  \\
y'_{4}  
\end{pmatrix}

=

\begin{pmatrix}
x_{1} & y_{1} & 1 & 0 & 0 & 0 & -x_{1}x'_{1} & -y_{1}x'_{1}\\
0 & 0 & 0 & x_{1} & y_{1} & 1 & -x_{1}y'_{1} & -y_{1}y'_{1}\\
x_{2} & y_{2} & 1 & 0 & 0 & 0 & -x_{2}x'_{2} & -y_{2}x'_{2}\\
0 & 0 & 0 & x_{2} & y_{2} & 1 & -x_{2}y'_{2} & -y_{2}y'_{2}\\
x_{3} & y_{3} & 1 & 0 & 0 & 0 & -x_{3}x'_{3} & -y_{3}x'_{3}\\
0 & 0 & 0 & x_{3} & y_{3} & 1 & -x_{3}y'_{3} & -y_{3}y'_{3}\\
x_{4} & y_{4} & 1 & 0 & 0 & 0 & -x_{4}x'_{4} & -y_{4}x'_{4}\\
0 & 0 & 0 & x_{4} & y_{4} & 1 & -x_{4}y'_{4} & -y_{4}y'_{4}\\
\end{pmatrix}


\begin{pmatrix}
h_{11}  \\
h_{12}  \\
h_{13}  \\
h_{21}  \\
h_{22}  \\
h_{23}  \\
h_{31}  \\
h_{32}
 \end{pmatrix}
```

上記の連立方程式を解析的に求めることも可能です。
通常は４点以上の点を取って再現する方がより再現度が上がることが自明です。

点が多くなってきたときの解法ですが、SVDを用いた最小２乗法による解決です。
SVDにより得られた解の行列Vの最後の行が今回の最小２乗法による解になります。
（点群は平均０、分散１で正規化された空間から観測した値を使用します。）

なぜそうなるか確認したい方は下記をご覧下さい。

[特異値分解](http://www.misojiro.t.u-tokyo.ac.jp/~murota/lect-kisosuri/singularvalue050129.pdf "特異値分解") 

ここでまとめておくと

・ホモグラフィー行列は画像を平面から斜め、また斜めから平面へ移してくれる行列
・行列の計算は４点を取れば計算可能
・４点以上取って計算する場合はSVDを使用する。

#アフィン変換

ここからはホモグラフィ行列を簡単化することで解ける問題を限定する代わりに必要な点を３点にし、３点にすることによる利点を活かせるアフィン変換について紹介していきます。

```math

\begin{pmatrix}
x'  \\
y'  \\
1
\end{pmatrix}

=

\begin{pmatrix}
a_{1} & a_{2} & t_{1}\\
a_{3} & a_{4} & t_{2}\\
0 & 0 & 1
\end{pmatrix}


\begin{pmatrix}
x  \\
y  \\
1
\end{pmatrix}
```

上記の式が意味していることはスケーリング１にしています。
つまり変換量を抑えているということです。

四角形から平行四辺形の変換や拡大、縮小、平行移動まで行えますが四角形から台形への変換は行えません。

```math
  x' = \frac{h_{11}x + h_{21}y + h_{31}}{h_{13}x + h_{23}y + 1}\\

  y' = \frac{h_{12}x + h_{22}y + h_{32}}{h_{13}x + h_{23}y + 1}

```

上記の式を確認するとx,yの値に応じてx',y'の値を共通に変換できています。つまり座標の位置によって全体の値の大きさを変えることができていたのですが、下記のように変わったので、単純な変換と平行移動しか扱えなくなりました。その代わり必要な点が３点で済みます。この利点は分割アフィンワーピングが使用できる点です。

```math
  x' = h_{11}x + h_{21}y + h_{31}\\

  y' = h_{12}x + h_{22}y + h_{32}

```

##画像を別の画像に埋め込む

画像を別の画像に埋め込むのに分割アフィンワーピングを使用します。

下記のような二つの画像を合わせたいときに通常であれば点群を取って、ホモグラフィ行列では点をサンプリングして計算していきますが、計算誤差などにより画像の頂点に合わせることは難しくなります。

しかしドロネー三角形分割法を用いたアフィン変換だと、画像の頂点を合わせることが可能になります。

<img width="750" alt="Screen Shot 2016-02-04 at 21.39.02.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/cbbeb343-b0d9-20c2-1886-7cab4c6331c4.png">

ではなぜドロネー三角形分割法を用いたアフィン変換を用いれば画像の頂点を合わせることが可能なのか説明していきます。

アフィン変換で重要な点は３点のみで良い点です。

そしてドロネー三角形分割法はある点群が得られたときにその点群を最小角度が最大となる三角形で結んでいく手法です。
三角形の最小角度が最大となる選び方については下記をご覧下さい。

[Delaunay 三角形分割](
http://jp.mathworks.com/help/matlab/math/delaunay-triangulation.html?requestedDomain=jp.mathworks.com&requestedDomain=jp.mathworks.com "Delaunay 三角形分割") 

例を下記に示します。

<img width="674" alt="Screen Shot 2016-02-04 at 21.44.21.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/f1fc1778-3dd0-111e-8060-13d28a3b8087.png">

[Programming Computer Vision with Python](
http://programmingcomputervision.com/downloads/ProgrammingComputerVision_CCdraft.pdf "Programming Computer Vision
with Python") 

つまり全ての点に対してアフィン変換がかかるので頂点を合わせた変換が可能になります。
計算方法はDLT法でアフィン変換も同様に計算します。

まとめると

・アフィン変換はホモグラフィ変換よりも変換できるバリエーションは減る
・その代わりに３点のみで変換可能
・３点で良いとドロネー三角形分割法が使用できるため、頂点を合わせた変換が可能

#画像の位置合わせ

画像の位置合わせの場合は画像自体は同様の画像を使用するため相似変換が使用されます。

```math

\begin{pmatrix}
x'  \\
y'  \\
1
\end{pmatrix}

=

\begin{pmatrix}
s\ cos(\theta) & -s\ sin(\theta) & t_{x}\\
s\ sin(\theta) & s\ cos(\theta) & t_{y}\\
0 & 0 & 1
\end{pmatrix}


\begin{pmatrix}
x  \\
y  \\
1
\end{pmatrix}
```

sは拡大率、θは回転率を表しています。
座標全体に対して拡大、縮小、回転をかけられるので、位置合わせには最適な変換となっています。
下記の画像の平均を単純に取ったものと位置合わせして平均を取ったものを比較してみます。
位置合わせをした方が精度よく画像を再現できていることが確認できます。

<img width="690" alt="Screen Shot 2016-02-04 at 22.15.25.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/925d7e42-33d0-f603-dc89-5aef9f14bedc.png">

[Programming Computer Vision with Python](
http://programmingcomputervision.com/downloads/ProgrammingComputerVision_CCdraft.pdf "Programming Computer Vision
with Python") 


<img width="497" alt="Screen Shot 2016-02-04 at 22.13.58.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/ed35d903-06e6-377d-950b-1308491eb339.png">

[Programming Computer Vision with Python](
http://programmingcomputervision.com/downloads/ProgrammingComputerVision_CCdraft.pdf "Programming Computer Vision
with Python") 


#RANSAC

DLT法しか紹介していないですが、DLT法はノイズには頑健ではないため、別の手法が提案されています。
この手法はノイズに対して頑健な代わりに時間がかかるアルゴリズムです。

具体的な手法については下記の資料を参考にしてください。

[Preemptive RANSAC by David Nister.](
http://www.slideshare.net/enddl22/preemptive-ransac-by-david-nister)

#画像のつなぎ合わせ

同じ場所から撮影された画像がある。
これを２つ目の図のように画像をつなぎ合わせたい場合

<img width="647" alt="Screen Shot 2016-02-04 at 22.36.12.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/f8dd9816-5c73-f6be-6a4a-e9d720b3b73d.png">

[Programming Computer Vision with Python](
http://programmingcomputervision.com/downloads/ProgrammingComputerVision_CCdraft.pdf "Programming Computer Vision
with Python") 


<img width="727" alt="Screen Shot 2016-02-04 at 22.34.33.png" src="https://qiita-image-store.s3.amazonaws.com/0/10496/27e8cc02-a96f-cd5f-39ac-5b4dbedfb7e9.png">

[Programming Computer Vision with Python](
http://programmingcomputervision.com/downloads/ProgrammingComputerVision_CCdraft.pdf "Programming Computer Vision
with Python") 

１：画像の特徴点をSHIFTによって取得（詳しくは下記を参照）
[画像処理の数式を見て石になった時のための、金の針](
http://qiita.com/icoxfog417/items/adbbf445d357c924b8fc "画像処理の数式を見て石になった時のための、金の針") 
２：画像間のホモグラフィー行列をRANSACを用いて計算
３：中央の画像を選択
４：中央の画像の左右に０を追加して、変形画像をつなぎ合わせる。
５：平行成分によって左右のどちらにつなげるかを判断している。

画像間の差が目立つので市販のソフトでは輝度を正規化して滑らかにする処理を入れています。

#ipython notebook

アフィン変換のみですがコードで記述しています。

https://github.com/SnowMasaya/Image_Processing_for_Wedding

#参考

[実践 コンピュータビジョン](https://www.oreilly.co.jp/books/9784873116075/)
本記事の素案

[ペンシルベニア州立大学 コンピューターサイエンス](http://www.cse.psu.edu/~rtc12/CSE486/lecture16.pdf)
３次元からの画像変換についてのレクチャーが記述しており、今回の記事のように２次元に限定されていないため分かりやすい

[ミシガン大学](http://wwweb.eecs.umich.edu/vision/teaching/EECS442_2011/lectures/discussion2.pdf)

[ミシガン大学のまとめ](http://wwweb.eecs.umich.edu/vision/teaching/EECS442_2012/class_schedule.html)

画像変換の基本的なところからの資料なので基本的なところを抑えられる。

[Preemptive RANSAC by David Nister.](http://www.slideshare.net/enddl22/preemptive-ransac-by-david-nister?qid=8a42dc4f-37d1-4ee5-82e5-a0d8bfeb0ba6&v=default&b=&from_search=1)

RANSACのアルゴリズムがすぐにわかる優れたスライド

[大人になってからの再学習：アフィン変換](http://d.hatena.ne.jp/Zellij/20120523/p1)

アフィン変換の分かりやすい説明