DIS(Dense Inverse Search)の仕組み・実行とパラメータ調整

## DISの仕組み
### DISの概要
DIS(Dense Inverse Search)はLucas-Kanade法とHorn-Schunck法を組み合わせることでDense(ピクセル単位のベクトルの計算)にOptical Flowを推論するアルゴリズムです。DISではLucas-Kanade法における計算負荷の軽減にあたって、Inverse Searchが行われます。

![DIS_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3936087/46fadddf-de89-22bf-cd75-65dc54bdbebc.png)
DIS論文 Algorithm.1

DISでは上記のようなアルゴリズムに基づいて密なOptical Flowの計算が行われます。ピラミッドのスケール毎に1)～5)の処理が実行されます。1)～5)の処理の詳細については次項以降で詳しく確認します。


### 1) Creation of a grid
DISにおけるスケールごとのイテレーションではまずスケールに対応したgridを生成します。このgirdに基づいて疎なOptical FlowのアルゴリズムであるInverse Search(Lucas-Kanade法)が実行されます。

gridの生成では、「ピラミッドのスケール」と「パッチの重複(overlap)の度合い」に基づいてLucas-Kanade法を実行する点が決定されます。DISの論文ではスケール$s$で生成されるLucas-Kanade法を計算する格子点の数を$N_{s}$のように定義されます。$N_{s}$は後ろの処理を表現するにあたってよく出てくるので抑えておくと良いです。


### 2) Initialization
次に1)によって決めた$N_{s}$個の格子点に対し、ベクトルの初期化を行います。

それぞれの格子点のインデックスを$i \in N_{s}$、$i$番目の格子点のOptical Flowを$\mathbf{u}_{i}$と表すと、$i$番目の格子点のOptical Flowは下記のような式で計算されます。

```math
\mathbf{u}_{i} = U_{s+1} \left( \frac{\mathbf{x}}{\theta_{sd}} \right) \cdot \theta_{sd} \in \mathbb{R}^{2}
```

上記の式は「$\mathbf{u}_{i}$の初期ベクトルが1つ粗いスケールの密なOptical Flowの値に基づいて決定されること」を意味します。

$\theta_{sd}$はイテレーションにおけるダウンダンプリングの比率に対応し、$\theta_{ss} \geq \theta_{sd} \geq \theta_{sf}$が成立しています。

また、スケールのイテレーションの一番初めの$\theta_{sf}$では$U_{s+1}$が存在しないので、$\mathbf{u}_{i}$にはゼロベクトルが与えられます。


### 3) Inverse Search
DISではLucas-Kanade法の改良版であるInverse Searchを用いて疎なOptical Flowの計算が行われます。以下、Inverse Searchについて詳しく確認します。

```math
\begin{align}
\Delta \mathbf{u} &= \mathrm{argmin}_{\Delta \mathbf{u}'} \sum_{\mathbf{x} \in X} \left( I_{t+1}(\mathbf{x} + \mathbf{u} + \Delta \mathbf{u}') - T(\mathbf{x}) \right)^{2} \quad (1) \\
\mathbf{u} & \longleftarrow \mathbf{u} + \Delta \mathbf{u}
\end{align}
```

Lucas-Kanade法では上記のような式に基づく繰り返し演算によってOptical Flowの$\mathbf{u}$の計算が行われます。上記の演算では$(1)$式を解くにあたってヘッセ行列を計算する必要がある一方で、この演算の計算量は他の演算に比較して計算量の多いものとなります。Lucas-Kanade法の導出の詳細については下記で詳しく取り扱ったので、下記をご確認ください。

https://qiita.com/Chi_corp_123/items/dab2e1c8206b451ddd00

「ヘッセ行列の計算量が大きい」という課題の解決にあたって、Inverse SearchではSource画像ではなくTarget画像を動かすことでヘッセ行列を固定します。

```math
\begin{align}
\Delta \mathbf{u} &= \mathrm{argmin}_{\Delta \mathbf{u}'} \sum_{\mathbf{x} \in X} \left( T(\mathbf{x} - \Delta \mathbf{u}') - I_{t+1}(\mathbf{x} + \mathbf{u})\right)^{2} \quad (2) \\
\mathbf{u} & \longleftarrow \mathbf{u} + \Delta \mathbf{u}
\end{align}
```

上記に基づいて$\Delta \mathbf{u}$を得る場合、繰り返し演算におけるヘッセ行列は全て同じものを用いることができ、計算が1度で済みます。$(2)$式ではTarget画像とSource画像の項の順序を入れ替えていることから、このような手法をInverse Searchと言います。

このような演算を行うことでLucas-Kanade法の高速化を実現することが可能になります。


### (4) Densification
DensificationではInverse Search(Lucas-Kanade法)によって計算された疎なOptical Flow(特定の位置に対応するベクトルのみを出力)の$\mathbf{u}_{i}$を元に密なOptical Flow(画像の全てのピクセルに対応するベクトルを出力)を得る処理を行います。

密なOptical Flowの$U_{s}(\mathbf{x})$は位置$\mathbf{x}$の周辺の疎なOptical Flowの重みづけ和に基づいて下記のような計算式によって計算されます。

```math
\begin{align}
U_{s}(\mathbf{x}) &= \frac{1}{Z} \sum_{i=1}^{N_s} \frac{\lambda_{i, \mathbf{x}}}{\max{(1, ||d_{i}(\mathbf{x})||^{2})}} \cdot \mathbf{u}_{i} \\
d_{i}(\mathbf{x}) &= I_{t+1}(\mathbf{x}+\mathbf{u}_{i}) - T(\mathbf{x}) \\
Z &= \sum_{i=1}^{N_s} \frac{\lambda_{i, \mathbf{x}}}{\max{(1, ||d_{i}(\mathbf{x})||^{2})}}
\end{align}
```

重みの計算にあたってはまず$\mathbf{x}$の周辺のOptical Flowのみを用いるにあたって$\lambda_{i, \mathbf{x}}$を用います。$\lambda_{i, \mathbf{x}}$は$i$番目の疎なOptical Flowが位置$\mathbf{x}$の周辺にある場合は$1$、ない場合は$0$を出力します。

また、$d_{i}(\mathbf{x})$によって重みの計算が行われます。$d_{i}(\mathbf{x})$は移動先のSourceの画素値とTargetの画素値が大きく異なる場合大きな値を出力するので、$d_{i}(\mathbf{x})$の二乗の逆数を重みとすることで画素値が近い$\mathbf{u}_{i}$を優先的に採用します。

$Z$は正規化に用いる項であり、この項を採用することでの$\lambda_{i, \mathbf{x}}=1$の項の数に関わらずベクトルを同じスケールで表現することができます。


### (5) Variational Refinement
DISのアルゴリズムによって生成されるピクセル毎のFlow(ベクトル)の集合を$U \in \mathbb{R}^{W \times H \times 2}$のようにおくとき、生成されたベクトルについて下記のようにenergyの$E(U)$を定義します。

```math
\begin{align}
E(U) &= \int_{\Omega} \sigma \Psi(E_{I}) + \gamma \Psi(E_{G}) + \alpha \Psi(E_{S}) d \mathbf{x} \quad (3) \\
\Psi(a^{2}) &= \sqrt{a^{2} + \epsilon^{2}}  \\
\mathbf{x} &= \left( \begin{array}{c} x \\ y \end{array} \right)
\end{align}
```

$\Omega$による積分によって1ピクセルに関する計算を画像全体に適用することが読み取れます。Variational Refinementでは上記が小さくなるように$U$の値を調整します。$(3)$式に出てくる$E_{I}$は下記のように定義される指標です。
```math
\begin{align}
E_{I} &= \mathbf{u}^{\mathrm{T}} \bar{\mathbf{J}}_{0} \mathbf{u} \in \mathbb{R} \\
\mathbf{J}_{0} &= \beta_{0} (\nabla_{3} I) (\nabla_{3}^{\mathrm{T}} I) \in \mathbb{R}^{3 \times 3} \\
\mathbf{u} &= \left( \begin{array}{c} u \\ v \\ 1 \end{array} \right), \quad \nabla_{3} = \left( \begin{array}{c} \displaystyle \frac{\partial}{\partial x} \\ \displaystyle \frac{\partial}{\partial y} \\ \displaystyle \frac{\partial}{\partial t} \end{array} \right)
\end{align}
```

上記はBrightness Constancy($(\nabla_{3} I) \mathbf{u}=0$)に近づくように定義される指標です。$I \in \mathbb{R}$は特定の1ピクセルにおけるグレースケールの画素値を表すと理解しておくと良いです。$E_{I}$と同様に$E_{G}$は下記のように定義されます。

```math
\begin{align}
E_{G} &= \mathbf{u}^{\mathrm{T}} \bar{\mathbf{J}}_{xy} \mathbf{u} \in \mathbb{R} \\
\bar{\mathbf{J}}_{xy} &= \beta_{x} (\nabla_{3} I_{dx}) (\nabla_{3}^{\mathrm{T}} I_{dx}) + \beta_{y} (\nabla_{3} I_{dy}) (\nabla_{3}^{\mathrm{T}} I_{dy}) \in \mathbb{R}^{3 \times 3}
\end{align}
```

また、Smoothnessに基づく項である$E_{S}$は下記のように定義されます。

```math
E_{S} = ||\nabla u||^{2} + ||\nabla v||^{2}
```

上記の$u, v$はそれぞれ特定の1ピクセルにおけるDISの出力結果に対応します。



## DISの実行とパラメータ調整
### DISの実行
DISの実行にあたってはOpenCVを用いることが多いです。OpenCVにはPythonバインディングも用意されており、Pythonからも使用することができます。以下ではPythonを用いた実行について確認します。

```python
import cv2 as cv

image1_gray = cv.cvtColor(image1, cv.COLOR_BGR2GRAY)
image2_gray = cv.cvtColor(image2, cv.COLOR_BGR2GRAY)

inst = cv.DISOpticalFlow.create(cv.DISOPTICAL_FLOW_PRESET_MEDIUM)
flow = inst.calc(image1_gray, image2_gray, None)
```

DISは上記のようなPython実装を用いることで実行することができます。`image1`と`image2`はRGBの3チャネルを持った画像を読み込んだものであり、DISのアルゴリズムによって`image1`から`image2`の変化を各ピクセル毎に(dense)縦横の2Dのベクトルで出力します。

ライブラリのインストールにあたっては`pip install opencv-python`を実行することで入手することができます。


### DISのパラメータ調整
`DISOpticalFlow`クラスには下記のようなパラメータの確認を行うメソッド(`get`)が用意されています。

| メソッド | 概要 |
|:-----------|:-----------|
|`getFinestScale()`|フローが計算されるガウシアンピラミッドの最も細かいレベル|
|`getPatchSize()`|マッチングのための画像パッチのサイズ（ピクセル単位）|
|`getPatchStride()`|隣接するパッチ間のストライド|
|`getGradientDescentIterations()`|Inverse Searchにおける、勾配降下法の最大反復回数|
|`getVariationalRefinementIterations()`|スケールごとのVariational Refinementの固定点反復の数|
|`getVariationalRefinementAlpha()`|滑らかさ項の重み、$(1)$式の$\alpha$に対応|
|`getVariationalRefinementDelta()`|色彩恒常性項の重み、$(1)$式の$\sigma$に対応|
|`getVariationalRefinementGamma()`|勾配恒常性項の重み、$(1)$式の$\gamma$に対応|
|`getUseMeanNormalization()`|パッチ距離を計算する際に、パッチの平均正規化を使用するかどうか|
|`getUseSpatialPropagation()`|良好なオプティカルフローベクトルの空間伝搬を使用するかどうか|

上記を実際に実行するにあたっては、下記のようなコードを実行すれば良いです。

```python
import cv2 as cv

inst = cv.DISOpticalFlow.create(cv.DISOPTICAL_FLOW_PRESET_MEDIUM)

print(inst.getPatchSize())
```

DISで用いられるパッチサイズはデフォルトでは$8 \times 8$であるので、特に指定のない場合は上記を実行した場合`8`が出力されます。パラメータのデフォルト値の確認にあたって、一覧にまとめたそれぞれのメソッドの実行結果を以下で確認しました。

```python
import cv2 as cv

inst = cv.DISOpticalFlow.create(cv.DISOPTICAL_FLOW_PRESET_MEDIUM)

print(inst.getFinestScale())
print(inst.getPatchSize())
print(inst.getPatchStride())
print(inst.getGradientDescentIterations())
print(inst.getVariationalRefinementIterations())
print(inst.getVariationalRefinementAlpha())
print(inst.getVariationalRefinementDelta())
print(inst.getVariationalRefinementGamma())
print(inst.getUseMeanNormalization())
print(inst.getUseSpatialPropagation())
```

・実行結果
```
1
8
3
25
5
20.0
5.0
10.0
True
True
```

ここまではパラメータの確認に用いる`get`について確認しましたが、パラメータの指定にあたっては`set`が用意されています。`set`のメソッドはたとえば下記のように実行することができます。

```python
import cv2 as cv

inst = cv.DISOpticalFlow.create(cv.DISOPTICAL_FLOW_PRESET_MEDIUM)

print(inst.getPatchSize())
inst.setPatchSize(2)
print(inst.getPatchSize())
```

・実行結果
```
8
2
```

DISにおけるパラメータの調整については上記のように`set`のメソッドを用いることで行うことができます。また、公式のドキュメントは下記より確認することができるので合わせてご確認ください。

https://docs.hsp.moe/OpenCV/453/cpp/ja/classcv_1_1_d_i_s_optical_flow.html