画像処理&機械学習で水の世界を鮮やかに！~Sea-thru法~

皆様こんにちは．18のhiraと申します．普段は水中ロボットを作るアクア研で活動し，情報通信系の平均点を底から支えています（引きずり下ろしているともいう）．

さて，海で泳いだことがある人も多いかと思いますが，水の中はどうも見えにくいですよね．特に「青緑がかっていて暗い」と感じるのではないでしょうか．そして肉眼だけでなく，カメラで撮った写真も同様に色味の無く暗い写真になってしまいます．

この画像の様に水中においては，それも遠ければ遠いほど急速に青緑がかり暗く写ってしまいます．せっかく極彩色で美しいサンゴ礁を撮っても，見直したくなる画像とはならないでしょう．

せっかくダイビングをしたというのに真っ暗な写真ばっかりでは面白くありません．理系が恥も外聞も耐え忍びキラキラした海に繰り出すというのなら，桃色のサンゴ，真っ黄色の魚，紫の貝殻などなど，ただ海を楽しみに来た至極真っ当な方々の目では見ることのできない総天然色の「本物の色」を見てみたいですよね？

では，なぜ海の中では物が「青緑がかっていて暗く」見えてしまうのでしょうか．実はこの眼（カメラ）に見えている暗い画像は，本来あるべき鮮明な画像が光学的に減衰し，水自体の青色が加わったものです．

水中画像は以下のような数式で表せます．

I_c = D_c + B_c

ここで I_c はカメラ(or眼．カメラも眼もRGB3チャンネルを用いるのでどっちでもいい．以下カメラで統一)見えている暗い画像の色（私たちが普通に今まで見てきた水中の世界）， D_c は本来の鮮やかな画像の色が水の影響で減衰した色， B_c は物体とカメラ間に挟まれた水自体の色です．添え字の ~_c は c=R,G,B の色の3原色それぞれを意味します．この I_c は最大値1のRGB値で，R=G=B=0なら黒，R=G=B=1なら白，R=1,G=B=0なら真っ赤となります．

もっと複雑には次のようになります．

I_c = J_c e^{-\beta_c^D(\vec{v}_D)z} + B_c^\infty(1-e^{-\beta_c^B(\vec{v}_B)z})\quad... 1

右辺第1項が D_c に，第2項が B_c に対応しています． J_c が求める鮮やかな元の「本物の色」です． z がカメラから物体までの距離です．第1項より \beta_c^D(\vec{v}_D)>0 とすると鮮明な画像( J_c )が距離 z の指数関数により減少していることが分かります．これにより元の鮮明な画像が高いRGBの値を持っていても，距離が遠いと指数関数の影響により減少し，暗くなってしまいます．そして第2項は水の色を表します．先ほどとは反対に物体との距離が遠いほど B_c という値が顕著になります．この B_c こそが物体とカメラの間の空間にある水そのものの色です．間に水がたくさんあるほど，水の色が目立ってくるのだと分かります．

ここから鮮明な画像 J_c を求めるにはどうすればよいでしょうか．この方法がDerya Akkaynak博士とTali Treibitz博士が開発した”Sea-thru法”です（参考文献参照）．手順としてはまず各ピクセルの B_c を求め， I_c 画像から除き，画像から水の色を除きます．そして D_c となった画像に対し，また各ピクセルごとの指数関数の逆数 e^{\beta_c^D(\vec{v}_D)z} を求め，これを掛けることで本来の色 J_c を求めます．下の画像が博士らが実際に処理を施した画像のビフォーアフターです．左側が撮影した画像 I_c ，右側が本来の鮮明さを取り戻した J_c に相当します．

この方法では，鮮明化したい暗い画像に加え，その画像における画面上の物体までの距離をピクセルごとに求めた深度マップを用います．あるピクセルに対し，その色と距離から真の色を求めていきます．「深度」とありますが潜水深度ではなく，カメラから被写体までの水平方向距離のことを指します．そしてこのSea-thru法に特徴的なことは鮮明化したい暗い画像とその深度マップ以外には画像や諸情報（水深や何かしらのセンサーの値），人による手動の調整が要らないことです．画像と深度マップを与えるだけで完結します（有れば有ればでより正確にできる）．

では実際にSea-thru法をやってみましょう．ここで使用する画像はSea-thru法の発表時に公開（論文中に保存先URLが記載）された画像と深度マップを使用します．この中から下の4910というナンバリングの画像を用います．

この時深度としてはカメラから0.9m~3.2m程度の範囲です．下の図が深度マップを図示したもので，手前側の一番近い青色が0.9m，奥の緑っぽい所が3.2m程度，黄色が遠すぎて無視する深度無限大の点です．

まず B_c こと水そのものの色を求め，除いてみます．しかし与えられたのは画像と深度マップだけです．海や湖ごとに大きく異なるであろう澄み具合や水の微妙な色の情報は与えられていません． I_c には D_c も含まれているのに，どのように B_c だけを見抜くかが問題です．

ここで使うのが「下位1%の暗いピクセル」です．この暗いピクセルは主に画像の特に暗い点や影にあたります．RGBそれぞれではなく、一旦画像を白黒グレースケールとし、グレースケール中で最も暗い点を選びます。ここではそもそも色情報が存在しなかったといえるので， D_c\approx 0 とみなせ， I_c \approx B_c と言えます．すなわち水の色そのものの色値がのっているだけのピクセルと考えられます．そしてこの下位1%の暗いピクセルを， z の範囲を10等分したそれぞれの距離の階層の中で求めます．例えば距離が0m~10mの幅がある画像あれば，0m~1mの範囲,1m~2mの範囲…,9m~10mの範囲の各範囲における暗い点を求め，距離のばらつきを持たせた暗い点( I_c \approx B_c の点)を求めます．求まった点の色の値，距離の情報から B_c の推定値 \hat{B_c} を以下の式から求めます．

\hat{B_c} = B_c^\infty(1-e^{-\beta_c^Bz})+J'_ce^{-\beta_c^{D'}z}\quad...2

これは B_c の式とよく似ています（正しくは I_c=B_c+D_c 似ているが， D_c っぽい第2項は暗いピクセルでも残った D_c の残余）．上で求めた暗い点を各 ~_c=R,G,B ごとに z と色の値 B_c をプロットすると2式のような形になります．このプロットに近似できる B_c^\infty ， \beta_c^B ， J'_c ， \beta_c を求めます．それぞれ[0,1]，[0,5]，[0,1]，[0,5]の範囲に収めます．MATLABならfit関数とfittype関数を用いることで与えた点プロットに対応する，求める形の数式（この場合2式の形）の B_c^\infty 等各係数を簡単に求められます．この大量の点をフィッティングできる1つの最適な曲線を探すのが機械学習です．あっさりです．最新技術なんてそんなもんです．RGB3チャンネルそれぞれの係数の組み合わせを求めます．c=G(緑色)に関するグレースケール下位1%のピクセルの色値と近似曲線は次のようになりました．

c=R,B も同様に求め，これにより2式からある距離 z における水の色 B_c を求めることが可能になりました．これと z の深度マップを組み合わせることで，各ピクセルごとの水の色の推定値が分かります．推定される水の色を抽出した B_c 画像は以下のようになりました．

下位1%の暗いピクセルを分析することで画像における距離ごとの水の色の振る舞いを推定できるようになりました．全ピクセル分求め，これを I_c 画像から引くことで1全ピクセルの D_c が求まります．元画像 I_c から上の水の色画像 B_c を除いた D_c だけの画像は以下のようになりました．ちなみに背景が真っ黒の場所がありますが，これは距離が求められない(遠すぎて計測不可)の場所です．これ以降の計算では無視します．

次はこの D_c に対し e^{\beta_c^D(\vec{v}_D)z} を求めます。そのためにはまずlocal illuminant mapというものを求めます。

local illuminantr mapは(も)よく分かりませんが、画像をめっちゃボヤけさせて浮かび上がる光の濃淡のマッピングみたいです。画像をボヤけさせると言いましたが、方法としては隣り合ったピクセルの平均を取り、それを何回も繰り返しどんどん画像を平滑化させていきます．local illuminant mapは次の処理で作成します．

\begin{aligned} a_c'(x,y) &= \frac{1}{N_e}\sum_{N_e}a_c(x',y')\quad\mathrm{but}\quad N_e(x',y')=(x',y')\mathrm{with}||z(x,y)-z(x',y')||\leq\epsilon\\a(x,y)&=D_c(x,y)p+a_c'(x,y)(1-p)\end{aligned}

a_c(x,y) は今までの I_c のようなピクセル (x,y) における色値です．1行目ではピクセル (x,y) に隣接する上下左右4つのピクセル (x',y') の値 a(x',y') の平均値を新たな (x,y) の仮値 a_c'(x,y) としています．この時，隣接する4つすべての点を平均するとは限らず， (x,y) との深度距離の差が \epsilon 以下の物のみを選びます．つまり物体が途切れ，段差になっている境界同士は平均化しません．2行目でこの a_c'(x,y) の値に係数 1-p を掛けたものと， D_c の (x,y) に係数pを掛けたものの和を新しい a_c(x,y) とし，再び隣接するピクセルとの平均処理を繰り返します．この図においては真ん中のサンゴの塊と周辺の水底の部分は大きな段差となっているので，この境界ををまたいで平均化することはありません．

やってみたところ，この処理はトンでもなく時間がかかりました．google colaboratoyを使っても200万回ほど繰り返し200時間くらいかかりました(合ってる気がしない…?)．値としては p=1/100000~1/1000000 ぐらいとめっちゃ小さくするとボヤけやすいみたいです．また \epsilon にあたる距離制限は0.1(10cm)程度としました．

修正もしながら遠隔授業のバックグラウンドで働かせること1ヵ月，出来上がったボヤけ画像のlocal illuminat mapがこれです．

それでは \beta_c^D(z) を求めます． \beta_c^D(z) は

\beta_c^D(z) = a e^{bz}+c e^{dz}

とあらわされ， a ， b ， c ， d はそれぞれ下限 [0,-\infty,0,-\infty] ，上限 [\infty,0,\infty,0] で与えられます．ここで変数 f=2 をもちいて \hat{E_c}=fa_c とすると \beta_c^D(z) は以下の関係を持ちます．

\hat{z}=-\frac{\log(\hat{E_c})}{B_c^D(z)}

求める \hat{B_c^D(z)} は上式において

\min||z-\hat{z}||

をもたらすものです．ここでは上の定義を少し改めます．あるピクセルに対して深度マップの値 z と a_c(x,y) が判明しているので，

\hat{\beta_c^D(z)}=-\frac{\log(\hat{E_c})}{z}\quad...3

を満たす \hat{\beta_c^D(z)} = a e^{bz}+c e^{dz} を求めます．求め方は \hat{\beta_c} 同様MATLABのfit関数を用いれば簡単です．あっさり機械学習です．ただ，この a_c(x,y) ですが図のように必ずしも深度 z で決まった値をとるわけではなく，ムラがあります．特にこの画像では同じ程度の z の場所でも中央のサンゴでは緑が強く，周辺の水底では青が強いです．理想的にはもっとボカしを繰り返して a_c が可能な限り均一になるべきなのかもしれませんが，繰り返すのも時間がかかるので上図の a_c を使いました．緑部と青部の混在してもフィッティングする線が変に歪んでしまう気もするので，画像両端の水底の青部だけを抜き出し，さらにある距離 z での平均の a_c を求めこの平均値を使います．両端の水底青色領域での a_c (平均)と距離の関係は次のグラフのようになります．

上から赤緑青の距離による強度変化です．カメラに写っている最小深度距離が1mなのでそこを0としてmm単位の距離でプロットされています．緑と青が怖いくらい似ていますが大丈夫なのか…?理想としては赤色みたく距離に比例して小さくなってくれることですが，緑青ではサンゴの足元から徐々ににじみ出た互いの色で盛り上がりができています．

構わず続けて3式の値でフィッティングすると \hat{\beta_c^D(z)} の値は次のようになりました．

うーん，こちらも理想的には赤色みたくスロープ状に小さくなるのが良いんですが…ともあれ求まってしまったので，これを用いて計算していきます．ちなみにこの \hat{\beta_c^D(z)} を用いて \hat{z} を求めてみたのが次のプロットです．

横軸が深度マップの(正しい)深度 z ，縦軸がフィッティングを用いて得られた \hat{z} です．理想的には y=x の関係になります．赤色は良さそうですが緑と青がまだ上手く行っていなかったようです．

難点はありそうですが，いよいよ目標である J_c を求めます． J_c は以下の式から求めます．

J_c = D_ce^{\hat{\beta_c^D(z)}z}\quad...4

D_c 全てのピクセルにこれまで求めた \hat{\beta_c^D(z)} 用いて4式の指数関数を掛けた結果，次の画像となりました．

水底は石っぽい灰色で，サンゴっぽい何かは赤みを帯びています．ただ，画像右上は水色となっています．これは右上の D_c に赤色成分が少なく，緑と青が残っているからです．

上の図は D_c の画像から1枚目が赤色成分のみを，2枚目が緑色成分のみを取り出したものです．ただでさえ赤色は減衰し易いのですが水色となった右上では特に少ないです．そのため赤色の e^{\hat{\beta_c^D(z)}z} を掛けてもほとんど0のままで，緑，青の成分のみが残り水色となります．水色となった理由としては D_c=I_c-B_c より D_c を求めた際，僕は0未満になってしまった値を0としていました．ただ，これだと「わずかでも色味があった」という事実を全く無視してしまうことになります．そのためいくら指数関数を掛けようと0のままで，当然緑と青しか残りません．0とするのではなく，「極小だけど残っている」として元の値（どうせ小さい）の何割かの微小値として残しておくと良かったかもしれません．

さて，いろいろ改善点などありますが，元の青緑に比べればそれっぽい結果にはなったのではと思います．ただ，この方法ではまだまだ時間がかかります．本記事で使用した深度マップは中央のサンゴの写真を周囲から何十枚も撮り，それらを画像処理で解析して深度距離を算出したのもです．もちろんリアルタイムに得られるマップではないので，これに代わる深度マップの取得方法が必要です．また a_c の作成には膨大な数の平均化処理が必要です．まさか本記事みたいに1フレーム200時間かけるわけにはいかないので，高速な平均化方法が望まれます．幸い博士らは現在リアルタイムでこのSet-thru法を行う方法を開発中だそうです．そう遠くない未来には画像どころか動画で美しい真の色で泳ぐ魚を見られるかもしれません．

僕の次の目標としては画像処理とは打って変わって，FPGAを用いたソナーをつくろうと思っています．やっぱり水中といえばソナー．見えない世界を見るのがソナー．海底マッピングや自動制御など水中ロボットの醍醐味はソナーから．年末年始ははんだの熱で暖まり，HDLプログラミング御事始めとなるかもしれません．目先の成績は視界0ですが．

 

参考文献，資料

・Derya Akkaynak Tali Treibitz「Sea-thru: A Method For Removing Water From Underwater Images」

・Derya Akkaynak「Sea-thru: A Method For Removing Water From Underwater Images @Israel Computer Vision Day 2019」