TOP > 国内特許検索 > 画像再構成装置、画像再構成方法、画像再構成プログラム、CT装置 > 明細書

明細書 :画像再構成装置、画像再構成方法、画像再構成プログラム、CT装置

発行国 日本国特許庁(JP)
公報種別 特許公報(B2)
特許番号 特許第5190825号 (P5190825)
登録日 平成25年2月8日(2013.2.8)
発行日 平成25年4月24日(2013.4.24)
発明の名称または考案の名称 画像再構成装置、画像再構成方法、画像再構成プログラム、CT装置
国際特許分類 A61B   6/03        (2006.01)
FI A61B 6/03 350N
請求項の数または発明の数 14
全頁数 28
出願番号 特願2008-544217 (P2008-544217)
出願日 平成19年11月13日(2007.11.13)
国際出願番号 PCT/JP2007/072339
国際公開番号 WO2008/059982
国際公開日 平成20年5月22日(2008.5.22)
優先権出願番号 2006307058
優先日 平成18年11月13日(2006.11.13)
優先権主張国 日本国(JP)
審査請求日 平成22年11月15日(2010.11.15)
特許権者または実用新案権者 【識別番号】504255685
【氏名又は名称】国立大学法人京都工芸繊維大学
発明者または考案者 【氏名】西川 幸宏
個別代理人の代理人 【識別番号】110001427、【氏名又は名称】特許業務法人前田特許事務所
【識別番号】100077931、【弁理士】、【氏名又は名称】前田 弘
【識別番号】100110939、【弁理士】、【氏名又は名称】竹内 宏
【識別番号】100110940、【弁理士】、【氏名又は名称】嶋田 高久
【識別番号】100113262、【弁理士】、【氏名又は名称】竹内 祐二
【識別番号】100115059、【弁理士】、【氏名又は名称】今江 克実
【識別番号】100115691、【弁理士】、【氏名又は名称】藤田 篤史
【識別番号】100117581、【弁理士】、【氏名又は名称】二宮 克也
【識別番号】100117710、【弁理士】、【氏名又は名称】原田 智雄
【識別番号】100121728、【弁理士】、【氏名又は名称】井関 勝守
【識別番号】100124671、【弁理士】、【氏名又は名称】関 啓
【識別番号】100131060、【弁理士】、【氏名又は名称】杉浦 靖也
審査官 【審査官】小田倉 直人
参考文献・文献 特開2006-17568(JP,A)
特開2001-286463(JP,A)
特表2001-509400(JP,A)
調査した分野 A61B 6/03
特許請求の範囲 【請求項1】
対象物に放射線を照射して得られた投影像(以下「放射線投影像」と呼ぶ)から前記対象物の断面像を求める装置であって、
前記対象物の現在の推定断面像から演算により得られた投影像と前記放射線投影像との差を含む評価関数(以下「エネルギー」と呼ぶ)(E)を求める手段(a)と、
前記現在の推定断面像の一部を改変する手段(b)と、
前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差を含むエネルギー(E)を求める手段(c)と、
前記エネルギー(E)と前記エネルギー(E)との差(ΔE)を求める手段(d)と、
前記改変を受理するか否かを前記差(ΔE)と受理確率を制御する温度パラメータ(T)とを用いた受理関数に基づいて判定し、判定結果を前記現在の推定断面像に反映させる手段(e)と、
前記手段(a)~(e)による一連の処理の繰り返しが所定回数に達するごとに前記温度パラメータ(T)の値を変更する手段(f)と、
を備えることを特徴とする画像再構成装置。
【請求項2】
請求項1において、
前記手段(a)は、
前記現在の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記現在の推定断面像の局所領域の標準偏差との和を含むエネルギー(E)を求め、
前記手段(c)は、
前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記改変後の推定断面像の局所領域の標準偏差との和を含むエネルギー(E)を求める、
ことを特徴とする画像再構成装置。
【請求項3】
請求項1において、
前記手段(a)は、
前記現在の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記現在の推定断面像の局所領域のエントロピーとの和を含むエネルギー(E)を求め、
前記手段(c)は、
前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記改変後の推定断面像の局所領域のエントロピーとの和を含むエネルギー(E)を求める、
ことを特徴とする画像再構成装置。
【請求項4】
請求項1において、
前記手段(a)は、
前記現在の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記現在の推定断面像の局所領域の標準偏差と、前記現在の推定断面像の局所領域のエントロピーとの和を含むエネルギー(E)を求め、
前記手段(c)は、
前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記改変後の推定断面像の局所領域の標準偏差と、前記改変後の推定断面像の局所領域のエントロピーとの和を含むエネルギー(E)を求める、
ことを特徴とする画像再構成装置。
【請求項5】
請求項1において、
前記手段(a),(c),(d)に代えて、
[数1]によりΔHを算出し、算出したΔHを構成要素として含むΔE(ΔE=ΔH+…)を求める手段(h)を備える、
【数1】
JP0005190825B2_000019t.gif
ここで、前記対象物の現在の推定断面像をf(x,y)、前記手段(b)による改変部分をΔf(x,y)とすると、Δf(x,y)は座標(x,y)においてのみΔμの値を持ち、それ以外は0であるような断面像である。p(r,θ)は、前記対象物の現在の推定断面像から演算により得られた投影像である。p(r,θ)は、前記対象物の放射線投影像である。rは、投影像を撮影した1次元検出器のチャンネル位置である。θは、投影角度である。r(θ)=xcosθ+ysinθである。
ことを特徴とする画像再構成装置。
【請求項6】
請求項1において、
前記手段(a),(c),(d)に代えて、
[数2]によりΔHを算出し、算出したΔHを構成要素として含むΔE(ΔE=ΔH+…)を求める手段(h)を備える、
【数2】
JP0005190825B2_000020t.gif
ここで、前記対象物の現在の推定断面像をf(x,y)、前記手段(b)による改変部分をΔf(x,y)とすると、Δf(x,y)は座標(x,y)においてのみΔμの値を持ち、それ以外は0であるような断面像である。p(r,θ)は、前記対象物の現在の推定断面像から演算により得られた投影像である。p(r,θ)は、前記対象物の放射線投影像である。rは、投影像を撮影した1次元検出器のチャンネル位置である。θは、投影角度である。r(θ)=xcosθ+ysinθである。Mは、投影角度の数である。
ことを特徴とする画像再構成装置。
【請求項7】
請求項5または6において、
前記手段(h)は、
[数3]によりΔσを算出し、算出したΔσに係数cを掛けたもの(cΔσ)と前記ΔHとの和を構成要素として含むΔE(ΔE=ΔH+cΔσ+…)を求める、
【数3】
JP0005190825B2_000021t.gif
なお、σは座標(x,y)の周囲d×d画素の輝度値の標準偏差であり[数4]により表される。f,fはそれぞれ前記手段(b)による改変前後のf(x,y)の値である。
【数4】
JP0005190825B2_000022t.gif
ただし、
【数5】
JP0005190825B2_000023t.gif
【数6】
JP0005190825B2_000024t.gif
である。
ことを特徴とする画像再構成装置。
【請求項8】
請求項5または6において、
前記手段(h)は、
[数7]によりΔSを算出し、算出したΔSに温度パラメータ(T)を掛けたもの(-TΔS)と前記ΔHとの和を構成要素として含むΔE(ΔE=ΔH-TΔS+…)を求める、
【数7】
JP0005190825B2_000025t.gif
なお、Sは座標(x,y)の周囲d×d画素の局所領域画像のエントロピーであり[数8]により表される。
【数8】
JP0005190825B2_000026t.gif
また、
N:前記局所領域画像中の画素の総数、
:画素値がデジタル値でiである画素の総数、
:画素値がデジタル値でjである画素の総数、
k:定数、
であり、前記手段(b)による改変により画素値がデジタル値でiからjに変更されるものとする。
ことを特徴とする画像再構成装置。
【請求項9】
請求項1において、
前記手段(e),(f)に代えて、
前記改変を受理するか否かを前記差(ΔE)と受理確率を制御する温度パラメータ(T)とを用いた受理関数に基づいて判定し、判定結果を前記現在の推定断面像に反映させることを予約する手段(e1)と、
前記手段(a)~(d)および(e1)による一連の処理の繰り返しが所定回数に達するごとに、前記手段(e1)における予約を前記現在の推定断面像に反映させ、前記温度パラメータ(T)の値を変更する手段(f1)と、
を備えることを特徴とする画像再構成装置。
【請求項10】
対象物に放射線を照射して得られた投影像から前記対象物の断面像を求める装置であって、
前記対象物の放射線投影像p(r,θ)の逆投影像g(x,y)をフィルタなしの逆投影操作により計算する手段(m1)と、
前記対象物の現在の推定断面像f(x,y)から投影像p(r,θ)を計算し、さらに当該投影像p(r,θ)の逆投影像g(x,y)をフィルタなしの逆投影操作により計算する手段(m2)と、
前記対象物の現在の推定断面像f(x,y)に対する変更値を画素値とする画像Δμ(x,y)を生成する手段(m3)と、
各画素値について[数9]を適用して、画像ΔH(x,y)を生成する手段(m4)と、
【数9】
JP0005190825B2_000027t.gif
ここで、
M:投影角度の数、
である。
前記ΔH(x,y)を用いてΔE(x,y)を計算する手段(m5)と、
ここで、
ΔE(x,y)は、評価関数E(x,y)とE(x,y)との差である。
(x,y)は、前記推定断面像f(x,y)から演算により得られる投影像p(r,θ)と前記放射線投影像p(r,θ)との差を含む評価関数である。
(x,y)は、前記推定断面像f(x,y)に前記手段(m3)により得られた画像Δμ(x,y)を足し合わせたもの{f(x,y)+Δμ(x,y)}から演算により得られる投影像{p(r,θ)+Δp(r,θ)}と前記放射線投影像p(r,θ)との差を含む評価関数である。
前記ΔEが正である座標(x,y)について、前記Δμ(x,y)を0にする手段(m6)と、
前記推定断面像f(x,y)に前記手段(m6)により得られた画像Δμ(x,y)を足し合わせたものを新たな推定断面像f(x,y)として前記手段(m2)~(m6)による処理を繰り返す手段(m7)と、
を備えることを特徴とする画像再構成装置。
【請求項11】
対象物に放射線を照射して得られた投影像から前記対象物の断面像を求める装置であって、
前記対象物の放射線投影像p(r,θ)の逆投影像g(x,y)を[数10]により求める手段(m1)と、
【数10】
JP0005190825B2_000028t.gif
ここで、
r:投影像を撮影した1次元検出器のチャンネル位置、
θ:投影角度、
である。
前記対象物の現在の推定断面像f(x,y)から演算により投影像p(r,θ)を求め、さらに当該投影像p(r,θ)の逆投影像g(x,y)を[数11]により求める手段(m2)と、
【数11】
JP0005190825B2_000029t.gif
前記対象物の現在の推定断面像f(x,y)に対する変更値を画素値とする画像Δμ(x,y)を生成する手段(m3)と、
各画素値について[数12]を適用して、画像ΔH(x,y)を生成する手段(m4)と、
【数12】
JP0005190825B2_000030t.gif
ここで、
M:投影角度の数、
である。
各画素値について[数13]を適用して、画像Δσ(x,y)を生成する手段(m5)と、
【数13】
JP0005190825B2_000031t.gif
なお、σは座標(x,y)の周囲d×d画素の輝度値の標準偏差であり[数14]により表される。f,fはそれぞれ変更前後のf(x,y)の値である。
【数14】
JP0005190825B2_000032t.gif
ただし、
【数15】
JP0005190825B2_000033t.gif
【数16】
JP0005190825B2_000034t.gif
である。
各画素値について[数17]を適用して、画像ΔS(x,y)を生成する手段(m6)と、
【数17】
JP0005190825B2_000035t.gif
なお、Sは座標(x,y)の周囲d×d画素の局所領域画像のエントロピーであり[数18]により表される。
【数18】
JP0005190825B2_000036t.gif
また、
N:前記局所領域画像中の画素の総数、
:画素値がデジタル値でiである画素の総数、
:画素値がデジタル値でjである画素の総数、
k:定数、
であり、画素値がデジタル値でiからjに変更されるものとする。
[数19]に基づいてΔE(x,y)を計算する手段(m7)と、
【数19】
JP0005190825B2_000037t.gif
ここで、
c:係数
T:仮想温度(温度パラメータ)
である。
前記ΔEが正である座標(x,y)について、前記Δμ(x,y)を0にする手段(m8)と、
前記推定断面像f(x,y)に前記手段(m8)により得られた画像Δμ(x,y)を足したものを新たな推定断面像f(x,y)とする手段(m9)と、
前記Tをα倍(α<1)し、前記手段(m2)~(m9)による処理を繰り返す手段(m10)と、
を備えることを特徴とする画像再構成装置。
【請求項12】
対象物に放射線を照射して得られた投影像(以下「放射線投影像」と呼ぶ)から前記対象物の断面像を求める方法であって、
前記対象物の現在の推定断面像から演算により得られた投影像と前記放射線投影像との差を含む評価関数(以下「エネルギー」と呼ぶ)(E)を求めるステップ(a)と、
前記現在の推定断面像の一部を改変するステップ(b)と、
前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差を含むエネルギー(E)を求めるステップ(c)と、
前記エネルギー(E)と前記エネルギー(E)との差(ΔE)を求めるステップ(d)と、
前記改変を受理するか否かを前記差(ΔE)と受理確率を制御する温度パラメータ(T)とを用いた受理関数に基づいて判定するステップ(e)と、
前記判定結果を前記現在の推定断面像に反映させて前記ステップ(a)に戻るステップ(f)と、
前記ステップ(a)~(f)の繰り返しが所定回数に達するごとに前記温度パラメータ(T)の値を変更するステップ(g)と、
前記ステップ(e)における判定結果が所定の停止条件を満たしているか否かを判定し、満たしていれば処理を終了するステップ(h)と、
を備えることを特徴とする画像再構成方法。
【請求項13】
対象物に放射線を照射して得られた投影像(以下「放射線投影像」と呼ぶ)から前記対象物の断面像を求めるためのコンピュータプログラムであって、
コンピュータに、
前記対象物の現在の推定断面像から演算により得られた投影像と前記放射線投影像との差を含む評価関数(以下「エネルギー」と呼ぶ)(E)を求めるステップ(a)、
前記現在の推定断面像の一部を改変するステップ(b)、
前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差を含むエネルギー(E)を求めるステップ(c)、
前記エネルギー(E)と前記エネルギー(E)との差(ΔE)を求めるステップ(d)、
前記改変を受理するか否かを前記差(ΔE)と受理確率を制御する温度パラメータ(T)とを用いた受理関数に基づいて判定するステップ(e)、
前記判定結果を前記現在の推定断面像に反映させて前記ステップ(a)に戻るステップ(f)、
前記ステップ(a)~(f)の繰り返しが所定回数に達するごとに前記温度パラメータ(T)の値を変更するステップ(g)、
を実行させるための画像再構成プログラム。
【請求項14】
対象物に放射線を照射して投影像を得る手段(A)と、
前記投影像から前記対象物の断面像を求める手段(B)とを備えており、
前記手段(B)は、
前記対象物の現在の推定断面像から演算により得られた投影像と前記対象物に放射線を照射して得られた投影像(以下「放射線投影像」と呼ぶ)との差を含む評価関数(以下「エネルギー」と呼ぶ)(E)を求める手段(b1)と、
前記現在の推定断面像の一部を改変する手段(b2)と、
前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差を含むエネルギー(E)を求める手段(b3)と、
前記エネルギー(E)と前記エネルギー(E)との差(ΔE)を求める手段(b4)と、
前記改変を受理するか否かを前記差(ΔE)と受理確率を制御する温度パラメータ(T)とを用いた受理関数に基づいて判定し、判定結果を前記現在の推定断面像に反映させる手段(b5)と、
前記手段(b1)~(b5)による一連の処理の繰り返しが所定回数に達するごとに前記温度パラメータ(T)の値を変更する手段(b6)とを備える、
ことを特徴とするCT装置。
発明の詳細な説明 【技術分野】
【0001】
本発明は、対象物の放射線投影像から断面像を再構成する技術に関する。
【背景技術】
【0002】
計算機トモグラフィー(CT)は、X線などの光線を使って、様々な角度からの物体の投影像を撮影し、その後、計算によって断面像を得る技術である。図1に典型的なX線CT装置の模式図を示す。
【0003】
典型的なCT装置では、X線光源が回転移動することにより、測定対象の様々な角度からの投影像を取得する仕組みになっている。こうして得られた投影像から計算操作によって断面像を得るのがCTである。通常は、Filtered Back-projection法(FBP)によって投影像から断面像を変換操作によって求める。FBPでは、投影像にフィルタ(基本的には微分フィルタ)を施した後、もともとの投影方向に「逆投影」という操作を施すことで断面像を得る。このとき、微分フィルタはノイズや誤差を増幅し、アーティファクト(本来存在しない誤差や虚像)を生みやすいという特徴がある。さらに、逆投影操作はそのアーティファクトを断面像全体に伝播させる効果がある。そのため、CTにおけるアーティファクトは問題のある部分だけに留まらず、断面像全体を損なうという点で致命的になることが多い。
【0004】
アーティファクトの多くは、FBPに含まれるフィルタ操作や逆投影操作が原因なのだから、FBPを用いなければ、アーティファクトのない断面像が得られるはずである。FBP以外で断面像を計算する方法として、Algebraic Reconstruction Technique(ART)が歴史的に重要である。ARTはFBPが考案される前から利用されるほど歴史のある方法である。ARTでは断面像の計算をフィッティング問題と捉え、断面像をパラメータ、投影像をフィッティング対象として、断面像から計算した投影像(p)が、実験で求めた投影像(p)に一致するように逐次的に断面像を修正する。ARTでは(p-p)が0になるように断面像を漸近的に改良してゆく点が特徴となっている。ARTはFBPを全く用いないかわりに、断面像を求める計算に時間がかかるため、今では特殊な用途(地震波の解析など)にしか用いられない。ARTはFBPほど極端なアーティファクトは現れないが、得られる断面像はFBPの方が自然であることが多い。
【0005】
さらに、フィルタ操作や逆操作に由来しないアーティファクトの要因として、投影像におけるデータの欠落・過少も挙げられる。データが少なければ結果として得られる断面像の画質も低下するのは当たり前である。FBPではデータの欠落・過少も致命的なアーティファクトを生み出すことが知られており、大きな問題と認識されている。フィッティングに基づくCT(ARTを含む)はデータの欠落や過少に対してFBPよりは強いとされている。しかしながら、CTにおけるデータの欠落は、極端に「条件が悪い」問題とされており、ARTなどを用いても改善が得られにくい。ARTでは(p-p)を誤差として考えるのであるが、条件が悪いときはフィッティングに失敗しやすい。いわゆる二乗誤差(p-p)の最小化を行った方がフィッティングとしては安定する。二乗誤差を最小化する方法としては、最小二乗法が最も標準的である。最小二乗法では、一辺がパラメータの個数だけある正方行列の逆行列を求める。CTにおけるパラメータは断面像を構成するそれぞれの画素値であるため、パラメータの個数は莫大になる。1000×1000画素の断面像だと、パラメータの数は100万個になり、行列の要素数は1兆個と莫大である。従って、正攻法の最小二乗法では行列が大きくなりすぎで破綻してしまう。そこで、正攻法の最小二乗法を避け、Simultaneous iterative reconstruction technique (SIRT)やIterative least squaretechnique (ILST)が考案された。ARTと同様に断面像の計算をフィッティング問題として捉えるのであるが、先に述べたように、正攻法の最小二乗法を避けなければならないため、SIRTもILSTも計算の途中に投影操作の逆演算としてFBPが利用される。そのため、FBPに起因するフィルタ操作と逆投影にまつわる問題を根本的には解決しない。おそらくそれが理由となり、SIRTやILSTでは各種アーティファクトの「低減」効果が報告されているに留まる。

【非特許文献1】Yazdi M, 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.
【発明の開示】
【発明が解決しようとする課題】
【0006】
以上の議論から、アーティファクトの無い断面像を得るには、FBPを用いない、データの欠落に対して強い、という性質が重要であることがわかる。そこで、ARTと同様にFBPを用いず、フィッティングの評価関数として二乗誤差を利用しようというのは、誰でも思いつくことである。にもかかわらず、CTにおけるアーティファクトの問題が根本的に解決していないのは、この単純な考えの実現は簡単ではないからである。第一の問題は、断面像の計算はフィッティングの問題としては大規模すぎて、計算に時間がかかることである。適切な高速化が必須である。第二の問題は、最小二乗法はフィッティングアルゴリズムとしては「条件が悪い」問題に対して弱いことである。そのため、最小二乗法に基づく既存のアルゴリズム(ILSTなど)の延長では、アーティファクトの問題を解決できないと思われる。第三の問題は、SIRTやILSTがそうであるように、既存のアルゴリズムではFBPを完全に排除することができていないことである。
【課題を解決するための手段】
【0007】
本発明の第一の新規性は、フィッティング方法としてSimulated Annealing法(SA)を採用したことである。SAは公知の技術であるが、フィッティングに時間がかかることが知られており、CTと本質的に相性が悪い。しかしながら、敢えてSAをCTに適用することで、FBPを全く用いずとも、二乗誤差の最小化による断面像の計算ができるようになった。また、SAはデータの欠落など条件の悪いフィッティングに対してかなり安定であるという性質があり、この点においてもアーティファクトの根本的な解決に有利である。以上の事柄を鑑み、SAをCTに適用することで、アーティファクトの根本的解決に効果を見出した点が、本発明の重要な新規性である。単純にSAのCTへの適用を考えたとき、断面像から投影像を求める計算を無数に繰り返すことになる。まともに計算すると、断面像から投影像を求める計算は、FBPとほぼ同じであり、SAをCTに適用すると、FBPに比べて数百万倍の時間がかかってしまうことになる。現在の高速な計算機を以ってしても、「年」単位の計算時間を要してしまう。そこで、本発明では、式変形によって計算量を大幅に低減することで、この問題を解決した。この工夫は、SAをCTに適用する際に必要な技術であり、本発明の重要な要素である。
【0008】
第二の新規性は、フィッティングの評価関数として、二乗誤差だけでなく、アーティファクトを積極的に破壊するスムージング項とエントロピー項を導入した点である。従来のCTでは、pとpとの差あるいは二乗誤差だけを対象としていた。この場合、アーティファクトを御座なりにした方が良いフィッティング結果(小さな誤差)が達成できるという可能性が常にあり、現実の多くのケースではそのような条件になっていた。すなわち、アーティファクトを別のアーティファクトで相殺する状態になっていたのである。そのため、アーティファクトはなかなか消えないのである。本発明においても二乗誤差だけを評価関数とした場合には、アーティファクトは完全には消えなかった。この事実は、アーティファクトのない断面像を得るには、二乗誤差以外の要請が望ましいことを意味している。本発明では、統計力学の基本にのっとり、エントロピー項やスムージング項を考案した。これらの項は、断面像は滑らかで自然な濃淡画像であるべきだ、という当たり前の要請を数式の形で表現したものである。エントロピー項は、アーティファクトを破壊し、断面像全体の画質を均質にするように、フィッティングを誘導する。スムージング項は、エントロピー項で誘導される断面像のザラツキを押さえる効果がある。これらの項を導入することで、アーティファクトが消え、かつ自然な断面像が得られるようになった。エントロピー項とスムージング項は、それぞれ単独でもアーティファクトを低減する効果があるが、2つを組み合わせるとさらに効果的である。
【0009】
なお、本発明においては、一般の定義よりも広い概念として、X線・可視光・電波を含む電磁波、電子や荷電粒子からなる粒子線、媒体の振動である音波等を全て含めて「放射線」と呼ぶことにする。
【発明の効果】
【0010】
本発明の第一の効果として、データの欠落が要因になっているアーティファクトの低減効果が挙げられる。データの欠落が問題になる場合として、観察対象に不透明部位がある場合、投影角度に制限がある場合、投影角度の刻みが不均一である場合、コーンビームやヘリカルスキャンなどの三次元CTの場合が挙げられる。
【0011】
特に効果があることがわかっているのは、観察対象に不透明部位がある場合に現れるメタルアーティファクトの低減効果である。メタルアーティファクトとは、観察対象中にX線に対して不透明な部位(多くの場合、金属部分)が存在する場合、CTで得られる断面像全体(不透明部位だけでない)が破壊的に乱されることを指す。メタルアーティファクトの原因は、不透明部分で投影像の輝度が不連続に変化すること、不透明部分の情報が欠落することである。FBPにおいて微分フィルタが適用されると、輝度の不連続な変化は異常値となり、逆投影操作によって、不透明部位を中心とした筋状のアーティファクトを生み出す。さらに、情報の欠落によって不透明部位とは直接関係のない部分に、予期しない明暗が生じる。FBPを用いず、情報の欠落に対して安定なSAを利用した本発明がメタルアーティファクトの除去に有効であろうことは想像に難くない。
【0012】
実用上、重要なのは、コーンビームやヘリカルスキャンへの応用であろう。両者は3次元CTと呼ばれ、現在急速に普及しつつある技術である。しかしながら、3次元CTには特有のアーティファクトが現れることが知られており、その原因はわかっているのだが解決には至っていない。その原因は、コーンビームの場合には、データの欠落である。コーンビームでは完全な断面像を得るための条件を満たすことができず、データの欠落に由来するアーティファクトが現れる。ヘリカルスキャンでのアーティファクトの根本的な原因は、逆投影操作にある。ヘリカルスキャンでは系の幾何学的異方性(らせん)がフィルタ操作・逆投影操作に影響し、風車アーティファクトが現れる。本発明は、フィルタ操作・逆投影操作を必要とせず、データの欠落にも強いので、3次元CTにおける諸問題を解決できる。
【0013】
投影角度が不均一である場合として挙げられるのは、地震波による地球内部のCTや、人工衛星からの電波を用いた大気状態のCTによる解析である。これらはFBPが利用できない典型として知られている。本発明を利用することで、解析精度の向上を期待できる。
【0014】
本発明の第二の効果は、投影像撮影の迅速化、X線被爆量の低減である。SAはデータの欠落だけでなく、過少にも安定であるという特徴があり、本発明も同様である。CTの場合、データの過少に相当するのは、投影角度の数が少ない場合である。すなわち、本発明を利用すると、従来のCTより投影角度の数を少なくすることができる。投影角度の数とは放射線による投影像の撮影枚数である。撮影枚数は、撮影時間と被爆量に比例するので、撮影枚数が少なければ、撮影時間の短縮、被爆量の低減につながる。
【0015】
また、データの過少の別のパターンとして、投影像の画質が不良(S/N比が低い)な場合が挙げられる。本発明はそのような場合にも比較的安定であることがわかっている。投影像の画質の低下を許容できるとなると、これも撮影時間の短縮・被爆量の低減につながる。あるいは、S/N比が極端に低いSPECT,PETにおいて本発明は画質の向上に寄与できるだろう。
【0016】
本発明の第三の効果は、断面像の輝度値の決定精度が高いことである。本発明では測定した投影像をかなり忠実に再現する断面像が得られる。再現の精度はFBPより2桁程度高い。これは二乗誤差を最小化するというフィッティングアルゴリズムの恩恵である。輝度値の決定精度が高いことで、輝度値の定量性が保証され、輝度値を用いた密度測定が可能になる。これは骨密度測定などにおける測定精度の向上につながる。また、異常のある部位(腫瘍のある臓器など)の検出精度の向上にも寄与するだろう。
【0017】
本発明の第四の効果は、本発明で得られる断面像のコントラストがFBPよりも高いことである。第三の効果で述べたように、本発明は輝度値の決定精度が高いのであるが、その副次的な効果として、断面像のコントラストが高くなるのである。コントラストが高いと見掛けの空間分解能も高くなる傾向がある。その結果、本発明で得られる断面像は、従来のものより高画質になる。この効果は、特殊な条件・用途のCTではなく、普通のCTにも本発明が役立つという点で、特筆すべきものである。
【発明を実施するための最良の形態】
【0018】
(第1の実施形態)
本発明の実施形態によるX線CT装置の概略構成を図1に示す。このX線CT装置は、撮影部とコンピュータとを備えている。撮影部は、X線を測定対象に照射するための光源と、測定対象を透過したX線を検出する検出器とを備えており、測定対象の周囲を多数の方向からX線により投影して投影像を取得する。コンピュータは、X線CT装置全体を制御する制御部と、撮影部により得られたX線投影像に基づいて測定対象の関心領域の断面像を生成する画像再構成処理部とを備えている。なお、図1に示した構成はこれ以降の各実施形態に共通のものであり、各実施形態では画像再構成処理部において行われる処理がそれぞれ異なっている。
【0019】
本実施形態の画像再構成処理部は、投影像から断面像を求めるためのフィッティング方法としてSimulated Annealing法(SA)を採用している。そこで、まず、Simulated Annealing法(SA)のフレームワークを説明する。SAはMonte Carlo法から派生したフィッティングアルゴリズムで、乱数に基づいてフィッティングを進める点と、熱力学を模して仮想エネルギーや仮想温度を取り扱う点、が特徴となっている。SA自体は公知の技術であり以下のステップ(i)~(vi)により行われる。
(i)パラメータを乱数に基づいて1つ選択し、さらに乱数に基づいてそのパラメータを変更する(乱数)。
(ii)変更後の評価を行う。評価関数として仮想的なエネルギーEを考える。通常のSAではEは二乗誤差の総和である。変更前後のEの変化をΔEとする(評価)。
(iii)評価が改善される(ΔE<0)なら、その変更を受け入れる(変更)。
(iv)exp(-ΔE/T)の確率で改悪を受け入れる。
(v)Tを少しずつ減じる。
(vi)(i)から繰り返す。
【0020】
SAでは、(iv)にあるように、ボルツマン統計に従って改悪を受け入れることで、局所解から脱出できる可能性を担保している。そのため、局所解にとらわれることなく、大域解に辿りつくことができ、条件の悪いフィッティングに対しても安定に機能する。また、Tを徐々に減じることで、大域解にソフトランディングするようになっている。(i)から(v)まで一通り実行することを1モンテカルロステップと呼ぶ。SAではモンテカルロステップを無数に繰り返すことで、フィッティングが進行する。計算に必要な時間は、1モンテカルロステップに必要な時間×必要なモンテカルロステップ数である。必要なモンテカルロステップ数は、パラメータの数・自由度に比例する。
【0021】
次に、請求項に沿った形で本実施形態の説明を行う。CTをフィッティング問題と考えた場合、フィッティングパラメータは断面像(f(x,y))である。フィッティング対象となるデータは、測定した投影像(p(r,θ))である。rは投影像を撮影した1次元検出器のチャンネル位置、θは投影角度である。座標の定義を図2に示す。p(r,θ)は、角度θを変えながら投影像を測定することで得られるデータセットであり、横軸にr、縦軸にθを取ると、2次元の画像とみなすことができる。このようなデータをsinogramと呼ぶ。典型的なsinogramを図3に示す。
【0022】
極論すれば、CTはsinogramから断面像を求める技術である。本実施形態では、暫定の断面像f(x,y)と測定した投影像p(r,θ)との二乗誤差を考えるが、二乗誤差を計算するために、数101を使って、f(x,y)から計算で投影像(p(r,θ))を求める。
【数101】
JP0005190825B2_000002t.gif

【0023】
数101は図2において、sの方向に和をとっており、f(x,y)のs方向の投影像の計算であることがわかる。こうして求めたp(r,θ)を用いて、二乗誤差Hは以下のように計算できる。
【数102】
JP0005190825B2_000003t.gif

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

【0025】
Tは仮想温度(温度パラメータ)、Sはエントロピー、σは画素値の標準偏差、cはスムージング項の強さを表す係数である。TSがエントロピー項、cσがスムージング項である。Sやσの定義・計算法は後述する。本実施形態では、これらの数式を使い、図4に示すような手順で断面像を計算する。
・ステップ(a):何らかの形で求めた暫定の断面像f(x,y)に対し、評価関数(以下「エネルギー」と呼ぶ)(E)を求める(ST10およびST20)。
・ステップ(b):乱数で(x,y)およびΔμを選択し、断面像f(x,y)の一部を改変する(ST30)。
・ステップ(c):改変後の断面像f(x,y)に対して、評価関数Eを求める(ST40およびST50)。
・ステップ(d):前記エネルギー(E)と前記エネルギー(E)との差(ΔE)を求める(ST60)。
・ステップ(e):前記改変を採用するか否かを前記差(ΔE)および温度パラメータ(T)を用いた受理関数(典型的にはボルツマン統計)により判定する(ST70~ST110)。
・ステップ(f):前記ステップ(a)に戻る(ST120)。
・ステップ(g):前記ステップ(a)~(f)の繰り返しが所定回数に達するごとに仮想温度(T)の値を変更する(ST130)。
・ステップ(h):前記ステップ(e)における判定結果が所定の停止条件を満たしているか否かを判定し、満たしていれば処理を終了する。ここでは、ST80,ST100が実行された場合を「成功」、ST110が実行された場合を「失敗」とし、成功の確率が10%(この値は適宜調整可能)を下回ると終了する。そして、終了時点において採用されている推定断面像を測定対象の断面像とし、これをコンピュータのディスプレイに表示したり記録媒体に記録したりする。
【0026】
ステップ(a)からステップ(h)の一連の操作が請求項12に対応し、本実施形態ではこれら一連の操作を図1の画像再構成処理部が行う。
【0027】
SAの基本手順(i)~(vi)との対応を考えると、ステップ(b)が(i)、ステップ(a)(c)(d)が(ii)、ステップ(e)に(iii)と(iv)が含まれ、ステップ(g)が(v)、ステップ(f)が(vi)にそれぞれ対応する。このことから、本実施形態はCTにSAを忠実に当てはめたものであることがわかる。
【0028】
なお、上記ステップ(h)における処理終了の判定ついては、図1の画像再構成処理部において行うのではなく、コンピュータのディスプレイに推定断面像を逐次表示し、これをユーザが目で見て処理の終了をコンピュータに指示できるようにしてもよい。
【0029】
(第2の実施形態)
数103に示されたEを計算するためには、数101と数102を経ることになり、s,r,θに関して級数和を実行することになる。3重の積分(級数和)を計算するので、大変時間がかかる。すなわち、1モンテカルロステップに必要な計算時間が長いことを意味する。さらに、CTではパラメータの数が莫大であり、その結果、SAの実施に必要な総計算時間は、現在の計算機でも「年単位」になってしまう。
【0030】
そこで、本実施形態では、Eを計算するのではなく、変更を行ったときの変化分ΔEだけを主に計算する。
【数104】
JP0005190825B2_000005t.gif

【0031】
今、暫定断面像f(x,y)に関する変更をΔf(x,y)とする。Δf(x,y)は座標(x,y)においてのみΔμの値を持ち、それ以外は0であるような断面像である。そのΔf(x,y)の投影像Δp(r,θ)は数101と同じ方法で計算できる。このΔp(r,θ)を用いて、ΔHは以下のように計算できる。
【数105】
JP0005190825B2_000006t.gif

【0032】
これを書き下して、
【数106】
JP0005190825B2_000007t.gif
となる。今、Δf(x,y)は(x,y)においてのみ値を持つので、Δp(r,θ)はr(θ)=xcosθ+ysinθにおいてのみ値Δμを持ち、それ以外は0である。従って、数106の{}の中身も、r(θ)=xcosθ+ysinθにおいてのみ値を持つ。そのため、級数和はrとθの両方に対して行う必要はなく、以下のように表すことができる。
【数107】
JP0005190825B2_000008t.gif

【0033】
数107ではθに関する級数和だけになっている点が重要である。これにより、計算量を大幅に削減することができる。pやpはデジタル画像であるため、数107の計算にあたっては、通常は、r(θ)に関する内挿を行う必要がある。そのため、数107の{ }内をさらに展開することはできない。しかしながら、誤差を容認して数107を展開すると次式が得られる。
【数108】
JP0005190825B2_000009t.gif

【0034】
ただし、Mは投影角度の数である。数107の代わりに数108を用いるとΔHの値に1%程度の誤差が生じる。しかしながら、数108は数107より高速であり、本発明での利用価値は高い。
【0035】
次にΔσの計算を説明する。σは座標(x,y)付近の輝度値の標準偏差である。(x,y)付近として、(x,y)の周囲d×d画素を考える(実施例ではd=5を用いている)。それらの画素の標準偏差は次式で求めることができる。
【数109】
JP0005190825B2_000010t.gif
ただし、
【数110】
JP0005190825B2_000011t.gif
【数111】
JP0005190825B2_000012t.gif
である。数109~111を用いると、Δσは次のように計算することができる。
【数112】
JP0005190825B2_000013t.gif
ただし、f,fは変更前および後のf(x,y)の値である。
【0036】
次にΔSの計算であるが、その前にエントロピーSの定義を行う。一般に計算機が取り扱う画像はデジタル画像であり、画素の座標(x,y)だけでなく、画素の値もデジタルの値となっている。そこで、1つの画素を1つの量子と考え、画素値を量子状態とみなすことにする。すると画像は複数の量子からなる系であると考えることができる。統計力学に則ると、そのような系に対するエントロピーは次のように定義される。
【数113】
JP0005190825B2_000014t.gif

【0037】
Nは画像中の画素の総数であり、Nは画素値がデジタル値でiである画素の総数である。kは通常の物理ではボルツマン定数であるが、本発明は物理とは直接関係がないため任意である。本発明では、σと同様にSも(x,y)の周囲d×d画素で定義されるとした。
【0038】
数113で定義されたSに対してΔSを考える。変更によって、画素値はデジタル値でiからjに変更されるとする。すると、ΔSは次式のようになる。
【数114】
JP0005190825B2_000015t.gif
すなわち、変更によって、Nが1減って、Nが1増えるのである。数114を展開すると、非常に簡単な式が得られる。
【数115】
JP0005190825B2_000016t.gif

【0039】
以上をまとめて、本実施形態の手順を書き下すと次のようになる。
(I)乱数を使って(x,y)およびΔμを選択する。
(II)数104、108(または107)、112、115を使ってΔEを計算する。
(III)評価が改善される(ΔE<0)なら、f(x,y)にΔμを足す。
(IV)ΔE>0の時にも、exp(-ΔE/T)の確率で、f(x,y)にΔμを足す。
(V)Tを少しずつ減じる。
(VI)(I)から繰り返す。
(VII)所定の停止条件を満たしているか否かを判定し、満たしていれば処理を終了する。ここでは、(III),(IV)においてΔμが足された場合を「成功」、それ以外を「失敗」とし、成功の確率が所定値(例えば10%、この値は適宜調整可能)を下回ると終了する。そして、終了時点において採用されている推定断面像を測定対象の断面像とし、これをコンピュータのディスプレイに表示したり記録媒体に記録したりする。
【0040】
手順(I)~(VI)はSAの基本手順(i)~(vi)を忠実に踏襲していることがわかる。なお、手順(I)~(VII)の処理は図1の画像再構成処理部において行われる。
【0041】
なお、上記手順(VII)における処理終了の判定ついては、図1の画像再構成処理部において行うのではなく、コンピュータのディスプレイに推定断面像を逐次表示し、これをユーザが目で見て処理の終了をコンピュータに指示できるようにしてもよい。
【0042】
(第3の実施形態)
次に、上で述べた方法の重要な変形について述べる。上で述べた方法では数107あるいは数108を利用し、θに関する級数和を求めている。そのため、画素1つの変更を検討するのにM回の繰り返し演算を要する。本実施形態はこの部分の計算量をさらに減らすものである。
【0043】
まず、p(r,θ)の逆投影像g(x,y)を考える。
【数116】
JP0005190825B2_000017t.gif
数116では、フィルタ操作をしていないため、g(x,y)は断面像f(x,y)に似ているが、ぼやけたような画像になる。g(x,y)を用いると、数108は次のようになる。
【数117】
JP0005190825B2_000018t.gif
ただし、g(x,y)はp(r,θ)の逆投影像で、数116と同様の方法で計算される。数117では繰り返し演算がないという点で数108より優れている。g(x,y)は計算過程で全く変化しないため、予め計算しておくことができる。一方、g(x,y)はf(x,y)の1点を変更する度にわずかずつ変化するため、厳密には(I)~(IV)の手順を1回経るごとに再計算が必要である。しかしながら、f(x,y)の変更に伴うg(x,y)の変化は少ないという仮定をすると、別の手順を適用することができ、それが本実施形態である。以下に手順を示す。なお、以下の手順(1)~(9)の処理は図1の画像再構成処理部において行われる。
(1)p(r,θ)からg(x,y)を求める。
(2)f(x,y)からp(r,θ)を求め、さらにp(r,θ)からg(x,y)を求める。
(3)f(x,y)に対する変更値を画素値とする画像Δμ(x,y)を乱数で生成する。
(4)各画素値について数117を適用し、ΔH(x,y)という画像を計算する。
(5)同様に各画素値について数112、数115を適用し、Δσ(x,y)、ΔS(x,y)という画像を計算する。
(6)数104に基づいてΔE(x,y)を計算する。
(7)ΔEが正である座標(x,y)について、Δμ(x,y)を0にする。
(8)f(x,y)にΔμ(x,y)を足す。
(9)Tをα倍(α<1)し、(2)から繰り返す。
(10)所定の停止条件を満たしているか否かを判定し、満たしていれば処理を終了する。ここでは、(7)においてΔμ(x,y)が0にされた場合を「失敗」、されなかった場合を「成功」とし、成功の確率が所定値(例えば10%、この値は適宜調整可能)を下回ると終了する。そして、終了時点において採用されている推定断面像を測定対象の断面像とし、これをコンピュータのディスプレイに表示したり記録媒体に記録したりする。
【0044】
なお、上記手順(10)における処理終了の判定ついては、図1の画像再構成処理部において行うのではなく、コンピュータのディスプレイに推定断面像を逐次表示し、これをユーザが目で見て処理の終了をコンピュータに指示できるようにしてもよい。
【0045】
手順(I)~(VII)と、手順(1)~(10)を図5に模式的に示した。図5ではパラメータの値(本件では画素の値)を軸とし、仮想エネルギーEを濃淡で表している(Eが小さいほど濃い)。図示の都合により、パラメータの数(画素の数)は2個の場合を取り扱っている。図5(a)に示す手順(I)~(VII)では、ジグザグに折れ曲がりながら、時にでたらめな方向に進みながら、最終的にはEの最小値を探索してゆく。一方、図5(b)に示した手順(1)~(10)では、Eの勾配を考慮しながら、Eの最小値を探索してゆく。図5(b)の方が効率が良いということが直感的に理解できるだろう。数117による計算量の低減及び、図5(b)に示した探索効率の向上により、本実施形態は計算速度的に有利である。手順(1)~(10)はSAの基本アルゴリズム(i)~(vi)を忠実に踏襲するものではないが、(I)~(VII)の手順の自然な発展形である。
【0046】
(その他の実施形態)
上述の各実施形態で説明した処理を実行する画像再構成処理部は、これらの処理をコンピュータに実行させるためのプログラム、当該プログラムがインストールされたコンピュータ、さらには、これらの処理を実行する専用LSIなどを使って実現可能である。
【0047】
(実施例)
実施例としては、シミュレーションを用いたメタルアーティファクト除去効果を取り扱う。まず、比較のために、[非特許文献1]に掲載されている手法での結果を図6(非特許文献のfig.5)に示す。
【0048】
図6において、左のカラム(a)(d)(g)はMetal Artifactがないオリジナルのファントム画像である。中央(b)(e)(h)は、オリジナルのファントムの所定の位置に不可視領域を設定し、通常のFBP法により再構成した画像である。右(c)(f)(i)は、上記と同様に不可視領域を設定し[非特許文献1]の手法により再構成した画像であり、従来技術で最もうまくMetal Artifactが除去できている例の1つである。
【0049】
そこで、[非特許文献1]から図6(a)の画像を取り込み、文献と同じ位置に不可視領域を設定し、シミュレーションで投影像を計算した。結果を図7に示す。図7が本実施例におけるp(r,θ)に対応する。図7について、本発明によって再構成した結果および、比較のための図6(c)の拡大図を図8に示した。効果は一目瞭然である。本実施形態(b)ではMetal Artifactの特徴を見つけ出すことすら困難である。
【0050】
ただし本実施形態(b)ではエッジ付近が若干ボケている。これは「断面像を滑らかにする因子」(スムージング項cσ)がやや強いためである。Metal Artifactを低減するには、その因子をやや強めに設定する必要があることがわかっている。このように本実施形態によれば、仮想エネルギーEを計算する際の各因子のバランスには改良の余地があるものの、Metal Artifact除去の効果は非常に高い。
【0051】
なお、参考までに、スムージング項cσ、エントロピー項TSを設定せずに本実施形態のアルゴリズムを実行した場合の結果を図9に示す。図9(a)はスムージング項cσなし(エントロピー項TSのみ)の場合の結果、図9(b)はエントロピー項TSなし(スムージング項cσのみ)の結果、図9(c)はエントロピー項TS、スムージング項cσともになしの場合の結果を示す。図9と図8(b)とを比較すると、本実施形態においては、エネルギーEの因子としてスムージング項cσとエントロピー項TSの両方を設定したほうが、いずれか一方のみを設定した場合、ともに設定しない場合よりもMetal Artifact除去の効果が高いことが分かる。スムージング項だけでは、不透明を中心に放射状に広がる縞模様がアーティファクトして残るという特徴がある。一方、エントロピー項のみでは、断面像がざらつきS/Nが低いという特徴がある。また、最終的な画質は、徐冷過程におけるスムージング項の係数cと仮想温度Tの比が重要であるようである。
【0052】
次に、回転角度制限がある場合の実施例を図10に示す。図10(a)は、FBPのみ(通常のCT)による再構成結果を示す。図10(b)は、ILST(最小二乗法に基づく画像復元法)を適用した結果を示す。図10(c)は、本実施形態による再構成結果を示す。これらを比較すると、角度制限時に現れるアーティファクトに対して本実施形態が有効であることが分かる。角度制限によるアーティファクトには、円状の領域が両端の尖ったアーモンド状の領域になり、アーモンド状領域の尖っていない胴回りの輝度が反転するという特徴がある。本実施形態によると、それら両方の特徴が低減される。
【産業上の利用可能性】
【0053】
本発明は、Metal Artifactに対して顕著な効果があるので、Metal Artifactが深刻な分野、たとえば歯のCT、金属製インプラントを含むCTに特に有用である。
【0054】
また、本発明は情報欠落のある投影像のセットに対して一般的に有効である。情報欠落が顕著な場合の一つに、投影角度制限がある場合を挙げることができる。投影角度制限が問題になる分野は、三次元電子顕微鏡法、マンモグラフィーのCT版、平行移動CT法(特開2006-71472号公報)などである。情報欠落が問題になる別の場合として、コーンビームCT、ヘリカルスキャンCTなど3次元CTが挙げられる。本発明は、3次元CTに現れるアーティファクトを除去する目的にも有効である。
【0055】
さらには、情報が極端に少ない系に関しても利用可能である。例えば、蛍光X線CT、地震波による地球内部のCTなどに有用である。
【0056】
また、本発明の副次的な効果として、断面像における輝度値(X線では線吸収係数に対応する)を従来法よりも高い精度で決定できることが挙げられる。この効果は骨密度測定等における精度改善に応用可能である。
【0057】
本発明を用いると、従来の方法よりコントラストの高い再構成像を得ることができる。そのため、アーティファクト等が問題とならない通常のX線CTにおいても、本発明の利用価値は高い。また、本発明はデータの過少に対しても安定であることから、投影像測定の時間短縮、ひいてはX線被爆量の低減に有効である。これらをあわせて考えると、本発明は、既存のすべてのX線CT技術を置き換える可能性を秘めている。
【図面の簡単な説明】
【0058】
【図1】X線CT装置の模式図である。
【図2】断面像f(x,y)、投影像p(r,θ)の関係を示す模式図である。
【図3】典型的なp(r,θ)の例(横軸に検出器のチャンネル位置、縦軸に投影角度をとった画像)である。
【図4】本発明の基本手順を示すフローチャートである。
【図5】(a)は、手順(I)~(VI)のフィッティング経過の模式図である。(b)は、手順(1)~(9)のフィッティング経過の模式図である。
【図6】[非特許文献1]に掲載されている手法での結果を示す図である。
【図7】図6(a)からシミュレーションによって求めたsinogramである(白い部分は不可視領域)。
【図8】(a)は、図6(c)をさらに拡大した図である。(b)は、本発明での結果を示す図である。
【図9】仮想エネルギーEの各項の効果を表した模式図である。
【図10】角度制限がある場合のアーティファクトと本発明でのアーティファクト低減効果を示す図である。
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8
【図10】
9