(57)【要約】  コンポーネントを含むオブジェクトの画像データを処理する撮像システムである。撮像システムは、画像データを取得する撮像デバイスとプロセッサとを含む。プロセッサは、撮像デバイスから画像データを受信して、コンポーネントのコンポーネントモデルを得て、撮像デバイスの撮像デバイスモデルを得て、コンポーネントモデルと撮像デバイスモデルとに基づいて、非制限的な対象物の関数を構築して、非制限的な対象物の関数と画像データとに基づいて、コンポーネントを含むオブジェクトのモデルを構築し、さらに、モデルに基づいて、コンポーネントを含むオブジェクトの画像を表示する表示デバイスも含まれる。【選択図】図１
2014509037.tif 000002
【特許請求の範囲】
【請求項１】
コンポーネントを含むオブジェクトの画像データを処理する撮像システムであって、
画像データを取得する撮像デバイスと、
前記撮像デバイスから前記画像データを受信して、前記コンポーネントのコンポーネントモデルを得て、前記撮像デバイスの撮像デバイスモデルを得て、前記コンポーネントモデルと前記撮像デバイスモデルとに基づいて、非制限的な対象物の関数を構築して、前記非制限的な対象物の関数と前記画像データとに基づいて、前記コンポーネントを含む前記オブジェクトのモデルを構築するプロセッサと、
前記モデルに基づいて、前記コンポーネントを含む前記オブジェクトの画像を表示する表示デバイスと
を備える撮像システム。
【請求項２】
前記プロセッサは、前記コンポーネントを含む前記オブジェクトの前記モデルの構築を、前記非制限的な対象物の関数において、登録パラメータ及び画像パラメータを変化させることに基づいて、前記画像データの前記非制限的な対象物の関数を解くことに基づいて行い、
前記コンポーネントを含む前記オブジェクトの前記モデルの構築は、解かれた前記非制限的な対象物の関数についての前記画像パラメータ及び前記登録パラメータに基づいて行われる、請求項１に記載の撮像システム。
【請求項３】
前記プロセッサは、前記非制限的な対象物の関数を解くことを、
前記画像パラメータを一定に保ちつつ、前記登録パラメータを変化させることに基づいて、前記画像データの前記非制限的な対象物の関数を最適化することと、
最適化された前記登録パラメータを一定に保ちつつ、前記画像パラメータを変化させることに基づいて、前記画像データの前記非制限的な対象物の関数を最適化することに基づいて行う、請求項２に記載の撮像システム。
【請求項４】
前記登録パラメータは、前記コンポーネントの位置及び向きを特定する、請求項２または３に記載の撮像システム。
【請求項５】
前記画像パラメータは、前記コンポーネントを含む前記オブジェクトのパラメータ化された特徴を特定する、請求項２から４のいずれか一項に記載の撮像システム。
【請求項６】
前記コンポーネントモデルは、前記コンポーネントのパラメータ化された特徴と登録された特徴とを定義する、請求項１から５のいずれか一項に記載の撮像システム。
【請求項７】
前記撮像デバイスモデルは、前記撮像デバイスの特徴に基づいて、前記コンポーネントを含む前記オブジェクトの前記モデルを構築するための特徴を定義する、請求項１から６のいずれか一項に記載の撮像システム。
【請求項８】
前記非制限的な対象物の関数の構築は、前記コンポーネントモデルと、前記コンポーネントを含む前記オブジェクト内の第２のコンポーネントの第２のコンポーネントモデル、及び、前記撮像デバイスモデルに基づいて行われる、請求項１から７のいずれか一項に記載の撮像システム。
【請求項９】
コンポーネントを含むオブジェクトの画像データを処理する方法であって、
撮像デバイスから前記画像データを受信する段階と、
前記コンポーネントのコンポーネントモデルを得る段階と、
前記撮像デバイスの撮像デバイスモデルを得る段階と、
前記コンポーネントモデルと前記撮像デバイスモデルとに基づいて、非制限的な対象物の関数を構築する段階と、
前記非制限的な対象物の関数と前記画像データとに基づいて、前記コンポーネントを含む前記オブジェクトのモデルを構築する段階と
を備える方法。
【請求項１０】
前記コンポーネントを含む前記オブジェクトの前記モデルを構築する段階は、
前記非制限的な対象物の関数において、登録パラメータ及び画像パラメータを変化させることに基づいて、前記画像データの前記非制限的な対象物の関数を解く段階を含み、
前記コンポーネントを含む前記オブジェクトの前記モデルを構築する段階は、
解かれた前記非制限的な対象物の関数についての前記画像パラメータ及び前記登録パラメータに基づいて行われる、請求項９に記載の方法。
【請求項１１】
前記非制限的な対象物の関数を解く段階は、
前記画像パラメータを一定に保ちつつ、前記登録パラメータを変化させることに基づいて、前記画像データの前記非制限的な対象物の関数を最適化する段階と、
最適化された前記登録パラメータを一定に保ちつつ、前記画像パラメータを変化させることに基づいて、前記画像データの前記非制限的な対象物の関数を最適化する段階と
を含む、請求項１０に記載の方法。
【請求項１２】
前記登録パラメータは、前記コンポーネントの位置及び向きを特定する、請求項１０または１１に記載の方法。
【請求項１３】
前記画像パラメータは、前記コンポーネントを含む前記オブジェクトのパラメータ化された特徴を特定する、請求項１０から１２のいずれか一項に記載の方法。
【請求項１４】
前記コンポーネントモデルは、前記コンポーネントのパラメータ化された特徴と登録された特徴とを定義する、請求項１０から１３のいずれか一項に記載の方法。
【請求項１５】
前記撮像デバイスモデルは、前記撮像デバイスの特徴に基づいて、前記コンポーネントを含む前記オブジェクトの前記モデルを構築するための特徴を定義する、請求項１０から１４のいずれか一項に記載の方法。
【請求項１６】
前記非制限的な対象物の関数の構築は、前記コンポーネントモデルと、前記コンポーネントを含む前記オブジェクト内の第２のコンポーネントの第２のコンポーネントモデル、及び、前記撮像デバイスモデルに基づいて行われる、請求項１０から１５のいずれか一項に記載の方法。
【請求項１７】
コンピューティングプラットフォームに実行されると、前記コンピューティングプラットフォームに、コンポーネントを含むオブジェクトの画像データを処理する方法を含む処理を実行させるプログラムであって、
前記方法は、
撮像デバイスから画像データを受信する段階と、
前記コンポーネントのコンポーネントモデルを得る段階と、
前記撮像デバイスの撮像デバイスモデルを得る段階と、
前記コンポーネントモデルと前記撮像デバイスモデルとに基づいて、非制限的な対象物の関数を構築する段階と、
前記非制限的な対象物の関数と前記画像データとに基づいて、前記コンポーネントを含む前記オブジェクトのモデルを構築する段階と
を含む、プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
［関連出願の相互参照］
本願は、２０１１年３月２４日提出の米国仮出願第６１／４６７，１８７号の優先権を主張しており、この内容全体をここに参照として組み込む。
【０００２】
本発明は、政府の援助を受けてなされている（アメリカ国立衛生研究所による許可番号ＣＡ−１２７４４４）。アメリカ政府も本発明に一定の権利を有している。
【０００３】
本発明は、画像データの処理に関しており、より詳しくは、画像データのモデルに基づく処理に関する。
【背景技術】
【０００４】
本明細書で言及される文献、特許出願公開公報、及び特許の全ての内容を、ここに参照として組み込むこととする。
【０００５】
断層撮影においては、多くの場合、画像全体の一部が先験的に分かっている。たとえば整形外科においては、脊柱手術の場合の脊椎頸ネジ及び脊椎頸用ロッド、人口関節のための膝及び臀部のインプラント、及び、外傷の場合の固定用プレート及びネジなどのコンポーネントについては、既に確かな知識がある。画像の補助を受けて行う手術においては、外科手術用ツールを画像視野内に配置する場合も多い。これらコンポーネントが金属製である場合には、該当元素が突出している測定値は、光子が不足することによって信号対雑音比が低くなることにより支障が出る場合がある。同様に、高度に減衰されているコンポーネントの再構成には、ゼロに近い突出した値の数学的反転を行う場合があるので、アルゴリズムがバイアスに対して敏感に反応しすぎるようになる場合がある（多エネルギー効果（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）。
【００１１】
より近年の研究では、ＰＩＣＣＳ等のアルゴリズムで、生体組織の事前スキャンを組み込んだ非常に具体的な画像の予備知識（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）。しかしこれら方法で撮像した場合であってもまだ撮像品質が低い。
【００１２】
したがって、画像データの処理を向上させる必要性がある。
【発明の概要】
【００１３】
本発明の一実施形態におけるコンポーネントを含むオブジェクトの画像データを処理する撮像システムは、画像データを取得する撮像デバイスとプロセッサとを含む。プロセッサは、撮像デバイスから画像データを受信して、コンポーネントのコンポーネントモデルを得て、撮像デバイスの撮像デバイスモデルを得て、コンポーネントモデルと撮像デバイスモデルとに基づいて、非制限的な対象物の関数を構築して、非制限的な対象物の関数と画像データとに基づいて、コンポーネントを含むオブジェクトのモデルを構築し、さらに、モデルに基づいて、コンポーネントを含むオブジェクトの画像を表示する表示デバイスも含まれる。
【００１４】
本発明の一実施形態におけるコンポーネントを含むオブジェクトの画像データを処理する方法は、コンポーネントを含むオブジェクトの画像データを処理する方法であって、  撮像デバイスから画像データを受信する段階と、コンポーネントのコンポーネントモデルを得る段階と、撮像デバイスの撮像デバイスモデルを得る段階と、コンポーネントモデルと撮像デバイスモデルとに基づいて、非制限的な対象物の関数を構築する段階と、非制限的な対象物の関数と画像データとに基づいて、コンポーネントを含むオブジェクトのモデルを構築する段階とを含む。
【００１５】
記載を考慮することでさらなる目的及び利点が明らかになる。
【図面の簡単な説明】
【００１６】
【図１】本発明の一実施形態におけるシステムのブロック図である。
【００１７】
【図２】本発明の一実施形態におけるコンポーネントを含むオブジェクトのモデルを構築するための処理フローチャートの一例を示す。
【００１８】
【図３】本発明の一実施形態における登録パラメータ及び画像パラメータに基づいてコンポーネントを含むオブジェクトのモデルを構築する処理フローチャートの一例を示す。
【００１９】
【図４】本発明の一実施形態における２次元のカーネルベースの補間の例を示す。
【００２０】
【図５】本発明の一実施形態におけるネジと背景の量とのモデルの例を示す。
【００２１】
【図６】本発明の一実施形態における、一連の推定された軸方向の断面（axial slice）及びインプラントの姿勢を例示する。
【００２２】
【図７】本発明の一実施形態における収束グラフを示す。
【００２３】
【図８】本発明の一実施形態における異なる推定値に基づくネジと背景の量を、軸方向画像及び冠状画像の間で比較した様子である。
【００２４】
【図９】本発明の一実施形態における減衰画像を例示する。
【００２５】
【図１０】本発明の一実施形態における、一方向の単一コンポーネントからなるネジの減衰画像を例示する。
【発明を実施するための形態】
【００２６】
本発明の一実施形態を以下に詳述する。実施形態を説明するにあたり、明瞭に記載するために一定の用語を利用する。しかし本発明は、選択されたこの特定の用語に限定されることを意図はしていない。当業者であれば、本発明の広義の構想を逸脱せずに、他の均等物であるコンポーネントを利用したり、他の方法を開発したりすることもできることを想到するだろう。本明細書で利用する参考文献は、全てをそれぞれ参照として組み込む。
【００２７】
図１は、本発明の一実施形態におけるシステム１００のブロック図である。システム１００は、撮像デバイス１０２と、撮像デバイスモデルソース１０４、コンポーネントモデルソース１０６、プロセッサ１０８、及び表示デバイス１１０を含んでよい。
【００２８】
撮像デバイス１０２は、オブジェクトの画像データを取得するデバイスである。画像データは、該オブジェクトのそれぞれ異なる断面画像に関するデータであってよい。異なる断面画像は、ｘ線をそれぞれ異なる焦点面に投影させることで得られたものであってよい。たとえば撮像デバイス１０２は、コンピュータ断層撮影（ＣＴ）マシンであってよく、撮像されるオブジェクトは、内部にコンポーネントを埋め込まれたヒトであってよく、画像データは、ＣＴマシンからの投影画像データであってよい。コンポーネントは、金属ネジ等の任意のオブジェクトであってよい。
【００２９】
撮像デバイスモデルソース１０４は、撮像デバイスモデルの提供元であってよい。撮像デバイスモデルは、コンポーネントを含むオブジェクトの画像データに基づいて、コンポーネントを含むオブジェクトのモデルを構築するために利用されてよい撮像デバイスのモデルであってよい。撮像デバイスモデルは、撮像デバイス１０２の特徴に基づいて、コンポーネントを含むオブジェクトのモデルを構築する特徴を定義するものであってよい。たとえば特徴によって、モデルが、撮像デバイス１０２が撮像する焦点面、撮像デバイス１０２が得る画像の雑音等を考慮に入れて構築されるべきである、と指定されてよい。撮像デバイスモデルソース１０４は、撮像デバイスモデルを格納するデバイスまたは撮像デバイスモデルを生成するデバイスであってよい。
【００３０】
コンポーネントモデルソース１０６は、コンポーネントモデルの提供元であってよい。コンポーネントモデルはコンポーネントのモデルであってよい。コンポーネントモデルは、コンポーネントのパラメータ化された特徴及び登録された特徴を定義するものであってよい。パラメータ化された特徴とは、該コンポーネントの減衰特徴を特定する特徴であってよい。たとえばパラメータ化された特徴には、減衰係数、ハウンスフィールド値、材料の組成、または材料の密度が含まれていてよい。パラメータ化された特徴は、コンポーネントの組成が異種であるか同種であるかを示してよい。登録された特徴は、コンポーネントの形状、及びコンポーネントが変形可能である場合、コンポーネントが変化することができる形状、許される位置、向き等を示す特徴であってよい。たとえばコンポーネントモデルは、コンポーネントが画像データに持つ影響をモデル化するために利用可能なコンポーネントの材料組成及び構造を特定するコンピュータ支援設計（ＣＡＤ）モデルであってよい。コンポーネントモデルソース１０６は、撮像デバイスモデルを格納するデバイスまたは撮像デバイスモデルを生成するデバイスであってよい。
【００３１】
プロセッサ１０８は、撮像デバイス１０２から画像データを受信して、撮像デバイスモデルを撮像デバイスモデルソース１０４から取得して、コンポーネントモデルソース１０６からコンポーネントモデルを取得する処理ユニットであってよい。処理ユニットは、コンピュータ等のコンピューティングデバイスであってよい。コンポーネントが画像データに関する画像の視野にあることがわかっている場合には、プロセッサ１０８は、画像データ、撮像デバイスモデル、及びコンポーネントモデルソースを利用して、コンポーネントを含むオブジェクトのモデルを構築してよい。たとえば、コンポーネントとしてネジが、撮像対象のヒトの内部にあることがわかっている場合には、コンポーネントを含むオブジェクトは、ヒト及びネジになり、プロセッサ１０８は、ネジを含むヒトのモデルを構築するために、ネジのコンポーネントモデルを利用することができる。
【００３２】
画像データ、撮像デバイスモデル、及びコンポーネントモデルソースに基づいてコンポーネントを含むオブジェクトのモデルを構築することで、画質が上がり、これによりデータの忠実さの低いデータが取得されてもよくなる。コンポーネントを含むオブジェクトのモデルを構築するにあたり、プロセッサ１０８は、さらに、コンポーネントの登録パラメータを計算してよい。登録パラメータは、コンポーネントの位置、向き、及び形状を特定するものであってよい。プロセッサ１０８は、さらに、コンポーネントを含むオブジェクトのモデルに基づいてコンポーネントを含むオブジェクトの画像を生成してもよい。
【００３３】
表示デバイス１１０は、プロセッサ１０８から生成された画像を受信して、生成された画像を表示するデバイスであってよい。たとえば表示デバイス１１０は、液晶表示（ＬＣＤ）デバイスであってよい。表示デバイス１１０は、さらに、プロセッサ１０８から登録パラメータを受信して、登録パラメータを表示してもよい。
【００３４】
撮像デバイス１０２、撮像デバイスモデルソース１０４、コンポーネントモデルソース１０６、プロセッサ１０８、及び表示デバイス１１０は、それぞれが別のデバイスであっても、１つの統合デバイスであっても、別のデバイス及び統合デバイスの組み合わせであってもよい。撮像デバイス１０２、撮像デバイスモデルソース１０４、コンポーネントモデルソース１０６、プロセッサ１０８、及び表示デバイス１１０は、ネットワークを介して（たとえばイントラネットまたはインターネット）または１以上の統合デバイス内の回路を介して通信できてよい。
【００３５】
図２は、本発明の一実施形態におけるコンポーネントを含むオブジェクトのモデルを構築するための処理フローチャート２００の一例を示す。最初に、プロセッサ１０８は、撮像デバイス１０２から、コンポーネントを含むオブジェクトの画像データを受信して（ブロック２０２）、コンポーネントモデルソース１０６からコンポーネントのコンポーネントモデルを取得して（ブロック２０４）、撮像デバイスモデルソース１０４から撮像デバイス１０２の撮像デバイスモデルを取得してよい（ブロック２０６）。
【００３６】
上述したように、プロセッサ１０８は、撮像デバイスモデルソース１０４及びコンポーネントモデルソース１０６のうち１以上とは別であっても、これらに統合されていてもよい。プロセッサ１０８が別の実体である場合には、プロセッサは、関連するソースとネットワーク経由で通信することにより、関連するモデルを取得してよい。プロセッサ１０８が統合されている場合には、プロセッサは、関連するソースとして機能して、モデルを生成することで、関連するモデルを得てよい。
【００３７】
コンポーネントモデルと撮像デバイスモデルとを利用することで、プロセッサ１０８は、コンポーネントモデルと撮像デバイスモデルとに基づいて、非制限的な対象物の関数（unconstrained objective function）を構築してよい（ブロック２０８）。非制限的な対象物の関数は、撮像デバイスモデル（システマチックな雑音を含む）とコンポーネントモデルとに基づいて、尤度関数から導出されてよい。コンポーネントモデルは、マスキング処理と加算とを利用してオブジェクトに統合されることで、非制限的な対象物の関数を生成してよい。特に、対象物の関数は、オブジェクトの背景を特定する画像パラメータが、最適化処理において任意の値を自由にとることができるために、非制限的と称されてよい。
【００３８】
プロセッサはさらに、画像データに関連付けられている１以上の画像の視野にあることがわかっている１以上の追加のコンポーネントのための１以上の追加のコンポーネントモデルに基づいて、非制限的な対象物の関数を構築してよい。
【００３９】
非制限的な対象物の関数を利用して、プロセッサ１０８は、非制限的な対象物の関数及び画像データに基づいて、コンポーネントを含むオブジェクトのモデルを構築してよい（ブロック２１０）。
【００４０】
図３は、本発明の一実施形態における登録パラメータ及び画像パラメータに基づいてコンポーネントを含むオブジェクトのモデルを構築する処理フローチャート３００の一例を示す。図３は、プロセッサ１０８が、どのように、非制限的な対象物の関数及び画像データに基づいてコンポーネントを含むオブジェクトのモデルを構築することができるかを示すより詳しい実施形態を示していてよい（ブロック２１０）。
【００４１】
プロセッサ１０８は、画像パラメータを初期値に設定し、登録パラメータを初期値に設定してよい（ブロック３０２）。たとえば登録パラメータの初期値は、コンポーネントを、オブジェクトの最長の軸に合わせられたオブジェクトの中心に配置するものであってよい。画像パラメータは、コンポーネントを含むオブジェクトのパラメータ化された特徴を特定するパラメータであってよい。たとえば、画像パラメータは、コンポーネントを含むオブジェクトの背景の減衰を特定してよい。登録パラメータは、コンポーネントの位置及び向きを特定してよい。
【００４２】
初期値を利用して、画像パラメータは一定に保ちつつ、登録パラメータを変化させることで、非制限的な対象物の関数を画像データ向けに最適化することができる（ブロック３０４）。たとえばコンポーネントの推定位置を移動させてよい。非制限的な対象物の関数を増やすと、新たな位置を維持することができる。非制限的な対象物の関数を減らすと、コンポーネントの推定位置を前の位置に戻すことができる。
【００４３】
初期値及び更新された登録パラメータを利用することで、最適化された登録パラメータを一定に保ちつつ画像パラメータを変化させることで、非制限的な対象物の関数を画像データに最適化することができる（ブロック３０６）。たとえば、画像の背景の減衰は、画像パラメータを変更させることで推定されてよい。
【００４４】
更新登録パラメータと画像パラメータとの間の最適化プロセスは、所望の回数繰り返し実行することができる。
【００４５】
最終的な登録パラメータ及び画像パラメータを利用することで、プロセッサ１０８は、コンポーネントを含むオブジェクトのモデルを構築することができる（ブロック３０８）。たとえばプロセッサ１０８は、最終的な登録パラメータを利用して、コンポーネントの位置及び向きを決定することができ、最終的な画像パラメータを利用して、画像の各部分からの減衰を決定することができる。
【００４６】
＜例＞
製造されたコンポーネント（外科ツール、インプラント等）が断層の視野内にある確率は、増え続けている。この理由の１つは、老齢化が進み、人口装具の普及が進んだために、診断的な撮像をうける人が、インプラント（特に臀部及び膝のインプラント）を埋め込まれている確率が高くなっていることにある。別の理由としては、外科手術の補助手段として、手術中に撮像を行うケースが増えており（コーンビームＣＴ）、手術道具及びデバイス（外科用ネジ及びプレート）が対象の生体組織内またはその付近に配置される確率が上がっているからである。これらコンポーネントが金属を含んでいる場合、再構成された量は、コンポーネントの付近またはコンポーネントから離れた組織の画質に悪影響を与える、激しいアーチファクトを含む可能性がある。これらコンポーネントの物理的モデルが存在していることから、この知識を再構成アルゴリズムに統合して、これらアーチファクトを低減させる、またとないチャンスである。
【００４７】
以下に、コンポーネントの配置及び組成に関する既知の情報を明示的に組み込む、モデルに基づくペナルティ付き尤度推定法を説明する。この方法は、既知のコンポーネントそれぞれの生体組織ならびに位置及び姿勢を共同推定する代替最大化法（alternating maximization method）を利用する。提案する方法は、実質的に光子が不足していても、シミュレーションされた脊椎の脊椎頸ネジの再構成において金属のインプラントの境界付近であっても、アーチファクトのほとんどない画像を生成することができることが証明されている。デバイスの姿勢を同時に推定することができるので、品質保証及び処置の検証に有用でありうるデバイスの定量的情報を提供することができる。
【００４８】
手術道具、固定ハードウェア、インプラントの物理的モデルが利用可能なことから、再構成ルーチンに非常に具体的な予備知識を組み込むことができ、さらなる利点を生じる可能性がある。これらコンポーネントは、製造されたオブジェクトを表していることから、これらの材料組成及び構造を完全特定できるＣＡＤモデルが利用可能になる。下記で、このような既知の物理的モデルを再構成プロセスに統合するアルゴリズムを説明する。具体的には、オブジェクトの量のモデル自身は、１）（未知の）背景の生体組織の量、及び、２）撮像視野にあることがわかっているコンポーネント（１または複数）の組み合わせである。
【００４９】
コンポーネントの形状及び減衰の分布は分かっているものの（たとえば、デバイスの形状及び材料の内容を特定するＣＡＤモデルから導出される）、位置及び姿勢は分かっていない。したがって、各コンポーネント登録に関するパラメータも、オブジェクトのモデルの一部である。結果生成される再構成スキームは、画像パラメータと登録パラメータとの両方の関数である対象物を有することになり、これら２つのパラメータのセットは共同して推定される。この方法は、共同で画像再構成及び登録を行おうとする他の対象物の関数といくらか類似している点がある（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）。しかし以下に要点をまとめたこの方法は、とりわけ、非制限的な対象物の関数を利用しており、姿勢がわからない任意の数の既知のコンポーネントについて生成されている。
【００５１】
この方法は以下のようにまとめられる。セクションＩで、尤度ベースの対象物の関数と、画像の量を、未知の背景の生体組織内の任意の数の既知のコンポーネント（姿勢は分かっていない）の組み合わせとしてモデル化する再構成アルゴリズムとを記載している。セクションＩＩでは、この方法の実行を説明しており、従来の分析及び繰り返し法比較している。
【００５２】
＜Ａ．フォワード・モデル＞
透過型断層撮影システムの以下の計測モデルを例にとる。
【数１】
2014509037.tif 000003
【００５３】
平均投影計測値
【数Ａ】
2014509037.tif 000004
は、ベールの法則により、オブジェクトの線積分
【数Ｂ】
2014509037.tif 000005
に関連している。各計測値は、この光線の経路の減衰されていないｘ線のフルエンスと、この検出器部材の利得とに関連しているスカラー値
【数Ｃ】
2014509037.tif 000006
を有している。離散化されたオブジェクト（たとえばボクセルベースを利用して）については、線積分を以下に表すことができる。
【数２】
2014509037.tif 000007
【００５４】
ここでaijは、ｊ番目のボクセル（または）他のベース）μｊのｉ番目の線積分に対する貢献を示している。このモデルは、以下のように、行列ベクトル形式にコンパクトに書き換えることができる。
【数３】
2014509037.tif 000008
【００５５】
ここでＤ｛・｝は、ベクトル引数からの対角行列を形成する演算子を表しており、システム行列Ａは、aij全ての集合（collection）を表している。通常、数３は、オブジェクトと計測値との間の完全な関係を表しており、この数３から、対象物の関数を導出することができる。さらに、オブジェクトを、減衰分布が既知である任意の数のコンポーネントμＩ（ｎ）と、未知の生体組織の背景量μ＊との組み合わせとしてパラメータ化するよう選択することができる。
【数４】
2014509037.tif 000009
ここでＷ（λ）は、ベクトルがパラメータ化する変換演算子を表している。この変換は一般的なものであるが、ここでは、３次元（３Ｄ）の剛体変換（rigid transformation）に着目している。この場合には、ベクトルλは、コンポーネントの向きを定義する６つの要素（変換及び回転値）からなる。数４の第２項は、Ｎ個の既知の減衰量μＩ（ｎ）の減衰値（それぞれが、そのパラメータベクトルλ（ｎ）に従って（剛体）変換されている）の合計を表している。数４の第１項は、各既知のコンポーネントに対応している、１組の変換されたマスクｓ（ｎ）と乗算して修正された背景の生体組織の貢献μ＊を表している。具体的には、各マスクは、各コンポーネントのサポート領域を表している。典型的なマスクは、端部における部分容積効果を除いて、ほとんどがバイナリである。数４及び本明細書全体で利用されている積の演算子Пは、行列またはベクトルのオペランドの、要素ごとの積を表す。
【００５６】
＜Ｂ．尤度ベースの対象物の関数＞  数３及び数４は、平均計測値とパラメータ化された量との間の関係を表している。雑音モデルを選択することで、尤度関数が導き出される。ポアソン・モデルを呼び出して、以下の対数尤度を導出することができる。
【数５】
2014509037.tif 000010
【００５７】
ここで、Ｎ個のコンポーネント全ての登録パラメータのセットの集合は、以下のように表される。
【数６】
2014509037.tif 000011
【００５８】
以下の暗示的な形式で、ペナルティ付きの尤度推定量（penalized-likelihood estimator）を利用するよう選択される。
【数７】
2014509037.tif 000012
【００５９】
ここでＲ（・）は、非常に雑音の多い画像を抑制するペナルティ付き関数を示す。数７における暗示的な推定量は、「既知のコンポーネント再構成」（ＫＣＲ）法と称すことができる。この推定量は、生体組織の背景、及び、各既知のコンポーネントに関する登録パラメータのセットの両方を共同して近似する。正則化された項は一般的であるが、対の二次ペナルティ（pairwise quadratic penalty）を利用してよい。
【数８】
2014509037.tif 000013
【００６０】
ここで、Ｎｊは、ボクセルｊの周りの一次的周辺ボクセル（first-order neighborhood of voxels）を示す。
【００６１】
＜Ｃ．変換演算子＞
数７の共同推定量は、数４の変換演算子Ｗを特定する必要がある。数４の公式化は、一般的であり、変換は、ワーピング演算及び他の形態のパラメータ化された変形を表していてよい。同様に、Ｗは、旋回軸または中心点等の拘束された動きを持つコンポーネントの変換を表していてよい。この例は、剛体変換に着目している。
【００６２】
「移動」画像から「剛体」画像への登録に関する変換演算子は、通常、１）移動画像の格子化された部分を、変換画像のこれらの点の位置へマッピングすること、２）移動画像の各点の周りのボクセルの周辺部分を利用して、変換画像のボクセル値を生成すること、という２つの部分に分解することができる。３Ｄの剛体変換において、移動画像を変換画像への点のマッピングは、以下のように示されてよい。
【数９】
2014509037.tif 000014
【００６３】
［ｘｊ  ｙｊ  ｚｊ］Ｔは異動画像のボクセル位置を表しており、［ｘ'ｊ  ｙ'ｊ  ｚ'ｊ］Ｔは変換画像の位置である。ベクトル［ｘｃ  ｙｃ  ｚｃ］Ｔは、画像の量の中心の位置であり、［ｘＴ  ｙＴ  ｚＴ］Ｔは、変換パラメータである。各軸の周りの回転行列は、以下のように定義される。
【数１０】
2014509037.tif 000015
【００６４】
θ、Ψ、及びφが、各軸周りの回転パラメータを示す。したがって変換のための６つのエレメントのパラメータベクトルは、λ＝［ｘＴ  ｙＴ  ｚＴ  θ  Ψ  φ］Ｔと明示的に定義することができる。数９は、［ｘｊ  ｙｊ  ｚｊ］Ｔから［ｘ'ｊ  ｙ'ｊ  ｚ'ｊ］Ｔへのｊ番目のボクセル位置のマッピングを示しているが、マッピング全体は、３つの「移動画像」ベクトルｘ、ｙ、及びｚと、移動画像及び変換画像のすべての格子点の集合である３つの変換画像ｘ'、ｙ'、及びｚ'で表される。
【００６５】
点のマッピングが定義されている場合には、補間段階も定義されている。各次元で分離可能なカーネルベースの補間法を考慮することができる。
【００６６】
図４は、本発明の一実施形態における２次元のカーネルベースの補間の例を示す。左上の図は、移動画像（ｖ）の変換画像（ｕ）へのマッピングを示す。右上の図は、移動画像の値の周辺部分に基づいて、１つの変換画像の点を計算することを示している。下の図は変換画像の位置の関数であるカーネルが、各次元沿いに連続して適用されて、この点における変換画像の値を生成することを示している。
【００６７】
変換画像の各位置ｊにおいて、移動画像の各点の周辺部分にカーネルを適用して、変換画像の値を近似する。一次元（１Ｄ）カーネルは、周辺領域の各次元沿いに連続して適用され、合計されることで、各パスが、１つ次元を下げる。たとえば、２Ｄの場合には、２つの１Ｄのカーネルを適用して、１つのボクセル値を得る。３Ｄの場合には、この演算は、元の変換に戻ることを意味しており、数学的には以下のようになる。
【数１１】
2014509037.tif 000016
【００６８】
ｕ及びｖは、それぞれ変換画像及び移動画像を示す。Ｋ個の演算子は、全ての位置において各方向に１Ｄのカーネルを利用することを表しており、この点のマッピングの関数である。多くのカーネルの選択があるが、以下のように定義されるＢスプラインカーネル（P. Thevenaz, et al., "Interpolation revisited," IEEE Trans Med Imaging, vol. 19, pp. 739-58, Jul 2000）を利用する。
【数１２】
2014509037.tif 000017
【００６９】
このカーネルは、１）変換画像の負ではない性質を維持して（負ではない「移動」画像）、負の減衰係数が推定値に入らないようにする、２）ｋ（ｔ）が、微分可能であり、勾配ベースの方法を数７で提案する目的の最適化で利用することができる。
【００７０】
以下のセクションでは、分析勾配情報（analytic gradient information）を利用して、数７の解を近似する。勾配情報は、必ず、推定するパラメータの変換演算子の導関数を含むので、これら等式はここから導出される。具体的には、λの要素に対するＷの演算子の導関数を導出する。連鎖律を利用して、以下のように記述することができる。
【数１３】
2014509037.tif 000018
【００７１】
ここで、λのｍ番目の導出関数である変換演算子は、
【数Ｄ】
2014509037.tif 000019
と表す。
【数Ｅ】
2014509037.tif 000020
演算子は、導関数カーネル（derivative kernel）を適用する以外は、元のカーネル演算子に類似している。数１２のＢスプラインカーネルについては、導関数カーネルは以下のように計算されてよい。
【数１４】
2014509037.tif 000021
【００７２】
数１３の残りの偏導関数（remaining partial derivatives）は、数９をλｍで微分して計算することができる。具体的には以下のように計算することができる。
【数１５】
2014509037.tif 000022
【数１６】
2014509037.tif 000023
【００７３】
【数Ｆ】
2014509037.tif 000024
は、数１０に示す回転行列の要素ごとの導関数を表している。
【００７４】
＜Ｄ．最適化法＞
前のサブセクションでは、フォワード・モデルと、数７に示す推定値に含まれるすべての要素を説明した。この推定値は、対象物の関数の暗示的なマキシマイザ（maximizer）であるので、最適化アルゴリズムは、解を繰り返し近似するよう適合されていてよい。前のセクションの展開から、対象物の関数は微分可能であり、推定パラメータ全てに対する導関数の分析的表現が可能であることがわかるだろう。
【数１７】
2014509037.tif 000025
【数１８】
2014509037.tif 000026
【００７５】
対象物の関数の分析勾配を前提にすると、パラメータ
【数Ｇ】
2014509037.tif 000027
について共同で解くためには、いずれの勾配ベースのアルゴリズムであっても利用できると素直に考えることができるが、これは、断層撮像再構成専用のアルゴリズムの長年に及ぶ開発の歴史を無視することになる。したがって、画像の量のパラメータμ＊が、断層に特化した画像更新を利用して更新され、登録パラメータΛが従来の勾配ベースの方法によって更新される、代替最小化方法（alternating minimization approach）を利用することを提案する。
【００７６】
Λは、６Ｎ個の要素しか含まず比較的小さいので（剛体変換されたコンポーネント）、ＢＦＧＳ（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）。（Ｎが大きい場合には、限定されたメモリの変形を利用すればいいだろう（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）を利用するとよい。Ｅｒｄｏｇａｎ及びＦｅｓｓｌｅｒが考えた導出法（H. Erdogan and J. A. Fessler, "Monotonic algorithms for transmission tomography," IEEE Trans Med Imaging, vol. 18, pp. 801-14, Sep 1999）では、数５を、限界対数尤度の合計として、以下のような表している。
【数１９】
2014509037.tif 000028
【００７８】
以下の定義を行うことができる。
【数２０】
2014509037.tif 000029
【００７９】
Ａが負ではない要素をもつ投影演算子を表しており、Ｗが、負ではない性質を維持するＢスプライン近似を利用して選択されており、ｓ（ｎ）、μ＊、及びμｉ（ｎ）も負ではないので、ｔｉ≧０及びｄｉ≧０である点に留意されたい（行列Ｂは、修正されたシステム行列の一種を表している）。ｔｉ≧０及びｌｉ＝ｔｉ＋ｄｉ≧０であるので、ｇｉ（ｔｉ）が、（H. Erdogan and J. A. Fessler, "Monotonic algorithms for transmission tomography," IEEE Trans Med Imaging, vol. 18, pp. 801-14, Sep 1999）の理論１の条件を満たしており、最適な曲率をもつｇｉ（ｔｉ）の放物線サロゲート法を以下のように書き換えることができる。
【数２１】
2014509037.tif 000030
【００８０】
ここで、上付き文字［ｋ］は、ｋ番目の繰り返しの値を示し、ｃｉは、放物線サロゲート法の曲率を示しており、［・］+は、負ではないようにする関数を示している（ゼロで切り捨てる）。ｇの第１及び第２の導関数はそれぞれ
【数Ｈ】
2014509037.tif 000031
を表し、数１９及び数２０から直接計算される。数２１のサロゲート法を利用して、対数尤度のサロゲート法全体を形成することができるが、（H. Erdogan and J. A. Fessler, "Ordered subsets algorithms for transmission tomography," Phys Med Biol, vol. 44, pp. 2 -51 , Nov 1999）でも利用されたＤｅ  Ｐｉｅｒｒｏのコンベクシティ・トリック（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）を利用する。
【数２２】
2014509037.tif 000032
【００８１】
分離可能なサロゲートを得るためには、
【数２３】
2014509037.tif 000033
とする。
【００８２】
【数Ｉ】
2014509037.tif 000034
を選択することができる。したがって正則化されていない対象物（unregularized objective）については、以下の更新を行うことができる。
【数２４】
2014509037.tif 000035
【００８３】
これら更新を（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）を適用して収束率を向上させ、予め計算された曲率を利用することで、計算時の負荷を低減させる。
【００８４】
＜Ｅ．実装＞
上述のようにして得た方法を調査するべく、全てのルーチンを呼び出しとともにＭａｔｌａｂ（マサチューセッツ州のＮａｔｉｃｋのＭａｔｈｗｏｒｋｓ，Ｉｎｃ社）に実装して、ＫＣＲ方法のより計算の多い側面の外部ライブラリを作成した（custom）。具体的には、マッチングされたＳｉｄｄｏｎ（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）プロジェクタ及びバックプロジェクタを、ＣＵＤＡラインブラリを利用してＣ／Ｃ＋＋に実装して、ＧＰＵで実行できるようにした。同じ方法は、ＧＰＵによる３Ｄカーネルベースの補間についても行うことができる。
【００８５】
表１に、ＫＣＲの代替最大化法全体の擬似コードを示す。コードは、登録パラメータ更新を含む１つ目のものと、順序付けられたサブセットベースの量の更新を含む２つ目のものとからなる、２つの分離したブロックからなる。この擬似コードでは、一般化された繰り返しの上付き文字[ｋ、ｐ、ｍ]を利用して、主要な繰り返し数ｋ、登録ブロックのサブ繰り返し数ｐ、及び、画像更新ブロックサブ繰り返し数ｍを示す。シンボル[・]は、具体的なインデックスに依存しない値に利用される。
【００８６】
表１は、アルゴリズムの幾つかの自由なパラメータ選択を示す。つまり、登録サブ繰り返し数Ｐ、順序付けられたサブセット数Ｍ、及び、繰り返し総数を示す。スケーリング係数σは、最初の逆ヘシアン推定値Ｈ[０，０，・]の類推に利用され、対象物の関数のスケールに基づいて、または、適切なステップサイズについて大まかなスケールを見つけるという初期ライン検索を利用して、経験的に選択する必要がある。この最適化の他の目立った特徴には、Ｈを主要なループの繰り返しの間で再利用することと、二次の概算ペナルティ（
【数Ｊ】
2014509037.tif 000036
でそれぞれ表される）の一次及び二次導関数を利用して画像更新に正規化を含めることとが含まれる。
【表１】
2014509037.tif 000037
表１：ＫＣＲアルゴリズムが利用する代替最大化の概略
【００８７】
Ｍ＝１の場合、画像更新は、局所的最大値に単調に収束するはずであり、Ｗｏｌｆｅ条件を満たすライン検索を利用する場合、登録更新は、局所的最大値に単調に収束するはずである。したがってこれらの条件下では、アルゴリズム全体が局所的最大値に収束することも予期される。通常の順序付けをされたサブセット法が厳密には収束しないのと同様に、Ｍ＞１の場合には真の収束は予期されない。他のアルゴリズムの置換も可能である（主要ループ外の曲率を予め計算して、順序付けされた個々のサブセットの画像更新間で登録更新をミキシングすることを含む）。
【００８８】
＜セクションＩＩ：結果＞
＜Ａ．シミュレーション研究＞
ＫＣＲ法のパフォーマンスを調査するために、シミュレーションされたコーンビーム計算断層撮影（ＣＢＣＴ）システムを利用して複数の研究を行った。この作業は、脊柱手術における脊椎頸ネジ配置をやりやすくするための撮像に着目して行われた（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）。
【００８９】
図５は、本発明の一実施形態におけるネジと背景の量とのモデルの例を示す。（Ａ）は、多軸の脊椎頸ネジのＣＡＤモデルを示す。モデルは、旋回軸の頭部及びネジ山を設けられたネジの２つの主要部品からなる。生体組織の背景の量については、高品質の従来のＣＴスキャンをカスタムのＲＡＮＤＯファントム（RANDO phantom：Ｐｈａｎｔｏｍ  Ｌａｂｏｒａｔｏｒｙ，Ｉｎｃ．社製）の腰静脈領域に行うことで得られるデジタルのファントム（digital phantom）を利用することができる。このファントムの多面断層（multiplanar slices）を図５に示す。（Ｂ）は、軸方向の断層を示し、（Ｃ）は矢状方向の断層を示し、（Ｄ）は、冠状の断層を示す。
【００９０】
すべての調査について、３００ｘ３００ｘ１２８ボクセルの画像量（１立方ミリメートル）が考慮された。脊椎頸ネジも同様にボクセル化され、（４）のオブジェクトモデルを利用して画像の量に配置された。（４）ではさらに、ＣＡＤモデルから導出されたボクセル化されたコンポーネントマスクｓ（ｎ）が必要となる。脊椎頸ネジは、０．３ｍｍ−１の線形減衰係数の同質チタンであると仮定する。比較のために、デジタルファントムの柔らかい組織の背景が０．０２ｍｍ−１であった。システムの形状は、ソースから検出器までの距離が１２００ｍｍであり、ソースから軸までの距離が６００ｍｍであるフラットパネルベースのＣ−アームシステムをエミュレートする目的で選択された（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）。プロジェクションは、（１．５５２ｘ１．５５２）ｍｍ２で、３６０角度を利用して３６０度においてシミュレーションされ、３６０ｘ１５０画素であった。シミュレーションされたデータは、単一エネルギーのフォワード・モデル及びポアソン雑音を利用して生成された。露光レベルは、減衰されていないビームにつき検出素子１つについて、均一的に１０４個の光子に固定された。以下のセクションで示されるように、この露光レベルでは、従来の方法を利用した場合、光子不足によるアーチファクトが顕著に出る。
【００９１】
３つの別個のネジ配置シナリオが調査された。１つ目が、一方向のモノリシック（単一コンポーネントの）脊椎頸ネジ、２つ目が、双方向のモノリシックネジ、３つ目が、一方向の２コンポーネントの多軸ネジである。登録パラメータΛを、腰椎への１つにネジ（１または複数）の適切な配置をシミュレーションするために選択した。
【００９２】
＜Ｂ．収束特性＞
ＫＣＲ法の収束特性を調査するために、上述した一方向の単一コンポーネントのネジ配置シナリオからシミュレーションしたデータを再構成した。背景の生体組織の最初の類推は、０．０３ｍｍ−１の上部減衰で切り捨てられたフィルタリングされた逆投影法（FBP）画像について行われた（ネジは背景の量μ＊には存在すべきではないので、生体組織から脊椎頸ネジを取り除くための大まかな試みにおいて）。登録パラメータについては初めの類推が複数回行われた。解決法から十分遠い複数の初めの類推において、ネジの位置の推定が局所的な最大値で失われた「キャプチャー問題」が観察された。（もっともよく行われるのは、再構成された視野を完全に出て、問題が従来のペナルティ月尤度推定にまで減るようにする）。しかし、投影データの真の位置にあるコンポーネントが少なくとも大体重なっている初期化において、真の姿勢に向かって方法が収束した（雑音による誤りを除いて）。キャプチャー範囲内の最初の類推については、姿勢の推定における残りの誤りが、１ボクセル内に、０．５度（回転パラメータにおいて）未満で１０４個の光子の雑音レベルにおいて一貫して見られた。
【００９３】
図６は、本発明の一実施形態における、一連の推定された軸方向の断面及びインプラントの姿勢を例示する。本図は、それぞれ登録及び画像サブ繰り返しについてＰ＝１０及びＭ＝６０を利用して、最初の類推からの始めたサンプルの繰り返しを示しており、最初の４回の繰り返し（Iteration）を示している。最初の類推においては、コンポーネントを任意に配置してよい。アルゴリズムの共同的な性質が明らかであり、外部ループの繰り返しの間に背景画像の縞状アーチファクトの低減とコンポーネントの姿勢両方に改善が見られた。図面に示すように、登録の更新と画像の更新とは、繰り返しを連続する間に同時に実行されてよい。
【００９４】
図７は、本発明の一実施形態における収束グラフを示す。図７の（ａ）は、６０個のサブセットを利用した画像更新と、４つの▽、６つの○、及び１０個の□のＢＦＧＳの繰り返しを利用した登録の更新とを交互に適用した、ＫＣＲの収束プロットを示している。図７の（ｂ）は、１０ＢＦＧＳの繰り返しと、１つの▽、５つの○、１０個の□、２０個の◇、４０個の△、及び６０個の☆の順序付けされたサブセットベースの画像の更新とを交互に繰り返し、且つ、動的な数のサブセット
【数Ｋ】
2014509037.tif 000038
を利用する方法を利用した場合の収束プロットを示している。対象物の関数の差は、登録ブロック（特定シンボル）と画像更新（・）の後にプロットされている。
【００９５】
複数のサブ繰り返しを様々に選択した場合の相対効果をよりよく理解するために、図７の収束プロットを計算した。２つのセットのプロットは、固定数の順序付けされたサブセットの更新に応じて、登録更新の数を変化させた場合の収束レートを表している。具体的には、ｋ番目の繰り返しと、解Φ＊（後者の値は、Ｐ＝０の場合に、動的な数のサブセットを利用して、ＫＣＲを１００回繰り返して近似されているが、これに関しては後述する）における対象物の関数の差をプロットしている。
【００９６】
図７の（ａ）は、６０のサブセットで４、６、及び１０の登録更新を利用して３つの方法を行った最初の２５回の繰り返しを示す場合を示す。対象物の関数の差は、登録更新（整数ｋの値同士の中間）、及び、画像更新（シンボル・で表されている）の両方の後に示されている。収束レートは、さらなるＢＦＧＳの更新がある度に増加しているように見えるが、Ｐ＝１０を超える数になると、実質的に同じ収束プロットになり、登録更新に重きをおいても得られるものが少ないことが示唆されている。さらには、画像更新を行うと、登録更新よりも対象物の関数の増加が大きくなることもわかった。
【００９７】
図７の（ｂ）では、登録更新数がＰ＝１０に固定されており、Ｍを１、５、１０、２０、４０、及び６０サブセットの間で変化させた。これらのケースにおいて、収束レートは、サブセット数に対する依存が強かった（パフォーマンスは従来の順序付けされたサブセット法と類似していた）。すべての方法が対象物を単調に増加させたかにおもえたが、Ｍ＝６０を超えてから平坦部が見られる（順序付けされたサブセットは厳密には収束しないために予期されることではある）。この平坦部を避けて、良好な収束レートを維持するために、サブセット数を動的に変化させる方法を実装する。具体的には、１２０個のサブセットから始めて、シーケンス９０、６０、４０、３０、２０、１０、３、２、及び１について、４回外部繰り返しを行うたびにサブセット数を低減させる。この方法は図７の（ｂ）の一番下のプロットに表されている。動的なサブセット法は最終的には１つのサブセットになるので、この方法では真の収束が得られることが予期される。サブセット数を変更した場合、４回繰り返すごとに関数値が段階的に向上する点に留意されたい。
【００９８】
＜Ｃ．画質の実験＞
画質の実験について、３つの脊椎頸ネジンプラントのシナリオと３つの異なる再構成法の全てを考慮した。具体的には、Ｐ＝１０にしてＫＣＲ法を２００回繰り返して、動的なサブセット法を適用した（２０回繰り返すごとに上述したシーケンスでＭが低減する）。これに比較して、ＦＢＰ法でも、従来のペナルティ付き尤度（ＰＬ）推定値を利用してデータを再構成した。ＦＢＰ処理では、データは、対数変換の前にゼロより少し高いデータを閾値として前処理された（ゼロのカウント計測値をなくすために）。（メジアンフィルタを利用してこれらゼロの値に補間を行う、ここに示されない代替的な方法を利用した場合であっても、定性的にほぼ同じ結果が生じた）。
【００９９】
ＰＬ法及びＫＣＲの両方で、正則化パラメータ（β）の同じ値の二次ペナルティが利用された。対象物同士は非常に類似しているので、パラメータが整合した場合には、これら２つの方法によって同様の解像度特性の画像が生成される。ペナルティ尤度法（ＫＣＲを含む）の解像度特性は、一般的に空間変形例であり（space-variant）、オブジェクトに依存している特性によって、正確な解像度のマッチングは難しくなる。βは、ＦＢＰ画像同士の空間解像度を定性的に整合させるように選択される。
【０１００】
図８は、本発明の一実施形態における異なる推定値に基づくネジと背景の量を、軸方向画像及び冠状画像の間で比較した様子である。図８は、ＫＣＲの画質を、ＦＢＰ及び従来のＰＬ推定の画質と比較したものである。各方法及び各脊椎頸ネジンプラントについて軸方向画像及び冠状画像が提示されている。一番上は、一方向の単一コンポーネントのネジであり、中間が、二方向の単一コンポーネントのネジであり、一番下が、一方向の２コンポーネントの多軸ネジである。すべての画像の窓及びレベルは、それぞれ５０及び１５０ＨＵである。真の画像（左）のオーバレイは、コンポーネントの真の姿勢を示しており、ＫＣＲ再構成画像（右）のものは、登録／再構成同時プロセスにより推定される姿勢を示している。
【０１０１】
図８では、結果は、左の列に示す新のデジタルファントム量とともに示されている。前述と同様に、コンポーネントの位置もオーバレイを利用して示されている。このオーバレイは、真の量の画像及びＫＣＲ画像に行われている（しかし、ＦＢＰまたは従来のＰＬでは、コンポーネント登録モデルを含まず、画像モデルしか含まないために、行われない）。
【０１０２】
図９は、本発明の一実施形態における減衰画像を例示する。単一コンポーネントの一方向のネジの場合の減衰全体（μ）の図では、（左）が真のファントムであり、（右）がＫＣＲである。ＫＣＲは、背景の減衰（μ＊）と登録パラメータΛの両方を推定するので、結果は従来の減衰画像としても（エラー！、参照元が見つかりません（右））、オーバレイとしても（エラー！、参照元が見つかりません８（右））、示すことができ、表示の好み及び動的な範囲に合致するほうを選択することができる。すべての画像の窓及びレベルは、それぞれ８００及び３００ＨＵである。
【０１０３】
そのほうが望ましければ、いつでも従来の減衰画像を、（４）オーバレイのないＫＣＲを選択することで示すことができる（たとえば、図９に示す単一コンポーネントの一方向のネジの場合）。各再構成技術及び各ネジ配置シナリオについて軸方向の画像及び冠状画像が示されている。冠状の断層は、茎までのネジの配置を示すために選択されている。この表示は、ネジが皮質骨までに達しているかを評価する際の診断においては特に重要である（横方向に侵害している場合（次善のネジの配置）、または医療的な侵害の場合（脊柱管に侵入している場合））。
【０１０４】
光子の露光レベルが１０４の場合、ネジの減衰が高いことによる光子の不足によるアーチファクトは、ＦＢＰ再構成において非常に顕著である。これらのアーチファクトは、軸方向の画像の周辺領域の中低度のオントラストの詳細を極度に失わせ、双方向のネジの配置シナリオにおいて最悪の影響がでる（このシナリオでは、金属量が増加することで光子の不足量が最大になる）。これと比較して、ＰＬ画像では大きな向上がみられ、縞状のアーチファクトの殆どが消えている。しかし、まだ特にネジの境界付近では残余アーチファクトが顕著であり、これにより椎骨の多くがあいまいであり、冠状画像に冠状の侵入があっても評価し難い。ＫＣＲ法は、殆どアーチファクトのない画像を生成している。残余アーチファクトはいくらかあるものの（そのほとんどはネジのヘッドあたりにある）、アーチファクトのコントラストは、ＦＢＰ及びＰＬ両方と比較して、非常に低い。さらに、冠状画像からは、ネジの配置において冠状の侵入がないことがよくわかり、真の量の画像によく合致している。ＫＣＲ法の画質は、ネジの境界付近であっても高く維持されている。
【０１０５】
図１０は、本発明の一実施形態における、一方向の単一コンポーネントからなるネジの減衰画像を例示する。ＫＣＲの正則化パラメータを変化させた効果が示されている。明らかな雑音対解像度のトレードオフ関係が示されているが（βが大きくなると、雑音が低くなり、空間解像度が犠牲になる）、画像全てにおいて、概ね縞状のアーチファクトがない。画像全ての窓及びレベルは、それぞれ５００及び１５０ＨＵである。
【０１０６】
図１０は、正則化パラメータの範囲内の、一方向の単一コンポーネントの例を示す。具体的には、前の研究と同じ正則化パラメータを用いて、βの値を１０倍大きくしたり１０分の１にしたりした。正則化のレベルを変化させた結果、βを小さくすると雑音が大きくなると予期され、βを大きくすると空間解像度が低下すると予期される。しかし画像全てにおいて、ネジによる光子不足による縞状のアーチファクトは基本的になかった。このことは、縞状のアーチファクトがなくなったのが単に正則化を行ったからではなく、ネジの構造及び組成に関する予備知識の統合が、実際に再構成を向上させる、ということを示唆している。正則化の値が低いと、よりゆっくり収束するように思われ、量の推定に雑音が多くなった場合、ＫＣＲの対象物関数は、局所的な最少値（特に登録パラメータに関して）の影響を受けやすいことがわかるだろう。変換推定の誤りは、図１０の、最低から最高のβで順序付けた３つの正則化のケースについて、０．３１、０．１３、及び０．０５ボクセルであった。同様に、平均絶対角度誤りは、０．２４、０．０９、及び０．０４度であった。したがって、βのこの特定の範囲においては、変換推定は、大きな正則化を行ったことの恩恵を受けている。
【０１０７】
＜セクションＩＩＩ：考察＞
この例は、オブジェクトのモデルに既知のコンポーネントを組み込むことを可能とするＫＣＲを参照する新規な断層撮影再構成法を説明する。コンポーネントの位置及び姿勢が、背景の減衰と同時に推定される。この技術の利用例は、脊柱手術において金属の脊椎頸ネジを含む画像に示したが、この技術は、形状及び組成がわかっているいずれの種類のコンポーネント（異質な組成のオブジェクトも含む）にも利用可能であるという意味で汎用性がある。シミュレーションされた脊椎頸ネジの再構成では、ＫＣＲは、光子不足に関するアーチファクトを大幅に低減させ、インプラントの境界までの生体組織をクリアに視覚化しており、従来のＦＢＰ及びＰＬ法が実質的にアーチファクトを含んでいるのと対照的である。したがい、ＫＣＲ法は、介入的なガイド（たとえば脊柱手術等）及び外科手術用工具が３Ｄ画質を曖昧にする他の介入のためのコーンビームＣＴ（ＣＢＣＴ）ベースの撮像分野において、及び、大きな将来性が見込めると思われる（[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）。さらに、ＫＣＲは、光子の不足するような場合において手術を成功させる可能性があるために、線量（dose）を低減させるための重要な手段となるだろう。ＫＣＲは、臀部及び膝のインプラント及びＣＴベースの生検針の案内等の他の用途においても重要だろう。端部を保存するペナルティ／予備知識を利用するＫＣＲ方法を単純に拡大させる方法（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）も利用することができ、これにより、雑音と解像度のトレードオフ関係がさらに向上する。
【０１０８】
上述した基本的なＫＣＲの枠組みは、ｘ線のスペクトルに関する多エネルギービーム及び物理学（たとえば分散）を利用しないフォワード・モデルを利用している。しかし、多エネルギーのフォワード・モデル（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）をフォワード・モデルに組み込むことは簡単である。
【０１０９】
一定のコンポーネントの知識が不正確である場合には（たとえば、外科手術時に脊柱の固定に利用される湾曲する金属ロッド、損耗が著しいので手術で直すことが必要な膝または臀部のインプラント、または、挿入時に撓む針の案内）、このＫＣＲの枠組みに対する実装に限界があることが想定される。しかし、（９）のマッピングによって、変形可能な登録モデルに簡単な修正を加えて、これらのシナリオのいくつかを対応可能とすることができるだろうことも付け加えておく（特に、差が、比較的低いパラメータの空間の次元で表されうるような場合）。
【０１１０】
上でまとめた結果は、ＫＣＲから生じる画質に重きをおいて説明されたが、推定値からは、画像量と登録パラメータという２つの出力がある。画像の量が一次産物ではあるが、姿勢の値の正確な推定自身も、体内に導入されたコンポーネントの位置を知るための医療的な価値があるだろう。たとえばＫＣＲ法の姿勢の出力によって、実際の道具を配置するための術前のプランを比較したり、外科手術の正確さを定量化したりする、といった有用性があるだろう。
【図１】
2014509037.tif 000039
【図２】
2014509037.tif 000040
【図３】
2014509037.tif 000041
【図４】
2014509037.tif 000042
【図５】
2014509037.tif 000043
【図６】
2014509037.tif 000044
【図７】
2014509037.tif 000045
【図８】
2014509037.tif 000046
【図９】
2014509037.tif 000047
【図１０】
2014509037.tif 000048
【国際調査報告】
2014509037.tif 000049
2014509037.tif 000050
2014509037.tif 000051