日本ロボット学会誌 Vol. xx

解No. xx, pp.1〜6, 200x1

説ディジタルフロー解析に基づく画像解析に関する研究

本多健二∗1
東京海洋大学 ∗1

1. は

じ め に

近年，生体イメージング技術の進歩により，細胞レベルでの様々な生体現象の観測が可能となってきている．これ
により, 顕微鏡による細胞の観察から, 計算機によって細胞画像の定量的データを得る, バイオイメージインフォマティ
クスと呼ばれる研究分野が盛んである. 細胞数や細胞核のD ཎ⏬ീサイズ，細胞核の形態的な特徴といった解析もその一つでE ኱ὠࡢᡭἲ

図 1 細胞画像あり，従来の観測者の目視による定性的な計測を計算機が行うことにより，定量的な計測値に基づく，客観的な判断メータの選定など試行錯誤的な手続きが必要となる [4] [5]．

基準を確立することが期待されている．
細胞画像を解析するためには，まず対象となる細胞を背このように閾値による分類では，画像の統計的性質によっ景から切り出す必要がある．画像二値化は対象の切り出してうまく対象を切り出すことが困難な場合があることが分に使われるもっとも基本的な処理の一つであり，二値化のかる．

精度がその後の結果に大きな影響を与えることが多い．こ

一方，画像の統計的性質によらず安定した画像二値化をれまで画像二値化は様々な手法が提案されているが，従来，行う方法の一つとして，幾何学的な構造から画像を捉える

画像の二値化では，閾値によって画素を分類する手法が多方法が考えられる．画像の幾何学的な解析は古くから提案く用いられている．Sezgin らは，これまで提案されてきたされている. 画像パターン上の点 (x, y) における濃淡の強閾値に基づく二値化手法を，その生成アルゴリズムに基づ度を f (x, y) と表すと，画像パターンは 2 変数関数 f (x, y)いて分類し，細胞画像を含む，NDT(NonDestructiveによって構成される曲面とみなすことができる．HaralickTesting) 画像，文書画像に対して各手法の比較を行っている [1]．例えば，図 1(a) に示す細胞画像 [2] に対して，細胞らは，ヘシアンの固有値の符号に基づき曲面の凹凸構造を内の核を切り出すことを考えてみる．画像二値化の代表的4 種類のクラスに分類し，各画素を記述する方法を提案しており [6]，この分類に基づいて画像の二値化を行うことがな手法である大津の手法は，クラス間の分散とクラス内の考えられる．しかし，従来の計算機において，幾何学的な分散との比によって分離度を定義し，分離度が最大となる

解析は時間のかかる処理であり，画像二値化への応用は検閾値を求めることにより，自動的に二値化を行う手法であ討されていない. 近年の計算機処理能力の飛躍的な向上にる [3]．しかし，大津の手法は濃淡値のヒストグラムが双峰

よってこれまで困難とされていた処理が容易になり，幾何性の分布をなすことを前提とした手法であり，二値化の結学的な手法の可能性も再検討する価値がある.

果は画像の統計的な性質に依存する．図 1(a) に示す細胞画

そこで，本論文では画像の幾何学的な構造に着目し，曲像は，細胞が密に存在する部分と疎な部分が存在し，濃淡面の凹凸構造に基づいて細胞画像の二値化を行う手法を提値の統計的性質が一様ではない. そのため，図 1(b) に示す案する. Haralick らの方法による記述は画素がどのクラスように，不安定な結果となる．濃度値の分布を調べる範囲に属するかという局所的な情報のみが扱われており，そのを，画像全体からより局所的に調べることで，二値化の精凹凸構造が，より大きな幾何学的構造の中でどのような形度を向上させることができることもあるが，適切な結果を態をとっているかといった，大域的な視点での解析はなさ得るためには局所領域の大きさや統計的処理のためのパラれていない．提案手法では，任意の曲面が，その勾配に基づく流れの場をつくるという性質に基づき，曲面の大域的なキーワード：画像二値化，ディジタルフロー，細胞画像

日本ロボット学会誌 xx 巻 xx 号

幾何学的構造をよく表すことのできる流線を定義する．こ—1—200x 年 xx 月

本多健二

2

図6
図 2 画素 x 0 とその 8 近傍

最少被覆フロー集合

ᴟ኱್ 凸凸 x1 < x2x1

凸 x2 < x3x2x3

凹 x3 > x4x4凹 ᴟᑠ್

図 7 双極フローの凹凸構造をもとめる．u = f (x, y) とすると u = f (x, y) の 1 階全

図 3 4x4 画像

微分

du =

∂u
∂u
dx +dy∂x∂y（1 ）∂uにおいて，∂u∂x , ∂y が微分可能ならば，x, y に関する du の全微分として d2 u が得られる．

d2 u = d(du) =図 4 被覆フロー集合∂2u∂ 2u 2∂2u
dxdy + 2 dy 2 （2）dx +22∂x∂x∂y∂y

ここで，u の x に関する微分を fx ，fx の y に関する微分を fxy ，u の y に関する微分を fy ，勾配とヘシアンを 
fxfxx∇f =H=fyfxyfxyfyy

とすると，式 (1) は，fxx fx2 + 2fxy fx fy + fyy fy2 = (∇f, H∇f )（3 ）となる．
本論文では，流線の変曲点がなす幾何学的構造，すなわ
ち (∇f, H∇f )=0 となる画素の集合を，(∇f, H ∇f ) 零交
図 5 DF 集合

差線と呼ぶ．

2. 2 ディジタルフローの流線は，微分幾何学的性質より，極大点から流れ出て極小点へと流れ込む．そして，流線上の凹凸構造を解析する

ディジタル画像において流線を定義する．以下では，そのための準備を行う．

ことによって細胞画像の二値化を行う. 更に，切り出され〔定義 1〕(近傍) 画像中の任意の画素 x 0 に関して，図 2 にた結果を従来の二値化手法と比較し，幾何学的構造にもとおける要素の集合 {x 1 ,x 2 ,· · · x 8 } を 8 近傍，{x 1 x 3 ,x 5 ,x 7 }ずく提案手法の有効性を検証する.

を 4 近傍と呼ぶ．ここで添字の集合 {1,2,· · · ,8} を S ，2. 準{1,3,5,7} を S1 で表す．
〔定義 2〕(極大点，極小点) 画素 x 0 = (x, y) における濃度値を f (x 0 ) とする．画素 x 0 が次の条件を満たすとき，x 0 を極大点と呼ぶ．

備

2. 1 (∇f, H ∇f ) 零交差線流線の凹凸がなす幾何学的構造について調べる．
画素 (x, y) の濃度値を f (x, y) とし，f (x, y) の 2 階微分

JRSJ Vol. xx No. xx—2—xx, 200x日本ロボット学会誌 LATEX 2ε クラスファイルの使い方

͐ᴟᑠⅬ

͐ᴟ኱Ⅼ

図 9 領域スクリーニング

図 8 細胞画像の凹凸構造3

図 10 スクリーニング画像と DF からなる領域

3. ディジタルフロー解析に基づいた
画像二値化
㞽஺ᕪ⥺

3. 1 双極フローの凹凸構造
DF は曲面の大域的な幾何学的構造をよく表すことができる．ここでは，双極フローの凹凸構造を解析するため，以㞽஺ᕪ⥺

下の定義を行う．
D 㡿ᇦ$ࡢቃ⏺ࣃࢱ࣮ࣥ͐ᴟ኱್

〔定義 5〕(被覆フロー集合，最小被覆フロー集合) DF のE 㡿ᇦ%ࡢቃ⏺ࣃࢱ࣮ࣥ͐ᴟᑠ್

集合であって，画像中の全ての画素が含まれているものを，

図 11 境界のパターン全ての k(k ∈ S) に対して，f (x k ) ≤ f (x 0 )被覆フロー集合という．特に，被覆フロー集合の中で，元の数が最小のものを最小被覆フロー集合と呼ぶ．
（4 ）各画素について DF が存在するため，図 3 に示す画像のまた，画素 x 0 が次の条件を満たすとき，x 0 を極小点と場合，全部で 16 本の DF が得られる．ここで，同じ系列呼ぶ．の DF は 1 つの DF とみなす．図 4 に図 3 の画像における被覆フロー集合を示す．
また，以下の定義も行う．

全ての k(k ∈ S) に対して，f (x k ) ≤ f (x 0 )〔定義 6〕(DF 集合) 画素 x 0 に流れる DF の集合を，画（5 ）〔定義 3〕(最大勾配画素，最小勾配画素) 画素 x 0 と近傍画素 {x k :k ∈ S} との濃度値の傾きを勾配と呼び，近傍画素 {x k :k ∈ S} の内，最大勾配を与える近傍画素 x k をx 0 の最大勾配画素と呼ぶ．
また，近傍画素 {x k :k ∈ S} の内，最小勾配を与える近傍画素 x k を x 0 の最小勾配画素と呼ぶ．
〔定義 4〕(双極フロー) 各画素 x i において，最大勾配画素をたどって極大点 a に至る系列と，最小勾配画素をたどって極小点 b に至る系列を合わせた点の系列

素 x 0 の DF 集合と呼ぶ．
図 5 に図 3 の画像における DF 集合を示す．被覆フロー集合と各画素の DF 集合を用いて，以下の手順に従い最小被覆フローが得られる．

(step1) 各画素の DF 集合を参照し，DF が 1 本しか通っていない画素が存在した場合，その DF は選択される．

(step2) 被覆フロー集合を参照し，(step1) において選択された DF が通らない画素を通る DF がある場合は，そのを双極フローの長さと呼ぶ．以下，双極フローを DF(DipoleDF を選択する．もしこの DF が 2 つ以上ある場合は，フローの長さ M の短い方を選択する．
図 6 に図 3 の画像における最小被覆フロー集合を示す．
最小被覆フロー集合の各 DF の凹凸を逐次判定することで画像の凹凸構造が得られる．DF の凹凸の判定は，DF を追跡する際に，図 7 に示すように，凸である極大点から，凹である極小点へと，DF に沿って，前後の濃度値の変化量Flow) とする．の増減を比較することで決定する．

x 1 (= a),x 2 , · · · x i , · · · x M (= b)（6 ）が存在する．この点の系列を双極フローと呼ぶ．また，M

日本ロボット学会誌 xx 巻 xx 号

—3—200x 年 xx 月

本多健二

411111111111111111

1
111111111222
2 22 222112

1111

111111111

1111111222111111111111111

11112222222222222222222 21111121122222222222

111111111112211111111111111111

2111222112

221111111122222 2112 212221212222 2212121222212222122222
122 222122 21222222122222
12222 2222
12222 222
22 2 22222
2
12222222 2212222222 2
122222211221221222222 11222222
2
1222
11222 21222
1
112
1212 22
11 2
1
122
1222
22111112 2
12 22
1
12 2 2
1
21 1
2 2 2
1
112 2 2
1
1 1
222
1
1 12
22
21
1
1
2 2 2
1
21 1 2
21
2
22
1
222
22
2 2
11
1122
22
2
2 2 21111
1111112
2
222
112 1
22
11
22 2 12211 1
2 21111111222
12 2
21211111
2222
2 2
12 2
12 2222 211111
222
11222122222
222222222222222
22
222 2
222222
222
22222222222221
1111111111111

11

1111111

111

111

11111

11111111

11111211

1122212222222

112122222222222

22

22

1111
11
1
11122222

211112

2222
2

22

222
2 2222222222222
2
22222222222
222
2

1
1111
11
1

112
2
1
1
2 2
2 22222222
222
22
2 2
222 2
222222222
2222222222222111111

111

111111111111111111111111111111111

112111111111111111111

111

111111221
111112222221

11121

1

11221222

11111111111

122222221111111122222222

22
12 22
122 22 221222122212122
222
1222222222
222
112 22222222211222211222
1
11122
112211222222
111112222111112 2222
111122211 11111 1222
1 11
2221112 22211221 121112
1 1122 22
1 112211112
1 11 2
121
22
22
1
22
111
1222
22111111 2
11122222111
2212222221111111111
2222 2
22
2
22222
2 2
2
22
112 22222
2112
22
2 2 2
22
22222
22222
2222222222222Ȫ 㸪U 

11111111111111222222
2211111111222222222222

111

11

12212222222

11111111122222222222111

1121112

222

2
22222
2

2222222222
2

2222222
211111111111
11221111112
2222
2 2
22 2
2 222222
22222
22222222
22222
2222222222

22222

Ȫ 㸪U 

11121112Ȫ 㸪U 

1
11111
1122222111222211111111111111111111111111111111111111111111111

11111111111111111111111111111111111111111111111111111111111111111111111

111111111111111111111

111

1111111111111
111111111112Ȫ 㸪U 図 12 各スケールにおける境界クラス

図 1 の細胞画像における，DF の凹凸構造の解析結果をための客観的な指標が必要となる．そこで，スケールを選択

図 8 に示す．図 8 に示すように，DF の凹凸構造によって，する指標として，領域スクリーニングによって得られた画像の凹凸構造の境界に着目する．この境界は，(∇f, H ∇f )画像の幾何学的構造を捉えていることが確認できる．

零交差線であり, 以下では，(∇f, H ∇f ) 零交差線が，安定

3. 2 領域スクリーニング
DF 解析による凹凸構造は，画像の幾何学的な特徴をよく捉えることができる．しかし，DF による原始的な凹凸構造は，核以外の微細な構造も捉えている．画像の幾何学的な構造をより正確に捉えるためには，曲面の大域的な凹凸構造と微細な凹凸構造を判別する必要がある．このために，凹凸構造の中に含まれる極値に着目する．凹凸構造に極値を重ね合わせた様子を図 8 に示す．凹凸構造に含まれる極値に着目すると，図 8 に示すように，凹凸構造の中に極値を含む連結領域と含まない連結領域に分けることができる．
そこで，極小値をもたない凹連結領域を凸連結領域へ，極大値を含まない凸連結領域を凹連結領域へ変換することで領域スクリーニングを行う．図 9 に，図 1 に示す細胞画像のスクリーニング後の画像を示す．スクリーニングにより，画像の大域的な構造が捉えられている．
次に，凹凸構造の境界線の構造に着目する．画像の凹凸構造の境界線は，単純で安定な構造と，入り組んで不安定な (∇f, H ∇f ) 零交差線と不安定な (∇f, H ∇f ) 零交差線に分類できることを示し，安定な (∇f, H ∇f ) 零交差線で囲まれるスケールを探すことによって対象を切り出す．
まず，微分幾何学的な性質により，ある極大点からある極小点へ流れ込む DF の集合は，一つの領域を形成する．
図 10 に，図 8 の各極大点から各極小点へ流れ込む DF がつくる領域を示す．図 10 の各領域を観察すると，各領域における (∇f, H ∇f ) 零交差線の構造は，図 11 に示す 2 つのパターンに分類される．図 11(a) と (b) はそれぞれ，図10 の斜線で示した領域 A と B の (∇f, H ∇f ) 零交差線の様子を示している．
図 11 に示す零交差線の数に着目すると，DF 上の零交差線は以下の 2 つのパターンに分類される．

(パターン 1) DF 上に，零交差線を 1 つ持つ場合．(図11(a))(パターン 2) DF 上に，零交差線を 2 つ以上含む場合．
(図 11(b))
そこで，零交差線を分類するため，以下を定義する．

な構造が混在していることが確認できる．DF 解析による曲〔定義 10〕(境界クラス) 境界クラス C f は，面の凹凸構造のような微分幾何学的な構造は，ガウシアンフィルタのスケールパラメータによって変化する性質があ⎧⎨1, (DF 上に零交差線を 1 つ持つとき．
)
Cf =⎩2, (DF 上に零交差線を 2 つ以上含むとき．
)

る．そして，曲面の凹凸構造の変化によって境界線の形態も変化する．そこで，ガウシアンフィルタのスケールパラメータを調整して，DF 解析による記述によって境界線の形（7 ）態変化を捉えながら適切なスケールパラメータを選択することで二値化を行う．以降，ガウシアンフィルタのスケールパラメータをスケールと呼ぶ．次項では，自動的に安定な構造が得られるような指標を決め，それによりスケールを調整する手法を検討する．

図 1 の細胞画像の画像中央の核に着目し，スケールを段階的に変化させて, クラス分けした結果を図 12 に示す．図3. 3 境界クラス

12 は，画像中央の核が全て安定境界で囲まれるまでのい適切なスケールを選択するためには，それを自動化する

JRSJ Vol. xx No. xx

以下，境界クラス 1 の境界を安定境界，境界クラス 2 の境界を不安定境界と呼ぶ．

くつかのスケールと，境界に対して安定境界が占める割合—4—xx, 200x

日本ロボット学会誌 LATEX 2ε クラスファイルの使い方

5⣽⬊⏬ീ ⏬ീ

.LWWOHUࡢᡭἲ.DSXUࡢᡭἲ

<HQࡢᡭἲ

*URXQGWUXWK኱ὠࡢᡭἲ1,EODFNࡢᡭἲᥦ᱌ᡭἲ

図 13 評価用データセットの例と各手法による二値化画像r を示す．図 12 に示すように，スケールが増加するに従い，核が安定領域で覆われてくる．図 1 の細胞画像の場合，σ = 2.2 において初めて対象領域の境界全てが安定な境界となる．

4. 実

験

分割されていれば Over-Segmentation(以下，OS)，2 つの細胞が 1 つに融合していれば Under-Segmentation(以下，US)，画像には存在する細胞が二値化結果に含まれていない
場合を Under-Ditection(以下，UD) とし，CS，OS，US，UD と判定された細胞数を求め，Ground Truth に示された全体の細胞数に対するそれぞれの比率を求めた．また，提案手法の有効性を示すため，従来の二値化手法との比本来画像にない細胞が二値化結果に存在した場合を Over-較を行う. 実験に用いたデータセットは，Broad BioimageDitection(以下，OD) として，OD と判定された細胞数と全体の細胞数における比率を求めた．その結果を図??に示す．図??より，Niblack 法，提案手法を除く手法では，OS，US がともに多く検出される．これは統計的性質が一様でない画像に対し，大域的に閾値を決定することによる．OD

Benchmark Collection datasets [2] を用いた．データセッ
トは Drosophila Kc167 cells に含まれる 4 種類のグレースケール画像である．一部のデータセットと，人間の判断による Ground truth を図 13 の 1 列目に示す．
比較対象には，Sezgin らの実験において高い精度の結果は，入力画像に細胞状のものが観測できるものの，実験にが得られた手法のうち，Kittler 法 [7]，Kapur 法 [8]，Yen用いた Ground Truth では切り出されていない場合に計測法 [9] の 3 手法と，画像二値化の代表的な手法として大津される．
また，切り出しの精度を確認するため，全体の細胞数に法，Niblack 法の 2 手法で計 5 種類の手法を用いた.
スケールに関して，スケールの初期値，増加分，上限を実対して CS と判定された細胞数の比率を各手法ごとに比較

験的に求めそれぞれ σ0 = 1.4，Δσ = 0.2，σn = 3.0 とし，した結果を図 14 に示す．図 14 に示すように，Niblack の切り出しの条件である安定境界が占める比率を 95 ％以上

手法及び，提案手法を除く各手法は，画像によって切り出とした．各手法との比較は同一のスケールで比較する．ス

しの精度にばらつきがある．局所的な領域を設定し解析すケールは実験的に求め，σ = 2.2 とした．

る手法は，対象のおおまかな構造をまず捉え，必要に応じ二値化手法は ImageJ [10] での実装である．一部のデータて一部分に注意を向け詳しく観察を行うという，人間のもセットについて，各手法による二値化結果を図 13 の 2 列目かのの見方に近い手法であると言える．実験において，局所ら 4 列目に示す. 提案手法の各画像 (512pixel×512pixel) に的な手法の Niblack の手法は，安定して細胞を切り出しておける 1 スケール当たりの平均処理時間は 1.22 秒であった．

いる．また，提案手法においても，局所的な二値化手法とまず，細胞がどのような状態で切り出されたかを確認す同等の結果が得られている．

るため，データセットに含まれる Ground Truth に対し，次に，領域の一致度を評価する. 評価は Sezgin らが領域二値化によって得られた細胞が，正しく切り出されていればの一致度として用いた指標に準じて行う. Sezgin らは領域Correctly-Segmentation(以下，CS)，1 つの細胞が 2 つにの一致度として RAE(Relative foreground area error) を日本ロボット学会誌 xx 巻 xx 号

—5—200x 年 xx 月

本多健二

6120

いて，AR2 が 1 に近いことから，Niblack の手法，提案手100法において Ground truth における領域をほぼ切り出している．AR1 より，Niblack 法において領域を大き目に切り80出す傾向があり，提案手法が Ground truth との一致度が⏬ീ1⏬ീ260高い．実験に用いた Ground Truth は人間の目による判断⏬ീ340⏬ീ4で対象と背景を切り分けているため，人間の心理的な影響も少なからず影響していると考えられるが，提案手法が人20間の判断に近い基準で切り出せているといえる．

0KapurKilerotsu

yenniblack

proposed

以上の結果より，提案手法が，細胞画像の二値化におい

図 14 細胞の計数結果の比較

Aoて有効に機能し，従来手法との比較においても同等もしくは高い精度で細胞を切り出せる二値化手法であるといえる．

AT

また，これは DF の変曲点による切り出しが精度の点で安定に機能することを示しており，提案手法の有効性を確認

Ao ∩ ATR1 =Ao ∩ ATATできた．

R2 =参

Ao ∩ ATAo

図 15 評価指標の定義
1.2

10.8

kapur
kiler

AR2 0.6

otsu
yen

0.4

niblack
proposed

0.2

000.20.40.6
AR10.8

図 16領域の一致度

11.2用いている. 本論文では，Sezgin らの方法に準じて，前述の実験において CS と判定されたそれぞれの細胞に対し式(8) に示すように評価する．
T
R1 = A0A∩ATT
R2 = A0A∩A0考文

献[ 1 ] M. Sezgin and B. Sankur:
Survey over Image Thresholding
Techniques and Quantitative Performance Evaluation , Journalof Electronic Imaging, 13(1) (2004).
[ 2 ] BBBC(Broad Bioimage Benchmark Collection):
http://www.broadinstitute.org/bbbc/.
[ 3 ] 大津展之，”判別および最小 2 乗規準に基づく自動しきい値選定法，”
信学論 (D)，Vol.63，No.4，pp.349-356，1980.
[ 4 ] W. Niblack, An Introduction to Digital Image Processing,
pp.115?116, 1986.
[ 5 ] J. Sauvola, M. Pietika, Adaptive document image binarization, Pattern Recognition 33(2), pp.225-236, 2000.
[ 6 ] M Haralick, LT Watson, T. J. Laﬀey, ”The Topographic PrimalSketch,” The International Journal of Robotics，2, 1, pp.195217, 1983.
[ 7 ] Kittler, J, Illingworth, J (1986), ”Minimum error thresholding”,
Pattern Recognition 19 41-47[ 8 ] Kapur, JN; Sahoo, PK, Wong, ACK (1985), ”A New Method forGray-Level Picture Thresholding Using the Entropy of the Histogram”, Graphical Models and Image Processing 29(3): 273285[ 9 ] Yen JC, Chang FJ, Chang S (1995), ”A New Criterion for Automatic Multilevel Thresholding”, IEEE Trans. on Image Processing 4 (3): 370-378[10] ImageJ: http://imagej.nih.gov/ij/.

（8 ）ここで，A0 は Ground truth，AT は二値化によって切り出された領域の面積を表す. 図 15 に各領域の様子を示す．各画像の R1，R2 の平均をそれぞれ AR1，AR2 とし，4 種類の細胞画像に対しての各手法ごとの AR1，AR2 を図16 に示す．図 16 に示すように，Niblack 法と提案手法を除く各手法は分散が大きい. これは，細胞数の計数結果と同様に，画像によって統計的性質が一様でないことが，切り出す形状にも影響していることによる. これに対し，局所的な手法である Niblack 法と提案手法の分散は小さく，同程度となっている. 一方，切り出される領域の大きさにつJRSJ Vol. xx No. xx—6—xx, 200x