撮像装置および画像処理方法
【課題】位相イメージングにおいて、ノイズが低減された高品位の位相像あるいは微分位相像を生成可能な技術を提供する。
【解決手段】干渉パターンのデータからｘ方向とｙ方向の微分位相データを抽出し、抽出されたｘ方向の微分位相データをｘ方向に関して微分しｘ方向の二階微分位相データを算出し、ｙ方向の微分位相データをｙ方向に関して微分しｙ方向の二階微分位相データを算出する。そして、ｘ方向とｙ方向の二階微分位相データを関数として含む二階微分方程式を解くことによって、被検体の位相像のデータを算出する。
【解決手段】干渉パターンのデータからｘ方向とｙ方向の微分位相データを抽出し、抽出されたｘ方向の微分位相データをｘ方向に関して微分しｘ方向の二階微分位相データを算出し、ｙ方向の微分位相データをｙ方向に関して微分しｙ方向の二階微分位相データを算出する。そして、ｘ方向とｙ方向の二階微分位相データを関数として含む二階微分方程式を解くことによって、被検体の位相像のデータを算出する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明はシアリング干渉計により得られる干渉パターンから位相像や微分位相像のデータを生成するための技術に関する。
【背景技術】
【０００２】
シアリング干渉計とは、同一光源から照射された光線などのコヒーレントな光を分割し、一方に被検体による波面の歪曲を生成し、もう一方をわずかにずらすことで干渉縞を形成し、この干渉縞のずれから被検体による光の位相変化を可視化する手段である。このようなシアリング干渉計は、例えば位相差顕微鏡などの分野で被検体の視認性を向上させる目的などに利用されている。
【０００３】
シアリング干渉計の一つとして、トールボット干渉法を利用したものが知られている。トールボット干渉法では、光源からの光を被検体に照射し、被検体を透過することで位相が変化した光を回折格子に導く。回折格子により回折された光は、回折格子からトールボット距離と呼ばれる所定の距離だけ離れた位置に干渉パターンを形成する。Ｘ線を用いた場合、この干渉パターンの周期は非常に小さいため、干渉パターンを直接検出することが難しいことがある。そこで、干渉パターンが形成される位置に吸収格子を配置し、吸収格子によって干渉パターンの一部を遮ることでモアレを形成し、このモアレ像を検出器で検出する方法が提案されている。このモアレ像を解析することにより、被検体の位相の情報を得ることができる。このように光の位相を利用したイメージングは一般に位相イメージング（Ｘ線を用いる場合はＸ線位相イメージング）と呼ばれている。
【０００４】
ただし、干渉パターンによって直接解析できるのは、あくまで被検体を透過することによって変化した光の位相の空間的な微分量（微分位相像、シア像とも呼ばれる）である。空間的な微分はエッジを強調する効果があるため、微分位相像は光学顕微鏡などの分野では被検体像を強調するために用いられる。また、この微分位相像から元の位相を計算し、位相像を得ることも可能である。
【０００５】
非特許文献１には、二次元トールボット干渉法を用いたＸ線位相イメージング手法が開示されている。具体的には、二次元の回折格子により形成された干渉パターンからＸ方向の微分位相像とＹ方向の微分位相像を取得し、それらに対しフーリエ積分を行うことによって、位相像を計算している。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】C. Kottler, C. David, F. Pfeiffer and O. Bunk, “A two-directional approach for grating based differential phase contrast imaging using hard x-rays,” Opt. Express 15 (2007) 1175-1181.
【発明の概要】
【発明が解決しようとする課題】
【０００７】
シアリング干渉計を利用した位相イメージングは、画像診断や各種の非破壊検査への応用が可能である。中でもＸ線位相イメージングは、ＣＴ（コンピュータ・トモグラフィ）などＸ線吸収量を画像化する従来の手法に比べて被ばく量を大幅に低減できる、筋肉や内蔵などの軟組織も画像化できるといった利点から、新たな画像診断手法として期待されている。
【０００８】
この種の位相イメージングにおいては、ノイズやアーチファクトの低減による画質向上が基本的な課題として存在する。特に、光源としてＸ線を用いる場合には、量子ノイズが不可避であり、しかも被ばく量低減のために線量を下げるほどノイズが増える傾向にあるため、ノイズ対策は重要となる。また、Ｘ線ではなく可視光を用いる装置であっても、検査対象によっては（例えば眼球）低い光量しか用いることができない場合があり、ＳＮ比の悪化が問題となる。
【０００９】
また、位相イメージングにおいては「位相折りたたみ（位相ラッピング）」という問題があることが知られている。位相折りたたみは、干渉パターンから得られる微分位相の値が２πの範囲に折りたたまれて（ラップされて）しまう現象である。例えば−πからπの間で折りたたまれる場合、実際には３／２πという位相が正しい値であったとしても、それは−１／２πとして測定される。非特許文献１においてはこのような位相折りたたみの問題は無視されている。そのため、非特許文献１の手法で生成した位相像には、本来何もない、すなわち位相が一定になる空間に、位相のゆがみやアーチファクトが出現する場合がある。
【００１０】
位相折りたたみへの対策として、「位相接続（位相のアンラッピング）」と呼ばれる処理が知られている。位相接続とは、位相が折りたたまれた箇所を検出し、位相を元に戻す操作である。位相接続に関しては様々な手法が提案されている。最も基本的な操作は隣接ピクセルとの位相の差を計算する手法である。その差の絶対値が一定以上（例えば３／２π以上）であるならば、そこはもともとスムーズにつながっていた位相が折りたたまれたものと疑うことができ、位相を元に戻す操作をする。具体的には位相飛びが生じているピクセル以降のピクセル全ての位相に全て２πか、−２πを足すことで位相隣接ピクセルとの位相差を３／２π以下にする操作を行う。従来の位相接続は、積分とは独立した操作として、微分位相像全体に対して適用する必要があるため、位相像を得るまでに要する処理時間が長くなるという問題があった。
【００１１】
本発明は上記実情に鑑みてなされたものであり、位相イメージングにおいて、ノイズが低減された高品位の位相像あるいは微分位相像を生成可能な技術を提供することを目的とする。
【００１２】
本発明のさらなる目的は、位相イメージングにおいて、位相像あるいは微分位相像におけるノイズの低減、並びに、位相折りたたみによる位相のゆがみやアーチファクトの抑制を、簡易な処理で実現するための技術を提供することにある。
【課題を解決するための手段】
【００１３】
本発明の第１態様は、シアリング干渉計により得られる干渉パターンに基づいて位相像のデータを生成する撮像装置であって、被検体を透過した電磁波により形成された干渉パターンのデータから、第１の方向に関する位相の変化を表す第１の微分位相データと、前記第１の方向に交差する第２の方向に関する位相の変化を表す第２の微分位相データとを抽出する微分位相データ抽出手段と、抽出された前記第１の微分位相データを前記第１の方向に関して微分することによって第１の二階微分位相データを算出すると共に、抽出された前記第２の微分位相データを前記第２の方向に関して微分することによって第２の二階微分位相データを算出する二階微分位相データ算出手段と、前記第１および第２の二階微分位相データを関数として含む二階微分方程式を解くことによって、前記被検体の位相像のデータを算出する位相データ算出手段と、を有する撮像装置である。
【００１４】
本発明の第２態様は、シアリング干渉計により得られる干渉パターンに基づいて位相像のデータを生成する画像処理方法であって、情報処理装置が、被検体を透過した電磁波により形成された干渉パターンのデータから、第１の方向に関する位相の変化を表す第１の
微分位相データと、前記第１の方向に交差する第２の方向に関する位相の変化を表す第２の微分位相データとを抽出する工程と、情報処理装置が、抽出された前記第１の微分位相データを前記第１の方向に関して微分することによって第１の二階微分位相データを算出すると共に、抽出された前記第２の微分位相データを前記第２の方向に関して微分することによって第２の二階微分位相データを算出する工程と、情報処理装置が、前記第１および第２の二階微分位相データを関数として含む二階微分方程式を解くことによって、前記被検体の位相像のデータを算出する工程と、を有する画像処理方法である。
【００１５】
本発明の第３態様は、上記画像処理方法の各工程を情報処理装置に実行させるプログラムである。
【発明の効果】
【００１６】
本発明によれば、位相イメージングにおいて、ノイズが低減された高品位の位相像あるいは微分位相像を生成することができる。また、位相イメージングにおいて、位相像あるいは微分位相像におけるノイズの低減、並びに、位相折りたたみによる位相のゆがみやアーチファクトの抑制を、簡易な処理で実現することも可能となる。
【図面の簡単な説明】
【００１７】
【図１】撮像装置の構成を示す図である。
【図２】回折格子、干渉パターン、吸収格子の構成例を示す図である。
【図３】情報処理装置の機能構成を示す図である。
【図４】干渉パターンから微分位相像を抽出する手法を示す図である。
【図５】従来例における微分位相像から位相像を取得する処理を示す図である。
【図６】実施例における微分位相像から位相像を取得する処理を示す図である。
【図７】シミュレーションで用いた被検体と微分位相像を示す図である。
【図８】従来例の手法で得られた位相像を示す図である。
【図９】実施例の手法で得られた位相像を示す図である。
【図１０】実施例の手法で得られた位相像から生成した微分位相像を示す図である。
【発明を実施するための形態】
【００１８】
以下、本発明の好ましい実施の形態について添付の図面に基づいて説明する。なお、各図において、同一の部材については同一の参照番号を付し、重複する説明は省略する。
【００１９】
本実施形態では微分位相像を得る手段としてＸ線トールボット干渉計を使用した例について説明する。ただし、本発明が適用可能な微分位相像の取得手法はこれに限られるわけではなく、被検体（被検物）の微分位相像を得ることができれば、如何なる形態のシアリング干渉計を用いてもよい。また測定に用いる光はＸ線に限らず、如何なる波長の電磁波を用いてもよい。なお本明細書においてＸ線とはエネルギーが２以上２００ｋｅＶ以下の光を指す。
【００２０】
（撮像装置の構成）
図１は本実施形態における撮像装置の構成を示した図である。図１に示した撮像装置１は、概略、干渉計１０と情報処理装置１１とから構成される。干渉計１０は、２次元Ｘ線トールボット干渉計であり、Ｘ線を発生させるＸ線源（光源）１１０、Ｘ線を回折する回折格子３１０、Ｘ線の一部を遮る吸収格子４１０、Ｘ線を検出する検出器５１０を備えている。情報処理装置１１は、演算部（コンピュータ本体）６１０と画像表示装置７１０を備えている。情報処理装置１１は干渉計１０にケーブルで接続されており、干渉計１０の各部の制御を行う。また干渉計１０の検出器５１０で得られた干渉パターンのデータは情報処理装置１１に伝送され、演算部６１０において後述する各種の演算処理に供される。演算部６１０による演算結果に基づいた画像は画像表示装置７１０に出力される。
【００２１】
測定対象である被検体２１０は、Ｘ線源１１０と回折格子３１０の間に配置される。Ｘ線源１１０から照射され被検体２１０を透過したＸ線は回折格子３１０により回折され、トールボット距離と呼ばれる所定の距離Ｌの面上に明部と暗部が配列方向に並んだ干渉パターン８１０を形成する。本明細書では、Ｘ線（光）の強度が大きい所を明部、小さい所を暗部とする。本実施形態に用いられる回折格子３１０は位相型の回折格子である。回折格子として振幅型の回折格子を用いることもできるが、位相型の回折格子の方がＸ線量（光量）の損失が少ないので有利である。
【００２２】
図２（ａ）は位相型の回折格子３１０（以下単に位相格子３１０と呼ぶ）の平面上の構成の一例を表した図である。３１１は位相の基準となる領域、３１２は基準の領域３１１に対してπほどＸ線の位相が変化する領域である。図２（ｂ）は位相格子３１０によってトールボット距離に生じる干渉パターン８１０の明部と暗部を表している。これらの組み合わせには何種類かあるが、本発明では公知のものを含めて如何なる構成のものを用いてもよいため、ここでは詳しい説明を割愛する。
【００２３】
この干渉パターン８１０から微分位相像を抽出する方法には大きく分けて二種類考えられる。一つは、吸収格子４１０を少しずつずらしながら複数回撮像し、その結果から微分位相像を計算する縞走査法である。もう一つは、自己像としての干渉パターン８１０の周期と少しずらした吸収格子４１０を利用して撮像する方法で、フーリエ変換法と呼ばれる。以下ではフーリエ変換法に関して説明するが、どちらの手法でも微分位相像を取得できる。なお吸収格子は、遮蔽格子または振幅格子とも呼ばれることもある。
【００２４】
フーリエ変換法では図２（ｃ）で示すような吸収格子４１０を用いる。吸収格子４１０には、干渉パターン８１０の明部に対応する透過部４１１と、暗部に対応する遮蔽部４１２とが設けられている。ただし、透過部４１１の周期（ピッチ）は干渉パターン８１０の明部の周期とは異なる。このため、吸収格子４１０を透過した光像を実際に撮像すると、干渉パターン８１０と吸収格子４１０の周期の違いによるモアレが観測される。
【００２５】
通常、位相格子３１０による干渉パターン８１０の周期は数μｍから十数μｍ程度であるが、吸収格子４１０によりモアレを作ることにより、その周期を数十μｍ以上に拡大することができる。これによって数十μｍ程度の分解能の検出器５１０で干渉パターンを撮像することができる。しかしながら、例えば、検出器５１０の分解能が十分高い場合や、位相格子３１０と検出器５１０の距離が大きい場合など、検出器５１０で干渉パターン８１０そのものを解像可能な場合には、吸収格子４１０を設けなくてもよい。
【００２６】
（情報処理装置の機能）
図３は、情報処理装置１１が有する機能を示すブロック図である。情報処理装置１１は、概略、干渉計１０の制御（例えばＸ線の照射、検出器５１０からのデータ取り込みなど）を行う制御部３０と、干渉計１０から取り込まれたデータに基づいて位相像のデータを生成するデータ処理部３１とを有する。データ処理部３１は、微分位相データ抽出部３２、二階微分位相データ算出部３３、及び、位相データ算出部３４から構成されている。微分位相データ抽出部３２は、干渉パターン（モアレ像）のデータを数値的に解析して微分位相データを抽出する機能であり、二階微分位相データ算出部３３は、微分位相データを微分することによって二階微分位相データを算出する機能である。位相データ算出部３４は、二階微分位相データを関数として含む二階微分方程式を解くことによって、被検体の位相像のデータを算出する機能である。これらの機能の詳細については後述する。
【００２７】
情報処理装置１１の演算部６１０は、ＣＰＵ（中央演算処理装置）、主記憶装置（ＲＡＭなど）、補助記憶装置（ＨＤＤ、ＳＳＤなど）、各種Ｉ／Ｆを有するコンピュータで構
成することができる。図３に示した機能は、補助記憶装置に格納されたプログラムが主記憶装置にロードされ、ＣＰＵにより実行されることで実現されるものである。もちろん、この構成はあくまで一例であり、本発明の範囲を限定するものではない。
【００２８】
（撮像装置の動作）
以下、本実施形態の撮像装置１において実行される画像処理方法、並びに従来例との差異を説明する。ただしここでは、干渉パターンのデータとして、干渉計１０により実際に被検物を計測して得たデータではなく、コンピュータシミュレーションによって生成したデータを用いる。干渉パターンのデータから位相像のデータを生成する一連の処理については、情報処理装置１１のデータ処理部３１によって行った。
【００２９】
シミュレーションでは、Ｘ線として１７．５ｋｅＶの単色光を仮定し、位相格子として図２（ａ）に示す構造の位相格子３１０を仮定した。位相格子３１０には、Ｘ線を十分透過させる基盤部３１１に対して柱状の構造物３１２がピッチｐ１で並んでいると仮定する。構造物３１２を透過したＸ線は、基盤部３１１を透過したＸ線に対して、位相がπ変化するようにした。このような構造の位相格子３１０は例えば構造物３１２がシリコンで作成された場合、２２．６μｍ程度の奥行きにすることによって実現される。
【００３０】
さらに図２（ａ）においてピッチｐ１を８μｍと仮定すると図１における位相格子３１０と吸収格子４１０の距離Ｌを約１１．３ｃｍに設定することにより、吸収格子４１０の直前で干渉パターン８１０が形成される。この干渉パターン８１０の周期ｐ２は４μｍである。吸収格子４１０によりモアレを形成することで、干渉パターン８１０を拡大し、検出器５１０でその強度を検出する。簡単のため、検出器５１０に含まれる画素素子の大きさは４μｍとして計算した。
【００３１】
トールボットのような干渉計の場合、入射光には高いコヒーレンス性が必要とされる。そのような高いコヒーレンス性を実験室や医療装置のような大きさがせいぜい数ｍ程度までの施設で用いる場合、微小な点光源を用いてコヒーレンス性を出す手法が用いられる。このような場合、入射光は平行波ではなく球面波になる。しかし、この違いは本発明の本質を変えるものではないため、本発明は平行光でも非平行光（球面波など）でも適用可能である。以下の実施例では平行光を仮定してシミュレーションを行った。
【００３２】
被検体としては二種類用意した。一つ目は図７（ａ）に示すような直径４００μｍのカルシウム球を想定した。もう一つは図７（ｂ）に示すような直径２００μｍのカルシウム球が４つ重なっている被検体である。計算の範囲は５１２μｍ×５１２μｍの範囲となる。
【００３３】
微分位相像の抽出方法としては前述の通りフーリエ変換法を使用した。そのために、吸収格子４１０の透過部４１１の周期ｐ３（図２（ｃ）参照）は４．５μｍとした。これによりほぼ４画素で一周期のモアレが実現することとなる。
【００３４】
図４（ａ）から図４（ｄ）はフーリエ変換法により干渉パターンから微分位相像を抽出する方法、すなわち微分位相データ抽出部３２の処理内容を概念的に示した図である。図２（ｂ）のような干渉パターン８１０を拡大した結果、図４（ａ）のようなモアレ像が取得できたとする。フーリエ変換法ではこのモアレ像をフーリエ変換する。これによってフーリエ空間内では９つのスペクトルに分離される。この９つのスペクトルのうち１次スペクトルにあたるスペクトル９１１、９１２、９１３、９１４のいずれかの周辺を窓関数で切り取って原点に戻し（図４（ｃ））、これを逆フーリエ変換することで位相微分像を得ることができる。図４（ｄ）の左側は、スペクトル９１４から生成された微分位相像の例であり、ｘ方向（第１の方向）に関する位相の変化を表している。図４（ｄ）の右側は、
スペクトル９１１から生成された微分位相像の例であり、ｙ方向（第２の方向）に関する位相の変化を表している。
【００３５】
図７（ｃ）は、図７（ａ）の被検体から得られたｘ方向、ｙ方向それぞれの微分位相像を示す。また図７（ｄ）は、図７（ｂ）の被検体から得られたｘ方向、ｙ方向それぞれの微分位相像を示す。球体のエッジ部分は最も急激に位相の微分が変化する個所となるため、この部分で位相飛びが発生しやすい。また、フーリエ変換法では、前述の窓関数の取り方によっては微分位相像にアーチファクト（本来は存在しない像）が発生する場合がある。図７（ｃ）、および図７（ｄ）の微分位相像でも若干のアーチファクトが見られる。
【００３６】
二次元の場合はｘ方向とｙ方向の微分位相像がそれぞれ独立に取得できる。この方向の違いは、図４（ｂ）のフーリエ空間においてどのスペクトルを切り取るかによって制御できる。例えばスペクトル９１１又は９１２を切り取ればｙ方向の微分位相像が得られ、スペクトル９１３又は９１４を切り取ればｘ方向の微分位相像が得られる。もちろんｘ、ｙの直交する２方向ではなく、互いに交差する任意の２方向についての微分位相像を得ることも可能である。以下では取得されたｘ方向微分位相像をφｘ（ｘ，ｙ）、ｙ方向微分位相像をφｙ（ｘ，ｙ）とおく。（ｘ，ｙ）は座標である。これらのｘ方向微分位相像とｙ方向微分位相像から、積分された位相像を求める。
【００３７】
（１）従来の位相像取得方法
従来例では、ｘ方向微分位相像とｙ方向微分位相像を別々に積分するか、フーリエ変換などの手段を用いて積分する。以下ではより誤差の少ない位相が求められるフーリエ変換による積分法（フーリエ積分法）を従来例として述べる。
【００３８】
図５はフーリエ積分法による位相像取得方法の手順を示すフローチャートである。まず前述のように、ｘ方向微分位相像φｘ（ｘ，ｙ）とｙ方向微分位相像φｙ（ｘ，ｙ）を取得する（ステップＳ５０）。ｇ（ｘ，ｙ）を次の式で定義する（ステップＳ５１）。
【数１】
ここで、ｉは虚数単位を示す。
【００３９】
ｇ（ｘ，ｙ）をフーリエ変換することでＧ（ｋｘ，ｋｙ）を得て（ステップＳ５２）、次の（式３）によって位相Ｐ（ｘ，ｙ）を取得する（ステップＳ５３、Ｓ５４）。
【数２】
ここで、Ｆ［・・・］は括弧の中をフーリエ変換することを示す演算子であり、Ｆ−１［・・・］は逆フーリエ変換を示す演算子である。ｋｘとｋｙはそれぞれｘ方向、ｙ方向の波数を表す。（式３）はすなわち（式１）で定義されたｇをフーリエ変換し積分演算子１／（２πｉ×（ｋｘ＋ｋｙ））を作用させて逆フーリエ変換をかけることと一致する。こうして得たＰ（ｘ，ｙ）は被検体の位相像を表す（ステップＳ５５）。
【００４０】
図８（ａ）は、図７（ｃ）で示した微分位相像より、図５に示した従来例による方法で取得した位相像の例である。カルシウム球の位相が球形に回復されているのがわかる。しかし、位相の折りたたみにより被検体の外周部で位相が不連続になっており、取得した位相像にゆがみが生じている。特に何もない真空部（一定）であるはずの被検体の外周部に
、本来ありえない位相のゆがみが生じていることがわかる。
【００４１】
図８（ｂ）は、図７（ｄ）で示した微分位相像より、図５に示した従来例による方法で取得した位相像の図である。４つのカルシウム球の位相が復元されていることがわかる。しかし、先の例と同様に本来は真空であるはずの領域に位相のゆがみができている。また、アーチファクトの影響で本来は球形である位相がジャギーを持っているように見える。
【００４２】
（２）実施例の位相像取得方法
図６は、本発明の実施例における位相像取得方法の手順を示すフローチャートである。従来例では、微分位相像を直接積分することで位相像を得ていたのに対し、本実施例では、微分位相像をｘ，ｙそれぞれの方向について微分し二階微分位相像を生成した後、二階微分位相像からなる二階微分方程式を解いて位相像を得る点に特徴がある。
【００４３】
まず、微分位相データ抽出部３２が、干渉パターンのデータ（モアレ像のデータ）から、ｘ方向微分位相像φｘ（ｘ，ｙ）とｙ方向微分位相像φｙ（ｘ，ｙ）を取得する（ステップＳ６０）。微分位相データ抽出部３２が生成したｘ方向微分位相像のデータ（第１の微分位相データ）およびｙ方向微分位相像のデータ（第２の微分位相データ）は、二階微分位相データ算出部３３に引き渡される。
【００４４】
二階微分位相データ算出部３３は、ｘ方向微分位相像をｘ方向に関して微分することによってｘ方向二階微分位相像を算出すると共に、ｙ方向微分位相像をｙ方向に関して微分することによってｙ方向二階微分位相像を算出する（ステップＳ６１）。ｘ方向二階微分位相像ψｘ（ｘ，ｙ）、ｙ方向二階微分位相像ψｙ（ｘ，ｙ）は、下記式で表される。
【数３】
【００４５】
ｘ方向二階微分位相像のデータ（第１の二階微分位相データ）およびｙ方向二階微分位相像のデータ（第２の二階微分位相データ）は、位相データ算出部３４に引き渡される。
【００４６】
位相データ算出部３４は、既知関数としてｘ方向、ｙ方向それぞれの二階微分位相像を含む二階微分方程式を解くことで、未知関数であるｘ方向、ｙ方向それぞれの位相像を求める（ステップＳ６２〜Ｓ６６）。そのアルゴリズムについて以下詳しく述べる。
【００４７】
まず位相データ算出部３４は、ｘ方向とｙ方向の二次微分位相像を用いて、ρ（ｘ，ｙ）を次のように定義する（ステップＳ６２）。
【数４】
ここで、Ｗ［・・・］は、［・・・］の値のうち、所定の範囲を外れた値をその所定の範囲内に収める処理を行う演算子である。微分位相像φｘ、φｙの値は基本的に連続的に変化するため、二階微分位相像ψｘ、ψｙの値は０を中心とする狭い範囲（せいぜい−π〜πの範囲）の内に収まる。しかし、位相の折りたたみが生じると、微分位相像φｘ、φｙの値が−πからπ、又はπから−πへと不連続に変化するため、二階微分位相像ψｘ、ψｙの値が２π又は−２πとなる。すなわち、位相の折りたたみが生じた箇所は、二階微分位相像ψｘ、ψｙにおいてスパイク状の異常値となって現れる。Ｗ［・・・］の操作は
、この異常値を除去（修正）する操作、すなわち、位相の折りたたみが生じた状態を元の位相の状態に戻す操作（位相接続）に相当する。以下、Ｗ［・・・］をアンラッピング演算子と呼ぶ。
【００４８】
アンラッピング演算子にはさまざまなデザインが考えられるが、本実施例では、具体的かつもっとも簡単な演算子として、次のルールを持つ演算子Ｗ［・・・」を用いる。
（ａ）［・・・］の値がπ以上のものに関して、その値から２πを引く
（ｂ）［・・・］の値が−π以下のものに関して、その値に２πを足す
以上のルールを適用することで元関数φｘ（ｘ，ｙ）やφｙ（ｘ，ｙ）に位相折りたたみが生じていても、位相接続がなされたのと同じ効果を持つ。
【００４９】
実際の位相接続では、上記のみならずノイズによる特異点の処理などを含めた様々な手法が存在する。Ｗ［・・・］もそれらの処理方法を応用した様々な演算アルゴリズムが考えられるが、本発明はそれらの手法に本質的には依存せず様々な方法を適用可能である。
【００５０】
位相データ算出部３４は、（式６）のρ（ｘ，ｙ）から次のポワソン方程式を解くことによって位相Ｐ（ｘ，ｙ）を計算する。
【数５】
【００５１】
このポワソン方程式を解く具体的アルゴリズムも様々な手法が考えられるがその一つとしては、フーリエ変換による二階の積分法が考えられる。すなわち、位相データ算出部３４は、ρ（ｘ，ｙ）をフーリエ変換することでＨ（ｋｘ，ｋｙ）を得て（ステップＳ６３）、次の（式９）によって位相Ｐ（ｘ，ｙ）を算出する（ステップＳ６４、Ｓ６５）。
【数６】
（式９）では（式６）で定義されたρ（ｘ，ｙ）に対してフーリエ変換を求め、１／（４π２（ｋｘ２＋ｋｙ２））の二階の積分演算子を作用させた後に逆フーリエ変換して、取得した値を位相Ｐ（ｘ，ｙ）とする。以上の処理により、二階微分位相像の二階積分である位相像を得ることができる。
【００５２】
図９（ａ）は、図７（ｃ）で示した微分位相像より、図６に示した本実施例の手法によって取得した位相像の例である。従来例の図８（ａ）と比べると、位相折りたたみによる歪みが改善されていることがわかる。これは従来例では位相が不連続になっていた箇所が、ある程度修復された影響によるものであり、本実施例の手法の有効性を示している。
【００５３】
図９（ｂ）は、図７（ｄ）で示した微分位相像より、図６に示した本実施例の手法によって取得した位相像の例である。先の例と同様に、従来例の図８（ｂ）と比べると位相折りたたみによる歪みが改善されており、本実施例の手法の有効性を示している。また、従来例の図８（ｂ）では目立っているアーチファクトについても、図９（ｂ）では低減していることがわかる。
【００５４】
また本実施例の手法はノイズ耐性に関しても優れている。表１は、従来例の位相像取得方法と本実施例の位相像取得方法とで、同じノイズに対するノイズリダクションの機能を
比べた結果である。入力画像として平均０、分散が０．５のランダムパターンの画像２枚を微分位相像とみなして位相を計算し、その位相像にどれだけのノイズが存在するかを比較した表である。この実験では、純粋にノイズリダクションを評価するため、被検体を含まない画像を用いている。表１の数値は分散であり、値が小さいほどノイズが少ない（ノイズが低減されている）ことを示す。１０回の計算を行なった結果、すべての回において本実施例の手法でノイズが減少していた。
【表１】
【００５５】
このように本実施例の手法は従来の手法よりもノイズ低減の効果もあることが分かる。この実験では被検体を含まない画像を用いたので、このノイズ低減の効果は、前述したアンラッピング演算子Ｗによる位相接続による画像の改善とは、本質的に別の効果である。すなわち位相接続を行わなくとも（つまり（式６）において演算子Ｗを省略しても）、本実施例の手法が従来の手法に比べて有利な効果を有していることがわかる。
【００５６】
図１０（ａ）、図１０（ｂ）はノイズ低減効果の応用例を示している。図１０（ａ）は、図７（ａ）の被検体の像に対し分散０．２５のノイズを重畳した場合のｘ方向及びｙ方向の微分位相像である。図１０（ｂ）は、図１０（ａ）の微分位相像に対し本実施例の位相像取得方法を適用した後、位相像をｘ方向、ｙ方向それぞれに関し微分することで生成した微分位相像である。図１０（ａ）の元の微分位相像においてみられたアーチファクトやノイズが、図１０（ｂ）では改善されており、高品位な微分位相像が得られている。すなわち本実施例の手法は、微分位相像の画質を向上するという目的にも利用できる。画像診断や検査の内容によっては、位相像をみるよりも、エッジが強調された微分位相像のほうが有用な場合もある。よって、高画質の微分位相像を生成できる機能も実用上非常に有利である。構成としては、例えば、図３の位相データ算出部３４の後段に、位相像のデータを微分する微分位相データ算出部を設ければよい。
【００５７】
以上述べたように、本実施例の位相像取得方法によれば、位相イメージングにおいて、ノイズが低減された高品位の位相像および微分位相像を生成することが可能である。そして、位相接続処理を施したことにより、位相折りたたみに起因する位相のゆがみやアーチファクトを改善することもできる。しかも本実施例では、二階微分位相像に対して位相接続処理を施す構成としたので、従来（微分位相像に対して位相接続処理を施す）に比べて、位相接続処理が格段に簡単となり、また高速に処理できる。なぜなら、二階微分位相像では位相飛びがスパイク状の異常値となって現れるため、本実施例のように簡易なアルゴリズムで修正できると共に、近傍ピクセルの値を調べる必要がないため積分演算の中に演算子Ｗとして組み込むことも容易だからである。
【００５８】
なお、上記実施形態および実施例は本発明の一具体例を示したものにすぎず、本発明の
範囲を限定する趣旨のものではない。例えば、上記実施例では、直交する２方向について演算を行ったが、第１の方向及び第２の方向については互いに交差する限りにおいて任意に設定することができる。また上記実施例ではトールボット干渉計を例示したが、本発明は他のシアリング干渉計にも好適に適用することができる。
【符号の説明】
【００５９】
１：撮像装置、１０：干渉計、１１：情報処理装置、３２：微分位相データ抽出部、３３：二階微分位相データ算出部、３４：位相データ算出部、２１０：被検体、８１０：干渉パターン
【技術分野】
【０００１】
本発明はシアリング干渉計により得られる干渉パターンから位相像や微分位相像のデータを生成するための技術に関する。
【背景技術】
【０００２】
シアリング干渉計とは、同一光源から照射された光線などのコヒーレントな光を分割し、一方に被検体による波面の歪曲を生成し、もう一方をわずかにずらすことで干渉縞を形成し、この干渉縞のずれから被検体による光の位相変化を可視化する手段である。このようなシアリング干渉計は、例えば位相差顕微鏡などの分野で被検体の視認性を向上させる目的などに利用されている。
【０００３】
シアリング干渉計の一つとして、トールボット干渉法を利用したものが知られている。トールボット干渉法では、光源からの光を被検体に照射し、被検体を透過することで位相が変化した光を回折格子に導く。回折格子により回折された光は、回折格子からトールボット距離と呼ばれる所定の距離だけ離れた位置に干渉パターンを形成する。Ｘ線を用いた場合、この干渉パターンの周期は非常に小さいため、干渉パターンを直接検出することが難しいことがある。そこで、干渉パターンが形成される位置に吸収格子を配置し、吸収格子によって干渉パターンの一部を遮ることでモアレを形成し、このモアレ像を検出器で検出する方法が提案されている。このモアレ像を解析することにより、被検体の位相の情報を得ることができる。このように光の位相を利用したイメージングは一般に位相イメージング（Ｘ線を用いる場合はＸ線位相イメージング）と呼ばれている。
【０００４】
ただし、干渉パターンによって直接解析できるのは、あくまで被検体を透過することによって変化した光の位相の空間的な微分量（微分位相像、シア像とも呼ばれる）である。空間的な微分はエッジを強調する効果があるため、微分位相像は光学顕微鏡などの分野では被検体像を強調するために用いられる。また、この微分位相像から元の位相を計算し、位相像を得ることも可能である。
【０００５】
非特許文献１には、二次元トールボット干渉法を用いたＸ線位相イメージング手法が開示されている。具体的には、二次元の回折格子により形成された干渉パターンからＸ方向の微分位相像とＹ方向の微分位相像を取得し、それらに対しフーリエ積分を行うことによって、位相像を計算している。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】C. Kottler, C. David, F. Pfeiffer and O. Bunk, “A two-directional approach for grating based differential phase contrast imaging using hard x-rays,” Opt. Express 15 (2007) 1175-1181.
【発明の概要】
【発明が解決しようとする課題】
【０００７】
シアリング干渉計を利用した位相イメージングは、画像診断や各種の非破壊検査への応用が可能である。中でもＸ線位相イメージングは、ＣＴ（コンピュータ・トモグラフィ）などＸ線吸収量を画像化する従来の手法に比べて被ばく量を大幅に低減できる、筋肉や内蔵などの軟組織も画像化できるといった利点から、新たな画像診断手法として期待されている。
【０００８】
この種の位相イメージングにおいては、ノイズやアーチファクトの低減による画質向上が基本的な課題として存在する。特に、光源としてＸ線を用いる場合には、量子ノイズが不可避であり、しかも被ばく量低減のために線量を下げるほどノイズが増える傾向にあるため、ノイズ対策は重要となる。また、Ｘ線ではなく可視光を用いる装置であっても、検査対象によっては（例えば眼球）低い光量しか用いることができない場合があり、ＳＮ比の悪化が問題となる。
【０００９】
また、位相イメージングにおいては「位相折りたたみ（位相ラッピング）」という問題があることが知られている。位相折りたたみは、干渉パターンから得られる微分位相の値が２πの範囲に折りたたまれて（ラップされて）しまう現象である。例えば−πからπの間で折りたたまれる場合、実際には３／２πという位相が正しい値であったとしても、それは−１／２πとして測定される。非特許文献１においてはこのような位相折りたたみの問題は無視されている。そのため、非特許文献１の手法で生成した位相像には、本来何もない、すなわち位相が一定になる空間に、位相のゆがみやアーチファクトが出現する場合がある。
【００１０】
位相折りたたみへの対策として、「位相接続（位相のアンラッピング）」と呼ばれる処理が知られている。位相接続とは、位相が折りたたまれた箇所を検出し、位相を元に戻す操作である。位相接続に関しては様々な手法が提案されている。最も基本的な操作は隣接ピクセルとの位相の差を計算する手法である。その差の絶対値が一定以上（例えば３／２π以上）であるならば、そこはもともとスムーズにつながっていた位相が折りたたまれたものと疑うことができ、位相を元に戻す操作をする。具体的には位相飛びが生じているピクセル以降のピクセル全ての位相に全て２πか、−２πを足すことで位相隣接ピクセルとの位相差を３／２π以下にする操作を行う。従来の位相接続は、積分とは独立した操作として、微分位相像全体に対して適用する必要があるため、位相像を得るまでに要する処理時間が長くなるという問題があった。
【００１１】
本発明は上記実情に鑑みてなされたものであり、位相イメージングにおいて、ノイズが低減された高品位の位相像あるいは微分位相像を生成可能な技術を提供することを目的とする。
【００１２】
本発明のさらなる目的は、位相イメージングにおいて、位相像あるいは微分位相像におけるノイズの低減、並びに、位相折りたたみによる位相のゆがみやアーチファクトの抑制を、簡易な処理で実現するための技術を提供することにある。
【課題を解決するための手段】
【００１３】
本発明の第１態様は、シアリング干渉計により得られる干渉パターンに基づいて位相像のデータを生成する撮像装置であって、被検体を透過した電磁波により形成された干渉パターンのデータから、第１の方向に関する位相の変化を表す第１の微分位相データと、前記第１の方向に交差する第２の方向に関する位相の変化を表す第２の微分位相データとを抽出する微分位相データ抽出手段と、抽出された前記第１の微分位相データを前記第１の方向に関して微分することによって第１の二階微分位相データを算出すると共に、抽出された前記第２の微分位相データを前記第２の方向に関して微分することによって第２の二階微分位相データを算出する二階微分位相データ算出手段と、前記第１および第２の二階微分位相データを関数として含む二階微分方程式を解くことによって、前記被検体の位相像のデータを算出する位相データ算出手段と、を有する撮像装置である。
【００１４】
本発明の第２態様は、シアリング干渉計により得られる干渉パターンに基づいて位相像のデータを生成する画像処理方法であって、情報処理装置が、被検体を透過した電磁波により形成された干渉パターンのデータから、第１の方向に関する位相の変化を表す第１の
微分位相データと、前記第１の方向に交差する第２の方向に関する位相の変化を表す第２の微分位相データとを抽出する工程と、情報処理装置が、抽出された前記第１の微分位相データを前記第１の方向に関して微分することによって第１の二階微分位相データを算出すると共に、抽出された前記第２の微分位相データを前記第２の方向に関して微分することによって第２の二階微分位相データを算出する工程と、情報処理装置が、前記第１および第２の二階微分位相データを関数として含む二階微分方程式を解くことによって、前記被検体の位相像のデータを算出する工程と、を有する画像処理方法である。
【００１５】
本発明の第３態様は、上記画像処理方法の各工程を情報処理装置に実行させるプログラムである。
【発明の効果】
【００１６】
本発明によれば、位相イメージングにおいて、ノイズが低減された高品位の位相像あるいは微分位相像を生成することができる。また、位相イメージングにおいて、位相像あるいは微分位相像におけるノイズの低減、並びに、位相折りたたみによる位相のゆがみやアーチファクトの抑制を、簡易な処理で実現することも可能となる。
【図面の簡単な説明】
【００１７】
【図１】撮像装置の構成を示す図である。
【図２】回折格子、干渉パターン、吸収格子の構成例を示す図である。
【図３】情報処理装置の機能構成を示す図である。
【図４】干渉パターンから微分位相像を抽出する手法を示す図である。
【図５】従来例における微分位相像から位相像を取得する処理を示す図である。
【図６】実施例における微分位相像から位相像を取得する処理を示す図である。
【図７】シミュレーションで用いた被検体と微分位相像を示す図である。
【図８】従来例の手法で得られた位相像を示す図である。
【図９】実施例の手法で得られた位相像を示す図である。
【図１０】実施例の手法で得られた位相像から生成した微分位相像を示す図である。
【発明を実施するための形態】
【００１８】
以下、本発明の好ましい実施の形態について添付の図面に基づいて説明する。なお、各図において、同一の部材については同一の参照番号を付し、重複する説明は省略する。
【００１９】
本実施形態では微分位相像を得る手段としてＸ線トールボット干渉計を使用した例について説明する。ただし、本発明が適用可能な微分位相像の取得手法はこれに限られるわけではなく、被検体（被検物）の微分位相像を得ることができれば、如何なる形態のシアリング干渉計を用いてもよい。また測定に用いる光はＸ線に限らず、如何なる波長の電磁波を用いてもよい。なお本明細書においてＸ線とはエネルギーが２以上２００ｋｅＶ以下の光を指す。
【００２０】
（撮像装置の構成）
図１は本実施形態における撮像装置の構成を示した図である。図１に示した撮像装置１は、概略、干渉計１０と情報処理装置１１とから構成される。干渉計１０は、２次元Ｘ線トールボット干渉計であり、Ｘ線を発生させるＸ線源（光源）１１０、Ｘ線を回折する回折格子３１０、Ｘ線の一部を遮る吸収格子４１０、Ｘ線を検出する検出器５１０を備えている。情報処理装置１１は、演算部（コンピュータ本体）６１０と画像表示装置７１０を備えている。情報処理装置１１は干渉計１０にケーブルで接続されており、干渉計１０の各部の制御を行う。また干渉計１０の検出器５１０で得られた干渉パターンのデータは情報処理装置１１に伝送され、演算部６１０において後述する各種の演算処理に供される。演算部６１０による演算結果に基づいた画像は画像表示装置７１０に出力される。
【００２１】
測定対象である被検体２１０は、Ｘ線源１１０と回折格子３１０の間に配置される。Ｘ線源１１０から照射され被検体２１０を透過したＸ線は回折格子３１０により回折され、トールボット距離と呼ばれる所定の距離Ｌの面上に明部と暗部が配列方向に並んだ干渉パターン８１０を形成する。本明細書では、Ｘ線（光）の強度が大きい所を明部、小さい所を暗部とする。本実施形態に用いられる回折格子３１０は位相型の回折格子である。回折格子として振幅型の回折格子を用いることもできるが、位相型の回折格子の方がＸ線量（光量）の損失が少ないので有利である。
【００２２】
図２（ａ）は位相型の回折格子３１０（以下単に位相格子３１０と呼ぶ）の平面上の構成の一例を表した図である。３１１は位相の基準となる領域、３１２は基準の領域３１１に対してπほどＸ線の位相が変化する領域である。図２（ｂ）は位相格子３１０によってトールボット距離に生じる干渉パターン８１０の明部と暗部を表している。これらの組み合わせには何種類かあるが、本発明では公知のものを含めて如何なる構成のものを用いてもよいため、ここでは詳しい説明を割愛する。
【００２３】
この干渉パターン８１０から微分位相像を抽出する方法には大きく分けて二種類考えられる。一つは、吸収格子４１０を少しずつずらしながら複数回撮像し、その結果から微分位相像を計算する縞走査法である。もう一つは、自己像としての干渉パターン８１０の周期と少しずらした吸収格子４１０を利用して撮像する方法で、フーリエ変換法と呼ばれる。以下ではフーリエ変換法に関して説明するが、どちらの手法でも微分位相像を取得できる。なお吸収格子は、遮蔽格子または振幅格子とも呼ばれることもある。
【００２４】
フーリエ変換法では図２（ｃ）で示すような吸収格子４１０を用いる。吸収格子４１０には、干渉パターン８１０の明部に対応する透過部４１１と、暗部に対応する遮蔽部４１２とが設けられている。ただし、透過部４１１の周期（ピッチ）は干渉パターン８１０の明部の周期とは異なる。このため、吸収格子４１０を透過した光像を実際に撮像すると、干渉パターン８１０と吸収格子４１０の周期の違いによるモアレが観測される。
【００２５】
通常、位相格子３１０による干渉パターン８１０の周期は数μｍから十数μｍ程度であるが、吸収格子４１０によりモアレを作ることにより、その周期を数十μｍ以上に拡大することができる。これによって数十μｍ程度の分解能の検出器５１０で干渉パターンを撮像することができる。しかしながら、例えば、検出器５１０の分解能が十分高い場合や、位相格子３１０と検出器５１０の距離が大きい場合など、検出器５１０で干渉パターン８１０そのものを解像可能な場合には、吸収格子４１０を設けなくてもよい。
【００２６】
（情報処理装置の機能）
図３は、情報処理装置１１が有する機能を示すブロック図である。情報処理装置１１は、概略、干渉計１０の制御（例えばＸ線の照射、検出器５１０からのデータ取り込みなど）を行う制御部３０と、干渉計１０から取り込まれたデータに基づいて位相像のデータを生成するデータ処理部３１とを有する。データ処理部３１は、微分位相データ抽出部３２、二階微分位相データ算出部３３、及び、位相データ算出部３４から構成されている。微分位相データ抽出部３２は、干渉パターン（モアレ像）のデータを数値的に解析して微分位相データを抽出する機能であり、二階微分位相データ算出部３３は、微分位相データを微分することによって二階微分位相データを算出する機能である。位相データ算出部３４は、二階微分位相データを関数として含む二階微分方程式を解くことによって、被検体の位相像のデータを算出する機能である。これらの機能の詳細については後述する。
【００２７】
情報処理装置１１の演算部６１０は、ＣＰＵ（中央演算処理装置）、主記憶装置（ＲＡＭなど）、補助記憶装置（ＨＤＤ、ＳＳＤなど）、各種Ｉ／Ｆを有するコンピュータで構
成することができる。図３に示した機能は、補助記憶装置に格納されたプログラムが主記憶装置にロードされ、ＣＰＵにより実行されることで実現されるものである。もちろん、この構成はあくまで一例であり、本発明の範囲を限定するものではない。
【００２８】
（撮像装置の動作）
以下、本実施形態の撮像装置１において実行される画像処理方法、並びに従来例との差異を説明する。ただしここでは、干渉パターンのデータとして、干渉計１０により実際に被検物を計測して得たデータではなく、コンピュータシミュレーションによって生成したデータを用いる。干渉パターンのデータから位相像のデータを生成する一連の処理については、情報処理装置１１のデータ処理部３１によって行った。
【００２９】
シミュレーションでは、Ｘ線として１７．５ｋｅＶの単色光を仮定し、位相格子として図２（ａ）に示す構造の位相格子３１０を仮定した。位相格子３１０には、Ｘ線を十分透過させる基盤部３１１に対して柱状の構造物３１２がピッチｐ１で並んでいると仮定する。構造物３１２を透過したＸ線は、基盤部３１１を透過したＸ線に対して、位相がπ変化するようにした。このような構造の位相格子３１０は例えば構造物３１２がシリコンで作成された場合、２２．６μｍ程度の奥行きにすることによって実現される。
【００３０】
さらに図２（ａ）においてピッチｐ１を８μｍと仮定すると図１における位相格子３１０と吸収格子４１０の距離Ｌを約１１．３ｃｍに設定することにより、吸収格子４１０の直前で干渉パターン８１０が形成される。この干渉パターン８１０の周期ｐ２は４μｍである。吸収格子４１０によりモアレを形成することで、干渉パターン８１０を拡大し、検出器５１０でその強度を検出する。簡単のため、検出器５１０に含まれる画素素子の大きさは４μｍとして計算した。
【００３１】
トールボットのような干渉計の場合、入射光には高いコヒーレンス性が必要とされる。そのような高いコヒーレンス性を実験室や医療装置のような大きさがせいぜい数ｍ程度までの施設で用いる場合、微小な点光源を用いてコヒーレンス性を出す手法が用いられる。このような場合、入射光は平行波ではなく球面波になる。しかし、この違いは本発明の本質を変えるものではないため、本発明は平行光でも非平行光（球面波など）でも適用可能である。以下の実施例では平行光を仮定してシミュレーションを行った。
【００３２】
被検体としては二種類用意した。一つ目は図７（ａ）に示すような直径４００μｍのカルシウム球を想定した。もう一つは図７（ｂ）に示すような直径２００μｍのカルシウム球が４つ重なっている被検体である。計算の範囲は５１２μｍ×５１２μｍの範囲となる。
【００３３】
微分位相像の抽出方法としては前述の通りフーリエ変換法を使用した。そのために、吸収格子４１０の透過部４１１の周期ｐ３（図２（ｃ）参照）は４．５μｍとした。これによりほぼ４画素で一周期のモアレが実現することとなる。
【００３４】
図４（ａ）から図４（ｄ）はフーリエ変換法により干渉パターンから微分位相像を抽出する方法、すなわち微分位相データ抽出部３２の処理内容を概念的に示した図である。図２（ｂ）のような干渉パターン８１０を拡大した結果、図４（ａ）のようなモアレ像が取得できたとする。フーリエ変換法ではこのモアレ像をフーリエ変換する。これによってフーリエ空間内では９つのスペクトルに分離される。この９つのスペクトルのうち１次スペクトルにあたるスペクトル９１１、９１２、９１３、９１４のいずれかの周辺を窓関数で切り取って原点に戻し（図４（ｃ））、これを逆フーリエ変換することで位相微分像を得ることができる。図４（ｄ）の左側は、スペクトル９１４から生成された微分位相像の例であり、ｘ方向（第１の方向）に関する位相の変化を表している。図４（ｄ）の右側は、
スペクトル９１１から生成された微分位相像の例であり、ｙ方向（第２の方向）に関する位相の変化を表している。
【００３５】
図７（ｃ）は、図７（ａ）の被検体から得られたｘ方向、ｙ方向それぞれの微分位相像を示す。また図７（ｄ）は、図７（ｂ）の被検体から得られたｘ方向、ｙ方向それぞれの微分位相像を示す。球体のエッジ部分は最も急激に位相の微分が変化する個所となるため、この部分で位相飛びが発生しやすい。また、フーリエ変換法では、前述の窓関数の取り方によっては微分位相像にアーチファクト（本来は存在しない像）が発生する場合がある。図７（ｃ）、および図７（ｄ）の微分位相像でも若干のアーチファクトが見られる。
【００３６】
二次元の場合はｘ方向とｙ方向の微分位相像がそれぞれ独立に取得できる。この方向の違いは、図４（ｂ）のフーリエ空間においてどのスペクトルを切り取るかによって制御できる。例えばスペクトル９１１又は９１２を切り取ればｙ方向の微分位相像が得られ、スペクトル９１３又は９１４を切り取ればｘ方向の微分位相像が得られる。もちろんｘ、ｙの直交する２方向ではなく、互いに交差する任意の２方向についての微分位相像を得ることも可能である。以下では取得されたｘ方向微分位相像をφｘ（ｘ，ｙ）、ｙ方向微分位相像をφｙ（ｘ，ｙ）とおく。（ｘ，ｙ）は座標である。これらのｘ方向微分位相像とｙ方向微分位相像から、積分された位相像を求める。
【００３７】
（１）従来の位相像取得方法
従来例では、ｘ方向微分位相像とｙ方向微分位相像を別々に積分するか、フーリエ変換などの手段を用いて積分する。以下ではより誤差の少ない位相が求められるフーリエ変換による積分法（フーリエ積分法）を従来例として述べる。
【００３８】
図５はフーリエ積分法による位相像取得方法の手順を示すフローチャートである。まず前述のように、ｘ方向微分位相像φｘ（ｘ，ｙ）とｙ方向微分位相像φｙ（ｘ，ｙ）を取得する（ステップＳ５０）。ｇ（ｘ，ｙ）を次の式で定義する（ステップＳ５１）。
【数１】
ここで、ｉは虚数単位を示す。
【００３９】
ｇ（ｘ，ｙ）をフーリエ変換することでＧ（ｋｘ，ｋｙ）を得て（ステップＳ５２）、次の（式３）によって位相Ｐ（ｘ，ｙ）を取得する（ステップＳ５３、Ｓ５４）。
【数２】
ここで、Ｆ［・・・］は括弧の中をフーリエ変換することを示す演算子であり、Ｆ−１［・・・］は逆フーリエ変換を示す演算子である。ｋｘとｋｙはそれぞれｘ方向、ｙ方向の波数を表す。（式３）はすなわち（式１）で定義されたｇをフーリエ変換し積分演算子１／（２πｉ×（ｋｘ＋ｋｙ））を作用させて逆フーリエ変換をかけることと一致する。こうして得たＰ（ｘ，ｙ）は被検体の位相像を表す（ステップＳ５５）。
【００４０】
図８（ａ）は、図７（ｃ）で示した微分位相像より、図５に示した従来例による方法で取得した位相像の例である。カルシウム球の位相が球形に回復されているのがわかる。しかし、位相の折りたたみにより被検体の外周部で位相が不連続になっており、取得した位相像にゆがみが生じている。特に何もない真空部（一定）であるはずの被検体の外周部に
、本来ありえない位相のゆがみが生じていることがわかる。
【００４１】
図８（ｂ）は、図７（ｄ）で示した微分位相像より、図５に示した従来例による方法で取得した位相像の図である。４つのカルシウム球の位相が復元されていることがわかる。しかし、先の例と同様に本来は真空であるはずの領域に位相のゆがみができている。また、アーチファクトの影響で本来は球形である位相がジャギーを持っているように見える。
【００４２】
（２）実施例の位相像取得方法
図６は、本発明の実施例における位相像取得方法の手順を示すフローチャートである。従来例では、微分位相像を直接積分することで位相像を得ていたのに対し、本実施例では、微分位相像をｘ，ｙそれぞれの方向について微分し二階微分位相像を生成した後、二階微分位相像からなる二階微分方程式を解いて位相像を得る点に特徴がある。
【００４３】
まず、微分位相データ抽出部３２が、干渉パターンのデータ（モアレ像のデータ）から、ｘ方向微分位相像φｘ（ｘ，ｙ）とｙ方向微分位相像φｙ（ｘ，ｙ）を取得する（ステップＳ６０）。微分位相データ抽出部３２が生成したｘ方向微分位相像のデータ（第１の微分位相データ）およびｙ方向微分位相像のデータ（第２の微分位相データ）は、二階微分位相データ算出部３３に引き渡される。
【００４４】
二階微分位相データ算出部３３は、ｘ方向微分位相像をｘ方向に関して微分することによってｘ方向二階微分位相像を算出すると共に、ｙ方向微分位相像をｙ方向に関して微分することによってｙ方向二階微分位相像を算出する（ステップＳ６１）。ｘ方向二階微分位相像ψｘ（ｘ，ｙ）、ｙ方向二階微分位相像ψｙ（ｘ，ｙ）は、下記式で表される。
【数３】
【００４５】
ｘ方向二階微分位相像のデータ（第１の二階微分位相データ）およびｙ方向二階微分位相像のデータ（第２の二階微分位相データ）は、位相データ算出部３４に引き渡される。
【００４６】
位相データ算出部３４は、既知関数としてｘ方向、ｙ方向それぞれの二階微分位相像を含む二階微分方程式を解くことで、未知関数であるｘ方向、ｙ方向それぞれの位相像を求める（ステップＳ６２〜Ｓ６６）。そのアルゴリズムについて以下詳しく述べる。
【００４７】
まず位相データ算出部３４は、ｘ方向とｙ方向の二次微分位相像を用いて、ρ（ｘ，ｙ）を次のように定義する（ステップＳ６２）。
【数４】
ここで、Ｗ［・・・］は、［・・・］の値のうち、所定の範囲を外れた値をその所定の範囲内に収める処理を行う演算子である。微分位相像φｘ、φｙの値は基本的に連続的に変化するため、二階微分位相像ψｘ、ψｙの値は０を中心とする狭い範囲（せいぜい−π〜πの範囲）の内に収まる。しかし、位相の折りたたみが生じると、微分位相像φｘ、φｙの値が−πからπ、又はπから−πへと不連続に変化するため、二階微分位相像ψｘ、ψｙの値が２π又は−２πとなる。すなわち、位相の折りたたみが生じた箇所は、二階微分位相像ψｘ、ψｙにおいてスパイク状の異常値となって現れる。Ｗ［・・・］の操作は
、この異常値を除去（修正）する操作、すなわち、位相の折りたたみが生じた状態を元の位相の状態に戻す操作（位相接続）に相当する。以下、Ｗ［・・・］をアンラッピング演算子と呼ぶ。
【００４８】
アンラッピング演算子にはさまざまなデザインが考えられるが、本実施例では、具体的かつもっとも簡単な演算子として、次のルールを持つ演算子Ｗ［・・・」を用いる。
（ａ）［・・・］の値がπ以上のものに関して、その値から２πを引く
（ｂ）［・・・］の値が−π以下のものに関して、その値に２πを足す
以上のルールを適用することで元関数φｘ（ｘ，ｙ）やφｙ（ｘ，ｙ）に位相折りたたみが生じていても、位相接続がなされたのと同じ効果を持つ。
【００４９】
実際の位相接続では、上記のみならずノイズによる特異点の処理などを含めた様々な手法が存在する。Ｗ［・・・］もそれらの処理方法を応用した様々な演算アルゴリズムが考えられるが、本発明はそれらの手法に本質的には依存せず様々な方法を適用可能である。
【００５０】
位相データ算出部３４は、（式６）のρ（ｘ，ｙ）から次のポワソン方程式を解くことによって位相Ｐ（ｘ，ｙ）を計算する。
【数５】
【００５１】
このポワソン方程式を解く具体的アルゴリズムも様々な手法が考えられるがその一つとしては、フーリエ変換による二階の積分法が考えられる。すなわち、位相データ算出部３４は、ρ（ｘ，ｙ）をフーリエ変換することでＨ（ｋｘ，ｋｙ）を得て（ステップＳ６３）、次の（式９）によって位相Ｐ（ｘ，ｙ）を算出する（ステップＳ６４、Ｓ６５）。
【数６】
（式９）では（式６）で定義されたρ（ｘ，ｙ）に対してフーリエ変換を求め、１／（４π２（ｋｘ２＋ｋｙ２））の二階の積分演算子を作用させた後に逆フーリエ変換して、取得した値を位相Ｐ（ｘ，ｙ）とする。以上の処理により、二階微分位相像の二階積分である位相像を得ることができる。
【００５２】
図９（ａ）は、図７（ｃ）で示した微分位相像より、図６に示した本実施例の手法によって取得した位相像の例である。従来例の図８（ａ）と比べると、位相折りたたみによる歪みが改善されていることがわかる。これは従来例では位相が不連続になっていた箇所が、ある程度修復された影響によるものであり、本実施例の手法の有効性を示している。
【００５３】
図９（ｂ）は、図７（ｄ）で示した微分位相像より、図６に示した本実施例の手法によって取得した位相像の例である。先の例と同様に、従来例の図８（ｂ）と比べると位相折りたたみによる歪みが改善されており、本実施例の手法の有効性を示している。また、従来例の図８（ｂ）では目立っているアーチファクトについても、図９（ｂ）では低減していることがわかる。
【００５４】
また本実施例の手法はノイズ耐性に関しても優れている。表１は、従来例の位相像取得方法と本実施例の位相像取得方法とで、同じノイズに対するノイズリダクションの機能を
比べた結果である。入力画像として平均０、分散が０．５のランダムパターンの画像２枚を微分位相像とみなして位相を計算し、その位相像にどれだけのノイズが存在するかを比較した表である。この実験では、純粋にノイズリダクションを評価するため、被検体を含まない画像を用いている。表１の数値は分散であり、値が小さいほどノイズが少ない（ノイズが低減されている）ことを示す。１０回の計算を行なった結果、すべての回において本実施例の手法でノイズが減少していた。
【表１】
【００５５】
このように本実施例の手法は従来の手法よりもノイズ低減の効果もあることが分かる。この実験では被検体を含まない画像を用いたので、このノイズ低減の効果は、前述したアンラッピング演算子Ｗによる位相接続による画像の改善とは、本質的に別の効果である。すなわち位相接続を行わなくとも（つまり（式６）において演算子Ｗを省略しても）、本実施例の手法が従来の手法に比べて有利な効果を有していることがわかる。
【００５６】
図１０（ａ）、図１０（ｂ）はノイズ低減効果の応用例を示している。図１０（ａ）は、図７（ａ）の被検体の像に対し分散０．２５のノイズを重畳した場合のｘ方向及びｙ方向の微分位相像である。図１０（ｂ）は、図１０（ａ）の微分位相像に対し本実施例の位相像取得方法を適用した後、位相像をｘ方向、ｙ方向それぞれに関し微分することで生成した微分位相像である。図１０（ａ）の元の微分位相像においてみられたアーチファクトやノイズが、図１０（ｂ）では改善されており、高品位な微分位相像が得られている。すなわち本実施例の手法は、微分位相像の画質を向上するという目的にも利用できる。画像診断や検査の内容によっては、位相像をみるよりも、エッジが強調された微分位相像のほうが有用な場合もある。よって、高画質の微分位相像を生成できる機能も実用上非常に有利である。構成としては、例えば、図３の位相データ算出部３４の後段に、位相像のデータを微分する微分位相データ算出部を設ければよい。
【００５７】
以上述べたように、本実施例の位相像取得方法によれば、位相イメージングにおいて、ノイズが低減された高品位の位相像および微分位相像を生成することが可能である。そして、位相接続処理を施したことにより、位相折りたたみに起因する位相のゆがみやアーチファクトを改善することもできる。しかも本実施例では、二階微分位相像に対して位相接続処理を施す構成としたので、従来（微分位相像に対して位相接続処理を施す）に比べて、位相接続処理が格段に簡単となり、また高速に処理できる。なぜなら、二階微分位相像では位相飛びがスパイク状の異常値となって現れるため、本実施例のように簡易なアルゴリズムで修正できると共に、近傍ピクセルの値を調べる必要がないため積分演算の中に演算子Ｗとして組み込むことも容易だからである。
【００５８】
なお、上記実施形態および実施例は本発明の一具体例を示したものにすぎず、本発明の
範囲を限定する趣旨のものではない。例えば、上記実施例では、直交する２方向について演算を行ったが、第１の方向及び第２の方向については互いに交差する限りにおいて任意に設定することができる。また上記実施例ではトールボット干渉計を例示したが、本発明は他のシアリング干渉計にも好適に適用することができる。
【符号の説明】
【００５９】
１：撮像装置、１０：干渉計、１１：情報処理装置、３２：微分位相データ抽出部、３３：二階微分位相データ算出部、３４：位相データ算出部、２１０：被検体、８１０：干渉パターン
【特許請求の範囲】
【請求項１】
シアリング干渉計により得られる干渉パターンに基づいて位相像のデータを生成する撮像装置であって、
被検体を透過した電磁波により形成された干渉パターンのデータから、第１の方向に関する位相の変化を表す第１の微分位相データと、前記第１の方向に交差する第２の方向に関する位相の変化を表す第２の微分位相データとを抽出する微分位相データ抽出手段と、
抽出された前記第１の微分位相データを前記第１の方向に関して微分することによって第１の二階微分位相データを算出すると共に、抽出された前記第２の微分位相データを前記第２の方向に関して微分することによって第２の二階微分位相データを算出する二階微分位相データ算出手段と、
前記第１および第２の二階微分位相データを関数として含む二階微分方程式を解くことによって、前記被検体の位相像のデータを算出する位相データ算出手段と、
を有することを特徴とする撮像装置。
【請求項２】
前記位相データ算出手段は、位相像のデータを算出する際に、前記第１および第２の二階微分位相データに対して位相接続処理を施すことを特徴とする請求項１に記載の撮像装置。
【請求項３】
前記位相接続処理は、前記第１および第２の二階微分位相データの値のうち所定の範囲を外れた値を、前記所定の範囲内に収める処理であることを特徴とする請求項２に記載の撮像装置。
【請求項４】
前記二階微分方程式は、既知関数である前記第１の二階微分位相データと前記第２の二階微分位相データとの和が、未知関数である前記被検体の位相像のデータを前記第１の方向に二階微分したものと前記第２の方向に二階微分したものとの和に等しいとする方程式であることを特徴とする請求項１〜３のうちいずれか１項に記載の撮像装置。
【請求項５】
前記位相データ算出手段によって算出された前記被検体の位相像のデータを微分することによって、前記被検体の微分位相像のデータを生成する微分位相データ算出手段をさらに有することを特徴とする請求項１〜４のうちいずれか１項に記載の撮像装置。
【請求項６】
シアリング干渉計により得られる干渉パターンに基づいて位相像のデータを生成する画像処理方法であって、
情報処理装置が、被検体を透過した電磁波により形成された干渉パターンのデータから、第１の方向に関する位相の変化を表す第１の微分位相データと、前記第１の方向に交差する第２の方向に関する位相の変化を表す第２の微分位相データとを抽出する工程と、
情報処理装置が、抽出された前記第１の微分位相データを前記第１の方向に関して微分することによって第１の二階微分位相データを算出すると共に、抽出された前記第２の微分位相データを前記第２の方向に関して微分することによって第２の二階微分位相データを算出する工程と、
情報処理装置が、前記第１および第２の二階微分位相データを関数として含む二階微分方程式を解くことによって、前記被検体の位相像のデータを算出する工程と、
を有することを特徴とする画像処理方法。
【請求項７】
請求項６に記載の画像処理方法の各工程を情報処理装置に実行させることを特徴とするプログラム。
【請求項１】
シアリング干渉計により得られる干渉パターンに基づいて位相像のデータを生成する撮像装置であって、
被検体を透過した電磁波により形成された干渉パターンのデータから、第１の方向に関する位相の変化を表す第１の微分位相データと、前記第１の方向に交差する第２の方向に関する位相の変化を表す第２の微分位相データとを抽出する微分位相データ抽出手段と、
抽出された前記第１の微分位相データを前記第１の方向に関して微分することによって第１の二階微分位相データを算出すると共に、抽出された前記第２の微分位相データを前記第２の方向に関して微分することによって第２の二階微分位相データを算出する二階微分位相データ算出手段と、
前記第１および第２の二階微分位相データを関数として含む二階微分方程式を解くことによって、前記被検体の位相像のデータを算出する位相データ算出手段と、
を有することを特徴とする撮像装置。
【請求項２】
前記位相データ算出手段は、位相像のデータを算出する際に、前記第１および第２の二階微分位相データに対して位相接続処理を施すことを特徴とする請求項１に記載の撮像装置。
【請求項３】
前記位相接続処理は、前記第１および第２の二階微分位相データの値のうち所定の範囲を外れた値を、前記所定の範囲内に収める処理であることを特徴とする請求項２に記載の撮像装置。
【請求項４】
前記二階微分方程式は、既知関数である前記第１の二階微分位相データと前記第２の二階微分位相データとの和が、未知関数である前記被検体の位相像のデータを前記第１の方向に二階微分したものと前記第２の方向に二階微分したものとの和に等しいとする方程式であることを特徴とする請求項１〜３のうちいずれか１項に記載の撮像装置。
【請求項５】
前記位相データ算出手段によって算出された前記被検体の位相像のデータを微分することによって、前記被検体の微分位相像のデータを生成する微分位相データ算出手段をさらに有することを特徴とする請求項１〜４のうちいずれか１項に記載の撮像装置。
【請求項６】
シアリング干渉計により得られる干渉パターンに基づいて位相像のデータを生成する画像処理方法であって、
情報処理装置が、被検体を透過した電磁波により形成された干渉パターンのデータから、第１の方向に関する位相の変化を表す第１の微分位相データと、前記第１の方向に交差する第２の方向に関する位相の変化を表す第２の微分位相データとを抽出する工程と、
情報処理装置が、抽出された前記第１の微分位相データを前記第１の方向に関して微分することによって第１の二階微分位相データを算出すると共に、抽出された前記第２の微分位相データを前記第２の方向に関して微分することによって第２の二階微分位相データを算出する工程と、
情報処理装置が、前記第１および第２の二階微分位相データを関数として含む二階微分方程式を解くことによって、前記被検体の位相像のデータを算出する工程と、
を有することを特徴とする画像処理方法。
【請求項７】
請求項６に記載の画像処理方法の各工程を情報処理装置に実行させることを特徴とするプログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【公開番号】特開２０１３−１０２９５１（Ｐ２０１３−１０２９５１Ａ）
【公開日】平成２５年５月３０日（２０１３．５．３０）
【国際特許分類】
【出願番号】特願２０１１−２４８７１７（Ｐ２０１１−２４８７１７）
【出願日】平成２３年１１月１４日（２０１１．１１．１４）
【出願人】（０００００１００７）キヤノン株式会社 (59,756)
【出願人】（５０４１３３１１０）国立大学法人電気通信大学 (383)
【Ｆターム（参考）】
【公開日】平成２５年５月３０日（２０１３．５．３０）
【国際特許分類】
【出願日】平成２３年１１月１４日（２０１１．１１．１４）
【出願人】（０００００１００７）キヤノン株式会社 (59,756)
【出願人】（５０４１３３１１０）国立大学法人電気通信大学 (383)
【Ｆターム（参考）】
[ Back to top ]