JP5885405B2 - 撮像装置、干渉縞解析プログラム及び干渉縞解析方法 - Google Patents

撮像装置、干渉縞解析プログラム及び干渉縞解析方法 Download PDF

Info

Publication number
JP5885405B2
JP5885405B2 JP2011131498A JP2011131498A JP5885405B2 JP 5885405 B2 JP5885405 B2 JP 5885405B2 JP 2011131498 A JP2011131498 A JP 2011131498A JP 2011131498 A JP2011131498 A JP 2011131498A JP 5885405 B2 JP5885405 B2 JP 5885405B2
Authority
JP
Grant status
Grant
Patent type
Prior art keywords
imaging apparatus
phase
image
interference pattern
window
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.)
Active
Application number
JP2011131498A
Other languages
English (en)
Other versions
JP2013002845A (ja )
Inventor
長井 健太郎
健太郎 長井
Original Assignee
キヤノン株式会社
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
Grant date

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01TMEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/40Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment with arrangements for generating radiation specially adapted for radiation diagnosis
    • A61B6/4035Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment with arrangements for generating radiation specially adapted for radiation diagnosis the source being combined with a filter or grating
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/48Diagnostic techniques
    • A61B6/484Diagnostic techniques involving phase contrast X-ray imaging
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01BMEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B9/00Instruments as specified in the subgroups and characterised by the use of optical measuring means
    • G01B9/02Interferometers for determining dimensional properties of, or relations between, measurement objects
    • G01B9/02097Self-interferometers, i.e. the object beam interfering with a shifted version of itself
    • G01B9/02098Self-interferometers, i.e. the object beam interfering with a shifted version of itself shearing interferometers
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00Investigating or analysing materials by the use of wave or particle radiation not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/20Investigating or analysing materials by the use of wave or particle radiation 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/20075Investigating or analysing materials by the use of wave or particle radiation 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
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/42Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment with arrangements for detecting radiation specially adapted for radiation diagnosis
    • A61B6/4291Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment with arrangements for detecting radiation specially adapted for radiation diagnosis the detector being combined with a grid or grating
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00Investigating materials by wave or particle radiation
    • G01N2223/30Accessories, mechanical or electrical features
    • G01N2223/345Accessories, mechanical or electrical features mathematical transformations on beams or signals, e.g. Fourier

Description

本発明は撮像装置に関し、特にシアリング干渉計を用いて被検体の情報を得る撮像装置、撮像装置に用いられるプログラム、及び解析方法に関する。

X線を含む様々な波長の光の干渉を利用して被検体を撮像、計測する技術が知られている。

この技術について簡単に説明をする。

位相波面の揃った(コヒーレントな)光を被検体に対して照射すると、光は被検体の形状や組成によって波面が変化する。波面に変化が生じた光を何らかの方法を用いて干渉させて干渉パターン(干渉縞)を形成させ、この干渉パターンを解析して光の位相波面を回復すると、被検体の位相、散乱、吸収に関する情報を算出することができる。

シアリング干渉計は上述のような光の干渉を利用して光のシア波面を計測する干渉計であり、シアリング干渉計で検出される干渉パターンは被検体によって生じる光の波面変化を微分した情報を有している。

この技術の代表的な適用例としてはレンズなどの形状測定を行う波面測定技術がある。

また、その他の適用例として、X線を用いて被検体の微分位相像を得る技術もある。

この技術はX線を被写体に照射し、被検体の形状や組成によって生じたX線の位相差を計測する技術であり、被検体の内部構造に関する情報を有する微分位相像を算出することが可能である。

干渉によって得られた干渉パターンから被検体によって生じた光の波面変化を計算する方法は位相回復法と呼ばれる。

位相回復法には幾つかの種類があるが、その一つにフーリエ変換法と呼ばれる方法があり、なかでも非特許文献1に記載されているような、干渉パターンに対して窓関数を掛けてからフーリエ変換を行う方法を窓フーリエ変換法と呼ぶ。

窓フーリエ変換法は、窓関数を用いないフーリエ変換法と比較すると一般的にノイズに対して耐性があるといった特徴がある。

"Windowed Fourier transform method for demodulation of carrier fringes,¨Opt.Eng.43(7)1472−1473(July 2004)

窓フーリエ変換法に用いる窓関数の大きさ(目安としては半値全幅がよく用いられる)が小さいほど、干渉パターンを局所的に周波数成分で表現することができる。そのため、空間分解能が向上する。

しかし、窓関数の大きさが小さいほど、波数空間において隣接するスペクトル同士の重なりが大きくなり、周波数の分解能が低下するという課題が生じる。

このように窓フーリエ変換は、空間分解能と周波数分解能の一方を高くすると他方が低くなるという課題を有する。

そこで、本発明は隣接するスペクトルの重なりの影響を低減することができる計算方法、プログラムと、撮像装置を提供することを目的とする。

その目的を達成するために、本発明の一側面としての撮像装置は、シアリング干渉計と、前記シアリング干渉計により得られる干渉パターンから被検体の情報を算出する算出手段と、を備えた撮像装置であって、前記算出手段は、前記干渉パターンを窓フーリエ変換することにより得られる波数空間上の互いに異なる座標におけるフーリエ成分を表す式を3以上用い、前記フーリエ成分を表す3以上の式を連立方程式として解くことにより、前記被検体の情報を算出し、前記被検体の情報は、前記被検体の吸収像、散乱像、微分位相像、位相像、の少なくとも1つに関する情報であることを特徴とする。

本発明のその他の側面については、以下で説明する実施の形態で明らかにする。

本発明によると、窓フーリエ変換を用いた位相回復を行う際に隣接するスペクトルの重なりの影響を低減することができる計算方法、プログラムと、撮像装置を提供することができる。

実施形態の撮像装置の概略図。 1次元のトールボット干渉計に用いる格子と干渉パターンの例の模式図。 2次元のトールボット干渉計に用いる格子と干渉パターンの例の模式図。 実施形態における位相回復に用いた座標を示した波数空間の模式図。 実施例と比較例のシミュレーションに用いた被検体と検出されるモアレの模式図。 実施例により得られた128×128画素の微分位相像。 比較例により得られた128×128画素の微分位相像。 窓フーリエ変換法における逐次位相変換のイメージ図。

周波数分解能を保ちつつ空間分解能を向上させる又は空間分解能を保ちつつ周波数分解能を向上させるためには、波数空間において隣接するスペクトルの重なりの影響を考慮して位相回復を行えばよいことが本発明の発明者が鋭意検討した結果わかった。スペクトルの重なりの影響を考慮して位相回復を行う方法として、例えば、スペクトルのフィッティングにより隣接するスペクトルを分離して位相回復を行う方法が考えられる。

しかし、スペクトルのフィッティングを行うには多大なデータを扱うことになる。例えば、1000×1000画素の画像を得たとする。それぞれの画素を中心として窓関数を適用し、フーリエ変換を実施すると、元画像の1画素毎に1000×1000画素の波数空間のデータが得られる。この結果、データの合計は1000の4乗となり、スペクトルの分離と位相回復に多大な時間又はコンピュータのリソースを必要とする。

そこで、上述のようなスペクトルのフィッティングによりスペクトルを分離する方法よりも短時間又は低リソースで行うことができる位相回復方法を行う撮像装置について、添付の図面に基づいて以下で説明をする。なお、各図において、同一の部材については同一の参照番号を付し、重複する説明は省略する。

本実施形態ではシアリング干渉計としてトールボット干渉計を使用した撮像装置について説明する。ただし、本特許が適用可能なシアリング干渉計はこれに限られるわけではなく、さまざまな形態のシアリング干渉計に適用可能である。

図1は本実施形態における撮像装置の構成を示した図である。図1に示した撮像装置1は、トールボット干渉計2と、算出手段として計算機610と、を備えている。トールボット干渉計2は、光源としてのX線源110と、X線を回折する回折格子310、X線の一部を遮る遮蔽格子410、X線を検出する検出器510を有している。更に撮像装置1は計算機610による計算結果に基づいた画像を表示する画像表示装置710と接続されて撮像システムを構成している。

以下、各構成について説明をする。

X線源としては、連続X線を出射するX線源を用いても、特性X線を出射するX線源を用いても良く、平行X線(平行光)を出射するX線源を用いても、発散X線(発散光)を出射するX線源を用いても良い。但し、本明細書においてX線とはエネルギーが2以上100keV以下の光を指す。

また、X線源110からのX線は回折格子310で回折されることにより、干渉パターンを形成する必要があるため、X線源110からのX線には干渉パターンを形成できる程度の空間的可干渉性が求められる。

X線源110からのX線は回折格子310により回折され、トールボット距離と呼ばれる所定の距離をおいて明部と暗部が配列した干渉パターンを形成する。但し、本明細書では、X線(光)の強度が大きい所を明部、小さい所を暗部とする。

本実施形態に用いられる回折格子310は位相型の回折格子である。回折格子として振幅型の回折格子を用いることもできるが、位相型の回折格子の方がX線量(光量)の損失が少ないので有利である。

図2(a)は1次元の干渉パターンを形成する位相格子310aの平面上の構成の一例を表した図である。311は位相の基準となる部分、312は311に対してπ分X線の位相が変化する部分である。図2(b)は位相格子310aによって形成される干渉パターン810aの明部811と暗部812を表している。

図3(a)は2次元の干渉パターンを形成する位相格子310bの構成例である。311は位相の基準となる部分、312は311に対してX線の位相をπ変化させる部分である。図3(b)は位相格子310bによって形成される干渉パターン810bの明部811と暗部812を表している。

遮蔽格子410はX線を透過する透過部とX線を遮蔽する遮蔽部が配列した構造を有しており、回折格子310からトールボット距離はなれて配置される。これにより、遮蔽格子410を透過したX線はモアレを形成する。尚、遮蔽部は遮蔽格子を透過したX線がモアレを形成する程度にX線を遮ればよいため、X線を完全に遮らなくても良い。

光としてX線を用いたトールボット干渉計の場合、回折格子により形成される干渉パターンの周期は通常数μmから大きくとも数十μ程度であるのに対し、一般的に用いられるX線の検出器の解像度は数十μmから数百μm程度である。そのためこの干渉パターンを直接検出することは難しい。そこで、本実施形態のように遮蔽格子410を用いてモアレを形成し、このモアレを検出する方法が用いられることが多い。このようにモアレを形成する場合、遮蔽格子のピッチは干渉パターンのピッチと同じか、僅かに異なる値をとることができ、形成したいモアレのピッチによって決めることができる。尚、モアレのピッチは遮蔽格子の遮蔽部と透過部が配列した方向と、干渉パターンの明部と暗部が配列した方向のなす角度によっても変化する。モアレの周期は様々な値をとることができるが、一般的に検出器510の検出素子の3画素分以上が望ましいとされている。

図2(c)は図2(b)の干渉パターンに用いられる遮蔽格子410aの平面上の構成を表した図である。また、図3(c)は図3(b)の干渉パターンに用いられる遮蔽格子410bの平面上の構成を表した図である。図2(c)の遮蔽格子410a、図3(c)の遮蔽格子410b共に、透過部411と遮蔽部412が周期的に配列している。

図2と図3に示した回折格子と遮蔽格子の組み合わせは一例であり、他の組み合わせを用いることもできる。尚、本発明はこれらの格子の構成に依存しない。また、干渉パターンを直接検出する場合は遮蔽格子は不要である。

検出器510は、X線を検出することのできる検出素子(例えばCCD)を有しており、遮蔽格子を透過して形成されたモアレの強度分布を検出する。本実施形態の撮像装置ではモアレの強度分布を検出するが、干渉パターンの強度分布を直接検出し、解析しても良い。尚、本実施形態の場合は説明のために干渉パターンとモアレを区別して説明したが、モアレも干渉パターンの一種である。つまり、本実施形態ではモアレを検出し、検出したモアレを解析するのでモアレを用いて説明をするが、干渉パターンを直接検出する場合もモアレを検出する場合と同様に解析を行うことができる。

計算機610はトールボット干渉計2の検出器510の検出結果に基づいて被検体の微分位相像の情報を計算する。

計算機610によって行われる計算方法(位相回復方法)について説明するために、まず比較例としてスペクトルのフィッティングによりスペクトルを分離して微分位相像の情報を計算する位相回復方法について説明をする。

2次元の窓フーリエ変換は次の式によって定義される。

ここで、f(x、y)は元関数、g(x、y)は窓関数、(x、y)は座標、(u,v)は窓関数の中心、(k,k)は波数を表す。WF[・・・]は括弧内の関数に対して窓フーリエ変換を行うことを示す演算子である。あるモアレの強度分布I(x,y)に対して窓フーリエ変換を行うと、それぞれの窓関数の中心位置(u,v)毎に波数空間(k,k)が得られる。

例えば、離散化した1000×1000画素の画像(本実施形態ではモアレの強度分布)を窓フーリエ変換すると、画像の夫々の画素を窓関数の中心位置とした窓毎に1000×1000画素の波数空間が得られる。つまり、1000×1000画素の波数空間が1000×1000枚得られる。これは、1000×1000画素の画像の情報が、1000の4乗の情報に変換されることを意味する。

これを図にすると図8のようになる。図8(a)はモアレI(x、y)の模式図である。窓関数g(u,v)で切り取られる範囲900は任意の座標(u,v)を中心座標としている。この範囲900内にフーリエ変換を実行すると、図8(b)に示したような波数空間9000が得られる。この波数空間は、0次スペクトルの911、1次スペクトルの912、913、914、915といったスペクトルを有しており、これらのスペクトルから被検体によるX線の波面の位相変化やX線吸収量、X線散乱の情報を計算することができる。尚、1次スペクトルはモアレの周期に由来したスペクトルである。

このような波数空間は通常、窓関数の中心座標(u,v)毎に計算される。

つまり、範囲900から窓関数の中心座標を変化させた範囲901内にフーリエ変換を実行すると波数空間9001が得られる。同様に範囲902内にフーリエ変換を実行すると波数空間9002が、範囲903内にフーリエ変換を実行すると波数空間9003が、範囲904内にフーリエ変換を実行すると9004が得られる。

上述のように窓フーリエ変換において窓関数で切り取られる範囲900の半径を縮小すると、隣接するスペクトルが重なり合ってしまう可能性がある。

そこで、隣接するスペクトルを分離する。スペクトル911〜914は窓関数の形状でスペクトルをフィッティングできると考えられるため、本比較例ではこの方法によるフィッティングを用いてスペクトルを分離する。

例えば、窓関数としてガウシアン窓を使用した場合、そのフーリエ変換は同じガウシアン窓となるため、波数空間上のスペクトルもガウシアン窓でフィッティングすればよい。

しかし、ガウシアン窓でフィッティングするにはすべての(u,v)の組み合わせについて窓フーリエ変換を実行し、図8(b)のように波数空間の全座標におけるフーリエ成分を計算する必要がある。そのため、1000×1000画素の画像を用いて位相回復を行う場合、上述のようにすべての画素を窓関数の中心として窓フーリエ変換を実行し、得られた波数空間毎に「1000×1000画素に相当する波数空間」を用いて位相回復を行う必要がある。そのため位相回復には多大な時間がかかる。特に画像が巨大化すると指数関数的に位相回復に必要な計算時間または計算機のリソースが増す。

そこで、本実施形態では波数空間(k,k)のマップを作成せず、(式1)における(k,k)の組み合わせのうち数点のみのフーリエ成分を計算して位相回復を行うことで、計算量を減らす。

本実施形態の計算機610によって行われる位相回復方法について説明する。
まず、2次元の位相イメージングを仮定し、モアレが近似的に次のような形で記述できると仮定する。

ここで、a(x,y)は被検体による吸収の量を表し、b(x,y)はモアレの振幅を表す。また、P(x,y)およびP(x,y)は測定する位相を表す。これらは場所によって異なる値をとりうる。ωおよびωはそれぞれx方向、y方向のモアレの周期を表す。モアレの形状は(式2)で表される形状に限定されるわけではなく、あくまで一例であり、本実施形態はさまざまなモアレ(干渉パターン)に適用可能である。例えば、モアレの方向が画面のx軸方向、y軸方向に沿っていないモアレは(式2)よりも複雑な式で表される。詳しく説明はしないが、このような場合は回転移動変換などを行うことで最終的に(式2)に帰結することができる。

また、(式2)の第3項を0にすると(式2)は1次元のモアレを示す。以下の説明は1次元のモアレにおいても適用できる。

(式2)を(式1)に代入する。ただし、窓関数g(x,y)の幅は十分に小さく、上記のa(x,y)、b(x,y)、P(x,y)およびP(x,y)はその範囲内で一定値であると近似できるとする。そのため、以下の議論ではそれぞれa,b,P,Pと略記する。また、窓関数g(x,y)のフーリエ変換をG(x,y)と記述する。

(式2)を(式1)に代入すると下記(式3)が得られる。

図4は2次元の位相イメージングの場合に、ある(u,v)を中心座標として窓フーリエ変換を施した時に得られる波数空間(k,k)のマップを示した図である。上述のように本実施形態ではこのようなマップを作成せず、波数空間における数点のみのフーリエ成分を算出するが、ここでは本実施形態の説明のためにこのようなマップを用いる。

ここで、(0,0)は原点であり0次スペクトルのピークの位置を示し(ω,0)、(−ω,0)、(0,ω)及び(0,−ω)は2次元モアレの1次スペクトルのピークを示す。以下図4(a)に示すように(0,0)、(ω,0)と(−ω,0)を用いて位相回復を行う方法を説明する。

この3つの座標におけるフーリエ成分を表す式を用いて位相を回復する。各々の座標におけるフーリエ成分の値は(式3)から以下のよう表すことができる。

尚、(ω,0)と(−ω,0)は原点に対して点対称であり、複素共役の関係にあるので、(式6)は(式5)から求めることもできる。

(式4)(式5)(式6)を連立方程式として解く。
(式4)にexp(−iωu)G(ω,0)をかけて(式5)との差をとると、


となる。

同様に、(式4)にexp(+iωu)G(−ω,0)をかけて(式6)との差をとると、

となる。ただし、ここで上式を導出するためには窓関数のフーリエ変換の特徴として以下の特性を考慮する必要がある。
G(−ω,ω)=G(ω,ω)・・・(y成分も同様)
G(ω,ω)G(ω,ω)=G(ω+ω,ω)・・・(y成分も同様)
このため(式4)の第4項と第5項は(式5)および(式6)の第4項と第5項と互いに打ち消しあう。

このように、(0,0)、(ω,0)、(−ω,0)の3つの座標におけるフーリエ成分を表わす式から(式7)と(式8)が導出できる。

(式7)と(式8)に、検出結果から算出したフーリエ成分の値({WF[I(x,y)](u,v,0,0)}、{WF[I(x,y)](u,v,ω,0)}、{WF[I(x,y)](u,v,−ω,0)}と窓関数のフーリエ変換の値(G(ω,0)、G(0,0)、G(−ω,0)、G(2ω,0)、G(−2ω,0))を代入すれば連立方程式の形でbとPが計算できる。

ここで、(ω,0)と(−ω,0)は原点に対して点対称であり、複素共役の関係にあるため、2つの座標におけるフーリエ成分の値は等しい。よって、(0,0)と(ω,0)又は(0,0)と(−ω,0)の2つの座標のみのフーリエ成分の値を用いることでbとPを計算することが可能である。つまり、本実施形態では、波数空間上の座標におけるフーリエ成分を表す式を3つ用い、それらのフーリエ成分を表す式から導かれる連立方程式を、波数空間上の2つの座標におけるフーリエ成分の値を用いて解く。ここで、波数空間上の2つの座標とは、第1の座標(ここでは原点)と、第1の座標と異なり且つ、第1の座標と原点に対して点対称の関係にない第2の座標(ここでは(ω,0)又は(−ω,0))のことを指す。

図4(b)に示した、(0、0)、(0、ω)と(0、−ω)を用いて同様に計算を行うことでPが求められる。また、aに関しても、求めたbとP、Pを(式4)、(式5)(式6)のいずれかに代入することで得られる。このように、複素共役の関係を利用することで実質3つの座標におけるフーリエ成分の値を用いて(式3)で仮定したa、b、P、Pの4つの値が求められる。

これら、a、b、P、Pから、被検体の吸収像、散乱像、微分位相像を得ることができ、さらに微分位相像を積分すると位相像を得ることができる。

以上より、本実施形態を用いると、比較例のようにスペクトルのフィッティングを行うよりも高速又は低リソースで窓フーリエ変換を利用した位相回復法を行うことができる。

以上の例では0次スペクトルのピークと1次スペクトルのピークを用いた例について説明した。しかし、a、b、P、Pの計算に用いる波数空間上の座標(k,k)の組み合わせや用いる座標の数はこれに限定されない。本実施形態ではフーリエ成分を表わす式を5つ用い、これらの式から導かれる連立方程式を解くことでa、b、P、Pを計算したが、用いる座標の数が増えるほどより正確にa、b、P、Pを計算することができる。例えば、フーリエ成分を表わす式を用いて連立方程式を複数たててP1の値を複数求め、最小二乗法を用いて最終的なP1としても良い。但し、検出器の画素を単位とする窓関数の区間をRとすると、Rを超える数の式を用いても、算出できるa、b、P、Pの正確さは向上しない。窓フーリエ変換によって得られる波数空間は元の窓関数の範囲内の画素の情報のみを含むからである。一方、用いる座標の数が増えるほど計算量は多くなる。そこで、計算に用いる式の数は、5以上R以下が好ましい。尚、1次元のモアレの位相回復を行う場合や、1次元の微分位相像を得たい場合など、Pを求める必要がない場合はフーリエ成分を表わす式を3つ以上用いればよい。この場合も上記のように複素共役の関係を用いれば、2つの座標におけるフーリエ成分の値からa,b,P1の値を計算することができる。本明細書において1次元の微分位相像とは、位相像を1方向に微分して得られる像である。Pを求める必要がないときは、Rを超える数の座標における式を用いても、算出できるa、b、P、の正確さは向上しないため、座標におけるフーリエ成分を表わす式を3以上R以下用いることが好ましい。尚、ガウシアン窓を用いる場合は、ガウシアン窓の分散をσとしたとき、99%の情報が存在する範囲である±3σに含まれる画素を窓関数の区間とする。

また、モアレによっては0次や1次のスペクトルのみならず、さらに高次のスペクトルも存在する。この高次スペクトルのピークを用いた場合でも同様に連立方程式を立てて計算することが可能である。例えば図8(b)で示されたスペクトル916,917,918,919といった2次スペクトルなどを利用する方法も考えられる。

また、用いる座標はスペクトルのピークでなくても良い。但し、フーリエ成分の絶対値が大きい座標を用いる方がノイズの影響を受けにくいので好ましい。

また、用いる座標はX軸またはY軸上にある方が、XまたはY軸上にない座標を用いる場合と比べて計算が簡潔にできるため好ましい。

本実施形態では、計算機610による計算を簡潔にするために複素共役の関係を用いたが、複素共役の関係を用いなくても位相を回復することはできる。この場合、連立方程式に代入するフーリエ成分の値も3つ以上必要になる。

尚、窓フーリエ変換は畳み込みの定理により、以下のように表すことも可能である。窓関数が原点で対称な奇関数、すなわちg(x)=g(−x)であると仮定すれば、

となる。ここでF[・・・]は通常のフーリエ変換、F−1[・・・]は逆フーリエ変換を表す。(式9)から、元関数のフーリエ変換F[f(x,y)]に、波数空間上の窓関数F[g(u−x,v−y)exp[ik(u−x)+ik(v−y)]]を乗じ、その逆フーリエ変換を求めることが、元関数に窓関数を乗じてからフーリエ変換を実行することと同じであることがわかる。本実施形態では(式1)を用いて位相回復を行ったが、(式9)を用いて位相回復を行っても良い。

以上、計算機610による位相回復方法について説明をした。上述のような計算を計算機610で行うには、上述のような計算を実行させるプログラムを計算機610に組み込めばよい。

(実施例)
実施形態に記載した撮像装置を用いて、位相回復を行ったシミュレーション結果を実施例として示す。

シミュレーションは、回折格子として図3(a)に示した位相格子310b、遮蔽格子として図3(c)に示した遮蔽格子410b、検出器として128×128画素の検出器を備えた撮像装置1を用いた。また、被検体は図5(a)に示すような球状の被検体1001を用い、この被検体1001が検出器の検出範囲の中央に設置されているとしてシミュレーションを行った。

図5(b)は図5(a)の被検体と128×128画素の検出器を用いた時に検出器によって検出されるモアレである。この検出結果を用いて、上述した位相回復方法により得られた微分位相像を図6の(a)と(b)に示す。図6(a)がX方向の微分位相像で、図6(b)がY方向の微分位相像である。

同様のシミュレーションを256×256画素と512×512画素の検出結果を用いて行い、実際に計算にかかった時間について、表1にまとめた。計算時間はそれぞれ、0.5秒、0.7秒、1.5秒であった。

(比較例)
実施例と同様に、上述した比較例と同じ方法を用いて位相回復を行ったシミュレーション結果を示す。本比較例の撮像装置は計算機で行われる位相回復方法のみが実施例と異なり、そのほかの構成は実施例と同じである。

用いた計算機のリソース上、一度に窓フーリエ空間全体を計算することは不可能であるため、窓フーリエ変換をすべての(u,v)の組み合わせについて計算を行って逐次(u,v)おける微分位相を決定する方法で実行した。実施例同様、図5(b)に示した検出結果を用いて得られた微分位相像を図7(a)と(b)に示す。図7(a)がX方向の微分位相像で、図7(b)がY方向の微分位相像である。

図6と図7を比較すると、同じような微分位相像が得られていることがわかる。

また、実施例同様に同様のシミュレーションを256×256画素、512×512画素の検出結果を用いて行い、実際に計算にかかった時間について、表1にまとめた。計算時間はそれぞれ、105秒、1044秒、20938秒と画素数に従って指数的に長い時間を要することがわかる。
また、実際に計算にかかった時間を実施例と比較すると、実施例の方が計算時間が短縮されていることが分かる。

本発明の骨子は、窓フーリエ変換により得られるフーリエ成分の値を表わす式を用いて位相を回復することによって、波数空間の全座標におけるフーリエ成分を計算せずに位相回復を行うことである。これにより、窓フーリエ変換による位相回復法を短時間又は低リソースで実行することができる。したがって、その理念に基づく手法であるのであれば本発明の実施の形態は上記のみに限定されない。

1 シアリング干渉計
2 撮像装置
610 計算機

Claims (12)

  1. シアリング干渉計と、前記シアリング干渉計により得られる干渉パターンから被検体の情報を算出する算出手段と、を備えた撮像装置であって、
    前記算出手段は、
    前記干渉パターンを窓フーリエ変換することにより得られる波数空間上の互いに異なる座標におけるフーリエ成分を表す式を3以上用い、前記フーリエ成分を表す3以上の式を連立方程式として解くことにより、前記被検体の情報を算出し、
    前記被検体の情報は、
    前記被検体の吸収像、散乱像、微分位相像、位相像、の少なくとも1つに関する情報であることを特徴とする撮像装置。
  2. 前記算出手段は、
    前記フーリエ成分を表す式を3以上R以下用い、
    前記被検体の1次元の、吸収像、散乱像、微分位相像、位相像、の少なくとも1つに関する情報を算出することを特徴とする請求項1に記載の撮像装置。
    但し、Rとは検出器の画素を単位とする窓関数の区間を示す。
  3. 前記算出手段は、
    前記フーリエ成分を表す式を5以上R以下の座標における前記フーリエ成分を表わす式を用い
    前記算出手段は、
    前記被検体の2次元の、吸収像、散乱像、微分位相像、位相像、の少なくとも1つに関する情報を算出することを特徴とする請求項1に記載の撮像装置。
    但し、Rとは検出器の画素を単位とする窓関数の区間を示す。
  4. 前記フーリエ成分を表す3以上の式のうち少なくとも1つは、
    前記干渉パターンの周期に由来する前記フーリエ成分を表す式であることを特徴とする請求項1乃至のいずれか1項に記載の撮像装置。
  5. 前記フーリエ成分を表す3以上の式の1つが、
    前記波数空間の原点における前記フーリエ成分を表す式であることを特徴とする請求項1乃至のいずれか1項に記載の撮像装置。
  6. 前記算出手段は
    前記波数空間の2以上の座標における前記フーリエ成分を算出した値を用いて前記連立方程式を解くことを特徴とする請求項1乃至に記載の撮像装置。
  7. 前記フーリエ成分を表す式のうち少なくとも2つは、
    前記波数空間の原点に対して点対称の関係にある2つの座標における前記フーリエ成分を表す式であることを特徴とする請求項1乃至のいずれか1項に記載の撮像装置。
  8. 前記フーリエ成分を表す式は、
    前記波数空間のX軸又はY軸上の座標における前記フーリエ成分を表す式であることを特徴とする請求項1乃至のいずれか1項に記載の撮像装置。
  9. 前記フーリエ成分を表す3以上の式は、それぞれ、下記式に示す窓フーリエ変換を定義する式に前記干渉パターンを表す式を代入して得られる式に、前記波数空間上の座標を代入して得られる式であることを特徴とする請求項1乃至8のいずれか1項に記載の撮像装置。
  10. 前記算出手段は、
    前記座標におけるフーリエ成分の値と、
    前記座標における窓関数のフーリエ変換の値とを前記フーリエ成分を表す式に代入することにより、前記連立方程式を解くことを特徴とする請求項1乃至9のいずれか1項に記載の撮像装置。
  11. シアリング干渉計により得られる干渉パターンから被検体の情報を算出するプログラムであって、
    前記干渉パターンを窓フーリエ変換することにより得られる波数空間上の互いに異なる座標におけるフーリエ成分を表す式を3以上用い、前記フーリエ成分を表す3以上の式を連立方程式として解くことにより、前記被検体の情報を算出し、
    前記被検体の情報は、
    前記被検体の吸収像、散乱像、微分位相像、位相像、の少なくとも1つに関する情報であることを特徴とするプログラム。
  12. シアリング干渉計により得られる干渉パターンから被検体の情報を得る方法であって、前記干渉パターンを窓フーリエ変換することにより得られる波数空間上の互いに異なる座標におけるフーリエ成分を表す式を3以上用い、前記フーリエ成分を表す3以上の式を連立方程式として解くことにより、前記被検体の情報を算出し、
    前記被検体の情報は、
    前記被検体の吸収像、散乱像、微分位相像、位相像、の少なくとも1つに関する情報であることを特徴とする方法。
JP2011131498A 2011-06-13 2011-06-13 撮像装置、干渉縞解析プログラム及び干渉縞解析方法 Active JP5885405B2 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2011131498A JP5885405B2 (ja) 2011-06-13 2011-06-13 撮像装置、干渉縞解析プログラム及び干渉縞解析方法

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
JP2011131498A JP5885405B2 (ja) 2011-06-13 2011-06-13 撮像装置、干渉縞解析プログラム及び干渉縞解析方法
PCT/JP2012/064494 WO2012173017A1 (en) 2011-06-13 2012-05-30 Imaging apparatus and program and method for analyzing interference pattern
US14125060 US20140114615A1 (en) 2011-06-13 2012-05-30 Imaging apparatus and program and method for analyzing interference pattern
EP20120729751 EP2718699A1 (en) 2011-06-13 2012-05-30 Imaging apparatus and program and method for analyzing interference pattern

Publications (2)

Publication Number Publication Date
JP2013002845A true JP2013002845A (ja) 2013-01-07
JP5885405B2 true JP5885405B2 (ja) 2016-03-15

Family

ID=46354455

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2011131498A Active JP5885405B2 (ja) 2011-06-13 2011-06-13 撮像装置、干渉縞解析プログラム及び干渉縞解析方法

Country Status (4)

Country Link
US (1) US20140114615A1 (ja)
EP (1) EP2718699A1 (ja)
JP (1) JP5885405B2 (ja)
WO (1) WO2012173017A1 (ja)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101387951B1 (ko) * 2013-05-10 2014-04-22 한국기계연구원 싱글 필드 방식의 엔코더를 이용한 웹 이송 속도 측정 방법
JP2015190776A (ja) * 2014-03-27 2015-11-02 キヤノン株式会社 画像処理装置および撮像システム
JP2015205174A (ja) * 2014-04-10 2015-11-19 キヤノン株式会社 画像処理装置および画像処理装置の制御方法

Family Cites Families (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE69637015T2 (de) * 1995-03-16 2007-08-16 Fei Co., Hillsboro Methode zur rekonstruktion von teilchenwellen in einer teilchen-optischen vorrichtung
WO2005045529A3 (en) * 2003-11-04 2006-01-05 Henry A Hill Characterization and compensation of errors in multi-axis interferometry system
JP2005156403A (ja) * 2003-11-27 2005-06-16 Canon Inc シアリング干渉を利用した測定方法及び装置、それを利用した露光方法及び装置、並びに、デバイス製造方法
JP2007234685A (ja) * 2006-02-28 2007-09-13 Canon Inc 測定装置、当該測定装置を有する露光装置及びデバイス製造方法
WO2008126516A1 (ja) * 2007-04-10 2008-10-23 Naoki Suehiro 送信方法、送信装置、受信方法及び受信装置
JP2009277712A (ja) * 2008-05-12 2009-11-26 Canon Inc 測定装置および露光装置
WO2010050483A1 (ja) * 2008-10-29 2010-05-06 キヤノン株式会社 X線撮像装置およびx線撮像方法
US8195435B2 (en) * 2008-12-19 2012-06-05 Tokyo Electron Limited Hybrid diffraction modeling of diffracting structures
US7949095B2 (en) * 2009-03-02 2011-05-24 University Of Rochester Methods and apparatus for differential phase-contrast fan beam CT, cone-beam CT and hybrid cone-beam CT
JP5328437B2 (ja) * 2009-03-25 2013-10-30 キヤノン株式会社 透過波面測定方法、屈折率分布測定方法、光学素子の製造方法、及び透過波面測定装置
US20120116703A1 (en) * 2009-04-24 2012-05-10 Nicolas Pavillon Method and apparatus for enhanced spatial bandwidth wavefronts reconstructed from digital interferograms or holograms
JP2011131498A (ja) 2009-12-24 2011-07-07 Micro-Tec Co Ltd スクリーン印刷機及びスクリーン印刷方法
JP5538936B2 (ja) * 2010-02-10 2014-07-02 キヤノン株式会社 解析方法、プログラム、記憶媒体、x線位相イメージング装置
WO2011121962A4 (en) * 2010-03-31 2012-01-12 Canon Kabushiki Kaisha Optical coherence tomographic imaging apparatus and control apparatus therefor

Also Published As

Publication number Publication date Type
JP2013002845A (ja) 2013-01-07 application
WO2012173017A1 (en) 2012-12-20 application
US20140114615A1 (en) 2014-04-24 application
EP2718699A1 (en) 2014-04-16 application

Similar Documents

Publication Publication Date Title
Nugent Coherent methods in the X-ray sciences
Kemao Two-dimensional windowed Fourier transform for fringe pattern analysis: principles, applications and implementations
Marone et al. Regridding reconstruction algorithm for real-time tomographic imaging
Su et al. Dynamic 3-D shape measurement method based on FTP
US7034949B2 (en) Systems and methods for wavefront measurement
Wieneke Iterative reconstruction of volumetric particle distribution
WO2011139213A1 (en) Touch determination by tomographic reconstruction
Wunderlich et al. Image covariance and lesion detectability in direct fan-beam x-ray computed tomography
Pen et al. Detection of dark matter skewness in the virmos-descart survey: implications for ω0
WO2012050510A1 (en) Touch determination by tomographic reconstruction
Olivo et al. Modelling of a novel x-ray phase contrast imaging technique based on coded apertures
US7574249B2 (en) Device-less gating of physiological movement for improved image detection
WO2010050483A1 (ja) X線撮像装置およびx線撮像方法
US20100322380A1 (en) X-ray detector for phase contrast imaging
WO2013062471A2 (en) Touch determination by tomographic reconstruction
JP2007203074A (ja) 投影または断層撮影による位相コントラスト画像の作成方法
Masajada et al. Vortex points localization problem in optical vortices interferometry
Wielgus et al. Evaluation of amplitude encoded fringe patterns using the bidimensional empirical mode decomposition and the 2D Hilbert transform generalizations
Thüring et al. Non-linear regularized phase retrieval for unidirectional X-ray differential phase contrast radiography
Rosales-Guzmán et al. Experimental detection of transverse particle movement with structured light
Hadrava Disentangling of spectra of multiple stars
Pellegrini et al. Quantifying the contribution of arm postural tremor to the outcome of goal-directed pointing task by displacement measures
Trusiak et al. Two-shot fringe pattern phase-amplitude demodulation using Gram-Schmidt orthonormalization with Hilbert-Huang pre-filtering
US20150131777A1 (en) Medical imaging system and medical image processing apparatus
Trusiak et al. Hilbert-Huang processing for single-exposure two-dimensional grating interferometry

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20140612

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20150623

A521 Written amendment

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20150824

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: 20160112

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20160209