JP4887491B2 - 医用画像処理方法及びその装置、プログラム - Google Patents

Description

本発明は、被写体を輪切状に連続して撮影して得られる３次元表示可能な画像群を画像処理する画像処理方法に関し、特に、異なる時期の同一被写体の画像群を比較可能にする画像処理方法に関する。

２次元のみで利用される異なる時期の同一被写体の画像を比較可能にする技術が、特開２００５−１２２４８号公報に画像読影支援方法として開示されている。
この背景技術の画像読影支援方法は、互いに異なる撮影条件の下で同時期に同一の被写体を撮影して取得した２以上の画像からなる第１の画像群の各々の画像（第１の画像とする）について、互いに異なる撮影条件の下で第１の画像群の撮影とは異なる時期にその被写体を撮影して取得した１以上の画像からなる第２の画像群の各々の画像（第２の画像とする）との一致性を表す指標値を算出し、算出した指標値が所定の基準を満たす第１の画像と第２の画像との組合せを１つ抽出し、この組合せを構成する両画像に含まれる被写体の位置を合わせる位置合わせ処理を行い、さらに、位置合わせ処理後の第１の画像と第２の画像との差分による差分画像を生成する構成である。

この背景技術の画像読影支援方法によれば、第１の画像群及び第２の画像群の各画像は互いに異なる撮影条件の下で撮影され取得されたものであるから、被写体の姿勢や状態の異なる複数の画像を含んでいるため、撮影時に被写体の姿勢や状態が変化していたとしても、その変化による各々の状態で撮影された第１の画像群の各画像と第２の画像群の各画像との中からより適切な第１の画像と第２の画像との組合せを抽出することが可能になり、位置合わせ処理後の第１の画像と第２の画像とは、高い精度で位置合わせがなされたものとなり、診断者は、効率的に第１の画像と第２の画像との比較読影を行うことが可能になる。
特開２００５−１２２４８号公報

取り扱う画像は画像群であるものの、その画像群全体で被写体を立体的に表示することはできないものである。したがって、現在主流となりつつある立体表示可能な画像群に対して処理を実施することができないという課題を有する。

また、背景技術を用いてもっとも近似する画像同士で２次元空間上においてローカルマッチングを実施したとしても、被写体の姿勢の変化及び被写体の変形（被写体の臓器の変形を含む）がある場合には、経時差分画像上におけるアーチファクトが大きく生じるという課題を有する。

本発明は前記課題を解決するためになされたものであり、異なる時期の同一被写体の立体表示可能な画像群を比較可能にする画像処理装置を提供することを目的とする。また、アーチファクトの低減した経時差分画像を生成する画像処理装置を提供することも目的である。

本発明に係る医用画像処理方法は、第１の時期に被写体の断層をＣＴにより連続して撮影して立体表示可能な第１時期画像群と、第２の時期に被写体の断層をＣＴにより連続して撮影して立体表示可能な第２時期画像群とを対比するためのコンピュータを用いた医用画像処理方法であって、ＣＴ画像に撮影された前記被写体の注目部位画像の相互相関値を基準として、前記第１時期画像群を構成する第１時期ＣＴ画像と前記第２時期画像群を構成する第２時期ＣＴ画像との対応する組み合わせを求める組合せ演算工程と、前記組合せ演算工程で求めた前記第１時期ＣＴ画像と前記第２時期ＣＴ画像との対応する組み合わせにおいて、注目部位を基準とした画像上のズレをグローバルマッチングの第１シフトベクトルとして求める第１シフトベクトル演算工程と、前記第１時期画像群が形成する前記被写体の立体空間内で、所定の大きさの複数の空間領域をテンプレート空間領域として設定するテンプレート空間設定工程と、前記第２時期画像群が形成する前記被写体の立体空間内で、前記第１時期画像群における前記テンプレート空間領域と対応し、当該対応するテンプレート空間領域の中心座標と前記第１シフトベクトルから求まる、前記テンプレート空間領域より大きい探索空間領域を設定する探索空間設定工程と、前記探索空間領域内で対応する前記テンプレート空間領域を移動及び回転させながら、対応する前記テンプレート空間領域における前記第１時期画像群と第２時期画像群とのそれぞれのＣＴ値の相互相関値を求め、当該相互相関値が最も高い前記テンプレート空間領域の位置からローカルマッチングの第２シフトベクトルを求める第２シフトベクトル演算工程と、前記第１時期画像群内の第２シフトベクトルを求めていない座標を所定間隔で取り、当該座標の第２シフトベクトルを既に求めた第２シフトベクトルを用いて線形補間により求める線形補間工程と、前記第２シフトベクトルが演算された各位置において、前記第１時期画像群と前記第２時期画像群との相互相関値の負値を外部エネルギーとし、前記各位置において隣接する前記第２シフトベクトル間の滑らかさの度合いを内部エネルギーとし、前記各位置におけるエネルギーの総量が最小となるように、前記第２シフトベクトルを更新するシフトベクトル更新工程と、前記第１シフトベクトル及び前記第２シフトベクトルから前記第２時期画像群をワーピング処理するワーピング工程とを含み、前記第２シフトベクトル演算工程が、前記テンプレート空間領域におけるＸ軸、Ｙ軸及びＺ軸のそれぞれの方向の移動量及び回転角度を遺伝子とし、対応する前記第１時期画像群と第２時期画像群とのそれぞれのＣＴ値における相互相関値を適合度とする遺伝アルゴリズムを用いて第２シフトベクトルを演算するものである。

このように本発明によれば、立体表示可能な撮影時期の異なる画像群同士において被写体の注目部位画像の立体的構造を加味したローカルマッチングのシフトベクトルを求めることができ、被写体の姿勢の変化及び被写体の変形が立体的に生じたとしても第１時期画像群と第２時期画像群を使用者が適切に対比することができるという効果を奏する。
また、ローカルマッチングを求める場合にテンプレート空間の相互相関値だけでなく既に求められた近傍のシフトベクトルとの調和を加味するので、対比時のアーチファクトの低減を実現することができるという効果を奏する。
さらに、テンプレート領域又は仮想テンプレート領域の単純な平行移動だけではなく、テンプレート領域又は仮想テンプレート領域自体の回転をさせながら相関値を求めるので、より相関値の高いローカルマッチングのシフトベクトルを求めることができるという効果を有する。つまり、第１時期画像群と第２時期画像群をより合致させた条件で対比可能となり、不一致により生じるアーチファクトを低減することができるという効果を有する。平行移動の場合には対応する画素が常に存在する一方、回転を加えると回転角度により対応する画素が存在しない場合もあるが、この場合に近傍の画素の平均値を求め、その平均値を対応する画素の階調値として用いることで対応することができる。
さらにまた、移動パラメータ（前後、左右、上下）の他に回転パラメータ（ロール、ピッチ、ヨー）も加わり、ローカルマッチングのシフトベクトルの算出時に増加した計算量を、遺伝アルゴリズムを利用して大幅に低減することができるという効果を有する。遺伝子として回転パラメータだけでなく、移動パラメータも加えることができる。
対象テンプレート空間の中心座標を除く画素の全てについて、ローカルマッチングを求める場合に対象仮想テンプレート空間の相関値だけでなく既に求められた近傍のシフトベクトルとの調和を加味することもできるし、その一部の画素について適用することもできる。その一部の画素の典型例として、対象テンプレート空間の中心座標、すなわち、テンプレートＶＯＩの１次補間座標、２次補間座標、・・・、ｎ次補間座標を挙げることができる。ここで、対象テンプレート空間の中心座標においても適用することも可能である。
第１の時期が第２の時期よりも後でもよいし前であってもよい。後記する実施形態では第１の時期が第２の時期よりも後で説明している。

本発明に係る医用画像処理方法は必要に応じて、前記探索空間領域が使用者の要求に応じて変更可能であるものである。
このように本発明によれば、探索領域又は仮想探索領域を変更可能にしているので、探索領域又は仮想探索領域を拡大変更することでより相関値の高いローカルマッチングのシフトベクトルを求めることができ、他方、探索領域又は仮想探索領域を縮小変更することで処理に必要な計算量が低減することができるという効果を有する。

アーチファクトが生じている場合に、探索領域又は仮想探索領域の拡大を指示することでアーチファクトを低減する可能性がある。ＣＴスキャナーの装置間の相違、断層のスライス幅の相違、患者の姿勢の変化、病変部の変形により、最適な探索領域の大きさは異なる。最適な探索領域とは、アーチファクトの影響が可視レベルで殆どなくなり、且つ、できるだけ探索領域が小さい場合の探索領域である。

本発明に係る医用画像処理方法は必要に応じて、前記第１時期画像群とワーピング処理後の前記第２時期画像群とを同時に出力する工程を含むものである。
このように本発明によれば、ワーピング処理がなされた、つまり、第２時期画像群の画素移動が実行された後に、断層画像を適切に出力することができるという効果を有する。

本発明に係る医用画像処理方法は必要に応じて、前記第１時期画像群とワーピング処理後の前記第２時期画像群との差分画像を求めて出力する工程を含むものである。
このように本発明によれば、第１時期の断層画像との差分画像を求めるので、この差分画像から使用者が第１時期と第２時期との被写体自体の内部及び外部形状の変化を容易に認識することができるという効果を有する。

以上の発明は、方法以外に、装置及びプログラムとして把握することもできる。
これら前記の発明の概要は、本発明に必須となる特徴を列挙したものではなく、これら複数の特徴のサブコンビネーションも発明となり得る。

本発明の実施の形態を図面に基づいて詳細に説明する。本発明は多くの異なる形態で実施可能である。したがって、本実施形態の記載内容のみで解釈すべきではない。また、本実施の形態の全体を通して同じ要素には同じ符号を付けている。
本実施の形態では、主にシステムについて説明するが、当業者であれば明らかな通り、本発明はコンピュータで使用可能なプログラム又は方法としても実施できる。したがって、本発明は、ハードウェアとしての実施形態、ソフトウェアとしての実施形態またはソフトウェアとハードウェアとの組合せの実施形態をとることができる。プログラムは、ハードディスク、ＣＤ−ＲＯＭ、光記憶装置または磁気記憶装置等の任意のコンピュータ可読媒体に記録できる。

（本発明の第１の実施形態）
［１．システム構成］
図１は本実施形態に係る画像処理装置を含むシステム構成図である。
本実施形態に係るシステムは、被写体を高さ方向に所定間隔で撮影するＣＴスキャナー１０と、ＣＴスキャナー１０で生成された画像群を画像処理してディスプレイに出力するコンピュータ３０とからなる。ＣＴスキャナー１０とコンピュータ３０はＬＡＮ２０を介して接続している。

［２．ハードウェア構成］
図２は本実施形態に係る画像処理装置が構築されているコンピュータのハードウェア構成図である。
画像処理装置が構築されているコンピュータ３０は、ＣＰＵ(Central Processing Unit)３１１、ＲＡＭ(Random Access Memory)３１２、ＲＯＭ(Read Only Memory)３１３、フラッシュメモリ(Flash memory)３１４、外部記憶装置であるＨＤ(Hard disk)３１５、ＬＡＮ(Local Area Network)カード３１６、マウス３１７、キーボード３１８、ビデオカード３１９、このビデオカード３１９と電気的に接続する表示装置であるディスプレイ３１９ａ、３１９ｂ、１３９ｃ、サウンドカード３２０、このサウンドカード３２０と電気的に接続する音出力装置であるスピーカ３２０ａ及びフロッピーディスク（登録商標）、ＣＤ−ＲＯＭ、ＤＶＤ−ＲＯＭ等の記憶媒体を読み書きするドライブ３２１からなる。ここで、コンピュータ３０はディスプレイを３台有しており、１つのサウンドカード３２０と接続しているが、それぞれのディスプレイ毎にサウンドカードを有している構成であってもよい。

コンピュータ３０に画像処理プログラムがインストール、つまり、ＨＤ３１５に複製されメインメモリ上に読み出し可能な状態とされ、実際にメインメモリ上に画像処理プログラムが読み出されＣＰＵ３１１がその画像処理プログラムに従って動作することでコンピュータ３０は画像処理装置の動作を行う。

［３．画像処理装置］
図３は本実施形態に係る画像処理装置のブロック構成図である。
画像処理装置３０は、入力部３１、画像処理部３２及び出力部３３からなる。
入力部３１は、ＣＴスキャナー１０で生成された画像群を取り込む機能を有する。ある被写体に関してＣＴスキャナー１０で生成された画像群はＬＡＮ２０を介してコンピュータ３０に送信され、コンピュータ３０はＨＤ３１５に記録する。入力部１０は送信された画像群又はＨＤ３１５に記録された画像群から同一被写体に関する現在画像群及び過去画像群を取り込む。

コンピュータ３０にはＣＴスキャナー１０を制御する機能を有する構成にすることもでき、この場合には使用者がＣＴスキャナー１０に対して撮影命令を出力し、ＣＴスキャナー１０が指定された撮影条件で被写体を撮影し、画像群がコンピュータ３０に送信されることになる。
画像処理部３２は入力部３１により取り込まれた画像群を画像処理し、グローバルマッチング部３２ａ、ローカルマッチング部３２ｂ、ワーピング部３２ｃ及び差分画像生成部３２ｄからなる。

グローバルマッチング部３２ａは、肺野領域の相関を基準として現在画像群の現在画像と過去画像群の過去画像の組み合せ、現在画像の肺野領域を組み合わせられた過去画像内で移動させ、その中で最も相関が高かった位置を特定し、その位置を基準とした現在画像と対応する過去画像とのズレをグローバルマッチングのシフトベクトルとして求める機能を有する。現在又は過去に拘らず、画像群は同一被写体の所定間隔の断層画像からなる。撮影条件により現在画像群のｎ枚目が必ず過去画像群のｎ枚目に相当しないことから、現在画像群と過去画像群の中で対応する画像同士を見つける必要がある。ただし、このように現在画像と過去画像の組み合わせを見つける方法としては、使用者自ら指定する方法、最初の組み合わせを見つけてそれ以降はその組み合わせを基準にする方法をとることができる。

ローカルマッチング部３２ｂは、現在画像群からなる被写体の立体空間に所定間隔でテンプレートＶＯＩを複数設定し、このテンプレートＶＯＩを中心とする所定矩形空間（長方体、立方体）のテンプレート領域に対し、過去画像群からなる被写体の立体空間内のＣＴ値（画素値）に関して相関が高い領域を求め、領域同士のズレをローカルマッチングのシフトベクトルとして求める機能を有する。ローカルマッチング部３２ｂはグローバルマッチングのシフトベクトルを前提にしてローカルマッチングのシフトベクトルを求めている。

図４は本実施形態に係るローカルマッチングの説明図である。あるテンプレートＶＯＩを中心とするテンプレート領域が現在画像群からなる被写体の立体空間で設定される。また、このテンプレートＶＯＩにグローバルマッチングのシフトベクトルを加算して求めることができる対応する探索ＶＯＩを中心とするテンプレート領域よりも大きい探索領域も設定される。そして、設定されたテンプレート領域を探索領域内で移動させ、各地点の相互相関値を求め、移動させた中で最も高い相互相関値が求められた地点を求め、ローカルマッチングのシフトベクトルを求める。

この図４で説明した仕組みで各テンプレートＶＯＩのローカルマッチングのシフトベクトルを求めることができる。テンプレートＶＯＩ以外の点のローカルマッチングのシフトベクトルについては、シフトベクトルを用いて線形補間にて求めることができる。
ワーピング部３２ｃは、このテンプレートＶＯＩ以外の点のローカルマッチングのシフトベクトルを求め、グローバルマッチングのシフトベクトルとローカルマッチングのシフトベクトルから過去画像群を３次元にワーピングする機能を有する。

図５又は図６は本実施形態に係るelastic matching法による補間の説明図である。図５（ａ）は２次元に関するものであり、図５（ｂ）は３次元に関するものである。図５（ａ）は説明の便宜上図示している。テンプレートＶＯＩを元に１次補間点を求め、テンプレートＶＯＩ及び１次補間点を元に２次補間点を求めている。つまり、ｎ次元補間点は、テンプレートＶＯＩ及び（ｎ−１）次元以下の補間点から求めることができる。図５（ｂ）に示すように、３次元であっても同様である。なお、テンプレートＶＯＩから突出している白抜き矢印はテンプレートＶＯＩのローカルマッチングのシフトベクトルである。

図６（ａ）に示すように、本実施形態では注目部位の肺野領域内にテンプレートＶＯＩを設定している。前説したテンプレートＶＯＩ以外の点のローカルマッチングのシフトベクトルを求める場合には線形補間を使用するが、図６（ｂ）に示すように、線形補間を利用した場合には近傍であっても大きくベクトルの方向及び大きさが異なる。発明者の鋭意努力によりこの近傍点同士の相違とアーチファクトの発生が因果関係があることが判明している。これはローカルマッチングにおけるテンプレートＶＯＩのシフトベクトルから線形補間を行う場合、３次元空間における線形補間の柔軟性が欠如し、アーチファクトが生じることが原因の１つである。

したがって、アーチファクトの低減を図るため、全てのローカルマッチングのシフトベクトルに対するElastic matchingによる非線形補間を実行する。Elastic matchingにより、３次元過去−現在画像群間の相関、及び、ローカルシフトベクトル間の円滑さを同時に考慮でき、ローカルシフトベクトルの正確な補間が可能となる。

本実施形態では３次元情報を対象としてElastic matching法を適用し、シフトベクトルは３次元空間内に生じている変形の度合いを表している。
本実施形態で用いるElastic matching法では、画像間の相関を表す外部エネルギーと、シフトベクトル間の円滑さを表す内部エネルギーを、各補間点において求める。外部エネルギーは、補間を行うシフトベクトルにおける、現在画像と過去画像におけるテンプレートの相関値の負値で与える。したがって、相関値が高いほど、外部エネルギーは低い値をとる。一方、内部エネルギーは、シフトベクトルの１次、及び、２次微分値のノルムで与え、シフトベクトルが滑らかなほど、内部エネルギーは低い値をとることになる。以下に、外部エネルギー及び内部エネルギーの式をそれぞれ示す。

ただし、t(x, y, z)は過去画像における３次元テンプレートを、f(x, y, z)は現在画像における３次元テンプレートをそれぞれ示す。また、（１）式においてVはテンプレート画像の定義域を示す。さらに、（２）式においてt(dx, dy, dz)は過去画像におけるローカルシフトベクトルを示す。すなわち、外部エネルギーと内部エネルギーの重み付き線形和を、その補間点におけるローカルエネルギーとして、全ての補間点のローカルエネルギーの総量が最小となるように、補間シフトベクトルを決定する。
このようなElastic matching法を用いることで、図６（ｃ）に示すように、近傍で類似するベクトルを有するシフトベクトルを求めることができる。
差分画像生成部３２ｄは、現在画像とワーピングに得られた過去画像の差分をとり差分画像を生成する機能を有する。

［４．動作］
図７は本実施形態に係る画像処理装置の概要フローチャートである。
まず、現在画像群と過去画像群とを合わせるために過去画像群の拡大、縮小（ステップ１００）、グローバルマッチング（ステップ２００）、ローカルマッチング（ステップ３００）、ワーピング（ステップ４００）、差分処理（５００）を順に実行する。各処理について以下図面を用いて説明する。

図８は本実施形態に係るグローバルマッチングの詳細フローチャートである。
ＣＰＵ３１１はＧａｕｓｓｉａｎフィルタで現在画像群及び過去画像群をぼかす（ステップ２０１）。ＣＰＵ３１１は現在画像群から一の画像を取り出す（ステップ２１１）。ＣＰＵ３１１は対象現在画像から肺野を含むテンプレート画像を生成する（ステップ２２１）。ＣＰＵ３１１はテンプレート画像をクロージングする（ステップ２２２）。ここで、クロージング処理によってスムージングしている。

ＣＰＵ３１１は過去画像群から一の過去画像を取り出す（ステップ２３１）。ＣＰＵ３１１はテンプレート画像と対象過去画像から２次元相関値を求める（ステップ２４１）。ＣＰＵ３１１は対象過去画像のすべての領域で求めたか否かを判断し（ステップ２５１）、全て求めていない判断した場合にはＣＰＵ３１１は比較領域を変更（ステップ２６１）してステップ２４１に戻る。

前記ステップ２５１で全ての領域で求めたと判断した場合には、ＣＰＵ３１１は全ての過去画像を対象としたか否かを判断する（ステップ２７１）。全てを対象としていないと判断した場合にはステップ２３１に戻る。全てを対象としたと判断した場合には、ＣＰＵ３１１は最も相関値の高い過去画像を特定する（ステップ２８１）。ＣＰＵ３１１は対象現在画像と特定した過去画像を関連付けてそのときのシフトベクトルを記録する（ステップ２８２）。
ＣＰＵ３１１は全ての現在画像を処理したか否かを判断し（ステップ２９１）、全て処理していないと判断した場合にステップ２１１に戻り、全て処理したと判断した場合にはグローバルマッチングを終了する。

図９は本実施形態に係るローカルマッチングの詳細フローチャートである。
ＣＰＵ３１１は肺野領域の現画像群を取り出す（ステップ３０１）。ＣＰＵ３１１は肺野領域の現画像群に対してテンプレートＶＯＩを設定する（ステップ３０２）。
ＣＰＵ３１１は対象テンプレートＶＯＩを取り出す（ステップ３１１）。ＣＰＵ３１１は対象テンプレートＶＯＩにグローバルマッチングのシフトベクトルを用いて過去画像群での探索領域ＶＯＩを求める（ステップ３１２）。ＣＰＵ３１１は現在画像群の対象テンプレート領域と過去画像群の探索領域内の対象テンプレート対応領域との３次元相関値を求める（ステップ３１３）。ＣＰＵ３１１は探索領域内の全ての位置で３次元相関値を求めたか否かを判断する（ステップ３１４）。求めていないと判断した場合には、ＣＰＵ３１１は探索領域内で対象テンプレート領域をボクセル毎に移動させて（ステップ３１５）、ステップ３１３に戻る。

前記ステップ３１４で求めたと判断した場合には、ＣＰＵ３１１は現在対象としているテンプレート領域に関して求めた相関値の中で最も高い値であった位置を特定し、その位置からシフトベクトルを求める（ステップ３２１）。ＣＰＵ３１１は求めたシフトベクトルをテンプレートＶＯＩに関連付けてメインメモリ上に記録する（ステップ３２２）。
ＣＰＵ３１１は全てのテンプレートＶＯＩのシフトベクトルを求めたか否かを判断し（ステップ３３１）、求めていないと判断した場合にはステップ３１１に戻る。前記ステップ３３１で求めたと判断した場合には、ローカルマッチングを終了する。

図１０は本実施形態に係るワーピングの詳細フローチャートその１である。
ＣＰＵ３１１はテンプレートＶＯＩから１次補間点を求める（ステップ４０１）。ＣＰＵ３１１はテンプレートＶＯＩのシフトベクトルから線形補間法で１次補間点でのシフトベクトルを求める（ステップ４０２）。ＣＰＵ３１１はテンプレートＶＯＩ及び１次補間点を仮想テンプレートＶＯＩとする（ステップ４０３）。

ＣＰＵ３１１は仮想対象テンプレートＶＯＩを取り出す（ステップ４１１）。ＣＰＵ３１１は仮想対象テンプレートＶＯＩにグローバルマッチングのシフトベクトルを用いて過去画像群での仮想探索領域ＶＯＩを求める（ステップ４１２）。ＣＰＵ３１１は現在画像群の仮想対象テンプレート領域と過去画像群の仮想探索領域内の対象テンプレート対応領域とのエネルギー（外部エネルギー及び内部エネルギー）を求める（ステップ４１３）。ここで、内部エネルギーを求める場合には、近傍の仮想テンプレートＶＯＩのシフトベクトルを用いる。たとえば、仮想対象テンプレートＶＯＩの３軸方向に位置する仮想テンプレートＶＯＩのシフトベクトルを用いる。ＣＰＵ３１１は仮想探索領域内の全ての位置でエネルギーを求めたか否かを判断する（ステップ４１４）。求めていないと判断した場合には、ＣＰＵ３１１は仮想探索領域内で対象テンプレート領域をボクセル毎に移動させて（ステップ４１５）、ステップ４１３に戻る。

前記ステップ４１４で求めたと判断した場合には、ＣＰＵ３１１は現在対象としている仮想テンプレート領域に関して求めた相関値の中で最も高い値であった位置を特定し、その位置からシフトベクトルを求める（ステップ４２１）。ＣＰＵ３１１は求めたシフトベクトルを仮想テンプレートＶＯＩに関連付けてメインメモリ上に記録する（ステップ４２２）。

ＣＰＵ３１１は全ての仮想テンプレートＶＯＩのシフトベクトルを求めたか否かを判断し（ステップ４３１）、求めていないと判断した場合にはステップ４１１に戻る。求めていると判断した場合には、ＣＰＵ３１１は新たに更新したシフトベクトルをもつ仮想テンプレートＶＯＩ数が全ての仮想テンプレートの５[％]以下であるか否かを判断する（ステップ４３２）。以下でないと判断した場合にはステップ４１１に戻る。
ここまでの処理によりメインメモリ上には仮想テンプレートＶＯＩ、つまり、テンプレートＶＯＩ及び１次補間点のローカルシフトベクトルが記録されている。

図１１は本実施形態に係るワーピングの詳細フローチャートその２及び差分処理の詳細フローチャートである。テンプレートＶＯＩ及び１次補間点のローカルシフトベクトルは既に求まっているので、その他の点についてローカルシフトベクトルを求める。
ＣＰＵ３１１はテンプレートＶＯＩ及び（ｎ−１）次以下の補間点からｎ次補間点を求める（ステップ４４１）。ＣＰＵ３１１はテンプレートＶＯＩ及び（ｎ−１）次以下の補間点のローカルシフトベクトルから線形補間法でｎ次補間点のローカルシフトベクトルを求める（ステップ４４２）。ＣＰＵ３１１は求めたローカルシフトベクトルをｎ次補間点に関連付けてメインメモリに記録する（ステップ４４３）。ＣＰＵ３１１は現在画像群の全てのボクセルでシフトベクトルを求めたか否かを判断し（ステップ４４４）、求めていないと判断した場合にはＣＰＵ３１１はｎをインクリメント（ステップ４４５）してステップ４４１に戻る。

前記ステップ４４４で求めていると判断した場合にはＣＰＵ３１１は求めた各ボクセルにおけるグローバルシフトベクトル及びローカルシフトベクトルを用いて過去画像群をワーピングする（ステップ４５１）。
ＣＰＵ３１１は現在画像群及びワーピングした過去画像群から差分画像を求める（ステップ５０１）。ＣＰＵ３１１はディスプレイ３１９ａに差分画像を出力する（ステップ５１１）。各差分画像とともに現在画像及び過去画像をディスプレイ３１９ｂ、３１９ｃに出力することもできる。

［５．効果］
このように本実施形態に係る画像処理装置によれば、現在画像群を構成する現在画像と過去画像群を構成する過去画像とを単に最も相関する画像同士に対してグローバルマッチング、ローカルマッチング及びワーピングを実行して差分画像を出力するのではなく、現在画像群からなる立体物と過去画像群からなる立体物に対してローカルマッチングを行い、ワーピングして差分画像を出力するので、２軸ではなく３軸方向の動きに対応してアーチファクトを低減することができる。

（その他の実施形態）
［１．テンプレート領域の回転：Obtaining cross-correlation value employing rotation operator on the VOI+遺伝的アルゴリズム］
肺野領域内における，同一部位の経時的な位置ずれを補正するために、ローカルマッチングを行う。ローカルマッチングでは、具体的に肺野領域内に多数設置したＶＯＩ(Volume of Interest)における、シフトベクトルを求める。ここで、シフトベクトルとは、位置ずれを補正するための、過去画像の変形量(平行移動量)を示す。前記第１の実施形態では、経時画像間の同一部位の類似性を考慮するために、３次元相互相関値を計算し、シフトベクトルを決定する。始めに、現在肺野領域内にテンプレートＶＯＩを、過去肺野領域内に探索ＶＯＩを多数設置する。また隣り合うＶＯＩが半分ずつオーバーラップするように，等間隔にテンプレートＶＯＩを配置する。次に，テンプレートＶＯＩを用い、対応する探索ＶＯＩ内を１ボクセルずつ走査しながら３次元相互相関値を計測する。探索ＶＯＩの中心から相互相関値が最大となる点への移動量を、求めるシフトベクトルとする。このようにして前記第１の実施形態ではシフトベクトルを求めたが、テンプレートＶＯＩの走査時に、x，y，z軸に対するテンプレートの回転処理を行いながら相互相関値を計測することにより，より相互相関値の高いシフトベクトルの算出が可能となる（図１２参照）。

ここで、相互相関値の計測に回転処理を加えることにより，計算コストが格段に増加する。そこで，計算コストの低減を図るために，ＧＡ(Genetic Algorithm)を用いた探索法を提案する。ＧＡを用いた探索法のフローチャートを図１３に示す。遺伝子には(dx,dy,dz,dθ,dφ,dψ)を用い、適合度には３次元相互相関値を用いる。前記ステップ３１５、４１５の移動と同様に回転を加えて処理し、メインメモリの記録時にｘ、ｙ、ｚだけでなく、θ、φ、ψも記録する。回転を加えて３次元相互相関値を求める場合に回転後の位置に対応するＣＴ値がない場合があるが、その位置の周辺のＣＴ値を平均値を用いることとする。
具体的に、前記第１の実施形態に適用すると、ローカルマッチングへの適用の場合を図１４が示し、ワーピングへの適用の場合を図１５が示す。すなわち、探索領域に適用する場合と、仮想探索領域に適用する場合である。図１５では適合度にエネルギーを用いたが、３次元相互相関値であってもよい。

［２．探索領域変更：Expanded search region］
Elastic matching法において、求めるシフトベクトルのノルムの上限を決定する、探索領域のスケールの拡張を行う。探索領域のスケールを大きくすることにより、大きな変形を考慮できるが、計算コストが大きくなる。従来のスケールでは、左右肺野において、下葉部における経時的変形量が大きく、Elastic matching法により最適なシフトベクトルが得られなかった。そこで，従来の探索領域のスケールを拡張することにより（図１６参照）、この問題を解決する。ここで、探索領域だけでなく仮想探索領域を拡大してもよい。

図１７は前記第１の実施形態に適用した場合を示す。差分画像を出力した後に、探索領域を変更するか否かを判断する（ステップ８０１）。この前に利用者に対して探索領域の変更によりアーチファクトの低減が望める旨を報知することが望ましい。変更しないと判断した場合にはステップ８０１に戻る。変更すると判断した場合には現在の探索領域のサイズ、求めたローカルシフトベクトル、ワーピング後の過去画像及び差分画像をＨＤ３１５に記録する（ステップ８０２）。利用者により指定された探索領域のサイズで、前記ステップ３１３ないしステップ５０１を実行する（ステップ８１０）。求まった差分画像をステップ８０２で記録した差分画像と対比可能に出力し（ステップ８２１）、ステップ８０１に戻る。ここで、ステップ８０１に戻るのは対比結果によりアーチファクトの低減が達成できた場合にはさらに探索領域のサイズを変更して差分画像を求めるためである。

［３．適用範囲拡張：Expanded applying region］
Elastic matching法の適用範囲の拡張を行う。前記第１の実施形態では、テンプレートＶＯＩとその１次補間点において求められたシフトベクトルに対してのみ適用していたが、さらに２次補間点におけるシフトベクトルについても同様に処理を行うことにより、経時差分画像上に生じるアーチファクトの低減を図る。２次補間点以上のｎ次補間点までElastic matching法を適用することもできる。このとき、テンプレートＶＯＩからｎ次補間点全てにElastic matching法を適用してもよいし、部分的に適用してもよい。部分的にはとは、例えば、テンプレートＶＯＩ、３次補間点のみにElastic matching法を適用することである。

［４．分立体画像］
前記第１の実施形態においては差分画像を生成していたが、現在画像群からなる被写体の現在の立体空間と、ワーピング後の過去画像群からなる被写体の過去の立体空間とを画像ではなく直接空間同士として比較することもでき、複数の差分画像を差分画像に毎使用者が参照して経時的変化のある部位（異常組織）を発見するのではなく、立体的な視点で差分立体空間を参照することでより迅速に経時的変化のある部位を発見することができる。また、経時的変化のある部位を立体的に見ることができ、その大きさ、形状をより直感的に把握することができる。

そして、経時的変化のある部位を使用者が特定し、現在画像群での領域を顕示して表示することで、他部位と経時的変化のある部位とを比較して見ることができる。
また、現在と過去の２時点ではなく、他の複数の過去の時点の過去画像群を参照し、前記経時的変化のある部位をそれぞれの過去画像群での領域でも顕示して表示し、順次表示していくことで、経時的変化のある部位の経時的変化を把握することができる。

以上の前記各実施形態により本発明を説明したが、本発明の技術的範囲は実施形態に記載の範囲には限定されず、これら各実施形態に多様な変更又は改良を加えることが可能である。そして、かような変更又は改良を加えた実施の形態も本発明の技術的範囲に含まれる。このことは、特許請求の範囲及び課題を解決する手段からも明らかなことである。

まず、補間点における初期シフトベクトルを、ＶＯＩ上のシフトベクトルの線形補間により決定する。次に、初期シフトベクトルを用い、各補間点における初期ローカルエネルギーを求めておく。さらに各補間点においてシフトベクトルをＧｒｅｅｄｙアルゴリズム（D.J. William ans M.Shak, “A fast algorithm for active contours and curvature estimation”, Computer Vis. Graph, Image Process, Image Understand, 55, 14-26(1992)）により更新する。ここで、補間点におけるシフトベクトルは、３次元ベクトル(dx, dy, dz)と表すことができる。Ｇｒｅｅｄｙアルゴリズムにより、各補間点のシフトベクトルは、ローカルエネルギーが最小となるように、ベクトル検索領域（Ｎ×Ｎ×Ｎ）内のシフトベクトルに更新される。ここで、Ｎは計算コストとシフトベクトル値の変動を考慮し、実験的に求めた値である。
以上の処理を全ての補間点に適用し、シフトベクトルの更新を行う。更新する補間点数が、全補間点の５パーセント以下になるまでElastic matching法による補間を繰り返す。

最後に、ＶＯＩ及び１次補間点におけるシフトベクトルから、線形補間法にて過去ＣＴ画像全体のローカルシフトベクトルを算出する。得られた各ボクセルのローカルシフトベクトルとグローバルシフトベクトルを用い、過去画像を非線形に変形してワーピング画像を作成後、現在ＣＴ画像からワーピングされた過去ＣＴ画像を減算することにより経時的サブトラクションＣＴ画像を作成した。

図１８にローカルシフトベクトルの補間に線形補間を用いた場合と、Elastic matching法を用いた場合の経時的サブトラクションＣＴ画像の比較を示す。図１８において、左図は線形補間のみを用いた場合のサブトラクションＣＴ画像を、右図はElastic matching法により得られるサブトラクションＣＴ画像をそれぞれ示す。同図に示すように、Elstic matching法を用いることにより、血管部の位置ずれによるアーチファクトが軽減されていることがわかる。

［２．探索領域変更：Expanded search region］の結果を図１９に示す。左図は従来のスケールを示し、右図は新たに設定したスケールにより求められた差分画像を示す。丸印囲み部分を対比すると右図のがよりアーチファクトが低減していることが分かる。

［３．適用範囲拡張：Expanded applying region］の結果を図２０に示す。左図は１次補間点まで求めた場合を示し、右図は２次補間点まで求めた場合を示す。丸印囲み部分を対比すると右図がよりアーチファクトが低減していることが分かる。

本発明の第１の実施形態に係る画像処理装置を含むシステム構成図である。 本発明の第１の実施形態に係る画像処理装置が構築されているコンピュータのハードウェア構成図である。 本発明の第１の実施形態に係る画像処理装置のブロック構成図である。 本発明の第１の実施形態に係るローカルマッチングの説明図である。 本発明の第１の実施形態に係るelastic matching法による補間の説明図である。 本発明の第１の実施形態に係るelastic matching法による補間の説明図である。 本発明の第１の実施形態に係る画像処理装置の概要フローチャートである。 本発明の第１の実施形態に係るグローバルマッチングの詳細フローチャートである。 本発明の第１の実施形態に係るローカルマッチングの詳細フローチャートである。 本発明の第１の実施形態に係るワーピングの詳細フローチャートその１である。 本発明の第１の実施形態に係るワーピングの詳細フローチャートその２及び差分処理の詳細フローチャートである。 本発明のその他の実施形態に係るテンプレート領域の回転操作の説明図である。 本発明のその他の実施形態に係る遺伝アルゴリズムを適用したシフトベクトルの求め方の説明図である。 本発明のその他の実施形態に係る遺伝アルゴリズム適用時の詳細フローチャートである。 本発明のその他の実施形態に係る遺伝アルゴリズム適用時の詳細フローチャートである。 本発明のその他の実施形態に係る探索領域の拡大操作の説明図である。 本発明のその他の実施形態に係る探索領域拡大の詳細フローチャートである。 実施例の実験結果（第１の実施形態の構成）である。 実施例の実験結果（探索領域変更の構成）である。 実施例の実験結果（適用範囲拡張の構成）である。

符号の説明

１０ ＣＴスキャナー
２０ ＬＡＮ
３０ コンピュータ
３１ 入力部
３２ 画像処理部
３２ａ グローバルマッチング部
３２ｂ ローカルマッチング部
３２ｃ ワーピング部
３２ｄ 差分画像生成部
３３ 出力部
３１１ ＣＰＵ
３１２ ＲＡＭ
３１３ ＲＯＭ
３１４ フラッシュメモリ
３１５ ＨＤ
３１６ ＬＡＮカード
３１７ マウス
３１８ キーボード
３１９ ビデオカード
３１９ａ ディスプレイ
３１９ｂ ディスプレイ
３１９ｃ ディスプレイ
３２０ サウンドカード
３２０ａ スピーカ
３２１ ドライブ

Claims (6)

1. 第１の時期に被写体の断層をＣＴにより連続して撮影して立体表示可能な第１時期画像群と、第２の時期に被写体の断層をＣＴにより連続して撮影して立体表示可能な第２時期画像群とを対比するためのコンピュータを用いた医用画像処理方法であって、
   ＣＴ画像に撮影された前記被写体の注目部位画像の相互相関値を基準として、前記第１時期画像群を構成する第１時期ＣＴ画像と前記第２時期画像群を構成する第２時期ＣＴ画像との対応する組み合わせを求める組合せ演算工程と、
   前記組合せ演算工程で求めた前記第１時期ＣＴ画像と前記第２時期ＣＴ画像との対応する組み合わせにおいて、注目部位を基準とした画像上のズレをグローバルマッチングの第１シフトベクトルとして求める第１シフトベクトル演算工程と、
   前記第１時期画像群が形成する前記被写体の立体空間内で、所定の大きさの複数の空間領域をテンプレート空間領域として設定するテンプレート空間設定工程と、
   前記第２時期画像群が形成する前記被写体の立体空間内で、前記第１時期画像群における前記テンプレート空間領域と対応し、当該対応するテンプレート空間領域の中心座標と前記第１シフトベクトルから求まる、前記テンプレート空間領域より大きい探索空間領域を設定する探索空間設定工程と、
   前記探索空間領域内で対応する前記テンプレート空間領域を移動及び回転させながら、対応する前記テンプレート空間領域における前記第１時期画像群と第２時期画像群とのそれぞれのＣＴ値の相互相関値を求め、当該相互相関値が最も高い前記テンプレート空間領域の位置からローカルマッチングの第２シフトベクトルを求める第２シフトベクトル演算工程と、
   前記第１時期画像群内の第２シフトベクトルを求めていない座標を所定間隔で取り、当該座標の第２シフトベクトルを既に求めた第２シフトベクトルを用いて線形補間により求める線形補間工程と、
   前記第２シフトベクトルが演算された各位置において、前記第１時期画像群と前記第２時期画像群との相互相関値の負値を外部エネルギーとし、前記各位置において隣接する前記第２シフトベクトル間の滑らかさの度合いを内部エネルギーとし、前記各位置におけるエネルギーの総量が最小となるように、前記第２シフトベクトルを更新するシフトベクトル更新工程と、
   前記第１シフトベクトル及び前記第２シフトベクトルから前記第２時期画像群をワーピング処理するワーピング工程とを含み、
   前記第２シフトベクトル演算工程が、前記テンプレート空間領域におけるＸ軸、Ｙ軸及びＺ軸のそれぞれの方向の移動量及び回転角度を遺伝子とし、対応する前記第１時期画像群と第２時期画像群とのそれぞれのＣＴ値における相互相関値を適合度とする遺伝アルゴリズムを用いて第２シフトベクトルを演算することを特徴とする医用画像処理方法。
2. 前記探索空間領域が使用者の要求に応じて変更可能である
   前記請求項１に記載の医用画像処理方法。
3. 前記第１時期画像群とワーピング処理後の前記第２時期画像群とを同時に出力する工程を含む
   前記請求項１または２に記載の医用画像処理方法。
4. 前記第１時期画像群とワーピング処理後の前記第２時期画像群との差分画像を求めて出力する工程を含む
   前記請求項１または２に記載の医用画像処理方法。
5. 第１の時期に被写体の断層を連続して撮影して立体表示可能な第１時期画像群と、第２の時期に被写体の断層を連続して撮影して立体表示可能な第２時期画像群とを対比するための医用画像処理装置であって、
   ＣＴ画像に撮影された前記被写体の注目部位画像の相互相関値を基準として、前記第１時期画像群を構成する第１時期ＣＴ画像と前記第２時期画像群を構成する第２時期ＣＴ画像との対応する組み合わせを求める組合せ演算手段と、
   前記組合せ演算手段で求めた前記第１時期ＣＴ画像と前記第２時期ＣＴ画像との対応する組み合わせにおいて、注目部位を基準とした画像上のズレをグローバルマッチングの第１シフトベクトルとして求める第１シフトベクトル演算手段と、
   前記第１時期画像群が形成する前記被写体の立体空間内で、所定の大きさの複数の空間領域をテンプレート空間領域として設定するテンプレート空間設定手段と、
   前記第２時期画像群が形成する前記被写体の立体空間内で、前記第１時期画像群における前記テンプレート空間領域と対応し、当該対応するテンプレート空間領域の中心座標と前記第１シフトベクトルから求まる、前記テンプレート空間領域より大きい探索空間領域を設定する探索空間設定手段と、
   前記探索空間領域内で対応する前記テンプレート空間領域を移動及び回転させながら、対応する前記テンプレート空間領域における前記第１時期画像群と第２時期画像群とのそれぞれのＣＴ値の相互相関値を求め、当該相互相関値が最も高い前記テンプレート空間領域の位置からローカルマッチングの第２シフトベクトルを求める第２シフトベクトル演算手段と、
   前記第１時期画像群内の第２シフトベクトルを求めていない座標を所定間隔で取り、当該座標の第２シフトベクトルを既に求めた第２シフトベクトルを用いて線形補間により求める線形補間手段と、
   前記第２シフトベクトルが演算された各位置において、前記第１時期画像群と前記第２時期画像群との相互相関値の負値を外部エネルギーとし、前記各位置において隣接する前記第２シフトベクトル間の滑らかさの度合いを内部エネルギーとし、前記各位置におけるエネルギーの総量が最小となるように、前記第２シフトベクトルを更新するシフトベクトル更新手段と、
   前記第１シフトベクトル及び前記第２シフトベクトルから前記第２時期画像群をワーピング処理するワーピング手段とを備え、
   前記第２シフトベクトル演算手段が、前記テンプレート空間領域におけるＸ軸、Ｙ軸及びＺ軸のそれぞれの方向の移動量及び回転角度を遺伝子とし、対応する前記第１時期画像群と第２時期画像群とのそれぞれのＣＴ値における相互相関値を適合度とする遺伝アルゴリズムを用いて第２シフトベクトルを演算することを特徴とする医用画像処理装置。
6. 第１の時期に被写体の断層を連続して撮影して立体表示可能な第１時期画像群と、第２の時期に被写体の断層を連続して撮影して立体表示可能な第２時期画像群とを対比するための医用画像処理プログラムであって、
   コンピュータを、
   ＣＴ画像に撮影された前記被写体の注目部位画像の相互相関値を基準として、前記第１時期画像群を構成する第１時期ＣＴ画像と前記第２時期画像群を構成する第２時期ＣＴ画像との対応する組み合わせを求める組合せ演算手段と、
   前記組合せ演算手段で求めた前記第１時期ＣＴ画像と前記第２時期ＣＴ画像との対応する組み合わせにおいて、注目部位を基準とした画像上のズレをグローバルマッチングの第１シフトベクトルとして求める第１シフトベクトル演算手段と、
   前記第１時期画像群が形成する前記被写体の立体空間内で、所定の大きさの複数の空間領域をテンプレート空間領域として設定するテンプレート空間設定手段と、
   前記第２時期画像群が形成する前記被写体の立体空間内で、前記第１時期画像群における前記テンプレート空間領域と対応し、当該対応するテンプレート空間領域の中心座標と前記第１シフトベクトルから求まる、前記テンプレート空間領域より大きい探索空間領域を設定する探索空間設定手段と、
   前記探索空間領域内で対応する前記テンプレート空間領域を移動及び回転させながら、対応する前記テンプレート空間領域における前記第１時期画像群と第２時期画像群とのそれぞれのＣＴ値の相互相関値を求め、当該相互相関値が最も高い前記テンプレート空間領域の位置からローカルマッチングの第２シフトベクトルを求める第２シフトベクトル演算手段と、
   前記第１時期画像群内の第２シフトベクトルを求めていない座標を所定間隔で取り、当該座標の第２シフトベクトルを既に求めた第２シフトベクトルを用いて線形補間により求める線形補間手段と、
   前記第２シフトベクトルが演算された各位置において、前記第１時期画像群と前記第２時期画像群との相互相関値の負値を外部エネルギーとし、前記各位置において隣接する前記第２シフトベクトル間の滑らかさの度合いを内部エネルギーとし、前記各位置におけるエネルギーの総量が最小となるように、前記第２シフトベクトルを更新するシフトベクトル更新手段と、
   前記第１シフトベクトル及び前記第２シフトベクトルから前記第２時期画像群をワーピング処理するワーピング手段として機能させ、
   前記第２シフトベクトル演算手段が、前記テンプレート空間領域におけるＸ軸、Ｙ軸及びＺ軸のそれぞれの方向の移動量及び回転角度を遺伝子とし、対応する前記第１時期画像群と第２時期画像群とのそれぞれのＣＴ値における相互相関値を適合度とする遺伝アルゴリズムを用いて第２シフトベクトルを演算することを特徴とする医用画像処理プログラム。