(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2022-10-07
(45)【発行日】2022-10-18
(54)【発明の名称】インテリアCT画像生成方法
(51)【国際特許分類】
A61B 6/03 20060101AFI20221011BHJP
G06T 1/00 20060101ALI20221011BHJP
【FI】
A61B6/03 350N
G06T1/00 290B
A61B6/03 ZDM
(21)【出願番号】P 2019508707
(86)(22)【出願日】2018-02-09
(86)【国際出願番号】 JP2018004603
(87)【国際公開番号】W WO2018179905
(87)【国際公開日】2018-10-04
【審査請求日】2020-12-24
(31)【優先権主張番号】P 2017061390
(32)【優先日】2017-03-27
(33)【優先権主張国・地域又は機関】JP
【国等の委託研究の成果に係る記載事項】(出願人による申告)平成28年度、国立研究開発法人科学技術振興機構、戦略的創造研究推進事業「位相画像解析研究」委託研究、産業技術力強化法第19条の適用を受ける特許出願
【前置審査】
(73)【特許権者】
【識別番号】504171134
【氏名又は名称】国立大学法人 筑波大学
(74)【代理人】
【識別番号】100165179
【氏名又は名称】田▲崎▼ 聡
(74)【代理人】
【識別番号】100188558
【氏名又は名称】飯田 雅人
(74)【代理人】
【識別番号】100175824
【氏名又は名称】小林 淳一
(74)【代理人】
【識別番号】100152272
【氏名又は名称】川越 雄一郎
(74)【代理人】
【識別番号】100181722
【氏名又は名称】春田 洋孝
(72)【発明者】
【氏名】工藤 博幸
【審査官】遠藤 直恵
(56)【参考文献】
【文献】特表2003-502130(JP,A)
【文献】特開昭52-110581(JP,A)
【文献】特開平11-019078(JP,A)
【文献】米国特許出願公開第2011/0142316(US,A1)
(58)【調査した分野】(Int.Cl.,DB名)
A61B 6/00-6/14
G06T 1/00
(57)【特許請求の範囲】
【請求項1】
部分的に測定した対象物全体をカバーする全体投影データを対象物に関する先験情報に利用して厳密な画像再構成を行うインテリアCTの画像生成方法であって、
CT測定の幾何学系において撮像対象の内部の関心領域(ROI)を通過する量子ビームを測定してインテリアCT投影データを取得する工程と、
前記CT測定の幾何学系において前記撮像対象の外部にある
任意のセグメントからの前記撮像対象の全体を通過する
ROIを通らないビームを測定して部分的な全体投影データを取得する工程を含んでおり、更に、
取得した前記インテリアCT投影データ
から、前記任意のセグメントからの前記ROIを通らない部分的な全体投影データ
を前記先験情報として利用して、前記ROIを厳密に再構成する処理工程を含んでい
ることを特徴とするインテリアCTの画像生成方法。
【請求項2】
請求項1のインテリアCTの画像生成方法において、前記セグメントは、前記撮像対象を取り囲む曲線に対応する複数の点を含んだ奇数個(1、3、5、…)のセグメントであることを特徴とするインテリアCTの画像生成方法。
【請求項3】
請求項2のインテリアCTの画像生成方法において、前記インテリアCT投影データ取得工程を、360度円軌道ファンビームで行い、前記部分的な全体投影データの取得を、前記インテリアCT投影データを取得する前記360度円軌道ファンビームの円軌道内に含まれている前記セグメントからファンビームで行うことを特徴とするインテリアCTの画像生成方法。
【請求項4】
請求項2のインテリアCTの画像生成方法において、前記インテリアCT投影データ取得工程を、ファンビームショートスキャンで行い、前記部分的な全体投影データの取得を、前記インテリアCT投影データを取得する前記ファンビームショートスキャンの円弧軌道内に含まれている前記セグメントからファンビームで行うことを特徴とするインテリアCTの画像生成方法。
【請求項5】
請求項2のインテリアCTの画像生成方法において、前記インテリアCT投影データ取得工程を、180度平行ビームスキャンで行い、前記部分的な全体投影データの取得を、前記インテリアCT投影データを取得する前記180度平行ビームスキャンの軌道内に含まれている前記セグメントから平行ビームで行うことを特徴とするインテリアCTの画像生成方法。
【請求項6】
請求項3、4、5の何れか一項のインテリアCTの画像生成方法において、前記セグメントはその角度が、単数の場合には30度以上に、3個又は5個の場合には、10度以上になるように設定されていることを特徴とするインテリアCTの画像生成方法。
【請求項7】
請求項1のインテリアCTの画像生成方法において、前記ROIの再構成処理工程は、解析的画像再構成法、逐次近似画像再構成法、統計的画像再構成法の何れか、又は、それらの組み合わせにより実行されることを特徴とするインテリアCTの画像生成方法。
【請求項8】
請求項1のインテリアCTの画像生成方法において、アクティブコリメータを利用して、前記インテリアCT投影データと前記部分的な全体投影データを取得する量子ビームの開口角度を制御することを特徴とするインテリアCTの画像生成方法。
【請求項9】
請求項1のインテリアCTの画像生成方法において、前記部分的な全体投影データとして、位置合わせのためのScout-Viewスキャン投影データを利用することを特徴とするインテリアCTの画像生成方法。
【請求項10】
請求項1のインテリアCTの画像生成方法において、前記部分的な全体投影データを、前記インテリアCT投影データよりも低い分解能の他の量子ビームを照射する測定を利用して取得することを特徴とするインテリアCTの画像生成方法。
【請求項11】
請求項1のインテリアCTの画像生成方法において、量子ビームを対象物に照射して測定される前記インテリアCT投影データ及び前記全体投影データが、量子ビームと物体の相互作用により発生する、吸収、位相シフト、散乱、回折、屈折の少なくとも何れか一つの空間的な物理量分布の、前記量子ビームが通過する直線上の線積分になっていることを特徴とするインテリアCTの画像生成方法。
【請求項12】
部分的に測定した対象物全体をカバーする全体投影データを対象物に関する先験情報に利用して厳密な画像再構成を行うインテリアCTの画像生成方法であって、
CT測定の幾何学系において撮像対象の内部の関心領域(ROI)を通過する量子ビームを測定してトランケーションを含んだインテリアCT投影データを取得する工程と、
前記CT測定の幾何学系において前記撮像対象の外部にある
任意のセグメントからの前記撮像対象の全体を通過する
ROIを通らない量子ビームを測定してトランケーションを含まない部分的な全体投影データを取得する工程を含んでおり、更に、
取得したトランケーションを含んだ前記インテリアCT投影データ
から、トランケーションを含まない前記
任意のセグメントからの前記ROIを通らない部分的な全体投影データ
を前記先験情報として利用して、前記ROIを厳密に再構成する処理工程を含んで
いることを特徴とするインテリアCTの画像生成方法。
【請求項13】
請求項12のインテリアCTの画像生成方法において、前記ROIの再構成処理工程は、解析的画像再構成法、逐次近似画像再構成法、統計的画像再構成法の何れか、又は、それらの組み合わせにより実行されることを特徴とするインテリアCTの画像生成方法。
【請求項14】
請求項12のインテリアCTの画像生成方法において、前記トランケーションを含んだ前記インテリアCT投影データにおけるヒルベルト変換の計算不可能性を、前記トランケーションを含まない前記部分的な全体投影データを用いて計算可能にすることにより、前記ROIを再構成することを特徴とするインテリアCTの画像生成方法。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、物体内部における物理量分布の線積分値を測定してデータ処理により物理量分布を画像生成するCTの画像生成方法に関し、特に、インテリアCTの画像生成方法に関する。
【背景技術】
【0002】
まず、インテリアCT(コンピュータトモグラフィー)の原理について説明する。なお、以下では説明を分かりやすくするためCT測定に使用する量子ビームがX線の場合を想定しているが、他の全ての量子ビームを用いたCTについても同様に適用可能なことは当業者であれば明らかである。CTにおいて、その目的から、例えば心臓や乳房など小さな検査の関心領域(ROI: Region of Interest)の画像のみで十分であり、ROI外部の画像は不要である状況が想定される。通常のCTでは、このような状況であっても、
図1(a)に示すように、ROIを通過するX線のみではなく、対象物全体をカバーするX線を照射して投影データを測定している。しかし、直感的にROIを通らない直線上の投影データはROIの情報を持っていないため、不要なことが予想される。そこで、このような無駄な測定をなくすため、
図1(b)に示すように、ROIを通過するX線のみを照射して投影データを測定してROIのみの画像を生成する方式のCTが、インテリアCTと呼ばれるものである。このようなインテリアCTには、(1)ROI外部の被曝量の低減、(2)検出器やビーム幅の削減、(3)視野に収まらない大きな物体の撮影が可能になること、などの長所がある。また、マイクロ(ナノ)CTなどの小視野を拡大して撮影するCT装置では、試料が視野をはみ出してしまうインテリアCTの状況が頻繁に発生する。
【0003】
このように、インテリアCTにおいては、測定した投影データから画像を生成する画像再構成の問題は、以下のように定式化できる。
図2に示すように、実線で示すROI Ωを通る投影データは全て測定し、一方、点線で示すROI Ωを通らない投影データは全て測定しないものとする。そして、このような不完全な投影データから、如何にしてROI Ωの画像を数学的に厳密に再構成するかが、インテリアCTにおける画像再構成問題である。
【0004】
一方、インテリアCTに関する従来技術においては、まず、非特許文献[1]において、NattererはインテリアCTの画像再構成は解が一意に定まらず数学的に厳密な画像再構成は不可能であることを証明した。そこで、長年様々な近似解法が研究されてきたが、画像の低周波成分に誤差が発生して実用に至らなかった。例えば、代表的な画像再構成法として、各方向投影データの欠損部分を滑らかな関数で外挿してからフィルタ補正逆投影(FBP: Filtered Backprojection)法により画像再構成を行う手法、不完全な投影データから代数的再構成(ART: Algebraic Reconstruction Technique)法や同時反復型逐次近似再構成(SIRT: Simultaneous Iterative Reconstruction Technique)法などの逐次近似法を適用する手法、などが研究された(非特許文献[2]~[4])。
【0005】
図3に、インテリアCTにおいて発生する典型的なアーティファクトの例を示す。画像の低周波成分が歪むシェーディングアーティファクトやROI周辺部の画像値が盛り上がるカッピングアーティファクトがよく発生する。これらの先行研究に対して、2007年になりようやく数学的に厳密な画像再構成法が発見され、インテリアCTにおける画像再構成の新しい方向性が生まれた。これらの厳密解法の重要なキーポイントを一言で要約すると、『インテリアCTの投影データのみでは厳密な画像再構成は不可能だが、対象物に関するごく僅かな先験情報があれば数学的に厳密な画像再構成が可能になる』と述べることができる。
【0006】
更に、上述した厳密解法の発展について述べると、まず、2007年のYeらの論文(非特許文献[5])、2008年のKudoらの論文(非特許文献[6])、2008年のWangらの特許(特許文献[1])では、
図4(a)に示すように、ROI Ω内部の任意小領域B(先験情報領域)における画像値が既知という先験情報があれば、ROI Ωの厳密な画像再構成が可能であることが証明された。次に、2009年のYuらの論文(非特許文献[7])と2014年のWangらの特許(特許文献[2])では、
図4(b)に示すようにROI Ω全体で画像値が区分的一様(完全な一定値を持つ有限個の領域で構成されていること)であれば、ROI Ωの厳密な画像再構成が可能であることが証明された。更に、2010年のYangらの論文(非特許文献[8])では、Yuらの手法を拡張して、
図4(b)に示すように、ROI Ω全体で画像値が区分的多項式(有限次数の多項式で表される変化を持つ有限個の領域で構成されていること)であればROI Ωの厳密な画像再構成が可能であることが証明された。最後に、2017年の工藤らの現在出願中の特許出願(PCT出願)やKudoの論文(非特許文献[9])では、Yuら、Wangら、Yangらの手法において必要な先験情報を大幅に削減することに成功して、ROI Ω内部の任意小領域B(先験情報領域)において画像値が区分的一様または区分的多項式であれば、ROI Ωの厳密な画像再構成が可能であることが証明された。
【0007】
加えて、非特許文献[10]では、KatsevichのFBP法と呼ばれるFBP型の構造を持つ画像再構成法が、非特許文献[11]では、KatsevichのFBP法におけるフィルタ補正投影データが対称性を持つことが、そして、非特許文献[12]では、位相CTにも大きく異なる2つの原理が存在することが示されている。
【先行技術文献】
【特許文献】
【0008】
【文献】米国特許第7,697,658号
【文献】米国特許第8,811,700号
【非特許文献】
【0009】
【文献】Natterer F: The Mathematics of Computerized Tomography. Wiley, 1986
【文献】Ogawa K, Nakajima M, Yuta S: A reconstruction algorithm from truncated projections. IEEE Transactions on Medical Imaging 3: 34-40, 1984
【文献】Ohnesorge B, Flohr T, Schwarz K, Heiken JP, Bae KT: Efficient correction for CT image artifacts caused by objects extending outside the scan field of view. Medical Physics 27: 39-46, 2000
【文献】Ohyama N, Shiraishi A, Honda T, Tsujiuchi J: Analysis and improvement in region-of-interest tomography. Applied Optics 23: 4105-4110, 1984
【文献】Ye Y, Yu H, Wei Y, Wang G: A general local reconstruction approach based on a truncated Hilbert transform. International Journal of Biomedical Imaging 2007: Article ID 63634, 2007
【文献】Kudo H, Courdurier M, Noo F, Defrise M: Tiny a priori knowledge solves the interior problem in computed tomography. Physics in Medicine and Biology 53: 2207-2231, 2008
【文献】Yu H, Wang G: Compressed sensing based interior tomography. Physics in Medicine and Biology 54: 2791-2805, 2009
【文献】Yang J, Yu H, Jiang M, Wang G: High order total variation minimization for interior tomography. Inverse Problems, 26: Article ID 35013, 2010
【文献】Kudo H: Practical interior tomography. Proceedings of International Forum on Medical Imaging in Asia 2017, Paper No. O7-I-2, 2017
【文献】Katsevich A, Analysis of an exact inversion algorithm for spiral cone-beam CT. Physics in Medicine and Biology 47: 2583-2598, 2002
【文献】Noo F, Defrise M, Clackdoyle R, Kudo H: Image reconstruction from fan-beam projections on less than a short scan. Physics in Medicine and Biology 47: 2525-2546, 2002
【文献】百生敦: 位相コントラストX線イメージング. 放射光 10: 273-285, 1997
【発明の概要】
【発明が解決しようとする課題】
【0010】
上述した全ての厳密解法においては、数学的に厳密な画像再構成を可能にするためには、対象物に関する値が既知、区分的一様なことが既知、区分的多項式なことが既知、などの先験情報が必要になる。
【0011】
一方、現実のCTイメージングにおいて先験情報を獲得する方法としては、例えば、(1)同一患者の以前のCT画像や他モダリティで撮影した画像からの推定、(2)他の患者の同じ部位を撮影したCT画像からの推定、(3)上記の特許出願(PCT出願)で提案された近似解法で再構成したアーティファクトを含む画像からユーザが手動または画像解析で推定するなど、多様な方法が考えられる。そして、これらの先験情報が得られるイメージング状況であれば、厳密解法により従来の近似解法よりはるかに高精度な画像再構成が可能である。しかし、撮影前に対象物に関する先験情報が既知であるという状況は必ずしも多くなく、一般性に欠けるものと言っても過言ではない。
【0012】
そこで、本発明は、上述した従来技術における課題、即ち、先験情報が利用できない場合においても、厳密解法により高精度な画像再構成が可能となるより実用的で汎用性にも優れたインテリアCT画像生成方法を提供することをその目的とする。
【課題を解決するための手段】
【0013】
上記の目的を達成するため、本発明によれば、まず、インテリアCTの画像生成方法であって、CT測定の幾何学系において撮像対象の内部の関心領域(ROI)を通過する量子ビームを測定してインテリアCT投影データを取得する工程と、前記CT測定の幾何学系において前記撮像対象の外部にあるセグメントからの前記撮像対象の全体を通過する量子ビームを測定して部分的な全体投影データを取得する工程を含んでおり、更に、取得した前記インテリアCT投影データと前記部分的な全体投影データに基づいて前記ROIを厳密に再構成する処理工程を含んでいることを特徴とするインテリアCTの画像生成方法が提供される。
【0014】
本発明では、上記に記載のインテリアCTの画像生成方法において、前記セグメントは、前記撮像対象を取り囲む曲線に対応する複数の点を含んだ奇数個(1、3、5、…)のセグメントであることが好ましく、更には、前記インテリアCT投影データ取得工程を、360度円軌道ファンビームで行い、前記部分的な全体投影データの取得を、前記インテリアCT投影データを取得する前記360度円軌道ファンビームの円軌道内に含まれている前記セグメントからファンビームで行うことが、又は、前記インテリアCT投影データ取得工程を、ファンビームショートスキャンで行い、前記部分的な全体投影データの取得を、前記インテリアCT投影データを取得する前記ファンビームショートスキャンの円弧軌道内に含まれている前記セグメントからファンビームで行うことが、又は、前記インテリアCT投影データ取得工程を、180度平行ビームスキャンで行い、前記部分的な全体投影データの取得を、前記インテリアCT投影データを取得する前記180度平行ビームスキャンの軌道内に含まれている前記セグメントから平行ビームで行うことが好ましい。
【0015】
また、本発明では、上記に記載のインテリアCTの画像生成方法において、前記セグメントはその角度が、最低でも一つ以上の投影データが含まるように設定されており、更に、上記に記載のインテリアCTの画像生成方法において、前記ROIの再構成処理工程は、解析的画像再構成法、逐次近似画像再構成法、統計的画像再構成法の何れか、又は、それらの組み合わせにより実行されることが好ましく、更に、アクティブコリメータを利用して、前記インテリアCT投影データと前記部分的な全体投影データを取得する量子ビームの開口角度を制御すること、又は、前記部分的な全体投影データとして、位置合わせのためのScout-Viewスキャン投影データを利用すること、又は、前記部分的な全体投影データを、前記インテリアCT投影データよりも低い分解能の他の量子ビームを照射する測定を利用して取得することが好ましく、更に、量子ビームを対象物に照射して測定される前記インテリアCT投影データ及び前記全体投影データが、量子ビームと物体の相互作用により発生する、吸収、位相シフト、散乱、回折、屈折の少なくとも何れか一つの空間的な物理量分布の、前記量子ビームが通過する直線上の線積分になっていることが好ましい。
【0016】
加えて、本発明によれば、上記の目的を達成するため、インテリアCTの画像生成方法であって、CT測定の幾何学系において撮像対象の内部の関心領域(ROI)を通過する量子ビームを測定してトランケーションを含んだインテリアCT投影データを取得する工程と、前記CT測定の幾何学系において前記撮像対象の外部にあるセグメントからの前記撮像対象の全体を通過する量子ビームを測定してトランケーションを含まない部分的な全体投影データを取得する工程を含んでおり、更に、取得したトランケーションを含んだ前記インテリアCT投影データとトランケーションを含まない前記部分的な全体投影データに基づいて前記ROIを厳密に再構成する処理工程を含んでいるインテリアCTの画像生成方法が提供される。
【0017】
また、本発明では、上記に記載のインテリアCTの画像生成方法において、前記ROIの再構成処理工程は、解析的画像再構成法、逐次近似画像再構成法、統計的画像再構成法の何れか、又は、それらの組み合わせにより実行されることが好ましく、又は、前記トランケーションを含んだ前記インテリアCT投影データにおけるヒルベルト変換の計算不可能性を、前記トランケーションを含まない前記部分的な全体投影データを用いて計算可能にすることにより、前記ROIを再構成することが好ましい。
【発明の効果】
【0018】
上述した本発明によれば、先験情報が利用できない場合においても厳密解法により高精度な画像再構成が可能となる、実用性や汎用性においてより優れたインテリアCT画像生成方法を提供するという優れた効果を発揮する。
【図面の簡単な説明】
【0019】
【
図1】通常のCTとの比較により本発明が関わるインテリアCTについて説明する図である。
【
図2】インテリアCTにおける画像再構成問題の定義を説明する図である。
【
図3】インテリアCTにおいて発生する典型的なアーティファクトの例を示す図である。
【
図4】画像値の先験情報が既知である、又は、画像値が区分的一様又は区分的多項式であればROIの厳密な画像再構成が可能であることを説明する図である。
【
図5】本発明の画像再構成方法である先験情報を利用することなくインテリアCTの数学的に厳密な画像再構成を実現する手法を説明する図である。
【
図6】本発明の画像再構成方法を一般化する場合の考え方を説明する図である。
【
図7】上記画像再構成方法を、180度平行ビーム、ファンビームショートスキャン、多角形軌道ファンビームスキャンに適用した場合の説明図である。
【
図8】通常のCTとインテリアCTの場合における360度円軌道ファンビームCTのデータ収集の様子を示す図である。
【
図9】非インテリアCTとインテリアCTの投影データに画像再構成法を適用した際の、違いや計算の破綻について検討した結果を示す図である。
【
図10】2つの対向位置の投影データ座標に対する、KatsevichのFBP法におけるフィルタ補正投影データの対称性を説明する図である。
【
図11】全体投影データの角度範囲Eが大きい場合の状況を説明する図である。
【
図12】全体投影データの角度範囲Eが任意の小さな円弧セグメントの場合の状況を説明する図である。
【
図13】X線源軌道、物体サポート、ROI Ωが原点を中心とする円でない一般的な場合を考えた系1と系2の設定を説明する図である。
【
図14】本発明の有効性を示すために行った数値シミュレーション実験の結果を示す図である。
【
図15】本発明の有効性を示すために行った数値シミュレーション実験の結果を示す図である。
【
図16】本発明の有効性を示すために行った追加の数値シミュレーション実験の結果を示す図である。
【
図17】全体投影データを測定するセグメントを複数の小円弧セグメントに分割して設定する原理を説明する図である。
【
図18】回折格子を用いたX線位相顕微鏡やTalbot型干渉計による方法における幾何学的関係を示す図である。
【
図19】本発明の画像再構成方法を利用した装置の一例である、一般的なX線CT装置の外観構成を示す図である。
【
図20】上記X線CT装置の内部構成の一例を示す図である。
【発明を実施するための形態】
【0020】
以下、本発明の実施の形態について詳細に説明するが、それに先立ち、本発明者等が提案するインテリアCT厳密解法について述べる。
【0021】
<本発明のインテリアCT厳密解法>
本発明では、上述したように従来技術における対象物に関する先験情報を利用することなく、インテリアCTの数学的に厳密な画像再構成を実現する手法が提案されており、これについて以下に詳述する。
【0022】
インテリアCTにおいて解の一意性(画像再構成問題の解が一つに定まること)が成り立ち、厳密な画像再構成を可能にするためには、インテリアCT投影データに加えて、何らかの付加的な情報が必要である。従来の厳密解法では、既述のように、対象物に関する先験情報を利用して解が一意に定まるようにしているが、本発明ではこれに代わるものとして、インテリアCTでは通常は測定しない投影データ、即ち、ROI Ωを通らない余分な投影データを測定して、この補足データを利用する。即ち、提案する厳密解法の重要な特徴は、ROI Ωを通らない必要最小限の余分な投影データを測定することにある。
【0023】
なお、上記の手法はCTの多様な幾何学系(X線源の軌道や光源からのX線の出し方等)にも適用できる一般性を有するが、以下では、最も実用的であるX線源が円軌道上を360度動くファンビームCTについてその原理を説明する。
【0024】
図5(a)に示すように、360度の円軌道のファンビームCTにおいて、ROI Ωを通過するX線のみを照射して投影データを測定するインテリアCTの状況を考える。もちろん、インテリアCT投影データだけでは数学的に厳密なROI Ωの画像再構成は不可能であり、
図5(b)に示すように、付加的な情報として、円軌道の一部の円弧セグメントEから対象物全体をカバーするX線を照射して、全体の投影データを測定する。この部分的な全体投影データ(以下、単に「全体投影データ」とも言う)をインテリアCT投影データに加味して、厳密なROI Ωの再構成ができるようにする。この場合に問題になるのは、大きく以下の2点である。
【0025】
(1)全体投影データを測定する円弧セグメントEは、どの程度短くできるか(もしも広範囲のEが必要であれば、非インテリアCTに近づいてしまい、提案手法の長所は小さいと思われる)。
(2)円弧セグメントEの全体投影データがあれば、厳密なROI Ωの画像再構成が可能であることが数学的に証明できるか。
【0026】
本発明では、上記(1)の問題に対しては、驚くことに、Eは任意の小さな円弧セグメントでよいことが証明できた。
【0027】
また、上記(2)の問題に対しては、過去の対象物に関する先験情報を用いたインテリアCT厳密解法の研究で解の一意性を示すツールとして使用された微分逆投影(DBP: Differentiated Backprojection)法とトランケーションヒルベルト変換を組み合わせた画像再構成法(非特許文献[5]-[9])では証明できないが、FBP法におけるフィルタリング処理にトランケーションヒルベルト変換を導入した新しい画像再構成法を用いることで証明できた(詳細は後述する)。ここで、本発明で明らかになった360度円軌道ファンビームCTにおけるインテリアCT画像再構成問題の解の一意性をまとめると、次のようになる。
【0028】
[解の一意性]
インテリアCT投影データに加えて任意の円弧セグメントE(いくら小さくともよい)の全体スキャン(左右のトランケーションなし)投影データがあれば、インテリアCTの画像再構成の解は一意である。
【0029】
この新しいインテリアCTの厳密解法について、以下において、まず証明することに成功した解の一意性の結果を述べ、次に、その証明を示す。
【0030】
<得られたインテリアCTにおける解の一意性>
ここでは、後に説明するFBP法のフィルタリング処理にトランケーションヒルベルト変換を導入した新しい画像再構成法に基づき証明することに成功した、インテリアCT画像再構成における解の一意性の結果をまとめて述べる。なお、記号の定義として、対象画像(物体)をf(x,y)、ROIをΩで表す。
【0031】
(1)360度円軌道ファンビームCTの場合のインテリアCTの解の一意性について
図5(a)に示す360度円軌道ファンビームCTにおいて、ROI Ωを通過するX線のみを照射して投影データ測定を行うインテリアCTの状況を考える。なお、インテリアCT投影データのみでは画像再構成問題の解は一意に定まらないので、
図5(b)に示すように、円軌道上の任意の小円弧セグメントEから対象物の全体をカバーするX線を照射して全体投影データを測定する。このとき、次の解の一意性が成り立つ。
【0032】
[解の一意性(360度円軌道ファンビームスキャン)]
ROI Ωを通る全ての投影データに加えて、任意の小円弧セグメントE(いくら小さくともよい)から全体投影データを測定すれば、ROI Ωで画像f(x,y)は一意に定まり、Ωの厳密な再構成が可能である。
【0033】
なお、この一意性の問題点は、ここでは、360度円軌道ファンビームCTに特化しているため、他の幾何学系(例えば、180度平行ビームスキャン、ファンビームショートスキャン、非円軌道を用いたファンビームスキャン)を用いたCTに適用できないことである。そこで、任意の幾何学系の場合におけるインテリアCTの解の一意性について、以下に述べる。
【0034】
(2)任意の幾何学系の場合のインテリアCTの解の一意性
そこで、360度円軌道ファンビームCTの場合の解の一意性を、任意の幾何学系で測定した投影データに適用可能なように一般化する。
【0035】
まず、
図6に基づいて、一般化する際の考え方を述べる。任意の幾何学系を用いた場合においても、360度円軌道ファンビームCTの場合においても、物体の直線上の線積分値の集合を測定している点は共通であり、両者の投影データの間には座標変換(投影データ間の座標変換は、CT分野では、リビニングと呼ばれる)の関係が成立する。そこで、任意の幾何学系で測定した投影データを、一旦、
図6(a)に示すように、仮想的な360度円軌道ファンビームCTの投影データに座標変換して、そのデータが360度円軌道ファンビームCTの場合における一意性の条件を満足しているか調べ、満足していれば解の一意性が成立すると判断するのが、基本となる考え方である。この考え方で任意の幾何学系で測定した投影データに適用可能な解の一意性を導出すると、最終的に、次の結論が得られる。
【0036】
[解の一意性(任意の幾何学系)]
以下の2つの両方の条件が満足されるように投影データが測定されていれば、ROI Ωで画像f(x,y)は一意に定まり、Ωの厳密な再構成が可能である。
【0037】
(条件1)ROI Ωを通る全ての投影データが測定されていること(インテリアCTの測定条件)。
【0038】
(条件2)物体外部にある任意の(物体を内部に含むある円に対応する)小円弧セグメントEを通る全ての直線上の投影データ(部分的な全体投影データ)が測定されていること(解を一意に定めるために余分に測定する投影データの条件)。
【0039】
ただし、小円弧セグメントEの長さはいくら短くともよく、どの場所にあってもよい(但し、物体を内部に含むある円に対応するセグメントであることは必要)。
【0040】
上述した2つの幾何学的条件の意味を、
図6(b)に示す。
【0041】
更に、上述の一般化した解の一意性が有効な3つの事例を以下に述べる。いずれも、360度円軌道ファンビームCTの場合の解の一意性では、厳密なROI Ωの再構成可能性が証明できない例になっている。
【0042】
(a)180度平行ビームスキャンの場合
CTのデータ収集法の中で、最も基本的な180度平行ビームスキャンでインテリアCTを実施する場合を考える。動径をr、角度をθとして投影データをp(r,θ)(投影角度範囲-π/2≦θ<π/2)で表す。いま、-ε≦θ≦ε(εは小角度)の角度範囲では(トランケーションなしの)全体投影データが測定され、それ以外ではインテリアCT投影データしか測定されないとする。このとき、小円弧セグメントEを
図7(a)に示すようにとれば上述の解の一意性の条件を満足していることが分かり、ROI再構成の解は一意である。
【0043】
[解の一意性(180度平行ビーム)]
-π/2≦θ<π/2の角度範囲のインテリアCT投影データに加え任意の小角度範囲E(いくら小さくともよい)で全体投影データを測定すれば、ROI Ωで画像f(x,y)は一意に定まり、Ωの厳密な再構成が可能である。
【0044】
なお、この一意性は、360度円軌道ファンビームCTの場合の解の一意性において円軌道の半径を無限にした極限と解釈することができ、180度平行ビームスキャンにおいても360度円軌道ファンビームCTと同様な解の一意性が成立する。なお、この結果は任意の幾何学系に一般化した一意性を用いて初めて証明できる。
【0045】
(b)ファンビームショートスキャンの場合
図7(b)に示すファンビームショートスキャンの場合を考える。円軌道上のX線源の位置をβ∈[-π/2-α
max,π/2+α
max)(α
maxはショートスキャンの条件から決まるオーバースキャン角度、非特許文献[11])、直線検出器上の座標をuとしてファンビーム投影データをg(u,β)で表す。いま、-ε≦β≦ε(εは小角度)の角度範囲で(トランケーションなしの)全体投影データが測定され、それ以外ではインテリアCT投影データしか測定されないとする。このとき、小円弧セグメントEを
図7(b)に示すようにとれば、上述の解の一意性の条件を満足していることが分かり、ROI再構成の解は一意である。
【0046】
[解の一意性(ファンビームショートスキャン)]
-π/2-αmax≦β<π/2+αmaxの角度範囲のインテリアCT投影データに加えて、任意の小角度範囲E(いくら小さくともよい)で全体投影データを測定すれば、ROI Ωで画像f(x,y)は一意に定まり、Ωの厳密な再構成が可能である。
【0047】
(c)多角形軌道を用いたファンビームスキャン:
図7(c)に示す正五角形軌道を用いたファンビームスキャンを考える。正五角形軌道上のX線源の位置をβ∈[-π,π)(βは正五角形の中心から軌道上の点を見た方位角)、直線検出器上の座標をuとしてファンビーム投影データをg(u,β)で表す。今、-ε≦β≦ε(εは小角度で正五角形軌道の一辺がEになるように決定)の角度範囲では(トランケーションなしの)全体投影データが測定され、それ以外ではインテリアCT投影データしか測定されないとする。
【0048】
このとき、小円弧セグメントEを
図7(c)に示すようにとれば、上述の解の一意性の条件を満足していることが分かり、ROI再構成の解は一意である。
【0049】
[解の一意性(多角形軌道を用いたファンビームスキャン)]
-π≦β<πの角度範囲のインテリアCT投影データに加えて、任意の小角度範囲E(いくら小さくともよい)で全体投影データを測定すれば、ROI Ωで画像f(x,y)は一意に定まり、Ωの厳密な再構成が可能である。
【0050】
もちろん、上述した3つの事例は、任意の幾何学系を用いた解の一意性が有効であることを示す一端であり、本発明の手法により解の一意性が証明できるCTのデータ収集法の範囲はかなり広いことは、当業者であれば明らかであろう。
【0051】
(3)測定の実施例
上述したように、本発明の特徴は、通常のインテリアCT投影データに加えて部分軌道Eから全体投影データを測定して、厳密な画像再構成を可能にする点にあり、このような測定を実現する実施例としては、各CT装置のハードウェア構成やイメージング状況に応じて、多様なものが考えられる。以下では、典型的な5つの実施例(a)~(e)について述べる。
【0052】
(a)純粋に2回の測定を行う
1回目の投影データ測定は小角度範囲Eから物体全てをカバーするX線を用いて行い、2回目の測定はROIを通過するX線のみを照射してインテリアCT投影データ全てが測定されるように行う。
【0053】
(b)位置合わせのScout-Viewスキャン投影データを全体投影データに利用する
マイクロ(ナノ)CT装置などでは、撮影を行う前に見たいROIを上手く視野内に収める位置決めの目的で、物体全てが視野内に入るような低分解能のScout-Viewスキャンが行われる。このScout-Viewスキャンの投影データを、部分軌道Eの全体投影データとして利用することができる。
【0054】
(c)複数X線源を用いた装置
近年、複数のX線源を用いた医療用X線CT装置やマイクロ(ナノ)CT装置が実用されている。これらの装置では、複数のX線源のうち1つを低分解能で全体投影データを測定するように、他のX線源でROI Ωを高分解能で撮影するように設定すれば、低分解能投影データを部分軌道Eの全体投影データとして利用することができる。
【0055】
(d)アクティブコリメータを利用する方法
開口角度を動的に変化可能なアクティブコリメータをX線源の前に設置して、部分軌道Eではトランケーションがないように開口角度を大きく、それ以外では開口角度を小さくしてROIを通過するX線のみを照射するように制御し、投影データ測定を行う。
【0056】
(e)以前に撮影した画像から計算で全体投影データを推定する方法
同一の患者の以前に撮影したCT画像や他の患者の同じ部位を撮影したCT画像が利用できる状況であれば、これらの画像から部分軌道Eの投影データのうちトランケーションにより、未観測のデータを計算で推定することができる。
【0057】
<具体的な画像再構成法>
上述したように、ある幾何学系で測定した投影データが解の一意性の条件を満足していれば、具体的な画像再構成法としては(数学的に正しい方法であれば)どのような手法を用いても正しく画像再構成が可能である。したがって、本発明に対応する画像再構成法の選択肢は無数にある。具体的な例としては、以下の代表的な画像再構成法の全てが利用できる。
【0058】
(1)解析的画像再構成法
まず、360度円軌道ファンビームCTの場合には、後述する解の一意性を証明するのに用いたFBP法におけるフィルタリング処理にトランケーションヒルベルト変換を導入した新しい画像再構成法を用いることができる。360度円軌道ファンビームCT以外の幾何学系の場合には、上述の画像再構成法を測定に用いた各々の幾何学系に拡張した直接再構成法、または、一旦360度ファンビーム投影データに座標変換(リビニング)した後に、上述の画像再構成法を適用するリビニング法を用いることができる。もちろん、本発明で証明に利用した画像再構成法と異なる解析的再構成法が将来発見されれば、それで置き換えることも可能であろう。
【0059】
(2)逐次近似画像再構成法
画像の画素値を一列に並べた画像ベクトルと測定したインテリアCT投影データを一列に並べた投影データベクトルを、各々、x, bとすると、画像再構成問題は線形方程式Ax=bを解く問題と定式化できる。この線形方程式をART法、SIRT法、同時反復型ART(SART: Simultaneous Algebraic Reconstruction Technique)法などのCT画像再構成分野における代表的な反復法を用いて解く。
【0060】
(3)統計的画像再構成法
画像xから計算される予測投影データAxと実測投影データbの距離distance(Ax,b)を雑音の統計的性質に基づいて定義しておき、この距離を凸最適化の手法を用いて最小化する。代表的な距離としては、重み付き最小2乗(WLS: Weighted Least Squares)法、放射型CTのPoisson分布対数尤度、透過型CTのPoisson分布対数尤度などがある。距離の最小化に用いる凸最適化の反復法も多数のものが存在する。
【0061】
即ち、解の一意性が保証されているインテリアCT測定データであれば、具体的な画像再構成法の選択には制限はなく(数学的に正しい方法であれば)、何れでもよいということである。
【0062】
<360度ファンビームCTの場合の解の一意性の証明>
以降では、本発明の基礎となっており、その実用可能性が最も高いと思われる360度円軌道ファンビームCTの場合の解の一意性を証明する。なお、任意の幾何学系に一般化する方法は既に説明している。なお、証明は、解の一意性の条件を満足する際に具体的に厳密な画像再構成が可能な手法を、2002年に非特許文献[10]で発見されたKatsevichのFBP法と呼ばれるFBP型の構造を持つ画像再構成法に基づき示すことで行う。
【0063】
(1)基礎事項
まず、証明で考える問題設定と使用する用語・記号の定義を行い、証明に使用する画像再構成法であるKatsevichのFBP法について説明する。
図8に、通常のCTとインテリアCTの場合における360度円軌道ファンビームCTのデータ収集の様子を示す。画像化の対象となる物体をf(x,y)で表す。そして、物体のサポート(対象物が存在する領域)は原点を中心とする半径aの円、ROI Ωは原点を中心とする半径dの円の状況を考える(ただし、a>d)。また、X線源が動く軌道は原点を中心とする半径Rの円a (β)=(Rcos β,Rsin β)(β∈[-π,π))とする(ただし、βはX線源の位置を表す変数でR>a>d)。そして、各X線源の位置βから測定されたファンビーム投影データを直線検出器上の座標uを用いて表す。検出器は対象物の外部に設置するが、画像再構成の研究でよく行われるように、実検出器と平行な原点を通る仮想検出器を考え、その座標uを用いて投影データを表す。即ち、ファンビーム投影データを(u,β)の2変数関数としてg(u,β)で表す。ただし、ROI Ω、物体サポート、X線源軌道の全てが原点を中心とする円と仮定するのは、数式の論述を複雑にしないためであり、X線源軌道の内部に物体サポートが含まれ、物体サポートの内部にROI Ωが含まれる幾何学的配置であれば、以下の証明は積分範囲を読み替えるなどの少ない変更で容易に拡張できることを述べておく(この点は証明の後に再度述べる)。
【0064】
通常のCTでは、物体全てをカバーするようにX線を照射するため、投影データはg(u,β)(-A≦u≦A, β∈[-π,π))であるが(ただし、A=aR/(R2-a2)1/2)、インテリアCTではROI Ωを通過するX線のみが照射されるため、投影データは検出器座標uに関してより狭いg(u,β)(-D≦u≦D, β∈[-π,π))の範囲(ただし、D= dR/(R2-d2)1/2)でしか測定されない(d<aのとき明らかにD<Aである)。この各方向投影データg(u,β)の左右が欠損することをトランケーションと呼び、インテリアCTの画像再構成を困難にしている要因である。
【0065】
本発明では、通常の非インテリアCTを対象として2002年に非特許文献[10]でKatsevichにより発見されたFBP法を証明に用いる。この画像再構成法の具体的な処理手順を、以下にまとめて示す。
【0066】
【0067】
【0068】
【0069】
ただし、p.vは積分のCauchyの主値を表し、式(2)で表される積分変換はヒルベルト変換(Hilbert Transform)と呼ばれる。通常のFBP法では、再構成画像のぼけを補正するフィルタリングはランプフィルタにより一回の畳み込み積分で行われるが、KatsevichのFBP法では、ランプフィルタを微分とヒルベルト変換に分解してStep 1で微分をStep 2でヒルベルト変換を計算している点が異なる。このような特殊なFBP法を証明に用いる理由は、追って明らかになるが、ヒルベルト変換の巧妙な性質を利用するのが目的である。
【0070】
(2)解の一意性の証明
KatsevichのFBP法は、通常の非インテリアCTを想定した画像再構成法である。そこで、非インテリアCTとインテリアCTの投影データにこの画像再構成法を適用した際に、どのような違いがあり、インテリアCTではどこで計算が破綻するかを考えてみる。その結果をまとめて
図9に表の形で示す。
【0071】
図9における(a)欄及び(b)欄においては、微分は近傍データのみで計算できる局所演算であり、積分は積分範囲の全データが存在して初めて計算できることに注意して、上記の式(1)-(3)により計算を行ったときに各関数が計算できる範囲を示している。式(1)の微分は局所演算であるため、非インテリアCTでも、インテリアCTでも、何れの場合でも問題なく計算できる。しかし、式(2)のヒルベルト変換は、-A≦u≦Aの範囲でg(u,β)が測定されている非インテリアCTでは問題なく計算できるが、投影データの左右にトランケーションがあり-D≦u≦Dの範囲でしかg(u,β)が測定されていないインテリアCTでは計算できない。即ち、インテリアCTでは、ヒルベルト変換の箇所でFBP法の計算が破綻して続行できなくなる。
【0072】
そこで、本発明では、以下に述べるように、小さなX線源の角度範囲E={β|-ε≦β≦ε}(εは小角度)において(トランケーションなしの)全体投影データが測定されていれば、上述の計算不可能なヒルベルト変換の箇所を正しく計算してフィルタ補正された投影データgF(u,β)(-D≦u≦D, β∈[-π,π))を求め、逆投影してROI Ωの数学的に厳密な画像再構成が可能なことを示す。
【0073】
図9における(c)欄に、その新手法で計算できる各関数の範囲を示す。なお、以降では簡単のため、全体投影データが測定されるX線源の角度範囲Eは-ε≦β≦εとしているが、Eが他の角度範囲でも証明は全く同じである。さて、それを証明する。まず、証明に使用するヒルベルト変換の重要な2つの性質を述べる。一つ目は、『フィルタ補正投影データの対称性』と呼ばれる性質である。式(1),(2)で微分フィルタとヒルベルト変換フィルタを作用させたg
F(u,β)(いわばランプフィルタを作用させた投影データ)は、次の補題1に示す対称性を持つ。
【0074】
[補題1]
図10(a)に示すように、(u
0,β
0)と(u
1,β
1)を2つの対向位置の投影データ座標とするとき、KatsevichのFBP法におけるフィルタ補正投影データg
F(u,β)は次の対称性を持つ。
【数4】
【0075】
補題1の証明はNooらの非特許文献[11]を参照されたい。式(4)のような対称性が測定された元投影データg(u,β)について成立するのは自明だが、KatsevichのFBP法におけるフィルタリング後の投影データg
F(u,β)についても成立するのが、この驚きの性質である。次に、二つ目の性質は『トランケーションヒルベルト変換の計算可能性』と呼ばれるものである。いま、
図10(b)に示すように、サポート(関数のゼロでない区間)[-A,A]に含まれる一部区間[-D,D]でしか観測されていない関数g(u)があり、その次式で定義されるヒルベルト変換を計算したい状況を考える。
【数5】
【0076】
もちろん、通常式(5)を計算するには全区間[-A,A]のg(u)が必要であるが、一部区間[-D,D]でしか観測されていなくとも、ヒルベルト変換した結果Hg(u')が[-D,D]に含まれる任意の小区間u'∈P⊂[-D,D](先験情報区間)で既知であれば、Hg(u')は[-D,D]全体で一意に定まり正しく計算できるというのがこの性質である。それを以下の補題2にまとめておく。
【0077】
[補題2]
サポートが[-A,A]の関数g(u)が[-D,D]⊂[-A,A]のみで観測され、かつ、g(u)のヒルベルト変換Hg(u')が[-D,D]に含まれる任意の小区間u'∈P⊂[-D,D]で既知であれば、Hg(u')は観測区間[-D,D]全体で一意に定まり観測データg(u)(-D≦u≦D)から正しく計算可能である。
【0078】
補題2の証明とその具体的な計算方法は、Kudoらの非特許文献[6]を参照されたい。非特許文献[6]では、凸射影法と呼ばれる先験情報を拘束条件として用いる反復法により、上述のトランケーションヒルベルト変換を計算している。別の非反復な計算方法として、行列の特異値分解を用いた擬似逆行列計算でも計算可能である。以上の2つの性質を用いて、目的の『X線源の小角度範囲-ε≦β≦ε(εは小角度)において(トランケーションなしの)全体投影データが測定されていれば、上述の計算不可能なヒルベルト変換の箇所を正しく計算してフィルタ補正された投影データgF(u,β)(-D≦u≦D, β∈[-π,π))を一意に求めることができる』事実を示す。証明は以下に述べる2段階の手順で行う。まず、第1ステップでは、全体投影データの角度範囲-ε≦β≦εが十分に大きい(εが大きい)場合について証明を行い、次に、第2ステップで、全体投影データの角度範囲-ε≦β≦εはいくらでも小さくできることを証明する。
【0079】
(a)第1ステップ(全体投影データの角度範囲Eが大きい場合)
まず、
図11に示すような状況を考える。
図11において、実線は全体投影データg(u,β)(-A≦u≦A)が測定されているX線源軌道の角度範囲E={β| -ε≦β≦ε}を表し、点線はインテリアCT投影データg(u,β)(-D≦u≦D)のみが測定されている軌道の角度範囲N=[-π,π)\Eを表す(ただし、小角度εはROI Ωの境界に接する垂直線が軌道から切り取る円弧セグメントから決まる)。この場合、角度範囲Eのg(u,β)は全体投影データであるため、式(1),(2)で微分とヒルベルト変換したg
F(u,β)(-A≦u≦A)が問題なく計算できる。しかし、インテリアCT投影データしか測定されていない角度範囲Nでは、微分は計算できるがヒルベルト変換は計算できない状況になる。これに対して、補題1と補題2を上手く組み合わせて用いると、次の手順によりg
F(u,β)(-D≦u≦D)が正しく計算できる。
【0080】
各β0∈Nのトランケーション投影データに対して、以下のStep 1とStep 2を行う。
【0081】
[Step 1](ヒルベルト変換先験情報構築)
β
0に対向する
図11中のB-Cの範囲の微分とヒルベルト変換が計算されているg
F(u,β)をβ
0側に折り返して集め補題1の対称性を用いると、g
F(u,β
0)の先験情報が構築できる(即ち、小区間u∈P⊂[-D,D]においてg
F(u,β
0)が既知という状況が構築できる)。
【0082】
[Step 2](トランケーションヒルベルト変換)
Step 1で構築したgF(u,β0)(u∈P)を先験情報に用いて、補題2のトランケーションヒルベルト変換の手法により微分フィルタをかけたデータgD(u,β0)(-D≦u≦D)から、フィルタ補正されたデータgF(u,β0)(-D≦u≦D)を計算する。
【0083】
全体投影データの角度範囲Eが大きい
図11に示す設定の場合には、(空集合でない)先験情報区間Pを全てのβ∈Nについて構築でき、この方法で全てのβ∈Nについてg
F(u,β) (-D≦u≦D)を正しく計算できる。よって、フィルタ補正された投影データg
F(u,β)(-D≦u≦D, β∈[-π,π))を式(3)により逆投影して、ROI Ωの数学的に厳密な画像再構成が可能である。
【0084】
(b)第2ステップ(全体投影データの角度範囲Eが任意の小さな円弧セグメントの場合)
第1ステップの証明だけでは、解の一意性の結論で述べた全体投影データの角度範囲が任意の小円弧セグメントEでよいことは示せない。このはるかに強い結果を示すには、もう一工夫が必要である。どのようにしてそれを証明するかを、
図12に示す。
【0085】
第1ステップで述べた対向データを利用してヒルベルト変換可能なX線源軌道の角度範囲を拡大する操作を、複数回繰り返すのが基本アイデアである。
【0086】
図12の例では、元々ヒルベルト変換可能な角度範囲Eからスタートして、一回、第1ステップの操作を行うとヒルベルト変換可能な角度範囲E
1は真ん中の図まで増大し、もう一回、第1ステップの操作を行うと全てのβ∈[-π,π)についてヒルベルト変換可能なようになる。この例では2回の操作で十分だが、一般にEが空集合でない限り、無限回の操作で必ず360度角度範囲β∈[-π,π)についてヒルベルト変換できる状態に到達することが幾何学的に証明できる。よって、式(3)の逆投影によりROI Ωの画像を生成するのに十分なフィルタ補正投影データg
F(u,β)(-D≦u≦D, β∈[-π,π))が計算できるので、ROI Ωを厳密に正しく再構成できる。
【0087】
最終的に得られた解の一意性に関する結論を定理の形でまとめると、以下のようになる。
【0088】
[定理]
図8(b)の設定(X線源軌道は原点中心の半径Rの円、物体サポートは原点中心の半径aの円、ROI Ωは原点中心の半径dの円)の360度円軌道ファンビームCTにおいて、ROI Ωを通る全ての投影データに加えて任意の小円弧セグメントE(いくら小さくともよい)から全体投影データを測定すれば、ROI Ωで画像f(x,y)は一意に定まり、Ωの厳密な再構成が可能である。
【0089】
上述の定理では、証明に現れる数式が複雑になり過ぎないように、X線源軌道、物体サポート、ROI Ωが原点を中心とする円の場合を考えたが、この証明を拡張して以下の2つの系を証明することは容易である。具体的には、証明に用いる画像再構成法の積分範囲などが多少変わってくるが、証明の本質と手順は同じである。
【0090】
[系1]
図13(a)の設定(X線源軌道は半径Rの円、物体サポートは円軌道の内部に含まれる任意の円、ROI Ωは物体サポートの内部に含まれる任意の円)の360度円軌道ファンビームCTにおいて、ROI Ωを通る全ての投影データに加えて任意の小円弧セグメントE(いくら小さくともよい)から全体投影データを測定すれば、ROI Ωで画像f(x,y)は一意に定まり、Ωの厳密な再構成が可能である。
【0091】
[系2]
図13(b)の設定(X線源軌道は半径Rの円、物体サポートは円軌道の内部に含まれる任意の凸形状領域、ROI Ωは物体サポートの内部に含まれる任意の凸形状領域)の360度円軌道ファンビームCTにおいて、ROI Ωを通る全ての投影データに加えて任意の小円弧セグメントE(いくら小さくともよい)から全体投影データを測定すれば、ROI Ωで画像f(x,y)は一意に定まり、Ωの厳密な再構成が可能である。
【0092】
<シミュレーション実験結果>
図14及び
図15に、本発明の有効性を示すために行った数値シミュレーション実験の結果を示す。実験では、人体頭部のCT画像を模擬した数値ファントムを用いて画像再構成を行った。360度あたりのファンビーム投影データ数(View数)は1,000であり、全体投影データを測定する小円弧セグメントE内の投影データ数MがM=300,M=100,M=33の3パターンについて、画像再構成を行った。なお、ここで、M=300はEの角度範囲が108度、M=100はEの角度範囲が36度、M=33はEの角度範囲が12度の場合に相当する。また、上述のセグメントEが1個の場合に加えて、セグメントEが3個の場合についても実験を行った。Eが3個のセグメントの場合は、M=94×3個(34度×3個)とM=25×3個(9度×3個)の設定とした。
図15(a)は、小円弧セグメントEが1個の場合の再構成画像であり、
図15(b)は小円弧セグメントEが3個の場合の再構成画像である。
【0093】
画像再構成に用いた手法は、解の一意性を証明するのに使用したKatsevichのFBP法におけるフィルタリング処理にトランケーションヒルベルト変換を導入して厳密に画像再構成できるようにした厳密解法である。また、比較のため、全ての投影データが全体投影データの場合の通常のFBP法、全ての投影データがインテリアCT投影データの場合に最もよく使用されるローカルFBP法(投影データの欠損部分を滑らかな関数で外挿してからFBP法で画像再構成)についても実験を行った。
【0094】
これらの図からも明らかなように、ローカルFBP法では、ROI周辺部の画像値が盛り上がるカッピングアーティファクトが強く発生しているのに対して、本発明のセグメントEの個数が1個の手法ではM=300とM=100についてはほぼ完全にROIを画像再構成することができている。なお、Eの角度範囲が小さなM=33の場合については、M=300やM=100と比較すると低周波数成分のアーティファクトは多少増大しているが、ローカルFBP法と比較すると、はるかに誤差が少ないことが分かる。また、セグメントEの個数が3個の場合の方が同等な角度範囲のEを用いる1個の場合と比較して、よりアーティファクトが少ない画像再構成が実現できる傾向が見られた。
【0095】
<追加シミュレーション実験結果>
図16に、本発明の有効性を示すために行った追加の数値シミュレーション実験の結果を示す。実験では、医療用X線CT装置で測定した人体頭部のCT実画像を数値ファントムとして画像再構成を行った。180度平行ビームスキャンを想定して、180度あたりの投影データ数(View数)は800と設定して、上述のセグメントEが1個、全体投影データを測定する小円弧セグメントE内の投影データ数MがM=1の場合(即ち全体投影データの数が僅か1個のみ)という最も難しい場合について、画像再構成を行った。なお、ここで、M=1はEの角度範囲が0.225度の場合に相当する。
図16(a)は原画像、
図16(b)は全体投影データを用いない場合のローカルFBP法による再構成画像、
図16(c)は全体投影データを用いない場合の逐次近似法による再構成画像、
図16(d)は1個の全体投影データを用いる本発明の手法による再構成画像である。
【0096】
画像再構成に用いた手法は、段落[0059]で述べた逐次近似画像再構成法である。具体的には、画像の画素値を一列に並べた画像ベクトルと測定したインテリアCT投影データと1個の全体投影データを一列に並べた投影データベクトルを、各々、x, bとして、画像再構成問題を線形方程式Ax=bを解く問題に定式化して、この線形方程式をCT画像再構成分野における代表的な反復法であるART法を用いて解いた。ただし、実装にあたっては、反復が進むにつれて計算誤差が増大する影響を抑制するため、トータルバリエーション(TV: Total Variation)による正則化(非特許文献[7])を弱く施した。
【0097】
この図から明らかなように、全体投影データを用いない通常のインテリアCTにおいては、ローカルFBPでも逐次近似法でもROI周辺部の画素値が盛り上がるカッピング効果と画像の直流成分シフトによるアーティファクトが強く発生して脳内の構造が見えないほど画質が劣化しているのに対して、1個のセグメントEから僅か1個の全体投影データを測定して補足的に用いる本発明の手法ではROIをほぼ完全に画像再構成することができている。本実験では、全体投影データに誤差や雑音が全くない場合を想定したが、全体投影データの測定精度が高い場合には、本実験結果に見られるように1個のセグメントEから測定した僅か1個の全体投影データのみを補足して、インテリアCTにおいて発生するアーティファクトをほぼ完全に除去することができると結論づけられる。
【0098】
<複数の全体投影データを測定する円弧セグメントを用いた方法>
これまでの議論では、画像再構成問題の解を一意に定めるために全体投影データを測定する円弧セグメントEが単一の場合のみを考えてきた。この場合でも、上述の証明の第2ステップで行ったヒルベルト変換可能な角度範囲を拡大する操作を複数回逐次的に繰り返すことで、全体投影データを測定するセグメントEをいくらでも小さくできることは数学的に証明できるが、この第2ステップの操作回数が多いと画像再構成が複雑になり安定性(解が一意に定まりやすく雑音や計算誤差に敏感に影響されないこと)が悪くなることが危惧され、上述のシミュレーション実験でもEが小さいとアーティファクトが増大することが確認された。
【0099】
そこで、この問題点を緩和するため、ヒルベルト変換可能な角度範囲を拡大する操作が一回だけ(証明の第1ステップだけ)でも全体投影データを測定するセグメントEを十分に小さくできる方法として、Eを複数の小円弧セグメントに分割して設定する方法を考案した。
【0100】
図17に、本手法の原理を示す。全体投影データを測定する円弧セグメントEが単一の場合には、ヒルベルト変換可能な角度範囲を拡大する操作を一回で済ませるには、(証明で述べたように)
図17(a)に示すようにEはかなり大きくとることが必要である。これに対して、
図17(b),(c)に示すように、全体投影データを測定する角度範囲Eを奇数n個(n=3,5,…)の等間隔に配置した小円弧セグメントに分割して設定すれば、どの方向のインテリアCT投影データβ∈Nも対向データから一部が分かっている先験情報区間u∈Pが存在する状況を構築することがはるかに容易になり、トータルのEの角度範囲を大幅に小さくすることが可能である(
図17(a)と
図17(b),(c)の比較から明らかである)。即ち、全体投影データを測定するセグメントEを分割して配置することが可能であれば、より小さなEでより安定に画像再構成を行うことができる。
【0101】
<位相X線CTへの拡張>
これまでは、主に物体内部のX線吸収係数分布を画像化する吸収CTを想定して説明を行ってきたが、近年X線が物体を通過する際の位相シフト量分布を画像化する位相CTの実用化が進んでいる。本発明で述べた解の一意性は、かかる位相CTの場合にも全く同じように成立することが示される。以下にその理由を簡潔に述べる。
【0102】
位相CTにも以下の大きく異なる2つの原理が存在する(非特許文献[12])。
【0103】
(1)単色X線を用いたBonse-Hart型干渉計による方法
この方法では、測定した干渉縞パターンから位相復元処理を行って得られるデータは、位相シフト量空間分布の線積分値になる。よって、任意の幾何学系(γ,η)で測定した投影データp(γ,η)を仮想360度円軌道ファンビームCTの投影データg(u,β)に座標変換することができ、これまでの議論と同様に変換されたデータg(u,β)が解の一意性を満足することが厳密にROI再構成可能な条件となる。したがって、解の一意性は全く同じであり、本発明の全ての結論が成立する。
【0104】
(2)回折格子を用いたX線位相顕微鏡やTalbot型干渉計による方法
この方法では、測定した干渉縞パターンやモアレ縞パターンから位相復元処理を行って得られる投影データは、位相シフト量空間分布の線積分値を各光線に垂直な方向に一次微分したデータ(微分位相投影データと呼ばれる)になる。いま、任意の幾何学系(γ,η)で測定した微分位相投影データをp
D(γ,η)で表し、各データが微分される方向の単位法線ベクトルをn(γ,η)とすると、各データp
D(γ,η)と同じ直線に対応するKatsevichのFBP法で現れる式(1)の微分フィルタ投影データg
D(u,β)の間には、次の関係式が成り立つ。
【数6】
【0105】
ただし、a'(β)≡(d/dβ)a(β)=(-Rsin β,Rcos β)である。式(6)の証明は、Nooらの非特許文献[11]を参照されたい。なお、
図18には、式(6)の幾何学的関係を示す。
【0106】
したがって、式(6)を用いて任意の幾何学系で測定した微分位相投影データpD(γ,η)を仮想360度円軌道ファンビームCTにおけるKatsevichのFBP法の微分投影データgD(u,β)に変換することができ、これまでの議論と同様に、変換されたデータgD(u,β)が解の一意性を満足することが厳密にROI再構成可能な条件となる。したがって、解の一意性は全く同じであり本発明の全ての結論が成立する。
【0107】
以上に詳細に述べたように、本発明は、先験情報が既知であるという状況でない場合でも厳密解法によって従来の近似解法よりはるかに高精度な画像再構成が可能な、より一般性のあるインテリアCT画像生成方法であり、当該方法は、物体内部における物理量分布の線積分値を測定してデータ処理により物理量分布を画像生成する原理に基づくあらゆるCTイメージング装置においても適用可能である。
【0108】
なお、CTとは、通常、X線吸収係数分布を画像生成する吸収X線CTのことを指す場合が多く、そこで、以下には、本発明の画像生成方法を適用する装置の最も一般的な例として、X線を用いて被検体の内部の断面像を得るX線CT装置について、図面を参照しながら概略を説明する。
【0109】
図19は、上述した画像生成方法を利用して物体内部における物理量分布の線積分値を測定してデータ処理により物理量分布を画像生成する一般的なX線CT装置の全体外観構成を示す。即ち、X線CT装置は、以下にも述べるX線照射部等の構成要素を収納すると共に、その中央部には被検体が位置される略円筒状の空洞部を備えたガントリー部1と、その上面に被検体を載置する天板(クレードル)4を備えた基台部2と、そして、データ処理装置であるコンピュータ(ここでは図示せず)、得られた画像などを表示するためのディスプレイ装置5や必要な入力を行うためのキーボード6等を含むコンソール部3を備えている。
【0110】
ガントリー部1のハウジング内部やコンソール部3内には、
図20にも示すように、X線CT装置を構成する構成要素が設けられている。その一例として、図示のように、ハウジングの内部には、サンプルに対してX線を扇状に照射するX線発生装置10と、当該装置から照射されて被検体を透過したX線を検出する円弧状のX線検出装置20が、ここでは図示しないが、例えばリング状のフレーム上に取り付けられている。
【0111】
一方、X線発生装置10とX線検出装置20との間の空間には、その上面に被検体を載置(セッティング)する天板30(
図19の符号4に対応)が設けられている。なお、その一部にX線発生装置10とX線検出装置20が取り付けられる部材は、回転駆動部50を介して、例えば、上記ガントリー部1の内部に設けられたモータ等の回転駆動機構により、所定の方向に所定の回転速度で回転する(図の矢印を参照)。一方、被検体を載置するための天板30は、上記X線発生装置10とX線検出装置20の回転面の略中央部の円筒状の空間に対向して配置され、サンプル載置台移動部60により移動される。更に、ガントリー部1には、上記X線発生装置10に対して高電圧を発生して供給するためのX線高電圧部40、モータの回転制御によりX線発生装置10とX線検出装置20が取り付けられる部材を回転駆動する回転駆動部50などが設けられている。
【0112】
上述したX線検出装置20からの検出信号は、データ収集部70に入力されて画像データとして収集され、更に、画像再生部75において、サンプルの内部の断面画像または3次元画像等として再生される。なお、図中の符号76は、画像再生部75においてサンプルの内部の断面画像または3次元画像を再生する際に使用される記憶装置(画像メモリ)である。また、当該画像再生部75により再生されたサンプルの内部の断面画像または3次元画像は、例えば、液晶表示装置等により構成された画像表示部80(
図19の符号5に対応)において表示される。なお、この画像表示部80に、所謂、タッチパネル(図示せず)を組み込むことによれば、当該画像表示部80により、装置の操作に必要な入力を行うことが可能となる。しかしながら、本発明はこれに限定されることなく、当該装置は、当該タッチパネルに代え、キーボード(
図19の符号6に対応)やテンキーやマウス等を備えてもよい。
【0113】
図中の符号90は、上述したX線CT装置を構成する各部の動作を制御するための制御部(
図19の符号3に対応)を示している。より具体的には、例えば、中央演算処理装置(CPU)、そして、RAMやROM等の記憶装置(メモリ)、更には、HDD等の外部記憶装置等により構成されており、そして、記憶装置内に格納した各部の動作を制御するためのソフトウェアやファームウェアに基づいて必要な制御を実行する。
【0114】
そして、上述した本発明であるインテリアCTの画像生成方法は、上記X線CT装置を構成する画像再生部75において、例えば、ソフトウェアとしてRAMやROM等の記憶装置(メモリ)内に格納され、中央演算処理装置(CPU)によって実行される。
【0115】
なお、本発明が応用できる技術はこれに止まらず、例えば、X線を照射した際の位相シフト分布の線積分データから位相シフト分布を画像生成する位相X線CT、体内に投与した放射性薬剤分布の画像を生成する核医学イメージング装置であるPET(ポジトロンエミッションCT)やSPECT(単光子放射型CT)、超音波・マイクロ波・音波・地震波などの波動を用いたCT、電子線CT、投影データからの画像再構成を利用したMRI(磁気共鳴イメージング)などにも応用可能である。即ち、本発明において「物体」または「画像」とは、画像化する物理量の空間分布のことを指し、「投影データ」とは、その直線上の線積分値を表す測定データのことを指す。
【0116】
また、位相シフト、量子ビームの位相シフト、回折、又は回折を含む投影データの数値は、光学素子の追加あるいはその位置変更により取得した単数又は複数の投影データのセットから抽出され、当該抽出された前記量子ビームの位相シフト、回折、又は回折を含む前記投影データの数値を用いて画像を再構成することも可能である。
【0117】
以上には、本発明の種々の実施例になるインテリアCTの画像生成方法について詳細に述べた。しかしながら、本発明は、上述した実施例のみに限定されるものではなく、様々な変形例が含まれることは明らかである。例えば、上記した実施例は本発明を分かりやすく説明するためにシステム全体を詳細に説明したものであり、必ずしも説明した全ての構成を備えるものに限定されるものではない。また、ある実施例の構成の一部を他の実施例の構成に置き換えることが可能であり、また、ある実施例の構成に他の実施例の構成を加えることも可能である。また、各実施例の構成の一部について、他の構成の追加・削除・置換をすることが可能であろう。
【産業上の利用性】
【0118】
本発明は、物体内部における物理量分布の線積分値を測定してデータ処理により物理量分布を画像生成する画像再構成方法、特に、インテリアCTの画像再構成方法を提供する。
【符号の説明】
【0119】
1…ガントリー部、3…コンソール部、4、30…天板、10…X線発生装置、20…X線検出装置、90…制御部