JP2016061608A - 画像処理方法、画像処理装置、撮像システム - Google Patents

Abstract

【課題】微分干渉計で得られた干渉像のデータから、簡単な処理で、位相の畳み込みの影響が少ない被検体画像を取得するための技術を提供する。
【解決手段】画像処理装置が、微分干渉計により得られた干渉像のデータを取得し、前記干渉像の波数空間上のスペクトルに対し、前記スペクトルの中心に関して反対称な関数もしくは該関数の項を含む多項式関数で表現されるフィルタをかけ、フィルタされた結果を用いて画像データを生成する。
【選択図】図８

Description

本発明は、微分干渉計により得られた干渉像の画像処理技術に関する。

従来、精密に物質を測定する手段として位相を用いた測定法が用いられてきた。位相を用いた測定法は、波面のそろった（コヒーレントな）入射光に対して干渉を起こし、その干渉縞を計測することで、波長の数分の一から数十分の一の位相差による入射光波面（位相）の変化を測定する。このような測定法を利用した干渉計は、例えばレンズの表面のわずかな凹凸を測定するのに好適な手段である。

さらに、干渉を用いた波面計測手法の中でも近年注目を集めているのが、数十ナノメール以下の波長の光（電磁波）、すなわちＸ線を利用したＸ線位相イメージングである。Ｘ線位相イメージングは、従来の被検体の吸収によるコントラストを画像化するＸ線吸収像とは異なり、被検体に対するＸ線の透過時に生じる位相の変化を検出する手法である。位相変化を検出することで従来の吸収像では検出が困難であった生体軟組織などの被検体内部における吸収係数の低い箇所を観察することができる。

Ｘ線位相イメージングの例として、Ｘ線を用いたトールボット干渉法について述べる（特許文献１）。Ｘ線トールボット干渉計では、光源からのＸ線が被検体を透過し、それに伴って光の入射位相が変化する。被検体を透過した光は、回折格子と呼ばれる周期的パターンを持った格子で回折されることによって、回折格子から所定の距離だけ離れた位置に第一の干渉パターンを形成する。この第一の干渉パターンの変化を解析することで前述の入射光波面の変化を測定する。

前述の周期的なパターンを持った回折格子のパターン周期は、装置長や入射光の波長等の条件によって変化する。一般的なＸ線の場合、その周期は数μｍのオーダーとなる。また、それによって生じる第一の干渉パターンも同様に数μｍのオーダーの周期となることが知られている。このような場合、一般的に用いられている検出器では、分解能がせいぜい数十μｍであるため、第一の干渉パターンを検出することは不可能である。そのため、干渉パターンが形成される位置に第一の干渉パターンとほぼ同じ周期の遮蔽格子を配置する。遮蔽格子で第一の干渉パターンの一部を遮ることにより周期が数百μｍ程度の第二の干渉パターン、すなわちモアレを形成し、このモアレを検出器で検出することによって干渉パターンの変化を間接的に測定することができる。モアレの形成法には周期を調節した遮蔽格子を第一の格子による第一の干渉パターンと同じ向きにアライメントして出す場合（拡大モアレ）と、格子を回転させてモアレの周期と向きを調節する場合（回転モアレ）がある。

トールボット干渉計は微分干渉計であり、干渉像（第一の干渉パターン又は第二の干渉パターン）から取得される一次情報は、隣接した二点の位相差、すなわち微分位相の情報である。前述の回折格子や遮蔽格子として、例えばストライプ形状の一次元の格子を用いた場合には、一軸方向に沿った方向の微分位相情報が得られる。一方、回折格子や遮蔽格子として、例えば市松格子状の二次元の格子を用いた場合には、例えばｘ軸方向とｙ軸方向の二次元の微分位相情報を得ることができる。位相の定量性を計算するにあたり、微分情報を積分する操作を通して元の被検体による入射光波面の変化、すなわち位相変化量を表す像を回復する必要がある。二次元格子を利用することで、少なくとも二つの軸に沿った二次元の微分情報を積分に使用することができるため、ノイズの少ない積分位相像が得られることが知られている。

二次元微分干渉計において、微分位相から位相を計算する手段としては、二次元の微分位相を複素空間上に投影し、フーリエ空間内で積分演算子となる作用素を作用させる手法が知られている。また、特許文献１では、微分位相をさらに微分して得た二次微分位相からポワソン方程式におけるソース項（特許文献１の式６におけるρ）を導きだし、二重積分作用素によって位相を求めることでより低ノイズの位相を求める手法について開示している。

これらの位相回復手法は、一旦微分位相情報を実空間上に投影する操作を必要とする。言い換えると、微分位相は複素関数の偏角という形で取得されるために、微分位相を求める際に複素関数の偏角を求める操作が必要になるということである。複素関数の偏角という形で数値を求める場合、位相の畳みこみ（ラッピング）という問題が生じる。これは実際に得られる偏角が、−πからπまでの２πの範囲内に限定され、この範囲外の数値はすべて、−πからπまでの数値に投影（畳み込み）されてしまい、例えば０と２πの区別がつかなくなることを示している。位相の畳み込みは取得される微分位相の定量性を悪化させるほか、画像に表示した場合にアーチファクトとなって定性的な印象にも影響を与える。

この問題に対処するため、正しい位相を求める位相接続（アンラッピング）に関する研究も広く行われており、様々な手法が提案されている。しかし、この位相接続処理アルゴリズムは数学的に複雑な計算を必要とし、いずれのアルゴリズムも処理時間がかかることが知られている。それゆえ、微分干渉計を利用したシステムでは、撮像（干渉像の取得）から微分位相像や位相像を得るまで、さらにはそれらの画像から特徴抽出を行うまでに、ある程度の処理時間を要してしまうという課題がある。特に、二次元微分干渉計のようにｘ軸方向の微分位相とｙ軸方向の微分位相を同時に求める必要がある場合、その処理コストは二倍となる。そのため、リアルタイム性が要求されるアプリケーション（例えば、医療や検査など）に本手法が適さないという問題があった。

特開２０１３−１０２９５１号公報

本発明は上記実情に鑑みなされたものであって、微分干渉計で得られた干渉像のデータから、簡単な処理で、位相の畳み込みの影響が少ない被検体画像を取得するための技術を提供することを目的とする。

本発明の第一態様は、コンピュータが、微分干渉計により得られた干渉像のデータを取得する取得ステップと、コンピュータが、前記干渉像の波数空間上のスペクトルに対し、前記スペクトルの中心に関して反対称な関数もしくは該関数の項を含む多項式関数で表現されるフィルタをかけるフィルタステップと、コンピュータが、フィルタされた結果を用いて画像データを生成する生成ステップと、を含むことを特徴とする画像処理方法である。

本発明の第二態様は、本発明に係る画像処理方法の各ステップをコンピュータに実行させることを特徴とするプログラムである。

本発明の第三態様は、微分干渉計により得られた干渉像のデータを取得する取得手段と
、前記干渉像の波数空間上のスペクトルに対し、前記スペクトルの中心に関して反対称な関数もしくは該関数の項を含む多項式関数で表現されるフィルタをかけるフィルタ手段と、フィルタされた結果を用いて画像データを生成する生成手段と、を有することを特徴とする画像処理装置である。

本発明の第四態様は、微分干渉計からなる撮像装置と、前記撮像装置により得られた干渉像のデータに基づき画像データを生成する、本発明に係る画像処理装置と、を有することを特徴とする撮像システムである。

本発明によれば、微分干渉計で得られた干渉像のデータから、簡単な処理で、位相の畳み込みの影響が少ない被検体画像を取得することができる。

Ｘ線位相イメージングの撮像システムの構成を模式的に示す図。 被検体、フリンジ像、ｘ軸方向とｙ軸方向の微分位相像の例を示す図。 図２の微分位相像を直接従来方法に適用した例を示す図。 本発明の実施例１に関する図。 本発明の実施例２に関する図。 本発明の実施例３に関する図。 本発明の実施例４に関する図。 本発明の実施形態に係る画像処理方法の流れを示す図。

以下、本発明の好ましい実施の形態について添付の図面に基づいて説明する。本発明は、微分干渉計により得られた干渉像のデータから、簡単な処理で（従来のような位相接続の処理を行わずに）、位相の畳み込みの影響が少ない被検体画像を取得する技術に関する。以下に述べる実施形態では、二次元微分干渉計の一つである二次元トールボットＸ線位相イメージング装置を例に挙げて説明する。ただし、本発明の適用範囲にこれに限られず、干渉像に含まれる位相変化から被検体の情報を取得する装置全般に本発明を好ましく適用することができる。もちろん二次元でなく一次元の微分干渉計にも適用可能である。また測定に用いる光はＸ線に限らず、如何なる波長の電磁波を用いてもよい。なお本明細書においてＸ線とはエネルギーが２以上２００ｋｅＶ以下の光を指す。

図１はトールボットＸ線位相イメージングの撮像システムの構成を示した図である。撮像システムは、Ｘ線トールボット干渉計からなる撮像装置１０と、撮像装置１０で取得したＸ線画像を処理する画像処理装置１１とから構成される。撮像装置１０は、Ｘ線を発生させるＸ線源（光源）１１０と、Ｘ線を回折する回折格子１３０と、Ｘ線の一部を遮る遮蔽格子１４０と、Ｘ線を検出する検出器１５０とを備えている。画像処理装置１１は、撮像装置１０に接続された演算部（コンピュータ）１６０と、演算部１６０による演算結果に基づいた画像を表示する画像表示装置１７０とを有している。

演算部１６０は、ＣＰＵ（中央演算処理装置）、ＲＡＭ（ランダムアクセスメモリ）、補助記憶装置などのハードウェア資源を備えた汎用のコンピュータにより構成できる。後述する画像処理、各種演算、および制御は、補助記憶装置に格納されたプログラムをＣＰＵが読み込み実行することで実現されるものである。なお、演算部１６０の機能のうちの一部又は全部をＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）のような回路で構成することもできる。

Ｘ線源１１０からのＸ線は回折格子１３０により回折され、トールボット距離と呼ばれ
る所定の距離をおいて明部と暗部が配列方向に並んだ干渉パターン１８０を形成する。

通常、回折格子１３０による第一の干渉パターン１８０の周期は数μｍから十数μｍ程度である。そこで、第一の干渉パターン１８０と同じかわずかに異なる周期をもつ遮蔽格子１４０を、第一の干渉パターン１８０が形成される位置に配置する。そうすると、第一の干渉パターン１８０と遮蔽格子１４０によりモアレが形成され、干渉パターンの周期を数十μｍ以上あるいは無限（遮蔽格子上に形成される干渉パターンと遮蔽格子とのピッチとは配列方向とが一致する場合）に拡大することができる。モアレの周期は、用いる位相回復方法と検出器の画素サイズを考慮して適宜決めることができるが、本実施形態においては、モアレの周期が画素サイズの２倍以上、検出器の撮像範囲以下であることが好ましい。このモアレ（空間的な周期性をもつパターン）を二次元イメージセンサである検出器１５０により画像化し、二次元画像を得る。このような仕組みにより、数十μｍ平方程度の分解能の検出器１５０で、数μｍから十数μｍの周期をもつ干渉パターンのイメージングを可能にしている。ただし、検出器１５０の空間分解能が十分に高い場合には、遮蔽格子１４０を省略し、第一の干渉パターン１８０をそのまま撮像してもよい。以下、検出器１５０で得られた周期的なパターンをもつ画像を、干渉像又はフリンジ像とも呼ぶ。

回折格子１３０には、周期的にＸ線の位相を変調する位相型の回折格子（位相格子）でも、周期的にＸ線の振幅を変調する振幅型の回折格子（遮蔽格子）でもよいが、Ｘ線の損失が少ないため位相格子が用いられることが多い。遮蔽格子１４０には、Ｘ線透過部とＸ線遮蔽部とが配列された格子が用いられることが多い。

測定時には、被検体１２０を回折格子１３０の前又は回折格子１３０と遮蔽格子１４０との間に設置する。Ｘ線は一般に透過性が高いために被検体１２０を透過するが、その際に透過した物質の分子構造と密度に応じた位相の変化が生じる。この位相の変化は第一の干渉パターン１８０の配置に影響を与える。そのため遮蔽格子１４０によるモアレにも歪みを生じさせる。この歪みは演算部１６０で計算することで位相の微分情報として取得できる。以上がＸ線位相イメージング装置の概要である。

次に演算部１６０により実行される位相回復法について説明する。位相回復法には大きく分けてフーリエ変換法と縞走査法がある。フーリエ変換法は一枚のフリンジ像から、縞走査法は複数枚のフリンジ像から計算によって微分位相を計算する手法である。フーリエ変換法ではフリンジ像をフーリエ変換し、微分位相情報を含むスペクトルを取り出し、そのスペクトルを逆フーリエ変換し、偏角をとることで微分位相を求める。また、縞走査法では複数のフリンジ像に複素定数を乗算した値を加算し、その値の偏角をとることで微分位相を求める。いずれの手法においても偏角を取る前の複素位相情報を一旦取得することを特徴とする。

例えば、フーリエ変換法の場合を説明する。実空間の座標を（ｘ，ｙ）、フーリエ空間（波数空間）の座標を（ｋ_ｘ，ｋ_ｙ）、フリンジ像をａ（ｘ，ｙ）、そのフーリエ変換をＡ（ｋ_ｘ，ｋ_ｙ）と表記すると、ｘ軸方向の微分位相情報を求めるには、

により、ｘ軸方向の複素位相情報Ｐ_ｘ（ｘ，ｙ）を求める。ここで、Ｇ（ｋ_ｘ，ｋ_ｙ）はフーリエ空間上の窓関数であり、求めたいｘ軸方向の微分位相情報のスペクトル中心に相当する（ｑ_ｘ，ｑ_ｙ）の点を中心とした所定範囲の成分を取り出すために使用される。また、Ｆ^−１は逆フーリエ変換を示す。つまり、式１の処理は、フリンジ像のフーリエ変換から、窓関数によりｘ軸方向の微分位相情報に対応するスペクトルのみ抽出し、抽出した
スペクトルの中心を原点（ｋ_ｘ，ｋ_ｙ）＝（０，０）に平行移動した後、逆フーリエ変換する、というものである。ｙ軸方向の複素位相情報Ｐ_ｙ（ｘ，ｙ）に関しても同様に求める（ただし窓関数が異なる）。

また、縞走査法では回折格子１３０と遮蔽格子１４０の相対位置を変えつつ、複数のフリンジ像ａ^ｎ（ｘ，ｙ）を取得したのちに、以下の式でｘ軸方向の複素位相情報Ｐ_ｘ（ｘ，ｙ）を求める。

ここで、ｒ_ｘｘ ^ｎとｒ_ｘｙ ^ｎは複素定数であり、回折格子１３０と遮蔽格子１４０の相対位置によって決定される各画素におけるフリンジの位相を表す定数である。ｉは虚数単位である。ｙ軸方向の複素位相情報Ｐ_ｙ（ｘ，ｙ）に関しても同様に求める（ただし複素定数が異なる）。

いずれの手法にせよ、計算の過程で複素位相情報Ｐ_ｘ（ｘ，ｙ）、Ｐ_ｙ（ｘ，ｙ）を取得したのちに、
Ａｒｇ［Ｐ_ｘ（ｘ，ｙ）］、Ａｒｇ［Ｐ_ｙ（ｘ，ｙ）］ （式３）
により偏角を計算することで、ｘ軸方向とｙ軸方向の微分位相が求められる。Ａｒｇは偏角をとる操作を意味する。このような手法では、偏角をとる際に位相の畳みこみが生じるため、微分位相の値は−πからπの間で限定されてしまう。このような畳みこみは取得される微分位相の定量性を悪化させるほか、画像に表示した場合にアーチファクトとなって定性的な印象にも影響を与える。さらにこうして取得された画像から積分演算や特許文献１に示すような画像処理を行うと、その影響がさらに拡大される。

これら従来の技術による結果について、コンピュータ・シミュレーションを用いて示したものが図２（ａ）〜図２（ｃ）、および図３である。図２（ａ）はシミュレーションで用いた被検体であり、球状のプラスチック球を想定している。図２（ｂ）は図１の装置をシミュレートした結果取得されたフリンジ像（モアレ像）である。このフリンジ像をフーリエ変換法で解析すると、図２（ｃ）のような二次元の微分位相像が取得される（左側がｘ軸方向の微分位相像、右側がｙ軸方向の微分位相像である）。しかし、これらの微分位相像では、被検体の端部など大きく値が変わる部分において、符号２１１で示したような位相の畳みこみが発生していることがわかる。

図２（ｃ）の微分位相像を用いて特許文献１に開示されている手法を用いて作成した画像が図３になる。図３は特許文献１の式６のρにあたる画像である。特に被検体の上下左右の端部においてアーチファクト３０１が生じていることが分かる。

このような位相畳みこみを元の位相の値に戻す操作が位相接続である。特許文献１では簡易的な位相接続の方法として、二次微分を行いその値がπ以上になっているものに対して２πを引くもしくは足すという簡易的な操作を行っている。このような操作は当然、ｘ軸方向の微分位相、ｙ軸方向の微分位相の双方に対して行う必要があり、処理コストがかかる。このような簡易的な方法では解決できない畳み込みもあり、そのような場合にはさらに複雑なアルゴリズムを用いて位相接続を行う必要がある。いずれにせよ、特許文献１に開示された手法は位相畳みこみの影響を受けやすく、それを避けるために位相接続を行うと画像を取得するために時間を要する。それゆえ、例えば空港における手荷物検査や製品の出荷検査、医療検査など即時性を要する画像処理においては特許文献１に開示の手法は適用が難しい。

そこで、本発明の実施形態に係る撮像システムでは、図８に示す画像処理方法を用いて被検体画像の再構成を行う。まず演算部１６０が、微分干渉計である撮像装置１０で得られたフリンジ像（干渉像）のデータを取得する（ステップ８００）。そして、演算部１６０が、フリンジ像のデータを基に波数空間における所望の範囲のスペクトルを取得し（ステップ８０１）、取得されたスペクトルに対し反対称性をもつフィルタをかける（ステップ８０２）。その後、フィルタされた結果を用いて被検体画像データの生成を行う（ステップ８０３）。各ステップは、独立して行う必要はなく、同時に行っても良い。

ステップ８０１の処理は、フリンジ像のデータから微分位相情報に対応するスペクトルを抽出する操作である。二次元の場合は、ｘ軸方向の微分位相情報に対応する第１のスペクトルとｙ軸方向の微分位相情報に対応する第２のスペクトルの両方をそれぞれ取得する。例えば、フリンジ像のデータから前述した式１又は式２で示される複素位相情報Ｐ_ｘ（ｘ，ｙ）とＰ_ｙ（ｘ，ｙ）を計算し、それらをフーリエ変換することで、第１のスペクトルと第２のスペクトルを得ることができる。あるいは、式１から分かるように、位相回復処理にフーリエ変換法を用いた場合には、演算の途中でＰ_ｘ（ｘ、ｙ）及びＰ_ｙ（ｘ、ｙ）のフーリエ変換を取得するプロセス（つまりステップ８００と８０１と同じ処理）が含まれる。よって、逆フーリエ変換を行う前のデータ（式１の右辺の［］内のデータ）に対し、ステップ８０２の操作を行うことで、本方法と同じ効果が得られる。他にも、フーリエ変換法、縞走査法、又はその他なんらかの手法で求めた位相像から再びＰ_ｘ（ｘ、ｙ）とＰ_ｙ（ｘ、ｙ）を求め、それらのフーリエ変換を計算してもよい。

ステップ８０２で用いる「反対称性をもつフィルタ」とは、ステップ８０１で抽出したスペクトルの中心（スペクトルピークの頂点位置）に関して反対称な関数で表現されるフィルタ、又は、そのような関数の項を含む多項式関数で表現されるフィルタである。二次元の場合には、第１のスペクトルに対しては、第１のスペクトルの中心を通り、かつ、ｋ_ｘ軸に垂直な平面に関して反対称な第１の関数もしくは第１の関数の項を含む多項式関数で表現される第１のフィルタを適用する。また、第２のスペクトルに対しては、第２のスペクトルの中心を通り、かつ、ｋ_ｙ軸に垂直な平面に関して反対称な第２の関数もしくは第２の関数の項を含む多項式関数で表現される第２のフィルタを適用する。第１のフィルタ（第１の関数）及び第２のフィルタ（第２の関数）の具体例は実施例１〜実施例４にて後述する。なお、第１の関数及び第２の関数は実施例１〜実施例３で示すような実関数でもよいし、実施例４で示すような複素関数でもよい。複素関数の場合は、実部あるいは虚部に反対称性があればよい。尚、フーリエ変換法を行う場合、フリンジ像をフーリエ変換して得られる波数空間から微分位相情報のスペクトルを取り出す際に、反対称性をもつフィルタをかけることでステップ８０１と８０２とを同時に行うことができる。

ステップ８０３では、例えば、フィルタ後のスペクトルを逆フーリエ変換等により波数空間から実空間に変換することで、被検体画像を生成する。二次元の場合には、フィルタ後の第１のスペクトルと第２のスペクトルのそれぞれの実部を合成した後、逆フーリエ変換等を行うとよい。このようにして得られた被検体画像は、ステップ８０２で用いるフィルタの選び方しだいで、様々な特徴をもつ画像となる。詳細は後述するが、例えば、エッジが強調された微分位相像や位相接続された微分位相像に相当する画像が得られる。いずれの場合でも位相の畳み込みがほとんどない高品位な被検体画像が得られる。ステップ８０３では、さらにこれらの画像を用いて位相像を生成したり、位相像を再び微分して微分位相像を生成したりしてもよい。

以上述べた画像処理方法によれば、位相接続処理を行わなくても、波数空間上でのフィルタ処理を行うだけで、位相の畳み込みの影響が少ない被検体画像を取得することができる。したがって、従来方法に比べて処理コストの低減ならびに処理時間の大幅な短縮を図
ることができる。なお、本実施形態の画像処理方法は、演算にフーリエ変換や逆フーリエ変換を含むが、高速フーリエ変換などの手法を用いることで短時間での処理が可能である。

次に、本実施形態の画像処理方法の具体例を実施例１〜実施例４に示す。
［実施例１］
図４（ａ）は、実施例１の画像処理方法の流れを示している。まず演算部１６０が、撮像装置１０からフリンジ像のデータを取得する（ステップ４００）。そして、式１又は式２で示されるｘ軸方向の複素位相情報Ｐ_ｘ（ｘ，ｙ）とｙ軸方向の複素位相情報Ｐ_ｙ（ｘ，ｙ）を計算し（ステップ４０１）、そのフーリエ変換を計算する（ステップ４０２）。ステップ４０１と４０２の処理が、ｘ軸方向の微分位相情報に対応する第１のスペクトル及びｙ軸方向の微分位相情報に対応する第２のスペクトルを取得する処理（図８のステップ８０１）に対応する。続いて、演算部１６０が、第１のスペクトル（Ｐ_ｘ（ｘ，ｙ）のフーリエ変換）に対し第１のフィルタをかけるとともに、第２のスペクトル（Ｐ_ｙ（ｘ，ｙ）のフーリエ変換）に対し第２のフィルタをかける（ステップ４０３）。そして、第１のフィルタ後の第１のスペクトルの実部と第２のフィルタ後の第２のスペクトルの実部とを合成し（ステップ４０４）、その合成結果を逆フーリエ変換して被検体画像を生成する（ステップ４０５）。なお、図４（ａ）では、ステップ４０１〜４０３においてｘ軸に関する演算とｙ軸に関する演算とを並列処理しているが、それらをシリアルに演算してもよい。

実施例１では、第１のスペクトルに作用させる第１のフィルタとして、式４の第１の関数ｆ_ｘ（ｋ_ｘ，ｋ_ｙ）からなるフィルタを用い、

第２のスペクトルに作用させる第２のフィルタとして、式５の第２の関数ｆ_ｙ（ｋ_ｘ，ｋ_ｙ）

からなるフィルタを用いる。ただし、第１のスペクトルの中心、第２のスペクトルの中心はいずれも波数空間の原点（ｋ_ｘ，ｋ_ｙ）＝（０，０）にあるものとする。

実施例１の方法で取得された被検体画像を図４（ｂ）に示す。実施例１のフィルタ（式４、式５）を用いると、被検体画像として、微分位相のエッジが強調された画像が得られる。この画像は特許文献１の式６のρにあたる画像（図３参照）に近い特徴を有しているが、図３に比べて位相の折り畳みによるアーチファクトが消失していることが分かる。また、実施例１の被検体画像（図４（ｂ））を特許文献１の式６のρの代わりに用いて、特許文献１の式７のポワソン方程式を解くことで、位相を計算することも可能であり、このようにして得られた画像が図４（ｃ）である。図２（ｃ）に比べるとアーチファクトが消失しており画質が向上していることが分かる。さらに、図４（ｂ）もしくは図４（ｃ）の情報から積分像を求めることも可能である。

［実施例２］
以下では実施例２について説明する。実施例２の画像処理方法の流れを図５（ａ）に示す。ステップ５００〜５０５の処理は図４（ａ）のステップ４００〜４０５の処理と同じであり、実施例１との違いはステップ５０３で用いるフィルタが違う点のみであるため、処理の詳細については説明を割愛する。

実施例２では、第１のスペクトルに作用させる第１のフィルタとして、式６の第１の関数ｆ_ｘ（ｋ_ｘ，ｋ_ｙ）からなるフィルタを用い、

第２のスペクトルに作用させる第２のフィルタとして、式７の第２の関数ｆ_ｙ（ｋ_ｘ，ｋ_ｙ）

からなるフィルタを用いる。

実施例２の方法で取得された被検体画像を図５（ｂ）に示す。実施例２のフィルタ（式６、式７）を用いると、実施例１のフィルタと同様、被検体画像として、微分位相のエッジが強調された画像が得られる。この画像も特許文献１の式６のρにあたる画像（図３参照）に近い特徴を有しているが、図３に比べて位相の折り畳みによるアーチファクトが消失していることが分かる。また、実施例２の被検体画像（図５（ｂ））を特許文献１の式６のρの代わりに用いて、特許文献１の式７のポワソン方程式を解くことで、位相を計算することも可能であり、このようにして得られた画像が図５（ｃ）である。図２（ｃ）に比べるとアーチファクトが消失しており画質が向上していることが分かる。さらに、実施例１と同様に、図５（ｂ）もしくは図５（ｃ）の情報から積分像を求めることも可能である。

［実施例３］
以下では実施例３について説明する。実施例３の画像処理方法の流れを図６（ａ）に示す。ステップ６００〜６０５の処理は図４（ａ）のステップ４００〜４０５の処理と同じであり、実施例１との違いはステップ６０３で用いるフィルタが違う点のみであるため、処理の詳細については説明を割愛する。

実施例３では、第１のスペクトルに作用させる第１のフィルタとして、式８の第１の関数ｆ_ｘ（ｋ_ｘ，ｋ_ｙ）からなるフィルタを用い、

第２のスペクトルに作用させる第２のフィルタとして、式９の第２の関数ｆ_ｙ（ｋ_ｘ，ｋ_ｙ）

からなるフィルタを用いる。ここで、Ｓｇｎとは、関数の符号に応じて−１（負の場合）、＋１（正の場合）、０（０の場合）を返す符号関数である。

実施例３の方法で取得された被検体画像を図６（ｂ）に示す。実施例１と実施例２の結果と同様に微分位相のエッジが強調された画像となる。また、図３の画像と比べると、アーチファクトがほとんどないことが分かる。

また、図６（ｂ）の被検体画像をフーリエ変換し、その結果に

を乗じて逆フーリエ変換し、それぞれの実部、虚部を図示したものが図６（ｃ）である。これは位相像に相当する。位相像もまた図２（ｃ）で示したアーチファクトが消えており画質が改善したことが分かる。また、実施例１および２と同様に図６（ｂ）もしくは図６（ｃ）の情報から積分像を求めることも可能である。

［実施例４］
実施例４ではフィルタで用いられる反対称な関数が複素関数でもよいことを示す。実施例４の画像処理方法の流れを図７（ａ）に示す。ステップ７００〜７０５の処理は図４（ａ）のステップ４００〜４０５の処理と同じであり、実施例１との違いはステップ７０３で用いるフィルタが違う点のみであるため、処理の詳細については説明を割愛する。

実施例４では、第１のスペクトルに作用させる第１のフィルタとして、式１１の第１の関数ｆ_ｘ（ｋ_ｘ，ｋ_ｙ）からなるフィルタを用い、

第２のスペクトルに作用させる第２のフィルタとして、式１２の第２の関数ｆ_ｙ（ｋ_ｘ，ｋ_ｙ）

からなるフィルタを用いる。

実施例４の方法で取得された被検体画像を図７（ｂ）に示す。図３の画像と比べると、アーチファクトがほとんどないことが分かる。図７（ｂ）の画像から実施例３と同様の手法で位相像に相当する画像を求めたものが図７（ｃ）である。これによると実施例１から実施例３の方法で得られる位相像よりもエッジが強調された擬似位相像が取得されている。図２（ｃ）と比べると位相畳みこみによるアーチファクトも存在しない。このような画像は例えば被検体エッジの位置と方向を同時に表示するような画像処理に使用可能であり、より可視性の高い特徴像を提供可能である。

以上より、本実施例の画像処理方法により、Ｘ線位相イメージング装置において可視性の高い被検体の特徴的画像をより簡便に取得可能である。

（他の実施例）
以上述べた実施例は本発明の一具体例を示したものにすぎず、本発明の範囲をそれらの具体例に限定する趣旨のものではない。

例えば、上記実施例では、二次元微分干渉計で得られた二次元情報に対する処理を説明
したが、本発明は一次元微分干渉計、あるいは二次元微分干渉計で取得した情報の一軸方向のみの情報を用いても同様の操作を行うことが可能である。また、Ｘ線以外の様々な波長、種類の電磁波を用いた干渉計（例えば、光学干渉計）などにも適用可能である。また、本発明が適用される撮像装置は、Ｘ線源又は電磁波源と別体に構成され、Ｘ線源又は電磁波源と組み合わせられることで撮像が可能になる構成でも良い。尚、本発明及び本明細書において撮像装置とは、周期性をもつパターンを撮像する装置であれば良く、被検体の情報を画像化するものに限定されない。

上記実施例で用いたフィルタは一例であり、反対称性をもつフィルタであれば、どのような関数のフィルタであっても上記実施例と同様の効果を期待できる。例えば、実施例１では、第１のフィルタとしてｆ_ｘ（ｋ_ｘ，ｋ_ｙ）＝ｋ_ｘなる関数（原点を通る傾きπ／２の平面）を用いたが、ｆ_ｘ（ｋ_ｘ，ｋ_ｙ）＝ａ×ｋ_ｘ（ａ：定数）のようにスケール（傾き）を変更してもよい。あるいは、ｆ_ｘ（ｋ_ｘ，ｋ_ｙ）＝ｋ_ｘ＋ｂ（ｂ：定数）のようにオフセット（切片）を変更してもよい。また、ｋ_ｘやａ×ｋ_ｘを一つの項として含む多項式関数を第１のフィルタとして用いることもできる。第２のフィルタについても同様である。また、実施例２〜４で用いたフィルタについても同様である。

上述した画像処理装置の具体的な実装は、ソフトウェア（プログラム）による実装と、ハードウェアによる実装のいずれも可能である。例えば、画像処理装置に内蔵されたコンピュータ（マイコン、ＣＰＵ、ＭＰＵ、ＦＰＧＡ等）のメモリにコンピュータプログラムを格納し、当該コンピュータプログラムをコンピュータに実行させて、各処理を実現させてもよい。また、本発明の全部または一部の処理を論理回路により実現するＡＳＩＣ等の専用プロセッサを設けることも好ましい。また、本発明は、クラウド環境におけるサーバーにも適用可能である。

また、例えば、記憶装置に記録されたプログラムを読み込み実行することで前述した実施形態の機能を実現するシステムや装置のコンピュータによって実行されるステップからなる方法によっても、本発明を実施することができる。この目的のために、上記プログラムは、例えば、ネットワークを通じて、又は、上記記憶装置となり得る様々なタイプの記録媒体（つまり、非一時的にデータを保持するコンピュータ読取可能な記録媒体）から、上記コンピュータに提供される。よって、上記コンピュータ（ＣＰＵ、ＭＰＵ等のデバイスを含む）、上記方法、上記プログラム（プログラムコード、プログラムプロダクトを含む）、上記プログラムを非一時的に保持するコンピュータ読取可能な記録媒体は、いずれも本発明の範疇に含まれる。

１１０：Ｘ線源、１２０：被検体、１３０：回折格子、１４０：遮蔽格子、１５０：検出器、１６０：演算部、１７０：画像表示装置、１８０：干渉パターン

Claims (14)

1. コンピュータが、微分干渉計により得られた干渉像のデータを取得する取得ステップと、
   コンピュータが、前記干渉像の波数空間上のスペクトルに対し、前記スペクトルの中心に関して反対称な関数もしくは該関数の項を含む多項式関数で表現されるフィルタをかけるフィルタステップと、
   コンピュータが、フィルタされた結果を用いて画像データを生成する生成ステップと、を含むことを特徴とする画像処理方法。
2. 前記干渉像のデータは、微分位相情報を含むデータであり、
   前記フィルタステップでは、前記干渉像のデータから前記微分位相情報に対応するスペクトルを取得し、前記取得されたスペクトルに対し前記フィルタをかける
   ことを特徴とする請求項１に記載の画像処理方法。
3. 前記干渉像のデータは、ｘ軸方向の微分位相情報とｙ軸方向の微分位相情報とを含むデータであり、
   前記フィルタステップでは、前記干渉像のデータから前記ｘ軸方向の微分位相情報に対応する第１のスペクトルと前記ｙ軸方向の微分位相情報に対応する第２のスペクトルとをそれぞれ取得し、前記第１のスペクトルに対し第１のフィルタをかけ、前記第２のスペクトルに対し第２のフィルタをかける
   ことを特徴とする請求項２に記載の画像処理方法。
4. 前記第１のフィルタは、前記第１のスペクトルの中心を通り、かつ、ｋ_ｘ軸（ｋ_ｘはｘ軸方向の波数）に垂直な平面に関して反対称な第１の関数もしくは該第１の関数の項を含む多項式関数で表現されるフィルタであり、
   前記第２のフィルタは、前記第２のスペクトルの中心を通り、かつ、ｋ_ｙ軸（ｋ_ｙはｙ軸方向の波数）に垂直な平面に関して反対称な第２の関数もしくは該第２の関数の項を含む多項式関数で表現されるフィルタである
   ことを特徴とする請求項３に記載の画像処理方法。
5. 前記第１のスペクトルの中心をｋ_ｘ＝０とした場合に、前記第１の関数は、
   

   

   であり、
   前記第２のスペクトルの中心をｋ_ｙ＝０とした場合に、前記第２の関数は、
   

   

   である
   ことを特徴とする請求項４に記載の画像処理方法。
6. 前記第１のスペクトルの中心をｋ_ｘ＝０とした場合に、前記第１の関数は、
   

   

   であり、
   前記第２のスペクトルの中心をｋ_ｙ＝０とした場合に、前記第２の関数は、
   

   

   である
   ことを特徴とする請求項４に記載の画像処理方法。
7. 前記第１のスペクトルの中心をｋ_ｘ＝０とした場合に、前記第１の関数は、
   

   

   であり、
   前記第２のスペクトルの中心をｋ_ｙ＝０とした場合に、前記第２の関数は、
   

   

   である
   ことを特徴とする請求項４に記載の画像処理方法。
8. 前記第１のスペクトルの中心を（ｋ_ｘ，ｋ_ｙ）＝（０，０）とした場合に、前記第１の関数は、
   

   

   であり、
   前記第２のスペクトルの中心を（ｋ_ｘ，ｋ_ｙ）＝（０，０）とした場合に、前記第２の関数は、
   

   

   である
   （ただし、Ａｒｇは偏角をとる操作であり、ｉは虚数単位である）
   ことを特徴とする請求項４に記載の画像処理方法。
9. 前記生成ステップでは、第１のフィルタ後の第１のスペクトルの実部と第２のフィルタ後の第２のスペクトルの実部とを合成した結果を波数空間から実空間に変換することにより、前記画像データを生成する
   ことを特徴とする請求項３〜８のうちいずれか１項に記載の画像処理方法。
10. コンピュータが、前記画像データから位相像又は微分位相像のデータを生成するステップをさらに有する
    ことを特徴とする請求項１〜９のうちいずれか１項に記載の画像処理方法。
11. 請求項１〜１０のうちいずれか１項に記載の画像処理方法の各ステップをコンピュータに実行させることを特徴とするプログラム。
12. 微分干渉計により得られた干渉像のデータを取得する取得手段と、
    前記干渉像の波数空間上のスペクトルに対し、前記スペクトルの中心に関して反対称な関数もしくは該関数の項を含む多項式関数で表現されるフィルタをかけるフィルタ手段と、
    フィルタされた結果を用いて画像データを生成する生成手段と、
    を有することを特徴とする画像処理装置。
13. 微分干渉計からなる撮像装置と、
    前記撮像装置により得られた干渉像のデータに基づき画像データを生成する、請求項１２に記載の画像処理装置と、
    を有することを特徴とする撮像システム。
14. 前記微分干渉計がＸ線トールボット干渉計であることを特徴とする請求項１３に記載の撮像システム。

Images