JP6495945B2 - 微分位相コントラストイメージングにおけるスペクトル位相アンラッピングのためのバイアスフリー正則化 - Google Patents

Description

本発明は、測定された位相シフトデータの処理方法、信号処理装置、コンピュータプログラム要素、及びコンピュータ読み取り可能媒体に関する。

医療分野を含む技術的応用において、位相コントラストイメージングは、軟組織におけるコントラストが優れているため、吸収のみによるイメージング方法よりも優れていると認識されている。幾つかの応用では、位相コントラストイメージングは、干渉計をイメージング装置に搭載することを必要とする。イメージングされるオブジェクトにより生じる所望の屈折及びそれによる位相シフトに関係することが分かった干渉パターンの一定の位相シフトを、干渉計の配置により、抽出することができる。残念ながら、干渉計により取得される位相シフトデータは不明確なことがある。これは、干渉計の角度範囲が限定されているからである。その範囲を超えるどんな測定も、その範囲に「ラップ（ｗｒａｐｐｅｄ）」され、実際、観察されるのは位相シフト自体ではなく、「その範囲を法とする（ｍｏｄｕｌｏ ｔｈｅ ｒａｎｇｅ）」測定位相シフト値である。位相ラッピング（ｐｈａｓｅ ｗｒａｐｐｉｎｇ）は、位相コントラスト画像に可視アーティファクトを生じることがあるので、望ましくない。この不明確さに対抗するため、過去、多くの「位相アンラッピング（ｐｈａｓｅ ｕｎ−ｗｒａｐｐｉｎｇ）」アルゴリズムが提案されている。しかし、既存のアルゴリズムにより得られるアンラッピングされた位相データの精度は、時として高くない。

それゆえ、測定された位相シフトデータを処理する時に精度を向上する必要があるだろう。

本発明の目的は、独立請求項の主題により達成され、さらに別の実施形態は従属請求項に記載されている。

留意点として、以下に説明する本発明の態様は、信号、信号処理装置、コンピュータプログラム要素、及びコンピュータ読み取り可能媒体に等しく当てはまる。

本発明の第１の態様による方法は、位相シフトデータをアンラッピングする方法であって、
−検出器が測定した位相シフトデータを受け取るステップと、
−正則化項を有しない非凸状又は非凹状の目的関数に、パラメータとして、前記位相シフトデータを含めるステップと、
−前記目的関数の臨界点を見つけるステップと、
−前記臨界点のみの中から、正則化目的関数のターゲット臨界点を見つけるステップであって、前記正則化目的関数は前記目的関数と正則化項とにより校正されるステップと、
−アンラッピングされた位相シフトデータの推定値として、前記ターゲット臨界点を出力するステップとを含む。

第２の態様による方法は、位相シフトデータをアンラッピングする別の方法であって、
−検出器が測定した位相シフトデータを受け取るステップと、
−非凸状又は非凹状の目的関数と正則化項とにより構成された正則化目的関数に、パラメータとして、前記位相シフトデータを含めるステップと、
−前記正則化目的関数の臨界点を見つけるステップであって、前記臨界点はリード臨界点を含むステップと、
−前記正則化項を有しない前記目的関数の非正則化臨界点を見つけるステップと、
−見つかった前記非正則化臨界点のうちから、少なくとも１つの他の非正則化臨界点のうち少なくとも１つより前記リード臨界点に近いターゲットポイントを選択するステップと、
−アンラッピングされた位相シフトデータの推定値として、選択された前記ターゲット点を出力するステップとを含む。

換言すると、上記方法は、入力位相シフトデータの位相ラッピングアーティファクトを補正し、そう補正したデータを、アンラッピングされた位相シフトデータとして出力する。

ここで用いる目的関数は、最適化理論の分野で既知のものである。言い換えると、目的関数は、前記関数の領域の要素を実数にマッピングする実数値関数である。目的関数は、場合に応じて、効用又はコスト関数と呼ばれることもある。これらの関数の臨界点は、その目的関数が極大又は極小（もしくは最大又は最小）となる、その関数の領域中の偏角又は位置である。それゆえ臨界点は場合によっては最大値又はマキシマイザ（ｍａｘｉｍｉｚｅｒ）又はミニマイザ（ｍｉｎｉｍｉｚｅｒ）とも呼ばれる。ここでは目的関数の最大化又は最小化のどちらも想定されている。いずれかの方法で臨界点を探索するステップは、一実施形態では、既知の数値アルゴリズムにより、目的関数を解いて、前記臨界点を求めることになる。

提案の方法により、正則化項に存在により生じるバイアス効果を回避又は少なくとも緩和することができる。同様に正則化項（ｒｅｇｕｌａｒｉｚｅｒ）又は正則化（ｒｅｇｕｌａｒｉｚｅｒ）は、通常は単に目的関数の領域における要素の関数である実数値関数である。これにより、例えば、探索する臨界点の平滑性又は「小ささ」など所望の特性の強制が可能となる。観察されているところでは、目的関数が複数の極小値を有するとき（すなわち、目的関数が凸状でも凹状でもないとき）、正則化項は、目的関数の最小値を領域の正しい範囲にシフトするようにバイアスをかける。これは望ましいが、目的関数は、その極値に、前記範囲に向かう望ましくないバイアス（例えば、絶対値が小さい）を生じる。探しているアンラッピングされた位相データはこれら極値のうちにあるので、このバイアスは、正則化された最適化に基づく多くの位相アンラッピングアルゴリズムの出力に生き残る。

このバイアス効果を解消するため、ここで２つの代替的方法を提案する。一方法では、まず探索されるのは正則化しない目的関数であり、次に結果が正則化された関数に適用される。別の実施形態では、まず正則化された関数が探索され、次に結果が非正則化目的関数に適用される。換言すると、探索ステップの順序は逆である。両方の方法において、出力（すなわち、アンラッピングされた位相シフトデータ）は結局のところ非正則化目的関数の点である。言い換えると、正則化された目的関数の臨界点は、アンラッピングされた位相シフトデータとしては返されない。むしろ、正則化された目的関数の臨界点は、非正則化目的関数の臨界点から選択するガイド手段として用いられる。この方法で、正則化項の存在により生じる探索ステップにおけるバイアスを、回避又は少なくとも緩和することができるが、正則化要件を考慮することもできる。さらに言い換えると、本方法により、正則化する必要性と、前記正則化項により生じるバイアスを回避又は緩和することとを両立させることができる。

一実施形態では、前記目的関数は前記測定データについては、尤度関数である。尤度関数の土台となる確率密度は、ラッピングされたガウシアン密度、ミーゼス型密度、又はその他の好適な確率密度であって、角度データなどの「方向」データの確率モデル化に適したもののうちどれかである。

しかし、目的関数は確率的なものでなければならないというのではない。提案の方法は、目的関数の最適化を用いるどの位相アンラッピングアルゴリズムにも適用可能であり、目的関数がどう得られたか、統計的／確率的か、またはその他の信号モデリングアプローチかにはよらない。

第２の態様による方法では、一実施形態では、ターゲットポイントは、（好適な距離尺度で測定した時）リードポイントに最も近い。しかし、これは必須ではない。ターゲットポイントがリードポイントの所定近傍内にあるだけで十分であり、すべての非正則的臨界点を考慮したとき、リードポイントに「最も近い」ものではない別の実施形態も想定されるからである。

一実施形態では、正則化項はチホノフ（Ｔｉｋｈｏｎｏｖ）正則化である。これにより、検出器の空間的解像度が限定されているため、測定データがアンダーサンプリングされていてもよいが説明できる。

一実施形態では、前記方法のどちらでも、あるピクセルの測定データは、隣接ピクセルの（の一部又は全部）で測定されたデータから独立に処理される。言い換えると、前記方法は、あるピクセルで測定されたデータに適用されたとき、検出器Ｄの前記あるピクセルの隣接ピクセルで測定されたデータに依存しない。

「リード（ｌｅａｄ）」臨界点は、（正則化された）目的関数が、その目的関数により評価されたときに少なくとも１つの他の臨界点のうち少なくとも１つが返す値より「ベターな（ｂｅｔｔｅｒ）」値を返す臨界点である。修飾語「ベター」は手元の最適化タスクに依存する。例えば、最小化問題では、値は小さいほどベターである。逆に、最大化問題では、値は大きいほどベターである。 特に、一実施形態では、リード臨界点は、目的関数が他のすべての臨界点のうち最大又は最小となる臨界点であってもよい。しかしこれは必須ではない。リードポイントが、単に所定閾値よりも良い値を返せば十分である、又は所定数の臨界点が返すよりも良い値を返す、別の実施形態も想定しているからである。

ここで「臨界点（ｃｒｉｔｉｃａｌ ｐｏｉｎｔ）」は、目的関数が局所又は極大（又は最小又は最大）になる点である。ここで、どの最適化問題の場合でも、至る所で「臨界点」は極小又は極大どちらかである。ある最適化について、臨界点が一旦は極小点であり次に極大点であるような混合状態は、ここでは想定していない。

スペクトル微分位相コントラストイメージングのためのイメージング装置を示す図である。 位相ラッピングにより生じる画像ノイズ又はアーティファクトの一例を示す図である。 異なる実施形態による信号プロセッサを示すブロック図である。 正則化及び非正則化目的関数を示すグラフである。 位相シフトデータ処理の異なる実施形態による方法を示すフローチャートである。

図１は、位相コントラストイメージング機能を有する、特に微分位相コントラストイメージング（ＤＰＣＩ）を有するイメージングシステムＩＭの基本的コンポーネントを示す。Ｘ線放射波ＸＢを発生するＸ線源ＸＲがある。Ｘ線放射波は検査領域を通った後、検出器Ｄの検出器ピクセルｐ_ｘにより検出可能である。位相コントラストイメージング機能は、Ｘ線源ＸＲと放射線感知検出器Ｄとの間に配置された干渉計により行われる。

干渉計（非限定的な一実施形態では、タルボットタイプ又はタルボット・ラウタイプのものである）は、２つの（タルボットタイプ）回折格子Ｇ_１、Ｇ_２、又はより好ましくは３つの（タルボット・ラウタイプ）回折格子Ｇ_０、Ｇ_１及びＧ_２を含む。Ｘ線源側にある第１の減衰回折格子Ｇ_０は、そろった周期ｐ_０を有し、Ｘ線源ＸＲで放射されたＸ線放射の波面の空間的コヒーレンスを生じる。

吸収回折格子Ｇ_１（周期ｐ_１を有する）は、Ｘ線源から距離Ｄのところに配置され、さらに下流で周期ｐ_２の干渉パターンを生じる。前記干渉パターンは、検出器Ｄにより検出できる。ここで、（イメージングされる）サンプルがＸ線源と検出器との間の検査領域に導入されると、干渉パターンの位相がシフトする。この干渉パターンシフトΔφは、サンプルを通る各光路に沿って集積された屈折による位相シフトΦΔの傾きに比例する（そのため、ＤＣＰＩと呼ばれる）。言い換えると、干渉の位相変化を測定するとき、これにより、サンプルにおける屈折により生じる位相シフトのシフト（すなわち傾き）を抽出することができる。

残念ながら、干渉パターンの位相シフトは、一般的に、直接的に空間的分解をするには小さすぎる。ほとんどのＸ線検出器の分解能ではこれをできない。それゆえ、この干渉パターン位相シフトを「サンプリング」するため、干渉パターンと同じ周期ｐ_２を有する第２の減衰回折格子Ｇ_２が、回折格子Ｇ_１から距離ｌのところに配置される。干渉パターン位相シフトの実際の抽出（及びサンプルにより生じる位相勾配の抽出）は、ここに想定されるすべての異なる実施形態による多数の異なる方法で実現できる。

一般的に、微分位相抽出に必要なものは、検出器Ｄと少なくとも１つの回折格子との間の相対的な運動である。これの実現は、一実施形態では、例えばアナライザ回折格子Ｇ_２を異なる離散的な回折格子位置にわたり横方向に（すなわち、回折格子に直交するｘ方向に沿って）動かすアクチュエータを用い、各回折格子位置において各ピクセルｐｘでの強度を測定することにより行うことができる。各ピクセルにおける強度は正弦波状に振動することが分かる。言い換えると、各ピクセルは、アナライザ回折格子Ｇ_２が動く間に、時間の関数として（或いは異なる回折格子位置の関数として）、（各ピクセルにおける）異なる強度の時系列を記録する。このアプローチ（「位相ステッピング」）は、Ｆ． Ｐｆｅｉｆｆｅｒ ｅｔ ａｌ著「Ｐｈａｓｅ ｒｅｔｒｉｅｖａｌ ａｎｄ ｄｉｆｆｅｒｅｎｔｉａｌ ｐｈａｓｅ−ｃｏｎｔｒａｓｔ ｉｍａｇｉｎｇ ｗｉｔｈ ｌｏｗ−ｂｒｉｌｌｉａｎｃｅ Ｘ−ｒａｙ ｓｏｕｒｃｅｓ」（Ｎａｔｕｒｅ Ｐｈｙｓ． Ｌｅｔｔ． ２， ２５８−２６１ （２００６））に記載されている。

例えば、前に参照した文献Ｅｐｐｌｅの第２９１０４頁の式（１ａ）、（１ｂ）に記載されているように、各ピクセルｐ_ｘにおける振動する強度信号は、他の量と共に、干渉パターンΔφ（〜ΦΔ）の所望の位相シフト（位相シフトは以下、単にφと記す）を「エンコード（ｅｎｃｏｄｅｓ）」する。

干渉パターンシフトΔφの取得は、特にＡ／Ｄ変換回路を含むデータ取得回路（図示せず）によりピクセルｐ_ｘごとの検出器信号を処理して数値形式にすることにより行える。数字の処理は、フーリエ分析モジュール（図示せず）により、又は好適な曲線フィッティング演算をしてピクセルごとの測定位相シフトΔφを取得するように構成されたモジュールにより行われる。これらのパラメータの収集は、処理ユニットにより処理できる。

図１のイメージングシステムＩＭ中の検出器Ｄは、光子計数タイプのものである。検出器Ｄは、検出装置の光子計数及びエネルギー分解機能を提供する回路ＡＳＩＣを含む。この機能により、スペクトルイメージングが可能である。各エネルギー応答を、画像化されるオブジェクトが構成されていると考えられる具体的な基材（ｂａｓｅ ｍａｔｅｒｉａｌ）（骨、水、脂質、造影剤など）と関連付けできるからである。どれかのピクセルｐ_ｘでピックアップされたパルスは、異なるエネルギーレベルＥ_ｊ＝Ｅ_１、Ｅ_２、Ｅ_３を表すエネルギー閾値（「ビン」とも呼ぶ）と比較される。パルスが一定の閾値を越えると、前記エネルギー閾値のカウンタがセットされ、例えばフーリエ分析器により信号が処理され、エネルギーＥ_１、Ｅ_２又はＥ_３の位相シフトφ_ｊになる。

言い換えると、検出器の光子エネルギー分解をする実施形態では、異なるエネルギーレベルＥ_ｊで測定された位相シフトデータφ_ｊが、異なるエネルギーチャネルＥ_１乃至Ｅ_３を通して提供される（図にはかかるエネルギーレベルＥ_１、Ｅ_２及びＥ_３のみを示したが、これは単に例示を目的としたものであり、２より大きい数のエネルギーレベルもここで想定される）。

各エネルギーチャネルの回折パターンの位相シフト角度φ_ｊは、以下にさらに詳しく説明するように、信号プロセッサＳＰによりピクセルごとに処理され、位相の不明確さをアンラップする。シングルプロセッサＳＰは、ワークステーションＤＰＵなどのデジタル処理ユニット上で実行できる。アンラップされた位相角は、好適なイメージングソフトウェアにより、位相コントラスト画像にレンダリングされてもよく、スクリーンＭＴに表示され、又はさもなければ処理又は格納されてもよい。

前にも触れたように、位相ラッピングは、なかんずく、異なる位相コントラストイメージングにおいて生じるかも知れない現象である。干渉計のダイナミックレンジは通常、ある角度範囲に、例えば−πないしπの間に限定されているからである。言い換えると、干渉計は、３つ（又は少なくとも２つ）の回折格子Ｇ_０ないしＧ_２のシステムであり、干渉パターンの位相シフトの測定は実際にはできず、「２πを法とする（ｍｏｄｕｌｏ ２π）」真の位相シフトのみを測定できる：範囲外位相シフトへの干渉計の応答は、長さが許すだけ多くの回数、単位円の周りに（ラジアン単位で長さを表した）真の位相シフトを「ラッピング（ｗｒａｐｐｉｎｇ）」して、ラッピングされた長さの終端がある単位円の円周上の角度位置を読むことに対応する。巻き数に関する情報は失われる。干渉計の応答として得られるのは、終端の角度位置のみである。位相アンラッピングアルゴリズムは、ラッピングされた位相及び追加情報又はモデル仮定から、失われた情報を回復することを目的としている。この情報の喪失は、位相コントラスト画像データに「人工的な」不連続性を生じ、それが画像ノイズになり、画像中の重要な細部を潜在的に不明確にするので、望ましくない。図２のペインＡ）とＢ）の例示した位相コントラスト画像は、何がうまくいかなくなるかを示す。左側の図２Ａ）は、ファントムボディを示す（この場合、２５ｋｅＶでイメージングした人の足であり、その構造は干渉計のダイナミックレンジを越える屈折勾配を誘発する）。ペインＡ）は従来の吸収Ｘ線画像を示す。右のペインＢ）は対応する位相コントラスト画像を示す。２本の矢印は、ノイズが集積した、又は位相ラップ不連続により生じたアーティファクトを有する画像部分を示す。画像Ｂ）のその他の画像部分は、十分に滑らかであり、位相ラッピングはない。再構成などの後処理中に、又は位相ラッピングされたデータを吸収コントラスト画像データと合成する時に、位相ラッピングによりさらに嫌なアーティファクトが生じる。

それゆえ、位相シフトに関する不明確さを解消して、測定された（角度）位相シフトデータをアンラップすることが望ましい。シャノンのサンプリング条件により、多くの検出器ベースのシステムにおいて、微分位相コントラスト画像を実際にサンプリングできるとは保証できないことが分かっている。検出器の空間分解能（すなわち、隣接ピクセル間のスペース（ｉｎｔｅｒ−ｓｐａｃｅ））は、通常、（空間的）サンプリング定理を満たすほど小さくない。この事実は、標準的な位相アンラッピングアルゴリズムを使う場合、結果の信頼性を著しく損なうことが、出願人には分かった。

図３を参照して、位相の不明確性を解消する、すなわち「位相アンラッパー」として動作するように構成された信号プロセッサＳＰの２つの実施形態Ａ、Ｂが示されている。より具体的には、ここに提案のシグナルプロセッサは、最適化問題を解くことを含む既存の位相アンラッピングアルゴリズムと組み合わせることを意図したものである。

かかる位相アンラッピングアルゴリズムの一例は、前出の引例Ｅｐｐｌｅの第２９１０５頁の式（５）に記載されている。他の位相アンラッピングアルゴリズムのように、Ｅｐｐｌｅは、ある目的関数ｆの最大化又は最小化するものを見つけることを記載している。目的関数は一般的に、所定の許容範囲（関数の「領域」であり、この場合はラジアン単位で測定される角度値である）の要素を数にマッピングするコスト関数又は効用関数である。目的関数は、一般的に、多数の極小又は極大を含む２次元の曲線又は３次元以上の曲面を確定するスカラー値の関数である。最適化で探している点は、極小又は極大となる点（集合的に「極値」という）である。これらの点は、（最大化問題の場合）マキシマイザ又は（最小化問題の場合）ミニマイザ又は集合的に「臨界点」と呼ばれるが、ある最適化問題の場合には、臨界点は極小又は極大のいずれかであり、両方ではない。点が「臨界」点であるか否かは、関数が何をするかのみならず、具体的な領域、すなわち関数が適用される又は評価される要素の範囲にも依存する。領域は、連続（例えば、実数）であってもよいが、前記関数を離散的な部分領域（ｓｕｂ−ｄｏｍａｉｎ）（例えば、すべての整数）に制約する場合、どの点が臨界点であるかないかに関する結果は異なることもある。臨界点の必要ではないが十分な条件は、目的関数が微分可能であり、目的関数の勾配すなわち導関数が０になることである。目的関数の近傍の単純なサンプリングにより、見つかった臨界点がミニマイザ又はマキシマイザであるか決定できる。

以下、最大化問題を考える時（すなわち、マキシマイザを探す、探索する、又は解く時）、言うまでもなく、述べることはすべて逆の問題、すなわち最小化問題に等しく当てはまる。各最大化問題は最小化問題として書き換えることができ、逆も真だからである。言い換えると、ここに提案する方法とシグナルプロセッサは、最大化又は最小化のいずれかを想定している。

位相アンラッピングアルゴリズムの具体的な場合に戻り、測定され潜在的にラッピングされた位相角が、考えているタイプの目的関数に適用され、数値アルゴリズムが臨界点を探す。臨界点は、アンラッピングされた位相として（場合によっては、何らかのリスケーリング後に）出力されてもよく、ＣＴ画像再構成又はその他の画像処理バックエンドのため、再構成器に送られることもできる。

より具体的に、Ｅｐｐｌｅに記載された例として、提案のシグナルプロセッサＳＰが演算する目的関数は、信号モデリングアプローチに、例えば最尤（Ｍａｘｉｍｕｍ−Ｌｉｋｅｌｉｈｏｏｄ （ＭＬ））法その他に基づいても良い。

目的関数は、信号モデリングアプローチに応じて定式化され、異なるエネルギーレベルＥ_１、Ｅ_２、Ｅ_３における測定される（場合によってはラッピングされている）位相シフトは、目的関数にパラメータとして適用される。対応するアンラッピングされた位相シフトデータは、好適な数値アルゴリズム又は手法により解く又は「探索」したい自由変数として、目的関数中に含まれる。言い換えると、目的関数は、測定された、場合によってはラッピングされている位相シフトデータを、（未知の）アンラッピングされた位相シフトデータに関連付ける。信号処理装置ＳＰは、図３及び図５を参照して詳しく説明するように、好適な数値アルゴリズムを適用することにより、アンラッピングされた位相シフトデータを計算するように動作する。

上記の拡張として、（Ｅｐｐｌｅのような）多数の既知の位相アンラッピングアルゴリズムは、目的関数の臨界点を解くだけでなく、正則化された目的関数の臨界点を解く。この正則化は正則化項Ｒを用いて行う。正則化項は、時としてペナルティ関数と呼ばれることもあるが、考慮中の目的関数と算術的に合成されて、正則化目的関数を構成する項である。ペナルティ又は正則化項を用いて、最適化問題に対して、解（すなわち、臨界点）の望ましい特性を強制する。例えば、正則化項が無い目的関数を用いたとすれば、可能性のある他の解を排除する一定の平滑性条件が満たされることを望んでもよい。他の実施形態では、最適化問題に対して、大きな解よりも「小さい」解を好んでも良い。小ささはペナルティ項に組み込まれたＬ_１又はＬ_２などの好適なノルムにより測定される。

目的関数ｆに対する正則化項Ｒの効果が、最小化問題を考えた場合の図４に図示した。図４は、異なる極小値及び極大値を有する、正則化されていない目的関数ｆを示すグラフである。言い換えると、関数ｆは非凸関数である。ここに提案の方法とシグナルプロセッサは、かかる非凸／非凹関数に具体的に応用されていることが分かる。見て分かるように、関数ｆの極小は、（この場合位相角ｑ＝φの範囲となるように取られる）関数の領域の要素で生じる。（例えば、２項の関数和を構成して）この関数を正則化項Ｒと組み合わせて、正則化目的関数ｆ＋Ｒを構成する場合、そのグラフ／曲面は限定に向けて実際上押しつけられる。正則化項Ｒを用いる意図は（最小化の場合に）、原点からの距離に比例して、極小値に「ペナルティを課す（ｐｅｎａｌｉｚｅ）」ことである。言い換えると、「小さい」最小値（すなわち、原点に近いもの）を有利にする。しかし、残念ながら、他の意図しない効果もある：正則化項Ｒの適用により、臨界点自体が原点にシフトされるというバイアスｂが生じる。このバイアスは時として不正確さを引き起こすことが観察されている。ここに提案の図３に示した信号プロセッサは、このバイアスを除去し、正則化関数を使う利点を保持するが、元の問題より、返される解がより正確であるという付加的な利点を有するように動作する。再び、最小化について図４を参照して述べたことはすべて、最大化及び非凹関数に準用可能であり、これらの実施形態もここで想定している。

簡単に、図３から分かるように、信号プロセッサＳＲの実施形態Ａ）、Ｂ）は、好適なインターフェースＩＮとＯＵＴを含み、ＩＮは、ピクセルごとに、（場合によっては位相ラッピング破損（ｐｈａｓｅ ｗｒａｐ ｃｏｒｒｕｐｔｅｄ）している）測定された位相角を受け取り、ＯＵＴは、一部又は全部の測定値の位相ラッピング破損が無くなった、又は緩和された（前記ピクセルの）アンラッピングされた位相角推定値を出力する。

信号プロセッサＳＰは少なくとも２つの最適化モジュールＲとＮＲとを含む。一方のモジュールＲは、正則化された目的関数の臨界点を探すオプティマイザとして動作し、他方のモジュールＮＲは、正則化されていない目的関数の臨界点を解く又は探す。

大まかに、図３Ａ）、３Ｂ）に模式的に示したように、処理パイプラインにおける２つの最適化モジュールＮＲとＲの順序は逆にすることができる。言い換えると、図Ａの実施形態では、位相角φ_ｊは、最初に、正則化されていない目的関数に適用され、さらに下流において、この第１段階の最適化の結果が、目的関数が正則化されている第２の段階でさらに処理される。図３の実施形態Ｂ）では、この処理パイプラインは反転され、受け取られた位相角が最初に、正則化された目的関数に適用され、次いでそれからの出力が、目的関数が正則化されていない第２のモジュールに入力される。要するに、ここで提案しているものは、２段階信号処理モジュールであって、２つのオプティマイザを含み、１つは正則化項が無い目的関数に関する最適化をし、１つは正則化項がある目的関数に関する最適化をする。２つのオプティマイザモジュールはどちらの順序で適用することもできる。図３の信号プロセッサＳＰの動作を、最初に図５のフローチャートＡ）を参照して、ここで説明する。

ステップＳ３０５において、（潜在的に位相ラッピングされた）エネルギー分解位相信号φ_ｊが、イメージングシステムＩＭの干渉イメージング機器Ｇ_０−Ｇ_２を通して検出器Ｄで「見た」データとして受け取られる。

ステップＳ３１０において、そのように受け取られた位相シフトデータは、正則化項を有しない目的関数Ｆに適用される。すなわち、パラメータのプレースホルダーに測定データが追加される。前に示唆したとおり、正則化項の存在により生じる臨界点の原点へのシフトにともなう混乱は、目的関数の曲面に２以上の最小値がある場合のみにあてはまる。それゆえ、以下の発明は、非凸状（又は、オプティマイザが最大化を行うか、最小化を行うかに応じて、非凹状））の目的関数に適用するものである。

図５のフローチャートＡ）に戻り、ステップＳ３１５において、信号プロセッサＳＰの第１のオプティマイザモジュールＮＲは前記（非凸状又は非凹状の）目的関数の臨界点を探す。特に、非正則化項Ｒは、処理パイプラインのこの段階に含まれる。Ｓ３１５における探索又は解決動作は、これらの臨界点を効果的に解決する既存のどの数値的最小化又は最大化アルゴリズムによっても実現することができる。ここで想定している例は、下降シンプレックス法（ｄｏｗｎｈｉｌｌ ｓｉｍｐｌｅｘ ｍｅｔｈｏｄ）、（非線形）共役勾配法（（ｎｏｎｌｉｎｅａｒ） ｃｏｎｊｕｇａｔｅ ｇｒａｄｉｅｎｔ ｍｅｔｈｏｄ）、ニュートン・ラフソン法（Ｎｅｗｔｏｎ−Ｒａｐｈｓｏｎ ｍｅｔｈｏｄ）などである。好適なものを選択して、手元の具体的な目的関数に関する不適切性に合わせることができる。目的関数の構造は、一般的に、選択された信号モデルにより決定される。実施例は、（好適なタイプの確率分布が与えられた時の）ＭＬアプローチの統計的アプローチである。この場合、計算上の理由のため、通常は、尤度関数の対数をマイナスにしたものを最適化する。一実施形態では、今の検出器の場合には、（ラッピングされた）ガウシアン又はミーゼス型密度が信号モデルにモデル化される。

こうして見つかった臨界点は、段階Ｓ３２０で第２のオプティマイザＲに転送される。そこで前記臨界点は、段階Ｓ３１５におけるモジュールＮＲの前出の目的関数と、正則化項とにより構成される、合成された目的関数、すなわち正則化された目的関数、に関して処理される。しかし、今回は、段階Ｓ３１０とは反対に、前出の連続領域は離散的になり、すなわちステップＳ３１５で見つかった臨界点のセットの離散空間に制約される。言い換えると、ステップＳ３２０において、臨界点のみで探索を行い、正則化目的関数の目標臨界点を見つける。ステップＳ３２０で新しい正則化された最適化をする一方法は、臨界点の離散的有限集合においてその正則化目的関数を評価して、目的関数が最小値となる又は最小値を「返す」ものを見つける、又は目的関数が少なくとも所定の許容できる閾値より小さい値となるものを見つけることである。この条件を満たすより多くの点があるとき、すなわち２以上の臨界点が前記閾値より小さくなる正則化目的関数の値を返すとき、前記点のうちの１つの点を返すランダム実験を行い、その点をユーザインターフェースに適格点として表示して人間ユーザに最終的に選択させても良い。例えば、前の測定された位相シフト値を、的確な臨界点の各々で置き換えることにより、適格な臨界点の各々に対してサンプル位相コントラスト画像を生成することができ、ユーザはサンプル画像の視覚的平滑性を調べて最終的な選択をサポートすることができる。

ステップＳ３２５において、具体的な最適化アルゴリズム及び／又は用いられる目的関数に応じて、直接、又は何らかの変換もしくはスケーリングの後に、関連ピクセルの目標点が、アンラッピングされた位相シフトデータとして出力される。

一実施形態によると、上記のステップＳ３０５−Ｓ３２５は各単一ピクセルデータに対して独立に（しかし、場合によっては平行して、又は順次的に）処理される。言い換えると、どのピクセル又はプロセスの隣接データの評価も、本アプローチには含まれない。これは、ピクセルのアンラッピングされた位相値に到達するために、前記ピクセルの隣接ピクセル値の評価を要する出願人の国際出願公開第ＷＯ２０１２／０３８８５７に開示されたようなその他のアプローチとは区別される。

ここで、図５のフローチャートＢ）を参照する。これは、図３の実施形態ｂ）による信号処理ユニットにより実装可能な方法に対応する。

この方法は、実施形態Ａ）の前の方法に似ているが、処理パイプラインにおけるステップＳ３１５とＳ３２０の順序が基本的に逆になっている。

より詳細には、ステップＳ４０５において、任意のピクセルの位相シフトデータを受け取る。

測定データを、ステップＳ４１０において、モジュールＲで実装された正則化目的関数に適用する。正則化された目的関数は、ステップＳ３１０を参照して上で説明したものと同様に、非凸状目的関数と正則化項とにより構成される。

ステップＳ３１５の変形として、ステップＳ４１５において、処理パイプライン中のモジュールＲを用いて正則化目的関数の臨界点を見つける探索を行う。見つかった臨界点には、「リード（ｌｅａｄ）」臨界点、すなわち正則化された目的関数を「最も良く」満たすもの、又は少なくとも見つかったその他の臨界点の一部又は全部よりよいものがある。例えば、リード臨界点は、前記リード臨界点で正則化目的関数が評価されたときに最小値を返す、又はそう評価されたときに最大値を返す。言い換えると、リード臨界点は２つのシーンを「リード」する。それは、最適化タスクがどのように定式化されているかに応じて、コストが最小になる、又は効用が最大になるという点で優れている。しかし、前記リード点は、次の２つのステップＳ４２０、Ｓ４２５を参照して今から説明するように、探索を所望のアンラッピングされた位相シフトデータに「リード」する又は誘導する。幾つかの実施形態では、リード点は、グローバルにベストではないかも知れないが、単にベターであればよい。すなわち、残りの臨界点の大部分より正則化された目的関数により好ましい応答をかえせばよい、又は所定の閾値よりもベターな応答を返せばよい。

ステップＳ４２０において、処理フローはモジュールＮＲに移り、「非正則化」臨界点、すなわち非正則化目的関数の臨界点を見つける。すなわち、ステップＳ４１５の目的関数が、正則化項無しで、それだけで最適化される。

次いで、ステップＳ４２５において、ステップＳ４２０で見つかった非正則化臨界点のうちから、正則化最適化問題を解く時にステップＳ４１５で見つかったリード臨界点に最も近いものを選択することにより、非正則化最適化問題を解く。Ｌ_１やＬ_２などの幾つかの好適なノルムにより、近さを測定する。幾つかの実施形態では、２つ以上の「最も近い」点があるかも知れない。非正則化臨界点は、リード臨界点に、所定の距離閾値より近いことが要求されているだけだからである。後者の場合、ランダムアルゴリズムを実行して、アルゴリズムを結論的にしてもよい。あるいは、実施形態Ａ）について上で説明したように、最終的な選択のため、ユーザインターフェースを提示してもよい。

次いでステップＳ４３０において、場合によっては好適な変換又はスケーリング操作の適用後、アンラッピングされた位相シフトデータとして、ステップＳ４２５で選択された点が出力される。

両方の場合に、２つの実施形態から分かるように、出力点は常に非正則化最適化問題のものである。しかし、正則化最適化問題は、正則化要件を満たすために、まだ行われる。しかし、目的関数の正則化バージョンを最適化することにより見つかる臨界点は、非正則化目的関数に関して探索された非正則化臨界点から選択するガイド手段又は「キュー（ｃｕｅｓ）」としてのみ用いられる。これにより、ステップＳ３２５又はＳ４３０で出力される最終結果が、それぞれ正則化項の影響によるバイアスを受けないことが保証される。

上で触れたように、正則化項と目的関数の具体的な形式は、具体的な位相アンラッピングアルゴリズム及び／又はその土台として用いられる信号モデルに依存する。それゆえ、本方法は、非凸状又は非凹状の目的関数に関する最適化問題を含む既存の位相アンラッピングアルゴリズムへのアドオンと考えても良い。

また、臨界点を見つける２つのモジュールＲとＮＲで用いられる数値的手法は、同じであってもなくてもよい。

例えば、一実施形態では、Ｅｐｐｌｅに記載されたものと同様に、正則化項Ｒを含む目的関数は、次の通りであってもよい：

ここで、
ｗは異なるエネルギービンを示し、
φ_ｍは測定された、潜在的にラッピングされた位相シフトであり、
Ｍは推定されたアンラッピングされた位相シフト
（外１）

であり、
Ｅ_ｗはエネルギービンｗの有効エネルギーであり、
（外２）

は測定された位相シフトφ_ｍの標準偏差であり、
Ｒは正則化関数であり、
λは正則化関数Ｒの強さを制御する正則化パラメータである。

モデルパラメータ
（外３）

、Ｅ_ｗ、λは、校正測定により取得してもよく、λはテスト測定でベストな結果になるように調整される。

Ｅｐｐｌｅに記載されているように、上記のモデルは単色放射で用いることもできるが、式（１）の土台となる信号モデルは、アンラッピング問題の位相勾配の周知のエネルギー依存性を利用する。

式（１）による上記の目的関数ｌ（Ｍ＝φ）は、最大尤度アプローチから得られる。アンラッピングされた位相シフトデータは、最大尤度関数を用いた推定位相ｔである。より具体的に、位相勾配データは、前記最大尤度関数の対数の最小化により推定される。

位相データ測定検出器はミーゼス型密度（ｖｏｎ Ｍｉｓｅｓ ｄｅｎｓｉｔｙ）により決まると仮定した。しかし、他の密度をここに含めてもよい。例えば、ラッピングされたガウシアン密度を仮定してもよい。Ｔ Ｗｅｂｅｒその他は、「Ｍｅａｓｕｒｅｍｅｎｔｓ ａｎｄ ｓｉｍｕｌａｔｉｏｎｓ ａｎａｌｙｚｉｎｇ ｔｈｅ ｎｏｉｓｅ ｂｅｈａｖｉｏｒ ｏｆ ｇｒａｔｉｎｇ−ｂａｓｅｄ Ｘ−ｒａｙ ｐｈａｓｅ−ｃｏｎｔｒａｓｔ ｉｍａｇｉｎｇ」（Ｎｕｃｌｅａｒ Ｉｎｓｔｒｕｍｅｎｔｓ ａｎｄ Ｍｅｔｈｏｄｓ ｉｎ Ｐｈｙｓｉｃｓ Ｒｅｓｅａｒｃｈ Ａ６４８（２０１１）， ｐｐ Ｓ２７３−Ｓ２７５）において、ミーゼスが位相シフト測定値を十分な精度で記述することを確認し（Ｗｅｂｅｒの図４を参照）している。また、Ｔ Ｗｅｂｅｒその他は、「Ｎｏｉｓｅ ｉｎ ｘ−ｒａｙ ｇｒａｔｉｎｇ−ｂａｓｅｄ ｐｈａｓｅ−ｃｏｎｔｒａｓｔ ｉｍａｇｉｎｇ」（Ｍｅｄ． Ｐｈｙｓ． ３８ （７）， Ｊｕｌｙ ２０１１， ｐｐ ４１３３−４１４０）において、例えば上記の式（１）で
（外４）

を推定するのに使えるミーゼス分布（ｖｏｎ Ｍｉｓｅｓ ｄｉｓｔｒｉｂｕｔｉｏｎ）の変形例を説明している。好適な密度のその他のオプションは、Ｈ Ｇｕｄｂｊａｒｔｓｓｏｎその他により、「Ｔｈｅ Ｒｉｃｉａｎ Ｄｉｓｔｒｉｂｕｔｉｏｎ ｏｆ Ｎｏｉｓｙ ＭＲＩ Ｄａｔａ」（Ｍａｇｎ． Ｒｅｓｏｎ． Ｍｅｄ．， １９９５， Ｄｅｃ， ３４（６）， ｐｐ ９１０−９１４）の式（５）に記載されている。

言うまでもなく、測定データの統計をモデル化するのに用いられる密度又は分布のタイプに応じて、尤度関数は、式（１）とは異なる構造を有し、臨界点を見つける他の数値的手法を必要とするかも知れない。しかし、当業者には言うまでもなく、上記の二段階の正則化及び非正則化のアプローチは、どんな目的関数にも、それが最大尤度アプローチから得られたものであろうとなかろうと、適用可能である。例えば、他の数値的方法に基づく目的関数もここでは想定している。

見て分かるように、正則化項は、検討している具体的なサンプル点Ｍだけの関数である。一実施形態では、正則化項はチホノフ（Ｔｉｋｈｏｎｏｖ）正則化である。これにより、位相勾配の全位相コントラスト画像（図２のペインＢ）を参照）にエイリアシングが生じ、画像領域における平滑化制約による正則化が不都合になることを説明できる。その替わりに、Ｍの二乗にペナルティを課し、小さな値の位相勾配Ｍに有利に働くチホノフ正則化が提案されている。

一例として、目的関数（１）に図３の方法Ａ）のステップＳ３１５を適用すると、非正則対数尤度関数：

の物理的自由変数Ｍ（＝φ）の許容範囲／領域における（この場合には）極小を見つける必要がある。

特に、ステップＳ３１５には正則化はない。

数値的に見つかった極小値の位置をＭ１、Ｍ２、・・・Ｍｎと記す。式（２）には正則化は適用されていないので、これらの値におけるＭの真値の最良の推測は、バイアスの影響を受けない。正則化は後のステップＳ３２０で行われ、ここでは離散的最小化問題：

から、妥当な「ターゲット」極小が得られる。

この最適化問題は、非正則化対数尤度関数のバイアスフリー極小値Ｍ１、Ｍ２、・・・Ｍｎのみを考える点を除けば、元の正則化問題（２）に非常に似ている。換言すると、式（２）の最適化で考えられる連続領域は、部分集合｛Ｍ１，Ｍ２，…，Ｍｎ｝で定義される離散領域に制約されない。このように、位相アンラッピングの最終値Ｍｋは、正則化により誘発されるバイアスの影響を受けないか、又は影響が小さい。

あるいは、図３のＢ）の代替的方法を適用する場合、臨界点の連続領域において式（１）を解き、これらの臨界点において式（１）を評価することによりリード臨界点（ｌｅａｄ ｃｒｉｔｉｃａｌ ｐｏｉｎｔ）を求める（ｅｓｔａｂｌｉｓｈ）。次いで、連続領域で再び式（２）を解いて、非正則化臨界点を求める。次いで、出力として、前記の非正則化臨界点から、リード臨界点に最も近い又は少なくとも（距離閾値と比べて）十分近い点を選択する。

一実施形態では、画像データ処理装置ＳＰは、Ｍａｔｌａｂ（Ｒ）などの好適な科学計算プラットフォームでプログラムされ、計算システム（イメージャのワークステーションＤＰＵなど）で実行するのに適したＣ＋＋又はＣのルーチンに変換されてもよい。

本発明の他の一実施形態では、適切なシステムにおいて、上記の実施形態の一つによる方法の方法ステップを実行するように較正されたことを特徴とするコンピュータプログラムまたはコンピュータプログラム要素が提供される。

コンピュータプログラム要素は、コンピュータユニットに記憶されてもよい。コンピュータユニットも本発明の一実施形態の一部であってもよい。このコンピューティングユニットは、上記の方法のステップを実行するまたは実行を誘起するように構成され得る。さらに、上記の装置のコンポーネントを動作させるように構成されていてもよい。コンピューティングユニットは、自動的に動作し、及び／またはユーザの命令を実行するように構成されている。コンピュータプログラムはデータプロセッサのワーキングメモリにロードされる。データプロセッサは、本発明の方法を実行するように構成されている。

本発明のこの実施形態は、初めから本発明を用いるコンピュータプログラムと、アップデートにより本発明を用いるプログラムになる既存のプログラムとの両方をカバーする。

さらに、コンピュータプログラム要素は、上記の方法の実施形態の手順を満たす必要なすべてのステップを提供できる。

本発明のさらに別の一実施形態によると、ＣＤ−ＲＯＭなどのコンピュータ読み取り可能媒体が提供され、そのコンピュータ読み取り可能媒体は、前のセクションで説明したコンピュータプログラム要素を記憶したものである。

コンピュータプログラムは、光記憶媒体や他のハードウェアとともに、またはその一部として供給される固体媒体などの適切な媒体に記憶及び／または配布することができ、インターネットや有線または無線の電気通信システムなどを介して他の形式で配信することもできる。

しかし、コンピュータプログラムは、ワールドワイドウェブ等のネットワーク上で提供されてもよく、そのようなネットワークからデータプロセッサのワーキングメモリにダウンロードされてもよい。本発明のさらにべつの実施形態では、コンピュータプログラム要素をダウンロードできるようにする媒体が提供され、そのコンピュータプログラム要素は本発明の上記の実施形態の一つによる方法を実行するように構成されている。

留意すべき点として、本発明の実施形態を、異なる主題を参照して説明する。具体的に、一部の実施形態を方法の請求項を参照して説明し、他の一部の実施形態を装置の請求項を参照して説明する。しかし、本技術分野の当業者は、上記の説明と以下の説明から、特に断らないかぎり、一種類の主題に属する特徴の任意の組み合わせに加えて、異なる複数の主題に関係する特徴の間の任意の組み合わせも本出願で開示されていると考えられることが分かるであろう。しかし、すべての特徴は組み合わせて、特徴の単なる和以上のシナジー効果を提供することができる。

図面と上記の説明に詳しく示し本発明を説明したが、かかる例示と説明は例であり限定ではない。本発明は開示した実施形態には限定されない。請求項に記載した発明を実施する際、図面、本開示、及び従属項を研究して、開示した実施形態のその他のバリエーションを、当業者は理解して実施することができるであろう。

請求項において、「有する（ｃｏｍｐｒｉｓｉｎｇ）」という用語は他の要素やステップを排除するものではなく、「１つの（“ａ” ｏｒ “ａｎ”）」という表現は複数ある場合を排除するものではない。単一のプロセッサまたはその他のアイテムが請求項に記載した複数のユニットの機能を満たすこともできる。相異なる従属クレームに手段が記載されているからといって、その手段を組み合わせて有利に使用することができないということではない。請求項に含まれる参照符号は、その請求項の範囲を限定するものと解してはならない。

Claims (16)

1. 位相シフトデータをアンラッピングする方法であって、
   検出器が測定した位相シフトデータを受け取るステップと、
   正則化項を有しない非凸状又は非凹状の目的関数に、パラメータとして、前記位相シフトデータを含めるステップと、
   前記目的関数の臨界点を見つけるステップと、
   前記臨界点のみの中から、正則化目的関数のターゲット臨界点を見つけるステップであって、前記正則化目的関数は前記目的関数と正則化項とにより校正されるステップと、
   アンラッピングされた位相シフトデータの推定値として、前記ターゲット臨界点を出力するステップとを含む、方法。
2. 位相シフトデータをアンラッピングする方法であって、
   検出器が測定した位相シフトデータを受け取るステップと、
   非凸状又は非凹状の目的関数と正則化項とにより構成された正則化目的関数に、パラメータとして、前記位相シフトデータを含めるステップと、
   前記正則化目的関数の臨界点を見つけるステップであって、前記臨界点はリード臨界点を含むステップと、
   前記正則化項を有しない前記目的関数の非正則化臨界点を見つけるステップと、
   見つかった前記非正則化臨界点のうちから、少なくとも１つの他の非正則化臨界点のうち少なくとも１つより前記リード臨界点に近いターゲット臨界点を選択するステップと、
   アンラッピングされた位相シフトデータの推定値として、選択された前記ターゲット臨界点を出力するステップとを含む、方法。
3. 前記目的関数は、測定された前記位相シフトデータに関する尤度関数から得られる、
   請求項１または２に記載の方法。
4. 前記尤度関数の土台となる確率密度はラッピングされたガウシアン密度である、
   請求項３に記載の方法。
5. 前記尤度関数の土台となる確率密度はミーゼス型密度である、
   請求項３に記載の方法。
6. 前記正則化項はチホノフ正則化である、請求項１乃至５いずれか一項に記載の方法。
7. 前記方法において、あるピクセルにおいて測定された位相シフトデータは、隣接するピクセルにおいて測定された位相シフトデータとは独立して処理される、
   請求項１乃至６いずれか一項に記載の方法。
8. 位相シフトデータをアンラッピングする信号処理装置であって、
   検出器が測定した位相シフトデータを受け取る手段と、
   正則化項を有しない非凸状又は非凹状の目的関数に、パラメータとして、前記位相シフトデータを含める手段と、
   前記目的関数の臨界点を見つける手段と、
   前記臨界点のみの中から、正則化目的関数のターゲット臨界点を見つけるステップであって、前記正則化目的関数は前記目的関数と正則化項とにより校正される手段と、
   アンラッピングされた位相シフトデータの推定値として、前記ターゲット臨界点を出力する手段とを有する、
   信号処理装置。
9. 位相シフトデータをアンラッピングする信号処理装置であって、
   検出器が測定した位相シフトデータを受け取る手段と、
   非凸状又は非凹状の目的関数と正則化項とにより構成された正則化目的関数に、パラメータとして、前記位相シフトデータを含める手段と、
   前記正則化目的関数の臨界点を見つける手段であって、前記臨界点はリード臨界点を含む手段と、
   前記正則化項を有しない前記目的関数の非正則化臨界点を見つける手段と、
   見つかった前記非正則化臨界点のうちから、少なくとも１つの他の非正則化臨界点のうち少なくとも１つより前記リード臨界点に近いターゲット臨界点を選択する手段と、
   アンラッピングされた位相シフトデータの推定値として、選択された前記ターゲット臨界点を出力する手段とを有する、
   信号処理装置。
10. 前記目的関数を、測定された前記位相シフトデータに関する尤度関数から求めるように構成された、請求項８または９に記載の信号処理装置。
11. 前記尤度関数の土台となる確率密度はラッピングされたガウシアン密度である、請求項１０に記載の信号処理装置。
12. 前記正則化項としてチホノフ正則化を用いるように構成された、
    請求項８乃至１１いずれか一項に記載の信号処理装置。
13. あるピクセルにおいて測定された位相シフトデータを、隣接するピクセルで測定された位相シフトデータとは独立して処理するように構成された、
    請求項８乃至１２いずれか一項に記載の信号処理装置。
14. 前記検出器はエネルギー分解的である、
    請求項８乃至１３いずれか一項に記載の信号処理装置。
15. 請求項８乃至１４いずれか一項に記載の信号処理装置を制御するコンピュータプログラム要素であって、データ処理ユニットにより実行されると、請求項１乃至７いずれか一項に記載の方法の各ステップを実行させるコンピュータプログラム要素。
16. 請求項１５に記載のプログラム要素を記憶したコンピュータ読み取り可能媒体。

Images