説明

位相情報の解析方法、該位相情報の解析プログラム、記憶媒体、X線位相イメージング装置

【課題】窓フーリエ変換法を用いた解析において、解像度の更なる向上を図ることが可能となる位相情報の解析方法等を提供する。
【解決手段】光あるいはX線を含む波長の波の干渉によるモアレ像の周期的パターンを解析し、位相波面を含む位相情報を取得する位相情報の解析方法であって、
前記モアレ像の周期的なパターンの少なくとも一部を、窓関数によって窓フーリエ変換する工程と、
前記窓フーリエ変換されたモアレ像における、位相情報を担うスペクトルの情報と、前記位相情報を担うスペクトルの情報に重畳している位相情報を担っていないスペクトルの情報とを解析的に計算する工程と、
前記位相情報を担っていないスペクトルの情報を、前記位相情報を担うスペクトルの情報から分離し、前記位相波面を含む位相情報を取得する工程と、を有する構成とする。

【発明の詳細な説明】
【技術分野】
【0001】
本発明は、位相情報の解析方法、該位相情報の解析プログラム、記憶媒体、X線位相イメージング装置に関する。
特に、任意の位相波面をもつ光などの入射波を干渉させて作成したモアレ(干渉縞もしくは強度パターン)等の周期的なパターンからもとの入射波の位相波面もしくは位相波面の微分を計算する技術に関するものである。
【背景技術】
【0002】
光やX線を含む様々な波長の波を用いて干渉を生じさせて被検体の形状計測に用いる技術が知られている。
これらの計測技術では位相波面の揃った(コヒーレントな)入射光を被検体に対し照射し反射もしくは透過させる。
この反射光もしくは透過光は被検体の形状や組成によって波面が変化する事が知られている。
そこで、この変化を何らかの方法を用いて干渉させることにより、モアレ像(これは干渉縞とも称されるが、以下ではこれをモアレと記す。)に変換し、そのパターンを解析する。これによって変化した位相波面もしくはその位相波面の微分像(位相微分像)を計算することが可能となる。
これらの技術の適用例として代表的な例としてはレンズなどの形状測定を行う波面測定技術などがある。
また、近年は医療分野でもX線を用いた技術として、X線位相イメージング手法が知られている。
この手法は入射X線を被写体に透過させた際に、被写体の屈折率の違いによる位相差をモアレ縞などの周期的パターンを利用して取得する技術である。
被写体内の構成物質はそれぞれ異なる屈折率を有するため、波面の変化にそれに応じた特徴が現れる。そこでその位相波面を干渉などによって検出する。
【0003】
これら干渉によって得られた強度パターンから元の波面もしくは入射光の位相波面の変化を計算する技術が位相回復法である。
位相回復法には幾つかの種類があるが、その一つに窓フーリエ変換法がある(非特許文献1)。
この手法によれば、強度パターンに対して窓関数を掛けてフーリエ変換を行う窓フーリエ変換法を用いてパターンを解析することによって、位相波面の形状を計算することができる。
【先行技術文献】
【非特許文献】
【0004】
【非特許文献1】“Windowed Fourier transform method for demodulation of carrier fringes,¨ Opt.Eng.43(7) 1472−1473(July 2004)
【発明の概要】
【発明が解決しようとする課題】
【0005】
上記した非特許文献1における窓フーリエ変換方法は、旧来行なわれていた通常のフーリエ変換法に比べて、ノイズに対する画像の耐性を向上することができるが、つぎのような課題を有している。
窓フーリエ変換法では、基本的に位相波面の形状変化が比較的滑らかな被検知体に対して有効であるが、使用する窓関数の大きさ(目安としては半値全幅が良く用いられる)によっては、取得した位相波面の画像が乱れる場合がある。
特に、窓関数の半値全幅を小さくとると一定のパターンが波面回復画像に重畳し、正確な位相波面像が取得できなくなる場合がある。
また、一方で窓関数の半値全幅を大きくとると画像は乱れにくくなるが、全体的な解像度が犠牲となる。
そのため、被検知体の形状によっては正確な位相波面の微細な形状が取得できなくなるという問題を有している。
これらは、窓フーリエ変換法の原理的な問題であり、そのため仮に元画像のノイズが無視できたとしても、解像度には限界が生じるという根本的な問題である。
【0006】
本発明は、上記課題に鑑み、窓フーリエ変換法を用いた解析において、解像度の更なる向上を図ることが可能となる位相情報の解析方法等の提供を目的とする。
【課題を解決するための手段】
【0007】
本発明は、光あるいはX線を含む波長の波の干渉によるモアレ像の周期的パターンを解析し、位相波面を含む位相情報を取得する位相情報の解析方法であって、
前記モアレ像の周期的なパターンの少なくとも一部を、窓関数によって窓フーリエ変換する工程と、
前記窓フーリエ変換されたモアレ像における、位相情報を担うスペクトルの情報と、前記位相情報を担うスペクトルの情報に重畳している位相情報を担っていないスペクトルの情報とを解析的に計算する工程と、
前記位相情報を担っていないスペクトルの情報を、前記位相情報を担うスペクトルの情報から分離し、前記位相波面を含む位相情報を取得する工程と、
を有することを特徴とする。
【発明の効果】
【0008】
本発明によれば、窓フーリエ変換法を用いた解析において、解像度の更なる向上を図ることが可能となる位相情報の解析方法等を実現することができる。
【図面の簡単な説明】
【0009】
【図1】本発明の実施形態を説明するモアレ像から波面の変化を計算する過程を表したフローチャートである。
【図2】本発明の実施形態の説明に用いたタルボ型干渉計を示した図である。
【図3】窓フーリエ変換によるモアレパターンのスペクトについて説明する概念図である。
【図4】本発明の実施例1で使用した非検知物の構造を表した図である。
【図5】本発明の実施例で用いた位相格子の構造を表した図である。図5(a)は実施例1で用いた縞格子、図5(b)は実施例2で用いた市松格子である。
【図6】本発明の実施例1の説明に用いたモアレ像を表した図である。
【図7】本発明の実施例1と従来例を比較したモアレ像から波面回復を行った図である。図7(a)は従来例の結果を示した図、図7(b)は実施例1の結果を示した図である。
【図8】本発明の実施例2の説明に用いたモアレ像を表した図である。
【図9】本発明の実施例2と従来例を比較したY軸方向の位相波面微分像を示した図である。図9(a)は従来例を示した図、図9(b)は実施例2を示した図である。
【図10】本発明の実施例2と従来例を比較したX軸方向の位相波面微分像を示した図である。図10(a)は従来例を示した図、図10(b)は実施例2を示した図である。
【発明を実施するための形態】
【0010】
本発明に係る位相情報の解析方法は、窓フーリエ変換法によりモアレ像の周期的パターンを解析するに際し、位相情報を担うスペクトルの情報から、位相情報を担っていないスペクトルの情報を解析的に分離する。
ここで、解析的とは、2点以上のデータから0次成分と1次成分のスペクトルデータを方程式に基づく計算によって算出することをいう。
すなわち、本発明の解析手法では、所定の窓関数を用いているため、フーリエ変換後のスペクトルの形状を予測することができる。そのため、0次成分と1次成分のスペクトルデータを分離する場合においても、各スペクトルデータの形状を方程式に基づいて算出することができる。
例えば、窓関数にガウシアンを使った場合、モアレ像はガウシアンによる変換によって、0次スペクトルと1次のスペクトルは近似的にガウシアンが重なった形になる。
そこで、それぞれのスペクトルがガウシアンであると仮定することにより、0次スペクトルと1次スペクトルについて解析的に分離を行うことが可能となる。
そして、分離した1次スペクトルから波面形状を計算することにより、解析的にピーク分離を行わなかった場合と比較して、より精細な波面のデータを取得することが可能となる。
【0011】
このように、本実施形態の構成によれば、解像度の更なる向上を図ることが可能となる。一方、従来の窓フーリエ変換法では、使用する窓関数の大きさによっては、取得した位相波面の画像が乱れる場合がある。以下に、これらについて更に説明する。
まず、その前提として、窓フーリエ変換法の概略について説明すると、窓フーリエ変換法では窓フーリエ変換によって切り取られた箇所のフーリエ変換は、バックグラウンドの0次スペクトルとモアレパターンによる1次スペクトルに分割される。
この1次スペクトルから窓関数によって切り取られた範囲における位相波面情報が取得でき、これら窓関数の位置をずらしていくことで、各位置の位相波面情報を繋げることができる。これによって取得した画面における位相波面形状を形成することができる。
このような、窓フーリエ変換法において、解像度を上昇させる方法の一つとして、窓関数の切り取り半径を小さく取る方法がある。
しかし、窓関数を小さく取ると、前記窓フーリエ変換による波数空間内での0次スペクトルと1次スペクトルが互いに重なり、位相波面情報に0次スペクトルの影響が出る。そのため回復された波面形状に歪みが出る。
【0012】
図3は窓関数によって窓フーリエ変換されたモアレパターンを模式的に示した図である。
図3において、30は0次スペクトル、31は1次スペクトルを示す。
図3(a)では窓関数の大きさを大きめに取った場合を示し、図3(b)は窓関数の大きさを小さめに取った場合である。
図3(a)では中央にある0次スペクトルとその両脇にある1次スペクトルはほぼ独立したスペクトルになっているため、1次スペクトルの情報をそのスペクトルの値とすればよい。また、このようにスペクトル情報を取得すれば、回復した位相波面像が乱れにくい。
一方、図3(b)では0次スペクトルと1次スペクトルがお互いに裾を広げた形で干渉しあっている。
そのため、0次と1次のスペクトルデータが重なり合い1次スペクトルの情報を単独で取得することが難しくなっている。このため単純に1次スペクトルの値を切り出すだけでは正確な位相波面形状を取得できなくなり、0次スペクトルと1次スペクトルが重畳したデータとなって取得される。
これに対して、本実施形態の構成によれば、使用する窓関数の大きさにかかわらず、窓フーリエ成分から0次スペクトルと1次スペクトルを解析的に分離することによって、0次スペクトルからの影響を排除して解像度の更なる向上を図ることが可能となる。
また、本実施形態の構成として、このような位相情報の解析方法を、コンピュータに実行させる位相情報の解析プログラムを構成することができる。
また、この位相情報の解析プログラムを記憶したコンピュータが読み取り可能の記憶媒体を構成することができる。
【0013】
つぎに、本実施形態の位相情報の解析方法について主として位相波面の情報を計算する箇所について説明する。
非特許文献1では、干渉縞を計算するために窓フーリエ変換法を用いたCarrier Fringeと言う手法を紹介している。
本実施形態ではこのような従来例のモアレ像の周期的なパターンの一部を、窓関数によって切り出してフーリエ変換し、そのスペクトルのデータから逐次位相を決定して行く、窓フーリエ変換法における位相波面の情報を計算する箇所を改善したものである。
図1に、本実施形態における従来例の窓フーリエ変換法による計算手法の手順を改善したフローチャートを示す。
図1に示すように、最初に、ステップ11において、モアレ画像(干渉縞)を取得する。
次に、ステップ12において、取得したモアレ画像に対し窓フーリエ変換を実施する。
窓フーリエ変換に使用する窓関数は様々な種類が使用可能である。
次に、ステップ13において、この窓フーリエ変換から特に1次スペクトル、すなわちモアレの周波数に一致したスペクトルのデータを抽出する。
その際、従来例ではこの取得された1次スペクトルに対応する箇所のデータの値を、そのまま使用していたのに対して、本実施形態では0次スペクトルと1次スペクトルを解析的に分離し、1次スペクトルから0次スペクトルの影響を排除する。
そのため、1次スペクトルのデータに0次スペクトルのデータが重畳していることを前提としてその差分を計算する。
この差分を計算のために、0次スペクトルに対応するスペクトルのデータも取得し、1次スペクトルに重畳した0次スペクトルの情報を解析的に求める手順を増やした。
その際に、高速化のために0次スペクトル、1次スペクトル共に形状がガウシアンで近似できるとし、フィッティングによって2つのスペクトルを分離する手順を用いた。
次に、ステップ14において、この抽出したデータから位相角を求めることによって、位相波面の変化量を求める。
上記のように求められた位相角は−πからπまでの間にラップ(畳み込み)されたデータになるため、次に、ステップ15において、その不連続点を解析して補正を行うアンラップ(位相接続)を行う。
以上のようにして、窓フーリエ成分から0次スペクトルと1次スペクトルを解析的に分離することによって、0次スペクトルからの影響を排除して得られた像が、波面の変化、もしくはその微分を表す情報となる。
微分の情報の場合はさらに積分することによって、位相波面の変化を求めることができる。
【0014】
つぎに、図2を用いて本実施形態のX線位相イメージング装置の構成例について説明する。
以下で説明する本実施例では、干渉システムとしてX線位相イメージング装置、特にその中でもタルボ干渉計を用いた構成例について説明する。
X線位相イメージング装置は近年医療応用などで注目を集めている技術である。医療応用においては被検体が人体であり、その微細な構造を精度良く取得する技術が欠かせないためである。
その中でも、タルボ干渉計は医療用X線位相イメージングの候補として現在盛んに研究が行なわれている。
但し、本発明は、タルボ干渉計やX線位相イメージング装置に限定されるものではなく、モアレや周期パターンを利用した計測技術一般に応用可能である。
タルボ干渉計を用いたX線位相イメージング装置については、A.Momose,et al.,Jpn.J.Appl.Phys.42,L866(2003)に詳細に記載されている。
【0015】
図2に、タルボ干渉計を用いたX線位相イメージング装置の構成例を示す。
図2において、210はX線源、220は被検知物、230は位相格子、240は吸収格子、250は検出装置、260は演算装置、261はCPUである。
上記X線位相イメージング装置を用いた場合において、X線を発生してから被検知物を透過させ元波面を得るまでの流れを説明する。
位相格子230はX線源からのX線が被検知物を透過する際、該X線の位相もしくは強度を変調する手段を構成する。
検出装置250は前記位相格子を透過したX線によるタルボ効果によって生じる強度情報を検出する手段を構成する。
演算手段260は、X線検出手段で取得された強度情報から、前記位相格子に入射したX線の位相情報を取得する手段を構成しており、上記した本発明の位相情報の解析方法をコンピュータに実行させるコンピュータシステムを備えている。
【0016】
これらの構成による作用について説明すると、まず、放射線発生部であるX線源210から発生したX線が被検知物220を透過する。
X線は、該検知物220を透過する際に、被検知物220の形状等に応じた波面の変化及び吸収を生じる。
被検知物220を透過したX線は、位相格子230を透過するとモアレを形成する。
X線はモアレ像が形成される位置に配された吸収格子240を透過することによって撮像装置の解像度に合うように強調される。
吸収格子240を透過したX線のモアレ縞の強度情報は、検出部である検出装置250によって検出される。
検出装置250は、放射線の干渉縞の強度情報を検出することのできる素子のことであり、例えばCCD(Charge Coupled Device)等の撮像装置が挙げられる。
検出装置250で検出された干渉縞の強度情報は、前述の解析方法における各工程の演算を行う演算装置260によって解析され、位相微分情報、即ち波面を特定の軸方向に微分した像となる。
演算装置260はCPU(Central Processing Unit)261を有している。
【実施例】
【0017】
以下に、本実施例について説明する。
[実施例1]
実施例では、コンピュータシミュレーションによる計算例について説明する。シミュレーションで使用したパラメータは以下のとおりである。
まず、X線源210より照射されるX線は17.7keVすなわち波長0.7ÅのコヒーレントなX線、即ち位相波面の揃ったX線が入射すると仮定した。
この入射したX線は被検知物220によって位相波面が変化するが、本実施例で用いた非検知物は図4のように200μm直径のリン酸カルシウム球41が4つ重なっているものを仮定した。
また、ここでは前述の位相格子として4μm縞型π格子(縞格子)を用いた。
ここで4μm縞型π格子とは、図5(a)に示す通り入射したX線の位相をπ変化させる箇所501と変化させない箇所502が1:1の縞模様で構成されており、一組の縞模様の幅が4μm周期であることを示す。
このとき、検出装置250で検出されるモアレ像は例えば図6のようになる。
【0018】
このモアレ像から波面回復を行うと図7のようになる。
図7では比較のため、図7(a)に従来例の結果を示し、図7(b)に実施例の結果を示している。
なお、従来例では非特許文献1に基づいた結果を示している。
また、窓関数にはガウシアンを用いた。窓関数の半値全幅の大きさは画像上で2ピクセルとした。
【0019】
本実施例と従来例の違いは、図1で示した波面情報を計算する手順におけるステップ13である。
従来例では参考データとして1次スペクトルに対応する箇所のデータを取得する際に、単純にデータの値をそのまま使用されるが、0次スペクトルと1次スペクトルを解析的に分離する手順が加えられている。これについては、実施形態で説明したとおりであり、重複となるので説明は省略する。
【0020】
図7は従来例と実施例で回復した位相波面の微分像にどのような違いが出るかを示した図である。
従来例では画面に横縞の模様が入っている。これは窓フーリエ変換した際の1次スペクトル像に0次スペクトルの像が重畳している為に発生した偽の像である。それに対して、本実施例では0次スペクトルが分離されているためにこのような縞模様は見られない。
このことは本発明の有効性を示している。この本発明の効果は被検知体の微細構造を再現するために窓関数を小さく設定するほど有効である。
【0021】
[実施例2]
実施例2では、位相格子として縞格子を用いた実施例1と異なり、4μmの市松π格子(市松格子)を用いた。
ここで、4μm市松π格子とは位相がπ変化する部分511としない部分512が市松格子状に交互にあらわれる形状で図5(b)に示される形状である。
窓関数の半値全幅の大きさは実施例1の場合と同様に画像上に対して2ピクセルを用いた。
このとき検出装置250で検出されるモアレ像は例えば図8で示されるようなものになり二次元の構造を持つ。
【0022】
図9及び図10は、従来例と実施例とを比較して示した回復した位相波面の微分像である。
本実施例でも、実施例1と同様に波面を計算する手順におけるステップ13の中で、0次スペクトルと1次スペクトルを解析的に分離する手順を加えた。
そのため、このような分離を行っていない従来例では縞模様が重畳している。
これに対して、本実施例の場合は余計な縞模様が重畳しない鮮明な画像がX、Y両方の位相微分像において取得できた。
すなわち、本発明はモアレの構造が一次元か二次元に拠らず有効であることが分かった。
モアレの形状変化から波面の変化ないしは位相の情報を解析する際に利用可能である。
本発明は、以上で説明した実施例1、2で使用したX線やタルボ装置などの装置に限定されるものではなく、可視光等のX線より長波長領域の電磁波を用いた一般のモアレ像解析にも利用可能である。すなわち、光あるいはX線を含む波長の波の干渉によるモアレ像の解析に利用できる。
また、以上で説明した実施例では窓フーリエ変換における窓関数をガウシアンにして解析を行ったが、これらはあくまで例示であり本発明における窓関数の形状はこれらに制限されることを意味しない。これら窓関数は任意の形状のものとそれに対応した解析手法を用いることができる。
【符号の説明】
【0023】
210:X線源
220:被検知物
230:位相格子
240:吸収格子
250:検出装置、
260:演算装置
261:CPU

【特許請求の範囲】
【請求項1】
光あるいはX線を含む波長の波の干渉によるモアレ像の周期的パターンを解析し、位相波面を含む位相情報を取得する位相情報の解析方法であって、
前記モアレ像の周期的なパターンの少なくとも一部を、窓関数によって窓フーリエ変換する工程と、
前記窓フーリエ変換されたモアレ像における、位相情報を担うスペクトルの情報と、前記位相情報を担うスペクトルの情報に重畳している位相情報を担っていないスペクトルの情報とを解析的に計算する工程と、
前記位相情報を担っていないスペクトルの情報を、前記位相情報を担うスペクトルの情報から分離し、前記位相波面を含む位相情報を取得する工程と、
を有することを特徴とする位相情報の解析方法。
【請求項2】
前記窓関数としてガウシアンを用い、前記位相情報を担うスペクトルおよび前記位相情報を担っていないスペクトルのそれぞれが、ガウシアンの形状を有していると仮定してこれらのスペクトルの情報を解析的に計算することを特徴とする請求項1に記載の位相情報の解析方法。
【請求項3】
請求項1または請求項2に記載の位相情報の解析方法を、コンピュータに実行させることを特徴とする位相情報の解析プログラム。
【請求項4】
請求項3に記載の位相情報の解析プログラムを記憶したことを特徴とするコンピュータが読み取り可能の記憶媒体。
【請求項5】
X線源と、
前記X線源からのX線が被検知物を透過する際、該X線の位相もしくは強度を変調する位相格子と、
前記位相格子を透過したX線によるタルボ効果によって生じる強度情報を検出する検出手段と、
前記検出手段で取得された強度情報から、前記位相格子に入射したX線の位相情報を取得する演算手段と、を有するX線位相イメージング装置であって、
前記演算手段が、請求項1または、請求項2に記載の位相情報の解析方法をコンピュータに実行させるコンピュータシステムを備えていることを特徴とするX線位相イメージング装置。
【請求項6】
前記位相格子が、縞格子または市松格子で構成されていることを特徴とする請求項5に記載のX線位相イメージング装置。

【図1】
image rotate

【図2】
image rotate

【図3】
image rotate

【図4】
image rotate

【図5】
image rotate

【図6】
image rotate

【図7】
image rotate

【図8】
image rotate

【図9】
image rotate

【図10】
image rotate


【公開番号】特開2011−163937(P2011−163937A)
【公開日】平成23年8月25日(2011.8.25)
【国際特許分類】
【出願番号】特願2010−27214(P2010−27214)
【出願日】平成22年2月10日(2010.2.10)
【出願人】(000001007)キヤノン株式会社 (59,756)
【Fターム(参考)】