(57)【特許請求の範囲】
【請求項１】
制約付き最適化を使用して微分画像から積分画像を再生する方法において、
該方法は、
【数１】
0005818909.tif 000022
と記される、制約
0005818909.tif 000023
を条件としてコスト関数
0005818909.tif 000024
を最小化する画像を識別するステップを含んでおり、
ただしｆは画像ベクトル、Ｔは変換演算子行列、Ｄxは微分演算子行列、Ｗは雑音標準偏差の逆数を含む対角行列、φは測定した画像ベクトル、εは雑音パワーについての境界であり、
0005818909.tif 000025
は、ベクトルのpノルムを示しており、ただしpは正の整数である、
ことを特徴とする方法。
【請求項２】
請求項１に記載の方法において、
１つのサンプルから量的ｘ線画像を取得するためにｘ線用装置、殊に硬ｘ線用装置によって前記微分画像を取得し、
該装置は、
ａ．Ｘ線源と、
ｂ．Ｇ０，Ｇ１およびＧ２と称される３つの回折格子、または、Ｇ１およびＧ２と称される少なくとも２つの回折格子と、
ｃ．複数の個別ピクセルを有しかつ空間変調検出感度を備えた位置感応型検出器と、
ｄ．当該検出器の画像を記録するための手段と、
ｅ．個別ピクセル毎に、吸収が主であるピクセル、および／または、微分位相コントラストが主であるピクセル、および／または、ｘ線散乱が主であるピクセルとして対象体の特性を識別するため、一連の画像においてピクセル毎の強度を評価する手段とを備え、
ｆ．ただし０からπまたは２πまで連続的または段階的に前記サンプルを回転させるか、または当該サンプルに対して前記装置および前記線源を回転させることによって前記一連の画像を収集する、
ことを特徴とする方法。
【請求項３】
請求項１または２に記載の方法において、
当該方法をいわゆる「近接方式」または「Talbot方式」のいずれかで動作させる、
ことを特徴とする方法。
【請求項４】
請求項１から３までのいずれか１項に記載の方法において、
Ｇ１は、吸収回折格子または位相回折格子であるライン型回折格子（Ｇ１）であり、前記位相回折格子は、低吸収型回折格子であるが大きなＸ線位相シフトを形成し、該Ｘ線位相シフトはπ／２またはその倍数である、
ことを特徴とする方法。
【請求項５】
請求項２に記載の方法、または請求項２を引用する請求項３または４に記載の方法において、
Ｇ２は、Ｇ１の自己画像の周期と同じ周期を有する高Ｘ線吸収コントラストを有するライン型回折格子であり、Ｇ２は、ラインがＧ１のラインと平行になっている前記検出器の前に配置される、
ことを特徴とする方法。
【請求項６】
請求項１から５までのいずれか１項に記載の方法において、
平行ビーム、準平行、扇形ビームまたは円錐形ビームモードになるように動作を選択し、Ｇ０，Ｇ１およびＧ２は対応して平坦、円筒形または円錐形の形状をそれぞれ有する、
ことを特徴とする方法。
【請求項７】
請求項１から６までのいずれか１項に記載の方法において、
２次元回折格子による全視野モードかまたは１次元回折格子による走査モードになるように動作を選択する、
ことを特徴とする方法。
【請求項８】
請求項２記載の方法、または請求項２を引用する請求項３から７までのいずれか１項記載の方法において、
近接方式動作に対し、前記回折格子間の間隔を当該近接方式内で自由に選択し、かつ、Talbot方式に対し、
【数２】
0005818909.tif 000026
ｌ＝１，２，３，…，Ｄnは平行Ｘ線ビームが使用される場合の奇数の分数Talbot距離であり、これに対してＤn,sphは扇形または円錐形のＸ線ビームが使用される場合の奇数の分数Talbot距離であり、Ｌは線源とＧ１との間の間隔である、
ことを特徴とする方法。
【請求項９】
請求項１から８までのいずれか１項に記載の方法において、
１つの回折格子Ｇ０，Ｇ１またはＧ２を、別の回折格子に対して機械式にシフトすることによって位相ステッピングを行う、
ことを特徴とする方法。
【請求項１０】
請求項２に記載の方法、または請求項２を引用する請求項３から９までのいずれか１項に記載の方法において、
前記回折格子構造をプレーナ技術によって作製する、
ことを特徴とする方法。
【請求項１１】
請求項９に記載の方法において、
前記演算子Ｄxは、前記回折格子の前記位相ステッピング方向における任意の次数の微分演算子である、
ことを特徴とする方法。
【請求項１２】
請求項９または１１に記載の方法において、
前記変換演算子Ｔは、前記回折格子の前記ステッピング方向に対して垂直な方向における任意の次数の微分演算子である、
ことを特徴とする方法。
【請求項１３】
請求項１から１２までのいずれか１項に記載の方法において、
前記重み付け演算子Ｗは、各ピクセルにおけるＤＰＣ画像の標準偏差の逆数１／σDPCを含む対角重み付け行列である、
ことを特徴とする方法。
【請求項１４】
請求項１から１３までのいずれか１項に記載の方法において、
前記制約付き最適化問題を２次錐計画（ＳＯＣＰ）に計算しなおすことによって当該制約付き最適化問題を解く、
ことを特徴とする方法。
【請求項１５】
請求項１から１４までのいずれか１項に記載の方法において、
前記制約付き最適化問題を、例えば
【数３】
0005818909.tif 000027
にしたがって制約のない形にすることによって解く、
ことを特徴とする方法。
【請求項１６】
請求項１５に記載の方法において、
勾配降下法または（非線形）共役勾配アルゴリズムまたは別の方法によって前記制約のない最適化問題を解く、
ことを特徴とする方法。
【請求項１７】
制約付き最適化を使用して微分画像から積分画像を再生するシステムにおいて、
該システムは、
ａ）  
【数４】
0005818909.tif 000028
と記される、制約
0005818909.tif 000029
を条件としてコスト関数
0005818909.tif 000030
を最小化する画像を識別する手段を含んでおり、
ｂ）ｆは画像ベクトル、Ｔは変換演算子行列、Ｄxは微分演算子行列、Ｗは、雑音標準偏差の逆数を含む対角行列、φは測定した画像ベクトル、εは、雑音パワーについての境界であり、
0005818909.tif 000031
は、ベクトルのpノルムを示しており、ただしpは正の整数である、
ことを特徴とするシステム。
【請求項１８】
請求項１７に記載のシステムにおいて、
前記微分画像は、１つのサンプルから量的ｘ線画像を得るためにｘ線用装置によって、殊に硬ｘ線用装置によって取得され、
該装置は、
ａ．Ｘ線源と、
ｂ．Ｇ０，Ｇ１およびＧ２と称される３つの回折格子、または、Ｇ１およびＧ２と称される少なくとも２つの回折格子と、
ｃ．複数の個別ピクセルを有しかつ空間変調検出感度を備えた位置感応型検出器と、
ｄ．当該検出器の画像を記録するための手段と、
ｅ．個別ピクセル毎に、吸収が主であるピクセル、および／または、微分位相コントラストが主であるピクセル、および／または、ｘ線散乱が主であるピクセルとして対象体の特性を識別するため、一連の画像においてピクセル毎の強度を評価する手段とを備え、
ｆ．  ただし、前記一連の画像は、０からπまたは２πまで連続的または段階的に前記サンプルを回転させるか、または当該サンプルに対して前記装置および前記線源を回転させることによって収集される、
ことを特徴とするシステム。
【請求項１９】
請求項１７または１８に記載のシステムにおいて、
該システムは、いわゆる「近接方式」または「Talbot方式」のいずれかで動作される、
ことを特徴とするシステム。
【請求項２０】
請求項１７から１９までのいずれか１項に記載のシステムにおいて、
Ｇ１は、吸収回折格子または位相回折格子であるライン型回折格子（Ｇ１）であり、前記位相回折格子は、低吸収型回折格子であるが大きなＸ線位相シフトを形成し、該Ｘ線位相シフトはπ／２またはその倍数である、
ことを特徴とするシステム。
【請求項２１】
請求項１７から２０までのいずれか１項に記載のシステムにおいて、
Ｇ２は、Ｇ１の自己画像の周期と同じ周期を有する高Ｘ線吸収コントラストを有するライン型回折格子であり、Ｇ２は、ラインがＧ１のラインと平行になっている検出器の前に配置される、
ことを特徴とするシステム。
【請求項２２】
請求項１７から２１までのいずれか１項に記載のシステムにおいて、
平行ビーム、準平行、扇形ビームまたは円錐形ビームモードになるように動作が選択され、Ｇ０，Ｇ１およびＧ２は相応に平坦、円筒形または円錐形の形状をそれぞれ有する、
ことを特徴とするシステム。
【請求項２３】
請求項１７から２２までのいずれか１項に記載のシステムにおいて、
２次元回折格子による全視野モードかまたは１次元回折格子による走査モードになるように動作が選択される、
ことを特徴とするシステム。
【請求項２４】
請求項１８記載のシステム、又は請求項１８を引用する請求項１９から２３までのいずれか1項記載のシステムにおいて、
近接方式動作に対し、前記回折格子間の間隔を当該近接方式内で自由に選択し、かつ、
Talbot方式に対し、
【数５】
0005818909.tif 000032
ｌ＝１，２，３，…，Ｄnは平行Ｘ線ビームが使用される場合の奇数の分数Talbot距離であり、これに対してＤn,sphは扇形または円錐形のＸ線ビームが使用される場合の奇数の分数Talbot距離であり、Ｌは線源とＧ１との間の間隔である、
ことを特徴とするシステム。
【請求項２５】
請求項１７から２４までのいずれか１項に記載のシステムにおいて、
１つの回折格子Ｇ０，Ｇ１またはＧ２を、別の回折格子に対して機械式にシフトすることによって位相ステッピングが行われる、
ことを特徴とするシステム。
【請求項２６】
請求項１８記載のシステム、又は請求項１８を引用する請求項１９から２５までのいずれか1項記載のシステムにおいて、
前記回折格子構造はプレーナ技術によって作製されている、
ことを特徴とするシステム。
【請求項２７】
請求項２５に記載のシステムにおいて、
前記演算子Ｄxは、前記回折格子の前記位相ステッピング方向における任意の次数の微分演算子である、
ことを特徴とするシステム。
【請求項２８】
請求項２５または２７に記載のシステムにおいて、
前記変換演算子Ｔは、前記回折格子の前記ステッピング方向に対して垂直な方向における任意の次数の微分演算子である、
ことを特徴とするシステム。
【請求項２９】
請求項１７から２８までのいずれか１項に記載のシステムにおいて、
前記重み付け演算子Ｗは、各ピクセルにおけるＤＰＣ画像の標準偏差の逆数１／σDPCを含む対角重み付け行列である、
ことを特徴とするシステム。
【請求項３０】
請求項１７から２９までのいずれか１項に記載のシステムにおいて、
前記制約付き最適化問題は、当該制約付き最適化問題を２次錐計画（ＳＯＣＰ）に計算しなおすことによって解かれる、
ことを特徴とするシステム。
【請求項３１】
請求項１７から３０までのいずれか１項に記載のシステムにおいて、
前記制約付き最適化問題は、例えば
【数６】
0005818909.tif 000033
にしたがって制約のない形にされることによって解かれる、
ことを特徴とするシステム。
【請求項３２】
請求項３１に記載のシステムにおいて、
勾配降下法または（非線形）共役勾配アルゴリズムによって前記制約のない最適化問題が解かれる、
ことを特徴とするシステム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、制約付き最適化を使用して微分画像から積分画像を再生する方法およびシステムに関する。
【０００２】
慣用の可視光光学系とは異なり、Ｘ線光学系における屈折率は１に極めて近くかつ１よりも小さいことはよく知られている。１次近似において、媒体における小さくかつ無視可能な非等方性に対し、組織の光学特性を特徴付ける屈折率は（Ｘ線吸収を含めて）その複素数の形態で、すなわちｎ＝１−δ−ｉβで表すことができる。ただしδは、位相シフト特性を表す屈折率の実部の減少であるのに対し、虚部βは、サンプルの吸収特性を表す。慣用の吸収ベースのラジオグラフィでは、Ｘ線位相シフト情報は一般的に、画像再構成に直接利用されない。しかしながら10keVを上回る光子エネルギにおいて、また（低原子番号素子から構成される）軽材料に対しては、δは一般的にβよりも３桁大きいため、上記の位相シフト項は、減衰項よりも顕著に作用する。この結果、位相コントラストモダリティは、慣用の吸収ベースのイメージングに比べて格段に大きい画像コントラストを形成できるのである。
【０００３】
さらに吸収エッジから遠ざかると、δはＸ線エネルギの二乗に反比例するのに対してβは、エネルギの４乗で減少する。このメカニズムの重要な結果は、吸収に比べてはるかに低い線量堆積で位相信号が得られることであり、生体試料または生体系などにおいて放射ダメージを考慮しなければならない場合には極めて重要なことである。
【０００４】
上記の位相信号を記録するためにいくつかのアプローチが開発されている。これらのアプローチは、（結晶による）干渉法、位相伝搬法、アナライザ結晶に基づく技術、またはｘ線回折格子に分類可能である。ここで説明する発明は、殊に後者の技術に関連するものである。
【０００５】
回折格子ベースのｘ線イメージングセットアップにおいては実質的に、対象体におけるｘ線の偏向が検出される。このような偏向は、対象体における位相シフトグラジエントでの反射によって発生して微分位相コントラスト（ＤＰＣ differential phase contrast）になるか、またはサンプルにおける不均一性での散乱によっていわゆる暗視野法（ＤＦＩ dark field image）になり得る。ＤＰＣ画像信号は、イメージ処理ルーチンによって位相コントラスト（ＰＣ phase contrast）画像を得るために使用可能である。
【図面の簡単な説明】
【０００６】
【図１】制約付き最適化を使用して微分画像から積分画像を再生するシステムを示す図である。
【図２】modified Shepp-Loganファントムによって形成された、雑音の多いＤＰＣ画像の積分を示す図である。
【０００７】
図１に示したように、２つ回折格子（Ｇ１およびＧ２）または３つの回折格子（Ｇ０，Ｇ１およびＧ２）を有するセットアップを適用してｘ線の偏向を記録することができる。２回折格子セットアップの場合、ｘ線源はその空間コヒーレンスについて所定の要求を満たす必要があるが、３回折格子セットアップでは空間コヒーレンスは不要である。Ｇ０回折格子が必要になるのは、Ｘ線源のサイズがp2＊ｌ／ｄよりも大きい場合であり、ただしp2は、Ｇ２の周期であり、ｌは、Ｘ線源とＧ１との間の間隔であり、ｄは、Ｇ１とＧ２との間の間隔である。したがって３回折格子セットアップは、インコヒーレントｘ線源に使用するのに適しており、殊にｘ線管に使用するのに適している。
【０００８】
慣用の減衰コントラスト（ＡＣ attenuation contrast）と、ＤＰＣおよびＤＦＩコントラストとを分離するために位相ステッピングアプローチを適用する。複数の回折格子のうちの１つは、入射ビームに対して横方向にずらされ、その間に複数の画像が取得される。検出器面内の各ピクセルにおける輝度信号は、変位に応じて振動する。この振動の平均値は、上記の減衰コントラスト（ＡＣ）を表す。この振動の位相は、波面位相プロフィールの１次微分に直接結びつけることができ、したがってＤＰＣ信号に結びつけることができる。上記の振動の振幅は、対象体におけるｘ線の散乱に依存し、したがってＤＦＩ信号を生ずる。
【０００９】
（２つまたは３つの）回折格子に対して、いくつかのアプローチが提案されて適用されている。回折格子Ｇ０（必要な場合）は、Ｘ線源に最も近い回折格子である。この回折格子は一般的に、周期p0の吸収ラインの透過回折格子からなる。これは、同じ周期のラインだけから放射線を発する線源によって置き換えることができる。回折格子Ｇ１は、ｘ線源のさらに下流に配置される。この回折格子は、周期p1のラインから構成される。回折格子Ｇ２は、このセットアップの最も下流にある回折格子である。この回折格子は一般的に、周期p2の吸収ラインの透過回折格子からなる。この回折格子は択一的に、回折格子に類似した感度を有する同じ周期の検出システムによって置換することが可能である。
【００１０】
セットアップの２つの方式は、いわゆる「近接方式」と「Talbot方式」とに区別することができる。「近接方式」では回折格子周期p、回折格子間隔ｄおよびｘ線波長λは、回折作用が無視できるように選択される。この場合、すべての回折格子を吸収ラインから構成しなければならない。「Talbot方式」では、回折構造における回折が重要である。正確な判定基準は、回折構造のデューティサイクルに依存し、かつ、この回折格子が吸収型であるか位相シフト型であるかに依存するためこれらの２つの方式をはっきりと区別することは容易ではない。例えば、吸収ラインおよび0.5のデューティサイクルを有する回折格子に対し、「近接方式」に対する条件はｄ≧p2/2λである。ここでＧ１は、吸収型であるか、または有利には位相シフト型である回折格子ラインから構成すべきである。いくつかの位相シフト量が考えられ、有利にはπ／２またはその数倍である。回折格子の周期は、回折光子間の相対的な間隔に整合しなければならない。「Talbot方式」のセットアップの場合、良好なコントラストを得るためにTalbot効果を考慮しなければならない。回折格子周期および間隔に対する式は、従来技術において公知である。
【００１１】
サンプルは多くの場合にＧ１とＧ０と間に（または２回折格子セットアップの場合にはＧ１の上流に）配置されるが、Ｇ１とＧ２との間にサンプルを配置することも有利になり得る。
【００１２】
本発明は、上述したすべてのケースに関連しており、すなわち、２回折格子および３回折格子のケース、「近接方式」のケースおよび「Talbot方式」に関連しており、またＧ１の上流または下流に配置されるサンプルに対するものである。
【００１３】
さらにここに示した発明はまた、スキャニングベースのシステムと組み合わせて、平坦な回折格子を使用した平行および準平行な幾何学形状に対し、または円筒形または球形に湾曲した回折格子を使用した小型の扇型ビームまたは円錐型ビームに対しても動作する。
【００１４】
上述のシステムを用いたイメージング実験は、Ｘ線屈折の影響を受け易い。（光軸に対する）屈折角αと、上記のシステムによって測定した微分位相シフトとの間の関係は、
【数１】
0005818909.tif 000002
によって得られ、ただしｄは回折格子間の間隔であり、p2はＧ２の周期である。ｘ軸は、位相ステッピング方向に相応する。さらにサンプルの後における波の位相プロフィールφ(ｘ)とα(ｘ)との間の関係は、［１２］によって得られ、
【数２】
0005818909.tif 000003
である。
式（２）を式（１）に代入することにより、
【数３】
0005818909.tif 000004
が得られる。
したがってΔφは、φ(ｘ)の１次微分に比例する。この式によれば、φ(ｘ)の再構成にはｘ方向におけるＤＰＣ測定値を積分が必要であり、
【数４】
0005818909.tif 000005
である。
【００１５】
一般的に雑音の多い信号を積分することにより、雑音エラーおよび雑音分散が累積される。その結果、積分の方向に振幅が増大する水平方向のストライプ状のアーチファクトが発生する。図２には、modified Shepp-Loganファントムによって形成された、雑音の多いＤＰＣ画像の積分が示されており、図２ｃではストライプ状のアーチファクトを明瞭に視認することができる。
【００１６】
したがって本発明の目的は、積分された位相コントラスト画像の品質を改善するシステムおよび方法を提供することである。
【００１７】
この目的は、本発明により、独立請求項１および１９に記載された特徴的構成によって達成される。本発明の好適な実施形態は、請求項２から１８、および、２０から３６にそれぞれ示されている。
【００１８】
制約付き最適化を使用することによって微分画像から積分画像を再生する本発明によるシステムおよび本発明による方法には、
ａ）
【数５】
0005818909.tif 000006
と記される、制約
0005818909.tif 000007
の条件下でコスト関数
0005818909.tif 000008
を最小化する画像を識別するための手段が含まれており、
ｂ）  ただしｆは画像ベクトル、Ｔは変換演算子行列、Ｄxは微分演算子行列、Ｗは、雑音標準偏差の逆数を含む対角行列、φは測定した画像ベクトル、εは、雑音パワーについての境界である。
0005818909.tif 000009
は、ベクトルのpノルムを示しており、ただしpは正数、有利には区間［０，２］内の整数である。
【００１９】
本発明の有利な実施形態は、サンプルから量的ｘ線画像を取得するためにｘ線用の装置、殊に硬ｘ線用の装置から、微分画像が取得される場合に得られ、ここでこの装置は、
ａ．Ｘ線源と、
ｂ．Ｇ０，Ｇ１およびＧ２と称される３つの回折格子、または、Ｇ１およびＧ２と称される少なくとも２つの回折格子と、
ｃ．複数の個別ピクセルを有しかつ空間変調検出感度を備えた位置感応型検出器と、
ｄ．この検出器の画像を記録するための手段と、
ｅ．個別ピクセル毎に、吸収が主であるピクセル、および／または、微分位相コントラストが主であるピクセル、および／または、ｘ線散乱が主であるピクセルとして対象体の特性を識別するため、一連の画像においてピクセル毎の強度を評価するための手段と、
ｆ．ここで上記の一連の画像は、０からπまたは２πで連続的または段階的に上記のサンプルを回転させるか、またはこのサンプルに対して上記の装置および線源を回転させることによって収集される。
【００２０】
上記のシステムおよび方法は有利にも、いわゆる「近接方式」または「Talbot方式」のいずれかで動作させることが可能である。
【００２１】
回折格子Ｇ１の有利な実施形態では、Ｇ１は、吸収回折格子または位相回折格子であるライン型回折格子であり、位相回折格子は、低吸収型回折格子であるが大きなＸ線位相シフトを形成し、このＸ線位相シフトは有利にはπ／２またはその倍数である。
【００２２】
したがって回折格子Ｇ２は、Ｇ１の自己画像の周期と同じ周期を有する高Ｘ線吸収コントラストを有するライン型回折格子として実現することができ、Ｇ２は、ラインがＧ１のラインと平行になっている検出器の前に配置される。
【００２３】
上記のシステムおよび方法では、動作に所定の自由が許容され、動作は、平行ビーム、準平行、扇形ビームまたは円錐形ビームモードになるように選択でき、またＧ０，Ｇ１およびＧ２は相応に平坦、円筒形または円錐形の形状をそれぞれ有する。したがって上記の動作は、２次元回折格子による全視野モードかまたは１次元回折格子による走査モードになるように選択可能である。
【００２４】
上記のシステムおよび方法をセットアップするため、２つのタイプの動作をつぎのように定めることができる。すなわちａ）近接方式動作に対し、上記の回折格子間の間隔は、この近接方式内では自由に選択され、かつｂ）Talbot方式に対し
【数６】
0005818909.tif 000010
ｌ＝１，３，５，…であり、Ｄnは平行Ｘ線ビームが使用される場合の奇数の分数Talbot距離であり、これに対してＤn,sphは扇形または円錐形のＸ線ビームが使用される場合の奇数の分数のTalbot距離であり、Ｌは線源とＧ１との間の間隔である。
【００２５】
微分位相コントラストアプローチの利点のフルに活用するため、上記のシステムおよび方法は、上記の位相ステッピングが、１つの回折格子Ｇ０，Ｇ１またはＧ２を他の回折格子に対して機械式にシフトさせることによって行われるようにして実行される。
【００２６】
回折格子構造に関しては、この回折格子構造は有利にはプレーナ技術によって作製される。
【００２７】
ＣＴ，ＭＲＩまたは超音波装置から得られる医用画像などのような画像において微分位相情報を分離するため、欧州特許明細書EP 10167569.2にしたがってこの微分位相情報を取得することができる。この欧州特許明細書は、参照によってこの明細書に含まれるものとする。
【００２８】
本発明の別の有利な実施形態において、Ｇ１とＧ２との間の位相関係は、上記の強度曲線を１次のテーラー級数によって展開することのできる値に正確に一致することができ、かつ、上記の微分位相情報は有利には国際出願明細書PCT/EP2010/051291（WO 2010/089319）にしたがって取得することができる。この国際出願明細書は参照によってこの明細書に含まれるものとする。
【００２９】
数学的なコスト関数を解く間、演算子Ｄxは、回折格子の位相ステッピング方向における任意の次数の微分演算子として設計可能である。これに相応して上記の変換演算子Ｔは、回折格子のステッピング方向に対して垂直な方向における任意の次数の微分演算子になるように選択可能である。さらに上記の重み付け演算子Ｗは、各ピクセルにおけるＤＰＣ画像の標準偏差の逆数１／σDPCを含む対角重み付け行列となるように選択可能である。
【００３０】
有利には上記の制約付き最適化問題は、これを２次錐計画（ＳＯＣＰ）に計算しなおすことによって解くことができる。択一的にはこの制約付き最適化問題は、例えばこれを
【数７】
0005818909.tif 000011
にしたがって制約のない形にすることによって解くことができる。
【００３１】
上記の制約のない最適化問題は択一的に、勾配降下法または（非線形）共役勾配アルゴリズムによって解くことができる。
【００３２】
この結果、ここで提案する新しい発明により、雑音の多いＤＰＣ画像を直接積分することによって生じるストライプ状のアーチファクトの問題が扱われる。基本的なアイデアは、制約付きの最適化問題を解くことにより、垂直の画像方向における変動を抑止することである。上記の測定データとの一貫性を維持しながら、コスト関数を最小化することによって上記の位相画像を取り出す。雑音の多いＤＰＣ測定値に適用すれば、これにより、上記の画像において上記の積分によって強制的に一層小さな変動が形成されるようになり、したがって画像品質が改善されるのである。
本発明の好適な実施例を以下、詳しく説明する。
【００３３】
以下では、医用画像、または構造解析または材料試験などから得られる画像などの画像をベクトルとして表す。画像ベクトルｆ(ｉ)は、画像Ｉ(ｘ,ｙ)
【数８】
0005818909.tif 000012
において列毎にピクセル値を抽出することによって得られ、ここでｎx×ｎyは画像サイズである。画像ベクトルのサイズは１×ｎ、ただしｎ＝ｎx・ｎyである。画像変換は、画像ベクトルに適用可能な演算子（行列）によって表される。演算子行列のサイズはｍ×ｎであり、ｍは変換されるベクトルのサイズである。多くの場合、ｍ＝ｎが成り立つ。
【００３４】
上記の新しい方法は、標準回帰モデルに基づいており、ＤＰＣ画像の測定値は、
φ＝Ｄx・φ＋ｗ                                    （６）
によって得られる。
【００３５】
Ｄxはφの微分測定をモデル化する測定値演算子であり、またｗは、ＤＰＣ画像における雑音をモデル化するランダムベクトルである。Ｄxは、ｘ方向におけるφの有限差分変換によって実現可能である。
【００３６】
雑音モデルについては、ＤＰＣ画像の１つのピクセルにおける雑音分散は、
【数９】
0005818909.tif 000013
によって得られる。ただしＮはフォトンの個数であり、Ｖはこのピクセルにおける平均しま可視度である。ＮおよびＶは、上記の測定のＡＣおよびＤＦＩ画像を使用して計算することができる。
【００３７】
ＤＰＣ測定値から位相を取り出すため、制約付き最適化問題を定める。すなわち、
【数１０】
0005818909.tif 000014
である。
【００３８】
この最適化問題においてｆは画像ベクトルであり、Ｔは変換演算子であり、φは測定したＤＰＣ画像ベクトルであり、またεは上記の雑音パワーの境界である。Ｗは、この雑音モデルを考慮するための対角重み付け行列であり、ただしｗi,i＝１／σDPC,iである。
0005818909.tif 000015
は、ベクトルのpノルムを示しており、これは
【数１１】
0005818909.tif 000016
と定義される。
【００３９】
（８）に示した上記の問題は、データ整合性の制約
【数１２】
0005818909.tif 000017
を満たす任意のｆに対してベクトルＴｆの最小pノルムを求めようとするものである。
【００４０】
（８）の最小化項において上記の画像は、（例えばフリーエ変換、ウェーブレット、有限差分など）行列Ｔによって表される任意の線形変換領域に変換可能である。 これにより、データ制約の利用が極めてフレキシブルになる。ＤＰＣ測定の場合、有限差分変換
Ｔ＝Ｄy      （１０）
が有利である。それは、上記の積分画像は、大きな強度変動（水平方向のストライプ）により、位相ステッピングに対して垂直な方向に歪まされるからである。
【００４１】
問題（８）を解くのには多種多様な可能性がある。凸の場合（p≧１）、これは計算しなおすことができ、２次錐計画として解くことができる。択一的な選択肢は、ラグランジュ乗算子を使用して
【数１３】
0005818909.tif 000018
または一般的に任意の量のデータ制約
【数１４】
0005818909.tif 000019
を用いて（８）を制約のない形に計算し直すことである。
【００４２】
問題のこのような定式化は正則化として知られており、これは主に不良設定逆問題に適用されている。ラグランジュ乗算子（または正則化パラメタ）λiは、データ整合項
0005818909.tif 000020
と比較した正則化項
0005818909.tif 000021
の重み付けを制御する。
【００４３】
問題（１２）は、任意個の正則化項を許容する強力な特性を有し、このことは、上記の対象体について一層アプリオリな知識が利用可能な場合に殊に有益である。
【００４４】
上記の正則化項におけるノルムパラメタｐの選択は、使用される変換演算子Ｔに依存する。一般的な選択はｐ＝２である。それは、この場合に問題（１１）／（１２）は線形になり、明示的な解が存在するからである。Ｔ＝Ｄyの場合、ｌ2ノルムは不鮮明に結び付くことがあり、したがって画像の解像度が低減される。その一方でｌiノルム（ｐ＝１）の最小化は、画像においてエッジを保持することが知られており、したがってｌ2ノルムよりも一層良好な画像の解像度が得られることが期待される［１］。ｌ1ノルム最小化により、明示的な解が存在しない非線形最適化問題が得られる。
【００４５】
ここではこの非線形最適化問題を解くために反復式のアルゴリズムが使用される。これは、非線形共役勾配（ＮＬＣＧ non-linear Conjugate Gradients）法をベースにしており、これは、大規模線形システムの逆問題に対して高速に収束するという特徴を有する［２］。
【００４６】
参考文献
［１］L. Rudin，S. OsherおよびE. Fatemi，"Nonlinear total variation based noise removal algorithms"，Phys. D Nonlinear Phenom. 60, 259-268 (1992)
［２］J. NocedalおよびS. Wright, Numerical optimization (Springer verlag, 1999)
【図１】
0005818909.tif 000034
【図２】
0005818909.tif 000035