(57)【要約】真のサポートが正確にわかっていない場合においても、高精度に画像を回復することができる画像再構成方法。本方法では、まず、初期画像Ｉをｇｉｎｉｔｉａｌとする（Ｓ１３００）。そして、測定されたサポートに膨張処理を施し（Ｓ１４００）、サポートＤを表す画像ｄを生成する（Ｓ１５００）。そして、画像ｄにスネークスを適用し（Ｓ１７００）、抽出した物体Ｄ’を新たなサポートＤとする（Ｓ１８００）。そして、得られたサポートＤと測定された原画像のフーリエ振幅｜Ｆ｜を用いて、ｇｉｎｉｔｉａｌに対してＭ回ＥＲアルゴリズムを適用し、出力画像ｇｎを得る（Ｓ１９００）。そして、得られたｇｎをそれぞれｇｉｎｉｔｉａｌ、ｄとする（Ｓ２０００、Ｓ２１００）。そして、Ｓ１７００〜Ｓ２１００を所定の繰り返し回数（Ｎ回）実行し、画像を再構成する。Ｎ回後に得られる出力画像ｇＮが再構成画像である。【特許請求の範囲】
【請求項１】
原画像のフーリエ振幅を入力するステップと、
原画像が占めるサポートを入力するステップと、
入力したサポートに膨張処理を施すステップと、
入力したフーリエ振幅および膨張処理を施したサポートを用いて、位相回復アルゴリズムおよびサポートを収縮させる機能を有する物体抽出アルゴリズムを併用してサポート条件を更新することにより、画像を再構成するステップと、
を有する画像再構成方法。
【請求項２】
前記位相回復アルゴリズムはエラーリダクションアルゴリズム（ＥＲ）であり、前記物体抽出アルゴリズムはスネークスアルゴリズム（Ｓｎａｋｅｓ）である請求項１記載の画像再構成方法。
【請求項３】
前記サポート条件の更新は、前記位相回復アルゴリズムを一または複数回繰り返すことにより得られる出力画像に前記物体抽出アルゴリズムを適用して画像を抽出し、抽出した画像を新たなサポート条件とすることによって行われる請求項１記載の画像再構成方法。
【請求項４】
前記位相回復アルゴリズムおよび前記物体抽出アルゴリズムを適用する順番はあらかじめ決められている請求項１記載の画像再構成方法。
【請求項５】
前記膨張処理の後に適用する最初のアルゴリズムは前記物体抽出アルゴリズムである請求項４記載の画像再構成方法。
【請求項６】
終了時のアルゴリズムは前記位相回復アルゴリズムである請求項４記載の画像再構成方法。
【請求項７】
原画像のフーリエ振幅を入力するステップと、
原画像が占めるサポートを入力するステップと、
入力したサポートに膨張処理を施すステップと、
入力したフーリエ振幅および膨張処理を施したサポート領域を用いて、位相回復アルゴリズムおよびサポートを収縮させる機能を有する物体抽出アルゴリズムを併用してサポート条件を更新することにより、画像を再構成するステップと、
をコンピュータに実行させる画像再構成プログラム。
【請求項８】
最適化アルゴリズムを用いて、特定の周波数において構戒される領域の位相の組み合わせが原画像と一致する初期画像を探索するステップと、
探索した初期画像を用いて、位相回復アルゴリズムにより画像を再構成するステップと、
を有する画像再構成方法。
【請求項９】
前記最適化アルゴリズムは遺伝的アルゴリズム（ＧＡ）であり、前記位相回復アルゴリズムはエラーリダクションアルゴリズム（ＥＲ）である請求項８記載の画像再構成方法。
【請求項１０】
前記探索ステップの前に、入力したフーリエ振幅にローパスフィルタ処理を施し、フーリエ振幅の要素数を減らすステップをさらに有する請求項８記載の画像再構成方法。
【請求項１１】
最適化アルゴリズムを用いて、特定の周波数において構成される領域の位相の組み合わせが原画像と一致する初期画像を探索するステップと、
探索した初期画像を用いて、位相回復アルゴリズムにより画像を再構成するステップと、
をコンピュータに実行させる画像再構成プログラム。
【請求項１２】
複数の時間フレームで構成される動画像を入力するステップと、
前記複数の時間フレームに対して、連続する２つの時間フレームの間の動き情報を算出するステップと、
算出した各時間フレームに対応する動き情報を用いて、各時間フレームのフーリエ振幅およびサポートを更新するステップと、
更新した各時間フレームのフーリエ振幅およびサポートをそれぞれ入力データとして請求項１記載の画像再構成方法に適用することにより、各時間フレームの再構成画像を導出するステップと、
を有する画像再構成方法。
【請求項１３】
前記動き情報は、平行移動、回転、拡大、および縮小を行列の形で表す動き関数Ｆによって表現され、
各時間フレームに対応する動き情報は、第（ｋ−１）フレーム（ここで、ｋ＝１，２，‥，Ｋ：Ｋは整数）の再構成画像から得られるサポートＤｋ−１を動き関数Ｆに適用して得られるサポートＦ（Ｄｋ−１）と、第ｋフレームの測定サポートＤｋとの差の絶対値の二乗が最小になるときの動き関数ＦＯＰＴとして算出される、
請求項１２記載の画像再構成方法。
【請求項１４】
前記動き情報は、平行移動、回転、拡大、および縮小を行列の形で表す動き関数Ｆによって表現され、
各時間フレームに対応する動き情報は、
第（ｋ−１）フレーム（ここで、ｋ＝１，２，‥，Ｋ：Ｋは整数）の再構成画像から得られるサポートＤｋ−１を動き関数Ｆに適用して得られるサポートＦ（Ｄｋ−１）と、第ｋフレームの測定サポートＤｋとの差の絶対値の二乗が最小になる動き関数ＦＯＰＴを算出するステップと、
算出した関数ＦＯＰＴに対して平行移動、回転、拡大、および縮小の各パラメータの値を変更して得られる関数Ｆ（ｎ）（ここで、ｎ＝１，２，‥，Ｎ：Ｎは整数）を算出するステップと、
算出した関数ＦＯＰＴおよび関数Ｆ（ｎ）の（Ｎ＋１）個の関数の中から、フーリエ誤差が最小になる関数Ｆ（ｉ）を選択するステップと、
により算出される、
請求項１２記載の画像再構成方法。
【請求項１５】
複数の時間フレームで構成される動画像を入力するステップと、
前記複数の時間フレームに対して、連続する２つの時間フレームの間の動き情報を算出するステップと、
算出した各時間フレームに対応する動き情報を用いて、各時間フレームのフーリエ振幅およびサポートを更新するステップと、
更新した各時間フレームのフーリエ振幅およびサポートをそれぞれ入力データとして請求項１記載の画像再構成方法に適用することにより、各時間フレームの再構成画像を導出するステップと、
をコンピュータに実行させる画像再構成プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】  本発明は、画像再構成方法に関し、特に、位相回復アルゴリズムを用いた画像再構成の高精度化のための画像領域の拘束条件の適応的設定に関する。
【背景技術】
【０００２】  電子顕微鏡や天体観測、Ｘ線結晶学などの様々な分野で位相回復法が必要とされている。位相回復法とは、観測で得られるフーリエ振幅を用いて、失われた位相を回復することにより、画像を再構成する方法であり、様々な手法が提案されている（非特許文献１〜非特許文献４）。なお、位相回復法におけるアルゴリズムは、一般に、「位相回復アルゴリズム」と呼ばれる。【０００３】  その代表的な手法としては、主にＦｉｅｎｕｐによって提案されたフーリエ反復法がある（非特許文献４）。この手法においては、フーリエ領域の拘束条件と画像領域の拘束条件の両方を満足させるように、フーリエ変換および逆フーリエ変換を繰り返し適用し、画像を再構成する（図１参照）。フーリエ反復法における最も基本的なアルゴリズムは、エラーリダクション（ＥＲ：Ｅｒｒｏｒ  Ｒｅｄｕｃｔｉｏｎ）アルゴリズムである（非特許文献４）。【０００４】  フーリエ反復法の拘束条件の一つであるフーリエ領域の拘束条件には、測定された原画像のフーリエ振幅が用いられる。また、フーリエ反復法のもう一つの条件である画像領域の拘束条件には、原画像は非負の実数値画像であるとみなす非負条件と、原画像が占める領域以外では値が０であるという条件（サポート条件と呼ばれる）とが用いられる。サポート条件には、測定されるサポートが用いられる。ここで、「サポート」とは、原画像が占める領域のことである。また、「原画像」とは、観察対象画像のことである。
【非特許文献１】Ｒ．Ｗ．Ｇｅｒｃｈｂｅｒｇ，″Ｓｕｐｅｒ−ｒｅｓｏｌｕｔｉｏｎ  ｔｈｒｏｕｇｈ  ｅｒｒｏｒ  ｅｎｅｒｇｙ  ｒｅｄｕｃｔｉｏｎ，″Ｏｐｔ．Ａｃｔａ，ｖｏｌ．２１，ｐｐ．７０９−７２０，１９７４
【非特許文献２】Ｂ．Ｊ．Ｈｏｅｎｄｅｒｓ，″Ｏｎ  ｔｈｅ  ｓｏｌｕｔｉｏｎ  ｏｆ  ｔｈｅ  ｐｈａｓｅ  ｒｅｔｒｉｅｖａｌ  ｐｒｏｂｌｅｍ，″Ｊ．Ｍａｔｈ．Ｐｈｙｓ．，ｖｏｌ．１６，ｐｐ．１７１９−１７２５，１９７５
【非特許文献３】Ｒ．Ｅ．Ｂｕｒｇｅ，Ｍ．Ａ．Ｆｉｄｄｙ，Ａ．Ｈ．Ｇｒｅｅｎａｗａｙ，ａｎｄ  Ｇ．Ｒｏｓｓ，″Ｔｈｅ  ｐｈａｓｅ  ｐｒｏｂｌｅｍ，″Ｐｒｏｃ．Ｒ．Ｓｏｃ．Ｌｏｎｄ，ｖｏｌ．３５０，ｐｐ．１９１−２１２，１９７６
【非特許文献４】Ｊ．Ｒ．Ｆｉｅｎｕｐ，″Ｐｈａｓｅ  ｒｅｔｒｉｅｖａｌ  ａｌｇｏｒｉｔｈｍｓ：ａ  ｃｏｍｐａｒｉｓｏｎ，″Ａｐｐｌ．Ｏｐｔ．，ｖｏｌ．２１，ｐｐ．２７５８−２７６９，１９８２
【発明の開示】
【発明が解決しようとする課題】
【０００５】  しかしながら、従来のフーリエ反復法においては、次のような問題がある。すなわち、上記のように、画像領域の拘束条件としては、非負条件とサポート条件が用いられ、後者のサポート条件としては、測定されるサポートが用いられる。ところが、実際に測定されるサポートは、左右にずれたり、回転したり、膨張や収縮の影響を受けたりするなど、真のサポートからずれて測定される場合がある。そのため、後で詳述するように、場合によっては、再構成画像が原画像以外の画像に収束して停滞するという問題が生じる。なお、ここで、「真のサポート」とは、原画像の輪郭およびその内部から成る領域のことである。【０００６】  本発明は、かかる点に鑑みてなされたものであり、真のサポートが正確にわかっていない場合においても、位相回復アルゴリズムを用いて高精度に画像を回復することができる画像再構成方法を提供することを目的とする。
【課題を解決するための手段】
【０００７】  本発明は、原画像のフーリエ振幅を入力するステップと、原画像が占めるサポートを入力するステップと、入力したサポートに膨張処理を施すステップと、入力したフーリエ振幅および膨張処理を施したサポートを用いて、位相回復アルゴリズムおよびサポートを収縮させる機能を有する物体抽出アルゴリズムを併用してサポート条件を更新することにより、画像を再構成するステップとを有するようにした。
【発明の効果】
【０００８】  本発明によれば、真のサポートが正確にわかっていない場合においても、位相回復アルゴリズムを用いて高精度に画像を回復することができる画像再構成方法を得ることができる。
【図面の簡単な説明】
【０００９】
【図１】フーリエ反復法の原理を示すブロック図
【図２】原画像を示す図
【図３】図２の原画像のフーリエ振幅を示す図
【図４】初期画像を示す図
【図５】図４の初期画像を用いて得られた再構成画像を示す図
【図６】繰り返しに対するフーリエ誤差の変化を示すグラフ
【図７】図７Ａは、測定されたサポートが真のサポートから上に３０画素ずれた場合を示す図、図７Ｂは、測定されたサポートが真のサポートから下に３０画素ずれた場合を示す図、図７Ｃは、測定されたサポートが真のサポートから左に３０画素ずれた場合を示す図、図７Ｄは、測定されたサポートが真のサポートから右に３０画素ずれた場合を示す図
【図８】図８Ａは、測定されたサポートが真のサポートから上に３０画素ずれた場合におけるＥＲアルゴリズムによる再構成画像を示す図、図８Ｂは、その場合のフーリエ誤差の推移を示す図
【図９】図９Ａは、測定されたサポートが真のサポートから下に３０画素ずれた場合におけるＥＲアルゴリズムによる再構成画像を示す図、図９Ｂは、その場合のフーリエ誤差の推移を示す図
【図１０】図１０Ａは、測定されたサポートが真のサポートから左に３０画素ずれた場合におけるＥＲアルゴリズムによる再構成画像を示す図、図１０Ｂは、その場合のフーリエ誤差の推移を示す図
【図１１】図１１Ａは、測定されたサポートが真のサポートから右に３０画素ずれた場合におけるＥＲアルゴリズムによる再構成画像を示す図、図１１Ｂは、その場合のフーリエ誤差の推移を示す図
【図１２】図１２Ａは、測定されたサポートが真のサポートに対して５度回転した場合を示す図、図１２Ｂは、測定されたサポートが真のサポートに対して１５度回転した場合を示す図、図１２Ｃは、測定されたサポートが真のサポートに対して３０度回転した場合を示す図、図１２Ｄは、測定されたサポートが真のサポートに対して４５度回転した場合を示す図
【図１３】図１３Ａは、測定されたサポートが真のサポートに対して５度回転した場合におけるＥＲアルゴリズムによる再構成画像を示す図、図１３Ｂは、その場合のフーリエ誤差の推移を示す図
【図１４】図１４Ａは、測定されたサポートが真のサポートに対して１５度回転した場合におけるＥＲアルゴリズムによる再構成画像を示す図、図１４Ｂは、その場合のフーリエ誤差の推移を示す図
【図１５】図１５Ａは、測定されたサポートが真のサポートに対して３０度回転した場合におけるＥＲアルゴリズムによる再構成画像を示す図、図１５Ｂは、その場合のフーリエ誤差の推移を示す図
【図１６】図１６Ａは、測定されたサポートが真のサポートに対して４５度回転した場合におけるＥＲアルゴリズムによる再構成画像を示す図、図１６Ｂは、その場合のフーリエ誤差の推移を示す図
【図１７】真のサポートを示す図
【図１８】図１８Ａは、フィルタ窓３×３を用いて図１７の画像に膨張処理を施した場合を示す図、図１８Ｂは、フィルタ窓７×７を用いて図１７の画像に膨張処理を施した場合を示す図
【図１９】図１９Ａは、３×３のフィルタ窓を用いて膨張処理を施した場合におけるＥＲアルゴリズムによる再構成画像を示す図、図１９Ｂは、その場合のフーリエ誤差の推移を示す図
【図２０】図２０Ａは、７×７のフィルタ窓を用いて膨張処理を施した場合におけるＥＲアルゴリズムによる再構成画像を示す図、図２０Ｂは、その場合のフーリエ誤差の推移を示す図
【図２１】図２１Ａは、フィルタ窓３×３を用いて図１７の画像に収縮処理を施した場合を示す図、図２１Ｂは、フィルタ窓７×７を用いて図１７の画像に収縮処理を施した場合を示す図
【図２２】図２２Ａは、３×３のフィルタ窓を用いて収縮処理を施した場合におけるＥＲアルゴリズムによる再構成画像を示す図、図２２Ｂは、その場合のフーリエ誤差の推移を示す図
【図２３】図２３Ａは、７×７のフィルタ窓を用いて収縮処理を施した場合におけるＥＲアルゴリズムによる再構成画像を示す図、図２３Ｂは、その場合のフーリエ誤差の推移を示す図
【図２４】膨張および収縮が混在したサポートを含む画像を示す図
【図２５】図２５Ａは、膨張および収縮が混在したサポートを用いた場合におけるＥＲアルゴリズムによる再構成画像を示す図、図２５Ｂは、その場合のフーリエ誤差の推移を示す図
【図２６】大きいフィルタ窓（１７×１７）を用いて膨張処理を施した場合におけるＥＲアルゴリズムによる再構成画像を示す図
【図２７】スネークスにおける初期輪郭の設定の例を示す図
【図２８】本発明の実施の形態１に係る画像再構成方法のアルゴリズムの手順を示すフローチャート
【図２９】実施の形態１の方法を実現するための全体構成を示す概略図
【図３０】図３０Ａは、図２４のサポート条件を用いた場合においてＥＲアルゴリズムを適用して得られた再構成画像を示す図、図３０Ｂは、図２４のサポート条件を用いた場合において実施の形態１の方法を適用して得られた再構成画像を示す図
【図３１】図３１Ａは、四角のサポートを示す図、図３１Ｂは、図３１Ａのサポートを用いてＥＲアルゴリズムを適用して得られた再構成画像を示す図
【図３２】図３２Ａは、測定されたサポートが収縮の影響を受けている場合において実施の形態１の方法を適用して得られた再構成画像を示す図、図３２Ｂは、５０回後の出力画像を示す図、図３２Ｃは、１００回後の出力画像を示す図、図３２Ｄは、１５０回後の出力画像を示す図、図３２Ｅは、２００回後の出力画像を示す図、図３２Ｆは、２５０回後の出力画像を示す図
【図３３】図３３Ａは、測定されたサポートが膨張の影響を受けている場合において実施の形態１の方法を適用して得られた再構成画像を示す図、図３３Ｂは、５０回後の出力画像を示す図、図３３Ｃは、１００回後の出力画像を示す図、図３３Ｄは、１５０回後の出力画像を示す図、図３３Ｅは、２００回後の出力画像を示す図、図３３Ｆは、２５０回後の出力画像を示す図
【図３４】図３４Ａは、測定されたサポートが回転１度の影響を受けている場合において実施の形態１の方法を適用して得られた再構成画像を示す図、図３４Ｂは、５０回後の出力画像を示す図、図３４Ｃは、１００回後の出力画像を示す図、図３４Ｄは、１５０回後の出力画像を示す図、図３４Ｅは、２００回後の出力画像を示す図、図３４Ｆは、２５０回後の出力画像を示す図
【図３５】図３５Ａは、原画像を示す図、図３５Ｂは、収縮の影響を受けた測定サポートを示す図、図３５Ｃは、図３５Ｂのサポートを用いて従来のフーリエ反復法を適用して得られた再構成画像を示す図、図３５Ｄは、図３５Ｂのサポートを用いて実施の形態２の方法を適用して得られた再構成画像を示す図
【図３６】図３６Ａは、原画像を示す図、図３６は、回転の影響を受けた測定サポートを示す図、図３６Ｃは、図３６Ｂのサポートを用いて従来のフーリエ反復法を適用して得られた再構成画像を示す図、図３６Ｄは、図３６Ｂのサポートを用いて実施の形態２の方法を適用して得られた再構成画像を示す図
【図３７】電子顕微鏡の一般的な構成を示すブロック図
【図３８】本発明の実施の形態２に係る画像再構成方法のアルゴリズムの手順を示すフローチャート
【図３９】図３９Ａは、原画像を示す図、図３９Ｂは、初期画像を示す図、図３９Ｃは、図３９Ｂの初期画像から従来のフーリエ反復法を適用して得られた再構成画像を示す図、図３９Ｄは、図３９Ｂの初期画像から実施の形態２の方法を適用して得られた再構成画像を示す図
【図４０】本発明の実施の形態３に係る画像再構成方法のアルゴリズムの手順を示すフローチャート
【図４１】図４０のフーリエ振幅／サポート更新処理の内容を示すフローチャート
【図４２】連続する４つの時間フレームにおける測定サポートの一例を示す図
【図４３】連続する２つの時間フレームにおける動き情報の定義を説明するための図であり、特に、図４３Ａは、移動を示す図、図４３Ｂは、回転を示す図、図４３Ｃは、拡大・縮小を示す図、図４３Ｄは、移動を示す図
【図４４】図４４Ａは、真のサポートを示す図、図４４Ｂは、図４４Ａの真のサポートに対して＋７画素膨張している測定サポートを示す図、図４４Ｃは、図４４Ａの真のサポートに対して−７画素収縮している測定サポートを示す図
【図４５】図４５Ａは、真のサポートを示す図、図４５Ｂは、図４５Ａの真のサポートに対して４５度回転している測定サポートを示す図
【図４６】図４６Ａは、関数ＦＯＰＴが最適でない場合において当該関数ＦＯＰＴにより更新されたサポートを示す図、図４６Ｂは、図４６Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像のフーリエ誤差の推移を示す図
【図４７】図４７Ａは、関数ＦＯＰＴと関数Ｆ（ｎ）（ｎ＝１，２，‥，Ｎ）の（Ｎ＋１）個の関数のうち、関数ＦＯＰＴ以外の選択された最適な関数Ｆ（ｉ）により更新されたサポートを示す図、図４７Ｂは、図４７Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像のフーリエ誤差の推移を示す図
【図４８】図４８Ａは、関数ＦＯＰＴと関数Ｆ（ｎ）（ｎ＝１，２，‥，Ｎ）の（Ｎ＋１）個の関数のうち、関数ＦＯＰＴ最適な関数Ｆ（ｉ）以外のある関数により更新されたサポートを示す図、図４８Ｂは、図４８Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像のフーリエ誤差の推移を示す図
【図４９】図４９Ａは、測定結果のサポートを示す図、図４９Ｂは、図４９Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像のフーリエ誤差の推移を示す図
【図５０】第０フレームの処理を説明するための図であり、特に、図５０Ａは、測定された第０フレームのフーリエ振幅を示す図、図５０Ｂは、測定された第０フレームのサポートを示す図、図５０Ｃは、図５０Ａのフーリエ振幅および図５０Ｂのサポートに実施の形態１の方法を適用して得られた再構成画像を示す図
【図５１】第ｋフレームの処理を説明するための図であり、特に、図５１Ａは、測定された第ｋフレームのフーリエ振幅を示す図、図５１Ｂは、測定された第ｋフレームのサポートを示す図、図５１Ｃは、図５１Ａのフーリエ振幅および図５１Ｂのサポートに実施の形態３の方法を適用して得られた再構成画像を示す図
【図５２】図５２Ａは、フィルタ窓９×９を用いて図４４Ｃの測定サポートに膨張処理を施して得られた初期サポートを示す図、図５２Ｂは、図５２Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像を示す図
【図５３】図５３Ａは、フィルタ窓１１×１１を用いて図４４Ｃの測定サポートに膨張処理を施して得られた初期サポートを示す図、図５３Ｂは、図５３Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像を示す図
【図５４】図５４Ａは、フィルタ窓１３×１３を用いて図４４Ｃの測定サポートに膨張処理を施して得られた初期サポートを示す図、図５４Ｂは、図５４Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像を示す図
【図５５】図５５Ａは、フィルタ窓１７×１７を用いて図４４Ｃの測定サポートに膨張処理を施して得られた初期サポートを示す図、図５５Ｂは、図５５Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像を示す図
【図５６】図５６Ａは、フィルタ窓２５×２５を用いて図４４Ｃの測定サポートに膨張処理を施して得られた初期サポートを示す図、図５６Ｂは、図５６Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像を示す図
【図５７】図５７Ａは、フィルタ窓３１×３１を用いて図４４Ｃの測定サポートに膨張処理を施して得られた初期サポートを示す図、図５７Ｂは、図５７Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像を示す図
【図５８】図５８Ａは、フィルタ窓３３×３３を用いて図４４Ｃの測定サポートに膨張処理を施して得られた初期サポートを示す図、図５８Ｂは、図５８Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像を示す図
【発明を実施するための最良の形態】
【００１０】  以下、本発明の実施の形態について、図面を参照して詳細に説明する。【００１１】  （実施の形態１）
本発明者は、フーリエ反復法における問題点を検討し、例えば、ＥＲアルゴリズムによる再構成画像が、（１）回転、（２）収縮、（３）膨張と収縮の混在、の場合に収束しないことを見いだした。そして、この場合においても高精度な画像再構成を実現するフーリエ反復法を考案する必要があることを見いだした。さらに、そのためには、サポート条件を真のサポートに近づくように更新していくことが必要であることを見いだした。【００１２】  本発明は、測定されたサポートに膨張処理を施した後、測定されたフーリエ振幅および膨張処理を施したサポートを用いて、位相回復アルゴリズム（例えば、フーリエ反復法におけるＥＲアルゴリズム）およびサポートを収縮させる機能を有する物体抽出アルゴリズム（例えば、スネークス）を併用してサポート条件を更新することにより、画像を再構成する。【００１３】  まず、フーリエ反復法において最も基本的なアルゴリズムであるＥＲアルゴリズムについて説明する。【００１４】  図１は、フーリエ反復法の原理を示すブロック図である。ＥＲアルゴリズムは、フーリエ変換（図１の（ｉ））と逆フーリエ変換（図１の（ｉｉｉ））を通して、フーリエ領域の拘束条件（図１の（ｉｉ））と画像領域の拘束条件（図１の（ｉｖ））を繰り返し用い、画像を解に収束させるアルゴリズムである。ｍ（ｍ＝１，２，‥）回目の反復を例にとると、その処理手順は、以下のように表される。【００１５】  （ｉ）入力画像のフーリエ変換
入力画像ｇｍ−１（ｘ，ｙ）をフーリエ変換して、フーリエ係数Ｇｍ−１（ｕ，ｖ）を次の（式１）のように算出する。
【数１】
2005114568.tif 000002
ただし、θｍ−１はｍ−１回目の回復位相である。【００１６】  （ｉｉ）フーリエ領域の拘束条件の適用
（式１）で算出されたフーリエ振幅｜Ｇｍ−１（ｕ，ｖ）｜を測定された原画像のフーリエ振幅｜Ｆ（ｕ，ｖ）｜に置き換え、次の（式２）に示すようにフーリエ係数Ｇ’ｍ−１（ｕ，ｖ）を算出する。
【数２】
2005114568.tif 000003
【００１７】  （ｉｉｉ）逆フーリエ変換による出力画像の生成
（式２）で算出されたフーリエ係数Ｇ’ｍ−１（ｕ，ｖ）を逆フーリエ変換して、画像ｇ’ｍ−１（ｘ，ｙ）を次の（式３）により得る。
【数３】
2005114568.tif 000004
【００１８】  （ｉｖ）画像領域の拘束条件の適用
（式３）で得られたｇ’ｍ−１（ｘ，ｙ）に対して、次の（式４）に示す画像領域の拘束条件を適用し、ｍ回目の入力画像ｇｍ（ｘ，ｙ）を算出する。
【数４】
2005114568.tif 000005
ただし、Ｄはｇ’ｍ−１（ｘ，ｙ）の値が画像領域の拘束条件を満たすサンプル点の集合である。なお、画像領域の拘束条件とは、上記のように、原画像は非負の実数値画像であるとみなす非負条件と、原画像が占める領域以外では値が０であるというサポート条件のことである。【００１９】  そして、上記一連の処理を所定の繰り返し回数だけ行うか、またはエラー関数が所定の値以下となるまで繰り返される。また、初期入力画像ｇ０には、非負の乱数値を与えた画像が、一般に用いられる。【００２０】  次に、実際の画像を用いて、ＥＲアルゴリズムの性能を示す。図２は、原画像を示す図、図３は、図２の原画像のフーリエ振幅を示す図、図４は、初期画像を示す図、図５は、図４の初期画像からＥＲアルゴリズムにより再構成された画像を示す図である。なお、「原画像」とは、図２に示すように再構成されるべき画像であり、「初期画像」とは、図４に示すように最初に与えられるランダムな値を持つ画像である。【００２１】  ここで、フーリエ領域の拘束条件で用いられるフーリエ振幅としては、原画像（図２）から算出されたフーリエ振幅（図３）を用いる。また、画像領域の拘束条件であるサポート条件においては、サポートは正確にわかっているものとする。初期画像には、図４に示すように、非負の乱数値を与えた画像を用いる。この初期画像を用いてＥＲアルゴリズムを５００回繰り返し適用した場合の再構成画像を、図５に示す。図５から、再構成画像は正しい画像（図２の原画像）を再構成しているといえる。【００２２】  また、再構成画像を数値的に評価するため、次の（式５）を用いる。ここで、ＥＦｍはフーリエ誤差（Ｆｏｕｒｉｅｒ  ｅｒｒｏｒ）である。
【数５】
2005114568.tif 000006
図６は、フーリエ誤差の変化を示している。図６から、再構成画像は少ない繰り返し数（ｉｔｅｒａｔｉｏｎ）で正しい画像（原画像）に収束しているといえる。【００２３】  以上のことから、フーリエ反復法においては、正確なフーリエ領域の拘束条件および画像領域の拘束条件を用いれば、精度良く再構成画像を原画像に収束させることが可能であることがわかる。【００２４】  このように、フーリエ反復法は、フーリエ領域の拘束条件と画像領域の拘束条件の両方を満足する画像に収束させる手法であり、画像領域の拘束条件には、一般に、サポートが正確にわかっているものとしている。しかし、実際に測定されるサポートは、左右にずれたり、回転したり、膨張や収縮の影響を受けるなど、真のサポートからずれて測定される場合がある。そのため、場合によっては、再構成画像が原画像以外の画像に収束してしまうという問題が生じる。【００２５】  本発明者は、この点を明確にすべく、以下の検討を行った。すなわち、ここでは、実際の画像を用いて、測定されたサポートが真のサポートからずれて測定された場合のＥＲアルゴリズムへの影響を示す。ここでは、測定されたサポートが真のサポートからずれた場合として、次の３つのケースを想定する。（ケース１）測定されたサポートが真のサポートから上下左右にずれて測定された場合。（ケース２）測定されたサポートが真のサポートから回転して測定された場合。（ケース３）測定されたサポートが真のサポートと比較して膨張／収縮して測定された場合。ただし、フーリエ領域の拘束条件において、測定される振幅としては、原画像の振幅を用いる。また、ＥＲアルゴリズムの繰り返しは５００回とする。【００２６】  （ケース１）
ここでは、図７に示すように、測定されたサポート（白）が真のサポート（灰色）と比較してそれぞれ上下左右にずれた場合を考える。そして、このサポート条件を用いてＥＲアルゴリズムを適用した際に得られた再構成画像およびフーリエ誤差の推移を、図８〜図１１のＡ、Ｂ２つの図面にそれぞれ示す。なお、ここでは、実際の測定において生じるずれの最大値を３０画素と想定し、３０画素上下左右にずれた場合のサポートを用いた。すなわち、図７Ａおよび図８は、上に３０画素ずれた場合であり、図７Ｂおよび図９は、下に３０画素ずれた場合であり、図７Ｃおよび図１０は、左に３０画素ずれた場合であり、図７Ｄおよび図１１は、右に３０画素ずれた場合である。【００２７】  図８〜図１１のＡ、Ｂ２つの図面から、測定結果が上下左右にずれたサポートを用いても、再構成画像は正しい画像（図２）に収束することがわかる。【００２８】  （ケース２）
ここでは、図１２に示すように、測定されたサポート（白）が真のサポート（灰色）と比較して回転した場合を考える。そして、このサポート条件を用いてＥＲアルゴリズムを適用した際に得られた再構成画像およびフーリエ誤差の推移を、図１３〜図１６のＡ、Ｂ２つの図面にそれぞれ示す。なお、ここでは、実際の測定において生じる回転のずれの最大値を４５度と想定し、例として５度、１５度、３０度、４５度回転した場合のサポートの測定結果を用いた。すなわち、図１２Ａおよび図１３は、５度回転した場合であり、図１２Ｂおよび図１４は、１５度回転した場合であり、図１２Ｃおよび図１５は、３０度回転した場合であり、図１２Ｄおよび図１６は、４５度回転した場合である。【００２９】  図１３〜図１６のＡ、Ｂ２つの図面から、サポートの測定結果が真のサポートから回転した場合、再構成画像は正しい画像（図２）に収束しないことがわかる。【００３０】  （ケース３）
まず、真のサポート（図１７）に膨張処理を施した場合を図１８に示す。ここでは、実際の測定において生じる膨張の程度は、最大７×７のフィルタ窓程度と想定し、例としてフィルタ窓３×３、７×７を用いた（それぞれ図１８Ａ、図１８Ｂ）。また、図１８Ａおよび図１８Ｂをそれぞれ測定されたサポートとして用いてＥＲアルゴリズムを適用した場合の再構成画像およびフーリエ誤差の推移を、図１９、図２０のＡ、Ｂ２つの図面にそれぞれ示す。【００３１】  図１９〜図２０のＡ、Ｂ２つの図面から、７×７以下のフィルタ窓程度の膨張処理を施した場合のサポートを用いても、再構成画像は正しい画像（図２）に収束することがわかる。【００３２】  次に、真のサポート（図１７）に収縮処理を施した場合を図２１に示す。ここでも、上記のように、実際の測定において生じる収縮の程度は、最大７×７程度のフィルタ窓と想定し、例としてフィルタ窓３×３、７×７を用いた（それぞれ図２１Ａ、図２１Ｂ）。また、図２１Ａおよび図２１Ｂをそれぞれ測定されたサポートとして用いてＥＲアルゴリズムを適用した場合の再構成画像およびフーリエ誤差の推移を、図２２、図２３のＡ、Ｂ２つの図面にそれぞれ示す。【００３３】  図２２〜図２３のＡ、Ｂ２つの図面から、収縮処理を施した場合のサポートを用いた場合、再構成画像は正しい画像（図２）に収束しないことがわかる。【００３４】  ところで、実際の測定では、上記（ケース１）〜（ケース３）が混在するサポートが測定される。そこで、ここでは、例として真のサポート（図１７）に図２４に示すように膨張および収縮が混在したサポートを示す画像を用いて、ＥＲアルゴリズムを適用した際の再構成画像およびフーリエ誤差の推移を、図２５Ａおよび図２５Ｂにそれぞれに示す。【００３５】  図２５Ａおよび図２５Ｂから、膨張および収縮が混在したサポートを用いた場合、収縮の影響により、再構成画像は正しい画像（図２）に収束しないことがわかる。【００３６】  上記の実験から、（１）おおよその輪郭（真実に近いもの）が得られているが、そのサポートが回転した形で与えられた場合、（２）おおよその輪郭（真実に近いもの）が得られているが、そのサポートが収縮された形で与えられた場合、（３）膨張と収縮の混在の影響を受けたサポートが与えられた場合、において、ＥＲアルゴリズムによる再構成画像は原画像に収束せず停滞してしまうことがわかった。これは、本発明者が実験により得た新しい知見である。【００３７】  そこで、本発明者は、上記新しい知見に基づき、サポートが正確にわかっていない場合、特に、サポートが回転、収縮、膨張と収縮の混在した場合であっても、高精度に画像を再構成することができる新しいフーリエ反復法を考案した。【００３８】  すなわち、上記のように、測定されたサポートが膨張（または上下左右のずれ）のみの影響を受けている場合、再構成画像は正しい画像に収束する。そこで、本実施の形態では、サポートの測定結果に膨張処理を施すことにより、たとえ回転、収縮、膨張と収縮の混在したサポートであったとしても、真のサポートを含むように新たなサポートを生成する。このとき、膨張されている部分もさらに膨張されるため、図２６（フィルタ窓１７×１７の場合���に示すように、膨張処理のフィルタ窓が大きすぎる場合、再構成画像が正しい画像に収束しないという問題が生じる。そこで、本実施の形態では、サポート条件を真のサポートに近づくように更新していくことで、この問題を回避する。サポート条件の更新には、ＥＲアルゴリズムを数回繰り返すことにより得られた出力画像に「スネークス（Ｓｎａｋｅｓ）」と呼ばれるアルゴリズムを用いて画像を抽出し、その抽出画像を新たなサポート条件とする。このようにスネークスを用いることで、膨張した画像を真のサポートに近づけることが可能になる。したがって、上記の処理を繰り返し行うことにより、高精度な再構成画像を実現することが可能になる。【００３９】  要するに、本実施の形態では、スネークスを適用する前に膨張処理を施すとともに、ＥＲアルゴリズムとスネークスを併用してサポート条件の更新を行うようにしている。【００４０】  ここで、スネークスとは、対象物体抽出手法の１つであり、動的輪郭モデルとも呼ばれている。スネークスの詳細は、非特許文献５（Ｍ．Ｋａｓｓ，Ａ．Ｗｉｔｋｉｎ，ａｎｄ  Ｄ．Ｔｅｒｚｏｐｏｕｌｏｓ，“Ｓｎａｋｅｓ：Ａｃｔｉｖｅ  ｃｏｎｔｏｕｒ  ｍｏｄｅｌｓ，”Ｉｎｔ．Ｊ．Ｃｏｍｐｕｔ．Ｖｉｓ．，ｖｏｌ．１，ｐｐ．３１２−３３１，１９８８）に示されている。スネークスは、真のサポートが得られていない場合に、任意（適宜）に初期輪郭を設定し、それを、スネークスの理論で定義されているエネルギー関数に基づいて収縮させることにより、真のサポートを提示するという機能を有する。より具体的には、ＥＲアルゴリズムである程度回復された画像にスネークスを適用し、それを新たなサポートとすることで、真のサポートからずれて測定されたサポート（回転、収縮、膨張と収縮の混在）を徐々に真のサポートに近づけることができる。【００４１】  なお、スネークスの初期輪郭は、抽出対象物体が含まれるように設定される。図２７は、抽出対象物体である目標（Ｔａｒｇｅｔ）に対し、スネークスにおける初期輪郭（Ｉｎｉｔｉａｌ  ｃｏｎｔｏｕｒ）の設定の例を示している。【００４２】  また、エネルギー関数は、濃淡画像における画素の濃淡値を高度とみなして３次元表現し、観察対象を包含するように弾性体を配置することに相当する。弾性体は、（力を取り去ると）収縮し、濃淡値差が大きいエッジ部分でその収縮を停止し、エッジ抽出が可能になる。「スネークス」という用語は、上記のアルゴリズムを指す場合と、上記の弾性体を指す場合の両方がある。そこで、ここでは、前者を「スネークス（Ｓｎａｋｅｓ）アルゴリズム」とも表記する。【００４３】  一方、ＥＲアルゴリズムは、上記のように、フーリエ振幅とサポートを用いて、失われた位相を回復することにより、原画像を再構成するという機能を有する。【００４４】  ここで、スネークスアルゴリズムについて説明する。スネークスは、画像平面（ｘ，ｙ）上にｓを媒介変数として表現された閉曲線ｖ（ｓ）＝（ｘ（ｓ），ｙ（ｓ））（０≦ｓ≦１）であり、次の（式６）で定義されるエネルギー関数（「スネークスのエネルギー」と呼ばれる）を最小化するように変形する。
【数６】
2005114568.tif 000007
本実施の形態では、離散画像を対象とするため、ｖ（ｓ）を離散点ｖｉ＝（ｘｉ，ｙｉ）（ｉ＝１，２，‥，Ｎ）に置き換えて表現する。したがって、（式６）は次の（式７）で書き換えられる。
【数７】
2005114568.tif 000008
（式７）において、Ｅｉｎｔ（ｖｉ）は内部エネルギーと呼ばれ、輪郭モデルの滑らかさを表し、次の（式８）で定義される。
【数８】
2005114568.tif 000009
ただし、αｉ、βｉは重み係数を表す。このＥｉｎｔ（ｖｉ）の働きにより、スネークスは滑らかに収縮する力を受ける。また、Ｅｉｍａｇｅ（ｖｉ）は画像エネルギーと呼ばれ、しばしば、次の（式９）で定義されるＥｅｄｇｅ（ｖｉ）が用いられる。
【数９】
2005114568.tif 000010
ただし、ｗｅｄｇｅは重み係数、Ｉ（ｖｉ）は画像の位置ｖｉにおける輝度値を表す。このＥｅｄｇｅ（ｖｉ）の働きにより、スネークスは画像の特徴であるエッジなどに引き寄せられる力を受ける。さらに、Ｅｃｏｎ（ｖｉ）は外部エネルギーと呼ばれ、画像の特徴に注目して、多様な定義が可能であることが非特許文献５に記載されている。そこで、本実施の形態では、凹形状の画像を抽出対象にしているため、凹形状の輪郭の抽出が可能になるように、外部エネルギーとして、非特許文献６（坂口  俊文，大山  公一，″面積項をもつスネーク，″信学春季全大，Ｄ−５５５，１９９１）において次の（式１０）で定義されたＥａｒｅａを用いる。
【数１０】
2005114568.tif 000011
（式１０）のＥａｒｅａは面積項と呼ばれ、ｖｉ（ｉ＝１，２，‥，Ｎ）が囲む面積を表している。この面積項をスネークスのエネルギーに加えることで、輪郭モデルによって囲まれた面積に比例したエネルギーが増加し、エネルギーの最小化の過程で面積を減少させる力が働き、凹形状の輪郭の抽出が可能になる。【００４５】  図２８は、本実施の形態に係るスネークスを用いたフーリエ反復法のアルゴリズムの手順を示すフローチャートである。【００４６】  まず、ステップＳ１０００では、測定された原画像のフーリエ振幅｜Ｆ（ｕ，ｖ）｜を入力する。【００４７】  そして、ステップＳ１１００では、測定されたサポートＤを入力する。【００４８】  そして、ステップＳ１２００では、初期画像Ｉ（ｘ，ｙ）を入力する。【００４９】  なお、ステップＳ１０００〜ステップ１２００の順序は、これに限定されず、自由に変更可能である。【００５０】  そして、ステップＳ１３００では、ステップＳ１２００で入力した初期画像Ｉ（ｘ，ｙ）をｇｉｎｉｔｉａｌ（ｘ，ｙ）とする。ここで、ｇｉｎｉｔｉａｌ（ｘ，ｙ）は、ステップＳ１９００の処理の入力画像である。【００５１】  そして、ステップＳ１４００では、ステップＳ１１００で入力したサポートＤに膨張処理を施し、この膨張処理を施した結果を新たなサポートＤとする。ここで、膨張処理は、当該膨張処理後の新たなサポートが真のサポートを含むように施される。【００５２】  そして、ステップＳ１５００では、サポート以外では値が０である画像ｄ（ｘ，ｙ）を生成する。【００５３】  そして、ステップＳ１６００では、カウンタｎを１に設定する。【００５４】  そして、ステップＳ１７００では、現在の画像ｄ（ｘ，ｙ）にスネークスアルゴリズムを適用し、物体Ｄ’を抽出する。【００５５】  そして、ステップＳ１８００では、ステップＳ１７００で抽出した物体Ｄ’を新たなサポートＤとする。【００５６】  そして、ステップＳ１９００では、ステップＳ１０００で入力したフーリエ振幅｜Ｆ（ｕ，ｖ）｜と現在のサポートＤを用いて、ステップＳ１３００で得られたｇｉｎｉｔｉａｌ（ｘ，ｙ）の画像に対して、Ｍ回ＥＲアルゴリズムを適用し、出力画像ｇｎ（ｘ，ｙ）を得る。【００５７】  そして、ステップＳ２０００では、ステップＳ１９００で得た出力画像ｇｎ（ｘ，ｙ）を新たな画像ｇｉｎｉｔｉａｌ（ｘ，ｙ）とする。【００５８】  そして、ステップＳ２１００では、ステップＳ１９００で得た出力画像ｇｎ（ｘ，ｙ）を新たな画像ｄ（ｘ，ｙ）とする。【００５９】  そして、ステップＳ２２００では、カウンタｎを１だけインクリメントする。【００６０】  そして、ステップＳ２３００では、カウンタｎの値があらかじめ設定された繰り返し回数（Ｎ）よりも大きいか否かを判断する。この判断の結果として、ｎ≦Ｎの場合は（Ｓ２３００：ＮＯ）、ステップＳ１７００に戻り、ｎ＞Ｎの場合は（Ｓ２３００：ＹＥＳ）、ステップＳ２４００に進む。【００６１】  ステップＳ２４００では、ステップＳ１７００〜ステップＳ２１００の処理をＮ回繰り返して得られる出力画像ｇＮ（ｘ，ｙ）を、本手法により最終的に得られる再構成画像として出力し、上記一連の処理を終了する。【００６２】  要するに、本手法では、手順１として、初期画像Ｉ（ｘ，ｙ）をｇｉｎｉｔｉａｌ（ｘ，ｙ）とする。そして、手順２として、測定されたサポートに膨張処理を施す。そして、手順３として、サポートＤを表す画像ｄ（ｘ，ｙ）を生成する。そして、手順４として、画像ｄ（ｘ，ｙ）にスネークスアルゴリズムを適用し、物体Ｄ’を抽出する。そして、手順５として、手順４で抽出した物体Ｄ’を新たなサポートＤとする。そして、手順６として、手順５で得られたサポートＤおよび測定された原画像のフーリエ振幅｜Ｆ（ｕ，ｖ）｜を用いて、ｇｉｎｉｔｉａｌ（ｘ，ｙ）に対して、Ｍ回ＥＲアルゴリズムを適用し、出力画像ｇｎ（ｘ，ｙ）を得る。そして、手順７として、得られた出力画像ｇｎ（ｘ，ｙ）を画像ｇｉｎｉｔｉａｌ（ｘ，ｙ）に入力する。そして、手順８として、得られた出力画像ｇｎ（ｘ，ｙ）を画像ｄ（ｘ，ｙ）に入力する。そして、手順９として、手順４〜手順８をあらかじめ設定された繰り返し回数（Ｎ回）実行することにより、画像を再構成する。Ｎ回後に得られる出力画像ｇＮ（ｘ，ｙ）が、本手法により最終的に得られる再構成画像である。【００６３】  なお、本実施の形態では、ＥＲアルゴリズムとスネークスアルゴリズムの併用に関して、適用の順番として、スネークスアルゴリズムから開始し、ＥＲアルゴリズムで終了するようにしているが、これに限定されるわけでない。すなわち、必ずしもスネークスアルゴリズムから開始する必要はなく、また、ＥＲアルゴリズムで終了する必要もない。ただし、以下の理由から、本実施の形態では、スネークスアルゴリズムから開始し、ＥＲアルゴリズムで終了するようにしている。【００６４】  まず、スネークスアルゴリズムから開始する理由は、次のとおりである。ＥＲアルゴリズムでは、サポート以外の部分は画像が回復されないため、スネークスアルゴリズムを適用する前に、設定する初期スネークスが、完全に真のサポートを包含するする必要がある。したがって、真のサポートに収縮を施したものがサポートとして与えられているおそれがある場合は、それに膨張処理を施した後にスネークスアルゴリズムを適用するのが好ましい。そのため、収縮の影響を受けたサポートの場合は、はじめに膨張処理を施した後にスネークスアルゴリズムから開始する。ただし、与えられたサポートが真のサポートを完全に包含している場合、例えば、回転や膨張の場合は、必ずしもスネークスアルゴリズムから開始する必要はない。【００６５】  また、ＥＲアルゴリズムで終了する理由は、次のとおりである。例えば、スネークスアルゴリズムで終了するとした場合、スネークスアルゴリズムがｍ回目に終了したとすると、その結果は、ｍ＋１回目のＥＲアルゴリズムの適用の際に用いるサポートとなる。そのため、Ｍ回目でＥＲアルゴリズムを終了する場合、Ｍ回目のＥＲアルゴリズムの後にスネークスアルゴリズムを適用してサポートを更新しても、その更新されたサポートは使用されないため、その処理が無駄になってしまう。そのため、本実施の形態では、スネークスアルゴリズムではなく、ＥＲアルゴリズムで終了するようにしている。【００６６】  図２９は、本手法を実現するための全体構成を示す概略図である。【００６７】  まず、ステッ��Ｓ３０００では、原画像のフーリエ振幅を測定する。この測定は、ソース（入射波の源）からコヒーレントな波を試料に入射し、散乱された波を検出器によって回折像として記録することによって行われる。【００６８】  そして、ステップＳ３１００では、サポートの測定を行う。この測定は、物理的なレンズを用いて検出器によって粗い実像を取得することによって行われる。なお、ステップＳ３０００とステップＳ３１００の順序は、もちろん、逆であってもよい。【００６９】  そして、ステップＳ３２００では、再構成画像の生成を行う。具体的には、計算機において、ステップＳ３０００で測定したフーリエ振幅とステップＳ３１００で測定したサポートを用いて、図２８に示す本手法（アルゴリズム）を適用して再構成画像ｇＮを取得する。【００７０】  次に、実験例を示す。ここでは、実際の画像を用いて、測定されたサポートが真のサポートからずれた場合の本手法の有効性を示す。【００７１】  条件は、次のとおりである。画像領域の拘束条件であるサポート条件において、測定されたサポートとして、真のサポート（図１７）に膨張および収縮が混在したサポートを示す画像（図２４）を用いる。また、フーリエ領域の拘束条件において、測定される振幅として原画像の振幅（図３）を用いる。【００７２】  上記の条件を用いて、初期画像（図４）にＥＲアルゴリズムおよび本手法を適用して再構成した画像を、それぞれ、図３０Ａおよび図３０Ｂに示す。すなわち、図３０Ａは、ＥＲアルゴリズムを適用して得られた再構成画像を示し、図３０Ｂは、本手法を適用して得られた再構成画像を示している。図３０Ａからは、ＥＲアルゴリズムを用いた場合、再構成画像は原画像（図２）に収束していないことがわかる。一方、図３０Ｂからは、本手法を用いた場合、再構成画像は原画像（図２）に収束していること、つまり、本手法は、サポートが正確にわかっていない場合においても、高精度に画像を回復することが可能であることがわかる。【００７３】  図３１は、サポートを四角にした場合のＥＲアルゴリズムの結果を示している。すなわち、図３１Ａは、四角のサポートを示し、図３１Ｂは、図３１Ａのサポートと原画像の振幅（図３）を用いて、ＥＲアルゴリズムを適用して得られた再構成画像を示している。ここで、原画像は、図２に示す画像である。この場合、再構成画像は原画像（図２）に収束していないことがわかる。【００７４】  図３２〜図３４は、それぞれ、収縮、膨張、回転の影響を受けたサポートを用いた場合の本手法による再構成画像を示している。すなわち、図３２は、測定されたサポートが収束の影響を受けている場合の本手法における再構成画像であり、図３３は、測定されたサポートが膨張の影響を受けている場合の本手法における再構成画像であり、図３４は、測定されたサポートが回転１度の影響を受けている場合の本手法における再構成画像である。スネークスに基づく本手法で用いるパラメータの値は、次のとおりである。本実験では、スネークスアルゴリズムの（式８）および（式９）における３つのパラメータαｉ、βｉ、ｗｅｄｇｅは、それぞれ、αｉ＝０．４、βｉ＝０．４、ｗｅｄｇｅ＝１．５に設定した。このスネークスアルゴリズムで用いら���る３つのパラメータは、抽出対象物体の形状によっては、異なる値に設定される。また、図２８のステップＳ１９００におけるＭ（ＥＲアルゴリズムの適用回数）は、Ｍ＝５０に、ステップＳ２３００におけるＮ（手順の繰り返し回数）は、Ｎ＝５にそれぞれ設定した。図３２〜図３４からは、本手法を用いた場合、再構成画像は原画像（図２）に収束していることがわかる。なお、図３２〜図３４のＡの各図面に示す初期画像において、白で示す画像は、測定されたサポートに膨張処理を施して得られた新たなサポートである。【００７５】  次に、本手法と従来のフーリエ反復法との比較例を示す。【００７６】  図３５は、測定されたサポートが収縮の影響を受けている場合である。すなわち、図３５Ａは、原画像を示す図、図３５Ｂは、収縮の影響を受けた測定サポートを示す図、図３５Ｃは、図３５Ｂのサポートを用いて従来のフーリエ反復法を適用して得られた再構成画像を示す図、図３５Ｄは、図３５Ｂのサポートを用いて本手法を適用して得られた再構成画像を示す図である。従来のフーリエ反復法を用いた場合は、図３５Ｃに示すように、再構成画像は原画像（図３５Ａ）に収束していないが、本手法を用いた場合は、図３５Ｄに示すように、再構成画像は原画像（図３５Ａ）に収束している。【００７７】  また、図３６は、測定されたサポートが回転の影響を受けている場合である。すなわち、図３６Ａは、原画像を示す図、図３６Ｂは、回転の影響を受けた測定サポートを示す図、図３６Ｃは、図３６Ｂのサポートを用いて従来のフーリエ反復法を適用して得られた再構成画像を示す図、図３６Ｄは、図３６Ｂのサポートを用いて本手法を適用して得られた再構成画像を示す図である。従来のフーリエ反復法を用いた場合は、図３６Ｃに示すように、再構成画像は原画像（図３６Ａ）に収束していないが、本手法を用いた場合は、図３６Ｄに示すように、再構成画像は原画像（図３６Ａ）に収束している。【００７８】  このように、本実施の形態によれば、測定されたサポートに対し、スネークスアルゴリズムを適用する前に膨張処理を施し、ＥＲアルゴリズムとスネークスアルゴリズムを併用してサポート条件の更新を行うため、サポートが正確にわかっていない場合においても、高精度に画像を回復することができる。【００７９】  なお、本実施の形態では、スネークスを適用するが、この理由は次のとおりである。膨張の影響を大きく受けたサポートを用いてＥＲアルゴリズムを実行すると、図２６に示すように正しい原画像に収束しない。一方、真のサポートを用いた場合、図５に示すように、ＥＲアルゴリズムによる再構成画像は原画像に収束する。実験から、膨張処理を施したサポートを用いて繰り返しＥＲアルゴリズムを適用すると、再構成画像は徐々に真のサポートに近づいていくことがわかった。ただし、膨張の程度が大きい場合は、真のサポートの収束の速度は遅く、場合によっては、完全な収束を見ることができなかった。そこで、ある回数ごとに、再構成した画像をスネークアルゴリズムで抽出し、サポートを収縮し、徐々に真のサポートに近づけるようにした。その結果、徐々に再構成画像が原画像に収束する。なお、スネークスはエッジが明確でない対象物を抽出する手法として有効であることが知られているため、完全に収束していない再構成画像のエッジをおおよそ抽出し、サポートを収縮させることに有効であると考えた。【００８０】  したがって、「収縮」が可能であれば、スネークス以外の方法でも応用は可能である。具体的には、例えば、次の３つが考えられる。【００８１】  第１は、エッジ抽出アルゴリズムと領域抽出などの画像処理との組み合わせである。すなわち、エッジを抽出しておき、そのエッジで構成される領域を選び出し、選ばれた領域（またはその膨張領域）をサポートにする処理である。【００８２】  第２は、テクスチャ抽出を用いたサポートの抽出である。テクスチャ抽出手法とは、例えば、画像の中から、格子模様の領域を選び出したり、または、水玉模様の領域を選び出したりする方法である。サポートが未知の場合または誤差を含む領域が与えられている場合には、ＥＲアルゴリズムの再構成画像は真実の画像ではなく、観察対象物体の位置に縞模様が現れる。したがって、テクスチャ抽出の一般的な方法でも適用は可能である。【００８３】  第３は、最適化手法である。例えば、曖昧なサポート内部に物体が存在していると仮定して、最適化手法であるＧＡ（遺伝的アルゴリズム：Ｇｅｎｅｔｉｃ  Ａｌｇｏｒｉｔｈｍ）で考えられるすべての物体形状を探索範囲として、最適な物を求める方法である。【００８４】  また、本発明の適用分野は、原画像のフーリエ振幅と原画像が占める領域（サポート）とが既知情報として得られる装置であれば、どのような分野でも適用可能である。例えば、一般に電子顕微鏡や天体観測、Ｘ線結晶学では、上記２つの情報が測定によって得られるため、適用が可能である。装置の全体構成としては、図２９に示すように、まず原画像のフーリエ振幅とサポートを所定の装置で測定し、この２つの情報を用いて計算機で本発明の位相反復アルゴリズムを用いて再構成画像を取得する。【００８５】  図３７は、電子顕微鏡の一般的な構成を示すブロック図である。この電子顕微鏡１００において、電子源１１０で発生した電子線は、平行照射用レンズ系１２０によって平行な電子ビームとなり、試料１３０に照射され、散乱される。散乱された回折波は、対物レンズ１４０を通して検出器１５０によって記録される。ここでは、原画像のフーリエ振幅と原画像が占める領域（サポート）は、検出器１５０によって測定される。電子源１１０、平行照射用レンズ系１２０、試料１３０、対物レンズ１４０、および検出器１５０は、真空の筐体１６０に収納されている。【００８６】  検出器１５０はコンピュータ１７０に接続されている。コンピュータ１７０は、ソフトウエアとして本発明の位相反復アルゴリズムを格納しており、検出器１５０で測定される原画像のフーリエ振幅およびサポートを用いて、再構成画像を生成する。【００８７】  上記のように、本発明は、原画像のフーリエ振幅と原画像が占める領域（サポート）が既知情報として得られていれば、適用分野を問わず、そのまま適用可能である。ただし、本発明は、サポートに測定誤差が生じた場合を考慮しているが、測定されたフーリエ振幅にノイズが生じた場合は、何らかの方法でノイズを除去する必要がある。この問題を解決する手法としては、例えば、非特許文献７（高城  洋明，高橋  徹，野津  剛，″ノイズを含んだフーリエ振幅からの物体回復アルゴリズム，″光学，ｖｏｌ．２５，ｐｐ．２２３−２２７，１９９６）や非特許文献８（高城  洋明，伊藤  克彦，藤崎  敏朗，″位相回復に用いるフーリエ反復法の改善，″光学，ｖｏｌ．３２，ｐｐ．３９−４５，２００３）がある。したがって、測定されたフーリエ振幅にノイズが生じた場合は、ＥＲアルゴリズムの部分をそれらの手法に置き換えることで解決される。【００８８】  （実施の形態２）
上記ＥＲアルゴリズムは、必要とする計算量が少なく、２次元物体にも容易に適用できるなど、多くの長所を有している。しかし、ＥＲアルゴリズムは、フーリエ振幅とサポートが正確にわかっていたとしても、初期画像によっては、正しい画像（解）に到達せず、それ以外の画像に停滞してしまうことが知られている。非特許文献９（高城  洋明，高橋  徹，寺鳥  誠範，永野  直広，″繰り返し型位相回復アルゴリズムにおける停滞問題の数値シミュレーションによる考察，″光学，ｖｏｌ．２１，ｐｐ．１１９−１２７，１９９２）参照。【００８９】  非特許文献９によれば、ＥＲアルゴリズムは、２×２Ｌ型物体の場合において、原画像と初期画像の３つの最高周波数成分の位相が一致しない場合、完全なる停滞または不完全なる停滞を引き起こすという性質を有する。ここで、完全なる停滞とは、原画像と異なる収束画像を有する場合であり、不完全なる停滞とは、途中で停滞を起こすが、いずれ原画像に収束する場合である。具体的には、２×２Ｌ型物体の場合、任意の初期画像が３つの最高周波数成分の位相の違いによって４つの領域に分割され、そのうち、原画像に収束できる領域は、完全なる停滞の場合は１つ、不完全なる停滞の場合は２つ存在する。また、４つの領域を構成する３つの境界は、任意の初期画像における３つの最高周波数成分のうちの１つのフーリエ振幅がゼロ（０）となる画像の位置に生じる。【００９０】  本発明者は、２×２Ｌ型物体におけるＥＲアルゴリズムの上記性質を検討した結果、完全なる停滞を防ぎ、再構成画像を原画像に収束させるためには、任意の初期画像の特定の周波数（ここで、特定の周波数とは、例えば、２×２Ｌ型物体の場合、３つの最高周波数成分のことである）から生じる４つの領域から、原画像に収束する領域を見つける必要があることを見いだした。すなわち、この領域は、３つの最高周波数成分の位相から構成されているため、再構成画像を原画像に収束させるためには、初期画像の最高周波数成分の位相の組み合わせを原画像の位相の組み合わせに一致させる必要があることを見いだした。【００９１】  また、非特許文献９では、２×２Ｌ型の小規模物体を例として数値シミュレーションを行っており、再構成画像内部の濃淡値が複雑な物体（以下「大規模物体」という）において上記停滞問題の原因が完全に解決されたとは言いがたいことを見いだした。本発明では、大規模物体も対象にしている。大規模物体の場合には、３つの最高周波数成分のフーリエ振幅がゼロとなる物体の位置に生じる境界の他に、ある特定の周波数のフーリエ振幅がゼロとなり、その位置に境界が多数生じるため、任意の初期画像から生じる領域が増え、それに伴って、ある特定の周波数から構成される領域の位相の組み合わせも膨大な数になる。膨大な位相の組み合わせからなる領域から、原画像が存在する領域を見つけるのは、困難である。【００９２】  本発明は、ある特定の周波数によって構成される領域の位相に注目し、とり得る位相の組み合わせの中から最適化手法を用いて、任意の初期画像から特定の周波数において構成される領域の位相の組み合わせとそれに対応する周波数における原画像の位相の組み合わせとを一致させるものである。本実施の形態では、���適化手法として遺伝的アルゴリズム（ＧＡ：Ｇｅｎｅｔｉｃ  Ａｌｇｏｒｉｔｈｍ）を用いる。【００９３】  また、遺伝的アルゴリズムにおける探索を容易にするため、遺伝的アルゴリズムを適用する前に、フーリエ振幅｜Ｆ｜にローパスフィルタ処理を施し、フーリエ振幅の要素数を減らす。すなわち、フーリエ振幅の要素数を減らすことにより、｜Ｆ｜のサイズを小さくすることができ、遺伝的アルゴリズムによって探索すべき位相の組み合わせを少なくすることができる。これにより、大規模物体においても、短い計算時間で正しい画像を得ることが可能になる。【００９４】  ここで、遺伝的アルゴリズムは、組み合わせ最適化問題の最適解を探索する手法であり、生物集団が自然淘汰を繰り返しながら進化していく現象をモデル化したものである。すなわち、まず、固体つまり解候補を複数個生成し、それら固体の環境への適応度に応じて、選択淘汰、交叉、突然変異といった遺伝的操作を繰り返していくことにより、最適解を求める手法である。遺伝的アルゴリズムは、探索空間が広く、従来の手法では探索が困難であった最適化問題に対し実用解を与えること、解を評価するための評価関数が柔軟に設定できること、探索過程において評価関数の微分値を必要としないことなどの長所を有する。【００９５】  図３８は、本実施の形態に係る遺伝的アルゴリズムの手順を示すフローチャートである。【００９６】  まず、ステップＳ４０００では、遺伝子型の決定を行う。すなわち、遺伝的アルゴリズムの課題となる事象のモデル化を行う。具体的には、対象とする問題の解候補を「０」、「１」の決まった長さの文字列にコード化する。各文字列が１つの個体に相当する。本手法における各個体は、「０」または「π」からなる領域を構成する位相を表し、遺伝子型はこの位相を表現する「０」または「１」の文字列を遺伝子として持つ。【００９７】  そして、ステップＳ４１００では、初期集団の生成を行う。すなわち、ステップＳ４０００で決定した遺伝子型に従って、複数の生物個体を生成する。具体的には、とり得る解候補をランダムに生成する。【００９８】  そして、ステップＳ４２００では、固体の評価を行う。すなわち、適応環境における適応度の算定を行う。具体的には、各固体を、次の評価関数により評価する。ここで、適応度の評価関数としては、次の２つが考えられる。【００９９】  第１は、フーリエ誤差を示す上記（式５）を評価関数とする場合である。この場合、評価関数の値が小さい個体ほど、適応度は高い。この評価関数を使う場合、ランク戦略を使って個体を選択する（個体の選択はステップＳ４３００で行う）。ここで、ランク戦略とは、適応度によって各個体をランク付けし、あらかじめ各ランクに対して決められた確率で子孫を残すようにするものである。各個体の選択確率は適応度によらずにランクに依存する。【０１００】  第２は、フーリエ誤差の逆数、つまり、次の（式１１）を評価関数とする場合である。
【数１１】
2005114568.tif 000012
この場合、原画像においてある繰り返し回数ｍのときのＥＦｍがゼロになるため、評価関数が無限大になる。この場合の解決法としては、２つある。１つ目は、無限大になるのを防ぐため、ＥＦｍが０に限りなく近い値となるように、微小な正の値を加算することである。２つ目は、ランク戦略を使って、数字を割り当て、個体を選択することである（個体の選択はステップＳ４３００で行う）。ランク戦略を使えば、同じ評価値の場合は同一順序でよいため、評価関数の値が無限大になっても問題はない。【０１０１】  そして、ステップＳ４３００では、遺伝子的操作を行う。遺伝子的操作は、選択・淘汰、交叉、および突然変異からなる。選択・淘汰では、各個体の評価に基づいて個体の選択・淘汰を行う。交叉では、遺伝子型として定義されている遺伝子情報を個体間で入れ替え、新しい個体を生成する。突然変異では、遺伝子型として定義されている遺伝子情報の記号列を、特定の確率で変化させ、新しい個体を生成する。突然変異は、一定の確率で個体の遺伝子を「０」から「１」に、または、「１」から「０」に置き換えることにより遂行される。【０１０２】  そして、ステップＳ４４００では、終了条件を満たすか否かを判断する。ここで、終了条件とは、探索を終了する条件のことであって、具体的には、例えば、総世代交代数があらかじめ設定された値に達したとき、または、適応度に基づいた条件を設定しその条件を満たしたときに、探索を終了する。この判断の結果として、終了条件を満たしている場合は（Ｓ４４００：ＹＥＳ）、上記一連の処理を終了し、終了条件を満たしていない場合は（Ｓ４４００：ＮＯ）、ステップＳ４２００に戻って、ステップＳ４２００〜ステップＳ４３００の処理を繰り返す。【０１０３】  以下、具体例を用いて説明する。まず、原画像として次の（式１２）に示す２×２の小規模物体でＬ型のサポートを含む画像ｆ（４×４のサイズ）を用いる。
【数１２】
2005114568.tif 000013
ただし、ａ、ｂ、ｃは正の実数である。次に、初期画像ｇ０を次の（式１３）とし、
【数１３】
2005114568.tif 000014
ｍ回目の回復画像ｇｍを次の（式１４）で表す。
【数１４】
2005114568.tif 000015
ただし、ａ０、ｂ０、ｃ０は正の実数である。ここで、（式１４）をフーリエ変換すると次の（式１５）が導かれる。
【数１５】
2005114568.tif 000016
【０１０４】  ２×２のサポートをもつ画像の場合、最高周波数成分のとり得る位相の組み合わせ（（ａｍ−ｂｍ−ｃｍ）の位相，（ａｍ＋ｂｍ−ｃｍ）の位相，（ａｍ−ｂｍ＋ｃｍ）の位相）は、（０，０，０）、（π，０，π）、（π，π，０）、（π，０，０）の計４つのうちのいずれかである。ここで、最高周波数成分とは、（式１４）にフーリエ変換を施して得られたフーリエ係数が実数となる周波数、つまり、（式１５）においてフーリエ係数が実数となる周波数のことである。上記のように、ＥＲアルゴリズムによって再構成画像を原画像に収束させるためには、原画像の最高周波数成分の位相と一致している初期画像からＥＲアルゴリズムを開始しなければならない。そこで、上記４つの位相の組み合わせの全体を探索空間として、各個体は、そのうちの１つの位相を表現する遺伝的アルゴリズムを設計する。個体の評価は、その個体が表現する位相と一致する初期画像にＥＲアルゴリズムを施した場合に得られる再構成画像のフーリエ誤差を用いる。このような遺伝的アルゴリズム処理を行えば、フーリエ誤差が最小となる個体がその探索により得られ、再構成画像は、最良、つまり、真実の画像となる。【０１０５】  本手法が対象とする大規模物体の場合には、図３８のように遺伝的アルゴリズムを適用する前に、得られた原画像のフーリエ振幅｜Ｆ（ｕ，ｖ）｜（ｕ＝１，‥，Ｐ、ｖ＝１，‥，Ｑ）にローパスフィルタ処理を施したフーリエ振幅｜Ｆ（ｕ’，ｖ’）｜（ｕ’＝１，‥，Ｐ’、ｖ’＝１，‥，Ｑ’）を算出する。次に、算出されたフーリエ振幅｜Ｆ（ｕ’，ｖ’）｜を用いて、図３８に示す遺伝的アルゴリズムにより画像を再構成する。再構成された画像は、サイズがＰ’＝Ｐ／Ｌ１、Ｑ’＝Ｑ／Ｌの縮小画像となる。すなわち、得られた縮小画像が正しく収束したとすると、この縮小画像は、真実の画像（Ｐ画素×Ｑ画素サイズの画像）に対して、ローパスフィルタ処理を施した後、ダウンサンプリングを行って得られる縮小された画像（Ｐ／Ｌ１画素×Ｑ／Ｌ２画素サイズの画像）と一致する。画像のダウンサンプリングを施した縮小画像は、縮小する前の画像の輪郭線の詳細を知ることはできないが、大体の形状を得るためには有効な情報を与える。【０１０６】  そこで、本実施の形態では、｜Ｆ（ｕ，ｖ）｜（ｕ＝１，‥，Ｐ、ｖ＝１，‥，Ｑ）から｜Ｆ（ｕ’，ｖ’）｜（ｕ’＝１，‥，Ｐ’、ｖ’＝１，‥，Ｑ’）を得ることにより、フーリエ振幅から構成される領域の位相の組み合わせの最適値を高速かつ高精度に算出しておき、その最適値を用いて、真実の画像の輪郭線の大体の形状を得ることができる。そして、得られた大体の形状に基づいて、スネークスの初期値を設定することにより、高精度な画像再構成を実現することができる。この結果、収束速度の高速化を図ることができるとともに、大規模物体（つまり、再構成画像内部の濃淡値が複雑な物体）への適用が可能になる。【０１０７】  ｜Ｆ（ｕ，ｖ）｜（ｕ＝１，‥，Ｐ、ｖ＝１，‥，Ｑ）から｜Ｆ（ｕ’，ｖ’）｜（ｕ’＝１，‥，Ｐ’、ｖ’＝１，‥，Ｑ’）を得る具体的な処理としては、原画像のフーリエ振幅｜Ｆ（ｕ，ｖ）｜（ｕ＝１，‥，Ｐ、ｖ＝１，‥，Ｑ）に対して、ローパスフィルタ処理を施し、フーリエ振幅｜Ｆ（ｕ’，ｖ’）｜（ｕ’＝１，‥，Ｐ’、ｖ’＝１，‥，Ｑ’）を得る。ここで、Ｐ、Ｑは正の整数であり、Ｐ’、Ｑ’は、Ｐ’＜Ｐ、Ｑ’＜Ｑを満たす。通常は、Ｐ’＝Ｐ／Ｌ１、Ｑ’＝Ｑ／Ｌ２が整数となるような整数Ｌ１、Ｌ２を用いる。この場合、ローパスフィルタ処理を施すことにより、フーリエ振幅の低周波領域が得られ、そのサイズは真実のサイズの縮小（ここでは、縦サイズ１／Ｌ１、横サイズ１／Ｌ２）となる。【０１０８】  このようにサイズを縮小することにより、遺伝的アルゴリズムによって探索すべきフーリエ振幅のサイズを小さくすることができ、探索を容易に行うことが可能になる。すなわち、｜Ｆ（ｕ，ｖ）｜（ｕ＝１，‥，Ｐ、ｖ＝１，‥，Ｑ）の最適値を求めるよりは、｜Ｆ（ｕ’，ｖ’）｜（ｕ’＝１，‥，Ｐ’、ｖ’＝１，‥，Ｑ’）の最適値を求めるほうが、探索は容易である。要素数が減少することにより、その組み合わせ総数が減少するためである。【０１０９】  ここで、実際の画像を用いて本手法（遺伝的アルゴリズムを用いた位相回復手法）の有効性を示す。【０１１０】  図３９Ａは、原画像を示す図、図３９Ｂは、初期画像を示す図、図３９Ｃは、図３９Ｂの初期画像から従来のフーリエ反復法を適用して得られた再構成画像を示す図、図３９Ｄは、図３９Ｂの初期画像から本手法を適用して得られた再構成画像を示す図である。従来のフーリエ反復法を用いた場合は、図３９Ｃに示すように、再構成画像は原画像（図３９Ａ）に収束していないが、本手法を用いた場合は、図３９Ｄに示すように、再構成画像は原画像（図３９Ａ）に収束している。【０１１１】  このように、本実施の形態によれば、とり得る位相の組み合わせの中から遺伝的アルゴリズムを用いて初期画像と原画像のある特定の周波数の位相の組み合わせを一致させるため、常に再構成画像を原画像に収束させることができ、特に、正確なフーリエ振幅とサポートがわかっている場合、停滞問題を引き起こさず、高精度に画像を回復させることができる。【０１１２】  また、遺伝的アルゴリズムを適用する前に、フーリエ振幅にローパスフィルタ処理を施し、フーリエ振幅の要素数を減らす、つまり、画像をダウンサンプリングし、画像サイズを縮小するため、短い計算時間で正しい画像を得ることができ、大規模画像においても同様の効果を得ることができる。【０１１３】  なお、本実施の形態では、最適化手法として遺伝的アルゴリズムを用いているが、これに限定されるわけではない。例えば、グローバルサーチ（広域探索）である遺伝的アルゴリズム（ＧＡ）以外に、ローカルサーチ（近傍探索または局所探索）であるＳＡ（Ｓｉｍｕｌａｔｅｄ  Ａｎｎｅａｌｉｎｇ：焼きなまし法）やトレリスサーチ（Ｔｒｅｌｌｉｓ  Ｓｅａｒｃｈ）の適用も可能である。また、ＧＡ＋ＳＡなどの組み合わせも可能である。【０１１４】  また、各実施の形態における手法は併用可能である。すなわち、実施の形態１における手法は、サポートが正確にわかっていない場合の解決策に関する手法であり、実施の形態２における手法は、ＥＲアルゴリズムの問題点（たとえフーリエ振幅とサポートが正確にわかっていても正しい解には収束しない場合がある）を改善した再構成手法であるため、必要に応じて併用することができる。【０１１５】  より具体的には、実施の形態１の手法の再構成手法としては、フーリエ振幅とサポートが正確にわかっている場合は、正しい解に収束するということを前提に、ＥＲアルゴリズムを用いている。しかしながら、ＥＲアルゴリズムは、フーリエ振幅とサポートが正確にわかっていても、正しい解に収束しない場合が生じるため、実施の形態１の手法のＥＲアルゴリズムの部分を実施の形態２の再構成手法に置き換えることが可能である。これにより、より高精度な再構成手法を実現することができる。【０１１６】  例えば、両者を併用する場合は、実施の形態１の手法における画像を再構成する部分を、ＥＲアルゴリズムから実施の形態２の再構成手法に置き換える。具体的には、図２８のステップＳ１９００におけるＥＲアルゴリズムの部分を実施の形態２の再構成手法に置き換える。【０１１７】  また、併用しない場合も考えられる。【０１１８】  例えば、実施の形態２の手法は、フーリエ振幅とサポートが正確にわかっている場合に用いられる。ただし、実際に測定されるサポートは、真のサポートとずれて測定される。これは、実施の形態１の手法を用いることにより解決される。【０１１９】  また、実施の形態１の手法は、測定されるサポートが正確にわかっていない場合に用いられる。一般に、測定されるサポートは真のサポートとずれているため、この手法を用いることにより正しい画像に再構成することができる。ただし、画像を再構成する際、ＥＲアルゴリズムを用いるため、正しい画像に収束しない場合がある。これは、実施の形態２の手法をＥＲアルゴリズムの代わりに用いることにより、解決される。【０１２０】  （実施の形態３）
実施の形態３は、本発明を動画像に応用した場合である。【０１２１】  実施の形態１、２は、静止画像を再構成する場合を想定しており、特に、実施の形態１の画像再構成方法（スネークスを用いた位相回復手法）を単純に動画像に適用することはできない。その理由は、次のとおりである。【０１２２】  図２８に示す実施の形態１の方法は、測定されたフーリエ振幅とサポート（細部情報である原画像の輪郭情報）を用いて画像を再構成する、換言すれば、フーリエ振幅と原画像の輪郭情報から、原画像の内部の細部情報を回復させるものである。したがって、膨張の度合を大きくし過ぎると、最外郭の輪郭形状が変化し、細部情報の誤差が大きくなり、画像の再構成に失敗する場合がある。すなわち、初期設定のサポートが真のサポートを含んでいない場合には、原画像の細部情報が欠落し、真のサポートに収束しないため、画像の再構成に失敗する場合がある。なお、初期設定のサポートが四角の場合にも、原画像の細部情報が大幅に欠落してしまい、また、余分な領域が広いため、画像の再構成に失敗する場合がある。【０１２３】  この点、動画像では、撮像対象が動くため、測定されたサポート（ひいては初期設定のサポート）が真のサポートを含んでいない場合がある。この場合、上記のように、原画像の細部情報が欠落するため、動画像の一部の再構成に失敗し、動画像の回復の精度が低下する場合がある。【０１２４】  本発明者は、静止画像に適用される実施の形態１の手法を動画像に応用するためには、複数の時間フレームの動画像における動き（後述するように、関数で表現）を実施の形態１の手法に導入する必要があることを見いだした。また、そのためには、ある時間フレームで得られた高精度の再構成画像におけるサポートを用いて次の時間フレームにおけるフーリエ振幅およびサポートを更新する必要があることを見いだしたのである。【０１２５】  本発明は、動画像を構成する複数の時間フレームにおいて、各時間フレームのフーリエ振幅およびサポートから動き関数（動きを表す関数）を導出し、導出した動き関数により各時間フレームのフーリエ振幅およびサポートを更新し、更新した各時間フレームのフーリエ振幅およびサポートを実施の形態１の手法に適用して各時間フレームの再構成画像を導出することにより、動画像を再構成するものである。【０１２６】  なお、以下では、適宜、測定により得られるフーリエ振幅およびサポートのことをそれぞれ「測定フーリエ振幅」および「測定サポート」と呼び、更新により得られるフーリエ振幅およびサポートのことをそれぞれ「更新フーリエ振幅」および「更新サポート」と呼び、再構成画像から得られるサポートのことを「計算サポート」と呼ぶことにする。【０１２７】  図４０は、本発明の実施の形態３に係る画像再構成方法のアルゴリズムの手順を示すフローチャートである。また、図４１は、図４０のフーリエ振幅／サポート更新処理の内容を示すフローチャートである。【０１２８】  まず、ステップＳ５０００では、動画像を構成する（Ｋ＋１）個の時間フレームを取得する。ここで、Ｋは所定の整数であり、また、時間フレームの間隔はあらかじめ適切に設定されている。【０１２９】  例えば、動きのある被撮像物体（細胞や結晶など）を撮像して観察する場合、撮像中、被撮像物体は、移動（平行移動、回転）と変形（伸縮、部分形状の変化）を繰り返し、図４２に示すように、平行移動、回転、拡大、縮小の動きを示す。図４２の例では、任意の４つの連続する時間フレームｋ、（ｋ＋１）、（ｋ＋２）、（ｋ＋３）を示している。【０１３０】  そして、ステップＳ５１００では、ステップＳ５０００で取得した各時間フレームにおけるフーリエ振幅｜Ｆ０（ｕ，ｖ）｜〜｜ＦＫ（ｕ，ｖ）｜およびサポートＤ０〜ＤＫをそれぞれ測定する。後述するように、測定サポートを用いて、あるサポートから次のサポートへの動き情報（平行移動、回転、拡大、縮小）を表現する関数を得ることができる。なお、以下、便宜上、ｋ（ｋ＝０，１，‥，Ｋ）番目の時間フレーム（以下「第ｋフレーム」という）におけるフーリエ振幅およびサポートを、それぞれ、｜Ｆｋ｜、Ｄｋで表すことにする。【０１３１】  そして、ステップＳ５２００では、カウンタｋの値を０に設定する。【０１３２】  そして、ステップＳ５３００では、第ｋフレームに対して、実施の形態１における図２８に示す一連の処理（スネークスを用いた位相回復アルゴリズム）を実行する。ただし、ｋ＝１，２，‥，Ｋの場合、ステップＳ１０００で入力する測定フーリエ振幅およびステップＳ１１００で入力する測定サポートを、後述するステップＳ５８００で更新したフーリエ振幅およびサポートにそれぞれ読み替えるものとする。【０１３３】  そして、ステップＳ５４００では、ステップＳ５３００で得た出力画像（つまり、ステップＳ２４００で得た出力画像ｇＮ（ｘ，ｙ））を第ｋフレームの再構成画像ｇｋＮ（ｘ，ｙ）として保存する。【０１３４】  そして、ステップＳ５５００では、ステップＳ５４００で得た再構成画像ｇｋＮ（ｘ，ｙ）から計算サポートＤｋを取得し、保存する。【０１３５】  そして、ステップＳ５６００では、カウンタｋの値を１だけインクリメントする。【０１３６】  そして、ステップＳ５７００では、カウンタｋの値がＫよりも大きいか否かを判断する。この判断の結果として、ｋ≦Ｋの場合は（Ｓ５７００：ＮＯ）、ステップＳ５８００に進み、ｋ＞Ｋの場合は（Ｓ５７００：ＹＥＳ）、最後の第Ｋフレームまで処理が終了したものと判断して、ステップＳ５９００に進む。【０１３７】  ステップＳ５８００では、第ｋフレームに対して（ここで、ｋ＝１，２，‥，Ｋ）、フーリエ振幅およびサポートの更新処理を行う。具体的には、１つ前の第（ｋ−１）フレームの計算サポートＤｋ−１と当該第ｋフレームの測定サポートＤｋとから動き関数を導出し、導出した動き関数により当該第ｋフレームのフーリエ振幅およびサポートを更新した後、ステップ３５３００に戻る。【０１３８】  本ステップの具体的な処理手順は、例えば、図４１に示すとおりである。【０１３９】  まず、ステップＳ５８１０では、ステップＳ５５００で取得し保存した第（ｋ−１）フレームの計算サポートＤｋ−１を入力する。【０１４０】  そして、ステップＳ５８２０では、ステップＳ５１００で測定した第ｋフレームの測定サポートＤｋを入力する。【０１４１】  なお、ステップＳ５８１０とステップＳ５８２０の２つの処理は、逆の順番であってもよい。【０１４２】  そして、ステップＳ５８３０では、サポートＤｋ−１からサポートＤｋへの動きを表す関数として関数ＦＯＰＴを算出する。【０１４３】  ここで、本実施の形態で用いる動き関数Ｆを定義しておく。サポートＤｋ−１からサポートＤｋへの動き関数Ｆは、Ｄｋ＝Ｆ（Ｄｋ−１）によって定義される。この関係は、第（ｋ−１）フレームのサポートをＤｋ−１（ｘ，ｙ）、第ｋフレームのサポートをＤｋ（Ｘ，Ｙ）とすると、平行移動（例えば、ｘ０、ｙ０）、回転（例えば、θ）、拡大・縮小（例えば、ａ、ｅ）のパラメータを用いて、次の（式１６）、つまり、
【数１６】
2005114568.tif 000017
で表現することができる。すなわち、関数Ｆは、平行移動、回転、拡大・縮小のアフィン変換を表している。【０１４４】  例えば、図４３に示すように、連続する２つの時間フレームにおいて、破線の四角形から実線の四角形への動きに関して、図４３Ａは、（−ｘ０，−ｙ０）移動を示し、図４３Ｂは、θ回転を示し、図４３Ｃは、ａ倍、ｅ倍の拡大・縮小を示し、図４３Ｄは、（ｘ０，ｙ０）移動を示している。【０１４５】  なお、動き関数Ｆは、（式１６）では３×３行列で表現されているが、次の（式１７）のように、
【数１７】
2005114568.tif 000018
２×２行列で表現することもできる。【０１４６】  また、これらの動き関数Ｆのパラメータの値域、つまり、被撮像物体の移動（平行移動、回転）および変形（伸縮、部分形状の変化）の度合（変化量の上限）は、事前の知識によりおおよそ得ることができる。例えば、一例として、膨張・収縮の場合は、真のサポート（図４４Ａ）に対して±７画素であり（図４４Ｂは膨張の場合、図４４Ｃは収縮の場合）、回転の場合は、真のサポート（図４５Ａ）に対して４５度である（図４５Ｂ）。【０１４７】  このような動き関数Ｆを用いて、関数Ｆｏｐｔを次のように定義する。すなわち、本実施の形態では、第（ｋ−１）フレームの計算サポートＤｋ−１を動き関数Ｆに適用して得られるサポートＦ（Ｄｋ−１）と、第ｋフレームの測定サポートＤｋとの差の絶対値の二乗、つまり、｜Ｆ（Ｄｋ−１）−Ｄｋ｜２が最小になるときの動き関数Ｆを、関数ＦＯＰＴと定義する。ここで、Ｄｋは、Ｄｋ（ｘ，ｙ）のことであり、第ｋフレームの測定サポートの画素値を表している。また、Ｆ（Ｄｋ−１）は、Ｄｋ−１（ｘ，ｙ）、つまり、第（ｋ−１）フレームの計算サポートに、関数Ｆで表したアフィン変換を施した後の画素値を表している。したがって、｜Ｆ（Ｄｋ−１）−Ｄｋ｜２は、これらの画素値の差の絶対値を二乗したものである。【０１４８】  そして、ステップＳ５８４０では、上記した被撮像物体の移動（平行移動、回転）および変形（伸縮、部分形状の変化）の度合（変化量の上限）の範囲内において、平行移動、回転、拡大・縮小の各パラメータの値を変更して関数ＦＯＰＴを変形させることにより、つまり、関数ＦＯＰＴに対して平行移動、回転、拡大・縮小の各パラメータの値を変更することにより、Ｎ個の関数Ｆ（ｎ）（ｎ＝１，２，‥，Ｎ）を算出する。ここで、Ｎは所定の整数であり、Ｎはあらかじめ適切に設定されている。【０１４９】  そして、ステップＳ５８５０では、ステップＳ５８３０で得た関数ＦＯＰＴとステップＳ５８４０で得たＮ個の関数Ｆ（ｎ）（ｎ＝１，２，‥，Ｎ）を用いて、対応する測定フーリエ振幅｜Ｆｋ｜および測定サポートＤｋをそれぞれ更新し、更新後の（Ｎ＋１）組のフーリエ振幅およびサポートから、（式５）を用いて、対応するフーリエ誤差をそれぞれ算出する。【０１５０】  ここで、測定サポートＤｋの更新は、（式１６）または（式１７）によって行われる。また、測定フーリエ振幅｜Ｆｋ｜の更新は、（式１６）または（式１７）における関数Ｆが示す動き情報のうち、回転のパラメータθのみを用いて、次の（式１８）、つまり、
【数１８】
2005114568.tif 000019
によって行われる。ここで、（ｕ，ｖ）は、更新前のフーリエ振幅｜Ｆｋ（ｕ，ｖ）｜、つまり、ステップＳ５１００で測定したフーリエ振幅の成分であり、（Ｕ，Ｖ）は、更新後のフーリエ振幅｜Ｆｋ（Ｕ，Ｖ）｜の成分である。このように、フーリエ振幅の更新は、サポートの更新に伴って一義的に行われる。【０１５１】  そして、ステップＳ５８６０では、関数Ｆｏｐｔと関数Ｆ（ｎ）（ｎ＝１，２，‥，Ｎ）の（Ｎ＋１）個の関数の中から、ステップＳ５８５０で算出したフーリエ誤差が最小となる関数を選択して関数Ｆ（ｉ）とする。【０１５２】  そして、ステップＳ５８７０では、ステップＳ５８６０で選択した関数Ｆ（ｉ）により第ｋフレームのフーリエ振幅｜Ｆｋ｜およびサポートＤｋを更新する。【０１５３】  なお、上記のように、フーリエ振幅の更新は、動き情報として回転のみを用い、サポートの更新は、動き情報として、回転に加えて、平行移動、拡大、縮小を用いる。したがって、フーリエ振幅は、回転の情報を表す２×２行列によって更新が行われ（（式１８）参照）、サポートは、通常、回転に加えて平行移動、拡大、縮小の情報を表す３×３行列によって更新が行われる（（式１６）参照）。ただし、上記のように、（式１６）で示す３×３行列は、（式１７）のように２×２行列でも表すことができる。【０１５４】  このように、本手法によれば、関数ＦＯＰＴのみを用いてフーリエ振幅｜Ｆｋ｜およびサポートＤｋを更新する場合に生じる問題を解決することができる。【０１５５】  すなわち、具体的には、関数ＦＯＰＴのみを用いて測定フーリエ振幅｜Ｆｋ｜および測定サポートＤｋを更新する場合には、上記のように、｜Ｆ（Ｄｋ−１）−Ｄｋ｜２が最小になる動き関数ＦＯＰＴを算出し、その関数ＦＯＰＴにより各時間フレームの測定フーリエ振幅および測定サポートを更新し、更新後のフーリエ振幅およびサポートにより動画を再構成する。しかし、関数ＦＯＰＴは、元々は測定サポートから算出するため、測定誤差が生じた場合、関数ＦＯＰＴは最適にはならず、高精度な再構成画像を得ることができないという問題がある。【０１５６】  そこで、この問題を解決するために、本手法では、サポートに測定誤差が生じた場合においても、最適な関数Ｆ（ｉ）を算出することにより、高精度な動画の再構成を実現するようにしている。【０１５７】  ここで、実際の画像を用いて具体的に説明する。【０１５８】  図４６Ａは、関数ＦＯＰＴが最適でない場合において当該関数ＦＯＰＴにより更新されたサポート（Ｄｋ＝ＦＯＰＴ（Ｄｋ−１））を示し、図４６Ｂは、図４６Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像のフーリエ誤差の推移を示している。また、図４７Ａは、関数ＦＯＰＴと関数Ｆ（ｎ）（ｎ＝１，２，‥，Ｎ）の（Ｎ＋１）個の関数のうち、関数ＦＯＰＴ以外の選択された最適な関数Ｆ（ｉ）（関数Ｆ１と表記する）により更新されたサポート（Ｄｋ＝Ｆ１（Ｄｋ−１））を示し、図４７Ｂは、図４７Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像のフーリエ誤差の推移を示している。また、図４８Ａは、関数Ｆｏｐｔと関数Ｆ（ｎ）（ｎ＝１，２，‥，Ｎ）の（Ｎ＋１）個の関数のうち、関数Ｆｏｐｔと最適な関数Ｆ（ｉ）以外のある関数（関数Ｆ２と表記する）により更新されたサポート（Ｄｋ＝Ｆ２（Ｄｋ−１））を示し、図４８Ｂは、図４８Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像のフーリエ誤差の推移を示している。また、図４９Ａは、測定結果のサポートを示し、図４９Ｂは、図４９Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像のフーリエ誤差の推移を示している。ここで、関数Ｆ１、Ｆ２は、測定サポート形状から類推されたもの、つまり、上記被撮像物体の移動（平行移動、回転）および変形（伸縮、部分形状の変化）の度合の範囲から類推されたものとする。【０１５９】  図４６〜図４９の例では、例えば、サポートに測定誤差が生じたため、ＦＯＰＴではなく関数Ｆ１が最適な場合であり、関数Ｆ１により更新したサポートを用いて得られた再構成画像のフーリエ誤差が最も小さく、ベストであることがわかる。【０１６０】  このように、本手法によれば、サポートに測定誤差が生じたため、最適なＦＯＰＴを見つけることができなかった場合でも、測定フーリエ振幅および測定サポートを更新するのに最適な関数Ｆ１を算出することができる。【０１６１】  最後に、図４０のステップ３５９００では、ステップＳ５４００で順次保存した各時間フレームの再構成画像ｇ０Ｎ（ｘ，ｙ）〜ｇＫＮ（ｘ，ｙ）から動画像ファイルを生成して、上記一連の処理を終了する。【０１６２】  要するに、本手法では、まず、第０フレームの測定フーリエ振幅｜Ｆ０｜（図５０Ａ）および測定サポートＤ０（図５０Ｂ）に図２８に示す実施の形態１の方法（スネークスを用いた位相回復手法）を適用して第０フレームの再構成画像ｇ０（図５０Ｃ）を求め、得られた再構成画像ｇ０から計算サポートＤ０を取得する。そして、第０フレームの計算サポートＤ０と第１フレームの測定サポートＤ１を用いて、｜Ｆ（Ｄ０）−Ｄ１｜２が最小になる動き関数ＦＯＰＴを求める。そして、平行移動、回転、拡大・縮小の各パラメータの値を変更して関数ＦＯＰＴを変形させることによりＮ個の関数Ｆ（ｎ）（ｎ＝１，２，‥，Ｎ）を求め、得られた関数ＦＯＰＴと関数Ｆ（ｎ）（ｎ＝１，２，‥，Ｎ）の（Ｎ＋１）個の関数の中から、フーリエ誤差が最小になる最適な関数Ｆ（ｉ）を選択し、選択した関数Ｆ（ｉ）により、第１フレームの測定フーリエ振幅｜Ｆ１｜および測定サポートＤ１を更新する。そして、更新後のフーリエ振幅｜Ｆ１｜およびサポートＤ１を用いて第１フレームの再構成画像ｇ１を求め、得られた再構成画像ｇ１から計算サポートＤ１を取得する。そして、以下同様にして、第２フレームから第Ｋフレームまでのすべてのフレームについて、第（ｋ−１）フレームの計算サポートＤｋ−１と第ｋフレームの計算サポートＤｋを用いて（ここでは、ｋ＝２，３，‥，Ｋ）、最適な関数Ｆ（ｉ）を選択し、選択した関数Ｆ（ｉ）により、第ｋフレームの測定フーリエ振幅｜Ｆｋ｜（図５１Ａ）および測定サポートＤｋ（図５１Ｂ）を更新し、更新後のフーリエ振幅｜Ｆｋ｜およびサポートＤｋを用いて第ｋフレームの再構成画像ｇｋ（図５１Ｃ）を求め、得られた再構成画像ｇｋから計算サポートＤｋを取得する。このようにして、順次、測定フー���エ振幅｜Ｆ１｜〜｜ＦＫ｜および測定サポートＤ１〜ＤＫを更新し、これら更新した各時間フレームのフーリエ振幅｜Ｆ１｜〜｜ＦＫ｜およびサポートＤ１〜ＤＫを、図２８に示す実施の形態１の方法に適用し、各時間フレームの再構成画像を導出することにより動画を再構成する。【０１６３】  ここで、実施の形態１の画像再構成手法（位相回復手法）を単純に動画像に適用できない理由について、図面を用いて再度説明する。【０１６４】  上記のように、一般的な撮像物体における、膨張・収縮時のサポートの観測誤差は、真のサポートに対して±７画素であり、回転時のサポートの観測誤差は、真のサポートに対して４５度である。【０１６５】  図５２Ａは、フィルタ窓９×９を用いて図４４Ｃの測定サポート（真のサポートに対して−７画素収縮している場合）に膨張処理を施して得られた初期サポートを示し、図５２Ｂは、図５２Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像を示している。また、図５３Ａは、フィルタ窓１１×１１を用いて図４４Ｃの測定サポートに膨張処理を施して得られた初期サポートを示し、図５３Ｂは、図５３Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像を示している。また、図５４Ａは、フィルタ窓１３×１３を用いて図４４Ｃの測定サポートに膨張処理を施して得られた初期サポートを示し、図５４Ｂは、図５４Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像を示している。また、図５５Ａは、フィルタ窓１７×１７を用いて図４４Ｃの測定サポートに膨張処理を施して得られた初期サポートを示し、図５５Ｂは、図５５Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像を示している。また、図５６Ａは、フィルタ窓２５×２５を用いて図４４Ｃの測定サポートに膨張処理を施して得られた初期サポートを示し、図５６Ｂは、図５６Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像を示している。また、図５７Ａは、フィルタ窓３１×３１を用いて図４４Ｃの測定サポートに膨張処理を施して得られた初期サポートを示し、図５７Ｂは、図５７Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像を示している。また、図５８Ａは、フィルタ窓３３×３３を用いて図４４Ｃの測定サポートに膨張処理を施して得られた初期サポートを示し、図５８Ｂは、図５８Ａのサポートを用いて実施の形態１の方法を適用して得られた再構成画像を示している。【０１６６】  図５２〜図５７から、膨張の度合がフィルタ窓９×９からフィルタ窓３１×３１である場合は、再構成画像は正しい画像に収束することがわかる。これに対し、図５８から、膨張の度合がフィルタ窓３３×３３である場合は、画像の再構成に失敗していることがわかる。この画像再構成の失敗は、膨張による細部情報欠落に起因するものである。【０１６７】  すなわち、実施の形態１の方法は、上記のように、フーリエ振幅とサポート（細部情報である原画像の輪郭情報）を用いて画像を再構成するものであり、より具体的には、フーリエ振幅と原画像の輪郭情報から、原画像の内部の細部情報を回復させるものである。したがって、膨張の度合を大きくし過ぎると、最外郭の��郭形状が変化し、細部情報の誤差が大きくなり、画像の再構成に失敗する。【０１６８】  しかし、実施の形態３の方法によれば、前フレームの高精度の再構成画像から更新したサポートにより画像を再構成するため、膨張による細部情報欠落による画像再構成の失敗を回避することができる。【０１６９】  このように、実施の形態によれば、被撮像物体が移動・変形する場合においても、動画像を高精度に再構成することができる。【０１７０】  本明細書は、２００４年５月２４日出願の特願２００４−１５４００６に基づく。この内容はすべてここに含めておく。
【産業上の利用可能性】
【０１７１】  本発明に係る画像再構成方法は、真のサポートが正確にわかっていない場合においても、位相回復アルゴリズムを用いて高精度に画像を回復することができる画像再構成方法として有用である。また、大規模画像においても、正確なフーリエ振幅とサポートがわかっている場合、停滞問題を引き起こさず、高精度に画像を回復させることができる位相回復アルゴリズムにおける画像の再構成方法として有用である。【０１７２】  さらには、被撮像物体が移動・変形する場合においても、動画像を高精度に再構成することができる画像再構成方法として有用である。
【図１】
2005114568.tif 000020
【図２】
2005114568.tif 000021
【図３】
2005114568.tif 000022
【図４】
2005114568.tif 000023
【図５】
2005114568.tif 000024
【図６】
2005114568.tif 000025
【図７】
2005114568.tif 000026
【図８】
2005114568.tif 000027
【図９】
2005114568.tif 000028
【図１０】
2005114568.tif 000029
【図１１】
2005114568.tif 000030
【図１２】
2005114568.tif 000031
【図１３】
2005114568.tif 000032
【図１４】
2005114568.tif 000033
【図１５】
2005114568.tif 000034
【図１６】
2005114568.tif 000035
【図１７】
2005114568.tif 000036
【図１８】
2005114568.tif 000037
【図１９】
2005114568.tif 000038
【図２０】
2005114568.tif 000039
【図２１】
2005114568.tif 000040
【図２２】
2005114568.tif 000041
【図２３】
2005114568.tif 000042
【図２４】
2005114568.tif 000043
【図２５】
2005114568.tif 000044
【図２６】
2005114568.tif 000045
【図２７】
2005114568.tif 000046
【図２８】
2005114568.tif 000047
【図２９】
2005114568.tif 000048
【図３０】
2005114568.tif 000049
【図３１】
2005114568.tif 000050
【図３２】
2005114568.tif 000051
【図３３】
2005114568.tif 000052
【図３４】
2005114568.tif 000053
【図３５】
2005114568.tif 000054
【図３６】
2005114568.tif 000055
【図３７】
2005114568.tif 000056
【図３８】
2005114568.tif 000057
【図３９】
2005114568.tif 000058
【図４０】
2005114568.tif 000059
【図４１】
2005114568.tif 000060
【図４２】
2005114568.tif 000061
【図４３】
2005114568.tif 000062
【図４４】
2005114568.tif 000063
【図４５】
2005114568.tif 000064
【図４６】
2005114568.tif 000065
【図４７】
2005114568.tif 000066
【図４８】
2005114568.tif 000067
【図４９】
2005114568.tif 000068
【図５０】
2005114568.tif 000069
【図５１】
2005114568.tif 000070
【図５２】
2005114568.tif 000071
【図５３】
2005114568.tif 000072
【図５４】
2005114568.tif 000073
【図５５】
2005114568.tif 000074
【図５６】
2005114568.tif 000075
【図５７】
2005114568.tif 000076
【図５８】
2005114568.tif 000077
【国際調査報告】
2005114568.tif 000078
2005114568.tif 000079
2005114568.tif 000080
2005114568.tif 000081
2005114568.tif 000082
2005114568.tif 000083
2005114568.tif 000084