ファジーハフ変換によるロバスト画像処理

九州大学学術情報リポジトリ

Kyushu University Institutional Repository

ファジーハフ変換によるロバスト画像処理

岡田, 正之

九州芸術工科大学

https://doi.org/10.11501/3168352

出版情報：Kyushu Institute of Design, 1999, 博士（工学）, 課程博士

の

ファジーハフ変換による ロバスト画像処理

Robust Image Processing by Fuzzy Hough Transform

1999年12月

岡 田 正 之

Masayuki OKADA

目次

第1章 序論 1

1.1 研究の目的と背景 . • • • • • • • • • • • • • • • • • • • • • • 1 1.2 論文の概要と構成

.

• • • • • • • • • • • • • • • • • • • • • • 4

第2章 口パス卜回帰による画像処理 ⁷

2.1 まえがき . ^{• • • • • • • • • • • • • • • • • • • • • • • • • • • 7}

2.2 GNC法による非線形回帰 • • • • • • • • • • • • • • • • • 9 2.2.1

GNC法による回帰アルゴリズム .

• • • • • • • • • • 9 2.2.2 例題 . ^{• • • • • • • • • • • • • • • • • • • • • • • • • 13}

2.2.3 むすび . ^{• • • • • • • • • • • • • • • • • • • • • • • • 13}

2.3 最小2乗トリム回帰のアナログ解法 . • • • • • • • • • • • • 15 2.3.1 LTS推定 • • • • • • • • • • • • • • • • • • • • • • • 15 2.3.2 アナログ解法

.

• • • • • • • • • • • • • • • • • • • • 16 2.3.3 実験例

.

^{• • • • • • • • • • • • • • • • • • • • • • • • 19}

2.3.4 むすび . ^{• • • • • • • • • • • • • • • • • • • • • • • • 19}

2.4 ロバストエネルギー最小化によるエッジ保存平滑化 . • • • 25 2.4.1 エネルギー最小化による平滑化

.

• • • • • • • • • • 25 2.4.2 最小2乗エネルギー • • • • • • • • • • • • • • • • • 25 2.4.3

M推定エネルギー .

• • • • • • • • • • • • • • • • • • 26 2.4.4

LMedS推定エネルギー

. • • • • • • • • • • • • • • • 27 2.4.5 実験例

.

^{• • • • • • • • • • • • • • • • • • • • • • • • 29}

2.4.6 むすび

.

^{• • • • • • • • • • • • • • • • • • • • • • • • 30}

第3章 ファジーハフ変換 38

3.1 まえがき

.

^{• • • • • • • • • • • • • • • • • • • • • • • • • • • 38}

3.2 可変型テンプレートマッチング

.

• • • • • • • • • • • • • • 38 3.3 ファジーハフ変換

.

• • • • • • • • • • • • • • • • • • • • • • 41 3.4 むすび ^{• • • • • • • • • • • • • • • • • • • • • • • • • • • • 44}

第4章 ファジ-ハフ変換による画像処理 49

4.1 まえがき .

_{• • • • • • • • • •} • • • • • • • • • • • • • • • • •

49 4.2 ファジーハフ変換による関数回帰 .

• • • • • • • • • • • • •

51 4.3 代表点による補間 .

_{• • • • •} • • • • • • • • • • • • • • • • •

60 4.4 画像の領域分割 .

_{• • • • • •} • • • • • • • • • • • • • • • • •

64 4.5 ランダムドットステレオネ見.

• • • • • • • • • • • • • • • • •

66 4.6

連想記憶

.

^{• • • •} • • • • • • • • • • • • • • • • • • • • • • •

74 4.7

むすび

.

• • • • • • • • • • • • • • • • • • • • • • • • • • • •

75

第5章 ア二一リングファジ-ハフ変換

5.1 まえがき .

5.2 アニーリングファジーハフ変換

5.2.1 ロバスト推定としての表式 .

_{• • •} • • • • • • • • • •

82 5.2.2 アニーリング .

_{• • • • • •} • • • • • • • • • • • • • •

83 5.2.3 オプテイカルフロー推定への応用 .

• • • • • • • • •

87 5.3 ニューラルネットによるファジーハフ変換 .

• • • • • • • •

91 5.4

むすび . • • • • • • • • • • • • • • • • • • • • • • • • • • • •

94 78

78 81

第6章 フアジハフ変換による運動検出

6.1 まえがき .

6.2 運動知覚モデル . 6.2.1 局所運動の検出

氏リ ハb FO

円i

ou ny ny QU

6.2.2 ファジーハフ変換による最尤推定 .

• • • • • • • • •

98 6.2.3 協調による統合 .

• • • • • • • • • • • • • • • • • • •

105 6.2.4

むすび

.

• • • • • • • • • • • • • • • • • • • • • • • •

110 6.3

動き 検出

.

• • • • • • • • • • • • • • • • • • • • • • • • • • •

117 6.3.1

ファジーハフ変換による運動検出

.

_{• • •} • • • • • •

117 6.3.2

逐次運動抽出 • • • • • • • • • • • • • • • • • • • •

119 6.3.3

実験

.

• • • • • • • • • • • • • • • • • • • • • • • • •

120

6.3.4

むすび

. 121

第7章 モードフィルタ

7.1 まえがき .

128

128

7.2

モードフィルタ . • • • • • • • • • • • • • • • • • • • • • • •

128

7.2.1

最大事後確率推定によるフィルタリング

.

_{• • • • •}

129

7.2.2 モードフィルタ .

• • • • • • • • • • • • • • • • • • •

130

」

7.2.3 カラー座標について ^{- ・ ・} . 131

7.2.4 簡略モードフィルタ . 131

7.3 むすぴ ^{- ・ ・ ・ ・ ・} 134

第8章 結論 137

」

第1章 序論

1.1 研究の目的と背景

コンビュータや通信技術の発達により3 現在社会に浸透したインター ネットもマルチメディア化が進み， WWW (World Wide Web)で，

世界各地からのカラー静止画像や動画像を見ることができる. 例えば， 火 星探査機から米国NASAに送られた火星表面のカラー画像をWWWで 見ることができる. これは3 大量の画像データを扱う画像処理技術や画 像圧縮技術に支えられている. さらに動画像を扱うためには動画像処理 技術や動画像圧縮技術が必須である.

これらの静止画像や動画像のディジタル化において3 多少とも画像に は色々な種類の歪を生じる. 例えばカメラ等のような光学系機械を用い て得られる画像は， カメラレンズの口径によって帯域が限定されたり， さ らにレンズの収差による歪を受ける. またカメラの動揺によるボケが生 じることもある. 走査画像においても3 走査線の拡がりの幅， サンプリ ング間隔， 量子化レベルなどによる歪を生じる. その他に通信回線上で も3 画像の圧縮， 復元などによりノイズが付加される. このように光学 系 ， 検出機器， 増幅器3 量子化， 各種計算機処理過程においてノイズが 混入する. このようなノイズは， 画像の有効利用の面から大きな障害に なることがある.

画像データ自体が不連続な変化成分を含む非定常， 非ガウス分布信号 であり， これに上記のノイズ等が加わって生成画像データとなる. このノ イズ成分にはしばしば突発的な大振幅のものが含まれることがある. そ の場合， 線形でシフト ・ インバリアントな(場所に依存しなしì )処理で は， この突発的ノイズに対応できず， 処理出力がその突発的なノイズに 大きく影響されて不自然なものとなる.

これは， ノイズの統計的性質を見誤ったために生じるもの であると解 釈できる. 突発的なノイズもある確率分布に従うはず で， その確率分布 の推定が不適当であるためにノイズが除去できず信号成分も歪ませてし

まうのである. しかし現実的には， ノイズの統計的性質を完全に把握す ることはできないので， ノイズが定常ガウス成分の他にもこのような非 定常， 非ガウスノイズも含むような場合でも， ノイズの統計的性質を正 確には知らなくても安定して対処できるような画像処理技術が必要とな る. 本論文では， このような問題に対処する方法として， ロバスト性が 高い非線形な回帰アルゴリズムを提案し， 簡単な実験を行ってその有効 性を示すとともに種々の画像処理への応用を試みる. 本論文で扱う画像 処理技術はノイズ除去， 画像平滑化， 画像セグメンテーション， 動画像 処理などである.

まず，画像平滑化については，線形な信号処理では除去することが難し い雑音(例えばインパルス性ノイズ(図2.3を参照))なとを対象に，各 種の順序統計に基づく非線形処理が提案されて一定の成果を収めている [1] . その原型はメジアンフィルタであり，線形処理では不可能であった信 号の急激な変化点(エッジ)保存とインパルス性ノイズ除去を同時に達 成できるものである. メジアンフィルタは処理が単純な割にはこのよう に優れた特長を持つが，次のような欠点もある. 1つ目は，非インパルス 性ノイズ3 例えばガウス性白色ノイズなどの除去能力が線形処理に比べ て劣っていることで， 2つ自は， ステップ状のエッジ信号は保存するが，

細かく変化する信号を劣化させてしまうことである. そこで， 信号の復 元精度を向上させるために，単純なメジアンフィルタを重み付きメジア ンフィルタに拡張したり[2]， メジアン以外の順序統計フィルタョ 例えば α-トリム平均フィルタやDW-MTM(Double Window Modified Trimmed Mean)[2]などが提案されている. これらの改良によってノイズ除去能力 とともに画像の細部保存性も向上してきている. しかしながら，これら の非線形フィルタでも非定常， 非ガウスノイズの頻度が高くなると性能 が劣化してくる. 例えばメジアンフィルタは全データの半分まで非定常，

非ガウスノイズが含まれでも対処できるが， 非定常， 非ガウスノイズが 半分以上の割合になってくると急速に性能が劣化してくる. また， メジ アンフィルタに重みを付けたファイルタでも， ノイズが全データの半分 以上になるとメジアンフィルタと同程度の出力になる. そこで本論文で は従来の非線形フィルタより もロバスト性が更に高いフィルタを提案す る.

次に画(象のセグメンテーションは， コンビュータビジョンやパターン 認識の前処理段階として非常に重要な画像処理技術である. すなわち各 物体を構成する境界線や面を抽出することが，画像認識をするために必

」

要となる. 対象とする物体が背景とともに他の物体と混在する場合には セグメンテーションは不可欠な処理であり， 2次元， 3次元にかかわらず3 セグメンテーシヨンの善し悪しがその後の処理の性能に大きく影響する.

画像をいくつかの領域に分割するためには各物体の面を一つの領域とし て取り出す離散的な処理と， その領域の濃淡値やテクスチャを推定する 連続的な処理との両方が同時に行われなければならない. すなわち， 領 域の分割は， 一つの面のなかでの連続性と他の面との不連続性が同時に 問われる. 一般に複雑な背景から一つの図を浮かび上がらせるというこ とは， その図の連続性と他のパターンとの不連続性とを認識することで あり3 そこでは連続的な処理と離散的な処理の両方が行われることが不 可欠である. このような連続と不連続とが混在した画像復元のモデルと しては， 標準正則化理論に基づくRBFニューラルネット

[

³

]

^{やマルコフ}

ランダム場モデル

[

⁴

]

など最適化に基づく手法がこれまで多く提案されて いる. これらのモデルは元々は線形であるが， ラインプロセスと呼ばれ る離散変数を導入することにより非線形化して不連続も抽出できるよう にされた

[

⁵

]

. さらにアウトライアプロセスと呼ばれる離散変数も導入す ることによってインパルスノイズも除去できるような拡張も行われてい る

[ ⁶ ] ^.

^{しかしその結果，} 最適化問題は非凸になり， 解を求めるのが非常 に困難となる. そこで本論文ではこれらの手法とは異なるアプローチの 方法を提案する. 従来の手法が画像の濃淡値を画像平面の関数として表 していたのに対し3 本論文では画像平面と濃淡値とを合わせた空間での 確率分布として表現する. これによって不連続や隠蔽などが自然に扱え るようになる.

最後に動画像でも静止画像と同じような処理すなわちノイズ除去やセ グメンテーションが要求される. 静止画像にはない動画像に特有の画像 特徴は物体やカメラの動きであり， 連続した複数のフレームから画像中 の各物体の動きを検出する処理が必要となる. すなわちここでも上記の セグメンテーションと同様な処理が必要となり3 物体の動きがかかわっ てくることにより問題が複雑になってくる. 運動の検出法は大きくオプ テイカルフロー

[

⁷

]

^{とマッチング法}

[

⁸

]

とに分けられる. オプテイカルフ ローでは平滑化によってノイズが除去される. また， マッチング法では ブロックマッチングによってノイズの影響が平均化される. しかし， こ れらの方法でも領域分割と各領域の運動推定とを両立させることは非常 に難しく， これまでにいろいろな改良が加えられてきている. 動き領域 の分割法ではハフ変換を使って領域抽出を行う方法

[

⁹

]

^{やM推定などのロ}

」

バスト推定を使う方法

[

¹⁰

]

などが提案されている. またマッチング法で は 1画素マッチング法やアフィン動き推定と領域分割とを交互に行う反 復法

[

¹¹

]

などが提案されている. 本論文でもオプテイカルフローとマッ チング法の両方についてロバストな推定法を提案する.

1.2 論文の概要と構成

第1章では， 本研究の背景と扱っている問題を示し，あわせて論文 の概要について述べる.

第2章では，本論文の基礎となる非線形な関数回帰の大域的最適近似解 を求める手法について論じる. 非線形な関数回帰のほとんどは非凸な最 適化問題であるので多数の局所解を持ち3 最適解を求めるのが困難であ る. そこで 2.2節では，非凸な画像復元の最適解を求めるアナログ的な手 法であるGNC法

[

¹²

]

を代表的なロバスト推定法であるM推定とLMedS 推定に応用し，メジアン選択をアナログな非線形計画問題として定式化 してこれらの推定解を求める. 2.3節では，M推定よりもスケール変動に 関してロバストな最小2乗トリム推定解を求めることを試みる. この非 線形回帰法は基本的には混合整数計画問題であり，非常に多くの極小解 を持つ， そこでホップフィールドネットワークの解法に似たアナログ解法 を用い3 決定論的アニーリング法で近似解を求める方法を提案する. 2.4 節では，エッジを保存する平滑化法について， M推定に含まれる関値の

ようなパラメータを含まない方法としてLMedS推定に基づく平滑化法を 提案し， 簡単な実験を通して最小2乗法，M推定，メジアンフィルタな どと比較検討する.

第3章は，この章以降の基礎となるロバスト回帰の1種であるファジー ハフ変換を提案する. まず最初に3 ハフ変換がテンプレートマッチング やロバスト回帰と等価であることを最適化問題による定式化から導いて3 この最適化問題の目的関数にエントロビーに似た関数を付け加えること によりファジーハフ変換が得ら れることを示す. また3 テンプレートが 高次の関数で表される場合のように， データ点とテンプレートとの距離 が解析的に表せないような場合に対する近似解法も導く.

第4章では，ファジーハフ変換による関数回帰の有効性を示す

.

関数回 帰は，画像復元， 物体検出， 領域分割などの広範な画像処理の基礎とな る手法である. 画像処理やパターン認識での関数回帰では， ノイズの平

滑化とエッジなどに対応する関数値の不連続の保存が重要である. 最小 2乗法や正則化マルコフランダム場モデルは線形処理なのでエッジ等が なまる. また， 1価関数しか扱えず透明視などが扱えない. そこで， ロ バスト性の高い非線形のファジーハフ変換で関数回帰を行う. ここでは 関数回帰をテンプレートマッチングとして捉え， 独立変数と従属変数と を合わせた空間を画像平面と同一視し3 関数をこの空間でのテンプレー トとみなして陰関数として扱う. この手法を， 疎らなデータからの関数 回帰に基づくデータ点の削減， クラスタリングによる領域分割，ランダ ムドットステレオ視ヘ応用し透明視も再現できることを示す.

ファジーハフ変換においてもファジ一度の設定や初期値設定など実際 の適用において難しい問題がある. ファジーハフ変換のこのような難点 を解決するために，第5章ではアニーリング法を導入する. アニーリング 法は基本的に整数計画問題に導入されるものである. そこで， ここでは ファジーハフ変換を整数計画問題として定式化してアニーリングを導入 する. ファジーハフ変換にアニーリングを導入することによりファジ一度 の設定や初期値の設定問題の解決だけでなく， 複数個のパターンを抽出 することが容易になる. ファジーハフ変換でのテンプレート抽出の基本 的なステップは， 大域最適解すなわちパラメータ空間の分布の最大ピー クを1つ取り出すことである. 複数個のパターンを抽出するには，取り 出されたテンプレートに所属するデータを除去した後に次の大域最適解 を取り出すことを繰り返す. この手法により抽出の各段階では大域最適 解を1つ求めればよいことになる.

第6章では，ファジーハフ変換を動画像からの運動検出に応用する. 6.2 節では， 運動知覚のニューラルネットモデルを提案し， パーパーポール

錯視の例をシミュレーションしてJ�\理実験結果を再現する. 本モデルは 局所運動の検出に単純な相関だけを使う. シミュレーションにおける仮 定は， 1) 2値画像を扱う， 2)時間を離散化して局所運動を相関で検 出する， 3)運動情報は局所的なニユーロン結合により空間的に伝播す る， 4)窓枠による隠蔽の効果も取り入れる，等である. 6.3節では，ファ ジーハフ変換に基づいて動画像から主要な動きを逐次に抽出する方法を 提案する. ハフ変換はパラメータ空間における分布のモードを取り出す ため3 主要な動きから順に逐次取り出すことができる. 逐次抽出法では 抽出の停止条件を与えるだけでよく， 抽出された動きの主要順位もわか る. 複数個の運動パターンを抽出するにはアニーリングハフ変換で各パ ターンごとの大域最適解を求め， それに属するデータを除去して取り出

せばよい. 各領域の運動は， アフィン変換動きモデルで表す. 本方法は 大域的動き領域分割を与え， 運動分割と運動推定の2つの問題を同時に 解決する. 簡単な動画像の例で本方法の有効性を示す.

第7章では， カラー画像のエッジを保存する非線形平滑化フィルタと してモードフィルタを提案する. カラー画像では画素は3つの値を持つ ベクトルとして表され，これらの要素には互いに相関がある.このような 信号にフィルタ操作を行う場合， 特に非線形な処理を施す場合には， ベ クトルの各要素に対して個別に処理を施してもうまくいかない. そこで ここでは， ユークリッド距離に基づくカラー画像の非線形平滑化法とし てモードフィルタを提示し，カラー値をCIE-L*♂V座標に変換してフィ ルタ処理をする. また，計算量を削減するためにモードフィルタを簡略 化したフィルタがカラー画像のある種の情報圧縮であるスケッチ画像を 得るのに有用であることを示す.

最後に第8章では， 本研究のまとめと今後の課題を述べる.

第2章 口バスト回帰による画像 処理

2.1 まえがき

我々が手にし得る画像には大なり小なり， 必ずノイズが含まれてい る. たとえば，インコヒーレント照明された物体はすでにその表面にお いて波面の統計的ゆらぎに起因するノイズを含む. 光学系，検出機器3 増 幅器，量子化， 各種計算機処理の過程，また通信回線上においてもノイ ズが混入することは明らかである. このようなノイズは画像の有効利用 の面から大きな障害になる. 画像ノイズを除去する手法は色々研究され てきている. 微細なノイズ(図2.2)，たとえばノイズの発生確率がガウ ス分布に沿って発生し分布長が短いノイズなら線形処理法で除去できる が3 ノイズにはガウス分布ノイズのようなものだけでなくノイズの発生 位置(時間)や大きさが不規則なインパルス性ノイズ(図2.3)も含まれ る(図2.4). これらのノイズを線形処理法で除去することは非常に難し い. 特に，画像データなどはエッジ保存とインパルス性ノイズの除去を 同時に行なわなければならない. そこでこの章では，画像データにおけ るガウス分布ノイズやインパルス性ノイズ等を平滑化し，エッジ保存に も対応したロバスト性が高い非線形な回帰アルゴリズムを提案する.

2.2節では，非線形な関数回帰法であるM推定とL11edS推定について 最適解に近い解を求める方法としてGNC(Graduated Non-Convexity)法 によるアルゴリズムを提案する[13].

2.3節では，ロバスト推定法の1種である最小2乗トリム推定を混合整 数計画問題として定式化し，その近似解を決定論的アニーリングによる アナログ解法で求める方法を提案する[14] .

2.4節では，最小2乗法に基づく平滑化法である正則化法のインパルス ノイズの除去能力とエッジの保存性とを高めた方法として，LMedS推定 に基づく方法を提案する[15].

図2.1:信号 図2.2:ガウス分布ノイズ

図2.4:ガウス分布ノイズとイン 図2.3:インパルス性ノイズ

パルス性ノイズの混合

2.2 GNC法による非線形回帰

観測データに関数を当てはめる関数回帰は統計処理や信号処理の基 礎である. 最も基本的な最小2乗法は線形解法で解けるので計算が容易 であるが，外れデータに敏感である. そこでロバスト性が高い非線形な 回帰アルゴリズムが種々考案されている[16]，[17]

.

しかしこれらの非線 形な回帰のほとんどは非凸な最適化問題であるので多数の局所最適解を もち，最適解を求めるのが困難である. これまで個々の問題についてディ ジタルアルゴリズムが工夫されている[18]，[19]が，任意の次元の回帰関

数に適用できるような一般的なアルゴリズムはまだ構成されていない.

本節では非凸な画像復元の最適解を求めるアナログ的な手法であるGNC 法[12]を代表的なロバスト推定法であるM推定とLMedS推定に応用し た非線形回帰アルゴリズムを提案する. 本方法の特長は特にディジタル 的な処理であるメジアン選択をアナログな非線形計画問題として定式化 し ，GNC法によって大域最適解に近い 解を求めることである.

2.2.1

GNC法による回帰アルゴリズム

観測データ (Xl，

dr)，…，(Xη， dn)ぅXiξRP，diξRP に関数 j(x，B)，xε

RP，BεRT，jεRqを当てはめるとする.

最も基本的な最小2乗法では Oを次の最適化問題の解として 求める.

mln 玄Ilj(xi，

の- dil12 〆'z，、、 qL

11ム

、、El，J

ここで11.11 はユークリッドノルムである. これの解は T θj(xi，B)

乞[j(

X i， B) - dJ l

� J ':; �' _{� } =} 0

θ0 (2.2)

の解として得られる. 通常fはOに関して線形であるから式(2.2)は0 に関する連立l次方程式になり容易に解ける. しかしユークリッドノル ムでは外れデータほど大きなウエイトを持つので ，一つでも外れデータ があるとゆが正解から大きくずれる.

この難点、の解消法としてまず思いつくことは次式のようにノルムを変 えて外れデータの影響を弱めることである.

mln 乞ゆ[f(Xil B) - di] (2.3)

M推定では ゆは例えば次のような関数である.

、11ノ 、11/

α

<一> α

Uu u

fFss、、 〆'z，、、、

つL 円L

uu α /|〈11 一一 、、1『FノUu r'a，、、 AV (2.4)

ここでαはf から α以上離れた外れデータを無視するしきい値であり，

問題ごとに適切に設定する必要がある. 問題(2.3)は非凸な最適化問題 であり複数個の極小解を持つ. そこでここでは(2.3)の解をGNC法で 求めることにする. まず(2.3)を次のような極限として表す.

min Z ^{J Wglf(川)} - ^{di] (2.5)}

ここでの(y)はg竺Oでνの凸な関数となりg↑∞で式(2.4)のゆに なるような関数， 例えば

ゆ9

(y) = - ^{ln [1}

⁺

^{e仰_y2)] (2.6)}

である. GNC法では 式(2.5)の mlnとlimの順序を入れ替えて(局所 最適解を求めるのであればこの交換は許される)

min 2二4>g[f(Xil B) - di] (2.7)

の解を， 9 �

0から出発して g�∞まで連続的に追跡する.

9 � 0で は( 2.7 )の解は最小 2乗解であり，

9↑∞で(2.3)の局所最適解が得ら

れる. この解追跡法は決定論的アニーリング法の一種であり最適解に近

い解が得られる. 式( 2.7 )の解は

T θf(XilB)

� 4>�[f(Xil B) - d i]T V J ^��'

^{V J}

= ^{0 (2.8)}

」

の解として得られる. ここで%はのの導関数でありロバスト推定では 影響関数と呼ばれる[16]. ここでは式(2.6 )の 微分

め ^{( U)=} ₁

₊

_eg(y2-α) ^{U (2.9)}

を使うことにする. M推定は最小2乗法よりもロバストであるが， α

の値を適切に設定しなければならないという難点がある. 通常αの 適切な 値の見積りは困難である. そのようなパラメータを含まない推定法とし てLMedS推定[18]がある. M推定が 最小2乗(Least Mean of Squares略 して

L

MS)法の2乗(square )を関数のに変えたものであったのに対し 平均(Mean)をメジアン( Median)に変えたもの，すなわち式(2.1 )を

min med {llf (Xi， B) - di112} (2.10)

に変えたものがLMedS (Least ìv1edian of Squares)推定である. Z = med { 11 f ( X i， B) -di 112 }とおくとzは

� sgn (z - Ilf(Xi， B) -di112) = 0 〆'『t、、 qL

11ム11i

、、EE，，f

の解である. [20]. ここでsgn (ν)はUくOで-1 ， ^{Y >} 0で1となる 関 数である. sgn (ν) = limg→∞ ta凶(gy)であるから

�

tanh(g[zg -Ilf(九B)-di 112) ⁼ 0 (2.12 )

の解をZgとすると 式( 2.11 )の解

zはZ limg→∞Zgである.

すなわ ち式(2.10)の解はZ - limg→∞mlnZgとして得られる. mln Zgの解は θZg/θO二Oの解として得られる. そこでθZg/δ0を求めよう. 式(2.12)

を微分すると

土 sech2(仇-llf(川) -di 112) { 包 ^{- (f ( X i} ， B) -di)

^{T f}}

f(川 � ^}二O

θB \J \" ""' 2 ， θ0

(2.13)

を得る. 従ってθZg/θ8=0の解は

主 sech2(山一IIf(川)-ぬ112)(f(川) _ di)Tåf(xi，8 _。。 L ^{= 0 . (2.14)}

から求まる.通常f(x，のは0に関して線形であるから f(x，8)= f[(x)8+

fo(x)と書き表せ， θf/θ8=fl(X)となる.これを式(2.14)に代入して0 について形式的に解くと

。 = [喜悦州仇- ^Ilf(Xi， 8)一仇112])flf[(f(Xi，8)ーム)T ] ^-l

乞sech2(g[Zg- Ilf(Xi， 8) -diI12])[diーん(Xi)] (2.15)

となる.この式の右辺はO によるので反復法で解かねばならない.以上 よりGNC法によるLMedS推定は， 例えば go= 0.1，α= l.l，ξ1 -ε2 10-Ro。を任意の値として次のような手順になる.

( i )

⁹

= go

(込) 8二80

( iii ) 式(2.12)をニュートン法で解いてらを求める

(か)

仇=式(2.15) の右辺

(

^v

) 18η- 81

^{<ε1 ?}

yes =今(札)ヘ行く

η0=今O二Oη，(iii)ヘ行く (札) 18n-801くと2

^?

yes二今終了

間二今g=α*go，80=8η，(込)ヘ行く

」

2.2.2

例題

簡単な例題としてzも dも1次元として1次式すなわちf(x，の二 αx+b

の当てはめをしてみた. データはOから1までの区間に9個等間 隔にぬをとり，このうち5個のdiはdi二0.5Xiとし3 残りの4個にはそ

の値にノイズを加えた. まず最小2乗法はノイズに強く影響された. ノイ ズがすべてプラスの値のときは当然ながらbがプラスになり，正負混ざっ ているノイズのときでも一つでも大きなノイズがあるとαとbの値は正解 0.5，0から大きくずれた. 次にM推定ではαの値が最小のノイズよりも 小さければαとbは0.5と0になったが，αがこれより大きいとこの値か らずれ，十分大きなαでは最小2乗法と同じ値になった. 最後にLMedS 推定ではノイズをいろいろ変えてシュミレーションしてみたすべての場 合において正解に収束した. 収束の様子の一例を図2.5に示す. 9竺Oの ときのαとbの値は最小2乗法であり，正解から大きくずれている. なお gを大きな値から始めるとαとbの最初の値は最小2乗解からずれ，最終 結果がαとbの初期値によって変わるようになり，初期値によっては正解 でない値に収束した. 従って十分小さなgから始めなければならない.

2.2.3

むすび

代表的な非線形関数回帰法である日推定とLMedS推定の求解法とし てGNC法に基づくアナログ解法を提案し，簡単な実験によりGNC法の 有効性とLMedS推定のロバスト性を確認した.

」

0.8 0.6 0.4

\ \ \ \ \ \ \ \ 、\ \ \ \ \ \ \

a ^一一一

b

^一一一

0.2

0

←一一一一一…一一一一一一一一---/---- 一一 ーで九一一一一一一一一一一一一一一^\

\ \

ー0.2 -0.4 -0.6 -0.8

0.1 10

9

図2.5:αとbの収束の様子

2.3 最小2乗トリム回帰のアナログ解法

外れデータ(outlier)と呼ばれる大きな誤差を持つデータが頻繁に含ま れるようなニューラルネット学習[21，22]や画像処理[17]などでは，線形 な回帰法である最小2乗法に代わるロバスト推定による非線形な回帰法 が使われる. 最もよく使われているのはM推定であるが， スケール変動 に関してもロバストな方法としてRousseeuw[23]によって最小2乗メジア ン(Least Median of Squares)法や更にそれよりも漸近的な性質が優れて いる最小2乗トリム(Least Trimmed Squares)法などが提案され，これら はまだM推定ほと広くは知られていないが画像処理への応用などが試み られている[19，24]

.

しかしこれらの非線形回帰法は基本的には混合整数 計画問題であり， 非常に多くの極小解を持つ[2 5] ので， 実用的な計算時 間で近似解を求める確率解法などが提案されている[23]. しかしこれらの 近似解法でも回帰パラメータの数が増えると急速に求解が困難になって くる. そこでここではこれらのディジタル解法とは異なるアプローチと して， ホップフィールドニューラルネットの解法に似たアナログ解法を提 案する. まず最小2乗トリム(以下LTSと略す )法を混合整数計画問題と して定式化し， その近似解を決定論的アニーリング法で求める方法を提 案し，この方法によって良い近似解が得られることを理論的に傍証する.

2.3.1

LTS推定

点、Xiでの値がYiであるデータがn個与えられているとする. XiもYiも 一般にベクトルである. このデータに関数f(x，のを当てはめるとする. () は求めるべきパラメータベクトルである. このときLTS推定は次のよう に定義される[23].

m j ⁿ � ^{(γ2)i:n (2.16)}

(γ2)川はη個の残差r-;=(Yi-f(九()))2

(i

⁼

1， ・川)のうちの小さい方

からt番めのものである. すなわち残差の小さいものh個の総和を最小に するOをパラメータの推定値とする方法である.九の値は適宜設定する.

最もよく使われるのはηを奇数としてh=(η+1)/2 である.九=ηの場 合が最小2乗法である. 一般にLTS推定はNP困難な問題であり， 確率 解法[23]などで近似解が求められる. ここでは後述する解法との比較の

ため，まず次のような単純な解法を考える.

」

81: eを最小2乗法の値に初期設定する.

82:残差を求めてソーティングし， 小さいもの九個を 選ぶ.

83:選んだh個の残差の総和を最小にするOを求める.

84: eの収束をチェックしてまだ収束していなければ 82ヘ戻る.

ステップ83は九個の残差による最小2乗推定であるから， この解法は 最小2乗法とソーティングの繰り返しであり， 単純である反面局所最適

解に捕われ易いと予想される.

2.3.2

アナログ角手法

前節の解法はソーティングを含むディジタル解法であるが， ここでは それに代わるアナログ解法を考える.そのためにまず式(2.16)をそれと

等価な次の式に書き換える.

TT E ^{P バ} |

> ^(2.17)

subj.to 玄 _{Pi二九， Piε{O， 1}} I

この混合整数計画問題をアナログ解法で解くために， ニューラルネット による組み合わせ最適化の手法[26]に倣って， 式(2.17)の目的関数にエ ントロビーg=ε[pi1npi +

(1 - Pi)ln(l -

Pi)]を付け加えて次のような非 線形計画問題に書き換える.

mm Z pt ヴ +Tg

p，o t=1

subj.to 乞 Pi = h， 0

_�

Pi

_�

1 i

⁼

l

(2.18)

式(2.18)の解は温度TがOのとき式(2.17)の解になる. 式(2.18)の解は 次のラグランジュ関数の停留点である.

すなわち

L二 乞 Pir ; + Tg +入( 乞 Pi -h)

i=l

^1ニ1

δL

円

札

一一 ニ ザ+T

ln

^{-_ -'_" _} +入=0

θ'Pi

^{. � ， -}

--- 1 - Pi

(2.19)

(2.20)

ハU ハU

h 2 2

一

j 一

針 一 釘

仇

p

nZ凶 η乞凶 紅一似弘一MW

(2.21)

(2.22)

から求まる. まず式(2.20)から

Pi二五百刊)/T

l (2.23)

となる.これを式(2.21)に代入して入について解けば入は 仔(i= 1，. . .，η) の関数となり， それを式(2.23)に代入したものを式(2.22)に代入すれば，

式(2.22)はOだけに関する式になるのでこれを解けばOが求まる. なお式 (2.21)や(2.22)は解析的には解けないのでニュートン法などの数値計算で 解を求めなければならない. この解を大きなTから出発してTを徐々に 減少させながら追跡しTこOで式(2.17)の解を得るのが決定論的アニー リング法[26]と呼ばれる方法である. アニーリング法は前節の解法のよ うに直接(2.17)を解くよりも計算時間がかかるが， より良い解が得られ ることが経験的に知られている. ここではそのことの理論的な傍証を与

えておく.

[性質1 ]決定論的アニーリンク、法においてTを単調に減少させれば�Pi仔 も単調に減少する

.

(証明)f = �Pirl， 9二�[Pilnpi+ (1 - Pi)ln(l -Pi)]とおいて， 変数pと0 とを合わせてq(すなわちq = [p，

B])とすると式(2.19)は

L

=

f(q)十Tg(q) +入(Aq -

h) ^(2.24)

と書き表される.ここで

A

は[1γ・.，1，0，・リ0]という行列である.まず 一般的に次の 補助定理が成り立つ.

[補助定理]任意のベクトルsについて次の等式が成り立つ.

d f θf

_{A T} \T

dq

- =

(

^{ー +}

^A'l's)

^'1'一三

dT

^{\ 8a ' - - -/}

dT (2.25)

(補助定理の 証明)式(2.21)すなわちAq

- h

^{= 0を Tで微分すると} Ad q/dT

=

0となるので任意のベクトルsについてSTAdq/dT = 0と なり， この式とdf/dT

=

(θf/θq )T dq/dTとから式(2.25)が成り立つ.

」

この補助定理を使って性質1が以下のように証明される. まず式(2.20) と(2.22)は まとめて

θf θ9

.T'

�+T一三+

Al 入=0

θq θq

と書けるからこれを Tで微分すると

( θ

_一一

2 /

+T

^{θ2g\ d q}

一

一

一

)

一一+一一+AT--=

^θ9

^，

^， I T d入 O åq2 δq2/ dT θq dT

となる. これにがL/θq2=θ2f/θq2 +Tθ2g/θq2を代入すると

d q θ2 L

\ _ 1 I

åg

_.. T'

d入

一

=-( -) - 1(

i +

ポー)

dT θq2 θq dT

(2.26)

(2.27)

(2.28)

を得る. 式(2.26)から得られる θg/θq=-(θf/δq+AT入)/Tを式(2.28) に代入すると

dq θ2 L\ _ l θf

^I

A T ，

fT1

Æ T d入

否 ニ(

石

^�)-l(�

� ^{+A入 -}

^{T A1}

^��)/T

^(2.29)

となる. 上記の補助定理から

df

一一=

I θ (

一

一+Al入 f

， .T _\ _

T A1

rn.T'

��)1 d入 、T'd � � q d T

\θq ^{， -}-

..

- -

- dT / dT

であるからこれに式(2.29)を代入すると

df T θ2 L \ _ l

-=zl(

ー

す)

-

1 z/T

dT - \ åa 2

(2.30)

(2.31)

を得る. ただしz - θf/θq+AT入-TATd入

/dTである.

各Tで qはL の極小点、で、あるからがL/θq2は正定値行列でありまた Tは正であるから 式(2.31)よりdf/d T ^> 0 である. 従ってTが減ると f も単調に減少する.

(証明終)

この単調減少 性によって決定論的 アニーリングによって f= �Piイが小 さい解が得られる. なおTが十分 大きいときにはPi二九/n(i=l，. リη) となるのでOは最小2乗解となる. 従ってアニーリング法では初期値に よらずいつで も最小2乗解から出発することになる. 従ってアニーリン グ法では初期値の設定に関する問題はない. 一方前節の解法では得られ る結果が初期値によって大きく変わる. 適切な初期値は一般には分から ないので前節では中立的な最小2乗解を初期値とすることにした. これ は アニーリング法との比較に関しても公平である.

2.3.3

実験例

前節のアニーリング法の性能を検証するために簡単なデータで実験し てみた.Xi， Yiとも1次元とし，直線f =αx+bを当てはめた.すな わちe= [α， b]である.九二(n + 1)/2 とした.アニーリング法ではT は最初100 とし，0.9倍ずつ10-5まで減少させた.まず図2. 6 のような 11点のデータでやってみたところ同図に示すような回帰直線がそれぞれ 得られた.式(2.16)の残差は2.3.1節の方法では0.0042 であり，一方ア ニーリング法ではOとなり最適解の直線が得られた.性質1の確認のた めこのときのアニーリング法でT減少とともにf = �Pirlが減少する様 子を図2.7に示す.次に文献[22]に似たデータとしてy二0.5xという 直線に平均値0 ，標準偏差0.05のガウスノイズと平均値0.2

，標準偏差

0.1のガウスノイズとを2:1の割合で加えた図2.8のデータで回帰直線を 求めたところ同図に示す結果となり，図2. 6同様アニーリング法が最も 誤差の小さい解を与えた.式(2.16)の残差は2.3.1節の方法で0.041，ァ ニーリング法で0.0016 である.また画像処理への適用例として画像強度 の時空間微分によるオプテイカルフローの推定に応用してみる.時刻tで の第(i， j)画素でのグレイレベルをILと記す.このグレイレベルの勾配 DX;j =

(Il+1，j

-

1i�

_1，j)/2， DYÔ =

(1[，;+1

- 1;，j_1)/2 ， Dη ⁼ 勾十

1

_-4に

よる第(i ， j)画素での速度ベクトル川=(りふ吟)に関する拘束は

DXLuZ+D只;υむ+DTi� = 0

(2.32)

である.2 個の変数υふ坊に対して式(2.32)1個だけでは式が足りないの で，近傍の画素での勾配も用いてυふ坊が推定される.このとき最小2 乗法だとノイズや物体境界などで推定値が歪むのでここではLTS法を適 用してみる

.この場合式(2.16)のη(今の場合 rij)は式(2.32)の左辺であ

り，式(2.16)の2はウインドウ内の画素に対してとられる.図2.9(a)は 最小2乗法で得られたフローの例であり，一方本方法では図2.9(b)のよ うに物体の境界付近での誤差が図2.9(a)よりも小さい結果が得られた.

2.3.4

むすび

ロバスト推定法の1種であるLTS推定の近似解法としてアニーリング によるアナログ解法を提案し，誤差の小さい解が得られることを理論的

-・ーー-

に傍証し， 簡単な実験によって確認した. ニューラルネットの学習や画 像処理への応用および計算の高速化などが今後の課題である.

-・ーーー'ー

0.6

0.5

。 0.4

0.3

0.2 ---ーー<5^- ()'

0.1 .0

。

。 0.2 0.4

.0' _----ーミ> 0

0.6 0.8

。

三〉

('j

1.2

図2.6:回帰直線(実線:最小2乗法?破線:2.3.1節の方法?点線:アニーリン グ法)

-・-

0.14 0.12 0.1 0.08 0.06 0.04 0.02

0

1 e-05 0.0001 0.001 0.01 0.1 10 100 T

図2.7: T減少によりfが減る様子

[

^-・・ー-

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

0 -0.1 -0.2

。

。 。

。

0.2 0.4

。 。 。

。 。 。

。 。 。

0.6 0.8

図2.8:回帰直線(1点鎖線:元の直線?実線:最小2乗法?破線:2節の方法?点 線:アニーリング法)

lVlVlVJVJVJVJVlvjvlVlVlVlVlVJVV lvlVIVJVJVJVJVJVJVJVIVIVIVIVJVJV

l vlψlvvvvv l vlVlV l V l VlVlwuv v

方

-小

々 ペ J々々々々↓lvlvlv

JV

W/小介外介介〉hvivlVえぐ小ぶ寸ふソ市介爪小小小小小冷iγ入ヤl J ャ 匂/ア『/\小小小A小小余

l |

ヤト\人'ψ中fAド引不小小小小小引小」ヘ\ト\〈山vふ

|

爪

|

爪lF小小小小小小小ふヘ\ぺ

〉

ルーかハ庁 山V

一小小小小小

A

r

小六い\べ，vv

寸

K

/ト

一 六

_小小小

小小小爪

l

F々

M

V

l

吋\\一元

小小小小小

小↑ボホ

『 V WV併￠同小小小小小小↑司みJJV

求

、小引JVJVJVJV

l

vル

l

A寸

J V

必ヤァ小 〉 ↓VJVIVIVlVlV寸 寸ヲJ JV l羽

lV

lv

vv

vψ

ψ

ψ

lVlV lV

1V

↓ J

ヤlVlVIVしVLVLVヤヤヤヤlvlVlVIVLVしV

\レ\レ\レ\レ吋/ψ\レ\レ\レ\レ\レ\レヤ\レ吋/ψ吋/吋/

\レヤ\レ吋/ψ吋/'\レ\レ\レ\レ\レ\レ\レ\レ吋/吋/吋/吋/

(p)姉よJM潮路

lV

lV lvlVlVWJVJVJVJVJV JVIV IVlVlvv lVlVlVJVJVJVJVJVJVJVIVlVlVlvιv lv vv JVJVlvjv

lv

lVlVlV lv v

lVlV lv vvv

lv lv

lvlVlVlV

lV

↓ ル

オペlVUA寸寸vvlvlvlvlVlvムムLJV/弘hJAぺAIlVlvlVlvφ|〉\ノゾ

匂/、/ A山 小小 小 小 小小爪l爪li〉〉〉

小

/ 苛 ス小小小小小小爪 1 小 7 〉〉 AYA小小小小小介吟AIAlA1爪l小 〉 爪Y小|小小 Aド小小引↑爪l爪l爪l爪| 〉 爪Y小|冷小小ArA下小ホーホlAーベ爪l 〉

bl介IR小小小小小小小↑-\ホ

l 〉 ル\ ιr

f小小小爪小小↑ポー

ホl 〉 LW\除JK\K 小小小小 爪lAー↑↑ 〉〉 J

本

\

糸

l 、v iv j v l v l V l V泳 | fv弘〉 ー ザ\ l v>し 山 v i v J V l v l V l V l V l V l V? J L切ヤlvvvvψψψlVIVlVlV↓ヤ

lVlVlVLVLVLVLVヤヤヤヤlVlVlVLV

\レヤ\レ吋/吋/吋/'\レ\レ\レ\レ\レ\レ\レ\レ\レ\レ\レ\レ

\レ\レヤ吋/、↓\レ\レ\レ\レ\レ\レ\レ\レヤ吋/ψ\レ\レ (σ)円L叶∞路

図M・。

リ川口lニ、lい\

tv A

\レ\レヤ刈/ψ吋/'\レ\レ\レ\レ\レ\レ--.v--.v\V\V吋/'\レ

圃園田ーー-

2.4 ロバストエネルギ一最小化によるエッジ保存 平滑化

画像復元などではエッジを保存しつつノイズを除去することが要求さ れる. エッジを保存してインパルスノイズも取り除けるフィルタとして メジアンフィルタが有名であるが， ガウスノイズの平滑化能力が低いの で種々の改良が加えられている[18]. メジアンフィルタは平均フィルタを ロバスト化したものであるが， 別の流れとして正則化フィルタをロバス ト化する研究がある. 正則化フィルタはガウスフィルタとほとんど同じ 特性を持つ線形フィルタであり， インパルスノイズの除去やエッジの保 存に難がある[27]. そこで正則化にラインプロセスを導入し非線形化して エッジの保存性を高める方法が提案されており [5]， さらに同様な0，1変 数の追加によりインパルスノイズの除去も行えるように拡張されている [6]. このような拡張は正則化の基礎となる最小2乗法をロバスト推定の 1種であるM推定に拡張したことと等価である[4]・ しかしこの方法には パラメータの設定が難しいという難点がある. そこでここではM推定に 含まれるしきい値のようなパラメータを含まない方法としてLMedS推定

[23] に基づく平滑化法を提案する.

2.4.1 エネルギ一最小化による平滑化

簡単のため1次元データについて記述する. 第t番目のデータをdi， そ れからの復元値をfiとする. ここではfiをあるエネルギーの最小解とし て求める方法を考え， エネルギーとしてまず最初に最も基本的な最小2 乗法に基づくもの， 次にM推定に基づくもの3 そして最後にLMedS推定 に基づくものを考える.

2.4.2 最小2乗エネルギー

まず最初にhを次の式から決める方法を考える.

ffiln

α乞(ん- di)2

+

ßL(λ- fi+I)2 (2.33)

このエネルギーは下に凸であるから最適解は唯1個であり， 次の連立1 次方程式を解けば得られる.

」

^-ー-

α(fi - di)

+

ß (2fi - fi-l - fi+l) = ₀ (2.34)

これは最も単純な標準正則化法[27] である.

2.4.3

M推定エネルギー

式(2.33)の2乗エネルギーは変位が大きくなるほど復元力が強くなる.

従ってインパルスノイズ付近では復元値が大きく歪み3またエッジもなま る.これを改善する一つの方法は2乗関数を途中で飽和させて式(2.33)を mln α乞σα(fi- di) + ß '2二σβ(h - fi+l) (2.35) と変えるものである.特に

、、B，，〆 、、Bl''

p p・

4FU く ，TL

> 一 Z Z

J''BE、 ，，SEt、、 つ mノ つ- pa q ・4 4'u /EEEノ、BEEt、 一一 、、E，，JZ P σ (2.36)

のとき，ロバスト推定の1種類であるM推定，そのなかでもHinich-Talwar

推定[28]に対応する.らは正定数である.エネルギーの第1J頁をσαにす る効果はらよりも大きいインパルスノイズを取り除けることである.ま た第2:t頁をσpにすることによりジャンプ幅がtβよりも大きいエッジが保 存される.式(2.35)のエネルギーは式(2.33)のエネルギーにラインプロ セス[5] と呼ばれる0，1変数を付け加えて3その変数を消去しでも得られ る

[

4

].

なおエネルギーの第1項に付加されるOぅ1変数はノイズプロセス とかサンプリングプロセスとかいろいろな呼び方がある[6].式(2.35)の エネルギーは非凸なので多くの局所最適解を持ち， よい解を得るのは容

易ではないが， 一つの方法としてアニーリングがある. 例えば

仰)= 十 ^[1

⁺

eg(t�-x2)]

₊

_{t; 川}

という関数はg:::Oでx2と等価となり， 9↑∞で式(2.36)ののになる.

そこで式(2.35)のσαとσpとを式(2.37)のふとおとに置き換えた式の 最適解をg:::Oからg↑∞まで追跡する.この最適解は

αι(fi - di)

+

ß[f，ß(fi - fi-l)

+

f，ß(h - fi+l)]

⁼

0 (2.38)

の解として得られる.9竺Oのときこの解は式(2.33)の解すなわち最小2 乗解であり， 9↑∞で式(2.35)の解すなわちM推定の解の一つに収束す る.このような解追跡法は決定論的アニーリングと呼ばれ比較的よい解 が得られることが経験的に知られている.

瓦

^{-田園圃ーー-}

2.4.4

LMedS推定エネルギー

M推定による平滑化ではらやtßなどのパラメータを適切な値に設定す ればインパルスノイズも取り除けエッジも保存される. しかしこれらの 適切な値はインパルスノイズやエッジの高さに依存し， またαやFの値 によっても変わるので一般的に設定するのは困難である. そこでここで はロバスト推定の別の方法であるLMedS(Least Median _of Squares)推定 に基づくエネルギーを考える. M推定が最小2乗法の2乗を飽和関数に 変えたものであったのに対し， 最小2乗平均の平均を中央値(メジアン) に変えたものがLMedS推定である[23]. そこでまず式(2.33)をそれと等 価な次の式に変形する.

m1n α 乞 ^m

+ß � ^mean{

^{(fj -ん+1)2}

^{I j}

⁼

^{i - 1， i， i}

⁺

^1}

(2.39) ここでmeanは平均すなわち単純な総和を表す. なお簡単のためウインド

ウ幅は3とした. LMedS推定ではこれを

m1n α � ^{Pi+ ß} 乞 ^{qi ， (2.40)}

Pi二med{(fj - dj)2

I ^j

=

i - 1ぅi，i+1}

^(2.41) qi二med{(fj-fj+1)2

I ^j

=

i - 1， i， i

+

1}

(2.42) と変更する. medは 中央値を表す. 式(2.40)のエネルギーは各ウインド

ウ内の三つの変位のうち二つが小さくなることを要求し残る一つはいく ら大きくても構わないということであるから， 第1項によりインパルス ノイズが平滑化され， 第2項によりエッジが保存される. 式(2.40)もま た非凸な最 適化問題であり解くのが困難である. そこでここでもアニー

リングを用いることにする. 例えば

�

^ゆ(Pi

-(ん-dj)2)=0 J=i-1

ゆ(x)⁼ tanh(gx)

(2.43)

(2.44)

の解Piはg竺Oのとき式(2.39)のmeanすなわち平均となりg↑∞のと き式(2.41)のn1edすなわちメジアンとなる[20]. そこでgを0から徐々

-ー�

に増加させて式(2.40)をこのようにファジー化した解を追跡することに する.式(2.40)をファジーイじした解は

α _j=-iー1θfi 芝

^生

^{三 +ß} 芝 ^包二

^O

' ，-

jニi-2θfi (2.45)

の解 として得ら れる.式(2.45) のなかの偏微分は式(2.43)やqiに関する 同様な式をfi で微分することにより

弘二l二叫i-1，i)(fi - di) /θfi

む

^(仇-l，j)

(2.46) 2+1

8Pi = 2グ(九)(h - di)/乞グ(Pi，j)

θfi

(2.47)

寄 ^{= 州山+ 叫 +l，i)(fi υω岬i) 川川 )}

^{(川(仏λト日川川一イイ叫dιω}

^句

ⁱ⁾

i+2 ^(2.48)

と求まる.ここで

空訪乍云rZ

= 2M州州ゆグ釧州内

f

句

'(qiωω( qi-2，

仇ιいいいZ は叫一-ト引2 2，υ ，，ユtト川叫J-1)以)(fi - fi-I) 一→

/乞グ(qi-υ)

j=-i-3

計上=収(q日-l)(fi - fi-1) +グ(qi-1，i)(fi - fi+1)]/乞グ(qi-1，j)

j=-i-2

� ⁼ 2[cþ'(qi，i-I)(h - h-1) i+1

+引か)(fi

-

h+1)]/乞グ(qi，j )

j=-i-1

23T =W(qt+131)(ftーん) i+2

/乞グ(qi+1，j

)=-2

(2.49)

(2.50)

(2.51)

(2.52)

Pi，j二Pi

-

(ん- j)2 d (2.53)

qi，j二qi

- (ん-

fj+1)2 (2.54)

である.これらを式(2.45)に代入してfiに関して形式的に解くと

l'

_αγ1di + ß(γ2fi-1 +γ3h+1)

ん 一 川 、

(2.55)