はじめに
以前投稿した記事で、埋め込み境界法で用いたい境界点を得るために、レベルセット法による画像のセグメンテーションを用いました。
このときはとりあえず使えればいいやということで、公開されているプログラムを単に実行させてレベルセット関数を得ていました。
それでもまあ問題ないと思うのですが、中身や理屈がわかっていないプログラムをそのまま使うのはあまりよろしくないと思いますので、多少勉強してきました。
本記事では、以下の事を行ないました。
- Chan Vese Level Set Method による画像のセグメンテーションの簡単な説明
- 鴨川の衛星写真にこの手法を用いて飛び石の境界を判定
- VP法(Volume Penalization法)の簡単な説明・レベルセット法の利用
- 浅水方程式への実装
- 飛び石周りの流れ・水位場の数値計算
- Mayaviで可視化
あとは適当に考察とかもしてみます。
Chan Vese Level Set Method による画像のセグメンテーション
そもそも、画像のセグメンテーションとは、与えられた画像を分割し注目したい部分を切り出すことです。最近では、機械学習なんかで非常にすごいことになっていそうな気がします。医療や自動運転への応用など、どうやらいろいろと役に立っているようです。
この画像のセグメンテーション手法の一つが、レベルセット法によるセグメンテーションです。レベルセット法とは、レベルセット関数という、境界面からの符号付き距離関数を用いて物体の内外を表現する方法です。この曲線と画像の値(例えば輝度)をもちいた評価関数(後述)を最小化させるような曲線$C= \{ (x,y) : \phi(x,y)=0 \}$を求めて、物体の境界を取り出します。
レベルセット法の符号は、とりあえず今は Chan & Vese (2001) に従います。
(前回用いたPythonコードでは逆。計算の符号が逆になるだけで問題ない)
|
| レベルセット関数 $\phi$ の概念図 |
エネルギー汎関数
$\Omega \subset \mathbb{R}^2 \cdots$有界な開集合、$\partial \Omega\cdots \Omega$の境界、$u_0: \bar{\Omega} \rightarrow \mathbb{R}\cdots$画像、$C$は$\Omega$内の曲線とする。Chan-Vese レベルセット法では、次のようなエネルギー汎関数を最小化することを目指します。
\newcommand{\pdv}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\bm}[1]{ \boldsymbol{#1} }
\newcommand{\dd}{ \mathrm{d} }
\newcommand{\ve}{ \varepsilon }
\newcommand{\odv}[2]{\frac{\dd #1}{\dd #2}}
\begin{align}
F(c_1,c_2,C)=&\mu\cdot\mathrm{Length}(C)+\nu\cdot\mathrm{Area}(\mathrm{Inside}(C)) \\
&+\lambda_1 \int_{\mathrm{Inside}(C)}|u_0(x,y)-c_1|^2\dd x \dd y \\
&+\lambda_2 \int_{\mathrm{Outside}(C)}|u_0(x,y)-c_2|^2\dd x \dd y
\end{align}
ここで、$\mu\geq0, \nu\geq 0, \lambda_1,\lambda_2>0$は定数のパラメータ(論文中ではとりあえず$\lambda_1=\lambda_2=1,\nu=0$としている)。$c_1$は曲線内部の輝度の平均値、$c_2$は外部の輝度平均値。
この汎関数の右辺各項の意味は
- 第1項 ... 曲線の長さ
- 第2項 ... 内部の面積
- 第3項&第4項 ... 曲線内部・外部の値とその平均との差分(分散のようなもの)
で、これを最小化することで輪郭を抽出します。
第3項と第4項による効果は、簡単なケース、例えば輝度値が二値の画像を考えると理解しやすいです。
この項を最小化する曲線が、輝度値を二分する境界ではないときには、どちらかの項が正になってしまいます。この項を最小化する曲線は、輝度値の境界線そのものとなります。
レベルセット関数による記述
各領域をレベルセット関数で表すと、
\begin{align}
C &= \{ (x,y)\in\Omega : \phi(x,y)=0 \} \\
\mathrm{Inside}(C) &= \{ (x,y)\in\Omega : \phi(x,y)>0 \} \\
\mathrm{Outside}(C) &= \{ (x,y)\in\Omega : \phi(x,y)<0 \}
\end{align}
のように書ける。これを使って、$C$を$\phi$に置き換えた汎関数は以下のようになります。
\begin{align}
F(c_1,c_2,\phi) =& \mu\int_\Omega \delta(\phi(x,y))|\nabla\phi(x,y)|\dd x \dd y \\
&+\nu\int_\Omega H(\phi(x,y))\dd x \dd y \\
&+\lambda_1 \int_\Omega |u_0(x,y)-c_1|^2 H(\phi(x,y)) \dd x \dd y \\
&+\lambda_2 \int_\Omega |u_0(x,y)-c_2|^2 (1-H(\phi(x,y)))\dd x \dd y
\end{align}
ここで、$H$はヘヴィサイド関数で
H(z)=\left\{\begin{array}{ll}
1, & \mathrm{if} \: z\geq 0, \\
0, & \mathrm{if} \: z < 0,
\end{array}\right. \quad \delta(z)=\frac{\dd}{\dd z} H(z)
です。
オイラーラグランジュ方程式の導出
未知関数$\phi$を求めるための オイラーラグランジュ方程式を導出します。まず、扱いやすいようにヘヴィサイド関数とデルタ関数を正則化しておきます。
$H_\varepsilon$を$H$の$C^2(\bar{\Omega})$正則化とし、$\delta_\varepsilon=H_\varepsilon'$とします。それぞれを正則化したものに変えた汎関数を $F_\varepsilon$としておきます。
この汎関数 $F_\varepsilon$ を最小化(停留化)するような式を求めるために、解析力学を勉強したときを思い出して汎関数の変分を考えます。$c_1,c_2$をfixして、$\phi$に関した変分は以下のようになります。
\begin{align}
\delta F_\ve[\phi] &= F_\ve[\phi+\delta\phi]-F_\ve[\phi] \\
& = \int_\Omega [\mu \left\{ \delta_\ve(\phi+\delta\phi)|\nabla(\phi+\delta\phi)|-\delta_\ve(\phi)|\nabla\phi| \right\} \\
& \quad \quad \quad + \nu \{H_\ve(\phi+\delta\phi)-H_\ve(\phi)\} \\
& \quad \quad \quad + \lambda_1 |u_0 - c_1|^2 \{H_\ve(\phi+\delta\phi)-H_\ve(\phi)\} \\
& \quad \quad \quad - \lambda_2 |u_0 - c_2|^2 \{H_\ve(\phi+\delta\phi)-H_\ve(\phi)\} ] \dd x \dd y
\end{align}
積分中の第2,3,4項は単純に$H_\ve$の微分すなわち$\delta_\ve(\phi)$に係数(と$\delta\phi$)をかけているだけ。第1項はすこし面倒で、以下のように変換しておきます。$\delta$の二次以上の項は省いています。
\begin{align}
& \delta_\ve(\phi+\delta\phi) \simeq \delta_\ve(\phi) + \odv{\delta_\ve}{z}(\phi) \delta\phi \\
& |\nabla(\phi+\delta\phi)| \simeq |\nabla\phi| + \nabla |\nabla\phi| \cdot \delta\nabla\phi = |\nabla\phi| + \frac{\nabla\phi}{|\nabla\phi|}\cdot \nabla\delta\phi \\
\end{align}
より
\begin{align}
&\int_\Omega \mu \left\{ \delta_\ve(\phi+\delta\phi)|\nabla(\phi+\delta\phi)|-\delta_\ve(\phi)|\nabla\phi| \right\}\dd x \dd y \\
&\simeq \int_\Omega \mu \left\{ \left( \delta_\ve(\phi) + \odv{\delta_\ve}{z}(\phi) \delta\phi \right) \left( |\nabla\phi| + \frac{\nabla\phi}{|\nabla\phi|}\cdot \nabla\delta\phi \right) -\delta_\ve(\phi)|\nabla\phi| \right\}\dd x \dd y \\
&= \int_\Omega \mu \left\{ |\nabla\phi| \odv{\delta_\ve}{z}(\phi) \delta\phi + \delta_\ve(\phi) \frac{\nabla\phi}{|\nabla\phi|}\cdot \nabla\delta\phi \right\}\dd x \dd y \\
\end{align}
ここで、
\begin{align}
&\int_\Omega \left\{ \delta_\ve(\phi) \frac{\nabla\phi}{|\nabla\phi|}\cdot \nabla\delta\phi \right\}\dd x \dd y \\
& = \int_\Omega \nabla \cdot \left\{ \delta_\ve(\phi) \frac{\nabla\phi}{|\nabla\phi|} \delta\phi \right\}\dd x \dd y - \int_\Omega \nabla \cdot \left\{ \delta_\ve(\phi) \frac{\nabla\phi}{|\nabla\phi|} \right\} \delta\phi \dd x \dd y \\
& = \int_{\partial \Omega} \frac{\delta_\ve(\phi)}{|\nabla\phi|} \pdv{\phi}{n} \delta\phi \dd S
- \int_\Omega \delta_\ve(\phi) \nabla \cdot \left( \frac{\nabla\phi}{|\nabla\phi|} \right) \delta\phi \dd x \dd y
- \int_\Omega \nabla\delta_\ve(\phi) \cdot \frac{\nabla\phi}{|\nabla\phi|} \delta\phi \dd x \dd y \\
& = \int_{\partial \Omega} \frac{\delta_\ve(\phi)}{|\nabla\phi|} \pdv{\phi}{n} \delta\phi \dd S
- \int_\Omega \delta_\ve(\phi) \nabla \cdot \left( \frac{\nabla\phi}{|\nabla\phi|} \right) \delta\phi \dd x \dd y
- \int_\Omega \odv{\delta_\ve}{z}\nabla\phi \cdot \frac{\nabla\phi}{|\nabla\phi|} \delta\phi \dd x \dd y \\
& = \int_{\partial \Omega} \frac{\delta_\ve(\phi)}{|\nabla\phi|} \pdv{\phi}{n} \delta\phi \dd S
- \int_\Omega \delta_\ve(\phi) \nabla \cdot \left( \frac{\nabla\phi}{|\nabla\phi|} \right) \delta\phi \dd x \dd y
- \int_\Omega \odv{\delta_\ve}{z}|\nabla\phi| \delta\phi \dd x \dd y
\end{align}
なので、
\begin{align}
&\int_\Omega \left\{ \delta_\ve(\phi+\delta\phi)|\nabla(\phi+\delta\phi)|-\delta_\ve(\phi)|\nabla\phi| \right\}\dd x \dd y \\
&\simeq \int_{\partial \Omega} \frac{\delta_\ve(\phi)}{|\nabla\phi|} \pdv{\phi}{n} \delta\phi \dd S
- \int_\Omega \delta_\ve(\phi) \nabla \cdot \left( \frac{\nabla\phi}{|\nabla\phi|} \right) \delta\phi \dd x \dd y
\end{align}
である。$\delta\phi$の項がゼロとなることを要請すると、以下の式が得られる。
\begin{gather}
-\delta_\ve(\phi) \left[ \mu\nabla \cdot \left( \frac{\nabla\phi}{|\nabla\phi|} \right) -\nu -\lambda_1 |u_0 - c_1|^2 + \lambda_2 |u_0 - c_2|^2 \right] = 0 \:\: \mathrm{in} \:\: \Omega \\
\frac{\delta_\ve(\phi)}{|\nabla\phi|} \pdv{\phi}{n} = 0 \:\: \mathrm{on} \:\: \partial\Omega
\end{gather}
最急降下法的に時間発展方程式にすると、$t\geq 0$として
\begin{gather}
\pdv{\phi}{t}=\delta_\ve(\phi) \left[ \mu\nabla \cdot \left( \frac{\nabla\phi}{|\nabla\phi|} \right) -\nu -\lambda_1 |u_0 - c_1|^2 + \lambda_2 |u_0 - c_2|^2 \right] = 0 \:\: \mathrm{in} \:\: (0,\infty)\times \Omega \\
\frac{\delta_\ve(\phi)}{|\nabla\phi|} \pdv{\phi}{n} = 0 \:\: \mathrm{on} \:\: \partial\Omega
\end{gather}
適当に初期条件のレベルセット関数を決めて($\phi(0,x,y)=\phi_0(x,y) \in \Omega$)、これをうまいこと解いていくことで、評価関数$F$を極小にするような輪郭を得る手法がこの Chan-Vese レベルセット法である。
参考プログラムは、細かいところは多少違うが、おおむねこの方程式を解いて輪郭を抽出しています。レベルセット関数の再初期化とかもとりあえずやってます。
航空写真に適用
このセグメンテーションのプログラムも自分で書こうと思ったのですが、とりあえず一旦は前回と同じように参考文献のコードで計算を実行しました。
|
| 今回使った画像(google mapより) |
今回は、たくさんの障害物まわりの流れを解こうと思って、上のような画像をとってきました。画像は適当なソフトで傾きと大きさを調整しています。これに、以下のような初期条件を与えて、セグメンテーションを行ないました。
|
| 初期条件(左：ゼロレベルセット、右：レベルセット関数) |
結果は以下のようになりました。ローカルミニマムを求めるので左下の鳥の石は抽出できていませんし、影や歩行者の関係で石の輪郭をちゃんと抽出できていなかったりしていますが、とりあえずよしとします。鳥の石とかは初期条件をもう少し調節すれば抽出できると思います。
|
| 結果 (左：ゼロレベルセット、右：レベルセット関数) |
この周りの流れを計算したいと思います。
VP法 (Volume Penalization 法)
複雑境界の周りの流れを計算する手法として、前回は Direct forcing による埋め込み境界法を用いました。この手法では、境界を表現する点を選定してくる必要があったり、補間の時に隣接2点の距離が必要だったりと、ぼちぼち面倒です。
そこで、今回は VP法 (Volume Penalization 法, Angot et al., 1999 など) という手法を用います。これは、簡単に説明すると、固体部分をマスク関数(VOF関数のようなもの)で表して、固体部分では多孔質として流れを解く手法です。
運動量保存の式は、以下のようになります。
\pdv{\bm{u}}{t} + \bm{u}\cdot\nabla\bm{u}= -\nabla P + \nu \Delta \bm{u} - \frac{1}{\eta}\chi(\bm{x}) (\bm{u}-\bm{u}_O)
$\chi$はマスク関数で、固体中で1, それ以外で0とします。$\bm{u}_O$は固体の速度です。この最後の項が固体の影響を表します。
今、物体からの距離であるレベルセット関数を持っているので、こいつをうまいこと0~1の関数に変換します。今回は以下の関数で変換します。
\chi = \frac{1}{2} \left\{ 1+\tanh\left( \frac{\phi}{2\epsilon} \right) \right\}
ここで、$\Lambda<\chi<1-\Lambda$となる幅を$\Delta$とすると、$\epsilon$は
\epsilon = \frac{\Delta}{4} \frac{1}{\tanh^{-1}(1-2\Lambda)}
となるパラメータとしています。適当にパラメータを決めて変換したのが以下のマスク関数です。
|
| 今回使うマスク関数 |
VP法の浅水方程式への適用
目指しているのは川の流れを計算することなので、単純な2次元非圧縮よりも浅水方程式で解きたいですよね。そこで、VP法を浅水方程式に適用している研究を少し探してみました。すると、海洋モデルに適用した研究(Reckinger et al. 2012 など)が見つかりました。どうやらこれは圧縮性流体にVP法を用いた研究(Liu & Vasilyev, 2007 など)を浅水方程式に拡張しているようです。本記事ではこれらを参考にします。支配方程式は
\begin{align}
\pdv{h}{t} &= -\left[ 1 + \left(\frac{1}{\Phi}-1\right)\chi \right] \nabla\cdot(h\bm{u}) \\
\pdv{\bm{u}}{t} + \bm{u}\cdot\nabla\bm{u} &= -g\nabla h + \nu \Delta \bm{u} - \frac{1}{\eta}\chi (\bm{u}-0)
\end{align}
$\Phi$は空隙率、$\eta$はペナルティパラメータで、今回は$\Phi=10^{-2}, \eta=10^{-4}$としておきます。
河床勾配や河底の地形、底面摩擦などは一旦考えないようにします。渦粘性係数も今は定数($\nu=5\times10^{-3}$)としておきます。
数値計算
計算条件
- 今回こしらえた画像が、230×300ピクセルなので、格子点数も$N_x,N_y=230,300$としました。
- $\Delta x = \Delta y = 0.05\mathrm{m}$、すなわち$L_x=11.5\mathrm{m}, L_y=15.0\mathrm{m}$としました。
|
| 計算領域の模式図 |
- 境界条件
- 流入境界：$uh=Q$(定数), $v=0$, $\partial h/\partial x=0$
- 流出境界：$\partial(u,v,h)/\partial x=0$
- 壁境界(上下):$\partial h/\partial y=0$, slip境界
- 初期条件：$u=U_0, h=H_0$ (定数, $Q=U_0H_0$), $v=0$
- $\Delta t = 0.001$
- CFL条件：
- 移流：$\Delta t \leq \Delta x / U \simeq 0.01\mathrm{s}$
- 粘性：$\Delta t \leq \Delta x^2 / 2\nu \simeq 0.25 \mathrm{s}$
- 浅水波の伝搬：多孔質中の浅水波の位相速度 $c_{sw}=\sqrt{u^2+\frac{1}{\Phi}(gh-u^2)}\simeq\sqrt{gh/\Phi}$
より、$\Delta t \leq \Delta x / \sqrt{gh/\Phi} \simeq 0.05/\sqrt{10\times 1*100}\simeq 0.0016\mathrm{s}$
- CFL条件：
- 運動量保存式のペナルティ項は硬いため、$\exp(-\chi \Delta t/\eta)$ を掛けて時間発展
- 離散化：Arakawa C グリッド
- 差分スキーム：
- 移流：とりあえず1次風上差分
- 粘性・圧力勾配：2次中心差分
- 時間積分：SSP-RK3
結果
$Q=U_0H_0=0.4*0.25$として$t=150$まで計算しました。水位$h$の最大・最小値を見ると定常になってそうなので良しとしましょう。
|
| $h$の最大値と最小値の時間変化 |
計算結果のgifはこのようになりました。亀石を避けて流れ場ができています。また、流出までの距離が短いせいか、渦粘性や数値粘性のせいか、飛び石の配置のせいかわかりませんが、思いのほか流れは定常になっています。カルマン渦列的なのはあんまり見えません。
|
| 色：水位、矢印：流速 |
また、せっかく水位を解いているので、3Dプロットをしたくなります。亀石と水位を同時にプロットしたいのですが、pyplotだと両者の交わりがうまく描画されなかったり、plotlyだとアスペクト比を整えられなさそうなど、面倒な事がありました。これらの手法に別にこだわる必要はないので、今回はpythonで3次元プロットを行なうライブラリである mayavi を用いました。
水面を描いて、そこに色を塗る(たとえば流速の大きさ)スクリプトをとりあえず載せておきます。
import numpy as np
from mayavi import mlab
# 各パラメータを読み込む
xc = np.load("xc.npy") # セルセンターでのy座標
yc = np.load("yc.npy") # セルセンターでのy座標
Ts = np.load("Ts.npy") # 時間
vof = np.load("vof.npy") # vof関数(マスク関数)
xv, yv = np.meshgrid(xc, yc)
mlab.figure(bgcolor=(1,1,1),size=(800, 600))
# 描画する配列の配列を先に読み込んで用意
eta_data = [np.load(f"eta{t}.npy") for t in range(301)]
vel_data = [np.load(f"vel{t}.npy") for t in range(301)]
s=mlab.mesh(xv, yv, vel, scalars=vel,colormap="rainbow", opacity=0.6)
vmin, vmax = -30, 30 # カラーバーの最大、最小
s.module_manager.scalar_lut_manager.use_default_range = False
s.module_manager.scalar_lut_manager.data_range = (vmin, vmax)
v=mlab.mesh(xv, yv, vof*0.4, color=(0.8,0.8,0.8)) # マスク関数を適当に水位付近に
# 文字
title = mlab.title('t = '+str(round(Ts[0],1)), size=0.5, color=(0,0,0))
# カラーバー表示
cb = mlab.colorbar(s, title='Velocity', orientation='vertical')
cb.label_text_property.font_size = 12
cb.label_text_property.color = (0,0,0) # ラベルの色を黒に(これをしないと白い文字で書かれるので見えない)
cb.title_text_property.color = (0,0,0) # カラーバーのタイトルの色
@mlab.animate(delay=500)
def anim():
for t in range(301):
title.text = 't = '+str(round(Ts[t],1)) # タイトルを更新
s.mlab_source.set(z=eta_data[t], scalars=vel_data[t])
mlab.savefig(r"/kame_img/vort"+str(t).zfill(3)+".png")
yield
a=anim()
mlab.show()
このコードを適宜変えて動かしたらアニメーションが出てきます。連番画像を保存してffmpegでgifにしたものが以下になります。やはり三次元でのプロットは美しくてよいですね。
|
| 色：流速 |
|
| 色：渦度 |
|
| 色：水位 |
考察
水位を見ると、亀石の前で37cmくらい、後ろで10cm弱程度になっています。
|
| y方向平均水位 |
計算結果を眺めてるだけで満足してはまだ不十分かもしれません。計算結果はあくまでも方程式の数値解ですから、大事なのはそれが現実の現象をどれだけ再現できているかということです。これは川の流れという身近な現象ですから、実際に測れるものです。
ということで、鴨川にきました。この日はすぐそこの出町桝形商店街で七夕まつりをやってたせいか人がたくさんいます。猛暑の晴れの日ということで、流量も落ち着いていますし、涼を求めて鴨川に吸い込まれたのでしょうか。トンビも虎視眈々と食料を狙いながら舞っています。
|
| 夏の鴨川と夕涼み |
川に入ってこそこそと水深を測ってきました。人目をはばかりながら測ったところ、水深は、
- 上流側：15cm~30cm強
- 下流側：5cm~15cm弱
- 石の裏側は特に浅かった(砂が堆積)
といった印象でした。流速は、水面の気泡の速度を動画で撮って追跡したら上流側で0.45m/sくらいでした。数値実験の初期条件・境界条件はこの値を参考にしています。数値計算結果を見てみると、上流側で現実より深くなっています。下流側ではそこまで差はないですが、数値解のほうが気持ち浅めになっている気がします。
あと、当たり前なのですが、河底は平らなわけがなく、自然のもの(砂の堆積)や人工物でめちゃくちゃデコボコしています。このあたりも本当はモデルに組み込む必要があるかもしれません。
結果としましては、大まかには再現できているがまだ不十分(実験設定？計算条件？)といった印象です。いろいろ改善していくところがありそうです。
改善のアイデアとしては、実験設定の観点では
- 摩擦項 (底面、側面)
- 河床地形(傾斜+でこぼこ)の項
- 流量の見積もりの改善
などが考えられます。計算条件では
- 移流スキームの改善(今は1次風上)
- 渦粘性係数(今は定数)
- VP法のパラメータ
- レベルセット法で抽出した輪郭をさらに高精度に
など、まだまだ改善・調査できる部分が多くあります。今後の課題です。
おわりに
今回はレベルセット法でセグメンテーションした境界を用いて物体周りの浅水方程式を解きました。
空撮→レベルセット関数→VP法で計算という流れは割と便利そうだと思います。空から撮影できればいろいろ応用もできるかもしれません。暇なときに遊んでいきます。
これは単なる趣味の数値計算ですので、詰めが甘い部分もあると思います。何か間違い等ございましたらお伝えください。また、改善・応用のアイデアとかがあれば教えてくれると嬉しいです。
参考文献
- Chan-Vese level-sets in Python: https://github.com/kevin-keraudren/chanvese
- Chan, T. F., & Vese, L. A. (2001). Active contours without edges. IEEE Transactions on image processing, 10(2), 266-277.
- Angot, P., Bruneau, C. H., & Fabrie, P. (1999). A penalization method to take into account obstacles in incompressible viscous flows. Numerische Mathematik, 81(4), 497-520.
- Reckinger, S. M., Vasilyev, O. V., & Fox-Kemper, B. (2012). Adaptive volume penalization for ocean modeling. Ocean Dynamics, 62, 1201-1215.
- Liu, Q., & Vasilyev, O. V. (2007). A Brinkman penalization method for compressible flows in complex geometries. Journal of Computational Physics, 227(2), 946-966.