【LucasKanade法とは】最小二乗法とガウシアンピラミッドによる推定の原理・計算式

オプティカルフロー推定におけるLucasKanadeは？最小二乗法とガウシアンピラミッドによる推定の原理・計算式について解説します。

LucasKanade法とは
【最小二乗法】オプティカルフロー推定
【問題点】最小二乗法だけでは計算精度が悪い
【山登り法】大きい移動量も精度良く推定
LucasKanade法の流れまとめ
Python版OpenCVによる実装例
関連ページ

LucasKanade法とは

LucasKanade法は1981年にLucas、金出さんらによって作られたオプティカルフローを推定するアルゴリズムの1つです。LucasKanade法では、「最小二乗法」と「山登り法」を組み合わせることでより良いオプティカルフロー $(u, v)$ の値を推定します。

–	説明
①	最小二乗法でオプティカルフロー(移動物体の動きベクトル)を推定
②	最小二乗法で求めた推定解は実際に使うと真値と大きくズレてしまう
③	山登り法を使って推定解を真値に近づける

前回記事

【オプティカルフローとは】推定の原理・特徴・計算式

オプティカルフローとは？移動量の推定方法、原理、計算式についてまとめました。

【最小二乗法】オプティカルフロー推定

1画素あたりの拘束方程式は1つであるため、このままではオプティカルフロー $(u, v)$ を一意に求めることが出来ません。
そこで、LucasKanade法では「近い画素も同じ動きをする」という仮定を用いることで、方程式の個数を増やし、解 $(u,v)$ を推定します。

例

画像I上の隣接し合う2点 $p_1, p_2$ が同じ動きをすると仮定すれば、下式のようにその2点の画素に対する拘束式の $(u, v)$ は同じになります。

(1) $\begin{eqnarray*} \left{ \begin{array}{r@{\,}c@{\,}r@{\,}c@{\,}r@{\;\leq\;}r} I_{x1}(p_1)u+I_{y1}(p_1)v&=&-\Delta I_1\\ I_{x2}(p_2)u+I_{y2}(p_2)v&=&-\Delta I_2 \end{array}% \right. \end{eqnarray*}$

よって、この2式を連立して解けば、未知数が2個、方程式も2個なので $(u, v)$ を一意に求めることができます。
これがLucasKanade法の基本的な考え方です。
ところが、実際には隣接する2点だけでなく、もっと広い領域で動きが同じであると仮定します。
例えば $3\times 3(=9)$ の画素の動きがすべて同じだとすれば，下記のような9つの方程式が得られます。

(2) $\begin{eqnarray*} \left{ \begin{array}{r@{\,}c@{\,}r@{\,}c@{\,}r@{\;\leq\;}r} I_{x1}(p_1)u+I_{y1}(p_1)v&=&-\Delta I_1\\ I_{x2}(p_2)u+I_{y2}(p_2)v&=&-\Delta I_2\\ \vdots\\ I_{x9}(p_9)u+I_{y9}(p_9)v&=&-\Delta I_9 \end{array}% \right. \end{eqnarray*}$

この場合，未知数2個、方程式9個となります。よって、9つの方程式を満たす2つの未知数 $(u, v)$ は一意に求めることができません。
これをグラフに表すと下画像のようになります。

未知数が求まらないというのは、上のグラフで説明すると9つの点を通る直線は存在しないということです。
しかし、9つの点は全く散らばっているわけではありません。（近い点同士なので画素の動きは似ている）
よって，9つの方程式を近似的に満足する未知数uとvは求めることができるはずです。
この事を先程のグラフで書き表すと下のような近似直線を求めることになります。

この近似直線を求める作業が、最初に述べた「最小二乗法」です。
このように「最小二乗法」を利用してオプティカルフローを推定することが「LucasKanade法」の基本的な考え方の1つです。
次の項目では、この流れを数式（線形代数）でまとめます。
先程の最小二乗法によってオプティカルフロー $(u, v)$ を推定する一連の流れを線形代数の式で表していきます。
$M$ 個の画素が同じ動きをすると仮定すれば、つぎの $M$ 個のオプティカルフロー拘束式が得られます。

(3) $\begin{eqnarray*} \left{ \begin{array}{r@{\,}c@{\,}r@{\,}c@{\,}r@{\;\leq\;}r} I_{x1}u+I_{y1}v&=&-\Delta I_1\\ I_{x2}u+I_{y2}v&=&-\Delta I_2\\ \vdots\\ I_{xM}u+I_{yM}v&=&-\Delta I_M \end{array}% \right. \end{eqnarray*}$

ただし

(4) $\begin{eqnarray*} \left{ \begin{array}{r@{\,}c@{\,}r@{\,}c@{\,}r@{\;\leq\;}r} I_{xi}=I_{x}(p_i)\\ I_{yi}=I_{y}(p_i)\\ \Delta I_i=\Delta I(p_i) \end{array}% \right. \end{eqnarray*}$

M個の方程式を1つの線形方程式にまとめると次のようになります。

(5) $\begin{eqnarray*} AX=B \end{eqnarray*}$

ただし

(6) $\begin{eqnarray*} A= \begin{bmatrix} I_{x1} & I_{y1} \\ I_{x2} & I_{y2} \\ \vdots & \vdots \\ I_{xM} & I_{yM} \end{bmatrix} , X= \begin{bmatrix} u \\ v \\ \end{bmatrix} , B= \begin{bmatrix} \Delta I_1 \\ \Delta I_2 \\ \vdots \\ \Delta I_M \end{bmatrix} \end{eqnarray*}$

このとき，以下の定理を用いればオプティカルフローを推定することが出来ます。

定理

(5)式に対して次式を解けば $(u,v)$ の近似値を求めることができます。

(7) $\begin{eqnarray*} A^TAX=A^TB \end{eqnarray*}$

ただし

(8) $\begin{eqnarray*} |A^TA|\neq0 \end{eqnarray*}$

つまり、 $A^TA$ の逆行列が存在すれば、最小二乗解(オプティカルフローの推定値)が一意に求まります。
また、線形代数の特性より、以下のようにして得られた解が良いか悪いかを調べることが出来ます。

–	説明
①	$\|A^TA\|$ が大きい
②	$A^TA$ の2つの固有値 $\lambda 1, \lambda 2$ が微小でない
③	$\lambda 1/\lambda 2$ が大きすぎない( $\lambda 1>\lambda 2$ )

$A^TA$ を見てやるだけで簡単に解の安定性を調べることが出来ます。それにより、その画素が追跡しやすいか否かがわかり、特徴追跡をおこなうときに役立ちます。( $|A^TA|$ の大きい点，つまり特徴点だけを抽出)
物体追跡で一番大事なことは、追跡対象を選択することです。よって $|A^TA|$ の大きい点を対象に選べば
安定的に追跡できることになります。
線形代数で表現することで、現代制御理論のように簡単に安定判別ができます。

計算例

隣接しあう3つの画素の動きが全く同じであるとし、以下の3つのオプティカルフロー拘束式が得られたとします。

(9) $\begin{eqnarray*} \left{ \begin{array}{r@{\,}c@{\,}r@{\,}c@{\,}r@{\;\leq\;}r} I_{x1}u+I_{y1}v&=&-\Delta I_1\\ I_{x2}u+I_{y2}v&=&-\Delta I_2\\ I_{x3}u+I_{y3}v&=&-\Delta I_3 \end{array}% \right. \end{eqnarray*}$

ただし， $I_{x1}=1, I_{y1}=-1, \Delta I_1=1, I_{x2}=1, I_{y2}=-2, \Delta I_2=0.1$ ,
$I_{x3}=3, I_{y3}=-3, \Delta I_3=9.8$ とします。
このとき，線形方程式の定数行列 $A, B$ と行列 $X$ はつぎのようになる

(10) $\begin{eqnarray*} A= \begin{bmatrix} 1 & -1 \\ 1 & -2 \\ 3 & 4 \end{bmatrix} , X= \begin{bmatrix} u \ v \ \end{bmatrix} , B= \begin{bmatrix} 1 \ 0.1 \ 9.8 \end{bmatrix} \end{eqnarray*}$

$A^TA\neq0$ となるので，(1)式で最小二乗解Xを求めることができます。

(11) $\begin{eqnarray*} X=(A^TA)^{-1}A^TB\approx \begin{bmatrix} 1.99 \\ 0.95 \\ \end{bmatrix} \end{eqnarray*}$

よって画素の移動量は $u=1.99, v=0.95$ となります。

【問題点】最小二乗法だけでは計算精度が悪い

オプティカルフロー拘束式は、以下の仮定①②の下で成立するものでした。

–	説明
仮定①	移動前後で追跡したい点の画素値は変化しない
仮定②	追跡したい点の移動量は小さい（<1画素）

しかし、当然ながら現実の世界では，これらの仮定はほとんど成立しません。
特に仮定②は成立しない場合がほとんどでしょう。（移動物体は1フレームで1画素以上移動することが多い）
仮定②が成立しない中でも最小二乗法でオプティカルフローを推定することができますが、その解は真値と大きくズレてしまいます。
そこで、最小二乗法だけでなく「山登り法」という方法を用いて仮定②が成立しない場合でもオプティカルフローを正確に推定します。
「山登り法」については次の項目で紹介します。

【山登り法】大きい移動量も精度良く推定

先程の項目でも紹介しましたが、オプティカルフロー拘束式は以下の2つの仮定で成立します。

–	説明
仮定①	移動前後で追跡したい点の画素値は変化しない
仮定②	追跡したい点の移動量 $(u, v)$ は小さい（<1画素）

しかし、 $u,v$ が大きい場合は仮定②が崩れてしまい、推定解が真値と大きくズレてしまいます。
ここで、もし $u,v$ の近似解 $u',v'$ が得られ、かつその精度が良ければ万々歳です。これを数式で表すと次のようになります。

(12) $\begin{eqnarray*} \left{ \begin{array}{r@{\,}c@{\,}r@{\,}c@{\,}r@{\;\leq\;}r} |u-u'|<0.5\\ |v-v'|<0.5 \end{array}% \right. \end{eqnarray*}$

ここで

(13) $\begin{eqnarray*} \left{ \begin{array}{r@{\,}c@{\,}r@{\,}c@{\,}r@{\;\leq\;}r} \Delta u=u-u'\\ \Delta v=v-v' \end{array}% \right. \end{eqnarray*}$

とおくと、u,vは

(14) $\begin{eqnarray*} \left{ \begin{array}{r@{\,}c@{\,}r@{\,}c@{\,}r@{\;\leq\;}r} u=u'+\Delta u\\ v=v'+\Delta v \end{array}% \right. \end{eqnarray*}$

と表せます。このとき、仮定1)より

(15) $\begin{eqnarray*} I(x+u,y+v)=H(x,y)\\ I(x+u'+\Delta u,y+v'+\Delta v)=H(x,y)\\ I(x+u'+\Delta u,y+v'+\Delta v)-H(x,y)=0 \end{eqnarray*}$

となるので、この式を1次項までテイラー展開してやると

(16) $\begin{eqnarray*} \frac{\partial I}{\partial x}\Delta u+\frac{\partial I}{\partial y}\Delta v+I(x+u',y+v')-H(x,y)=0 \end{eqnarray*}$

という方程式が得られます。この式は下式のオプティカルフロー拘束式とよく似ています。

(17) $\begin{eqnarray*} \frac{\partial I}{\partial x}\Delta u+\frac{\partial I}{\partial y}\Delta v+\Delta I=0 \end{eqnarray*}$

よって、近似解 $u',v'$ が得られれば、最小二乗法で $\Delta u,\Delta v$ が推定でき、真値に近い $u,v$ が計算できます。

あとは近似解 $u',v'$ をどのように求めるかが課題となります。そこで登場するのが「山登り法」です。
「山登り法」では、近似解 $u',v'$ を求めるために、下図のような多重解像度画像のピラミッド（通称ガウンシアンピラミッド）を用いられます。

このピラミッドは画像 $H,I$ を元の解像度(一番下)から1/2ずつ小さくした画像を重ねています。
このように解像度を落とすことで、画素を物理的に大きくし、強引に仮定②を満足させます。
例えば $6\times 6$ の画像ついて考えると下図のようになります。
※赤色：移動物体

■元の解像度・・・移動量は2[px]

■解像度1/2・・・移動量は1[px]

このように解像度を落とし続ければ、いずれ画素の移動量が1以下になり仮定2を満たせるので、拘束式を解いてオプティカルフローが求まります。このオプティカルフローが近似解 $u',v'$ となります。
この近似解を使って、解像度が一段上の画像でまた近似解を求めます。この作業を元画像の解像度まで繰り返し行うことで、最終的に真値に近いオプティカルフローを求めることが出来ます。これが山登り法です。

LucasKanade法の流れまとめ

LucasKanade法の処理の流れをまとめると以下のようになります。

–	説明
1	画像 $H,I$ に平滑フィルタをかけ、間引きにより解像度が1/2(画素数は1/4)になる画像 $H_1とI_1$ を作ります。この処理を繰り返して $H_2,I_2,H_3, I_3...H_n, I_n$ を作成します。そして、画素の大きさが画像中の画素の最大移動量より大きくなるまで画像 $H, I$ の画像ピラミッドを作成します。
2	最初に, すべての画素における移動量の近似解を $u=0, v=0$ として画像のレベルを $i=n$ とします。
3.1	$H_i, I_i$ における各画素の $\Delta u ,\Delta v$ を(16)式で計算します。
3.2	各画素の移動量 $(u,v)=(u'+\Delta u, v'+\Delta v)$ を(5)式で計算する.
3.3	解像度のレベルが一段高い画像 $H_{i-1}, I_{i-1}$ の画像の近似的な移動量 $u',v'$ を計算します。 ① $x,y$ がともに偶数値の画素の近似移動量 $(u', v')=I$ の $(x/2, y/2$ )画素の $(2u,2v)$ ②それ以外の画素に対して、隣(上下左右)の $(x, y)$ ともに偶数値の画素の移動量の平均値を $(u',v')$ とします。 ③カウント変数を更新します。 $(i=i-1)$
4	$i=0$ になるまで3.1～3.3の処理を繰り返します。