三木成彦＊藤田徹也＊＊金久正弘＊＊＊ （昭和59年4月14日受理）

水平方向の運動による流れ画像の点拡がり関数同定について（1）

三木成彦＊藤田徹也＊＊金久正弘＊＊＊

（昭和59年4月14日受理）

en An ldentification of Point Spread Function of lmages Blurred by Object Motion in Horizontal Direction 一 Part 1

Shigehiko MIKI， Tetsuya FuJITA and Seiko KANEKu

（Received April 14， 1984）

 This paper deals with an identjfication of point spread function （PSF） of noisy images degraded by object motion in horizontal direction． First， assuming a spatial correlation of exponential fQrm and separable variables，

the image signal is modeled by a stationary process at the output of a discrete−time linear dynamical system randomiy excited．

 Next， the mot ion blur is represented by a PSF and a Kalman est imator for nQisy irnages of known PSF is derived． Finally， by examining the relation between an assumed PSF and the theoretical improvement index of the estimator， it is revealed that the gain estimation is most important in PSF identification． This result is very useful for finding a PSF identification method．

1． ま え が き

 本稿は，加法性雑音ととも運動ぼけにより劣化した画像 の点拡がり関数（Point Spread F囎ction， P S F）同定に ついて考察したものである。

 最近，加法性雑音のみならず，運動ぼけをも伴うような 流れ画像の復元問題に関する研究が盛んに行なわれてい る1）、5）。PS．F同定問題についてもいくつかの研究報告が なされている6）〜8）。しかしながら，原画像信号が既知とい

う前提があったり6），その同定法と画像推定法との関係が 明らかでない7）など，いまだ不十分である。これは，本論 文で扱うような流れ画像処理の分野における同定問題が制 御工学の分野におけるそれとはかなり異なった性質を有す るからである9）。

 筆者らは以前に，PSF未知の場合への一つの試みと して，等加速度直線運動という還動の種類はわかってい るが，加速度が未知の場合についての考察結果を報告し

た10）。ここではこれを一歩進め，直線運動の方向はわかっ ているが1），PSF未知なる場合を取り扱うことにする。

特に直線運動の方向を水平方向とする。本稿では，流れ画 像の復元を目的とした同定を扱っているため，同定の良さ を表わす指標として，その同定されたPSFを用いたとき の推定画像の改善度を用いる。そこでまず，画像信号と観 測信号のモデルについて述べ，次に画像信号を逐次的に推 定するカルマン・フィルタ10）を導出する。最後に，仮定さ れたPSFとそれを用いてカルマン・フィルタにより信号 推定を行なった時に得られる理論的改善度の関係を調べ，

PSF同定における最も重要なことはPSFのゲイン推定 であることを明らかにする。

＊ 津山工業高等専門学校電気工学科

＊＊読売テレビ放送株式会社

．＊＊＊福岡工業大学

  2．モデルの設定およびカルマン・フィルタ

2．1 画像信号モデル11）

対象とする．2次元画像をMXNの画素にディジタル化 し，i行ブ列聖の画素の明るさ（画像信号）をfi，ゴとすれ ば（Fig．1参照）

ノi．ゴ＝Xi，ゴ＋μ， （1）

津山高専紀要第22号（1984）

12：．⁝二・⁝−−ーー−凹

z 2 一一一 j 一一一一一一一一一一 r 一一一 N

fLj

Fig．1 Discrete representation of two−dimensional    image（ルf＝ハr＝嵩10）．

となる。ここで，μは信号ノ垣の平均値であり，毎ゴは次 のような平均値と自己相関関数を持つ広義定常確率過程と

する12）。

E ［xi， i］ ＝ O

R（為，1）＝E［陶， 露 ＋lt， ＋」］

       lil lkl

   ＝a2x ph pv ， O〈ph， pv〈1

（2）

（3）

 ここで，E［・］は集合平均を表わし，σ2xは毎ゴの分

翫．oあ． o。． H；r−do．Ain．7k llZ ilrtfl…亜商凸「片rfi1 lr9（S ．． t・ ur N幽 Iu、r   昌 rr■    v  ■v− 、 ・■ U  、 「1 Vrぜ  L 》 「 →ラ  山Pり 1−4 ■−1回  ノ ，｝Tua

距離に対する相関係数である。又，この式はテレビ画像な どについて実験的にもよく検証されたものである。このよ うに本稿では胸部X線画像処理や心臓画像処理などのよう に前もって対象画像が定まっている場合ではなく，上記の 統計的性質以外の情報が前もって与えられていない場合を 取り扱うことにする。

 このとき，ノ珂をそのi行目の要素として持つM次元列 ベクトルF とすれば（Fig．1参照），次のような画像モデ ル式を得る。

Xl・＋1＝ph Xj＋Vlrj， 1 ＝O，・・・…，N−1 Fゴ＝xノ十Ω，ノ＝1，……，ハr

（4）

（5）

 ここで，XノおよびΩはそれぞれ毎」およびμをi行目 の要素として持つM次元列ベクトル，Wjはステップノに おけるいわゆる システム雑音 を示すM次元列ベクトル である。又，Xゴの初期値XOは平均値OM×1（M×1零行 列），共分散行列crx2 Bなるランダム変数Wiは平均値 eM・1，共分散行列σu2 Bなる白色ガゥス雑音であり， XO とWjは無相関である。更に， BはM次元正方行列で， B の（r，t） bl素br，tは次式で示されるものである。

br，t＝pvir−tl ^（6）

σ麗2二（1−Ph2）σx2 ^（7）

 このように，このモデルでは画像信号の水平および垂直 方向の相関情報が，それぞれ（4）式のXノの係数およびWi の共分散行列の中に含まれている。

 なお，一般性に固執するならば，PSFと共に上記の画 像に関する統計量およびパラメータはすべて未知13）である との前提に立つべきであろうが，それはいたずらに問題を 複雑にして本質を見失う恐れなしとしない。従ってここで は，PSF以外はすべて既知であると仮定しておく。

 2．2観測モデル

 いま，簡単化のため，カメラのイメージングシステムと 対象物体間の相対運動により水平方向にのみ生じた空間不 変の劣化を考え，（∫，の要素の観測値夕珂は次式で表わ

されると仮定する。

  」 i，ゴ＝Σht fi，ゴ＋1÷砂ゴぜ，1≦ゴ〈M，

     1＝一m

       m一』1≦ニノ≦；N−n       （8）

 ここで，％ゴは加法性ランダム雑音でv平均値0，分散 σv2なるガウス雑音とする。又， htは空間的な劣化を表わ す点拡がり関数（PSF）である。更に，m，π＞0はぼ けの範囲に対応し，例えば物体に対してカメラが左から右 レ舌hいチ7十巨A〜斗初＝nr7s ｝IA  S・i「Lつア善』紳｝一ジ…一AL          ULVJソ，！国〃^冒一訓冒 ワ1：9LSYJV／Vd・ 御mU（」

    」吻口Vd噛−v なる。

 錫ゴおよび 毎をi行目の要素として持つM次元列ベク トルをそれぞれyノおよびγノとすれば（8）式からy，は次 のようになる。

     の

  】rゴ＝ Σ hlF ＋1十V， ，  m＋1≦（ブ≦1＞一n   （9）

    1・＝一Lm

 ここで，玲は平均値◎M、1，共分散行列σv2 IMの白色 ガゥス雑音で，WゴおよびXoとは無相関とする。但し， IM はM次元単位行列である。

 2．3 カルマン・フィルタ

 1（4）式，（5）式の信号モデルおよび（9）式の観測モデル をカルマン・フィルタ適用可能な形に状態拡張すれば，次 式のようになる。

Sゴ＋1＝ΦS 十σゴ」 m≦ニブ≦コV−n      （10）

yノ＝H（S，十C）十γン，  〃z十1≦望ブ〈ム「一〃  （11）

但し，

Sl ＝［X」LmT Xf；一tn＋IT Xs Xs ＋．T］T

H＝［h−mlM h−m＋11M ho lM hn Inf］

ひ芦［01・MO1×M…◎1・MWゴT］T c＝［nT aT ．Tn］T

C12）

（13）

（14）

（15）

更に又，σu2は次式により与えられる。

一2 一

￠ ＝＝

O・M

IM  Oル「一一一一一一一一一一…一一一一一一一◎M      ロへ             サロ      ま

◎M  IM  ＼＼    i

         、、       1   、、        、、      、、         1

      ＼、 ＼＼  ◎M          ㊥M   IM

eM一一一一一一一一一一一一一一一一一一一一一一一一一eM phJM

（16）

 上月中，Tおよび㊥Mはそれぞれ行列の転置およびM 次元正方零行列を表わす。又 σゴは白色であり玲とは無 相関である。更に，Sゴの初期値Smについては

 E［Sm］＝＝OM（m＋n＋Dxl （i7）

      B

E［SmSmT］ ＝ a．2

    PhB一一一一一一一一一phm＋nB

［絶・認＼！

 F  x．． B phB

phm＋nB一一一一一一一一一一一一一一一一］phB B

（18）

 更に又，（10）式，（11）式のモデル式は可制御性および可 観測性を満足している。従って，（10）式，（11）式にカルマ

ン・フィルタを適用することができる。

 さて，このカルマン・フィルタは，信号モデル式が画像 信号をよく表わしている限りMSE（Mean Square Error）

最小の意味で最適である。しかしながら，Mは通常かなり 大きな値になるので，上式のままでは計算時間や記憶容量 が大きくなりすぎ実用的でない。従って，Mが1の任意の 倍数であると仮定して画像を1行N列のストリップ毎に処 理する1次元ストリップ処理法（1次元ベクトル処理法）

が報告されている14）〜16）。このストリップ処理法は1＝1の 場合は特にズカラ処理法と呼ばれる。このストリップ処理 法は垂直方向の相関情報の損失を伴う。この損失を償う方 法として，Overlap処理法15）， Co皿posite filtering法16）な どが報告されている。

 なお，このMSE基準における最適が必ずしも画像処理 における最適ではない17）ことは事実であるが，統計的性質 のみが与えられている場合，このMES以外に適当な基準 がなく，この基準によっても一定の改善度が得られること も事実である10）。

3．仮のPSFと理論的改善度

 本章では，前章で示したカルマン・フィルタを用いて画 像を復元する場合に対する仮のPSFと理論的改善度の関 係について考察する。但し，（9）式のmおよびnは既知と する。又簡単化のため，以下ではスカラ処理法を採用 し，n＝Oとして議論を行なう。このとき，観測値ッ垣に 対する処理モデルは，正確なモデル（10）式〜（18）式に対応

してそれぞれ次のようになる。

Sノ＋1＝ΦSノ十σノ，  m≦＜ブ≦：1＞

y房躍H（sゴ＋c）＋ 毎

Sl 一一［xi，」・一m xi，」・一m＋1．．． xi，1・］T

H＝［h＿mh＿僻1…ho］

σゴー［00…OUi，ゴ］T

C＝［ lb IL 一・ pb ］T

rp ＝＝

Q 1． 0一一一一一一一一一一一・一一・一9

  0  t． SNN

  o

O  l

O一一一一一一一一一一一一一一〇 Ph

E ［Sml ＝＝e（m＋1） xl

E ［Sm SmT］ ＝ cr2x

1 Ph．…一一一一……一一一一一一一Phm

    、       1     ヒ      ぞ      モへ      

Rh 1  ＼、    i  1＼＼＼ ＼『  l

i＼ミ肥

汁解一一一一一一一一一二初1

（19）

（20）

（21）

（22）

（23）

（24）

（25）

（26）

（27）

 但し，賜ゴは平均値0，分散σu2なる白色ガウス雑音で 賜ゴとは無相関である。

 Hが未知なる場合の上記モデルに対する逐次推定式は次 のようになる。

 ム         へ

 Sゴゆゴ＝Sゴ1，ヴ＿1十K／Aゴ， 吻十1≦1ブ≦＝N    （28）

 Kゴ需Pゴ1ゴー1HT［HPゴーゴーIHT十σv2］一1，m十1≦＝ゴ≦；N（29）

Pゴ1ゴ＝Pゴげ一1−KブHPゴ・i−1，

P ヴー1＝Φ1）ゴー1，ノー1ΦT十Φu，

       ム Aゴ＝（夕i，」一HC）一HSゴli，ノ＿1，

ハ      へ sノ［i，ノー1＝・ΦSJLI I i，ノー1

Ou＝

0…一一一一一一一一…一一一一一・一…O I＼、      i

I    、、      ［ 1      、、       ［ 1         、、      I l      、、      L I       、       「 1       、       1        ヘコ       ロ

1      0   0 i

O一一…一一一一一一〇 σ 2

m＋1s｛ jfE｛ N m≦ニノ≦＝N

m＋1≦（ノ≦N

（30）

（31）

（32）

（33）

（34）

     ム      ハ

 ここで，Sゴi ，ゴおよびSゴr房＿1はそれぞれッ房を観測 したときのSブのフィルタリング推定値および一段予測値 を意味している。又 馬 はカルマン・ゲインと呼ばれる ものであり，Pゴげはフィルタリング推定値の誤差共分散 行列，Pノげ＿1は一段予測推定値の誤差共分散行列である。

 なお，（29）式〜（32）式において，普通のカルマン・フィ

津山高専紀要第22号（1984）

ルタの場合のHの代りに仮のPSFH＝［h＿m h＿m＋1……ho］

が用いられているのは，Hが未知であるからである。

 さて，Xi，ゴの推定値2i，ゴの平均2乗誤差MSEは特定の 有効桁数に対し，離散時間to後にある一定の値に収束す る。そこで，この収束状態での改善度をしらべるために次 のようなパラメータを定義する。

   Var｛h＿〃3ノンーm十……十ho∫ゴ｝

         σ2S／N＝

       ハMSEx＝E［（κ，ヴーXi，ゴ）2］， 1≦9 i≦9M，  0〈ノ≦N MSEア＝E［ ヴ2］＝σv2， 1〈i≦＝M， to≦9ゴ≦91＞

η＝1010910（MSEx／MSEン） （dB）

G＝h＿m＋………十ho

（35）

（36）

（37）

（38）

（39）

 ここで，Var｛・｝およびηはそれぞれ分散および理論的 改善度を表わし，GはPSFのゲインである。いま， h＿1

＝＝O．3， ho＝O．8， S／N＝100， Ph＝O．974， p＝s8．5， ax2

＝440，有効桁数7に対してto＝9とする。これらのデー タの場合の，仮定したPSFと信号推定の理論的改善度の 関係をFig．2に示す。 Fig．2を見れば，処理結果がPSF

10

G−1，1 G−1，125 G冨1．075

       バ

Aゴ＝（HSj−HSノ「 ヴー1）十（G−G）μ十 ヴ

但し

G＝h−m＋h−m＋1十 ＋ho

5

︵曽︶置

G冨1．15 G−1，05

（h−i G−ho）

   O O，5 ・ 1，0 一

       ho

Fig．2 Assumed PSF ho， h．1 vs． theoretical improvernent    index n （h−1＝O． 3，ho＝O． 8， S／N＝100， Ph＝O． 974，

   p・ ＝＝ 58． 5， ax2＝＝ 440， to＝9 for seven signifieant    figures of MSE）．

のゲイン推定の誤りに大きく影響されていることがわか る。従って，PSF同定における最も重要なことはPSF のゲイン推定であることがわかる。いま，カルマン・フィ ルタの式からこのことを確かめてみよう。（28）式より，

      バ       ム

フィルタリング推定二一S川，ノは一段予測値Si［垣＿1に観 測値ア財による修正項K／Aノが加わったものである。こ のK／Aノのうち，ア耐を直接含んでいるのはAゴであるの で，Aノを調べてみよう。（20），（22），（24），（32）式より

（40）

（41）

 Sノの各要素の標準偏差が21であるのに対しμが58，5で あること，および上式の第1項が真のPSFHと仮定し たPSFHの個々の値に関係しているのに対し，第2項 はHとHの個々の値ではなく，それらのゲインGとGの差 に関係していることより，第2項の方が第1項よりAノに 対する影響が大きく，かつGがAゴに大きく影響すること がわかる。即ち，Fig．2の結果がカルマン・フィルタの式 からも納得のいくことがわかる。一般的に平均値μはFig．2 の場合のようにかなり大きな値になるので，Fig．2から得 られたこの結果には一般性があるものと思われる。ちなみ に，Fig．3の原画像から計算：されたμおよびax2は58．5お よび440である。なお，ここではスカラ処理の場合を扱っ たが，これを多次元ベクトル処理にすれば，ηはさらに向 上するものと考えられる。

Fig．3 Original image（104×130）．

4．あ と が き

 本稿では，加法性雑音とともに水平方向の運動ぼけによ り劣化した画像のPSF同定について述べた。

 まず，信号モデルおよび観測モデルを示した。次に，画 像信号推定のためのカルマン・フィルタを導出した。最後 に，仮定したPSFに対するカルマン・フィルタの理論的 改善度を調べることにより，PSFのゲインを正しく推定 することが，PSF同定には最も重要であることを明らか にした。この結果に基づいたPSF同定法18）や2次元方向 の流れ画像のPSF同定19）については別の機会に発表する 予定である。

 最後に，日頃いろいろ御討議頂く京都工芸繊維大学樋口 助教授大阪大学田村助教授および大阪電気通信大学富永 助教授に深謝致します。

一4一

         文     献

1 ） A． O． Abeutalib and L． M． Silverman ； IEEE Trans．，

  CAS−22−3（1975−3），278．

2） A． O． Aboutalib， M． S． Murphy and L． M． Silverman；

  IEEE Trans．， AC−22−3（1977−6），294．

3 ） J． W． Woods and V． K． lngle ； IEEE Trans．， ASSP−29−2   （1981−4），188．

4）下野哲雄 外3名；電子通信学会論文誌（D），」63一

  互声一11 （昭55−11），939。

5）片山徹；システムと制御1，26−8（昭57一一8），533．

6 ） M． P． Ekstrom ； IEEE Trans．， C−22−4（1973一一4）， 322．

7） M．S． Murphy and L． M． Silverrnan； Proc． 1978 IEEE   Conf． on Decision and Control （San Diego， CA），

  （1979−1），783．

8）片山徹，田中雅博；システムと制御，25−11（昭56−

  1Z），709．

9）三木成彦外2名；電子通信学会技術研究報告，PRL

  80−77 （1981−Ol），63．

10）三木成彦外2名；電子通信学会論文誌（D），」65・一1）一   9（昭57−9），1137．

11）金久正弘．三木成彦；日経エレクトロニクス，211（昭   54−4），106．

12） L． E． Franks；Bell Syst．Tech． J．， 45−4（1966−4），609．

13）三木成彦外3名；電子通信学会論文誌（D），」67−D−4   （昭59−4），512．

14）N．E． Nahi and C． A． Franco；】：EEE Trans。， CO〕畷一21   （1973−4），305．

15） J． W． Woods and C． H． Radewan ； EEE Trans．， IT一一23−

  4（1977−7），473．

16） S． Miki and S．Kaneku ； Proc． 2nd Scandinavian Conf．

  on lrnage Analysis（Helsinki）， （1981−6）， 248．

17） K S． Fu and A． Rosenfeld ； IEEE Trans．， C−25−12   （1976−12），1336．

18）三木成彦外2名；電子通信学会技術研究報告，IE 81−

  84（1981−10），69．

19）三木成彦 外2名；電子通信学会情報システム部門全   国大会，No．128（1983）．