& 3 3 ' ' (., (Pixel), (Light Intensity) (Random Variable). (Joint Probability). V., V = {,,, V }. i x i x = (x, x,, x V ) T. x i i (State Variable),

全文

(1)Vol.2009-CVIM-167 No.36 2009/6/10. 情報処理学会研究報告 IPSJ SIG Technical Report. 中に含まれているとは限らない場合に有効に機能することが期待されている. 不完全性を伴うデータを扱う情報処理において, その不完全性からくる不確実性を数式を用いて表現する. 大規模確率場と確率的画像処理の深化と展開. 伝統的数学の道具としての確率・統計をもちいるというわけである. 通常, あつかわれる自然画像, 風景画像は空間的に一様な性質を持っている場合は少なく, 従来の空間フィルター. 田中和之†1. だけでは対応しきれない場合がある. 信号処理においてはこの問題を解決する手段のひとつとして適応フィルターとよばれるものが提案されている. これに対して確率的画像処理の場合は, 生成モデルの性質は空間的に一様な統計的性質をもつものを仮定する. 生成される画. 確率的画像処理はマルコフ確率場とベイズ推定をもとにして定式化される. その具体的なアルゴリズムは確率伝搬法により構成される. マルコフ確率場と確率伝搬法は物理モデルおよび物理的計算技法と深く関係したアプローチである. 本稿ではこれらの理論的枠組みとその基礎となる情報統計力学的視点について解説する.. 像はランダムネスを伴うため, その生成モデルの統計的性質が強く反映されている部分と必ずしもそうとはいえない部分が共存する画像が生成される. このランダムネスをともなう統計的性質を系統的に制御する鍵がベイズの公式 (Bayes Formula) である. 確率的画像処理では画像生成モデルとベイズの公式から構成される確率モデルの統計量. Deeping and Expansion of Large-Scale Random Fields and Probabilistic Image Processing. を計算することが主な手順となる. しかし, この計算には通常は画像の大きさに対して指数関数的に計算量が増加してしまう計算の複雑さ (Compulational Complexity) の問題がついてまわる. つい最近まで, あまり有効な方法として認知されることはなかった. しかし, 確率伝搬法 (Belief Propagation) とよばれる方法が 1980 年代半ばに人工知能の分野. Kazuyuki Tanaka†1. で提案され, その後, 実用的な確率的画像処理のアルゴリズムが提案されるようになる. また, 画像生成モデル, 確率伝搬法の数理構造のなかに物理モデルと物理的計算技法との関係. The mathematical frameworks of probabilistic image processing are formulated by means of Markov random fields and Bayesian inference. Some practical algorithms are constructed by applying belief propagation methods. Markov random fields and belief propagation methods can be regarded as one of statistical-mechanical approaches to probabilistic information processing. In the present paper, we review some fundamental frameworks and statisticalmechanical approaches of probabilistic image processing.. が明らかになり, 情報統計力学 (Statistical-Mechanical Informatics) という物理学における新しい研究分野の主要な研究テーマのひとつとしても認知されつつある.1)–3) 本稿では, もっとも基本的な問題設定であるノイズ除去という用途において確率的画像処理についての説明を進めることにする. ノイズ除去はノイズの加えられてしまった画像からノイズの加えられる前のもとの画像を推定する問題である. これを信号処理の言葉になおせばノイズがデータ, すなわち観測信号であり, もとの画像が原信号ということになる. つまり原信号をデータから推定するという問題である. 確率的画像処理の場合はもとの画像. 1. はじめに. の統計的性質とランダムネスを伴うノイズを確率を用いてあらわし, ベイズの公式 (Bayes. 確率的画像処理とは画像の生成モデル (Genrrative Model) を確率をつかって仮定す. Formula) を用いて, データが与えられたときに, そのデータを生成した画像はもともとはど. ることを基本とする方法である. その用途はノイズ除去, 領域分割, エッジ検出などの不完全. のような画像であったかをあらわす確率を導き, 推定がおこなわれる. ベイズの公式をつかっ. データを伴う処理にあり, 要請される処理を進める上で必要とする情報のすべてがデータの. たこのようなアプローチはベイズ推定 (Bayesian Inference) などとよばれている. ベイズ推定においてもとの画像の統計的性質をあらわす確率を事前確率 (Prior Probability), データが与えられたときに, そのデータを生成した画像はもともとはどのような画像であっ. †1 東北大学大学院情報科学研究科 Graduate School of Information Sciences, Tohoku University. たかをあらわす確率を事後確率 (Posterior Probability) とよんでいる.. 1. c 2009 Information Processing Society of Japan °.

(2) Vol.2009-CVIM-167 No.36 2009/6/10. 情報処理学会研究報告 IPSJ SIG Technical Report. により与えられるものとする.. 2. 画像生成モデルと大規模確率場. V と E で定義されるグラフ (V, E) 上で画像 x =. (x1 , x2 , · · ·, x12 )T が与えられる. 図 1 はこのグラフ (V, E) 上で定義される確率 Pr{X = x}. 画像処理において確率を導入する場合に, それぞれの画素 (Pixel) での光の強度, すなわち. のグラフ表現である.. Lmn?o pqb & 21 &/3 5 4 1 3 '64 ' 1 ' 3 ( 4 (

(3) 12&/3 ) 4 1 3 * 4 1 ' 3 + 4 1 ( 3 , 4 ) 1 ) 3* * 1 * 3+ + 1 + 3, , 1 ) 3E- 4 4 1 * 3F&G. 44 1 + 3F&?& 4 4 1 , 3H& 54 - 798 :<;%=?> &/. 7@;A= :<;B;%>C&0& 7@;B;A:<;BD0> & # $ !!" #%$ IKJMLON5PFQRP<STPVURPFWXPZYOPZ[\PH]\PF^XP_N6`\P_NaN5PN5Qb cdJeLXf N5PFQghP@fiQRP<SjghP!fkSTPVUghP!fiWXPHYRg P@fBY\PZ[Xg P@fB[\PZ]Xg P@fB^\PN6`Xg P@f Nl`\PNmNOghP@f N6`OPNmNOghP f NmNlP_N5QghP@f NlPFWTghP!f2QRPHYRghP!fkSTPH[Rg P@f2URPZ]Xg P@fBW\PH^TghP@fAYOPN6`Xg P@fA[OPNmNOghP!fB]\P_N5Qg0b. 輝度値 (Light Intensity) に対して確率変数 (Random Variable) に対応させることになる. つまり画素の総数と同じ次元をもつ結合確率 (Joint Probability) をつかってあらわされることになる. |V | 個の画素が正方格子上にならべられている場合を想定する. そのひとつひとつの画素に通し番号をつけて, すべての画素からなる集合を V = {1, 2, · · ·, |V |} とあらわすことにする. i 番目の画素の輝度値を xi とするとひとつの画像は x = (x1 , x2 , · · ·, x|V | )T. という縦ベクトルであらわされることになる. このとき xi は画素 i における光の強度の状態変数 (State Variable), 縦ベクトル x は画像に対する状態ベクトル (State Vector) とよぶことができる. 各画素 i の輝度値に対する確率変数 (Random Variable) を Xi と定義すると, この画像全体には X = (X1 , X2 , · · ·, X|V | )T という確率変数のベクトル, すなわ. 図 1 もとの画像 x に対する式 (1) の確率 Pr{X = x} のグラフ表現である.. ち確率ベクトル (Random Vector) がわりあてられることになる. つまり画像 x の確率. は Pr{X = x} によりあらわされる. そこでユーザーがノイズ除去において最終的に得たい画像の性質を確率 Pr{X = x} を. この確率は何をあらわすかを 2 値画像 (Q = 2) で説明しよう. 式 (1) の確率分布がどのよ. つかってどのように表現するかについて考えてみよう. たとえば一枚の紙の上に黒い点がひ. うな画像を生成するかについて考えてみる. まずもっとも簡単な = 1, 2 , E = {1, Ã V! Ã 2}! 0 x1 を並べ上げると , によるグラフ (V, E) でまず考える. 可能なすべての状態 x = 0 x2. ©. とつのっていたらこれはゴミとみなして手で払ってのけようとするであろう. ノイズ除去はこの操作をアルゴリズム化することが基本となる. このことは「局所的にあるひとつの画素. Ã ! Ã ! Ã !. に着目し, その近傍画素と階調値ができるだけ似た値になる画像がノイズが加えられる以前. 1. の画像である」と仮定した操作とみなすことができる. ここまでは空間フィルターの考え方. 0. と同じである. しかし, 確率をつかう場合は上記の操作を選択する場合と選択しない場合のを用いて, この操作を実現するひとつの選択肢は次の表式である.. ª. Pr X = x ∝. Y. {i,j}∈E. 0 1. Pr{X =. 2 種類が自動的・適応的に可能であるようにデザインすることができる. 確率 Pr{X = x}. ©. ,. 1. ,. Ã 1! 1 1. ª. ©. ª. の 4 つである. この 4 つの状態の確率の大小関係は. } = Pr{X =. Ã ! 0 0. } > Pr{X =. Ã ! 1. } = Pr{X =. 0. Ã ! 0 1. }. (3). と与えられる (図 2 参照). つまり隣の画素の状態が決まれば自分の画素の状態が決まること. ¡. ¢ 1 exp − α(xi − xj )2 2. になる.. (1). . . . . . E はすべての最近接画素対 {i, j} の集合である. 画素の集合が V = 1, 2, · · ·, 12 が縦 3 行, . 横 4 列の正方格子の上に並べられている場合を考えよう. 最近接画素対の集合は. ©. E = {1, 2}, {2, 3}, {3, 4}, {5, 6}, {6, 7}, {7, 8}, {9, 10}, {10, 11}, {11, 12},. ª. {1, 5}, {2, 6}, {3, 7}, {4, 8}, {5, 9}, {6, 10}, {7, 11}, {8, 12}. 図2. . ˘ ¯ V = {1, 2}, E = {1, 2} で与えられる場合に x = (0, 0)T , (0, 1)T , (1, 0)T , (1, 1)T の式 (1) の確率 Pr{X = x} の値の大小関係.. (2). 2. c 2009 Information Processing Society of Japan °.

(4) Vol.2009-CVIM-167 No.36 2009/6/10. 情報処理学会研究報告 IPSJ SIG Technical Report . . 次に. V = {1, 2, 3, 4, 5, 6, 7, 8, 9}. ©. {1, 4}, {2, 5}, {3, 6}, {4, 7}, {5, 8}, {6, 9}. ª.

(5). . . . . . . (5) . . . . . . で与えられるグラフ (V, E) において画素 5 に着目するという場合を考える. 画素 5 以外の画素 i の輝度値がすべて 1, すなわち xi = 1(∀i∈V \{5}) であるとき画素 5 の輝度値は x5 = 0. . (4). E = {1, 2}, {2, 3}, {4, 5}, {5, 6}, {7, 8}, {8, 9},. . . .

(6). . . . . . . .

(7). . . . . . . . . .

(8). . .

(9). . . . . . と x5 = 1 で比べれば確率 Pr{X = x} の値が大きくなるのは x5 = 1 のときである (図 3 参照). . 図 4 V = {1, 2, · · ·, 12}, E が式 (2) で与えられる場合に画素 6 と画素 7 に着目し, それ以外の画素の輝度値が!xi = 1(i∈{1, 2,!5, 9, 10}), ! ! xi != 0(i∈{3, 4, 8, 11, 12}) であるとき画素 6 と 7 の輝度値 x6 0 1 0 1 は = , , , で比べれば確率 Pr{X = x} の値が最も大きくなるのは x7 0 0 1 1 ! ! x6 1 = のときである. x7 0. . . . . . . ©. . だきたい. ただし, Pr X = x . . . . さいころをふれば常に x1 = x2 = · · · = x|V | を満たす画像 x がでるかといえばそうではな. . . い. むしろ, まず E に属する {i, j} のなかで xi = xj を満たす {i, j} の本数が多い (xi 6= xj. . . . の値に比例して x がでやすくなっている. もちろん,. x1 = x2 = · · · = x|V | のときに一番確率が大きくなり, その面が最もでやすくなる. しかし,. . ª. . を満たす最近接画素対 {i, j} の本数が少ない) 画像 x ほどでやすくなると考えていただきた. ©. 図 3 V = {1, 2, · · ·, 9}, E が式 (5) で与えられる場合に画素 5 に着目し, それ以外の画素の輝度値が xi = 1(∀i∈V \{5}) であるとき画素 5 の輝度値は x5 = 0 と x5 = 1 で比べれば式 (1) の確率 Pr{X = x} の値が大きくなるのは x5 = 1 のときである.. ª. い. 式 (1) の確率分布 Pr X = x をもとにマルコフ連鎖モンテカルロ (Markov Chain. Monte Carlo: MCMC) 法4),5) による画像生成プログラムを構成し, これにより生成された画像 x = (x1 , x2 , · · ·x|V | )T を図 5 に示す.. . もう少し大きなグラフとして V = {1, 2, · · ·, 12} で E が式 (2) で与えられる場合を考. .

(10) . . . えてみよう. このグラフ (V, E) 上で画素 6 と画素 7 に着目し, それ以外の画素の輝度値が xi != 1(i∈{1, xi = 4, 8, 11, 12}) であるとき画素 6 と 7 の輝度値 Ã Ã ! 2,Ã5, 9,!10}), Ã ! Ã 0(i∈{3, !. x6 x7. るのは. =. Ã. x6. 0. ,. 0. !. =. 1 0. Ã ! 1. ,. 0 1. ,. 1 1. を比べれば確率 Pr{X = x} の値が最も大きくな ˘ ¯ 図 5 式 (1) の確率分布 Pr X = x をもとに構成されたマルコフ連鎖モンテカルロ法により生成された画像 T x = (x1 , x2 , · · ·x|V | ) .. のときである. つまり, 確率 (1) に従えば, 画素 6, 7 は周りの画素. x7 0 の状態をみながら自分の画素の状態が決められていることを示している.. もっと画素の数が多い場合を考えてみよう. まず, 2|V | 面体のさいころを考えていた. 図 5 で α の小さな時に本来は確率の低いはずの画像が生成されている. このことを説明. 3. c 2009 Information Processing Society of Japan °.

(11) Vol.2009-CVIM-167 No.36 2009/6/10. 情報処理学会研究報告 IPSJ SIG Technical Report. するために, xi 6= xj を満たす最近接画素対 {i, j} の本数. K(x) ≡. X. 2. {i,j}∈E. (xi − xj ) =. X. {i,j}∈E. (1 − δxi ,xj ). µi ≡. (6). ©. ª. いうことができる. ところがこの本数 K(x) が全最近接画素対のなかの 10% 以下の画像 x が出. ¯ ≡ {x|K(x)/|E| > 0.1} の確率 Pr{A} ¯ ≡ とその余事象 A. X. {x|K(x)/|E|≤0.1}. X. {x|K(x)/|E|>0.1}. 考えると状況は変わってくる. α がある程度以上に大きいときは不等式 ¯ Pr{A} > Pr{A}. ©. Pr X = x. ©. ª. 1 X. xi Pr{X = x}. (10). x|V | =0. 1.7627... において最大値をいることがわかる. つまりゆらぎが最大となるということであ. ª. る. 図 5 で α = 2 のときの画像はこのゆらぎが最大となる点の付近で生成される画像と見. .

(12) . Pr X = x を. (7). が成り立ち, 上で述べてきたことと同じである, しかし, α がある程度より小さくなると不等号の向きが逆転し,. ¯ Pr{A} < Pr{A}. x1 =0x2 =0. ···. |V |→ + ∞ における式 (9) の共分散 Cov[Xi , Xj ] の α 依存性は図 6 の通りであり, α =. という量を導入する. 与えられた画像 x は K(x) が小さいほど確率 Pr X = x の高い画像と現するという事象 A = {x|K(x)/|E| ≤ 0.1} の確率 Pr{A} ≡. 1 1 X X. (8). という不等式が成り立ってしまう. 後者の場合は事象 A に属するひとつひとつの画像の確率 Pr{X = x} は大きいが, その事象の標本点の個数 |A| が少ないため, 総和でみたときに. 図6. 確率 Pr{A} が小さくなってしまう. α が大きいほど xi = xj ({i, j}∈E) を満たす画像 x が. . . . . |V | → +∞ における式 (9) の共分散 Cov[Xi , Xj ].. でやすくなるが, 小さければ xi = xj ({i, j}∈E) を満たさない画像 x のなかのひとつがでやすくなる. 以上が図 5 で α の大きさによりさまざまの性質をもつ画像 x が生成されること. なすことができる. 通常の自然画像を 2 値化した画像と図 5 の α = 2 のときの生成画像を. の解釈である.. 比較すると, 全体としてみれば似ているとは言い難いが部分的なパターンを取り出してみる. 不等式 (7) から不等式 (8) に切り替わる α の値がどこかにあるはずである. しかもその切. とそこに類似性が見いだされる (図 7 参照). この類似性を画像処理に利用するというのがベ. り替わりの点でどのようなことがおこっているのかも興味深い. 図 5 で α = 2 のときの画. イズ推定による画像処理の発想の基本である. Q = 256 の場合の確率分布 (1) におけるマ. 像が α の小さいときの性質と α の大きいときの性質をあわせもつ画像であることに気がつ. ルコフ連鎖モンテカルロ法による生成画像 x は図 8 の通りである. ちょっと違う視点から図 3 と図 4 の例を説明しよう. 図 3 の例で V \{5} の画素の状態が. くであろう. このときなにかのゆらぎが大きくなっていることが予想される. どのようなゆらぎが大きくなっているのであろうか. 最近接画素対の輝度値の相関に対するゆらぎと考え. 固定されている時の画素 5 の状態に対する確率は条件付き確率をつかって次のようにあら. るのが自然であろう. そのゆらぎをあらわす量のひとつとして共分散. わされる.. Cov[Xi , Xj ] ≡. 1 1 X X. x1 =0x2 =0. ···. 1 X. (xi − µi )(xj − µj )Pr{X = x} (∀{i, j}∈E). Pr{X5 = x5 |X1 = 1, X2 = 1, X3 = 1, X4 = 1, X6 = 1, X7 = 1, X8 = 1, X9 = 1} Pr{X1 = 1, X2 = 1, X3 = 1, X4 = 1, X5 = x5 , X6 = 1, X7 = 1, X8 = 1, X9 = 1} = Pr{X1 = 1, X2 = 1, X3 = 1, X4 = 1, X6 = 1, X7 = 1, X8 = 1, X9 = 1} (11). (9). x|V | =0. があげられる. µi は Xi の期待値である.. 4. c 2009 Information Processing Society of Japan °.

(13) Vol.2009-CVIM-167 No.36 2009/6/10. 情報処理学会研究報告 IPSJ SIG Technical Report. . . 図 7 通常の自然画像を 2 値化した画像と図 5 の α = 2 のときの生成画像の類似性.. . . . . . . .   X1 x1      X2     x2     x2     X2     X3     x3       X3     x3   X4     x4      X    x     4    4 XV =  X5  , XV \{5} =   , xV =  x5  , xV \{5} =      X6     x6   X6     x6       X7     x7  X    x     7    7        X8   x8   X8   x8 . . X1. X9. X9.

(14) . x1. x9. (13). x9. という記号を導入する. この記号をつかうと式 (11)-(12) は画素 5 以外の画素の輝度値が. xi (∀i∈V \{5}) であるときの画素 5 の輝度値に対する確率は次のように与えられる. Pr{X5 = x5 |XV \{5} = xV \{5} } = Pr{XV \{5} = xV \{5} } =. 図 8 式 (1) の確率分布におけるマルコフ連鎖モンテカルロ (Markov Chain Monte Carlo) 法による生成画像 x.. 1 X. x5 =0. Pr{XV = xV }. 式 (14)-(15) に式 (1) を代入すると次の式が得られる. Y ¢ ¡ 1 Pr{X5 = x5 |XV \{5} = xV \{5} } ∝ exp − α(x5 − xj )2 2. Pr{X1 = 1, X2 = 1, X3 = 1, X4 = 1, X6 = 1, X7 = 1, X8 = 1, X9 = 1} =. 1 X. Pr{XV = xV } Pr{XV \{5} = xV \{5} }. Pr{X1 = 1, X2 = 1, X3 = 1, X4 = 1, X5 = x5 , X6 = 1, X7 = 1, X8 = 1, X9 = 1}. (14). (15). (16). j∈∂5. x5 =0. (12). ∂5 ≡ {2, 4, 6, 8}. (17). Pr{X5 = x5 |XV \{5} = xV \{5} } = Pr{X5 = x5 |Xk = xk , ∀k∈∂5}. (18). ∂5 は画素 5 のすべての最近接画素の集合をあらわしている. 式 (16). 式が長いので. が成り立つことを示している.. 同様にして, 図 4 の例で V \{6, 7} の画素の状態が固定されている場合に. 5. c 2009 Information Processing Society of Japan °.

(15) Vol.2009-CVIM-167 No.36 2009/6/10. 情報処理学会研究報告 IPSJ SIG Technical Report. . . X1. . . . x1. . 式 (22) と式 (22) はいずれも着目画素以外の画素の状態を固定したとき, 着目画素の状態. . . はその最近接画素の状態のみによって決定されることをあらわしている. この性質をもつ確. X  x  X1 x1  2  2     X     X3     2  x2     x3       X4     X3   x3     x4           X4   x4   X5   x5          X  x   X6   x6   5  5     XV =   , xV =   (19)  , XV \{6,7} =   , xV \{6,7} =   X8   x8   X7   x7           X9   x9   X8   x8               X9   x9  X x     10 10         X       X x x 11 11   10   10       X11   x11  X12 x12 X12. 率場 XV をマルコフ確率場 (Markov Random Field: MRF) とよばれ, 次の等式としてまとめられる.. Pr{Xi = xi |XV \{i} = xV \{i} } = Pr{Xi = xi |Xk = xk , ∀k∈∂{i}} (∀i∈V ). の正方格子からなるグラフ (V, E) に対して書き下すと次のような表式としてまとめられる. Y ¡ 1 ¢ exp − α(xi − xk )2 (26) Pr{Xi = xi |XV \{i} = xV \{i} } ∝ 2 k∈∂i. Pr{Xi = xi , Xj = xj |XV \{i,j} = xV \{i,j} } ³ Y ¡ 1 ¢´ ∝ exp − α(xi − xk )2 2 k∈∂i\{j}. x12. ¡. ×exp −. という記号を導入することで, 画素 6 と 7 の状態に対する確率は条件付き確率をつかって次のようにあらわされる.. Pr{X6 = x6 , X7 = x7 |XV \{6,7} Pr{XV \{6,7} = xV \{6,7} } =. Pr{XV = xV } = xV \{6,7} } = Pr{XV \{6,7} = xV \{6,7} }. 1 1 X X. x6 =0x7 =0. Pr{XV = xV }. ¢ 1 ×exp − α(x6 − x7 )2 2. ³ Y. k∈∂j\{i}. ¡. exp −. ´. ¢ 1 α(xj − xk )2 2. (27). (20). マルコフ確率場の性質をもつ確率分布 (1) は画像のさまざまの性質をあらわすことが可能であるが, それだけでは個々のデータに含まれる情報をとりこむことができない. この. (21). データをとりこむメカニズムが必要である. ノイズ除去の場合にはデータはもともとの画像. x にノイズが加えられることで生成されたと考える. 最も基本的はノイズは加法的白色ガ. Pr{X6 = x6 , X7 = x7 |XV \{6,7} = xV \{6,7} } ³ Y ¡ 1 ¢´ exp − α(x6 − xk )2 ∝ 2. ¡. ¢ 1 α(xi − xj )2 2. 3. ノイズ生成モデルとベイズ推定. ウスノイズ (Additive White Gaussian Noise) である. 各画素ごとに独立に発生した. 式 (20)-(21) に式 (1) を代入すると次の式が得られる.. k∈∂6\{7}. (25). ∂i は画素 i のすべての最近接画素の集合をあらわしている. また, 式 (22) と式 (22) を一般. ³ Y. k∈∂7\{6}. ∂{6} ≡ {2, 5, 7, 10}, ∂{7} ≡ {3, 6, 8, 11}. ni (i∈V ) という値による構成されるノイズ n = (n1 , n2 , · · ·, n|V | ) が画像に加算されて生成. されたのがデータ y = (y1 , y2 , · · ·, y|V | ) であると仮定する.. ¡. ¢ 1 exp − α(x7 − xk )2 2. ´. y =x+n. (28). データ y = (y1 , y2 , · · ·, y|V | ) の確率ベクトルを Y = (Y1 , Y2 , · · ·, Y|V | ) とあらわすことにす. (22). ると加法的白色ノイズの確率は次の式で与えられる (図 9 参照). Y ¡ ¢ 1 Pr{Y = y|X = x} ∝ exp − 2 (yi − xi )2 2σ. (23). やはり. (29). i∈V. Pr{X6 = x6 , X7 = x7 |XV \{6,7} = xV \{6,7} }. = Pr{X6 = x6 , X7 = x7 |Xk = xk , ∀k∈∂{6, 7}}. Pr{Y = y|X = x} はもとの画像 x が与えられたときのデータ y が生成される条件付. き確率であるが今ほしいのは逆で, データ y が与えられたときのもとの画像 x に対する条. (24). が成り立つこととなる.. 件付き確率 Pr{X = x|Y = y} である. これは式 (1) と式 (29) からベイズの公式 (Bayes. 6. c 2009 Information Processing Society of Japan °.

(16) 情報処理学会研究報告 IPSJ SIG Technical Report. -. 0. -. . . Vol.2009-CVIM-167 No.36 2009/6/10. . .. / /. 1 2 3 0 1 2 3 4 7- 5 +- - 4 -65 -+- - . . 3. 3 4 4 5 5 6 7 8 9 6 7 8 9 : 3>; 303 3 : 3<; 3=3 3 . . .

(17). % & "!#!#!$ %'& ()*% +*& !#!#!,+ %'& . 図9. . 式 (29) の条件付き確率分布のグラフ表現..

(18)

(19) # " ! %$', & '- &#()(*(+& ,- 2 $./&, 01&- #()()(+&0 ,- 2. 図 11 式 (30) の事後確率のグラフ表現.. Pr{X = x|Y = y} について, 各確率変数ごとのモーメントや共分散などの統計量を計算する問題に帰着される. 例えば確率変数 Xi の期待値 (1 次のモーメント) を求めようとする場合, 定義に従えば. XX. Q−1 Q−1. E[Xi ] =. x1 =0x2 =0. X. Q−1. ···. xi Pr{X = x|Y = y}. (31). x|V | =0. という操作をおこなうことになる. これは一般には O(exp(|V |)) の計算量を必要としてし. . 図 10. . まう. この統計量を (数値計算の誤差を除いたとしても) 厳密なものではないが問題によっ. . ては厳密計算の値に近い値が得られるように設計することができる近似計算技法として確. 式 (29) の条件付き確率分布によるデータ y の生成過程.. 率伝搬法 (Belief Propagation) が最近注目をあつめている. 確率伝搬法は各最近接格子点対 {i, j} に対して Mi→j (xj ) と Mj→i (xi ) の 2 種類のメッ. Formula) をつかって次のように導かれる. Pr{Y = y|X = x}Pr{X = x} Pr{Y = y} ³Y ¡ ¢´³ Y ¡ 1 ¢´ 1 ∝ exp − α(xi − xj )2 exp − 2 (yi − xi )2 2σ 2. セージ (Message) とよばれる量が導入される. 式 (30) で与えられた事後確率に対しては. Pr{X = x|Y = y} =. i∈V. このメッセージは次のメッセージ伝搬規則 (Message Propagation Rule) によるくりかえし計算を収束するまで実行することで計算される (図 12 参照).. (30). {i,j}∈E. 率 (Posterior Probability) とよばれる. また, 確率 Pr{X = x} は事前確率 (Prior. Probability) とよばれる. 式 (1) の事前確率 Pr{X = x} と式 (30) の事後確率 Pr{X = x|Y = y} はそれぞれ図 1 と図 9 の正方格子によるグラフ表現であたえられる.. X. Q−1. Mj→i (xi ) ← Constant×. ベイズの公式により導出された式 (30) の条件付き確率 Pr{X = x|Y = y} は事後確. z=0. ¡. exp −. ¢ 1 1 α(xi − z)2 − 2 (z − yj )2 2 2σ. ×. Y. Mk→j (z) (∀{i, j}∈E). (32). k∈∂j \{i}. 具体的な計算は結合確率. 7. c 2009 Information Processing Society of Japan °.

(20) Vol.2009-CVIM-167 No.36 2009/6/10. 情報処理学会研究報告 IPSJ SIG Technical Report. X. Q−1. Mi→j (xj ) ← Constant×. z=0. ¡. exp −. ¢ 1 1 α(xj − z)2 − 2 (z − yi )2 2 2σ. ×. Y. Mk→j (z) (∀{i, j}∈E). いただきたい. 確率伝搬法は人工知能の確率推論におけるアルゴリズムとして 1980 年代に考案され, その後, 誤り訂正符号における低密度パリティ検査 (Low Density Parity Check:. LDPC) 符号, 移動体通信の符号分割多重接続 (Code Division Multiple Access: CDMA) (33). 方式, 組み合わせ最適化における充足可能性 (Satisfiability: SAT) 問題に対するキーアプリ. ∂j は画素 j の最近接画素の集合, ∂j \{i} は画素 j の最近接画素の集合から画素 i をのぞい. 献3),6)–8) を参照していただきたい. さらに確率伝搬法は物理学における平均場理論との類. k∈∂i \{i}. ケーションとしても新しい計算の概念を生み出しつつある. その最近の展開については文. た集合を意味する. メッセージ伝搬規則 (32)-(33) にあらわれる. Q. k∈∂j \{i}. 似性が指摘され，情報統計力学の研究テーマとしても研究が進められている3),6),9) ．. は, 高々, 各画素. の最近接画素数なので計算量は O(1) となる. つまり, くりかえし計算によるアルゴリズム. 4. 情報統計力学と統計的性能評価. の 1 回の更新ごとの計算量は O(|V |) というわけである. モーメントや共分散はメッセージ. {Mi→j (x), Mj→i (x)|∀{i, j}∈E} をもちいてあらわすことができるので, 得られたメッセー. 前節で説明したベイズ推定による画像処理システムの性能の評価をするひとつの方法として. ジを代入することで求められる. 図 13 に示した出力結果はこの確率伝搬法により得られた. . . . . . . . .

(21). . たくさんのサンプルを生成して標本平均 (Sample Average) を計算するという方法がまず考えられる. これは, まず事前確率 (1) に従ってもとの画像 x を 100 個生成し, それぞれに対. . . . して各画像 x ごとに条件付き確率 (29) に従って 100 個ずつのデータ y を生成する. そうして, ~ Y ~ =~ 得られた 10000 個のサンプルデータ y から Pr{X| y } に従ってもとの画像 x の推定値. . x b=x b(y) をそれぞれに求めた上で x と x b との平均二乗誤差 d(x, x b) =. 1 |V |. P. i∈V. (xi − x bi )2. についての標本平均を計算すると, この事後確率に対する画素あたりの誤差を実験的に評価できる. しかし, これではいずれの場合も 10000 回のテストを行わなければならないことになる. また, より信頼性の高い結果を得るためには 1000000 回, 100000000 回とテスト回数. 図 12 確率伝搬法におけるメッセージ伝搬規則のグラフ表現.. を増やす必要があるかもしれない. ところが, このプロセスを式で表すと. XX. Q−1 Q−1. ものである.. x1 =0x2 =0.

(22) . X ³Z. Q−1. ···. x|V | =0. +∞ −∞. Z. +∞. −∞. ···. Z. +∞. −∞. ¡. ¢. d x, x b(y). ´. ×Pr{Y = y|X = x}dy1 dy2 · · ·dy|V | Pr{X = x}. (34). という量を計算することに相当し, これを計算できてしまえばテストを 10000 回行わなくてもよいことになる. 実は式 (34) であらわされる量は 20 世紀の後半に多くの物理学者がスピングラス理論 (Spin Glass Theory) とよばれる分野のなかで挑戦してきた量, すなわち配位平均 (Configuration Average) の初等的な場合に相当する. このことは, 20 世紀後半に行われたスピングラス理論に関する多くの理論的成果が確率的情報処理の問題の理解図 13. 確率伝搬法によるノイズ除去例.. に役立てられる可能性を示しており, 実際に画像処理はもとより, 誤り訂正符号, 移動体通信における CDMA 復調方式, 組み合わせ最適化などの様々な問題で情報統計力学によるアプ. 確率伝搬法を用いた画像処理の基本的なアルゴリズム更正法の詳細は文献2) を参照して. ローチが新しい研究成果をあげている3),10)–14) .. 8. c 2009 Information Processing Society of Japan °.

(23) Vol.2009-CVIM-167 No.36 2009/6/10. 情報処理学会研究報告 IPSJ SIG Technical Report. 7) C. M. Bishop：Pattern Recognition and Machine Learning, Springer, 2006. 8) M. J. Wainwright and M. I. Jordan: Graphical Models, Exponential Families, and Variational Inference, now Publishing Inc, 2008. 9) M. Opper and D. Saad (eds): Advanced Mean Field Methods —Theory and Practice, MIT Press, 2001. 10) 西森秀稔: スピングラス理論と情報統計力学, 岩波書店, 1999. 11) H. Nishimori: Statistical Physics of Spin Glasses and Information Processing, —An Introduction, Oxford University Press, 2001. 12) 樺島祥介: 物理の世界/学習と情報の平均場理論, 岩波書店, 2002. 13) Y. Kabashima and D. Saad: Statistial Mechanics of Low-Density Parity-Check Codes (Topical Review), Journal of Physics A: Mathematical and General, Vol. 37, No.6, pp.R1-R43, 2004. 14) M. M´ezard and A. Montanari: Information, Physics, and Computation, Oxford University Press, 2009. 15) A. S. Willsky: Multiresolution Markov models for signal and image processing, Proceedings of IEEE, vol.90, no.8, pp.1396-1458, 2002.. 5. まとめ本稿ではノイズ除去を例にとって確率的画像処理について解説した. 事前確率 (1) を領域分割, エッジ検出などの用途にあわせて適切に変更することにより, さまざまの確率的画像処理フィルターを設計することができる. そのさまざまの用途に対する最近の展開については文献15) で概観することができる. 確率的画像処理の情報統計力学的展開についての詳細は文献1),11) を参照していただきたい. 謝辞本研究の一部は文部科学省科学研究費補助金 (No.18079002, No.18079008) および東北大学電気情報系グローバル COE プログラム「情報エレクトロニクスシステム教育研究拠点」の補助のもとでおこなわれたものである.. 参考. 文. 献. (平成 21 年 5 月 9 日受付). 1) K. Tanaka: Statistical-Mechanical Approach to Image Processing (Topical Review), Journal of Physics A: Mathematical and General, Vol.35, No.37, pp.R81R150, 2002. 2) 田中和之: 確率モデルによる画像処理技術入門, 森北出版, 2006. 3) 田中和之 (編著): 臨時別冊・数理科学 SGC ライブラリ「確率的情報処理と統計力学— 様々なアプローチとそのチュートリアル」, サイエンス社, 2006. 4) 伊庭幸人: ベイズ統計と統計物理, 岩波書店 (2003). 5) 伊庭幸人, 種村正美: 統計科学のフロンティア計算統計 II —マルコフ連鎖モンテカルロ法とその周辺, 岩波書店, 2005. 6) 汪金芳, 田栗正章, 手塚集, 樺島祥介, 上田修功: 統計科学のフロンティア計算統計 I —確率計算の新しい手法, 岩波書店, 2003.. (平成. 年. 月. 日採録). 田中和之昭和 36 年生. 平成元年東北大学大学院工学研究科電子工学専攻博士課程修了. 同年東北大学工学部助手. 平成 6 年室蘭工業大学工学部助教授. 平成 11 年東北大学大学院情報科学研究科助教授. 平成 19 年より東北大学大学院情報科学研究科教授. 確率的情報処理と情報統計力学に関する研究に従事. 工学博士. 日本物理学会, 電子情報通信学会, 人工知能学会, 情報理論とその応用学会, IEEE-CIS 各会員.. 9. c 2009 Information Processing Society of Japan °.

(24)