撮像装置、干渉縞解析プログラム及び干渉縞解析方法
【課題】 本発明は隣接するスペクトルの重なりの影響を低減することができ、且つ、波数空間上のスペクトルをフィッティングして位相回復を行うよりも少ない計算量で位相回復を行うことができる計算方法、プログラムと、撮像装置を提供すること。
【解決手段】 撮像装置は、シアリング干渉計と、シアリング干渉計により得られる干渉パターンから被検体の情報を算出する算出手段を備える。
算出手段は、干渉パターンを窓フーリエ変換することにより得られる波数空間上の座標におけるフーリエ成分を表す式を３以上用い、それらのフーリエ成分を表す式を連立方程式として解くことで被検体の情報を算出する。
【解決手段】 撮像装置は、シアリング干渉計と、シアリング干渉計により得られる干渉パターンから被検体の情報を算出する算出手段を備える。
算出手段は、干渉パターンを窓フーリエ変換することにより得られる波数空間上の座標におけるフーリエ成分を表す式を３以上用い、それらのフーリエ成分を表す式を連立方程式として解くことで被検体の情報を算出する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は撮像装置に関し、特にシアリング干渉計を用いて被検体の情報を得る撮像装置、撮像装置に用いられるプログラム、及び解析方法に関する。
【背景技術】
【０００２】
Ｘ線を含む様々な波長の光の干渉を利用して被検体を撮像、計測する技術が知られている。
【０００３】
この技術について簡単に説明をする。
【０００４】
位相波面の揃った（コヒーレントな）光を被検体に対して照射すると、光は被検体の形状や組成によって波面が変化する。波面に変化が生じた光を何らかの方法を用いて干渉させて干渉パターン（干渉縞）を形成させ、この干渉パターンを解析して光の位相波面を回復すると、被検体の位相、散乱、吸収に関する情報を算出することができる。
【０００５】
シアリング干渉計は上述のような光の干渉を利用して光のシア波面を計測する干渉計であり、シアリング干渉計で検出される干渉パターンは被検体によって生じる光の波面変化を微分した情報を有している。
【０００６】
この技術の代表的な適用例としてはレンズなどの形状測定を行う波面測定技術がある。
【０００７】
また、その他の適用例として、Ｘ線を用いて被検体の微分位相像を得る技術もある。
【０００８】
この技術はＸ線を被写体に照射し、被検体の形状や組成によって生じたＸ線の位相差を計測する技術であり、被検体の内部構造に関する情報を有する微分位相像を算出することが可能である。
【０００９】
干渉によって得られた干渉パターンから被検体によって生じた光の波面変化を計算する方法は位相回復法と呼ばれる。
【００１０】
位相回復法には幾つかの種類があるが、その一つにフーリエ変換法と呼ばれる方法があり、なかでも非特許文献１に記載されているような、干渉パターンに対して窓関数を掛けてからフーリエ変換を行う方法を窓フーリエ変換法と呼ぶ。
【００１１】
窓フーリエ変換法は、窓関数を用いないフーリエ変換法と比較すると一般的にノイズに対して耐性があるといった特徴がある。
【先行技術文献】
【非特許文献】
【００１２】
【非特許文献１】“Ｗｉｎｄｏｗｅｄ Ｆｏｕｒｉｅｒ ｔｒａｎｓｆｏｒｍ ｍｅｔｈｏｄ ｆｏｒ ｄｅｍｏｄｕｌａｔｉｏｎ ｏｆ ｃａｒｒｉｅｒ ｆｒｉｎｇｅｓ，¨Ｏｐｔ．Ｅｎｇ．４３（７）１４７２−１４７３（Ｊｕｌｙ ２００４）
【発明の概要】
【発明が解決しようとする課題】
【００１３】
窓フーリエ変換法に用いる窓関数の大きさ（目安としては半値全幅がよく用いられる）が小さいほど、干渉パターンを局所的に周波数成分で表現することができる。そのため、空間分解能が向上する。
【００１４】
しかし、窓関数の大きさが小さいほど、波数空間において隣接するスペクトル同士の重なりが大きくなり、周波数の分解能が低下するという課題が生じる。
【００１５】
このように窓フーリエ変換は、空間分解能と周波数分解能の一方を高くすると他方が低くなるという課題を有する。
【００１６】
そこで、本発明は隣接するスペクトルの重なりの影響を低減することができる計算方法、プログラムと、撮像装置を提供することを目的とする。
【課題を解決するための手段】
【００１７】
その目的を達成するために、本発明の一側面としての撮像装置は、シアリング干渉計と、前記シアリング干渉計により得られる干渉パターンから被検体の情報を算出する算出手段と、を備えた撮像装置であって、前記算出手段は、前記干渉パターンを窓フーリエ変換することにより得られる波数空間上座標におけるフーリエ成分を表す式を３以上用い、前記フーリエ成分を表す式を連立方程式として解くことを特徴とする。
【００１８】
本発明のその他の側面については、以下で説明する実施の形態で明らかにする。
【発明の効果】
【００１９】
本発明によると、窓フーリエ変換を用いた位相回復を行う際に隣接するスペクトルの重なりの影響を低減することができる計算方法、プログラムと、撮像装置を提供することができる。
【図面の簡単な説明】
【００２０】
【図１】実施形態の撮像装置の概略図。
【図２】１次元のトールボット干渉計に用いる格子と干渉パターンの例の模式図。
【図３】２次元のトールボット干渉計に用いる格子と干渉パターンの例の模式図。
【図４】実施形態における位相回復に用いた座標を示した波数空間の模式図。
【図５】実施例と比較例のシミュレーションに用いた被検体と検出されるモアレの模式図。
【図６】実施例により得られた１２８×１２８画素の微分位相像。
【図７】比較例により得られた１２８×１２８画素の微分位相像。
【図８】窓フーリエ変換法における逐次位相変換のイメージ図。
【発明を実施するための形態】
【００２１】
周波数分解能を保ちつつ空間分解能を向上させる又は空間分解能を保ちつつ周波数分解能を向上させるためには、波数空間において隣接するスペクトルの重なりの影響を考慮して位相回復を行えばよいことが本発明の発明者が鋭意検討した結果わかった。スペクトルの重なりの影響を考慮して位相回復を行う方法として、例えば、スペクトルのフィッティングにより隣接するスペクトルを分離して位相回復を行う方法が考えられる。
【００２２】
しかし、スペクトルのフィッティングを行うには多大なデータを扱うことになる。例えば、１０００×１０００画素の画像を得たとする。それぞれの画素を中心として窓関数を適用し、フーリエ変換を実施すると、元画像の１画素毎に１０００×１０００画素の波数空間のデータが得られる。この結果、データの合計は１０００の４乗となり、スペクトルの分離と位相回復に多大な時間又はコンピュータのリソースを必要とする。
【００２３】
そこで、上述のようなスペクトルのフィッティングによりスペクトルを分離する方法よりも短時間又は低リソースで行うことができる位相回復方法を行う撮像装置について、添付の図面に基づいて以下で説明をする。なお、各図において、同一の部材については同一の参照番号を付し、重複する説明は省略する。
【００２４】
本実施形態ではシアリング干渉計としてトールボット干渉計を使用した撮像装置について説明する。ただし、本特許が適用可能なシアリング干渉計はこれに限られるわけではなく、さまざまな形態のシアリング干渉計に適用可能である。
【００２５】
図１は本実施形態における撮像装置の構成を示した図である。図１に示した撮像装置１は、トールボット干渉計２と、算出手段として計算機６１０と、を備えている。トールボット干渉計２は、光源としてのＸ線源１１０と、Ｘ線を回折する回折格子３１０、Ｘ線の一部を遮る遮蔽格子４１０、Ｘ線を検出する検出器５１０を有している。更に撮像装置１は計算機６１０による計算結果に基づいた画像を表示する画像表示装置７１０と接続されて撮像システムを構成している。
【００２６】
以下、各構成について説明をする。
【００２７】
Ｘ線源としては、連続Ｘ線を出射するＸ線源を用いても、特性Ｘ線を出射するＸ線源を用いても良く、平行Ｘ線（平行光）を出射するＸ線源を用いても、発散Ｘ線（発散光）を出射するＸ線源を用いても良い。但し、本明細書においてＸ線とはエネルギーが２以上１００ｋｅＶ以下の光を指す。
【００２８】
また、Ｘ線源１１０からのＸ線は回折格子３１０で回折されることにより、干渉パターンを形成する必要があるため、Ｘ線源１１０からのＸ線には干渉パターンを形成できる程度の空間的可干渉性が求められる。
【００２９】
Ｘ線源１１０からのＸ線は回折格子３１０により回折され、トールボット距離と呼ばれる所定の距離をおいて明部と暗部が配列した干渉パターンを形成する。但し、本明細書では、Ｘ線（光）の強度が大きい所を明部、小さい所を暗部とする。
【００３０】
本実施形態に用いられる回折格子３１０は位相型の回折格子である。回折格子として振幅型の回折格子を用いることもできるが、位相型の回折格子の方がＸ線量（光量）の損失が少ないので有利である。
【００３１】
図２（ａ）は１次元の干渉パターンを形成する位相格子３１０ａの平面上の構成の一例を表した図である。３１１は位相の基準となる部分、３１２は３１１に対してπ分Ｘ線の位相が変化する部分である。図２（ｂ）は位相格子３１０ａによって形成される干渉パターン８１０ａの明部８１１と暗部８１２を表している。
【００３２】
図３（ａ）は２次元の干渉パターンを形成する位相格子３１０ｂの構成例である。３１１は位相の基準となる部分、３１２は３１１に対してＸ線の位相をπ変化させる部分である。図３（ｂ）は位相格子３１０ｂによって形成される干渉パターン８１０ｂの明部８１１と暗部８１２を表している。
【００３３】
遮蔽格子４１０はＸ線を透過する透過部とＸ線を遮蔽する遮蔽部が配列した構造を有しており、回折格子３１０からトールボット距離はなれて配置される。これにより、遮蔽格子４１０を透過したＸ線はモアレを形成する。尚、遮蔽部は遮蔽格子を透過したＸ線がモアレを形成する程度にＸ線を遮ればよいため、Ｘ線を完全に遮らなくても良い。
【００３４】
光としてＸ線を用いたトールボット干渉計の場合、回折格子により形成される干渉パターンの周期は通常数μｍから大きくとも数十μ程度であるのに対し、一般的に用いられるＸ線の検出器の解像度は数十μｍから数百μｍ程度である。そのためこの干渉パターンを直接検出することは難しい。そこで、本実施形態のように遮蔽格子４１０を用いてモアレを形成し、このモアレを検出する方法が用いられることが多い。このようにモアレを形成する場合、遮蔽格子のピッチは干渉パターンのピッチと同じか、僅かに異なる値をとることができ、形成したいモアレのピッチによって決めることができる。尚、モアレのピッチは遮蔽格子の遮蔽部と透過部が配列した方向と、干渉パターンの明部と暗部が配列した方向のなす角度によっても変化する。モアレの周期は様々な値をとることができるが、一般的に検出器５１０の検出素子の３画素分以上が望ましいとされている。
【００３５】
図２（ｃ）は図２（ｂ）の干渉パターンに用いられる遮蔽格子４１０ａの平面上の構成を表した図である。また、図３（ｃ）は図３（ｂ）の干渉パターンに用いられる遮蔽格子４１０ｂの平面上の構成を表した図である。図２（ｃ）の遮蔽格子４１０ａ、図３（ｃ）の遮蔽格子４１０ｂ共に、透過部４１１と遮蔽部４１２が周期的に配列している。
【００３６】
図２と図３に示した回折格子と遮蔽格子の組み合わせは一例であり、他の組み合わせを用いることもできる。尚、本発明はこれらの格子の構成に依存しない。また、干渉パターンを直接検出する場合は遮蔽格子は不要である。
【００３７】
検出器５１０は、Ｘ線を検出することのできる検出素子（例えばＣＣＤ）を有しており、遮蔽格子を透過して形成されたモアレの強度分布を検出する。本実施形態の撮像装置ではモアレの強度分布を検出するが、干渉パターンの強度分布を直接検出し、解析しても良い。尚、本実施形態の場合は説明のために干渉パターンとモアレを区別して説明したが、モアレも干渉パターンの一種である。つまり、本実施形態ではモアレを検出し、検出したモアレを解析するのでモアレを用いて説明をするが、干渉パターンを直接検出する場合もモアレを検出する場合と同様に解析を行うことができる。
【００３８】
計算機６１０はトールボット干渉計２の検出器５１０の検出結果に基づいて被検体の微分位相像の情報を計算する。
【００３９】
計算機６１０によって行われる計算方法（位相回復方法）について説明するために、まず比較例としてスペクトルのフィッティングによりスペクトルを分離して微分位相像の情報を計算する位相回復方法について説明をする。
【００４０】
２次元の窓フーリエ変換は次の式によって定義される。
【００４１】
【数１】
【００４２】
ここで、ｆ（ｘ、ｙ）は元関数、ｇ（ｘ、ｙ）は窓関数、（ｘ、ｙ）は座標、（ｕ，ｖ）は窓関数の中心、（ｋｘ，ｋｙ）は波数を表す。ＷＦ［・・・］は括弧内の関数に対して窓フーリエ変換を行うことを示す演算子である。あるモアレの強度分布Ｉ（ｘ，ｙ）に対して窓フーリエ変換を行うと、それぞれの窓関数の中心位置（ｕ，ｖ）毎に波数空間（ｋｘ，ｋｙ）が得られる。
【００４３】
例えば、離散化した１０００×１０００画素の画像（本実施形態ではモアレの強度分布）を窓フーリエ変換すると、画像の夫々の画素を窓関数の中心位置とした窓毎に１０００×１０００画素の波数空間が得られる。つまり、１０００×１０００画素の波数空間が１０００×１０００枚得られる。これは、１０００×１０００画素の画像の情報が、１０００の４乗の情報に変換されることを意味する。
【００４４】
これを図にすると図８のようになる。図８（ａ）はモアレＩ（ｘ、ｙ）の模式図である。窓関数ｇ（ｕ，ｖ）で切り取られる範囲９００は任意の座標（ｕ，ｖ）を中心座標としている。この範囲９００内にフーリエ変換を実行すると、図８（ｂ）に示したような波数空間９０００が得られる。この波数空間は、０次スペクトルの９１１、１次スペクトルの９１２、９１３、９１４、９１５といったスペクトルを有しており、これらのスペクトルから被検体によるＸ線の波面の位相変化やＸ線吸収量、Ｘ線散乱の情報を計算することができる。尚、１次スペクトルはモアレの周期に由来したスペクトルである。
【００４５】
このような波数空間は通常、窓関数の中心座標（ｕ，ｖ）毎に計算される。
【００４６】
つまり、範囲９００から窓関数の中心座標を変化させた範囲９０１内にフーリエ変換を実行すると波数空間９００１が得られる。同様に範囲９０２内にフーリエ変換を実行すると波数空間９００２が、範囲９０３内にフーリエ変換を実行すると波数空間９００３が、範囲９０４内にフーリエ変換を実行すると９００４が得られる。
【００４７】
上述のように窓フーリエ変換において窓関数で切り取られる範囲９００の半径を縮小すると、隣接するスペクトルが重なり合ってしまう可能性がある。
【００４８】
そこで、隣接するスペクトルを分離する。スペクトル９１１〜９１４は窓関数の形状でスペクトルをフィッティングできると考えられるため、本比較例ではこの方法によるフィッティングを用いてスペクトルを分離する。
【００４９】
例えば、窓関数としてガウシアン窓を使用した場合、そのフーリエ変換は同じガウシアン窓となるため、波数空間上のスペクトルもガウシアン窓でフィッティングすればよい。
【００５０】
しかし、ガウシアン窓でフィッティングするにはすべての（ｕ，ｖ）の組み合わせについて窓フーリエ変換を実行し、図８（ｂ）のように波数空間の全座標におけるフーリエ成分を計算する必要がある。そのため、１０００×１０００画素の画像を用いて位相回復を行う場合、上述のようにすべての画素を窓関数の中心として窓フーリエ変換を実行し、得られた波数空間毎に「１０００×１０００画素に相当する波数空間」を用いて位相回復を行う必要がある。そのため位相回復には多大な時間がかかる。特に画像が巨大化すると指数関数的に位相回復に必要な計算時間または計算機のリソースが増す。
【００５１】
そこで、本実施形態では波数空間（ｋｘ，ｋｙ）のマップを作成せず、（式１）における（ｋｘ，ｋｙ）の組み合わせのうち数点のみのフーリエ成分を計算して位相回復を行うことで、計算量を減らす。
【００５２】
本実施形態の計算機６１０によって行われる位相回復方法について説明する。
まず、２次元の位相イメージングを仮定し、モアレが近似的に次のような形で記述できると仮定する。
【００５３】
【数２】
【００５４】
ここで、ａ（ｘ，ｙ）は被検体による吸収の量を表し、ｂ（ｘ，ｙ）はモアレの振幅を表す。また、Ｐ１（ｘ，ｙ）およびＰ２（ｘ，ｙ）は測定する位相を表す。これらは場所によって異なる値をとりうる。ω１およびω２はそれぞれｘ方向、ｙ方向のモアレの周期を表す。モアレの形状は（式２）で表される形状に限定されるわけではなく、あくまで一例であり、本実施形態はさまざまなモアレ（干渉パターン）に適用可能である。例えば、モアレの方向が画面のｘ軸方向、ｙ軸方向に沿っていないモアレは（式２）よりも複雑な式で表される。詳しく説明はしないが、このような場合は回転移動変換などを行うことで最終的に（式２）に帰結することができる。
【００５５】
また、（式２）の第３項を０にすると（式２）は１次元のモアレを示す。以下の説明は１次元のモアレにおいても適用できる。
【００５６】
（式２）を（式１）に代入する。ただし、窓関数ｇ（ｘ，ｙ）の幅は十分に小さく、上記のａ（ｘ，ｙ）、ｂ（ｘ，ｙ）、Ｐ１（ｘ，ｙ）およびＰ２（ｘ，ｙ）はその範囲内で一定値であると近似できるとする。そのため、以下の議論ではそれぞれａ，ｂ，Ｐ１，Ｐ２と略記する。また、窓関数ｇ（ｘ，ｙ）のフーリエ変換をＧ（ｘ，ｙ）と記述する。
【００５７】
（式２）を（式１）に代入すると下記（式３）が得られる。
【００５８】
【数３】
【００５９】
図４は２次元の位相イメージングの場合に、ある（ｕ，ｖ）を中心座標として窓フーリエ変換を施した時に得られる波数空間（ｋｘ，ｋｙ）のマップを示した図である。上述のように本実施形態ではこのようなマップを作成せず、波数空間における数点のみのフーリエ成分を算出するが、ここでは本実施形態の説明のためにこのようなマップを用いる。
【００６０】
ここで、（０，０）は原点であり０次スペクトルのピークの位置を示し（ω１，０）、（−ω１，０）、（０，ω２）及び（０，−ω２）は２次元モアレの１次スペクトルのピークを示す。以下図４（ａ）に示すように（０，０）、（ω１，０）と（−ω１，０）を用いて位相回復を行う方法を説明する。
【００６１】
この３つの座標におけるフーリエ成分を表す式を用いて位相を回復する。各々の座標におけるフーリエ成分の値は（式３）から以下のよう表すことができる。
【００６２】
【数４】
【００６３】
【数５】
【００６４】
【数６】
【００６５】
尚、（ω１，０）と（−ω１，０）は原点に対して点対称であり、複素共役の関係にあるので、（式６）は（式５）から求めることもできる。
【００６６】
（式４）（式５）（式６）を連立方程式として解く。
（式４）にｅｘｐ（−ｉω１ｕ）Ｇ（ω１，０）をかけて（式５）との差をとると、
【００６７】
【数７】
【００６８】
となる。
【００６９】
同様に、（式４）にｅｘｐ（＋ｉω１ｕ）Ｇ（−ω１，０）をかけて（式６）との差をとると、
【００７０】
【数８】
【００７１】
となる。ただし、ここで上式を導出するためには窓関数のフーリエ変換の特徴として以下の特性を考慮する必要がある。
Ｇ（−ωａ，ωｂ）＝Ｇ（ωａ，ωｂ）・・・（ｙ成分も同様）
Ｇ（ωａ，ωｃ）Ｇ（ωｂ，ωｃ）＝Ｇ（ωａ＋ωｂ，ωｃ）・・・（ｙ成分も同様）
このため（式４）の第４項と第５項は（式５）および（式６）の第４項と第５項と互いに打ち消しあう。
【００７２】
このように、（０，０）、（ω１，０）、（−ω１，０）の３つの座標におけるフーリエ成分を表わす式から（式７）と（式８）が導出できる。
【００７３】
（式７）と（式８）に、検出結果から算出したフーリエ成分の値（｛ＷＦ［Ｉ（ｘ，ｙ）］（ｕ，ｖ，０，０）｝、｛ＷＦ［Ｉ（ｘ，ｙ）］（ｕ，ｖ，ω１，０）｝、｛ＷＦ［Ｉ（ｘ，ｙ）］（ｕ，ｖ，−ω１，０）｝と窓関数のフーリエ変換の値（Ｇ（ω１，０）、Ｇ（０，０）、Ｇ（−ω１，０）、Ｇ（２ω１，０）、Ｇ（−２ω１，０））を代入すれば連立方程式の形でｂとＰ１が計算できる。
【００７４】
ここで、（ω１，０）と（−ω１，０）は原点に対して点対称であり、複素共役の関係にあるため、２つの座標におけるフーリエ成分の値は等しい。よって、（０，０）と（ω１，０）又は（０，０）と（−ω１，０）の２つの座標のみのフーリエ成分の値を用いることでｂとＰ１を計算することが可能である。つまり、本実施形態では、波数空間上の座標におけるフーリエ成分を表す式を３つ用い、それらのフーリエ成分を表す式から導かれる連立方程式を、波数空間上の２つの座標におけるフーリエ成分の値を用いて解く。ここで、波数空間上の２つの座標とは、第１の座標（ここでは原点）と、第１の座標と異なり且つ、第１の座標と原点に対して点対称の関係にない第２の座標（ここでは（ω１，０）又は（−ω１，０））のことを指す。
【００７５】
図４（ｂ）に示した、（０、０）、（０、ω２）と（０、−ω２）を用いて同様に計算を行うことでＰ２が求められる。また、ａに関しても、求めたｂとＰ１、Ｐ２を（式４）、（式５）（式６）のいずれかに代入することで得られる。このように、複素共役の関係を利用することで実質３つの座標におけるフーリエ成分の値を用いて（式３）で仮定したａ、ｂ、Ｐ１、Ｐ２の４つの値が求められる。
【００７６】
これら、ａ、ｂ、Ｐ１、Ｐ２から、被検体の吸収像、散乱像、微分位相像を得ることができ、さらに微分位相像を積分すると位相像を得ることができる。
【００７７】
以上より、本実施形態を用いると、比較例のようにスペクトルのフィッティングを行うよりも高速又は低リソースで窓フーリエ変換を利用した位相回復法を行うことができる。
【００７８】
以上の例では０次スペクトルのピークと１次スペクトルのピークを用いた例について説明した。しかし、ａ、ｂ、Ｐ１、Ｐ２の計算に用いる波数空間上の座標（ｋｘ，ｋｙ）の組み合わせや用いる座標の数はこれに限定されない。本実施形態ではフーリエ成分を表わす式を５つ用い、これらの式から導かれる連立方程式を解くことでａ、ｂ、Ｐ１、Ｐ２を計算したが、用いる座標の数が増えるほどより正確にａ、ｂ、Ｐ１、Ｐ２を計算することができる。例えば、フーリエ成分を表わす式を用いて連立方程式を複数たててＰ１の値を複数求め、最小二乗法を用いて最終的なＰ１としても良い。但し、検出器の画素を単位とする窓関数の区間をＲとすると、Ｒ２を超える数の式を用いても、算出できるａ、ｂ、Ｐ１、Ｐ２の正確さは向上しない。窓フーリエ変換によって得られる波数空間は元の窓関数の範囲内の画素の情報のみを含むからである。一方、用いる座標の数が増えるほど計算量は多くなる。そこで、計算に用いる式の数は、５以上Ｒ２以下が好ましい。尚、１次元のモアレの位相回復を行う場合や、１次元の微分位相像を得たい場合など、Ｐ２を求める必要がない場合はフーリエ成分を表わす式を３つ以上用いればよい。この場合も上記のように複素共役の関係を用いれば、２つの座標におけるフーリエ成分の値からａ，ｂ，Ｐ１の値を計算することができる。本明細書において１次元の微分位相像とは、位相像を１方向に微分して得られる像である。Ｐ２を求める必要がないときは、Ｒを超える数の座標における式を用いても、算出できるａ、ｂ、Ｐ１、の正確さは向上しないため、座標におけるフーリエ成分を表わす式を３以上Ｒ以下用いることが好ましい。尚、ガウシアン窓を用いる場合は、ガウシアン窓の分散をσとしたとき、９９％の情報が存在する範囲である±３σに含まれる画素を窓関数の区間とする。
【００７９】
また、モアレによっては０次や１次のスペクトルのみならず、さらに高次のスペクトルも存在する。この高次スペクトルのピークを用いた場合でも同様に連立方程式を立てて計算することが可能である。例えば図８（ｂ）で示されたスペクトル９１６，９１７，９１８，９１９といった２次スペクトルなどを利用する方法も考えられる。
【００８０】
また、用いる座標はスペクトルのピークでなくても良い。但し、フーリエ成分の絶対値が大きい座標を用いる方がノイズの影響を受けにくいので好ましい。
【００８１】
また、用いる座標はＸ軸またはＹ軸上にある方が、ＸまたはＹ軸上にない座標を用いる場合と比べて計算が簡潔にできるため好ましい。
【００８２】
本実施形態では、計算機６１０による計算を簡潔にするために複素共役の関係を用いたが、複素共役の関係を用いなくても位相を回復することはできる。この場合、連立方程式に代入するフーリエ成分の値も３つ以上必要になる。
【００８３】
尚、窓フーリエ変換は畳み込みの定理により、以下のように表すことも可能である。窓関数が原点で対称な奇関数、すなわちｇ（ｘ）＝ｇ（−ｘ）であると仮定すれば、
【００８４】
【数９】
【００８５】
となる。ここでＦ［・・・］は通常のフーリエ変換、Ｆ−１［・・・］は逆フーリエ変換を表す。（式９）から、元関数のフーリエ変換Ｆ［ｆ（ｘ，ｙ）］に、波数空間上の窓関数Ｆ［ｇ（ｕ−ｘ，ｖ−ｙ）ｅｘｐ［ｉｋｘ（ｕ−ｘ）＋ｉｋｙ（ｖ−ｙ）］］を乗じ、その逆フーリエ変換を求めることが、元関数に窓関数を乗じてからフーリエ変換を実行することと同じであることがわかる。本実施形態では（式１）を用いて位相回復を行ったが、（式９）を用いて位相回復を行っても良い。
【００８６】
以上、計算機６１０による位相回復方法について説明をした。上述のような計算を計算機６１０で行うには、上述のような計算を実行させるプログラムを計算機６１０に組み込めばよい。
【００８７】
（実施例）
実施形態に記載した撮像装置を用いて、位相回復を行ったシミュレーション結果を実施例として示す。
【００８８】
シミュレーションは、回折格子として図３（ａ）に示した位相格子３１０ｂ、遮蔽格子として図３（ｃ）に示した遮蔽格子４１０ｂ、検出器として１２８×１２８画素の検出器を備えた撮像装置１を用いた。また、被検体は図５（ａ）に示すような球状の被検体１００１を用い、この被検体１００１が検出器の検出範囲の中央に設置されているとしてシミュレーションを行った。
【００８９】
図５（ｂ）は図５（ａ）の被検体と１２８×１２８画素の検出器を用いた時に検出器によって検出されるモアレである。この検出結果を用いて、上述した位相回復方法により得られた微分位相像を図６の（ａ）と（ｂ）に示す。図６（ａ）がＸ方向の微分位相像で、図６（ｂ）がＹ方向の微分位相像である。
【００９０】
同様のシミュレーションを２５６×２５６画素と５１２×５１２画素の検出結果を用いて行い、実際に計算にかかった時間について、表１にまとめた。計算時間はそれぞれ、０．５秒、０．７秒、１．５秒であった。
【００９１】
（比較例）
実施例と同様に、上述した比較例と同じ方法を用いて位相回復を行ったシミュレーション結果を示す。本比較例の撮像装置は計算機で行われる位相回復方法のみが実施例と異なり、そのほかの構成は実施例と同じである。
【００９２】
用いた計算機のリソース上、一度に窓フーリエ空間全体を計算することは不可能であるため、窓フーリエ変換をすべての（ｕ，ｖ）の組み合わせについて計算を行って逐次（ｕ，ｖ）おける微分位相を決定する方法で実行した。実施例同様、図５（ｂ）に示した検出結果を用いて得られた微分位相像を図７（ａ）と（ｂ）に示す。図７（ａ）がＸ方向の微分位相像で、図７（ｂ）がＹ方向の微分位相像である。
【００９３】
図６と図７を比較すると、同じような微分位相像が得られていることがわかる。
【００９４】
また、実施例同様に同様のシミュレーションを２５６×２５６画素、５１２×５１２画素の検出結果を用いて行い、実際に計算にかかった時間について、表１にまとめた。計算時間はそれぞれ、１０５秒、１０４４秒、２０９３８秒と画素数に従って指数的に長い時間を要することがわかる。
また、実際に計算にかかった時間を実施例と比較すると、実施例の方が計算時間が短縮されていることが分かる。
【００９５】
本発明の骨子は、窓フーリエ変換により得られるフーリエ成分の値を表わす式を用いて位相を回復することによって、波数空間の全座標におけるフーリエ成分を計算せずに位相回復を行うことである。これにより、窓フーリエ変換による位相回復法を短時間又は低リソースで実行することができる。したがって、その理念に基づく手法であるのであれば本発明の実施の形態は上記のみに限定されない。
【００９６】
【表１】
【符号の説明】
【００９７】
１ シアリング干渉計
２ 撮像装置
６１０ 計算機
【技術分野】
【０００１】
本発明は撮像装置に関し、特にシアリング干渉計を用いて被検体の情報を得る撮像装置、撮像装置に用いられるプログラム、及び解析方法に関する。
【背景技術】
【０００２】
Ｘ線を含む様々な波長の光の干渉を利用して被検体を撮像、計測する技術が知られている。
【０００３】
この技術について簡単に説明をする。
【０００４】
位相波面の揃った（コヒーレントな）光を被検体に対して照射すると、光は被検体の形状や組成によって波面が変化する。波面に変化が生じた光を何らかの方法を用いて干渉させて干渉パターン（干渉縞）を形成させ、この干渉パターンを解析して光の位相波面を回復すると、被検体の位相、散乱、吸収に関する情報を算出することができる。
【０００５】
シアリング干渉計は上述のような光の干渉を利用して光のシア波面を計測する干渉計であり、シアリング干渉計で検出される干渉パターンは被検体によって生じる光の波面変化を微分した情報を有している。
【０００６】
この技術の代表的な適用例としてはレンズなどの形状測定を行う波面測定技術がある。
【０００７】
また、その他の適用例として、Ｘ線を用いて被検体の微分位相像を得る技術もある。
【０００８】
この技術はＸ線を被写体に照射し、被検体の形状や組成によって生じたＸ線の位相差を計測する技術であり、被検体の内部構造に関する情報を有する微分位相像を算出することが可能である。
【０００９】
干渉によって得られた干渉パターンから被検体によって生じた光の波面変化を計算する方法は位相回復法と呼ばれる。
【００１０】
位相回復法には幾つかの種類があるが、その一つにフーリエ変換法と呼ばれる方法があり、なかでも非特許文献１に記載されているような、干渉パターンに対して窓関数を掛けてからフーリエ変換を行う方法を窓フーリエ変換法と呼ぶ。
【００１１】
窓フーリエ変換法は、窓関数を用いないフーリエ変換法と比較すると一般的にノイズに対して耐性があるといった特徴がある。
【先行技術文献】
【非特許文献】
【００１２】
【非特許文献１】“Ｗｉｎｄｏｗｅｄ Ｆｏｕｒｉｅｒ ｔｒａｎｓｆｏｒｍ ｍｅｔｈｏｄ ｆｏｒ ｄｅｍｏｄｕｌａｔｉｏｎ ｏｆ ｃａｒｒｉｅｒ ｆｒｉｎｇｅｓ，¨Ｏｐｔ．Ｅｎｇ．４３（７）１４７２−１４７３（Ｊｕｌｙ ２００４）
【発明の概要】
【発明が解決しようとする課題】
【００１３】
窓フーリエ変換法に用いる窓関数の大きさ（目安としては半値全幅がよく用いられる）が小さいほど、干渉パターンを局所的に周波数成分で表現することができる。そのため、空間分解能が向上する。
【００１４】
しかし、窓関数の大きさが小さいほど、波数空間において隣接するスペクトル同士の重なりが大きくなり、周波数の分解能が低下するという課題が生じる。
【００１５】
このように窓フーリエ変換は、空間分解能と周波数分解能の一方を高くすると他方が低くなるという課題を有する。
【００１６】
そこで、本発明は隣接するスペクトルの重なりの影響を低減することができる計算方法、プログラムと、撮像装置を提供することを目的とする。
【課題を解決するための手段】
【００１７】
その目的を達成するために、本発明の一側面としての撮像装置は、シアリング干渉計と、前記シアリング干渉計により得られる干渉パターンから被検体の情報を算出する算出手段と、を備えた撮像装置であって、前記算出手段は、前記干渉パターンを窓フーリエ変換することにより得られる波数空間上座標におけるフーリエ成分を表す式を３以上用い、前記フーリエ成分を表す式を連立方程式として解くことを特徴とする。
【００１８】
本発明のその他の側面については、以下で説明する実施の形態で明らかにする。
【発明の効果】
【００１９】
本発明によると、窓フーリエ変換を用いた位相回復を行う際に隣接するスペクトルの重なりの影響を低減することができる計算方法、プログラムと、撮像装置を提供することができる。
【図面の簡単な説明】
【００２０】
【図１】実施形態の撮像装置の概略図。
【図２】１次元のトールボット干渉計に用いる格子と干渉パターンの例の模式図。
【図３】２次元のトールボット干渉計に用いる格子と干渉パターンの例の模式図。
【図４】実施形態における位相回復に用いた座標を示した波数空間の模式図。
【図５】実施例と比較例のシミュレーションに用いた被検体と検出されるモアレの模式図。
【図６】実施例により得られた１２８×１２８画素の微分位相像。
【図７】比較例により得られた１２８×１２８画素の微分位相像。
【図８】窓フーリエ変換法における逐次位相変換のイメージ図。
【発明を実施するための形態】
【００２１】
周波数分解能を保ちつつ空間分解能を向上させる又は空間分解能を保ちつつ周波数分解能を向上させるためには、波数空間において隣接するスペクトルの重なりの影響を考慮して位相回復を行えばよいことが本発明の発明者が鋭意検討した結果わかった。スペクトルの重なりの影響を考慮して位相回復を行う方法として、例えば、スペクトルのフィッティングにより隣接するスペクトルを分離して位相回復を行う方法が考えられる。
【００２２】
しかし、スペクトルのフィッティングを行うには多大なデータを扱うことになる。例えば、１０００×１０００画素の画像を得たとする。それぞれの画素を中心として窓関数を適用し、フーリエ変換を実施すると、元画像の１画素毎に１０００×１０００画素の波数空間のデータが得られる。この結果、データの合計は１０００の４乗となり、スペクトルの分離と位相回復に多大な時間又はコンピュータのリソースを必要とする。
【００２３】
そこで、上述のようなスペクトルのフィッティングによりスペクトルを分離する方法よりも短時間又は低リソースで行うことができる位相回復方法を行う撮像装置について、添付の図面に基づいて以下で説明をする。なお、各図において、同一の部材については同一の参照番号を付し、重複する説明は省略する。
【００２４】
本実施形態ではシアリング干渉計としてトールボット干渉計を使用した撮像装置について説明する。ただし、本特許が適用可能なシアリング干渉計はこれに限られるわけではなく、さまざまな形態のシアリング干渉計に適用可能である。
【００２５】
図１は本実施形態における撮像装置の構成を示した図である。図１に示した撮像装置１は、トールボット干渉計２と、算出手段として計算機６１０と、を備えている。トールボット干渉計２は、光源としてのＸ線源１１０と、Ｘ線を回折する回折格子３１０、Ｘ線の一部を遮る遮蔽格子４１０、Ｘ線を検出する検出器５１０を有している。更に撮像装置１は計算機６１０による計算結果に基づいた画像を表示する画像表示装置７１０と接続されて撮像システムを構成している。
【００２６】
以下、各構成について説明をする。
【００２７】
Ｘ線源としては、連続Ｘ線を出射するＸ線源を用いても、特性Ｘ線を出射するＸ線源を用いても良く、平行Ｘ線（平行光）を出射するＸ線源を用いても、発散Ｘ線（発散光）を出射するＸ線源を用いても良い。但し、本明細書においてＸ線とはエネルギーが２以上１００ｋｅＶ以下の光を指す。
【００２８】
また、Ｘ線源１１０からのＸ線は回折格子３１０で回折されることにより、干渉パターンを形成する必要があるため、Ｘ線源１１０からのＸ線には干渉パターンを形成できる程度の空間的可干渉性が求められる。
【００２９】
Ｘ線源１１０からのＸ線は回折格子３１０により回折され、トールボット距離と呼ばれる所定の距離をおいて明部と暗部が配列した干渉パターンを形成する。但し、本明細書では、Ｘ線（光）の強度が大きい所を明部、小さい所を暗部とする。
【００３０】
本実施形態に用いられる回折格子３１０は位相型の回折格子である。回折格子として振幅型の回折格子を用いることもできるが、位相型の回折格子の方がＸ線量（光量）の損失が少ないので有利である。
【００３１】
図２（ａ）は１次元の干渉パターンを形成する位相格子３１０ａの平面上の構成の一例を表した図である。３１１は位相の基準となる部分、３１２は３１１に対してπ分Ｘ線の位相が変化する部分である。図２（ｂ）は位相格子３１０ａによって形成される干渉パターン８１０ａの明部８１１と暗部８１２を表している。
【００３２】
図３（ａ）は２次元の干渉パターンを形成する位相格子３１０ｂの構成例である。３１１は位相の基準となる部分、３１２は３１１に対してＸ線の位相をπ変化させる部分である。図３（ｂ）は位相格子３１０ｂによって形成される干渉パターン８１０ｂの明部８１１と暗部８１２を表している。
【００３３】
遮蔽格子４１０はＸ線を透過する透過部とＸ線を遮蔽する遮蔽部が配列した構造を有しており、回折格子３１０からトールボット距離はなれて配置される。これにより、遮蔽格子４１０を透過したＸ線はモアレを形成する。尚、遮蔽部は遮蔽格子を透過したＸ線がモアレを形成する程度にＸ線を遮ればよいため、Ｘ線を完全に遮らなくても良い。
【００３４】
光としてＸ線を用いたトールボット干渉計の場合、回折格子により形成される干渉パターンの周期は通常数μｍから大きくとも数十μ程度であるのに対し、一般的に用いられるＸ線の検出器の解像度は数十μｍから数百μｍ程度である。そのためこの干渉パターンを直接検出することは難しい。そこで、本実施形態のように遮蔽格子４１０を用いてモアレを形成し、このモアレを検出する方法が用いられることが多い。このようにモアレを形成する場合、遮蔽格子のピッチは干渉パターンのピッチと同じか、僅かに異なる値をとることができ、形成したいモアレのピッチによって決めることができる。尚、モアレのピッチは遮蔽格子の遮蔽部と透過部が配列した方向と、干渉パターンの明部と暗部が配列した方向のなす角度によっても変化する。モアレの周期は様々な値をとることができるが、一般的に検出器５１０の検出素子の３画素分以上が望ましいとされている。
【００３５】
図２（ｃ）は図２（ｂ）の干渉パターンに用いられる遮蔽格子４１０ａの平面上の構成を表した図である。また、図３（ｃ）は図３（ｂ）の干渉パターンに用いられる遮蔽格子４１０ｂの平面上の構成を表した図である。図２（ｃ）の遮蔽格子４１０ａ、図３（ｃ）の遮蔽格子４１０ｂ共に、透過部４１１と遮蔽部４１２が周期的に配列している。
【００３６】
図２と図３に示した回折格子と遮蔽格子の組み合わせは一例であり、他の組み合わせを用いることもできる。尚、本発明はこれらの格子の構成に依存しない。また、干渉パターンを直接検出する場合は遮蔽格子は不要である。
【００３７】
検出器５１０は、Ｘ線を検出することのできる検出素子（例えばＣＣＤ）を有しており、遮蔽格子を透過して形成されたモアレの強度分布を検出する。本実施形態の撮像装置ではモアレの強度分布を検出するが、干渉パターンの強度分布を直接検出し、解析しても良い。尚、本実施形態の場合は説明のために干渉パターンとモアレを区別して説明したが、モアレも干渉パターンの一種である。つまり、本実施形態ではモアレを検出し、検出したモアレを解析するのでモアレを用いて説明をするが、干渉パターンを直接検出する場合もモアレを検出する場合と同様に解析を行うことができる。
【００３８】
計算機６１０はトールボット干渉計２の検出器５１０の検出結果に基づいて被検体の微分位相像の情報を計算する。
【００３９】
計算機６１０によって行われる計算方法（位相回復方法）について説明するために、まず比較例としてスペクトルのフィッティングによりスペクトルを分離して微分位相像の情報を計算する位相回復方法について説明をする。
【００４０】
２次元の窓フーリエ変換は次の式によって定義される。
【００４１】
【数１】
【００４２】
ここで、ｆ（ｘ、ｙ）は元関数、ｇ（ｘ、ｙ）は窓関数、（ｘ、ｙ）は座標、（ｕ，ｖ）は窓関数の中心、（ｋｘ，ｋｙ）は波数を表す。ＷＦ［・・・］は括弧内の関数に対して窓フーリエ変換を行うことを示す演算子である。あるモアレの強度分布Ｉ（ｘ，ｙ）に対して窓フーリエ変換を行うと、それぞれの窓関数の中心位置（ｕ，ｖ）毎に波数空間（ｋｘ，ｋｙ）が得られる。
【００４３】
例えば、離散化した１０００×１０００画素の画像（本実施形態ではモアレの強度分布）を窓フーリエ変換すると、画像の夫々の画素を窓関数の中心位置とした窓毎に１０００×１０００画素の波数空間が得られる。つまり、１０００×１０００画素の波数空間が１０００×１０００枚得られる。これは、１０００×１０００画素の画像の情報が、１０００の４乗の情報に変換されることを意味する。
【００４４】
これを図にすると図８のようになる。図８（ａ）はモアレＩ（ｘ、ｙ）の模式図である。窓関数ｇ（ｕ，ｖ）で切り取られる範囲９００は任意の座標（ｕ，ｖ）を中心座標としている。この範囲９００内にフーリエ変換を実行すると、図８（ｂ）に示したような波数空間９０００が得られる。この波数空間は、０次スペクトルの９１１、１次スペクトルの９１２、９１３、９１４、９１５といったスペクトルを有しており、これらのスペクトルから被検体によるＸ線の波面の位相変化やＸ線吸収量、Ｘ線散乱の情報を計算することができる。尚、１次スペクトルはモアレの周期に由来したスペクトルである。
【００４５】
このような波数空間は通常、窓関数の中心座標（ｕ，ｖ）毎に計算される。
【００４６】
つまり、範囲９００から窓関数の中心座標を変化させた範囲９０１内にフーリエ変換を実行すると波数空間９００１が得られる。同様に範囲９０２内にフーリエ変換を実行すると波数空間９００２が、範囲９０３内にフーリエ変換を実行すると波数空間９００３が、範囲９０４内にフーリエ変換を実行すると９００４が得られる。
【００４７】
上述のように窓フーリエ変換において窓関数で切り取られる範囲９００の半径を縮小すると、隣接するスペクトルが重なり合ってしまう可能性がある。
【００４８】
そこで、隣接するスペクトルを分離する。スペクトル９１１〜９１４は窓関数の形状でスペクトルをフィッティングできると考えられるため、本比較例ではこの方法によるフィッティングを用いてスペクトルを分離する。
【００４９】
例えば、窓関数としてガウシアン窓を使用した場合、そのフーリエ変換は同じガウシアン窓となるため、波数空間上のスペクトルもガウシアン窓でフィッティングすればよい。
【００５０】
しかし、ガウシアン窓でフィッティングするにはすべての（ｕ，ｖ）の組み合わせについて窓フーリエ変換を実行し、図８（ｂ）のように波数空間の全座標におけるフーリエ成分を計算する必要がある。そのため、１０００×１０００画素の画像を用いて位相回復を行う場合、上述のようにすべての画素を窓関数の中心として窓フーリエ変換を実行し、得られた波数空間毎に「１０００×１０００画素に相当する波数空間」を用いて位相回復を行う必要がある。そのため位相回復には多大な時間がかかる。特に画像が巨大化すると指数関数的に位相回復に必要な計算時間または計算機のリソースが増す。
【００５１】
そこで、本実施形態では波数空間（ｋｘ，ｋｙ）のマップを作成せず、（式１）における（ｋｘ，ｋｙ）の組み合わせのうち数点のみのフーリエ成分を計算して位相回復を行うことで、計算量を減らす。
【００５２】
本実施形態の計算機６１０によって行われる位相回復方法について説明する。
まず、２次元の位相イメージングを仮定し、モアレが近似的に次のような形で記述できると仮定する。
【００５３】
【数２】
【００５４】
ここで、ａ（ｘ，ｙ）は被検体による吸収の量を表し、ｂ（ｘ，ｙ）はモアレの振幅を表す。また、Ｐ１（ｘ，ｙ）およびＰ２（ｘ，ｙ）は測定する位相を表す。これらは場所によって異なる値をとりうる。ω１およびω２はそれぞれｘ方向、ｙ方向のモアレの周期を表す。モアレの形状は（式２）で表される形状に限定されるわけではなく、あくまで一例であり、本実施形態はさまざまなモアレ（干渉パターン）に適用可能である。例えば、モアレの方向が画面のｘ軸方向、ｙ軸方向に沿っていないモアレは（式２）よりも複雑な式で表される。詳しく説明はしないが、このような場合は回転移動変換などを行うことで最終的に（式２）に帰結することができる。
【００５５】
また、（式２）の第３項を０にすると（式２）は１次元のモアレを示す。以下の説明は１次元のモアレにおいても適用できる。
【００５６】
（式２）を（式１）に代入する。ただし、窓関数ｇ（ｘ，ｙ）の幅は十分に小さく、上記のａ（ｘ，ｙ）、ｂ（ｘ，ｙ）、Ｐ１（ｘ，ｙ）およびＰ２（ｘ，ｙ）はその範囲内で一定値であると近似できるとする。そのため、以下の議論ではそれぞれａ，ｂ，Ｐ１，Ｐ２と略記する。また、窓関数ｇ（ｘ，ｙ）のフーリエ変換をＧ（ｘ，ｙ）と記述する。
【００５７】
（式２）を（式１）に代入すると下記（式３）が得られる。
【００５８】
【数３】
【００５９】
図４は２次元の位相イメージングの場合に、ある（ｕ，ｖ）を中心座標として窓フーリエ変換を施した時に得られる波数空間（ｋｘ，ｋｙ）のマップを示した図である。上述のように本実施形態ではこのようなマップを作成せず、波数空間における数点のみのフーリエ成分を算出するが、ここでは本実施形態の説明のためにこのようなマップを用いる。
【００６０】
ここで、（０，０）は原点であり０次スペクトルのピークの位置を示し（ω１，０）、（−ω１，０）、（０，ω２）及び（０，−ω２）は２次元モアレの１次スペクトルのピークを示す。以下図４（ａ）に示すように（０，０）、（ω１，０）と（−ω１，０）を用いて位相回復を行う方法を説明する。
【００６１】
この３つの座標におけるフーリエ成分を表す式を用いて位相を回復する。各々の座標におけるフーリエ成分の値は（式３）から以下のよう表すことができる。
【００６２】
【数４】
【００６３】
【数５】
【００６４】
【数６】
【００６５】
尚、（ω１，０）と（−ω１，０）は原点に対して点対称であり、複素共役の関係にあるので、（式６）は（式５）から求めることもできる。
【００６６】
（式４）（式５）（式６）を連立方程式として解く。
（式４）にｅｘｐ（−ｉω１ｕ）Ｇ（ω１，０）をかけて（式５）との差をとると、
【００６７】
【数７】
【００６８】
となる。
【００６９】
同様に、（式４）にｅｘｐ（＋ｉω１ｕ）Ｇ（−ω１，０）をかけて（式６）との差をとると、
【００７０】
【数８】
【００７１】
となる。ただし、ここで上式を導出するためには窓関数のフーリエ変換の特徴として以下の特性を考慮する必要がある。
Ｇ（−ωａ，ωｂ）＝Ｇ（ωａ，ωｂ）・・・（ｙ成分も同様）
Ｇ（ωａ，ωｃ）Ｇ（ωｂ，ωｃ）＝Ｇ（ωａ＋ωｂ，ωｃ）・・・（ｙ成分も同様）
このため（式４）の第４項と第５項は（式５）および（式６）の第４項と第５項と互いに打ち消しあう。
【００７２】
このように、（０，０）、（ω１，０）、（−ω１，０）の３つの座標におけるフーリエ成分を表わす式から（式７）と（式８）が導出できる。
【００７３】
（式７）と（式８）に、検出結果から算出したフーリエ成分の値（｛ＷＦ［Ｉ（ｘ，ｙ）］（ｕ，ｖ，０，０）｝、｛ＷＦ［Ｉ（ｘ，ｙ）］（ｕ，ｖ，ω１，０）｝、｛ＷＦ［Ｉ（ｘ，ｙ）］（ｕ，ｖ，−ω１，０）｝と窓関数のフーリエ変換の値（Ｇ（ω１，０）、Ｇ（０，０）、Ｇ（−ω１，０）、Ｇ（２ω１，０）、Ｇ（−２ω１，０））を代入すれば連立方程式の形でｂとＰ１が計算できる。
【００７４】
ここで、（ω１，０）と（−ω１，０）は原点に対して点対称であり、複素共役の関係にあるため、２つの座標におけるフーリエ成分の値は等しい。よって、（０，０）と（ω１，０）又は（０，０）と（−ω１，０）の２つの座標のみのフーリエ成分の値を用いることでｂとＰ１を計算することが可能である。つまり、本実施形態では、波数空間上の座標におけるフーリエ成分を表す式を３つ用い、それらのフーリエ成分を表す式から導かれる連立方程式を、波数空間上の２つの座標におけるフーリエ成分の値を用いて解く。ここで、波数空間上の２つの座標とは、第１の座標（ここでは原点）と、第１の座標と異なり且つ、第１の座標と原点に対して点対称の関係にない第２の座標（ここでは（ω１，０）又は（−ω１，０））のことを指す。
【００７５】
図４（ｂ）に示した、（０、０）、（０、ω２）と（０、−ω２）を用いて同様に計算を行うことでＰ２が求められる。また、ａに関しても、求めたｂとＰ１、Ｐ２を（式４）、（式５）（式６）のいずれかに代入することで得られる。このように、複素共役の関係を利用することで実質３つの座標におけるフーリエ成分の値を用いて（式３）で仮定したａ、ｂ、Ｐ１、Ｐ２の４つの値が求められる。
【００７６】
これら、ａ、ｂ、Ｐ１、Ｐ２から、被検体の吸収像、散乱像、微分位相像を得ることができ、さらに微分位相像を積分すると位相像を得ることができる。
【００７７】
以上より、本実施形態を用いると、比較例のようにスペクトルのフィッティングを行うよりも高速又は低リソースで窓フーリエ変換を利用した位相回復法を行うことができる。
【００７８】
以上の例では０次スペクトルのピークと１次スペクトルのピークを用いた例について説明した。しかし、ａ、ｂ、Ｐ１、Ｐ２の計算に用いる波数空間上の座標（ｋｘ，ｋｙ）の組み合わせや用いる座標の数はこれに限定されない。本実施形態ではフーリエ成分を表わす式を５つ用い、これらの式から導かれる連立方程式を解くことでａ、ｂ、Ｐ１、Ｐ２を計算したが、用いる座標の数が増えるほどより正確にａ、ｂ、Ｐ１、Ｐ２を計算することができる。例えば、フーリエ成分を表わす式を用いて連立方程式を複数たててＰ１の値を複数求め、最小二乗法を用いて最終的なＰ１としても良い。但し、検出器の画素を単位とする窓関数の区間をＲとすると、Ｒ２を超える数の式を用いても、算出できるａ、ｂ、Ｐ１、Ｐ２の正確さは向上しない。窓フーリエ変換によって得られる波数空間は元の窓関数の範囲内の画素の情報のみを含むからである。一方、用いる座標の数が増えるほど計算量は多くなる。そこで、計算に用いる式の数は、５以上Ｒ２以下が好ましい。尚、１次元のモアレの位相回復を行う場合や、１次元の微分位相像を得たい場合など、Ｐ２を求める必要がない場合はフーリエ成分を表わす式を３つ以上用いればよい。この場合も上記のように複素共役の関係を用いれば、２つの座標におけるフーリエ成分の値からａ，ｂ，Ｐ１の値を計算することができる。本明細書において１次元の微分位相像とは、位相像を１方向に微分して得られる像である。Ｐ２を求める必要がないときは、Ｒを超える数の座標における式を用いても、算出できるａ、ｂ、Ｐ１、の正確さは向上しないため、座標におけるフーリエ成分を表わす式を３以上Ｒ以下用いることが好ましい。尚、ガウシアン窓を用いる場合は、ガウシアン窓の分散をσとしたとき、９９％の情報が存在する範囲である±３σに含まれる画素を窓関数の区間とする。
【００７９】
また、モアレによっては０次や１次のスペクトルのみならず、さらに高次のスペクトルも存在する。この高次スペクトルのピークを用いた場合でも同様に連立方程式を立てて計算することが可能である。例えば図８（ｂ）で示されたスペクトル９１６，９１７，９１８，９１９といった２次スペクトルなどを利用する方法も考えられる。
【００８０】
また、用いる座標はスペクトルのピークでなくても良い。但し、フーリエ成分の絶対値が大きい座標を用いる方がノイズの影響を受けにくいので好ましい。
【００８１】
また、用いる座標はＸ軸またはＹ軸上にある方が、ＸまたはＹ軸上にない座標を用いる場合と比べて計算が簡潔にできるため好ましい。
【００８２】
本実施形態では、計算機６１０による計算を簡潔にするために複素共役の関係を用いたが、複素共役の関係を用いなくても位相を回復することはできる。この場合、連立方程式に代入するフーリエ成分の値も３つ以上必要になる。
【００８３】
尚、窓フーリエ変換は畳み込みの定理により、以下のように表すことも可能である。窓関数が原点で対称な奇関数、すなわちｇ（ｘ）＝ｇ（−ｘ）であると仮定すれば、
【００８４】
【数９】
【００８５】
となる。ここでＦ［・・・］は通常のフーリエ変換、Ｆ−１［・・・］は逆フーリエ変換を表す。（式９）から、元関数のフーリエ変換Ｆ［ｆ（ｘ，ｙ）］に、波数空間上の窓関数Ｆ［ｇ（ｕ−ｘ，ｖ−ｙ）ｅｘｐ［ｉｋｘ（ｕ−ｘ）＋ｉｋｙ（ｖ−ｙ）］］を乗じ、その逆フーリエ変換を求めることが、元関数に窓関数を乗じてからフーリエ変換を実行することと同じであることがわかる。本実施形態では（式１）を用いて位相回復を行ったが、（式９）を用いて位相回復を行っても良い。
【００８６】
以上、計算機６１０による位相回復方法について説明をした。上述のような計算を計算機６１０で行うには、上述のような計算を実行させるプログラムを計算機６１０に組み込めばよい。
【００８７】
（実施例）
実施形態に記載した撮像装置を用いて、位相回復を行ったシミュレーション結果を実施例として示す。
【００８８】
シミュレーションは、回折格子として図３（ａ）に示した位相格子３１０ｂ、遮蔽格子として図３（ｃ）に示した遮蔽格子４１０ｂ、検出器として１２８×１２８画素の検出器を備えた撮像装置１を用いた。また、被検体は図５（ａ）に示すような球状の被検体１００１を用い、この被検体１００１が検出器の検出範囲の中央に設置されているとしてシミュレーションを行った。
【００８９】
図５（ｂ）は図５（ａ）の被検体と１２８×１２８画素の検出器を用いた時に検出器によって検出されるモアレである。この検出結果を用いて、上述した位相回復方法により得られた微分位相像を図６の（ａ）と（ｂ）に示す。図６（ａ）がＸ方向の微分位相像で、図６（ｂ）がＹ方向の微分位相像である。
【００９０】
同様のシミュレーションを２５６×２５６画素と５１２×５１２画素の検出結果を用いて行い、実際に計算にかかった時間について、表１にまとめた。計算時間はそれぞれ、０．５秒、０．７秒、１．５秒であった。
【００９１】
（比較例）
実施例と同様に、上述した比較例と同じ方法を用いて位相回復を行ったシミュレーション結果を示す。本比較例の撮像装置は計算機で行われる位相回復方法のみが実施例と異なり、そのほかの構成は実施例と同じである。
【００９２】
用いた計算機のリソース上、一度に窓フーリエ空間全体を計算することは不可能であるため、窓フーリエ変換をすべての（ｕ，ｖ）の組み合わせについて計算を行って逐次（ｕ，ｖ）おける微分位相を決定する方法で実行した。実施例同様、図５（ｂ）に示した検出結果を用いて得られた微分位相像を図７（ａ）と（ｂ）に示す。図７（ａ）がＸ方向の微分位相像で、図７（ｂ）がＹ方向の微分位相像である。
【００９３】
図６と図７を比較すると、同じような微分位相像が得られていることがわかる。
【００９４】
また、実施例同様に同様のシミュレーションを２５６×２５６画素、５１２×５１２画素の検出結果を用いて行い、実際に計算にかかった時間について、表１にまとめた。計算時間はそれぞれ、１０５秒、１０４４秒、２０９３８秒と画素数に従って指数的に長い時間を要することがわかる。
また、実際に計算にかかった時間を実施例と比較すると、実施例の方が計算時間が短縮されていることが分かる。
【００９５】
本発明の骨子は、窓フーリエ変換により得られるフーリエ成分の値を表わす式を用いて位相を回復することによって、波数空間の全座標におけるフーリエ成分を計算せずに位相回復を行うことである。これにより、窓フーリエ変換による位相回復法を短時間又は低リソースで実行することができる。したがって、その理念に基づく手法であるのであれば本発明の実施の形態は上記のみに限定されない。
【００９６】
【表１】
【符号の説明】
【００９７】
１ シアリング干渉計
２ 撮像装置
６１０ 計算機
【特許請求の範囲】
【請求項１】
シアリング干渉計と、前記シアリング干渉計により得られる干渉パターンから被検体の情報を算出する算出手段と、を備えた撮像装置であって、
前記算出手段は、
前記干渉パターンを窓フーリエ変換することにより得られる波数空間上の座標におけるフーリエ成分を表す式を３以上用い、前記フーリエ成分を表す式を連立方程式として解くことを特徴とする撮像装置。
【請求項２】
前記被検体の情報は、
前記被検体の吸収像、散乱像、微分位相像、位相像、の少なくとも１つに関する情報であることを特徴とする請求項１に記載の撮像装置。
【請求項３】
前記算出手段は、
前記フーリエ成分を表す式を３以上Ｒ以下用い、
前記被検体の１次元の、吸収像、散乱像、微分位相像、位相像、の少なくとも１つに関する情報を算出することを特徴とする請求項１に記載の撮像装置。
但し、Ｒとは検出器の画素を単位とする窓関数の区間を示す。
【請求項４】
前記算出手段は、
前記フーリエ成分を表す式を５以上Ｒ２以下の座標における前記フーリエ成分を表わす式であり、
前記算出手段は、
前記被検体の２次元の、吸収像、散乱像、微分位相像、位相像、の少なくとも１つに関する情報を算出することを特徴とする請求項１に記載の撮像装置。
但し、Ｒとは検出器の画素を単位とする窓関数の区間を示す。
【請求項５】
前記フーリエ成分を表す式のうち少なくとも１つは、
前記干渉パターンの周期に由来する前記フーリエ成分を表す式であることを特徴とする請求項１乃至４のいずれか１項に記載の撮像装置。
【請求項６】
前記フーリエ成分を表す式の１つが、
前記波数空間の原点における前記フーリエ成分を表す式であることを特徴とする請求項１乃至５のいずれか１項に記載の撮像装置。
【請求項７】
前記連立方程式は、
前記波数空間の２以上の座標における前記フーリエ成分を算出した値を用いて解くことを特徴とする請求項１乃至６に記載の撮像装置。
【請求項８】
前記フーリエ成分を表す式のうち少なくとも２つは、
前記波数空間の原点に対して点対称の関係にある２つの座標における前記フーリエ成分を表す式であることを特徴とする請求項１乃至７のいずれか１項に記載の撮像装置。
【請求項９】
前記フーリエ成分を表す式は、
前記波数空間のＸ軸又はＹ軸上の座標における前記フーリエ成分を表す式であることを特徴とする請求項１乃至８のいずれか１項に記載の撮像装置。
【請求項１０】
シアリング干渉計により得られる干渉パターンから被検体の情報を算出するプログラムであって、
前記干渉パターンを窓フーリエ変換することにより得られる波数空間上の座標におけるフーリエ成分を表す式を３以上用い、前記フーリエ成分を表す式を連立方程式として解くことを特徴とするプログラム。
【請求項１１】
シアリング干渉計により得られる干渉パターンから被検体の情報を得る方法であって、
前記干渉パターンを窓フーリエ変換することにより得られる波数空間上の座標におけるフーリエ成分を表す式を３以上用い、前記フーリエ成分を表す式を連立方程式として解くことを特徴とする方法。
【請求項１】
シアリング干渉計と、前記シアリング干渉計により得られる干渉パターンから被検体の情報を算出する算出手段と、を備えた撮像装置であって、
前記算出手段は、
前記干渉パターンを窓フーリエ変換することにより得られる波数空間上の座標におけるフーリエ成分を表す式を３以上用い、前記フーリエ成分を表す式を連立方程式として解くことを特徴とする撮像装置。
【請求項２】
前記被検体の情報は、
前記被検体の吸収像、散乱像、微分位相像、位相像、の少なくとも１つに関する情報であることを特徴とする請求項１に記載の撮像装置。
【請求項３】
前記算出手段は、
前記フーリエ成分を表す式を３以上Ｒ以下用い、
前記被検体の１次元の、吸収像、散乱像、微分位相像、位相像、の少なくとも１つに関する情報を算出することを特徴とする請求項１に記載の撮像装置。
但し、Ｒとは検出器の画素を単位とする窓関数の区間を示す。
【請求項４】
前記算出手段は、
前記フーリエ成分を表す式を５以上Ｒ２以下の座標における前記フーリエ成分を表わす式であり、
前記算出手段は、
前記被検体の２次元の、吸収像、散乱像、微分位相像、位相像、の少なくとも１つに関する情報を算出することを特徴とする請求項１に記載の撮像装置。
但し、Ｒとは検出器の画素を単位とする窓関数の区間を示す。
【請求項５】
前記フーリエ成分を表す式のうち少なくとも１つは、
前記干渉パターンの周期に由来する前記フーリエ成分を表す式であることを特徴とする請求項１乃至４のいずれか１項に記載の撮像装置。
【請求項６】
前記フーリエ成分を表す式の１つが、
前記波数空間の原点における前記フーリエ成分を表す式であることを特徴とする請求項１乃至５のいずれか１項に記載の撮像装置。
【請求項７】
前記連立方程式は、
前記波数空間の２以上の座標における前記フーリエ成分を算出した値を用いて解くことを特徴とする請求項１乃至６に記載の撮像装置。
【請求項８】
前記フーリエ成分を表す式のうち少なくとも２つは、
前記波数空間の原点に対して点対称の関係にある２つの座標における前記フーリエ成分を表す式であることを特徴とする請求項１乃至７のいずれか１項に記載の撮像装置。
【請求項９】
前記フーリエ成分を表す式は、
前記波数空間のＸ軸又はＹ軸上の座標における前記フーリエ成分を表す式であることを特徴とする請求項１乃至８のいずれか１項に記載の撮像装置。
【請求項１０】
シアリング干渉計により得られる干渉パターンから被検体の情報を算出するプログラムであって、
前記干渉パターンを窓フーリエ変換することにより得られる波数空間上の座標におけるフーリエ成分を表す式を３以上用い、前記フーリエ成分を表す式を連立方程式として解くことを特徴とするプログラム。
【請求項１１】
シアリング干渉計により得られる干渉パターンから被検体の情報を得る方法であって、
前記干渉パターンを窓フーリエ変換することにより得られる波数空間上の座標におけるフーリエ成分を表す式を３以上用い、前記フーリエ成分を表す式を連立方程式として解くことを特徴とする方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【公開番号】特開２０１３−２８４５（Ｐ２０１３−２８４５Ａ）
【公開日】平成２５年１月７日（２０１３．１．７）
【国際特許分類】
【出願番号】特願２０１１−１３１４９８（Ｐ２０１１−１３１４９８）
【出願日】平成２３年６月１３日（２０１１．６．１３）
【出願人】（０００００１００７）キヤノン株式会社 (59,756)
【Ｆターム（参考）】
【公開日】平成２５年１月７日（２０１３．１．７）
【国際特許分類】
【出願日】平成２３年６月１３日（２０１１．６．１３）
【出願人】（０００００１００７）キヤノン株式会社 (59,756)
【Ｆターム（参考）】
[ Back to top ]