JP4568730B2 - 劣化情報復元方法と復元装置 - Google Patents

Description

本発明は劣化情報の復元に関する。詳しくは、伝達することによって劣化した情報から原情報を復元する方法と装置に関する。

伝達することによって劣化した情報から原情報を正確に推定したいという要望が存在する。このような場合には、一般に伝達系の伝達関数を用いた推定が用いられる。伝達系の伝達関数が既知であって、原情報の全ての周波数領域に対して伝達関数が０でない値を持つ場合には、伝達関数の逆フィルタを用いて、伝達後の情報から原情報を完全に復元することができる。

上記した方法による復元は、伝達関数がある周波数領域で０となる場合には、そのまま適用することができない。上記の逆フィルタは、伝達関数が０となる周波数領域において定義することができないためである。従って、このような伝達系では、伝達関数が０となる周波数領域の情報が、伝達の過程で失われる。すなわち、情報を伝達する場合、その情報はもとのままの状態で伝わることは無く、伝達の過程で劣化する。例えば光学系の像を伝達する場合には、その像は収差や機器の誤差に起因して細かな微細部分つまり空間周波数の高い部分がカットされて劣化し、劣化した状態で認識される。

上記のように、伝達の過程で特定の周波数領域が除去されて劣化した情報に関して、伝達系の伝達特性と劣化後の情報に基づいて復元する技術が提案されている。劣化した情報を復元する技術は、主に画像処理の分野において提案されている。劣化した画像を復元する技術としては、例えば非特許文献１や２に記載の、Ｒｉｃｈａｒｄｓｏｎ―Ｌｕｃｙアルゴリズムを用いる方法が知られている。
Ｒｉｃｈａｒｄｓｏｎ−Ｌｕｃｙアルゴリズムを用いる方法では、画像における光の結像を一つの事象として捉え、確率統計の技術分野で用いられている手法を用いて、原画像を復元する。Ｒｉｃｈａｒｄｓｏｎ−Ｌｕｃｙアルゴリズムを用いる方法では、原画像における照度の分布に関して、その分布を正規化して、原画像における光の結像という事象の確率密度関数の分布として捉える。また劣化画像における照度の分布に関して、その分布を正規化して、劣化画像における光の結像という事象の確率密度関数の分布として捉える。上記のように考えると、光学系の伝達特性である点像強度分布関数（ＰＳＦ；Point Spread Function）は、原画像に光の点が結像しているという条件の下で劣化画像に光が結像する確率の分布を示す、条件付確率の確率密度関数の分布として捉えることができる。Ｒｉｃｈａｒｄｓｏｎ−Ｌｕｃｙアルゴリズムを用いる方法では、劣化画像の分布とＰＳＦの分布のそれぞれを用いて、Ｂａｙｅｓの理論に基づき、劣化画像の分布を実現する原画像の分布のうちで最大に尤もらしい分布を、反復計算によって推定する。ＰＳＦの分布は、例えば光学系のパラメータから計算によって算出することもできるし、実際に点像を伝達した場合の画像の分布を実験より求めて算出することもできる。
上記のＲｉｃｈａｒｄｓｏｎ−Ｌｕｃｙアルゴリズムを用いる方法は劣化した画像を復元する方法であるが、同様の方法を用いることによって、例えば電位の時間履歴といった他の種類の情報についても、劣化後の情報から原情報を復元することが可能である。
Ｗ．Ｈ．リチャードソン（W. H. Richardson），「ベイズ理論に基づく画像復元の反復計算方法」（"Bayesian-based iterative method of image restoration"），アメリカ光学会誌（Journal of Optical Society of America），（米国），１９７２年，６２巻，ｐ５５−５９ Ｌ．Ｂ．ルーシー（L. B. Lucy），「観測される分布を修正する反復計算手法」（"An iterative technique for the rectification of observed distributions"），アストロノミカル・ジャーナル（Astronomical Journal），（米国），１９７４年，７９巻，ｐ７４５−７５４

Ｒｉｃｈａｒｄｓｏｎ−Ｌｕｃｙアルゴリズムを用いると、劣化した画像はかなり良く復元される。原画像に文字が表示されており、劣化画像ではその文字をまったく認識できないほど画像が劣化していても、Ｒｉｃｈａｒｄｓｏｎ−Ｌｕｃｙアルゴリズムを用いて推定された原画像ではなんとかその文字を判読できる程度には、復元することができる。

上記のように劣化画像から原画像を復元するためには、光学系の伝達特性を正確に評価することが必要とされる。原画像を復元する過程では、光学系の伝達特性に基づいた推定を実施するため、光学系の伝達特性が不正確である場合には、復元される原画像も不正確なものとなってしまう。劣化画像に基づく原画像の復元を行う場合、光学系の伝達特性に関する正確な評価が必要不可欠である。
しかしながら、光学系の伝達特性を正確に評価することが困難な場合がある。光学系の伝達特性は、レンズの開口数、照明波長、収差、デフォーカス等のパラメータに基づいて算出される。これらのパラメータのうち幾つかが不明な場合や、幾つかが不正確な場合には、算出される光学系の伝達特性は不正確なものとなってしまう。このように光学系の伝達特性が不正確な場合には、上述したような従来の画像復元方法を用いて正確に原画像を復元することはできない。原画像の復元に先立ち光学系の伝達特性の正確な評価が必要となる。
また、光学系の伝達特性がまったく判明せず、劣化画像のみに基づいて原画像を復元したい場合がある。このような場合、従来の画像復元方法を用いることができず、原画像を復元することができなかった。

本発明では上記課題を解決する。本発明では劣化情報の分布と原情報の分布とに基づいて、反復計算を実施することによって、伝達系の伝達特性を正確に推定することが可能な技術を提供する。
また本発明では劣化情報の分布のみに基づいて、反復計算を実施することによって、伝達系の伝達特性と、原情報の分布の双方を正確に推定することが可能な技術を提供する。

本発明の一つの方法は、原情報と劣化情報から伝達系の周波数応答の分布を推定する方法である。この方法は、劣化情報の分布を特定する工程と、原情報の分布のスペクトル分布を特定する工程と、伝達系のインパルス応答の最初の推定分布を特定する工程を備えている。さらに、（１）インパルス応答の推定分布をフーリエ変換して第１の関数を得る工程と、（２）前記第１の関数に前記原情報の分布のスペクトル分布を乗じて第２の関数を得る工程と、（３）前記第２の関数を逆フーリエ変換して第３の関数を得る工程と、（４）前記劣化情報の分布を前記第３の関数で除して第４の関数を得る工程と、（５）前記第４の関数をフーリエ変換して第５の関数を得る工程と、（６）前記第５の関数に前記原情報の分布のスペクトル分布の反転関数を乗じて第６の関数を得る工程と、（７）前記第６の関数を逆フーリエ変換して第７の関数を得る工程と、（８）前記インパルス応答の推定分布に前記第７の関数を乗じて、インパルス応答の次の推定分布を得る工程を１サイクルとする工程群を備えている。さらに、前記のサイクルの工程（８）で得られたインパルス応答の次の推定分布をインパルス応答の推定分布に置き換えて、前記（１）から（８）の工程を繰返す工程と、インパルス応答の推定分布をフーリエ変換して、伝達系の周波数応答の分布を得て出力する工程とを備えている。
本明細書でインパルス応答とは、原情報の分布が単位インパルス関数で与えられる場合の劣化情報の分布を示す。情報として画像を扱う場合には、インパルス応答とはＰＳＦの分布のことを指す。
本明細書で周波数応答とは、原情報の分布のスペクトル分布に対する、劣化情報の分布のスペクトル分布の比を示す。情報として画像を扱う場合には、周波数応答とは光学的伝達関数（ＯＴＦ；Optical Transfer Function）の分布のことを指す。

本発明の方法で扱う情報としては、画像に限らず、電気信号の時間履歴などについても適用が可能であるが、以下では情報として画像を扱う場合を取り上げて、本発明の方法の原理について説明する。
本発明の方法を説明するために、まず最初にＢａｙｅｓの理論に基づく画像の復元方法について説明する。以下では、原画像は白黒の画像であり、ある光学系を介して原画像が伝達して、白黒の劣化画像が形成される場合について説明する。以下では、原画像と劣化画像はそれぞれの大きさが同一であり、画像内の点の位置を座標（ｘ、ｙ）で表現することが可能であって、原画像の照度分布はｆ_ｒ（ｘ、ｙ）、劣化画像の照度分布はｇ_ｒ（ｘ、ｙ）で表現されるものとする。

上記の原画像の分布ｆ_ｒ（ｘ、ｙ）とｇ_ｒ（ｘ、ｙ）は、以下に示す二次元フーリエ変換を実施することによって、ｘに関する空間周波数ｓと、ｙに関する空間周波数ｔについてのスペクトルを得ることができる。

上記の光学系のＯＴＦは、原画像の分布ｆ_ｒ（ｘ、ｙ）の空間スペクトルＦ_ｒ（ｓ、ｔ）、劣化画像の分布ｇ_ｒ（ｘ、ｙ）の空間スペクトルＧ_ｒ（ｓ、ｔ）に対して、次の関係を満たす複素関数Ｈ_ｒ（ｓ、ｔ）のことをいう。

一般的な光学系の場合、どのような原画像の分布ｆ_ｒ（ｘ、ｙ）に対しても、光学系と撮影条件によって決まるＯＴＦであるＨ_ｒ（ｓ、ｔ）を用いて、劣化後の画像ｇ_ｒ（ｘ、ｙ）を得ることができる。

複素関数であるＯＴＦは、複素振幅の大きさを表す振幅伝達関数（ＭＴＦ；Modulation Transfer Function）Ｍ_ｔ（ｓ、ｔ）と、位相ズレを表す位相伝達関数（ＰＴＦ；Phase Transfer Function）Ｐ_ｔ（ｓ、ｔ）を用いて、次式で表現される。

上記のＯＴＦを用いることによって、光学系の位相に関する特性についても、正確に評価することができる。上記のＯＴＦは、光学系の特性パラメータから計算することができる。

本発明の方法では、原画像の分布と劣化画像の分布を確率密度関数として扱い、Ｂａｙｅｓの理論に基づいて原画像の推定を行う。上記で設定された原画像分布ｆ_ｒ（ｘ、ｙ）と劣化画像分布ｇ_ｒ（ｘ、ｙ）は、下記の正規化を行うことによって、確率密度関数として扱うことが可能になる。

上記の正規化に合わせて、光学的伝達関数Ｈ_ｒ（ｓ、ｔ）についても正規化する。Ｈ_ｒ（ｓ、ｔ）は、空間周波数が０の点における値を用いて正規化する。

正規化された原画像の分布ｆ（ｘ、ｙ）と、劣化画像の分布ｇ（ｘ、ｙ）は、非負の関数であり、定義された領域での積分値が１であるから、確率密度関数として扱うことができる。上記の場合において、ｆ（ｘ、ｙ）は原画像の座標（ｘ、ｙ）における結像という事象の確率密度関数である。ｇ（ｘ、ｙ）は劣化画像の座標（ｘ、ｙ）における結像という事象の確率密度関数である。

原画像および劣化画像の分布を、確率密度関数と見なすことが可能な場合、Ｂａｙｅｓの理論に基づいて、劣化画像の分布から、その劣化画像を生じさせている原画像の分布を推定することができる。
原画像の座標（ｘ_１、ｙ_１）において点光源が存在する事象をＶ（ｘ_１、ｙ_１）、劣化画像の座標（ｘ_２、ｙ_２）において点像が結像する事象をＡ（ｘ_２、ｙ_２）とした場合、それぞれの事象の確率Ｐ（Ｖ）およびＰ（Ａ）は以下で表現される。

また、原画像の座標（ｘ_１、ｙ_１）に点光源が存在している場合に、劣化画像の座標（ｘ_２、ｙ_２）に結像する確率は、事象Ｖ（ｘ_１、ｙ_１）の生起を条件とする事象Ａ（ｘ_２、ｙ_２）の生起確率である。上記の確率は光学系のＰＳＦであるｈ（ｘ、ｙ）を用いて以下で表現される。

Ｂａｙｅｓの理論に基づくと、劣化画像内の点（ｘ_２、ｙ_２）に点像を結像させる場合の原画像の分布Ｐ（Ｖ（ｘ、ｙ）｜Ａ（ｘ_２、ｙ_２））は、以下で推定される。

上式の右辺のＰ（Ｖ（ｘ、ｙ））、Ｐ（Ａ（ｘ_２、ｙ_２）｜Ｖ（ｘ、ｙ））に、式（８）および式（１０）を代入すると、次式を得る。

上式の左辺は、劣化画像において点像が結像している場合に、推定される原画像の分布を表している。上式に劣化画像の分布ｇ（ｘ、ｙ）を乗じて積分することによって、劣化画像の分布ｇ（ｘ、ｙ）を実現するための原画像の分布ｆ（ｘ、ｙ）を得ることができる。
上式の両辺に、Ｐ（Ａ（ｘ_２、ｙ_２））＝ｇ（ｘ_２、ｙ_２）を乗じて、すべての（ｘ_２、ｙ_２）に関して積分すると、次式が得られる。

上式の左辺に式（９）を代入すると、その積分の結果はＰ（Ｖ（ｘ、ｙ））であり、ｆ（ｘ、ｙ）に等しい。従って、Ｂａｙｅｓの理論に基づくと、以下の関係が成り立つ。

上記の関係は、分布ｆ（ｘ、ｙ）が真の原画像の分布である場合に成立すると考えられる。すなわち、上式を満たす分布ｆ（ｘ、ｙ）を算出することが、劣化画像の復元に相当する。
上記の関係を満たす分布ｆ（ｘ、ｙ）は、式（１４）の右辺のｆ（ｘ、ｙ）をｆ_ｋ（ｘ、ｙ）とし、式（１４）の左辺のｆ（ｘ、ｙ）をｆ_ｋ＋１（ｘ、ｙ）として、ｆ_ｋ（ｘ、ｙ）に関する反復計算を実施し、ｆ_ｋ（ｘ、ｙ）の収束値を求めることで算出することができる。上記の反復計算によって求まるｆ_ｋ（ｘ、ｙ）の収束値は、Ｂａｙｅｓの理論に基づく原画像の推定分布に相当する。

上記では、原画像の分布ｆ（ｘ、ｙ）や劣化画像の分布ｇ（ｘ、ｙ）を正規化した場合について説明したが、実際の反復計算を実施する上では、これらの分布は正規化することなく、そのまま使用することができる。

上記の反復計算においては、反復計算を実施する前に、原画像の最初の推定分布ｆ_０（ｘ、ｙ）を設定しておく。最初の推定分布ｆ_０（ｘ、ｙ）としては、任意の分布を設定することができる。一般的には、劣化画像の分布ｇ（ｘ、ｙ）は原画像の分布ｆ（ｘ、ｙ）から大きく異なることはないため、最初の推定分布ｆ_０（ｘ、ｙ）としては、劣化画像の分布ｇ（ｘ、ｙ）を用いることが好ましい。

式（１４）の右辺は、ＰＳＦであるｈ（ｘ、ｙ）をもちいた畳み込み積分を備えている。一般に、光学系の位相特性まで含めてＰＳＦを正確に評価することは困難であり、上記の反復計算を正確な位相特性を含めて実施することは困難である。正確な位相特性を含まないＰＳＦを用いた反復計算は、誤った収束の結果をもたらすため、原画像の正確な復元の妨げとなる。
そこで、ＰＳＦの代わりに正確な位相特性を含ませることが容易であるＯＴＦを用いることで、より正確に原画像を復元することが可能となる。また、復元の過程で位相特性を正確に評価するために、原画像の推定分布ｆ_ｋ（ｘ、ｙ）（ｋ＝０，１，２，・・・）を複素関数に拡張し、その実部が画像分布を表現するものとして取扱う。上式の右辺に、フーリエ変換と逆フーリエ変換を用いることで、ＯＴＦを用いた形式に変更することができる。フーリエ変換をＦＴ（）、逆フーリエ変換をＦＴ^−１（）で表現すると、式（１４）は以下で表現される。

上記の反復計算を、ｆ_ｋが収束するまで繰り返し実施することによって、原画像を推定することができる。ｆ_ｋが収束したか否かの判定は、例えば反復計算の回数を予め設定しておいて、反復計算の回数によって判定してもよいし、ｆ_ｋとｆ_ｋ＋１の差分を算出して、算出される差分の絶対値の総和があるしきい値以下となるか否かで判断してもよい。

上記した原理を用いる復元方法は、以下に示す工程を順次実施していくことで実現される。以下では繰り返し計算によって推定される原画像の分布をｆ_ｋ（ｋ＝０、１、・・・）とする。原画像の推定分布ｆ_ｋは、復元の過程で位相に関する特性を正確に評価するために、複素関数として扱う。
まず伝達系の特性に基づいて、ＯＴＦであるＨ（ｓ、ｔ）を特定する。
次に原画像の最初の推定分布として、ｆ_０（ｘ、ｙ）の実部をｇ（ｘ、ｙ）として、ｆ_０（ｘ、ｙ）の虚部を０とする。
次に以下に示す演算を、ｆ_ｋ（ｘ、ｙ）が収束するまで繰り返し実施する。ここでＦＴ（）は二次元のフーリエ変換を表し、ＦＴ^−１（）は二次元の逆フーリエ変換を表す。またＨ^＃（ｓ、ｔ）は、Ｈ（ｓ、ｔ）の反転関数であり、Ｈ^＃（ｓ、ｔ）＝Ｈ（−ｓ、−ｔ）である。

上記の反復計算を繰返し実施し、最終的な原画像の推定分布ｆ_ｋを得る。最終的に得られる原画像の推定分布ｆ_ｋ（ｘ、ｙ）の実部が、原画像の復元画像ｆ（ｘ、ｙ）に相当する。

上記の方法では、原画像から劣化画像への伝達特性として、正確な位相特性を考慮したＯＴＦを用いて反復計算を実施することができる。従って、Ｒｉｃｈａｒｄｓｏｎ−Ｌｕｃｙアルゴリズムを用いる場合にくらべて、より正確な推定を実施することが可能と考えられる。

上記の方法では、畳み込み積分を用いることなく、フーリエ変換と、逆フーリエ変換と、四則演算を用いて、原画像を推定することができる。このため、Ｒｉｃｈａｒｄｓｏｎ−Ｌｕｃｙアルゴリズムを用いる場合にくらべて、処理に要する時間を大幅に短縮することが可能となる。

上記した方法を用いることで、既知である劣化画像分布とＯＴＦ分布とに基づいて、未知である原画像分布を復元することができる。上記と同様の原理によって、ＯＴＦ分布が未知である場合に、既知である劣化画像分布と原画像分布とに基づいて、ＯＴＦ分布を復元することが可能となる。

式（１２）において、ｘ_２−ｘを新たにｘ’とし、ｙ_２−ｙを新たにｙ’とすると、次式が得られる。

上式の両辺にＰ（Ａ（ｘ_２、ｙ_２））を乗じる。

上式の左辺は、Ｂａｙｅｓの理論に基づくと、Ｐ（Ａ（ｘ_２、ｙ_２）｜Ｖ（ｘ_２−ｘ’、ｙ_２−ｙ’））×Ｐ（Ｖ（ｘ_２−ｘ’、ｙ_２−ｙ’））に等しい。

上式に、式（８）、式（９）および式（１０）を代入する。

上式の両辺を、（ｘ_２、ｙ_２）に関して積分する。

上式左辺の積分は１に等しいため、結局、Ｂａｙｅｓの理論に基づくと、以下の関係が成立する。

上式の関係を満たすｈ（ｘ、ｙ）が算出されると、算出されたｈ（ｘ、ｙ）をフーリエ変換することによって、ＯＴＦ分布Ｈ（ｓ、ｔ）を推定することができる。

本発明に係る方法では、式（１４）の右辺のｈ（ｘ、ｙ）をｈ_ｋ（ｘ、ｙ）とし、式（２４）の左辺のｈ（ｘ、ｙ）をｈ_ｋ＋１（ｘ、ｙ）として、ｈ_ｋ（ｘ、ｙ）に関する反復計算を実施し、ｈ_ｋ（ｘ、ｙ）の収束値を求める。上記の反復計算によって求まるｈ_ｋ（ｘ、ｙ）の収束値をフーリエ変換して、Ｂａｙｅｓの理論に基づくＯＴＦの推定分布Ｈ_ｋ（ｘ、ｙ）を得ることができる。

上記した原理を用いるＯＴＦの推定方法は、以下に示す工程を順次実施していくことで実現される。以下では繰り返し計算によって推定されるＰＳＦの分布をｈ_ｋ（ｋ＝０、１、・・・）とする。ＰＳＦの推定分布ｈ_ｋは、反復計算の過程で位相に関する特性を正確に評価するために、複素関数として扱う。
まず劣化画像の分布ｇ（ｘ、ｙ）と、原画像の分布ｆ（ｘ、ｙ）を特定する。本発明の方法では、ｇ（ｘ、ｙ）の実部を劣化画像における照度の分布として、ｇ（ｘ、ｙ）の虚部を０とする。またｆ（ｘ、ｙ）の実部を原画像における照度の分布として、ｆ（ｘ、ｙ）の虚部を０とする。
次に原画像の分布ｆ（ｘ、ｙ）をフーリエ変換して、原画像の分布のスペクトル分布Ｆ（ｓ、ｔ）を算出する。
次に、ＰＳＦの最初の推定分布ｈ_０（ｘ、ｙ）を特定する。ＰＳＦの最初の推定分布ｈ_０（ｘ、ｙ）は、どのような分布としてもよい。例えば、ＰＳＦの最初の推定分布ｈ_０（ｘ、ｙ）は、実部を１として、虚部を０とする。
次に以下に示す演算を、ｈ_ｋ（ｘ、ｙ）が収束するまで繰り返し実施する。ここでＦＴ（）は二次元のフーリエ変換を表し、ＦＴ^−１（）は二次元の逆フーリエ変換を表す。またＦ^＃（ｓ、ｔ）は、Ｆ（ｓ、ｔ）の反転関数であり、Ｆ^＃（ｓ、ｔ）＝Ｆ（−ｓ、−ｔ）である。

ｈ_ｋ（ｘ、ｙ）の収束の判定は、例えば反復計算の回数が予め定められた回数に到達したか否かを基準に判別してもよいし、ｈ_ｋ（ｘ、ｙ）とｈ_ｋ＋１（ｘ、ｙ）の差分の分布を算出して、その差分の絶対値が全ての（ｘ、ｙ）に対してあるしきい値以下となるか否かを基準に判別してもよいし、ｆ_ｋとｆ_ｋ＋１の差分を算出して、算出される差分の絶対値を（ｘ、ｙ）に関して積分した値があるしきい値以下となるか否かで判断してもよい。
ｈ_ｋ（ｘ、ｙ）が収束したら、それをフーリエ変換して、ＯＴＦの推定分布Ｈ_ｋ（ｓ、ｔ）を算出する。

上記の方法を用いることによって、伝達系のＯＴＦを好適に推定することができる。本発明の方法を用いる場合、原画像分布が特定の周波数帯を含まない場合であっても、全ての周波数帯に対するＯＴＦを推定することができる。

上記したように、劣化画像とＯＴＦとが既知であれば原画像を復元することが可能であり、劣化画像と原画像が既知であればＯＴＦを推定することが可能である。
これらの方法を組み合わせることによって、劣化画像のみに基づいて、ＯＴＦを推定し、原画像を復元することが可能となる。

本発明の他の一つの方法は、劣化情報から原情報を復元する方法である。その方法は、劣化情報の分布を特定する工程と、原情報の最初の推定分布を特定する工程と、伝達系のインパルス応答の最初の推定分布を特定する工程と、（Ａ）前記劣化情報の分布と、原情報の推定分布と、インパルス応答の推定分布に基づいて、インパルス応答の推定分布を更新する工程と、（Ｂ）前記劣化情報の分布と、原情報の推定分布と、インパルス応答の推定分布に基づいて、原情報の推定分布を更新する工程と、上記（Ａ）と（Ｂ）の工程を交互に繰返す工程と、原情報の推定分布に基づいて、原情報を出力する工程とを備えている。

上記の方法において前記（Ａ）のインパルス応答の推定分布を更新する工程は、（Ａ１）原情報の推定分布をフーリエ変換して、原情報の推定分布のスペクトル分布を得る工程と、（Ａ２）インパルス応答の推定分布をフーリエ変換して第１の関数を得る工程と、（Ａ３）前記第１の関数に前記原情報の推定分布のスペクトル分布を乗じて第２の関数を得る工程と、（Ａ４）前記第２の関数を逆フーリエ変換して第３の関数を得る工程と、（Ａ５）前記劣化情報の分布を前記第３の関数で除して第４の関数を得る工程と、（Ａ６）前記第４の関数をフーリエ変換して第５の関数を得る工程と、（Ａ７）前記第５の関数に前記原情報の推定分布のスペクトル分布の反転関数を乗じて第６の関数を得る工程と、（Ａ８）前記第６の関数を逆フーリエ変換して第７の関数を得る工程と、（Ａ９）前記インパルス応答の推定分布に前記第７の関数を乗じて、インパルス応答の次の推定分布を得る工程と、（Ａ１０）インパルス応答の次の推定分布をインパルス応答の推定分布に置き換える工程と、を備えている。

上記の方法において前記（Ｂ）の原情報の推定分布を更新する工程は、（Ｂ１）インパルス応答の推定分布をフーリエ変換して、伝達系の周波数応答の推定分布を得る工程と、（Ｂ２）原情報の推定分布をフーリエ変換して第１の関数を得る工程と、（Ｂ３）前記第１の関数に前記周波数応答の推定分布を乗じて第２の関数を得る工程と、（Ｂ４）前記第２の関数を逆フーリエ変換して第３の関数を得る工程と、（Ｂ５）前記劣化情報の分布を前記第３の関数で除して第４の関数を得る工程と、（Ｂ６）前記第４の関数をフーリエ変換して第５の関数を得る工程と、（Ｂ７）前記第５の関数に前記周波数応答の推定分布の反転関数を乗じて第６の関数を得る工程と、（Ｂ８）前記第６の関数を逆フーリエ変換して第７の関数を得る工程と、（Ｂ９）原情報の推定分布に前記第７の関数を乗じて、原情報の次の推定分布を得る工程と、（Ｂ１０）原情報の次の推定分布を原情報の推定分布に置き換える工程と、を備えている。

上記の方法では、原画像とＰＳＦのそれぞれについて、適当な分布を仮定しておき、上記したＰＳＦの推定（すなわち、ＯＴＦの推定）と原画像の推定とを交互に繰返し実施していくことによって、原画像を復元する。
ＯＴＦに関しては、劣化画像の分布と、原画像の推定分布を用いて、Ｂａｙｅｓの理論に基づく計算によって、より改善された推定分布を得ることができる。ＯＴＦの推定に用いる原画像の推定分布が、真の原画像の分布に近いほど、真のＯＴＦの分布に近い推定分布を得ることができる。
原画像に関しては、劣化画像の分布と、ＯＴＦの推定分布を用いて、Ｂａｙｅｓの理論に基づく計算によって、より改善された推定分布を得ることができる。原画像の推定に用いるＯＴＦの推定分布が、真のＯＴＦの分布に近いほど、真の原画像の分布に近い推定分布を得ることができる。
従って、上記したＯＴＦの推定と原画像の推定を交互に繰返し実施していくことによって、原画像の推定分布はより真の原画像の分布に近づいていき、ＯＴＦの推定分布はより真のＯＴＦの分布に近づいていき、結果として良好に復元された原画像と良好に推定されたＯＴＦの双方を得ることができる。

本発明の方法では、以下の工程を順に実施していくことによって、劣化画像のみから原画像を復元することができる。
まず原画像の最初の推定分布と、ＰＳＦの最初の推定分布を設定する。原画像の最初の推定分布は、どのような分布を用いてもよい。一般に、劣化画像の分布は原画像の分布から大きく異なることはないため、原画像の最初の推定分布としては劣化画像の分布を用いることが好ましい。ＰＳＦの最初の推定分布は、どのような分布を用いてもよい。
次に、式（２７）〜（２９）の計算を実施して、より改善されたＰＳＦの推定分布ｈ_ｋ（ｓ、ｔ）を取得して、得られたｈ_ｋ（ｓ、ｔ）をフーリエ変換してＯＴＦの推定分布Ｈ_ｋ（ｓ、ｔ）を取得する。上記の推定を実施することによって、より真のＯＴＦ分布に近いＯＴＦの推定分布を得ることができる。
次に、式（１６）〜（１８）の計算を実施して、より改善された原画像の推定分布ｆ_ｋ（ｘ、ｙ）を取得する。上記の推定を実施することによって、より真の原画像分布に近い原画像の推定分布を得ることができる。
上記したＯＴＦの推定と原画像の推定を、交互に繰返し実施していくことによって、推定される原画像の分布は真の原画像の分布に近づいていき、推定されるＯＴＦの分布は真のＯＴＦの分布に近づいていく。従って、上記の反復計算を実施することによって、原画像を復元することが可能となる。

本発明者が上記の方法を用いて原画像の復元を実施したところ、劣化画像のみから原画像を好適に復元することができた。

上記した劣化情報の復元方法においては、前記（Ａ）のインパルス応答の推定分布を更新する工程は、前記（Ａ２）から（Ａ１０）の工程を繰返し実行し、前記（Ｂ）の原情報の推定分布を更新する工程は、前記（Ｂ２）から（Ｂ１０）の工程を繰返し実行する、ことが好ましい。
上記のインパルス応答の推定分布の更新は、繰返し実施することによって、より尤もらしいインパルス応答の推定分布を得ることができる。また、上記の原情報の推定分布の更新は、繰返し実施することによって、より尤もらしい推定分布を得ることができる。上記のように各推定分布の更新において処理を繰返し実施することによって、より尤もらしい推定分布を用いて後の工程を実施することが可能となるため、原情報を復元するまでに要する反復計算の回数を低減することができると考えられる。

上述した周波数応答の推定方法あるいは劣化情報の復元方法は、各工程をコンピュータに実行させるプログラムとして具現化することができる。コンピュータのハードウェア構成の例を図１２に示す。

上述した周波数応答の分布の推定方法は、下記のような装置を用いることによって、好適に実施することができる。図８に本発明の装置１０００の機能ブロック図を例示する。本発明の装置１０００は、原情報と劣化情報から伝達系の周波数応答の分布を推定する装置である。その装置１０００は、劣化情報の分布を特定する手段１００２と、原情報の分布のスペクトル分布を特定する手段１００４と、伝達系のインパルス応答の最初の推定分布を特定する手段１００６と、（１）インパルス応答の推定分布をフーリエ変換して第１の関数を得る手段１００８と、（２）前記第１の関数に前記原情報の分布のスペクトル分布を乗じて第２の関数を得る手段１０１０と、（３）前記第２の関数を逆フーリエ変換して第３の関数を得る手段１０１２と、（４）前記劣化情報の分布を前記第３の関数で除して第４の関数を得る手段１０１４と、（５）前記第４の関数をフーリエ変換して第５の関数を得る手段１０１６と、（６）前記第５の関数に前記原情報の分布のスペクトル分布の反転関数を乗じて第６の関数を得る手段１０１８と、（７）前記第６の関数を逆フーリエ変換して第７の関数を得る手段１０２０と、（８）前記インパルス応答の推定分布に前記第７の関数を乗じて、インパルス応答の次の推定分布を得る手段１０２２と、インパルス応答の次の推定分布をインパルス応答の推定分布に置き換えて、前記（１）から（８）の手段に繰返し処理を実施させる手段１０２４と、インパルス応答の推定分布をフーリエ変換して、伝達系の周波数応答の分布を得て出力する手段１０２６を備えている。

上述した劣化情報の復元方法は、下記のような装置を用いることによって、好適に実施することができる。図９、図１０および図１１に本発明の装置１１００の機能ブロック図を例示する。本発明の装置１１００は、劣化情報から原情報を復元する装置である。その装置１１００は、劣化情報の分布を特定する手段１１０２と、原情報の最初の推定分布を特定する手段１１０４と、伝達系のインパルス応答の最初の推定分布を特定する手段１１０６と、（Ａ）前記劣化情報の分布と、原情報の推定分布と、インパルス応答の推定分布に基づいて、インパルス応答の推定分布を更新する手段１１０８と、（Ｂ）前記劣化情報の分布と、原情報の推定分布と、インパルス応答の推定分布に基づいて、原情報の推定分布を更新する手段１１１０と、上記（Ａ）と（Ｂ）の手段に交互に繰返し処理を実施させる手段１１１２と、原情報の推定分布に基づいて、原情報を出力する手段１１１４を備えている。
図１０に前記（Ａ）のインパルス応答の推定分布を更新する手段１１０８の機能ブロック図を例示する。前記（Ａ）のインパルス応答の推定分布を更新する手段１１０８は、（Ａ１）原情報の推定分布をフーリエ変換して、原情報の推定分布のスペクトル分布を得る手段１１１６と、（Ａ２）インパルス応答の推定分布をフーリエ変換して第１の関数を得る手段１１１８と、（Ａ３）前記第１の関数に前記原情報の推定分布のスペクトル分布を乗じて第２の関数を得る手段１１２０と、（Ａ４）前記第２の関数を逆フーリエ変換して第３の関数を得る手段１１２２と、（Ａ５）前記劣化情報の分布を前記第３の関数で除して第４の関数を得る手段１１２４と、（Ａ６）前記第４の関数をフーリエ変換して第５の関数を得る手段１１２６と、（Ａ７）前記第５の関数に前記原情報の推定分布のスペクトル分布の反転関数を乗じて第６の関数を得る手段１１２８と、（Ａ８）前記第６の関数を逆フーリエ変換して第７の関数を得る手段１１３０と、（Ａ９）前記インパルス応答の推定分布に前記第７の関数を乗じて、インパルス応答の次の推定分布を得る手段１１３２と、（Ａ１０）インパルス応答の次の推定分布をインパルス応答の推定分布に置き換える手段１１３２を備えている。
図１１に前記（Ｂ）の原情報の推定分布を更新する手段１１１０の機能ブロック図を例示する。前記（Ｂ）の原情報の推定分布を更新する手段１１１０は、（Ｂ１）インパルス応答の推定分布をフーリエ変換して、伝達系の周波数応答の推定分布を得る手段１１３４と、（Ｂ２）原情報の推定分布をフーリエ変換して第１の関数を得る手段１１３６と、（Ｂ３）前記第１の関数に前記周波数応答の推定分布を乗じて第２の関数を得る手段１１３８と、（Ｂ４）前記第２の関数を逆フーリエ変換して第３の関数を得る手段１１４０と、（Ｂ５）前記劣化情報の分布を前記第３の関数で除して第４の関数を得る手段１１４２と、（Ｂ６）前記第４の関数をフーリエ変換して第５の関数を得る手段１１４４と、（Ｂ７）前記第５の関数に前記周波数応答の推定分布の反転関数を乗じて第６の関数を得る手段１１４６と、（Ｂ８）前記第６の関数を逆フーリエ変換して第７の関数を得る手段１１４８と、（Ｂ９）原情報の推定分布に前記第７の関数を乗じて、原情報の次の推定分布を得る手段１１５０と、（Ｂ１０）原情報の次の推定分布を原情報の推定分布に置き換える手段１１５０を備えている。

図１０および図１１に示すように、上記した劣化情報の復元装置１１００において、前記（Ａ）のインパルス応答の推定分布を更新する手段１１０８は、前記（Ａ２）から（Ａ１０）の手段に繰返し処理を実施させる手段１１５２を備え、前記（Ｂ）の原情報の推定分布を更新する手段１１１０は、前記（Ｂ２）から（Ｂ１０）の手段に繰返し処理を実施させる手段１１５４を備えることが好ましい。

本発明の周波数応答を推定する方法、プログラムあるいは装置を用いることで、劣化情報の分布と原情報の分布とに基づいて、伝達系の周波数応答を推定することができる。上記の技術を用いる場合、原情報の分布が特定の周波数帯を備えていない場合であっても、全ての周波数に対する周波数応答を取得することができる。
また本発明の原情報を復元する方法、プログラムあるいは装置を用いることで、劣化情報の分布のみに基づいて、周波数応答の推定と原情報の復元をすることができる。
上記した方法、プログラムあるいは装置は、ペプチド、タンパクを始めとする生体高分子のナノレベルの構造の動態可視化、天文学分野における画像解析への応用が可能である。

本発明の第１実施例に係る方法のフローチャートを示す図である。 本発明の第２実施例に係る方法のフローチャートを示す図である。 本発明の第２実施例に係る方法におけるＰＳＦの推定の工程のフローチャートを示す図である。 本発明の第２実施例に係る方法における原画像の推定の工程のフローチャートを示す図である。 本発明の第１実施例に係る方法の概要を示す図である。 本発明の第２実施例に係る方法の概要を示す図である。 原画像から劣化画像への光学系を介した伝達を模式的に示す図である。 本発明の周波数応答の分布の推定装置１０００の機能ブロック図である。 本発明の劣化情報の復元装置１１００の機能ブロック図である。 復元装置１１００の（Ａ）インパルス応答の推定分布を更新する手段１１０８の機能ブロック図である。 復元装置１１００の（Ｂ）原情報の推定分布を更新する手段１１１０の機能ブロック図である。 コンピュータのハードウェア構成例を示す図である。

符号の説明

１０・・・原画像
１２・・・光学系
１４・・・劣化画像
５２・・・ＰＳＦの実部の推定分布（繰返し数２回）
５４・・・ＰＳＦの実部の推定分布（繰返し数５０回）
５６・・・ＰＳＦの実部の推定分布（繰返し数２５０回）
６２・・・原画像の推定分布
６４・・・ＰＳＦの実部の推定分布
５０２、５０４、５０６、５０８、５１０、５１２、５１４、５１６、５１８、５２０、５２２、５２４、５２６、５２８、５３０・・・第１実施例に係る方法の各ステップ
６０２、６０４、６０６、６０８、６１０、６１２、６１４、６１６・・・第２実施例に係る方法の各ステップ
７０２、７０４、７０６、７０８、７１０、７１２、７１４、７１６、７１８、７２０、
７２２、７２４、７２６・・・ＯＴＦの推定の工程の各ステップ
８０２、８０４、８０６、８０８、８１０、８１２、８１４、８１６、８１８、８２０、８２２、８２４・・・原画像の推定の工程の各ステップ
１０００・・・周波数応答の分布の推定装置
１００２、１００４、１００６、１００８、１０１０、１０１２、１０１４、１０１６、１０１８、１０２０、１０２２、１０２４、１０２６・・・周波数応答の分布の推定装置を構成する各手段
１１００・・・劣化情報の復元装置
１１０２、１１０４、１１０６、１１０８、１１１０、１１１２、１１１４、１１１６、１１１８、１１２０、１１２２、１１２４、１１２６、１１２８、１１３０、１１３２、１１３４、１１３６、１１３８、１１４０、１１４２、１１４４、１１４６、１１４８、１１５０、１１５２、１１５４・・・劣化情報の復元装置を構成する各手段

以下、本発明を具現化した実施例について図面を参照して説明する。

（第１実施例）
図面を参照しながら、本実施例の方法について説明する。図１は本実施例の方法を説明するフローチャートである。
本実施例では、図７に示すように、白黒の原画像１０が、光学系１２を介して伝達された結果、白黒の劣化画像１４となった場合について、原画像１０と劣化画像１４から光学系１２のＯＴＦを推定する方法を扱う。上記の場合において、本実施例の方法は、白黒の原画像１０と、白黒の劣化画像１４とを用いて反復計算を実施し、光学系１２のＯＴＦを推定する。原画像１０と劣化画像１４は、それぞれの画像の大きさが同一であり、画像上の点の位置を（ｘ、ｙ）で表記することができる。
本実施例では、原画像１０を記述する分布ｆ（ｘ、ｙ）と、劣化画像１４を記述する分布ｇ（ｘ、ｙ）と、反復計算の過程で推定するＰＳＦの分布ｈ_ｋ（ｘ、ｙ）（ｋ＝０，１，２，・・・）は、複素関数として扱い、それぞれの位相特性を考慮する。本実施例では、分布ｇ（ｘ、ｙ）の実部を劣化画像における照度の分布とし、分布ｆ（ｘ、ｙ）の実部を原画像における照度の分布とする。

まずステップ５０２では、劣化画像の分布ｇ（ｘ、ｙ）を特定する。本実施例では、ｇ（ｘ、ｙ）の位相特性は不明なため、ｇ（ｘ、ｙ）の実部を劣化画像における照度の分布とし、ｇ（ｘ、ｙ）の虚部は全て０とする。

ステップ５０４では、原画像の分布ｆ（ｘ、ｙ）を特定する。本実施例では、ｆ（ｘ、ｙ）の位相特性は不明なため、ｆ（ｘ、ｙ）の実部を原画像における照度の分布とし、ｆ（ｘ、ｙ）の虚部は全て０とする。

ステップ５０６では、光学系１２のＰＳＦの最初の推定分布ｈ_０（ｘ、ｙ）を設定する。ＰＳＦの最初の推定分布ｈ_０（ｘ、ｙ）としては、どのような分布を設定してもよい。本実施例の方法では、すべての（ｘ、ｙ）に対してｈ_０（ｘ、ｙ）＝１となるように設定する。またステップ５０６では、反復計算の繰返し数ｋを０に設定する。

ステップ５０８では、原画像の分布ｆ（ｘ、ｙ）をフーリエ変換して、その結果をＦ（ｓ、ｔ）に設定する。ｓはｘ方向の空間周波数であり、ｔはｙ方向の空間周波数である。上記のフーリエ変換は、二次元平面内での空間周波数に関するものであり、次式で定義される。

上式および図１のＦＴ（）は、二次元フーリエ変換を示す。
上記したフーリエ変換は、高速フーリエ変換を用いることで好適に実施することができる。
上記によって計算される関数Ｆ（ｓ、ｔ）は、原画像分布ｆ（ｘ、ｙ）の空間スペクトルを表す。

ステップ５１０では、推定されるＰＳＦ分布ｈ_ｋ（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ_１（ｓ、ｔ）に設定する。関数Ｋ_１（ｓ、ｔ）は、第１の関数に相当する。ｋは非負の整数であり、反復計算の回数に応じて後のステップ５２６で増加していく。

ステップ５１２では、ステップ５１０で算出された関数Ｋ_１（ｓ、ｔ）に、ステップ５０８で算出された原画像のスペクトル分布Ｆ（ｓ、ｔ）を乗じて、関数Ｋ_２（ｓ、ｔ）を算出する。関数Ｋ_２（ｓ、ｔ）は、第２の関数に相当する。

ステップ５１４では、ステップ５１２で算出された関数Ｋ_２（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ_３（ｘ、ｙ）に設定する。関数Ｌ_３（ｘ、ｙ）は、第３の関数に相当する。上記の逆フーリエ変換は、二次元平面内での空間周波数に関するスペクトルから、実空間の分布を算出するものであり、次式で定義される。

上式および図１のＦＴ^−１（）は、二次元の逆フーリエ変換を示す。

ステップ５１６では、ステップ５０２で特定された劣化画像分布ｇ（ｘ、ｙ）を、ステップ５１４で算出された関数Ｌ_３（ｘ、ｙ）で除して、関数Ｌ_４（ｘ、ｙ）を算出する。関数Ｌ_４（ｘ、ｙ）は、第４の関数に相当する。

ステップ５１８では、ステップ５１６で算出された関数Ｌ_４（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ_５（ｓ、ｔ）に設定する。関数Ｋ_５（ｓ、ｔ）は、第５の関数に相当する。

ステップ５２０では、ステップ５１８で算出された関数Ｋ_５（ｓ、ｔ）に、Ｆ^＃（ｓ、ｔ）を乗じて、関数Ｋ_６（ｓ、ｔ）を算出する。関数Ｋ_６（ｓ、ｔ）は、第６の関数に相当する。Ｆ^＃（ｓ、ｔ）は、ステップ５０８で算出された原画像のスペクトル分布Ｆ（ｓ、ｔ）の反転関数であって、Ｆ^＃（ｓ、ｔ）＝Ｆ（−ｓ、−ｔ）の関係を満たす。

ステップ５２２では、ステップ５２０で算出した関数Ｋ_６（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ_７（ｘ、ｙ）に設定する。関数Ｌ_７（ｘ、ｙ）は、第７の関数に相当する。

ステップ５２４では、ＰＳＦの推定分布ｈ_ｋ（ｘ、ｙ）に、ステップ５２２で算出された関数Ｌ_７（ｘ、ｙ）を乗じて、改善されたＰＳＦの推定分布ｈ_k+1（ｘ、ｙ）を算出する。前記関数Ｌ７は一般に虚部を備える複素関数となるから、上記の方法によってＰＳＦの推定分布を位相特性も含めて改善することができる。

ステップ５２６では、改善されたＰＳＦの推定分布ｈ_k+1（ｘ、ｙ）とＰＳＦの推定分布ｈ_k（ｘ、ｙ）の差分を算出し、その差分の絶対値がすべての（ｘ、ｙ）に対してあるしきい値εを下回るか否かを判断する。前記差分の絶対値がある（ｘ、ｙ）に対してしきい値ε以上となる場合（ステップ５２６でＮＯの場合）、改善されたＰＳＦの推定分布ｈ_k+1（ｘ、ｙ）は未だ収束に達していないと判断して、処理はステップ５２８へ進む。前記差分の絶対値がすべての（ｘ、ｙ）に対してしきい値εを下回る場合（ステップ５２６でＹＥＳの場合）、改善されたＰＳＦの推定分布ｈ_k+1（ｘ、ｙ）は収束に達したと判断して、処理はステップ５３０へ進む。

ステップ５２８では、反復計算の繰り返し数ｋを１増加させる。処理はステップ５１０へ進み、ステップ５１０からステップ５２４までの処理を再度実施する。

ステップ５３０では、反復計算の結果得られたＰＳＦの推定分布ｈ_k+1（ｘ、ｙ）をフーリエ変換して、光学系１２のＯＴＦの推定分布Ｈ（ｓ、ｔ）を算出する。
上記の方法によって算出されたＨ（ｓ、ｔ）は、ディスプレーに表示してもよいし、印刷装置を用いて紙に印刷してもよいし、ハードディスクなどの記憶装置に記憶してもよいし、通信回線を経由して他のコンピュータに送信してもよい。

図５に、上記の方法の方法によって劣化画像１４と原画像１０から光学系のＯＴＦ分布を推定する様子を示す。図５には、推定されるＯＴＦ分布に対応するＰＳＦ分布の実部を示している。ＰＳＦの推定分布５２は、上記の反復計算を２回繰返した場合の結果を示し、ＰＳＦの推定分布５４は、上記の反復計算を５０回繰返した場合の結果を示し、ＰＳＦの推定分布５６は、上記の反復計算を２５０回繰返した場合を示している。上記の反復計算の繰返し数を増やすことによって、推定されるＰＳＦ分布は収束していく。すなわち、推定されるＯＴＦ分布は収束していく。

（第２実施例）
本発明に係る他の一つの実施例を、図面を参照しながら説明する。
本実施例では、図７に示すように、白黒の原画像１０が、光学系１２を介して伝達された結果、白黒の劣化画像１４となった場合について、劣化画像１４のみに基づいて、原画像１０の復元と、光学系１２のＯＴＦの推定を行う。上記の場合において、本実施例の方法は、原画像の分布とＯＴＦの分布のそれぞれについて、適当な分布を仮定しておく。次に仮定した原画像の分布と劣化画像の分布とを用いた反復計算によって、改善された光学系１２のＯＴＦ分布を推定する。次に、推定されたＯＴＦ分布と劣化画像の分布とを用いた反復計算によって、改善された原画像の分布を推定する。その後、改善された原画像の分布と劣化画像の分布を用いた反復計算を実施して、さらに改善されたＯＴＦ分布を推定する。上記のように、原画像分布の推定とＯＴＦ分布の推定を交互に繰返し実施していき、最終的にＢａｙｅｓ理論に基づいた原画像の推定分布とＯＴＦの推定分布を取得する。

原画像１０と劣化画像１４は、それぞれの画像の大きさが同一であり、画像上の点の位置を（ｘ、ｙ）で表記することができる。
本実施例では、劣化画像１４を記述する分布ｇ（ｘ、ｙ）と、反復計算の過程で推定される原画像１０を記述する分布ｆ_ｋ（ｘ、ｙ）およびｆ_ｋ ^ｎ（ｘ、ｙ）と、反復計算の過程で推定されるＰＳＦの分布ｈ_ｋ（ｘ、ｙ）およびｈ_ｋ ^ｍ（ｓ、ｙ）は、全て複素関数として扱い、それぞれの位相特定を考慮する。本実施例では、分布ｇ（ｘ、ｙ）の実部を劣化画像における照度の分布とし、分布ｆ_ｋ（ｘ、ｙ）およびｆ_ｋ ^ｎ（ｘ、ｙ）の実部を原画像における照度の推定分布とする。

図２に本実施例の画像復元方法のフローチャートを示す。
ステップ６０２では、劣化画像の分布ｇ（ｘ、ｙ）を特定する。本実施例の方法では、ｇ（ｘ、ｙ）の実部を劣化画像における照度の分布とし、ｇ（ｘ、ｙ）の虚部を全て０とする。

ステップ６０４では、原画像の最初の推定分布ｆ_０（ｘ、ｙ）を設定する。原画像の最初の推定分布ｆ_０（ｘ、ｙ）としては、任意の分布を設定することができる。本実施例の方法では、原画像の最初の推定分布ｆ_０（ｘ、ｙ）として、劣化画像の分布ｇ（ｘ、ｙ）を用いる。一般的に、原画像の分布と劣化画像の分布は大きく異ならないと考えられるため、上記のように原画像の最初の推定分布ｆ_０（ｘ、ｙ）を設定することによって、原画像の復元に伴う反復計算の回数を低減することができる。
またステップ６０４では、ＰＳＦの最初の推定分布ｈ_０（ｘ、ｙ）を設定する。ＰＳＦの最初の推定分布ｈ_０（ｘ、ｙ）としては、任意の分布を設定することができる。本実施例の方法では、ＰＳＦの最初の推定分布ｈ_０（ｘ、ｙ）として、全ての（ｘ、ｙ）に対して実部が１であり、虚部が０である関数を用いる。
さらにステップ６０４では、反復計算の繰返し数ｋを０に設定する。

ステップ６０６では、反復計算の繰返し数ｋを１増加させる。

ステップ６０８では、ＰＳＦの推定分布ｈ_ｋ（ｘ、ｙ）を算出する。この計算については後述する。

ステップ６１０では、ＰＳＦの推定分布ｈ_ｋ（ｘ、ｙ）をフーリエ変換して、ＯＴＦの推定分布Ｈ_ｋ（ｓ、ｔ）を算出する。

ステップ６１２では、原画像の推定分布ｆ_ｋ（ｘ、ｙ）を算出する。この計算については後述する。

ステップ６１４では、ステップ６１２で更新された原画像の推定分布ｆ_ｋ（ｘ、ｙ）と、更新される前の原画像の推定分布ｆ_ｋ−１との差分を算出し、その差分がすべての（ｘ、ｙ）に対してあるしきい値εを下回るか否かを判断する。前記差分がある（ｘ、ｙ）に対してしきい値ε以上の場合（ステップ６１４でＮＯの場合）、原画像の推定分布ｆ_ｋ（ｘ、ｙ）は未だ収束に達していないと判断して、処理はステップ６０６へ進み、ステップ６０６からステップ６１２までの処理を再度繰返す。前記差分がすべての（ｘ、ｙ）に対してしきい値εを下回る場合（ステップ６１４でＹＥＳの場合）、原画像の推定分布ｆ_ｋ（ｘ、ｙ）は収束に達したと判断して、処理はステップ６１６へ進む。

ステップ６１６では、反復計算の結果得られた原画像の推定分布ｆ_ｋ（ｘ、ｙ）を、復元された原画像の分布ｆ（ｘ、ｙ）に設定する。本実施例の方法では、反復計算の過程で画像を記述する分布を複素関数として扱っているため、上記で設定される原画像の分布ｆ（ｘ、ｙ）は一般に複素関数となる。本実施例の方法では、復元された原画像の分布ｆ（ｘ、ｙ）の実部が、原画像における照度の分布を表す。
上記の方法によって算出されたｆ（ｘ、ｙ）は、ディスプレーに表示してもよいし、印刷装置を用いて紙に印刷してもよいし、ハードディスクなどの記憶装置に記憶してもよいし、通信回線を経由して他のコンピュータに送信してもよい。

以下ではステップ６０８のｈ_ｋ（ｘ、ｙ）の推定について詳細に説明する。本実施例の方法では、ｈ_ｋ（ｘ、ｙ）の推定として、図３のフローチャートに示すステップ７０２〜７２６の工程を実施する。

ステップ７０２では、ｈ_ｋ（ｘ、ｙ）の最初の推定分布ｈ_ｋ ^０（ｘ、ｙ）を設定する。本実施例の方法では、すでに取得されているｈ_ｋ−１（ｘ、ｙ）を、ｈ_ｋ（ｘ、ｙ）の最初の推定ｈ_ｋ ^０（ｘ、ｙ）とする。また、ステップ７０２では、ｈ_ｋ（ｘ、ｙ）の推定に係る反復計算の繰返し数ｍを０に設定する。

ステップ７０４では、原画像のスペクトルの推定分布Ｆ（ｓ、ｔ）を設定する。本実施例の方法では、すでに取得されている原画像の推定分布ｆ_ｋ−１（ｘ、ｙ）をフーリエ変換して、その結果をＦ（ｓ、ｔ）とする。

ステップ７０６では、ｈ_ｋ（ｘ、ｙ）の推定分布ｈ_ｋ ^ｍ（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ_１（ｓ、ｔ）に設定する。関数Ｋ_１（ｓ、ｔ）は、第１の関数に相当する。

ステップ７０８では、ステップ７０６で算出された関数Ｋ_１（ｓ、ｔ）に、ステップ７０４で算出された原画像のスペクトルの推定分布Ｆ（ｓ、ｔ）を乗じて、関数Ｋ_２（ｓ、ｔ）を算出する。関数Ｋ_２（ｓ、ｔ）は、第２の関数に相当する。

ステップ７１０では、ステップ７０８で算出された関数Ｋ_２（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ_３（ｘ、ｙ）に設定する。関数Ｌ_３（ｘ、ｙ）は、第３の関数に相当する。

ステップ７１２では、図２のステップ６０４で特定された劣化画像分布ｇ（ｘ、ｙ）を、図３のステップ７１０で算出された関数Ｌ_３（ｘ、ｙ）で除して、関数Ｌ_４（ｘ、ｙ）を算出する。関数Ｌ_４（ｘ、ｙ）は、第４の関数に相当する。

ステップ７１４では、ステップ７１２で算出された関数Ｌ_４（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ_５（ｓ、ｔ）に設定する。関数Ｋ_５（ｓ、ｔ）は、第５の関数に相当する。

ステップ７１６では、ステップ７１４で算出された関数Ｋ_５（ｓ、ｔ）に、Ｆ^＃（ｓ、ｔ）を乗じて、関数Ｋ_６（ｓ、ｔ）を算出する。関数Ｋ_６（ｓ、ｔ）は、第６の関数に相当する。Ｆ^＃（ｓ、ｔ）は、ステップ７０４で算出された原画像のスペクトルの推定分布Ｆ（ｓ、ｔ）の反転関数であり、Ｆ^＃（ｓ、ｔ）＝Ｆ（−ｓ、−ｔ）の関係を満たす。

ステップ７１８では、ステップ７１６で算出された関数Ｋ_６（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ_７（ｘ、ｙ）に設定する。関数Ｌ_７（ｘ、ｙ）は、第７の関数に相当する。

ステップ７２０では、ｈ_ｋ（ｘ、ｙ）の推定分布ｈ_ｋ ^ｍ（ｘ、ｙ）に、ステップ７１８で設定された関数Ｌ_７（ｘ、ｙ）を乗じて、その結果をｈ_ｋ（ｘ、ｙ）の改善された推定分布ｈ_ｋ ^ｍ＋１（ｘ、ｙ）に設定する。

ステップ７２２では、ｈ_ｋ（ｘ、ｙ）の推定のための反復計算の繰返し数ｍを１増加させる。

ステップ７２４では、ｈ_ｋ（ｘ、ｙ）の推定のための反復計算の繰返し数ｍが、しきい値以上になったか否かを判断する。本実施例では、ｈ_ｋ（ｘ、ｙ）の推定のために実施する反復計算の回数を５回としている。繰り返し数ｍが５未満の場合（ステップ７２４でＮＯの場合）、処理はステップ７０６へ進み、ステップ７０６からステップ７２２までの処理を再度実施する。繰り返し数ｍが５以上の場合（ステップ７２４でＹＥＳの場合）、処理はステップ７２６へ進む。

ステップ７２６では、上記の反復計算の結果得られる分布ｈ_ｋ ^５（ｘ、ｙ）を、ｈ_ｋ（ｘ、ｙ）の推定分布として設定する。
上記の反復計算は、劣化画像分布ｇ（ｘ、ｙ）と、原画像分布ｆ_ｋ−１（ｘ、ｙ）から、ＰＳＦの分布ｈ_ｋ（ｘ、ｙ）を推定することに相当する。原画像分布ｆ_ｋ−１（ｘ、ｙ）が、真の原画像分布ｆ（ｘ、ｙ）に近いほど、上記の反復計算によって推定されるｈ_ｋ（ｘ、ｙ）は、真のＰＳＦの分布ｈ（ｘ、ｙ）に近づく。

以下では図２のステップ６１２のｆ_ｋ（ｘ、ｙ）の推定について詳細に説明する。本実施例の方法では、ｆ_ｋ（ｘ、ｙ）の推定として、図４のフローチャートに示すステップ８０２〜８２４の工程を実施する。

ステップ８０２では、ｆ_ｋ（ｘ、ｙ）の最初の推定分布ｆ_ｋ ^０（ｘ、ｙ）を設定する。本実施例の方法では、既に取得されているｆ_ｋ−１（ｘ、ｙ）を、ｆ_ｋ（ｘ、ｙ）の最初の推定分布ｆ_ｋ ^０（ｘ、ｙ）として設定する。また、ステップ８０２では、ｆ_ｋ（ｘ、ｙ）の推定のための反復計算の繰返し数ｎを０に設定する。

ステップ８０４では、ｆ_ｋ（ｘ、ｙ）の推定分布ｆ_ｋ ^ｎ（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ_１（ｓ、ｔ）に設定する。関数Ｋ_１（ｓ、ｔ）は、第１の関数に相当する。

ステップ８０６では、ステップ８０４で算出された関数Ｋ_１（ｓ、ｔ）に、図２のステップ６１０で算出されたＯＴＦの推定分布Ｈ_ｋ（ｓ、ｔ）を乗じて、関数Ｋ_２（ｓ、ｔ）を算出する。関数Ｋ_２（ｓ、ｔ）は、第２の関数に相当する。

ステップ８０８では、ステップ８０６で算出された関数Ｋ_２（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ_３（ｘ、ｙ）に設定する。関数Ｌ_３（ｘ、ｙ）は、第３の関数に相当する。

ステップ８１０では、図２のステップ６０２で特定された劣化画像分布ｇ（ｘ、ｙ）を、図４のステップ８０８で算出された関数Ｌ_３（ｘ、ｙ）で除して、関数Ｌ_４（ｘ、ｙ）を算出する。関数Ｌ_４（ｘ、ｙ）は、第４の関数に相当する。

ステップ８１２では、ステップ８１０で算出された関数Ｌ_４（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ_５（ｓ、ｔ）に設定する。関数Ｋ_５（ｓ、ｔ）は、第５の関数に相当する。

ステップ８１４では、ステップ８１２で算出された関数Ｋ_５（ｓ、ｔ）に、Ｈ_ｋ ^＃（ｓ、ｔ）を乗じて、関数Ｋ_６（ｓ、ｔ）を算出する。Ｈ_ｋ ^＃（ｓ、ｔ）は、図２のステップ６１０で算出されたＯＴＦの推定分布Ｈ_ｋ（ｓ、ｔ）の反転関数であり、Ｈ_ｋ ^＃（ｓ、ｔ）＝Ｈ_ｋ（−ｓ、−ｔ）の関係を満たす。

ステップ８１６では、ステップ８１４で算出された関数Ｋ_６（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ_７（ｘ、ｙ）に設定する。関数Ｌ_７（ｘ、ｙ）は、第７の関数に相当する。

ステップ８１８では、ｆ_ｋ（ｘ、ｙ）の推定分布ｆ_ｋ ^ｎ（ｘ、ｙ）に、ステップ８１６で算出された関数Ｌ_７（ｘ、ｙ）を乗じて、ｆ_ｋ（ｘ、ｙ）の改善された推定分布ｆ_ｋ ^ｎ＋１（ｘ、ｙ）を算出する。

ステップ８２０では、ｆ_ｋ（ｘ、ｙ）の推定のための反復計算の繰返し数ｎを１増加する。

ステップ８２２では、ｆ_ｋ（ｘ、ｙ）の推定のための反復計算の繰返し数ｎが、しきい値以上になったか否かを判断する。本実施例では、ｆ_ｋ（ｘ、ｙ）の推定のために実施する反復計算の回数を５回としている。繰り返し数ｎが５未満の場合（ステップ８２２でＮＯの場合）、処理はステップ８０４へ進み、ステップ８０４からステップ８２０までの処理を再度実施する。繰り返し数ｎが５以上の場合（ステップ８２２でＹＥＳの場合）、処理はステップ８２４へ進む。

ステップ８２４では、上記の反復計算の結果得られる分布ｆ_ｋ ^５（ｘ、ｙ）を、ｆ_ｋ（ｘ、ｙ）の推定分布として設定する。
上記の反復計算は、劣化画像分布ｇ（ｘ、ｙ）と、ＯＴＦ分布Ｈ_ｋ（ｓ、ｔ）とに基づいて、原画像分布ｆ_ｋ（ｘ、ｙ）を推定することに相当する。ＯＴＦ分布Ｈ_ｋ（ｓ、ｔ）が、真のＯＴＦ分布Ｈ（ｓ、ｔ）に近いほど、上記の反復計算によって推定されるｆ_ｋ（ｘ、ｙ）は、真の原画像の分布ｆ（ｘ、ｙ）に近づく。

図６に、上記の方法を用いて、劣化画像１４のみに基づいて反復計算を実施して、原画像６２を復元した結果を示す。劣化画像１４では輪郭がぼやけてほとんど認識できない文字が、復元された原画像６２でははっきりと文字として認識することができる。また、図６には原画像６２の復元と並行して推定される光学系１２のＯＴＦに対応するＰＳＦの実部の分布６４を示す。

上記の実施例では、白黒の原画像を復元する場合を説明したが、原画像が色を備えている場合にも、同様の手法を用いて復元することができる。原画像が色を備えている場合には、原画像におけるＲＧＢそれぞれの色成分の照度の分布ｆｒ（ｘ、ｙ）、ｆｇ（ｘ、ｙ）、ｆｂ（ｘ、ｙ）について、劣化画像におけるＲＧＢそれぞれの色成分の照度の分布ｇｒ（ｘ、ｙ）、ｇｇ（ｘ、ｙ）、ｇｂ（ｘ、ｙ）に基づいて、上記の実施例の方法を用いて個別に推定することができる。例えば、原画像におけるＲの色成分の照度の分布ｆｒ（ｘ、ｙ）は、劣化画像におけるＲの色成分の照度の分布ｇｒ（ｘ、ｙ）を用いた反復計算によって、推定することができる。同様にして、原画像におけるＧの色成分の照度の分布ｆｇ（ｘ、ｙ）、Ｂの色成分の照度の分布ｆｂ（ｘ、ｙ）についても、推定することが可能である。上記によって推定される原画像のＲＧＢそれぞれの照度分布から、原画像を復元することができる。

上記の実施例では、反復計算の過程において、先にＯＴＦの推定分布を更新して、次に原画像の推定分布を更新する例を説明した。ＯＴＦの推定分布の更新と、原画像の推定分布の更新は、交互に繰返し実施すればよいため、先に原画像の推定分布を更新して、次にＯＴＦの推定分布を更新してもよい。

以上、本発明の実施形態について詳細に説明したが、これらは例示に過ぎず、特許請求の範囲を限定するものではない。特許請求の範囲に記載の技術には、以上に例示した具体例を様々に変形、変更したものが含まれる。
また、本明細書または図面に説明した技術要素は、単独であるいは各種の組み合わせによって技術的有用性を発揮するものであり、出願時請求項記載の組み合わせに限定されるものではない。また、本明細書または図面に例示した技術は複数目的を同時に達成するものであり、そのうちの一つの目的を達成すること自体で技術的有用性を持つものである。

Claims (9)

1. 原画像と劣化画像から光学系の位相特性を反映した光学的伝達関数（ＯＴＦ；Optical Transfer Function）の分布を推定する方法であって、
   劣化画像の分布を特定する工程と、
   原画像の分布のスペクトル分布を特定する工程と、
   光学系の点像強度分布関数（ＰＳＦ；Point Spread Function）の最初の推定分布を特定する工程と、
   （１）ＰＳＦの推定分布をフーリエ変換して第１の関数を得る工程と、
   （２）前記第１の関数に前記原画像の分布のスペクトル分布を乗じて第２の関数を得る工程と、
   （３）前記第２の関数を逆フーリエ変換して第３の関数を得る工程と、
   （４）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る工程と、
   （５）前記第４の関数をフーリエ変換して第５の関数を得る工程と、
   （６）前記第５の関数に前記原画像の分布のスペクトル分布の反転関数を乗じて第６の関数を得る工程と、
   （７）前記第６の関数を逆フーリエ変換して第７の関数を得る工程と、
   （８）前記ＰＳＦの推定分布に前記第７の関数を乗じて、ＰＳＦの次の推定分布を得る工程と、
   ＰＳＦの次の推定分布をＰＳＦの推定分布に置き換えて、前記（１）から（８）の工程を繰返す工程と、
   ＰＳＦの推定分布をフーリエ変換して、ＯＴＦの分布を得て出力する工程とを備えることを特徴とするＯＴＦの分布の推定方法。
2. 劣化画像から原画像を復元する方法であって、
   劣化画像の分布を特定する工程と、
   原画像の最初の推定分布を特定する工程と、
   光学系の点像強度分布関数（ＰＳＦ；Point Spread Function）の最初の推定分布を特定する工程と、
   （Ａ）前記劣化画像の分布と、原画像の推定分布と、ＰＳＦの推定分布に基づいて、ＰＳＦの推定分布を更新する工程と、
   （Ｂ）前記劣化画像の分布と、原画像の推定分布と、ＰＳＦの推定分布に基づいて、原画像の推定分布を更新する工程と、
   上記（Ａ）と（Ｂ）の工程を交互に繰返す工程と、
   原画像の推定分布に基づいて、原画像を出力する工程と
   を備え、
   前記（Ａ）のＰＳＦの推定分布を更新する工程は、
   （Ａ１）原画像の推定分布をフーリエ変換して、原画像の推定分布のスペクトル分布を得る工程と、
   （Ａ２）ＰＳＦの推定分布をフーリエ変換して第１の関数を得る工程と、
   （Ａ３）前記第１の関数に前記原画像の推定分布のスペクトル分布を乗じて第２の関数を得る工程と、
   （Ａ４）前記第２の関数を逆フーリエ変換して第３の関数を得る工程と、
   （Ａ５）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る工程と、
   （Ａ６）前記第４の関数をフーリエ変換して第５の関数を得る工程と、
   （Ａ７）前記第５の関数に前記原画像の推定分布のスペクトル分布の反転関数を乗じて第６の関数を得る工程と、
   （Ａ８）前記第６の関数を逆フーリエ変換して第７の関数を得る工程と、
   （Ａ９）前記ＰＳＦの推定分布に前記第７の関数を乗じて、ＰＳＦの次の推定分布を得る工程と、
   （Ａ１０）ＰＳＦの次の推定分布をＰＳＦの推定分布に置き換える工程と、
   を備え、
   前記（Ｂ）の原画像の推定分布を更新する工程は、
   （Ｂ１）ＰＳＦの推定分布をフーリエ変換して、光学系の位相特性を反映した光学的伝達関数（ＯＴＦ；Optical Transfer Function）の推定分布を得る工程と、
   （Ｂ２）原画像の推定分布をフーリエ変換して第１の関数を得る工程と、
   （Ｂ３）前記第１の関数に前記ＯＴＦの推定分布を乗じて第２の関数を得る工程と、
   （Ｂ４）前記第２の関数を逆フーリエ変換して第３の関数を得る工程と、
   （Ｂ５）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る工程と、
   （Ｂ６）前記第４の関数をフーリエ変換して第５の関数を得る工程と、
   （Ｂ７）前記第５の関数に前記ＯＴＦの推定分布の反転関数を乗じて第６の関数を得る工程と、
   （Ｂ８）前記第６の関数を逆フーリエ変換して第７の関数を得る工程と、
   （Ｂ９）原画像の推定分布に前記第７の関数を乗じて、原画像の次の推定分布を得る工程と、
   （Ｂ１０）原画像の次の推定分布を原画像の推定分布に置き換える工程とを備えることを特徴とする、劣化画像の復元方法。
3. 前記（Ａ）のＰＳＦの推定分布を更新する工程は、
   前記（Ａ２）から（Ａ１０）の工程を繰返し実行し、
   前記（Ｂ）の原画像の推定分布を更新する工程は、
   前記（Ｂ２）から（Ｂ１０）の工程を繰返し実行する、
   ことを特徴とする請求項２の劣化画像の復元方法。
4. 請求項１の方法における各工程をコンピュータに実行させるためのプログラム。
5. 請求項２の方法における各工程をコンピュータに実行させるためのプログラム。
6. 請求項３の方法における各工程をコンピュータに実行させるためのプログラム。
7. 原画像と劣化画像から光学系の位相特性を反映した光学的伝達関数（ＯＴＦ；Optical Transfer Function）の分布を推定する装置であって、
   劣化画像の分布を特定する手段と、
   原画像の分布のスペクトル分布を特定する手段と、
   光学系の点像強度分布関数（ＰＳＦ；Point Spread Function）の最初の推定分布を特定する手段と、
   （１）ＰＳＦの推定分布をフーリエ変換して第１の関数を得る手段と、
   （２）前記第１の関数に前記原画像の分布のスペクトル分布を乗じて第２の関数を得る手段と、
   （３）前記第２の関数を逆フーリエ変換して第３の関数を得る手段と、
   （４）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る手段と、
   （５）前記第４の関数をフーリエ変換して第５の関数を得る手段と、
   （６）前記第５の関数に前記原画像の分布のスペクトル分布の反転関数を乗じて第６の関数を得る手段と、
   （７）前記第６の関数を逆フーリエ変換して第７の関数を得る手段と、
   （８）前記ＰＳＦの推定分布に前記第７の関数を乗じて、ＰＳＦの次の推定分布を得る手段と、
   ＰＳＦの次の推定分布をＰＳＦの推定分布に置き換えて、前記（１）から（８）の手段に繰返し処理を実施させる手段と、
   ＰＳＦの推定分布をフーリエ変換して、ＯＴＦの分布を得て出力する手段とを備えることを特徴とするＯＴＦの分布の推定装置。
8. 劣化画像から原画像を復元する装置であって、
   劣化画像の分布を特定する手段と、
   原画像の最初の推定分布を特定する手段と、
   光学系の点像強度分布関数（ＰＳＦ；Point Spread Function）の最初の推定分布を特定する手段と、
   （Ａ）前記劣化画像の分布と、原画像の推定分布と、ＰＳＦの推定分布に基づいて、ＰＳＦの推定分布を更新する手段と、
   （Ｂ）前記劣化画像の分布と、原画像の推定分布と、ＰＳＦの推定分布に基づいて、原画像の推定分布を更新する手段と、
   上記（Ａ）と（Ｂ）の手段に交互に繰返し処理を実施させる手段と、
   原画像の推定分布に基づいて、原画像を出力する手段と
   を備え、
   前記（Ａ）のＰＳＦの推定分布を更新する手段は、
   （Ａ１）原画像の推定分布をフーリエ変換して、原画像の推定分布のスペクトル分布を得る手段と、
   （Ａ２）ＰＳＦの推定分布をフーリエ変換して第１の関数を得る手段と、
   （Ａ３）前記第１の関数に前記原画像の推定分布のスペクトル分布を乗じて第２の関数を得る手段と、
   （Ａ４）前記第２の関数を逆フーリエ変換して第３の関数を得る手段と、
   （Ａ５）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る手段と、
   （Ａ６）前記第４の関数をフーリエ変換して第５の関数を得る手段と、
   （Ａ７）前記第５の関数に前記原画像の推定分布のスペクトル分布の反転関数を乗じて第６の関数を得る手段と、
   （Ａ８）前記第６の関数を逆フーリエ変換して第７の関数を得る手段と、
   （Ａ９）前記ＰＳＦの推定分布に前記第７の関数を乗じて、ＰＳＦの次の推定分布を得る手段と、
   （Ａ１０）ＰＳＦの次の推定分布をＰＳＦの推定分布に置き換える手段と、
   を備え、
   前記（Ｂ）の原画像の推定分布を更新する手段は、
   （Ｂ１）ＰＳＦの推定分布をフーリエ変換して、光学系の位相特性を反映した光学的伝達関数（ＯＴＦ；Optical Transfer Function）の推定分布を得る手段と、
   （Ｂ２）原画像の推定分布をフーリエ変換して第１の関数を得る手段と、
   （Ｂ３）前記第１の関数に前記ＯＴＦの推定分布を乗じて第２の関数を得る手段と、
   （Ｂ４）前記第２の関数を逆フーリエ変換して第３の関数を得る手段と、
   （Ｂ５）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る手段と、
   （Ｂ６）前記第４の関数をフーリエ変換して第５の関数を得る手段と、
   （Ｂ７）前記第５の関数に前記ＯＴＦの推定分布の反転関数を乗じて第６の関数を得る手段と、
   （Ｂ８）前記第６の関数を逆フーリエ変換して第７の関数を得る手段と、
   （Ｂ９）原画像の推定分布に前記第７の関数を乗じて、原画像の次の推定分布を得る手段と、
   （Ｂ１０）原画像の次の推定分布を原画像の推定分布に置き換える手段とを備えることを特徴とする、劣化画像の復元装置。
9. 前記（Ａ）のＰＳＦの推定分布を更新する手段は、
   前記（Ａ２）から（Ａ１０）の手段に繰返し処理を実施させる手段を備え、
   前記（Ｂ）の原画像の推定分布を更新する手段は、
   前記（Ｂ２）から（Ｂ１０）の手段に繰返し処理を実施させる手段を備える、
   ことを特徴とする請求項８の劣化画像の復元装置。

Images