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

Abstract

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

Description

本発明は、位相情報の解析方法、該位相情報の解析プログラム、記憶媒体、Ｘ線位相イメージング装置に関する。
特に、任意の位相波面をもつ光などの入射波を干渉させて作成したモアレ（干渉縞もしくは強度パターン）等の周期的なパターンからもとの入射波の位相波面もしくは位相波面の微分を計算する技術に関するものである。

光やＸ線を含む様々な波長の波を用いて干渉を生じさせて被検体の形状計測に用いる技術が知られている。
これらの計測技術では位相波面の揃った（コヒーレントな）入射光を被検体に対し照射し反射もしくは透過させる。
この反射光もしくは透過光は被検体の形状や組成によって波面が変化する事が知られている。
そこで、この変化を何らかの方法を用いて干渉させることにより、モアレ像（これは干渉縞とも称されるが、以下ではこれをモアレと記す。）に変換し、そのパターンを解析する。これによって変化した位相波面もしくはその位相波面の微分像（位相微分像）を計算することが可能となる。
これらの技術の適用例として代表的な例としてはレンズなどの形状測定を行う波面測定技術などがある。
また、近年は医療分野でもＸ線を用いた技術として、Ｘ線位相イメージング手法が知られている。
この手法は入射Ｘ線を被写体に透過させた際に、被写体の屈折率の違いによる位相差をモアレ縞などの周期的パターンを利用して取得する技術である。
被写体内の構成物質はそれぞれ異なる屈折率を有するため、波面の変化にそれに応じた特徴が現れる。そこでその位相波面を干渉などによって検出する。

これら干渉によって得られた強度パターンから元の波面もしくは入射光の位相波面の変化を計算する技術が位相回復法である。
位相回復法には幾つかの種類があるが、その一つに窓フーリエ変換法がある（非特許文献１）。
この手法によれば、強度パターンに対して窓関数を掛けてフーリエ変換を行う窓フーリエ変換法を用いてパターンを解析することによって、位相波面の形状を計算することができる。

"Ｗｉｎｄｏｗｅｄ Ｆｏｕｒｉｅｒ ｔｒａｎｓｆｏｒｍ ｍｅｔｈｏｄ ｆｏｒ ｄｅｍｏｄｕｌａｔｉｏｎ ｏｆ ｃａｒｒｉｅｒ ｆｒｉｎｇｅｓ，¨ Ｏｐｔ．Ｅｎｇ．４３（７） １４７２−１４７３（Ｊｕｌｙ ２００４）

上記した非特許文献１における窓フーリエ変換方法は、旧来行なわれていた通常のフーリエ変換法に比べて、ノイズに対する画像の耐性を向上することができるが、つぎのような課題を有している。
窓フーリエ変換法では、基本的に位相波面の形状変化が比較的滑らかな被検知体に対して有効であるが、使用する窓関数の大きさ（目安としては半値全幅が良く用いられる）によっては、取得した位相波面の画像が乱れる場合がある。
特に、窓関数の半値全幅を小さくとると一定のパターンが波面回復画像に重畳し、正確な位相波面像が取得できなくなる場合がある。
また、一方で窓関数の半値全幅を大きくとると画像は乱れにくくなるが、全体的な解像度が犠牲となる。
そのため、被検知体の形状によっては正確な位相波面の微細な形状が取得できなくなるという問題を有している。
これらは、窓フーリエ変換法の原理的な問題であり、そのため仮に元画像のノイズが無視できたとしても、解像度には限界が生じるという根本的な問題である。

本発明は、上記課題に鑑み、窓フーリエ変換法を用いた解析において、解像度の更なる向上を図ることが可能となる位相情報の解析方法等の提供を目的とする。

本発明は、光あるいはＸ線を含む波長の波の干渉によるモアレ像の周期的パターンを解析し、位相波面を含む位相情報を取得する位相情報の解析方法であって、
前記モアレ像の周期的なパターンの少なくとも一部を、窓関数によって窓フーリエ変換する工程と、
前記窓フーリエ変換されたモアレ像における、位相情報を担うスペクトルの情報と、前記位相情報を担うスペクトルの情報に重畳している位相情報を担っていないスペクトルの情報とを解析的に計算する工程と、
前記位相情報を担っていないスペクトルの情報を、前記位相情報を担うスペクトルの情報から分離し、前記位相波面を含む位相情報を取得する工程と、
を有することを特徴とする。

本発明によれば、窓フーリエ変換法を用いた解析において、解像度の更なる向上を図ることが可能となる位相情報の解析方法等を実現することができる。

本発明の実施形態を説明するモアレ像から波面の変化を計算する過程を表したフローチャートである。 本発明の実施形態の説明に用いたタルボ型干渉計を示した図である。 窓フーリエ変換によるモアレパターンのスペクトについて説明する概念図である。 本発明の実施例１で使用した非検知物の構造を表した図である。 本発明の実施例で用いた位相格子の構造を表した図である。図５（ａ）は実施例１で用いた縞格子、図５（ｂ）は実施例２で用いた市松格子である。 本発明の実施例１の説明に用いたモアレ像を表した図である。 本発明の実施例１と従来例を比較したモアレ像から波面回復を行った図である。図７（ａ）は従来例の結果を示した図、図７（ｂ）は実施例１の結果を示した図である。 本発明の実施例２の説明に用いたモアレ像を表した図である。 本発明の実施例２と従来例を比較したＹ軸方向の位相波面微分像を示した図である。図９（ａ）は従来例を示した図、図９（ｂ）は実施例２を示した図である。 本発明の実施例２と従来例を比較したＸ軸方向の位相波面微分像を示した図である。図１０（ａ）は従来例を示した図、図１０（ｂ）は実施例２を示した図である。

本発明に係る位相情報の解析方法は、窓フーリエ変換法によりモアレ像の周期的パターンを解析するに際し、位相情報を担うスペクトルの情報から、位相情報を担っていないスペクトルの情報を解析的に分離する。
ここで、解析的とは、２点以上のデータから０次成分と１次成分のスペクトルデータを方程式に基づく計算によって算出することをいう。
すなわち、本発明の解析手法では、所定の窓関数を用いているため、フーリエ変換後のスペクトルの形状を予測することができる。そのため、０次成分と１次成分のスペクトルデータを分離する場合においても、各スペクトルデータの形状を方程式に基づいて算出することができる。
例えば、窓関数にガウシアンを使った場合、モアレ像はガウシアンによる変換によって、０次スペクトルと１次のスペクトルは近似的にガウシアンが重なった形になる。
そこで、それぞれのスペクトルがガウシアンであると仮定することにより、０次スペクトルと１次スペクトルについて解析的に分離を行うことが可能となる。
そして、分離した１次スペクトルから波面形状を計算することにより、解析的にピーク分離を行わなかった場合と比較して、より精細な波面のデータを取得することが可能となる。

このように、本実施形態の構成によれば、解像度の更なる向上を図ることが可能となる。一方、従来の窓フーリエ変換法では、使用する窓関数の大きさによっては、取得した位相波面の画像が乱れる場合がある。以下に、これらについて更に説明する。
まず、その前提として、窓フーリエ変換法の概略について説明すると、窓フーリエ変換法では窓フーリエ変換によって切り取られた箇所のフーリエ変換は、バックグラウンドの０次スペクトルとモアレパターンによる１次スペクトルに分割される。
この１次スペクトルから窓関数によって切り取られた範囲における位相波面情報が取得でき、これら窓関数の位置をずらしていくことで、各位置の位相波面情報を繋げることができる。これによって取得した画面における位相波面形状を形成することができる。
このような、窓フーリエ変換法において、解像度を上昇させる方法の一つとして、窓関数の切り取り半径を小さく取る方法がある。
しかし、窓関数を小さく取ると、前記窓フーリエ変換による波数空間内での０次スペクトルと１次スペクトルが互いに重なり、位相波面情報に０次スペクトルの影響が出る。そのため回復された波面形状に歪みが出る。

図３は窓関数によって窓フーリエ変換されたモアレパターンを模式的に示した図である。
図３において、３０は０次スペクトル、３１は１次スペクトルを示す。
図３（ａ）では窓関数の大きさを大きめに取った場合を示し、図３（ｂ）は窓関数の大きさを小さめに取った場合である。
図３（ａ）では中央にある０次スペクトルとその両脇にある１次スペクトルはほぼ独立したスペクトルになっているため、１次スペクトルの情報をそのスペクトルの値とすればよい。また、このようにスペクトル情報を取得すれば、回復した位相波面像が乱れにくい。
一方、図３（ｂ）では０次スペクトルと１次スペクトルがお互いに裾を広げた形で干渉しあっている。
そのため、０次と１次のスペクトルデータが重なり合い１次スペクトルの情報を単独で取得することが難しくなっている。このため単純に１次スペクトルの値を切り出すだけでは正確な位相波面形状を取得できなくなり、０次スペクトルと１次スペクトルが重畳したデータとなって取得される。
これに対して、本実施形態の構成によれば、使用する窓関数の大きさにかかわらず、窓フーリエ成分から０次スペクトルと１次スペクトルを解析的に分離することによって、０次スペクトルからの影響を排除して解像度の更なる向上を図ることが可能となる。
また、本実施形態の構成として、このような位相情報の解析方法を、コンピュータに実行させる位相情報の解析プログラムを構成することができる。
また、この位相情報の解析プログラムを記憶したコンピュータが読み取り可能の記憶媒体を構成することができる。

つぎに、本実施形態の位相情報の解析方法について主として位相波面の情報を計算する箇所について説明する。
非特許文献１では、干渉縞を計算するために窓フーリエ変換法を用いたＣａｒｒｉｅｒ Ｆｒｉｎｇｅと言う手法を紹介している。
本実施形態ではこのような従来例のモアレ像の周期的なパターンの一部を、窓関数によって切り出してフーリエ変換し、そのスペクトルのデータから逐次位相を決定して行く、窓フーリエ変換法における位相波面の情報を計算する箇所を改善したものである。
図１に、本実施形態における従来例の窓フーリエ変換法による計算手法の手順を改善したフローチャートを示す。
図１に示すように、最初に、ステップ１１において、モアレ画像（干渉縞）を取得する。
次に、ステップ１２において、取得したモアレ画像に対し窓フーリエ変換を実施する。
窓フーリエ変換に使用する窓関数は様々な種類が使用可能である。
次に、ステップ１３において、この窓フーリエ変換から特に１次スペクトル、すなわちモアレの周波数に一致したスペクトルのデータを抽出する。
その際、従来例ではこの取得された１次スペクトルに対応する箇所のデータの値を、そのまま使用していたのに対して、本実施形態では０次スペクトルと１次スペクトルを解析的に分離し、１次スペクトルから０次スペクトルの影響を排除する。
そのため、１次スペクトルのデータに０次スペクトルのデータが重畳していることを前提としてその差分を計算する。
この差分を計算のために、０次スペクトルに対応するスペクトルのデータも取得し、１次スペクトルに重畳した０次スペクトルの情報を解析的に求める手順を増やした。
その際に、高速化のために０次スペクトル、１次スペクトル共に形状がガウシアンで近似できるとし、フィッティングによって２つのスペクトルを分離する手順を用いた。
次に、ステップ１４において、この抽出したデータから位相角を求めることによって、位相波面の変化量を求める。
上記のように求められた位相角は−πからπまでの間にラップ（畳み込み）されたデータになるため、次に、ステップ１５において、その不連続点を解析して補正を行うアンラップ（位相接続）を行う。
以上のようにして、窓フーリエ成分から０次スペクトルと１次スペクトルを解析的に分離することによって、０次スペクトルからの影響を排除して得られた像が、波面の変化、もしくはその微分を表す情報となる。
微分の情報の場合はさらに積分することによって、位相波面の変化を求めることができる。

つぎに、図２を用いて本実施形態のＸ線位相イメージング装置の構成例について説明する。
以下で説明する本実施例では、干渉システムとしてＸ線位相イメージング装置、特にその中でもタルボ干渉計を用いた構成例について説明する。
Ｘ線位相イメージング装置は近年医療応用などで注目を集めている技術である。医療応用においては被検体が人体であり、その微細な構造を精度良く取得する技術が欠かせないためである。
その中でも、タルボ干渉計は医療用Ｘ線位相イメージングの候補として現在盛んに研究が行なわれている。
但し、本発明は、タルボ干渉計やＸ線位相イメージング装置に限定されるものではなく、モアレや周期パターンを利用した計測技術一般に応用可能である。
タルボ干渉計を用いたＸ線位相イメージング装置については、Ａ．Ｍｏｍｏｓｅ，ｅｔ ａｌ．，Ｊｐｎ．Ｊ．Ａｐｐｌ．Ｐｈｙｓ．４２，Ｌ８６６（２００３）に詳細に記載されている。

図２に、タルボ干渉計を用いたＸ線位相イメージング装置の構成例を示す。
図２において、２１０はＸ線源、２２０は被検知物、２３０は位相格子、２４０は吸収格子、２５０は検出装置、２６０は演算装置、２６１はＣＰＵである。
上記Ｘ線位相イメージング装置を用いた場合において、Ｘ線を発生してから被検知物を透過させ元波面を得るまでの流れを説明する。
位相格子２３０はＸ線源からのＸ線が被検知物を透過する際、該Ｘ線の位相もしくは強度を変調する手段を構成する。
検出装置２５０は前記位相格子を透過したＸ線によるタルボ効果によって生じる強度情報を検出する手段を構成する。
演算手段２６０は、Ｘ線検出手段で取得された強度情報から、前記位相格子に入射したＸ線の位相情報を取得する手段を構成しており、上記した本発明の位相情報の解析方法をコンピュータに実行させるコンピュータシステムを備えている。

これらの構成による作用について説明すると、まず、放射線発生部であるＸ線源２１０から発生したＸ線が被検知物２２０を透過する。
Ｘ線は、該検知物２２０を透過する際に、被検知物２２０の形状等に応じた波面の変化及び吸収を生じる。
被検知物２２０を透過したＸ線は、位相格子２３０を透過するとモアレを形成する。
Ｘ線はモアレ像が形成される位置に配された吸収格子２４０を透過することによって撮像装置の解像度に合うように強調される。
吸収格子２４０を透過したＸ線のモアレ縞の強度情報は、検出部である検出装置２５０によって検出される。
検出装置２５０は、放射線の干渉縞の強度情報を検出することのできる素子のことであり、例えばＣＣＤ（Ｃｈａｒｇｅ Ｃｏｕｐｌｅｄ Ｄｅｖｉｃｅ）等の撮像装置が挙げられる。
検出装置２５０で検出された干渉縞の強度情報は、前述の解析方法における各工程の演算を行う演算装置２６０によって解析され、位相微分情報、即ち波面を特定の軸方向に微分した像となる。
演算装置２６０はＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）２６１を有している。

以下に、本実施例について説明する。
［実施例１］
実施例では、コンピュータシミュレーションによる計算例について説明する。シミュレーションで使用したパラメータは以下のとおりである。
まず、Ｘ線源２１０より照射されるＸ線は１７．７ｋｅＶすなわち波長０．７ÅのコヒーレントなＸ線、即ち位相波面の揃ったＸ線が入射すると仮定した。
この入射したＸ線は被検知物２２０によって位相波面が変化するが、本実施例で用いた非検知物は図４のように２００μｍ直径のリン酸カルシウム球４１が４つ重なっているものを仮定した。
また、ここでは前述の位相格子として４μｍ縞型π格子（縞格子）を用いた。
ここで４μｍ縞型π格子とは、図５（ａ）に示す通り入射したＸ線の位相をπ変化させる箇所５０１と変化させない箇所５０２が１：１の縞模様で構成されており、一組の縞模様の幅が４μｍ周期であることを示す。
このとき、検出装置２５０で検出されるモアレ像は例えば図６のようになる。

このモアレ像から波面回復を行うと図７のようになる。
図７では比較のため、図７（ａ）に従来例の結果を示し、図７（ｂ）に実施例の結果を示している。
なお、従来例では非特許文献１に基づいた結果を示している。
また、窓関数にはガウシアンを用いた。窓関数の半値全幅の大きさは画像上で２ピクセルとした。

本実施例と従来例の違いは、図１で示した波面情報を計算する手順におけるステップ１３である。
従来例では参考データとして１次スペクトルに対応する箇所のデータを取得する際に、単純にデータの値をそのまま使用されるが、０次スペクトルと１次スペクトルを解析的に分離する手順が加えられている。これについては、実施形態で説明したとおりであり、重複となるので説明は省略する。

図７は従来例と実施例で回復した位相波面の微分像にどのような違いが出るかを示した図である。
従来例では画面に横縞の模様が入っている。これは窓フーリエ変換した際の１次スペクトル像に０次スペクトルの像が重畳している為に発生した偽の像である。それに対して、本実施例では０次スペクトルが分離されているためにこのような縞模様は見られない。
このことは本発明の有効性を示している。この本発明の効果は被検知体の微細構造を再現するために窓関数を小さく設定するほど有効である。

［実施例２］
実施例２では、位相格子として縞格子を用いた実施例１と異なり、４μｍの市松π格子（市松格子）を用いた。
ここで、４μｍ市松π格子とは位相がπ変化する部分５１１としない部分５１２が市松格子状に交互にあらわれる形状で図５（ｂ）に示される形状である。
窓関数の半値全幅の大きさは実施例１の場合と同様に画像上に対して２ピクセルを用いた。
このとき検出装置２５０で検出されるモアレ像は例えば図８で示されるようなものになり二次元の構造を持つ。

図９及び図１０は、従来例と実施例とを比較して示した回復した位相波面の微分像である。
本実施例でも、実施例１と同様に波面を計算する手順におけるステップ１３の中で、０次スペクトルと１次スペクトルを解析的に分離する手順を加えた。
そのため、このような分離を行っていない従来例では縞模様が重畳している。
これに対して、本実施例の場合は余計な縞模様が重畳しない鮮明な画像がＸ、Ｙ両方の位相微分像において取得できた。
すなわち、本発明はモアレの構造が一次元か二次元に拠らず有効であることが分かった。
モアレの形状変化から波面の変化ないしは位相の情報を解析する際に利用可能である。
本発明は、以上で説明した実施例１、２で使用したＸ線やタルボ装置などの装置に限定されるものではなく、可視光等のＸ線より長波長領域の電磁波を用いた一般のモアレ像解析にも利用可能である。すなわち、光あるいはＸ線を含む波長の波の干渉によるモアレ像の解析に利用できる。
また、以上で説明した実施例では窓フーリエ変換における窓関数をガウシアンにして解析を行ったが、これらはあくまで例示であり本発明における窓関数の形状はこれらに制限されることを意味しない。これら窓関数は任意の形状のものとそれに対応した解析手法を用いることができる。

２１０：Ｘ線源
２２０：被検知物
２３０：位相格子
２４０：吸収格子
２５０：検出装置、
２６０：演算装置
２６１：ＣＰＵ

Claims (6)

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

Images