特許 6569733 画像再構成処理方法、画像再構成処理プログラム並びにそれを搭載した断層撮影装置

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6569733

(24)【登録日】2019年8月16日

(45)【発行日】2019年9月4日

(54)【発明の名称】画像再構成処理方法、画像再構成処理プログラム並びにそれを搭載した断層撮影装置

(51)【国際特許分類】

   A61B 6/03 20060101AFI20190826BHJP

【ＦＩ】

   A61B6/03 350X

【請求項の数】7

【全頁数】15

(21)【出願番号】特願2017-535173(P2017-535173)

(86)(22)【出願日】2015年8月17日

(86)【国際出願番号】JP2015073057

(87)【国際公開番号】WO2017029702

(87)【国際公開日】20170223

【審査請求日】2018年1月30日

(73)【特許権者】

【識別番号】000001993

【氏名又は名称】株式会社島津製作所

(74)【代理人】

【識別番号】100093056

【弁理士】

【氏名又は名称】杉谷 勉

(74)【代理人】

【識別番号】100142930

【弁理士】

【氏名又は名称】戸高 弘幸

(74)【代理人】

【識別番号】100175020

【弁理士】

【氏名又は名称】杉谷 知彦

(74)【代理人】

【識別番号】100180596

【弁理士】

【氏名又は名称】栗原 要

(74)【代理人】

【識別番号】100195349

【弁理士】

【氏名又は名称】青野 信喜

(72)【発明者】

【氏名】知野見 健太

【審査官】 松岡 智也

(56)【参考文献】

【文献】 特開２０１１−２０３１６０（ＪＰ，Ａ）

【文献】 特開２０１１−１５６３０２（ＪＰ，Ａ）

【文献】 特開２０１４−２２６５６２（ＪＰ，Ａ）

【文献】 米国特許出願公開第２０１３／０１５６１６３（ＵＳ，Ａ１）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ａ６１Ｂ ６／００ − ６／１４

Ｇ０１Ｎ ２３／００ − ２３／２２７６

(57)【特許請求の範囲】

【請求項1】

再構成処理を行う画像再構成処理方法であって、
逐次近似法により画像を更新する画像更新工程と、
当該画像更新工程での画像更新毎，一定間隔毎，一定の基準を満たすタイミングあるいは任意のタイミングに再構成画像から、当該再構成画像に基づく画像の更新で生成される再構成画像のヒストグラム上の幅を持った画素値分布を急峻なピークへと近づける時に作用させる物質制約値を推定する物質制約値推定工程と
を備え、
当該物質制約値推定工程で推定された前記物質制約値を用いて前記画像更新工程での逐次近似法により画像を更新して再構成処理を行うことを特徴とする画像再構成処理方法。

【請求項2】

請求項１に記載の画像再構成処理方法において、
（ａ）既知の構成物質の数に基づいて、前記物質制約値推定工程における前記物質制約値の推定がなされることを特徴とする画像再構成処理方法。

【請求項3】

請求項１に記載の画像再構成処理方法において、
（ｂ）パラメータとして与えられた構成物質の数に基づいて、前記物質制約値推定工程における前記物質制約値の推定がなされることを特徴とする画像再構成処理方法。

【請求項4】

請求項１から請求項３のいずれかに記載の画像再構成処理方法において、
（ｃ）パラメータとして与えられた物質制約値に基づいて、前記物質制約値推定工程における前記物質制約値の推定がなされることを特徴とする画像再構成処理方法。

【請求項5】

請求項１から請求項４のいずれかに記載の画像再構成処理方法において、
再構成画像のヒストグラムあるいは再構成画像のクラスタリング結果に基づいて、前記物質制約値推定工程における前記物質制約値の推定がなされることを特徴とする画像再構成処理方法。

【請求項6】

請求項１から請求項５のいずれかに記載の画像再構成処理方法をコンピュータに実行させることを特徴とする画像再構成処理プログラム。

【請求項7】

請求項６に記載の画像再構成処理プログラムを搭載した断層撮影装置であって、
当該画像再構成処理プログラムを実行する演算手段を備えることを特徴とする断層撮影装置。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、画像再構成処理方法、画像再構成処理プログラム並びにそれを搭載した断層撮影装置における逐次近似法による再構成アーティファクトの低減手法に関する。

【背景技術】

【0002】

断層撮影装置としてＸ線ＣＴ(Computed Tomography)装置を例に採って説明する。これまで、Ｘ線ＣＴにおける標準的な画像再構成手法としてフィルター補正逆投影法(FBP: Filtered Back Projection)が用いられてきた。近年では、計算機の性能向上も相まって、逐次近似法を用いた画像再構成の研究・実用化が進んでいる。Ｘ線ＣＴでは各種要因によるアーティファクトの発生が長年の課題となっているが、逐次近似法は、アーティファクトを低減するために複雑な物理モデルや事前情報（事前知識）などを反映することができる点が特徴となっており、これまでに様々な手法が提案されている（例えば、特許文献１、２および非特許文献１参照）。

【0003】

このうち、特許文献２：特開２０１１−１５６３０２号公報および非特許文献１で提案されている手法は、ベイズ理論に基づく事後確率最大化による推定（ＭＡＰ推定）を利用するもので、撮影サンプルの構成物質に関する情報（物質情報）を事前確率として与え、それによってより良い解を得ようとする。つまり、再構成画素が事前に指定された物質の画素値（Ｘ線吸収係数を表現したもの）を持つような効力を与え、それによってアーティファクトを低減しようとするアプローチである。

【0004】

また、特許文献１：米国特許第８，１７５，３６１号公報で提案されている手法は、遂次近似法において物質情報を用いた正則化（制約）を与えており、その正則化項を画像のヒストグラムから求めている。つまり、画像のヒストグラム（特許文献１のFIG. 4を参照）に基づいて物質情報を求め、その物質情報を用いた正則化項（特許文献１のregularization term R(X)を参照）を遂次近似式（特許文献１の（３）式を参照）の右辺の第２項に組み込んで、それによってアーティファクトを低減しようとするアプローチである。

【0005】

この物質情報の効果を画像ヒストグラムの観点から説明する。なお、説明に使用する図１０〜図１２のヒストグラムの縦軸は、画素数の最大値で正規化して表現しており、ヒストグラムの横軸は右に向かうに従って画素値が高くなる。例として、互いに異なるＸ線吸収係数を持つ４種類の材質で構成される撮影サンプルを考える。材質が純物質で、雑音がないなど理想的な状況を想定すると、図１０に示すように、再構成画像のヒストグラムには４個のピークが存在することになる。

【0006】

しかし、実際には様々な要因によってアーティファクトが発生し、ヒストグラム中の各ピークは、図１１に示すように幅を持った分布となる。一方、物質情報は再構成画素が取り得る画素値集合（４つの物質制約値）として与えられ、各画素値分布の中心は１つの物質制約値と対応することになる。物質情報が働くことで、図１２に示すように分布周辺の画素値は分布中心に向かって引き寄せられる。その結果、幅を持った画素値分布は急峻なピークへと近づき、理想的な画像、つまりアーティファクトが低減された画像が得られる。

【先行技術文献】

【特許文献】

【0007】

【特許文献1】米国特許第８，１７５，３６１号公報

【特許文献2】特開２０１１−１５６３０２号公報

【非特許文献】

【0008】

【非特許文献1】C. Lemmens: Suppression of Metal Artifacts in CT Using a Reconstruction Procedure That Combines MAP and Projection Completion, IEEE Transactions on Medical Imaging, Volume:28 Issue:2(2009)

【発明の概要】

【発明が解決しようとする課題】

【0009】

しかしながら、特許文献２：特開２０１１−１５６３０２号公報および非特許文献１で提案されている従来技術では、撮影サンプルの構成物質が既知である必要があるので、構成物質が未知の撮影サンプルに適用することができないという問題点がある。また、構成物質が（例えば合金などの）純物質でない場合には、混合比率によってＸ線吸収係数は異なってくるという問題点がある。このように、構成物質が既知だとしても与えるべき物質情報を正確に知ることができないケースもあり、やはり従来技術を適用することができない。

【0010】

一方、特許文献１：米国特許第８，１７５，３６１号公報で提案されている従来技術では、画像のヒストグラムに基づいて物質情報を推定するので、構成物質が未知の撮影サンプルにも適用することができる。しかし、遂次近似法を用いた画像再構成では、繰り返し回数（反復回数）が少ない時点では画像が十分に再構成できていない（アーティファクトが多く、ぼけている）。そのため、その時点で推定された物質情報も信頼できる値ではない。

【0011】

本発明は、このような事情に鑑みてなされたものであって、撮影サンプルの構成物質が既知・未知に関わらず適用可能で、信頼度が高い物質情報を推定して、それを用いてアーティファクトを低減することができる画像再構成処理方法、画像再構成処理プログラム並びにそれを搭載した断層撮影装置を提供することを目的とする。

【課題を解決するための手段】

【0012】

本発明は、このような目的を達成するために、次のような構成をとる。
すなわち、本発明の画像再構成処理方法は、再構成処理を行う画像再構成処理方法であって、 逐次近似法により画像を更新する画像更新工程と、当該画像更新工程での画像更新毎，一定間隔毎，一定の基準を満たすタイミングあるいは任意のタイミングに再構成画像から、当該再構成画像に基づく画像の更新で生成される再構成画像のヒストグラム上の幅を持った画素値分布を急峻なピークへと近づける時に作用させる物質制約値を推定する物質制約値推定工程とを備え、当該物質制約値推定工程で推定された物質制約値を用いて画像更新工程での逐次近似法により画像を更新して再構成処理を行うことを特徴とするものである。

【0013】

本発明の画像再構成処理方法によれば、物質情報推定工程では、画像更新工程での画像更新毎，一定間隔毎，一定の基準を満たすタイミングあるいは任意のタイミングに再構成画像から物質情報を推定し、物質情報推定工程で推定された物質情報を用いて、以降の逐次近似法での逐次近似式による計算において画像を更新する。このように、再構成画像から物質情報を推定することにより、本発明は、撮影サンプルの構成物質が既知・未知に関わらず適用可能である。また、特許文献１：米国特許第８，１７５，３６１号公報では推定された物質情報が固定であるのに対して、本発明では画像更新工程での画像更新毎，一定間隔毎，一定の基準を満たすタイミングあるいは任意のタイミングに再構成画像から物質情報を推定（更新）するので、例えば逐次近似式での繰り返し回数（反復回数）が少ない時点で推定された物質情報を、繰り返し回数が多い時点でも使い続けるというような問題を回避して、信頼度が高い物質情報を推定することができる。したがって、信頼度が高い物質情報を用いてアーティファクトを低減することができる。

【0014】

上述した本発明の画像再構成処理方法の一例は、（ａ）既知の構成物質の数に基づいて、物質制約値推定工程における物質制約値を推定することである。また、画像再構成処理方法の他の一例は、（ｂ）パラメータとして与えられた構成物質の数に基づいて、物質制約値推定工程における物質制約値を推定することである。さらに、画像再構成処理方法のさらなる他の一例は、（ｃ）パラメータとして与えられた物質制約値に基づいて、物質制約値推定工程における物質制約値を推定することである。

【0015】

特に、上記（ｂ）の場合には、撮影サンプルを構成する実際の物質数に関わらず構成物質の数がパラメータとして指定された場合、このパラメータに基づき物質情報を推定する。上記（ｃ）の場合には、構成物質が既知・未知に関わらず物質制約値がパラメータとして与えられた場合、このパラメータに基づき物質情報を推定する。上記（ａ）〜上記（ｃ）によって、誤った物質情報が推定される可能性が抑えられる。上記（ａ）〜上記（ｃ）のいずれかで物質情報を推定してもよいし、上記（ａ）〜上記（ｃ）を複数組み合わせて物質情報を推定してもよい。例えば、上記（ａ）および上記（ｃ）を組み合わせて物質情報を推定してもよいし、上記（ｂ）および上記（ｃ）を組み合わせて物質情報を推定してもよい。

【0016】

上述したこれらの発明の画像再構成処理方法において、再構成画像のヒストグラムに基づいて、物質制約値推定工程における物質制約値を推定してもよいし、再構成画像のクラスタリング結果に基づいて、物質制約値推定工程における物質制約値を推定してもよい。

【0017】

また、本発明の画像再構成処理プログラムは、これらの発明の画像再構成処理方法をコンピュータに実行させることを特徴とするものである。

【0018】

本発明の画像再構成処理プログラムによれば、これらの発明の画像再構成処理方法をコンピュータに実行させることによって、撮影サンプルの構成物質が既知・未知に関わらず適用可能で、信頼度が高い物質情報を推定して、それを用いてアーティファクトを低減することができる。

【0019】

また、本発明の断層撮影装置は、本発明の画像再構成処理プログラムを搭載した断層撮影装置であって、当該画像再構成処理プログラムを実行する演算手段を備えることを特徴とするものである。

【0020】

本発明の断層撮影装置によれば、画像再構成処理プログラムを実行する演算手段を備えることによって、撮影サンプルの構成物質が既知・未知に関わらず適用可能で、信頼度が高い物質情報を推定して、それを用いてアーティファクトを低減することができる。

【発明の効果】

【0021】

本発明に係る画像再構成処理方法によれば、再構成画像から物質情報を推定することにより、本発明は、撮影サンプルの構成物質が既知・未知に関わらず適用可能である。また、画像更新工程での画像更新毎，一定間隔毎，一定の基準を満たすタイミングあるいは任意のタイミングに再構成画像から物質情報を推定（更新）するので、例えば逐次近似式での繰り返し回数（反復回数）が少ない時点で推定された物質情報を、繰り返し回数が多い時点でも使い続けるというような問題を回避して、信頼度が高い物質情報を推定することができる。したがって、信頼度が高い物質情報を用いてアーティファクトを低減することができる。
また、本発明の画像再構成処理プログラムによれば、これらの発明の画像再構成処理方法をコンピュータに実行させることによって、撮影サンプルの構成物質が既知・未知に関わらず適用可能で、信頼度が高い物質情報を推定して、それを用いてアーティファクトを低減することができる。
また、本発明の断層撮影装置によれば、画像再構成処理プログラムを実行する演算手段を備えることによって、撮影サンプルの構成物質が既知・未知に関わらず適用可能で、信頼度が高い物質情報を推定して、それを用いてアーティファクトを低減することができる。

【図面の簡単な説明】

【0022】

【図1】実施例に係るＸ線ＣＴ装置の概略図およびブロック図である。

【図2】実施例に係る画像再構成処理のフローチャートである。

【図3】区分的ガウス関数を用いた目的関数の妥当性項に関する模式図である。

【図4】ヒストグラムを用いた物質情報推定の説明に供する模式図である。

【図5】強度係数ｅ^−αｄの説明に供する模式図である。

【図6】物質情報がないときの再構成結果である。

【図7】物質情報を更新しないときの再構成結果である。

【図8】物質情報を更新したときの再構成結果である。

【図9】変形例に係る画像再構成処理のフローチャートである。

【図10】理想的なヒストグラムの模式図である。

【図11】実際のヒストグラムの模式図である。

【図12】物質情報の効果によるヒストグラムの模式図である。

【発明を実施するための形態】

【実施例】

【0023】

以下、図面を参照して本発明の実施例を説明する。図１は、実施例に係るＸ線ＣＴ装置の概略図およびブロック図である。本実施例では、断層撮影装置としてＸ線ＣＴ装置を例に採って説明する。

【0024】

図１に示すように、本実施例に係るＸ線ＣＴ装置１は、対象物Ｏを撮像する撮像部２と、対象物Ｏを載置するステージ３と、そのステージ３を駆動するステージ駆動部４と、撮像部２を駆動する撮像駆動部５と、撮像部２のＸ線管２１に管電流や管電圧を与えるために高電圧を発生する高電圧発生部６と、撮像部２のＸ線検出器２２によって得られた投影データに対して再構成処理を行う再構成処理部７とを備えている。再構成処理部７は、本発明における演算手段に相当する。

【0025】

撮像部２は、対象物ＯにＸ線を照射するＸ線管２１と、Ｘ線管２１から照射され対象物Ｏを透過したＸ線を検出するＸ線検出器２２とを備えている。Ｘ線検出器２２については、イメージインテンシファイア（Ｉ．Ｉ）やフラットパネル型Ｘ線検出器（FPD: Flat Panel Detector）などに例示されるように、特に限定されない。本実施例では、Ｘ線検出器２２としてフラットパネル型Ｘ線検出器（ＦＰＤ）を例に採って説明する。

【0026】

ＦＰＤは、画素に対応して縦横に並べられた複数の検出素子からなり、Ｘ線を検出素子が検出して、検出されたＸ線のデータ（電荷信号）をＸ線検出信号として出力する。このようにして、Ｘ線管２１からＸ線を対象物Ｏに向けて照射し、ＦＰＤからなるＸ線検出器２２がＸ線を検出してＸ線検出信号を出力する。そして、Ｘ線検出信号に基づく画素値を画素（検出素子）に対応付けて並べることで投影データを取得する。

【0027】

ステージ駆動部４は、図示を省略するモータや駆動軸などから構成され、ステージ３を図中のＺ軸心周りに水平面内で回転させる。ステージ３の水平面内での回転によって対象物ＯもＺ軸心周りに水平面内で回転し、複数の投影データを取得する。

【0028】

撮像駆動部５は、ステージ駆動部４と同様に、図示を省略するモータや駆動軸などから構成されている。Ｘ線検出器２２にＸ線管２１が対向するようにそれぞれを移動させてＸ線ＣＴ撮影を行う。また、Ｘ線管２１またはＸ線検出器２２を水平方向（図中のＸ方向）に移動させて、Ｘ線ＣＴ撮影における拡大率を変更することも可能である。また、Ｘ線管２１またはＸ線検出器２２をＸ軸に対して傾斜させて、対象物Ｏの斜め方向から撮像することも可能である。

【0029】

高電圧発生部６は、高電圧を発生させて管電流や管電圧をＸ線管２１に与えることで、Ｘ線管２１からＸ線が発生して、対象物Ｏに対してＸ線を照射する。再構成処理部７は、後述する画像再構成処理プログラム８Ａを実行することで、対象物Ｏに関する再構成画像を取得する。再構成処理部７の具体的な機能については詳しく後述する。

【0030】

その他に、Ｘ線ＣＴ装置１は、メモリ部８と入力部９と出力部１０とコントローラ１１とを備えている。

【0031】

メモリ部８は、コントローラ１１を介して、Ｘ線検出器２２で得られた投影データや再構成処理部７で得られた再構成画像などのデータを書き込んで記憶し、適宜必要に応じて読み出して、コントローラ１１を介して、投影データや再構成画像を出力部１０に送り込んで出力する。メモリ部８は、ＲＯＭ（Read-only Memory）やＲＡＭ（Random-Access Memory）やハードディスクなどに代表される記憶媒体で構成されている。

【0032】

本実施例では、投影データや更新された再構成画像（「推定画像」とも呼ばれる）のヒストグラムをメモリ部８から読み出して、コントローラ１１を介して再構成処理部７に送り込んで、逐次近似法による画像更新や物質情報推定などの画像再構成処理（図２のフローチャートを参照）を行う。また、画像再構成処理プログラム８Ａがメモリ部８に記憶されており、コントローラ１１を介して、画像再構成処理プログラム８Ａをメモリ部８から再構成処理部７に読み出して、画像再構成処理プログラム８Ａを再構成処理部７が実行することで、図２のフローチャートに示す画像再構成処理を行う。画像再構成処理プログラム８Ａは、本発明における画像再構成処理プログラムに相当する。

【0033】

入力部９は、オペレータが入力したデータや命令をコントローラ１１に送り込む。入力部９は、キーボードと、マウスやジョイスティックやトラックボールやタッチパネルなどに代表されるポインティングデバイスで構成されている。

【0034】

出力部１０は、モニタなどに代表される表示部やプリンタなどで構成されている。本実施例では、投影データや再構成画像を出力部１０のモニタに表示する。

【0035】

コントローラ１１は、Ｘ線ＣＴ装置１を構成する各部分を統括制御する。Ｘ線検出器２２で得られた投影データや再構成処理部７で得られた再構成画像などのデータを、コントローラ１１を介して、メモリ部８に書き込んで記憶、あるいは出力部１０に送り込んで出力する。出力部１０が表示部の場合には出力表示し、出力部１０がプリンタの場合には出力印刷する。

【0036】

本実施例では、再構成処理部７やコントローラ１１は、中央演算処理装置（ＣＰＵ）などで構成されている。なお、再構成処理部７については、ＧＰＵ(Graphics Processing Unit) などで構成されてもよい。

【0037】

次に、再構成処理部７（図１を参照）の具体的な機能について、図２〜図５を参照して説明する。図２は、実施例に係る画像再構成処理のフローチャートであり、図３は、区分的ガウス関数を用いた目的関数の妥当性項に関する模式図であり、図４は、ヒストグラムを用いた物質情報推定の説明に供する模式図であり、図５は、強度係数ｅ^−αｄの説明に供する模式図である。

【0038】

（ステップＳ１）逐次近似法による画像更新
各種の逐次近似法により画像を更新する。なお、Ｘ線管２１（図１を参照）やＸ線検出器２２（図１を参照）の物理的特性（例えばビームハードニング・散乱等）を補正する処理を含んでいるのが好ましい。本実施例では当該処理を行っているが、これらの特性が無視できるような場合には、この処理を行わなくてもよい。その他、物理的特性の補正の有無や順序を適宜変更する場合も、本発明に含まれる。

【0039】

一般的に、目的関数最大化に基づく逐次近似法では、下記（１）式で表される目的関数Ｆを最大化する。なお、実際の計算では、関数の傾き（一階微分）のみから、関数の最小値を探索する勾配法のアルゴリズム（「最急降下法」とも呼ばれる）や、Newton法などの最適化アルゴリズムを使用する。また、局所解に陥ることを回避するために、遺伝的アルゴリズムや焼きなまし法などの、組み合わせ最適化法が組み込まれていてもよい。

【0040】

Ｆ（μ，ｙ）＝Ｄ（μ，ｙ）＋βＲ（μ） …（１）
ただし、上記（１）式中のμは再構成画像ベクトルであり、ｙは投影データである。Ｄは「データ項」と呼ばれ、実測データとの適合度を表し、実測投影（Ｘ線検出器２２で得られた実測の投影データ）および推定パラメータ（上記（１）式で推定された推定画像）から算出される尤度などで定義される。なお、μ，ｙはベクトルであるので、実際には太字で表記されることに留意されたい。

【0041】

Ｒは一般的に「ペナルティ項」などと呼ばれ、推定パラメータ（推定画像）の妥当性を反映する。本明細書では、以降Ｒを便宜上「妥当性項」と呼ぶ。本発明で利用している、ステップＳ２以降で詳しく述べる物質情報（物質制約値ｍ［１］，…，ｍ［Ｎ］）は、この妥当性項に反映される。なお、βは妥当性項Ｒの強さをコントロールする係数である。

【0042】

妥当性項の具体例としては、図３のような区分的ガウス関数などが挙げられる。この場合、推定された複数の物質制約値ｍ［１］，…，ｍ［Ｎ］を中心としたガウス分布を連結して妥当性項とする。推定画素値がガウス分布の中心に近いほど目的関数は大きくなり、推定画素値は物質制約値に近い値を取る作用が働くことになる。

【0043】

なお、各分布の平均および分散がパラメータとなるので、与えられるパラメータの数は（物質の数）×２となるが、繰り返し毎に更新するのは平均値のみである。つまり、後述するステップＳ６の「物質制約値の更新」とは、繰り返し毎に再構成画像のヒストグラムを参照して、各ガウス分布の中心位置を適切な位置に補正する（ずらす）ことを意味する。

【0044】

各ガウス分布の高さや幅（分散）については経験則に基づいて設定する。また、各ガウス分布の高さや幅（分散）を個々に設定してもよいし、互いに隣接するガウス分布の切り換え位置を適宜に変更して設定してもよい。このステップＳ１は、本発明における画像更新工程に相当する。

【0045】

（ステップＳ２）ヒストグラム生成
ステップＳ１で更新された再構成画像（推定画像）からヒストグラムを生成する。具体的には、図４に示すように、縦軸を画素数の最大値で正規化して、横軸を画素値相当のビン幅ｗとしたヒストグラムを生成する。上述したようにヒストグラムの横軸は右に向かうに従って画素値が高くなる。

【0046】

（ステップＳ３）ピーク検出
ステップＳ２で生成されたヒストグラムから極値をピークとして検出する。具体的には、図４に示すように、ヒストグラムのｋ番目のビンの高さをｈ［ｋ］とした時、ｈ［ｋ−１］＜ｈ［ｋ］＞ｈ［ｋ＋１］を満たす全てのピークを極値とみなす（図４では「●」で図示）。

【0047】

（ステップＳ４）ピーク評価
次に、ステップＳ３で検出された全てのピークについて、構成物質らしさを表す評価値を算出する。具体的には、図４に示すように、ピークとして検出されたｋ番目のビンの高さｈ［ｋ］（「ヒストグラム値」または「ピーク」とも呼ばれる）に対して、評価値ｅ［ｋ］を与える。この評価値ｅ［ｋ］が大きいほど、そのピークが撮影サンプルの構成物質を表していると考える。例えば、ピークとして検出されたヒストグラム値ｈ［ｋ］そのものを評価値ｅ［ｋ］とする方法が考えられる。

【0048】

その他に、例えば、撮影サンプルの構成物質の一部が既知の場合、その物質制約値を基準物質制約値とする。あるいは、（ｃ）のように物質制約値がパラメータとして与えられている場合、その物質制約値を基準物質制約値とする。図５に示すように、評価値ｅ［ｋ］の対象となる画素値と、基準物質制約値との間の距離（図５では「基準物質制約値からの距離」で表記）をｄ（ただしｄは非負）とすると、基準物質制約値からの距離ｄおよびヒストグラム値ｈ［ｋ］に基づいて評価値ｅ［ｋ］を定義する方法も考えられる。

【0049】

これは、推定される物質制約値は、既知またはパラメータとして与えられた物質制約値に近い値を取る、との考えに基づいている。つまり、基準物質制約値（図５を参照）から遠ざかれば評価値ｅ［ｋ］が小さくなると考えられる。そこで、自然対数の底（ネイピア数）ｅおよび基準物質制約値からの距離ｄを用いて、評価値ｅに関する強度係数をｅ^−αｄ（αは定数）とする。これにより基準物質制約値での画素値から遠ざかれば基準物質制約値からの距離ｄが大きくなる（長くなる）ので、強度係数、さらには評価値ｅ［ｋ］が小さく定義される。

【0050】

なお、基準物質制約値は不変である必要はない。このパラメータとして与えられた物質制約値や既知である物質制約値も更新対象としてもよい。

【0051】

以上により、評価値ｅ［ｋ］は下記（２）式で表される。
ｅ［ｋ］＝ｈ［ｋ］×ｅ^−αｄ …（２）
以上のように、ピークとして検出されたヒストグラム値ｈ［ｋ］そのものを評価値ｅ［ｋ］としてもよいし、上記（２）式のように基準物質制約値からの距離ｄおよびヒストグラム値ｈ［ｋ］に基づいて評価値ｅ［ｋ］を定義してもよい。

【0052】

（ステップＳ５）ピーク選定
次に、ステップＳ４で与えられた評価値ｅ［ｋ］、および必要であれば与えられたパラメータに基づいて、構成物質と思われるピークを選定する。例えば、撮影サンプルがＮ種類の物質で構成されているとすると、評価値が大きい順にＮ個のピークを抽出し、これらを選定されたピーク（選定ピーク）とする。構成物質の数が未知の場合（例えば御影石や自然鉱物（天然鉱物）の場合や、合金が完全に混じり合わない場合）には、一定の評価値以上のピークを選定ピークとする。

【0053】

その他に、撮影サンプルを構成する実際の物質数（例えばＮ個）に関わらず構成物質の数がパラメータ（Ｎ個と異なる数）として指定された場合、このパラメータに基づいて評価値が大きい順にピークを抽出し、これらを選定ピークとしてもよい。なお、選定されたピークがどの物質に対応するのかが判っている必要はない。また、選定されなかったピークは、この時点で破棄される。

【0054】

これにより、Ｎ個の物質情報を示すピーク（ピーク番号：ｋ_１，…，ｋ_Ｎ）が、図４に示すように選定されたことになる。選定されたピーク番号ｋ_１，…，ｋ_Ｎは連続番号となっていないので、選定されたピーク番号が、１，…，Ｎと連続番号になるように、ｋ_１＝１，…，ｋ_Ｎ＝Ｎと並び替える。そして、ｎ番目のヒストグラムのビンと対応する画素値をｖ［ｋ_ｎ］＝ｗ×（ｋ_ｎ−０．５）とすると、物質情報の具体的な値である物質制約値は、ｖ［ｋ_１］＝ｍ［１］，…，ｖ［ｋ_Ｎ］＝ｍ［Ｎ］として求まる。

【0055】

上記式（ｖ［ｋ_ｎ］＝ｗ×（ｋ_ｎ−０．５））では、ピーク番号ｋ_ｎからビン幅ｗの中央にずらすためにピーク番号ｋ_ｎから０．５を減算したが、物質制約値を求める式については上記式に限定されない。ピーク番号ｋ_ｎを変数とした線形関数により物質制約値を求めてもよい。

【0056】

（ステップＳ６）物質制約値の更新
以上の処理により、Ｎ個の物質情報（物質制約値＝各物質のＸ線吸収係数を表す画素値）であるｍ［１］，…，ｍ［Ｎ］が推定されて更新されたことになる。つまり、ステップＳ５で選定されたピーク（選定ピーク）に対応する画素値を、次のステップＳ１に戻る際のステップＳ１の逐次近似計算で使用する物質制約値として設定する。以上のように、ステップＳ２〜Ｓ６は、本発明における物質情報推定工程に相当する。

【0057】

（ステップＳ７）反復回数のカウンタ変数のインクリメント
逐次近似式での繰り返し回数（反復回数）のカウンタ変数をインクリメントする。

【0058】

（ステップＳ８）画像更新の終了？
逐次近似法による画像更新を終了する反復回数をＮ_ｉｔｅｒとし、カウンタ変数が反復回数Ｎ_ｉｔｅｒに達したか？否かを判断する。なお、反復回数Ｎ_ｉｔｅｒについてはオペレータが予め設定すればよい。カウンタ変数がＮ_ｉｔｅｒ以下の場合にはステップＳ１〜Ｓ６を継続するために、ステップＳ１に戻る。カウンタ変数がＮ_ｉｔｅｒを超えた場合には一連の計算を終了する。

【0059】

このようにして得られた推定画像を再構成画像として取得する。また、反復回数Ｎ_ｉｔｅｒを設定せずに、更新の度に得られた推定画像をオペレータが観察して、観察結果に基づいて一連の計算をオペレータが中断して、そのときに得られた推定画像を再構成画像として取得してもよい。または、何らかの収束評価値（目的関数の値など）が判定基準値を上回った、あるいは下回ったかによって判定してもよい。

【0060】

本実施例に係る画像再構成処理方法によれば、物質情報推定工程（図２ではステップＳ２〜Ｓ６）では、画像更新工程（図２ではステップＳ１）での画像更新毎に再構成画像から物質情報（本実施例では物質制約値ｍ［１］，…，ｍ［Ｎ］）を推定し、物質情報推定工程（ステップＳ２〜Ｓ６）で推定された物質情報（物質制約値ｍ［１］，…，ｍ［Ｎ］）を用いて、以降の逐次近似法での逐次近似式による計算において画像を更新する。このように、再構成画像から物質情報（物質制約値ｍ［１］，…，ｍ［Ｎ］）を推定することにより、本発明は、撮影サンプルの構成物質が既知・未知に関わらず適用可能である。

【0061】

また、「課題を解決するための手段」の欄でも述べたように、特許文献１：米国特許第８，１７５，３６１号公報では推定された物質情報が固定であるのに対して、本実施例では画像更新工程（ステップＳ１）での画像更新毎に再構成画像から物質情報（物質制約値ｍ［１］，…，ｍ［Ｎ］）を推定（更新）するので、例えば逐次近似式での繰り返し回数（反復回数）が少ない時点で推定された物質情報を、繰り返し回数が多い時点でも使い続けるというような問題を回避して、信頼度が高い物質情報（物質制約値ｍ［１］，…，ｍ［Ｎ］）を推定することができる。したがって、信頼度が高い物質情報（物質制約値ｍ［１］，…，ｍ［Ｎ］）を用いてアーティファクトを低減することができる。

【0062】

また、（ａ）既知の構成物質の数に基づいて、推定すべき物質情報を推定してもよいし、（ｂ）パラメータとして与えられた構成物質の数に基づいて、推定すべき物質情報を推定してもよいし、（ｃ）パラメータとして与えられた物質制約値に基づいて、推定すべき物質情報を推定してもよい。上記（ａ）の場合には、既知の構成物質の数（例えばＮ個）に基づいて、例えば、ステップＳ５でも述べたように評価値が大きい順にＮ個のピークを抽出し、これらを選定ピークとして、物質情報（物質制約値ｍ［１］，…，ｍ［Ｎ］）を推定する。

【0063】

また、上記（ｂ）の場合には、撮影サンプルを構成する実際の物質数（例えばＮ個）に関わらず構成物質の数（Ｎ個と異なる数）がパラメータとして指定された場合、このパラメータに基づき物質情報を推定する。例えば、ステップＳ５でも述べたように実際の物質数であるＮ個と異なる数がパラメータとして指定された場合、このパラメータに基づいて評価値が大きい順にピークを抽出し、これらを選定ピークとして、物質制約値を推定する。

【0064】

また、上記（ｃ）の場合には、構成物質が既知・未知に関わらず物質制約値がパラメータとして与えられた場合、このパラメータに基づき物質情報を推定する。例えば、ステップＳ４でも述べたように物質制約値がパラメータとして与えられた場合、その物質制約値を基準物質制約値とする。そして、基準物質制約値に基づいて評価値ｅ［ｋ］を求めて、ステップＳ５でも述べたように評価値ｅ［ｋ］に基づいてピークを選定して、物質制約値を推定する。

【0065】

上記（ａ）〜上記（ｃ）によって、誤った物質情報が推定される可能性が抑えられる。「課題を解決するための手段」の欄でも述べたように、上記（ａ）〜上記（ｃ）のいずれかで物質情報を推定してもよいし、上記（ａ）〜上記（ｃ）を複数組み合わせて物質情報を推定してもよい。例えば、上記（ａ）および上記（ｃ）を組み合わせて物質情報を推定してもよいし、上記（ｂ）および上記（ｃ）を組み合わせて物質情報を推定してもよい。

【0066】

本実施例では、図２のフローチャートに示すように、ステップＳ２の再構成画像のヒストグラムに基づいて、推定すべき物質情報（物質制約値ｍ［１］，…，ｍ［Ｎ］）を推定している。

【0067】

本実施例に係る画像再構成処理プログラム８Ａ（図１を参照）によれば、本実施例に係る画像再構成処理方法（図２のフローチャートを参照）をコンピュータ（本実施例では図１に示す再構成処理部７を構成するＣＰＵまたはＧＰＵ）に実行させることによって、撮影サンプルの構成物質が既知・未知に関わらず適用可能で、信頼度が高い物質情報（物質制約値ｍ［１］，…，ｍ［Ｎ］）を推定して、それを用いてアーティファクトを低減することができる。

【0068】

本実施例に係る断層撮影装置（本実施例ではＸ線ＣＴ装置）によれば、画像再構成処理プログラム８Ａを実行する演算手段（本実施例では図１に示す再構成処理部７を構成するＣＰＵまたはＧＰＵ）を備えることによって、撮影サンプルの構成物質が既知・未知に関わらず適用可能で、信頼度が高い物質情報（物質制約値ｍ［１］，…，ｍ［Ｎ］）を推定して、それを用いてアーティファクトを低減することができる。

【0069】

［再構成結果］
再構成結果について、図６〜図８を参照して説明する。図６は、物質情報がないときの再構成結果であり、図７は、物質情報を更新しないときの再構成結果であり、図８は、物質情報を更新したときの再構成結果である。図６や図７では、斜め方向のスジ状アーティファクト（ストリークアーティファクト）が観察される。これらに対して本発明のように物質情報を更新した場合には、図８ではアーティファクトがなくなり画質が改善したのが確認された。

【0070】

本発明は、上記実施形態に限られることはなく、下記のように変形実施することができる。

【0071】

（１）上述した実施例では、断層撮影装置としてＸ線ＣＴ装置を例に採って説明したが、逐次近似法による再構成処理を行う断層撮影装置であれば、特に限定されない。磁気共鳴診断装置(MRI: Magnetic Resonance Imaging) や光ＣＴ装置やＸ線以外の放射線（α線、β線、γ線など）断層撮影装置などに適用してもよい。

【0072】

（２）上述した実施例では、図１に示すような工業用や産業用の検査装置に適用したが、被検体を人体または小動物とした医用装置に適用してもよい。

【0073】

（３）単一波長Ｘ線（単色Ｘ線）や複数の波長からなるＸ線（多色Ｘ線）などに例示されるように、適用するＸ線の種類については特に限定されない。

【0074】

（４）上述した実施例では、図１に示す撮影態様であったが、トモシンセシス(tomosynthesis) などに例示されるように、断層撮影に関する撮影態様については特に限定されない。

【0075】

（５）上述した実施例では、画像更新工程（図２ではステップＳ１）での画像更新毎に再構成画像から物質情報（実施例では物質制約値ｍ［１］，…，ｍ［Ｎ］）を推定したが、一定間隔毎，一定の基準を満たすタイミングあるいは任意のタイミングに再構成画像から物質情報（物質制約値ｍ［１］，…，ｍ［Ｎ］）を推定してもよい。

【0076】

（６）上述した実施例でのヒストグラムについては、再構成画像の全画素から生成する場合に限らず、任意の画素集合、あるいはサンプリング周波数を下げる変換処理であるダウンサンプリング処理などの何らかの画像処理を適用した画像から生成することも可能である。

【0077】

（７）上述した実施例でのピーク検出（図２のフローチャートのステップＳ３を参照）については、極値周辺の面積を考慮して検出することも可能である。

【0078】

（８）上述した実施例では、図２のフローチャートに示すように、ステップＳ２の再構成画像のヒストグラムに基づいて、推定すべき物質情報（実施例では物質制約値ｍ［１］，…，ｍ［Ｎ］）を推定したが、図９のフローチャートに示すように、再構成画像のクラスタリング結果に基づいて、推定すべき物質情報（物質制約値ｍ［１］，…，ｍ［Ｎ］）を推定してもよい。具体的には、クラスタの平均を用い、与えられたクラスタ数ｋ個に分類するk-means法（「k平均法」とも呼ばれる）などに例示されるクラスタリング処理を再構成画素に対して適用し（ステップＴ２）、ステップＴ３で算出されたクラスタ重心を選定し（ステップＴ４）、選定されたクラスタ重心の位置の画素値を物質制約値とする（ステップＴ５）。図９のフローチャートを実行することで、実施例に係る図２のフローチャートと同等の効果が得られる。ここでのｋは、実施例でのヒストグラムのビン番号を表すｋと意味が相違することに留意されたい。

【産業上の利用可能性】

【0079】

以上のように、本発明は、Ｘ線ＣＴ装置（例えばトモシンセシス装置）やＭＲＩ装置や光ＣＴ装置などのような工業用や産業用の検査装置や医用装置などに適している。

【符号の説明】

【0080】

１ … Ｘ線ＣＴ装置
７ … 再構成処理部
８Ａ … 画像再構成処理プログラム
ｍ［１］，…，ｍ［Ｎ］ … 物質制約値
ｈ［ｋ］ … ヒストグラムのｋ番目のビンの高さ（ヒストグラム値）（ピーク）
ｅ［ｋ］ … 評価値

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】