画像復元方法
【課題】本発明は、十分に復元された復元像を得ることができる画像復元方法を提供する。
【解決手段】ステップ１０１〜１０４でｎ次元の理想像モデルの初期値を設定すると共に、ｎ次元のＰＳＦモデルの初期値として周波数空間でデルタ関数となる分布を設定し、ステップ１０５で制限条件と終了条件を与えた後、ステップ１０６〜１１４で理想像モデルを実空間あるいは周波数空間で像とＰＳＦモデルにより修正し、実空間の理想像モデルに対して実空間の制限を加え、ＰＳＦモデルを実空間あるいは周波数空間で像と理想像モデルにより修正し、周波数空間のＰＳＦモデルに対して周波数空間の制限を加えた後実空間に変換し、実空間のＰＳＦモデルに対して実空間の制限を加える処理を、ステップ１１５で所定の終了条件を満たすまで繰返し実行し、最終的な理想像モデルを復元像として生成する。
【解決手段】ステップ１０１〜１０４でｎ次元の理想像モデルの初期値を設定すると共に、ｎ次元のＰＳＦモデルの初期値として周波数空間でデルタ関数となる分布を設定し、ステップ１０５で制限条件と終了条件を与えた後、ステップ１０６〜１１４で理想像モデルを実空間あるいは周波数空間で像とＰＳＦモデルにより修正し、実空間の理想像モデルに対して実空間の制限を加え、ＰＳＦモデルを実空間あるいは周波数空間で像と理想像モデルにより修正し、周波数空間のＰＳＦモデルに対して周波数空間の制限を加えた後実空間に変換し、実空間のＰＳＦモデルに対して実空間の制限を加える処理を、ステップ１１５で所定の終了条件を満たすまで繰返し実行し、最終的な理想像モデルを復元像として生成する。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】本発明は、顕微鏡あるいは共焦点顕微鏡などの光学機器により取得した画像の画質を改善する画像復元方法に関するものである。
【０００２】
【従来の技術】顕微鏡あるいは共焦点顕微鏡などの光学機器により物体を観察した場合、観察画像は、元の物体に比べて光学機器のボケ分だけ劣化することが知られている。つまり、このときの観察画像の強度分布ｇは、元の物体の輝度分布ｆと光学機器の点像分布関数（ＰＳＦ）ｈの畳込み（ｃｏｎｖｏｌｕｔｉｏｎ）にノイズｎが加わった、
【０００３】
【数１】
で表現することができる。
【０００４】そして、このようなノイズｎが加わった画像から数値計算を用いて光学機器によるボケを取り除き、理想的な像を得る技術は、像の「復元」あるいは「デコンボリューション（ｄｅｃｏｎｖｏｌｕｔｉｏｎ）」と呼ばれている。
【０００５】ところで、光学顕微鏡では、深さ方向を等間隔で変位させつつ、試料像を撮像し、３次元画像（＝積層画像）を得る光学切片法（ｏｐｔｉｃａｌ ｓｅｃｔｉｏｎｉｎｇ）が用いられている。ところが、光学顕微鏡のＰＳＦは水平方向に比べて深さ方向に大きく広がっているので、積層画像の中の各画像は試料の輝度分布を正確に反映した断面像にはならず、その上下から漏れてくるボケが重畳している。
【０００６】そこで、従来、光学顕微鏡のＰＳＦを基に積層画像からボケを取り除き、理想的な３次元像を「復元」する計算処理が行われてきた。ここで用いられるＰＳＦは、理論値あるいは測定値である。
【０００７】しかし、光学顕微鏡の場合、実際のＰＳＦと理論値は必ずしも一致せず、また正確なＰＳＦの実測も困難である。そこで、「ブラインドデコンボリューション」という像復元の手法が検討されている。この手法は、１つの画像から理想像とＰＳＦの両方を復元するもので、ＰＳＦが未知であっても理想像を復元できるという長所がある。
【０００８】ここで、ブラインドデコンボリューションは、交互法（ａｌｔｅｒｎａｔｉｖｅ ｍｅｔｈｏｄ）というアルゴリズムで実行される。つまり、このような交互法では、初めに、理想像モデルｆとＰＳＦモデルｈに初期値を与え、次に、２つのモデルの畳込みが測定された実際の像ｇになるように両者を交互に修正するとともに、これらの過程を繰り返す。この過程で、モデルには適当な制限条件が加えられ、この交互修正を最適化ループとして何度か繰り返すうちに、２つのモデルは最適化されていく。
【０００９】図４は、従来の交互法を、さらに説明するためのもので、まず、ステップ４０１で、光学機器で測定した実際の像ｇを読み込む。次に、４０２で、実際の像ｇを、実空間から周波数空間へ変換する。これは畳込みの計算を、周波数空間で積として効率よく計算するためである。変換は一般に高速フーリエ変換（ＦＦＴ）で行われる。以下、実空間のものは、小文字（ｇ，ｆ，ｈ）で記し、周波数のものは、大文字（Ｇ，Ｆ，Ｈ）で記す。
【００１０】次に、４０３で、理想像モデルＦの初期値として、測定された実際の像Ｇを設定し、ステップ４０４で、ＰＳＦモデルｈの初期値として、インパルス即ち実空間のＤｉｒａｃのデルタ関数を設定する。そして、ステップ４０５で、制限条件と終了条件を与える。この場合、制限条件は、ＰＳＦの空間周波数には上限、即ちバンドリミットがあるので、このバンドリミットより与えることが可能であるが、現実は、その値は不明なので像のスペクトルから推定して値を設定する。
【００１１】図５（ａ）（ｂ）（ｃ）は、バンドリミットの決め方を説明するもので、ここでは、物体、ＰＳＦおよび像のスペクトルの関係を簡単のために１次元で表している。物体のスペクトル５０１は通常、測定に使用する光学機器のＰＳＦのスペクトル５０２より高周波の成分を含んでいる。よって、測定された像のスペクトル５０４は、ＰＳＦのスペクトル５０２のバンドリミット５０３と同じバンドリミット５０５をもつと考えられるから、この時の像のスペクトル５０４から推定して設定する。なお、３次元像の場合、３軸のバンドリミットを設定しておく。このバンドリミットは後述するステップ４１２で実行される。また、終了条件は一般に、最適化ループの繰り返し回数である。これは後述のステップ４１５で判断される。
【００１２】図４に戻って、ステップ４０６で、ＰＳＦモデルｈを、実空間から周波数空間へ変換する。次いで、ステップ４０７で、ＰＳＦモデルＨとの積が実際の像Ｇになるように、理想像モデルＦを修正する。この場合、理想像モデルＦの修正には、これまでに以下のような更新式が提案されている。・差分
【００１３】
【数２】
・ウィーナフィルタ付き差分
【００１４】
【数３】
・逆フィルタ
【００１５】
【数４】
・ウィーナフィルタ
【００１６】
【数５】
【００１７】次に、ステップ４０８で、理想像モデルＦを、周波数空間から実空間へ変換する。これは、ステップ４０９で実空間の非負制限を理想像モデルｆに適用するためである。変換は一般に高速逆フーリエ変換（ＩＦＦＴ）で行われる。
【００１８】ステップ４０９で、理想像モデルｆに、非負制限を加える。つまり、負の要素値を０で置き換え、それ以外の要素値はそのまま残す。理想像は光の強度分布であるから負の値をとらないからである。さらに、理想像モデルｆの総和が変化しないように補正する例もある。
【００１９】その後、ステップ４１０〜４１４で、理想像モデルとＰＳＦモデルを入れ変えて上述したステップ４０６〜４０９と同様の処理をすることで、ＰＳＦモデルを修正する。
【００２０】なお、ステップ４１２でのＰＳＦモデルＨに対する周波数空間の制限条件は、バンドリミットである。つまり、各軸のバンドリミットより周波数の大きい成分を０で置き換え、それ以外の成分はそのまま残す。また、ステップ４１４での実空間のＰＳＦモデルｈに対する制限条件としては、非負制限の他に、要素値の和を１に規格化する例もある。これは、像形成におけるエネルギー保存則を成立させるためである。
【００２１】その後、ステップ４１５で、最適化ループであるステップ４０６〜４１４の繰り返し回数がステップ４０５で設定した回数に達したら、処理を終了し、ステップ４１６で、最適化ループであるステップ４０６〜４１４による修正を経た理想像モデルｆが、実際の像ｇからボケを除去した復元像として出力される。
【００２２】以上の例ではモデルを周波数空間で修正していたが、実空間で修正する方法もある。その場合、処理の順序が若干変わって図６に示すようになる。この場合、ステップ６１０とステップ６１１は、図４のステップ４１０とステップ４１１ｓの順序が反対になる点が異なるだけでステップ６０１〜６１５での個々の処理は、上述した図４と同じなので、ここでの説明は省略する。なお、このような実空間での修正では、以下のような更新式が提案されている。
・積（multiplicative method)
【００２３】
【数６】
・最尤法（maximum likelihood estimation)（ref.７）
【００２４】
【数７】
【００２５】
【発明が解決しようとする課題】ところが、このようなブラインドデコンボリューションは不適切（ｉｌｌ−ｐｏｓｅｄ）な逆問題で、解の唯一性は保証されていない。つまり、１つの実際の像になる理想像モデルとＰＳＦモデルの組み合わせが無数に存在し、最終的に復元される２つのモデルは、モデルの初期値にも影響される。
【００２６】このため、上述したステップ４０４に示すように、ＰＳＦモデルの初期値として実空間のデルタ関数を取ると、理想像モデルの復元が十分に進行しない傾向が見られた。
【００２７】これを図７を参照して説明すると、実空間のデルタ関数は、図７（ａ）に示すように、周波数空間ですべての成分が１に等しいフラットな分布７０１になる。ここからＰＳＦモデルを出発すると、最適化ループを経るごとにＰＳＦモデルは真のＰＳＦ７０３に対して上から（各周波数成分が大きい方から）近づいてゆく傾向がある。このため、最終的なＰＳＦモデル、即ち復元ＰＳＦ７０２の各周波数成分は、図７（ｂ）に示すように対応する真のＰＳＦ７０３より大きいまま収束することが多い。つまり、理想像モデルはＰＳＦモデルとの畳込みが実際の像になるように最適化されてゆくので、図７R>７（ｃ）に示すように、最終的な理想像モデル、即ち復元像７０４の各周波数成分は、対応する真の物体７０５より小さくなり、実空間で見れば、細部が十分に復元されていない復元像になるという問題があった。本発明は上記事情に鑑みてなされたもので、十分に復元された復元像を得ることができる画像復元方法を提供することを目的とする。
【００２８】
【課題を解決するための手段】請求項１記載の発明は、物体と点像分布関数（ＰＳＦ）の畳込みとして取得されたｎ次元の像から、ＰＳＦによるボケを除去したｎ次元の復元像を得るための画像復元方法において、ｎ次元の理想像モデルに初期値として前記ｎ次元の像を設定するとともに、ｎ次元のＰＳＦモデルに初期値として周波数空間でデルタ関数となる分布を設定した後、前記理想像モデルを実空間あるいは周波数空間の像とＰＳＦモデルとにより修正し、実空間の前記理想像モデルに対して実空間の制限を加え、前記ＰＳＦモデルを実空間あるいは周波数空間の像と理想像モデルとにより修正し、周波数空間の前記ＰＳＦモデルに対して周波数空間の制限を加えた後、実空間に変換し、実空間の前記ＰＳＦモデルに対して実空間の制限を加えるような処理を所定の終了条件を満たすまで繰り返し実行し、最終的な理想像モデルを復元像として生成するようにしている。
【００２９】請求項２記載の発明は、請求項１記載の発明において、前記理想像モデルの初期値として、実際の像を設定している。請求項３記載の発明は、請求項１記載の発明において、前記理想像モデルと前記ＰＳＦモデルに対する実空間の制限が非負制限である。
【００３０】この結果、本発明によれば、周波数空間でデルタ関数となる分布をＰＳＦモデルの初期値にしてブラインドデコンボリューションを実施するが、この場合、周波数空間のデルタ関数は、実空間でフラットな分布、つまり完全なボケに対応したものになるので、この周波数空間のデルタ関数からＰＳＦモデルを出発すると、最適化ループを経るごとにＰＳＦモデルは、真のＰＳＦに対して各周波数成分が小さい方から近づくようになり、対応する復元像の各周波数成分は、対応する真の物体より小さくならず、実空間で見ると、像の復元を細部まで十分に進行させることができる。
【００３１】
【発明の実施の形態】以下、本発明の実施の形態を図面に従い説明する。
（第１の実施の形態）この第１の実施の形態は、周波数空間でモデルを更新する例を示すもので、ここでは、落射蛍光共焦点レーザ走査顕微鏡（ＥＦＣＬＳＭ）の画像復元に適用した例を示している。
【００３２】図１は、本発明による画像復元方法の処理の流れを示しており、ここでは、ＥＦＣＬＳＭで光学切片法により測定した３次元画像（積層画像）ファイルを読み込み、それを画像処理し、ボケを除去した復元像を表示する、あるいは画像ファイルとして出力するプログラムである。ただし、画像処理の部分以外は従来の技術でも説明した一般的な内容であるから説明を省略する。
【００３３】まず、ステップ１０１で、３次元画像ファイルを読み込み、実際の像ｇの３次元複素数データ配列の実部に代入する。ここでは、虚部に０を代入する。つまり、後で高速フーリエ変換（ＦＦＴ）を行うため、ｘ，ｙ，ｚ各軸のデータ点数が２の累乗になっていないときは、足りないデータ点に０詰めして、点数を２の累乗にしておく。また、実際の像ｇの全要素値の和Σｇを求めておき、後のステップ１０９で使用する。
【００３４】次に、ステップ１０２で、実空間の実際の像ｇを、ＦＦＴで周波数空間の実際の像Ｇに変換する。次いで、１０３で、理想像モデルＦを、実際の像Ｇと同じ大きさの３次元複素数データ配列として用意し、その初期値として周波数空間の実際の像Ｇの値を代入し、ステップ１０４で、ＰＳＦモデルＨを、実際の像Ｇと同じ大きさの３次元複素数データ配列として用意し、その初期値として、周波数空間のデルタ関数を代入する。
【００３５】
【数８】
【００３６】ここで、実空間のＰＳＦモデルｈの全要素値の和は、自動的に１になる。次に、ステップ１０５で、制限条件と終了条件を与える。このうちの制限条件は、実際の像Ｇから推定して、ＰＳＦモデルＨに対する３軸のバンドリミットを設定する。また、終了条件として、最適化ループの繰り返し回数Ｎを設定する。
【００３７】ステップ１０６で、実空間のＰＳＦモデルｈをＦＦＴで周波数空間のＰＳＦモデルＨに変換する。ただし、初期値を周波数空間で与えているので、１回目は必要ない。次に、ステップ１０７で、修正されたＰＳＦモデルＨとの積が実際の像Ｇになるように、理想像モデルＦの各成分を次式の差分による更新式で修正する。
【００３８】
【数９】
【００３９】そして、ステップ１０８で、周波数空間の理想像モデルＦを、ＩＦＦＴで実空間の理想像モデルｆに変換し、ステップ１０９で、実空間の理想像モデルｆに非負制限を加える。この場合、まず理想像モデルｆの、負の実部を持つ要素値を０で置き換え、それ以外の要素値は実部を複素数の絶対値で置き換え、虚部を０にする。次に、理想像モデルｆの総和が変化しないように、ステップ１０１で求めたΣｇを使って次式のような規格化を行う。
【００４０】
【数１０】
【００４１】次に、ステップ１１０で、実空間の理想像モデルｆをＦＦＴで周波数空間の理想像モデルＦに変換する。次いで、ステップ１１１で、修正された理想像モデルＦとの積が実際の像Ｇになるように、ＰＳＦモデルＨの各成分を次式の差分による更新式で修正する。
【００４２】
【数１１】
【００４３】次に、ステップ１１２で、周波数空間のＰＳＦモデルＨにバンドリミットを加える。即ち、ステップ１０５で設定した各軸のバンドリミットより周波数の大きい成分を０で置き換え、それ以外の成分はそのまま残す。
【００４４】次に、ステップ１１３で、周波数空間のＰＳＦモデルＨを、ＩＦＦＴで実空間のＰＳＦモデルｈに変換し、ステップ１１４で、実空間のＰＳＦモデルｈに非負制限を加える。この場合、まずＰＳＦモデルｈの、負の実部を持つ要素値を０で置き換え、それ以外の要素値は実部を複素数の絶対値で置き換え、虚部を０にする。次に、ＰＳＦモデルｈの全要素値の和が１になるように、次式のような規格化を行う。
【００４５】
【数１２】
【００４６】その後、最適化ループであるステップ１０６〜１１４の繰り返し回数がステップ１０５で設定した回数Ｎに達したら、ステップ１１５で、処理を終了し、ステップ１１６で、最適化ループであるステップ１０６〜１１４により繰り返し修正された理想像モデルｆの実部が、実際の像ｇからボケを除去した復元像として出力され、復元の画像処理を終了する。
【００４７】従って、このようにすれば、周波数空間でデルタ関数となる分布をＰＳＦモデルの初期値にしたブラインドデコンボリューションが実施されるが、この場合、周波数空間のデルタ関数は、実空間でフラットな分布、つまり完全なボケに対応したものになるので、そこから出発したＰＳＦモデルがシャープになり過ぎることはない。これにより、図２（ａ）に示すように、ＰＳＦモデルを周波数空間のデルタ関数２０１から出発すると、最適化ループを経るごとにＰＳＦモデルは、図２（ｂ）に示すように真のＰＳＦ２０３に対して基本的に下から（各周波数成分が小さい方から）近づくようになり、復元ＰＳＦ２０２の各周波数成分も、対応する真のＰＳＦ２０３より小さい方向で収束する傾向になる。この結果、図２（ｃ）に示すように、対応する復元像２０４の各周波数成分は、対応する真の物体２０５より小さくならないので、実空間で見ると、像の復元を細部まで十分に進行させることができる。つまり、周波数空間でデルタ関数となる分布をＰＳＦモデルの初期値に設定したのち、周波数空間でモデルを更新するブラインドデコンボリューションを実施することにより、像の復元を細部まで進行させることができるようになり、細部まで十分に復元された復元像を生成できる。
【００４８】なお、本発明は、更新式や他の制限条件、終了条件などに関係なく、すべての交互法によるブラインドデコンボリューションに有効である。例えば、本実施例は、ステップ１０７，ステップ１１１の更新式に（式２）を用いたが、それを（式３）〜（式５）などに置き換えてもよい。なお、（式５）のような更新式を使用した場合は、理想像モデルの初期値を設定する必要はない。
（第２の実施の形態）第１の実施の形態では、周波数空間でモデルを更新する方法を述べたが、この第２の実施の形態では、実空間でモデルを更新する方法の場合である。
【００４９】図３は、画像復元方法の処理の流れを示しており、この場合、ステップ３０４にＰＳＦモデルｈの初期値として、周波数空間でデルタ関数となるような実空間の分布
【００５０】
【数１３】
を設定する。また、上述した図１のステップ１０２、ステップ１０７、１０８に対応するステップ３０２、ステップ３０７、３０８は除かれ、さらに、ステップ３１０とステップ３１１の順序が反対になるものの、その他のステップ３０１、ステップ３０３〜３０６、ステップ３０９〜３１５での個々の処理は、上述した図１の対応するステップと同じなので、ここでの説明は省略する。また、この場合も、上述した（式６）〜（式７）の更新式を用いる。
【００５１】このようにしても、上述したのと同様な効果が期待できる。なお、上述した第１および第２の実施の形態では、落射蛍光共焦点レーザ走査顕微鏡に適用した例を述べたが、（式１）のような畳込みで形成されるすべての画像に適用できることは言うまでもない。また復元する画像は、上記のような３次元像に限らず、他の次元でも可能である。例えば２次元の画像を復元するには、２次元の実際の像ｇに対して、２次元の理想像モデルとＰＳＦモデルを上記と同様の方法で最適化すればよい。
【００５２】なお、本発明には、以下の発明も含まれる。
（１）請求項１記載の画像復元方法において、ＰＳＦモデルに対する周波数空間の制限がバンドリミットである。
（２）請求項１記載の画像復元方法において、繰り返しの処理が所定の回数に達したことを終了条件とする。
（３）請求項１記載の画像復元方法において、像を取得する装置が共焦点顕微鏡である。
（４）請求項１記載の画像復元方法において、像を取得する装置が顕微鏡である。
【００５３】
【発明の効果】以上述べたように、本発明によれば、周波数空間でデルタ関数となる分布をＰＳＦモデルの初期値に設定したのち、ブラインドデコンボリューションを実施することにより、像の復元を細部まで進行させることができ、細部まで十分に復元された復元像を生成することができる。
【図面の簡単な説明】
【図１】本発明の第１の実施の形態を説明するためのフローチャート。
【図２】第１の実施の形態での復元ＰＳＦと復元像を説明するための図。
【図３】本発明の第２の実施の形態を説明するためのフローチャート。
【図４】従来の画像復元方法の一例を説明するためのフローチャート。
【図５】従来の画像復元方法でのバンドリミットの決め方を説明する図。
【図６】従来の画像復元方法の他の例を説明するためのフローチャート。
【図７】従来例での復元ＰＳＦと復元像を説明するための図。
【符号の説明】
２０１…デルタ関数、
２０２…復元ＰＳＦ、
２０３…真のＰＳＦ、
２０４…復元像、
２０５…真の物体。
【０００１】
【発明の属する技術分野】本発明は、顕微鏡あるいは共焦点顕微鏡などの光学機器により取得した画像の画質を改善する画像復元方法に関するものである。
【０００２】
【従来の技術】顕微鏡あるいは共焦点顕微鏡などの光学機器により物体を観察した場合、観察画像は、元の物体に比べて光学機器のボケ分だけ劣化することが知られている。つまり、このときの観察画像の強度分布ｇは、元の物体の輝度分布ｆと光学機器の点像分布関数（ＰＳＦ）ｈの畳込み（ｃｏｎｖｏｌｕｔｉｏｎ）にノイズｎが加わった、
【０００３】
【数１】
で表現することができる。
【０００４】そして、このようなノイズｎが加わった画像から数値計算を用いて光学機器によるボケを取り除き、理想的な像を得る技術は、像の「復元」あるいは「デコンボリューション（ｄｅｃｏｎｖｏｌｕｔｉｏｎ）」と呼ばれている。
【０００５】ところで、光学顕微鏡では、深さ方向を等間隔で変位させつつ、試料像を撮像し、３次元画像（＝積層画像）を得る光学切片法（ｏｐｔｉｃａｌ ｓｅｃｔｉｏｎｉｎｇ）が用いられている。ところが、光学顕微鏡のＰＳＦは水平方向に比べて深さ方向に大きく広がっているので、積層画像の中の各画像は試料の輝度分布を正確に反映した断面像にはならず、その上下から漏れてくるボケが重畳している。
【０００６】そこで、従来、光学顕微鏡のＰＳＦを基に積層画像からボケを取り除き、理想的な３次元像を「復元」する計算処理が行われてきた。ここで用いられるＰＳＦは、理論値あるいは測定値である。
【０００７】しかし、光学顕微鏡の場合、実際のＰＳＦと理論値は必ずしも一致せず、また正確なＰＳＦの実測も困難である。そこで、「ブラインドデコンボリューション」という像復元の手法が検討されている。この手法は、１つの画像から理想像とＰＳＦの両方を復元するもので、ＰＳＦが未知であっても理想像を復元できるという長所がある。
【０００８】ここで、ブラインドデコンボリューションは、交互法（ａｌｔｅｒｎａｔｉｖｅ ｍｅｔｈｏｄ）というアルゴリズムで実行される。つまり、このような交互法では、初めに、理想像モデルｆとＰＳＦモデルｈに初期値を与え、次に、２つのモデルの畳込みが測定された実際の像ｇになるように両者を交互に修正するとともに、これらの過程を繰り返す。この過程で、モデルには適当な制限条件が加えられ、この交互修正を最適化ループとして何度か繰り返すうちに、２つのモデルは最適化されていく。
【０００９】図４は、従来の交互法を、さらに説明するためのもので、まず、ステップ４０１で、光学機器で測定した実際の像ｇを読み込む。次に、４０２で、実際の像ｇを、実空間から周波数空間へ変換する。これは畳込みの計算を、周波数空間で積として効率よく計算するためである。変換は一般に高速フーリエ変換（ＦＦＴ）で行われる。以下、実空間のものは、小文字（ｇ，ｆ，ｈ）で記し、周波数のものは、大文字（Ｇ，Ｆ，Ｈ）で記す。
【００１０】次に、４０３で、理想像モデルＦの初期値として、測定された実際の像Ｇを設定し、ステップ４０４で、ＰＳＦモデルｈの初期値として、インパルス即ち実空間のＤｉｒａｃのデルタ関数を設定する。そして、ステップ４０５で、制限条件と終了条件を与える。この場合、制限条件は、ＰＳＦの空間周波数には上限、即ちバンドリミットがあるので、このバンドリミットより与えることが可能であるが、現実は、その値は不明なので像のスペクトルから推定して値を設定する。
【００１１】図５（ａ）（ｂ）（ｃ）は、バンドリミットの決め方を説明するもので、ここでは、物体、ＰＳＦおよび像のスペクトルの関係を簡単のために１次元で表している。物体のスペクトル５０１は通常、測定に使用する光学機器のＰＳＦのスペクトル５０２より高周波の成分を含んでいる。よって、測定された像のスペクトル５０４は、ＰＳＦのスペクトル５０２のバンドリミット５０３と同じバンドリミット５０５をもつと考えられるから、この時の像のスペクトル５０４から推定して設定する。なお、３次元像の場合、３軸のバンドリミットを設定しておく。このバンドリミットは後述するステップ４１２で実行される。また、終了条件は一般に、最適化ループの繰り返し回数である。これは後述のステップ４１５で判断される。
【００１２】図４に戻って、ステップ４０６で、ＰＳＦモデルｈを、実空間から周波数空間へ変換する。次いで、ステップ４０７で、ＰＳＦモデルＨとの積が実際の像Ｇになるように、理想像モデルＦを修正する。この場合、理想像モデルＦの修正には、これまでに以下のような更新式が提案されている。・差分
【００１３】
【数２】
・ウィーナフィルタ付き差分
【００１４】
【数３】
・逆フィルタ
【００１５】
【数４】
・ウィーナフィルタ
【００１６】
【数５】
【００１７】次に、ステップ４０８で、理想像モデルＦを、周波数空間から実空間へ変換する。これは、ステップ４０９で実空間の非負制限を理想像モデルｆに適用するためである。変換は一般に高速逆フーリエ変換（ＩＦＦＴ）で行われる。
【００１８】ステップ４０９で、理想像モデルｆに、非負制限を加える。つまり、負の要素値を０で置き換え、それ以外の要素値はそのまま残す。理想像は光の強度分布であるから負の値をとらないからである。さらに、理想像モデルｆの総和が変化しないように補正する例もある。
【００１９】その後、ステップ４１０〜４１４で、理想像モデルとＰＳＦモデルを入れ変えて上述したステップ４０６〜４０９と同様の処理をすることで、ＰＳＦモデルを修正する。
【００２０】なお、ステップ４１２でのＰＳＦモデルＨに対する周波数空間の制限条件は、バンドリミットである。つまり、各軸のバンドリミットより周波数の大きい成分を０で置き換え、それ以外の成分はそのまま残す。また、ステップ４１４での実空間のＰＳＦモデルｈに対する制限条件としては、非負制限の他に、要素値の和を１に規格化する例もある。これは、像形成におけるエネルギー保存則を成立させるためである。
【００２１】その後、ステップ４１５で、最適化ループであるステップ４０６〜４１４の繰り返し回数がステップ４０５で設定した回数に達したら、処理を終了し、ステップ４１６で、最適化ループであるステップ４０６〜４１４による修正を経た理想像モデルｆが、実際の像ｇからボケを除去した復元像として出力される。
【００２２】以上の例ではモデルを周波数空間で修正していたが、実空間で修正する方法もある。その場合、処理の順序が若干変わって図６に示すようになる。この場合、ステップ６１０とステップ６１１は、図４のステップ４１０とステップ４１１ｓの順序が反対になる点が異なるだけでステップ６０１〜６１５での個々の処理は、上述した図４と同じなので、ここでの説明は省略する。なお、このような実空間での修正では、以下のような更新式が提案されている。
・積（multiplicative method)
【００２３】
【数６】
・最尤法（maximum likelihood estimation)（ref.７）
【００２４】
【数７】
【００２５】
【発明が解決しようとする課題】ところが、このようなブラインドデコンボリューションは不適切（ｉｌｌ−ｐｏｓｅｄ）な逆問題で、解の唯一性は保証されていない。つまり、１つの実際の像になる理想像モデルとＰＳＦモデルの組み合わせが無数に存在し、最終的に復元される２つのモデルは、モデルの初期値にも影響される。
【００２６】このため、上述したステップ４０４に示すように、ＰＳＦモデルの初期値として実空間のデルタ関数を取ると、理想像モデルの復元が十分に進行しない傾向が見られた。
【００２７】これを図７を参照して説明すると、実空間のデルタ関数は、図７（ａ）に示すように、周波数空間ですべての成分が１に等しいフラットな分布７０１になる。ここからＰＳＦモデルを出発すると、最適化ループを経るごとにＰＳＦモデルは真のＰＳＦ７０３に対して上から（各周波数成分が大きい方から）近づいてゆく傾向がある。このため、最終的なＰＳＦモデル、即ち復元ＰＳＦ７０２の各周波数成分は、図７（ｂ）に示すように対応する真のＰＳＦ７０３より大きいまま収束することが多い。つまり、理想像モデルはＰＳＦモデルとの畳込みが実際の像になるように最適化されてゆくので、図７R>７（ｃ）に示すように、最終的な理想像モデル、即ち復元像７０４の各周波数成分は、対応する真の物体７０５より小さくなり、実空間で見れば、細部が十分に復元されていない復元像になるという問題があった。本発明は上記事情に鑑みてなされたもので、十分に復元された復元像を得ることができる画像復元方法を提供することを目的とする。
【００２８】
【課題を解決するための手段】請求項１記載の発明は、物体と点像分布関数（ＰＳＦ）の畳込みとして取得されたｎ次元の像から、ＰＳＦによるボケを除去したｎ次元の復元像を得るための画像復元方法において、ｎ次元の理想像モデルに初期値として前記ｎ次元の像を設定するとともに、ｎ次元のＰＳＦモデルに初期値として周波数空間でデルタ関数となる分布を設定した後、前記理想像モデルを実空間あるいは周波数空間の像とＰＳＦモデルとにより修正し、実空間の前記理想像モデルに対して実空間の制限を加え、前記ＰＳＦモデルを実空間あるいは周波数空間の像と理想像モデルとにより修正し、周波数空間の前記ＰＳＦモデルに対して周波数空間の制限を加えた後、実空間に変換し、実空間の前記ＰＳＦモデルに対して実空間の制限を加えるような処理を所定の終了条件を満たすまで繰り返し実行し、最終的な理想像モデルを復元像として生成するようにしている。
【００２９】請求項２記載の発明は、請求項１記載の発明において、前記理想像モデルの初期値として、実際の像を設定している。請求項３記載の発明は、請求項１記載の発明において、前記理想像モデルと前記ＰＳＦモデルに対する実空間の制限が非負制限である。
【００３０】この結果、本発明によれば、周波数空間でデルタ関数となる分布をＰＳＦモデルの初期値にしてブラインドデコンボリューションを実施するが、この場合、周波数空間のデルタ関数は、実空間でフラットな分布、つまり完全なボケに対応したものになるので、この周波数空間のデルタ関数からＰＳＦモデルを出発すると、最適化ループを経るごとにＰＳＦモデルは、真のＰＳＦに対して各周波数成分が小さい方から近づくようになり、対応する復元像の各周波数成分は、対応する真の物体より小さくならず、実空間で見ると、像の復元を細部まで十分に進行させることができる。
【００３１】
【発明の実施の形態】以下、本発明の実施の形態を図面に従い説明する。
（第１の実施の形態）この第１の実施の形態は、周波数空間でモデルを更新する例を示すもので、ここでは、落射蛍光共焦点レーザ走査顕微鏡（ＥＦＣＬＳＭ）の画像復元に適用した例を示している。
【００３２】図１は、本発明による画像復元方法の処理の流れを示しており、ここでは、ＥＦＣＬＳＭで光学切片法により測定した３次元画像（積層画像）ファイルを読み込み、それを画像処理し、ボケを除去した復元像を表示する、あるいは画像ファイルとして出力するプログラムである。ただし、画像処理の部分以外は従来の技術でも説明した一般的な内容であるから説明を省略する。
【００３３】まず、ステップ１０１で、３次元画像ファイルを読み込み、実際の像ｇの３次元複素数データ配列の実部に代入する。ここでは、虚部に０を代入する。つまり、後で高速フーリエ変換（ＦＦＴ）を行うため、ｘ，ｙ，ｚ各軸のデータ点数が２の累乗になっていないときは、足りないデータ点に０詰めして、点数を２の累乗にしておく。また、実際の像ｇの全要素値の和Σｇを求めておき、後のステップ１０９で使用する。
【００３４】次に、ステップ１０２で、実空間の実際の像ｇを、ＦＦＴで周波数空間の実際の像Ｇに変換する。次いで、１０３で、理想像モデルＦを、実際の像Ｇと同じ大きさの３次元複素数データ配列として用意し、その初期値として周波数空間の実際の像Ｇの値を代入し、ステップ１０４で、ＰＳＦモデルＨを、実際の像Ｇと同じ大きさの３次元複素数データ配列として用意し、その初期値として、周波数空間のデルタ関数を代入する。
【００３５】
【数８】
【００３６】ここで、実空間のＰＳＦモデルｈの全要素値の和は、自動的に１になる。次に、ステップ１０５で、制限条件と終了条件を与える。このうちの制限条件は、実際の像Ｇから推定して、ＰＳＦモデルＨに対する３軸のバンドリミットを設定する。また、終了条件として、最適化ループの繰り返し回数Ｎを設定する。
【００３７】ステップ１０６で、実空間のＰＳＦモデルｈをＦＦＴで周波数空間のＰＳＦモデルＨに変換する。ただし、初期値を周波数空間で与えているので、１回目は必要ない。次に、ステップ１０７で、修正されたＰＳＦモデルＨとの積が実際の像Ｇになるように、理想像モデルＦの各成分を次式の差分による更新式で修正する。
【００３８】
【数９】
【００３９】そして、ステップ１０８で、周波数空間の理想像モデルＦを、ＩＦＦＴで実空間の理想像モデルｆに変換し、ステップ１０９で、実空間の理想像モデルｆに非負制限を加える。この場合、まず理想像モデルｆの、負の実部を持つ要素値を０で置き換え、それ以外の要素値は実部を複素数の絶対値で置き換え、虚部を０にする。次に、理想像モデルｆの総和が変化しないように、ステップ１０１で求めたΣｇを使って次式のような規格化を行う。
【００４０】
【数１０】
【００４１】次に、ステップ１１０で、実空間の理想像モデルｆをＦＦＴで周波数空間の理想像モデルＦに変換する。次いで、ステップ１１１で、修正された理想像モデルＦとの積が実際の像Ｇになるように、ＰＳＦモデルＨの各成分を次式の差分による更新式で修正する。
【００４２】
【数１１】
【００４３】次に、ステップ１１２で、周波数空間のＰＳＦモデルＨにバンドリミットを加える。即ち、ステップ１０５で設定した各軸のバンドリミットより周波数の大きい成分を０で置き換え、それ以外の成分はそのまま残す。
【００４４】次に、ステップ１１３で、周波数空間のＰＳＦモデルＨを、ＩＦＦＴで実空間のＰＳＦモデルｈに変換し、ステップ１１４で、実空間のＰＳＦモデルｈに非負制限を加える。この場合、まずＰＳＦモデルｈの、負の実部を持つ要素値を０で置き換え、それ以外の要素値は実部を複素数の絶対値で置き換え、虚部を０にする。次に、ＰＳＦモデルｈの全要素値の和が１になるように、次式のような規格化を行う。
【００４５】
【数１２】
【００４６】その後、最適化ループであるステップ１０６〜１１４の繰り返し回数がステップ１０５で設定した回数Ｎに達したら、ステップ１１５で、処理を終了し、ステップ１１６で、最適化ループであるステップ１０６〜１１４により繰り返し修正された理想像モデルｆの実部が、実際の像ｇからボケを除去した復元像として出力され、復元の画像処理を終了する。
【００４７】従って、このようにすれば、周波数空間でデルタ関数となる分布をＰＳＦモデルの初期値にしたブラインドデコンボリューションが実施されるが、この場合、周波数空間のデルタ関数は、実空間でフラットな分布、つまり完全なボケに対応したものになるので、そこから出発したＰＳＦモデルがシャープになり過ぎることはない。これにより、図２（ａ）に示すように、ＰＳＦモデルを周波数空間のデルタ関数２０１から出発すると、最適化ループを経るごとにＰＳＦモデルは、図２（ｂ）に示すように真のＰＳＦ２０３に対して基本的に下から（各周波数成分が小さい方から）近づくようになり、復元ＰＳＦ２０２の各周波数成分も、対応する真のＰＳＦ２０３より小さい方向で収束する傾向になる。この結果、図２（ｃ）に示すように、対応する復元像２０４の各周波数成分は、対応する真の物体２０５より小さくならないので、実空間で見ると、像の復元を細部まで十分に進行させることができる。つまり、周波数空間でデルタ関数となる分布をＰＳＦモデルの初期値に設定したのち、周波数空間でモデルを更新するブラインドデコンボリューションを実施することにより、像の復元を細部まで進行させることができるようになり、細部まで十分に復元された復元像を生成できる。
【００４８】なお、本発明は、更新式や他の制限条件、終了条件などに関係なく、すべての交互法によるブラインドデコンボリューションに有効である。例えば、本実施例は、ステップ１０７，ステップ１１１の更新式に（式２）を用いたが、それを（式３）〜（式５）などに置き換えてもよい。なお、（式５）のような更新式を使用した場合は、理想像モデルの初期値を設定する必要はない。
（第２の実施の形態）第１の実施の形態では、周波数空間でモデルを更新する方法を述べたが、この第２の実施の形態では、実空間でモデルを更新する方法の場合である。
【００４９】図３は、画像復元方法の処理の流れを示しており、この場合、ステップ３０４にＰＳＦモデルｈの初期値として、周波数空間でデルタ関数となるような実空間の分布
【００５０】
【数１３】
を設定する。また、上述した図１のステップ１０２、ステップ１０７、１０８に対応するステップ３０２、ステップ３０７、３０８は除かれ、さらに、ステップ３１０とステップ３１１の順序が反対になるものの、その他のステップ３０１、ステップ３０３〜３０６、ステップ３０９〜３１５での個々の処理は、上述した図１の対応するステップと同じなので、ここでの説明は省略する。また、この場合も、上述した（式６）〜（式７）の更新式を用いる。
【００５１】このようにしても、上述したのと同様な効果が期待できる。なお、上述した第１および第２の実施の形態では、落射蛍光共焦点レーザ走査顕微鏡に適用した例を述べたが、（式１）のような畳込みで形成されるすべての画像に適用できることは言うまでもない。また復元する画像は、上記のような３次元像に限らず、他の次元でも可能である。例えば２次元の画像を復元するには、２次元の実際の像ｇに対して、２次元の理想像モデルとＰＳＦモデルを上記と同様の方法で最適化すればよい。
【００５２】なお、本発明には、以下の発明も含まれる。
（１）請求項１記載の画像復元方法において、ＰＳＦモデルに対する周波数空間の制限がバンドリミットである。
（２）請求項１記載の画像復元方法において、繰り返しの処理が所定の回数に達したことを終了条件とする。
（３）請求項１記載の画像復元方法において、像を取得する装置が共焦点顕微鏡である。
（４）請求項１記載の画像復元方法において、像を取得する装置が顕微鏡である。
【００５３】
【発明の効果】以上述べたように、本発明によれば、周波数空間でデルタ関数となる分布をＰＳＦモデルの初期値に設定したのち、ブラインドデコンボリューションを実施することにより、像の復元を細部まで進行させることができ、細部まで十分に復元された復元像を生成することができる。
【図面の簡単な説明】
【図１】本発明の第１の実施の形態を説明するためのフローチャート。
【図２】第１の実施の形態での復元ＰＳＦと復元像を説明するための図。
【図３】本発明の第２の実施の形態を説明するためのフローチャート。
【図４】従来の画像復元方法の一例を説明するためのフローチャート。
【図５】従来の画像復元方法でのバンドリミットの決め方を説明する図。
【図６】従来の画像復元方法の他の例を説明するためのフローチャート。
【図７】従来例での復元ＰＳＦと復元像を説明するための図。
【符号の説明】
２０１…デルタ関数、
２０２…復元ＰＳＦ、
２０３…真のＰＳＦ、
２０４…復元像、
２０５…真の物体。
【特許請求の範囲】
【請求項１】 物体と点像分布関数（ＰＳＦ）の畳込みとして取得されたｎ次元の像から、ＰＳＦによるボケを除去したｎ次元の復元像を得るための画像復元方法において、ｎ次元の理想像モデルに初期値として前記ｎ次元の像を設定するとともに、ｎ次元のＰＳＦモデルに初期値として周波数空間でデルタ関数となる分布を設定した後、前記理想像モデルを実空間あるいは周波数空間の像とＰＳＦモデルとにより修正し、実空間の前記理想像モデルに対して実空間の制限を加え、前記ＰＳＦモデルを実空間あるいは周波数空間の像と理想像モデルとにより修正し、周波数空間の前記ＰＳＦモデルに対して周波数空間の制限を加えた後、実空間に変換し、実空間の前記ＰＳＦモデルに対して実空間の制限を加えるような処理を所定の終了条件を満たすまで繰り返し実行し、最終的な理想像モデルを復元像として生成することを特徴とする画像復元方法。
【請求項２】 前記理想像モデルの初期値として、実際の像を設定することを特徴とする請求項１記載の画像復元方法。
【請求項３】 前記理想像モデルと前記ＰＳＦモデルに対する実空間の制限が非負制限であることを特徴とする請求項１記載の画像復元方法。
【請求項１】 物体と点像分布関数（ＰＳＦ）の畳込みとして取得されたｎ次元の像から、ＰＳＦによるボケを除去したｎ次元の復元像を得るための画像復元方法において、ｎ次元の理想像モデルに初期値として前記ｎ次元の像を設定するとともに、ｎ次元のＰＳＦモデルに初期値として周波数空間でデルタ関数となる分布を設定した後、前記理想像モデルを実空間あるいは周波数空間の像とＰＳＦモデルとにより修正し、実空間の前記理想像モデルに対して実空間の制限を加え、前記ＰＳＦモデルを実空間あるいは周波数空間の像と理想像モデルとにより修正し、周波数空間の前記ＰＳＦモデルに対して周波数空間の制限を加えた後、実空間に変換し、実空間の前記ＰＳＦモデルに対して実空間の制限を加えるような処理を所定の終了条件を満たすまで繰り返し実行し、最終的な理想像モデルを復元像として生成することを特徴とする画像復元方法。
【請求項２】 前記理想像モデルの初期値として、実際の像を設定することを特徴とする請求項１記載の画像復元方法。
【請求項３】 前記理想像モデルと前記ＰＳＦモデルに対する実空間の制限が非負制限であることを特徴とする請求項１記載の画像復元方法。
【図１】
【図２】
【図７】
【図３】
【図４】
【図５】
【図６】
【図２】
【図７】
【図３】
【図４】
【図５】
【図６】
【公開番号】特開２０００−４３６３（Ｐ２０００−４３６３Ａ）
【公開日】平成１２年１月７日（２０００．１．７）
【国際特許分類】
【出願番号】特願平１０−１７００４４
【出願日】平成１０年６月１７日（１９９８．６．１７）
【出願人】（００００００３７６）オリンパス光学工業株式会社 (11,466)
【Ｆターム（参考）】
【公開日】平成１２年１月７日（２０００．１．７）
【国際特許分類】
【出願日】平成１０年６月１７日（１９９８．６．１７）
【出願人】（００００００３７６）オリンパス光学工業株式会社 (11,466)
【Ｆターム（参考）】
[ Back to top ]