劣化した画像を復元する装置、方法およびプログラム
【課題】 撮像時のぶれやぼけの影響を表現する点像分布関数を適切に取得して、劣化した画像から劣化のない画像を復元することが可能な技術を提供する。
【解決手段】 本発明の撮像装置は、被写体を撮像して被写体画像データを出力する第１撮像手段と、点光源を撮像して点像画像データを出力する第２撮像手段と、被写体画像データと点像画像データに基いて改善された被写体画像データを生成する画像処理手段を備えている。第１撮像手段と第２撮像手段は、共通する支持台によって支持されており、互いに同期して撮像するものである。被写体と点光源は、共通する載置台上に載置されている。画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と被写体画像データに基いて改善された被写体画像データを生成する。
【解決手段】 本発明の撮像装置は、被写体を撮像して被写体画像データを出力する第１撮像手段と、点光源を撮像して点像画像データを出力する第２撮像手段と、被写体画像データと点像画像データに基いて改善された被写体画像データを生成する画像処理手段を備えている。第１撮像手段と第２撮像手段は、共通する支持台によって支持されており、互いに同期して撮像するものである。被写体と点光源は、共通する載置台上に載置されている。画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と被写体画像データに基いて改善された被写体画像データを生成する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は撮像時のぶれやぼけの影響によって劣化した画像から、劣化のない画像を復元する技術に関する。
【背景技術】
【０００２】
被写体を撮像する際に、解像度が不足したり、画像がぶれてしまったり、ぼけてしまったりすることがある。このように劣化した画像から劣化のない理想的な画像を復元する技術の開発が進められている。
【０００３】
本発明者らは、これまでにも、劣化した画像から劣化のない画像を復元する技術を開発している。この技術は、例えば特許文献１や特許文献２に記載されている。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】国際公開第２００６／０４１１２６号
【特許文献２】国際公開第２００６／０４１１２７号
【発明の概要】
【発明が解決しようとする課題】
【０００５】
劣化した画像から劣化のない画像を復元するためには、撮像時のぶれやぼけの影響を表現する点像分布関数が必要とされる。適切な点像分布関数を取得することができれば、劣化した画像から劣化のない画像を正確かつ速やかに復元することが可能となる。
【０００６】
本発明は上記課題を解決する。本発明では、撮像時のぶれやぼけの影響を表現する点像分布関数を適切に取得して、劣化した画像から劣化のない画像を復元することが可能な技術を提供する。
【課題を解決するための手段】
【０００７】
本発明は撮像装置として具現化される。その撮像装置は、被写体を撮像して被写体画像データを出力する第１撮像手段と、点光源を撮像して点像画像データを出力する第２撮像手段と、被写体画像データと点像画像データに基いて改善された被写体画像データを生成する画像処理手段を備えている。その撮像装置において、第１撮像手段と第２撮像手段は、共通する支持台によって支持されており、互いに同期して撮像するものである。その撮像装置において、被写体と点光源は、共通する載置台上に載置されている。その撮像装置において、画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と被写体画像データに基いて改善された被写体画像データを生成する。
【０００８】
上記の撮像装置の動作内容を方法の発明として具現化することもできる。本発明の撮像方法は、第１撮像手段によって被写体を撮像して被写体画像データを出力する工程と、第２撮像手段によって点光源を撮像して点像画像データを出力する工程と、画像処理手段によって被写体画像データと点像画像データに基いて改善された被写体画像データを生成する工程を備えている。その撮像方法において、第１撮像手段と第２撮像手段は、共通する支持台によって支持されており、互いに同期して撮像するものである。その撮像方法において、被写体と点光源は、共通する載置台上に載置されている。その撮像方法において、画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と被写体画像データに基いて改善された被写体画像データを生成する。
なお本発明は上記の撮像方法における各工程をコンピュータに実行させるためのプログラムとしても具現化することができる。
【０００９】
第１撮像手段で撮像される被写体の画像は、ぶれやぼけによって劣化していることがある。上記の撮像装置では、第１撮像手段と同様の条件で撮像する第２撮像手段によって点光源を撮像して、第２撮像手段で撮像された画像から点像分布関数（ＰＳＦ：Point Spread Function）を取得する。そして、第１撮像手段で撮像された画像とＰＳＦを用いて、ぶれやぼけのない改善された被写体の画像を生成する。
【００１０】
以下では撮像時のぶれやぼけによって劣化している画像を劣化画像と呼び、ぶれやぼけのない理想的な画像を原画像と呼ぶ。劣化画像とＰＳＦから原画像を復元する原理について説明する。原画像と劣化画像はそれぞれの大きさが同一であり、いずれも画像内の点の位置を座標（ｘ、ｙ）で表現することが可能であって、原画像の照度分布はｆｒ（ｘ、ｙ）、劣化画像の照度分布はｇｒ（ｘ、ｙ）で表現されるものとする。またＰＳＦの分布をｈｒ（ｘ、ｙ）とする。
【００１１】
原画像の分布ｆｒ（ｘ、ｙ）と、劣化画像の分布ｇｒ（ｘ、ｙ）と、ＰＳＦの分布ｈｒ（ｘ、ｙ）の間には、次の関係が成り立つ。
【００１２】
【数１】
【００１３】
原画像の分布ｆｒ（ｘ、ｙ）と劣化画像の分布ｇｒ（ｘ、ｙ）は、以下に示す二次元フーリエ変換を実施することによって、ｘ軸に関する空間周波数ｓと、ｙ軸に関する空間周波数ｔについてのスペクトルをそれぞれ得ることができる。
【００１４】
【数２】
【数３】
【００１５】
ここでＦＴ（）は二次元フーリエ変換を示している。またＰＳＦの分布ｈｒ（ｘ、ｙ）について二次元フーリエ変換を実施することによって、周波数空間における伝達関数である光学的伝達関数（ＯＴＦ：Optical Transfer Function）の分布Ｈｒ（ｓ、ｔ）を得ることができる。
【００１６】
【数４】
【００１７】
式（１）の両辺についてフーリエ変換することで、以下に示す劣化画像のスペクトル分布Ｇｒ（ｓ、ｔ）と、原画像のスペクトル分布Ｆｒ（ｓ、ｔ）と、ＯＴＦ分布Ｈｒ（ｓ、ｔ）の関係が得られる。
【００１８】
【数５】
【００１９】
従って、劣化画像のスペクトル分布Ｇｒ（ｓ、ｔ）とＯＴＦ分布Ｈｒ（ｓ、ｔ）が既知であれば、次式により原画像のスペクトル分布Ｆｒ（ｓ、ｔ）を得ることができる。
【００２０】
【数６】
【００２１】
上式によって原画像のスペクトル分布Ｆｒ（ｓ、ｔ）が得られれば、二次元の逆フーリエ変換を施すことによって、原画像の分布ｆｒ（ｘ、ｙ）を得ることができる。
【００２２】
【数７】
【００２３】
ここでＦＴ−１（）は二次元の逆フーリエ変換を示している。以上のように、本発明の画像取得装置によれば、第２撮像手段で撮像された画像からＰＳＦの分布を取得し、第１撮像手段で撮像された画像とＰＳＦの分布を用いて、ぶれやぼけのない被写体の画像を得ることができる。
【００２４】
上記の撮像装置においては、前記画像処理手段が、点像分布関数と被写体画像データに基づいて改善された被写体画像データを生成する際に、被写体画像データから劣化画像の分布を特定する処理と、点像分布関数から光学的伝達関数を特定する処理と、原画像の最初の推定分布を特定する処理と、（１）原画像の推定分布をフーリエ変換して第１の関数を得る処理と、（２）前記第１の関数に前記光学的伝達関数を乗じて第２の関数を得る処理と、（３）前記第２の関数を逆フーリエ変換して第３の関数を得る処理と、（４）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る処理と、（５）前記第４の関数をフーリエ変換して第５の関数を得る処理と、（６）前記第５の関数に前記光学的伝達関数の反転関数を乗じて第６の関数を得る処理と、（７）前記第６の関数を逆フーリエ変換して第７の関数を得る処理と、（８）原画像の推定分布に前記第７の関数を乗じて、原画像の次の推定分布を得る処理と、原画像の次の推定分布を原画像の推定分布に置き換えて、前記（１）から（８）の処理を繰返す処理と、原画像の推定分布に基づいて、改善された被写体画像データを生成する処理を実行することが好ましい。
【００２５】
劣化画像から原画像を復元する際に、ＯＴＦ分布Ｈ（ｓ、ｔ）が０となる周波数領域が存在する場合には、数式（６）の演算を行うことができず、原画像を完全に復元することができない。そこで上記の撮像装置では、原画像の分布ｆｒ（ｘ、ｙ）と劣化画像の分布ｇｒ（ｘ、ｙ）を確率密度関数として扱い、Ｂａｙｅｓの理論に基づいて原画像の推定を行う。原画像分布ｆｒ（ｘ、ｙ）と劣化画像分布ｇｒ（ｘ、ｙ）は、下記の正規化を行うことによって、確率密度関数として扱うことが可能になる。
【００２６】
【数８】
【数９】
【００２７】
上記の正規化に合わせて、ＯＴＦ分布Ｈｒ（ｓ、ｔ）についても正規化する。Ｈｒ（ｓ、ｔ）は、空間周波数の原点での値を用いて正規化する。
【００２８】
【数１０】
【００２９】
正規化された原画像の分布ｆ（ｘ、ｙ）と、劣化画像の分布ｇ（ｘ、ｙ）は、非負の関数であり、定義された領域での積分値が１であるから、確率密度関数として扱うことができる。上記の場合において、ｆ（ｘ、ｙ）は原画像の座標（ｘ、ｙ）における結像という事象の確率密度関数である。ｇ（ｘ、ｙ）は劣化画像の座標（ｘ、ｙ）における結像という事象の確率密度関数である。
【００３０】
原画像および劣化画像の分布を、確率密度関数と見なすことが可能な場合、Ｂａｙｅｓの理論に基づいて、劣化画像の分布から、その劣化画像を生じさせている原画像の分布を推定することができる。
原画像の座標（ｘ、ｙ）に光が結像する事象は、その座標（ｘ、ｙ）に点光源が存在する事象として扱うことができる。原画像の座標（ｘ１、ｙ１）において点光源が存在する事象をＶ（ｘ１、ｙ１）、劣化画像の座標（ｘ２、ｙ２）において点像が結像する事象をＡ（ｘ２、ｙ２）とした場合、それぞれの事象の確率Ｐ（Ｖ）およびＰ（Ａ）は以下で表現される。
【００３１】
【数１１】
【数１２】
【００３２】
また、原画像の座標（ｘ１、ｙ１）に点光源が存在している場合に、劣化画像の座標（ｘ２、ｙ２）に結像する確率は、事象Ｖ（ｘ１、ｙ１）の生起を条件とする事象Ａ（ｘ２、ｙ２）の生起確率である。上記の確率はＰＳＦ分布ｈ（ｘ、ｙ）を用いて以下で表現される。
【００３３】
【数１３】
【００３４】
Ｂａｙｅｓの理論に基づくと、劣化画像内の点（ｘ２、ｙ２）に点像を結像させる場合の原画像の分布Ｐ（Ｖ（ｘ、ｙ）｜Ａ（ｘ２、ｙ２））は、以下で推定される。
【００３５】
【数１４】
【００３６】
上式の右辺のＰ（Ｖ（ｘ、ｙ））、Ｐ（Ａ（ｘ２、ｙ２）｜Ｖ（ｘ、ｙ））に、数式（１１）および数式（１３）を代入すると、次式を得る。
【００３７】
【数１５】
【００３８】
上式の左辺は、劣化画像において点像が結像している場合に、推定される原画像の分布を表している。上式に劣化画像の分布ｇ（ｘ、ｙ）を乗じて積分することによって、劣化画像の分布ｇ（ｘ、ｙ）を実現するための原画像の分布ｆ（ｘ、ｙ）を得ることができる。
上式の両辺に、Ｐ（Ａ（ｘ２、ｙ２））＝ｇ（ｘ２、ｙ２）を乗じて、すべての（ｘ２、ｙ２）に関して積分すると、次式が得られる。
【００３９】
【数１６】
【００４０】
上式の左辺に数式（１２）を代入すると、その積分の結果はＰ（Ｖ（ｘ、ｙ））であり、ｆ（ｘ、ｙ）に等しい。従って、Ｂａｙｅｓの理論に基づくと、以下の関係が成り立つ。
【００４１】
【数１７】
【００４２】
上記の関係は、分布ｆ（ｘ、ｙ）が真の原画像の分布である場合に成立すると考えられる。すなわち、上式を満たす分布ｆ（ｘ、ｙ）を算出することが、劣化画像の復元に相当する。
上記の関係を満たす分布ｆ（ｘ、ｙ）は、数式（１７）の右辺のｆ（ｘ、ｙ）をｆｋ（ｘ、ｙ）とし、数式（１７）の左辺のｆ（ｘ、ｙ）をｆｋ＋１（ｘ、ｙ）として、ｆｋ（ｘ、ｙ）に関する反復計算を実施し、ｆｋ（ｘ、ｙ）の収束値を求めることで算出することができる。上記の反復計算によって求まるｆｋ（ｘ、ｙ）の収束値は、Ｂａｙｅｓの理論に基づく原画像の推定分布に相当する。
【００４３】
上記では、原画像の分布ｆ（ｘ、ｙ）や劣化画像の分布ｇ（ｘ、ｙ）を正規化した場合について説明したが、実際の反復計算を実施する上では、これらの分布は正規化することなく、そのまま使用することができる。
【００４４】
上記の反復計算においては、反復計算を実施する前に、原画像の最初の推定分布ｆ０（ｘ、ｙ）を設定しておく。最初の推定分布ｆ０（ｘ、ｙ）としては、任意の分布を設定することができる。一般的には、劣化画像の分布ｇ（ｘ、ｙ）は原画像の分布ｆ（ｘ、ｙ）から大きく異なることはないため、最初の推定分布ｆ０（ｘ、ｙ）としては、劣化画像の分布ｇ（ｘ、ｙ）を用いることが好ましい。
【００４５】
数式（１７）の右辺は、ＰＳＦ分布ｈ（ｘ、ｙ）を用いた畳み込み積分を備えている。一般に、光学系の位相特性まで含めてＰＳＦを正確に評価することは困難であり、上記の反復計算を正確な位相特性を含めて実施することは困難である。正確な位相特性を含まないＰＳＦを用いた反復計算は、誤った収束の結果をもたらすため、原画像の正確な復元の妨げとなる。
そこで、ＰＳＦの代わりに正確な位相特性を含ませることが容易であるＯＴＦを用いることで、より正確に原画像を復元することが可能となる。また、復元の過程で位相特性を正確に評価するために、原画像の推定分布ｆｋ（ｘ、ｙ）（ｋ＝０，１，２，・・・）を複素関数に拡張し、その実部が画像分布を表現するものとして取扱う。上式の右辺に、フーリエ変換と逆フーリエ変換を用いることで、ＯＴＦを用いた形式に変更することができる。フーリエ変換をＦＴ（）、逆フーリエ変換をＦＴ−１（）で表現すると、数式（１７）は以下で表現される。
【００４６】
【数１８】
【００４７】
上記の反復計算を、ｆｋが収束するまで繰り返し実施することによって、原画像を推定することができる。ｆｋが収束したか否かの判定は、例えば反復計算の回数を予め設定しておいて、反復計算の回数によって判定してもよいし、ｆｋとｆｋ＋１の差分を算出して、算出される差分の絶対値の総和があるしきい値以下となるか否かで判断してもよい。
【００４８】
上記した原理を用いる画像復元処理は、以下に示す工程を順次実施していくことで実現される。以下では繰り返し計算によって推定される原画像の分布をｆｋ（ｋ＝０、１、・・・）とする。原画像の推定分布ｆｋは、復元の過程で位相に関する特性を正確に評価するために、複素関数として扱う。
まずＰＳＦ分布ｈ（ｘ、ｙ）を二次元フーリエ変換して、ＯＴＦ分布Ｈ（ｓ、ｔ）を特定する。
次に原画像の最初の推定分布として、ｆ０（ｘ、ｙ）の実部をｇ（ｘ、ｙ）として、ｆ０（ｘ、ｙ）の虚部を０とする。
次に以下に示す演算を、ｆｋ（ｘ、ｙ）が収束するまで繰り返し実施する。ここでＦＴ（）は二次元のフーリエ変換を表し、ＦＴ−１（）は二次元の逆フーリエ変換を表す。またＨ＃（ｓ、ｔ）は、Ｈ（ｓ、ｔ）の反転関数であり、Ｈ＃（ｓ、ｔ）＝Ｈ（−ｓ、−ｔ）である。
【００４９】
【数１９】
【数２０】
【数２１】
【００５０】
上記の反復計算を繰返し実施し、最終的な原画像の推定分布ｆｋを得る。最終的に得られる原画像の推定分布ｆｋ（ｘ、ｙ）が、原画像の復元画像ｆ（ｘ、ｙ）に相当する。
【００５１】
以上のように、上記の画像取得装置によれば、劣化画像の分布ｇ（ｘ、ｙ）とＰＳＦの分布ｈ（ｘ、ｙ）を用いて、Ｂａｙｅｓの理論に基づき、劣化画像の分布ｇ（ｘ、ｙ）を実現する原画像の分布ｆ（ｘ、ｙ）のうちで最大に尤もらしい分布を、反復計算によって推定することができる。この画像取得装置によれば、ＯＴＦの分布Ｈ（ｓ、ｔ）がある周波数領域で０となる場合、すなわち数式（６）による画像復元ができない場合であっても、原画像を復元することができる。ぼけやぶれのない被写体の画像を得ることができる。
【００５２】
あるいは、上記の撮像装置は、前記画像処理手段が、点像分布関数と被写体画像データに基づいて改善された被写体画像データを生成する際に、被写体画像データから劣化画像の分布を特定する処理と、原画像の最初の推定分布を特定する処理と、点像画像データに基いて取得された点像分布関数を、点像分布関数の最初の推定分布に設定する処理と、（Ａ）前記劣化画像の分布と、原画像の推定分布と、点像分布関数の推定分布に基づいて、点像分布関数の推定分布を更新する処理と、（Ｂ）前記劣化画像の分布と、原画像の推定分布と、点像分布関数の推定分布に基づいて、原画像の推定分布を更新する処理と、上記（Ａ）と（Ｂ）の工程を交互に繰返す処理と、原画像の推定分布に基づいて、改善された被写体画像データを生成する処理を実行するものであり、前記（Ａ）の点像分布関数の推定分布を更新する処理は、（Ａ１）原画像の推定分布をフーリエ変換して、原画像の推定分布のスペクトル分布を得る処理と、（Ａ２）点像分布関数の推定分布をフーリエ変換して第１の関数を得る処理と、（Ａ３）前記第１の関数に前記原画像の推定分布のスペクトル分布を乗じて第２の関数を得る処理と、（Ａ４）前記第２の関数を逆フーリエ変換して第３の関数を得る処理と、（Ａ５）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る処理と、（Ａ６）前記第４の関数をフーリエ変換して第５の関数を得る処理と、（Ａ７）前記第５の関数に前記原画像の推定分布のスペクトル分布の反転関数を乗じて第６の関数を得る処理と、（Ａ８）前記第６の関数を逆フーリエ変換して第７の関数を得る処理と、（Ａ９）前記点像分布関数の推定分布に前記第７の関数を乗じて、点像分布関数の次の推定分布を得る処理と、（Ａ１０）点像分布関数の次の推定分布を点像分布関数の推定分布に置き換える処理とを備え、前記（Ｂ）の原画像の推定分布を更新する処理は、（Ｂ１）点像分布関数の推定分布をフーリエ変換して、光学的伝達関数の推定分布を得る処理と、（Ｂ２）原画像の推定分布をフーリエ変換して第１の関数を得る処理と、（Ｂ３）前記第１の関数に前記光学的伝達関数の推定分布を乗じて第２の関数を得る処理と、（Ｂ４）前記第２の関数を逆フーリエ変換して第３の関数を得る処理と、（Ｂ５）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る処理と、（Ｂ６）前記第４の関数をフーリエ変換して第５の関数を得る処理と、（Ｂ７）前記第５の関数に前記光学的伝達関数の推定分布の反転関数を乗じて第６の関数を得る処理と、（Ｂ８）前記第６の関数を逆フーリエ変換して第７の関数を得る処理と、（Ｂ９）原画像の推定分布に前記第７の関数を乗じて、原画像の次の推定分布を得る処理と、（Ｂ１０）原画像の次の推定分布を原画像の推定分布に置き換える処理とを備えることが好ましい。
【００５３】
第２撮像手段で撮像された点像画像データから取得されるＰＳＦ分布は、第２撮像手段による撮像の際の光学的特性を反映したものであり、第１撮像手段による撮像の際の光学的特性と完全に一致している保証はない。従って、ＰＳＦ分布のわずかな相違によって、原画像を完全に復元することが困難な場合がある。そこで、上記の撮像装置では、第２撮像手段で撮像された点像画像データから取得されるＰＳＦ分布を、第１撮像手段による撮像の際の真のＰＳＦ分布の最初の推定値として扱い、原画像の分布とＰＳＦ分布の双方について、Ｂａｙｅｓの理論に基いた推定を行う。
【００５４】
上述したように、劣化画像分布ｇ（ｘ、ｙ）とＯＴＦ分布Ｈ（ｓ、ｔ）が既知であれば、未知である原画像分布ｆ（ｘ、ｙ）を復元することができる。上記と同様の原理によって、ＯＴＦ分布Ｈ（ｓ、ｔ）が未知である場合に、劣化画像分布ｇ（ｘ、ｙ）と原画像分布ｆ（ｘ、ｙ）が既知であれば、ＯＴＦ分布Ｈ（ｓ、ｔ）を復元することが可能となる。
【００５５】
数式（１５）において、ｘ２−ｘを新たにｘ’とし、ｙ２−ｙを新たにｙ’とすると、次式が得られる。
【００５６】
【数２２】
【００５７】
上式の両辺にＰ（Ａ（ｘ２、ｙ２））を乗じると、次式が得られる。
【００５８】
【数２３】
【００５９】
上式の左辺は、Ｂａｙｅｓの理論に基づくと、Ｐ（Ａ（ｘ２、ｙ２）｜Ｖ（ｘ２−ｘ’、ｙ２−ｙ’））×Ｐ（Ｖ（ｘ２−ｘ’、ｙ２−ｙ’））に等しいため、次式となる。
【００６０】
【数２４】
【００６１】
上式に、数式（１１）、数式（１２）および数式（１３）を代入すると、次式が得られる。
【００６２】
【数２５】
【００６３】
上式の両辺を、（ｘ２、ｙ２）に関して積分すると、次式が得られる。
【００６４】
【数２６】
【００６５】
上式左辺の積分は１に等しいため、結局、Ｂａｙｅｓの理論に基づくと、以下の関係が成立する。
【００６６】
【数２７】
【００６７】
上式の関係を満たすＰＳＦ分布ｈ（ｘ、ｙ）が算出されると、算出されたｈ（ｘ、ｙ）をフーリエ変換することによって、ＯＴＦ分布Ｈ（ｓ、ｔ）を推定することができる。
【００６８】
数式（２７）の右辺のｈ（ｘ、ｙ）をｈｋ（ｘ、ｙ）とし、数式（２７）の左辺のｈ（ｘ、ｙ）をｈｋ＋１（ｘ、ｙ）として、ｈｋ（ｘ、ｙ）に関する反復計算を実施し、ｈｋ（ｘ、ｙ）の収束値を求める。この反復計算は、次式で表される。
【００６９】
【数２８】
【００７０】
上記の反復計算によって求まるｈｋ（ｘ、ｙ）の収束値をフーリエ変換して、Ｂａｙｅｓの理論に基づくＯＴＦの推定分布Ｈｋ（ｘ、ｙ）を得ることができる。この関係は次式で表される。
【００７１】
【数２９】
【００７２】
上記した原理を用いるＯＴＦの推定処理は、以下に示す工程を順次実施していくことで実現される。以下では繰り返し計算によって推定されるＰＳＦの分布をｈｋ（ｋ＝０、１、・・・）とする。ＰＳＦの推定分布ｈｋは、反復計算の過程で位相に関する特性を正確に評価するために、複素関数として扱う。
まず劣化画像の分布ｇ（ｘ、ｙ）と、原画像の分布ｆ（ｘ、ｙ）を特定する。本発明の方法では、ｇ（ｘ、ｙ）の実部を劣化画像における照度の分布として、ｇ（ｘ、ｙ）の虚部を０とする。またｆ（ｘ、ｙ）の実部を原画像における照度の分布として、ｆ（ｘ、ｙ）の虚部を０とする。
次に原画像の分布ｆ（ｘ、ｙ）をフーリエ変換して、原画像の分布のスペクトル分布Ｆ（ｓ、ｔ）を算出する。
次に、ＰＳＦの最初の推定分布ｈ０（ｘ、ｙ）を特定する。
次に以下に示す演算を、ｈｋ（ｘ、ｙ）が収束するまで繰り返し実施する。ここでＦＴ（）は二次元のフーリエ変換を表し、ＦＴ−１（）は二次元の逆フーリエ変換を表す。またＦ＃（ｓ、ｔ）は、Ｆ（ｓ、ｔ）の反転関数であり、Ｆ＃（ｓ、ｔ）＝Ｆ（−ｓ、−ｔ）である。
【００７３】
【数３０】
【数３１】
【数３２】
【００７４】
ｈｋ（ｘ、ｙ）の収束の判定は、例えば反復計算の回数が予め定められた回数に到達したか否かを基準に判別してもよいし、ｈｋ（ｘ、ｙ）とｈｋ＋１（ｘ、ｙ）の差分の分布を算出して、その差分の絶対値が全ての（ｘ、ｙ）に対してあるしきい値以下となるか否かを基準に判別してもよいし、ｆｋとｆｋ＋１の差分を算出して、算出される差分の絶対値を（ｘ、ｙ）に関して積分した値があるしきい値以下となるか否かで判断してもよい。
ＰＳＦ分布ｈｋ（ｘ、ｙ）が収束したら、それをフーリエ変換して、ＯＴＦの推定分布Ｈｋ（ｓ、ｔ）を算出する。
【００７５】
上記したように、劣化画像とＯＴＦとが既知であれば原画像を復元することが可能であり、劣化画像と原画像が既知であればＯＴＦを推定することが可能である。
これらを組み合わせることによって、劣化画像のみに基づいて、ＯＴＦを推定し、原画像を復元することが可能となる。
すなわち、原画像とＰＳＦのそれぞれについて、適当な分布を仮定しておき、上記したＰＳＦの推定（すなわち、ＯＴＦの推定）と原画像の推定とを交互に繰返し実施していくことによって、原画像を復元する。
ＯＴＦに関しては、劣化画像の分布と、原画像の推定分布を用いて、Ｂａｙｅｓの理論に基づく計算によって、より改善された推定分布を得ることができる。ＯＴＦの推定に用いる原画像の推定分布が、真の原画像の分布に近いほど、真のＯＴＦの分布に近い推定分布を得ることができる。
原画像に関しては、劣化画像の分布と、ＯＴＦの推定分布を用いて、Ｂａｙｅｓの理論に基づく計算によって、より改善された推定分布を得ることができる。原画像の推定に用いるＯＴＦの推定分布が、真のＯＴＦの分布に近いほど、真の原画像の分布に近い推定分布を得ることができる。
従って、上記したＯＴＦの推定と原画像の推定を交互に繰返し実施していくことによって、原画像の推定分布はより真の原画像の分布に近づいていき、ＯＴＦの推定分布はより真のＯＴＦの分布に近づいていき、結果として良好に復元された原画像と良好に推定されたＯＴＦの双方を得ることができる。
【００７６】
すなわち、以下の工程を順に実施していくことによって、劣化画像のみから原画像を復元することができる。
まず原画像の最初の推定分布と、ＰＳＦの最初の推定分布を設定する。一般に、劣化画像の分布は原画像の分布から大きく異なることはないため、原画像の最初の推定分布としては劣化画像の分布を用いる。ＰＳＦの最初の推定分布としては、第２撮像手段によって撮像された点像画像データから取得されるＰＳＦ分布を用いる。
次に、数式（３０）〜（３２）の計算を実施して、より改善されたＰＳＦの推定分布ｈｋ（ｓ、ｔ）を取得して、得られたｈｋ（ｓ、ｔ）をフーリエ変換してＯＴＦの推定分布Ｈｋ（ｓ、ｔ）を取得する。上記の推定を実施することによって、より真のＯＴＦ分布に近いＯＴＦの推定分布を得ることができる。
次に、数式（１９）〜（２１）の計算を実施して、より改善された原画像の推定分布ｆｋ（ｘ、ｙ）を取得する。上記の推定を実施することによって、より真の原画像分布に近い原画像の推定分布を得ることができる。
上記したＯＴＦの推定と原画像の推定を、交互に繰返し実施していくことによって、推定される原画像の分布は真の原画像の分布に近づいていき、推定されるＯＴＦの分布は真のＯＴＦの分布に近づいていく。従って、上記の反復計算を実施することによって、原画像を復元することが可能となる。
【００７７】
上記の画像取得装置によれば、第１撮像手段による撮像時のＰＳＦと第２撮像手段による撮像時のＰＳＦが相違する場合であっても、ぶれのない被写体の画像を取得することができる。
【００７８】
本発明は天体撮像装置としても具現化することができる。その天体撮像装置は、天体を撮像して天体画像データを出力する第１撮像手段と、電離層に向けてレーザ光を照射するレーザ光照射手段と、レーザ光照射手段から照射されて電離層で反射したレーザ光を撮像して点像画像データを出力する第２撮像手段と、天体画像データと点像画像データに基いて改善された天体画像データを生成する画像処理手段を備えている。その天体撮像装置において、第１撮像手段と第２撮像手段は、互いに同期して撮像するものである。その天体撮像装置において、画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と天体画像データに基いて改善された天体画像データを生成する。
【００７９】
上記の天体撮像装置の動作内容を方法の発明として具現化することもできる。本発明の天体撮像方法は、第１撮像手段によって天体を撮像して天体画像データを出力する工程と、レーザ光照射手段によって電離層に向けてレーザ光を照射する工程と、レーザ光照射手段から照射されて電離層で反射したレーザ光を第２撮像手段によって撮像して点像画像データを出力する工程と、画像処理手段によって天体画像データと点像画像データに基いて改善された天体画像データを生成する工程を備えている。その天体撮像方法において、第１撮像手段と第２撮像手段は、互いに同期して撮像するものである。その天体撮像方法において、画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と天体画像データに基いて改善された天体画像データを生成する。
なお本発明は上記の天体撮像方法における各工程をコンピュータに実行させるためのプログラムとしても具現化することができる。
【００８０】
第１撮像手段で撮像される天体の画像は、大気のゆれによるシーイングの影響を受けて、ぼやけた画像となってしまう。本発明の天体撮像装置では、第１撮像手段で撮像された天体の画像と、第２撮像手段で撮像されたレーザ光の画像を用いて、シーイングの影響を除去して改善された天体の画像を取得する。第２撮像手段で撮像される画像は、レーザ光照射手段から電離層へ到達するまでの間にシーイングの影響を受け、さらに電離層で反射してから第２撮像手段へ到達するまでの間にシーイングの影響を受けた点光源の画像として扱うことができる。従って、第２撮像手段で撮像された点像画像データに基づいて、天体からの光が電離層を経由して第１撮像手段に到達するまでの間に受けるシーイングの影響を表現するＰＳＦ分布を取得することができる。第１撮像手段で撮像された天体画像データとＰＳＦ分布を用いて、シーイングの影響を排除した改善された天体画像データを得ることができる。
【００８１】
また、本発明者らは、ぶれ画像に特有の性質に着目して、ぶれ画像のみに基づいてぶれのない画像を復元する手法にも想到した。本発明のぶれ画像復元装置は、ぶれ画像の分布を特定する手段と、ぶれ画像の分布についての複素ケプストラムを算出する手段と、ぶれ画像の分布についての複素ケプストラムからパルス列を抽出する手段と、抽出されたパルス列をぶれ画像の分布についての複素ケプストラムから除去して、ぶれのない画像の分布についての複素ケプストラムを推定する手段と、ぶれのない画像の分布についての複素ケプストラムからぶれのない画像の分布を算出する手段を備えている。
【００８２】
上記のぶれ画像復元装置の動作内容を方法の発明として具現化することもできる。本発明のぶれ画像復元方法は、ぶれ画像の分布を特定する工程と、ぶれ画像の分布についての複素ケプストラムを算出する工程と、ぶれ画像の分布についての複素ケプストラムからパルス列を抽出する工程と、抽出されたパルス列をぶれ画像の分布についての複素ケプストラムから除去して、ぶれのない画像の分布についての複素ケプストラムを推定する工程と、ぶれのない画像の分布についての複素ケプストラムからぶれのない画像の分布を算出する工程を備えている。
なお本発明は上記のぶれ画像復元方法における各工程をコンピュータに実行させるためのプログラムとしても具現化することができる。
【００８３】
複素ケプストラムとは、ある関数について、フーリエ変換してスペクトル分布を取得し、そのスペクトル分布の複素対数関数を計算して、その結果を逆フーリエ変換したものをいう。これはフーリエ変換を続けて２回行ったもので代用することも出来る。以下ではぶれ画像を劣化画像と呼び、ぶれのない画像を原画像と呼ぶ。劣化画像の分布ｇ（ｘ、ｙ）、ＰＳＦの分布ｈ（ｘ、ｙ）、原画像の分布ｆ（ｘ、ｙ）には、以下の関係が成り立つ。
【００８４】
【数３３】
【００８５】
上式の両辺をフーリエ変換すると、劣化画像のスペクトル分布Ｇ（ｓ、ｔ）、ＯＴＦの分布Ｈ（ｓ、ｔ）、原画像のスペクトル分布Ｆ（ｓ、ｔ）の間で以下の関係が成り立つ。
【００８６】
【数３４】
【００８７】
上式の両辺について複素対数関数を計算すると、次式を得る。
【００８８】
【数３５】
【００８９】
上式の両辺について逆フーリエ変換することで、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）、ＰＳＦの複素ケプストラムｈｃ（ｘ、ｙ）、原画像の複素ケプストラムｆｃ（ｘ、ｙ）について、以下の関係式を得る。
【００９０】
【数３６】
【００９１】
一般のぶれ画像は、多数個のずれた画像が重なり合ったものとして扱うことができる。原画像の分布ｆ（ｘ、ｙ）に対して、全く同じ画像を、強度をｎ倍として、位置を（α、β）だけずらして重ね合わせた場合、劣化画像の分布ｇ（ｘ、ｙ）は次式で与えられる。
【００９２】
【数３７】
【００９３】
上式の両辺をフーリエ変換すると、劣化画像の分布のスペクトル分布Ｇ（ｓ、ｔ）と原画像の分布のスペクトル分布Ｆ（ｘ、ｙ）の間に、次の関係が成り立つ。
【００９４】
【数３８】
【００９５】
上式の両辺の複素対数関数を計算すると、次式が得られる。
【００９６】
【数３９】
【００９７】
ここで、ｎは１に比べて十分に小さいと仮定している。上式の両辺について逆フーリエ変換すると、次式を得る。
【００９８】
【数４０】
【００９９】
上述のように、上式の左辺ｇｃ（ｘ、ｙ）は、劣化画像の分布ｇ（ｘ、ｙ）についての複素ケプストラムである。上式の右辺第１項ｆｃ（ｘ、ｙ）は、原画像の分布ｆ（ｘ、ｙ）についての複素ケプストラムである。上式の右辺第２項以降は、ＰＳＦ分布ｈ（ｘ、ｙ）についての複素ケプストラムである。
上式から明らかなように、劣化画像の分布ｇ（ｘ、ｙ）の複素ケプストラムｇｃ（ｘ、ｙ）は、原画像の分布ｆ（ｘ、ｙ）の複素ケプストラムｆｃ（ｘ、ｙ）に、ＰＳＦ分布ｈ（ｘ、ｙ）の複素ケプストラムを重ね合わせたものとなっている。画像の劣化がぶれによるものの場合には、ＰＳＦ分布ｈ（ｘ、ｙ）の複素ケプストラムは、ずれた画像の強度ｎと、ずれ量（α、β）で定まるデルタ関数の系列（すなわちパルス列）によって表現される。従って、劣化画像の分布ｇ（ｘ、ｙ）の複素ケプストラムｇｃ（ｘ、ｙ）からこのパルス列を抽出して取り除くことができれば、原画像の分布ｆ（ｘ、ｙ）の複素ケプストラムｆｃ（ｘ、ｙ）を推定することができる。
【０１００】
劣化画像の分布ｇ（ｘ、ｙ）の複素ケプストラムｇｃ（ｘ、ｙ）からパルス列を抽出するには、例えば劣化画像の分布ｇ（ｘ、ｙ）の複素ケプストラムｇｃ（ｘ、ｙ）の振幅値がしきい値を超える個所をパルス列として識別する。パルス列が抽出されると、例えばノッチフィルタを用いることで、劣化画像の分布ｇ（ｘ、ｙ）の複素ケプストラムｇｃ（ｘ、ｙ）からパルス列を除去して、原画像の分布ｆ（ｘ、ｙ）の複素ケプストラムｆｃ（ｘ、ｙ）を推定することができる。
【０１０１】
推定された原画像の分布ｆ（ｘ、ｙ）の複素ケプストラムｆｃ（ｘ、ｙ）に関して、フーリエ変換を実施し、その結果の複素指数関数を計算することで、原画像のスペクトル分布Ｆ（ｓ、ｔ）が得られる。その後に原画像のスペクトル分布Ｆ（ｓ、ｔ）を逆フーリエ変換することで、原画像の分布ｆ（ｘ、ｙ）を復元することができる。
【０１０２】
なお上記では説明を簡単にするために２つの画像が重なり合ったぶれ画像の場合を例として説明したが、多数個の画像が重なり合っている一般的なぶれ画像の場合にも、同様の手法によって原画像の分布ｆ（ｘ、ｙ）を復元することができる。
【０１０３】
本発明の他の一つのぶれ画像復元装置は、ぶれ画像の分布を特定する手段と、ぶれ画像の分布についての複素ケプストラムを算出する手段と、ぶれ画像の分布についての複素ケプストラムからパルス列を抽出する手段と、抽出されたパルス列に基づいて点像分布関数を取得する手段と、ぶれ画像の分布と点像分布関数に基づいてぶれのない画像の分布を推定する手段を備えている。
【０１０４】
上記のぶれ画像復元装置の動作内容を方法の発明として具現化することもできる。本発明の他の一つのぶれ画像復元方法は、ぶれ画像の分布を特定する工程と、ぶれ画像の分布についての複素ケプストラムを算出する工程と、ぶれ画像の分布についての複素ケプストラムからパルス列を抽出する工程と、抽出されたパルス列に基づいて点像分布関数を取得する工程と、ぶれ画像の分布と点像分布関数に基づいてぶれのない画像の分布を推定する工程を備えている。
なお本発明は上記のぶれ画像復元方法における各工程をコンピュータに実行させるためのプログラムとしても具現化することができる。
【０１０５】
ここでもう一度、原画像の分布ｆ（ｘ、ｙ）が撮像時のぶれによって劣化画像の分布ｇ（ｘ、ｙ）となっている場合を考える。原画像の分布ｆ（ｘ、ｙ）に対して、全く同じ画像を強度をｎ倍として、位置を（α、β）だけずらして重ね合わせた場合、劣化画像の分布ｇ（ｘ、ｙ）は次式で与えられる。
【０１０６】
【数４１】
【０１０７】
すなわち、この場合のＰＳＦ分布ｈ（ｘ、ｙ）は、次式で与えられる。
【０１０８】
【数４２】
【０１０９】
数式（４０）と数式（４２）を対比すると、数式（４２）のＰＳＦの分布ｈ（ｘ、ｙ）の右辺第２項は、数式（４０）の劣化画像の分布ｇ（ｘ、ｙ）の複素ケプストラムｇｃ（ｘ、ｙ）に現れるパルス列のうち最大のもの（数式（４０）の右辺第２項）と同様に、位置が（α、β）であって、大きさがｎのデルタ関数となっていることが分かる。従って、劣化画像の分布ｇ（ｘ、ｙ）の複素ケプストラムｇｃ（ｘ、ｙ）から抽出されたパルス列のうち最大のものに、強度が１で原点に位置するデルタ関数を重ね合わせることで、ぶれの影響を示すＰＳＦの分布ｈ（ｘ、ｙ）を直接的に取得することができる。
【０１１０】
あるいは、ｎが１より十分に小さければ、劣化画像の分布ｇ（ｘ、ｙ）の複素ケプストラムｇｃ（ｘ、ｙ）に現れるパルス列に、強度が１で原点に位置するデルタ関数を単に重ね合わせることで、ＰＳＦの分布ｈ（ｘ、ｙ）の近似解を取得することができる。
【０１１１】
ぶれの影響を表現するＰＳＦの分布ｈ（ｘ、ｙ）が取得されれば、劣化画像の分布ｇ（ｘ、ｙ）とＰＳＦの分布ｈ（ｘ、ｙ）に基づいて、原画像の分布ｆ（ｘ、ｙ）を推定することができる。ぶれ画像のみに基づいて、ぶれのない画像を取得することができる。
【０１１２】
なお上記では説明を簡単にするために２つの画像が重なり合ったぶれ画像の場合を例として説明したが、多数個の画像が重なり合っている一般的なぶれ画像の場合にも、同様の手法によって原画像の分布ｆ（ｘ、ｙ）を復元することができる。
【発明の効果】
【０１１３】
本発明の装置、方法およびプログラムによれば、撮像時のぶれやぼけの影響を表現する点像分布関数を適切に取得して、劣化した画像から劣化のない画像を復元することができる。
【図面の簡単な説明】
【０１１４】
【図１】図１は実施例１の撮像装置１００の構成を示す図である。
【図２】図２は実施例１の撮像装置１００の動作を説明するフローチャートである。
【図３】図３は実施例１の撮像装置１００によってぶれのない被写体画像を取得する様子を示している。
【図４】図４は実施例１の変形例１の撮像装置１５０における画像復元処理を説明するフローチャートである。
【図５】図５は実施例１の変形例２の撮像装置１６０における画像復元処理を説明するフローチャートである。
【図６】図６は実施例１の変形例２の撮像装置１６０におけるＰＳＦ分布の推定処理を説明するフローチャートである。
【図７】図７は実施例１の変形例２の撮像装置１６０における原画像分布の推定処理を説明するフローチャートである。
【図８】図８は実施例２の天体撮像装置８００の構成を示す図である。
【図９】図９は実施例２の天体撮像装置８００によってシーイングの影響を排除した天体画像を取得する様子を示す図である。
【図１０】図１０は実施例２の天体撮像装置８００の動作を説明するフローチャートである。
【図１１】図１１は実施例３の画像復元装置１１００の構成を示す図である。
【図１２】図１２は実施例３の画像復元装置１１００の処理を説明するフローチャートである。
【図１３】図１３は実施例３の画像復元装置１１００で計算されるぶれ画像の複素ケプストラム分布の概形を示す図である。
【図１４】図１４は実施例３の画像復元装置１１００の機能ブロック図である。
【図１５】図１５は実施例４の画像復元装置１１５０の動作を説明するフローチャートである。
【図１６】図１６は実施例４の画像復元装置１１５０の機能ブロック図である。
【発明を実施するための形態】
【０１１５】
以下、本発明を具現化した実施例について図面を参照して説明する。
【０１１６】
（実施例１）
図面を参照しながら、本実施例の撮像装置１００について説明する。
図１は本実施例の撮像装置１００の構成を示す概略図である。本実施例の撮像装置１００は、撮像台１１０に載置された被写体Ｗを撮像する第１イメージセンサ１０４と、撮像台１１０に固定されたレーザ光源１０８と、レーザ光源１０８を撮像する第２イメージセンサ１０２と、制御装置１１２を備えている。本実施例では被写体Ｗは表面に半導体回路が形成されたフォトマスクであって、撮像装置１００は被写体Ｗの表面に形成された回路の形状を撮像する。
【０１１７】
第１イメージセンサ１０４と第２イメージセンサ１０２は、撮像台１１０の上面に対して平行な撮像面を備えている。第１イメージセンサ１０４と第２イメージセンサ１０２は、それぞれ撮像面に複数の撮像素子を備えており、撮像面に形成される像を撮像素子を介して認識し、ＲＧＢの各色成分の照度の分布を示す電気信号に変換して画像処理装置１１２へ送信する。第１イメージセンサ１０４と第２イメージセンサ１０２は、両者が一体的に移動するように、センサ支持台１０６に搭載されている。第１イメージセンサ１０４と第２イメージセンサ１０２は、レンズ収差等の特性が同一である。また、第１イメージセンサ１０４と第２イメージセンサ１０２は、それぞれのシャッタが同期して開閉するように制御装置１１２によって制御される。従って、撮像台１１０とセンサ支持台１０６の間のわずかな振動によって撮像時に第１イメージセンサ１０４で撮像される被写体Ｗの画像がぶれてしまう場合には、第２イメージセンサ１０２で撮像されるレーザ光源１０８の画像も同じようにぶれることになる。
【０１１８】
第２イメージセンサ１０２で撮像されるレーザ光源１０８の画像は、撮像台１１０に固定された点光源を撮像した画像ということができる。従って、第２イメージセンサ１０２で撮像される画像は、第２イメージセンサ１０２で撮像する際の点像分布関数（ＰＳＦ：Point Spread Function）の分布として扱うことができる。第２イメージセンサ１０２と第１イメージセンサ１０４は同一の特性を有しており、それぞれのシャッタが同期して開閉するように撮像するので、第２イメージセンサ１０２で撮像する際のＰＳＦの分布は、第１イメージセンサ１０４で撮像する際のＰＳＦの分布として扱うことができる。レーザ光源１０８の発光は制御装置１１２によって制御される。
【０１１９】
制御装置１１２は、第１イメージセンサ１０４から送信されるデータと、第２イメージセンサ１０２から送信されるデータに基づいて、被写体Ｗのぶれのない画像を生成する。
【０１２０】
以下では図２のフローチャートを参照しながら、撮像装置１００によって被写体Ｗのぶれのない画像を取得する方法について説明する。以下の説明では、撮像時にぶれてしまった被写体Ｗの画像を劣化画像と呼び、ぶれのない被写体Ｗの画像を原画像と呼ぶ。
【０１２１】
ステップＳ２０２では、被写体Ｗを撮像台１１０に載置し、位置を固定する。
【０１２２】
ステップＳ２０４では、点光源であるレーザ光源１０８の発光を開始する。
【０１２３】
ステップＳ２０６では、第１イメージセンサ１０４と第２イメージセンサ１０２を同期させて撮像する。第１イメージセンサ１０４では、被写体Ｗを撮像したデータが取得される。第２イメージセンサ１０２では点光源であるレーザ光源１０８を撮像したデータが取得される。第１イメージセンサ１０４で撮像されたデータは、被写体画像データとして制御装置１１２に記録される。第２イメージセンサ１０２で撮像されたデータは、点像画像データとして制御装置１１２に記録される。
【０１２４】
被写体画像データと点像画像データは、いずれもＲＧＢそれぞれの色成分の照度の分布を表現している。本実施例では、ＲＧＢそれぞれの色成分に関して、ステップＳ２０８とステップＳ２１０の処理を実施して、それぞれの色成分に関して復元された照度の分布を組み合わせることで、原画像を生成する。
【０１２５】
ステップＳ２０８では、点像画像データから、ＰＳＦの分布ｈ（ｘ、ｙ）を取得する。本実施例では、点像画像データの照度の分布を正規化して、ＰＳＦの分布ｈ（ｘ、ｙ）とする。
【０１２６】
ステップＳ２１０では、被写体画像データにおける照度の分布を劣化画像の分布ｇ（ｘ、ｙ）とし、劣化画像の分布ｇ（ｘ、ｙ）と、ステップＳ２０８で取得されたＰＳＦの分布ｈ（ｘ、ｙ）を用いて、画像復元処理を行う。
【０１２７】
本実施例では、ステップＳ２１０において、以下のようにして画像復元処理を行う。
まず劣化画像の分布ｇ（ｘ、ｙ）をフーリエ変換して、劣化画像のスペクトル分布Ｇ（ｓ、ｔ）を取得する。
【０１２８】
そして、ステップＳ２０８で取得されたＰＳＦの分布ｈ（ｘ、ｙ）をフーリエ変換して、光学的伝達関数（ＯＴＦ：Optical Transfer Function）の分布Ｈ（ｓ、ｔ）を取得する。
【０１２９】
その後に、劣化画像のスペクトル分布Ｇ（ｓ、ｔ）をＯＴＦの分布Ｈ（ｓ、ｔ）で除して、原画像のスペクトル分布Ｆ（ｓ、ｔ）を得る。
【０１３０】
最後に、原画像のスペクトル分布Ｆ（ｓ、ｔ）を逆フーリエ変換することで、原画像の分布ｆ（ｘ、ｙ）が得られる。
【０１３１】
上述のように原画像の分布ｆ（ｘ、ｙ）をＲＧＢの色成分それぞれについて取得し、それぞれの色成分の照度の分布を組み合わせることで、劣化画像から原画像を復元することができる。すなわち、撮像時にぶれてしまった被写体Ｗの画像から、ぶれのない被写体Ｗの画像を得ることができる。ぶれのない被写体Ｗの画像は、例えばディスプレイに表示してもよいし、画像データとしてハードディスクなどの記憶装置に記憶してもよいし、通信回線を経由して他のコンピュータに送信してもよい。
【０１３２】
図３は本実施例の撮像装置１００によって被写体Ｗのぶれのない画像を取得する様子を示している。本実施例の撮像装置１００によれば、第１イメージセンサ１０４で撮像された被写体画像３０２と、第２イメージセンサ１０２で撮像された点像画像３０４を用いて、ぶれのない被写体画像３０６を取得することができる。撮像台１１０とセンサ支持台１０６の間でわずかな振動があっても、ぶれのない被写体画像３０６を取得することができる。撮像台１１０やセンサ支持台１０６に防振機構を設けることなく、被写体Ｗのぶれのない画像を得ることが可能となる。
【０１３３】
（実施例１の変形例１）
本実施例の撮像装置１５０は、図１に示す実施例１の撮像装置１００と同様のハードウェア構成を備えているが、制御装置１１２が図２のステップＳ２１０において行う画像復元処理のみが異なっている。以下では実施例１との相違点のみについて詳細に説明する。
【０１３４】
図４は本実施例の撮像装置１５０において図２のステップＳ２１０で行う画像復元処理の詳細を示すフローチャートである。
【０１３５】
ステップＳ４０２では、劣化画像の分布ｇ（ｘ、ｙ）を特定する。本実施例の劣化画像の分布ｇ（ｘ、ｙ）は、光の照度の分布を位相特性も含めて記述する複素関数である。本実施例では、ｇ（ｘ、ｙ）の実部を被写体画像データにおける照度の分布と一致させ、ｇ（ｘ、ｙ）の虚部を全て０とする。
【０１３６】
次にステップＳ４０４では、ＰＳＦの分布ｈ（ｘ、ｙ）をフーリエ変換して、光学的伝達関数（ＯＴＦ：Optical Transfer Function）の分布Ｈ（ｓ、ｔ）を特定する。
【０１３７】
ステップＳ４０６では、ステップＳ４０４で特定された劣化画像の分布ｇ（ｘ、ｙ）を原画像の最初の推定分布ｆ０（ｘ、ｙ）に設定する。また反復計算の繰り返し数ｋを０に設定する。
【０１３８】
ステップＳ４０８では、原画像の推定分布ｆｋ（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ１（ｓ、ｔ）に設定する。関数Ｋ１（ｓ、ｔ）は、第１の関数に相当する。ｋは非負の整数であり、反復計算の回数に応じて後のステップＳ４２６で増加していく。
【０１３９】
ステップＳ４１０では、ステップＳ４０８で設定した第１の関数Ｋ１（ｓ、ｔ）に、ステップＳ４０４で特定したＨ（ｓ、ｔ）を乗じて、関数Ｋ２（ｓ、ｔ）を算出する。関数Ｋ２（ｓ、ｔ）は、第２の関数に相当する。
【０１４０】
ステップＳ４１２では、ステップＳ４１０で算出した第２の関数Ｋ２（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ３（ｘ、ｙ）に設定する。関数Ｌ３（ｘ、ｙ）は、複素関数であって、第３の関数に相当する。
【０１４１】
ステップＳ４１４では、ステップＳ４０２で特定した劣化画像分布ｇ（ｘ、ｙ）を、ステップＳ４１２で設定した第３の関数Ｌ３（ｘ、ｙ）で除して、関数Ｌ４（ｘ、ｙ）を算出する。関数Ｌ４（ｘ、ｙ）は、複素関数であって、第４の関数に相当する。
【０１４２】
ステップＳ４１６では、ステップＳ４１４で算出した第４の関数Ｌ４（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ５（ｓ、ｔ）に設定する。関数Ｋ５（ｓ、ｔ）は、第５の関数に相当する。
【０１４３】
ステップＳ４１８では、ステップＳ４１６で設定した第５の関数Ｋ５（ｓ、ｔ）に、Ｈ＃（ｓ、ｔ）を乗じて、関数Ｋ６（ｓ、ｔ）を算出する。Ｈ＃（ｓ、ｔ）は、ステップＳ４０４で特定したＨ（ｓ、ｔ）の反転関数であって、Ｈ＃（ｓ、ｔ）＝Ｈ（−ｓ、−ｔ）の関係を満たす。関数Ｋ６（ｓ、ｔ）は、第６の関数に相当する。
【０１４４】
ステップＳ４２０では、ステップＳ４１８で設定した第６の関数Ｋ６（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ７（ｘ、ｙ）に設定する。関数Ｌ７（ｘ、ｙ）は、複素関数であって、第７の関数に相当する。
【０１４５】
ステップＳ４２２では、原画像の推定分布ｆｋ（ｘ、ｙ）に、ステップＳ４２０で設定した第７の関数Ｌ７（ｘ、ｙ）を乗じて、原画像の改善された推定分布ｆｋ＋１（ｘ、ｙ）を算出する。前記第７の関数は一般に虚部を備える複素関数となるから、上記の方法によって原画像の推定分布を位相特性も含めて改善することができる。
【０１４６】
ステップＳ４２４では、ステップＳ４２２で算出された原画像の改善された推定分布ｆｋ＋１（ｘ、ｙ）と、原画像の推定分布ｆｋ（ｘ、ｙ）の差分を算出し、その絶対値が全ての（ｘ、ｙ）に対してしきい値εを下回るか否かを判断する。前記差分の絶対値がある（ｘ、ｙ）に対してしきい値ε以上の場合（ステップＳ４２４でＮＯの場合）、原画像の改善された推定分布ｆｋ＋１（ｘ、ｙ）はまだ収束に達していないと判断し、処理はステップＳ４２６へ進む。前記差分の絶対値が全ての（ｘ、ｙ）に対してしきい値εを下回る場合（ステップＳ４２４でＹＥＳの場合）、原画像の改善された推定分布ｆｋ＋１（ｘ、ｙ）は収束に達したものと判断して、処理はステップＳ４２８へ進む。
【０１４７】
ステップＳ４２６では、反復計算の繰り返し数ｋを１増加させる。処理はステップＳ４０８へ戻り、ステップＳ４０８からステップＳ４２４までの処理を再度実施する。
【０１４８】
ステップＳ４２８では、反復計算の結果得られた原画像の推定分布ｆｋ＋１（ｘ、ｙ）を原画像の分布ｆ（ｘ、ｙ）として出力する。
【０１４９】
以上の原画像の復元処理をＲＧＢの色成分それぞれについて行うことで、被写体Ｗのぶれのない画像を取得することができる。取得されたぶれのない画像は、例えばディスプレイに表示してもよいし、画像データとしてハードディスクなどの記憶装置に記憶してもよいし、通信回線を経由して他のコンピュータに送信してもよい。
【０１５０】
本実施例の撮像装置１５０によれば、劣化画像の分布ｇ（ｘ、ｙ）とＰＳＦの分布ｈ（ｘ、ｙ）を用いて、Ｂａｙｅｓの理論に基づき、劣化画像の分布ｇ（ｘ、ｙ）を実現する原画像の分布ｆ（ｘ、ｙ）のうちで最大に尤もらしい分布を、反復計算によって推定することができる。この撮像装置１５０は、ＯＴＦの分布Ｈ（ｓ、ｔ）がある周波数領域で０となる場合、すなわち実施例１の撮像装置１００による画像復元ができない場合であっても、原画像を復元することができる。
【０１５１】
（実施例１の変形例２）
本実施例の撮像装置１６０は、図１に示す実施例１の撮像装置１００や実施例１の変形例１の撮像装置１５０と同様のハードウェア構成を備えているが、制御装置１１２が図２のステップＳ２１０において行う画像復元処理の手法のみが異なっている。以下では実施例１および実施例１の変形例１との相違点のみについて詳細に説明する。
【０１５２】
図５は本実施例の撮像装置１６０において図２のステップＳ２１０で行う画像復元処理の詳細を示すフローチャートである。
【０１５３】
ステップＳ５０２では、劣化画像の分布ｇ（ｘ、ｙ）を特定する。本実施例の劣化画像の分布ｇ（ｘ、ｙ）は、光の照度の分布を位相特性も含めて記述する複素関数である。本実施例では、ｇ（ｘ、ｙ）の実部を被写体画像データにおける照度の分布と一致させ、ｇ（ｘ、ｙ）の虚部を全て０とする。
【０１５４】
ステップＳ５０４では、原画像の最初の推定分布ｆ０（ｘ、ｙ）を設定する。本実施例のでは、原画像の最初の推定分布ｆ０（ｘ、ｙ）として、劣化画像の分布ｇ（ｘ、ｙ）を用いる。一般的に、原画像の分布と劣化画像の分布は大きく異ならないと考えられるため、上記のように原画像の最初の推定分布ｆ０（ｘ、ｙ）を設定することによって、原画像の復元に伴う反復計算の回数を低減することができる。
またステップＳ５０４では、ＰＳＦの最初の推定分布ｈ０（ｘ、ｙ）を設定する。本実施例では、ＰＳＦの最初の推定分布ｈ０（ｘ、ｙ）として、図２のステップＳ２０８で取得されているＰＳＦの分布ｈ（ｘ、ｙ）を用いる。第２イメージセンサ１０２で撮像する際のＰＳＦの分布は、第１イメージセンサ１０４で撮像する際のＰＳＦの分布と大きく異ならないと考えられるため、上記のようにＰＳＦの最初の推定分布ｈ０（ｘ、ｙ）を設定することによって、原画像の復元に伴う反復計算の回数を低減することができる。
さらにステップＳ５０４では、反復計算の繰返し数ｋを０に設定する。
【０１５５】
ステップＳ５０６では、反復計算の繰返し数ｋを１増加させる。
【０１５６】
ステップＳ５０８では、ＰＳＦの推定分布ｈｋ（ｘ、ｙ）を算出して、ＰＳＦの推定分布を更新する。この計算については後述する。
【０１５７】
ステップＳ５１０では、ＰＳＦの推定分布ｈｋ（ｘ、ｙ）をフーリエ変換して、ＯＴＦの推定分布Ｈｋ（ｓ、ｔ）を算出する。
【０１５８】
ステップＳ５１２では、原画像の推定分布ｆｋ（ｘ、ｙ）を算出して、原画像の推定分布を更新する。この計算については後述する。
【０１５９】
ステップＳ５１４では、ステップＳ５１２で更新された原画像の推定分布ｆｋ（ｘ、ｙ）と、更新される前の原画像の推定分布ｆｋ−１との差分を算出し、その差分がすべての（ｘ、ｙ）に対してあるしきい値εを下回るか否かを判断する。前記差分がある（ｘ、ｙ）に対してしきい値ε以上の場合（ステップＳ５１４でＮＯの場合）、原画像の推定分布ｆｋ（ｘ、ｙ）は未だ収束に達していないと判断して、処理はステップＳ５０６へ戻り、ステップＳ５０６からステップＳ５１４までの処理を再度繰返す。前記差分がすべての（ｘ、ｙ）に対してしきい値εを下回る場合（ステップＳ５１４でＹＥＳの場合）、原画像の推定分布ｆｋ（ｘ、ｙ）は収束に達したと判断して、処理はステップＳ５１６へ進む。
【０１６０】
ステップＳ５１６では、反復計算の結果得られた原画像の推定分布ｆｋ（ｘ、ｙ）を、復元された原画像の分布ｆ（ｘ、ｙ）に設定する。本実施例の方法では、反復計算の過程で画像を記述する分布を複素関数として扱っているため、上記で設定される原画像の分布ｆ（ｘ、ｙ）は一般に複素関数となる。本実施例では、復元された原画像の分布ｆ（ｘ、ｙ）の実部が、原画像における照度の分布を表す。
【０１６１】
以下では図５のステップＳ５０８のｈｋ（ｘ、ｙ）の推定について詳細に説明する。本実施例では、ｈｋ（ｘ、ｙ）の推定として、図６のフローチャートに示すステップＳ６０２〜Ｓ６２６の工程を実施する。
【０１６２】
ステップＳ６０２では、ｈｋ（ｘ、ｙ）の最初の推定分布ｈｋ０（ｘ、ｙ）を設定する。本実施例の方法では、すでに取得されているｈｋ−１（ｘ、ｙ）を、ｈｋ（ｘ、ｙ）の最初の推定ｈｋ０（ｘ、ｙ）とする。また、ステップＳ６０２では、ｈｋ（ｘ、ｙ）の推定に係る反復計算の繰返し数ｍを０に設定する。
【０１６３】
ステップＳ６０４では、原画像の推定スペクトル分布Ｆ（ｓ、ｔ）を設定する。本実施例では、すでに取得されている原画像の推定分布ｆｋ−１（ｘ、ｙ）をフーリエ変換して、その結果をＦ（ｓ、ｔ）とする。
【０１６４】
ステップＳ６０６では、ｈｋ（ｘ、ｙ）の推定分布ｈｋｍ（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ１（ｓ、ｔ）に設定する。関数Ｋ１（ｓ、ｔ）は、第１の関数に相当する。
【０１６５】
ステップＳ６０８では、ステップＳ６０６で算出された関数Ｋ１（ｓ、ｔ）に、ステップＳ６０４で算出された原画像の推定スペクトル分布Ｆ（ｓ、ｔ）を乗じて、関数Ｋ２（ｓ、ｔ）を算出する。関数Ｋ２（ｓ、ｔ）は、第２の関数に相当する。
【０１６６】
ステップＳ６１０では、ステップＳ６０８で算出された関数Ｋ２（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ３（ｘ、ｙ）に設定する。関数Ｌ３（ｘ、ｙ）は、第３の関数に相当する。
【０１６７】
ステップＳ６１２では、図５のステップＳ５０４で特定された劣化画像分布ｇ（ｘ、ｙ）を、図６のステップＳ６１０で算出された関数Ｌ３（ｘ、ｙ）で除して、関数Ｌ４（ｘ、ｙ）を算出する。関数Ｌ４（ｘ、ｙ）は、第４の関数に相当する。
【０１６８】
ステップＳ６１４では、ステップＳ６１２で算出された関数Ｌ４（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ５（ｓ、ｔ）に設定する。関数Ｋ５（ｓ、ｔ）は、第５の関数に相当する。
【０１６９】
ステップＳ６１６では、ステップＳ６１４で算出された関数Ｋ５（ｓ、ｔ）に、Ｆ＃（ｓ、ｔ）を乗じて、関数Ｋ６（ｓ、ｔ）を算出する。関数Ｋ６（ｓ、ｔ）は、第６の関数に相当する。Ｆ＃（ｓ、ｔ）は、ステップＳ６０４で算出された原画像の推定スペクトル分布Ｆ（ｓ、ｔ）の反転関数であり、Ｆ＃（ｓ、ｔ）＝Ｆ（−ｓ、−ｔ）の関係を満たす。
【０１７０】
ステップＳ６１８では、ステップＳ６１６で算出された関数Ｋ６（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ７（ｘ、ｙ）に設定する。関数Ｌ７（ｘ、ｙ）は、第７の関数に相当する。
【０１７１】
ステップＳ６２０では、ｈｋ（ｘ、ｙ）の推定分布ｈｋｍ（ｘ、ｙ）に、ステップＳ６１８で設定された関数Ｌ７（ｘ、ｙ）を乗じて、その結果をｈｋ（ｘ、ｙ）の改善された推定分布ｈｋｍ＋１（ｘ、ｙ）に設定する。
【０１７２】
ステップＳ６２２では、ｈｋ（ｘ、ｙ）の推定のための反復計算の繰返し数ｍを１増加させる。
【０１７３】
ステップＳ６２４では、ｈｋ（ｘ、ｙ）の推定のための反復計算の繰返し数ｍが、しきい値以上になったか否かを判断する。本実施例では、ｈｋ（ｘ、ｙ）の推定のために実施する反復計算の回数を５回としている。繰り返し数ｍが５未満の場合（ステップＳ６２４でＮＯの場合）、処理はステップＳ６０６へ進み、ステップＳ６０６からステップＳ６２２までの処理を再度実施する。繰り返し数ｍが５以上の場合（ステップＳ６２４でＹＥＳの場合）、処理はステップＳ６２６へ進む。
【０１７４】
ステップＳ６２６では、上記の反復計算の結果得られる分布ｈｋ５（ｘ、ｙ）を、ｈｋ（ｘ、ｙ）の推定分布として設定する。
上記の反復計算は、劣化画像分布ｇ（ｘ、ｙ）と、原画像分布ｆｋ−１（ｘ、ｙ）から、ＰＳＦの分布ｈｋ（ｘ、ｙ）を推定することに相当する。原画像分布ｆｋ−１（ｘ、ｙ）が、真の原画像分布ｆ（ｘ、ｙ）に近いほど、上記の反復計算によって推定されるｈｋ（ｘ、ｙ）は、真のＰＳＦの分布ｈ（ｘ、ｙ）に近づく。
【０１７５】
以下では図５のステップＳ５１４のｆｋ（ｘ、ｙ）の推定について詳細に説明する。本実施例の方法では、ｆｋ（ｘ、ｙ）の推定として、図７のフローチャートに示すステップＳ７０２〜Ｓ７２４の工程を実施する。
【０１７６】
ステップＳ７０２では、ｆｋ（ｘ、ｙ）の最初の推定分布ｆｋ０（ｘ、ｙ）を設定する。本実施例の方法では、既に取得されているｆｋ−１（ｘ、ｙ）を、ｆｋ（ｘ、ｙ）の最初の推定分布ｆｋ０（ｘ、ｙ）として設定する。また、ステップＳ７０２では、ｆｋ（ｘ、ｙ）の推定のための反復計算の繰返し数ｎを０に設定する。
【０１７７】
ステップＳ７０４では、ｆｋ（ｘ、ｙ）の推定分布ｆｋｎ（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ１（ｓ、ｔ）に設定する。関数Ｋ１（ｓ、ｔ）は、第１の関数に相当する。
【０１７８】
ステップＳ７０６では、ステップＳ７０４で算出された関数Ｋ１（ｓ、ｔ）に、図５のステップＳ５１０で算出されたＯＴＦの推定分布Ｈｋ（ｓ、ｔ）を乗じて、関数Ｋ２（ｓ、ｔ）を算出する。関数Ｋ２（ｓ、ｔ）は、第２の関数に相当する。
【０１７９】
ステップＳ７０８では、ステップＳ７０６で算出された関数Ｋ２（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ３（ｘ、ｙ）に設定する。関数Ｌ３（ｘ、ｙ）は、第３の関数に相当する。
【０１８０】
ステップＳ７１０では、図５のステップＳ５０４で特定された劣化画像分布ｇ（ｘ、ｙ）を、図７のステップＳ７０８で算出された関数Ｌ３（ｘ、ｙ）で除して、関数Ｌ４（ｘ、ｙ）を算出する。関数Ｌ４（ｘ、ｙ）は、第４の関数に相当する。
【０１８１】
ステップＳ７１２では、ステップＳ７１０で算出された関数Ｌ４（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ５（ｓ、ｔ）に設定する。関数Ｋ５（ｓ、ｔ）は、第５の関数に相当する。
【０１８２】
ステップＳ７１４では、ステップＳ７１２で算出された関数Ｋ５（ｓ、ｔ）に、Ｈｋ＃（ｓ、ｔ）を乗じて、関数Ｋ６（ｓ、ｔ）を算出する。Ｈｋ＃（ｓ、ｔ）は、図５のステップＳ５１２で算出されたＯＴＦの推定分布Ｈｋ（ｓ、ｔ）の反転関数であり、Ｈｋ＃（ｓ、ｔ）＝Ｈｋ（−ｓ、−ｔ）の関係を満たす。
【０１８３】
ステップＳ７１６では、ステップＳ７１４で算出された関数Ｋ６（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ７（ｘ、ｙ）に設定する。関数Ｌ７（ｘ、ｙ）は、第７の関数に相当する。
【０１８４】
ステップＳ７１８では、ｆｋ（ｘ、ｙ）の推定分布ｆｋｎ（ｘ、ｙ）に、ステップＳ７１６で算出された関数Ｌ７（ｘ、ｙ）を乗じて、ｆｋ（ｘ、ｙ）の改善された推定分布ｆｋｎ＋１（ｘ、ｙ）を算出する。
【０１８５】
ステップＳ７２０では、ｆｋ（ｘ、ｙ）の推定のための反復計算の繰返し数ｎを１増加する。
【０１８６】
ステップＳ７２２では、ｆｋ（ｘ、ｙ）の推定のための反復計算の繰返し数ｎが、しきい値以上になったか否かを判断する。本実施例では、ｆｋ（ｘ、ｙ）の推定のために実施する反復計算の回数を５回としている。繰り返し数ｎが５未満の場合（ステップＳ７２２でＮＯの場合）、処理はステップＳ７０４へ進み、ステップＳ７０４からステップＳ７２０までの処理を再度実施する。繰り返し数ｎが５以上の場合（ステップＳ７２２でＹＥＳの場合）、処理はステップＳ７２４へ進む。
【０１８７】
ステップＳ７２４では、上記の反復計算の結果得られる分布ｆｋ５（ｘ、ｙ）を、ｆｋ（ｘ、ｙ）の推定分布として設定する。
上記の反復計算は、劣化画像分布ｇ（ｘ、ｙ）と、ＯＴＦ分布Ｈｋ（ｓ、ｔ）とに基づいて、原画像分布ｆｋ（ｘ、ｙ）を推定することに相当する。ＯＴＦ分布Ｈｋ（ｓ、ｔ）が、真のＯＴＦ分布Ｈ（ｓ、ｔ）に近いほど、上記の反復計算によって推定されるｆｋ（ｘ、ｙ）は、真の原画像の分布ｆ（ｘ、ｙ）に近づく。
【０１８８】
上記の実施例では、反復計算の過程において、先にＯＴＦの推定分布を更新して、次に原画像の推定分布を更新する例を説明した。ＯＴＦの推定分布の更新と、原画像の推定分布の更新は、交互に繰返し実施すればよいため、先に原画像の推定分布を更新して、次にＯＴＦの推定分布を更新してもよい。
【０１８９】
以上の原画像の復元処理をＲＧＢの色成分それぞれについて行うことで、被写体Ｗのぶれのない画像を取得することができる。取得されたぶれのない画像は、例えばディスプレイに表示してもよいし、画像データとしてハードディスクなどの記憶装置に記憶してもよいし、通信回線を経由して他のコンピュータに送信してもよい。
【０１９０】
本実施例の撮像装置１６０によれば、劣化画像の分布ｇ（ｘ、ｙ）とＰＳＦの分布ｈ（ｘ、ｙ）を用いて、Ｂａｙｅｓの理論に基づき、劣化画像の分布ｇ（ｘ、ｙ）を実現する原画像の分布ｆ（ｘ、ｙ）のうちで最大に尤もらしい分布を、反復計算によって推定することができる。この撮像装置１６０は、ＯＴＦの分布Ｈ（ｓ、ｔ）がある周波数領域で０となる場合、すなわち実施例１の撮像装置１００による画像復元ができない場合であっても、原画像を復元することができる。
【０１９１】
本実施例の撮像装置１６０によれば、第２イメージセンサ１０２で撮像された点像画像データに基いて取得されたＰＳＦの分布ｈ（ｘ、ｙ）が、第１イメージセンサ１０４で被写体Ｗを撮像する際の真のＰＳＦの分布と相違する場合であっても、劣化画像の分布ｇ（ｘ、ｙ）を実現するＰＳＦの分布ｈ（ｘ、ｙ）のうちで最大に尤もらしい分布を、反復計算によって推定することができる。実施例１の変形例１の撮像装置１５０にくらべ、より正確に原画像を復元することができる。
【０１９２】
（実施例２）
図面を参照しながら、本実施例の天体撮像装置８００について説明する。
図８は本実施例の天体撮像装置８００の構成を示す概略図である。本実施例の天体撮像装置８００は、天体を撮像する第１カメラ８０２と、電離層に向けてレーザ光を照射するレーザ照射器８０６と、レーザ照射器８０６から照射されて電離層で反射したレーザ光を撮像する第２カメラ８０４と、制御装置８０８を備えている。
【０１９３】
第１カメラ８０２と第２カメラ８０４とレーザ照射器８０６は、同一の地域内で、地上に載置されている。
【０１９４】
第１カメラ８０２と第２カメラ８０４は、それぞれ撮像面に複数の撮像素子を備えており、撮像面に形成される像を撮像素子を介して認識し、ＲＧＢの各色成分の照度の分布を示す電気信号に変換して制御装置８０８へ送信する。第１カメラ８０２と第２カメラ８０４は、レンズ収差等の特性が同一である。また、第１カメラ８０２と第２カメラ８０４は、それぞれのシャッタが同期して開閉するように制御装置８０８によって制御される。
【０１９５】
レーザ照射器８０６は、制御装置８０８によって制御されており、電離層に向けてレーザ光を照射する。レーザ照射器８０６で照射されたレーザ光は、電離層で反射して、第２カメラ８０４で受光される。
【０１９６】
第１カメラ８０２が撮像する天体の画像は、大気のゆれによるシーイングの影響を受けて、ぼやけた画像となってしまう。本実施例の天体撮像装置８００では、第１カメラ８０２で撮像された天体の画像と、第２カメラ８０４で撮像されたレーザ光の画像を用いて、シーイングの影響を除去した天体の画像を取得する。第２カメラ８０４で撮像される画像は、レーザ照射器８０６から電離層へ到達するまでの間にシーイングの影響を受け、さらに電離層で反射してから第２カメラ８０４へ到達するまでの間にシーイングの影響を受けた点光源の画像として扱うことができる。従って、第２カメラ８０４で撮像される画像に基づいて、第２カメラ８０４で天体を撮像する際のＰＳＦの分布を取得することができる。第２カメラ８０４と第１カメラ８０２は同一の特性を有しており、それぞれのシャッタが同期して開閉するように撮像するので、第２カメラ８０４で撮像する際のＰＳＦの分布は、第１カメラ８０２で撮像する際のＰＳＦの分布として扱うことができる。
【０１９７】
図９は本実施例の天体撮像装置８００を用いて惑星を撮像した様子を示している。図９に示すように、本実施例の天体撮像装置８００によれば、大気のゆれによるシーイングの影響によってぼやけてしまっている天体画像９０２から、シーイングの影響を排除した鮮明な天体画像９０４を取得することができる。
【０１９８】
図１０のフローチャートを参照しながら、天体撮像装置８００の動作について説明する。
【０１９９】
ステップＳ１００２では、レーザ照射器８０６によるレーザ光の照射を開始する。レーザ照射器８０６によって照射されたレーザ光は、電離層で反射して、第２カメラ８０４まで到達する。
【０２００】
ステップＳ１００４では、第１カメラ８０２と第２カメラ８０４を同期させて撮像する。第１カメラ８０２では、天体を撮像したデータが取得される。第２カメラ８０４では、レーザ照射器８０６から照射されて電離層によって反射されたレーザ光を撮像したデータが取得される。第１カメラで撮像されたデータは天体画像データとして制御装置８０８に記録される。第２カメラ８０４で撮像されたデータは点像画像データとして制御装置８０８に記録される。
【０２０１】
天体画像データと点像画像データは、いずれもＲＧＢそれぞれの色成分の照度の分布を表現している。本実施例では、ＲＧＢそれぞれの色成分に関して、ステップＳ１００６とステップＳ１００８の処理を実施する。
【０２０２】
ステップＳ１００６では、点像画像データから、ＰＳＦの分布ｈ（ｘ、ｙ）を取得する。本実施例の場合、第２カメラ８０４で撮像される画像は、地上と電離層の間を往復する間にシーイングの影響を受けたレーザ光の画像である。従って、天体を撮像する際のＰＳＦの分布としては、点像画像データから得られるＰＳＦの分布を用いる。
【０２０３】
ステップＳ１００８では、天体画像データにおける照度の分布と、ステップＳ１００６で取得されたＰＳＦの分布ｈ（ｘ、ｙ）を用いて、画像復元処理を行う。ステップＳ１００６で、すでに天体を撮像する際のＰＳＦの分布が取得されているので、実施例１や、実施例１の変形例１や、実施例１の変形例２で詳述した画像復元処理を施すことによって、シーイングの影響を排除した天体画像を得ることができる。
【０２０４】
本実施例の天体撮像装置８００によれば、大気のゆれによるシーイングの影響を排除した鮮明な天体画像を取得することができる。
【０２０５】
（実施例３）
図１１、図１２および図１３を参照しながら、実施例３の画像復元装置１１００について説明する。本実施例の画像復元装置１１００は、撮像時にぶれてしまった画像のみに基づいて、ぶれの影響を表現するＰＳＦの分布を推定して、ぶれのない理想的な画像を復元する。以下では、撮像時にぶれてしまった画像を劣化画像と呼び、ぶれのない理想的な画像を原画像と呼ぶ。
【０２０６】
図１１は画像復元装置１１００のハードウェア構成を示している。画像復元装置１１００は、後述する処理を実施するプログラムがインストールされた汎用のコンピュータである。画像復元装置１１００は、制御部１１０２と、入力部１１０４と、出力部１１０６と、演算部１１０８と、記憶部１１１０を備えている。記憶部１１１０は、プログラム記憶部１１１２と、データ記憶部１１１４を備えている。
画像復元装置１１００では、入力部１１０４を介して劣化画像のデータが入力されて、データ記憶部１１１４に記憶される。その後、制御部１１０２はプログラム記憶部１１１２に記憶されたプログラムに従って演算部１１０８で処理を行い、原画像を復元して原画像のデータをデータ記憶部１１１４に記憶する。データ記憶部１１１４に記憶された原画像のデータは、出力部１１０６を介して出力される。
なお画像復元装置１１００は、上記のような汎用のコンピュータの代わりに、同様の処理を実施可能な専用ロジックやＩＣ、ＦＰＧＡなどであってもよい。
【０２０７】
図１２のフローチャートを参照しながら、本実施例の画像復元装置１１００が行う処理の詳細を説明する。なお本実施例の画像復元装置１１００では、データ記憶部１１１４に記憶された劣化画像データから、ＲＧＢそれぞれの色成分の照度の分布ｇｒ（ｘ、ｙ）、ｇｇ（ｘ、ｙ）、ｇｂ（ｘ、ｙ）を特定しておいて、ＲＧＢそれぞれの色成分について図１２の処理を実行する。
【０２０８】
ステップＳ１２０２では、劣化画像の分布ｇ（ｘ、ｙ）を特定する。劣化画像の分布ｇ（ｘ、ｙ）は、光の照度の分布を位相特性も含めて記述する複素関数である。本実施例では、ｇ（ｘ、ｙ）の実部を劣化画像データにおける照度の分布と一致させ、ｇ（ｘ、ｙ）の虚部を全て０とする。
【０２０９】
ステップＳ１２０４では、劣化画像の分布ｇ（ｘ、ｙ）をフーリエ変換して、劣化画像のスペクトル分布Ｇ（ｓ、ｔ）を算出する。
【０２１０】
ステップＳ１２０６では、劣化画像のスペクトル分布Ｇ（ｓ、ｔ）の複素対数関数であるｌｎ（Ｇ（ｓ、ｔ））を計算して、関数Ｇｃ（ｓ、ｔ）に設定する。
【０２１１】
ステップＳ１２０８では、関数Ｇｃ（ｓ、ｔ）を逆フーリエ変換して、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）を算出する。
【０２１２】
図１３は図１２のステップＳ１２０８で算出される劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）の概略形状を示している。図１３に示すように、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）は、原画像の複素ケプストラムｆｃ（ｘ、ｙ）に、ぶれの影響を表現するＰＳＦの複素ケプストラムであるパルス列が重ね合わさった形状となっている。
【０２１３】
ステップＳ１２１０では、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）からパルス列を抽出する。本実施例では、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）の振幅値が所定のしきい値を超える個所をパルスとして識別する。
【０２１４】
ステップＳ１２１２では、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）から抽出されたパルス列を劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）から減じて、原画像の推定複素ケプストラムｆｃ（ｘ、ｙ）を算出する。
【０２１５】
ステップＳ１２１４では、原画像の推定複素ケプストラムｆｃ（ｘ、ｙ）をフーリエ変換して、関数Ｆｃ（ｓ、ｔ）を算出する。
【０２１６】
ステップＳ１２１６では、関数Ｆｃ（ｓ、ｔ）の複素指数関数であるｅｘｐ（Ｆｃ（ｓ、ｔ））を計算して、原画像の推定スペクトル分布Ｆ（ｓ、ｔ）を算出する。
【０２１７】
ステップＳ１２１８では、原画像の推定スペクトル分布Ｆ（ｓ、ｔ）を逆フーリエ変換して、原画像の推定分布ｆ（ｘ、ｙ）を算出し、原画像の推定分布ｆ（ｘ、ｙ）を原画像の分布として出力する。
【０２１８】
以上の処理をＲＧＢの色成分それぞれについて行うことで、ぶれのない画像を取得することができる。ぶれのない理想的な画像は、例えばディスプレイに表示してもよいし、画像データとしてハードディスクなどの記憶装置に記憶してもよいし、通信回線を経由して他のコンピュータに送信してもよい。
【０２１９】
なお図１４に本実施例のぶれ画像復元装置１１００が実現する機能のブロック図１４００を示している。データ入力ブロック１４０２から入力された劣化画像の分布ｇ（ｘ、ｙ）について、フーリエ変換ブロック１４０４の処理によって劣化画像のスペクトル分布Ｇ（ｓ、ｔ）が計算される。複素対数ブロック１４０６の処理によって関数Ｇｃ（ｓ、ｔ）が計算される。逆フーリエ変換ブロック１４０８の処理によって劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）が計算される。パルス抽出ブロック１４１０の処理によって劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）からパルス列が抽出される。減算ブロック１４１２によって劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）からパルス列が除去されて、原画像の推定複素ケプストラムｆｃ（ｘ、ｙ）が計算される。フーリエ変換ブロック１４１４の処理によって関数Ｆｃ（ｓ、ｔ）が計算される。複素指数ブロック１４１６の処理によって原画像の推定スペクトル分布Ｆ（ｓ、ｔ）が計算される。逆フーリエ変換ブロック１４１８の処理によって原画像の推定分布ｆ（ｘ、ｙ）が計算されて、データ出力ブロック１４２０に出力される。
【０２２０】
本実施例の画像復元装置１１００によれば、撮像時にぶれてしまった画像のみに基づいて、ぶれのない理想的な画像を復元することができる。
【０２２１】
（実施例４）
図面を参照しながら、本実施例の画像復元装置１１５０について説明する。本実施例の画像復元装置１１５０は、図１１に示す実施例３の画像復元装置１１００と同様のハードウェア構成を備えているが、プログラム記憶部１１１２にインストールされているプログラムの内容が異なる。
【０２２２】
図１５のフローチャートを参照しながら、本実施例の画像復元装置１１５０が行う処理の詳細を説明する。なお本実施例の画像復元装置１１５０では、データ記憶部１１１４に記憶された劣化画像データから、ＲＧＢそれぞれの色成分の照度の分布ｇｒ（ｘ、ｙ）、ｇｇ（ｘ、ｙ）、ｇｂ（ｘ、ｙ）を特定しておいて、ＲＧＢそれぞれの色成分について図１５の処理を実行する。
【０２２３】
ステップＳ１５０２では、劣化画像の分布ｇ（ｘ、ｙ）を特定する。劣化画像の分布ｇ（ｘ、ｙ）は、光の照度の分布を位相特性も含めて記述する複素関数である。本実施例では、ｇ（ｘ、ｙ）の実部を劣化画像データにおける照度の分布と一致させ、ｇ（ｘ、ｙ）の虚部を全て０とする。
【０２２４】
ステップＳ１５０４では、劣化画像の分布ｇ（ｘ、ｙ）をフーリエ変換して、劣化画像のスペクトル分布Ｇ（ｓ、ｔ）を算出する。
【０２２５】
ステップＳ１５０６では、劣化画像のスペクトル分布Ｇ（ｓ、ｔ）の複素対数関数であるｌｎ（Ｇ（ｓ、ｔ））を計算して、関数Ｇｃ（ｓ、ｔ）に設定する。
【０２２６】
ステップＳ１５０８では、関数Ｇｃ（ｓ、ｔ）を逆フーリエ変換して、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）を算出する。
【０２２７】
ステップＳ１５１０では、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）からパルス列を抽出する。本実施例では、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）の振幅値が所定のしきい値を超える個所をパルスとして識別する。
【０２２８】
ステップＳ１５１２では、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）から抽出されたパルス列より、ＰＳＦの分布ｈ（ｘ、ｙ）を特定する。画像の劣化がぶれによるものの場合には、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）に現れるパルス列から、ＰＳＦの分布ｈ（ｘ、ｙ）を直接的に求めることができる。
【０２２９】
ステップＳ１５１４では、ステップＳ１５０２で特定された劣化画像の分布ｇ（ｘ、ｙ）と、ステップＳ１５１２で特定されたＰＳＦの分布ｈ（ｘ、ｙ）を用いて、原画像の分布ｆ（ｘ、ｙ）を復元する。この画像の復元処理は、種々の手法によって行うことができる。
【０２３０】
例えば実施例１で詳述したように、ＰＳＦ分布ｈ（ｘ、ｙ）をフーリエ変換してＯＴＦ分布Ｈ（ｓ、ｔ）を取得し、劣化画像のスペクトル分布Ｇ（ｓ、ｔ）をＯＴＦ分布Ｈ（ｓ、ｔ）で除して原画像のスペクトル分布Ｆ（ｓ、ｔ）を算出し、原画像のスペクトル分布Ｆ（ｓ、ｔ）を逆フーリエ変換することによって、原画像の分布ｆ（ｘ、ｙ）を取得することができる。
【０２３１】
あるいは実施例１の変形例１で詳述したように、ＰＳＦ分布ｈ（ｘ、ｙ）と劣化画像の分布ｇ（ｘ、ｙ）を用いて、図４に示す画像復元処理を行って、原画像の分布ｆ（ｘ、ｙ）を取得することもできる。
【０２３２】
あるいは実施例１の変形例２のように、ＰＳＦ分布ｈ（ｘ、ｙ）を真のＰＳＦ分布の最初の推定分布ｈ０（ｘ、ｙ）として、劣化画像の分布ｇ（ｘ、ｙ）を用いて、図５、図６および図７に示す画像復元処理を行って、原画像の分布ｆ（ｘ、ｙ）を取得することもできる。
【０２３３】
なお図１６は本実施例のぶれ画像復元装置１１５０が実現する機能のブロック図１６００を示している。データ入力ブロック１６０２から入力された劣化画像の分布ｇ（ｘ、ｙ）について、フーリエ変換ブロック１６０４の処理によって劣化画像のスペクトル分布Ｇ（ｓ、ｔ）が計算される。複素対数ブロック１６０６の処理によって関数Ｇｃ（ｘ、ｙ）が計算される。逆フーリエ変換ブロック１６０８の処理によって劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）が計算される。パルス抽出ブロック１６１０の処理によって劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）からパルス列が抽出される。ＰＳＦ特定ブロック１６１２の処理によってＰＳＦ分布ｈ（ｘ、ｙ）が計算される。画像復元ブロック１６１４の処理によって原画像の分布ｆ（ｘ、ｙ）が計算されて、データ出力ブロック１６１６に出力される。
【０２３４】
本実施例の画像復元装置１１５０によれば、撮像時にぶれてしまった画像のみに基づいて、ぶれのない理想的な画像を復元することができる。
【０２３５】
以上、本発明の実施形態について詳細に説明したが、これらは例示に過ぎず、特許請求の範囲を限定するものではない。特許請求の範囲に記載の技術には、以上に例示した具体例を様々に変形、変更したものが含まれる。例えば、劣化復元の基本式である数式（１８）の代わりに、数式（１７）を反復計算に直接用いることも可能である。この場合にはフーリエ変換という過程は不要である。
また、本明細書または図面に説明した技術要素は、単独であるいは各種の組み合わせによって技術的有用性を発揮するものであり、出願時請求項記載の組み合わせに限定されるものではない。また、本明細書または図面に例示した技術は複数目的を同時に達成するものであり、そのうちの一つの目的を達成すること自体で技術的有用性を持つものである。
【符号の説明】
【０２３６】
１００：撮像装置
１０２：第２イメージセンサ
１０４：第１イメージセンサ
１０６：センサ支持台
１０８：レーザ光源
１１０：撮像台
１１２：制御装置
１５０：撮像装置
１６０：撮像装置
３０２：被写体画像
３０４：点像画像
３０６：ぶれのない被写体画像
８００：天体撮像装置
８０２：第１カメラ
８０４：第２カメラ
８０６：レーザ照射器
８０８：制御装置
９０２：天体画像
９０４：シーイングの影響を排除した天体画像
１１００：画像復元装置
１１０２：制御部
１１０４：入力部
１１０６：出力部
１１０８：演算部
１１１０：記憶部
１１１２：プログラム記憶部
１１１４：データ記憶部
１１５０：ぶれ画像復元装置
１４００：ブロック図
１４０２：データ入力ブロック
１４０４：フーリエ変換ブロック
１４０６：複素対数ブロック
１４０８：逆フーリエ変換ブロック
１４１０：パルス抽出ブロック
１４１２：減算ブロック
１４１４：フーリエ変換ブロック
１４１６：複素指数ブロック
１４１８：逆フーリエ変換ブロック
１４２０：データ出力ブロック
１６００：ブロック図
１６０２：データ入力ブロック
１６０４：フーリエ変換ブロック
１６０６：複素対数ブロック
１６０８：逆フーリエ変換ブロック
１６１０：パルス抽出ブロック
１６１２：ＰＳＦ特定ブロック
１６１４：画像復元ブロック
１６１６：データ出力ブロック
【技術分野】
【０００１】
本発明は撮像時のぶれやぼけの影響によって劣化した画像から、劣化のない画像を復元する技術に関する。
【背景技術】
【０００２】
被写体を撮像する際に、解像度が不足したり、画像がぶれてしまったり、ぼけてしまったりすることがある。このように劣化した画像から劣化のない理想的な画像を復元する技術の開発が進められている。
【０００３】
本発明者らは、これまでにも、劣化した画像から劣化のない画像を復元する技術を開発している。この技術は、例えば特許文献１や特許文献２に記載されている。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】国際公開第２００６／０４１１２６号
【特許文献２】国際公開第２００６／０４１１２７号
【発明の概要】
【発明が解決しようとする課題】
【０００５】
劣化した画像から劣化のない画像を復元するためには、撮像時のぶれやぼけの影響を表現する点像分布関数が必要とされる。適切な点像分布関数を取得することができれば、劣化した画像から劣化のない画像を正確かつ速やかに復元することが可能となる。
【０００６】
本発明は上記課題を解決する。本発明では、撮像時のぶれやぼけの影響を表現する点像分布関数を適切に取得して、劣化した画像から劣化のない画像を復元することが可能な技術を提供する。
【課題を解決するための手段】
【０００７】
本発明は撮像装置として具現化される。その撮像装置は、被写体を撮像して被写体画像データを出力する第１撮像手段と、点光源を撮像して点像画像データを出力する第２撮像手段と、被写体画像データと点像画像データに基いて改善された被写体画像データを生成する画像処理手段を備えている。その撮像装置において、第１撮像手段と第２撮像手段は、共通する支持台によって支持されており、互いに同期して撮像するものである。その撮像装置において、被写体と点光源は、共通する載置台上に載置されている。その撮像装置において、画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と被写体画像データに基いて改善された被写体画像データを生成する。
【０００８】
上記の撮像装置の動作内容を方法の発明として具現化することもできる。本発明の撮像方法は、第１撮像手段によって被写体を撮像して被写体画像データを出力する工程と、第２撮像手段によって点光源を撮像して点像画像データを出力する工程と、画像処理手段によって被写体画像データと点像画像データに基いて改善された被写体画像データを生成する工程を備えている。その撮像方法において、第１撮像手段と第２撮像手段は、共通する支持台によって支持されており、互いに同期して撮像するものである。その撮像方法において、被写体と点光源は、共通する載置台上に載置されている。その撮像方法において、画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と被写体画像データに基いて改善された被写体画像データを生成する。
なお本発明は上記の撮像方法における各工程をコンピュータに実行させるためのプログラムとしても具現化することができる。
【０００９】
第１撮像手段で撮像される被写体の画像は、ぶれやぼけによって劣化していることがある。上記の撮像装置では、第１撮像手段と同様の条件で撮像する第２撮像手段によって点光源を撮像して、第２撮像手段で撮像された画像から点像分布関数（ＰＳＦ：Point Spread Function）を取得する。そして、第１撮像手段で撮像された画像とＰＳＦを用いて、ぶれやぼけのない改善された被写体の画像を生成する。
【００１０】
以下では撮像時のぶれやぼけによって劣化している画像を劣化画像と呼び、ぶれやぼけのない理想的な画像を原画像と呼ぶ。劣化画像とＰＳＦから原画像を復元する原理について説明する。原画像と劣化画像はそれぞれの大きさが同一であり、いずれも画像内の点の位置を座標（ｘ、ｙ）で表現することが可能であって、原画像の照度分布はｆｒ（ｘ、ｙ）、劣化画像の照度分布はｇｒ（ｘ、ｙ）で表現されるものとする。またＰＳＦの分布をｈｒ（ｘ、ｙ）とする。
【００１１】
原画像の分布ｆｒ（ｘ、ｙ）と、劣化画像の分布ｇｒ（ｘ、ｙ）と、ＰＳＦの分布ｈｒ（ｘ、ｙ）の間には、次の関係が成り立つ。
【００１２】
【数１】
【００１３】
原画像の分布ｆｒ（ｘ、ｙ）と劣化画像の分布ｇｒ（ｘ、ｙ）は、以下に示す二次元フーリエ変換を実施することによって、ｘ軸に関する空間周波数ｓと、ｙ軸に関する空間周波数ｔについてのスペクトルをそれぞれ得ることができる。
【００１４】
【数２】
【数３】
【００１５】
ここでＦＴ（）は二次元フーリエ変換を示している。またＰＳＦの分布ｈｒ（ｘ、ｙ）について二次元フーリエ変換を実施することによって、周波数空間における伝達関数である光学的伝達関数（ＯＴＦ：Optical Transfer Function）の分布Ｈｒ（ｓ、ｔ）を得ることができる。
【００１６】
【数４】
【００１７】
式（１）の両辺についてフーリエ変換することで、以下に示す劣化画像のスペクトル分布Ｇｒ（ｓ、ｔ）と、原画像のスペクトル分布Ｆｒ（ｓ、ｔ）と、ＯＴＦ分布Ｈｒ（ｓ、ｔ）の関係が得られる。
【００１８】
【数５】
【００１９】
従って、劣化画像のスペクトル分布Ｇｒ（ｓ、ｔ）とＯＴＦ分布Ｈｒ（ｓ、ｔ）が既知であれば、次式により原画像のスペクトル分布Ｆｒ（ｓ、ｔ）を得ることができる。
【００２０】
【数６】
【００２１】
上式によって原画像のスペクトル分布Ｆｒ（ｓ、ｔ）が得られれば、二次元の逆フーリエ変換を施すことによって、原画像の分布ｆｒ（ｘ、ｙ）を得ることができる。
【００２２】
【数７】
【００２３】
ここでＦＴ−１（）は二次元の逆フーリエ変換を示している。以上のように、本発明の画像取得装置によれば、第２撮像手段で撮像された画像からＰＳＦの分布を取得し、第１撮像手段で撮像された画像とＰＳＦの分布を用いて、ぶれやぼけのない被写体の画像を得ることができる。
【００２４】
上記の撮像装置においては、前記画像処理手段が、点像分布関数と被写体画像データに基づいて改善された被写体画像データを生成する際に、被写体画像データから劣化画像の分布を特定する処理と、点像分布関数から光学的伝達関数を特定する処理と、原画像の最初の推定分布を特定する処理と、（１）原画像の推定分布をフーリエ変換して第１の関数を得る処理と、（２）前記第１の関数に前記光学的伝達関数を乗じて第２の関数を得る処理と、（３）前記第２の関数を逆フーリエ変換して第３の関数を得る処理と、（４）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る処理と、（５）前記第４の関数をフーリエ変換して第５の関数を得る処理と、（６）前記第５の関数に前記光学的伝達関数の反転関数を乗じて第６の関数を得る処理と、（７）前記第６の関数を逆フーリエ変換して第７の関数を得る処理と、（８）原画像の推定分布に前記第７の関数を乗じて、原画像の次の推定分布を得る処理と、原画像の次の推定分布を原画像の推定分布に置き換えて、前記（１）から（８）の処理を繰返す処理と、原画像の推定分布に基づいて、改善された被写体画像データを生成する処理を実行することが好ましい。
【００２５】
劣化画像から原画像を復元する際に、ＯＴＦ分布Ｈ（ｓ、ｔ）が０となる周波数領域が存在する場合には、数式（６）の演算を行うことができず、原画像を完全に復元することができない。そこで上記の撮像装置では、原画像の分布ｆｒ（ｘ、ｙ）と劣化画像の分布ｇｒ（ｘ、ｙ）を確率密度関数として扱い、Ｂａｙｅｓの理論に基づいて原画像の推定を行う。原画像分布ｆｒ（ｘ、ｙ）と劣化画像分布ｇｒ（ｘ、ｙ）は、下記の正規化を行うことによって、確率密度関数として扱うことが可能になる。
【００２６】
【数８】
【数９】
【００２７】
上記の正規化に合わせて、ＯＴＦ分布Ｈｒ（ｓ、ｔ）についても正規化する。Ｈｒ（ｓ、ｔ）は、空間周波数の原点での値を用いて正規化する。
【００２８】
【数１０】
【００２９】
正規化された原画像の分布ｆ（ｘ、ｙ）と、劣化画像の分布ｇ（ｘ、ｙ）は、非負の関数であり、定義された領域での積分値が１であるから、確率密度関数として扱うことができる。上記の場合において、ｆ（ｘ、ｙ）は原画像の座標（ｘ、ｙ）における結像という事象の確率密度関数である。ｇ（ｘ、ｙ）は劣化画像の座標（ｘ、ｙ）における結像という事象の確率密度関数である。
【００３０】
原画像および劣化画像の分布を、確率密度関数と見なすことが可能な場合、Ｂａｙｅｓの理論に基づいて、劣化画像の分布から、その劣化画像を生じさせている原画像の分布を推定することができる。
原画像の座標（ｘ、ｙ）に光が結像する事象は、その座標（ｘ、ｙ）に点光源が存在する事象として扱うことができる。原画像の座標（ｘ１、ｙ１）において点光源が存在する事象をＶ（ｘ１、ｙ１）、劣化画像の座標（ｘ２、ｙ２）において点像が結像する事象をＡ（ｘ２、ｙ２）とした場合、それぞれの事象の確率Ｐ（Ｖ）およびＰ（Ａ）は以下で表現される。
【００３１】
【数１１】
【数１２】
【００３２】
また、原画像の座標（ｘ１、ｙ１）に点光源が存在している場合に、劣化画像の座標（ｘ２、ｙ２）に結像する確率は、事象Ｖ（ｘ１、ｙ１）の生起を条件とする事象Ａ（ｘ２、ｙ２）の生起確率である。上記の確率はＰＳＦ分布ｈ（ｘ、ｙ）を用いて以下で表現される。
【００３３】
【数１３】
【００３４】
Ｂａｙｅｓの理論に基づくと、劣化画像内の点（ｘ２、ｙ２）に点像を結像させる場合の原画像の分布Ｐ（Ｖ（ｘ、ｙ）｜Ａ（ｘ２、ｙ２））は、以下で推定される。
【００３５】
【数１４】
【００３６】
上式の右辺のＰ（Ｖ（ｘ、ｙ））、Ｐ（Ａ（ｘ２、ｙ２）｜Ｖ（ｘ、ｙ））に、数式（１１）および数式（１３）を代入すると、次式を得る。
【００３７】
【数１５】
【００３８】
上式の左辺は、劣化画像において点像が結像している場合に、推定される原画像の分布を表している。上式に劣化画像の分布ｇ（ｘ、ｙ）を乗じて積分することによって、劣化画像の分布ｇ（ｘ、ｙ）を実現するための原画像の分布ｆ（ｘ、ｙ）を得ることができる。
上式の両辺に、Ｐ（Ａ（ｘ２、ｙ２））＝ｇ（ｘ２、ｙ２）を乗じて、すべての（ｘ２、ｙ２）に関して積分すると、次式が得られる。
【００３９】
【数１６】
【００４０】
上式の左辺に数式（１２）を代入すると、その積分の結果はＰ（Ｖ（ｘ、ｙ））であり、ｆ（ｘ、ｙ）に等しい。従って、Ｂａｙｅｓの理論に基づくと、以下の関係が成り立つ。
【００４１】
【数１７】
【００４２】
上記の関係は、分布ｆ（ｘ、ｙ）が真の原画像の分布である場合に成立すると考えられる。すなわち、上式を満たす分布ｆ（ｘ、ｙ）を算出することが、劣化画像の復元に相当する。
上記の関係を満たす分布ｆ（ｘ、ｙ）は、数式（１７）の右辺のｆ（ｘ、ｙ）をｆｋ（ｘ、ｙ）とし、数式（１７）の左辺のｆ（ｘ、ｙ）をｆｋ＋１（ｘ、ｙ）として、ｆｋ（ｘ、ｙ）に関する反復計算を実施し、ｆｋ（ｘ、ｙ）の収束値を求めることで算出することができる。上記の反復計算によって求まるｆｋ（ｘ、ｙ）の収束値は、Ｂａｙｅｓの理論に基づく原画像の推定分布に相当する。
【００４３】
上記では、原画像の分布ｆ（ｘ、ｙ）や劣化画像の分布ｇ（ｘ、ｙ）を正規化した場合について説明したが、実際の反復計算を実施する上では、これらの分布は正規化することなく、そのまま使用することができる。
【００４４】
上記の反復計算においては、反復計算を実施する前に、原画像の最初の推定分布ｆ０（ｘ、ｙ）を設定しておく。最初の推定分布ｆ０（ｘ、ｙ）としては、任意の分布を設定することができる。一般的には、劣化画像の分布ｇ（ｘ、ｙ）は原画像の分布ｆ（ｘ、ｙ）から大きく異なることはないため、最初の推定分布ｆ０（ｘ、ｙ）としては、劣化画像の分布ｇ（ｘ、ｙ）を用いることが好ましい。
【００４５】
数式（１７）の右辺は、ＰＳＦ分布ｈ（ｘ、ｙ）を用いた畳み込み積分を備えている。一般に、光学系の位相特性まで含めてＰＳＦを正確に評価することは困難であり、上記の反復計算を正確な位相特性を含めて実施することは困難である。正確な位相特性を含まないＰＳＦを用いた反復計算は、誤った収束の結果をもたらすため、原画像の正確な復元の妨げとなる。
そこで、ＰＳＦの代わりに正確な位相特性を含ませることが容易であるＯＴＦを用いることで、より正確に原画像を復元することが可能となる。また、復元の過程で位相特性を正確に評価するために、原画像の推定分布ｆｋ（ｘ、ｙ）（ｋ＝０，１，２，・・・）を複素関数に拡張し、その実部が画像分布を表現するものとして取扱う。上式の右辺に、フーリエ変換と逆フーリエ変換を用いることで、ＯＴＦを用いた形式に変更することができる。フーリエ変換をＦＴ（）、逆フーリエ変換をＦＴ−１（）で表現すると、数式（１７）は以下で表現される。
【００４６】
【数１８】
【００４７】
上記の反復計算を、ｆｋが収束するまで繰り返し実施することによって、原画像を推定することができる。ｆｋが収束したか否かの判定は、例えば反復計算の回数を予め設定しておいて、反復計算の回数によって判定してもよいし、ｆｋとｆｋ＋１の差分を算出して、算出される差分の絶対値の総和があるしきい値以下となるか否かで判断してもよい。
【００４８】
上記した原理を用いる画像復元処理は、以下に示す工程を順次実施していくことで実現される。以下では繰り返し計算によって推定される原画像の分布をｆｋ（ｋ＝０、１、・・・）とする。原画像の推定分布ｆｋは、復元の過程で位相に関する特性を正確に評価するために、複素関数として扱う。
まずＰＳＦ分布ｈ（ｘ、ｙ）を二次元フーリエ変換して、ＯＴＦ分布Ｈ（ｓ、ｔ）を特定する。
次に原画像の最初の推定分布として、ｆ０（ｘ、ｙ）の実部をｇ（ｘ、ｙ）として、ｆ０（ｘ、ｙ）の虚部を０とする。
次に以下に示す演算を、ｆｋ（ｘ、ｙ）が収束するまで繰り返し実施する。ここでＦＴ（）は二次元のフーリエ変換を表し、ＦＴ−１（）は二次元の逆フーリエ変換を表す。またＨ＃（ｓ、ｔ）は、Ｈ（ｓ、ｔ）の反転関数であり、Ｈ＃（ｓ、ｔ）＝Ｈ（−ｓ、−ｔ）である。
【００４９】
【数１９】
【数２０】
【数２１】
【００５０】
上記の反復計算を繰返し実施し、最終的な原画像の推定分布ｆｋを得る。最終的に得られる原画像の推定分布ｆｋ（ｘ、ｙ）が、原画像の復元画像ｆ（ｘ、ｙ）に相当する。
【００５１】
以上のように、上記の画像取得装置によれば、劣化画像の分布ｇ（ｘ、ｙ）とＰＳＦの分布ｈ（ｘ、ｙ）を用いて、Ｂａｙｅｓの理論に基づき、劣化画像の分布ｇ（ｘ、ｙ）を実現する原画像の分布ｆ（ｘ、ｙ）のうちで最大に尤もらしい分布を、反復計算によって推定することができる。この画像取得装置によれば、ＯＴＦの分布Ｈ（ｓ、ｔ）がある周波数領域で０となる場合、すなわち数式（６）による画像復元ができない場合であっても、原画像を復元することができる。ぼけやぶれのない被写体の画像を得ることができる。
【００５２】
あるいは、上記の撮像装置は、前記画像処理手段が、点像分布関数と被写体画像データに基づいて改善された被写体画像データを生成する際に、被写体画像データから劣化画像の分布を特定する処理と、原画像の最初の推定分布を特定する処理と、点像画像データに基いて取得された点像分布関数を、点像分布関数の最初の推定分布に設定する処理と、（Ａ）前記劣化画像の分布と、原画像の推定分布と、点像分布関数の推定分布に基づいて、点像分布関数の推定分布を更新する処理と、（Ｂ）前記劣化画像の分布と、原画像の推定分布と、点像分布関数の推定分布に基づいて、原画像の推定分布を更新する処理と、上記（Ａ）と（Ｂ）の工程を交互に繰返す処理と、原画像の推定分布に基づいて、改善された被写体画像データを生成する処理を実行するものであり、前記（Ａ）の点像分布関数の推定分布を更新する処理は、（Ａ１）原画像の推定分布をフーリエ変換して、原画像の推定分布のスペクトル分布を得る処理と、（Ａ２）点像分布関数の推定分布をフーリエ変換して第１の関数を得る処理と、（Ａ３）前記第１の関数に前記原画像の推定分布のスペクトル分布を乗じて第２の関数を得る処理と、（Ａ４）前記第２の関数を逆フーリエ変換して第３の関数を得る処理と、（Ａ５）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る処理と、（Ａ６）前記第４の関数をフーリエ変換して第５の関数を得る処理と、（Ａ７）前記第５の関数に前記原画像の推定分布のスペクトル分布の反転関数を乗じて第６の関数を得る処理と、（Ａ８）前記第６の関数を逆フーリエ変換して第７の関数を得る処理と、（Ａ９）前記点像分布関数の推定分布に前記第７の関数を乗じて、点像分布関数の次の推定分布を得る処理と、（Ａ１０）点像分布関数の次の推定分布を点像分布関数の推定分布に置き換える処理とを備え、前記（Ｂ）の原画像の推定分布を更新する処理は、（Ｂ１）点像分布関数の推定分布をフーリエ変換して、光学的伝達関数の推定分布を得る処理と、（Ｂ２）原画像の推定分布をフーリエ変換して第１の関数を得る処理と、（Ｂ３）前記第１の関数に前記光学的伝達関数の推定分布を乗じて第２の関数を得る処理と、（Ｂ４）前記第２の関数を逆フーリエ変換して第３の関数を得る処理と、（Ｂ５）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る処理と、（Ｂ６）前記第４の関数をフーリエ変換して第５の関数を得る処理と、（Ｂ７）前記第５の関数に前記光学的伝達関数の推定分布の反転関数を乗じて第６の関数を得る処理と、（Ｂ８）前記第６の関数を逆フーリエ変換して第７の関数を得る処理と、（Ｂ９）原画像の推定分布に前記第７の関数を乗じて、原画像の次の推定分布を得る処理と、（Ｂ１０）原画像の次の推定分布を原画像の推定分布に置き換える処理とを備えることが好ましい。
【００５３】
第２撮像手段で撮像された点像画像データから取得されるＰＳＦ分布は、第２撮像手段による撮像の際の光学的特性を反映したものであり、第１撮像手段による撮像の際の光学的特性と完全に一致している保証はない。従って、ＰＳＦ分布のわずかな相違によって、原画像を完全に復元することが困難な場合がある。そこで、上記の撮像装置では、第２撮像手段で撮像された点像画像データから取得されるＰＳＦ分布を、第１撮像手段による撮像の際の真のＰＳＦ分布の最初の推定値として扱い、原画像の分布とＰＳＦ分布の双方について、Ｂａｙｅｓの理論に基いた推定を行う。
【００５４】
上述したように、劣化画像分布ｇ（ｘ、ｙ）とＯＴＦ分布Ｈ（ｓ、ｔ）が既知であれば、未知である原画像分布ｆ（ｘ、ｙ）を復元することができる。上記と同様の原理によって、ＯＴＦ分布Ｈ（ｓ、ｔ）が未知である場合に、劣化画像分布ｇ（ｘ、ｙ）と原画像分布ｆ（ｘ、ｙ）が既知であれば、ＯＴＦ分布Ｈ（ｓ、ｔ）を復元することが可能となる。
【００５５】
数式（１５）において、ｘ２−ｘを新たにｘ’とし、ｙ２−ｙを新たにｙ’とすると、次式が得られる。
【００５６】
【数２２】
【００５７】
上式の両辺にＰ（Ａ（ｘ２、ｙ２））を乗じると、次式が得られる。
【００５８】
【数２３】
【００５９】
上式の左辺は、Ｂａｙｅｓの理論に基づくと、Ｐ（Ａ（ｘ２、ｙ２）｜Ｖ（ｘ２−ｘ’、ｙ２−ｙ’））×Ｐ（Ｖ（ｘ２−ｘ’、ｙ２−ｙ’））に等しいため、次式となる。
【００６０】
【数２４】
【００６１】
上式に、数式（１１）、数式（１２）および数式（１３）を代入すると、次式が得られる。
【００６２】
【数２５】
【００６３】
上式の両辺を、（ｘ２、ｙ２）に関して積分すると、次式が得られる。
【００６４】
【数２６】
【００６５】
上式左辺の積分は１に等しいため、結局、Ｂａｙｅｓの理論に基づくと、以下の関係が成立する。
【００６６】
【数２７】
【００６７】
上式の関係を満たすＰＳＦ分布ｈ（ｘ、ｙ）が算出されると、算出されたｈ（ｘ、ｙ）をフーリエ変換することによって、ＯＴＦ分布Ｈ（ｓ、ｔ）を推定することができる。
【００６８】
数式（２７）の右辺のｈ（ｘ、ｙ）をｈｋ（ｘ、ｙ）とし、数式（２７）の左辺のｈ（ｘ、ｙ）をｈｋ＋１（ｘ、ｙ）として、ｈｋ（ｘ、ｙ）に関する反復計算を実施し、ｈｋ（ｘ、ｙ）の収束値を求める。この反復計算は、次式で表される。
【００６９】
【数２８】
【００７０】
上記の反復計算によって求まるｈｋ（ｘ、ｙ）の収束値をフーリエ変換して、Ｂａｙｅｓの理論に基づくＯＴＦの推定分布Ｈｋ（ｘ、ｙ）を得ることができる。この関係は次式で表される。
【００７１】
【数２９】
【００７２】
上記した原理を用いるＯＴＦの推定処理は、以下に示す工程を順次実施していくことで実現される。以下では繰り返し計算によって推定されるＰＳＦの分布をｈｋ（ｋ＝０、１、・・・）とする。ＰＳＦの推定分布ｈｋは、反復計算の過程で位相に関する特性を正確に評価するために、複素関数として扱う。
まず劣化画像の分布ｇ（ｘ、ｙ）と、原画像の分布ｆ（ｘ、ｙ）を特定する。本発明の方法では、ｇ（ｘ、ｙ）の実部を劣化画像における照度の分布として、ｇ（ｘ、ｙ）の虚部を０とする。またｆ（ｘ、ｙ）の実部を原画像における照度の分布として、ｆ（ｘ、ｙ）の虚部を０とする。
次に原画像の分布ｆ（ｘ、ｙ）をフーリエ変換して、原画像の分布のスペクトル分布Ｆ（ｓ、ｔ）を算出する。
次に、ＰＳＦの最初の推定分布ｈ０（ｘ、ｙ）を特定する。
次に以下に示す演算を、ｈｋ（ｘ、ｙ）が収束するまで繰り返し実施する。ここでＦＴ（）は二次元のフーリエ変換を表し、ＦＴ−１（）は二次元の逆フーリエ変換を表す。またＦ＃（ｓ、ｔ）は、Ｆ（ｓ、ｔ）の反転関数であり、Ｆ＃（ｓ、ｔ）＝Ｆ（−ｓ、−ｔ）である。
【００７３】
【数３０】
【数３１】
【数３２】
【００７４】
ｈｋ（ｘ、ｙ）の収束の判定は、例えば反復計算の回数が予め定められた回数に到達したか否かを基準に判別してもよいし、ｈｋ（ｘ、ｙ）とｈｋ＋１（ｘ、ｙ）の差分の分布を算出して、その差分の絶対値が全ての（ｘ、ｙ）に対してあるしきい値以下となるか否かを基準に判別してもよいし、ｆｋとｆｋ＋１の差分を算出して、算出される差分の絶対値を（ｘ、ｙ）に関して積分した値があるしきい値以下となるか否かで判断してもよい。
ＰＳＦ分布ｈｋ（ｘ、ｙ）が収束したら、それをフーリエ変換して、ＯＴＦの推定分布Ｈｋ（ｓ、ｔ）を算出する。
【００７５】
上記したように、劣化画像とＯＴＦとが既知であれば原画像を復元することが可能であり、劣化画像と原画像が既知であればＯＴＦを推定することが可能である。
これらを組み合わせることによって、劣化画像のみに基づいて、ＯＴＦを推定し、原画像を復元することが可能となる。
すなわち、原画像とＰＳＦのそれぞれについて、適当な分布を仮定しておき、上記したＰＳＦの推定（すなわち、ＯＴＦの推定）と原画像の推定とを交互に繰返し実施していくことによって、原画像を復元する。
ＯＴＦに関しては、劣化画像の分布と、原画像の推定分布を用いて、Ｂａｙｅｓの理論に基づく計算によって、より改善された推定分布を得ることができる。ＯＴＦの推定に用いる原画像の推定分布が、真の原画像の分布に近いほど、真のＯＴＦの分布に近い推定分布を得ることができる。
原画像に関しては、劣化画像の分布と、ＯＴＦの推定分布を用いて、Ｂａｙｅｓの理論に基づく計算によって、より改善された推定分布を得ることができる。原画像の推定に用いるＯＴＦの推定分布が、真のＯＴＦの分布に近いほど、真の原画像の分布に近い推定分布を得ることができる。
従って、上記したＯＴＦの推定と原画像の推定を交互に繰返し実施していくことによって、原画像の推定分布はより真の原画像の分布に近づいていき、ＯＴＦの推定分布はより真のＯＴＦの分布に近づいていき、結果として良好に復元された原画像と良好に推定されたＯＴＦの双方を得ることができる。
【００７６】
すなわち、以下の工程を順に実施していくことによって、劣化画像のみから原画像を復元することができる。
まず原画像の最初の推定分布と、ＰＳＦの最初の推定分布を設定する。一般に、劣化画像の分布は原画像の分布から大きく異なることはないため、原画像の最初の推定分布としては劣化画像の分布を用いる。ＰＳＦの最初の推定分布としては、第２撮像手段によって撮像された点像画像データから取得されるＰＳＦ分布を用いる。
次に、数式（３０）〜（３２）の計算を実施して、より改善されたＰＳＦの推定分布ｈｋ（ｓ、ｔ）を取得して、得られたｈｋ（ｓ、ｔ）をフーリエ変換してＯＴＦの推定分布Ｈｋ（ｓ、ｔ）を取得する。上記の推定を実施することによって、より真のＯＴＦ分布に近いＯＴＦの推定分布を得ることができる。
次に、数式（１９）〜（２１）の計算を実施して、より改善された原画像の推定分布ｆｋ（ｘ、ｙ）を取得する。上記の推定を実施することによって、より真の原画像分布に近い原画像の推定分布を得ることができる。
上記したＯＴＦの推定と原画像の推定を、交互に繰返し実施していくことによって、推定される原画像の分布は真の原画像の分布に近づいていき、推定されるＯＴＦの分布は真のＯＴＦの分布に近づいていく。従って、上記の反復計算を実施することによって、原画像を復元することが可能となる。
【００７７】
上記の画像取得装置によれば、第１撮像手段による撮像時のＰＳＦと第２撮像手段による撮像時のＰＳＦが相違する場合であっても、ぶれのない被写体の画像を取得することができる。
【００７８】
本発明は天体撮像装置としても具現化することができる。その天体撮像装置は、天体を撮像して天体画像データを出力する第１撮像手段と、電離層に向けてレーザ光を照射するレーザ光照射手段と、レーザ光照射手段から照射されて電離層で反射したレーザ光を撮像して点像画像データを出力する第２撮像手段と、天体画像データと点像画像データに基いて改善された天体画像データを生成する画像処理手段を備えている。その天体撮像装置において、第１撮像手段と第２撮像手段は、互いに同期して撮像するものである。その天体撮像装置において、画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と天体画像データに基いて改善された天体画像データを生成する。
【００７９】
上記の天体撮像装置の動作内容を方法の発明として具現化することもできる。本発明の天体撮像方法は、第１撮像手段によって天体を撮像して天体画像データを出力する工程と、レーザ光照射手段によって電離層に向けてレーザ光を照射する工程と、レーザ光照射手段から照射されて電離層で反射したレーザ光を第２撮像手段によって撮像して点像画像データを出力する工程と、画像処理手段によって天体画像データと点像画像データに基いて改善された天体画像データを生成する工程を備えている。その天体撮像方法において、第１撮像手段と第２撮像手段は、互いに同期して撮像するものである。その天体撮像方法において、画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と天体画像データに基いて改善された天体画像データを生成する。
なお本発明は上記の天体撮像方法における各工程をコンピュータに実行させるためのプログラムとしても具現化することができる。
【００８０】
第１撮像手段で撮像される天体の画像は、大気のゆれによるシーイングの影響を受けて、ぼやけた画像となってしまう。本発明の天体撮像装置では、第１撮像手段で撮像された天体の画像と、第２撮像手段で撮像されたレーザ光の画像を用いて、シーイングの影響を除去して改善された天体の画像を取得する。第２撮像手段で撮像される画像は、レーザ光照射手段から電離層へ到達するまでの間にシーイングの影響を受け、さらに電離層で反射してから第２撮像手段へ到達するまでの間にシーイングの影響を受けた点光源の画像として扱うことができる。従って、第２撮像手段で撮像された点像画像データに基づいて、天体からの光が電離層を経由して第１撮像手段に到達するまでの間に受けるシーイングの影響を表現するＰＳＦ分布を取得することができる。第１撮像手段で撮像された天体画像データとＰＳＦ分布を用いて、シーイングの影響を排除した改善された天体画像データを得ることができる。
【００８１】
また、本発明者らは、ぶれ画像に特有の性質に着目して、ぶれ画像のみに基づいてぶれのない画像を復元する手法にも想到した。本発明のぶれ画像復元装置は、ぶれ画像の分布を特定する手段と、ぶれ画像の分布についての複素ケプストラムを算出する手段と、ぶれ画像の分布についての複素ケプストラムからパルス列を抽出する手段と、抽出されたパルス列をぶれ画像の分布についての複素ケプストラムから除去して、ぶれのない画像の分布についての複素ケプストラムを推定する手段と、ぶれのない画像の分布についての複素ケプストラムからぶれのない画像の分布を算出する手段を備えている。
【００８２】
上記のぶれ画像復元装置の動作内容を方法の発明として具現化することもできる。本発明のぶれ画像復元方法は、ぶれ画像の分布を特定する工程と、ぶれ画像の分布についての複素ケプストラムを算出する工程と、ぶれ画像の分布についての複素ケプストラムからパルス列を抽出する工程と、抽出されたパルス列をぶれ画像の分布についての複素ケプストラムから除去して、ぶれのない画像の分布についての複素ケプストラムを推定する工程と、ぶれのない画像の分布についての複素ケプストラムからぶれのない画像の分布を算出する工程を備えている。
なお本発明は上記のぶれ画像復元方法における各工程をコンピュータに実行させるためのプログラムとしても具現化することができる。
【００８３】
複素ケプストラムとは、ある関数について、フーリエ変換してスペクトル分布を取得し、そのスペクトル分布の複素対数関数を計算して、その結果を逆フーリエ変換したものをいう。これはフーリエ変換を続けて２回行ったもので代用することも出来る。以下ではぶれ画像を劣化画像と呼び、ぶれのない画像を原画像と呼ぶ。劣化画像の分布ｇ（ｘ、ｙ）、ＰＳＦの分布ｈ（ｘ、ｙ）、原画像の分布ｆ（ｘ、ｙ）には、以下の関係が成り立つ。
【００８４】
【数３３】
【００８５】
上式の両辺をフーリエ変換すると、劣化画像のスペクトル分布Ｇ（ｓ、ｔ）、ＯＴＦの分布Ｈ（ｓ、ｔ）、原画像のスペクトル分布Ｆ（ｓ、ｔ）の間で以下の関係が成り立つ。
【００８６】
【数３４】
【００８７】
上式の両辺について複素対数関数を計算すると、次式を得る。
【００８８】
【数３５】
【００８９】
上式の両辺について逆フーリエ変換することで、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）、ＰＳＦの複素ケプストラムｈｃ（ｘ、ｙ）、原画像の複素ケプストラムｆｃ（ｘ、ｙ）について、以下の関係式を得る。
【００９０】
【数３６】
【００９１】
一般のぶれ画像は、多数個のずれた画像が重なり合ったものとして扱うことができる。原画像の分布ｆ（ｘ、ｙ）に対して、全く同じ画像を、強度をｎ倍として、位置を（α、β）だけずらして重ね合わせた場合、劣化画像の分布ｇ（ｘ、ｙ）は次式で与えられる。
【００９２】
【数３７】
【００９３】
上式の両辺をフーリエ変換すると、劣化画像の分布のスペクトル分布Ｇ（ｓ、ｔ）と原画像の分布のスペクトル分布Ｆ（ｘ、ｙ）の間に、次の関係が成り立つ。
【００９４】
【数３８】
【００９５】
上式の両辺の複素対数関数を計算すると、次式が得られる。
【００９６】
【数３９】
【００９７】
ここで、ｎは１に比べて十分に小さいと仮定している。上式の両辺について逆フーリエ変換すると、次式を得る。
【００９８】
【数４０】
【００９９】
上述のように、上式の左辺ｇｃ（ｘ、ｙ）は、劣化画像の分布ｇ（ｘ、ｙ）についての複素ケプストラムである。上式の右辺第１項ｆｃ（ｘ、ｙ）は、原画像の分布ｆ（ｘ、ｙ）についての複素ケプストラムである。上式の右辺第２項以降は、ＰＳＦ分布ｈ（ｘ、ｙ）についての複素ケプストラムである。
上式から明らかなように、劣化画像の分布ｇ（ｘ、ｙ）の複素ケプストラムｇｃ（ｘ、ｙ）は、原画像の分布ｆ（ｘ、ｙ）の複素ケプストラムｆｃ（ｘ、ｙ）に、ＰＳＦ分布ｈ（ｘ、ｙ）の複素ケプストラムを重ね合わせたものとなっている。画像の劣化がぶれによるものの場合には、ＰＳＦ分布ｈ（ｘ、ｙ）の複素ケプストラムは、ずれた画像の強度ｎと、ずれ量（α、β）で定まるデルタ関数の系列（すなわちパルス列）によって表現される。従って、劣化画像の分布ｇ（ｘ、ｙ）の複素ケプストラムｇｃ（ｘ、ｙ）からこのパルス列を抽出して取り除くことができれば、原画像の分布ｆ（ｘ、ｙ）の複素ケプストラムｆｃ（ｘ、ｙ）を推定することができる。
【０１００】
劣化画像の分布ｇ（ｘ、ｙ）の複素ケプストラムｇｃ（ｘ、ｙ）からパルス列を抽出するには、例えば劣化画像の分布ｇ（ｘ、ｙ）の複素ケプストラムｇｃ（ｘ、ｙ）の振幅値がしきい値を超える個所をパルス列として識別する。パルス列が抽出されると、例えばノッチフィルタを用いることで、劣化画像の分布ｇ（ｘ、ｙ）の複素ケプストラムｇｃ（ｘ、ｙ）からパルス列を除去して、原画像の分布ｆ（ｘ、ｙ）の複素ケプストラムｆｃ（ｘ、ｙ）を推定することができる。
【０１０１】
推定された原画像の分布ｆ（ｘ、ｙ）の複素ケプストラムｆｃ（ｘ、ｙ）に関して、フーリエ変換を実施し、その結果の複素指数関数を計算することで、原画像のスペクトル分布Ｆ（ｓ、ｔ）が得られる。その後に原画像のスペクトル分布Ｆ（ｓ、ｔ）を逆フーリエ変換することで、原画像の分布ｆ（ｘ、ｙ）を復元することができる。
【０１０２】
なお上記では説明を簡単にするために２つの画像が重なり合ったぶれ画像の場合を例として説明したが、多数個の画像が重なり合っている一般的なぶれ画像の場合にも、同様の手法によって原画像の分布ｆ（ｘ、ｙ）を復元することができる。
【０１０３】
本発明の他の一つのぶれ画像復元装置は、ぶれ画像の分布を特定する手段と、ぶれ画像の分布についての複素ケプストラムを算出する手段と、ぶれ画像の分布についての複素ケプストラムからパルス列を抽出する手段と、抽出されたパルス列に基づいて点像分布関数を取得する手段と、ぶれ画像の分布と点像分布関数に基づいてぶれのない画像の分布を推定する手段を備えている。
【０１０４】
上記のぶれ画像復元装置の動作内容を方法の発明として具現化することもできる。本発明の他の一つのぶれ画像復元方法は、ぶれ画像の分布を特定する工程と、ぶれ画像の分布についての複素ケプストラムを算出する工程と、ぶれ画像の分布についての複素ケプストラムからパルス列を抽出する工程と、抽出されたパルス列に基づいて点像分布関数を取得する工程と、ぶれ画像の分布と点像分布関数に基づいてぶれのない画像の分布を推定する工程を備えている。
なお本発明は上記のぶれ画像復元方法における各工程をコンピュータに実行させるためのプログラムとしても具現化することができる。
【０１０５】
ここでもう一度、原画像の分布ｆ（ｘ、ｙ）が撮像時のぶれによって劣化画像の分布ｇ（ｘ、ｙ）となっている場合を考える。原画像の分布ｆ（ｘ、ｙ）に対して、全く同じ画像を強度をｎ倍として、位置を（α、β）だけずらして重ね合わせた場合、劣化画像の分布ｇ（ｘ、ｙ）は次式で与えられる。
【０１０６】
【数４１】
【０１０７】
すなわち、この場合のＰＳＦ分布ｈ（ｘ、ｙ）は、次式で与えられる。
【０１０８】
【数４２】
【０１０９】
数式（４０）と数式（４２）を対比すると、数式（４２）のＰＳＦの分布ｈ（ｘ、ｙ）の右辺第２項は、数式（４０）の劣化画像の分布ｇ（ｘ、ｙ）の複素ケプストラムｇｃ（ｘ、ｙ）に現れるパルス列のうち最大のもの（数式（４０）の右辺第２項）と同様に、位置が（α、β）であって、大きさがｎのデルタ関数となっていることが分かる。従って、劣化画像の分布ｇ（ｘ、ｙ）の複素ケプストラムｇｃ（ｘ、ｙ）から抽出されたパルス列のうち最大のものに、強度が１で原点に位置するデルタ関数を重ね合わせることで、ぶれの影響を示すＰＳＦの分布ｈ（ｘ、ｙ）を直接的に取得することができる。
【０１１０】
あるいは、ｎが１より十分に小さければ、劣化画像の分布ｇ（ｘ、ｙ）の複素ケプストラムｇｃ（ｘ、ｙ）に現れるパルス列に、強度が１で原点に位置するデルタ関数を単に重ね合わせることで、ＰＳＦの分布ｈ（ｘ、ｙ）の近似解を取得することができる。
【０１１１】
ぶれの影響を表現するＰＳＦの分布ｈ（ｘ、ｙ）が取得されれば、劣化画像の分布ｇ（ｘ、ｙ）とＰＳＦの分布ｈ（ｘ、ｙ）に基づいて、原画像の分布ｆ（ｘ、ｙ）を推定することができる。ぶれ画像のみに基づいて、ぶれのない画像を取得することができる。
【０１１２】
なお上記では説明を簡単にするために２つの画像が重なり合ったぶれ画像の場合を例として説明したが、多数個の画像が重なり合っている一般的なぶれ画像の場合にも、同様の手法によって原画像の分布ｆ（ｘ、ｙ）を復元することができる。
【発明の効果】
【０１１３】
本発明の装置、方法およびプログラムによれば、撮像時のぶれやぼけの影響を表現する点像分布関数を適切に取得して、劣化した画像から劣化のない画像を復元することができる。
【図面の簡単な説明】
【０１１４】
【図１】図１は実施例１の撮像装置１００の構成を示す図である。
【図２】図２は実施例１の撮像装置１００の動作を説明するフローチャートである。
【図３】図３は実施例１の撮像装置１００によってぶれのない被写体画像を取得する様子を示している。
【図４】図４は実施例１の変形例１の撮像装置１５０における画像復元処理を説明するフローチャートである。
【図５】図５は実施例１の変形例２の撮像装置１６０における画像復元処理を説明するフローチャートである。
【図６】図６は実施例１の変形例２の撮像装置１６０におけるＰＳＦ分布の推定処理を説明するフローチャートである。
【図７】図７は実施例１の変形例２の撮像装置１６０における原画像分布の推定処理を説明するフローチャートである。
【図８】図８は実施例２の天体撮像装置８００の構成を示す図である。
【図９】図９は実施例２の天体撮像装置８００によってシーイングの影響を排除した天体画像を取得する様子を示す図である。
【図１０】図１０は実施例２の天体撮像装置８００の動作を説明するフローチャートである。
【図１１】図１１は実施例３の画像復元装置１１００の構成を示す図である。
【図１２】図１２は実施例３の画像復元装置１１００の処理を説明するフローチャートである。
【図１３】図１３は実施例３の画像復元装置１１００で計算されるぶれ画像の複素ケプストラム分布の概形を示す図である。
【図１４】図１４は実施例３の画像復元装置１１００の機能ブロック図である。
【図１５】図１５は実施例４の画像復元装置１１５０の動作を説明するフローチャートである。
【図１６】図１６は実施例４の画像復元装置１１５０の機能ブロック図である。
【発明を実施するための形態】
【０１１５】
以下、本発明を具現化した実施例について図面を参照して説明する。
【０１１６】
（実施例１）
図面を参照しながら、本実施例の撮像装置１００について説明する。
図１は本実施例の撮像装置１００の構成を示す概略図である。本実施例の撮像装置１００は、撮像台１１０に載置された被写体Ｗを撮像する第１イメージセンサ１０４と、撮像台１１０に固定されたレーザ光源１０８と、レーザ光源１０８を撮像する第２イメージセンサ１０２と、制御装置１１２を備えている。本実施例では被写体Ｗは表面に半導体回路が形成されたフォトマスクであって、撮像装置１００は被写体Ｗの表面に形成された回路の形状を撮像する。
【０１１７】
第１イメージセンサ１０４と第２イメージセンサ１０２は、撮像台１１０の上面に対して平行な撮像面を備えている。第１イメージセンサ１０４と第２イメージセンサ１０２は、それぞれ撮像面に複数の撮像素子を備えており、撮像面に形成される像を撮像素子を介して認識し、ＲＧＢの各色成分の照度の分布を示す電気信号に変換して画像処理装置１１２へ送信する。第１イメージセンサ１０４と第２イメージセンサ１０２は、両者が一体的に移動するように、センサ支持台１０６に搭載されている。第１イメージセンサ１０４と第２イメージセンサ１０２は、レンズ収差等の特性が同一である。また、第１イメージセンサ１０４と第２イメージセンサ１０２は、それぞれのシャッタが同期して開閉するように制御装置１１２によって制御される。従って、撮像台１１０とセンサ支持台１０６の間のわずかな振動によって撮像時に第１イメージセンサ１０４で撮像される被写体Ｗの画像がぶれてしまう場合には、第２イメージセンサ１０２で撮像されるレーザ光源１０８の画像も同じようにぶれることになる。
【０１１８】
第２イメージセンサ１０２で撮像されるレーザ光源１０８の画像は、撮像台１１０に固定された点光源を撮像した画像ということができる。従って、第２イメージセンサ１０２で撮像される画像は、第２イメージセンサ１０２で撮像する際の点像分布関数（ＰＳＦ：Point Spread Function）の分布として扱うことができる。第２イメージセンサ１０２と第１イメージセンサ１０４は同一の特性を有しており、それぞれのシャッタが同期して開閉するように撮像するので、第２イメージセンサ１０２で撮像する際のＰＳＦの分布は、第１イメージセンサ１０４で撮像する際のＰＳＦの分布として扱うことができる。レーザ光源１０８の発光は制御装置１１２によって制御される。
【０１１９】
制御装置１１２は、第１イメージセンサ１０４から送信されるデータと、第２イメージセンサ１０２から送信されるデータに基づいて、被写体Ｗのぶれのない画像を生成する。
【０１２０】
以下では図２のフローチャートを参照しながら、撮像装置１００によって被写体Ｗのぶれのない画像を取得する方法について説明する。以下の説明では、撮像時にぶれてしまった被写体Ｗの画像を劣化画像と呼び、ぶれのない被写体Ｗの画像を原画像と呼ぶ。
【０１２１】
ステップＳ２０２では、被写体Ｗを撮像台１１０に載置し、位置を固定する。
【０１２２】
ステップＳ２０４では、点光源であるレーザ光源１０８の発光を開始する。
【０１２３】
ステップＳ２０６では、第１イメージセンサ１０４と第２イメージセンサ１０２を同期させて撮像する。第１イメージセンサ１０４では、被写体Ｗを撮像したデータが取得される。第２イメージセンサ１０２では点光源であるレーザ光源１０８を撮像したデータが取得される。第１イメージセンサ１０４で撮像されたデータは、被写体画像データとして制御装置１１２に記録される。第２イメージセンサ１０２で撮像されたデータは、点像画像データとして制御装置１１２に記録される。
【０１２４】
被写体画像データと点像画像データは、いずれもＲＧＢそれぞれの色成分の照度の分布を表現している。本実施例では、ＲＧＢそれぞれの色成分に関して、ステップＳ２０８とステップＳ２１０の処理を実施して、それぞれの色成分に関して復元された照度の分布を組み合わせることで、原画像を生成する。
【０１２５】
ステップＳ２０８では、点像画像データから、ＰＳＦの分布ｈ（ｘ、ｙ）を取得する。本実施例では、点像画像データの照度の分布を正規化して、ＰＳＦの分布ｈ（ｘ、ｙ）とする。
【０１２６】
ステップＳ２１０では、被写体画像データにおける照度の分布を劣化画像の分布ｇ（ｘ、ｙ）とし、劣化画像の分布ｇ（ｘ、ｙ）と、ステップＳ２０８で取得されたＰＳＦの分布ｈ（ｘ、ｙ）を用いて、画像復元処理を行う。
【０１２７】
本実施例では、ステップＳ２１０において、以下のようにして画像復元処理を行う。
まず劣化画像の分布ｇ（ｘ、ｙ）をフーリエ変換して、劣化画像のスペクトル分布Ｇ（ｓ、ｔ）を取得する。
【０１２８】
そして、ステップＳ２０８で取得されたＰＳＦの分布ｈ（ｘ、ｙ）をフーリエ変換して、光学的伝達関数（ＯＴＦ：Optical Transfer Function）の分布Ｈ（ｓ、ｔ）を取得する。
【０１２９】
その後に、劣化画像のスペクトル分布Ｇ（ｓ、ｔ）をＯＴＦの分布Ｈ（ｓ、ｔ）で除して、原画像のスペクトル分布Ｆ（ｓ、ｔ）を得る。
【０１３０】
最後に、原画像のスペクトル分布Ｆ（ｓ、ｔ）を逆フーリエ変換することで、原画像の分布ｆ（ｘ、ｙ）が得られる。
【０１３１】
上述のように原画像の分布ｆ（ｘ、ｙ）をＲＧＢの色成分それぞれについて取得し、それぞれの色成分の照度の分布を組み合わせることで、劣化画像から原画像を復元することができる。すなわち、撮像時にぶれてしまった被写体Ｗの画像から、ぶれのない被写体Ｗの画像を得ることができる。ぶれのない被写体Ｗの画像は、例えばディスプレイに表示してもよいし、画像データとしてハードディスクなどの記憶装置に記憶してもよいし、通信回線を経由して他のコンピュータに送信してもよい。
【０１３２】
図３は本実施例の撮像装置１００によって被写体Ｗのぶれのない画像を取得する様子を示している。本実施例の撮像装置１００によれば、第１イメージセンサ１０４で撮像された被写体画像３０２と、第２イメージセンサ１０２で撮像された点像画像３０４を用いて、ぶれのない被写体画像３０６を取得することができる。撮像台１１０とセンサ支持台１０６の間でわずかな振動があっても、ぶれのない被写体画像３０６を取得することができる。撮像台１１０やセンサ支持台１０６に防振機構を設けることなく、被写体Ｗのぶれのない画像を得ることが可能となる。
【０１３３】
（実施例１の変形例１）
本実施例の撮像装置１５０は、図１に示す実施例１の撮像装置１００と同様のハードウェア構成を備えているが、制御装置１１２が図２のステップＳ２１０において行う画像復元処理のみが異なっている。以下では実施例１との相違点のみについて詳細に説明する。
【０１３４】
図４は本実施例の撮像装置１５０において図２のステップＳ２１０で行う画像復元処理の詳細を示すフローチャートである。
【０１３５】
ステップＳ４０２では、劣化画像の分布ｇ（ｘ、ｙ）を特定する。本実施例の劣化画像の分布ｇ（ｘ、ｙ）は、光の照度の分布を位相特性も含めて記述する複素関数である。本実施例では、ｇ（ｘ、ｙ）の実部を被写体画像データにおける照度の分布と一致させ、ｇ（ｘ、ｙ）の虚部を全て０とする。
【０１３６】
次にステップＳ４０４では、ＰＳＦの分布ｈ（ｘ、ｙ）をフーリエ変換して、光学的伝達関数（ＯＴＦ：Optical Transfer Function）の分布Ｈ（ｓ、ｔ）を特定する。
【０１３７】
ステップＳ４０６では、ステップＳ４０４で特定された劣化画像の分布ｇ（ｘ、ｙ）を原画像の最初の推定分布ｆ０（ｘ、ｙ）に設定する。また反復計算の繰り返し数ｋを０に設定する。
【０１３８】
ステップＳ４０８では、原画像の推定分布ｆｋ（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ１（ｓ、ｔ）に設定する。関数Ｋ１（ｓ、ｔ）は、第１の関数に相当する。ｋは非負の整数であり、反復計算の回数に応じて後のステップＳ４２６で増加していく。
【０１３９】
ステップＳ４１０では、ステップＳ４０８で設定した第１の関数Ｋ１（ｓ、ｔ）に、ステップＳ４０４で特定したＨ（ｓ、ｔ）を乗じて、関数Ｋ２（ｓ、ｔ）を算出する。関数Ｋ２（ｓ、ｔ）は、第２の関数に相当する。
【０１４０】
ステップＳ４１２では、ステップＳ４１０で算出した第２の関数Ｋ２（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ３（ｘ、ｙ）に設定する。関数Ｌ３（ｘ、ｙ）は、複素関数であって、第３の関数に相当する。
【０１４１】
ステップＳ４１４では、ステップＳ４０２で特定した劣化画像分布ｇ（ｘ、ｙ）を、ステップＳ４１２で設定した第３の関数Ｌ３（ｘ、ｙ）で除して、関数Ｌ４（ｘ、ｙ）を算出する。関数Ｌ４（ｘ、ｙ）は、複素関数であって、第４の関数に相当する。
【０１４２】
ステップＳ４１６では、ステップＳ４１４で算出した第４の関数Ｌ４（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ５（ｓ、ｔ）に設定する。関数Ｋ５（ｓ、ｔ）は、第５の関数に相当する。
【０１４３】
ステップＳ４１８では、ステップＳ４１６で設定した第５の関数Ｋ５（ｓ、ｔ）に、Ｈ＃（ｓ、ｔ）を乗じて、関数Ｋ６（ｓ、ｔ）を算出する。Ｈ＃（ｓ、ｔ）は、ステップＳ４０４で特定したＨ（ｓ、ｔ）の反転関数であって、Ｈ＃（ｓ、ｔ）＝Ｈ（−ｓ、−ｔ）の関係を満たす。関数Ｋ６（ｓ、ｔ）は、第６の関数に相当する。
【０１４４】
ステップＳ４２０では、ステップＳ４１８で設定した第６の関数Ｋ６（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ７（ｘ、ｙ）に設定する。関数Ｌ７（ｘ、ｙ）は、複素関数であって、第７の関数に相当する。
【０１４５】
ステップＳ４２２では、原画像の推定分布ｆｋ（ｘ、ｙ）に、ステップＳ４２０で設定した第７の関数Ｌ７（ｘ、ｙ）を乗じて、原画像の改善された推定分布ｆｋ＋１（ｘ、ｙ）を算出する。前記第７の関数は一般に虚部を備える複素関数となるから、上記の方法によって原画像の推定分布を位相特性も含めて改善することができる。
【０１４６】
ステップＳ４２４では、ステップＳ４２２で算出された原画像の改善された推定分布ｆｋ＋１（ｘ、ｙ）と、原画像の推定分布ｆｋ（ｘ、ｙ）の差分を算出し、その絶対値が全ての（ｘ、ｙ）に対してしきい値εを下回るか否かを判断する。前記差分の絶対値がある（ｘ、ｙ）に対してしきい値ε以上の場合（ステップＳ４２４でＮＯの場合）、原画像の改善された推定分布ｆｋ＋１（ｘ、ｙ）はまだ収束に達していないと判断し、処理はステップＳ４２６へ進む。前記差分の絶対値が全ての（ｘ、ｙ）に対してしきい値εを下回る場合（ステップＳ４２４でＹＥＳの場合）、原画像の改善された推定分布ｆｋ＋１（ｘ、ｙ）は収束に達したものと判断して、処理はステップＳ４２８へ進む。
【０１４７】
ステップＳ４２６では、反復計算の繰り返し数ｋを１増加させる。処理はステップＳ４０８へ戻り、ステップＳ４０８からステップＳ４２４までの処理を再度実施する。
【０１４８】
ステップＳ４２８では、反復計算の結果得られた原画像の推定分布ｆｋ＋１（ｘ、ｙ）を原画像の分布ｆ（ｘ、ｙ）として出力する。
【０１４９】
以上の原画像の復元処理をＲＧＢの色成分それぞれについて行うことで、被写体Ｗのぶれのない画像を取得することができる。取得されたぶれのない画像は、例えばディスプレイに表示してもよいし、画像データとしてハードディスクなどの記憶装置に記憶してもよいし、通信回線を経由して他のコンピュータに送信してもよい。
【０１５０】
本実施例の撮像装置１５０によれば、劣化画像の分布ｇ（ｘ、ｙ）とＰＳＦの分布ｈ（ｘ、ｙ）を用いて、Ｂａｙｅｓの理論に基づき、劣化画像の分布ｇ（ｘ、ｙ）を実現する原画像の分布ｆ（ｘ、ｙ）のうちで最大に尤もらしい分布を、反復計算によって推定することができる。この撮像装置１５０は、ＯＴＦの分布Ｈ（ｓ、ｔ）がある周波数領域で０となる場合、すなわち実施例１の撮像装置１００による画像復元ができない場合であっても、原画像を復元することができる。
【０１５１】
（実施例１の変形例２）
本実施例の撮像装置１６０は、図１に示す実施例１の撮像装置１００や実施例１の変形例１の撮像装置１５０と同様のハードウェア構成を備えているが、制御装置１１２が図２のステップＳ２１０において行う画像復元処理の手法のみが異なっている。以下では実施例１および実施例１の変形例１との相違点のみについて詳細に説明する。
【０１５２】
図５は本実施例の撮像装置１６０において図２のステップＳ２１０で行う画像復元処理の詳細を示すフローチャートである。
【０１５３】
ステップＳ５０２では、劣化画像の分布ｇ（ｘ、ｙ）を特定する。本実施例の劣化画像の分布ｇ（ｘ、ｙ）は、光の照度の分布を位相特性も含めて記述する複素関数である。本実施例では、ｇ（ｘ、ｙ）の実部を被写体画像データにおける照度の分布と一致させ、ｇ（ｘ、ｙ）の虚部を全て０とする。
【０１５４】
ステップＳ５０４では、原画像の最初の推定分布ｆ０（ｘ、ｙ）を設定する。本実施例のでは、原画像の最初の推定分布ｆ０（ｘ、ｙ）として、劣化画像の分布ｇ（ｘ、ｙ）を用いる。一般的に、原画像の分布と劣化画像の分布は大きく異ならないと考えられるため、上記のように原画像の最初の推定分布ｆ０（ｘ、ｙ）を設定することによって、原画像の復元に伴う反復計算の回数を低減することができる。
またステップＳ５０４では、ＰＳＦの最初の推定分布ｈ０（ｘ、ｙ）を設定する。本実施例では、ＰＳＦの最初の推定分布ｈ０（ｘ、ｙ）として、図２のステップＳ２０８で取得されているＰＳＦの分布ｈ（ｘ、ｙ）を用いる。第２イメージセンサ１０２で撮像する際のＰＳＦの分布は、第１イメージセンサ１０４で撮像する際のＰＳＦの分布と大きく異ならないと考えられるため、上記のようにＰＳＦの最初の推定分布ｈ０（ｘ、ｙ）を設定することによって、原画像の復元に伴う反復計算の回数を低減することができる。
さらにステップＳ５０４では、反復計算の繰返し数ｋを０に設定する。
【０１５５】
ステップＳ５０６では、反復計算の繰返し数ｋを１増加させる。
【０１５６】
ステップＳ５０８では、ＰＳＦの推定分布ｈｋ（ｘ、ｙ）を算出して、ＰＳＦの推定分布を更新する。この計算については後述する。
【０１５７】
ステップＳ５１０では、ＰＳＦの推定分布ｈｋ（ｘ、ｙ）をフーリエ変換して、ＯＴＦの推定分布Ｈｋ（ｓ、ｔ）を算出する。
【０１５８】
ステップＳ５１２では、原画像の推定分布ｆｋ（ｘ、ｙ）を算出して、原画像の推定分布を更新する。この計算については後述する。
【０１５９】
ステップＳ５１４では、ステップＳ５１２で更新された原画像の推定分布ｆｋ（ｘ、ｙ）と、更新される前の原画像の推定分布ｆｋ−１との差分を算出し、その差分がすべての（ｘ、ｙ）に対してあるしきい値εを下回るか否かを判断する。前記差分がある（ｘ、ｙ）に対してしきい値ε以上の場合（ステップＳ５１４でＮＯの場合）、原画像の推定分布ｆｋ（ｘ、ｙ）は未だ収束に達していないと判断して、処理はステップＳ５０６へ戻り、ステップＳ５０６からステップＳ５１４までの処理を再度繰返す。前記差分がすべての（ｘ、ｙ）に対してしきい値εを下回る場合（ステップＳ５１４でＹＥＳの場合）、原画像の推定分布ｆｋ（ｘ、ｙ）は収束に達したと判断して、処理はステップＳ５１６へ進む。
【０１６０】
ステップＳ５１６では、反復計算の結果得られた原画像の推定分布ｆｋ（ｘ、ｙ）を、復元された原画像の分布ｆ（ｘ、ｙ）に設定する。本実施例の方法では、反復計算の過程で画像を記述する分布を複素関数として扱っているため、上記で設定される原画像の分布ｆ（ｘ、ｙ）は一般に複素関数となる。本実施例では、復元された原画像の分布ｆ（ｘ、ｙ）の実部が、原画像における照度の分布を表す。
【０１６１】
以下では図５のステップＳ５０８のｈｋ（ｘ、ｙ）の推定について詳細に説明する。本実施例では、ｈｋ（ｘ、ｙ）の推定として、図６のフローチャートに示すステップＳ６０２〜Ｓ６２６の工程を実施する。
【０１６２】
ステップＳ６０２では、ｈｋ（ｘ、ｙ）の最初の推定分布ｈｋ０（ｘ、ｙ）を設定する。本実施例の方法では、すでに取得されているｈｋ−１（ｘ、ｙ）を、ｈｋ（ｘ、ｙ）の最初の推定ｈｋ０（ｘ、ｙ）とする。また、ステップＳ６０２では、ｈｋ（ｘ、ｙ）の推定に係る反復計算の繰返し数ｍを０に設定する。
【０１６３】
ステップＳ６０４では、原画像の推定スペクトル分布Ｆ（ｓ、ｔ）を設定する。本実施例では、すでに取得されている原画像の推定分布ｆｋ−１（ｘ、ｙ）をフーリエ変換して、その結果をＦ（ｓ、ｔ）とする。
【０１６４】
ステップＳ６０６では、ｈｋ（ｘ、ｙ）の推定分布ｈｋｍ（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ１（ｓ、ｔ）に設定する。関数Ｋ１（ｓ、ｔ）は、第１の関数に相当する。
【０１６５】
ステップＳ６０８では、ステップＳ６０６で算出された関数Ｋ１（ｓ、ｔ）に、ステップＳ６０４で算出された原画像の推定スペクトル分布Ｆ（ｓ、ｔ）を乗じて、関数Ｋ２（ｓ、ｔ）を算出する。関数Ｋ２（ｓ、ｔ）は、第２の関数に相当する。
【０１６６】
ステップＳ６１０では、ステップＳ６０８で算出された関数Ｋ２（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ３（ｘ、ｙ）に設定する。関数Ｌ３（ｘ、ｙ）は、第３の関数に相当する。
【０１６７】
ステップＳ６１２では、図５のステップＳ５０４で特定された劣化画像分布ｇ（ｘ、ｙ）を、図６のステップＳ６１０で算出された関数Ｌ３（ｘ、ｙ）で除して、関数Ｌ４（ｘ、ｙ）を算出する。関数Ｌ４（ｘ、ｙ）は、第４の関数に相当する。
【０１６８】
ステップＳ６１４では、ステップＳ６１２で算出された関数Ｌ４（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ５（ｓ、ｔ）に設定する。関数Ｋ５（ｓ、ｔ）は、第５の関数に相当する。
【０１６９】
ステップＳ６１６では、ステップＳ６１４で算出された関数Ｋ５（ｓ、ｔ）に、Ｆ＃（ｓ、ｔ）を乗じて、関数Ｋ６（ｓ、ｔ）を算出する。関数Ｋ６（ｓ、ｔ）は、第６の関数に相当する。Ｆ＃（ｓ、ｔ）は、ステップＳ６０４で算出された原画像の推定スペクトル分布Ｆ（ｓ、ｔ）の反転関数であり、Ｆ＃（ｓ、ｔ）＝Ｆ（−ｓ、−ｔ）の関係を満たす。
【０１７０】
ステップＳ６１８では、ステップＳ６１６で算出された関数Ｋ６（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ７（ｘ、ｙ）に設定する。関数Ｌ７（ｘ、ｙ）は、第７の関数に相当する。
【０１７１】
ステップＳ６２０では、ｈｋ（ｘ、ｙ）の推定分布ｈｋｍ（ｘ、ｙ）に、ステップＳ６１８で設定された関数Ｌ７（ｘ、ｙ）を乗じて、その結果をｈｋ（ｘ、ｙ）の改善された推定分布ｈｋｍ＋１（ｘ、ｙ）に設定する。
【０１７２】
ステップＳ６２２では、ｈｋ（ｘ、ｙ）の推定のための反復計算の繰返し数ｍを１増加させる。
【０１７３】
ステップＳ６２４では、ｈｋ（ｘ、ｙ）の推定のための反復計算の繰返し数ｍが、しきい値以上になったか否かを判断する。本実施例では、ｈｋ（ｘ、ｙ）の推定のために実施する反復計算の回数を５回としている。繰り返し数ｍが５未満の場合（ステップＳ６２４でＮＯの場合）、処理はステップＳ６０６へ進み、ステップＳ６０６からステップＳ６２２までの処理を再度実施する。繰り返し数ｍが５以上の場合（ステップＳ６２４でＹＥＳの場合）、処理はステップＳ６２６へ進む。
【０１７４】
ステップＳ６２６では、上記の反復計算の結果得られる分布ｈｋ５（ｘ、ｙ）を、ｈｋ（ｘ、ｙ）の推定分布として設定する。
上記の反復計算は、劣化画像分布ｇ（ｘ、ｙ）と、原画像分布ｆｋ−１（ｘ、ｙ）から、ＰＳＦの分布ｈｋ（ｘ、ｙ）を推定することに相当する。原画像分布ｆｋ−１（ｘ、ｙ）が、真の原画像分布ｆ（ｘ、ｙ）に近いほど、上記の反復計算によって推定されるｈｋ（ｘ、ｙ）は、真のＰＳＦの分布ｈ（ｘ、ｙ）に近づく。
【０１７５】
以下では図５のステップＳ５１４のｆｋ（ｘ、ｙ）の推定について詳細に説明する。本実施例の方法では、ｆｋ（ｘ、ｙ）の推定として、図７のフローチャートに示すステップＳ７０２〜Ｓ７２４の工程を実施する。
【０１７６】
ステップＳ７０２では、ｆｋ（ｘ、ｙ）の最初の推定分布ｆｋ０（ｘ、ｙ）を設定する。本実施例の方法では、既に取得されているｆｋ−１（ｘ、ｙ）を、ｆｋ（ｘ、ｙ）の最初の推定分布ｆｋ０（ｘ、ｙ）として設定する。また、ステップＳ７０２では、ｆｋ（ｘ、ｙ）の推定のための反復計算の繰返し数ｎを０に設定する。
【０１７７】
ステップＳ７０４では、ｆｋ（ｘ、ｙ）の推定分布ｆｋｎ（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ１（ｓ、ｔ）に設定する。関数Ｋ１（ｓ、ｔ）は、第１の関数に相当する。
【０１７８】
ステップＳ７０６では、ステップＳ７０４で算出された関数Ｋ１（ｓ、ｔ）に、図５のステップＳ５１０で算出されたＯＴＦの推定分布Ｈｋ（ｓ、ｔ）を乗じて、関数Ｋ２（ｓ、ｔ）を算出する。関数Ｋ２（ｓ、ｔ）は、第２の関数に相当する。
【０１７９】
ステップＳ７０８では、ステップＳ７０６で算出された関数Ｋ２（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ３（ｘ、ｙ）に設定する。関数Ｌ３（ｘ、ｙ）は、第３の関数に相当する。
【０１８０】
ステップＳ７１０では、図５のステップＳ５０４で特定された劣化画像分布ｇ（ｘ、ｙ）を、図７のステップＳ７０８で算出された関数Ｌ３（ｘ、ｙ）で除して、関数Ｌ４（ｘ、ｙ）を算出する。関数Ｌ４（ｘ、ｙ）は、第４の関数に相当する。
【０１８１】
ステップＳ７１２では、ステップＳ７１０で算出された関数Ｌ４（ｘ、ｙ）をフーリエ変換して、その結果を関数Ｋ５（ｓ、ｔ）に設定する。関数Ｋ５（ｓ、ｔ）は、第５の関数に相当する。
【０１８２】
ステップＳ７１４では、ステップＳ７１２で算出された関数Ｋ５（ｓ、ｔ）に、Ｈｋ＃（ｓ、ｔ）を乗じて、関数Ｋ６（ｓ、ｔ）を算出する。Ｈｋ＃（ｓ、ｔ）は、図５のステップＳ５１２で算出されたＯＴＦの推定分布Ｈｋ（ｓ、ｔ）の反転関数であり、Ｈｋ＃（ｓ、ｔ）＝Ｈｋ（−ｓ、−ｔ）の関係を満たす。
【０１８３】
ステップＳ７１６では、ステップＳ７１４で算出された関数Ｋ６（ｓ、ｔ）を逆フーリエ変換して、その結果を関数Ｌ７（ｘ、ｙ）に設定する。関数Ｌ７（ｘ、ｙ）は、第７の関数に相当する。
【０１８４】
ステップＳ７１８では、ｆｋ（ｘ、ｙ）の推定分布ｆｋｎ（ｘ、ｙ）に、ステップＳ７１６で算出された関数Ｌ７（ｘ、ｙ）を乗じて、ｆｋ（ｘ、ｙ）の改善された推定分布ｆｋｎ＋１（ｘ、ｙ）を算出する。
【０１８５】
ステップＳ７２０では、ｆｋ（ｘ、ｙ）の推定のための反復計算の繰返し数ｎを１増加する。
【０１８６】
ステップＳ７２２では、ｆｋ（ｘ、ｙ）の推定のための反復計算の繰返し数ｎが、しきい値以上になったか否かを判断する。本実施例では、ｆｋ（ｘ、ｙ）の推定のために実施する反復計算の回数を５回としている。繰り返し数ｎが５未満の場合（ステップＳ７２２でＮＯの場合）、処理はステップＳ７０４へ進み、ステップＳ７０４からステップＳ７２０までの処理を再度実施する。繰り返し数ｎが５以上の場合（ステップＳ７２２でＹＥＳの場合）、処理はステップＳ７２４へ進む。
【０１８７】
ステップＳ７２４では、上記の反復計算の結果得られる分布ｆｋ５（ｘ、ｙ）を、ｆｋ（ｘ、ｙ）の推定分布として設定する。
上記の反復計算は、劣化画像分布ｇ（ｘ、ｙ）と、ＯＴＦ分布Ｈｋ（ｓ、ｔ）とに基づいて、原画像分布ｆｋ（ｘ、ｙ）を推定することに相当する。ＯＴＦ分布Ｈｋ（ｓ、ｔ）が、真のＯＴＦ分布Ｈ（ｓ、ｔ）に近いほど、上記の反復計算によって推定されるｆｋ（ｘ、ｙ）は、真の原画像の分布ｆ（ｘ、ｙ）に近づく。
【０１８８】
上記の実施例では、反復計算の過程において、先にＯＴＦの推定分布を更新して、次に原画像の推定分布を更新する例を説明した。ＯＴＦの推定分布の更新と、原画像の推定分布の更新は、交互に繰返し実施すればよいため、先に原画像の推定分布を更新して、次にＯＴＦの推定分布を更新してもよい。
【０１８９】
以上の原画像の復元処理をＲＧＢの色成分それぞれについて行うことで、被写体Ｗのぶれのない画像を取得することができる。取得されたぶれのない画像は、例えばディスプレイに表示してもよいし、画像データとしてハードディスクなどの記憶装置に記憶してもよいし、通信回線を経由して他のコンピュータに送信してもよい。
【０１９０】
本実施例の撮像装置１６０によれば、劣化画像の分布ｇ（ｘ、ｙ）とＰＳＦの分布ｈ（ｘ、ｙ）を用いて、Ｂａｙｅｓの理論に基づき、劣化画像の分布ｇ（ｘ、ｙ）を実現する原画像の分布ｆ（ｘ、ｙ）のうちで最大に尤もらしい分布を、反復計算によって推定することができる。この撮像装置１６０は、ＯＴＦの分布Ｈ（ｓ、ｔ）がある周波数領域で０となる場合、すなわち実施例１の撮像装置１００による画像復元ができない場合であっても、原画像を復元することができる。
【０１９１】
本実施例の撮像装置１６０によれば、第２イメージセンサ１０２で撮像された点像画像データに基いて取得されたＰＳＦの分布ｈ（ｘ、ｙ）が、第１イメージセンサ１０４で被写体Ｗを撮像する際の真のＰＳＦの分布と相違する場合であっても、劣化画像の分布ｇ（ｘ、ｙ）を実現するＰＳＦの分布ｈ（ｘ、ｙ）のうちで最大に尤もらしい分布を、反復計算によって推定することができる。実施例１の変形例１の撮像装置１５０にくらべ、より正確に原画像を復元することができる。
【０１９２】
（実施例２）
図面を参照しながら、本実施例の天体撮像装置８００について説明する。
図８は本実施例の天体撮像装置８００の構成を示す概略図である。本実施例の天体撮像装置８００は、天体を撮像する第１カメラ８０２と、電離層に向けてレーザ光を照射するレーザ照射器８０６と、レーザ照射器８０６から照射されて電離層で反射したレーザ光を撮像する第２カメラ８０４と、制御装置８０８を備えている。
【０１９３】
第１カメラ８０２と第２カメラ８０４とレーザ照射器８０６は、同一の地域内で、地上に載置されている。
【０１９４】
第１カメラ８０２と第２カメラ８０４は、それぞれ撮像面に複数の撮像素子を備えており、撮像面に形成される像を撮像素子を介して認識し、ＲＧＢの各色成分の照度の分布を示す電気信号に変換して制御装置８０８へ送信する。第１カメラ８０２と第２カメラ８０４は、レンズ収差等の特性が同一である。また、第１カメラ８０２と第２カメラ８０４は、それぞれのシャッタが同期して開閉するように制御装置８０８によって制御される。
【０１９５】
レーザ照射器８０６は、制御装置８０８によって制御されており、電離層に向けてレーザ光を照射する。レーザ照射器８０６で照射されたレーザ光は、電離層で反射して、第２カメラ８０４で受光される。
【０１９６】
第１カメラ８０２が撮像する天体の画像は、大気のゆれによるシーイングの影響を受けて、ぼやけた画像となってしまう。本実施例の天体撮像装置８００では、第１カメラ８０２で撮像された天体の画像と、第２カメラ８０４で撮像されたレーザ光の画像を用いて、シーイングの影響を除去した天体の画像を取得する。第２カメラ８０４で撮像される画像は、レーザ照射器８０６から電離層へ到達するまでの間にシーイングの影響を受け、さらに電離層で反射してから第２カメラ８０４へ到達するまでの間にシーイングの影響を受けた点光源の画像として扱うことができる。従って、第２カメラ８０４で撮像される画像に基づいて、第２カメラ８０４で天体を撮像する際のＰＳＦの分布を取得することができる。第２カメラ８０４と第１カメラ８０２は同一の特性を有しており、それぞれのシャッタが同期して開閉するように撮像するので、第２カメラ８０４で撮像する際のＰＳＦの分布は、第１カメラ８０２で撮像する際のＰＳＦの分布として扱うことができる。
【０１９７】
図９は本実施例の天体撮像装置８００を用いて惑星を撮像した様子を示している。図９に示すように、本実施例の天体撮像装置８００によれば、大気のゆれによるシーイングの影響によってぼやけてしまっている天体画像９０２から、シーイングの影響を排除した鮮明な天体画像９０４を取得することができる。
【０１９８】
図１０のフローチャートを参照しながら、天体撮像装置８００の動作について説明する。
【０１９９】
ステップＳ１００２では、レーザ照射器８０６によるレーザ光の照射を開始する。レーザ照射器８０６によって照射されたレーザ光は、電離層で反射して、第２カメラ８０４まで到達する。
【０２００】
ステップＳ１００４では、第１カメラ８０２と第２カメラ８０４を同期させて撮像する。第１カメラ８０２では、天体を撮像したデータが取得される。第２カメラ８０４では、レーザ照射器８０６から照射されて電離層によって反射されたレーザ光を撮像したデータが取得される。第１カメラで撮像されたデータは天体画像データとして制御装置８０８に記録される。第２カメラ８０４で撮像されたデータは点像画像データとして制御装置８０８に記録される。
【０２０１】
天体画像データと点像画像データは、いずれもＲＧＢそれぞれの色成分の照度の分布を表現している。本実施例では、ＲＧＢそれぞれの色成分に関して、ステップＳ１００６とステップＳ１００８の処理を実施する。
【０２０２】
ステップＳ１００６では、点像画像データから、ＰＳＦの分布ｈ（ｘ、ｙ）を取得する。本実施例の場合、第２カメラ８０４で撮像される画像は、地上と電離層の間を往復する間にシーイングの影響を受けたレーザ光の画像である。従って、天体を撮像する際のＰＳＦの分布としては、点像画像データから得られるＰＳＦの分布を用いる。
【０２０３】
ステップＳ１００８では、天体画像データにおける照度の分布と、ステップＳ１００６で取得されたＰＳＦの分布ｈ（ｘ、ｙ）を用いて、画像復元処理を行う。ステップＳ１００６で、すでに天体を撮像する際のＰＳＦの分布が取得されているので、実施例１や、実施例１の変形例１や、実施例１の変形例２で詳述した画像復元処理を施すことによって、シーイングの影響を排除した天体画像を得ることができる。
【０２０４】
本実施例の天体撮像装置８００によれば、大気のゆれによるシーイングの影響を排除した鮮明な天体画像を取得することができる。
【０２０５】
（実施例３）
図１１、図１２および図１３を参照しながら、実施例３の画像復元装置１１００について説明する。本実施例の画像復元装置１１００は、撮像時にぶれてしまった画像のみに基づいて、ぶれの影響を表現するＰＳＦの分布を推定して、ぶれのない理想的な画像を復元する。以下では、撮像時にぶれてしまった画像を劣化画像と呼び、ぶれのない理想的な画像を原画像と呼ぶ。
【０２０６】
図１１は画像復元装置１１００のハードウェア構成を示している。画像復元装置１１００は、後述する処理を実施するプログラムがインストールされた汎用のコンピュータである。画像復元装置１１００は、制御部１１０２と、入力部１１０４と、出力部１１０６と、演算部１１０８と、記憶部１１１０を備えている。記憶部１１１０は、プログラム記憶部１１１２と、データ記憶部１１１４を備えている。
画像復元装置１１００では、入力部１１０４を介して劣化画像のデータが入力されて、データ記憶部１１１４に記憶される。その後、制御部１１０２はプログラム記憶部１１１２に記憶されたプログラムに従って演算部１１０８で処理を行い、原画像を復元して原画像のデータをデータ記憶部１１１４に記憶する。データ記憶部１１１４に記憶された原画像のデータは、出力部１１０６を介して出力される。
なお画像復元装置１１００は、上記のような汎用のコンピュータの代わりに、同様の処理を実施可能な専用ロジックやＩＣ、ＦＰＧＡなどであってもよい。
【０２０７】
図１２のフローチャートを参照しながら、本実施例の画像復元装置１１００が行う処理の詳細を説明する。なお本実施例の画像復元装置１１００では、データ記憶部１１１４に記憶された劣化画像データから、ＲＧＢそれぞれの色成分の照度の分布ｇｒ（ｘ、ｙ）、ｇｇ（ｘ、ｙ）、ｇｂ（ｘ、ｙ）を特定しておいて、ＲＧＢそれぞれの色成分について図１２の処理を実行する。
【０２０８】
ステップＳ１２０２では、劣化画像の分布ｇ（ｘ、ｙ）を特定する。劣化画像の分布ｇ（ｘ、ｙ）は、光の照度の分布を位相特性も含めて記述する複素関数である。本実施例では、ｇ（ｘ、ｙ）の実部を劣化画像データにおける照度の分布と一致させ、ｇ（ｘ、ｙ）の虚部を全て０とする。
【０２０９】
ステップＳ１２０４では、劣化画像の分布ｇ（ｘ、ｙ）をフーリエ変換して、劣化画像のスペクトル分布Ｇ（ｓ、ｔ）を算出する。
【０２１０】
ステップＳ１２０６では、劣化画像のスペクトル分布Ｇ（ｓ、ｔ）の複素対数関数であるｌｎ（Ｇ（ｓ、ｔ））を計算して、関数Ｇｃ（ｓ、ｔ）に設定する。
【０２１１】
ステップＳ１２０８では、関数Ｇｃ（ｓ、ｔ）を逆フーリエ変換して、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）を算出する。
【０２１２】
図１３は図１２のステップＳ１２０８で算出される劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）の概略形状を示している。図１３に示すように、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）は、原画像の複素ケプストラムｆｃ（ｘ、ｙ）に、ぶれの影響を表現するＰＳＦの複素ケプストラムであるパルス列が重ね合わさった形状となっている。
【０２１３】
ステップＳ１２１０では、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）からパルス列を抽出する。本実施例では、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）の振幅値が所定のしきい値を超える個所をパルスとして識別する。
【０２１４】
ステップＳ１２１２では、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）から抽出されたパルス列を劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）から減じて、原画像の推定複素ケプストラムｆｃ（ｘ、ｙ）を算出する。
【０２１５】
ステップＳ１２１４では、原画像の推定複素ケプストラムｆｃ（ｘ、ｙ）をフーリエ変換して、関数Ｆｃ（ｓ、ｔ）を算出する。
【０２１６】
ステップＳ１２１６では、関数Ｆｃ（ｓ、ｔ）の複素指数関数であるｅｘｐ（Ｆｃ（ｓ、ｔ））を計算して、原画像の推定スペクトル分布Ｆ（ｓ、ｔ）を算出する。
【０２１７】
ステップＳ１２１８では、原画像の推定スペクトル分布Ｆ（ｓ、ｔ）を逆フーリエ変換して、原画像の推定分布ｆ（ｘ、ｙ）を算出し、原画像の推定分布ｆ（ｘ、ｙ）を原画像の分布として出力する。
【０２１８】
以上の処理をＲＧＢの色成分それぞれについて行うことで、ぶれのない画像を取得することができる。ぶれのない理想的な画像は、例えばディスプレイに表示してもよいし、画像データとしてハードディスクなどの記憶装置に記憶してもよいし、通信回線を経由して他のコンピュータに送信してもよい。
【０２１９】
なお図１４に本実施例のぶれ画像復元装置１１００が実現する機能のブロック図１４００を示している。データ入力ブロック１４０２から入力された劣化画像の分布ｇ（ｘ、ｙ）について、フーリエ変換ブロック１４０４の処理によって劣化画像のスペクトル分布Ｇ（ｓ、ｔ）が計算される。複素対数ブロック１４０６の処理によって関数Ｇｃ（ｓ、ｔ）が計算される。逆フーリエ変換ブロック１４０８の処理によって劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）が計算される。パルス抽出ブロック１４１０の処理によって劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）からパルス列が抽出される。減算ブロック１４１２によって劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）からパルス列が除去されて、原画像の推定複素ケプストラムｆｃ（ｘ、ｙ）が計算される。フーリエ変換ブロック１４１４の処理によって関数Ｆｃ（ｓ、ｔ）が計算される。複素指数ブロック１４１６の処理によって原画像の推定スペクトル分布Ｆ（ｓ、ｔ）が計算される。逆フーリエ変換ブロック１４１８の処理によって原画像の推定分布ｆ（ｘ、ｙ）が計算されて、データ出力ブロック１４２０に出力される。
【０２２０】
本実施例の画像復元装置１１００によれば、撮像時にぶれてしまった画像のみに基づいて、ぶれのない理想的な画像を復元することができる。
【０２２１】
（実施例４）
図面を参照しながら、本実施例の画像復元装置１１５０について説明する。本実施例の画像復元装置１１５０は、図１１に示す実施例３の画像復元装置１１００と同様のハードウェア構成を備えているが、プログラム記憶部１１１２にインストールされているプログラムの内容が異なる。
【０２２２】
図１５のフローチャートを参照しながら、本実施例の画像復元装置１１５０が行う処理の詳細を説明する。なお本実施例の画像復元装置１１５０では、データ記憶部１１１４に記憶された劣化画像データから、ＲＧＢそれぞれの色成分の照度の分布ｇｒ（ｘ、ｙ）、ｇｇ（ｘ、ｙ）、ｇｂ（ｘ、ｙ）を特定しておいて、ＲＧＢそれぞれの色成分について図１５の処理を実行する。
【０２２３】
ステップＳ１５０２では、劣化画像の分布ｇ（ｘ、ｙ）を特定する。劣化画像の分布ｇ（ｘ、ｙ）は、光の照度の分布を位相特性も含めて記述する複素関数である。本実施例では、ｇ（ｘ、ｙ）の実部を劣化画像データにおける照度の分布と一致させ、ｇ（ｘ、ｙ）の虚部を全て０とする。
【０２２４】
ステップＳ１５０４では、劣化画像の分布ｇ（ｘ、ｙ）をフーリエ変換して、劣化画像のスペクトル分布Ｇ（ｓ、ｔ）を算出する。
【０２２５】
ステップＳ１５０６では、劣化画像のスペクトル分布Ｇ（ｓ、ｔ）の複素対数関数であるｌｎ（Ｇ（ｓ、ｔ））を計算して、関数Ｇｃ（ｓ、ｔ）に設定する。
【０２２６】
ステップＳ１５０８では、関数Ｇｃ（ｓ、ｔ）を逆フーリエ変換して、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）を算出する。
【０２２７】
ステップＳ１５１０では、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）からパルス列を抽出する。本実施例では、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）の振幅値が所定のしきい値を超える個所をパルスとして識別する。
【０２２８】
ステップＳ１５１２では、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）から抽出されたパルス列より、ＰＳＦの分布ｈ（ｘ、ｙ）を特定する。画像の劣化がぶれによるものの場合には、劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）に現れるパルス列から、ＰＳＦの分布ｈ（ｘ、ｙ）を直接的に求めることができる。
【０２２９】
ステップＳ１５１４では、ステップＳ１５０２で特定された劣化画像の分布ｇ（ｘ、ｙ）と、ステップＳ１５１２で特定されたＰＳＦの分布ｈ（ｘ、ｙ）を用いて、原画像の分布ｆ（ｘ、ｙ）を復元する。この画像の復元処理は、種々の手法によって行うことができる。
【０２３０】
例えば実施例１で詳述したように、ＰＳＦ分布ｈ（ｘ、ｙ）をフーリエ変換してＯＴＦ分布Ｈ（ｓ、ｔ）を取得し、劣化画像のスペクトル分布Ｇ（ｓ、ｔ）をＯＴＦ分布Ｈ（ｓ、ｔ）で除して原画像のスペクトル分布Ｆ（ｓ、ｔ）を算出し、原画像のスペクトル分布Ｆ（ｓ、ｔ）を逆フーリエ変換することによって、原画像の分布ｆ（ｘ、ｙ）を取得することができる。
【０２３１】
あるいは実施例１の変形例１で詳述したように、ＰＳＦ分布ｈ（ｘ、ｙ）と劣化画像の分布ｇ（ｘ、ｙ）を用いて、図４に示す画像復元処理を行って、原画像の分布ｆ（ｘ、ｙ）を取得することもできる。
【０２３２】
あるいは実施例１の変形例２のように、ＰＳＦ分布ｈ（ｘ、ｙ）を真のＰＳＦ分布の最初の推定分布ｈ０（ｘ、ｙ）として、劣化画像の分布ｇ（ｘ、ｙ）を用いて、図５、図６および図７に示す画像復元処理を行って、原画像の分布ｆ（ｘ、ｙ）を取得することもできる。
【０２３３】
なお図１６は本実施例のぶれ画像復元装置１１５０が実現する機能のブロック図１６００を示している。データ入力ブロック１６０２から入力された劣化画像の分布ｇ（ｘ、ｙ）について、フーリエ変換ブロック１６０４の処理によって劣化画像のスペクトル分布Ｇ（ｓ、ｔ）が計算される。複素対数ブロック１６０６の処理によって関数Ｇｃ（ｘ、ｙ）が計算される。逆フーリエ変換ブロック１６０８の処理によって劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）が計算される。パルス抽出ブロック１６１０の処理によって劣化画像の複素ケプストラムｇｃ（ｘ、ｙ）からパルス列が抽出される。ＰＳＦ特定ブロック１６１２の処理によってＰＳＦ分布ｈ（ｘ、ｙ）が計算される。画像復元ブロック１６１４の処理によって原画像の分布ｆ（ｘ、ｙ）が計算されて、データ出力ブロック１６１６に出力される。
【０２３４】
本実施例の画像復元装置１１５０によれば、撮像時にぶれてしまった画像のみに基づいて、ぶれのない理想的な画像を復元することができる。
【０２３５】
以上、本発明の実施形態について詳細に説明したが、これらは例示に過ぎず、特許請求の範囲を限定するものではない。特許請求の範囲に記載の技術には、以上に例示した具体例を様々に変形、変更したものが含まれる。例えば、劣化復元の基本式である数式（１８）の代わりに、数式（１７）を反復計算に直接用いることも可能である。この場合にはフーリエ変換という過程は不要である。
また、本明細書または図面に説明した技術要素は、単独であるいは各種の組み合わせによって技術的有用性を発揮するものであり、出願時請求項記載の組み合わせに限定されるものではない。また、本明細書または図面に例示した技術は複数目的を同時に達成するものであり、そのうちの一つの目的を達成すること自体で技術的有用性を持つものである。
【符号の説明】
【０２３６】
１００：撮像装置
１０２：第２イメージセンサ
１０４：第１イメージセンサ
１０６：センサ支持台
１０８：レーザ光源
１１０：撮像台
１１２：制御装置
１５０：撮像装置
１６０：撮像装置
３０２：被写体画像
３０４：点像画像
３０６：ぶれのない被写体画像
８００：天体撮像装置
８０２：第１カメラ
８０４：第２カメラ
８０６：レーザ照射器
８０８：制御装置
９０２：天体画像
９０４：シーイングの影響を排除した天体画像
１１００：画像復元装置
１１０２：制御部
１１０４：入力部
１１０６：出力部
１１０８：演算部
１１１０：記憶部
１１１２：プログラム記憶部
１１１４：データ記憶部
１１５０：ぶれ画像復元装置
１４００：ブロック図
１４０２：データ入力ブロック
１４０４：フーリエ変換ブロック
１４０６：複素対数ブロック
１４０８：逆フーリエ変換ブロック
１４１０：パルス抽出ブロック
１４１２：減算ブロック
１４１４：フーリエ変換ブロック
１４１６：複素指数ブロック
１４１８：逆フーリエ変換ブロック
１４２０：データ出力ブロック
１６００：ブロック図
１６０２：データ入力ブロック
１６０４：フーリエ変換ブロック
１６０６：複素対数ブロック
１６０８：逆フーリエ変換ブロック
１６１０：パルス抽出ブロック
１６１２：ＰＳＦ特定ブロック
１６１４：画像復元ブロック
１６１６：データ出力ブロック
【特許請求の範囲】
【請求項１】
被写体を撮像して被写体画像データを出力する第１撮像手段と、
点光源を撮像して点像画像データを出力する第２撮像手段と、
被写体画像データと点像画像データに基いて改善された被写体画像データを生成する画像処理手段を備えており、
第１撮像手段と第２撮像手段は、共通する支持台によって支持されており、互いに同期して撮像するものであり、
被写体と点光源は、共通する載置台上に載置されており、
画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と被写体画像データに基いて改善された被写体画像データを生成する、撮像装置。
【請求項２】
前記画像処理手段が、点像分布関数と被写体画像データに基づいて改善された被写体画像データを生成する際に、
被写体画像データから劣化画像の分布を特定する処理と、
点像分布関数から光学的伝達関数を特定する処理と、
原画像の最初の推定分布を特定する処理と、
（１）原画像の推定分布をフーリエ変換して第１の関数を得る処理と、
（２）前記第１の関数に前記光学的伝達関数を乗じて第２の関数を得る処理と、
（３）前記第２の関数を逆フーリエ変換して第３の関数を得る処理と、
（４）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る処理と、
（５）前記第４の関数をフーリエ変換して第５の関数を得る処理と、
（６）前記第５の関数に前記光学的伝達関数の反転関数を乗じて第６の関数を得る処理と、
（７）前記第６の関数を逆フーリエ変換して第７の関数を得る処理と、
（８）原画像の推定分布に前記第７の関数を乗じて、原画像の次の推定分布を得る処理と、
原画像の次の推定分布を原画像の推定分布に置き換えて、前記（１）から（８）の処理を繰返す処理と、
原画像の推定分布に基づいて、改善された被写体画像データを生成する処理を実行することを特徴とする、請求項１の撮像装置。
【請求項３】
前記画像処理手段が、点像分布関数と被写体画像データに基づいて改善された被写体画像データを生成する際に、
被写体画像データから劣化画像の分布を特定する処理と、
原画像の最初の推定分布を特定する処理と、
点像画像データに基いて取得された点像分布関数を、点像分布関数の最初の推定分布に設定する処理と、
（Ａ）前記劣化画像の分布と、原画像の推定分布と、点像分布関数の推定分布に基づいて、点像分布関数の推定分布を更新する処理と、
（Ｂ）前記劣化画像の分布と、原画像の推定分布と、点像分布関数の推定分布に基づいて、原画像の推定分布を更新する処理と、
上記（Ａ）と（Ｂ）の工程を交互に繰返す処理と、
原画像の推定分布に基づいて、改善された被写体画像データを生成する処理を実行するものであり、
前記（Ａ）の点像分布関数の推定分布を更新する処理は、
（Ａ１）原画像の推定分布をフーリエ変換して、原画像の推定分布のスペクトル分布を得る処理と、
（Ａ２）点像分布関数の推定分布をフーリエ変換して第１の関数を得る処理と、
（Ａ３）前記第１の関数に前記原画像の推定分布のスペクトル分布を乗じて第２の関数を得る処理と、
（Ａ４）前記第２の関数を逆フーリエ変換して第３の関数を得る処理と、
（Ａ５）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る処理と、
（Ａ６）前記第４の関数をフーリエ変換して第５の関数を得る処理と、
（Ａ７）前記第５の関数に前記原画像の推定分布のスペクトル分布の反転関数を乗じて第６の関数を得る処理と、
（Ａ８）前記第６の関数を逆フーリエ変換して第７の関数を得る処理と、
（Ａ９）前記点像分布関数の推定分布に前記第７の関数を乗じて、点像分布関数の次の推定分布を得る処理と、
（Ａ１０）点像分布関数の次の推定分布を点像分布関数の推定分布に置き換える処理とを備え、
前記（Ｂ）の原画像の推定分布を更新する処理は、
（Ｂ１）点像分布関数の推定分布をフーリエ変換して、光学的伝達関数の推定分布を得る処理と、
（Ｂ２）原画像の推定分布をフーリエ変換して第１の関数を得る処理と、
（Ｂ３）前記第１の関数に前記光学的伝達関数の推定分布を乗じて第２の関数を得る処理と、
（Ｂ４）前記第２の関数を逆フーリエ変換して第３の関数を得る処理と、
（Ｂ５）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る処理と、
（Ｂ６）前記第４の関数をフーリエ変換して第５の関数を得る処理と、
（Ｂ７）前記第５の関数に前記光学的伝達関数の推定分布の反転関数を乗じて第６の関数を得る処理と、
（Ｂ８）前記第６の関数を逆フーリエ変換して第７の関数を得る処理と、
（Ｂ９）原画像の推定分布に前記第７の関数を乗じて、原画像の次の推定分布を得る処理と、
（Ｂ１０）原画像の次の推定分布を原画像の推定分布に置き換える処理とを備えることを特徴とする、請求項１の撮像装置。
【請求項４】
天体を撮像して天体画像データを出力する第１撮像手段と、
電離層に向けてレーザ光を照射するレーザ光照射手段と、
レーザ光照射手段から照射されて電離層で反射したレーザ光を撮像して点像画像データを出力する第２撮像手段と、
天体画像データと点像画像データに基いて改善された天体画像データを生成する画像処理手段を備えており、
第１撮像手段と第２撮像手段は、互いに同期して撮像するものであり、
画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と天体画像データに基いて改善された天体画像データを生成する、天体撮像装置。
【請求項５】
ぶれ画像からぶれのない画像を復元する装置であって、
ぶれ画像の分布を特定する手段と、
ぶれ画像の分布についての複素ケプストラムを算出する手段と、
ぶれ画像の分布についての複素ケプストラムからパルス列を抽出する手段と、
抽出されたパルス列をぶれ画像の分布についての複素ケプストラムから除去して、ぶれのない画像の分布についての複素ケプストラムを推定する手段と、
ぶれのない画像の分布についての複素ケプストラムからぶれのない画像の分布を算出する手段を備える、ぶれ画像復元装置。
【請求項６】
ぶれ画像からぶれのない画像を復元する装置であって、
ぶれ画像の分布を特定する手段と、
ぶれ画像の分布についての複素ケプストラムを算出する手段と、
ぶれ画像の分布についての複素ケプストラムからパルス列を抽出する手段と、
抽出されたパルス列に基づいて点像分布関数を取得する手段と、
ぶれ画像の分布と点像分布関数に基づいてぶれのない画像の分布を推定する手段を備える、ぶれ画像復元装置。
【請求項７】
第１撮像手段によって被写体を撮像して被写体画像データを出力する工程と、
第２撮像手段によって点光源を撮像して点像画像データを出力する工程と、
画像処理手段によって被写体画像データと点像画像データに基いて改善された被写体画像データを生成する工程を備えており、
第１撮像手段と第２撮像手段は、共通する支持台によって支持されており、互いに同期して撮像するものであり、
被写体と点光源は、共通する載置台上に載置されており、
画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と被写体画像データに基いて改善された被写体画像データを生成する、撮像方法。
【請求項８】
第１撮像手段によって天体を撮像して天体画像データを出力する工程と、
レーザ光照射手段によって電離層に向けてレーザ光を照射する工程と、
レーザ光照射手段から照射されて電離層で反射したレーザ光を第２撮像手段によって撮像して点像画像データを出力する工程と、
画像処理手段によって天体画像データと点像画像データに基いて改善された天体画像データを生成する工程を備えており、
第１撮像手段と第２撮像手段は、互いに同期して撮像するものであり、
画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と天体画像データに基いて改善された天体画像データを生成する、天体撮像方法。
【請求項９】
ぶれ画像からぶれのない画像を復元する方法であって、
ぶれ画像の分布を特定する工程と、
ぶれ画像の分布についての複素ケプストラムを算出する工程と、
ぶれ画像の分布についての複素ケプストラムからパルス列を抽出する工程と、
抽出されたパルス列をぶれ画像の分布についての複素ケプストラムから除去して、ぶれのない画像の分布についての複素ケプストラムを推定する工程と、
ぶれのない画像の分布についての複素ケプストラムからぶれのない画像の分布を算出する工程を備える、ぶれ画像復元方法。
【請求項１０】
ぶれ画像からぶれのない画像を復元する方法であって、
ぶれ画像の分布を特定する工程と、
ぶれ画像の分布についての複素ケプストラムを算出する工程と、
ぶれ画像の分布についての複素ケプストラムからパルス列を抽出する工程と、
抽出されたパルス列に基づいて点像分布関数を取得する工程と、
ぶれ画像の分布と点像分布関数に基づいてぶれのない画像の分布を推定する工程を備える、ぶれ画像復元方法。
【請求項１１】
請求項７の方法における各工程をコンピュータに実行させるためのプログラム。
【請求項１２】
請求項８の方法における各工程をコンピュータに実行させるためのプログラム。
【請求項１３】
請求項９の方法における各工程をコンピュータに実行させるためのプログラム。
【請求項１４】
請求項１０の方法における各工程をコンピュータに実行させるためのプログラム。
【請求項１】
被写体を撮像して被写体画像データを出力する第１撮像手段と、
点光源を撮像して点像画像データを出力する第２撮像手段と、
被写体画像データと点像画像データに基いて改善された被写体画像データを生成する画像処理手段を備えており、
第１撮像手段と第２撮像手段は、共通する支持台によって支持されており、互いに同期して撮像するものであり、
被写体と点光源は、共通する載置台上に載置されており、
画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と被写体画像データに基いて改善された被写体画像データを生成する、撮像装置。
【請求項２】
前記画像処理手段が、点像分布関数と被写体画像データに基づいて改善された被写体画像データを生成する際に、
被写体画像データから劣化画像の分布を特定する処理と、
点像分布関数から光学的伝達関数を特定する処理と、
原画像の最初の推定分布を特定する処理と、
（１）原画像の推定分布をフーリエ変換して第１の関数を得る処理と、
（２）前記第１の関数に前記光学的伝達関数を乗じて第２の関数を得る処理と、
（３）前記第２の関数を逆フーリエ変換して第３の関数を得る処理と、
（４）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る処理と、
（５）前記第４の関数をフーリエ変換して第５の関数を得る処理と、
（６）前記第５の関数に前記光学的伝達関数の反転関数を乗じて第６の関数を得る処理と、
（７）前記第６の関数を逆フーリエ変換して第７の関数を得る処理と、
（８）原画像の推定分布に前記第７の関数を乗じて、原画像の次の推定分布を得る処理と、
原画像の次の推定分布を原画像の推定分布に置き換えて、前記（１）から（８）の処理を繰返す処理と、
原画像の推定分布に基づいて、改善された被写体画像データを生成する処理を実行することを特徴とする、請求項１の撮像装置。
【請求項３】
前記画像処理手段が、点像分布関数と被写体画像データに基づいて改善された被写体画像データを生成する際に、
被写体画像データから劣化画像の分布を特定する処理と、
原画像の最初の推定分布を特定する処理と、
点像画像データに基いて取得された点像分布関数を、点像分布関数の最初の推定分布に設定する処理と、
（Ａ）前記劣化画像の分布と、原画像の推定分布と、点像分布関数の推定分布に基づいて、点像分布関数の推定分布を更新する処理と、
（Ｂ）前記劣化画像の分布と、原画像の推定分布と、点像分布関数の推定分布に基づいて、原画像の推定分布を更新する処理と、
上記（Ａ）と（Ｂ）の工程を交互に繰返す処理と、
原画像の推定分布に基づいて、改善された被写体画像データを生成する処理を実行するものであり、
前記（Ａ）の点像分布関数の推定分布を更新する処理は、
（Ａ１）原画像の推定分布をフーリエ変換して、原画像の推定分布のスペクトル分布を得る処理と、
（Ａ２）点像分布関数の推定分布をフーリエ変換して第１の関数を得る処理と、
（Ａ３）前記第１の関数に前記原画像の推定分布のスペクトル分布を乗じて第２の関数を得る処理と、
（Ａ４）前記第２の関数を逆フーリエ変換して第３の関数を得る処理と、
（Ａ５）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る処理と、
（Ａ６）前記第４の関数をフーリエ変換して第５の関数を得る処理と、
（Ａ７）前記第５の関数に前記原画像の推定分布のスペクトル分布の反転関数を乗じて第６の関数を得る処理と、
（Ａ８）前記第６の関数を逆フーリエ変換して第７の関数を得る処理と、
（Ａ９）前記点像分布関数の推定分布に前記第７の関数を乗じて、点像分布関数の次の推定分布を得る処理と、
（Ａ１０）点像分布関数の次の推定分布を点像分布関数の推定分布に置き換える処理とを備え、
前記（Ｂ）の原画像の推定分布を更新する処理は、
（Ｂ１）点像分布関数の推定分布をフーリエ変換して、光学的伝達関数の推定分布を得る処理と、
（Ｂ２）原画像の推定分布をフーリエ変換して第１の関数を得る処理と、
（Ｂ３）前記第１の関数に前記光学的伝達関数の推定分布を乗じて第２の関数を得る処理と、
（Ｂ４）前記第２の関数を逆フーリエ変換して第３の関数を得る処理と、
（Ｂ５）前記劣化画像の分布を前記第３の関数で除して第４の関数を得る処理と、
（Ｂ６）前記第４の関数をフーリエ変換して第５の関数を得る処理と、
（Ｂ７）前記第５の関数に前記光学的伝達関数の推定分布の反転関数を乗じて第６の関数を得る処理と、
（Ｂ８）前記第６の関数を逆フーリエ変換して第７の関数を得る処理と、
（Ｂ９）原画像の推定分布に前記第７の関数を乗じて、原画像の次の推定分布を得る処理と、
（Ｂ１０）原画像の次の推定分布を原画像の推定分布に置き換える処理とを備えることを特徴とする、請求項１の撮像装置。
【請求項４】
天体を撮像して天体画像データを出力する第１撮像手段と、
電離層に向けてレーザ光を照射するレーザ光照射手段と、
レーザ光照射手段から照射されて電離層で反射したレーザ光を撮像して点像画像データを出力する第２撮像手段と、
天体画像データと点像画像データに基いて改善された天体画像データを生成する画像処理手段を備えており、
第１撮像手段と第２撮像手段は、互いに同期して撮像するものであり、
画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と天体画像データに基いて改善された天体画像データを生成する、天体撮像装置。
【請求項５】
ぶれ画像からぶれのない画像を復元する装置であって、
ぶれ画像の分布を特定する手段と、
ぶれ画像の分布についての複素ケプストラムを算出する手段と、
ぶれ画像の分布についての複素ケプストラムからパルス列を抽出する手段と、
抽出されたパルス列をぶれ画像の分布についての複素ケプストラムから除去して、ぶれのない画像の分布についての複素ケプストラムを推定する手段と、
ぶれのない画像の分布についての複素ケプストラムからぶれのない画像の分布を算出する手段を備える、ぶれ画像復元装置。
【請求項６】
ぶれ画像からぶれのない画像を復元する装置であって、
ぶれ画像の分布を特定する手段と、
ぶれ画像の分布についての複素ケプストラムを算出する手段と、
ぶれ画像の分布についての複素ケプストラムからパルス列を抽出する手段と、
抽出されたパルス列に基づいて点像分布関数を取得する手段と、
ぶれ画像の分布と点像分布関数に基づいてぶれのない画像の分布を推定する手段を備える、ぶれ画像復元装置。
【請求項７】
第１撮像手段によって被写体を撮像して被写体画像データを出力する工程と、
第２撮像手段によって点光源を撮像して点像画像データを出力する工程と、
画像処理手段によって被写体画像データと点像画像データに基いて改善された被写体画像データを生成する工程を備えており、
第１撮像手段と第２撮像手段は、共通する支持台によって支持されており、互いに同期して撮像するものであり、
被写体と点光源は、共通する載置台上に載置されており、
画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と被写体画像データに基いて改善された被写体画像データを生成する、撮像方法。
【請求項８】
第１撮像手段によって天体を撮像して天体画像データを出力する工程と、
レーザ光照射手段によって電離層に向けてレーザ光を照射する工程と、
レーザ光照射手段から照射されて電離層で反射したレーザ光を第２撮像手段によって撮像して点像画像データを出力する工程と、
画像処理手段によって天体画像データと点像画像データに基いて改善された天体画像データを生成する工程を備えており、
第１撮像手段と第２撮像手段は、互いに同期して撮像するものであり、
画像処理手段は、点像画像データに基づいて点像分布関数を取得し、点像分布関数と天体画像データに基いて改善された天体画像データを生成する、天体撮像方法。
【請求項９】
ぶれ画像からぶれのない画像を復元する方法であって、
ぶれ画像の分布を特定する工程と、
ぶれ画像の分布についての複素ケプストラムを算出する工程と、
ぶれ画像の分布についての複素ケプストラムからパルス列を抽出する工程と、
抽出されたパルス列をぶれ画像の分布についての複素ケプストラムから除去して、ぶれのない画像の分布についての複素ケプストラムを推定する工程と、
ぶれのない画像の分布についての複素ケプストラムからぶれのない画像の分布を算出する工程を備える、ぶれ画像復元方法。
【請求項１０】
ぶれ画像からぶれのない画像を復元する方法であって、
ぶれ画像の分布を特定する工程と、
ぶれ画像の分布についての複素ケプストラムを算出する工程と、
ぶれ画像の分布についての複素ケプストラムからパルス列を抽出する工程と、
抽出されたパルス列に基づいて点像分布関数を取得する工程と、
ぶれ画像の分布と点像分布関数に基づいてぶれのない画像の分布を推定する工程を備える、ぶれ画像復元方法。
【請求項１１】
請求項７の方法における各工程をコンピュータに実行させるためのプログラム。
【請求項１２】
請求項８の方法における各工程をコンピュータに実行させるためのプログラム。
【請求項１３】
請求項９の方法における各工程をコンピュータに実行させるためのプログラム。
【請求項１４】
請求項１０の方法における各工程をコンピュータに実行させるためのプログラム。
【図１】
【図２】
【図４】
【図５】
【図６】
【図７】
【図８】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図３】
【図９】
【図２】
【図４】
【図５】
【図６】
【図７】
【図８】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図３】
【図９】
【公開番号】特開２０１２−１７８１８３（Ｐ２０１２−１７８１８３Ａ）
【公開日】平成２４年９月１３日（２０１２．９．１３）
【国際特許分類】
【出願番号】特願２０１２−１２７３６９（Ｐ２０１２−１２７３６９）
【出願日】平成２４年６月４日（２０１２．６．４）
【分割の表示】特願２００８−２６３８（Ｐ２００８−２６３８）の分割
【原出願日】平成２０年１月９日（２００８．１．９）
【出願人】（５９１０５６２２６）ライトロン株式会社 (3)
【Ｆターム（参考）】
【公開日】平成２４年９月１３日（２０１２．９．１３）
【国際特許分類】
【出願日】平成２４年６月４日（２０１２．６．４）
【分割の表示】特願２００８−２６３８（Ｐ２００８−２６３８）の分割
【原出願日】平成２０年１月９日（２００８．１．９）
【出願人】（５９１０５６２２６）ライトロン株式会社 (3)
【Ｆターム（参考）】
[ Back to top ]