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

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

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

(11)【特許番号】6558491

(24)【登録日】2019年7月26日

(45)【発行日】2019年8月14日

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

(51)【国際特許分類】

   G01N 23/046 20180101AFI20190805BHJP

   A61B 6/03 20060101ALI20190805BHJP

   A61B 5/055 20060101ALI20190805BHJP

【ＦＩ】

   G01N23/046

   A61B6/03 360J

   A61B6/03 350X

   A61B5/055 376

【請求項の数】8

【全頁数】15

(21)【出願番号】特願2018-503977(P2018-503977)

(86)(22)【出願日】2016年3月11日

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

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

(87)【国際公開日】20170914

【審査請求日】2018年8月16日

(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】

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

【請求項8】

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

【発明の詳細な説明】

【技術分野】

【0001】

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

【背景技術】

【0002】

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

【0003】

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

【0004】

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

【0005】

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

【0006】

上記手法は、目的関数最大化に基づく逐次近似法として考えることができる。この手法は、下記（１）式で表される目的関数Ｆを最大化することで、再構成画像を求める。

【0007】

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

【0008】

Ｒは一般的に「ペナルティ項」などと呼ばれ、推定パラメータ（推定画像）の妥当性を反映する。本明細書では、以降Ｒを便宜上「妥当性項」と呼ぶ。上述した物質情報は、この妥当性項に反映されており、例えば図１１のような区分的ガウス関数が利用される。なお、βは妥当性項Ｒの強さをコントロールする係数であり、経験的に決められることが多い。

【0009】

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

【0010】

μ^ｎ＋１＝μ^ｎ＋α×∇Ｆ（μ，ｙ）
＝μ^ｎ＋α×∇Ｄ（μ，ｙ）＋α×β×∇Ｒ（μ） …（２）
ただし、上記（２）式中の∇は勾配であり、推定パラメータ（再構成画像）に関する偏微分である。

【0011】

ｊ番目の画素における更新式は、下記（３）式で表される。

【0012】

【数1】

【先行技術文献】

【特許文献】

【0013】

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

【特許文献2】米国特許第８，９５８，６６０号公報

【非特許文献】

【0014】

【非特許文献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)

【発明の概要】

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

【0015】

しかしながら、離散信号を処理する以上、上記手法では高画質な再構成画像が得られないという問題点がある。具体的には、離散信号を処理する以上、例えば図１２に示すように、物体境界に位置するなど、複数の物質が混合している画素（以下、この画素を「混合画素」と呼ぶ）が存在することになる。混合画素は、各物質の画素値が一定比率で足しあわされた画素値（中間画素値）を持つ。

【0016】

物質情報はアーティファクトが発生している非混合画素に対しては有効に作用するが、このような中間画素値を持つ混合画素に対しては、特定の物質の画素値を持つような作用が働いてしまうので、不適切な効果が得られる。具体的には、図１３（ａ）は物質情報なしのときであり、図１３（ｂ）は物質情報ありのときに、微細な構造が消失したり（図１３（ｂ）の物体中央部を参照）、本来滑らかであるべき物質情報が不自然にがたついたりする（図１３（ｂ）の物体輪郭部を参照）。

【0017】

従来技術では、非混合画素と混合画素とを区別せずに、物質情報による制約をどの画素に対しても一定の強度で与えるので、この現象を避けることが困難であった。このことを上記（１）式において考えてみると、混合画素の存在を想定していないという意味で、目的関数に適切な事前知識が反映されていないとみることができる。

【0018】

本発明は、このような事情に鑑みてなされたものであって、高画質な再構成画像を得ることができる画像再構成処理方法、画像再構成処理プログラム並びにそれを搭載した断層撮影装置を提供することを目的とする。

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

【0019】

本発明は、このような目的を達成するために、次のような構成をとる。
すなわち、本発明の画像再構成処理方法は、再構成処理を行う画像再構成処理方法であって、逐次近似法により画像を更新する画像更新工程と、当該画像更新工程での画像の更新によって得られた再構成画像から、事前知識に対する重み係数マップを生成して、当該重み係数マップに基づいて各画素に対する事前知識の重み係数を調整することで前記重み係数マップを更新する重み係数マップ更新工程とを備え、当該重み係数マップ更新工程で更新された前記重み係数マップを事前知識に適用して前記画像更新工程での逐次近似法により画像を更新して再構成処理を行うことを特徴とするものである。

【0020】

本発明の画像再構成処理方法によれば、画像の更新によって得られた（推定中の）再構成画像から、事前知識に対する重み係数マップを生成して、当該重み係数マップに基づいて各画素に対する事前知識の重み係数を調整することで、高画質な再構成画像が得られないという問題点を解決することができる。つまり、各画素に対する事前知識の重み係数を調整することで、再構成画像の画素に対して過度な制約を与えることを回避することができ、高画質な再構成画像を得ることができる。

【0021】

なお、特許文献２：米国特許第８，９５８，６６０号公報で提案されている手法は、ボクセル依存型スケーリング・ファクタを算出するための係数マップを目的関数の勾配に適用する（特許文献２のclaim 1, 4を参照）。つまり、特許文献２では目的関数全体から算出される更新量に対して重み係数マップを生成する点で本発明と相違する。

【0022】

また、特許文献２では目的関数全体に対して重み係数を調整しようとすると、データ項にまで重み係数が掛けられ、データ項の更新量までが抑制されてしまい、逐次近似法における再構成の処理スピードが遅くなる。その結果、逐次近似法での繰り返し回数（反復回数）を増やさなければならない。これに対して、本発明の場合には重み係数マップを事前知識のみに適用するので、逐次近似法における再構成の処理スピードが速くなり、逐次近似法での繰り返し回数（反復回数）を減らすことができる。

【0023】

上述した本発明の画像再構成処理方法において、（ａ）互いに異なる複数の事前知識に対して、同一の重み係数マップを適用してもよいし、（ｂ）互いに異なる複数の事前知識に対して、それぞれ異なる重み係数マップを適用してもよい。このように、（ａ）、（ｂ）によって、複数の事前知識に本発明を適用することが可能となる。もちろん、１種類の事前知識に本発明を適用してもよい。

【0024】

重み係数マップの一例は、複数の物質が混合した画素である混合画素の情報を反映したマップである。混合画素の情報を用いて当該情報を反映したマップを重み係数マップとすることで、混合画素に対する事前知識（物質情報）の効果が弱められる。その結果、物質情報によるアーティスト低減効果を維持しながらも、微細構造が保持された再構成画像（図４を参照）や、滑らかな境界を持つ再構成画像（図５および図６のプロファイルを参照）が得られる。

【0025】

具体的には、再構成画像のエッジ情報を用いて重み係数マップを生成する。例えば、エッジ強度が大きいほど、重み係数マップの値が小さくなるように設定する。もちろん、それ以外にも再構成画像を表示（モニタリング）して、表示された再構成画像の（境界と認識される）任意の場所（画素）をオペレータ（ユーザ）が指定することで、任意の場所に手動で入力された情報を用いて重み係数マップを生成してもよい。なお、微細なエッジを連結するなどの目的のために、エッジ情報としてエッジ画像に平滑化フィルターなどを適用してもよい。

【0026】

なお、境界を滑らかにすることを目的として再構成処理後に再構成画像そのものに平滑化を行うことも考えられる。重み係数マップを伴わない物質情報制約付きの逐次近似法では、微細構造が保持されずに再構成画像が得られてしまうので、かかる微細構造が保持されずに得られた再構成画像に対して平滑化を行ったとしても、微細構造を復元することができない。

【0027】

上述した本発明の画像再構成処理方法において、重み係数マップを更新するタイミングは、画像更新毎，一定間隔毎，一定の基準を満たすタイミングあるいは任意のタイミングである。

【0028】

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

【0029】

本発明の画像再構成処理プログラムによれば、これらの発明の画像再構成処理方法をコンピュータに実行させることによって、再構成画像の画素に対して過度な制約を与えることを回避することができ、高画質な再構成画像を得ることができる。

【0030】

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

【0031】

本発明の断層撮影装置によれば、画像再構成処理プログラムを実行する演算手段を備えることによって、再構成画像の画素に対して過度な制約を与えることを回避することができ、高画質な再構成画像を得ることができる。

【発明の効果】

【0032】

本発明に係る画像再構成処理方法によれば、画像の更新によって得られた再構成画像から、事前知識に対する重み係数マップを生成して、当該重み係数マップに基づいて各画素に対する事前知識の重み係数を調整することで、高画質な再構成画像が得られないという問題点を解決することができる。つまり、各画素に対する事前知識の重み係数を調整することで、再構成画像の画素に対して過度な制約を与えることを回避することができ、高画質な再構成画像を得ることができる。
また、本発明の画像再構成処理プログラムによれば、これらの発明の画像再構成処理方法をコンピュータに実行させることによって、再構成画像の画素に対して過度な制約を与えることを回避することができ、高画質な再構成画像を得ることができる。
また、本発明の断層撮影装置によれば、画像再構成処理プログラムを実行する演算手段を備えることによって、再構成画像の画素に対して過度な制約を与えることを回避することができ、高画質な再構成画像を得ることができる。

【図面の簡単な説明】

【0033】

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

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

【図3】各ステップにおける生成データであり、（ａ）はステップＳ２で得られた再構成画像μ、（ｂ）はステップＳ３で得られたエッジ画像Ｅ、（ｃ）はステップＳ４で得られた重み係数マップＷである。

【図4】微細構造保持に関する重み係数マップの適用結果であり、（ａ）および（ｂ）は物質情報による制約がない逐次近似法による再構成画像、（ｃ）および（ｄ）は重み係数マップなしの物質情報による制約がある逐次近似法による再構成画像、（ｅ）および（ｆ）は重み係数マップありの物質情報による制約がある逐次近似法による再構成画像である。

【図5】ガタツキ低減に関する重み係数マップの適用結果であり、（ａ）は重み係数マップなしの再構成画像、（ｂ）は重み係数マップありの再構成画像である。

【図6】物質境界のプロファイルである。

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

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

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

【図10】物質情報の効果における適用事例であり、（ａ）は物質情報による制約がない逐次近似法による再構成画像、（ｂ）は物質情報による制約がある逐次近似法による再構成画像である。

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

【図12】混合画素の模式図である。

【図13】混合画素に適用した場合の物質情報の副次的効果であり、（ａ）は物質情報なしのときの模式図、（ｂ）は物質情報ありのときの模式図である。

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

【実施例】

【0034】

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

【0035】

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

【0036】

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

【0037】

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

【0038】

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

【0039】

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

【0040】

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

【0041】

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

【0042】

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

【0043】

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

【0044】

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

【0045】

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

【0046】

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

【0047】

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

【0048】

次に、再構成処理部７（図１を参照）の具体的な機能について、図２〜図６を参照して説明する。図２は、実施例に係る画像再構成処理のフローチャートであり、図３は、各ステップにおける生成データであり、図３（ａ）は、ステップＳ２で得られた再構成画像μであり、図３（ｂ）は、ステップＳ３で得られたエッジ画像Ｅであり、図３（ｃ）は、ステップＳ４で得られた重み係数マップＷであり、図４は、微細構造保持に関する重み係数マップの適用結果であり、図４（ａ）および図４（ｂ）は物質情報による制約がない逐次近似法による再構成画像であり、図４（ｃ）および図４（ｄ）は重み係数マップなしの物質情報による制約がある逐次近似法による再構成画像であり、図４（ｅ）および図４（ｆ）は重み係数マップありの物質情報による制約がある逐次近似法による再構成画像であり、図５は、ガタツキ低減に関する重み係数マップの適用結果であり、図５（ａ）は重み係数マップなしの再構成画像であり、図５（ｂ）は重み係数マップありの再構成画像であり、図６は、物質境界のプロファイルである。

【0049】

（ステップＳ１）初期重み係数設定
事前知識の重み係数マップをＷとする。事前知識の重み係数マップＷの各要素を１．０で初期化する。なお、ｊ番目の画素に対応するマップ係数をＷ_ｊで表す。また、何らかの初期マップを用いて重み係数マップＷを初期化してもよい。ここで、事前知識として、例えば上述した物質情報や複数の画素値における滑らかさを示す平坦情報などがある。

【0050】

（ステップＳ２）逐次近似法による画像更新
再構成画像をμとする。各種の逐次近似法により再構成画像μを更新する（図３（ａ）を参照）。最適化アルゴリズムとして上述した最急降下法を使用した場合、ｎ回目の繰り返しにおける再構成画像μ中の各画素μ_ｊの更新式は、下記（４）式で表される。

【0051】

【数2】

【0052】

上記（４）式は上記（３）式のβにＷ_ｊを掛けた形となっており、これにより画素位置に応じて事前知識の重みに異なる値を持たせることが可能となる。逆に、Ｗ_ｊを１とすれば上記（３）式と等価となる。上記（３）式は上記（１）式を最大化（最適化）する更新式であるのに対し、上記（４）式は上記（１）式から直接導出していないので、厳密には上記（１）式を最大化（最適化）する更新式にはなっていない。つまり、この係数βの置き換えは、間接的に目的関数の事前知識を修正していることになる。なお、上記では最急降下法を例に説明したが、最適化アルゴリズムとしては上述したNewton法などを用いてもよく、上述した遺伝的アルゴリズムや焼きなまし法などの組み合わせ最適化法が適用されていてもよい。

【0053】

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

【0054】

（ステップＳ３）再構成画像からエッジ抽出
エッジ画像をＥとする。ステップＳ２での処理の画像更新によって得られた再構成画像μからエッジを抽出し、エッジ画像Ｅを生成する（図３（ｂ）を参照）。なお、ｊ番目の画素におけるエッジ強度をＥ_ｊで表す。エッジ抽出手法としては、SobelフィルターやLaplacianフィルター等を利用する。これらのフィルタリング処理によって、再構成画像μの画素値に依存したエッジ強度Ｅ_ｊ（具体的には一次微分値・二次微分値）が算出される。なお、微細なエッジを連結するなどの目的のために、エッジ画像Ｅに平滑化フィルターなどを適用してもよい。

【0055】

なお、「課題を解決するための手段」の欄でも述べたように、境界を滑らかにすることを目的として再構成処理後に再構成画像そのものに平滑化を行うことも考えられる。重み係数マップを伴わない物質情報制約付きの逐次近似法では、微細構造が保持されずに再構成画像が得られてしまうので、かかる微細構造が保持されずに得られた再構成画像に対して平滑化を行ったとしても、微細構造を復元することができない。

【0056】

（ステップＳ４）重み係数マップ更新
ステップＳ３での処理で抽出されたエッジ画像Ｅを用いて、重み係数マップＷを更新する（図３（ｃ）を参照）。エッジ強度Ｅ_ｊ（ｊ：画素番号）が大きいほど、重み係数マップの値Ｗ_ｊが小さくなるように、以下のような計算式を用いる。
計算式の例：Ｗ_ｊ＝ｅｘｐ（−γ×Ｅ_ｊ） （γは定数）
この処理によって再構成画像μのエッジ部分に対応するマップ係数Ｗ_ｊに適当な値が設定されることになる。以上のように、このステップＳ４は、本発明における重み係数マップ更新工程に相当する。次に、逐次近似式での繰り返し回数（反復回数）のカウンタ変数ｎをインクリメントする。

【0057】

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

【0058】

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

【0059】

以上のステップＳ１〜Ｓ５の操作により、混合画素に対する事前知識（物質情報）の効果が弱められる。その結果、物質情報によるアーティスト低減効果を維持しながらも、微細構造が保持された再構成画像（図４を参照）や、滑らかな境界を持つ再構成画像（図５および図６のプロファイルを参照）が得られる。

【0060】

図４（ａ），図４（ｃ）および図４（ｅ）は、同一の画像からそれぞれ得られた逐次近似法による各々の再構成画像であり、図４（ｂ），図４（ｄ）および図４（ｆ）は、同一の画像からそれぞれ得られた逐次近似法による各々の再構成画像である。物質情報による制約がない場合には、図４（ａ）では左斜め方向のアーティファクトが生じていることが確認され、重み係数マップなしの物質情報による制約がある場合には、図４（ｄ）では微細構造が保持されていないことが確認されている。これらに対して、本実施例のように、重み係数マップありの物質情報による制約がある場合には、図４（ｅ）および図４（ｆ）ともに物質情報によるアーティスト低減効果を維持しながらも、微細構造が保持された再構成画像が得られることが確認された。これにより、本微細な構造が消失することもなく再構成画像を得ることができる。

【0061】

図５（ａ）および図５（ｂ）は、同一の画像からそれぞれ得られた逐次近似法による各々の再構成画像である。重み係数マップがない場合には、図５（ａ）では境界においてガタツキが発生し、図６の点線のプロファイルでは境界において急峻な段差が確認されている。これらに対して、本実施例のように、重み係数マップがある場合には、図５（ｂ）では境界においてガタツキがなくなり、図６の実線のプロファイルでは境界において段差が滑らかになっていることが確認された。これにより、本来滑らかであるべき物質情報が不自然にがたつくこともなく再構成画像を得ることができる。

【0062】

本実施例に係る画像再構成処理方法によれば、画像更新工程（図２ではステップＳ２）では、逐次近似法により画像（再構成画像μ）を更新し、重み係数マップ更新工程（図２ではステップＳ４）では、当該画像更新工程（ステップＳ２）での画像の更新によって得られた再構成画像μから、事前知識に対する重み係数マップＷを生成して、当該重み係数マップＷに基づいて各画素ｊに対する事前知識の重み係数Ｗ_ｊを調整することで重み係数マップＷを更新する。このように、画像の更新によって得られた（推定中の）再構成画像μから、事前知識に対する重み係数マップＷを生成して、当該重み係数マップＷに基づいて各画素ｊに対する事前知識の重み係数Ｗ_ｊを調整することで、高画質な再構成画像が得られないという問題点を解決することができる。つまり、各画素ｊに対する事前知識の重み係数Ｗ_ｊを調整することで、再構成画像μの画素に対して過度な制約を与えることを回避することができ、高画質な再構成画像μを得ることができる。

【0063】

なお、「課題を解決するための手段」の欄でも述べたように、特許文献２：米国特許第８，９５８，６６０号公報で提案されている手法は、ボクセル依存型スケーリング・ファクタを算出するための係数マップを目的関数の勾配に適用する。つまり、特許文献２では目的関数全体から算出される更新量に対して重み係数マップを生成する点で本発明の実施例と相違する。

【0064】

また、特許文献２では目的関数全体に対して重み係数を調整しようとすると、データ項にまで重み係数が掛けられ、データ項の更新量までが抑制されてしまい、逐次近似法における再構成の処理スピードが遅くなる。その結果、逐次近似法での繰り返し回数（反復回数）を増やさなければならない。これに対して、本発明の実施例の場合には重み係数マップＷを事前知識のみ（上記（４）式の妥当性項Ｒを参照）に適用するので、逐次近似法における再構成の処理スピードが速くなり、逐次近似法での繰り返し回数（反復回数）ｎを減らすことができる。

【0065】

本実施例では、重み係数マップとして、複数の物質が混合した画素である混合画素の情報を反映したマップを用いている。混合画素の情報を用いて当該情報を反映したマップを重み係数マップとすることで、混合画素に対する事前知識（物質情報）の効果が弱められる。その結果、物質情報によるアーティスト低減効果を維持しながらも、微細構造が保持された再構成画像（図４（ｅ）および図４（ｆ）を参照）や、滑らかな境界を持つ再構成画像（図５（ｂ）および図６の実線のプロファイルを参照）が得られる。

【0066】

具体的には、本実施例では再構成画像μのエッジ情報を用いて重み係数マップＷを生成する（図２のステップＳ４を参照）。例えば、Ｗ_ｊ＝ｅｘｐ（−γ×Ｅ_ｊ）のように、エッジ強度Ｅ_ｊ（ｊ：画素番号）が大きいほど、重み係数マップの値Ｗ_ｊが小さくなるように設定する。また、本実施例では、図２のフローチャートに示すように重み係数マップを更新するタイミングを画像更新毎としている。

【0067】

本実施例に係る画像再構成処理プログラム８Ａ（図１を参照）によれば、本実施例に係る画像再構成処理方法（図２のフローチャートを参照）をコンピュータ（本実施例では図１に示す再構成処理部７を構成するＣＰＵまたはＧＰＵ）に実行させることによって、再構成画像μの画素に対して過度な制約を与えることを回避することができ、高画質な再構成画像μを得ることができる。

【0068】

本実施例に係る断層撮影装置（本実施例ではＸ線ＣＴ装置）によれば、画像再構成処理プログラム８Ａを実行する演算手段（本実施例では図１に示す再構成処理部７を構成するＣＰＵまたはＧＰＵ）を備えることによって、再構成画像μの画素に対して過度な制約を与えることを回避することができ、高画質な再構成画像μを得ることができる。

【0069】

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

【0070】

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

【0071】

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

【0072】

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

【0073】

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

【0074】

（５）上述した実施例では、上記（４）式では１種類の事前知識（妥当性項Ｒ）を使用したが、複数の妥当性項（Ｒ_１，Ｒ_２，…）を導入した場合には、各妥当性項に共通する重み係数マップＷを生成してもよいし、各妥当性項に対応する異なる重み係数マップ（Ｗ_１，Ｗ_２，…）を生成してもよい。すなわち、（ａ）互いに異なる複数の事前知識（妥当性項Ｒ_１，Ｒ_２，…）に対して、同一の重み係数マップＷを適用してもよいし、（ｂ）互いに異なる複数の事前知識（妥当性項Ｒ_１，Ｒ_２，…）に対して、それぞれ異なる重み係数マップ（Ｗ_１，Ｗ_２，…）を適用してもよい。例えば、（ｂ）において、２つの妥当性項Ｒ_１，Ｒ_２に対して、それぞれ異なる重み係数マップＷ_１，Ｗ_２を適用する場合には、ｎ回目の繰り返しにおける再構成画像μ中の各画素μ_ｊの更新式は、下記（５）式で表される。β_１，β_２は、上記（１）式での説明でも述べたように妥当性項Ｒ_１，Ｒ_２の強さをコントロールする係数であり、経験的に決められることが多い。

【0075】

【数3】

【0076】

（６）上述した実施例では、重み係数マップを更新するタイミングを画像更新毎としていたが、重み係数マップを更新するタイミングは、一定間隔毎，一定の基準を満たすタイミングあるいは任意のタイミングであってもよい。

【0077】

（７）上述した実施例では、再構成画像のエッジ情報を用いて重み係数マップを生成することで、重み係数マップの生成を自動的に行ったが、重み係数マップの生成を手動で行ってもよい。例えば、再構成画像を表示（モニタリング）して、表示された再構成画像の（境界と認識される）任意の場所（画素）をオペレータ（ユーザ）が指定することで、任意の場所に手動で入力された情報を用いて重み係数マップを生成してもよい。

【0078】

（８）上述した実施例では、事前知識として物質情報を用いたが、事前知識であれば特に限定されず、例えば複数の画素値における滑らかを示す平坦情報などでもよく、物質情報と平坦情報とを組み合わせて用いてもよい。

【産業上の利用可能性】

【0079】

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

【符号の説明】

【0080】

１ … Ｘ線ＣＴ装置
７ … 再構成処理部
８Ａ … 画像再構成処理プログラム
Ｗ … 重み係数マップ
μ … 再構成画像
Ｅ … エッジ画像
Ｒ … 妥当性項（ペナルティ項）

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】