画像処理方法及びＸ線コンピュータ断層撮影装置
【課題】並列計算にも適用可能なＴＶ最小化反復再構成アルゴリズムを実現する画像処理方法及びＸ線コンピュータ断層撮影装置を提供すること。
【解決手段】データ収集回路１０４は、投影データを収集する。再構成部１１４は、投影データに基づいて画像ボリュームを再構成する。再構成部１１４は、投影データを各々が複数のビューを含む複数のサブセットに区分する。再構成部１１４は、画像ボリュームを更新するために複数のサブセットの各々について、複数のビューにＯＳ−ＳＡＲＴを並列的に実行する。再構成部１１４は、ＯＳ−ＳＡＲＴの完了後において、所定の規則に従って勾配ステップ幅を決定する。再構成部１１４は、勾配ステップ幅を利用して全変動を適応的に最小化する。
【解決手段】データ収集回路１０４は、投影データを収集する。再構成部１１４は、投影データに基づいて画像ボリュームを再構成する。再構成部１１４は、投影データを各々が複数のビューを含む複数のサブセットに区分する。再構成部１１４は、画像ボリュームを更新するために複数のサブセットの各々について、複数のビューにＯＳ−ＳＡＲＴを並列的に実行する。再構成部１１４は、ＯＳ−ＳＡＲＴの完了後において、所定の規則に従って勾配ステップ幅を決定する。再構成部１１４は、勾配ステップ幅を利用して全変動を適応的に最小化する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明の実施形態は、画像処理方法及びＸ線コンピュータ断層撮影装置に関する。
【背景技術】
【０００２】
ボリューム画像再構成に関して、以下の３つの例に示すような様々なグループによって反復アルゴリズムが開発されている。パン（Pan）博士らによるシカゴ大学グループは、疎ビュー（sparse view）Ｘ線ＣＴ再構成と角度制限Ｘ線ＣＴ再構成とに対する適用についての全変動（ＴＶ：total variation）最小化反復再構成アルゴリズムを提案した。ワン（Wang）博士らによるバージニアテクノロジー（Virginia technology）グループは、２００９年に、多ビューにおいて切り捨てられた投影データを有する関心領域（ＲＯｌ）再構成を目的としたＴＶ最小化アルゴリズムを発表した（すなわち、内部再構成問題）。最後に、ウィスコンシン大学マディソン校（Chen博士ら）は、先行画像拘束圧縮検出（ＰＩＣＣＳ：prior image constrained compressed sensing）方法を提案した。これら３つの先行技術においてＴＶは、同一画像上の画像ノイズを平滑化しエッジを維持するために利用される。従って、ＴＶは、画像処理における上記の効果を最適化するために最小化される。
【０００３】
以下の擬似コードは、非特許文献１に開示されたような１つの実施態様を示す。
【表１】
【０００４】
先行技術の取り組みにもかかわらず、幾つかの問題が未解決のままであり、改善が必要とされる。例えば、シカゴ大学グループの技術は、凸射影法（ＰＯＣＳ：projection on convex set）が用いられるので、画像処理を並列計算で実行することができない。この制約は、計算を終えるのに長時間かかるので、アルゴリズムを多ビューに関する３次元コーンビーム投影データに適用する場合において重要である。
【０００５】
シカゴ大学の手法は、正値性制限（positivity constraint）と計算集中的プロセスとを必要とする。実際に、この手法は、各レイについて画像ボリュームを更新する。例えば、ビュー毎に１００本のレイを有する９０個のビューの場合、シカゴ大学の手法は、第１の勾配降下ステップを実行する前に画像を９０００回（９０×１００）更新する。さらに、このＴＶ最小化ステップにおいて、固定のステップ幅（step size）と目的関数の正規化された勾配とが探索方向について利用される。典型的には、ステップ幅が極めて小さい値に固定されるので、ＴＶ最小化ステップは、多数の反復を必要とする。シカゴ大学グループの技術は、実測の投影データが負データである場合において、直接的に適用できないという正値性制限の追加の制限を有する。
【０００６】
一方、バージニアテクノロジーグループは、内部断層撮影法での多ビューからの投影データを利用する並列計算における画像処理を実施した（非特許文献２）。並列計算は、オーダードサブセット―同時代数的再構成法（ＯＳ−ＳＡＲＴ：ordered subset simultaneous algebraic reconstruction technique）に基づく。ＳＡＲＴと画像更新ステップとは、各サブセットについて反復される。すなわち、ＳＡＲＴが第１のサブセットにおいてのみ実行された後、画像は、直ちに更新される。これら２つの連続したステップは、各サブセットについて反復される。より詳細には、このＳＡＲＴステップにおいて、先行技術は、システムマトリクスの行列要素をキャッシュ（記憶）する必要がある。
【０００７】
以下の擬似コードは、非特許文献３で開示された１つの実施態様を示す。
【表２】
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】Sidky and Pan, Phys. Med. Biol., Vol.53. pp.4777-4807. 2008
【非特許文献２】Ge Wang and Ming Jiang, J of X-Ray Science and Technology. Pp 169-177,12 (2004)
【非特許文献３】Heugyong Yu and Ge Wang, Phys. Med. Biol. Pp.2791-2805, 54 (2009).
【発明の概要】
【発明が解決しようとする課題】
【０００９】
２００９年のバージニアテクノロジーの手法は、内部再構成問題を目的としている。しかしながら、この問題に厳密な解決策がないことが示されている。従って、この手法は、内部問題の近似解であると考えられる。また、２００９年のバージニアテクノロジーの手法は、高コントラストの対象に対する一時しのぎの解決策であり、さらに低コントラストの対象の撮像では多くの問題を有すると考えられる。さらに、２００９年のバージニアテクノロジーの手法は、ＴＶ最小化アルゴリズムの本来の長所である疎ビュー再構成問題に適用できないと考えられる。
【００１０】
ウィスコンシン大学のグループによる手法は、不都合に限定されている。この手法は、多くの場合において利用可能でない１組の先行画像を必要とする。さらに、この手法は、本質的にＰＯＣＳに基づいているので、並列計算を実行することができない。
【００１１】
以上の事情を考慮して、並列計算にも適用可能なＴＶ最小化反復再構成アルゴリズムを実施するための実用的な解決策が望まれている。
【００１２】
目的は、並列計算にも適用可能なＴＶ最小化反復再構成アルゴリズムを実現する画像処理方法及びＸ線コンピュータ断層撮影装置を提供することにある。
【課題を解決するための手段】
【００１３】
本実施形態に係る画像処理方法は、データ収集回路からの投影データを複数のサブセットに区分し、前記複数のサブセットの各々は複数のビューを含み、前記投影データに基づく画像ボリュームを更新するために、前記複数のサブセットの各々について前記複数のビューにＯＳ−ＳＡＲＴを並列的に実行し、前記ＯＳ−ＳＡＲＴの実行後、所定の規則に従って勾配ステップ幅を決定し、前記決定された勾配ステップ幅を利用して全変動を適応的に最小化する、ことを具備する。
【図面の簡単な説明】
【００１４】
【図１】本実施形態に係るＸ線コンピュータ断層撮影装置の構成を示す図。
【図２】本実施形態に係る全変動反復再構成（ＴＶ−ＩＲ）の処理の典型的な流れを示す図。
【図３】図２のオーダードサブセット同時代数的再構成（ＯＳ−ＳＡＲＴ）の処理の典型的な流れを示す図。
【図４】図２の直線探索法の処理の典型的な流れを示す図。
【図５】本実施形態に係るＴＶ−ＩＲに関する擬似コードのリストの一例を示す図。
【図６】本実施形態に係るＯＳ−ＳＡＲＴに関する疑似コードのリストの一例を示す図。
【図７】本実施形態に係る直線探索法に関する擬似コードのリストの一例を示す図。
【発明を実施するための形態】
【００１５】
以下、図面を参照しながら本実施形態に係わる画像処理方法及びＸ線コンピュータ断層撮影装置（Ｘ線ＣＴ装置）を説明する。図１は、本実施形態に係るＸ線ＣＴ装置の構成を示す図である。図１に示すように、Ｘ線ＣＴ装置は、マルチスライス型である。Ｘ線ＣＴ装置は、ガントリ１００を含む。図１に示されるガントリ１００は、側面から示されている。ガントリ１００は、Ｘ線管１０１、環状形状を有するフレーム１０２、及び多列又は２次元アレイ型のＸ線検出器１０３を含む。Ｘ線管１０１とＸ線検出器１０３とは、フレーム１０２上の被検体Ｓを横切る対角線上に搭載され、回転軸ＲＡに回転可能に支持されている。回転部１０７は、フレーム１０２を０．４秒／回転等の高速で回転させる。被検体Ｓは、回転軸ＲＡに沿って、図１の紙面に対する奥行き方向又は手前方向に移動される。
【００１６】
Ｘ線ＣＴ装置は、さらに、高電圧発生部１０９を含む。高電圧発生部１０９は、Ｘ線管１０１にＸ線を発生させるために、スリップリング１０８を介して管電圧をＸ線管１０１に印加する。Ｘ線は、被検体Ｓに向けて照射される。被検体Ｓの断面領域は、円で表わされている。Ｘ線検出器１０３は、被検体Ｓを透過したＸ線を検出するために、被検体Ｓを挟んでＸ線管１０１に対向する位置に配置される。
【００１７】
図１に示すように、Ｘ線ＣＴ装置は、さらに、Ｘ線検出器１０３からの検出信号を処理するための他の構成要素を含む。データ収集回路（ＤＡＳ：data acquisition system）１０４は、Ｘ線検出器１０３から出力された検出信号をチャネル毎に電圧信号に変換し、増幅し、デジタルデータに変換する。データ収集回路により発生されたデジタルデータは、生データとも呼ばれている。Ｘ線検出器１０３とＤＡＳ１０４とは、最大で９００ＴＰＰＲ（１回転当たりの総投影数：total number of projections per rotation）、９００ＴＰＰＲ〜１８００ＴＰＰＲ、及び９００ＴＰＰＲ〜３６００ＴＰＰＲで既定のＴＰＰＲを処理するように構成される。
【００１８】
デジタルデータは、非接触データ伝送部１０５によって、ガントリ１００の外部に設けられたコンソール内の前処理部１０６に供給される。前処理部１０６は、デジタルデータに感度補正等の前処理を施す。前処理が施されたデジタルデータは、投影データと呼ばれている。記憶部１１２は、再構成処理の前段階において、投影データを記憶する。記憶部１１２は、データ／制御バスを介してシステム制御部１１０に接続されている。また、記憶部１１２は、データ／制御バスを介して再構成部１１４、表示部１１６、入力部１１５、及びスキャン計画支援部２００に接続されている。スキャン計画支援部２００は、技術者がスキャン計画をたてるのを支援するための機能を有する。
【００１９】
再構成部１１４は、投影データに基づいて画像ボリュームを再構成する。再構成部１１４は、種々のソフトウェア構成要素とハードウェア構成要素とを含む。再構成部１１４は、並列計算に適切な反復再構成技法を利用してＴＶを最小化することに優れている。本実施形態に係る再構成部１１４は、全ボリューム反復再構成（ＴＶＩＲ：total volume iterative reconstruction）アルゴリズムを実行する。ＴＶＩＲアルゴリズムは、投影データにＯＳ−ＳＡＲＴステップとＴＶ最小化ステップとを実行する。２つのステップは、反復数が予め設定されているメインループにおいて連続して実行される。
【００２０】
ＴＶ最小化ステップの前段階において投影データは、ＯＳ−ＳＡＲＴに供される。投影データは、Ｎ個のサブセットに区分される。各サブセットは、特定数のビューを有する。なおＮは、予め設定される。ＯＳ−ＳＡＲＴにおいて各サブセットは、連続して処理されてもよい。あるいは、多重の中央処理装置（ＣＰＵ）やグラフィック処理装置（ＧＰＵ）などの特定のマイクロプロセッサを利用することによって、複数のサブセットは、並列的に処理されてもよい。
【００２１】
ＯＳ−ＳＡＲＴにおいて再構成部１１４は、２つの主要な処理を実行する。すなわち、再構成部１１４は、Ｎ個のサブセットの各々について、画像ボリュームを再投影して計算上の投影データを生成する。次に再構成部１１４は、生成された計算上の投影データと実測の投影データとの正規化された差分を逆投影（バックプロジェクション）して更新画像ボリュームを再構成する。より詳細には、再構成部１１４は、レイトレーシング法を利用することによって画像ボリュームを再投影する。レイトレーシング法においてシステムマトリクスの行列要素は、キャッシュ（記憶）されない。再構成部１１４は、サブセット内の全てのレイを同時に再投影する。この処理は、必要に応じて並列的に実行される。逆投影において、再構成部１１４は、ピクセルドリブン（pixel - driven）技術を利用して、サブセット内の全ての正規化された差分投影データを逆投影して、所望の更新画像ボリュームを生成する。再構成部１１４がサブセット内の全てのレイサム（すなわち、差分投影データ）を逆投影して画像ボリュームを生成するので、この処理は、必要に応じて並列的に実行される。これらの処理は、単一のＯＳ−ＳＡＲＴステップを完了するためにＮ個全てのサブセットについて施される。
【００２２】
ＴＶ最小化ステップにおいて再構成部１１４は、現在の画像ボリュームの目的関数が過去の画像ボリュームより小さくなるように、直線探索手法を利用して正のステップ幅を探索する。この手法において、再構成部１１４は、従来技術による固定のステップ幅に比して、ＴＶ最小化ステップにおいて可変のステップ幅を生成する。また、先行技術のパンの手法において探索手順は、負の正規化された勾配である。しかしながら、本実施形態において探索手順は、必要に応じて正規化されなくてもよい。
【００２３】
再構成部１１４は、分解能とノイズとを釣り合わせるために「ＴＶステップ」と呼ばれるパラメータを調整する。この調整処理は、ＴＶ最小化ステップにおいて実行される。再構成部１１４は、ＴＶ最小化ステップをＸ回だけ反復する。なおＸは、ある程度の分解能を犠牲にしながらノイズを抑制するために予め定められたパラメータである。再構成部１１４は、分解能とこのパラメータのノイズレベルとの間に成立するトレードオフを決定するのに優れている。
【００２４】
次に図２を参照しながら、本実施形態に係るＴＶ−ＩＲの処理の典型的な流れを説明する。図２に示される処理は、２つのループを含む。各ループにおいて特定のステップが反復される。第１のループ（以下、外側ループと呼ぶことにする。）は、ステップＳ２０、Ｓ３０、及びＳ４０を含む。第２のループ（以下、内側ループと呼ぶことにする。）は、ステップＳ３０及びＳ４０を含む。これらのループは、ある特定の条件が満たされるまで反復されるのではなく、所定の回数だけ反復される。ステップＳ２０及びＳ３０は、内側ループにおいて、ステップＳ１０で所定のＴＶステップ数が初期化されるように第１の所定数だけ反復される。内側ループが終了した時、ステップＳ２０、Ｓ３０、及びＳ４０は、外側ループにおいて、ステップＳ１０で所定の反復数が初期化されるように第２の所定数だけ反復される。
【００２５】
初期化ステップＳ１０において画像Ｘ０は、０に初期化される。また、ステップＳ１０においてＴＶステップ数は、第１の所定数に初期化され、反復数は、第２の所定数に初期化される。より詳細には、これらの所定数は、典型的には、経験的手法により決定される。ステップＳ２０においてＯＳ−ＳＡＲＴは、ＴＶ−ＩＲの処理過程中、実測の投影データに対して実行される。ステップＳ２０において、中間画像Ｘ１が出力される。ＯＳ−ＳＡＲＴの詳細は、後で図３を参照しながら説明する。ステップＳ３０において、所定のＴＶステップ数に基づいてＴＶを最小にするように、直線探索法が繰り返し実行される。ステップＳ３０が終了すると、ステップＳ２０における中間画像Ｘ１から画像Ｘ２が生成される。ステップＳ３０において生成された画像は、外側ループが再びステップＳ２０から繰り返される前にステップＳ４０において中間画像ホルダに記憶される。または、ステップＳ３０において生成された画像は、ステップＳ５０において画像Ｘ２が出力される前に、ステップＳ４０において変数Ｘ１に置き換えられる。
【００２６】
図３に示すフローチャートは、図２のステップ２０におけるＯＳ−ＳＡＲＴの処理の典型的な流れを示す。図示されたＯＳ−ＳＡＲＴの処理は、特定のステップが反復して実行される１つのループを含む。ループは、ステップＳ２３、Ｓ２４、及びＳ２５を含み、ある特定の条件が満たされるまで反復されるのではなく、所定数Ｎに従って反復される。ループが終了した場合、ＯＳ−ＳＡＲＴは、中間画像として変数Ｘ１を出力する。
【００２７】
図３に示すように、ステップＳ２１において、図２のステップＳ２０に示される画像Ｘ０と実測の投影データｐとが取り込まれる。ステップＳ２２において、実測の投影データｐは、サブセットと呼ばれるＮsets個のグループに区分される。Ｎ個のサブセットは、Ｎsubrays＝Ｎrays／Ｎsets個のレイサムをそれぞれ含む。Ｎ個サブセットの各々について、ステップＳ２３において、画像Ｘ０は再投影され計算上の投影データが生成される。ステップＳ２４において、実測の投影データと計算上の投影データとの間の正規化差分は、逆投影される。より詳細には、ステップＳ２３において、システムマトリクスの行列要素がキャッシュされないレイトレーシング（光線追跡）技術が利用されることによって、画像ボリュームが再投影される。さらに、ステップＳ２３において、サブセット内の全てのレイが同時に再投影される。この処理は、必要に応じて並列的に実行される。ステップＳ２４において、ピクセルドリブン技術を利用して、サブセット内の全ての正規化差分投影データが逆投影され、所望の更新画像ボリュームが生成される。ステップＳ２４において、サブセット内の全てのレイサム、すなわち正規化差分投影データが逆投影されることによって画像ボリュームが生成される。従って、この処理は、必要に応じて並列的に実行される。これらの処理がＮ個全てのサブセットに適用されることによって、単一のＯＳ−ＳＡＲＴステップが完了する。最後に、ステップＳ２５は、ステップＳ２４からの逆投影画像を画像Ｘ０に追加することによって画像Ｘ０を更新する。ステップＳ２５において生成された画像は、ステップＳ２３に進むことによりループが繰り返される前に中間画像ホルダに記憶される、又はステップＳ２６において画像Ｘ１が出力される前にＸ１に置き換えられる。
【００２８】
図４に示すフローチャートは、図２のステップＳ３０における直線探索法の処理の典型的な流れを示す。この直線探索法は、ＴＶ−ＩＲの処理過程において実行される。図示された直線探索法に関するフローチャートは、特定のステップが反復して実行される１つのループを含む。このループは、ステップＳ３２、Ｓ３３、及びＳ３４を含む。このループは、ある特定の停止条件が満たされるまでではなく、所定回数だけ反復される。ループが終了した時、直線探索法は、変数Ｘ２として画像を出力し、ステップ幅として変数ａを出力する。
【００２９】
図４に示すように、ステップＳ３１において画像Ｘ１は、図２のステップＳ３０を介してステップＳ２０から取り込まれる。さらにステップＳ３１においてステップ幅変数ａは、１．０等の所定値に初期化される。ステップＳ３２において画像Ｘ２は、式（Ｘ１−ａ＊ｇｒａｄ（ＴＶ（Ｘ１）））に従って計算される。ここで、ａはステップ幅変数であり、ｇｒａｄ（・）は所定の勾配演算子であり、ＴＶ（・）は所定の全変動演算子（ＴＶ演算子）である。同様にステップＳ３２において、画像変数Ｘ２が新しく割り当てられた後、ＴＶ（Ｘ２）が計算される。ここで、ＴＶ（・）は、同一の所定のＴＶ演算子である。ステップＳ３３において、ＴＶ（Ｘ２）の値とＴＶ（Ｘ１）の値とが比較される。ＴＶ（Ｘ２）＜ＴＶ（Ｘ１）が真でない場合、ステップＳ３４において現行のステップ幅値が半分にされる。ステップＳ３４が終了するとステップＳ３２へ進む。一方、ＴＶ（Ｘ２）＜ＴＶ（Ｘ１）が真の場合、ステップＳ３５において現行のステップ幅値ａと画像Ｘ２とが出力される。換言すれば、前述のステップＳ３２、Ｓ３３、及びＳ３４は、ＴＶ（Ｘ２）＜ＴＶ（Ｘ１）が真になるまで繰り返される。
【００３０】
図５は、本実施形態に係るＴＶＩＲに関する擬似コードのリストの一例を示す。メインのプロシージャ（procedure）は、画像の外観に影響する、第３行〜第１０行の反復数と第５行〜第８行のＴＶステップ数とを決定する２つのパラメータを有する。ノイズレベルを抑制するために、所定数だけ反復されるＴＶステップが利用される。擬似コードにおいて、反復数はＮiterで示され、ＴＶステップ数はＮTVで示される。Ｎiterは、ビュー数に反比例し、ＮTVはノイズレベルに比例する。プログラムを終了するための明確な収束判断基準が定義されないことに注意されたい。その代わりに、プログラムは、ある特定の所定回数だけ繰り返された後で停止する。
【００３１】
本実施形態においては、例えば、ビュー数が９００以上であり、投影データがほぼノイズを含まない場合、例えば、Ｎiterは２０に設定され、ＮTVは１に設定される。一方、投影データが多くのノイズを含む場合、ＮTVは、例えば、３〜９の何れかに設定される。実際には、ＮTVは９を超える値に設定される場合もある。ＮTVが大きいほど再構成画像が滑らかになるが、空間分解能が低下する。
【００３２】
図５の擬似コードにおいて、以下の表記法を用いることにする。画像ｆは、以下の（１）式において、システムマトリクスＡと列ベクトルｐによって示される投影データとによって表される。
【００３３】
ｐ＝Ａｆ （１）
画像ｆは、２次元においてＹＤＩＭ×ＸＤＩＭ次元を有し、３次元においてＺＤＩＭ×ＹＤＩＭ×ＸＤＩＭ次元を有する。従って、画像は、２次元において２次元行列ｆYDIM,XDIMにより表され、３次元において３次元行列ｆZDIM,YDIM,XDIMにより表される。順投影プロシージャにおいて、画像ｆは、２次元ではＹＤＩＭ×ＸＤＩＭを有する列ベクトルに再配列され、３次元ではＺＤＩＭ×ＹＤＩＭ×ＸＤＩＭを有する列ベクトルに再配列される。Ｎpixelsは、画像内の総画素数を表す。すなわち、２次元においてＮpixels＝ＹＤＩＭ×ＸＤＩＭであり、３次元においてＮpixels＝ＺＤＩＭ×ＹＤＩＭ×ＸＤＩＭである。
【００３４】
投影データｐは、投影ビュー数ＶＩＥＷＳを有する。投影データは、２次元データの場合、チャンネル数（ＣＨＡＮＮＥＬＳ）に対応する数のレイを有する。３次元の場合、投影データは、ＣＨＡＮＮＥＬＳ×セグメント数（ＳＥＧＭＥＮＴＳ）に対応する数のレイを有する。ＳＥＧＭＥＮＴＳは、追加の次元を示す。従って、総レイ数が、Ｎraysによって示される場合、２次元投影データは、Ｎrays＝ビュー数（ＶＩＥＷＳ）×ＣＨＡＮＮＥＬＳを有し、３次元投影データは、Ｎrays＝ＶＩＥＷＳ×ＳＥＧＭＥＮＴＳ×ＣＨＡＮＮＥＬＳを有する。投影データは、Ｎrays個の要素を有する列ベクトルｐによって示される。列ベクトルｐのｉ番目の要素は、ｉ番目のレイサムであるｐiによって示される。
【００３５】
システムマトリクスＡの次元は、Ｎrays×Ｎpixelsである。その行列要素は、ａi,jによって示される。従って、ｉ番目のレイサムは、（２）式で以下のように表わせる。
【数１】
【００３６】
システムマトリクスＡ内の全ての行列要素を記憶することは、その膨大な次元数のために膨大なコストがかかる。例えば、各ビューの８９６個のチャネルにおいて９００個のビューが測定され、サイズ５１２×５１２の画像行列が再構成される場合、システムマトリクスＡは、８０６，４００×２６２，１４４次元を有する。システムマトリクスＡは疎であるが、全ての行列要素を記憶するにはコストが高い。システムマトリクスＡを記憶しないようにする１つの方法は、転置を必要としない場合において、行列要素をその場で計算することである。１つの先行技術の手法において、１つのレイの行列要素だけが計算され記憶され、レイ毎に画像ボリュームが更新される。本実施形態においてシステムマトリクスＡは記憶されないが、その転置を逆投影演算子であると解釈する。
【００３７】
目的関数（費用関数）は、以下に示すように、２次元の場合に（３）式で規定され、３次元の場合に（４）式で規定されるようなＴＶである。変数εは、後述の探索方向の式によりシグマ内の各ｕ（・）項が分割されるので、各ｕ（・）項がゼロになることを防止するための微小量である。（３）式と（４）式とは、一連の離散勾配のｌlノルムを表わす。
【数２】
【００３８】
【数３】
【００３９】
図６は、図２のステップＳ２０に関する擬似コードのリストの一例を示す。この擬似コードが示す処理は、本実施形態に係るＴＶ−ＩＲにおけるＯＳ−ＳＡＲＴに含まれる。よく研究されたＯＳ―ＳＡＲＴを検討してきたが、本実施形態に係るＯＳ−ＳＡＲＴの例示的な実施態様は、先行技術で示されたようなＰＯＣＳの変形である。本実施形態は、ある特定のＯＳ−ＳＡＲＴが適正なＰＯＣＳ処理であることを証明するものではない。ＯＳ−ＳＡＲＴの利点は、必要に応じて並列で実施されることである。本実施形態に係るＯＳ−ＳＡＲＴにおいてシステムマトリクスＡは、キャッシュされない。
【００４０】
図６に示すように、ＯＳ−ＳＡＲＴにおいて投影データｐは、Ｎsets個のサブセットに区分される。ｔ番目のサブセットは、Ｎsubrays＝Ｎrays／Ｎsetsのレイサムを含むＳtによって示される。数学的には、ＯＳ−ＳＡＲＴの更新式は、以下の（５）式によって示される。
【数４】
【００４１】
ここで、システムマトリクスＡの行列要素と投影データｐとは、サブセットＳｔ内のレイに対応する。前述の式においてシステムマトリクスＡの転置が必要とされるが、本実施形態においては、システムマトリクスＡの行列要素を記憶せず、転置処理を逆投影演算子により実行する。（５）式中の３つの項は、以下の（５ａ）式、（５ｂ）式、（５ｃ）式のようにそれぞれ表される。
【数５】
【００４２】
第１項は、画像ｆ(k,t)の再投影である。
【数６】
【００４３】
第２項は、ｌ番目のレイの長さである。
【数７】
【００４４】
第３項は、ｊ番目の画素が逆投影された回数として表されてもよい。
【００４５】
図７は、図２のステップＳ３０に関する擬似コードのリストの一例を示す。この擬似コードが示す処理は、本実施形態に係るＴＶ−ＩＲにおける直線探索法に関する処理に含まれる。直線探索法は、ステップ幅αを求めるために利用される。ステップ幅αは、以下の（６）式又は図５の擬似コードの７．１行に示される更新式で利用される。
【数８】
【００４６】
その結果、ＵTV（ｆ(k+l)）≦ＵTV（ｆ(k)）＋μαｄTsearch▽ＵTV（ｆ(k)）となる。ここで、μは、０＜μ＜１を満たすスカラである。例えば、μ＝０．００１に設定される。
【００４７】
探索方向は、（３）式と（４）式とに示したような目的関数を最小化するために、勾配降下法において利用される。探索方向ｄsearchは、以下の（７）式のように規定される。
【数９】
【００４８】
（７）式は、Ｎpixels個の要素の列ベクトルとして表わされる。これは、２次元の場合（８）式により、３次元の場合（９）式により示される。
【数１０】
【００４９】
【数１１】
【００５０】
２次元の場合、（３）式の３つの項、すなわちｕ（ｋ−１，ｌ）、ｕ（ｋ，ｌ−１）、及びｕ（ｋ，ｌ）だけが変数ｆk,lを含む。これらの項は、（８）式を得るために、ｆk,lについてのみ微分される。
【００５１】
本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。
【符号の説明】
【００５２】
１００…ガントリ、１０１…Ｘ線管、１０２…フレーム、１０３…Ｘ線検出器、１０４…データ取得回路（ＤＡＳ）、１０５…非接触データ伝送部、１０６…前処理部、１０７…回転部、１０８…スリップリング、１０９…高電圧発生部、１１０…システム制御部、１１２…記憶部、１１４…再構成部、１１５…入力部、１１６…表示部
【技術分野】
【０００１】
本発明の実施形態は、画像処理方法及びＸ線コンピュータ断層撮影装置に関する。
【背景技術】
【０００２】
ボリューム画像再構成に関して、以下の３つの例に示すような様々なグループによって反復アルゴリズムが開発されている。パン（Pan）博士らによるシカゴ大学グループは、疎ビュー（sparse view）Ｘ線ＣＴ再構成と角度制限Ｘ線ＣＴ再構成とに対する適用についての全変動（ＴＶ：total variation）最小化反復再構成アルゴリズムを提案した。ワン（Wang）博士らによるバージニアテクノロジー（Virginia technology）グループは、２００９年に、多ビューにおいて切り捨てられた投影データを有する関心領域（ＲＯｌ）再構成を目的としたＴＶ最小化アルゴリズムを発表した（すなわち、内部再構成問題）。最後に、ウィスコンシン大学マディソン校（Chen博士ら）は、先行画像拘束圧縮検出（ＰＩＣＣＳ：prior image constrained compressed sensing）方法を提案した。これら３つの先行技術においてＴＶは、同一画像上の画像ノイズを平滑化しエッジを維持するために利用される。従って、ＴＶは、画像処理における上記の効果を最適化するために最小化される。
【０００３】
以下の擬似コードは、非特許文献１に開示されたような１つの実施態様を示す。
【表１】
【０００４】
先行技術の取り組みにもかかわらず、幾つかの問題が未解決のままであり、改善が必要とされる。例えば、シカゴ大学グループの技術は、凸射影法（ＰＯＣＳ：projection on convex set）が用いられるので、画像処理を並列計算で実行することができない。この制約は、計算を終えるのに長時間かかるので、アルゴリズムを多ビューに関する３次元コーンビーム投影データに適用する場合において重要である。
【０００５】
シカゴ大学の手法は、正値性制限（positivity constraint）と計算集中的プロセスとを必要とする。実際に、この手法は、各レイについて画像ボリュームを更新する。例えば、ビュー毎に１００本のレイを有する９０個のビューの場合、シカゴ大学の手法は、第１の勾配降下ステップを実行する前に画像を９０００回（９０×１００）更新する。さらに、このＴＶ最小化ステップにおいて、固定のステップ幅（step size）と目的関数の正規化された勾配とが探索方向について利用される。典型的には、ステップ幅が極めて小さい値に固定されるので、ＴＶ最小化ステップは、多数の反復を必要とする。シカゴ大学グループの技術は、実測の投影データが負データである場合において、直接的に適用できないという正値性制限の追加の制限を有する。
【０００６】
一方、バージニアテクノロジーグループは、内部断層撮影法での多ビューからの投影データを利用する並列計算における画像処理を実施した（非特許文献２）。並列計算は、オーダードサブセット―同時代数的再構成法（ＯＳ−ＳＡＲＴ：ordered subset simultaneous algebraic reconstruction technique）に基づく。ＳＡＲＴと画像更新ステップとは、各サブセットについて反復される。すなわち、ＳＡＲＴが第１のサブセットにおいてのみ実行された後、画像は、直ちに更新される。これら２つの連続したステップは、各サブセットについて反復される。より詳細には、このＳＡＲＴステップにおいて、先行技術は、システムマトリクスの行列要素をキャッシュ（記憶）する必要がある。
【０００７】
以下の擬似コードは、非特許文献３で開示された１つの実施態様を示す。
【表２】
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】Sidky and Pan, Phys. Med. Biol., Vol.53. pp.4777-4807. 2008
【非特許文献２】Ge Wang and Ming Jiang, J of X-Ray Science and Technology. Pp 169-177,12 (2004)
【非特許文献３】Heugyong Yu and Ge Wang, Phys. Med. Biol. Pp.2791-2805, 54 (2009).
【発明の概要】
【発明が解決しようとする課題】
【０００９】
２００９年のバージニアテクノロジーの手法は、内部再構成問題を目的としている。しかしながら、この問題に厳密な解決策がないことが示されている。従って、この手法は、内部問題の近似解であると考えられる。また、２００９年のバージニアテクノロジーの手法は、高コントラストの対象に対する一時しのぎの解決策であり、さらに低コントラストの対象の撮像では多くの問題を有すると考えられる。さらに、２００９年のバージニアテクノロジーの手法は、ＴＶ最小化アルゴリズムの本来の長所である疎ビュー再構成問題に適用できないと考えられる。
【００１０】
ウィスコンシン大学のグループによる手法は、不都合に限定されている。この手法は、多くの場合において利用可能でない１組の先行画像を必要とする。さらに、この手法は、本質的にＰＯＣＳに基づいているので、並列計算を実行することができない。
【００１１】
以上の事情を考慮して、並列計算にも適用可能なＴＶ最小化反復再構成アルゴリズムを実施するための実用的な解決策が望まれている。
【００１２】
目的は、並列計算にも適用可能なＴＶ最小化反復再構成アルゴリズムを実現する画像処理方法及びＸ線コンピュータ断層撮影装置を提供することにある。
【課題を解決するための手段】
【００１３】
本実施形態に係る画像処理方法は、データ収集回路からの投影データを複数のサブセットに区分し、前記複数のサブセットの各々は複数のビューを含み、前記投影データに基づく画像ボリュームを更新するために、前記複数のサブセットの各々について前記複数のビューにＯＳ−ＳＡＲＴを並列的に実行し、前記ＯＳ−ＳＡＲＴの実行後、所定の規則に従って勾配ステップ幅を決定し、前記決定された勾配ステップ幅を利用して全変動を適応的に最小化する、ことを具備する。
【図面の簡単な説明】
【００１４】
【図１】本実施形態に係るＸ線コンピュータ断層撮影装置の構成を示す図。
【図２】本実施形態に係る全変動反復再構成（ＴＶ−ＩＲ）の処理の典型的な流れを示す図。
【図３】図２のオーダードサブセット同時代数的再構成（ＯＳ−ＳＡＲＴ）の処理の典型的な流れを示す図。
【図４】図２の直線探索法の処理の典型的な流れを示す図。
【図５】本実施形態に係るＴＶ−ＩＲに関する擬似コードのリストの一例を示す図。
【図６】本実施形態に係るＯＳ−ＳＡＲＴに関する疑似コードのリストの一例を示す図。
【図７】本実施形態に係る直線探索法に関する擬似コードのリストの一例を示す図。
【発明を実施するための形態】
【００１５】
以下、図面を参照しながら本実施形態に係わる画像処理方法及びＸ線コンピュータ断層撮影装置（Ｘ線ＣＴ装置）を説明する。図１は、本実施形態に係るＸ線ＣＴ装置の構成を示す図である。図１に示すように、Ｘ線ＣＴ装置は、マルチスライス型である。Ｘ線ＣＴ装置は、ガントリ１００を含む。図１に示されるガントリ１００は、側面から示されている。ガントリ１００は、Ｘ線管１０１、環状形状を有するフレーム１０２、及び多列又は２次元アレイ型のＸ線検出器１０３を含む。Ｘ線管１０１とＸ線検出器１０３とは、フレーム１０２上の被検体Ｓを横切る対角線上に搭載され、回転軸ＲＡに回転可能に支持されている。回転部１０７は、フレーム１０２を０．４秒／回転等の高速で回転させる。被検体Ｓは、回転軸ＲＡに沿って、図１の紙面に対する奥行き方向又は手前方向に移動される。
【００１６】
Ｘ線ＣＴ装置は、さらに、高電圧発生部１０９を含む。高電圧発生部１０９は、Ｘ線管１０１にＸ線を発生させるために、スリップリング１０８を介して管電圧をＸ線管１０１に印加する。Ｘ線は、被検体Ｓに向けて照射される。被検体Ｓの断面領域は、円で表わされている。Ｘ線検出器１０３は、被検体Ｓを透過したＸ線を検出するために、被検体Ｓを挟んでＸ線管１０１に対向する位置に配置される。
【００１７】
図１に示すように、Ｘ線ＣＴ装置は、さらに、Ｘ線検出器１０３からの検出信号を処理するための他の構成要素を含む。データ収集回路（ＤＡＳ：data acquisition system）１０４は、Ｘ線検出器１０３から出力された検出信号をチャネル毎に電圧信号に変換し、増幅し、デジタルデータに変換する。データ収集回路により発生されたデジタルデータは、生データとも呼ばれている。Ｘ線検出器１０３とＤＡＳ１０４とは、最大で９００ＴＰＰＲ（１回転当たりの総投影数：total number of projections per rotation）、９００ＴＰＰＲ〜１８００ＴＰＰＲ、及び９００ＴＰＰＲ〜３６００ＴＰＰＲで既定のＴＰＰＲを処理するように構成される。
【００１８】
デジタルデータは、非接触データ伝送部１０５によって、ガントリ１００の外部に設けられたコンソール内の前処理部１０６に供給される。前処理部１０６は、デジタルデータに感度補正等の前処理を施す。前処理が施されたデジタルデータは、投影データと呼ばれている。記憶部１１２は、再構成処理の前段階において、投影データを記憶する。記憶部１１２は、データ／制御バスを介してシステム制御部１１０に接続されている。また、記憶部１１２は、データ／制御バスを介して再構成部１１４、表示部１１６、入力部１１５、及びスキャン計画支援部２００に接続されている。スキャン計画支援部２００は、技術者がスキャン計画をたてるのを支援するための機能を有する。
【００１９】
再構成部１１４は、投影データに基づいて画像ボリュームを再構成する。再構成部１１４は、種々のソフトウェア構成要素とハードウェア構成要素とを含む。再構成部１１４は、並列計算に適切な反復再構成技法を利用してＴＶを最小化することに優れている。本実施形態に係る再構成部１１４は、全ボリューム反復再構成（ＴＶＩＲ：total volume iterative reconstruction）アルゴリズムを実行する。ＴＶＩＲアルゴリズムは、投影データにＯＳ−ＳＡＲＴステップとＴＶ最小化ステップとを実行する。２つのステップは、反復数が予め設定されているメインループにおいて連続して実行される。
【００２０】
ＴＶ最小化ステップの前段階において投影データは、ＯＳ−ＳＡＲＴに供される。投影データは、Ｎ個のサブセットに区分される。各サブセットは、特定数のビューを有する。なおＮは、予め設定される。ＯＳ−ＳＡＲＴにおいて各サブセットは、連続して処理されてもよい。あるいは、多重の中央処理装置（ＣＰＵ）やグラフィック処理装置（ＧＰＵ）などの特定のマイクロプロセッサを利用することによって、複数のサブセットは、並列的に処理されてもよい。
【００２１】
ＯＳ−ＳＡＲＴにおいて再構成部１１４は、２つの主要な処理を実行する。すなわち、再構成部１１４は、Ｎ個のサブセットの各々について、画像ボリュームを再投影して計算上の投影データを生成する。次に再構成部１１４は、生成された計算上の投影データと実測の投影データとの正規化された差分を逆投影（バックプロジェクション）して更新画像ボリュームを再構成する。より詳細には、再構成部１１４は、レイトレーシング法を利用することによって画像ボリュームを再投影する。レイトレーシング法においてシステムマトリクスの行列要素は、キャッシュ（記憶）されない。再構成部１１４は、サブセット内の全てのレイを同時に再投影する。この処理は、必要に応じて並列的に実行される。逆投影において、再構成部１１４は、ピクセルドリブン（pixel - driven）技術を利用して、サブセット内の全ての正規化された差分投影データを逆投影して、所望の更新画像ボリュームを生成する。再構成部１１４がサブセット内の全てのレイサム（すなわち、差分投影データ）を逆投影して画像ボリュームを生成するので、この処理は、必要に応じて並列的に実行される。これらの処理は、単一のＯＳ−ＳＡＲＴステップを完了するためにＮ個全てのサブセットについて施される。
【００２２】
ＴＶ最小化ステップにおいて再構成部１１４は、現在の画像ボリュームの目的関数が過去の画像ボリュームより小さくなるように、直線探索手法を利用して正のステップ幅を探索する。この手法において、再構成部１１４は、従来技術による固定のステップ幅に比して、ＴＶ最小化ステップにおいて可変のステップ幅を生成する。また、先行技術のパンの手法において探索手順は、負の正規化された勾配である。しかしながら、本実施形態において探索手順は、必要に応じて正規化されなくてもよい。
【００２３】
再構成部１１４は、分解能とノイズとを釣り合わせるために「ＴＶステップ」と呼ばれるパラメータを調整する。この調整処理は、ＴＶ最小化ステップにおいて実行される。再構成部１１４は、ＴＶ最小化ステップをＸ回だけ反復する。なおＸは、ある程度の分解能を犠牲にしながらノイズを抑制するために予め定められたパラメータである。再構成部１１４は、分解能とこのパラメータのノイズレベルとの間に成立するトレードオフを決定するのに優れている。
【００２４】
次に図２を参照しながら、本実施形態に係るＴＶ−ＩＲの処理の典型的な流れを説明する。図２に示される処理は、２つのループを含む。各ループにおいて特定のステップが反復される。第１のループ（以下、外側ループと呼ぶことにする。）は、ステップＳ２０、Ｓ３０、及びＳ４０を含む。第２のループ（以下、内側ループと呼ぶことにする。）は、ステップＳ３０及びＳ４０を含む。これらのループは、ある特定の条件が満たされるまで反復されるのではなく、所定の回数だけ反復される。ステップＳ２０及びＳ３０は、内側ループにおいて、ステップＳ１０で所定のＴＶステップ数が初期化されるように第１の所定数だけ反復される。内側ループが終了した時、ステップＳ２０、Ｓ３０、及びＳ４０は、外側ループにおいて、ステップＳ１０で所定の反復数が初期化されるように第２の所定数だけ反復される。
【００２５】
初期化ステップＳ１０において画像Ｘ０は、０に初期化される。また、ステップＳ１０においてＴＶステップ数は、第１の所定数に初期化され、反復数は、第２の所定数に初期化される。より詳細には、これらの所定数は、典型的には、経験的手法により決定される。ステップＳ２０においてＯＳ−ＳＡＲＴは、ＴＶ−ＩＲの処理過程中、実測の投影データに対して実行される。ステップＳ２０において、中間画像Ｘ１が出力される。ＯＳ−ＳＡＲＴの詳細は、後で図３を参照しながら説明する。ステップＳ３０において、所定のＴＶステップ数に基づいてＴＶを最小にするように、直線探索法が繰り返し実行される。ステップＳ３０が終了すると、ステップＳ２０における中間画像Ｘ１から画像Ｘ２が生成される。ステップＳ３０において生成された画像は、外側ループが再びステップＳ２０から繰り返される前にステップＳ４０において中間画像ホルダに記憶される。または、ステップＳ３０において生成された画像は、ステップＳ５０において画像Ｘ２が出力される前に、ステップＳ４０において変数Ｘ１に置き換えられる。
【００２６】
図３に示すフローチャートは、図２のステップ２０におけるＯＳ−ＳＡＲＴの処理の典型的な流れを示す。図示されたＯＳ−ＳＡＲＴの処理は、特定のステップが反復して実行される１つのループを含む。ループは、ステップＳ２３、Ｓ２４、及びＳ２５を含み、ある特定の条件が満たされるまで反復されるのではなく、所定数Ｎに従って反復される。ループが終了した場合、ＯＳ−ＳＡＲＴは、中間画像として変数Ｘ１を出力する。
【００２７】
図３に示すように、ステップＳ２１において、図２のステップＳ２０に示される画像Ｘ０と実測の投影データｐとが取り込まれる。ステップＳ２２において、実測の投影データｐは、サブセットと呼ばれるＮsets個のグループに区分される。Ｎ個のサブセットは、Ｎsubrays＝Ｎrays／Ｎsets個のレイサムをそれぞれ含む。Ｎ個サブセットの各々について、ステップＳ２３において、画像Ｘ０は再投影され計算上の投影データが生成される。ステップＳ２４において、実測の投影データと計算上の投影データとの間の正規化差分は、逆投影される。より詳細には、ステップＳ２３において、システムマトリクスの行列要素がキャッシュされないレイトレーシング（光線追跡）技術が利用されることによって、画像ボリュームが再投影される。さらに、ステップＳ２３において、サブセット内の全てのレイが同時に再投影される。この処理は、必要に応じて並列的に実行される。ステップＳ２４において、ピクセルドリブン技術を利用して、サブセット内の全ての正規化差分投影データが逆投影され、所望の更新画像ボリュームが生成される。ステップＳ２４において、サブセット内の全てのレイサム、すなわち正規化差分投影データが逆投影されることによって画像ボリュームが生成される。従って、この処理は、必要に応じて並列的に実行される。これらの処理がＮ個全てのサブセットに適用されることによって、単一のＯＳ−ＳＡＲＴステップが完了する。最後に、ステップＳ２５は、ステップＳ２４からの逆投影画像を画像Ｘ０に追加することによって画像Ｘ０を更新する。ステップＳ２５において生成された画像は、ステップＳ２３に進むことによりループが繰り返される前に中間画像ホルダに記憶される、又はステップＳ２６において画像Ｘ１が出力される前にＸ１に置き換えられる。
【００２８】
図４に示すフローチャートは、図２のステップＳ３０における直線探索法の処理の典型的な流れを示す。この直線探索法は、ＴＶ−ＩＲの処理過程において実行される。図示された直線探索法に関するフローチャートは、特定のステップが反復して実行される１つのループを含む。このループは、ステップＳ３２、Ｓ３３、及びＳ３４を含む。このループは、ある特定の停止条件が満たされるまでではなく、所定回数だけ反復される。ループが終了した時、直線探索法は、変数Ｘ２として画像を出力し、ステップ幅として変数ａを出力する。
【００２９】
図４に示すように、ステップＳ３１において画像Ｘ１は、図２のステップＳ３０を介してステップＳ２０から取り込まれる。さらにステップＳ３１においてステップ幅変数ａは、１．０等の所定値に初期化される。ステップＳ３２において画像Ｘ２は、式（Ｘ１−ａ＊ｇｒａｄ（ＴＶ（Ｘ１）））に従って計算される。ここで、ａはステップ幅変数であり、ｇｒａｄ（・）は所定の勾配演算子であり、ＴＶ（・）は所定の全変動演算子（ＴＶ演算子）である。同様にステップＳ３２において、画像変数Ｘ２が新しく割り当てられた後、ＴＶ（Ｘ２）が計算される。ここで、ＴＶ（・）は、同一の所定のＴＶ演算子である。ステップＳ３３において、ＴＶ（Ｘ２）の値とＴＶ（Ｘ１）の値とが比較される。ＴＶ（Ｘ２）＜ＴＶ（Ｘ１）が真でない場合、ステップＳ３４において現行のステップ幅値が半分にされる。ステップＳ３４が終了するとステップＳ３２へ進む。一方、ＴＶ（Ｘ２）＜ＴＶ（Ｘ１）が真の場合、ステップＳ３５において現行のステップ幅値ａと画像Ｘ２とが出力される。換言すれば、前述のステップＳ３２、Ｓ３３、及びＳ３４は、ＴＶ（Ｘ２）＜ＴＶ（Ｘ１）が真になるまで繰り返される。
【００３０】
図５は、本実施形態に係るＴＶＩＲに関する擬似コードのリストの一例を示す。メインのプロシージャ（procedure）は、画像の外観に影響する、第３行〜第１０行の反復数と第５行〜第８行のＴＶステップ数とを決定する２つのパラメータを有する。ノイズレベルを抑制するために、所定数だけ反復されるＴＶステップが利用される。擬似コードにおいて、反復数はＮiterで示され、ＴＶステップ数はＮTVで示される。Ｎiterは、ビュー数に反比例し、ＮTVはノイズレベルに比例する。プログラムを終了するための明確な収束判断基準が定義されないことに注意されたい。その代わりに、プログラムは、ある特定の所定回数だけ繰り返された後で停止する。
【００３１】
本実施形態においては、例えば、ビュー数が９００以上であり、投影データがほぼノイズを含まない場合、例えば、Ｎiterは２０に設定され、ＮTVは１に設定される。一方、投影データが多くのノイズを含む場合、ＮTVは、例えば、３〜９の何れかに設定される。実際には、ＮTVは９を超える値に設定される場合もある。ＮTVが大きいほど再構成画像が滑らかになるが、空間分解能が低下する。
【００３２】
図５の擬似コードにおいて、以下の表記法を用いることにする。画像ｆは、以下の（１）式において、システムマトリクスＡと列ベクトルｐによって示される投影データとによって表される。
【００３３】
ｐ＝Ａｆ （１）
画像ｆは、２次元においてＹＤＩＭ×ＸＤＩＭ次元を有し、３次元においてＺＤＩＭ×ＹＤＩＭ×ＸＤＩＭ次元を有する。従って、画像は、２次元において２次元行列ｆYDIM,XDIMにより表され、３次元において３次元行列ｆZDIM,YDIM,XDIMにより表される。順投影プロシージャにおいて、画像ｆは、２次元ではＹＤＩＭ×ＸＤＩＭを有する列ベクトルに再配列され、３次元ではＺＤＩＭ×ＹＤＩＭ×ＸＤＩＭを有する列ベクトルに再配列される。Ｎpixelsは、画像内の総画素数を表す。すなわち、２次元においてＮpixels＝ＹＤＩＭ×ＸＤＩＭであり、３次元においてＮpixels＝ＺＤＩＭ×ＹＤＩＭ×ＸＤＩＭである。
【００３４】
投影データｐは、投影ビュー数ＶＩＥＷＳを有する。投影データは、２次元データの場合、チャンネル数（ＣＨＡＮＮＥＬＳ）に対応する数のレイを有する。３次元の場合、投影データは、ＣＨＡＮＮＥＬＳ×セグメント数（ＳＥＧＭＥＮＴＳ）に対応する数のレイを有する。ＳＥＧＭＥＮＴＳは、追加の次元を示す。従って、総レイ数が、Ｎraysによって示される場合、２次元投影データは、Ｎrays＝ビュー数（ＶＩＥＷＳ）×ＣＨＡＮＮＥＬＳを有し、３次元投影データは、Ｎrays＝ＶＩＥＷＳ×ＳＥＧＭＥＮＴＳ×ＣＨＡＮＮＥＬＳを有する。投影データは、Ｎrays個の要素を有する列ベクトルｐによって示される。列ベクトルｐのｉ番目の要素は、ｉ番目のレイサムであるｐiによって示される。
【００３５】
システムマトリクスＡの次元は、Ｎrays×Ｎpixelsである。その行列要素は、ａi,jによって示される。従って、ｉ番目のレイサムは、（２）式で以下のように表わせる。
【数１】
【００３６】
システムマトリクスＡ内の全ての行列要素を記憶することは、その膨大な次元数のために膨大なコストがかかる。例えば、各ビューの８９６個のチャネルにおいて９００個のビューが測定され、サイズ５１２×５１２の画像行列が再構成される場合、システムマトリクスＡは、８０６，４００×２６２，１４４次元を有する。システムマトリクスＡは疎であるが、全ての行列要素を記憶するにはコストが高い。システムマトリクスＡを記憶しないようにする１つの方法は、転置を必要としない場合において、行列要素をその場で計算することである。１つの先行技術の手法において、１つのレイの行列要素だけが計算され記憶され、レイ毎に画像ボリュームが更新される。本実施形態においてシステムマトリクスＡは記憶されないが、その転置を逆投影演算子であると解釈する。
【００３７】
目的関数（費用関数）は、以下に示すように、２次元の場合に（３）式で規定され、３次元の場合に（４）式で規定されるようなＴＶである。変数εは、後述の探索方向の式によりシグマ内の各ｕ（・）項が分割されるので、各ｕ（・）項がゼロになることを防止するための微小量である。（３）式と（４）式とは、一連の離散勾配のｌlノルムを表わす。
【数２】
【００３８】
【数３】
【００３９】
図６は、図２のステップＳ２０に関する擬似コードのリストの一例を示す。この擬似コードが示す処理は、本実施形態に係るＴＶ−ＩＲにおけるＯＳ−ＳＡＲＴに含まれる。よく研究されたＯＳ―ＳＡＲＴを検討してきたが、本実施形態に係るＯＳ−ＳＡＲＴの例示的な実施態様は、先行技術で示されたようなＰＯＣＳの変形である。本実施形態は、ある特定のＯＳ−ＳＡＲＴが適正なＰＯＣＳ処理であることを証明するものではない。ＯＳ−ＳＡＲＴの利点は、必要に応じて並列で実施されることである。本実施形態に係るＯＳ−ＳＡＲＴにおいてシステムマトリクスＡは、キャッシュされない。
【００４０】
図６に示すように、ＯＳ−ＳＡＲＴにおいて投影データｐは、Ｎsets個のサブセットに区分される。ｔ番目のサブセットは、Ｎsubrays＝Ｎrays／Ｎsetsのレイサムを含むＳtによって示される。数学的には、ＯＳ−ＳＡＲＴの更新式は、以下の（５）式によって示される。
【数４】
【００４１】
ここで、システムマトリクスＡの行列要素と投影データｐとは、サブセットＳｔ内のレイに対応する。前述の式においてシステムマトリクスＡの転置が必要とされるが、本実施形態においては、システムマトリクスＡの行列要素を記憶せず、転置処理を逆投影演算子により実行する。（５）式中の３つの項は、以下の（５ａ）式、（５ｂ）式、（５ｃ）式のようにそれぞれ表される。
【数５】
【００４２】
第１項は、画像ｆ(k,t)の再投影である。
【数６】
【００４３】
第２項は、ｌ番目のレイの長さである。
【数７】
【００４４】
第３項は、ｊ番目の画素が逆投影された回数として表されてもよい。
【００４５】
図７は、図２のステップＳ３０に関する擬似コードのリストの一例を示す。この擬似コードが示す処理は、本実施形態に係るＴＶ−ＩＲにおける直線探索法に関する処理に含まれる。直線探索法は、ステップ幅αを求めるために利用される。ステップ幅αは、以下の（６）式又は図５の擬似コードの７．１行に示される更新式で利用される。
【数８】
【００４６】
その結果、ＵTV（ｆ(k+l)）≦ＵTV（ｆ(k)）＋μαｄTsearch▽ＵTV（ｆ(k)）となる。ここで、μは、０＜μ＜１を満たすスカラである。例えば、μ＝０．００１に設定される。
【００４７】
探索方向は、（３）式と（４）式とに示したような目的関数を最小化するために、勾配降下法において利用される。探索方向ｄsearchは、以下の（７）式のように規定される。
【数９】
【００４８】
（７）式は、Ｎpixels個の要素の列ベクトルとして表わされる。これは、２次元の場合（８）式により、３次元の場合（９）式により示される。
【数１０】
【００４９】
【数１１】
【００５０】
２次元の場合、（３）式の３つの項、すなわちｕ（ｋ−１，ｌ）、ｕ（ｋ，ｌ−１）、及びｕ（ｋ，ｌ）だけが変数ｆk,lを含む。これらの項は、（８）式を得るために、ｆk,lについてのみ微分される。
【００５１】
本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。
【符号の説明】
【００５２】
１００…ガントリ、１０１…Ｘ線管、１０２…フレーム、１０３…Ｘ線検出器、１０４…データ取得回路（ＤＡＳ）、１０５…非接触データ伝送部、１０６…前処理部、１０７…回転部、１０８…スリップリング、１０９…高電圧発生部、１１０…システム制御部、１１２…記憶部、１１４…再構成部、１１５…入力部、１１６…表示部
【特許請求の範囲】
【請求項１】
データ収集回路からの投影データを複数のサブセットに区分し、前記複数のサブセットの各々は複数のビューを含み、
前記投影データに基づく画像ボリュームを更新するために、前記複数のサブセットの各々について前記複数のビューにＯＳ−ＳＡＲＴを並列的に実行し、
前記ＯＳ−ＳＡＲＴの実行後、所定の規則に従って勾配ステップ幅を決定し、
前記決定された勾配ステップ幅を利用して全変動を適応的に最小化する、
ことを具備する画像処理方法。
【請求項２】
前記投影データは、多ビューを有する、請求項１記載の画像処理方法。
【請求項３】
前記投影データは、疎ビューを有する、請求項１記載の画像処理方法。
【請求項４】
前記画像ボリュームは、正規再構成により発生される、請求項２又は３記載の画像処理方法。
【請求項５】
前記画像ボリュームは、内部再構成により発生される、請求項２又は３記載の画像処理方法。
【請求項６】
前記勾配ステップ幅は、直線探索法に基づいて決定される、請求項１記載の画像処理方法。
【請求項７】
前記勾配ステップ幅は、反復して更新される前記画像ボリュームのうちの現在の画像ボリュームの目的関数が過去の画像ボリュームの目的関数よりも小さいことを保証する、請求項６記載の画像処理方法。
【請求項８】
前記ＯＳ−ＳＡＲＴは、グラフィック処理装置によって実行される、請求項１記載の画像処理方法。
【請求項９】
前記ＯＳ−ＳＡＲＴは、中央処理装置によって実行される、請求項１記載の画像処理方法。
【請求項１０】
前記ＯＳ−ＳＡＲＴを実行することにおいて、
前記複数のサブセットの各々について、前記画像ボリュームを再投影して計算上の投影データを生成し、
前記計算上の投影データと実測の投影データとの間の正規化差分を逆投影して、更新後の前記画像ボリュームを再構成する、
請求項１記載の画像処理方法。
【請求項１１】
投影データを収集する収集部と、
前記投影データに基づいて画像ボリュームを再構成する再構成部と、を具備するＸ線コンピュータ断層撮影装置であって、
前記再構成部は、前記投影データを各々が複数のビューを含む複数のサブセットに区分し、前記画像ボリュームを更新するために前記複数のサブセットの各々について前記複数のビューにＯＳ−ＳＡＲＴを並列的に実行し、前記ＯＳ−ＳＡＲＴの完了後において所定の規則に従って勾配ステップ幅を決定し、前記勾配ステップ幅を利用して全変動を適応的に最小化する、
ことを特徴とするＸ線コンピュータ断層撮影装置。
【請求項１２】
前記収集部は、前記投影データを多ビューで収集する、ことを特徴とする請求項１１記載のＸ線コンピュータ断層撮影装置。
【請求項１３】
前記収集部は、前記投影データを疎ビューで収集する、ことを特徴とする請求項１１記載のＸ線コンピュータ断層撮影装置。
【請求項１４】
前記収集部は、前記投影データを正規再構成のために収集する、ことを特徴とする請求項１２又は１３記載のＸ線コンピュータ断層撮影装置。
【請求項１５】
前記収集部は、前記投影データを内部再構成のために収集する、ことを特徴とする請求項１２又は１３記載のＸ線コンピュータ断層撮影装置。
【請求項１６】
前記再構成部は、前記勾配ステップ幅を直線探索法に基づいて決定する、ことを特徴とする請求項１１記載のＸ線コンピュータ断層撮影装置。
【請求項１７】
前記再構成部は、現在の画像ボリュームの目的関数が過去の画像ボリュームの目的関数より小さくなるように前記勾配ステップ幅を設定する、ことを特徴とする請求項１６記載のＸ線コンピュータ断層撮影装置。
【請求項１８】
前記再構成部は、グラフィック処理装置（ＧＰＵ）により構成される、ことを特徴とする請求項１１記載のＸ線コンピュータ断層撮影装置。
【請求項１９】
前記再構成部は、中央処理装置（ＣＰＵ）により構成される、ことを特徴とする請求項１１記載のＸ線コンピュータ断層撮影装置。
【請求項２０】
前記再構成部は、前記複数のサブセットの各々について、前記画像ボリュームを再投影して計算上の投影データを生成し、前記計算上の投影データと実測の投影データとの間の正規化差分を逆投影し、更新後の前記画像ボリュームを再構成する、ことを特徴とする請求項１１記載のＸ線コンピュータ断層撮影装置。
【請求項２１】
データ収集回路で収集された投影データを複数のサブセットに区分するステップであって、前記複数のサブセットの各々は、複数のビューを含む区分ステップと、
前記複数のサブセットのうちの１つのサブセットの前記複数のビューにＯＳ―ＳＡＲＴを並列的に実行するＯＳ―ＳＡＲＴステップと、
前記ＯＳ―ＳＡＲＴステップにおける画像ボリュームを更新する更新ステップと
前記複数のサブセットに前記ＯＳ―ＳＡＲＴステップと前記更新ステップとを反復実行する反復ステップと、
前記反復ステップの完了後、所定の規則に従って勾配ステップ幅を決定する決定ステップと、
前記決定された勾配ステップ幅を利用して全変動を適応的に最小化する最小化ステップと、
を具備する画像処理装置。
【請求項２２】
投影データを収集する収集部と、
前記投影データに基づいて画像ボリュームを再構成する再構成部と、を具備するＸ線コンピュータ断層撮影装置であって、
前記再構成部は、前記投影データを各々が複数のビューを含む複数のサブセットに区分し、前記複数のサブセットのうちの１つのサブセットの前記複数のビューにＯＳ―ＳＡＲＴを並列的に実行し、前記ＯＳ―ＳＡＲＴの完了後において画像ボリュームを更新し、前記複数のサブセットに前記ＯＳ―ＳＡＲＴと前記更新とを反復実行し、反復実行の完了後、所定の規則に従って勾配ステップ幅を決定し、前記決定された勾配ステップ幅を利用して全変動を適応的に最小化する、
ことを特徴とするＸ線コンピュータ断層撮影装置。
【請求項１】
データ収集回路からの投影データを複数のサブセットに区分し、前記複数のサブセットの各々は複数のビューを含み、
前記投影データに基づく画像ボリュームを更新するために、前記複数のサブセットの各々について前記複数のビューにＯＳ−ＳＡＲＴを並列的に実行し、
前記ＯＳ−ＳＡＲＴの実行後、所定の規則に従って勾配ステップ幅を決定し、
前記決定された勾配ステップ幅を利用して全変動を適応的に最小化する、
ことを具備する画像処理方法。
【請求項２】
前記投影データは、多ビューを有する、請求項１記載の画像処理方法。
【請求項３】
前記投影データは、疎ビューを有する、請求項１記載の画像処理方法。
【請求項４】
前記画像ボリュームは、正規再構成により発生される、請求項２又は３記載の画像処理方法。
【請求項５】
前記画像ボリュームは、内部再構成により発生される、請求項２又は３記載の画像処理方法。
【請求項６】
前記勾配ステップ幅は、直線探索法に基づいて決定される、請求項１記載の画像処理方法。
【請求項７】
前記勾配ステップ幅は、反復して更新される前記画像ボリュームのうちの現在の画像ボリュームの目的関数が過去の画像ボリュームの目的関数よりも小さいことを保証する、請求項６記載の画像処理方法。
【請求項８】
前記ＯＳ−ＳＡＲＴは、グラフィック処理装置によって実行される、請求項１記載の画像処理方法。
【請求項９】
前記ＯＳ−ＳＡＲＴは、中央処理装置によって実行される、請求項１記載の画像処理方法。
【請求項１０】
前記ＯＳ−ＳＡＲＴを実行することにおいて、
前記複数のサブセットの各々について、前記画像ボリュームを再投影して計算上の投影データを生成し、
前記計算上の投影データと実測の投影データとの間の正規化差分を逆投影して、更新後の前記画像ボリュームを再構成する、
請求項１記載の画像処理方法。
【請求項１１】
投影データを収集する収集部と、
前記投影データに基づいて画像ボリュームを再構成する再構成部と、を具備するＸ線コンピュータ断層撮影装置であって、
前記再構成部は、前記投影データを各々が複数のビューを含む複数のサブセットに区分し、前記画像ボリュームを更新するために前記複数のサブセットの各々について前記複数のビューにＯＳ−ＳＡＲＴを並列的に実行し、前記ＯＳ−ＳＡＲＴの完了後において所定の規則に従って勾配ステップ幅を決定し、前記勾配ステップ幅を利用して全変動を適応的に最小化する、
ことを特徴とするＸ線コンピュータ断層撮影装置。
【請求項１２】
前記収集部は、前記投影データを多ビューで収集する、ことを特徴とする請求項１１記載のＸ線コンピュータ断層撮影装置。
【請求項１３】
前記収集部は、前記投影データを疎ビューで収集する、ことを特徴とする請求項１１記載のＸ線コンピュータ断層撮影装置。
【請求項１４】
前記収集部は、前記投影データを正規再構成のために収集する、ことを特徴とする請求項１２又は１３記載のＸ線コンピュータ断層撮影装置。
【請求項１５】
前記収集部は、前記投影データを内部再構成のために収集する、ことを特徴とする請求項１２又は１３記載のＸ線コンピュータ断層撮影装置。
【請求項１６】
前記再構成部は、前記勾配ステップ幅を直線探索法に基づいて決定する、ことを特徴とする請求項１１記載のＸ線コンピュータ断層撮影装置。
【請求項１７】
前記再構成部は、現在の画像ボリュームの目的関数が過去の画像ボリュームの目的関数より小さくなるように前記勾配ステップ幅を設定する、ことを特徴とする請求項１６記載のＸ線コンピュータ断層撮影装置。
【請求項１８】
前記再構成部は、グラフィック処理装置（ＧＰＵ）により構成される、ことを特徴とする請求項１１記載のＸ線コンピュータ断層撮影装置。
【請求項１９】
前記再構成部は、中央処理装置（ＣＰＵ）により構成される、ことを特徴とする請求項１１記載のＸ線コンピュータ断層撮影装置。
【請求項２０】
前記再構成部は、前記複数のサブセットの各々について、前記画像ボリュームを再投影して計算上の投影データを生成し、前記計算上の投影データと実測の投影データとの間の正規化差分を逆投影し、更新後の前記画像ボリュームを再構成する、ことを特徴とする請求項１１記載のＸ線コンピュータ断層撮影装置。
【請求項２１】
データ収集回路で収集された投影データを複数のサブセットに区分するステップであって、前記複数のサブセットの各々は、複数のビューを含む区分ステップと、
前記複数のサブセットのうちの１つのサブセットの前記複数のビューにＯＳ―ＳＡＲＴを並列的に実行するＯＳ―ＳＡＲＴステップと、
前記ＯＳ―ＳＡＲＴステップにおける画像ボリュームを更新する更新ステップと
前記複数のサブセットに前記ＯＳ―ＳＡＲＴステップと前記更新ステップとを反復実行する反復ステップと、
前記反復ステップの完了後、所定の規則に従って勾配ステップ幅を決定する決定ステップと、
前記決定された勾配ステップ幅を利用して全変動を適応的に最小化する最小化ステップと、
を具備する画像処理装置。
【請求項２２】
投影データを収集する収集部と、
前記投影データに基づいて画像ボリュームを再構成する再構成部と、を具備するＸ線コンピュータ断層撮影装置であって、
前記再構成部は、前記投影データを各々が複数のビューを含む複数のサブセットに区分し、前記複数のサブセットのうちの１つのサブセットの前記複数のビューにＯＳ―ＳＡＲＴを並列的に実行し、前記ＯＳ―ＳＡＲＴの完了後において画像ボリュームを更新し、前記複数のサブセットに前記ＯＳ―ＳＡＲＴと前記更新とを反復実行し、反復実行の完了後、所定の規則に従って勾配ステップ幅を決定し、前記決定された勾配ステップ幅を利用して全変動を適応的に最小化する、
ことを特徴とするＸ線コンピュータ断層撮影装置。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【公開番号】特開２０１１−１３９８９４（Ｐ２０１１−１３９８９４Ａ）
【公開日】平成２３年７月２１日（２０１１．７．２１）
【国際特許分類】
【出願番号】特願２０１０−２５７６８８（Ｐ２０１０−２５７６８８）
【出願日】平成２２年１１月１８日（２０１０．１１．１８）
【出願人】（０００００３０７８）株式会社東芝 (54,554)
【出願人】（５９４１６４５４２）東芝メディカルシステムズ株式会社 (4,066)
【Ｆターム（参考）】
【公開日】平成２３年７月２１日（２０１１．７．２１）
【国際特許分類】
【出願日】平成２２年１１月１８日（２０１０．１１．１８）
【出願人】（０００００３０７８）株式会社東芝 (54,554)
【出願人】（５９４１６４５４２）東芝メディカルシステムズ株式会社 (4,066)
【Ｆターム（参考）】
[ Back to top ]