(57)【要約】【課題】効率的なパラメータの推定探索を実現することを課題とする。【解決手段】Ｘ線ＣＴ装置と、投影ベース法を用いることにより画像を再構成する演算装置と、を有し、演算装置は、各基底物質厚さと、基底物質厚さの尤度とを座標軸として有する座標空間を設定し、Ｘ線減弱応答に基づいて、座標空間に設定された推定厚さ入力値３０１を起点とし、尤度等高線３００の尾根方向と直角の方向において、尤度が最も大きくなる第１推定厚さ３０８を探索する第１サーチを行い、推定厚さ入力値３０１からシフトした位置に設定されたシフト始点３０４を起点として、尤度が最も大きくなる第２推定厚さ３０９を探索する第２サーチを行い、第１推定厚さ３０８及び第２推定厚さ３０９を結ぶライン上において、尤度が最も高くなる最尤推定量を探索する第３サーチを行うことで、推定厚さ出力値３１２を探索する第３サーチを行うことを特徴とする。【選択図】図４
2019033781.tif 000002
【特許請求の範囲】
【請求項１】
投影ベース法を用いることにより、Ｎ個（Ｎ≧２）の第１のパラメータについてのＸ線エネルギ分布を有する被撮像対象のＸ線減弱応答に含まれるＭ個（Ｍ≦Ｎ）の第２のパラメータに対する所定の値の推定値に基づいて画像を再構成する演算部が、
前記Ｍ個の第２のパラメータと、前記Ｍ個の第２のパラメータの尤度とを座標軸として有する座標空間を設定し、
Ｘ線減弱応答に基づいて、前記座標空間に所定の起点を設定し、前記起点から前記座標空間における所定の方向において、前記尤度が最も大きくなる第１の値を探索する第１のサーチを行う第１のステップと、
前記第１のサーチとは起点及び方向の少なくとも一方が異なり、前記尤度が最も大きくなる第２の値を探索する第２のサーチを行う第２のステップと、
前記第１の値及び前記第２の値を結ぶライン上において、前記尤度が最も大きくなる最尤推定量を探索する第３のサーチを行い、前記最尤推定量を前記Ｍ個の第２のパラメータの推定値とする第３のステップと、
を行うことを特徴とするパラメータ推定方法。
【請求項２】
前記尤度の分布における尾根のラインが、前記座標空間における前記Ｍ個の第２のパラメータを示す座標軸に対してなす角度が入力部を介して入力され、
前記第１のステップ及び前記第２のステップにおける前記第１のサーチ及び前記第２のサーチは、前記尾根のラインに対して垂直方向に行われる
ことを特徴とする請求項１に記載のパラメータ推定方法。
【請求項３】
前記第２のサーチの起点は、前記座標空間上において、前記第１のサーチの起点から、前記尾根のラインに対して平行な方向に所定距離離れた箇所に設定される
ことを特徴とする請求項２に記載のパラメータ推定方法。
【請求項４】
Ｎ個（Ｎ≧２）の第１のパラメータについてのＸ線エネルギ分布を有する被撮像対象のＸ線減弱応答を得ることで、Ｎ個の投影データを取得するＸ線投影データ取得部と、
投影ベース法を用いることにより、前記Ｘ線減弱応答に含まれるＭ個（Ｍ≦Ｎ）の第２のパラメータに対する所定の値の推定値に基づいて画像を再構成する演算部と、
を有し、
前記演算部は、
前記Ｍ個の第２のパラメータと、前記Ｍ個の第２のパラメータの尤度とを座標軸として有する座標空間を設定し、
Ｘ線減弱応答に基づいて、前記座標空間に所定の起点を設定し、前記起点から前記座標空間における所定の方向において、前記尤度が最も大きくなる第１の値を探索する第１のサーチを行う第１のステップと、
前記第１のサーチとは起点及び方向の少なくとも一方が異なり、前記尤度が最も大きくなる第２の値を探索する第２のサーチを行う第２のステップと、
前記第１の値及び前記第２の値を結ぶライン上において、前記尤度が最も大きくなる最尤推定量を探索する第３のサーチを行い、前記最尤推定量を前記Ｍ個の第２のパラメータの推定値とする第３のステップと、
を行うことを特徴とするＸ線ＣＴシステム。
【請求項５】
前記尤度の分布における尾根のラインが、前記座標空間における前記Ｍ個の第２のパラメータを示す座標軸に対してなす角度が入力部を介して入力され、
前記第１のステップ及び前記第２のステップにおける前記第１のサーチ及び前記第２のサーチは、前記尾根のラインに対して垂直方向に行われる
ことを特徴とする請求項４に記載のＸ線ＣＴシステム。
【請求項６】
前記第２のサーチの起点は、前記座標空間上において、前記第１のサーチの起点から、前記尾根のラインに対して平行な方向に所定距離離れた箇所に設定される
ことを特徴とする請求項５に記載のＸ線ＣＴシステム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、パラメータ推定方法及びＸ線ＣＴシステムの技術に関する。
【背景技術】
【０００２】
Ｘ線ＣＴ（Computer Tomography）撮像装置（以下、Ｘ線ＣＴ装置と適宜記載する）はエネルギ情報を得られない電流モードのＸ線検出器を用いて、Ｘ線管による連続（非単色）エネルギ分布（Ｘ線エネルギ分布：以下、スペクトルと称する）を有するＸ線光子群を検出する構成が一般的である。
ところで、Ｘ線減弱係数は原子番号ごとに固有のエネルギ依存性を有している。従って、異なるスペクトルによる複数の減弱係数から原子番号に関する情報を得ることができる。しかし、電流モードのＸ線検出器では、エネルギ情報を得られないため原子番号に関する情報を取得できない。
【０００３】
複数のスペクトルを有するＸ線群による情報を有効に利用するための手法として大きく以下の二つの手法が挙げられる。一つ目はＤｕａｌ  Ｅｎｅｒｇｙ  ＣＴである。この技術では、検出器は電流モードのままで、２種のＸ線管電圧による二つの連続スペクトルを用いるものである。二つ目はｐｈｏｔｏｎ  ｃｏｕｎｔｉｎｇ  ＣＴ  （ＰＣＣＴ）や、ｓｐｅｃｔｒａｌ  ＣＴ等と呼ばれる技術である。この技術は、エネルギ情報が得られるパルスモード検出器を用いるものである。
【０００４】
Ｘ線ＣＴ装置は物質のＸ線を遮る能力（減弱係数）の違いを見るものである。しかしながら、Ｘ線減弱係数のエネルギに対する依存性は元素（原子番号）によって異なる。これを利用すればＮ種のエネルギ情報が得られる場合に、異なる実効原子番号の任意の物質Ｍ個（Ｍ≦Ｎ）を基底物質とした物質分解が可能となる（例えば、特許文献１参照）。この物質分解を投影から画像に戻す（画像再構成時）前に行うか、後で行うかにより、投影ベース法（pre-reconstruction method）、画像ベース法（post-reconstruction method）の２つの画像再構成法が存在する（例えば、非特許文献１参照）。臨床ＣＴでの高フラックス条件のためＰＣＣＴは検出素子を従来のＸ線ＣＴよりも微細化している。その複雑なスペクトル応答を取り扱うため、投影ベース法が一般的に用いられている。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】米国特許出願公開第２００６／０１０９９４９号明細書
【非特許文献】
【０００６】
【非特許文献１】ASSIS、 V. C.; SALVADEO、 D. H. P.; MASCARENHAS、 N. D.; LEVADA、 A. L. M. Double Noise Filtering in CT: Pre- and Post-Reconstruction. In: CONFERENCE ON GRAPHICS、 PATTERNS AND IMAGES、 28. (SIBGRAPI)、 2015、 Salvador. Proceedings. Los Alamitos: IEEE Computer Society´s Conference Publishing Services、 2015. On-line.
【発明の概要】
【発明が解決しようとする課題】
【０００７】
一般的に、投影ベース法による画像再構成方法では、実測されたスペクトル（実測スペクトル）と基底物質厚さ推定値に対応する想定スペクトル（平均スペクトル）との２者から得られる尤度に関して基底物質厚さの最尤推定が行われる。なお、平均スペクトルの「平均」とは、基底物質厚さに対応した（＝過去に調べた）カウントの期待値のことである。また、想定スペクトルとは、ある基底物質厚さに対し過去のデータを基に、これが平均値であると想定したスペクトルのことである。基底物質厚さと、想定スペクトル（平均スペクトル）とを基に、基底物質厚さの最尤推定量を算出することができる。
しかしながら、想定スペクトル（平均スペクトル）を基に基底物質厚さの最尤推定量を算出することは、一般的に高い計算負荷を伴う（詳細は後記）。
【０００８】
また、最尤推定のサーチ手法は複数存在するが、一般に、このような最尤推定では、繰り返し計算を行うことで、最尤推定量が求められる。
しかしながら、このサーチ手法によっては、繰り返し計算数が多くなり、計算負荷が増大するという課題がある。
【０００９】
このような背景に鑑みて本発明がなされたのであり、本発明は、効率的なパラメータの推定探索を実現することを課題とする。
【課題を解決するための手段】
【００１０】
前記した課題を解決するため、本発明は、投影ベース法を用いることにより、Ｎ個（Ｎ≧２）の第１のパラメータについてのＸ線エネルギ分布を有する被撮像対象のＸ線減弱応答に含まれるＭ個（Ｍ≦Ｎ）の第２のパラメータに対する所定の値の推定値に基づいて画像を再構成する演算部が、前記Ｍ個の第２のパラメータと、前記Ｍ個の第２のパラメータの尤度とを座標軸として有する座標空間を設定し、Ｘ線減弱応答に基づいて、前記座標空間に所定の起点を設定し、前記起点から前記座標空間における所定の方向において、前記尤度が最も大きくなる第１の値を探索する第１のサーチを行う第１のステップと、前記第１のサーチとは起点及び方向の少なくとも一方が異なり、前記尤度が最も大きくなる第２の値を探索する第２のサーチを行う第２のステップと、前記第１の値及び前記第２の値を結ぶライン上において、前記尤度が最も大きくなる最尤推定量を探索する第３のサーチを行い、前記最尤推定量を前記Ｍ個の第２のパラメータの推定値とする第３のステップと、を行うことを特徴とする。
その他の解決手段については実施形態中において適宜記載する。
【発明の効果】
【００１１】
本発明によれば、効率的なパラメータの推定探索を実現することができる。
【図面の簡単な説明】
【００１２】
【図１】第１実施形態に係るＸ線ＣＴシステム１０の概要を示す図である。
【図２】本実施形態で用いられる演算装置２の構成例を示す機能ブロック図である。
【図３】本実施形態で行われるＣＴ画像再構成処理を示すフローチャートである。
【図４】本実施形態で用いられる最適化処理の手法を示す図である。
【図５】本実施形態で用いられる最適化手法の手順を示すフローチャートである。
【発明を実施するための形態】
【００１３】
次に、本発明を実施するための形態（「実施形態」という）について、適宜図面を参照しながら詳細に説明する。
【００１４】
（Ｘ線ＣＴシステム１０）
図１は、第１実施形態に係るＸ線ＣＴシステム１０の概要を示す図である。
Ｘ線ＣＴシステム１０は、Ｘ線ＣＴ装置１、ＰＣ（Personal Computer）等の演算装置２、入力装置３、表示装置４等を有する。また、Ｘ線ＣＴ装置１は、ガントリ１１、Ｘ線管１２、寝台１３、Ｘ線検出器１５を有する。
被撮像対象である被検体Ｐは寝台１３上で静止している。
Ｘ線管１２とＸ線検出器１５は寝台１３上で静止している被検体Ｐを挟んで、互いに対向して配置されている。そして、入力装置３から操作情報が入力されると、ガントリ１１の回転駆動部（不図示）により、Ｘ線管１２とＸ線検出器１５の組が被検体Ｐの体軸周りを回転する。このとき、Ｘ線管１２から照射され、被検体Ｐによる減弱を受けたＸ線群がＸ線検出器１５によって計測される。
【００１５】
ここで、Ｘ線検出器１５はエネルギ情報を得ることが可能なパルスモードで動作するものとする。そして、Ｘ線検出器１５によって、Ｎ種のエネルギに関するカウント（すなわち、スペクトル）情報が取得される。また、取得されたサイノグラムの投影経路（回転角度位置、検出器位置）ごとに検出されたスペクトル情報に対し、演算装置２が画像再構成処理を行う。そして、表示装置４が演算結果である断層像を表示する。
【００１６】
（演算装置２）
図２は、本実施形態で用いられる演算装置２の構成例を示す機能ブロック図である。
演算装置２は、メモリ２１、ＣＰＵ（Central Processing Unit）２２、ＨＤ（Hard Disk）等の記憶装置２３、送受信装置２４を有する。
メモリ２１には、記憶装置２３からプログラムがロードされる。そして、メモリ２１にロードされたプログラムがＣＰＵ２２によって実行されることで、投影ベース処理部２１１、画像再構成処理部２１３、表示処理部２１４が具現化している。
【００１７】
投影ベース処理部２１１は、投影ベース法による物質分解処理を行う。また、投影ベース処理部２１１は、最適化処理部２１２を有している。最適化処理部２１２は、後記する最適化処理を行う。
画像再構成処理部２１３は、投影ベース処理部２１１によって物質分解処理が行われたデータに対して画像再構成処理を行う。
表示処理部２１４は、画像再構成処理の結果を表示装置４（図１参照）に表示する。
【００１８】
記憶装置２３には、リッジ構造代表角度２３１及び実測スペクトル投影データ２３２が格納されている。リッジ代表角度については後記する。実測スペクトル投影データ２３２は、今回、Ｘ線ＣＴ装置１において撮像されたもの（スペクトル）である。
送受信装置２４は、Ｘ線ＣＴ装置１等との送受信を行う。
【００１９】
（全体処理）
図３は本実施形態で行われるＣＴ画像再構成処理を示すフローチャートである。適宜、図１及び図２を参照する。なお、図３の処理において、最適化処理（Ｓ２）以外の処理は、一般的に行われている処理であるため、各処理の詳細な説明は省略する。
まず、Ｘ線ＣＴ装置１が被検体Ｐに対する撮像を行う（Ｓ１）。この撮像の結果、サイノグラムの各投影経路における実測スペクトル投影データ２３２が得られる。
【００２０】
次に、投影ベース処理部２１１が投影ベース処理を行う（Ｓ２）。なお、投影ベース処理中において、最適化処理部２１２が最適化処理を行う。最適化処理については図４及び図５を参照して後で説明する。
【００２１】
理想的には、ある投影経路のあるエネルギウィンドウには、その投影経路とエネルギで被検体Ｐに入射したＸ線が減弱したもののみがカウントされることが望ましい。しかし、現実には被検体Ｐでの散乱、Ｘ線検出器１５内部での散乱等により、入射時と異なる投影経路、エネルギのＸ線によるカウントが重畳しうる。実測スペクトル投影データ２３２はこのような複雑なスペクトル応答の結果として得られるものである。しかし、この応答を基底物質として選択した任意の物質厚さに対して事前に把握できれば、該基底物質厚さ（パラメータ）の組み合わせに対して実測スペクトル投影データ２３２が実現する確率を計算することができる。なお、この確率は、基底物質をパラメータとした尤度とみなせる。すなわち、該基底物質厚さ組み合わせに対して、今回測定された実測スペクトル投影データ２３２が実現する確率と、基底物質をパラメータとした尤度とは数値として同一である。
【００２２】
基底物質厚さ組み合わせと、想定(平均)スペクトルとの対応の既知情報としてはＸ線挙動の計算機シミュレーション結果でも、Ｘ線検出器１５による実測結果でも、その複合でもよい。このとき、尤度を最大にするような基底物質厚さの組み合わせ（基底物質推定厚さ投影データ２４１）を各投影経路で求める方法が投影ベース法である。投影ベース法は、一般的には特別な初期値を必要とせず、適切なデフォルト値から開始する。例えば、画素値が一様な画像等が初期値として用いられる。しかし、なんらかの別方法により解に近い値が既知であれば、これを推定厚さ初期値投影データ２４３として用いてもよい。推定厚さ初期値投影データ２４３は、例えば、単なる円柱を撮像した場合に得られると推定される投影データ等を基に作成される。
【００２３】
投影ベース処理後、画像再構成処理部２１３が投影ベース処理の結果出力される基底物質推定厚さ投影データ（推定パラメータ）２４１を用いて画像再構成処理を行う（Ｓ３）。
画像再構成処理部２１３は、基底物質推定厚さ投影データ２４１をＦＢＰ（Filtered back projection）等の画像再構成法により基底物質推定濃度画像２４２に変換する。
【００２４】
（最適化処理）
図４は、本実施形態で用いられる最適化処理の手法を示す図である。また、図５は、本実施形態で用いられる最適化手法の手順を示すフローチャートである。ここでは、図４をメインとして、適宜図５を参照しつつ説明を行う。なお、ここでの説明において、ステップ番号は図５のステップ番号である。
ここでは単純化のため、各投影経路で独立に厚さ推定を行う場合を考える。
尤度等高線（評価関数）３００は基底物質が２個である場合における基底物質厚さの尤度分布を示している。なお、本実施形態では、基底物質２個としているが、３個以上としてもよい。この場合、図４の座標空間において基底物質の数に応じた座標軸が設定される。また、図４の座標空間において、紙面垂直方向は尤度の座標軸である。
【００２５】
ここで投影ベース法において、どのように尤度を最大化するかは任意である。ごく原始的には基底物質厚さ組み合わせごとに、網羅的に尤度を調べて一番高い尤度を示す点を選択してもよい。つまり、図４における座標空間における基底物質厚さの組み合わせのすべて（座標空間上において、とり得るすべての座標値）について尤度を調べてもよい。
ただし、組み合わせの刻み幅を細かくしていくと、座標軸数乗で計算回数が増大し現実的な手法ではない。座標軸数乗とは、例えば、座標軸が３つあれば３３という意味である。
【００２６】
そこで、一般的には、効率的にピークにたどり着くように工夫された様々な最大化（最適化）手法が用いられ、最適値が得られる。例えば、ある候補点ごとに勾配方向を求め、ラインサーチを繰り返す最急降下法や、過去のサーチベクトルと共役になる方向にラインサーチを繰り返す共役勾配法と呼ばれる最適化手法が知られている。しかし、最急降下法は評価関数の細長い谷（最大化の場合は尾根）状の箇所で繰り返し回数が非常に多くなることが知られている。また、共役勾配法では前記した評価関数の谷・尾根構造に対しても、それが二次関数であれば最小回数で収束するが、非二次関数に対しては収束が保証されないことが知られている。本実施形態では、エネルギ情報を得るＸ線ＣＴ分野において、基底物質厚さをパラメータとした実測スペクトル投影データ２３２に対する尤度の分布が尾根構造を有し、また必ずしも二次関数ではないことに着目した新しい最適化手法を提案する。
【００２７】
まず、最適化処理部２１２は、図４に示すような、それぞれの基底物質厚さと、尤度（対数尤度）を座標軸とする座標空間を設定する。
図４の例では、座標軸は横軸が第１基底物質厚さを示している。また、縦軸が第２基底物質厚さを示している。そして、紙面に垂直な方向が尤度を示している。第１基底物質は例えば水であり、第２基底物質は例えば骨である。
ここで、対数尤度は尤度の増加に対して単調増加であるため、最大尤度と最大対数尤度をとる厚さ組み合わせの位置は同一となる。すなわち、尤度の推定に対数尤度を用いても、対数尤度ではない尤度を用いても結果は同じである。本実施形態では、桁数節約等のために対数尤度を尤度の代用として使用する。
【００２８】
図４に示すように、推定される尤度等高線３００は細長い直線状の尾根状の分布を有すると推定される。なお、尤度等高線３００は未知であるが、経験上、図４に示すような形状を有するものと推測される。
また、尤度等高線３００の尾根構造は、二次関数でないことが分かっている。このため共役勾配法による最適化が適していない場合に用いられる最適化手法として、以下で説明するリッジ（尾根）サーチ法を使用する。
【００２９】
リッジサーチ法は、最初の２回のラインサーチ（以下、サーチと記載する）で尾根上の２点を見つけ出し、その２点の延長線上で第３のサーチを行う最適化手法である。
【００３０】
まず、最適化処理部２１２は実測スペクトル投影データ２３２を記憶装置２３から取得する（Ｓ２０１）。ここで、取得される実測スペクトル投影データ２３２は今回撮像された結果である。すなわち、図３のステップＳ１で撮像された結果である。
【００３１】
次に、実測スペクトル投影データ２３２を基に推定厚さ入力値（所定のパラメータ）３０１が入力される（Ｓ２０２）。
推定厚さ入力値３０１は、前記したように、画素値が一様な画像でもよいし、前記した推定厚さ初期値投影データ２４３でもよいし、なんらかの別方式による解に近い値でもよいし、逐次的に前回算出された最尤推定量を入力に戻した値でもよい。
【００３２】
最適化処理部２１２は、リッジ構造実角度３０２と同じ角度であるリッジ構造代表角度２３１を記憶装置２３から取得する（Ｓ２０３）。
リッジ構造代表角度２３１は、事前に記憶装置２３に保持されている。
ここで、リッジ構造実角度３０２は、Ｘ線エネルギが単色の場合、そのＸ線エネルギにおける基底物質のＸ線減弱係数の比で決まる値である。ここでは、リッジ構造実角度３０２を、第１基底物質厚さの座標軸と、尤度等高線３００の尾根がなす角度とする。
【００３３】
実際には、Ｘ線ＣＴでは多色Ｘ線（複数の波長を含むＸ線）が用いられている。またエネルギごとに減弱後のカウント残量が異なるため、リッジ構造実角度３０２は基底物質厚さの組み合わせが先に決まらなければ確定しない。しかし、大まかなリッジ構造実角度３０２は広い基底物質厚さ範囲で近い値を有する。すなわち、図４に示す座標空間上の、どの箇所でもリッジ構造実角度３０２は、ほぼ同じ値となる。
【００３４】
なお、リッジ構造実角度３０２は、最適化処理の最後まで未知であるが、これまでの経験からおおよその角度は知られている。つまり、リッジ構造実角度３０２（尤度等高線３００の尾根方向）は分かっているが、尤度等高線３００の尾根が座標空間上のどこに存在しているかは分からない。
次に、最適化処理部２１２は、第１基底物質厚さの座標軸となす角度がリッジ構造代表角度２３１であって（つまり、尤度等高線３００の尾根と平行）、かつ、推定厚さ入力値３０１を通るリッジ構造代表線３２１を設定する。
【００３５】
そして、最適化処理部２１２は、リッジ構造代表線３２１上に、推定厚さ入力値３０１と異なるシフト始点３０４を設定する（Ｓ２０４）。
【００３６】
さらに、最適化処理部２１２は、リッジ構造代表角度２３１を基に第１サーチ方向３０５ａ及び第２サーチ方向３０５ｂを設定する（Ｓ２０５）。図４に示すように、第１サーチ方向３０５ａ及び第２サーチ方向３０５ｂは、リッジ構造代表線３２１に対して垂直の方向である。ここで、第１サーチ方向３０５ａは後記する第１サーチが行われる方向である。同様に、第２サーチ方向３０５ｂは後記する第２サーチが行われる方向である。
【００３７】
本実施形態で用いられるリッジサーチ法におけるサーチには始点とサーチ方向が必要である。前記したように、推定厚さ入力値３０１は第１サーチの始点となる。なお、リッジ構造代表線３２１に垂直な方向は２方向存在する（図４の符号３０５ａ１，３０５ａ２の方向）。従って、最適化処理部２１２は、必要に応じて、各方向の尤度勾配を取得し、尤度が増大する方向を選択する。また、第１サーチにおける初回移動距離（サーチベクトル長さ）について一般的な最適値は存在しないため、予め定めた固定値を用いたり、系に応じてユーザが適切に選択したりするものとする。
【００３８】
次に、最適化処理部２１２は、推定厚さ入力値３０１を始点として、第１サーチ方向３０５ａの方向で最も尤度（対数尤度）の高くなる点（第１推定厚さ３０８）を探索する第１サーチ（第１のサーチ）を行う（Ｓ２１１）。ここでの「厚さ」とは基底物質厚さのことである。なお、第１サーチはシフト始点３０４の設定前に行われてもよい。
【００３９】
この結果、第１推定厚さ（第１の値）３０８が求まる。図４の例では、第１推定厚さ３０８は、第１基底物質厚さと、第２基底物質厚さの組み合わせとなる。第１サーチは、第１サーチ方向３０５ａの方向へ、定点ごとに尤度を求めてもよい。あるいは、最適化処理部２１２が、サーチのための点の幅を１，２，４・・・として点を設定し、前の点の尤度より、後の点の尤度が低い点が現れたら、前の点と、次の点との間を詳細にサーチするようにしてもよい。
ここで、第１推定厚さ３０８での尤度をＴとすると、Ｔは、以下（１）の式で示される。
【００４０】
【数１】
2019033781.tif 000003
【００４１】
ここで、ＴＭはＭ個の基底物質厚さ（Ｍ個の第２のパラメータ）である。つまり、ＴＭ＝（Ｔ１，Ｔ２，・・・，ＴＭ）である。後記する図４の例ではＴＭ＝（Ｔ１，Ｔ２）である。ここで、Ｔ１は第１基底物質厚さであり、Ｔ２は第２基底物質厚さである。また、Ｃｏｂｓ（ＥＮ）は今回撮像されたＮ個のエネルギ（Ｎ個の第１のパラメータ）を有する実測スペクトルである。つまり、ＥＮ＝（Ｅ１，Ｅ２，・・・，Ｅｉ，・・・ＥＮ）である。ＴＭはＥＮから分かるパラメータ（ＴＭはＥＮに含まれるパラメータ）である。そして、Ｃｅｓｔ（ＴＭ）は、基底物質厚さがＴＭ＝（Ｔ１，Ｔ２，・・・，ＴＭ）であるときの、推定スペクトルである。Ｃｅｓｔ（ＴＭ）は、過去のデータを基に既知である。そして、ｆ（Ｃｃｅｓｔ（ＴＭ）｜Ｃｏｂｓ（Ｅｉ））は、Ｃｏｂｓ（Ｅｉ）の下で単位時間にＴＭが生じる確率である。なお、この確率はポアソン分布を仮定している。ここで、ＴＭは変数、ＥＮは固定値（ゆえに、Ｃｏｂｓ（ＥＮ）も固定値）である。第１サーチでは、最適化処理部２１２は、推定厚さ入力値３０１を始点として、第１サーチ方向３０５ａの方向へＴＭの値を変化させながら、式（１）に示される尤度Ｌ（・）が最大となるＴＭをＴとする。
直感的には、撮像された実測スペクトル投影データ２３２が生じている条件下において、Ｍ個の基底物質厚さＴＭが生じる確率が最大となるＴＭが求められる。つまり、式（１）は、今回の撮像で得られたＥＮが得られる可能性の高い（確率が高い）ＴＭ（（Ｔ１，Ｔ２，・・・，ＴＭ）の組み合わせ）を求めている。
【００４２】
このように、リッジサーチ法におけるサーチには被検体Ｐに関わる実測データである実測スペクトル投影データ２３２が必要である。また、実測スペクトル投影データ２３２はサーチで用いる尤度計算の材料である尤度等高線の計算にも用いられる。
【００４３】
このように、最適化処理部２１２は、リッジ構造代表角度２３１が事前に準備され、そのリッジ構造代表角度２３１に垂直な方向に第１サーチを行う。このようにすることで、最初のサーチ（第１サーチ）をリッジ構造実角度３０２に対し、ほぼ垂直に行うことが可能となる。つまり、最適化処理部２１２は、最短距離で尤度等高線３００の尾根上の点である第１推定厚さ３０８に到達することができる。そして、最適化処理部２１２は第１推定厚さ３０８を出力する（Ｓ２１２）。
【００４４】
次に、最適化処理部２１２は、シフト始点３０４を始点として、第２サーチ方向３０５ｂで最も尤度の高くなる点（第２推定厚さ（第２の値）３０９）を探索する第２サーチ（第２のサーチ）を行う（Ｓ２１３）。図４の例では、第２推定厚さ３０９は、第１基底物質厚さと、第２基底物質厚さの組み合わせとなる。
この結果、第２推定厚さ３０９が出力される（Ｓ２１４）。第２サーチの手順は第１サーチと同様である。
【００４５】
このように、最適化処理部２１２は、推定厚さ入力値３０１からリッジ構造代表線３２１方向に若干シフトを与えたシフト始点３０４を設定する。シフトの量は、例えば、照射時のＸ線量から予測できる物質分解能程度の幅である。
そして、最適化処理部２１２は、シフト始点３０４を始点として、第１サーチ方向３０５ａと平行な第２サーチ方向３０５ｂの方向へ第２サーチを行う。これにより、第１サーチと同様、最適化処理部２１２は、最短距離で尤度等高線３００の尾根上の点である第２推定厚さ３０９を得ることができる。
なお、第２サーチは、第１サーチとは起点または方向が異なればよいが、図４に示すように推定厚さ入力値３０１からリッジ構造代表線３２１方向に若干シフトを与えたシフト始点３０４を起点とし、リッジ構造代表線３２１に対して垂直な方向へサーチされることが望ましい。
【００４６】
図４に示すように、第１推定厚さ３０８及び第２推定厚さ３０９の２点を結ぶ直線は、尤度等高線３００の直線状尾根構造に相当することが期待できる。最適化処理部２１２は、この方向を第３サーチ方向３１０として設定する（Ｓ２２１）。
そして、最適化処理部２１２は、第３サーチ方向３１０で最も尤度の高くなる点を探索する第３サーチ（第３のサーチ）を行う（Ｓ２２２）。
なお、第３サーチ方向３１０も２方向存在する（図４の符号３１０ａ，３１０ｂ）。しかし、第１推定厚さ３０８及び第２推定厚さ３０９での尤度の大小は分かっていることから、尤度の増大する方向を第３サーチ方向３１０として選択すればよい。また、最適化処理部２１２は、第１推定厚さ３０８と第２推定厚さ３０９のうち、尤度の大きい方（図４の場合は第２推定厚さ３０９）を第３サーチの始点として選択する。
前記したように、第３サーチによるサーチは、理想的には尤度等高線３００の尾根方向の座標空間における全空間の最尤推定量となる。最適化処理部２１２は、第３サーチ方向３１０で最も尤度の高くなる点を推定厚さ出力値（パラメータの推定値）３１２として出力する（Ｓ２２３）。図４の例では、推定厚さ出力値３１２は、第１基底物質厚さと、第２基底物質厚さの組み合わせとなる。
【００４７】
この後、投影ベース処理部２１１は出力された推定厚さ出力値３１２を基に、基底物質推定厚さ投影データ２４１（図３参照）を生成する。
【００４８】
推定厚さ入力値３０１が完全に未知の場合にはシフト始点３０４は大きくとってもよい。また、推定厚さ出力値３１２を新たな推定厚さ入力値３０１として、必要に応じてこの操作（リッジサーチ法）を繰り返してもよい（図５の破線矢印）。このように最適化処理を繰り返すことで逐次的な解の改善を図ることができる。
【００４９】
本実施形態では、尤度分布が直線構造をなす場合に対して特化したリッジサーチ法による最適化手法が用いられる。このようにすることで、高負荷な尤度計算の繰り返し回数を大きく削減することが可能となる。
また、本実施形態では、第１サーチ及び第２サーチが尤度等高線３００の平面視した尾根のラインに対して、垂直方向に行われている。このようにすることで、推定厚さ入力値３０１及びシフト始点３０４から最短距離（つまり、最も低負荷）で第１推定厚さ３０８及び第２推定厚さ３０９を探索することができる。
【００５０】
さらに、本実施形態では、最適化処理部２１２は、シフト始点３０４を、推定厚さ入力値３０１から尤度等高線３００の尾根のラインに沿って、所定距離離れた箇所に設定している。このようにすることで、第２サーチの探索距離を短くすることができ、計算負荷を軽減することができる。
【００５１】
なお、本実施形態では、Ｐｈｏｔｏｎ  Ｃｏｕｎｔｉｎｇ  ＣＴ  （ＰＣＣＴ）や、ｓｐｅｃｔｒａｌ  ＣＴを想定しているが、Ｄｕａｌ  Ｅｎｅｒｇｙ  ＣＴに対して本実施形態の手法が適用されてもよい。この場合、基底物質厚さの代わりに基底物質エネルギが用いられることとなる。
【００５２】
本発明は前記した実施形態に限定されるものではなく、様々な変形例が含まれる。例えば、前記した実施形態は本発明を分かりやすく説明するために詳細に説明したものであり、必ずしも説明したすべての構成を有するものに限定されるものではない。
【００５３】
また、前記した各構成、機能、各部２１１〜２１４、記憶装置２３等は、それらの一部またはすべてを、例えば集積回路で設計すること等によりハードウェアで実現してもよい。また、図２に示すように、前記した各構成、機能等は、ＣＰＵ２２等のプロセッサがそれぞれの機能を実現するプログラムを解釈し、実行することによりソフトウェアで実現してもよい。各機能を実現するプログラム、テーブル、ファイル等の情報は、ＨＤに格納すること以外に、メモリ２１や、ＳＳＤ（Solid State Drive）等の記録装置、または、ＩＣ（Integrated Circuit）カードや、ＳＤ（Secure Digital）カード、ＤＶＤ（Digital Versatile Disc）等の記録媒体に格納することができる。
また、各実施形態において、制御線や情報線は説明上必要と考えられるものを示しており、製品上必ずしもすべての制御線や情報線を示しているとは限らない。実際には、ほとんどすべての構成が相互に接続されていると考えてよい。
【符号の説明】
【００５４】
１      Ｘ線ＣＴ装置
２      演算装置（演算部）
３      入力装置
４      表示装置
１０    Ｘ線ＣＴシステム
２１１  投影ベース処理部
２１２  最適化処理部
２１３  画像再構成処理部
２１４  表示処理部
２３１  リッジ構造代表角度
２３２  実測スペクトル投影データ
３００  尤度等高線
３０１  推定厚さ入力値（所定のパラメータ）
３０２  リッジ構造実角度
３０４  シフト始点
３０５ａ  第１サーチ方向
３０５ｂ  第２サーチ方向
３０８  第１推定厚さ（第１の値）
３０９  第２推定厚さ（第２の値）
３１０  第３サーチ方向
３１２  推定厚さ出力値（パラメータの推定値）
３２１  リッジ構造代表線
【図１】
2019033781.tif 000004
【図２】
2019033781.tif 000005
【図３】
2019033781.tif 000006
【図４】
2019033781.tif 000007
【図５】
2019033781.tif 000008