WO2008059982A1

WO2008059982A1 - Dispositif de reconfiguration d'image, procédé de reconfiguration d'image, programme de reconfiguration d'image, appareil de tomodensitométrie

Info

Publication number: WO2008059982A1
Application number: PCT/JP2007/072339
Authority: WO
Inventors: Yukihiro Nishikawa
Original assignee: National University Corporation Kyoto Institute Of Technology
Priority date: 2006-11-13
Filing date: 2007-11-13
Publication date: 2008-05-22
Also published as: EP2092891A1; EP2092891B1; EP2092891A4; KR20090079994A; CN101553171A; US8090182B2; US20090279768A1; JPWO2008059982A1; CN101553171B; JP5190825B2

Description

明細書画像再構成装置、画像再構成方法、画像再構成プログラム、 CT装置技術分野

本発明は、対象物の放射線投影像から断面像を再構成する技術に関する。背景技術

計算機トモグラフィー (CT)は、 X線などの光線を使って、様々な角度からの物体の投影像を撮影し、その後、計算によって断面像を得る技術である。図 1に典型的な X線 CT 装置の模式図を示す。

典型的な CT装置では、 X線光源が回転移動することにより、測定対象の様々な角度からの投影像を取得する仕組みになっている。こうして得られた投影像から計算操作によって断面像を得るのが CTである。通常は、 Filtered Back-projection法 (FBP)によつて投影像から断面像を変換操作によって求める。 FBPでは、投影像にフィルタ (基本的には微分フィルタ)を施した後、もともとの投影方向に「逆投影」という操作を施すことで断面像を得る。このとき、微分フィルタはノイズや誤差を増幅し、アーティファクト (本来存在しない誤差や虚像)を生みやすいという特徴がある。さらに、逆投影操作はそのァ一ティファクトを断面像全体に伝播させる効果がある。そのため、 CTにおけるアーティファクトは問題のある部分だけに留まらず、断面像全体を損なうという点で致命的になることが多い。

アーティファクトの多くは、 FBPに含まれるフィルタ操作や逆投影操作が原因なのだから、 FBPを用いなければ、アーティファクトのない断面像が得られるはずである。 FBP 以外で断面像を計算する方法として、 Algebraic Reconstruction Technique(ART)が歴史的に重要である。 ARTは FBPが考案される前から利用されるほど歴史のある方法である。 ARTでは断面像の計算をフィッティング問題と捉え、断面像をパラメータ、投影像をフィッティング対象として、断面像から計算した投影像 ( が、実験で求めた投影像 (P 0)に一致するように逐次的に断面像を修正する。 ARTでは (p-Po)が 0になるように断面像を漸近的に改良してゆぐが特徴となっている。 ARTは FBPを全く用いないかわりに、断面像を求める計算に時間がかかるため、今では特殊な用途 (地震波の解析など)にしか用いられない。 ARTは FBPほど極端なアーティファクトは現れないが、得られる断面像は FBPの方が自然であることが多い。

さらに、フィルタ操作や逆操作に由来しないアーティファクトの要因として、投影像におけるデータの欠落'過少も挙げられる。データが少なければ結果として得られる断面像の画質も低下するのは当たり前である。 FBPではデータの欠落 ·過少も致命的なァ一ティファクトを生み出すことが知られており、大きな問題と認識されている。フイツティングに基づく CT(ARTを含む)はデータの欠落や過少に対して FBPよりは強し、とされている。しかしながら、 CTにおけるデータの欠落は、極端に「条件が悪い」問題とされておリ、 ARTなどを用いても改善が得られにくい。 ARTでは (p-po)を誤差として考えるのであるが、条件が悪し、ときはフィッティングに失敗しやすし、。いわゆる二乗誤差 (p— Po)² の最小化を行った方がフィッティングとしては安定する。二乗誤差を最小化する方法としては、最小二乗法が最も標準的である。最小二乗法では、一辺がパラメータの個数だけある正方行列の逆行列を求める。 CTにおけるパラメータは断面像を構成するそれぞれの画素値であるため、パラメータの個数は莫大になる。 1000 X 1000画素の断面像だと、パラメータの数は 100万個になり、行列の要素数は 1兆個と莫大である。従つて、正攻法の最小二乗法では行列が大きくなりすぎで破綻してしまう。そこで、正攻法の取小二乗法 ¾：避け、 Simultaneous iterative reconstruction technique (S1R 八 ~» Iterative least squaretechnique (ILST)が考案された。 ARTと同様に断面像の計算をフィッティング問題として捉えるのであるが、先に述べたように、正攻法の最小二乗法を避けなければならないため、 SIRTも ILSTも計算の途中に投影操作の逆演算として FBP が利用される。そのため、 FBPに起因するフィルタ操作と逆投影にまつわる問題を根本的には解決しない。おそらくそれが理由となり、 SIRTや 1LSTでは各種アーティファク卜の「低減」効果が報告されているに留まる。

[非特許文献 1 ]

Yazdi Gingras L, BeaulieuL An adaptive approach to metal artifact reduction in helical computedtomography for radiation therapy treatment planning: experimental and clinicalstudies. Int J Radiat Oncol Biol Phys.62: 1224 - 1231， 2005. 発明の開示

発明が解決しょうとする課題

以上の議論から、アーティファクトの無い断面像を得るには、 FBPを用いない、データの欠落に対して強い、という性質が重要であることがわかる。そこで、 ARTと同様に FBPを用いず、フィッティングの評価関数として二乗誤差を利用しょうというのは、誰でも思いつくことである。にもかかわらず、 CTにおけるアーティファクトの問題が根本的に解決していないのは、この単純な考えの実現は簡単ではないからである。第一の問題は、断面像の計算はフィッティングの問題としては大規模すぎて、計算に時間がかかることである。適切な高速化が必須である。第二の問題は、最小二乗法はフィッティングアルゴリズムとしては「条件が悪し、」問題に対して弱いことである。そのため、最小二乗法に基づく既存のアルゴリズム (ILSTなど)の延長では、アーティファクトの問題を解決できないと思われる。第三の問題は、 SIRTや ILST力《そうであるように、既存のァルゴリズムでは FBPを完全に排除することができていなし、ことである。

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

本発明の第一の新規性は、フィッティング方法として Simulated Annealing法 (SA)を採用したことである。 SAは公知の技術であるが、フィッティングに時間がかかることが知られており、 CTと本質的に相性が悪し、。しかしながら、敢えて SAを CTに適用することで、 FBPを全く用いずとも、二乗誤差の最小化による断面像の計算ができるようになつた。また、 SAはデータの欠落など条件の悪しワイツティングに対してかなり安定であるという性質があり、この点においてもアーティファクトの根本的な解決に有利である。以上の事柄を鑑み、 SAを CTに適用することで、アーティファクトの根本的解決に効果を見出した点が、本発明の重要な新規性である。単純に SAの CTへの適用を考えたとき、断面像から投影像を求める計算を無数に繰り返すことになる。まともに計算すると、断面像から投影像を求める計算は、 FBPとほぼ同じであり、 SAを CTに適用すると、 FBPに比べて数百万倍の時間がかかってしまうことになる。現在の高速な計算機を以つてしても、「年 J単位の計算時間を要してしまう。そこで、本発明では、式変形によって計算量を大幅に低減することで、この問題を解決した。この工夫は、 SAを CTに適用する際に必要な技術であり、本発明の重要な要素である。第二の新規性は、フィッティングの評価関数として、二乗誤差だけでなぐアーティファクトを積極的に破壊するスムージング項とエントロピ一項を導入した点である。従来の CTでは、 pと p。との差あるいは二乗誤差だけを対象としていた。この場合、アーティファクトを御座なリにした方が良いフィッティング結果 (小さな誤差)が達成できるという可能性が常にあり、現実の多くのケースではそのような条件になっていた。すなわち、ァ一ティファクトを別のアーティファクトで相殺する状態になっていたのである。そのため、アーティファクトはなかなか消えないのである。本発明においても二乗誤差だけを評価関数とした場合には、アーティファクトは完全には消えなかった。この事実は、アーティファクトのない断面像を得るには、二乗誤差以外の要請が望ましいことを意味している。本発明では、統計力学の基本にのっとり、エントロピ一項やス厶ージング項を考案した。これらの項は、断面像は滑らかで自然な濃淡画像であるべきだ、という当たり前の要請を数式の形で表現したものである。エントロピ一項は、アーティファクトを破壊し、断面像全体の画質を均質にするように、フィッティングを誘導する。スムージング項は、ェントロピー項で誘導される断面像のザラツキを押さえる効果がある。これらの項を導入することで、アーティファクトが消え、かつ自然な断面像が得られるようになった。ェントロピー項とスムージング項は、それぞれ単独でもアーティファクトを低減する効果があるが、 2つを組み合わせるとさらに効果的である。

なお、本発明においては、一般の定義よりも広い概念として、 X線 ·可視光 ·電波を含む電磁波、電子や荷電粒子からなる粒子線、媒体の振動である音波等を全て含めて「放射線」と呼ぶことにする。

発明の効果

本発明の第一の効果として、データの欠落が要因になっているアーティファクトの低減効果が挙げられる。データの欠落が問題になる場合として、観察対象に不透明部位がある場合、投影角度に制限がある場合、投影角度の刻みが不均一である場合、コーンビームやへリカルスキャンなどの三次元 CTの場合が挙げられる。

特に効果があることがわかっているのは、観察対象に不透明部位がある場合に現れるメタルアーティファクトの低減効果である。メタルアーティファクトとは、観察対象中に X 線に対して不透明な部位 (多くの場合、金属部分)が存在する場合、 CTで得られる断面像全体 (不透明部位だけでない)が破壊的に乱されることを指す。メタルアーティファクトの原因は、不透明部分で投影像の輝度が不連続に変化すること、不透明部分の情報が欠落することである。 FBPにおいて微分フィルタが適用されると、輝度の不連続な変化は異常値となり、逆投影操作によって、不透明部位を中心とした筋状のアーティファク卜を生み出す。さらに、情報の欠落によって不透明部位とは直接関係のない部分に、予期しない明暗が生じる。 FBPを用いず、情報の欠落に対して安定な SAを利用した本発明力《メタルアーティファクトの除去に有効であろうことは想像に難くない。

実用上、重要なのは、コーンビームやへリカルスキャンへの応用であろう。両者は 3 次元 CTと呼ばれ、現在急速に普及しつつある技術である。しかしながら、 3次元 CTには特有のアーティファクトが現れることが知られており、その原因はわかっているのだが解決には至っていない。その原因は、コーンビームの場合には、データの欠落である。コーンビームでは完全な断面像を得るための条件を満たすことができず、データの欠落に由来するアーティファクトが現れる。ヘリカルスキャンでのアーティファクトの根本的な原因は、逆投影操作にある。ヘリカルスキャンでは系の幾何学的異方性 (らせん）がフィルタ操作'逆投影操作に影響し、風車アーティファクトが現れる。本発明は、フィルタ操作'逆投影操作を必要とせず、データの欠落にも強いので、 3次元 CTにおける諸問題を解決できる。

投影角度力《不均一である場合として挙げられるのは、地震波による地球内部の CT や、人工衛星からの電波を用いた大気状態の CTによる解析である。これらは FBPが利用できない典型として知られている。本発明を利用することで、解析精度の向上を期待できる。

本発明の第二の効果は、投影像撮影の迅速化、 X線被爆量の低減である。 SAはデータの欠落だけでなぐ過少にも安定であるという特徴があり、本発明も同様である。 C 丁の場合、データの過少に相当するのは、投影角度の数が少ない場合である。すなわち、本発明を利用すると、従来の CTより投影角度の数を少なくすることができる。投影角度の数とは放射線による投影像の撮影枚数である。撮影枚数は、撮影時間と被爆量に比例するので、撮影枚数が少なければ、撮影時間の短縮、被爆量の低減につながる。

また、データの過少の別のパターンとして、投影像の画質が不良 (S/N比が低い)な場合が挙げられる。本発明はそのような場合にも比較的安定であることがわかっている。投影像の画質の低下を許容できるとなると、これも撮影時間の短縮'被爆量の低減につながる。あるいは、 S/N比が極端に低い SPECT, PETにおいて本発明は画質の向上に寄与できるだろう。

本発明の第三の効果は、断面像の輝度値の決定精度が高いことである。本発明では測定した投影像をかなり忠実に再現する断面像が得られる。再現の精度は FBPよリ 2桁程度高し、。これは二乗誤差を最小化するというフィッティングアルゴリズムの恩恵である。輝度値の決定精度が高いことで、輝度値の定量性が保証され、輝度値を用いた密度測定が可能になる。これは骨密度測定などにおける測定精度の向上につながる。また、異常のある部位 (腫瘍のある臓器など)の検出精度の向上にも寄与するだろラ。

本発明の第四の効果は、本発明で得られる断面像のコントラストが FBPよりも高いことである。第三の効果で述べたように、本発明は輝度値の決定精度が高いのであるが、その副次的な効果として、断面像のコントラストが高くなるのである。コントラストが高いと見掛けの空間分解能も高くなる傾向がある。その結果、本発明で得られる断面像は、従来のものより高画質になる。この効果は、特殊な条件 ·用途の CTではな普通の C Tにも本発明が役立つという点で、特筆すべきものである。図面の簡単な説明

図 1は、 X線 CT装置の模式図である。

図 2は、断面像 f(_X，y)、投影像 p(r， Θ )の関係を示す模式図である。

図 3は、典型的な p。(r, Θ )の例 (横軸に検出器のチャンネル位置、縦軸に投影角度をとつた画像)である。

図 4は、本発明の基本手順を示すフローチャートである。

図 5(a)は、手順（ I Γ(νΐ)のフィッティング経過の模式図である。図 5(b)は、手順 (1 )〜(9) のフィッティング経過の模式図である。

図 6は、 [非特許文献 1 ]に掲載されている手法での結果を示す図である。

図 7は、図 6(a)からシミュレーションによって求めた sinogramである (白い部分は不可視領域 )。

図 8(a)は、図 6(c)をさらに拡大した図である。図 8(b)は、本発明での結果を示す図である。

図 9は、仮想エネルギー Eの各項の効果を表した模式図である。

図 1 0は、角度制限がある場合のアーティファクトと本発明でのアーティファクト低減効果を示す図である。発明を実施するための最良の形態

(第 1の実施形態）

本発明の実施形態による X線 CT装置の概略構成を図 1に示す。この X線 CT装置は、撮影部とコンピュータとを備えている。撮影部は、 X線を測定対象に照射するための光源と、測定対象を透過した X線を検出する検出器とを備えており、測定対象の周囲を多数の方向から X線により投影して投影像を取得する。コンピュータは、 X線 CT装置全体を制御する制御部と、撮影部により得られた X線投影像に基づいて測定対象の関心領域の断面像を生成する画像再構成処理部とを備えている。なお、図 1に示した構成はこれ以降の各実施形態に共通のものであり、各実施形態では画像再構成処理部におし、て行われる処理がそれぞれ異なってし、る。

本実施形態の画像再構成処理部は、投影像から断面像を求めるためのフイツティング方法として Simulated Annealing法 (SA)を採用している。そこで、まず、 Simulated Annealing法 (SA)のフレームワークを説明する。 SAは Monte Carlo法から派生したフィッティングアルゴリズムで、乱数に基づいてフィッティングを進める点と、熱力学を模して仮想エネルギーや仮想温度を取り扱う点、が特徴となっている。 SA自体は公知の技術であり以下のステップ（ i )〜(vi )により行われる。

( i )パラメータを乱数に基づいて 1つ選択し、さらに乱数に基づいてそのパラメータを変更するほ L数)。

( ii )変更後の評価を行う。評価関数として仮想的なエネルギー Eを考える。通常の SA では Eは二乗誤差の総和である。変更前後の Eの変化を Δ Eとする (評価)。

( iii )評価が改善される (Δ Eく 0)なら、その変更を受け入れる (変更)。

( iv ) exp (- Δ E/T)の確率で改悪を受け入れる。

( V )Tを少しずつ減じる。

( vi )( i )から繰り返す。 SAでは、（iv)にあるように、ボルツマン統計に従って改悪を受け入れることで、局所解から脱出できる可能性を担保している。そのため、局所解にとらわれることな大域解に迪リつくことができ、条件の悪いフィッティングに対しても安定に機能する。また、 T を徐々に減じることで、大域解にソフトランディングするようになっている。（ i )から ( V ) まで一通り実行することを 1モンテカルロステップと呼ぶ。 SAではモンテカルロステップを無数に繰り返すことで、フィッティングが進行する。計算に必要な時間は、 1モンテ力ルロステップに必要な時間 X必要なモンテカルロステップ数である。必要なモンテ力ノレ口ステップ数は、パラメータの数'自由度に比例する。

次に、請求項に沿った形で本実施形態の説明を行う。 CTをフィッティング問題と考えた場合、フィッティングパラメータは断面像 (f0，y))である。フィッティング対象となるデータは、測定した投影像 (p。(r， Θ ))である。 rは投影像を撮影した 1次元検出器のチャンネル位置、 Θは投影角度である。座標の定義を図 3に示す。 p₀(r, Θ )は、角度 Θを変えながら投影像を測定することで得られるデータセットであり、横軸に r、縦軸に Θを取ると、 2次元の画像とみなすことができる。このようなデータを sinogramと呼ぶ。典型的な sinogram ¾：図 3に示す。

極論すれば、 CTは sinogramから断面像を求める技術である。本実施形態では、暫定の断面像 f ,y)と測定した投影像 P。(r, Θ )との二乗誤差を考えるが、二乗誤差を計算するために、数 1 01を使って、 f(_X,y)から計算で投影像 (p(r, Θ ))を求める。

[数 101 ] pyr,め = f(r cos θ - s sm 6, r sm e + s cos Θ)

s

数 101は図 2において、 sの方向に和をとつておリ、 f(x，y)の s方向の投影像の計算であることがわかる。こうして求めた p(r，0 )を用いて、二乗誤差 Hは以下のように計算でさる。

[数 102]

通常のフィッティングでは、仮想エネルギー Eとして、数 1 02をそのまま利用するが本実施形態では、 H以外にスムージング項とエントロピ一項を導入している点が特徴となっている。 Eの定義は次の通りである。

[数 103]

E = H - TS ^ ca

Tは仮想温度 (温度パラメータ)、 Sはエントロピー、 σは画素値の標準偏差、 cはスム一ジング項の強さを表す係数である。 TSがエントロピ一項、 cびがスムージング項である。 Sや σの定義'計算法は後述する。本実施形態では、これらの数式を使い、図 4に示すような手順で断面像を計算する。

'ステップ (a):何らかの形で求めた暫定の断面像 f(x,y)に対し、評価関数 (以下「エネルギー」と呼ぶ XE。；)を求める (ST1 0および ST20)。

■ステップ (b):乱数で (x。,y。：)および△〃を選択し、断面像 f(x,y)の一部を改変する (ST3 0)。

'ステップ (c):改変後の断面像 f y)に対して、評価関数を求める (ST40および ST5

0)。

■ステップ (d)：前記エネルギー (E。)と前記エネルギー (E との差 ( Δ E)を求める (ST60)。 'ステップ (e):前記改変を採用するか否かを前記差 (Δ Ε)および温度パラメータ (T)を用いた受理関数 (典型的にはボルツマン統計)により判定する (ST70〜ST1 10)。

■ステップ (f)：前記ステップ (a)に戻る (ST1 20)。

'ステップ (g):前記ステップ (a)〜(わの繰り返しが所定回数に達するごとに仮想温度 (T) の値を変更する (ST1 30)。

'ステップ (h):前記ステップ (e)における判定結果が所定の停止条件を満たしているか否かを判定し、満たしていれば処理を終了する。ここでは、 ST80.ST1 00が実行された場合を「成功」、 ST1 10が実行された場合を「失敗」とし、成功の確率が 1 0%にの値は適宜調整可能)を下回ると終了する。そして、終了時点において採用されている推定断面像を測定対象の断面像とし、これをコンピュータのディスプレイに表示したり記録媒体に記録したりする。

ステップ (a)からステップ (h)の一連の操作が請求項 1 2に対応し、本実施形態ではこれら一連の操作を図 1の画像再構成処理部が行う。 SAの基本手順（ i )〜(vi )との対応を考えると、ステップ (b)が（ i )、ステップ (aXcXd) が（ ϋ )、ステップ (e)に (iii )と ( iv)が含まれ、ステップ (g)が ( v )、ステップ (f)が (vi )にそれぞれ対応する。このこと力、ら、本実施形態は CTに SAを忠実に当てはめたものであることがわかる。

なお、上記ステップ (h)における処理終了の判定ついては、図 1の画像再構成処理部において行うのではなぐコンピュータのディスプレイに推定断面像を逐次表示し、これをユーザが目で見て処理の終了をコンピュータに指示できるようにしてもよい。

(第 2の実施形態）

数 1 03に示された Eを計算するためには、数 101と数 1 02を経ることになリ、 s, r, Θ に関して級数和を実行することになる。 3重の積分 (級数和)を計算するので、大変時間がかかる。すなわち、 1モンテカルロステップに必要な計算時間が長し、ことを意味する。さらに、 CTではパラメータの数が莫大であり、その結果、 SAの実施に必要な総計算時間は、現在の計算機でも「年単位」になってしまう。

そこで、本実施形態では、 Eを計算するのではなぐ変更を行ったときの変化分 Δ Εだけを主に計算する。

[数 1 04]

AE = AH + cAa - TAS

今、暫定断面像 f ,y)に関する変更を A f ，y)とする。 A fO，y)は座標 (x₀,y₀)においてのみ△ /の値を持ち、それ以外は 0であるような断面像である。その ΔΚχ,γ)の投影像△ p(r, Θ )は数 1 01と同じ方法で計算できる。この△ p(r, Θ )を用いて、△ Hは以下のように計算できる。

[数 105].

Δ / = X { , Θ) + Ap(r, Θ) - p₀ (r,め }²―∑ {p(r, θ) - p₀ (r , θ)}

τ,θ τβ これを書き下して、

[数 106]

となる。今、 Af(x，y)は (x₀,y₀)においてのみ値を持つので、 Δ ρ(ι·，は r( 0 )=x_ocos 0 + y。sin Θにおいてのみ値 Δ μを持ち、それ以外は 0である。従って、数 1 06の {}の中身も、 r( 0 )=x₀cos Θ +y₀sin Θにおいてのみ値を持つ。そのため、級数和は rと Θの両方に対して行う必要はな以下のように表すことができる。

[数 1 07]

数 1 07では Θに関する級数和だけになつている点が重要である。これにより、計算量を大幅に削減することができる。 Pや Poはデジタル画像であるため、数 107の計算にあたっては、通常は、 K に関する内挿を行う必要がある。そのため、数 107の {}内をさらに展開することはできない。しかしながら、誤差を容認して数 1 07を展開すると次式が得られる。

[数 108]

ΑΗ = ΜΑμ² + — 。

( , } ただし、 Mは投影角度の数である。数 1 07の代わりに数 1 08を用いると Δ Hの値に 1 %程度の誤差が生じる。しかしながら、数 1 08は数 107より高速であり、本発明での利用価値は高い。

次に△ σの計算を説明する。 σは座標 (_Xo，y₀)付近の輝度値の標準偏差である。 (x₀,y 。)付近として、（x₀,y₀)の周囲 d X d画素を考える (実施例では d=5を用いている)。それらの画素の標準偏差は次式で求めることができる。

[数 109]

[数 1 10]

〈/(。， o )²〉 = 。 + i, y₀ +ゾ) ²

[数 1 1 1 ]

(Λ^χο，ァ o )〉 = ∑ Λ^χο + 。 + j)

" i,j=-d l 2 である。数 109〜 1 1 1を用いると、厶 σは次のように計算することができる。

[数 1 12]

Δσ = / 。，ァ。 )²〉 7 ¹- - {〈ル。，ァ。 )〉— ~ ただし、は変更前および後の f(x₀,y₀)の値である。

次に ASの計算であるが、その前にエントロピー Sの定義を行う。一般に計算機が取リ扱う画像はデジタル画像であり、画素の座標 (x,y)だけでなぐ画素の値もデジタルの値となっている。そこで、 1つの画素を 1つの量子と考え、画素値を量子状態とみなすことにする。すると画像は複数の量子からなる系であると考えることができる。統計力学に則ると、そのような系に対するエントロピ一は次のように定義される。

[数 1 13]

N\

S 二 k\n ―

Nは画像中の画素の総数であり、は画素値がデジタル値で iである画素の総数である。 kは通常の物理ではボルツマン定数であるが、本発明は物理とは直接関係がないため任意である。本発明では、 σと同様に Sも (x₀，y₀)の周囲 d X d画素で定義されるとした。

数 1 13で定義された Sに対して A Sを考える。変更によって、画素値はデジタル値で i から jに変更されるとする。すると、 Δ Sは次式のようになる。

[数 114]

Ni NI

AS = k\n kin

N,!^!---^,. -l)---(N +l)---N„! N,!N₂!-"N … Nゾ! "'N,,! すなわち、変更によって、が 1減って、 Njが 1増えるのである。数 114を展開すると、非常に簡単な式が得られる。

[数 115]

AS = k\nN_i -kln(N_J +l) 以上をまとめて、本実施形態の手順を書き下すと次のようになる。

( I )乱数を使って (x₀,y₀)および△ を選択する。

(11)数104、108(または107)、"2、"5を使って厶巳を計算する。

(ΠΙ)評価が改善される (ΔΕく 0)なら、 f(x₀,y₀)に△ を足す。

(IV) ΔΕ>0の時にも、 exp (— ΔΕ )の確率で、 fx₀,y₀)に△ を足す。

(V) Tを少しずつ減じる。

(VI) (I)から繰り返す。

(W)所定の停止条件を満たしているか否かを判定し、満たしていれば処理を終了する。ここでは、 (Π),(ΐν)において△ が足された場合を「成功」、それ以外を「失敗」とし、成功の確率が所定値 (例えば 10ο/ο、この値は適宜調整可能)を下回ると終了する。そして、終了時点において採用されてし、る推定断面像を測定対象の断面像とし、これをコンピュータのディスプレイに表示したり記録媒体に記録したりする。

手順（ I )〜(VI)は SAの基本手順（ i )〜(vi)を忠実に踏襲してし、ることがわかる。なお、手順（ I )〜(VH)の処理は図 1の画像再構成処理部において行われる。

なお、上記手順 (W)における処理終了の判定ついては、図 1の画像再構成処理部において行うのではなコンピュータのディスプレイに推定断面像を逐次表示し、これをユーザが目で見て処理の終了をコンピュータに指示できるようにしてもよい。

(第 3の実施形態）

次に、上で述べた方法の重要な変形について述べる。上で述べた方法では数 107あるいは数 1 08を利用し、 Θに関する級数和を求めている。そのため、画素 1つの変更を検討するのに M回の繰り返し演算を要する。本実施形態はこの部分の計算量をさらに減らすものである。

まず、 p(r, Θ )の逆投影像 g(x，y)を考える。

[数 1 1 6] g(x, y t = > p(x cos θ + y sm θ, θ) 数 1 1 6では、フィルタ操作をしていないため、 g(x,y)は断面像 f(x,y)に似ているが、ぼやけたような画像になる。 g(x,y)を用いると、数 1 08は次のようになる。

[数 1 1 7]

Δ = ΜΑμ + Δ〃 {g(x₀ , 。)一 g₀、 ₀ , 0 )} ただし、 g。(x,y)は p₀(r, Θ )の逆投影像で、数 1 1 6と同様の方法で計算される。数 1 1 7では繰り返し演算がないという点で数 1 08より優れている。 g₀(x,y)は計算過程で全く変化しないため、予め計算しておくことができる。一方、 g(x，y)は f(x,y)の 1点を変更する度にわずかずつ変化するため、厳密には（ I )〜(IV)の手順を 1回経るごとに再計算が必要である。しかしながら、 f(_X,y)の変更に伴う g(x,y)の変化は少ないという仮定をすると、別の手順を適用することができ、それが本実施形態である。以下に手順を示す。なお、以下の手順 (1 )〜(9)の処理は図 1の画像再構成処理部において行われる。

("I ) _Po(r，から g₀(x，y)を求める。

(2) f(_x,_y)から p(r， Θ )を求め、さらに p(r, Θ )から g(x，y)を求める。

(3) f(_X,y)に対する変更値を画素値とする画像厶 μ (x,y)を乱数で生成する。

(4)各画素値について数 1 1 7を適用し、 A H(x，y)という画像を計算する。

(5)同様に各画素値について数 1 1 2、数 1 1 5を適用し、 Δ σ(χ，ν)、 A S(x，y)という画像を計算する。

(6)数 1 04に基づいて A E(x，y)を計算する。

(7) Δ Eが正である座標 (x,y)について、 Δ (x,y)を 0にする。

(8) f(_X,y)に A jU (x,y)を足す。

(9) Tをひ倍 ( く 1 )し、（2)から繰り返す。 (10)所定の停止条件を満たしているか否力、を判定し、満たしていれば処理を終了する。ここでは、（7)において△ (x,y)が 0にされた場合を「失敗」、されなかった場合を「成功」とし、成功の確率が所定値 (例えば 10%、この値は適宜調整可能)を下回ると終了する。そして、終了時点において採用されている推定断面像を測定対象の断面像とし、これをコンピュータのディスプレイに表示したり記録媒体に記録したりする。

なお、上記手順 (1 0)における処理終了の判定ついては、図 1の画像再構成処理部において行うのではなぐコンピュータのディスプレイに推定断面像を逐次表示し、これをユーザが目で見て処理の終了をコンピュータに指示できるようにしてもよい。

手順（ I と、手順 (1 )〜(1 0)を図 5に模式的に示した。図 5ではパラメータの値

(本件では画素の値)を軸とし、仮想エネルギー Eを濃淡で表してし、る (Eが小さいほど濃い)。図示の都合により、パラメータの数 (画素の数)は 2個の場合を取り扱っている。図 5(a)に示す手順（ I )〜(VE)では、ジグザグに折れ曲がりながら、時にでたらめな方向に進みながら、最終的には Eの最小値を探索してゆく。一方、図 5(b)に示した手順 (1 )〜 (1 0)では、 Eの勾配を考慮しながら、 Eの最小値を探索してゆく。図 5(b)の方が効率が良いということが直感的に理解できるだろう。数 1 1 7による計算量の低減及び、図 5(b) に示した探索効率の向上により、本実施形態は計算速度的に有利である。手順 (1 )〜 (1 0)は SAの基本アルゴリズム（ i )〜(vi )を忠実に踏襲するものではなし、が、 ( I )〜 (VII)の手順の自然な発展形である。

(その他の実施形態）

上述の各実施形態で説明した処理を実行する画像再構成処理部は、これらの処理をコンピュータに実行させるためのプログラム、当該プログラムがインストールされたコンピュータ、さらには、これらの処理を実行する専用 LSIなどを使って実現可能である。

(実施例）

実施例としては、シミュレーションを用いたメタルアーティファクト除去効果を取り扱う。まず、比較のために、 [非特許文献 1 ]に掲載されている手法での結果を図 6(非特許文献の fig.5)に示す。

図 6において、左のカラム (a)(d)(g)は Metal Artifactがないオリジナルのファントム画像である。中央 (b)(e)(h)は、オリジナルのファントムの所定の位置に不可視領域を設定し、通常の FBP法により再構成した画像である。右 (cXfXi)は、上記と同様に不可視領域を設定し [非特許文献 1 ]の手法により再構成した画像であり、従来技術で最もうまく Metal Artifactが除去できている例の 1つである。

そこで、 [非特許文献 1]から図 6(a)の画像を取り込み、文献と同じ位置に不可視領域を設定し、シミュレーションで投影像を計算した。結果を図 7に示す。図 7が本実施例における p。(r，0 )に対応する。図 7について、本発明によって再構成した結果および、比較のための図 6(c)の拡大図を図 8に示した。効果は一目瞭然である。本実施形態 (b)では Metal Artifactの特徴を見つけ出すことすら困難である。

ただし本実施形態 (b)ではエッジ付近が若干ボケている。これは「断面像を滑らかにする因子」 (スム一ジング項 οσ)がやや強いためである。 Metal Artifactを低減するには、その因子をやや強めに設定する必要があることがわかっている。このように本実施形態によれば、仮想エネルギー Eを計算する際の各因子のバランスには改良の余地があるものの、 Metal Artifact除去の効果は非常に高い。

なお、参考までに、スムージング項、エントロピ一項 TSを設定せずに本実施形態のアルゴリズムを実行した場合の結果を図 9に示す。図 9(a)はスムージング項 C CTなし (エントロピ一項 TSのみ)の場合の結果、図 9(b)はエントロピ一項 TSなし (スム一ジング項 C CTのみ)の結果、図 9(c)はエントロピ一項 TS、スムージング項ともになしの場合の結果を示す。図 9と図 8(b)とを比較すると、本実施形態においては、エネルギー Εの因子としてスムージング項 c σとエントロピ一項 TSの両方を設定したほうが、いずれか一方のみを設定した場合、ともに設定しない場合よりも Metal Artifact除去の効果が高いことが分かる。スムージング項だけでは、不透明を中心に放射状に広力《る縞模様がアーティファクトして残るという特徴がある。一方、エントロピ一項のみでは、断面像がざらつき S/Nが低いという特徴がある。また、最終的な画質は、 (徐冷過程におけるスムージング項の係数 cと仮想温度 Tの比が重要であるようである。

次に、回転角度制限がある場合の実施例を図 1 0に示す。図 1 0(a)は、 FBPのみ (通常の CT)による再構成結果を示す。図 1 0(b)は、 ILST (最小二乗法に基づく画像復元法) を適用した結果を示す。図 1 0(c)は、本実施形態による再構成結果を示す。これらを比較すると、角度制限時に現れるアーティファクトに対して本実施形態が有効であることが分かる。角度制限によるアーティファクトには、円状の領域が両端の尖ったァ一モンド状の領域になり、アーモンド状領域の尖っていない胴回りの輝度が反転するといぅ特徴がある。本実施形態によると、それら両方の特徴が低減される。産業上の利用可能性

本発明は、 Metal Artifactに対して顕著な効果があるので、 Metal Artifactが深刻な分野、たとえば歯の CT、金属製インプラン卜を含む CTに特に有用である。

また、本発明は情報欠落のある投影像のセットに対して一般的に有効である。情報欠落が顕著な場合の一つに、投影角度制限がある場合を挙げることができる。投影角度制限が問題になる分野は、三次元電子顕微鏡法、マンモグラフィーの CT版、平行移動 CT法 (特開 2006-71472号公報)などである。情報欠落が問題になる別の場合として、コーンビーム CT、ヘリカルスキャン CTなど 3次元 CTが挙げられる。本発明は、 3次元 C 丁に現れるアーティファクトを除去する目的にも有効である。

さらには、情報が極端に少ない系に関しても利用可能である。例えば、蛍光 X線 C丁、地震波による地球内部の CTなどに有用である。

また、本発明の副次的な効果として、断面像における輝度値 (X線では線吸収係数に対応する)を従来法よりも高い精度で決定できることが挙げられる。この効果は骨密度測定等における精度改善に応用可能である。

本発明を用いると、従来の方法よりコントラストの高い再構成像を得ることができる。そのため、アーティファクト等が問題とならない通常の X線 CTにおいても、本発明の利用価値は高い。また、本発明はデータの過少に対しても安定であることから、投影像測定の時間短縮、ひいては X線被爆量の低減に有効である。これらをあわせて考えると、本発明は、既存のすべての X線 CT技術を置き換える可能性を秘めている。

Claims

請求の範囲

1 . 対象物に放射線を照射して得られた投影像 (以下「放射線投影像」と呼ぶ)から前記対象物の断面像を求める装置であって、

前記対象物の現在の推定断面像から演算により得られた投影像と前記放射線投影像との差を含む評価関数 (以下「エネルギー」と呼ぶ )(E。)を求める手段 (a)と、

前記現在の推定断面像の一部を改変する手段 (b)と、

前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差を含むエネルギー ( )を求める手段 (c)と、

前記エネルギー (E。；)と前記エネルギー ( )との差 (Δ Ε)を求める手段 ( と、前記改変を受理するか否かを前記差 (Δ Ε)と受理確率を制御する温度パラメータ (T) とを用いた受理関数に基づいて判定し、判定結果を前記現在の推定断面像に反映させる手段 (e)と、

前記手段 (a)〜(e)による一連の処理の繰り返しが所定回数に達するごとに前記温度パラメータ (T)の値を変更する手段 (わと、

を備えることを特徴とする画像再構成装置。

2. 請求項 1において、

前記手段 (a)は、

前記現在の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記現在の推定断面像の局所領域の標準偏差との和を含むエネルギー (E。：)を求め、

前記手段 (c)は、

前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記改変後の推定断面像の局所領域の標準偏差との和を含むエネルギー (E，) を求める、

ことを特徴とする画像再構成装置。

3. 請求項 1において、

前記手段 (a)は、

前記現在の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記現在の推定断面像の局所領域のエントロピーとの和を含むエネルギー (E。：)を求め、

前記手段 (c)は、

前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記改変後の推定断面像の局所領域のエントロピーとの和を含むエネルギー (E ，)を求める、

ことを特徴とする画像再構成装置。

4. 請求項 1において、

前記手段 (a)は、

前記現在の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記現在の推定断面像の局所領域の標準偏差と、前記現在の推定断面像の局所領域のエントロピーとの和を含むエネルギー (E。)を求め、

前記手段 (c)は、

前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記改変後の推定断面像の局所領域の標準偏差と、前記改変後の推定断面像の局所領域のエントロピーとの和を含むエネルギー (Ε,)を求める、

ことを特徴とする画像再構成装置。

5. 請求項 1において、

前記手段 (a) , (c) , (d)に代えて、

[数 1]により Δ Ηを算出し、算出した Δ Ηを構成要素として含む Δ Ε(Δ Ε= Δ Η十…)を求める手段 (h)を備える、

[数 1 ]

ΑΗ ρ₀ ( ), θ)]}

ここで、前記対象物の現在の推定断面像を f(x,y)、前記手段 (b)による改変部分を△ f(x， y)とすると、 Af(x,y)は座標 (Xo，y₀)においてのみ厶 μの値を持ち、それ以外は 0であるような断面像である。 p(r, Θ )は、前記対象物の現在の推定断面像から演算により得られた投影像である。 P₀(r, Θ )は、前記対象物の放射線投影像である。 rは、投影像を撮影した 1次元検出器のチャンネル位置である。 Θは、投影角度である。 Θ )=x₀cos Θ + y₀sin Θである。

ことを特徴とする画像再構成装置。

6. 請求項 1において、

前記手段 (a), (c), (d)に代えて、

[数 2]により AHを算出し、算出した ΔΗを構成要素として含む ΔΕ(ΔΕ=ΔΗ十…)を求める手段 (h)を備える、

[数 2]

ΑΗ = ΜΑ ²+Αμ^{ρ(Γ(θ),θ)-ρΜ )}

Θ ここで、前記対象物の現在の推定断面像を f(x,y)、前記手段 (b)による改変部分を ΔΚχ, y)とすると、 Af(x,y)は座標 (x₀，y₀)においてのみ△ μの値を持ち、それ以外は 0であるような断面像である。 p(r, Θ )は、前記対象物の現在の推定断面像から演算により得られた投影像である。 P。(r， Θ )は、前記対象物の放射線投影像である。 rは、投影像を撮影した 1次元検出器のチャンネル位置である。 Θは、投影角度である。 Θ )=x₀cos Θ + y。sin Θである。 Mは、投影角度の数である。

ことを特徴とする画像再構成装置。

7. 請求項 5または 6において、

前記手段 (h)は、

[数 3]により△ σを算出し、算出した Δ σに係数 cを掛けたもの (〇△ σ)と前記 ΔΗとの和を構成要素として含む ΔΕ(ΔΕ=ΔΗ+(ϊΔσ + · · ·)を求める、

[数 3]

なお、 σは座標 (x。,y。)の周囲 dxd画素の輝度値の標準偏差であり [数 4]により表される。はそれぞれ前記手段 (b)による改変前後の f(x₀，y₀)の値である。

[数 4]

び ⁼ ( V 0， 0 )²〉—〈/( 。，ァ 0 )) ただし、

[数 5] 。 ⁺ ，ァ。 +力

[数 6] ( 。 ⁺ ，。+

である。

ことを特徴とする画像再構成装置。

8. 請求項 5または 6において、

前記手段 (h)は、

[数 7]により ASを算出し、算出した ASに温度パラメータ (T)を掛けたもの (一 TAS)と前記厶 Hとの和を構成要素として含む ΔΕ(ΔΕ=ΔΗ— TAS + " を求める、

[数 7]

AS = k\nN_i -k\n(N_j +l) なお、 Sは座標 (x。,y。：)の周囲 d xd画素の局所領域画像のエントロピーであり [数 8]によリ表される。

[数 8] S二 In ^:

N .N₂ \' ' . Ni、 ' ' N_n\ また、

N：前記局所領域画像中の画素の総数、

：画素値がデジタル値で iである画素の総数、

Nj：画素値がデジタル値で jである画素の総数、

k:定数、

であり、前記手段 (b)による改変により画素値がデジタル値で iから jに変更されるものとする。

ことを特徴とする画像再構成装置。

9. 請求項 1において、

前記手段 (e), (f)に代えて、

前記改変を受理するか否かを前記差 (Δ Ε)と受理確率を制御する温度パラメ一タ (T) とを用いた受理関数に基づいて判定し、判定結果を前記現在の推定断面像に反映させることを予約する手段 (e1 )と、

前記手段 (a)~(d)および (e1 )による一連の処理の繰り返しが所定回数に達するごとに、前記手段 (e1 )における予約を前記現在の推定断面像に反映させ、前記温度パラメータ (T)の値を変更する手段 (f1 )と、

を備えることを特徴とする画像再構成装置。

1 0. 対象物に放射線を照射して得られた投影像から前記対象物の断面像を求める装置であって、

前記対象物の放射線投影像 p。(r, Θ )の逆投影像 g₀(x,y)をフィルタなしの逆投影操作により計算する手段 (ml )と、

前記対象物の現在の推定断面像 f(_X,y)から投影像 p(r, Θ )を計算し、さらに当該投影像 p(r, Θ )の逆投影像 g(x，y)をフィルタなしの逆投影操作により計算する手段 (m2)と、前記対象物の現在の推定断面像 fO ,y)に対する変更値を画素値とする画像 Δ (x,y) を生成する手段 (m3)と、

各画素値について [数 9]を適用して、画像 A H(x，y)を生成する手段 (m4)と、

[数 9]

Δ = ΜΑμ ² + AM{g{x₀ , y₀ ) ~ go (。，。 )} ここで、

Μ :投影角度の数、

である。

前記△ H(x,y)を用いて△ E(x，y)を計算する手段 (m5)と、

ここで、

Δ E(x,y)は、評価関数 E。(x，y)と Ε^χ,γ)との差である。

E。(x,y)は、前記推定断面像 f x,y)から演算により得られる投影像 P(r, Θ )と前記放射線投影像 P。(r, Θ )との差を含む評価関数である。

^Ocy)は、前記推定断面像 f(_X,y)に前記手段 (m3)により得られた画像△ μ (x，y)を足し合わせたもの {f(_X，y)+ Δ μ (x,y)}から演算により得られる投影像 ip(r, Θ )+ Δ p(r, 0 )}と前記放射線投影像 P。(r, Θ )との差を含む評価関数である。

前記 Δ Εが正である座標 (x，y)について、前記△ μ (x,y)を 0にする手段 (m6)と、前記推定断面像 f ， y)に前記手段 (m6)により得られた画像△ (x，y)を足し合わせたものを新たな推定断面像 f(x，y)として前記手段 (m2)〜(m6)による処理を繰り返す手段

(m7)と、

を備えることを特徴とする画像再構成装置。

1 1 . 対象物に放射線を照射して得られた投影像から前記対象物の断面像を求める装置であって、

前記対象物の放射線投影像 P。， Θ )の逆投影像 g。(x,y)を [数 1 0]により求める手段 (ml )と、

[数 1 0] g₀ (x, y) = Z Po ^x cos <9 + sin θ, θ) ここで、

r:投影像を撮影した 1次元検出器のチャンネル位置、

0 :投影角度、

である。

前記対象物の現在の推定断面像 f(x,y)から演算により投影像 P(r, Θ )を求め、さらに当該投影像 P(r,0)の逆投影像 g(x，y)を [数 11]により求める手段 (m2)と、

[数 11]

g(x,y)= > p(x cos Θ + ysm θ, θ) 前記対象物の現在の推定断面像 f(x,y)に対する変更値を画素値とする画像 Δ μ (x,y) を生成する手段 (m3)と、

各画素値について [数 12]を適用して、画像 AH(x,y)を生成する手段 (m4)と、

[数 12]

Δ/ = ΜΑμ ² +AM{gyx₀,y₀)-go ( 。，。 )} ここで、

M:投影角度の数、

である。

各画素値について [数 13]を適用して、画像△び0^)を生成する手段^5)と、

[数 13]

△^{σ =}、〈/(。， 0).

なお、 σは座標 (x。,y。)の周囲 dxd画素の輝度値の標準偏差であり

[数 14]により表れる。はそれぞれ変更前後の f(x₀,y₀)の値である。

[数" I4]

[数 15] (f(x₀， 0 〉二：^ Σ , ( 。 ^{+ Ζ}·，。 +ゾ

" i,j=-d/2

[数 16]

= - Ι ∑/( 。 "，ァ。 +ゾ)

« ゾ=一である。

各画素値について [数 17]を適用して、画像 AS(x,y)を生成する手段 (m6)と、

[数 17]

なお、 Sは座標 (x。,y。：)の周囲 d d画素の局所領域画像のエントロピーであり [数 18]により表される。

[数 18]

S = k\n

Λ^! ₂！… !… ！また、

Ν:前記局所領域画像中の画素の総数、

Ni：画素値がデジタル値で iである画素の総数、

Nj：画素値がデジタル値で jである画素の総数、

k:定数、

であり、画素値がデジタル値で ίから jに変更されるものとする。

[数 19]に基づいて AE(x,y)を計算する手段 (m7)と、

[数 19]

AE = AH + cAa-TAS C :係数

T:仮想温度 (温度パラメータ）

である。

前記 Δ Εが正である座標 (x，y)について、前記△〃(x,y)を 0にする手段 (m8)と、前記推定断面像 f(x,y)に前記手段 (m8)により得られた画像厶 μ (x,y)を足したものを新たな推定断面像 f(x,y)とする手段 (m9)と、

前記 Tをひ倍 ( く 1 )し、前記手段 (m2)〜(m9)による処理を繰り返す手段 (ml 0)と、を備えることを特徴とする画像再構成装置。

1 2. 対象物に放射線を照射して得られた投影像 (以下「放射線投影像」と呼ぶ)から前記対象物の断面像を求める方法であって、

前記対象物の現在の推定断面像から演算により得られた投影像と前記放射線投影像との差を含む評価関数 (以下「エネルギー」と呼ぶ )(E₀)を求めるステップ (a)と、前記現在の推定断面像の一部を改変するステップ (b)と、

前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差を含むエネルギー ( )を求めるステップ (c)と、

前記エネルギー (E₀)と前記エネルギーとの差 ( Δ E)を求めるステップ ( と、前記改変を受理するか否かを前記差 (Δ Ε)と受理確率を制御する温度パラメータ (T) とを用いた受理関数に基づいて判定するステップ ( と、

前記判定結果を前記現在の推定断面像に反映させて前記ステップ (a)に戻るステップ (わと、

前記ステップ (a)〜(わの繰り返しが所定回数に達するごとに前記温度パラメータ (T)の値を変更するステップ ( と、

前記ステップ (e)における判定結果が所定の停止条件を満たしているか否かを判定し、満たしていれば処理を終了するステップ (h)と、

を備えることを特徴とする画像再構成方法。

1 3. 対象物に放射線を照射して得られた投影像 (以下「放射線投影像」と呼ぶ)から前記対象物の断面像を求めるためのコンピュータプログラムであって、コンピュータに、

前記対象物の現在の推定断面像から演算により得られた投影像と前記放射線投影像との差を含む評価関数 (以下「エネルギー」と呼ぶ )(E₀)を求めるステップ (a)、前記現在の推定断面像の一部を改変するステップ (b)、

前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差を含むエネルギー (E を求めるステップ (c：)、

前記エネルギー (E。：)と前記エネルギー (E との差 (Δ Ε)を求めるステップ (d)、前記改変を受理するか否かを前記差 (Δ Ε)と受理確率を制御する温度パラメータ (T) とを用いた受理関数に基づいて判定するステップ (e)、

前記判定結果を前記現在の推定断面像に反映させて前記ステップ (a)に戻るステップ (わ、

前記ステップ (a)〜(f)の繰り返しが所定回数に達するごとに前記温度パラメータ (T)の値を変更するステップ (g)、

を実行させるための画像再構成プログラム。

14. 対象物に放射線を照射して投影像を得る手段 (A)と、

前記投影像から前記対象物の断面像を求める手段 (B)とを備えており、 .

前記手段 (B)は、

前記対象物の現在の推定断面像から演算により得られた投影像と前記対象物に放射線を照射して得られた投影像 (以下「放射線投影像」と呼ぶ)との差を含む評価関数 (以下「エネルギー Jと呼ぶ )(E₀)を求める手段 (b1 )と、

前記現在の推定断面像の一部を改変する手段 (b2)と、

前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差を含むエネルギー ( )を求める手段 (b3)と、

前記エネルギー (E₀)と前記エネルギー (E との差 (Δ Ε)を求める手段 (b4)と、前記改変を受理するか否かを前記差 (Δ Ε)と受理確率を制御する温度パラメータ (T) とを用いた受理関数に基づいて判定し、判定結果を前記現在の推定断面像に反映させる手段 (b5)と、

前記手段 (b1 )〜(b5)による一連の処理の繰り返しが所定回数に達するごとに前記温度パラメータ (T)の値を変更する手段 (b6)とを備える、ことを特徴とする CT装置。