JP5739976B2 - 位相画像を再構成する方法、コンピュータプログラム、位相画像を再構成する演算装置 - Google Patents

位相画像を再構成する方法、コンピュータプログラム、位相画像を再構成する演算装置 Download PDF

Info

Publication number
JP5739976B2
JP5739976B2 JP2013266132A JP2013266132A JP5739976B2 JP 5739976 B2 JP5739976 B2 JP 5739976B2 JP 2013266132 A JP2013266132 A JP 2013266132A JP 2013266132 A JP2013266132 A JP 2013266132A JP 5739976 B2 JP5739976 B2 JP 5739976B2
Authority
JP
Japan
Prior art keywords
phase
image
interference fringe
dimensional
parameter
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
JP2013266132A
Other languages
English (en)
Other versions
JP2014121614A (ja
Inventor
ジェイムズ ボーン ドナルド
ジェイムズ ボーン ドナルド
ウォン サンプソン
ウォン サンプソン
キエラン　ジェラード　ラルキン
ジェラード ラルキン キエラン
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Canon Inc
Original Assignee
Canon Inc
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Canon Inc filed Critical Canon Inc
Publication of JP2014121614A publication Critical patent/JP2014121614A/ja
Application granted granted Critical
Publication of JP5739976B2 publication Critical patent/JP5739976B2/ja
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/02—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
    • G01N23/04—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material
    • G01N23/041—Phase-contrast imaging, e.g. using grating interferometers
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/20—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by using diffraction of the radiation by the materials, e.g. for investigating crystal structure; by using scattering of the radiation by the materials, e.g. for investigating non-crystalline materials; by using reflection of the radiation by the materials
    • G01N23/20075—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by using diffraction of the radiation by the materials, e.g. for investigating crystal structure; by using scattering of the radiation by the materials, e.g. for investigating non-crystalline materials; by using reflection of the radiation by the materials by measuring interferences of X-rays, e.g. Borrmann effect
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography

Description

本発明は、２次元の干渉縞インターフェログラム画像のセットからの位相画像を再構成する方法、該方法を実行可能なプログラムが記録されたコンピュータ可読記憶媒体、及び該方法を実行可能なコンピュータシステムを備える装置に関する。
過去約１００年において、撮影装置において（通常は不可視である）位相変化を可視にするために多くの技術が開発された。これらの技術は、ゼルニケ位相コントラスト法、ノマルスキー微分干渉コントラスト法、一般化位相コントラスト法、フーコーナイフエッジ法、シュリーレン法、シャドウグラフ法、暗視野法及びワイヤテスト法を含む。最近、上記の技術のうちのいくつかのＸ線撮影への適用が進んだ。Ｘ線ビームの集束及び撮影が困難であることが主な原因であるＸ線の特に扱いにくい性質に対して、多くの新規の位相コントラスト技術が開発された。これらの技術は、ＴＩＥ（強度輸送方程式）、位相コントラストイメージング及びタイコグラフィを含む。別のそのような技術は「Ｘ線タルボモアレ干渉法」（ＸＴＭＩ）として既知であり、これは、１つ以上の微分位相画像を符号化する中間画像を生成する。単純な線形格子を使用して実現される場合、ＸＴＭＩ法は１つの微分位相画像を生成する。２次元（交差）格子を用いて実現されるＸＴＭＩは、２つの（交差又は直交）微分位相画像を生成する。
位相シフト法は周知であり、干渉縞の解析に広く使用されてきた。位相シフト干渉法（ＰＳＩ）は、１９７４年に最初に導入された。ＰＳＩでは、参照鏡を少しずつ線形に増分させて移動し且つ各ステップにおけるインターフェログラム画像を格納する必要がある。インターフェログラム画像における各点の値及び人工的に導入される位相シフトの量は、当該点の対応する光路長変化を計算するために使用される。複数の正弦波成分を用いる位相シフト技術が更に提案されている。
位相シフト法は、Ｘ線タルボモアレ干渉法に更に適用されてきた。１つの傾斜方向に対してのみ感知する線形格子を用いる場合、高調波を無視できるならば、画像内の点における低周波数背景強度パラメータ、干渉縞変調振幅パラメータ及び位相パラメータを回復するためには少なくとも３つの画像が必要であることは既知である。しかし、このシステムは、格子に対して直交する微分挙動のみを感知する。格子を回転して露光を更に２回行うことにより、全部で５回のみの露光で直交方向における微分位相及び変調パラメータを取得できる。これは機械的に複雑であり、測定中に格子を９０度回転する必要がある。
交差又は２次元（２Ｄ）格子（以下まとめて交差格子と称する）により、システムは２つの直交方向における微分位相及び変調を同時に感知できる。位相シフトを用いる交差格子を有するＸ線タルボシステムに対する撮影処理の数学モデルは、シーケンスにおいてｎ番目の画像の画素における撮像された強度Ｚｎが以下の形式で表されると仮定する。
式中、ａは低周波数背景強度であり、第２項であるｍ１，０ｃｏｓ（ξ１，０＋Φ１，０，ｎ）は、変調強度ｍ１，０、物体により生じた位相ξ１，０及び与えられた位相ステップφ１，０，ｎを有する被変調正弦波である。この被変調正弦波は、交差格子を形成する格子対の一方の要素から発生する。第３項は、同一の一般的な形式を有するが、交差格子の第２の要素から発生する。第４項及び第５項は同様の形式であるように見えるが、第２項及び第３項の相互作用から発生する。そのため、位相ステップφ１，−１，ｎ及びφ−１，１，ｎは個別に与えられるのではなく、以下のように最初の２つの項から生じる。
この強度の公式は９つの未知のパラメータを有するため、９つの未知のパラメータを求めるためには、９つの式に関連する９つの強度値（Ｚ１〜Ｚ９）を生成する少なくとも９つの位相ステップ画像が必要とされる。位相ステップの数を現在の限界である９ステップより減少できる方法が必要とされる。
最近、他の研究者が、窓フーリエ変換（ＷＦＴ）に基づくＸ線タルボモアレインターフェログラム解析に対して少ないパラメータを用いるモデルを提案した。これは、位相シフト法及び解析に適応可能である。位相シフト法の場合、この単純化されたモデルは以下の形式をとる。
しかし、このモデルが有するパラメータは少ないが、２つの交差する余弦項の変調振幅ｂが同一にされる。これは、実際は誤りである。従って、実際のＸ線タルボモアレインターフェログラムの挙動により適切に当てはまる単純化されたモデルが必要とされる。
交差格子を用いるＸ線タルボモアレ干渉法に対する位相シフト法の別の問題は、位相ステップが観察に適合するために使用されるモデルにより与えられる正確な値を有する必要があることである。解析法が標準位相ステップを明確に仮定するため、解析が必要とする標準位相ステップから位相ステップがずれた場合、解析により回復される結果として得られるパラメータに誤差が生じる。しかし、位相ステップはシステムにおける非常に小さい機械移動により生成されるため、通常、結果として得られる位相ステップに著しい誤差が生じる。１次元（１Ｄ）干渉縞系を用いる干渉位相シフトに対してこれらの位相ステップ誤差を補正する既知の方法は存在するが、交差格子を用いるＸ線タルボモアレ干渉法に対して位相ステップ誤差を補正する既知の方法は存在しない。これらの位相誤差を補正できる上記システムに対する解析方法が必要とされる。
物体の干渉縞インターフェログラム画像のセットから位相画像を再構成する方法が開示される。該方法は、８つ以下の撮像された干渉縞インターフェログラム画像を含む撮像された干渉縞インターフェログラム画像のセットを使用して、吸収パラメータ（ａ）、２次元振幅変調パラメータ（ｍｘ，ｍｙ）及び２次元位相変調パラメータ（ξｘ，ξｙ）の推定値を適切な閉形式解から判定し、パラメータ推定値を使用して位相画像を再構成する。
判定ステップは、パラメータ推定値を修正するのが望ましい。各画素に対するパラメータ推定値は、当該画素に対するシミュレーションされた干渉縞強度のセットを生成し且つシミュレーションされた干渉縞強度のセットと当該画素における撮像された干渉縞強度のセットとを比較することにより修正されるのが好ましい。あるいは、各画素に対するパラメータ推定値は、当該画素に対するシミュレーションされた干渉縞強度のセットを生成し且つシミュレーションされた干渉縞強度のセット及び撮像された干渉縞強度のセットにおける対応する画素の誤差を最小にすることにより修正される。
特定の実現例において、判定ステップは、パラメータの初期推定値を取得することと、パラメータにおける位相ステップ誤差を補正するために位相推定法を使用することと、補正された位相ステップを使用してパラメータに対する最適解を反復して判定することとを含む。
特定の実現例において、セットは以下を含む。
−５つの撮像された干渉縞インターフェログラム画像を含み、交差回折格子の位相ステップは、φｘ：［０，２π／５，４π／５，６π／５，８π／５］及びφｙ：［０，６π／５，２π／５，８π／５，４π／５］である
−６つの撮像された干渉縞インターフェログラム画像を含み、交差回折格子の位相ステップは、φｘ：［０，π／２，π／２，π，π，３π／２］及びφｙ：［π／２，０，π，π／２，３π／２，π］である
−７つの撮像された干渉縞インターフェログラム画像を含み、交差回折格子の位相ステップは、φｘ：［０，π／２，π，３π／２，０，π／２，π］及びφｙ：［０，３π／２，π，π／２，π，π／２，０］である
−８つの撮像された干渉縞インターフェログラム画像を含み、交差回折格子の位相ステップは、φｘ：［０，π／２，π，３π／２，０，π／２，π，３π／２］及びφｙ：［０，３π／２，π，π／２，π，π／２，０，３π／２］である
別の態様において、交差回折格子を有するＸ線干渉計により撮像された物体の干渉縞インターフェログラム画像のセットから位相画像を再構成する方法が開示される。本方法は、撮像された画像のセットを提供することと、干渉縞に関連する交差成分と干渉縞に関連する主成分とを関連付けることにより８つ以下の画像から画像を再構成できる再構成技術を提供することと、提供された再構成技術を使用して、提供された撮像された画像のセットから位相画像を再構成することとを備える。
他の態様が開示される。
図１は、１つの実現例に係る干渉縞インターフェログラム画像に対する位相復調方法を概略的に示すフローチャートである。 図２は、更なる実現例に係る干渉縞インターフェログラム画像に対する位相復調方法を概略的に示すフローチャートである。 図３は、図２からの反復位相復調方法の詳細を概略的に示すフローチャートである。 図４は、説明する構成に関係するＸ線タルボ干渉計の実験設定を概略的に示す図である。 図５Ａは、図４において説明する格子及び自己像の例を示す図である。 図５Ｂは、図４において説明する格子及び自己像の例を示す図である。 図５Ｃは、図４において説明する格子及び自己像の例を示す図である。 図６は、図１及び図２において説明する画像撮像処理の詳細の例を概略的に示すフローチャートである。 図７は、９つのステップから成る画像撮像処理の位相シフトパターンを示す図である。 図８は、８つのステップから成る画像撮像処理の位相シフトパターンを示す図である。 図９は、別の８つのステップから成る画像撮像処理の位相シフトパターンを示す図である。 図１０は、７つのステップから成る画像撮像処理の位相シフトパターンを示す図である。 図１１は、６つのステップから成る画像撮像処理の位相シフトパターンを示す図である。 図１２は、５つのステップから成る画像撮像処理の位相シフトパターンを示す図である。 図１３Ａは、説明する構成が実現される汎用コンピュータシステムを概略的に示すブロック図である。 図１３Ｂは、説明する構成が実現される汎用コンピュータシステムを概略的に示すブロック図である。
以下の図面を参照して、本発明の少なくとも１つの実施形態を以下に説明する。
状況
本開示は、位相シフト技術を使用するＸ線タルボ干渉法システムから取得されたインターフェログラム画像のセットからの位相情報の抽出に関する。
Ｘ線タルボ干渉法システムは、物体に関する情報を回復するために、タルボ効果を使用する多くの格子を介して、物体により生成された干渉縞を撮像する。この種類のＸ線撮影は、吸収ではなく位相差を使用してコントラストを生成するため、通常のＸ線撮影より精度が高い。
Ｘ線タルボ干渉法（ＸＴＩ）システム４００を図４に示す。図中、位相格子Ｇ１（４１０）及び撮影の対象である物体４０２の双方がＸ線により照明される。タルボ効果により、自己像４２０が位相格子Ｇ１の後方の位置に形成される。Ｘ線波の位相が物体４０２により乱されるため、自己像４２０は変形される。変形された自己像を解析することにより、物体４０２の特徴が演繹される。図４において、吸収格子Ｇ２（４３０）は、モアレ干渉縞を生成するために自己像４２０の位置に配置される。画像センサ４４０は、吸収格子４３０に近接し、モアレ干渉縞を撮像するように構成される。吸収格子Ｇ２（４３０）の周期は自己像４２０の周期と同様であるように構成されるため、自己像４２０と吸収格子Ｇ２（４３０）とを重ね合わせることにより生成されるモアレパターンは自己像の強度変化を増幅し、そのため、画像センサ４４０はパターンを更に容易に分解できる。
Ｘ線干渉法システム４００において使用される格子Ｇ１（４１０）及びＧ２（４３０）は、図５Ａ−５Ｃに示す２Ｄ構造を有する。図５Ａにおいて、５１０は位相格子Ｇ１（４１０）の共通設計を示し、暗い部分は位相シフトπが与えられるＧ１（４１０）上の領域を表し、明るい部分は位相シフトが与えられないＧ１（４１０）上の領域を表す。２Ｄ構造を有する位相格子Ｇ１を用いることにより、２次元の（交差する２方向に周期を有する）自己像が形成される。図５Ｂにおいて、５２０は吸収格子Ｇ２（４３０）の共通設計を示し、暗い部分はＸ線エネルギの大部分を吸収するＧ２（４３０）上の領域を表し、明るい部分はＸ線エネルギを通過させるＧ２（４３０）上の領域を表す。示されるＧ２（４３０）の例５２０は、モアレパターンを生成するために若干傾斜される。このように、２Ｄ構造を有する位相格子Ｇ１及び吸収格子Ｇ２を用いることにより、２次元のモアレパターンが形成される。Ｘ線タルボ干渉法システムは、図４に示されている構成とは異なるように構成されてもよいが、通常は図４のシステムの変形を表す。例えば、Ｘ線源と位相格子Ｇ１との間に線源格子と呼ばれる格子を配置してＸ線源からのＸ線の空間干渉性を向上させても良い。図５Ｃにおいて、５３０は、物体４０２が存在しない場合のシステム４００の自己像４２０の例を示す。吸収格子Ｇ２（４３０）が自己像４２０と同様の構造を有し且つ小さい角度回転された状態で自己像４２０と吸収格子Ｇ２とを重ね合わせることにより、モアレパターンが生成される。図４の物体４０２の位相情報はモアレパターンに残存するが、自己像４２０と比較してはるかに低い周波数で存在する。画像センサ４４０が限られた分解能しか有さない場合、分解能は、誇張された自己像であるモアレパターンを分解し且つ物体の位相情報を回復できれば十分である。尚、本発明及び本明細書では、このモアレパターンのような干渉縞を画像センサで取得した画像を干渉縞インターフェログラム画像と呼ぶ。特に、モアレ又は干渉縞が２次元であることを明確にする場合、２次元の、モアレ又は干渉縞を画像センサで取得した画像を２次元の干渉縞インターフェログラム画像と呼ぶことがある。
図４で説明するＸＴＩシステムにおいて、位相シフト技術は、物体により生じた吸収、振幅変調及び位相変調を回復するために適用されてもよい。画像のセットは、光路に沿うサンプルがモアレパターンの空間変位を介して取得されるように、各撮像の間に格子を移動させて撮像される。本明細書において、位相シフト技術を以下に詳細に説明する。
図１３Ａ及び図１３Ｂは、説明する種々の構成が実現される汎用コンピュータシステム１３００を示す。
図１３Ａに示すように、コンピュータシステム１３００は、コンピュータモジュール１３０１と、キーボード１３０２、マウスポインタ装置１３０３、スキャナ１３２６、カメラ１３２７及びマイク１３８０等の入力装置と、プリンタ１３１５、表示装置１３１４及びスピーカ１３１７を含む出力装置とを含む。コンピュータモジュール１３０１は、接続１３２３を介して図４の構成に従って構成されたＸ線装置１３９９に結合するように更に示される。示されるように、Ｘ線装置１３９９は、図４の物体４０２を表す対象１３２４の胴部の撮影を行うように構成される。外部変復調装置（モデム）送受信装置１３１６は、接続１３２１を介して通信ネットワーク１３２０との間で通信するためにコンピュータモジュール１３０１により使用されてもよい。通信ネットワーク１３２０は、インターネット、移動通信ネットワーク、プライベートＷＡＮ等のワイドエリアネットワーク（ＷＡＮ）又はローカルエリアネットワーク（ＬＡＮ）であってもよい。接続１３２１が電話線である場合、モデム１３１６は従来の「ダイヤルアップ式」モデムであってもよい。あるいは、接続１３２１が大容量（例えば、ケーブル）接続である場合、モデム１３１６は広帯域モデムであってもよい。無線モデムが通信ネットワーク１３２０への無線接続のために更に使用されてもよい。
通常、コンピュータモジュール１３０１は、少なくとも１つの処理ユニット１３０５とメモリユニット１３０６とを含む。例えばメモリユニット１３０６は、半導体ランダムアクセスメモリ（ＲＡＭ）及び半導体読み出し専用メモリ（ＲＯＭ）を有してもよい。コンピュータモジュール１３０１は、ビデオディスプレイ１３１４、スピーカ１３１７及びマイク１３８０に結合するオーディオ／ビデオインタフェース１３０７と、キーボード１３０２、マウス１３０３、スキャナ１３２６、カメラ１３２７及びオプションでジョイスティック又は他のヒューマンインタフェース装置（不図示）に結合するＩ／Ｏインタフェース１３１３と、外部モデム１３１６及びプリンタ１３１５に対するインタフェース１３０８とを含む多くの入出力（Ｉ／Ｏ）インタフェースを更に含む。いくつかの実施例において、モデム１３１６は、例えばインタフェース１３０８内等のコンピュータモジュール１３０１内に組み込まれてもよい。コンピュータモジュール１３０１は、接続１３２３を介してコンピュータシステム１３００をＸ線装置１３９９に結合させるローカルインタフェース１３１１を更に有する。インタフェース１３０８及び１３１１は、ＥｔｈｅｒｎｅｔＴＭ、ＢｌｕｅｔｏｏｔｈＴＭ無線構成又はＩＥＥＥ８０２．１１無線構成等の各規格で動作するように構成されてもよいが、多くの他の種類のインタフェースが実現されてもよい。通常、接続１３２３は、コンピュータ１３０１の指示に基づく装置１３９９によるＸ線画像の自動撮像を可能にし且つメモリ１３０６又はＨＤＤ１３１０に格納するために撮像した画像データをコンピュータ１３０１に通信するために双方向接続である。
Ｉ／Ｏインタフェース１３０８及び１３１３は、直列接続性及び並列接続性のいずれか又は双方を提供してもよい。直列接続性は、通常、ユニバーサルシリアルバス（ＵＳＢ）規格に従って実現され、対応するＵＳＢコネクタ（不図示）を有する。記憶装置１３０９が提供され、これは通常はハードディスクドライブ（ＨＤＤ）１３１０を含む。フロッピディスクドライブ及び磁気テープドライブ（不図示）等の他の記憶装置が更に使用されてもよい。光ディスクドライブ１３１２は、通常、不揮発性データソースとして動作するために提供される。例えば光ディスク（例えば、ＣＤ−ＲＯＭ、ＤＶＤ、Ｂｌｕ−ｒａｙ ＤｉｓｃＴＭ）、ＵＳＢ−ＲＡＭ、ポータブル外部ハードドライブ及びフロッピディスク等のポータブルメモリ装置が、システム１３００に対する適切なデータソースとして使用されてもよい。
通常、コンピュータモジュール１３０１の構成要素１３０５〜１３１３は、結果として当業者には既知であるコンピュータシステム１３００の従来の動作モードになるように相互接続バス１３０４を介して通信する。例えばプロセッサ１３０５は、接続１３１８を使用してシステムバス１３０４に結合される。同様に、メモリ１３０６及び光ディスクドライブ１３１２は接続１３１９によりシステムバス１３０４に結合される。説明する構成が実現されるコンピュータの例は、ＩＢＭ−ＰＣ及び互換機、Ｓｕｎ Ｓｐａｒｃｓｔａｔｉｏｎｓ、Ａｐｐｌｅ ＭａｃＴＭ又は同様のコンピュータシステムを含む。
位相復調方法は、説明する図１〜図１２の処理がコンピュータシステム１３００内で実行可能な１つ以上のソフトウェアアプリケーションプログラム１３３３として実現されてもよいコンピュータシステム１３００を使用して実現されてもよい。特に、位相復調方法のステップは、コンピュータシステム１３００内で実行されるソフトウェア１３３３内の命令１３３１（図１３Ｂを参照）により実現される。ソフトウェア命令１３３１は、各々が１つ以上の特定のタスクを実行する１つ以上のコードモジュールとして形成されてもよい。ソフトウェアは２つの別個の部分に分割されてもよく、第１の部分及び対応するコードモジュールが位相復調方法を実行し、第２の部分及び対応するコードモジュールが第１の部分とユーザとの間のユーザインタフェースを管理する。
ソフトウェアは、例えば後述する記憶装置を含むコンピュータ可読媒体に格納されてもよい。ソフトウェアは、コンピュータ可読媒体からコンピュータシステム１３００にロードされた後にコンピュータシステム１３００により実行される。そのようなソフトウェア又はコンピュータプログラムが記録されたコンピュータ可読媒体はコンピュータプログラム製品である。コンピュータシステム１３００におけるコンピュータプログラム製品の使用は、Ｘ線撮影及び画像処理、特にインターフェログラムの干渉縞の位相復調のための有利な装置を実現するのに好ましい。
通常、ソフトウェア１３３３はＨＤＤ１３１０又はメモリ１３０６に格納される。ソフトウェアは、コンピュータ可読媒体からコンピュータシステム１３００にロードされ、コンピュータシステム１３００により実行される。従って、例えばソフトウェア１３３３は、光ディスクドライブ１３１２により読み出される光学可読ディスク記憶媒体（例えば、ＣＤ−ＲＯＭ）１３２５に格納されてもよい。そのようなソフトウェア又はコンピュータプログラムが記録されたコンピュータ可読媒体はコンピュータプログラム製品である。コンピュータシステム１３００におけるコンピュータプログラム製品の使用は、Ｘ線撮影及び画像処理、特にインターフェログラムの干渉縞の位相復調のための装置を実現するのに好ましい。
いくつかの例において、アプリケーションプログラム１３３３は、１つ以上のＣＤ−ＲＯＭ１３２５上で符号化されてユーザに供給され且つ対応するドライブ１３１２を介して読み出されてもよく、あるいはネットワーク１３２０からユーザにより読み出されてもよい。更に、ソフトウェアは、他のコンピュータ可読媒体からコンピュータシステム１３００にロード可能である。コンピュータ可読記憶媒体は、実行及び／又は処理のためにコンピュータシステム１３００に記録された命令及び／又はデータを提供する何らかの非一時的な有形記憶媒体を示す。そのような記憶媒体の例は、コンピュータモジュール１３０１の内部であるか又は外部であるかに関係なく、フロッピディスク、磁気テープ、ＣＤ−ＲＯＭ、ＤＶＤ、Ｂｌｕ−ｒａｙ ＤｉｓｃＴＭ、ハードディスクドライブ、ＲＯＭ又は集積回路、ＵＳＢメモリ、光磁気ディスク、あるいはＰＣＭＣＩＡカード等のコンピュータ可読カードを含む。コンピュータモジュール１３０１へのソフトウェア、アプリケーションプログラム、命令及び／又はデータの提供に同様に関係してもよい一時的又は無形コンピュータ可読伝送媒体の例は、無線又は赤外線送信チャネル、別のコンピュータ又はネットワーク装置へのネットワーク接続、並びに電子メール送信及びウェブサイト等に記録された情報を含むインターネット又はイントラネットを含む。
上述のアプリケーションプログラム１３３３の第２の部分及び対応するコードモジュールは、ディスプレイ１３１４上にレンダリングされるか又は表される１つ以上のグラフィカルユーザインタフェース（ＧＵＩ）を実現するために実行されてもよい。通常はキーボード１３０２及びマウス１３０３の操作を介して、コンピュータシステム１３００及びアプリケーションのユーザは、ＧＵＩに関連するアプリケーションに制御コマンド及び／又は入力を提供するために機能的に適合可能な方法でインタフェースを操作してもよい。スピーカ１３１７を介して出力される音声プロンプト及びマイク１３８０を介して入力されるユーザの音声コマンドを利用するオーディオインタフェース等の機能的に適合可能なユーザインタフェースの他の形態が更に実現されてもよい。
図１３Ｂは、プロセッサ１３０５及び「メモリ」１３３４の詳細な概略ブロック図である。メモリ１３３４は、図１３Ａのコンピュータモジュール１３０１によりアクセス可能な全てのメモリモジュール（ＨＤＤ１３０９及び半導体メモリ１３０６を含む）の論理集合を表す。
最初にコンピュータモジュール１３０１の電源が投入されると、ＰＯＳＴ（ｐｏｗｅｒ−ｏｎ ｓｅｌｆ−ｔｅｓｔ）プログラム１３５０が実行する。通常、ＰＯＳＴプログラム１３５０は図１３Ａの半導体メモリ１３０６のＲＯＭ１３４９に格納される。ソフトウェアを格納したＲＯＭ１３４９等のハードウェア装置をファームウェアと呼ぶ場合がある。ＰＯＳＴプログラム１３５０は、適切な機能を保証するためにコンピュータモジュール１３０１内のハードウェアを調べ、通常は、プロセッサ１３０５とメモリ１３３４（１３０９、１３０６）と通常はＲＯＭ１３４９に同様に格納される正常な動作のための基本入出力システムソフトウェア（ＢＩＯＳ）モジュール１３５１とをチェックする。ＰＯＳＴプログラム１３５０が正常に実行すると、ＢＩＯＳ１３５１は図１３Ａのハードディスクドライブ１３１０を起動する。ハードディスクドライブ１３１０が起動することにより、ハードディスクドライブ１３１０に常駐するブートストラップローダプログラム１３５２がプロセッサ１３０５を介して実行する。これにより、オペレーティングシステム１３５３がＲＡＭメモリ１３０６にロードされ、動作を開始する。オペレーティングシステム１３５３は、プロセッサ管理、メモリ管理、装置管理、記憶装置管理、ソフトウェアアプリケーションインタフェース及び汎用ユーザインタフェースを含む種々の高レベル機能を実現するためのプロセッサ１３０５により実行可能なシステムレベルアプリケーションである。
オペレーティングシステム１３５３は、コンピュータモジュール１３０１上で実行する各処理又はアプリケーションが別の処理に割り当てられたメモリと衝突せずに実行するのに十分なメモリを有することを保証するために、メモリ１３３４（１３０９、１３０６）を管理する。更に、図１３Ａのシステム１３００において使用可能な異なる種類のメモリは、各処理が効率よく実行できるように適切に使用される必要がある。従って、集約メモリ１３３４は、メモリの特定のセグメントが割り当てられる方法を示すのではなく（特に明記しない限り）、コンピュータシステム１３００によりアクセス可能なメモリの全体図及びそれが使用される方法を提供することを意図する。
図１３Ｂに示すように、プロセッサ１３０５は、制御部１３３９、論理演算装置（ＡＬＵ）１３４０、並びにキャッシュメモリと呼ばれる場合があるローカル又は内部メモリ１３４８を含む多くの機能モジュールを含む。通常、キャッシュメモリ１３４８は、レジスタ部内に多くの記憶レジスタ１３４４〜１３４６を含む。１つ以上の内部バス１３４１は、これらの機能モジュールを機能的に相互接続する。通常、プロセッサ１３０５は、接続１３１８を使用してシステムバス１３０４を介して外部装置と通信するための１つ以上のインタフェース１３４２を更に有する。メモリ１３３４は、接続１３１９を使用してバス１３０４に結合される。
アプリケーションプログラム１３３３は、条件付き分岐命令及びループ命令を含んでもよい一連の命令１３３１を含む。プログラム１３３３は、プログラム１３３３の実行時に使用されるデータ１３３２を更に含んでもよい。命令１３３１は記憶場所１３２８、１３２９、１３３０に格納され、データ１３３２は記憶場所１３３５、１３３６、１３３７に格納される。命令１３３１及び記憶場所１３２８〜１３３０の相対サイズに依存して、特定の命令は、記憶場所１３３０内に示す命令により示されるように、単一の記憶場所に格納されてもよい。あるいは、命令は、記憶場所１３２８及び１３２９内に示す命令セグメントにより示されるように、各々が別個の記憶場所に格納される多くの部分に分割されてもよい。
一般に、プロセッサ１３０５は、そこで実行される命令セットを与えられる。プロセッサ１３０５は後続の入力を待ち、別の命令セットを実行することによりその後続の入力に対応する。各入力は、入力装置１３０２、１３０３の１つ以上により生成されたデータ、ネットワーク１３２０、１３０２のうち１つを介して外部ソースから受信されたデータ、記憶装置１３０６、１３０９の一方から検索されたデータ、あるいは対応する読取装置１３１２に挿入された記憶媒体１３２５から検索されたデータを含む多くのソースのうちの１つ以上から提供されてもよい。これらは全て図１３Ａに示される。いくつかの例において、命令セットを実行した結果、データが出力されてもよい。実行は、メモリ１３３４にデータ又は変数を格納することを更に含んでもよい。
開示する画像処理構成は、対応する記憶場所１３５５、１３５６、１３５７においてメモリ１３３４に格納される入力変数１３５４を使用する。画像処理構成は、対応する記憶場所１３６２、１３６３、１３６４においてメモリ１３３４に格納される出力変数１３６１を生成する。中間変数１３５８は、記憶場所１３５９、１３６０、１３６６及び１３６７に格納されてもよい。
図１３Ｂのプロセッサ１３０５を参照すると、レジスタ１３４４、１３４５、１３４６、論理演算装置（ＡＬＵ）１３４０及び制御部１３３９は協働して、プログラム１３３３を形成する命令セット内の各命令に対する「取り出し、復号化及び実行」サイクルを実行するために必要な一連のマイクロ命令を実行する。各「取り出し、復号化及び実行」サイクルは以下を含む。
（ｉ） 記憶場所１３２８、１３２９、１３３０から命令１３３１を取り出すか又は読み出す取り出し動作
（ｉｉ） 制御部１３３９が取り出された命令を判定する復号化動作
（ｉｉｉ） 制御部１３３９及び／又はＡＬＵ１３４０が命令を実行する実行動作
その後、次の命令に対する更なる「取り出し、復号化及び実行」サイクルが実行されてもよい。同様に、制御部１３３９が記憶場所１３３２に値を格納するか又は書き込む格納サイクルが実行されてもよい。
図１〜図１２の処理における各ステップ又はサブステップは、プログラム１３３３の１つ以上のセグメントに関連し、前述したプログラム１３３３のセグメントに対する命令セット内の各命令に対する「取り出し、復号化及び実行」サイクルを実行するために協働するプロセッサ１３０５内のレジスタ部１３４４、１３４５、１３４７、ＡＬＵ１３４０及び制御部１３３９により実行される。
あるいは、画像処理方法は、画像処理の機能又はサブ機能を実行する１つ以上の集積回路等の専用ハードウェアを使用して全体的又は部分的に実現されてもよい。そのような専用ハードウェアは、グラフィックプロセッサ、デジタル信号プロセッサ又は１つ以上のマイクロプロセッサ及び関連するメモリを含んでもよい。
概要
交差格子を用いるＸ線タルボモアレ干渉法において使用される格子は、シリコンウエハに穴をエッチングすることにより作成される。製造工程における精度の制約は、これらの穴が図５に示すような理想的な正方形のプロファイルより丸くなる傾向があることを意味する。正方形のプロファイルを実際に達成することはできないが、解析に対する好都合な形状が生成される。
上述の問題を解決する時に、干渉の一次項と交差項との関係に強い制約を与え、制約された関係を使用することにより位相シフト処理に必要なフレーム数を減少することが望ましい。
数学モデル
図４のＸ線タルボシステムの位相シフトシーケンスに対するモアレインターフェログラム内の所定の画素における強度を表すために必要な自由パラメータの数を減少する第１のステップとして、以下の分離可能な形式が提案される。
式中、第１の方向は、Ｘ線装置１３９９を使用して形成されたシステム４００のｘ軸とラベル付けされ、第２の方向はｙ軸とラベル付けされ、ａは総吸収係数であり、ｍｘはｘ成分の変調であり、ｍｙはｙ成分の変調であり、ξｘはｘ成分の位相であり、ξｙはｙ成分の位相であり、φｘ，ｎはｘ成分における与えられた位相ステップであり、φｙ，ｎはｙ成分における与えられた位相ステップである。
位相シフト技術は、搬送信号の一連の位相シフトを用いて時間サンプルを提供し、結果として得られる干渉縞強度Ｚｎ、ｎ＝０…Ｍ−１を撮像する。位相ステップφｘ，ｎ及びφｙ，ｎ、ｎ＝０…Ｍ−１を適切に選択することにより、ａ、ｍｘ、ｍｙ、ξｘ及びξｙの値を判定できるような結果として得られる干渉縞強度のセットが明らかになる。
Ｘ線タルボシステム４００において使用される格子の性質は、結果として得られる画像にいくつかの高調波が存在することを意味する。しかし、センサ４４０の影響を含む光学系の変調変換関数（ＭＴＦ）のため、通常、高調波の強度は非常に弱い。従って、これらの高調波は式（４）において無視される。
当該連立方程式の解法
式（４）を解いてａ、ｍｘ、ｍｙ、ξｘ及びξｙを求めるためには、一連の画像を撮像する必要がある。撮像された強度が、式（４）のようにｘ方向及びｙ方向の各々における分離可能な変調により適切にモデル化される場合、５回と少ない撮像を使用して式（４）の連立方程式を解くことが理論上可能である。
位相ステップの値φｘ，ｎ及びφｙ，ｎに依存して、本開示は５回以上のどんな撮像回数も使用する解を提供する。位相ステップが慎重に選択され且つ実験により達成される場合、いくつかの連立方程式に対して閉形式解が存在する。他の例において、閉形式解は未知であるが、式は依然として反復して解くことができる。しかし、反復解は、計算コストが閉形式解よりはるかに高く、確実に収束するために解に対する適切な初期推定値を必要とする。従って、可能な場合は閉形式解を使用することが望ましい。実験により達成された位相ステップが必要な位相ステップと異なる場合、これらの閉形式解では回復されたパラメータに誤差が生じる。しかし、一般に、達成される解は真の解に近似し、既知の位相ステップを使用して解を改善するために使用可能な反復法に対する適切な開始点となる解を形成する。これは、真の位相ステップがいくつかの手段により判定される必要があることを意味する。閉形式解が真の解に近似するため、反復解に必要な反復は非常に少なくなる。従って、閉形式の初期推定値、並びに真の位相ステップの適切な推定値を与えられた場合の反復解の一般的な性質による位相ステップ誤差に対するロバスト性により、この手法は高速画像処理を提供する。
モデルの適合
実際のシステムにおいて使用される格子４１０、４３０が正方形の開口を有さず且つシステム４００の変調伝達関数が高空間周波数成分の振幅を減少するため、厳密には、Ｘ線タルボシステム４００により生成されるモアレインターフェログラムを式（４）のような分離可能な形式で表すことはできない。これらの影響を考慮するために、分離可能な表現を以下の形式に展開して適合する。
式中、交差項の変調は因数ｂにより変更される。これにより、非正方形開口及び光学系のＭＴＦの影響を補正できる（開口は円形であり且つＭＴＦは円対称であると仮定する）。当然、実際のシステムはこれらの仮定に厳密に忠実ではない。しかし、実際のシステムに対してこれらを仮定することにより生じる誤差は非常に小さい。
閉形式解
図１は、閉形式解を提供する方法１００を示す。方法１００は、Ｘ線装置１３９９と協働するコンピュータシステム１３００上で実行可能なソフトウェアとして実現されるのが好ましい。ステップ１１０において、干渉縞インターフェーフェログラム画像が撮像され、格納及び処理のためにコンピュータ１３０１に提供される。いくつかの適用例において、画像はリモートのＸ線装置により撮像され、ネットワーク１３２０を介してコンピュータ通信により画像が提供されてもよい。ステップ１２０において、コンピュータ１３０１はアプリケーションを実行し、画像内の各画素に対して閉形式解を使用することにより、ステップ１１０の撮像から取得された干渉縞インターフェログラム画像（モアレ画像データ）から５つのパラメータａ、ｍｘ、ｍｙ、ξｘ及びξｙを推定する。ステップ１３０において、推定されたパラメータを使用して図４の物体の光路長ψ（ｘ，ｙ）のｘ導関数及びｙ導関数を求めることにより、ステップ１２０からの推定されたパラメータを使用して物体４０２（対象１３２４）の位相画像が再構成される。
撮像ステップ１１０の詳細を図６で更に説明する。ステップ６１０において、第１の干渉縞インターフェログラム画像が撮像される。ステップ６２０において、所望の位相ステップを取得するために、少なくとも１つの格子を空間移動する。理論上、格子Ｇ１（４１０）及びＧ２（４３０）は、Ｘ線方向に対して垂直な平面において「平行移動」の動きで移動可能である。Ｇ１とＧ２との間の相対位置は、移動処理の間、常に維持される。すなわち、それらの間の距離及び角度（モアレを生成するため）は同一のままである。２つの格子（Ｇ１及びＧ２）は、正弦波関数の十分なサンプルを生成するために、複数の小さい規則的ステップ（本開示によると、５つ、６つ、７つ又は８つのステップ）で移動される。物体がＸ線に対して垂直な平面において同様の方法であるが反対方向に移動される場合、同一の効果が達成される。しかし、物体の形状が通常は不規則であるため、物体を移動することは困難である。また、実際に２つの格子を同時に全く同一の量移動することは非現実的である。物体４０２、１３２４がＧ１（４１０）及びＧ２（４３０）に対して平行な平面において移動するように見えるように、図４に不図示であるが図１３Ａに格子４５０として示す追加の格子Ｇ０をＸ線源４０４に近接して構成し且つＸ線に対して垂直な方向にステップ移動するのが好ましい。その後、ステップ６３０において、ＸＴＩシステムは別の干渉縞画像を撮像する。
ステップ６４０において、全ての位相シフトが行われたかがチェックされる。５つの位相シフト画像を撮像するために、格子は少なくとも４回実際に移動される。尚、最初のシフトは０であると考えれば、５つの画像を得るためには単純に５回の移動が必要である。全ての移動が完了すると、ステップ６５０において、パラメータの推定のために、撮像された干渉縞画像の全てと公称（nominal）位相ステップφｘ，ｎ及びφｙ，ｎとがステップ１２０に渡される。
図６に示す画像撮像処理１１０において、異なる撮像回数が適用可能である。ステップ１１０における撮像回数に基づいて、異なる公式が閉形式解を使用してパラメータを推定するためにステップ１２０において使用される。
例１：固定比率を用いる９回の位相シフトによる推定
本例において、全部で９回の撮像が必要とされる。ステップ６４０において、９回全ての位相シフトが行われて９つの位相シフト画像が取得されたかがチェックされる。これらの人工的に導入される９つの位相ステップは、正弦波信号の少なくとも１つのサイクルを範囲に含むために、何らかの特定の値であるように選択される。本例は、ｘ方向の位相ステップφｘ，ｎとｙ方向の位相ステップφｙ，ｎとの間に固定比率を使用する。すなわち、φｙ，ｎ＝３φｘ，ｎとする。これにより、ｘ方向において範囲に含まれるサイクル数が整数である限りｙ方向において範囲に含まれるサイクル数は整数であることが保証される。
本例の場合、式（４）は、積を展開し且つ三角関数の公式を使用することにより、以下のように書き換えられる。
ここで、主成分は、格子の主要な支配（rulings）に関連付けられる強い成分である。これらの主成分は、以下の項に関連付けられる。
式（４）の積を展開することにより、通常は主成分より弱く且つ主成分のベクトル周波数のベクトル和である空間周波数ドメイン内の周波数である交差成分が得られる。これらの交差成分は、以下の項に関連付けられる。
本例において、ｘ方向の位相ステップφｘ，ｎは以下のように選択される。
［０，２π／９，４π／９，２π／３，８π／９，１０π／９，４π／３，１４π／９，１６π／９］（ラジアン）
それに対応して、ｙ方向の位相ステップφｙ，ｎは以下の通りである。
［０，２π／３，４π／３，０，２π／３，４π／３，０，２π／３，４π／３］（ラジアン）
図７は、位相ステップ（φｘ，ｎ，φｙ，ｎ）、ｎ＝０…Ｍ−１を示す。この場合、Ｍは５以上のどんな値も有することができるが、本例ではＭ＝９を使用する。尚、正弦波関数の周期性のため、（φｘ，ｎ，φｙ，ｎ）に対する上記の値は２πにラップ（位相畳みこみ）される。ラップされない場合、この位相シフトパターンは（ｘ，ｙ）平面における直線で実現可能である。実際、これは、単純な機械設計、並びに位相シフト処理におけるヒステリシスにより引き起こされる機械雑音の減少を示す。
式（６）及び式（９）を再度参照すると、φｙ，ｎ＝３φｘ，ｎである場合、φｎ −＝２φｘ，ｎ及びφｎ ＋＝４φｘ，ｎである。これは、ｘ方向において範囲に含まれるサイクル数が整数である場合にφｎ −及びφｎ ＋の双方が整数のサイクル数を範囲に含むことを保証する。換言すると、（ξｙ−ξｘ）成分に対する位相ステップφｎ −は以下の通りである。
［０，４π／９，８π／９，１２π／９，１６π／９，２π／９，６π／９，１０π／９，１４π／９］
また、（ξｙ＋ξｘ）成分に対する位相ステップφｎ ＋は以下の通りである。
［０，８π／９，１６π／９，６π／９，１４π／９，４π／９，１２π／９，２π／９，１０π／９］
式（６）における４つの成分ξｘ、ξｙ、ξｙ−ξｘ及びξｙ＋ξｘの全てに対する位相ステップが整数のサイクル数を範囲に含み、すなわち、
であるため、以下の基底は直交基底である。
式（６）において説明される連立一次方程式の場合、撮像された値ｚｎは、互いに対して直交する９つのベクトルの線形結合である。この連立一次方程式は、撮像した値を直交基底に投影することにより解くことができる。
αｘ、βｘ、αｙ、βｙから以下を回復できる。
α−、β−、α＋、β＋から以下を回復できる。
式（１１）及び式（１２）は、パラメータに対する別個の推定値の２つのセットを提供する。これは、パラメータのノイズを推定するのに有用である。更に、位相値の差分及び和は、ξｘ及びξｙが異なるラッピングアーティファクトを有するためξｘ及びξｙの値をアンラッピングするのに有用である場合がある。
図４の物体４０２の光路長ψ（ｘ，ｙ）のｘ導関数及びｙ導関数は、以下を使用してステップ１３０において求められる。
式中、（ｘ，ｙ）は画素位置であり、ｐｘ及びｐｙはそれぞれ、ｘ方向及びｙ方向における既知のモアレ搬送周期である。
例２：固定比率を用いる８回の位相シフトによる推定
本例において、全部で８回の撮像が必要とされるため、ステップ６４０において、８回全ての位相シフトが行われたかがチェックされる。これらの人工的に導入される８つの位相ステップは、正弦波信号の少なくとも１つのサイクルを範囲に含むために、何らかの特定の値であるように選択される。本例は、ｘ方向の位相ステップφｘ，ｎとｙ方向の位相ステップφｙ，ｎとの間に固定比率を使用する。すなわち、φｙ，ｎ＝３φｘ，ｎとする。この固定比率により、ｘ方向において範囲に含まれるサイクル数が整数である限りｙ方向において範囲に含まれるサイクル数は整数であることが保証される。
式（６）を再度参照すると、φｙ，ｎ＝３φｘ，ｎであるため、φｎ −＝２φｘ，ｎ及びφｎ ＋＝４φｘ，ｎであり、位相ステップは以下の通りである。
φｘ：［０，π／４，π／２，３π／４，π，５π／４，３π／２，７π／４］
φｙ：［０，３π／４，３π／２，π／４，π，７π／４，π／２，５π／４］
φ−：［０，π／２，π，３π／２，０，π／２，π，３π／２］
φ＋：［０，π，０，π，０，π，０，π］
図８は、Ｍ＝８とする位相ステップ（φｘ，ｎ，φｙ，ｎ）、ｎ＝０…Ｍ−１を示す。尚、正弦波関数の周期性を与えられる場合、（φｘ，ｎ，φｙ，ｎ）に対する上記の値は２πにラップされる。ラップされない場合、この位相シフトパターンは（ｘ，ｙ）平面における直線で実現可能である。実際、これは、単純な機械設計、並びに位相シフト処理におけるヒステリシスにより引き起こされる機械雑音の減少を示す。
上記の位相ステップが式（６）において使用され、以下を得られる。
結果として得られた連立方程式（１４）を解くことにより、以下を得られる。
第１の例と同様に、式（１１）を使用してｍｘ、ｍｙ、ξｘ及びξｙに対する第１の推定値セットを判定し且つ式（１２）を使用して第２の推定値セットを判定できる。尚、φ＋に対する位相ステップは０又はπのいずれかであり、これにより全てのｎに対してｓｉｎ（φｎ ＋）＝０になるため、β＋は連立方程式（１４）に存在しない。しかし、連立方程式（１４）は８回の撮像を用いる冗長な連立方程式であるため、このことが最終的な解に及ぼす影響は非常に限定される。
例３：市松模様を用いる８回の位相シフトによる推定
本例において、全部で８回の撮像が必要とされ、この場合も、ステップ６４０において、８回全ての位相シフトが行われたかがチェックされる。本例において適用される位相ステップは、ｘ−ｙ平面において市松模様を形成する。図９は、Ｍ＝８とする位相ステップ（φｘ，ｎ，φｙ，ｎ）、ｎ＝０…Ｍ−１を示す。位相ステップは、以下の通りである。
φｘ：［０，π／２，π，３π／２，０，π／２，π，３π／２］
φｙ：［０，３π／２，π，π／２，π，π／２，０，３π／２］
尚、正弦波関数の周期性のため、（φｘ，ｎ，φｙ，ｎ）に対する上記の値は２πにラップされる。上記の位相ステップを式（６）に挿入することにより、以下を得られる。
式（７）と同様に、以下の代入を行う。
これにより、（１６）の式は以下のようになる。
（１７）の上記の連立方程式は、過剰決定の連立方程式である。雑音を特徴付けるために、推定されたパラメータａ、ｍｘ、ｍｙ、ξｘ、ξｙの２つの別個のセットが取得される。
ランダム誤差の影響を軽減するためにこれらの推定値を平均することにより、以下を得られる。
第１の例と同様に、ｍｘ、ｍｙ、ξｘ、ξｙの値は式（１１）を使用して導き出される。
例４：市松模様を用いる７回の位相シフトによる推定
本例において、全部で７回の撮像が必要とされる。ステップ６４０において、７回全ての位相シフトが行われたかがチェックされる。本例において適用される位相ステップは、ｘ−ｙ平面において市松模様を形成する。図１０は、Ｍ＝７とする位相ステップ（φｘ，ｎ，φｙ，ｎ）、ｎ＝０…Ｍ−１を示す。位相ステップは、以下の通りである。
φｘ：［０，π／２，π，３π／２，０，π／２，π］
φｙ：［０，３π／２，π，π／２，π，π／２，０］
尚、正弦波関数の周期性を与えられる場合、（φｘ，ｎ，φｙ，ｎ）に対する上記の値は２πにラップされる。上記の位相ステップを式（６）に挿入することにより、以下を得られる。
式（７）と同様に、以下の代入を行う。
これにより、（２０）の式は以下のようになる。
（２１）における上記の連立方程式を解くことにより、以下を得られる。
第１の例と同様に、ｍｘ、ｍｙ、ξｘ、ξｙの値は式（１１）を使用して導き出される。
例５：市松模様を用いる６回の位相シフトによる推定
本例において、全部で６回の撮像が必要とされる。ステップ６４０において、６回全ての位相シフトが行われたかがチェックされる。本例において適用される位相ステップは、ｘ−ｙ平面において市松模様を形成する。図１１は、Ｍ＝６とする位相ステップ（φｘ，ｎ，φｙ，ｎ）、ｎ＝０…Ｍ−１を示す。位相ステップは以下の通りである。
φｘ：［０，π／２，π／２，π，π，３π／２］
φｙ：［π／２，０，π，π／２，３π／２，π］
尚、正弦波関数の周期性のため、（φｘ，ｎ，φｙ，ｎ）に対する上記の値は２πにラップされる。上記の位相ステップを式（６）に挿入し、式（７）と同様に以下の代入を行う。
これにより、解かれる連立方程式は以下のようになる。
（２１）における上記の連立方程式を解くことにより、以下を得られる。
第１の例と同様に、ｍｘ、ｍｙ、ξｘ、ξｙの値は式（１１）を使用して導き出される。
分離可能な変調を用いる非線形解
式（４）におけるパラメータは、同様にコンピュータ１３０１により実現されてもよい図２に示す反復法２００を使用して更に推定可能である。図２を参照すると、方法２００はステップ１１０から開始し、干渉縞インターフェログラム画像が図１と同一の方法で撮像される。ステップ２２０において、式（４）におけるパラメータの推定値が反復解から判定される。
ステップ１３０において、図４の物体の光路長ψ（ｘ，ｙ）のｘ導関数及びｙ導関数が求められ、方法２００はステップ２９９において終了する。再構成は、回復されたデータ（各画素において個別に回復される）を画像として処理することにより実行されてもよい。搬送位相が除去されると物体の特徴が更に可視であるため、微分光路長画像は一般に使用されるものである。
ステップ２２０の詳細を図３に説明する。図３では、ステップ３１０において、位相及び変調パラメータａ、ｍｘ、ｍｙ、ξｘ及びξｙの初期推定値が取得される。５つのパラメータの初期推定値は、第１の例と同一であってもよく、あるいは、本例において実現される方法が反復法であるため、単にランダムな値であってもよい。各画素に対するパラメータ推定値は、当該画素に対するシミュレーションされた（すなわち、反復された）干渉縞強度のセットを生成し且つシミュレーションされた干渉縞強度のセットと当該画素における撮像された干渉縞強度のセットとを比較することにより修正される。ステップ３２０において、位相ステップ誤差を補正するために、位相推定法が使用される。位相ステップ誤差は、各モアレ画像の撮像における格子４１０及び４３０により形成される交差回折格子の機械的な移動又はステップ移動により発生する。ステップ３３０において、ステップ３２０からの補正された位相ステップ及びステップ３１０からの初期推定値を使用して、パラメータの最適解を反復して求める。
式（６）、式（７）及び式（８）を再度参照すると、次式を得られる。
以下の項を展開する。
最終的に、次式を得られる。
式中、ｎ＝０，．．．，Ｍ−１とする。
図２の方法２００のステップ１１０において使用される位相ステップ数に依存して、Ｍの値は５以上のどんな整数でもよい。
例７：Ｍ回の位相シフトを用いる非線形復調
本例において、全部で９つの画像がステップ１１０において撮像される。式（２５）は、以下のように更に簡潔に書くことができる。
式中、ベクトルCは推定されるパラメータを表す。
残差関数を定義するのが好都合である。
式（２６）において、Zn(C)は、与えられた位相ステップ値（φｘ，ｎ，φｙ，ｎ）、ｎ＝０，．．．，Ｍ−１を用いて式（２５）により説明されるモデル化された画素の強度を示し、ｚｎは、ステップ１１０の出力である撮像された画素の強度を表す。f(C)を最小にするために、ガウス・ニュートンアルゴリズムを適用して平均二乗残差を最小にする。
これにより、以下のような反復解を得られる。
式中、Jはf(C)のM×５ヤコビ行列である。すなわち、
Jの列は以下の通りである。
何らかの初期推定値Coから開始する場合、ガウス・ニュートン法は、推定値が収束するまで、各ステップにおいてパラメータベクトルの推定値を更新する（f(Ck)＜ε、式中、εは何らかの小さい数である）。
上述の反復ガウス・ニュートンアルゴリズムは計算コストが高い場合があり、あるいは、特定の状況において収束しない場合がある。従って、Cの適切な初期推定値を選択することが重要である。この初期推定値は、ヒューリスティックに選択可能である。閉形式解と異なり、式（２８）のヤコビ行列Jが正則である限り、特定の位相シフトパターンは不要である。
それに対して、９回の撮像の場合、例１で説明したように、閉形式解が存在する。その場合、ガウス・ニュートンアルゴリズムにおける初期推定値は、例１における閉形式解を介して取得可能である。一般に、閉形式解からの初期推定値は、ヒューリスティックに選択される初期値よりはるかに正確である。この場合、式（２８）における位相ステップは例１におけるパターンを有する。
更に、ステップ３２０において、式（２８）における導入された位相ステップ（φｘ，ｎ，φｙ，ｎ）は、公称位相ステップと実験における実際の位相ステップとの間の誤差を減少するために、位相推定法を使用して補正される。
ステップ３２０において、どんな位相推定法が位相ステップ（φｘ，ｎ，φｙ，ｎ）を補正するために使用されてもよい。例えば、位相シフト画像のグループにおける個々の搬送の搬送周波数は、フーリエ変換のピークの位置を測定することにより判定可能であり、各位相シフト画像の相対位相は、それらのピークの複素位相を測定することにより判定可能である。
例８：８回の位相シフトを用いる非線形復調
本例では、全部で８つの画像が１１０において撮像される。例７と同様に、平均二乗センス（sense）が最小であるという意味で最適な推定値を取得するために、ガウス・ニュートンアルゴリズムが使用される。反復アルゴリズムに対するこの初期推定値C0は、ヒューリスティックに選択可能である。閉形式解と異なり、式（２８）のヤコビ行列Jが正則である限り、特定の位相シフトパターンは不要である。
８回の撮像の場合、例２及び例３で説明したように、閉形式解が存在する。その場合、初期推定値C0は、例２又は例３における閉形式解を介して取得可能である。一般に、閉形式解からの初期推定値は、ヒューリスティックに選択される初期値よりはるかに正確である。この場合、式（２８）における位相ステップは例２又は例３におけるパターンを有する。
更に、ステップ３２０において、式（２８）における導入された位相ステップ（φｘ，ｎ，φｙ，ｎ）は、公称位相ステップと実験における実際の位相ステップとの間の誤差を減少するために、位相推定法を使用して補正される。
例９：７回の位相シフトを用いる非線形復調
本例では、全部で７つの画像がステップ１１０において撮像される。例７と同様に、平均二乗センス（sense）が最小であるという意味で最適な推定値を取得するために、ガウス・ニュートンアルゴリズムが使用される。反復アルゴリズムに対する初期推定値C0は、ヒューリスティックに選択可能である。閉形式解と異なり、式（２８）のヤコビ行列Jが正則である限り、特定の位相シフトパターンは不要である。
７回の撮像の場合、例４で説明したように、閉形式解が存在する。その場合、初期推定値C0は、例４における閉形式解を介して取得可能である。一般に、閉形式解からの初期推定値は、ヒューリスティックに選択される初期値よりはるかに正確である。この場合、式（２８）における位相ステップは例４におけるパターンを有する。
更に、ステップ３２０において、式（２８）における導入された位相ステップ（φｘ，ｎ，φｙ，ｎ）は、公称位相ステップと実験における実際の位相ステップとの間の誤差を減少するために、位相推定法を使用して補正される。
例１０：６回の位相シフトを用いる非線形復調
本例では、全部で６つの画像が１１０において撮像される。例７と同様に、平均二乗センス（sense）が最小であるという意味で最適な推定値を取得するために、ガウス・ニュートンアルゴリズムが使用される。反復アルゴリズムに対するこの初期推定値C0は、ヒューリスティックに選択可能である。閉形式解と異なり、式（２８）のヤコビ行列Jが正則である限り、特定の位相シフトパターンは不要である。
６回の撮像の場合、例５で説明したように、閉形式解が存在する。その場合、初期推定値C0は、例５における閉形式解を介して取得可能である。一般に、閉形式解からの初期推定値は、ヒューリスティックに選択される初期値よりはるかに正確である。この場合、式（２８）における位相ステップは例５におけるパターンを有する。
更に、ステップ３２０において、式（２８）における導入された位相ステップ（φｘ，ｎ，φｙ，ｎ）は、公称位相ステップと実験における実際の位相ステップとの間の誤差を減少するために、位相推定法を使用して補正される。
例１１：５回の位相シフトを用いる非線形復調
本例において、全部で５回の撮像が必要とされる。ステップ６４０において、５回全ての位相シフトが行われたかがチェックされる。これらの人工的に導入される５つの位相シフトは、正弦波信号の少なくとも１つのサイクルを範囲に含むために、何らかの特定の値であるように選択される。本例は、ｘ方向の位相ステップφｘ，ｎとｙ方向の位相ステップφｙ，ｎとの間に固定比率を使用する。すなわち、φｙ，ｎ＝３φｘ，ｎとする。これは、ｘ方向において範囲に含まれるサイクル数が整数である限りｙ方向において範囲に含まれるサイクル数が整数であることを保証する。
図１２は、Ｍ＝５とする位相ステップ（φｘ，ｎ，φｙ，ｎ）、ｎ＝０…Ｍ−１を示す。位相ステップは以下の通りである。
φｘ：［０，２π／５，４π／５，６π／５，８π／５］
φｙ：［０，６π／５，２π／５，８π／５，４π／５］
換言すると、位相ステップは以下のように記述可能である。
式中、ｎ＝０，．．．，４である。尚、（φｘ，ｎ，φｙ，ｎ）に対する上記の値は、正弦波関数の周期性を与えられる場合、２πにラップされる。ラップされない場合、この位相シフトパターンは（ｘ，ｙ）平面における直線で実現可能である。実際、これは単純な機械設計、並びに位相シフト処理におけるヒステリシスにより引き起こされる機械雑音の減少を示す。
式（２９）を使用して式（６）のｓｉｎ（φｎ −）、ｓｉｎ（φｎ ＋）、ｃｏｓ（φｎ −）及びｃｏｓ（φｎ ＋）を置換する。
（３０）の上記の連立方程式を解くことにより、［ａ，ａ（αｘ＋α＋），ａ（βｘ−β＋），ａ（αｙ＋α−），ａ（βｙ−β−）］の値が算出される。
式（７）及び式（８）における定義を使用する場合、［ａ，ｍｘ，ｍｙ，ξｘ，ξｙ］の値は非線形連立方程式を解くことにより回復可能である。この非線形連立方程式は、ガウス・ニュートンアルゴリズム等の種々の手法を用いて解くことができる。ガウス・ニュートンアルゴリズムに対する初期推定値は、ヒューリスティックに選択可能である。
分離できない変調を用いる非線形解
前述の例は、式（４）で説明される変調モデルに基づき、この場合、ｘ及びｙ方向における変調は分離可能である。図５の交差格子の異なる構成を用いる場合、干渉縞インターフェログラム画像に対して異なる変調モデルを使用する方が精度が高い場合がある。
例えば異なる格子設計は、式（４）の積を展開した場合の以下の交差項の強度を変更する。
例１２：変化する交差項を用いる非線形復調
（３１）における交差項の振幅がａｂであると仮定すると、式（５）で説明したように、次式を得られる。
式中、ｂは０〜１で変化する。式は、以下のように書き換えられてもよい。
例７と同様に、平均二乗センス（sense）が最小であるという意味で最適な推定値を取得するために、ガウス・ニュートンアルゴリズムが使用される。ガウス・ニュートンアルゴリズムに対する初期推定値は、ヒューリスティックに選択されてもよく、あるいは、例１〜例５で説明した閉形式解が使用されてもよい。
分離できない場合、新規のヤコビ行列Jはf(C)のＭ×６ヤコビ行列である。すなわち、
Jの列は以下の通りである。
図２のステップ１１０で使用される位相ステップ数に依存して、Ｍの値は５以上のどんな整数であってもよい。
例１３：固定された交差項を用いる非線形復調
特定の格子構成を用いる場合、式（３２）において交差項の振幅ｂを推定できるため、交差項の相対強度は既知であると仮定できる。本例において、推定されるパラメータは例７と同一のままである。例７と同様に、平均二乗センス（sense）が最小であるという意味で最適な推定値を取得するために、ガウス・ニュートンアルゴリズムが使用される。ガウス・ニュートンアルゴリズムに対する初期推定値は、ヒューリスティックに選択されてもよく、あるいは、例１〜例５で説明した閉形式解が使用されてもよい。
本例に対するヤコビ行列Jは、f(C)のＭ×５ヤコビ行列である。すなわち、
Jの列は以下の通りである。
図２のステップ１１０で使用される位相ステップ数に依存して、Ｍの値は５以上のどんな整数であってもよい。本例では、干渉縞インターフェログラム画像が両方の次元で変化する交差項を有する場合に、交差項の振幅を、格子構造に関するデータを用いて個別に推定することができ且つ決定した振幅を用いてパラメータ推定値を修正することを示した。
結論
例は、物体により生じる吸収パラメータ（ａ）、振幅変調パラメータ（ｍｘ，ｍｙ）及び位相変調パラメータ（ξｘ，ξｙ）である位相及び変調パラメータを推定するために、位相シフト技術を使用して取得される多くのモアレ干渉縞パターン画像の処理を提供する。これらのパラメータは、画像シーケンスにおける画素位置毎に個別に回復される。画像と見なされ且つ式（１３）に従ってモアレパターンの搬送周波数を補正される場合、これらは、物体（４０２、１３２４）の微分光路長画像（∂ψ／∂ｘ、∂ψ／∂ｙ）を提供し、個々のＸ線画像又は干渉縞パターン画像から識別できない物体の細部を明らかにする可能性を有する。
各パラメータは、特定の画像により対応する表現を提供する。特定の画像は、以下を含む。
（ｉ） 式（７）の項（ξｘ）により示されるモアレ画像の各（ｘ，y）位置における第１の微分位相である第１の位相画像
（ｉｉ） 式（７）の項（ξｙ）により示されるモアレ画像の各（ｘ，ｙ）位置における第２の微分位相である第２の位相画像
（ｉｉｉ） 式（７）の項（ｍｘ）により示されるモアレ画像の各（ｘ，ｙ）位置における第１の変調である第１の変調画像
（ｉｖ） 式（７）の項（ｍｙ）により示されるモアレ画像の各（ｘ，ｙ）位置における第２の変調である第２の変調画像
（ｖ） 式（６）の項（ａ）により示される物体のＸ線吸収を表す背景（吸収）画像
本開示によると、少なくとも５つであり且つ８つ以下であるモアレ画像がステップ１１０において撮像され、ステップ１０２において処理される。Ｘ線の応用の場合、撮像画像数の減少は、物体（生体である可能性がある）へのＸ線照射が減少し、また、処理される必要のあるデータセットの総数が減少することで一般に通信及び計算に関する時間及びコストが減少されるため、有利である。
特定の画像、特に位相画像は、従来のＸ線法で識別できない物体の細部を明らかにするために単独で使用されてもよいが、物体の位相画像を完全に再構成するために特定の画像のセットが一括して使用されてもよい。画像は、表示装置１３１４上に表示されてもよい。
説明した処理ステップの場合、迅速な解を導出するために閉形式手法を使用できるが、これは位相ステップ誤差が生じやすい。反復手法は非常に遅いが、その性質上、位相ステップ誤差を調整できる。２つの手法は組み合わされて使用されてもよい。その場合、最初に閉形式手法を使用して反復手法の開始点が迅速に取得された後、反復手法により、不必要な反復を回避して、改善された解が提供される。その結果、合理的な計算リソースを使用して特定の画像の最適なセットが取得される。
説明した構成は、コンピュータ及びデータ処理産業に適用可能であり、特に、物体のＸ線撮影に適用可能である。構成は、損傷等を識別するための軟部組織の撮影に特に適用可能であり、画像処理により補足される安価なＸ線機器を使用するＭＲＩ及びＣＡＴスキャンの代替手段を提供する可能性を有する。
上記の説明は本発明のいくつかの実施形態のみを示し、変更及び／又は変形は本発明の主旨の範囲から逸脱することなく行われる。実施形態は例示であって、限定するものではない。

Claims (22)

  1. ２次元の干渉縞インターフェログラム画像のセットから位相画像を再構成する方法であって、
    ８つ以下の撮像された２次元の干渉縞インターフェログラム画像を含む前記干渉縞インターフェログラム画像のセットを使用して、吸収パラメータ、２次元振幅変調パラメータ、及び２次元位相変調パラメータの推定値をパラメータ推定値として、閉形式解から取得するステップと、
    前記パラメータ推定値を修正するステップと、
    前記パラメータ推定値を使用して、前記位相画像を再構成するステップと
    を備えることを特徴とする方法。
  2. 各画素に対する前記パラメータ推定値は、該画素に対する干渉縞強度のシミュレーションされたセットを生成し且つ前記干渉縞強度のシミュレーションされたセットと該画素における撮像された干渉縞強度のセットとを比較することにより修正されることを特徴とする請求項１記載の方法。
  3. 各画素に対する前記パラメータ推定値は、該画素に対する干渉縞強度のシミュレーションされたセットを生成し且つ前記干渉縞強度のシミュレーションされたセット及び撮像された干渉縞強度のセットにおける対応する画素の誤差を最小化することにより修正されることを特徴とする請求項１記載の方法。
  4. ２次元の干渉縞インターフェログラム画像のセットから位相画像を再構成する方法であって、
    ８つ以下の撮像された２次元の干渉縞インターフェログラム画像を含む前記干渉縞インターフェログラム画像のセットを使用して、吸収パラメータ、２次元振幅変調パラメータ、及び２次元位相変調パラメータの推定値をパラメータ推定値として、閉形式解から取得するステップと、
    前記パラメータ推定値を使用して、前記位相画像を再構成するステップと
    を備え、
    前記取得するステップは、
    前記パラメータの初期推定値を取得することと、
    前記パラメータにおける位相ステップ誤差を補正するために位相推定法を使用することと、
    前記補正された位相ステップ及び前記初期推定値を使用して、反復法により前記パラメータに対する反復解を取得することと
    を含むことを特徴とする方法。
  5. ２次元の干渉縞インターフェログラム画像のセットから位相画像を再構成する方法であって、
    ５つの撮像された２次元の干渉縞インターフェログラム画像を含む前記干渉縞インターフェログラム画像のセットを使用して、吸収パラメータ、２次元振幅変調パラメータ、及び２次元位相変調パラメータの推定値をパラメータ推定値として、閉形式解から取得するステップと、
    前記パラメータ推定値を使用して、前記位相画像を再構成するステップと
    を備えることを特徴とする方法。
  6. 前記画像の位相ステップは、φｘ：［０，２π／５，４π／５，６π／５，８π／５］及びφｙ：［０，６π／５，２π／５，８π／５，４π／５］であることを特徴とする請求項５記載の方法。
  7. ２次元の干渉縞インターフェログラム画像のセットから位相画像を再構成する方法であって、
    ６つの撮像された２次元の干渉縞インターフェログラム画像を含む前記干渉縞インターフェログラム画像のセットを使用して、吸収パラメータ、２次元振幅変調パラメータ、及び２次元位相変調パラメータの推定値をパラメータ推定値として、閉形式解から取得するステップと、
    前記パラメータ推定値を使用して、前記位相画像を再構成するステップと
    を備えることを特徴とする方法。
  8. 前記画像の位相ステップは、φｘ：［０，π／２，π／２，π，π，３π／２］及びφｙ：［π／２，０，π，π／２，３π／２，π］であることを特徴とする請求項７記載の方法。
  9. ２次元の干渉縞インターフェログラム画像のセットから位相画像を再構成する方法であって、
    ８つ以下の撮像された２次元の干渉縞インターフェログラム画像を含む前記干渉縞インターフェログラム画像のセットを使用して、吸収パラメータ、２次元振幅変調パラメータ、及び２次元位相変調パラメータの推定値をパラメータ推定値として、閉形式解から取得するステップと、
    前記パラメータ推定値を使用して、前記位相画像を再構成するステップと
    を備え、
    前記セットは、２次元の干渉縞インターフェログラム画像を７つ含み、前記画像の位相ステップは、φｘ：［０，π／２，π，３π／２，０，π／２，π］及びφｙ：［０，３π／２，π，π／２，π，π／２，０］であることを特徴とする方法。
  10. ２次元の干渉縞インターフェログラム画像のセットから位相画像を再構成する方法であって、
    ８つ以下の撮像された２次元の干渉縞インターフェログラム画像を含む前記干渉縞インターフェログラム画像のセットを使用して、吸収パラメータ、２次元振幅変調パラメータ、及び２次元位相変調パラメータの推定値をパラメータ推定値として、閉形式解から取得するステップと、
    前記パラメータ推定値を使用して、前記位相画像を再構成するステップと
    を備え、
    前記セットは、２次元の干渉縞インターフェログラム画像８つを含み、前記画像の位相ステップは、φｘ：［０，π／２，π，３π／２，０，π／２，π，３π／２］及びφｙ：［０，３π／２，π，π／２，π，π／２，０，３π／２］であることを特徴とする方法。
  11. 前記干渉縞インターフェログラム画像は、両方の次元で変化する交差項を有することを特徴とする請求項１に記載の方法。
  12. 更に、
    前記交差項の振幅を、格子構造に関するデータを用いて個別に推定するステップを備え、
    決定した振幅は、前記パラメータ推定値を修正するために使用される
    ことを特徴とする請求項１１に記載の方法。
  13. ２次元の干渉縞インターフェログラム画像のセットから位相画像を再構成する方法であって、
    前記画像のセットを提供されるステップと、
    前記画像に関連する交差成分と前記画像に関連する主成分とを関連付けることで提供された、８つ以下の画像から前記画像を再構成可能な再構成技術を使用して、提供された前記画像のセットから前記位相画像を再構成するステップと
    を備えることを特徴とする方法。
  14. コンピュータに、請求項１乃至１３の何れか１項に記載の方法の各ステップを実行させるためのコンピュータプログラム。
  15. ２次元の干渉縞インターフェログラム画像のセットから位相画像を再構成する演算装置であって、
    ８つ以下の２次元の干渉縞インターフェログラム画像を含む干渉縞インターフェログラム画像のセットを使用して、吸収パラメータ、２次元振幅変調パラメータ、及び２次元位相変調パラメータの推定値をパラメータ推定値として、閉形式解から取得する手段と、
    前記パラメータ推定値を修正する手段と、
    前記パラメータ推定値を使用して、前記位相画像を再構成する手段と
    を備えることを特徴とする演算装置。
  16. ２次元の干渉縞インターフェログラム画像のセットから位相画像を再構成する演算装置であって、
    ８つ以下の撮像された２次元の干渉縞インターフェログラム画像を含む前記干渉縞インターフェログラム画像のセットを使用して、吸収パラメータ、２次元振幅変調パラメータ、及び２次元位相変調パラメータの推定値をパラメータ推定値として、閉形式解から取得する手段と、
    前記パラメータ推定値を使用して、前記位相画像を再構成する手段と
    を備え、
    前記取得する手段は、
    前記パラメータの初期推定値を取得し、
    前記パラメータにおける位相ステップ誤差を補正するために位相推定法を使用し、
    前記補正された位相ステップ及び前記初期推定値を使用して、反復法により前記パラメータに対する反復解を取得する
    ことを特徴とする演算装置。
  17. ２次元の干渉縞インターフェログラム画像のセットから位相画像を再構成する演算装置であって、
    前記画像に関連する交差成分と前記画像に関連する主成分とを関連付けることで提供された、８つ以下の画像から前記画像を再構成可能な再構成技術を使用して、提供された前記画像のセットから前記位相画像を再構成することを特徴とする演算装置。
  18. 請求項１５乃至１７の何れか１項に記載の演算装置と、交差回折格子を有するＸ線干渉計システムと、を備えることを特徴とする装置。
  19. 前記演算装置は、前記位相画像を表示可能な表示装置を備えることを特徴とする請求項１８に記載の装置。
  20. ２次元の干渉縞インターフェログラム画像のセットから位相画像を再構成する方法であって、
    撮像された２次元の干渉縞インターフェログラム画像を含む前記干渉縞インターフェログラム画像のセットを使用して、吸収パラメータ、２次元振幅変調パラメータ、及び２次元位相変調パラメータの推定値をパラメータ推定値として、閉形式解から取得するステップと、
    前記パラメータ推定値を修正するステップと、
    前記パラメータ推定値を使用して、前記位相画像を再構成するステップと
    を備えることを特徴とする方法。
  21. ２次元の干渉縞インターフェログラム画像のセットから位相画像を再構成する方法であって、
    撮像された２次元の干渉縞インターフェログラム画像を含む前記干渉縞インターフェログラム画像のセットを使用して、吸収パラメータ、２次元振幅変調パラメータ、及び２次元位相変調パラメータの推定値をパラメータ推定値として、閉形式解から取得するステップと、
    前記パラメータ推定値を使用して、前記位相画像を再構成するステップと
    を備え、
    前記取得するステップは、
    前記パラメータの初期推定値を取得することと、
    前記パラメータにおける位相ステップ誤差を補正するために位相推定法を使用することと、
    前記補正された位相ステップ及び前記初期推定値を使用して、反復法により前記パラメータに対する反復解を取得することと
    を含むことを特徴とする方法。
  22. ２次元の干渉縞インターフェログラム画像のセットから位相画像を再構成する方法であって、
    前記画像のセットを提供されるステップと、
    前記画像に関連する交差成分と前記画像に関連する主成分とを関連付けることで提供された、画像から前記画像を再構成可能な再構成技術を使用して、提供された前記画像のセットから前記位相画像を再構成するステップと
    を備えることを特徴とする方法。
JP2013266132A 2012-12-24 2013-12-24 位相画像を再構成する方法、コンピュータプログラム、位相画像を再構成する演算装置 Expired - Fee Related JP5739976B2 (ja)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
AU2012268876 2012-12-24
AU2012268876A AU2012268876A1 (en) 2012-12-24 2012-12-24 Non-linear solution for 2D phase shifting

Publications (2)

Publication Number Publication Date
JP2014121614A JP2014121614A (ja) 2014-07-03
JP5739976B2 true JP5739976B2 (ja) 2015-06-24

Family

ID=50064346

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2013266132A Expired - Fee Related JP5739976B2 (ja) 2012-12-24 2013-12-24 位相画像を再構成する方法、コンピュータプログラム、位相画像を再構成する演算装置

Country Status (4)

Country Link
US (1) US20140177790A1 (ja)
EP (1) EP2746754A2 (ja)
JP (1) JP5739976B2 (ja)
AU (1) AU2012268876A1 (ja)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015038793A1 (en) * 2013-09-12 2015-03-19 The United States Of America, As Represented By The Secretary, Department Of Health & Human Services Demodulation of intensity modulation in x-ray imaging
AU2013260650B2 (en) * 2013-11-20 2015-07-16 Canon Kabushiki Kaisha Rotational phase unwrapping
CN104501781A (zh) * 2014-11-18 2015-04-08 中国科学院大学 基于已知先验信息约束的叠层成像技术
JP2016106721A (ja) * 2014-12-03 2016-06-20 キヤノン株式会社 画像処理装置および画像処理方法
JP6418066B2 (ja) 2015-05-21 2018-11-07 コニカミノルタ株式会社 位相画像処理方法、位相画像処理装置、画像処理装置及び画像処理プログラム
WO2016207423A1 (en) * 2015-06-26 2016-12-29 Koninklijke Philips N.V. Robust reconstruction for dark-field and phase contrast ct
CN105629695B (zh) * 2016-01-13 2019-01-18 河北工程大学 一种基于相位层叠衍射的全息成像方法
CN106618623B (zh) * 2017-01-11 2019-08-30 合肥工业大学 一次曝光的硬x射线光栅干涉仪的成像方法
US10897560B2 (en) * 2017-06-07 2021-01-19 Rambus Inc. Imaging devices and methods for reducing image artifacts
CN109186497B (zh) * 2018-10-30 2019-11-26 北京理工大学 牛顿迭代的数字莫尔移相干涉法中频谱混叠噪声去除方法
WO2020133531A1 (zh) * 2018-12-29 2020-07-02 深圳大学 X射线双相位光栅相衬成像系统
JP7381827B2 (ja) * 2019-05-30 2023-11-16 株式会社トプコン 光干渉測定装置および光干渉測定方法

Family Cites Families (28)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7341348B2 (en) * 2003-03-25 2008-03-11 Bausch & Lomb Incorporated Moiré aberrometer
WO2010050483A1 (ja) * 2008-10-29 2010-05-06 キヤノン株式会社 Ｘ線撮像装置およびｘ線撮像方法
DE102009004702B4 (de) * 2009-01-15 2019-01-31 Paul Scherer Institut Anordnung und Verfahren zur projektiven und/oder tomographischen Phasenkontrastbildgebung mit Röntgenstrahlung
US9348067B2 (en) * 2009-06-16 2016-05-24 Koninklijke Philips N.V. Tilted gratings and method for production of tilted gratings
JP2011108696A (ja) * 2009-11-13 2011-06-02 Nikon Corp 波面計測方法及び装置、並びに露光方法及び装置
JP5789613B2 (ja) * 2009-12-10 2015-10-07 コーニンクレッカ フィリップス エヌ ヴェ オンザフライ位相ステッピングを備えた非平行な格子装置、ｘ線システム及び使用方法
JP5900324B2 (ja) * 2010-03-18 2016-04-06 コニカミノルタ株式会社 Ｘ線撮影システム
JP2012090944A (ja) * 2010-03-30 2012-05-17 Fujifilm Corp 放射線撮影システム及び放射線撮影方法
JP2012090945A (ja) * 2010-03-30 2012-05-17 Fujifilm Corp 放射線検出装置、放射線撮影装置、放射線撮影システム
DE102010017425A1 (de) * 2010-06-17 2011-12-22 Karlsruher Institut für Technologie Geneigte Phasengitterstrukturen
BR112013009248A2 (pt) * 2010-10-19 2019-09-24 Koninl Philips Electronics Nv difração de reticulação para a geração de imagens por contraste de fase diferencial de raios x, disposição de detector de um sistema de raios x para a geração de imagens por contraste de fase de um objeto, sistema de geração de umagens médica de raios x para a geração de imagens por contraste de fase diferencial, método para a geração de imagens por contraste de fase diferencial, elemento de programa de computador para o controle de um aparelho e meio legível por computador
JP5331940B2 (ja) * 2010-10-27 2013-10-30 富士フイルム株式会社 放射線撮影システム及び放射線画像生成方法
EP2633813B1 (en) * 2010-10-29 2015-02-25 FUJIFILM Corporation Phase contrast radiation imaging device
JP2012115576A (ja) * 2010-12-02 2012-06-21 Fujifilm Corp 放射線画像検出装置、放射線撮影装置、放射線撮影システム
JP2012143537A (ja) * 2010-12-24 2012-08-02 Fujifilm Corp 放射線撮影システム及び画像処理方法
JP5915645B2 (ja) * 2011-03-23 2016-05-11 コニカミノルタ株式会社 医用画像表示システム
JP2012200567A (ja) * 2011-03-28 2012-10-22 Fujifilm Corp 放射線撮影システム及び放射線撮影方法
WO2012144317A1 (ja) * 2011-04-20 2012-10-26 富士フイルム株式会社 放射線撮影装置及び画像処理方法
JP2012236005A (ja) * 2011-04-26 2012-12-06 Fujifilm Corp 放射線撮影装置
JP2013063099A (ja) * 2011-09-15 2013-04-11 Canon Inc Ｘ線撮像装置
JP5475737B2 (ja) * 2011-10-04 2014-04-16 富士フイルム株式会社 放射線撮影装置及び画像処理方法
CN104066375B (zh) * 2012-01-24 2017-08-11 皇家飞利浦有限公司 多方向相衬x射线成像
JP5857800B2 (ja) * 2012-03-01 2016-02-10 コニカミノルタ株式会社 関節撮影装置及び撮影対象固定ユニット
JP2014090967A (ja) * 2012-11-06 2014-05-19 Canon Inc Ｘ線撮像装置
AU2012268882B2 (en) * 2012-12-24 2015-07-09 Canon Kabushiki Kaisha Estimating phase for phase-stepping algorithms
US9439613B2 (en) * 2013-02-12 2016-09-13 The Johns Hopkins University System and method for phase-contrast X-ray imaging
AU2013260650B2 (en) * 2013-11-20 2015-07-16 Canon Kabushiki Kaisha Rotational phase unwrapping
AU2013273822A1 (en) * 2013-12-23 2015-07-09 Canon Kabushiki Kaisha Modulation guided phase unwrapping

Also Published As

Publication number Publication date
US20140177790A1 (en) 2014-06-26
AU2012268876A1 (en) 2014-07-10
EP2746754A2 (en) 2014-06-25
JP2014121614A (ja) 2014-07-03

Similar Documents

Publication Publication Date Title
JP5739976B2 (ja) 位相画像を再構成する方法、コンピュータプログラム、位相画像を再構成する演算装置
JP5766259B2 (ja) 逆リース変換による差分画像の合成
US20140369623A1 (en) Estimating phase for phase-stepping algorithms
JP2016077904A (ja) 撮像方法、画像処理装置、コンピュータ可読媒体、方法、装置およびシステム
JP2016123853A (ja) 磁気共鳴イメージング装置及び画像処理装置
JP2017506925A (ja) 微分位相コントラスト撮像からの位相回復
Dabbech et al. Cygnus A super-resolved via convex optimization from VLA data
AU2013260650B2 (en) Rotational phase unwrapping
KR101727836B1 (ko) 영상 복원 장치 및 방법 및 장치
Servin et al. Modern fringe pattern analysis in interferometry
Milanfar A model of the effect of image motion in the radon transform domain
US8538101B2 (en) Image reconstruction
EP2059904A1 (fr) Procede d'imagerie informatise permettant une reconstruction tridimensionnelle a partir d'images radiographiques bidimensionnelles; dispositif de mise en oeuvre
JP2016087473A (ja) ２次元のフリンジ模様の復調において軸外周波数を低減するための非線形処理の方法、記憶媒体、および撮影システム
Qin et al. Video superresolution reconstruction based on subpixel registration and iterative back projection
Louis Combining image reconstruction and image analysis with an application to two-dimensional tomography
Vargas et al. Error analysis of the principal component analysis demodulation algorithm
JP6495945B2 (ja) 微分位相コントラストイメージングにおけるスペクトル位相アンラッピングのためのバイアスフリー正則化
KR101803139B1 (ko) 3d 홀로그램 왜곡 보정 방법 및 장치
US20160203601A1 (en) Method and device for restoring scan image and recording medium for storing same
Salvado et al. Removal of local and biased global maxima in intensity-based registration
US20230153968A1 (en) Bias correction for dark-field imaging based on sliding window phase retrieval
Dogan et al. Localization of point sources for systems governed by the wave equation
Gong et al. An experimental comparison of super-resolution reconstruction for image sequences
JP2016061608A (ja) 画像処理方法、画像処理装置、撮像システム

Legal Events

Date Code Title Description
A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20141017

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20141201

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20150127

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20150327

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20150424

LAPS Cancellation because of no payment of annual fees