2019-03-19

機械学習で桜の開花日を予想してみた

機械学習で予想した東京の開花日は３月２５日でした。３月１０日までの気象情報から推測しています。
今年は比較的暖冬で例年（平均：３月２６日）より早めでしょうと言われているので、それを反映できていると思います。ちなみに、weathermap[1]の予想日は３月２３日でした。（３月１４日現在）

方法

気象庁公開の気象データ[2]、さくらの開花日データ[3]を使います。
気象データは任意の観測地、観測期間、観測値を指定してcsv形式で入手することが出来ます。今回は、以下の条件の観測データを学習に使用しました。

観測地点：東京、銚子、甲府、熊谷、水戸、長野、名古屋、静岡、津、宇都宮、横浜の11地点
期間：2000年～2018年までの18年分
観測値：平均気温、最高気温、降水量の合計、日照時間、最深積雪、降雪量の合計、平均風速、最大風速、最大瞬間風速、平均蒸気圧、平均湿度、最小相対湿度、平均現地気圧、平均海面気圧、最低海面気圧、平均雲量、10分間降水量の最大の17つ
観測間隔：毎日

さくら開花前年の４月１日～３月１０日までの345日分のデータを使って学習します。

１データあたりの特徴量数は、345（日数）ｘ17（観測値の数）＝5865となります。
全データ数は、11（観測地点数）ｘ18（観測年数）＝198です。

学習にはランダムフォレストを使用しました。これはランダムに特徴量を選択しデータの予測を行う決定木を作成する方法です。単純な決定木よりもノイズに強いという特徴があります。また、学習の副産物として、どの特徴量が決定木の構成によく使われたのか知ることが出来るので、学習後にデータの分析をするのに役立ちます。

データの前処理

観測値が空白の日について、１日後の値で補完しました。
さくらの開花日は、１月１日を０として0-364で数値化しました。
２月２９日はデータから除外しました。
各特徴量は学習データの最大値、最小値で0-1に正規化しました。

学習結果のバリデーション

方法の良し悪しをクロスバリデーションで評価します。
クロスバリデーションでは、学習に使用するデータを分割して学習・予測・評価することを繰り返し、平均的な性能を見ることが出来ます。評価はMAEで行います。

${ MAE = \Sigma {|実際の値 - 予測値|} \tag{1} }$

これは予測誤差日数の平均ということになります。
比較を行うために、予測結果を過去の平均値とした時のMAEを計算しました。

	MAE
過去の平均値	5.36
予測値	2.67

さくらの開花日を過去の平均値として計算する場合5.36日の誤差があるが、予測をすると誤差は2.67日に減少することが分かりました。

ランダムフォレストでは、どの特徴量が予測に有効であるか知ることが出来ます。5865の特徴量のうち上位10個を挙げます。

特徴量名　重要度　重要度の累積値

Ave(tmp)[cels]_day58_2/28	0.21 	0.21 
Ave(tmp)[cels]_day299_10/27	0.04 	0.25 
Ave(tmp)[cels]_day11_1/12	0.04 	0.29 
Ave(tmp)[cels]_day364_12/31	0.02 	0.32 
Max(tmp)[cels]_day10_1/11	0.02 	0.34 
Ave(tmp)[cels]_day60_3/2	0.02 	0.36 
Max(tmp)[cels]_day59_3/1	0.02 	0.38 
Max(tmp)[cels]_day363_12/30	0.02 	0.40 
Ave(humid)[％]_day39_2/9	0.02 	0.41 
Ave(tmp)[cels]_day59_3/1	0.01 	0.42

Ave(tmp)[cels]..平均気温
Max(tmp)[cels]..最高気温
Ave(humid)[％]..平均湿度

２月２８日の平均気温が最も有効な特徴量でした。１０分の９が気温に関する特徴量であり、さくらの開花予測には気温が有効であることが分かりました。また、１０月以降の気温が開花日に大きく影響するが分かりました。

結果

学習に使用した観測点について今年の開花日を予想しました。
月日は四捨五入で計算しています。東京は82.58->83として３月２５日になりました。

地点	開花日（0-364）	開花日（月日）
銚子	83.91	3/26
甲府	81.41	3/23
熊谷	82.42	3/24
水戸	83.9	3/26
長野	93.83	4/05
名古屋	80.97	3/23
静岡	80.44	3/22
東京	82.58	3/25
津	81.15	3/23
宇都宮	85.32	3/27
横浜	82.73	3/25

コード

pyCodes/mlPrediction at master · misakikobayashi1984/pyCodes · GitHub

参考

[1]さくら開花予想2019, https://sakura.weathermap.jp/
[2]過去の気象データ, https://www.data.jma.go.jp/obd/stats/etrn/
[3]さくらの開花日, https://www.data.jma.go.jp/sakura/data/sakura003_06.html

2019-03-19

最小二乗フィッティングで糖尿病の進行度予測

糖尿病のサンプルデータを使って最小二乗フィッティングしました。糖尿病のサンプルデータを使って最小二乗フィッティングしました。

理論

観測値 $x,y$ は $y=f(x)$ というモデルで表されるとします。モデルが $x_{i} (i=1,..n)$ の線形結合で表される時、 $f(x) = \Sigma a_{i} x_{i}$ と書けます。 $a_{i}$ は $x_{i}$ の係数です。
次のような $k$ 個の観測値の集合がある場合を考えます。

${ (x_{11}, x_{12}, \dots ,x_{1n}, y_{1}), (x_{21}, x_{22}, \dots ,x_{2n}, y_{2}), \dots, (x_{k1}, x_{k2}, \dots ,x_{kn}, y_{k}) \tag{1} }$

この集合からモデルに最も当てはまる係数 $a_{i}$ を求めるには、モデルから求めた理論値 $f(x)$ と観測値 $y$ を最小にします。理論値と観測値の二乗誤差の総和は、

${ J = \Sigma {(f(x)- y_{i})}^{2} = \Sigma {( \Sigma a_{i} x_{i} - y_{i})}^{2} \tag{2} }$

と書けます。Jを最小にすればよいということが分かります。
(2)の $\Sigma a_{i} x_{i} - y_{i}$ を行列で表すと以下のようになります。

${ \begin{pmatrix} x_{11} && x_{12} && \dots && x_{1n} \\ x_{21} && x_{22} && \dots && x_{2n} \\ && && \dots && \\ x_{k1} && x_{k2} && \dots && x_{kn} \\ \end{pmatrix} \begin{pmatrix} a_{1} \\ a_{2} \\ \vdots \\ a_{n} \\ \end{pmatrix} - \begin{pmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{k} \\ \end{pmatrix} \tag{3} }$

(2)の $J$ を最小にするということは、(3)の大きさを最小にすればいいですね。簡略化＆理想的な場合を考えると、(3) $= 0$ と置いてこれを解けばよいということが分かります。

実装

scikit-learnのLinearRegression()を使います。[1]の例では、１次元の特徴量に限定してフィッティングしていますが、全特徴量を使ってフィッティングしました。Variance scoreはフィッテイングの良さを示す値で、データのばらつきをもとに計算します。[1]では0.47でしたが、全特徴量を使った場合0.59となり、よりよくフィッティング出来ています。

出力

Coefficients:  [ 2.06069866e-04 -7.05675908e-02  4.14834780e-01  2.49574079e-01 -7.11907159e-01  4.82683572e-01  9.08289960e-02  1.50462309e-01  6.01524431e-01  6.48063127e-02]
Mean squared error: 0.02
Variance score: 0.59

学習データ、テストデータをそれぞれプロットし求めたフィッティングを描画してみました。各特徴量に対する観測値の値です。プロット点はデータ、線は予測結果を表しています。各プロットのタイトルには求めた係数を表示しています。係数の大きさが大きいほど重要な特徴量と考えられるのです。ここで最も大きいものは $s1$ という特徴量でしたが、プロット点を見ると観測値と相関がないように見えます。また、データの傾向に反して負の係数であるのが気になります。未解決です。

f:id:mikekochang:20190319135801j:plain
Fig1. 学習データのプロット、フィッテイング結果

f:id:mikekochang:20190319135611p:plain
Fig2. テストデータのプロット、フィッテイング結果

コード

GitHub - misakikobayashi1984/pyCodes

参考

[1]Linear Regression Example, https://scikit-learn.org/stable/auto_examples/linear_model/plot_ols.html#sphx-glr-auto-examples-linear-model-plot-ols-py

2019-03-07

カメラの基本行列を理解する、エピポーラ線を描く

２つの画像における $\boldsymbol{E}$ の投影点を $\boldsymbol{x}_1$ , $\boldsymbol{x}_2$ とし、それを対応付ける行列 $\boldsymbol{E}$ を求めます。 $\boldsymbol{E}$ を基本行列と呼び、以下が成り立ちます[1]。

${ \boldsymbol{x}_1^T \boldsymbol{E} \boldsymbol{x}_2 = 0 \tag{1} }$

$\boldsymbol{E}$ の表すもの

$X$ , $\boldsymbol{x}_1$ , $\boldsymbol{x}_2$ の位置関係を図示するとFig1のようになっています。
f:id:mikekochang:20190307115528j:plain
Fig1. $\boldsymbol{X}$ , $\boldsymbol{x}_1$ , $\boldsymbol{x}_2$ の位置関係

$\boldsymbol{x}_1$ は $\boldsymbol{x}_2$ を回転、平行移動することで表すことができるので、回転を $\boldsymbol{R}$ 、平行移動を $\boldsymbol{t}$ として

${ \boldsymbol{x}_1 = \boldsymbol{R} \boldsymbol{x}_2 + \boldsymbol{t} \tag{2} }$

と書けます。 $\boldsymbol{R}$ $\boldsymbol{x}_2$ , $\boldsymbol{t}$ は同一平面上にあるので、この外積 $\boldsymbol{t} \times \boldsymbol{R} \boldsymbol{x}_2$ は平面上の垂線ベクトルを表します。
また、 $\boldsymbol{x}_1$ も同様に同一平面上にあるので、垂線ベクトルと直交する（＝垂線ベクトルとの内積が0となる）ことから

${ \boldsymbol{x}_1 (\boldsymbol{t} \times \boldsymbol{R} \boldsymbol{x}_2)= 0 \tag{3} }$

となります。ここで、外積を行列表現します。
$\boldsymbol{a}=(a_{1}, a_{2}, a_{3})$ の時

${ \boldsymbol{a} \times \boldsymbol{b} = \begin{pmatrix} 0 && -a_{3} && a_{2} \\ a_{3} && 0 && -a_{1} \\ -a_{2} && a_{1} && 0 \\ \end{pmatrix} \boldsymbol{b} \tag{4} }$
と表せるので、(3)は

${ \boldsymbol{x}_1^T \boldsymbol{E} \boldsymbol{x}_2 = 0 \tag{1} }$
ここで、
${ \boldsymbol{E} = \boldsymbol{TR} \tag{5} }$

${ \boldsymbol{T} = \begin{pmatrix} 0 && -t_{3} && t_{2} \\ t_{3} && 0 && -t_{1} \\ -t_{2} && t_{1} && 0 \\ \end{pmatrix} \tag{6} }$

と表せます。 $\boldsymbol{E}$ は $\boldsymbol{TR}$ のことですね。
また、(1)の両辺の転置を取ると、
${ \boldsymbol{x}_2^T \boldsymbol{E}^T \boldsymbol{x}_1 = 0 \tag{7} }$
となるので、 $\boldsymbol{x}_1$ を $\boldsymbol{x}_2$ に対応付ける行列は $\boldsymbol{E}^T$ であると言えます。

$\boldsymbol{E}$ の性質

$\boldsymbol{E}$ はカメラの内部パラメータで正規化済みの点 $\boldsymbol{x}_1$ , $\boldsymbol{x}_2$ の対応を求めたものです。正規化をしていない点 $\boldsymbol{x}_1^{'}$ , $\boldsymbol{x}_2^{'}$ に関しては、 $\boldsymbol{K}$ で正規化を行うと $\boldsymbol{x}_1 = \boldsymbol{K}^{-1} \boldsymbol{x}_1^{'}$ , $\boldsymbol{x}_2 = \boldsymbol{K}^{-1} \boldsymbol{x}_2^{'}$ と表せるので、

${ \boldsymbol{x}_1^{T} E \boldsymbol{x}_2 = (\boldsymbol{K}^{-1} \boldsymbol{x}_1^{'} )^{T} E (\boldsymbol{K}^{-1} \boldsymbol{x}_2^{'} ) = {\boldsymbol{x}_1^{'}}^{T} {\boldsymbol{K}^{-1}}^{T} E \boldsymbol{K}^{-1} \boldsymbol{x}_2 \tag{8} }$

と表せ、 $F = {\boldsymbol{K}^{-1}}^{T} E \boldsymbol{K}^{-1}$ と置き、 $F$ を基礎行列と呼びます。

$\boldsymbol{E}$ のランクについて、 $rank(\boldsymbol{T})=2$ （ $\boldsymbol{T}$ を変形してランクを調べる操作をすると分かる）、 $rank(\boldsymbol{R})=3$ から、 $rank(\boldsymbol{E})=rank(\boldsymbol{TR})=rank(\boldsymbol{T})=2$ である[2]と言えます。（ $\boldsymbol{F}$ についても同様です。）

エピ極、エピポーラ線の求め方

画像1の光学中心 $\boldsymbol{o}_1$ を画像2の光学中心に投影した点がエピ極です(Fig1.の $\boldsymbol{e}_2$ ）。また、 $\boldsymbol{e}_2$ と $\boldsymbol{x}_2$ を結んだ線をエピポーラ線と呼びます。すべてのエピポーラ線はエピ極を通過します。これは、

${ \boldsymbol{x}_1^{T} \boldsymbol{F} \boldsymbol{e}_2 \forall \boldsymbol{x}_1 \tag{9} }$

と表すことが出来ます。（どんな $\boldsymbol{x}_1$ に対しても $\boldsymbol{F}$ によって $\boldsymbol{e}_2$ と対応づけることが出来る。）ここから、

${ \boldsymbol{F} \boldsymbol{e}_2 = \boldsymbol{0} \tag{10} }$

ということが分かり、これを解くことで $\boldsymbol{e}_2$ が求まります。（特異値分解で解ける）

$\boldsymbol{x}_1$ を通るエピポーラ線は、
${ \boldsymbol{l}= \boldsymbol{F} \boldsymbol{x}_2 \tag{11} }$

として求めることが出来ます。これは $\boldsymbol{x}_1^{T} \boldsymbol{l}= \boldsymbol{0}$ という直線を表し、 $\boldsymbol{l}$ は直線の係数ということになります。

基礎行列の求め方

８点アルゴリズムというシンプルな方法で解くことを考えます。(1)を要素で表し、

${ (x_{1}, y_{1}, 1) \begin{pmatrix} f_{11} && f_{12} && f_{13} \\ f_{21} && f_{22} && f_{23} \\ f_{31} && f_{32} && f_{33} \\ \end{pmatrix} \begin{pmatrix} x_{2} \\ y_{2} \\ 1 \\ \end{pmatrix} = \boldsymbol{0} \tag{9} }$

とします。これを書き換えると、

${ \boldsymbol{b}_i \boldsymbol{f} = \boldsymbol{0} \tag{12} }$

ここで、
${ \boldsymbol{b}_i = (x_{2i}x_{1i},y_{2i}x_{1i}, x_{1i}, x_{2i}y_{1i}, y_{2i}y_{1i}, y_{1i}, x_{2i}z_{1i}, y_{2i}z_{1i}, 1) \tag{13} }$

${ \boldsymbol{f} = (f_{11}, f_{12}, f_{13}, f_{21}, f_{22}, f_{23}, f_{31}, f_{32}, f_{33})^{T} \tag{14} }$

と表せます。対応点が８点ある場合、

${ \boldsymbol{Bf} = 0 \tag{15} }$

${ \boldsymbol{B} = \begin{pmatrix} \boldsymbol{b}_{1} \\ \vdots \\ \boldsymbol{b}_{8} \\ \end{pmatrix} \tag{16} }$

となり、これを解くことで $\boldsymbol{f}$ が求まります。（ $\boldsymbol{B}$ の特異値分解で解きます。）

$\boldsymbol{F}$ の性質として、 $rank(\boldsymbol{F})=2$ となるはずですが、対応点データが正確でない限り求めた解は $rank(\boldsymbol{F})=2$ になるとは限りません。よって、最終的な解のランクが２になるように事後調整など行います[3]。

実装

以下の実装に、エピポーラ線を描く部分を加えました。tukurutanoshi.hateblo.jp

f:id:mikekochang:20190307115542j:plain
３段目左はカメラ１から物体を見た画像（画像１）、右はカメラ２から物体を見た画像（画像２）です。
カメラ１はワールド座標系と同じ座標軸、カメラ２はｚ方向に4だけ移動したもので、画像１，２ともエピポーラ線は画像中心で交わっています。

コード

pyCodes/recon3d at master · misakikobayashi1984/pyCodes · GitHub

一部計算は[4]のサンプルを使用しています。

参考

[1]視覚の幾何学３　キャリブレーション http://www.wakayama-u.ac.jp/~wuhy/cv11.pdf

[2]正則行列と逆行列の諸性質,https://risalc.info/src/inverse-matrix.html#rank

[3]画像の三次元理解のための最適化計算[III], http://www.iim.cs.tut.ac.jp/~kanatani/papers/ieicefit3.pdf

[4]実践コンピュータビジョン, https://www.oreilly.co.jp/pub/9784873116075/index.html

2019-03-05

２つの画像の対応点からから３次元上の点を求める

理論

ある３次元上の点 $\boldsymbol{X}=(X,Y,Z,1)$ をカメラ内部パラメータ行列 $A$ とカメラ外部パラメータ行列 $[R|t]$ を用いて画像上に投影した点を $\boldsymbol{x}=(x,y,1)$ とすると、以下の関係が成り立ちます[1]。
${ s\boldsymbol{x}= A[R|t] \boldsymbol{X} \tag{1} }$
または、
${ \displaystyle s \begin{pmatrix} x_{1} \\ y_{1} \\ 1 \\ \end{pmatrix} = \begin{pmatrix} f_{x} && 0 && c_{x} \\ 0 && f_{y} && c_{y} \\ 0 && 0 && 1 \\ \end{pmatrix} \begin{pmatrix} r_{11} && r_{12 }&& r_{13} && t_{1}\\ r_{21} && r_{22 }&& r_{23} && t_{2}\\ r_{31} && r_{32 }&& r_{33} && t_{3}\\ \end{pmatrix} \begin{pmatrix} X \\ Y \\ Z \\ 1 \\ \end{pmatrix} \tag{2} }$
$s$ はスケール調整の定数です。カメラ内部行列の $f_{x}$ , $f_{y}$ はそれぞれピクセル単位の焦点距離、 $c_{x}$ , $c_{y}$ は画像中心を表します。カメラの固有パラメータという感じです。カメラ外部行列はカメラ姿勢を表し、３次元物体の座標をカメラ座標に変換する役割を持ちます。

同じカメラで別の視点から $\boldsymbol{X}$ を投影した画像が２枚（画像１、画像２）あり、それぞれのカメラ内部・外部行列を合わせたものを $M_{1}=A[R|t]_{1}$ , $M_{2}=A[R|t] _{2}$ 、画像上の投影点を $x_{1}$ , $x_{2}$ とすると、それぞれに対して(1)を適応して、

${ s\boldsymbol{x}_{1}= M_{1} \boldsymbol{X} \tag{3} }$

${ s\boldsymbol{x}_{2}= M_{2} \boldsymbol{X} \tag{4} }$

と表せます。
これを変形すると、

${ \begin{pmatrix} M_{1} && -\boldsymbol{x}_{1}&& 0\\ M_{1} && 0 && -\boldsymbol{x}_{2}\\ \end{pmatrix} \begin{pmatrix} \boldsymbol{X} \\ s_{1} \\ s_{2} \\ \end{pmatrix} = \boldsymbol{0} \tag{5} }$

と書けます。これを解くことで $\boldsymbol{X}$ が求まります。特異値分解で解くことが出来ます。

試してみた

３次元上の物体をカメラ行列 $M_{1}$ , $M_{2}$ を使ってそれぞれ画像1,画像2に変換し、その画像を３次元上に復元する例を作成しました。
カメラ行列のパラメータはそれぞれ変えることが出来ます。

f:id:mikekochang:20190305134838j:plain
Fig1. １段目左：元の３次元上の物体、右：復元した３次元物体
（２段目：元の物体のカメラ座標上の点。ここでは説明していない）
３段目左：画像１、右：画像２

元の３次元物体と復元結果が一致しています。

コード

pyCodes/recon3d at master · misakikobayashi1984/pyCodes · GitHub

一部計算は[2]のサンプルを使用しています。

参考

[1]キャリブレーションと３次元構成, http://opencv.jp/opencv-2svn/c/camera_calibration_and_3d_reconstruction.html

[2]実践コンピュータビジョン, https://www.oreilly.co.jp/pub/9784873116075/index.html

2019-03-01

輪郭を変えて顔を太らせる、痩せさせる

tukurutanoshi.hateblo.jp

と同じモチベーションで輪郭を変えてみました。太ります。痩せます。

f:id:mikekochang:20190301121313j:plain

Fig1. 左上：元画像、左下：太り、右下：痩せ

顔特徴点の抽出

顔パーツを選択して置換する - つくるって楽しい、と同じです。

輪郭の内側／外側に三角形を作る

顔の輪郭の特徴点とその近くの特徴点を三角形で結びます。輪郭の外側にも適当に点を配置し三角形を作ります。
f:id:mikekochang:20190301121420j:plain
Fig2. 三角形の配置

輪郭点の移動と変形

輪郭点を移動させます。ここでは顔の外側方向にに5pixcelとしました。三角形ごとに移動前後から対応するAffine変換を求めて変形させます。
f:id:mikekochang:20190301124942j:plain
Fig3. 左：移動前後の特徴点（青：移動前、赤移動後）、右：変形結果

コード

https://github.com/misakikobayashi1984/pyCodes/tree/master/moveFaceLandmark

2019-02-27

移動座標の対応点からAffine変換の係数を求める

３組の対応点からAffine変換を求める方法について。opencvだとcv2.getAffineTransform()に対応します。なぜ３点なのでしょう。

理論

Affine変換により $(x_{i}, y_{i})$ が $(x_{i}^{'}, y_{i}^{'})$ に変換される時次のように書けます。

${ \begin{pmatrix} x_{i}^{'} \\ y_{i}^{'} \\ \end{pmatrix} = \begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ \end{pmatrix} \begin{pmatrix} x_{i} \\ y_{i} \\ 1 \\ \end{pmatrix} \tag{1} }$

(1)を書き換えます。

${ \displaystyle \begin{pmatrix} x_{i} && y_{i} && 1 && 0 && 0 && 0 \\ 0 && 0 && 0 && x_{i} && y_{i} && 1 \\ \end{pmatrix} \begin{pmatrix} a_{11} \\ a_{12} \\ a_{13} \\ a_{21} \\ a_{22} \\ a_{23} \\ \end{pmatrix} = \begin{pmatrix} x_{i}^{'} \\ y_{i}^{'} \\ \end{pmatrix} \tag{2} }$

(2)は１対応点のみの場合ですが、N個の対応点がある場合は次のように拡張できます。

${ \displaystyle \begin{pmatrix} x_{1} && y_{1} && 1 && 0 && 0 && 0 \\ 0 && 0 && 0 && x_{1} && y_{1} && 1 \\ \vdots && && && && && \vdots \\ x_{N} && y_{N} && 1 && 0 && 0 && 0 \\ 0 && 0 && 0 && x_{N} && y_{N} && 1 \\ \end{pmatrix} \begin{pmatrix} a_{11} \\ a_{12} \\ a_{13} \\ a_{21} \\ a_{22} \\ a_{23} \\ \end{pmatrix} = \begin{pmatrix} x_{1}^{'} \\ y_{1}^{'} \\ \vdots \\ x_{N}^{'} \\ y_{N}^{'} \\ \end{pmatrix} \tag{3} }$

(3)を以下のように表記します。
${ X \boldsymbol{a} = X^{'} \tag{4} }$

これは最小二乗法の行列表現[1]ですね。 $X$ が正則ならば $\boldsymbol{a}$ は
${ \boldsymbol{a} = X^{-1} X^{'} \tag{5} }$
と求めることが出来ます。 $\boldsymbol{a}$ の未知数は６つなので、対応点が３つあれば(5)の方法で求めることが出来ます。

また、対応点が４つ以上の場合は $X$ が正則でないので逆行列が求められませんが、以下の式で解くことが出来ます。
${ \boldsymbol{a} = (X^{T} X)^{-1} X^{T} X^{'} \tag{6} }$

(6)は $(X^{T} X)^{-1}$ が正則な時の唯一つの解になるそうです。[1]

試してみた

時計回りに60度、x,y方向にそれぞれ5,3ずつ移動する変換Mを用意します。（(1)の形式の変換行列です。）

theta = 60
tx, ty = 5, 3
M = np.array([[np.cos(np.deg2rad(theta)), -np.sin(np.deg2rad(theta)), tx],
                 [np.sin(np.deg2rad(theta)), np.cos(np.deg2rad(theta)), ty]], float)

Mを出力。

M: 
[[ 0.5       -0.8660254  5.       ]
 [ 0.8660254  0.5        3.       ]]

３点の場合

座標変換前の点の集合。

X = np.array([[15, 5, 1],
              [25, 10, 1],
              [20, 15, 1]])

変換。

X_tfm = M.dot(X.T)

前後をプロット。赤：移動前、青：移動後。
f:id:mikekochang:20190227164507j:plain

(5)から変換行列を計算。Mと一致しています。

M_slv = np.linalg.inv(XX).dot(XX_tfm)

M_slv: 
[ 0.5       -0.8660254  5.         0.8660254  0.5        3.       ]

４点の場合

座標変換前の点の集合。変換。

X = np.array([[15, 5, 1],
              [25, 10, 1],
              [20, 15, 1],
              [18, 18, 1]])
X_tfm = M.dot(X.T)

(6)から変換行列を計算。Mと一致しています。

M_slv = np.linalg.inv(XX.T.dot(XX)).dot(XX.T).dot(XX_tfm)

M_slv: 
[ 0.5       -0.8660254  5.         0.8660254  0.5        3.       ]

コード

pyCodes/calcAffineLeastSquare at master · misakikobayashi1984/pyCodes · GitHub

参考

[1]最小二乗法の行列表現（単回帰，多変数，多項式）,https://mathtrain.jp/leastsquarematrix

2019-02-26

顔パーツを選択して置換する

検出した顔パーツをすべて含むよう部分をスワップする例[1]を参考に、目、口などパーツを選択して置換する例を作成しました。Fig1では、右目、左目、口を置換しています。眼が青く、唇ぷるんになっています。

f:id:mikekochang:20190301101345j:plain
Fig1.　左上：パーツ交換元、左下：交換先、右下：交換後

顔パーツの検出

[1]と同様にdlibという機械学習のライブラリの顔パーツ検出を使用しています。dlibの顔パーツ検出はOne Millisecond Face Alignment with an Ensemble of Regression Trees[2]という回帰で特徴点を検出する方法の実装だそうです。
顔の特徴点を68点で記述したデータセットiBUG 300-W datasetを学習した結果がdlibで提供されています。

f:id:mikekochang:20190301101700j:plain
Fig2.　特徴点検出結果

パーツを三角形ごとに変形

顔パーツを構成する特徴点から３点ずつ選びパーツを三角形の集合にします。左目なら４つの三角形になります。三角形ごとに対応するAffine変換を求めてパーツを変形します。三角形に分けることでAffine変換のみで画像が変形できます。

f:id:mikekochang:20190301101729j:plain
Fig3.　三角形の変形

パーツを合成

変形後のパーツを合成します。単純な置き換えではパーツのつなぎ目が目立ってしまうので、openCVのseamlessClone()を使うことで変形先画像に違和感なく画像を合成します。ポワソン方程式を画像合成に応用した手法Poisson Image Editing[3]の実装ということです。

f:id:mikekochang:20190301102048j:plain

Fig3.　左：単純な置き換え、右：つなぎ目を目立たないように合成

コード

pyCodes/facePartsSwap at master · misakikobayashi1984/pyCodes · GitHub

参考

[1]Face Swap using OpenCV ( C++ / Python ),<https://www.learnopencv.com/face-swap-using-opencv-c-python/>

[2]Vahid Kazemi and Josephine Sullivan(2014),One Millisecond Face Alignment with an Ensemble of Regression Trees,CVPR 2014
[3]Patrick Perez, Michel Gangnet, and Andrew Blake(2003),"Poisson Image Editing",SIGGRAPH2003