JP2000004363A

JP2000004363A - 画像復元方法

Info

Publication number: JP2000004363A
Application number: JP17004498A
Authority: JP
Inventors: Toshiaki Matsuzawa; 聡明松沢
Original assignee: Olympus Optical Co Ltd
Current assignee: Olympus Corp
Priority date: 1998-06-17
Filing date: 1998-06-17
Publication date: 2000-01-07

Abstract

(57)【要約】【課題】本発明は、十分に復元された復元像を得ること
ができる画像復元方法を提供する。【解決手段】ステップ１０１〜１０４でｎ次元の理想像
モデルの初期値を設定すると共に、ｎ次元のＰＳＦモデ
ルの初期値として周波数空間でデルタ関数となる分布を
設定し、ステップ１０５で制限条件と終了条件を与えた
後、ステップ１０６〜１１４で理想像モデルを実空間あ
るいは周波数空間で像とＰＳＦモデルにより修正し、実
空間の理想像モデルに対して実空間の制限を加え、ＰＳ
Ｆモデルを実空間あるいは周波数空間で像と理想像モデ
ルにより修正し、周波数空間のＰＳＦモデルに対して周
波数空間の制限を加えた後実空間に変換し、実空間のＰ
ＳＦモデルに対して実空間の制限を加える処理を、ステ
ップ１１５で所定の終了条件を満たすまで繰返し実行
し、最終的な理想像モデルを復元像として生成する。

Description

【発明の詳細な説明】

【０００１】

【発明の属する技術分野】本発明は、顕微鏡あるいは共
焦点顕微鏡などの光学機器により取得した画像の画質を
改善する画像復元方法に関するものである。

【０００２】

【従来の技術】顕微鏡あるいは共焦点顕微鏡などの光学
機器により物体を観察した場合、観察画像は、元の物体
に比べて光学機器のボケ分だけ劣化することが知られて
いる。つまり、このときの観察画像の強度分布ｇは、元
の物体の輝度分布ｆと光学機器の点像分布関数（ＰＳ
Ｆ）ｈの畳込み（ｃｏｎｖｏｌｕｔｉｏｎ）にノイズｎ
が加わった、

【０００３】

【数１】で表現することができる。

【０００４】そして、このようなノイズｎが加わった画
像から数値計算を用いて光学機器によるボケを取り除
き、理想的な像を得る技術は、像の「復元」あるいは
「デコンボリューション（ｄｅｃｏｎｖｏｌｕｔｉｏ
ｎ）」と呼ばれている。

【０００５】ところで、光学顕微鏡では、深さ方向を等
間隔で変位させつつ、試料像を撮像し、３次元画像（＝
積層画像）を得る光学切片法（ｏｐｔｉｃａｌｓｅｃ
ｔｉｏｎｉｎｇ）が用いられている。ところが、光学顕
微鏡のＰＳＦは水平方向に比べて深さ方向に大きく広が
っているので、積層画像の中の各画像は試料の輝度分布
を正確に反映した断面像にはならず、その上下から漏れ
てくるボケが重畳している。

【０００６】そこで、従来、光学顕微鏡のＰＳＦを基に
積層画像からボケを取り除き、理想的な３次元像を「復
元」する計算処理が行われてきた。ここで用いられるＰ
ＳＦは、理論値あるいは測定値である。

【０００７】しかし、光学顕微鏡の場合、実際のＰＳＦ
と理論値は必ずしも一致せず、また正確なＰＳＦの実測
も困難である。そこで、「ブラインドデコンボリューシ
ョン」という像復元の手法が検討されている。この手法
は、１つの画像から理想像とＰＳＦの両方を復元するも
ので、ＰＳＦが未知であっても理想像を復元できるとい
う長所がある。

【０００８】ここで、ブラインドデコンボリューション
は、交互法（ａｌｔｅｒｎａｔｉｖｅｍｅｔｈｏｄ）
というアルゴリズムで実行される。つまり、このような
交互法では、初めに、理想像モデルｆとＰＳＦモデルｈ
に初期値を与え、次に、２つのモデルの畳込みが測定さ
れた実際の像ｇになるように両者を交互に修正するとと
もに、これらの過程を繰り返す。この過程で、モデルに
は適当な制限条件が加えられ、この交互修正を最適化ル
ープとして何度か繰り返すうちに、２つのモデルは最適
化されていく。

【０００９】図４は、従来の交互法を、さらに説明する
ためのもので、まず、ステップ４０１で、光学機器で測
定した実際の像ｇを読み込む。次に、４０２で、実際の
像ｇを、実空間から周波数空間へ変換する。これは畳込
みの計算を、周波数空間で積として効率よく計算するた
めである。変換は一般に高速フーリエ変換（ＦＦＴ）で
行われる。以下、実空間のものは、小文字（ｇ，ｆ，
ｈ）で記し、周波数のものは、大文字（Ｇ，Ｆ，Ｈ）で
記す。

【００１０】次に、４０３で、理想像モデルＦの初期値
として、測定された実際の像Ｇを設定し、ステップ４０
４で、ＰＳＦモデルｈの初期値として、インパルス即ち
実空間のＤｉｒａｃのデルタ関数を設定する。そして、
ステップ４０５で、制限条件と終了条件を与える。この
場合、制限条件は、ＰＳＦの空間周波数には上限、即ち
バンドリミットがあるので、このバンドリミットより与
えることが可能であるが、現実は、その値は不明なので
像のスペクトルから推定して値を設定する。

【００１１】図５（ａ）（ｂ）（ｃ）は、バンドリミッ
トの決め方を説明するもので、ここでは、物体、ＰＳＦ
および像のスペクトルの関係を簡単のために１次元で表
している。物体のスペクトル５０１は通常、測定に使用
する光学機器のＰＳＦのスペクトル５０２より高周波の
成分を含んでいる。よって、測定された像のスペクトル
５０４は、ＰＳＦのスペクトル５０２のバンドリミット
５０３と同じバンドリミット５０５をもつと考えられる
から、この時の像のスペクトル５０４から推定して設定
する。なお、３次元像の場合、３軸のバンドリミットを
設定しておく。このバンドリミットは後述するステップ
４１２で実行される。また、終了条件は一般に、最適化
ループの繰り返し回数である。これは後述のステップ４
１５で判断される。

【００１２】図４に戻って、ステップ４０６で、ＰＳＦ
モデルｈを、実空間から周波数空間へ変換する。次い
で、ステップ４０７で、ＰＳＦモデルＨとの積が実際の
像Ｇになるように、理想像モデルＦを修正する。この場
合、理想像モデルＦの修正には、これまでに以下のよう
な更新式が提案されている。・差分

【００１３】

【数２】・ウィーナフィルタ付き差分

【００１４】

【数３】・逆フィルタ

【００１５】

【数４】・ウィーナフィルタ

【００１６】

【数５】

【００１７】次に、ステップ４０８で、理想像モデルＦ
を、周波数空間から実空間へ変換する。これは、ステッ
プ４０９で実空間の非負制限を理想像モデルｆに適用す
るためである。変換は一般に高速逆フーリエ変換（ＩＦ
ＦＴ）で行われる。

【００１８】ステップ４０９で、理想像モデルｆに、非
負制限を加える。つまり、負の要素値を０で置き換え、
それ以外の要素値はそのまま残す。理想像は光の強度分
布であるから負の値をとらないからである。さらに、理
想像モデルｆの総和が変化しないように補正する例もあ
る。

【００１９】その後、ステップ４１０〜４１４で、理想
像モデルとＰＳＦモデルを入れ変えて上述したステップ
４０６〜４０９と同様の処理をすることで、ＰＳＦモデ
ルを修正する。

【００２０】なお、ステップ４１２でのＰＳＦモデルＨ
に対する周波数空間の制限条件は、バンドリミットであ
る。つまり、各軸のバンドリミットより周波数の大きい
成分を０で置き換え、それ以外の成分はそのまま残す。
また、ステップ４１４での実空間のＰＳＦモデルｈに対
する制限条件としては、非負制限の他に、要素値の和を
１に規格化する例もある。これは、像形成におけるエネ
ルギー保存則を成立させるためである。

【００２１】その後、ステップ４１５で、最適化ループ
であるステップ４０６〜４１４の繰り返し回数がステッ
プ４０５で設定した回数に達したら、処理を終了し、ス
テップ４１６で、最適化ループであるステップ４０６〜
４１４による修正を経た理想像モデルｆが、実際の像ｇ
からボケを除去した復元像として出力される。

【００２２】以上の例ではモデルを周波数空間で修正し
ていたが、実空間で修正する方法もある。その場合、処
理の順序が若干変わって図６に示すようになる。この場
合、ステップ６１０とステップ６１１は、図４のステッ
プ４１０とステップ４１１ｓの順序が反対になる点が異
なるだけでステップ６０１〜６１５での個々の処理は、
上述した図４と同じなので、ここでの説明は省略する。
なお、このような実空間での修正では、以下のような更
新式が提案されている。・積（multiplicative method)

【００２３】

【数６】・最尤法（maximum likelihood estimation)（ref.７）

【００２４】

【数７】

【００２５】

【発明が解決しようとする課題】ところが、このような
ブラインドデコンボリューションは不適切（ｉｌｌ−ｐ
ｏｓｅｄ）な逆問題で、解の唯一性は保証されていな
い。つまり、１つの実際の像になる理想像モデルとＰＳ
Ｆモデルの組み合わせが無数に存在し、最終的に復元さ
れる２つのモデルは、モデルの初期値にも影響される。

【００２６】このため、上述したステップ４０４に示す
ように、ＰＳＦモデルの初期値として実空間のデルタ関
数を取ると、理想像モデルの復元が十分に進行しない傾
向が見られた。

【００２７】これを図７を参照して説明すると、実空間
のデルタ関数は、図７（ａ）に示すように、周波数空間
ですべての成分が１に等しいフラットな分布７０１にな
る。ここからＰＳＦモデルを出発すると、最適化ループ
を経るごとにＰＳＦモデルは真のＰＳＦ７０３に対して
上から（各周波数成分が大きい方から）近づいてゆく傾
向がある。このため、最終的なＰＳＦモデル、即ち復元
ＰＳＦ７０２の各周波数成分は、図７（ｂ）に示すよう
に対応する真のＰＳＦ７０３より大きいまま収束するこ
とが多い。つまり、理想像モデルはＰＳＦモデルとの畳
込みが実際の像になるように最適化されてゆくので、図
７（ｃ）に示すように、最終的な理想像モデル、即ち復
元像７０４の各周波数成分は、対応する真の物体７０５
より小さくなり、実空間で見れば、細部が十分に復元さ
れていない復元像になるという問題があった。本発明は
上記事情に鑑みてなされたもので、十分に復元された復
元像を得ることができる画像復元方法を提供することを
目的とする。

【００２８】

【課題を解決するための手段】請求項１記載の発明は、
物体と点像分布関数（ＰＳＦ）の畳込みとして取得され
たｎ次元の像から、ＰＳＦによるボケを除去したｎ次元
の復元像を得るための画像復元方法において、ｎ次元の
理想像モデルに初期値として前記ｎ次元の像を設定する
とともに、ｎ次元のＰＳＦモデルに初期値として周波数
空間でデルタ関数となる分布を設定した後、前記理想像
モデルを実空間あるいは周波数空間の像とＰＳＦモデル
とにより修正し、実空間の前記理想像モデルに対して実
空間の制限を加え、前記ＰＳＦモデルを実空間あるいは
周波数空間の像と理想像モデルとにより修正し、周波数
空間の前記ＰＳＦモデルに対して周波数空間の制限を加
えた後、実空間に変換し、実空間の前記ＰＳＦモデルに
対して実空間の制限を加えるような処理を所定の終了条
件を満たすまで繰り返し実行し、最終的な理想像モデル
を復元像として生成するようにしている。

【００２９】請求項２記載の発明は、請求項１記載の発
明において、前記理想像モデルの初期値として、実際の
像を設定している。請求項３記載の発明は、請求項１記
載の発明において、前記理想像モデルと前記ＰＳＦモデ
ルに対する実空間の制限が非負制限である。

【００３０】この結果、本発明によれば、周波数空間で
デルタ関数となる分布をＰＳＦモデルの初期値にしてブ
ラインドデコンボリューションを実施するが、この場
合、周波数空間のデルタ関数は、実空間でフラットな分
布、つまり完全なボケに対応したものになるので、この
周波数空間のデルタ関数からＰＳＦモデルを出発する
と、最適化ループを経るごとにＰＳＦモデルは、真のＰ
ＳＦに対して各周波数成分が小さい方から近づくように
なり、対応する復元像の各周波数成分は、対応する真の
物体より小さくならず、実空間で見ると、像の復元を細
部まで十分に進行させることができる。

【００３１】

【発明の実施の形態】以下、本発明の実施の形態を図面
に従い説明する。（第１の実施の形態）この第１の実施の形態は、周波数
空間でモデルを更新する例を示すもので、ここでは、落
射蛍光共焦点レーザ走査顕微鏡（ＥＦＣＬＳＭ）の画像
復元に適用した例を示している。

【００３２】図１は、本発明による画像復元方法の処理
の流れを示しており、ここでは、ＥＦＣＬＳＭで光学切
片法により測定した３次元画像（積層画像）ファイルを
読み込み、それを画像処理し、ボケを除去した復元像を
表示する、あるいは画像ファイルとして出力するプログ
ラムである。ただし、画像処理の部分以外は従来の技術
でも説明した一般的な内容であるから説明を省略する。

【００３３】まず、ステップ１０１で、３次元画像ファ
イルを読み込み、実際の像ｇの３次元複素数データ配列
の実部に代入する。ここでは、虚部に０を代入する。つ
まり、後で高速フーリエ変換（ＦＦＴ）を行うため、
ｘ，ｙ，ｚ各軸のデータ点数が２の累乗になっていない
ときは、足りないデータ点に０詰めして、点数を２の累
乗にしておく。また、実際の像ｇの全要素値の和Σｇを
求めておき、後のステップ１０９で使用する。

【００３４】次に、ステップ１０２で、実空間の実際の
像ｇを、ＦＦＴで周波数空間の実際の像Ｇに変換する。
次いで、１０３で、理想像モデルＦを、実際の像Ｇと同
じ大きさの３次元複素数データ配列として用意し、その
初期値として周波数空間の実際の像Ｇの値を代入し、ス
テップ１０４で、ＰＳＦモデルＨを、実際の像Ｇと同じ
大きさの３次元複素数データ配列として用意し、その初
期値として、周波数空間のデルタ関数を代入する。

【００３５】

【数８】

【００３６】ここで、実空間のＰＳＦモデルｈの全要素
値の和は、自動的に１になる。次に、ステップ１０５
で、制限条件と終了条件を与える。このうちの制限条件
は、実際の像Ｇから推定して、ＰＳＦモデルＨに対する
３軸のバンドリミットを設定する。また、終了条件とし
て、最適化ループの繰り返し回数Ｎを設定する。

【００３７】ステップ１０６で、実空間のＰＳＦモデル
ｈをＦＦＴで周波数空間のＰＳＦモデルＨに変換する。
ただし、初期値を周波数空間で与えているので、１回目
は必要ない。次に、ステップ１０７で、修正されたＰＳ
ＦモデルＨとの積が実際の像Ｇになるように、理想像モ
デルＦの各成分を次式の差分による更新式で修正する。

【００３８】

【数９】

【００３９】そして、ステップ１０８で、周波数空間の
理想像モデルＦを、ＩＦＦＴで実空間の理想像モデルｆ
に変換し、ステップ１０９で、実空間の理想像モデルｆ
に非負制限を加える。この場合、まず理想像モデルｆ
の、負の実部を持つ要素値を０で置き換え、それ以外の
要素値は実部を複素数の絶対値で置き換え、虚部を０に
する。次に、理想像モデルｆの総和が変化しないよう
に、ステップ１０１で求めたΣｇを使って次式のような
規格化を行う。

【００４０】

【数１０】

【００４１】次に、ステップ１１０で、実空間の理想像
モデルｆをＦＦＴで周波数空間の理想像モデルＦに変換
する。次いで、ステップ１１１で、修正された理想像モ
デルＦとの積が実際の像Ｇになるように、ＰＳＦモデル
Ｈの各成分を次式の差分による更新式で修正する。

【００４２】

【数１１】

【００４３】次に、ステップ１１２で、周波数空間のＰ
ＳＦモデルＨにバンドリミットを加える。即ち、ステッ
プ１０５で設定した各軸のバンドリミットより周波数の
大きい成分を０で置き換え、それ以外の成分はそのまま
残す。

【００４４】次に、ステップ１１３で、周波数空間のＰ
ＳＦモデルＨを、ＩＦＦＴで実空間のＰＳＦモデルｈに
変換し、ステップ１１４で、実空間のＰＳＦモデルｈに
非負制限を加える。この場合、まずＰＳＦモデルｈの、
負の実部を持つ要素値を０で置き換え、それ以外の要素
値は実部を複素数の絶対値で置き換え、虚部を０にす
る。次に、ＰＳＦモデルｈの全要素値の和が１になるよ
うに、次式のような規格化を行う。

【００４５】

【数１２】

【００４６】その後、最適化ループであるステップ１０
６〜１１４の繰り返し回数がステップ１０５で設定した
回数Ｎに達したら、ステップ１１５で、処理を終了し、
ステップ１１６で、最適化ループであるステップ１０６
〜１１４により繰り返し修正された理想像モデルｆの実
部が、実際の像ｇからボケを除去した復元像として出力
され、復元の画像処理を終了する。

【００４７】従って、このようにすれば、周波数空間で
デルタ関数となる分布をＰＳＦモデルの初期値にしたブ
ラインドデコンボリューションが実施されるが、この場
合、周波数空間のデルタ関数は、実空間でフラットな分
布、つまり完全なボケに対応したものになるので、そこ
から出発したＰＳＦモデルがシャープになり過ぎること
はない。これにより、図２（ａ）に示すように、ＰＳＦ
モデルを周波数空間のデルタ関数２０１から出発する
と、最適化ループを経るごとにＰＳＦモデルは、図２
（ｂ）に示すように真のＰＳＦ２０３に対して基本的に
下から（各周波数成分が小さい方から）近づくようにな
り、復元ＰＳＦ２０２の各周波数成分も、対応する真の
ＰＳＦ２０３より小さい方向で収束する傾向になる。こ
の結果、図２（ｃ）に示すように、対応する復元像２０
４の各周波数成分は、対応する真の物体２０５より小さ
くならないので、実空間で見ると、像の復元を細部まで
十分に進行させることができる。つまり、周波数空間で
デルタ関数となる分布をＰＳＦモデルの初期値に設定し
たのち、周波数空間でモデルを更新するブラインドデコ
ンボリューションを実施することにより、像の復元を細
部まで進行させることができるようになり、細部まで十
分に復元された復元像を生成できる。

【００４８】なお、本発明は、更新式や他の制限条件、
終了条件などに関係なく、すべての交互法によるブライ
ンドデコンボリューションに有効である。例えば、本実
施例は、ステップ１０７，ステップ１１１の更新式に
（式２）を用いたが、それを（式３）〜（式５）などに
置き換えてもよい。なお、（式５）のような更新式を使
用した場合は、理想像モデルの初期値を設定する必要は
ない。（第２の実施の形態）第１の実施の形態では、周波数空
間でモデルを更新する方法を述べたが、この第２の実施
の形態では、実空間でモデルを更新する方法の場合であ
る。

【００４９】図３は、画像復元方法の処理の流れを示し
ており、この場合、ステップ３０４にＰＳＦモデルｈの
初期値として、周波数空間でデルタ関数となるような実
空間の分布

【００５０】

【数１３】を設定する。また、上述した図１のステップ１０２、ス
テップ１０７、１０８に対応するステップ３０２、ステ
ップ３０７、３０８は除かれ、さらに、ステップ３１０
とステップ３１１の順序が反対になるものの、その他の
ステップ３０１、ステップ３０３〜３０６、ステップ３
０９〜３１５での個々の処理は、上述した図１の対応す
るステップと同じなので、ここでの説明は省略する。ま
た、この場合も、上述した（式６）〜（式７）の更新式
を用いる。

【００５１】このようにしても、上述したのと同様な効
果が期待できる。なお、上述した第１および第２の実施
の形態では、落射蛍光共焦点レーザ走査顕微鏡に適用し
た例を述べたが、（式１）のような畳込みで形成される
すべての画像に適用できることは言うまでもない。また
復元する画像は、上記のような３次元像に限らず、他の
次元でも可能である。例えば２次元の画像を復元するに
は、２次元の実際の像ｇに対して、２次元の理想像モデ
ルとＰＳＦモデルを上記と同様の方法で最適化すればよ
い。

【００５２】なお、本発明には、以下の発明も含まれ
る。（１）請求項１記載の画像復元方法において、ＰＳＦモ
デルに対する周波数空間の制限がバンドリミットであ
る。（２）請求項１記載の画像復元方法において、繰り返し
の処理が所定の回数に達したことを終了条件とする。（３）請求項１記載の画像復元方法において、像を取得
する装置が共焦点顕微鏡である。（４）請求項１記載の画像復元方法において、像を取得
する装置が顕微鏡である。

【００５３】

【発明の効果】以上述べたように、本発明によれば、周
波数空間でデルタ関数となる分布をＰＳＦモデルの初期
値に設定したのち、ブラインドデコンボリューションを
実施することにより、像の復元を細部まで進行させるこ
とができ、細部まで十分に復元された復元像を生成する
ことができる。

【図面の簡単な説明】

【図１】本発明の第１の実施の形態を説明するためのフ
ローチャート。

【図２】第１の実施の形態での復元ＰＳＦと復元像を説
明するための図。

【図３】本発明の第２の実施の形態を説明するためのフ
ローチャート。

【図４】従来の画像復元方法の一例を説明するためのフ
ローチャート。

【図５】従来の画像復元方法でのバンドリミットの決め
方を説明する図。

【図６】従来の画像復元方法の他の例を説明するための
フローチャート。

【図７】従来例での復元ＰＳＦと復元像を説明するため
の図。

【符号の説明】

２０１…デルタ関数、２０２…復元ＰＳＦ、２０３…真のＰＳＦ、２０４…復元像、２０５…真の物体。

Claims

【特許請求の範囲】

【請求項１】物体と点像分布関数（ＰＳＦ）の畳込み
として取得されたｎ次元の像から、ＰＳＦによるボケを
除去したｎ次元の復元像を得るための画像復元方法にお
いて、ｎ次元の理想像モデルに初期値として前記ｎ次元の像を
設定するとともに、ｎ次元のＰＳＦモデルに初期値とし
て周波数空間でデルタ関数となる分布を設定した後、前記理想像モデルを実空間あるいは周波数空間の像とＰ
ＳＦモデルとにより修正し、実空間の前記理想像モデルに対して実空間の制限を加
え、前記ＰＳＦモデルを実空間あるいは周波数空間の像と理
想像モデルとにより修正し、周波数空間の前記ＰＳＦモデルに対して周波数空間の制
限を加えた後、実空間に変換し、実空間の前記ＰＳＦモデルに対して実空間の制限を加え
るような処理を所定の終了条件を満たすまで繰り返し実
行し、最終的な理想像モデルを復元像として生成するこ
とを特徴とする画像復元方法。
【請求項２】前記理想像モデルの初期値として、実際
の像を設定することを特徴とする請求項１記載の画像復
元方法。
【請求項３】前記理想像モデルと前記ＰＳＦモデルに
対する実空間の制限が非負制限であることを特徴とする
請求項１記載の画像復元方法。