JP2014509037A - 画像データのモデルに基づく処理 - Google Patents

画像データのモデルに基づく処理 Download PDF

Info

Publication number
JP2014509037A
JP2014509037A JP2014501302A JP2014501302A JP2014509037A JP 2014509037 A JP2014509037 A JP 2014509037A JP 2014501302 A JP2014501302 A JP 2014501302A JP 2014501302 A JP2014501302 A JP 2014501302A JP 2014509037 A JP2014509037 A JP 2014509037A
Authority
JP
Japan
Prior art keywords
component
model
image
imaging device
image data
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
JP2014501302A
Other languages
English (en)
Inventor
ウェブスター ステイマン、ジョセフ
エイチ. シエヴェルドセン、ジェフリー
Original Assignee
ザ・ジョンズ・ホプキンス・ユニバーシティー
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by ザ・ジョンズ・ホプキンス・ユニバーシティー filed Critical ザ・ジョンズ・ホプキンス・ユニバーシティー
Publication of JP2014509037A publication Critical patent/JP2014509037A/ja
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/0002Inspection of images, e.g. flaw detection
    • G06T7/0012Biomedical image inspection
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • G06T11/006Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/70Determining position or orientation of objects or cameras
    • G06T7/73Determining position or orientation of objects or cameras using feature-based methods
    • G06T7/75Determining position or orientation of objects or cameras using feature-based methods involving models
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10081Computed x-ray tomography [CT]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing
    • G06T2207/30052Implant; Prosthesis
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00Image generation
    • G06T2211/40Computed tomography
    • G06T2211/424Iterative

Abstract

コンポーネントを含むオブジェクトの画像データを処理する撮像システムである。撮像システムは、画像データを取得する撮像デバイスとプロセッサとを含む。プロセッサは、撮像デバイスから画像データを受信して、コンポーネントのコンポーネントモデルを得て、撮像デバイスの撮像デバイスモデルを得て、コンポーネントモデルと撮像デバイスモデルとに基づいて、非制限的な対象物の関数を構築して、非制限的な対象物の関数と画像データとに基づいて、コンポーネントを含むオブジェクトのモデルを構築し、さらに、モデルに基づいて、コンポーネントを含むオブジェクトの画像を表示する表示デバイスも含まれる。
【選択図】図1

Description

[関連出願の相互参照]
本願は、2011年3月24日提出の米国仮出願第61/467,187号の優先権を主張しており、この内容全体をここに参照として組み込む。
本発明は、政府の援助を受けてなされている(アメリカ国立衛生研究所による許可番号CA−127444)。アメリカ政府も本発明に一定の権利を有している。
本発明は、画像データの処理に関しており、より詳しくは、画像データのモデルに基づく処理に関する。
本明細書で言及される文献、特許出願公開公報、及び特許の全ての内容を、ここに参照として組み込むこととする。
断層撮影においては、多くの場合、画像全体の一部が先験的に分かっている。たとえば整形外科においては、脊柱手術の場合の脊椎頸ネジ及び脊椎頸用ロッド、人口関節のための膝及び臀部のインプラント、及び、外傷の場合の固定用プレート及びネジなどのコンポーネントについては、既に確かな知識がある。画像の補助を受けて行う手術においては、外科手術用ツールを画像視野内に配置する場合も多い。これらコンポーネントが金属製である場合には、該当元素が突出している測定値は、光子が不足することによって信号対雑音比が低くなることにより支障が出る場合がある。同様に、高度に減衰されているコンポーネントの再構成には、ゼロに近い突出した値の数学的反転を行う場合があるので、アルゴリズムがバイアスに対して敏感に反応しすぎるようになる場合がある(多エネルギー効果(polyenergetic effects)によって)。これら効果は両方とも、再構成された画像に縞状のアーチファクト(streak artifacts)を生成する傾向がある(B. De Man, et al, "Metal streak artifacts in X-ray computed tomography: A simulation study," IEEE Trans Nuclear Science, vol. 46, pp. 691-696, 1999; J. F. Barrett and N. Keat, "Artifacts in CT: recognition and avoidance," Radiographics, vol. 24, pp. 1679-91, Nov-Dec 2004)。
こういったアーチファクトが特に厄介な理由は、診断対象のすぐ外の部分の領域である場合が多く、アーチファクトがこの部分で最も顕著に表れる箇所であるからである。金属コンポーネントの近隣の画質が重要になる状況の例には、沈下または骨溶解が示されるインプラントを目視したいような場合(S. D. Stulberg, et al, "Monitoring pelvic osteolysis following total hip replacement surgery: an algorithm for surveillance," J Bone Joint Surg Am, vol. 84-A Suppl 2, pp. 116-22, 2002)、脊柱内に重要な構造がないようにするための、脊椎頸ネジの配置の評価(L. T. Holly and K. T. Foley, "Three-dimensional fiuoroscopy-guided percutaneous thoracolumbar pedicle screw placement. Technical note," J Neurosurg, vol. 99, pp. 324-9, Oct 2003; M. Y. Wang, et al, "Reliability of three-dimensional fluoroscopy for detecting pedicle screw violations in the thoracic and lumbar spine," Neurosurgery, vol. 54, pp. 1138-42; discussion 1142-3, May 2004]、及び、 生検針の案内[B. Daly, et al, "Percutaneous abdominal and pelvic interventional procedures using CT fluoroscopy guidance," AJR Am J Roentgenol, vol. 173, pp. 637-44, Sep 1999)が挙げられる。
金属の縞状アーチファクトを軽減させるためにこれまでに様々な方法が開発された(B. De Man, et al, "Reduction of metal steak artifacts in X-ray computed tomography using a transmission maximum a posteriori algorithm," IEEE Trans Nuclear Science, vol. 47, pp. 977-981, 2000 2000; G. H. Glover and N. J. Pelc, "An algorithm for the reduction of metal clip artifacts in CT reconstructions," Med Phys, vol. 8, pp. 799-807, Nov-Dec 1981; W. A. Kalender, et al, "Reduction of CT artifacts caused by metallic implants," Radiology, vol. 164, pp. 576-7, Aug 1987; H. Li, et al, "Metal artifact suppression from reformatted projections in multislice helical CT using dual-front active contours," Med Phys, vol. 37, pp. 5155-64, Oct 2010; D. D. Robertson, et al, "Total hip prosthesis metal-artifact suppression using iterative deblurring reconstruct0ion," J Comput Assist Tomogr, vol. 21, pp. 293-8, Mar- Apr 1997; G. Wang, et al, "Iterative deblurring for CT metal artifact reduction," IEEE Trans Med Imaging, vol. 15, pp. 657-64, 1996; O. Watzke and W. A. Kalender, "A pragmatic approach to metal artifact reduction in CT: merging of metal artifact reduced images," Eur Radiol, vol. 14, pp. 849-56, May 2004; B. P. Medoff, et al, "Iterative Convolution Backprojection Algorithms for Image-Reconstruction from Limited Data," Journal of the Optical Society of America, vol. 73, pp. 1493-1500, 1983; J. Rinkel, et al, "Computed tomographic metal artifact reduction for the detection and quantitation of small features near large metallic implants: a comparison of published methods," J Comput Assist Tomogr, vol. 32, pp. 621-9, Jul- Aug 2008)。
多くの方法では、金属の計測値を、失われたデータとして考えている。失われたデータは、再構成アルゴリズムから単純に除去したり(B. P. Medoff et al., "Iterative Convolution Backprojection Algorithms for Image-Reconstruction from Limited Data," Journal of the Optical Society of America, vol. 73, pp. 1493-1500, 1983、または、失われたデータの近隣の部分に基づいた値を利用して、埋めたりすることができる(G. H. Glover and N. J. Pelc, "An algorithm for the reduction of metal clip artifacts in CT reconstructions," Med Phys, vol. 8, pp. 799-807, Nov-Dec 1981; W. A. Kalender, et al., "Reduction of CT artifacts caused by metallic implants," Radiology, vol. 164, pp. 576-7, Aug 1987)。しかし、金属成分の正確な知識を利用する場合は稀である。
断層撮影は、一般的に、予備知識を再構成アルゴリズムに組み込むと有利に行うことができる。サンプリング量が少ない状況、信号対雑音比が低い状況では特にそうである。金属の縞状アーチファクトの補正を試みる方法では、その量の空間的な位置の特定、または、投影画像における、金属インプラントの位置の特定が必要となる場合が多い。この位置特定は、通常、金属コンポーネントの減衰係数が高い、という知識を利用して行われる。実際には、これは予備知識の組み込み方として弱い。
ペナルティ付き尤度再構成スキームでは、画像の一般知識をギブスの先例(Gibbs priors)またはペナルティ関数(penalty functions)によって組み込むことができる(K. Lange, "Convergence of EM image reconstruction algorithms with Gibbs smoothing," IEEE Trans Med Imaging, vol. 9, pp. 439-46, 1990; T. Hebert and R. Leahy, "A generalized EM algorithm for 3-D Bayesian reconstruction from Poisson data using Gibbs priors," IEEE Trans Med Imaging, vol. 8, pp. 194-202, 1989; J. B. Thibault, et al., "A three-dimensional statistical approach to improved image quality for multislice helical CT," Med Phys, vol. 34, pp. 4526-44, Nov 2007; J. Wang, et al, "Iterative image reconstruction for CBCT using edge-preserving prior," Med Phys, vol. 36, pp. 252-60, Jan 2009)。
より近年の研究では、PICCS等のアルゴリズムで、生体組織の事前スキャンを組み込んだ非常に具体的な画像の予備知識(very specific image prior)(G. H. Chen, et al., "Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets," Med Phys, vol. 35, pp. 660-3, Feb 2008])、及び、修正されたペナルティ付き尤度法が利用されている(J. Stayman, et al., "Penalized-likelihood reconstruction for sparse data acquisitions with unregistered prior images and compressed sensing penalties," in SPIE Medical Imaging, 2011)。しかしこれら方法で撮像した場合であってもまだ撮像品質が低い。
したがって、画像データの処理を向上させる必要性がある。
本発明の一実施形態におけるコンポーネントを含むオブジェクトの画像データを処理する撮像システムは、画像データを取得する撮像デバイスとプロセッサとを含む。プロセッサは、撮像デバイスから画像データを受信して、コンポーネントのコンポーネントモデルを得て、撮像デバイスの撮像デバイスモデルを得て、コンポーネントモデルと撮像デバイスモデルとに基づいて、非制限的な対象物の関数を構築して、非制限的な対象物の関数と画像データとに基づいて、コンポーネントを含むオブジェクトのモデルを構築し、さらに、モデルに基づいて、コンポーネントを含むオブジェクトの画像を表示する表示デバイスも含まれる。
本発明の一実施形態におけるコンポーネントを含むオブジェクトの画像データを処理する方法は、コンポーネントを含むオブジェクトの画像データを処理する方法であって、 撮像デバイスから画像データを受信する段階と、コンポーネントのコンポーネントモデルを得る段階と、撮像デバイスの撮像デバイスモデルを得る段階と、コンポーネントモデルと撮像デバイスモデルとに基づいて、非制限的な対象物の関数を構築する段階と、非制限的な対象物の関数と画像データとに基づいて、コンポーネントを含むオブジェクトのモデルを構築する段階とを含む。
記載を考慮することでさらなる目的及び利点が明らかになる。
本発明の一実施形態におけるシステムのブロック図である。
本発明の一実施形態におけるコンポーネントを含むオブジェクトのモデルを構築するための処理フローチャートの一例を示す。
本発明の一実施形態における登録パラメータ及び画像パラメータに基づいてコンポーネントを含むオブジェクトのモデルを構築する処理フローチャートの一例を示す。
本発明の一実施形態における2次元のカーネルベースの補間の例を示す。
本発明の一実施形態におけるネジと背景の量とのモデルの例を示す。
本発明の一実施形態における、一連の推定された軸方向の断面(axial slice)及びインプラントの姿勢を例示する。
本発明の一実施形態における収束グラフを示す。
本発明の一実施形態における異なる推定値に基づくネジと背景の量を、軸方向画像及び冠状画像の間で比較した様子である。
本発明の一実施形態における減衰画像を例示する。
本発明の一実施形態における、一方向の単一コンポーネントからなるネジの減衰画像を例示する。
本発明の一実施形態を以下に詳述する。実施形態を説明するにあたり、明瞭に記載するために一定の用語を利用する。しかし本発明は、選択されたこの特定の用語に限定されることを意図はしていない。当業者であれば、本発明の広義の構想を逸脱せずに、他の均等物であるコンポーネントを利用したり、他の方法を開発したりすることもできることを想到するだろう。本明細書で利用する参考文献は、全てをそれぞれ参照として組み込む。
図1は、本発明の一実施形態におけるシステム100のブロック図である。システム100は、撮像デバイス102と、撮像デバイスモデルソース104、コンポーネントモデルソース106、プロセッサ108、及び表示デバイス110を含んでよい。
撮像デバイス102は、オブジェクトの画像データを取得するデバイスである。画像データは、該オブジェクトのそれぞれ異なる断面画像に関するデータであってよい。異なる断面画像は、x線をそれぞれ異なる焦点面に投影させることで得られたものであってよい。たとえば撮像デバイス102は、コンピュータ断層撮影(CT)マシンであってよく、撮像されるオブジェクトは、内部にコンポーネントを埋め込まれたヒトであってよく、画像データは、CTマシンからの投影画像データであってよい。コンポーネントは、金属ネジ等の任意のオブジェクトであってよい。
撮像デバイスモデルソース104は、撮像デバイスモデルの提供元であってよい。撮像デバイスモデルは、コンポーネントを含むオブジェクトの画像データに基づいて、コンポーネントを含むオブジェクトのモデルを構築するために利用されてよい撮像デバイスのモデルであってよい。撮像デバイスモデルは、撮像デバイス102の特徴に基づいて、コンポーネントを含むオブジェクトのモデルを構築する特徴を定義するものであってよい。たとえば特徴によって、モデルが、撮像デバイス102が撮像する焦点面、撮像デバイス102が得る画像の雑音等を考慮に入れて構築されるべきである、と指定されてよい。撮像デバイスモデルソース104は、撮像デバイスモデルを格納するデバイスまたは撮像デバイスモデルを生成するデバイスであってよい。
コンポーネントモデルソース106は、コンポーネントモデルの提供元であってよい。コンポーネントモデルはコンポーネントのモデルであってよい。コンポーネントモデルは、コンポーネントのパラメータ化された特徴及び登録された特徴を定義するものであってよい。パラメータ化された特徴とは、該コンポーネントの減衰特徴を特定する特徴であってよい。たとえばパラメータ化された特徴には、減衰係数、ハウンスフィールド値、材料の組成、または材料の密度が含まれていてよい。パラメータ化された特徴は、コンポーネントの組成が異種であるか同種であるかを示してよい。登録された特徴は、コンポーネントの形状、及びコンポーネントが変形可能である場合、コンポーネントが変化することができる形状、許される位置、向き等を示す特徴であってよい。たとえばコンポーネントモデルは、コンポーネントが画像データに持つ影響をモデル化するために利用可能なコンポーネントの材料組成及び構造を特定するコンピュータ支援設計(CAD)モデルであってよい。コンポーネントモデルソース106は、撮像デバイスモデルを格納するデバイスまたは撮像デバイスモデルを生成するデバイスであってよい。
プロセッサ108は、撮像デバイス102から画像データを受信して、撮像デバイスモデルを撮像デバイスモデルソース104から取得して、コンポーネントモデルソース106からコンポーネントモデルを取得する処理ユニットであってよい。処理ユニットは、コンピュータ等のコンピューティングデバイスであってよい。コンポーネントが画像データに関する画像の視野にあることがわかっている場合には、プロセッサ108は、画像データ、撮像デバイスモデル、及びコンポーネントモデルソースを利用して、コンポーネントを含むオブジェクトのモデルを構築してよい。たとえば、コンポーネントとしてネジが、撮像対象のヒトの内部にあることがわかっている場合には、コンポーネントを含むオブジェクトは、ヒト及びネジになり、プロセッサ108は、ネジを含むヒトのモデルを構築するために、ネジのコンポーネントモデルを利用することができる。
画像データ、撮像デバイスモデル、及びコンポーネントモデルソースに基づいてコンポーネントを含むオブジェクトのモデルを構築することで、画質が上がり、これによりデータの忠実さの低いデータが取得されてもよくなる。コンポーネントを含むオブジェクトのモデルを構築するにあたり、プロセッサ108は、さらに、コンポーネントの登録パラメータを計算してよい。登録パラメータは、コンポーネントの位置、向き、及び形状を特定するものであってよい。プロセッサ108は、さらに、コンポーネントを含むオブジェクトのモデルに基づいてコンポーネントを含むオブジェクトの画像を生成してもよい。
表示デバイス110は、プロセッサ108から生成された画像を受信して、生成された画像を表示するデバイスであってよい。たとえば表示デバイス110は、液晶表示(LCD)デバイスであってよい。表示デバイス110は、さらに、プロセッサ108から登録パラメータを受信して、登録パラメータを表示してもよい。
撮像デバイス102、撮像デバイスモデルソース104、コンポーネントモデルソース106、プロセッサ108、及び表示デバイス110は、それぞれが別のデバイスであっても、1つの統合デバイスであっても、別のデバイス及び統合デバイスの組み合わせであってもよい。撮像デバイス102、撮像デバイスモデルソース104、コンポーネントモデルソース106、プロセッサ108、及び表示デバイス110は、ネットワークを介して(たとえばイントラネットまたはインターネット)または1以上の統合デバイス内の回路を介して通信できてよい。
図2は、本発明の一実施形態におけるコンポーネントを含むオブジェクトのモデルを構築するための処理フローチャート200の一例を示す。最初に、プロセッサ108は、撮像デバイス102から、コンポーネントを含むオブジェクトの画像データを受信して(ブロック202)、コンポーネントモデルソース106からコンポーネントのコンポーネントモデルを取得して(ブロック204)、撮像デバイスモデルソース104から撮像デバイス102の撮像デバイスモデルを取得してよい(ブロック206)。
上述したように、プロセッサ108は、撮像デバイスモデルソース104及びコンポーネントモデルソース106のうち1以上とは別であっても、これらに統合されていてもよい。プロセッサ108が別の実体である場合には、プロセッサは、関連するソースとネットワーク経由で通信することにより、関連するモデルを取得してよい。プロセッサ108が統合されている場合には、プロセッサは、関連するソースとして機能して、モデルを生成することで、関連するモデルを得てよい。
コンポーネントモデルと撮像デバイスモデルとを利用することで、プロセッサ108は、コンポーネントモデルと撮像デバイスモデルとに基づいて、非制限的な対象物の関数(unconstrained objective function)を構築してよい(ブロック208)。非制限的な対象物の関数は、撮像デバイスモデル(システマチックな雑音を含む)とコンポーネントモデルとに基づいて、尤度関数から導出されてよい。コンポーネントモデルは、マスキング処理と加算とを利用してオブジェクトに統合されることで、非制限的な対象物の関数を生成してよい。特に、対象物の関数は、オブジェクトの背景を特定する画像パラメータが、最適化処理において任意の値を自由にとることができるために、非制限的と称されてよい。
プロセッサはさらに、画像データに関連付けられている1以上の画像の視野にあることがわかっている1以上の追加のコンポーネントのための1以上の追加のコンポーネントモデルに基づいて、非制限的な対象物の関数を構築してよい。
非制限的な対象物の関数を利用して、プロセッサ108は、非制限的な対象物の関数及び画像データに基づいて、コンポーネントを含むオブジェクトのモデルを構築してよい(ブロック210)。
図3は、本発明の一実施形態における登録パラメータ及び画像パラメータに基づいてコンポーネントを含むオブジェクトのモデルを構築する処理フローチャート300の一例を示す。図3は、プロセッサ108が、どのように、非制限的な対象物の関数及び画像データに基づいてコンポーネントを含むオブジェクトのモデルを構築することができるかを示すより詳しい実施形態を示していてよい(ブロック210)。
プロセッサ108は、画像パラメータを初期値に設定し、登録パラメータを初期値に設定してよい(ブロック302)。たとえば登録パラメータの初期値は、コンポーネントを、オブジェクトの最長の軸に合わせられたオブジェクトの中心に配置するものであってよい。画像パラメータは、コンポーネントを含むオブジェクトのパラメータ化された特徴を特定するパラメータであってよい。たとえば、画像パラメータは、コンポーネントを含むオブジェクトの背景の減衰を特定してよい。登録パラメータは、コンポーネントの位置及び向きを特定してよい。
初期値を利用して、画像パラメータは一定に保ちつつ、登録パラメータを変化させることで、非制限的な対象物の関数を画像データ向けに最適化することができる(ブロック304)。たとえばコンポーネントの推定位置を移動させてよい。非制限的な対象物の関数を増やすと、新たな位置を維持することができる。非制限的な対象物の関数を減らすと、コンポーネントの推定位置を前の位置に戻すことができる。
初期値及び更新された登録パラメータを利用することで、最適化された登録パラメータを一定に保ちつつ画像パラメータを変化させることで、非制限的な対象物の関数を画像データに最適化することができる(ブロック306)。たとえば、画像の背景の減衰は、画像パラメータを変更させることで推定されてよい。
更新登録パラメータと画像パラメータとの間の最適化プロセスは、所望の回数繰り返し実行することができる。
最終的な登録パラメータ及び画像パラメータを利用することで、プロセッサ108は、コンポーネントを含むオブジェクトのモデルを構築することができる(ブロック308)。たとえばプロセッサ108は、最終的な登録パラメータを利用して、コンポーネントの位置及び向きを決定することができ、最終的な画像パラメータを利用して、画像の各部分からの減衰を決定することができる。
<例>
製造されたコンポーネント(外科ツール、インプラント等)が断層の視野内にある確率は、増え続けている。この理由の1つは、老齢化が進み、人口装具の普及が進んだために、診断的な撮像をうける人が、インプラント(特に臀部及び膝のインプラント)を埋め込まれている確率が高くなっていることにある。別の理由としては、外科手術の補助手段として、手術中に撮像を行うケースが増えており(コーンビームCT)、手術道具及びデバイス(外科用ネジ及びプレート)が対象の生体組織内またはその付近に配置される確率が上がっているからである。これらコンポーネントが金属を含んでいる場合、再構成された量は、コンポーネントの付近またはコンポーネントから離れた組織の画質に悪影響を与える、激しいアーチファクトを含む可能性がある。これらコンポーネントの物理的モデルが存在していることから、この知識を再構成アルゴリズムに統合して、これらアーチファクトを低減させる、またとないチャンスである。
以下に、コンポーネントの配置及び組成に関する既知の情報を明示的に組み込む、モデルに基づくペナルティ付き尤度推定法を説明する。この方法は、既知のコンポーネントそれぞれの生体組織ならびに位置及び姿勢を共同推定する代替最大化法(alternating maximization method)を利用する。提案する方法は、実質的に光子が不足していても、シミュレーションされた脊椎の脊椎頸ネジの再構成において金属のインプラントの境界付近であっても、アーチファクトのほとんどない画像を生成することができることが証明されている。デバイスの姿勢を同時に推定することができるので、品質保証及び処置の検証に有用でありうるデバイスの定量的情報を提供することができる。
手術道具、固定ハードウェア、インプラントの物理的モデルが利用可能なことから、再構成ルーチンに非常に具体的な予備知識を組み込むことができ、さらなる利点を生じる可能性がある。これらコンポーネントは、製造されたオブジェクトを表していることから、これらの材料組成及び構造を完全特定できるCADモデルが利用可能になる。下記で、このような既知の物理的モデルを再構成プロセスに統合するアルゴリズムを説明する。具体的には、オブジェクトの量のモデル自身は、1)(未知の)背景の生体組織の量、及び、2)撮像視野にあることがわかっているコンポーネント(1または複数)の組み合わせである。
コンポーネントの形状及び減衰の分布は分かっているものの(たとえば、デバイスの形状及び材料の内容を特定するCADモデルから導出される)、位置及び姿勢は分かっていない。したがって、各コンポーネント登録に関するパラメータも、オブジェクトのモデルの一部である。結果生成される再構成スキームは、画像パラメータと登録パラメータとの両方の関数である対象物を有することになり、これら2つのパラメータのセットは共同して推定される。この方法は、共同で画像再構成及び登録を行おうとする他の対象物の関数といくらか類似している点がある(J. A. Fessler, "Optimization transfer approach to joint registration/reconstruction for motion-compensated image reconstruction," presented at the ISBI, 2010; M. Jacobson and J. A. Fessler, "Joint estimation of image and deformation parameters in motion-corrected PET," in Proc. IEEE Nuc. Sci. Symp. Med. Im. Conf., 2003, pp. 3290-3294; S. Y. Chun and J. A. Fessler, "Joint image reconstruction and nonrigid motion estimation with a simple penalty that encourages local invertibility.," in Proc. SPIE 7258, Medical Imaging 2009: Phys. Med. Im., 2009, p. 72580U)。
Snyder, Murphyらは、対象物の関数に対する制約によってコンポーネントの正確な知識を組み込む、モデルに基づく方法を開発した(D. L. Snyder, et al., "Deblurring subject to nonnegativity constraints when known functions are present with application to object-constrained computerized tomography," IEEE Trans Med Imaging, vol. 20, pp. 1009-17, Oct 2001; R. J. Murphy et al, "Pose estimation of known objects during transmission tomographic image reconstruction," IEEE Trans Med Imaging, vol. 25, pp. 1392-404, Oct 2006; J. F. Williamson, et al, "Prospects for quantitative computed tomography imaging in the presence of foreign metal bodies using statistical image reconstruction," Med Phys, vol. 29, pp. 2404-18, Oct 2002)。しかし以下に要点をまとめたこの方法は、とりわけ、非制限的な対象物の関数を利用しており、姿勢がわからない任意の数の既知のコンポーネントについて生成されている。
この方法は以下のようにまとめられる。セクションIで、尤度ベースの対象物の関数と、画像の量を、未知の背景の生体組織内の任意の数の既知のコンポーネント(姿勢は分かっていない)の組み合わせとしてモデル化する再構成アルゴリズムとを記載している。セクションIIでは、この方法の実行を説明しており、従来の分析及び繰り返し法比較している。
<A.フォワード・モデル>
透過型断層撮影システムの以下の計測モデルを例にとる。
Figure 2014509037
平均投影計測値
Figure 2014509037
は、ベールの法則により、オブジェクトの線積分
Figure 2014509037
に関連している。各計測値は、この光線の経路の減衰されていないx線のフルエンスと、この検出器部材の利得とに関連しているスカラー値
Figure 2014509037
を有している。離散化されたオブジェクト(たとえばボクセルベースを利用して)については、線積分を以下に表すことができる。
Figure 2014509037
ここでaijは、j番目のボクセル(または)他のベース)μのi番目の線積分に対する貢献を示している。このモデルは、以下のように、行列ベクトル形式にコンパクトに書き換えることができる。
Figure 2014509037
ここでD{・}は、ベクトル引数からの対角行列を形成する演算子を表しており、システム行列Aは、aij全ての集合(collection)を表している。通常、数3は、オブジェクトと計測値との間の完全な関係を表しており、この数3から、対象物の関数を導出することができる。さらに、オブジェクトを、減衰分布が既知である任意の数のコンポーネントμ (n)と、未知の生体組織の背景量μとの組み合わせとしてパラメータ化するよう選択することができる。
Figure 2014509037
ここでW(λ)は、ベクトルがパラメータ化する変換演算子を表している。この変換は一般的なものであるが、ここでは、3次元(3D)の剛体変換(rigid transformation)に着目している。この場合には、ベクトルλは、コンポーネントの向きを定義する6つの要素(変換及び回転値)からなる。数4の第2項は、N個の既知の減衰量μ (n)の減衰値(それぞれが、そのパラメータベクトルλ(n)に従って(剛体)変換されている)の合計を表している。数4の第1項は、各既知のコンポーネントに対応している、1組の変換されたマスクs(n)と乗算して修正された背景の生体組織の貢献μ*を表している。具体的には、各マスクは、各コンポーネントのサポート領域を表している。典型的なマスクは、端部における部分容積効果を除いて、ほとんどがバイナリである。数4及び本明細書全体で利用されている積の演算子Пは、行列またはベクトルのオペランドの、要素ごとの積を表す。
<B.尤度ベースの対象物の関数> 数3及び数4は、平均計測値とパラメータ化された量との間の関係を表している。雑音モデルを選択することで、尤度関数が導き出される。ポアソン・モデルを呼び出して、以下の対数尤度を導出することができる。
Figure 2014509037
ここで、N個のコンポーネント全ての登録パラメータのセットの集合は、以下のように表される。
Figure 2014509037
以下の暗示的な形式で、ペナルティ付きの尤度推定量(penalized-likelihood estimator)を利用するよう選択される。
Figure 2014509037
ここでR(・)は、非常に雑音の多い画像を抑制するペナルティ付き関数を示す。数7における暗示的な推定量は、「既知のコンポーネント再構成」(KCR)法と称すことができる。この推定量は、生体組織の背景、及び、各既知のコンポーネントに関する登録パラメータのセットの両方を共同して近似する。正則化された項は一般的であるが、対の二次ペナルティ(pairwise quadratic penalty)を利用してよい。
Figure 2014509037
ここで、Nは、ボクセルjの周りの一次的周辺ボクセル(first-order neighborhood of voxels)を示す。
<C.変換演算子>
数7の共同推定量は、数4の変換演算子Wを特定する必要がある。数4の公式化は、一般的であり、変換は、ワーピング演算及び他の形態のパラメータ化された変形を表していてよい。同様に、Wは、旋回軸または中心点等の拘束された動きを持つコンポーネントの変換を表していてよい。この例は、剛体変換に着目している。
「移動」画像から「剛体」画像への登録に関する変換演算子は、通常、1)移動画像の格子化された部分を、変換画像のこれらの点の位置へマッピングすること、2)移動画像の各点の周りのボクセルの周辺部分を利用して、変換画像のボクセル値を生成すること、という2つの部分に分解することができる。3Dの剛体変換において、移動画像を変換画像への点のマッピングは、以下のように示されてよい。
Figure 2014509037
[xは異動画像のボクセル位置を表しており、[x' y' z'は変換画像の位置である。ベクトル[xは、画像の量の中心の位置であり、[xは、変換パラメータである。各軸の周りの回転行列は、以下のように定義される。
Figure 2014509037
θ、Ψ、及びφが、各軸周りの回転パラメータを示す。したがって変換のための6つのエレメントのパラメータベクトルは、λ=[x θ Ψ φ]と明示的に定義することができる。数9は、[xから[x' y' z'へのj番目のボクセル位置のマッピングを示しているが、マッピング全体は、3つの「移動画像」ベクトルx、y、及びzと、移動画像及び変換画像のすべての格子点の集合である3つの変換画像x'、y'、及びz'で表される。
点のマッピングが定義されている場合には、補間段階も定義されている。各次元で分離可能なカーネルベースの補間法を考慮することができる。
図4は、本発明の一実施形態における2次元のカーネルベースの補間の例を示す。左上の図は、移動画像(v)の変換画像(u)へのマッピングを示す。右上の図は、移動画像の値の周辺部分に基づいて、1つの変換画像の点を計算することを示している。下の図は変換画像の位置の関数であるカーネルが、各次元沿いに連続して適用されて、この点における変換画像の値を生成することを示している。
変換画像の各位置jにおいて、移動画像の各点の周辺部分にカーネルを適用して、変換画像の値を近似する。一次元(1D)カーネルは、周辺領域の各次元沿いに連続して適用され、合計されることで、各パスが、1つ次元を下げる。たとえば、2Dの場合には、2つの1Dのカーネルを適用して、1つのボクセル値を得る。3Dの場合には、この演算は、元の変換に戻ることを意味しており、数学的には以下のようになる。
Figure 2014509037
u及びvは、それぞれ変換画像及び移動画像を示す。K個の演算子は、全ての位置において各方向に1Dのカーネルを利用することを表しており、この点のマッピングの関数である。多くのカーネルの選択があるが、以下のように定義されるBスプラインカーネル(P. Thevenaz, et al., "Interpolation revisited," IEEE Trans Med Imaging, vol. 19, pp. 739-58, Jul 2000)を利用する。
Figure 2014509037
このカーネルは、1)変換画像の負ではない性質を維持して(負ではない「移動」画像)、負の減衰係数が推定値に入らないようにする、2)k(t)が、微分可能であり、勾配ベースの方法を数7で提案する目的の最適化で利用することができる。
以下のセクションでは、分析勾配情報(analytic gradient information)を利用して、数7の解を近似する。勾配情報は、必ず、推定するパラメータの変換演算子の導関数を含むので、これら等式はここから導出される。具体的には、λの要素に対するWの演算子の導関数を導出する。連鎖律を利用して、以下のように記述することができる。
Figure 2014509037
ここで、λのm番目の導出関数である変換演算子は、
Figure 2014509037
と表す。
Figure 2014509037
演算子は、導関数カーネル(derivative kernel)を適用する以外は、元のカーネル演算子に類似している。数12のBスプラインカーネルについては、導関数カーネルは以下のように計算されてよい。
Figure 2014509037
数13の残りの偏導関数(remaining partial derivatives)は、数9をλmで微分して計算することができる。具体的には以下のように計算することができる。
Figure 2014509037
Figure 2014509037
Figure 2014509037
は、数10に示す回転行列の要素ごとの導関数を表している。
<D.最適化法>
前のサブセクションでは、フォワード・モデルと、数7に示す推定値に含まれるすべての要素を説明した。この推定値は、対象物の関数の暗示的なマキシマイザ(maximizer)であるので、最適化アルゴリズムは、解を繰り返し近似するよう適合されていてよい。前のセクションの展開から、対象物の関数は微分可能であり、推定パラメータ全てに対する導関数の分析的表現が可能であることがわかるだろう。
Figure 2014509037
Figure 2014509037
対象物の関数の分析勾配を前提にすると、パラメータ
Figure 2014509037
について共同で解くためには、いずれの勾配ベースのアルゴリズムであっても利用できると素直に考えることができるが、これは、断層撮像再構成専用のアルゴリズムの長年に及ぶ開発の歴史を無視することになる。したがって、画像の量のパラメータμが、断層に特化した画像更新を利用して更新され、登録パラメータΛが従来の勾配ベースの方法によって更新される、代替最小化方法(alternating minimization approach)を利用することを提案する。
Λは、6N個の要素しか含まず比較的小さいので(剛体変換されたコンポーネント)、BFGS(Broyden-Fletcher-Goldfarb-Shanno)更新を利用して準ニュートン法を利用すると現実的である(H. Erdogan and J. A. Fessler, "Monotonic algorithms for transmission tomography," IEEE Trans Med Imaging, vol. 18, pp. 801-14, Sep 1999)。(Nが大きい場合には、限定されたメモリの変形を利用すればいいだろう(D. C. Liu and J. Nocedal, "On the Limited Memory BFGS Method for Large-Scale Optimization," Mathematical Programming, vol. 45, pp. 503-528, Dec 1989]))。括弧ライン検索(bracketing line search)を含む方法を適用することができる。
画像更新には、放物線サロゲート法(paraboloidal surrogates approach)を利用するとよい。Erdogan及びFesslerが考えた導出法(H. Erdogan and J. A. Fessler, "Monotonic algorithms for transmission tomography," IEEE Trans Med Imaging, vol. 18, pp. 801-14, Sep 1999)では、数5を、限界対数尤度の合計として、以下のような表している。
Figure 2014509037
以下の定義を行うことができる。
Figure 2014509037
Aが負ではない要素をもつ投影演算子を表しており、Wが、負ではない性質を維持するBスプライン近似を利用して選択されており、s(n)、μ、及びμ (n)も負ではないので、t≧0及びd≧0である点に留意されたい(行列Bは、修正されたシステム行列の一種を表している)。t≧0及びl=t+d≧0であるので、g(t)が、(H. Erdogan and J. A. Fessler, "Monotonic algorithms for transmission tomography," IEEE Trans Med Imaging, vol. 18, pp. 801-14, Sep 1999)の理論1の条件を満たしており、最適な曲率をもつg(t)の放物線サロゲート法を以下のように書き換えることができる。
Figure 2014509037
ここで、上付き文字[k]は、k番目の繰り返しの値を示し、cは、放物線サロゲート法の曲率を示しており、[・]+は、負ではないようにする関数を示している(ゼロで切り捨てる)。gの第1及び第2の導関数はそれぞれ
Figure 2014509037
を表し、数19及び数20から直接計算される。数21のサロゲート法を利用して、対数尤度のサロゲート法全体を形成することができるが、(H. Erdogan and J. A. Fessler, "Ordered subsets algorithms for transmission tomography," Phys Med Biol, vol. 44, pp. 2 -51 , Nov 1999)でも利用されたDe Pierroのコンベクシティ・トリック(convesity trick)(A. R. De Pierro, "On the relation between the ISRA and the EM algorithm for positron emission tomography," IEEE Trans Med Imaging, vol. 12, pp. 328-33, 1993)を利用する。
Figure 2014509037
分離可能なサロゲートを得るためには、
Figure 2014509037
とする。
Figure 2014509037
を選択することができる。したがって正則化されていない対象物(unregularized objective)については、以下の更新を行うことができる。
Figure 2014509037
これら更新を(H. Erdogan and J. A. Fessler, "Ordered subsets algorithms for transmission tomography," Phys Med Biol, vol. 44, pp. 2835-51 , Nov 1999)の技術を利用して直接修正して、様々なペナルティ関数を含ませ、順序付けられたサブセットのサブ繰り返し(ordered-subsets subiterations)を適用して収束率を向上させ、予め計算された曲率を利用することで、計算時の負荷を低減させる。
<E.実装>
上述のようにして得た方法を調査するべく、全てのルーチンを呼び出しとともにMatlab(マサチューセッツ州のNatickのMathworks,Inc社)に実装して、KCR方法のより計算の多い側面の外部ライブラリを作成した(custom)。具体的には、マッチングされたSiddon(R. L. Siddon, "Fast calculation of the exact radiological path for a three-dimensional CT array," Med Phys, vol. 12, pp. 252-5, Mar-Apr 1985)及び分離可能なフットプリント(Y. Long, et ah, "3D forward and back-projection for X-ray CT using separable footprints," IEEE Trans Med Imaging, vol. 29, pp. 1839-50, Nov 2010)プロジェクタ及びバックプロジェクタを、CUDAラインブラリを利用してC/C++に実装して、GPUで実行できるようにした。同じ方法は、GPUによる3Dカーネルベースの補間についても行うことができる。
表1に、KCRの代替最大化法全体の擬似コードを示す。コードは、登録パラメータ更新を含む1つ目のものと、順序付けられたサブセットベースの量の更新を含む2つ目のものとからなる、2つの分離したブロックからなる。この擬似コードでは、一般化された繰り返しの上付き文字[k、p、m]を利用して、主要な繰り返し数k、登録ブロックのサブ繰り返し数p、及び、画像更新ブロックサブ繰り返し数mを示す。シンボル[・]は、具体的なインデックスに依存しない値に利用される。
表1は、アルゴリズムの幾つかの自由なパラメータ選択を示す。つまり、登録サブ繰り返し数P、順序付けられたサブセット数M、及び、繰り返し総数を示す。スケーリング係数σは、最初の逆ヘシアン推定値H[0,0,・]の類推に利用され、対象物の関数のスケールに基づいて、または、適切なステップサイズについて大まかなスケールを見つけるという初期ライン検索を利用して、経験的に選択する必要がある。この最適化の他の目立った特徴には、Hを主要なループの繰り返しの間で再利用することと、二次の概算ペナルティ(
Figure 2014509037
でそれぞれ表される)の一次及び二次導関数を利用して画像更新に正規化を含めることとが含まれる。
Figure 2014509037
表1:KCRアルゴリズムが利用する代替最大化の概略
M=1の場合、画像更新は、局所的最大値に単調に収束するはずであり、Wolfe条件を満たすライン検索を利用する場合、登録更新は、局所的最大値に単調に収束するはずである。したがってこれらの条件下では、アルゴリズム全体が局所的最大値に収束することも予期される。通常の順序付けをされたサブセット法が厳密には収束しないのと同様に、M>1の場合には真の収束は予期されない。他のアルゴリズムの置換も可能である(主要ループ外の曲率を予め計算して、順序付けされた個々のサブセットの画像更新間で登録更新をミキシングすることを含む)。
<セクションII:結果>
<A.シミュレーション研究>
KCR法のパフォーマンスを調査するために、シミュレーションされたコーンビーム計算断層撮影(CBCT)システムを利用して複数の研究を行った。この作業は、脊柱手術における脊椎頸ネジ配置をやりやすくするための撮像に着目して行われた(J. H. Siewerdsen, et al., "Volume CT with a flat-panel detector on a mobile, isocentric C-arm: Preclinical investigation in guidance of minimally invasive surgery," Medical Physics, vol. 32, pp. 241-254, Jan 2005)。
図5は、本発明の一実施形態におけるネジと背景の量とのモデルの例を示す。(A)は、多軸の脊椎頸ネジのCADモデルを示す。モデルは、旋回軸の頭部及びネジ山を設けられたネジの2つの主要部品からなる。生体組織の背景の量については、高品質の従来のCTスキャンをカスタムのRANDOファントム(RANDO phantom:Phantom Laboratory,Inc.社製)の腰静脈領域に行うことで得られるデジタルのファントム(digital phantom)を利用することができる。このファントムの多面断層(multiplanar slices)を図5に示す。(B)は、軸方向の断層を示し、(C)は矢状方向の断層を示し、(D)は、冠状の断層を示す。
すべての調査について、300x300x128ボクセルの画像量(1立方ミリメートル)が考慮された。脊椎頸ネジも同様にボクセル化され、(4)のオブジェクトモデルを利用して画像の量に配置された。(4)ではさらに、CADモデルから導出されたボクセル化されたコンポーネントマスクs(n)が必要となる。脊椎頸ネジは、0.3mm−1の線形減衰係数の同質チタンであると仮定する。比較のために、デジタルファントムの柔らかい組織の背景が0.02mm−1であった。システムの形状は、ソースから検出器までの距離が1200mmであり、ソースから軸までの距離が600mmであるフラットパネルベースのC−アームシステムをエミュレートする目的で選択された(S. Schafer, et al., "Mobile C-arm cone-beam CT for guidance of spine surgery: Image quality, radiation dose, and integration with interventional guidance," Medical Physics, vol. 38, pp. 4563-4574, Aug 2011)。プロジェクションは、(1.552x1.552)mmで、360角度を利用して360度においてシミュレーションされ、360x150画素であった。シミュレーションされたデータは、単一エネルギーのフォワード・モデル及びポアソン雑音を利用して生成された。露光レベルは、減衰されていないビームにつき検出素子1つについて、均一的に104個の光子に固定された。以下のセクションで示されるように、この露光レベルでは、従来の方法を利用した場合、光子不足によるアーチファクトが顕著に出る。
3つの別個のネジ配置シナリオが調査された。1つ目が、一方向のモノリシック(単一コンポーネントの)脊椎頸ネジ、2つ目が、双方向のモノリシックネジ、3つ目が、一方向の2コンポーネントの多軸ネジである。登録パラメータΛを、腰椎への1つにネジ(1または複数)の適切な配置をシミュレーションするために選択した。
<B.収束特性>
KCR法の収束特性を調査するために、上述した一方向の単一コンポーネントのネジ配置シナリオからシミュレーションしたデータを再構成した。背景の生体組織の最初の類推は、0.03mm−1の上部減衰で切り捨てられたフィルタリングされた逆投影法(FBP)画像について行われた(ネジは背景の量μには存在すべきではないので、生体組織から脊椎頸ネジを取り除くための大まかな試みにおいて)。登録パラメータについては初めの類推が複数回行われた。解決法から十分遠い複数の初めの類推において、ネジの位置の推定が局所的な最大値で失われた「キャプチャー問題」が観察された。(もっともよく行われるのは、再構成された視野を完全に出て、問題が従来のペナルティ月尤度推定にまで減るようにする)。しかし、投影データの真の位置にあるコンポーネントが少なくとも大体重なっている初期化において、真の姿勢に向かって方法が収束した(雑音による誤りを除いて)。キャプチャー範囲内の最初の類推については、姿勢の推定における残りの誤りが、1ボクセル内に、0.5度(回転パラメータにおいて)未満で104個の光子の雑音レベルにおいて一貫して見られた。
図6は、本発明の一実施形態における、一連の推定された軸方向の断面及びインプラントの姿勢を例示する。本図は、それぞれ登録及び画像サブ繰り返しについてP=10及びM=60を利用して、最初の類推からの始めたサンプルの繰り返しを示しており、最初の4回の繰り返し(Iteration)を示している。最初の類推においては、コンポーネントを任意に配置してよい。アルゴリズムの共同的な性質が明らかであり、外部ループの繰り返しの間に背景画像の縞状アーチファクトの低減とコンポーネントの姿勢両方に改善が見られた。図面に示すように、登録の更新と画像の更新とは、繰り返しを連続する間に同時に実行されてよい。
図7は、本発明の一実施形態における収束グラフを示す。図7の(a)は、60個のサブセットを利用した画像更新と、4つの▽、6つの○、及び10個の□のBFGSの繰り返しを利用した登録の更新とを交互に適用した、KCRの収束プロットを示している。図7の(b)は、10BFGSの繰り返しと、1つの▽、5つの○、10個の□、20個の◇、40個の△、及び60個の☆の順序付けされたサブセットベースの画像の更新とを交互に繰り返し、且つ、動的な数のサブセット
Figure 2014509037
を利用する方法を利用した場合の収束プロットを示している。対象物の関数の差は、登録ブロック(特定シンボル)と画像更新(・)の後にプロットされている。
複数のサブ繰り返しを様々に選択した場合の相対効果をよりよく理解するために、図7の収束プロットを計算した。2つのセットのプロットは、固定数の順序付けされたサブセットの更新に応じて、登録更新の数を変化させた場合の収束レートを表している。具体的には、k番目の繰り返しと、解Φ(後者の値は、P=0の場合に、動的な数のサブセットを利用して、KCRを100回繰り返して近似されているが、これに関しては後述する)における対象物の関数の差をプロットしている。
図7の(a)は、60のサブセットで4、6、及び10の登録更新を利用して3つの方法を行った最初の25回の繰り返しを示す場合を示す。対象物の関数の差は、登録更新(整数kの値同士の中間)、及び、画像更新(シンボル・で表されている)の両方の後に示されている。収束レートは、さらなるBFGSの更新がある度に増加しているように見えるが、P=10を超える数になると、実質的に同じ収束プロットになり、登録更新に重きをおいても得られるものが少ないことが示唆されている。さらには、画像更新を行うと、登録更新よりも対象物の関数の増加が大きくなることもわかった。
図7の(b)では、登録更新数がP=10に固定されており、Mを1、5、10、20、40、及び60サブセットの間で変化させた。これらのケースにおいて、収束レートは、サブセット数に対する依存が強かった(パフォーマンスは従来の順序付けされたサブセット法と類似していた)。すべての方法が対象物を単調に増加させたかにおもえたが、M=60を超えてから平坦部が見られる(順序付けされたサブセットは厳密には収束しないために予期されることではある)。この平坦部を避けて、良好な収束レートを維持するために、サブセット数を動的に変化させる方法を実装する。具体的には、120個のサブセットから始めて、シーケンス90、60、40、30、20、10、3、2、及び1について、4回外部繰り返しを行うたびにサブセット数を低減させる。この方法は図7の(b)の一番下のプロットに表されている。動的なサブセット法は最終的には1つのサブセットになるので、この方法では真の収束が得られることが予期される。サブセット数を変更した場合、4回繰り返すごとに関数値が段階的に向上する点に留意されたい。
<C.画質の実験>
画質の実験について、3つの脊椎頸ネジンプラントのシナリオと3つの異なる再構成法の全てを考慮した。具体的には、P=10にしてKCR法を200回繰り返して、動的なサブセット法を適用した(20回繰り返すごとに上述したシーケンスでMが低減する)。これに比較して、FBP法でも、従来のペナルティ付き尤度(PL)推定値を利用してデータを再構成した。FBP処理では、データは、対数変換の前にゼロより少し高いデータを閾値として前処理された(ゼロのカウント計測値をなくすために)。(メジアンフィルタを利用してこれらゼロの値に補間を行う、ここに示されない代替的な方法を利用した場合であっても、定性的にほぼ同じ結果が生じた)。
PL法及びKCRの両方で、正則化パラメータ(β)の同じ値の二次ペナルティが利用された。対象物同士は非常に類似しているので、パラメータが整合した場合には、これら2つの方法によって同様の解像度特性の画像が生成される。ペナルティ尤度法(KCRを含む)の解像度特性は、一般的に空間変形例であり(space-variant)、オブジェクトに依存している特性によって、正確な解像度のマッチングは難しくなる。βは、FBP画像同士の空間解像度を定性的に整合させるように選択される。
図8は、本発明の一実施形態における異なる推定値に基づくネジと背景の量を、軸方向画像及び冠状画像の間で比較した様子である。図8は、KCRの画質を、FBP及び従来のPL推定の画質と比較したものである。各方法及び各脊椎頸ネジンプラントについて軸方向画像及び冠状画像が提示されている。一番上は、一方向の単一コンポーネントのネジであり、中間が、二方向の単一コンポーネントのネジであり、一番下が、一方向の2コンポーネントの多軸ネジである。すべての画像の窓及びレベルは、それぞれ50及び150HUである。真の画像(左)のオーバレイは、コンポーネントの真の姿勢を示しており、KCR再構成画像(右)のものは、登録/再構成同時プロセスにより推定される姿勢を示している。
図8では、結果は、左の列に示す新のデジタルファントム量とともに示されている。前述と同様に、コンポーネントの位置もオーバレイを利用して示されている。このオーバレイは、真の量の画像及びKCR画像に行われている(しかし、FBPまたは従来のPLでは、コンポーネント登録モデルを含まず、画像モデルしか含まないために、行われない)。
図9は、本発明の一実施形態における減衰画像を例示する。単一コンポーネントの一方向のネジの場合の減衰全体(μ)の図では、(左)が真のファントムであり、(右)がKCRである。KCRは、背景の減衰(μ)と登録パラメータΛの両方を推定するので、結果は従来の減衰画像としても(エラー!、参照元が見つかりません(右))、オーバレイとしても(エラー!、参照元が見つかりません8(右))、示すことができ、表示の好み及び動的な範囲に合致するほうを選択することができる。すべての画像の窓及びレベルは、それぞれ800及び300HUである。
そのほうが望ましければ、いつでも従来の減衰画像を、(4)オーバレイのないKCRを選択することで示すことができる(たとえば、図9に示す単一コンポーネントの一方向のネジの場合)。各再構成技術及び各ネジ配置シナリオについて軸方向の画像及び冠状画像が示されている。冠状の断層は、茎までのネジの配置を示すために選択されている。この表示は、ネジが皮質骨までに達しているかを評価する際の診断においては特に重要である(横方向に侵害している場合(次善のネジの配置)、または医療的な侵害の場合(脊柱管に侵入している場合))。
光子の露光レベルが10の場合、ネジの減衰が高いことによる光子の不足によるアーチファクトは、FBP再構成において非常に顕著である。これらのアーチファクトは、軸方向の画像の周辺領域の中低度のオントラストの詳細を極度に失わせ、双方向のネジの配置シナリオにおいて最悪の影響がでる(このシナリオでは、金属量が増加することで光子の不足量が最大になる)。これと比較して、PL画像では大きな向上がみられ、縞状のアーチファクトの殆どが消えている。しかし、まだ特にネジの境界付近では残余アーチファクトが顕著であり、これにより椎骨の多くがあいまいであり、冠状画像に冠状の侵入があっても評価し難い。KCR法は、殆どアーチファクトのない画像を生成している。残余アーチファクトはいくらかあるものの(そのほとんどはネジのヘッドあたりにある)、アーチファクトのコントラストは、FBP及びPL両方と比較して、非常に低い。さらに、冠状画像からは、ネジの配置において冠状の侵入がないことがよくわかり、真の量の画像によく合致している。KCR法の画質は、ネジの境界付近であっても高く維持されている。
図10は、本発明の一実施形態における、一方向の単一コンポーネントからなるネジの減衰画像を例示する。KCRの正則化パラメータを変化させた効果が示されている。明らかな雑音対解像度のトレードオフ関係が示されているが(βが大きくなると、雑音が低くなり、空間解像度が犠牲になる)、画像全てにおいて、概ね縞状のアーチファクトがない。画像全ての窓及びレベルは、それぞれ500及び150HUである。
図10は、正則化パラメータの範囲内の、一方向の単一コンポーネントの例を示す。具体的には、前の研究と同じ正則化パラメータを用いて、βの値を10倍大きくしたり10分の1にしたりした。正則化のレベルを変化させた結果、βを小さくすると雑音が大きくなると予期され、βを大きくすると空間解像度が低下すると予期される。しかし画像全てにおいて、ネジによる光子不足による縞状のアーチファクトは基本的になかった。このことは、縞状のアーチファクトがなくなったのが単に正則化を行ったからではなく、ネジの構造及び組成に関する予備知識の統合が、実際に再構成を向上させる、ということを示唆している。正則化の値が低いと、よりゆっくり収束するように思われ、量の推定に雑音が多くなった場合、KCRの対象物関数は、局所的な最少値(特に登録パラメータに関して)の影響を受けやすいことがわかるだろう。変換推定の誤りは、図10の、最低から最高のβで順序付けた3つの正則化のケースについて、0.31、0.13、及び0.05ボクセルであった。同様に、平均絶対角度誤りは、0.24、0.09、及び0.04度であった。したがって、βのこの特定の範囲においては、変換推定は、大きな正則化を行ったことの恩恵を受けている。
<セクションIII:考察>
この例は、オブジェクトのモデルに既知のコンポーネントを組み込むことを可能とするKCRを参照する新規な断層撮影再構成法を説明する。コンポーネントの位置及び姿勢が、背景の減衰と同時に推定される。この技術の利用例は、脊柱手術において金属の脊椎頸ネジを含む画像に示したが、この技術は、形状及び組成がわかっているいずれの種類のコンポーネント(異質な組成のオブジェクトも含む)にも利用可能であるという意味で汎用性がある。シミュレーションされた脊椎頸ネジの再構成では、KCRは、光子不足に関するアーチファクトを大幅に低減させ、インプラントの境界までの生体組織をクリアに視覚化しており、従来のFBP及びPL法が実質的にアーチファクトを含んでいるのと対照的である。したがい、KCR法は、介入的なガイド(たとえば脊柱手術等)及び外科手術用工具が3D画質を曖昧にする他の介入のためのコーンビームCT(CBCT)ベースの撮像分野において、及び、大きな将来性が見込めると思われる([J. H. Siewerdsen, et ah, "Volume CT with a flat-panel detector on a mobile, isocentric C-arm: Pre-clinical investigation in guidance of minimally invasive surgery," Medical Physics, vol. 32, pp. 241-254, Jan 2005)。さらに、KCRは、光子の不足するような場合において手術を成功させる可能性があるために、線量(dose)を低減させるための重要な手段となるだろう。KCRは、臀部及び膝のインプラント及びCTベースの生検針の案内等の他の用途においても重要だろう。端部を保存するペナルティ/予備知識を利用するKCR方法を単純に拡大させる方法(K. Lange, "Convergence of EM image reconstruction algorithms with Gibbs smoothing," IEEE Trans Med Imaging, vol. 9, pp. 439-46, 1990; C. R. Vogel and M. E. Oman, "Fast, robust total variation-based reconstruction of noisy, blurred images," IEEE Trans Image Process, vol. 7, pp. 813-24, 1998)も利用することができ、これにより、雑音と解像度のトレードオフ関係がさらに向上する。
上述した基本的なKCRの枠組みは、x線のスペクトルに関する多エネルギービーム及び物理学(たとえば分散)を利用しないフォワード・モデルを利用している。しかし、多エネルギーのフォワード・モデル(B. De Man, et al., "An iterative maximum-likelihood polychromatic algorithm for CT," IEEE Trans Med Imaging, vol. 20, pp. 999-1008, Oct 2001; I. A. Elbakri and J. A. Fessler, "Statistical image reconstruction for polyenergetic X-ray computed tomography," IEEE Trans Med Imaging, vol. 21, pp. 89-99, Feb 2002; I. A. Elbakri and J. A. Fessler, "Segmentation-free statistical image reconstruction for polyenergetic x-ray computed tomography with experimental validation," Phys Med Biol, vol. 48, pp. 2453-77, Aug 7 2003; S. Srivastava and J. A. Fessler, "Simplified statistical image reconstruction algorithm for polyenergetic X-ray CT," 2005 IEEE Nuclear Science Symposium Conference Record, Vols 1-5, pp. 1551-1555, 2005)を組み込むこともできる。この修正は、ビーム硬化が金属アーチファクトの主要な貢献要素であることから、実際の医療データに対する用途には重要だろう。同様に、公知の分散部(known scatter fraction)をフォワード・モデルに組み込むことは簡単である。
一定のコンポーネントの知識が不正確である場合には(たとえば、外科手術時に脊柱の固定に利用される湾曲する金属ロッド、損耗が著しいので手術で直すことが必要な膝または臀部のインプラント、または、挿入時に撓む針の案内)、このKCRの枠組みに対する実装に限界があることが想定される。しかし、(9)のマッピングによって、変形可能な登録モデルに簡単な修正を加えて、これらのシナリオのいくつかを対応可能とすることができるだろうことも付け加えておく(特に、差が、比較的低いパラメータの空間の次元で表されうるような場合)。
上でまとめた結果は、KCRから生じる画質に重きをおいて説明されたが、推定値からは、画像量と登録パラメータという2つの出力がある。画像の量が一次産物ではあるが、姿勢の値の正確な推定自身も、体内に導入されたコンポーネントの位置を知るための医療的な価値があるだろう。たとえばKCR法の姿勢の出力によって、実際の道具を配置するための術前のプランを比較したり、外科手術の正確さを定量化したりする、といった有用性があるだろう。

Claims (17)

  1. コンポーネントを含むオブジェクトの画像データを処理する撮像システムであって、
    画像データを取得する撮像デバイスと、
    前記撮像デバイスから前記画像データを受信して、前記コンポーネントのコンポーネントモデルを得て、前記撮像デバイスの撮像デバイスモデルを得て、前記コンポーネントモデルと前記撮像デバイスモデルとに基づいて、非制限的な対象物の関数を構築して、前記非制限的な対象物の関数と前記画像データとに基づいて、前記コンポーネントを含む前記オブジェクトのモデルを構築するプロセッサと、
    前記モデルに基づいて、前記コンポーネントを含む前記オブジェクトの画像を表示する表示デバイスと
    を備える撮像システム。
  2. 前記プロセッサは、前記コンポーネントを含む前記オブジェクトの前記モデルの構築を、前記非制限的な対象物の関数において、登録パラメータ及び画像パラメータを変化させることに基づいて、前記画像データの前記非制限的な対象物の関数を解くことに基づいて行い、
    前記コンポーネントを含む前記オブジェクトの前記モデルの構築は、解かれた前記非制限的な対象物の関数についての前記画像パラメータ及び前記登録パラメータに基づいて行われる、請求項1に記載の撮像システム。
  3. 前記プロセッサは、前記非制限的な対象物の関数を解くことを、
    前記画像パラメータを一定に保ちつつ、前記登録パラメータを変化させることに基づいて、前記画像データの前記非制限的な対象物の関数を最適化することと、
    最適化された前記登録パラメータを一定に保ちつつ、前記画像パラメータを変化させることに基づいて、前記画像データの前記非制限的な対象物の関数を最適化することに基づいて行う、請求項2に記載の撮像システム。
  4. 前記登録パラメータは、前記コンポーネントの位置及び向きを特定する、請求項2または3に記載の撮像システム。
  5. 前記画像パラメータは、前記コンポーネントを含む前記オブジェクトのパラメータ化された特徴を特定する、請求項2から4のいずれか一項に記載の撮像システム。
  6. 前記コンポーネントモデルは、前記コンポーネントのパラメータ化された特徴と登録された特徴とを定義する、請求項1から5のいずれか一項に記載の撮像システム。
  7. 前記撮像デバイスモデルは、前記撮像デバイスの特徴に基づいて、前記コンポーネントを含む前記オブジェクトの前記モデルを構築するための特徴を定義する、請求項1から6のいずれか一項に記載の撮像システム。
  8. 前記非制限的な対象物の関数の構築は、前記コンポーネントモデルと、前記コンポーネントを含む前記オブジェクト内の第2のコンポーネントの第2のコンポーネントモデル、及び、前記撮像デバイスモデルに基づいて行われる、請求項1から7のいずれか一項に記載の撮像システム。
  9. コンポーネントを含むオブジェクトの画像データを処理する方法であって、
    撮像デバイスから前記画像データを受信する段階と、
    前記コンポーネントのコンポーネントモデルを得る段階と、
    前記撮像デバイスの撮像デバイスモデルを得る段階と、
    前記コンポーネントモデルと前記撮像デバイスモデルとに基づいて、非制限的な対象物の関数を構築する段階と、
    前記非制限的な対象物の関数と前記画像データとに基づいて、前記コンポーネントを含む前記オブジェクトのモデルを構築する段階と
    を備える方法。
  10. 前記コンポーネントを含む前記オブジェクトの前記モデルを構築する段階は、
    前記非制限的な対象物の関数において、登録パラメータ及び画像パラメータを変化させることに基づいて、前記画像データの前記非制限的な対象物の関数を解く段階を含み、
    前記コンポーネントを含む前記オブジェクトの前記モデルを構築する段階は、
    解かれた前記非制限的な対象物の関数についての前記画像パラメータ及び前記登録パラメータに基づいて行われる、請求項9に記載の方法。
  11. 前記非制限的な対象物の関数を解く段階は、
    前記画像パラメータを一定に保ちつつ、前記登録パラメータを変化させることに基づいて、前記画像データの前記非制限的な対象物の関数を最適化する段階と、
    最適化された前記登録パラメータを一定に保ちつつ、前記画像パラメータを変化させることに基づいて、前記画像データの前記非制限的な対象物の関数を最適化する段階と
    を含む、請求項10に記載の方法。
  12. 前記登録パラメータは、前記コンポーネントの位置及び向きを特定する、請求項10または11に記載の方法。
  13. 前記画像パラメータは、前記コンポーネントを含む前記オブジェクトのパラメータ化された特徴を特定する、請求項10から12のいずれか一項に記載の方法。
  14. 前記コンポーネントモデルは、前記コンポーネントのパラメータ化された特徴と登録された特徴とを定義する、請求項10から13のいずれか一項に記載の方法。
  15. 前記撮像デバイスモデルは、前記撮像デバイスの特徴に基づいて、前記コンポーネントを含む前記オブジェクトの前記モデルを構築するための特徴を定義する、請求項10から14のいずれか一項に記載の方法。
  16. 前記非制限的な対象物の関数の構築は、前記コンポーネントモデルと、前記コンポーネントを含む前記オブジェクト内の第2のコンポーネントの第2のコンポーネントモデル、及び、前記撮像デバイスモデルに基づいて行われる、請求項10から15のいずれか一項に記載の方法。
  17. コンピューティングプラットフォームに実行されると、前記コンピューティングプラットフォームに、コンポーネントを含むオブジェクトの画像データを処理する方法を含む処理を実行させるプログラムであって、
    前記方法は、
    撮像デバイスから画像データを受信する段階と、
    前記コンポーネントのコンポーネントモデルを得る段階と、
    前記撮像デバイスの撮像デバイスモデルを得る段階と、
    前記コンポーネントモデルと前記撮像デバイスモデルとに基づいて、非制限的な対象物の関数を構築する段階と、
    前記非制限的な対象物の関数と前記画像データとに基づいて、前記コンポーネントを含む前記オブジェクトのモデルを構築する段階と
    を含む、プログラム。
JP2014501302A 2011-03-24 2012-03-26 画像データのモデルに基づく処理 Pending JP2014509037A (ja)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
US201161467187P 2011-03-24 2011-03-24
US61/467,187 2011-03-24
PCT/US2012/030578 WO2012129566A2 (en) 2011-03-24 2012-03-26 Model-based processing of image data

Publications (1)

Publication Number Publication Date
JP2014509037A true JP2014509037A (ja) 2014-04-10

Family

ID=46880088

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2014501302A Pending JP2014509037A (ja) 2011-03-24 2012-03-26 画像データのモデルに基づく処理

Country Status (5)

Country Link
US (1) US9177374B2 (ja)
EP (1) EP2689414A4 (ja)
JP (1) JP2014509037A (ja)
CA (1) CA2831079A1 (ja)
WO (1) WO2012129566A2 (ja)

Families Citing this family (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8897528B2 (en) * 2006-06-26 2014-11-25 General Electric Company System and method for iterative image reconstruction
JP6122269B2 (ja) * 2011-12-16 2017-04-26 キヤノン株式会社 画像処理装置、画像処理方法、及びプログラム
US9224205B2 (en) * 2012-06-14 2015-12-29 Qualcomm Incorporated Accelerated geometric shape detection and accurate pose tracking
US9610084B2 (en) 2012-09-12 2017-04-04 Peter Michael Sutherland Walker Method and apparatus for hip replacements
DE102013204552B4 (de) * 2013-03-15 2023-09-21 Siemens Healthcare Gmbh Verfahren zur artefaktfreien Wiedergabe von Metallteilen in dreidimensional rekonstruierten Bildern
WO2014160766A1 (en) * 2013-03-26 2014-10-02 The Johns Hopkins University Task-based source-detector trajectories for tomographic imaging
CN104517263B (zh) * 2013-09-30 2019-06-14 Ge医疗系统环球技术有限公司 减少计算机断层扫描图像重构中伪像的方法和装置
CN105374012B (zh) * 2014-08-27 2018-11-27 通用电气公司 用于消除由性能差异的探测器单元所致的条状伪影的方法
JP6413927B2 (ja) * 2015-05-25 2018-10-31 コニカミノルタ株式会社 動態解析装置及び動態解析システム
US10909732B2 (en) * 2016-01-29 2021-02-02 The General Hospital Corporation Systems and methods for joint image reconstruction and motion estimation in magnetic resonance imaging
US11403792B2 (en) 2017-06-19 2022-08-02 Washington University Deep learning-assisted image reconstruction for tomographic imaging
JP6862310B2 (ja) * 2017-08-10 2021-04-21 株式会社日立製作所 パラメータ推定方法及びx線ctシステム
FR3080700B1 (fr) * 2018-04-26 2021-04-02 Centre Nat Rech Scient Procede et dispositif de controle non-destructif d'une piece
DE102019004303A1 (de) 2019-06-18 2020-12-24 Ziehm Imaging Gmbh Verfahren und Vorrichtug der medizinischen Bildgebung zur Darstellung eines 3D-Volumens mit wenigstens einen eingebrachten Fremdobjekt
US11348259B2 (en) * 2020-05-23 2022-05-31 Ping An Technology (Shenzhen) Co., Ltd. Device and method for alignment of multi-modal clinical images using joint synthesis, segmentation, and registration
US11890124B2 (en) * 2021-02-01 2024-02-06 Medtronic Navigation, Inc. Systems and methods for low-dose AI-based imaging

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6898302B1 (en) 1999-05-21 2005-05-24 Emory University Systems, methods and computer program products for the display and visually driven definition of tomographic image planes in three-dimensional space
US7542791B2 (en) * 2003-01-30 2009-06-02 Medtronic Navigation, Inc. Method and apparatus for preplanning a surgical procedure
US7369695B2 (en) 2004-08-20 2008-05-06 General Electric Company Method and apparatus for metal artifact reduction in 3D X-ray image reconstruction using artifact spatial information
US7778392B1 (en) 2004-11-02 2010-08-17 Pme Ip Australia Pty Ltd Method of reconstructing computed tomography (CT) volumes suitable for execution on commodity central processing units (CPUs) and graphics processors, and apparatus operating in accord with those methods (rotational X-ray on GPUs)
US7840249B2 (en) 2004-11-24 2010-11-23 University Of Iowa Research Foundation Clinical micro-CT (CMCT) methods, techniques and apparatus
US20120150243A9 (en) * 2006-08-31 2012-06-14 Catholic Healthcare West (Chw) Computerized Planning Tool For Spine Surgery and Method and Device for Creating a Customized Guide for Implantations
JP4818959B2 (ja) * 2007-03-14 2011-11-16 富士フイルム株式会社 断層画像処理方法および装置ならびにプログラム
US8644568B1 (en) * 2008-07-25 2014-02-04 O.N.Diagnostics, LLC Automated patient-specific bone-implant biomechanical analysis
US8366719B2 (en) * 2009-03-18 2013-02-05 Integrated Spinal Concepts, Inc. Image-guided minimal-step placement of screw into bone
JP5611324B2 (ja) * 2009-04-25 2014-10-22 シーメンス アクチエンゲゼルシヤフトSiemens Aktiengesellschaft インプラントと生物の骨との相対位置を評価するシステム
US8842893B2 (en) * 2010-04-30 2014-09-23 Medtronic Navigation, Inc. Method and apparatus for image-based navigation

Also Published As

Publication number Publication date
EP2689414A4 (en) 2014-09-10
US9177374B2 (en) 2015-11-03
CA2831079A1 (en) 2012-09-27
EP2689414A2 (en) 2014-01-29
WO2012129566A3 (en) 2012-12-06
US20140010431A1 (en) 2014-01-09
WO2012129566A2 (en) 2012-09-27

Similar Documents

Publication Publication Date Title
Stayman et al. Model-based tomographic reconstruction of objects containing known components
US9177374B2 (en) Model-based processing of image data
Gong et al. Iterative PET image reconstruction using convolutional neural network representation
EP3683771B1 (en) Medical processing apparatus
JP7381224B2 (ja) 医用画像処理装置及びプログラム
EP3716214B1 (en) Medical image processing apparatus and method for acquiring training images
Schmidt et al. A spectral CT method to directly estimate basis material maps from experimental photon-counting data
Zhang et al. Regularization strategies in statistical image reconstruction of low‐dose x‐ray CT: A review
Sidky et al. Constrained ${\rm T} p {\rm V} $ Minimization for Enhanced Exploitation of Gradient Sparsity: Application to CT Image Reconstruction
CN106920246B (zh) 在存在金属伪影的情况下用于分割的不确定性图
Tilley et al. Penalized-likelihood reconstruction with high-fidelity measurement models for high-resolution cone-beam imaging
KR20190103227A (ko) 단층촬영 재구성에서 사용하기 위한 데이터의 딥러닝 기반 추정
Abdurahman et al. Beam hardening correction using cone beam consistency conditions
WO2014185078A1 (ja) X線ct画像処理方法,x線ct画像処理プログラム及びx線ct画像装置
JP2008528168A (ja) X線投影の補正又は拡張を行う装置及び方法
US20160206272A1 (en) Computed tomography having motion compensation
Riblett et al. Data‐driven respiratory motion compensation for four‐dimensional cone‐beam computed tomography (4D‐CBCT) using groupwise deformable registration
JP7341879B2 (ja) 医用画像処理装置、x線コンピュータ断層撮影装置及びプログラム
Pourmorteza et al. Reconstruction of difference in sequential CT studies using penalized likelihood estimation
Wu et al. Using uncertainty in deep learning reconstruction for cone-beam CT of the brain
EP3881289A1 (en) Artificial intelligence (ai)-based standardized uptake value (suv) correction and variation assessment for positron emission tomography (pet)
EP2742483B1 (en) Image processing method
Van Eyndhoven et al. Region-based iterative reconstruction of structurally changing objects in CT
EP3123446B1 (en) Image generation apparatus
Li et al. 3D coronary artery reconstruction by 2D motion compensation based on mutual information