JP4356831B2

JP4356831B2 - 複数画像間の非剛体レジストレーション方法

Info

Publication number: JP4356831B2
Application number: JP2003304747A
Authority: JP
Inventors: 潤桝本; 将吾畦元
Original assignee: Aze Ltd
Current assignee: Aze Ltd
Priority date: 2003-08-28
Filing date: 2003-08-28
Publication date: 2009-11-04
Anticipated expiration: 2023-08-28
Also published as: JP2005078176A

Description

本発明は主に、ＣＴ（computed tomography）・ＭＲＩ（magnetic resonance imaging）・核医学・ＣＲ（computed radiography）・ＤＳＡ（digital subtraction angiography）・ＤＲ（real-time radiography）等の放射線診断システムを用いて撮像した医療用の複数の画像間において、肝臓・膵臓・腎臓・脾臓・心臓等の臓器や肺・脳・血管あるいは腫瘍といった生体組織など、診断や観察の対象となる標的組織の位置合せ（各画像において標的組織の各部に対応する各画像点を各画像間で対応づける）を、標的組織の非剛体としての変形を考慮して行なう複数画像間の非剛体レジストレーション方法に関する。

従来、放射線診断システムを用いた標的組織の撮影は、外科的処置を伴わずに標的組織おける構造変化や機能障害、病変部の有無等を検知するために、臨床現場において広く行なわれている。また、３次元的に得られたＸ線ＣＴやＭＲＩ等による標的組織の画像データに基づき、標的組織の３次元表示画像や経時的に変化する動画像をコンピュータ・グラフィックスにより作成する技術は、近年著しい進歩を見せている。

例えば、腹部領域内の肝臓等において、造影剤を注入してから撮像されるまでの時間帯が互いに異なる（撮影時相の相異なる）複数の画像系列の３次元ＣＴ画像により得られた画像データに基づき肝臓等の領域を特定する技術が、本願発明者らにより提案されている（下記特許文献１参照）。このような技術を用いて作成される肝臓等に関する高精度な３次元表示画像や動画像は、肝臓等における疾病診断や外科手術計画に際し有効であり、臨床現場において広く用いられる傾向にある。

撮影時相の相異なる複数の画像系列は、１回の息止め期間によって得られない場合も多い。このような場合、相異なる画像系列間で互いに対応する画像の組において、標的組織の位置が正確に位置合せされている保証はなく、むしろ呼吸による体の動きにより各画像間で標的組織の位置ずれが起きている可能性が高い。このため、高精度な３次元表示画像や動画像を作成するためには、画像間での標的組織の位置合せが重要な課題となる。

従来、このような画像間での標的組織の位置合せを行なうための手法として、下記非特許文献１に記載されているような非剛体レジストレーションが知られている。この非剛体レジストレーションは、画像空間を所定の間隔ごとに区切る制御ポイントをパラメータとするレジストレーション評価関数を用いて、このレジストレーション評価関数が最大となるように制御ポイントの位置を制御することにより、画像空間の変形ベクトル場を推定するものであり、上記レジストレーション評価関数は、画像間の類似性の尺度に基づき画像間の対応づけをする項（以下「類似性尺度項」と称する）と、画像内の空間変動は滑らかであるという条件に基づき画像間の対応づけを制限する項（以下「平滑制限項」と称する）とから構成されている。

造影剤を用いたＣＴ・ＭＲＩ等により異なる時相間で撮像された生体組織の画像では、各生体組織の特性に応じて濃度値の分布（例えばコントラスト）がそれぞれに変化する。したがって、非剛体レジストレーションで用いられる上記類似性の尺度は、位置合せされる画像間の濃度値分布の相違に対処し得るものである必要がある。従来そのような類似性の尺度として、相互情報量（ＭＩ）や正規化相互情報（ＮＭＩ）に相当するエントロピー相関係数（ＥＣＣ）、ピアソンの積率相互相関（ＮＣＣ）等が知られている。これらの類似性の尺度は、相互に対応づけられる複数の画像の各々において各画像点が持つ濃度値を組にしてなる多次元ベクトルの同時分布に基づいて求められる。

特開２００２−３４５８０７号公報 Rueckert D Sonoda LI，Hayes C，et al：Nonrigid Registration Using Free-Form Deformations：Application to Breast MR Images，IEEE Transactions on Medical Imaging，vol.18，no.8，pp.712-721，1999. Forsey DR，Bartels RH：Hierarchical B-Spline refinement：ACM Transactions on Computer Graphics，vol.22，no.4，pp.205-212，1988.

しかし、上記ＭＩやＥＣＣ、ＮＣＣを類似性の尺度とする従来の非剛体レジストレーションは、画像間の位置合せ誤差が大きく、臨床現場において用いられる段階には達していない。このため、位置合せ誤差の少ない非剛体レジストレーションを可能とする新たな類似性の尺度の構築が待望されていた。

本発明は、このような事情に鑑みなされたものであり、画像間における標的組織の良好な位置合せを可能とする新たな類似性の尺度を用いた複数画像間の非剛体レジストレーション方法を提供することを第１の目的とする。

また、従来の非剛体レジストレーションでは、画像間の位置合せを各画像の撮像領域全体に対し実施するようになっている。このような手法は、上記同時分布の個人差が少なくて一般的な同時分布のモデルを構築し得る脳等の生体組織の場合には適用できるが、腹部における肝臓等の生体組織のように上記同時分布の個人差が大きくて一般的なモデルが構築しにくい場合は適用しにくい。

また、腹部領域を撮像した画像間の位置合せにおいては、隣接する複数の生体組織の境界で生じる生体組織相互間の位置ずれ（以下「組織滑動」と称する）が影響して、良好な位置合せが行なえないという問題がある。すなわち、上述したように非剛体レジストレーションで用いられるレジストレーション評価関数は、画像内の空間変動は滑らかであるという条件に基づき画像間の対応づけを制限する平滑制限項を備えているが、組織滑動が生じる領域では空間変動が不連続となる。このため、不連続な空間変動が生じる境界を横断するように平滑制限が適用されると、画像間の位置合せにおいて望ましくない挙動を示すという問題がある。

本発明は、このような事情に鑑み、組織滑動が生じる複数の組織を撮像した画像間において、標的組織を良好に位置合せ可能な複数画像間の非剛体レジストレーション方法を提供することを第２の目的とする。

上記第１の目的を達成するため、本発明に係る第１の複数画像間の非剛体レジストレーション方法は、
画像間で濃度値の分布が互いに異なるように所定の標的組織を撮像した複数の画像間において、前記標的組織の非剛体としての変形を考慮して該標的組織の位置合せを行なう複数画像間の非剛体レジストレーション方法において、
前記複数の画像の各々において前記標的組織に対応する各画像点が持つ濃度値を各標本値とする確率変数を前記画像別に設定するとともに、これら画像別の各確率変数を組としてなる多次元確率ベクトルの同時分布を推定し、
所定の局所領域内において前記各画像点を前記画像間で対応づける変換を設定し、
該変換によって前記各画像点が前記画像間で適正に対応づけられたとすれば、前記各画像点が前記各画像において持つ各濃度値を組としてなる多次元ベクトルは、前記同時分布に従うと仮定して、
前記多次元ベクトルが観測される尤度を示す、前記変換のパラメータを母数とする尤度関数を求め、
求められた前記尤度関数を前記各画像間の類似性の尺度として用いることを特徴とするものである。

また、上記第２の目的を達成するため、本発明に係る第２の複数画像間の非剛体レジストレーション方法は、
画像間で濃度値の分布が互いに異なるように所定の標的組織を他の組織と共に撮像した複数の画像間において、前記標的組織の非剛体としての変形を考慮して該標的組織の位置合せを行なう複数画像間の非剛体レジストレーション方法において、
前記複数の画像の各々において各画像点が持つ濃度値を各標本値とする確率変数を前記画像別に設定するとともに、これら画像別の各確率変数を組としてなる多次元確率ベクトルの同時分布を求め、
この同時分布に基づき、前記標的組織に略対応する対象分布範囲と他の非対象分布範囲とを選別し、
前記非対象分布範囲に対応する前記同時分布を均一化するとともに、
前記対象分布範囲に対応する前記同時分布を推定し、
所定の局所領域内において前記各画像点を前記画像間で対応づける変換を設定し、
該変換によって前記各画像点が前記画像間で適正に対応づけられたとすれば、前記各画像点が前記各画像において持つ各濃度値を組としてなる多次元ベクトルは、前記対象分布範囲における前記同時分布に従うと仮定して、
前記多次元ベクトルが観測される尤度を示す、前記変換のパラメータを母数とする尤度関数を求め、
求められた前記尤度関数を前記各画像間の類似性の尺度として用いることを特徴とするものである。

本発明において、前記複数の画像間で前記標的組織の位置合せを行なうことにより、該標的組織に対応する空間領域の変形ベクトル場を推定するようにしてもよい。

また、前記変換の前記パラメータは、Ｂ−スプラインの制御点の位置を示すものとすることができ、前記複数の画像は、Ｘ線ＣＴ撮影により得られたものとすることができる。

本発明に係る第１の複数画像間の非剛体レジストレーション方法によれば、上記構成を備えたことにより、複数画像間における標的組織の良好な位置合せが可能となる。

また、本発明に係る第２の複数画像間の非剛体レジストレーション方法によれば、上記構成を備えたことにより、画像間の位置合せにおける組織滑動の悪影響を低減でき、組織滑動が生じる複数の組織を撮像した画像間における標的組織の良好な位置合せが可能となる。

以下、本発明の一実施形態に係る複数画像間の非剛体レジストレーション方法（以下「本実施形態方法」と称する）について、図面を参照しながら説明する。なお、本実施形態方法は、造影剤を用いてＣＴ撮影された腹部領域における相異なる２時相の画像系列（造影剤を注入してから撮像されるまでの時間帯が相異なる画像系列）間において、肝臓を標的組織としてその位置合せを行なうものである。

本実施形態方法の基本的な手順では、まず上記２時相の画像系列で互いに対応する一対の画像Ｉ（χ)およびＪ（χ)（χは３次元座標ｘ，ｙ，ｚの組として表される位置ベクトルで、３次元画像において各格子点を構成するボクセル（ｖｏｘｅｌ）に相当する。χ＝（ｘ，ｙ，ｚ）^Ｔ、Ｔは転置を表す）において、各解剖学的構造（生体組織）に対応する各画像点の集合の族Γ＝｛γ_１，γ_２，γ_３，・・・，γ_ｎ｝を定義する。この一対の画像Ｉ（χ)およびＪ（χ)が正確に位置合せされている場合、所定の解剖学的構造に対応する各画像点の集合γ（γはΓに属する）は、画像Ｉ（χ)およびＪ（χ)において同時確率P(ｉ，ｊ｜γ)を伴って、それぞれ濃度値ｉ，ｊで現れると仮定する。換言すれば、一対の画像Ｉ（χ)およびＪ（χ)の各々において各解剖学的構造（γ）に対応する各画像点（χ）が持つ濃度値を各標本値とする確率変数（ＩまたはＪ）を画像別に設定するとともに、これら画像別の確率変数を組としてなる多次元確率ベクトル（本実施形態方法では２次元確率ベクトル（Ｉ，Ｊ））が所定の同時分布（Ｐ（Ｉ（χ），Ｊ（χ）｜γ））に従うと仮定する。図１にこのような同時分布の一例を示す。図１は一対の画像において各解剖学的構造に対応する同時分布の一例を示す図である。なお、図１に示す軸ｉは画像Ｉにおいて各画像が持つ濃度値を示し、軸ｊは画像Ｊにおいて各画像が持つ濃度値を示す。このことは他の図においても同様である。

次に、位置合せ（レジストレーション）されてない一対の画像Ｉ（χ)およびＪ（χ)について検討する。Ｊ（χ)の各画像点をＩ（χ)の各画像点に対応づけるレジストレーション変換Ｔ（χ)を設定し、この変換Ｔ（χ)によって上記各画像点が一対の画像Ｉ（χ)およびＪ（χ)間で適正に対応づけられたとすれば、前記各画像において持つ各濃度値を組としてなる２次元ベクトル（ｉ，ｊ）は上記同時分布に従うと仮定して、この２次元ベクトル（ｉ，ｊ）が観測される尤度を示す尤度関数を下式（１）により定義する。なお、このような形式の尤度関数は、統計解析の分野における尤度解析法、特に最尤推定法において一般的に用いられるものである。ただし最尤推定法では、母集団の分布における平均や共分散等を未知とし、これらを母数として尤度関数に組み込んでその最尤推定量を求めるのが一般的である。これに対し本発明では、母集団の分布は予め推定しておき（この推定に最尤推定法を適用することも可能）、レジストレーション変換Ｔ（χ)のパラメータ（本実施形態方法では後述するφ）を母数として尤度関数に組み込んでその最尤推定量を求める点において、最尤推定法における尤度関数の用い方とは異なっている。

また、最大対数尤度変換T_MLを、下式（２）により定義する。

脳のように、周辺組織との間の組織滑動が生じ難く、また同時分布の個人差が少ないため同時分布の一般的なモデルを構築し得る組織を標的組織とする場合、上式（１），（２）により定義される尤度関数は、画像間の位置合せを各画像の撮像領域全体に対し実施する非剛体レジストレーションにおける類似性の尺度として用いることが有効である。

しかし、本実施形態方法が標的組織とする肝臓のように、周辺組織との間の組織滑動が生じやすく、また同時分布の個人差が大きいため同時分布の一般的なモデルが構築しにくい場合は、次に概略的に示す２つの手順を必要とする。

１．位置合せされていない患者個々のＣＴ画像データから、標的組織である肝臓の同時分布を推定する。

２．他の組織の同時分布は類似性の尺度に影響を与えないように基本的に均一分布としてモデル化する。

このような２つの手順を実施することにより、肝臓に対応する画像領域内では位置合せが行なわれるが、他の周辺領域については位置合せが強制されない。この結果、肝臓と他の周辺組織との境界で発生する組織滑動による悪影響を低減することが可能となる。

以下、上記２つの手順の詳細を述べるが、その準備として、解剖学的構造の集合族Γは2つの組織、すなわち標的組織としての肝臓（Ｌ）と他の周辺組織（Ｏ）（Ｏは標的組織以外の全組織を示す）とからなるものとし、上式（１）を下式（３）のように変形しておく。

本実施形態方法で提案される画像間の類似性の尺度は、上式（３）の尤度関数においてＰ（ｉ，ｊ｜Ｌ），Ｐ（ｉ，ｊ｜Ｏ），Ｐ（Ｌ）およびＰ（Ｏ）を定義することによって得られる。

〈排他的条件〉
Ｐ（ｉ，ｊ｜Ｌ），Ｐ（ｉ，ｊ｜Ｏ），Ｐ（Ｌ）およびＰ（Ｏ）を定義する前に、本実施形態方法において導入する排他的条件について説明する。図２は排他的条件を示す図で、同図（ａ）は組織γ_３が排他的条件を満たさない場合を、同図（ｂ）は組織γ_３が排他的条件を満たさす場合をそれぞれ例示している。

組織γ_ｓに対応する２次元ベクトル（ｉ，ｊ）の同時分布Ｐ（ｉ，ｊ｜γ_ｓ）のｉ軸およびｊ軸への投影であるＰ（ｉ｜γ_ｓ）およびＰ（ｊ｜γ_ｓ）が、他の全ての組織γ_ｔ（ｔとｓは異なる）についての同時分布Ｐ（ｉ｜γ_ｔ）およびＰ（ｊ｜γ_ｔ）とそれぞれ無視できる程度しか重ならない場合に、Ｐ（ｉ，ｊ｜γ_ｓ）は排他的であると定義する。図２（ａ）において組織γ_３は、他の組織｛γ_４，γ_５，γ_６，γ_７，γ_８｝の分布領域が自身に関連した対象分布領域（斜線を付けた十字状の領域）と重なっているために排他的とはみなされないが、図２（ｂ）において組織γ_３は排他的であるとみなされる。数学的には、Ｐ（ｉ，ｊ｜γ_ｓ）は下式（４）の条件を満たすと排他的とみなされる。

ただし、εおよびρはそれぞれ０および１に十分近い値とする。

本実施形態方法では、標的組織すなわち肝臓（Ｌ）の同時分布が排他的条件を満たすと仮定する。実際の状況で排他的条件を課すのは制約が強すぎると思われるかもしれないが、非剛体レジストレーションについては、画像全体において排他的条件を満たすことは不要である。非剛体レジストレーションにおけるレジストレーション評価関数の計算は、画像全体ではなく各局所領域で行なわれるからである。したがって、この条件は局所領域のみで満たされればよく、これは実際の状況において満たすのに妥当な条件とみなすことが可能である。

図３は局所的排他的条件を示す図で、同図（ａ）は局所領域（Ａｒｅａ１〜５）を規定するための同時分布（同時ヒストグラム）を示し、同図（ｂ）は画像全域の同時分布を示している。また、同図（ａ）で規定された局所領域（Ａｒｅａ１〜５）の同時分布を同図（ｃ）〜（ｇ）にそれぞれ示す。組織γ_３は、画像全域の同時分布においては排他的とみなされないが、図３（ｄ）以外の各局所領域においては排他的であるとみなせる。排他的条件は、一対の画像Ｉ（χ)およびＪ（χ)において、標的組織（肝臓）の分布が周辺組織の分布からよく分離されることを保証する。これは、一対の画像Ｉ（χ)およびＪ（χ)の双方において標的および周辺組織領域の間の画像境界が検出可能であることを意味する。

以下、上式（３）におけるＰ（ｉ，ｊ｜Ｌ），Ｐ（ｉ，ｊ｜Ｏ），Ｐ（Ｌ）およびＰ（Ｏ）を推定する手順を説明する。まず、Ｐ（ｉ，ｊ｜Ｌ）を非剛体レジストレーションが行なわれる前の画像対から推定する手順について説明し、次にＰ（ｉ，ｊ｜Ｏ）をＰ（ｉ，ｊ｜Ｌ）が排他的条件を満たすように導出する手順について説明する。最後に、Ｐ（ｉ，ｊ｜Ｌ）およびＰ（ｉ，ｊ｜Ｏ）に基づき、Ｐ（Ｌ）およびＰ（Ｏ）を推定する手順について説明する。

〈Ｐ（ｉ，ｊ｜Ｌ）の推定〉
まず、標的組織（肝臓）の同時分布Ｐ（ｉ，ｊ｜Ｌ）は、患者個々のＣＴ画像データ間で大きな変動が見られるが、下式（５）で表される多変量正規分布関数（本実施形態方法では２変量正規分布関数）でよく近似されるものと仮定する。

次に、ＣＴ撮影により異なる時相で得た位置合せされていない２つの画像において、腹部ＣＴスキャンのための視野（ＦＯＶ）を脊椎位置に基づいて設定するとともに、肝臓組織に大部分占有されるように関心領域（ＶＯＩ）を設定する（図４（ａ）参照）。図４は肝臓領域の同時分布の推定手順を示す図で、同図（ａ）は関心領域を示し、同図（ｂ）は関心領域における同時分布を示している。また、同図（ｃ）は上式（５）に基づき推定されたＰ（ｉ，ｊ｜Ｌ）を示し、同図（ｄ）は後述の如く均一化された非対象領域の同時分布Ｐ（ｉ，ｊ｜Ｏ）を示している。なお、後述する大阪大学病院で得られた腹部ＣＴ画像データにおいては、各例で脊椎に対する肝臓の位置が大きく異ならないため、各患者についてＶＯＩの位置を固定することができることが確認された。

次いで、上式（５）における同時分布Ｐ（ｉ，ｊ｜Ｌ）の平均値および共分散行列を、２つの領域のＶＯＩから得た同時ヒストグラムを分析することによって推定する。平均値および共分散行列は、中心が同時ヒストグラムの最頻値にあって水平および垂直方向の幅がそれぞれｉ軸およびｊ軸に投影された１次元ヒストグラムの半値全幅（ＦＷＨＭ）の３倍に等しいヒストグラム領域から推定される。２つの領域はこの時点では位置合せされないが、同時分布の良好な近似を与える。図４（ｂ）および（ｃ）に、ＶＯＩから得た同時ヒストグラムおよび肝臓の推定同時分布をそれぞれ示す。なお、このような手法は、上記特許文献１において詳しく開示されている。

〈Ｐ（ｉ，ｊ｜Ｏ）の導出〉
周辺組織の同時分布Ｐ（ｉ，ｊ｜Ｏ）は、肝臓の同時分布Ｐ（ｉ，ｊ｜Ｌ）が排他的条件を満たすという前提にたって、基本的に均一な分布としてモデル化する。Ｐ（ｉ，ｊ｜Ｏ）は図２の斜線領域を除いて均一にモデル化すべきである。本実施形態方法では、Ｐ（ｉ，ｊ｜Ｌ）のｉ軸およびｊ軸それぞれへの投影Ｐ（ｉ｜Ｌ）およびＰ（ｊ｜Ｌ）に基づいて、図２の斜線領域の分布をモデル化する。これは、Ｐ（ｉ｜Ｌ）およびＰ（ｊ｜Ｌ）の正規化後、均一分布からＰ（ｉ｜Ｌ）およびＰ（ｊ｜Ｌ）の逆投影を引き算することによって実現される。

すなわち、まずＧ（ｉ｜Ｌ）およびＧ（ｊ｜Ｌ）をそれぞれ下式（６）および（７）により定義する。

このＧ（ｉ｜Ｌ）およびＧ（ｊ｜Ｌ）を用いてＰ（ｉ，ｊ｜Ｏ）は、下式（８）により均一にモデル化される。

〈Ｐ（Ｌ）およびＰ（Ｏ）の推定〉
組織ＬおよびＯの事前確率は基本的にそれらの体積に依存する。非剛体レジストレーションにおいて事前確率は尤度が計算される局所領域におけるそれらの比に依存すべきである。この比を正確に推定するのは容易ではないため、以下の経験的な方法を用いて事前確率を決定する。Ｐ（Ｌ）およびＰ（Ｏ）をそれぞれα（αは０以上１以下の数）および１−αとする。Ｐ（Ｌ）・Ｐ（ｉ，ｊ｜Ｌ）の最大確率がＰ（Ｏ）・Ｐ（ｉ，ｊ｜Ｏ）の最大確率に等しいとすると、αへの制限が下式（１０）によって得られる。

この制限から、αが下式（１１）により導き出される。

〈非剛体レジストレーション手順への類似性尺度の組み込み〉
上述した手順により定義された上式（３）の尤度関数を類似性の尺度として、下式（１２）のように定義されるレジストレーション評価関数（目的関数）中に組み込む。

ただし、Ｃ_{ｓｉｍｉｌａｒｉｔｙ}(φ)は上述の尤度として定義された類似性尺度項を表し、Ｃ_{ｓｍｏｏｔｈ}(φ)は平滑制限項を表し、λは２つの項のバランスを取る重みづけのパラメータである。また、Φは下式（１３）で与えられるレジストレーション変換を記述するパラメータを表す。

ただし、δ（ｘ；Φ）はＢ−スプライン（曲面）で記述される自由曲面の変形（ＦＦＤ）を表し、ΦはＢ−スプライン制御点の全体集合を表し、φ（φはΦに属する）は各局所領域に関係する制御点の部分集合を表す。

Ｃ_{ｓｉｍｉｌａｒｉｔｙ}(φ)およびＣ_{ｓｍｏｏｔｈ}(φ)はそれぞれ各局所領域における類似性尺度項および平滑制限項であり、下式（１４）および（１５）によって与えられる。

ただし、Ｖはφが関係する局所領域を表す。

上式（１２）で定義されるレジストレーション評価関数を最小化するため、本実施形態方法では階層グリッドを用いる最急降下アルゴリズムを使用する。このアルゴリズム手順は以下の通りである。
〈１〉制御点Φ^ｍ（ｍ＝０）を初期化する。
〈２〉上式(１２)の勾配ベクトルを下式（１６）に基づいて計算する。

〈３〉‖∇Ｃ‖＞ε（εは小さな正の数）の間において、制御点Φを下式（１７）に基づいて更新し、勾配ベクトル▽Ｃを再計算する。

〈４〉制御点Φの階層を細かくする（Φ^ｍ→Φ^ｍ+１）
〈５〉上記手順〈２〉〜〈４〉を最高階層度に達するまで繰り返す。
ただし、Φ^ｍはｍ番目の階層レベルのＢ−スプライン制御点グリッドを示す。

以下、本実施形態方法を、臨床現場で得られたＣＴ画像データに対して実施した結果について説明する。適用対象となったＣＴ画像データは、大阪大学病院（大阪府吹田市）および国立癌センター（東京都）において肝臓のダイナミックＣＴ撮影により得られた８患者のデータである。撮影条件を表１に示す。なお、上記アルゴリズム手順におけるΦ^ｍ→Φ^ｍ+１の演算は、本明細書の背景技術の欄に記載した非特許文献２において提案されている手法を用いている。また、階層グリッドは３つのレベルからなり、グリッド間隔は４２mm、２１mm、１０．５mmとした。

位置合せされる相異なる２つの時相の画像系列は、一方（第１時相）が造影剤注入前の時相あるいは動脈相第１相（造影剤の効果が小さい）のものであり、他方（第２時相）が門脈相のものである（造影剤の効果が大きい）。このような２つの時相における３次元ＣＴ画像データは１回の息止め期間内で得ることができず、呼吸運動により２つの時相の画像間において各組織の変形がかなり認められる。なお、得られた３次元ＣＴ画像データのサイズは、ボクセル単位で５１２×５１２×１５０〜２００である。

位置合せ結果を評価するため、専門家により両時相の画像から抽出された肝臓の輪郭線を用いた。そして、位置合せ誤差を、位置合せ前の画像Ｉ（χ）における肝臓の輪郭線と位置合せ後の画像Ｊ（Ｔ（χ））における肝臓の輪郭線との距離の平均値によって定義した。また、位置合せにより推定された、肝臓に対応する空間領域の変形ベクトル場の妥当性を視覚的に評価した。

図５に８つのＣＴ画像データに対する位置合せの評価結果を示す。図５に示すように、本実施形態方法における位置合せ誤差は、エントロピー相関係数（ＥＣＣ）およびピアソンの積率相互相関（ＮＣＣ）を類似性の尺度として用いた場合に比べ少ないことが確かめられた。

図６に本実施形態方法およびＮＣＣを用いる従来手法によりそれぞれ得られた変形ベクトル場を示す。図６は推定された変形ベクトル場を示す図で、同図（ａ）は従来手法により得られた変形ベクトル場を示し、同図（ｂ）は本実施形態方法により得られた変形ベクトル場を示している。このケースでは肋骨と肝臓との間で組織滑動が生じていた。すなわち、図７に示すように、呼吸運動により左側の肋骨が上向きに変位し右側の肋骨が動かない状態なのに対し、肝臓は下向きに変位していた。図７は腹部内の各生体組織の変位を示す図である。

図６（ａ）に示すように、従来手法では肋骨および肝臓を共に位置合せしようとするので、肝臓の左側領域が左側の肋骨の上向きの動きによる悪影響を受け、肝臓の右下側の領域が動かない右側の肋骨の悪影響を受けている。図６（ｂ）に示すように、本実施形態方法では肋骨（骨組織）に関しては位置合せされないため、腹部における組織滑動による影響を受けないことが確認された。

以上、本発明の実施形態について説明したが、本発明は上記実施形態に限定されるものではなく、種々に態様を変更することが可能である。

例えば上記実施形態では、造影剤を用いたＣＴ撮影により得られた腹部領域の画像データに基づき、肝臓を標的組織としてその位置合せを実施しているが、撮影手法はＣＴに限られずまた標的組織は肝臓に限られない。撮影手法としてはＣＴ以外に、ＭＲＩ・核医学・ＣＲ・ＤＳＡ・ＤＲ等の種々の放射線診断システムを用いることが可能である。他の撮影手法で得られた画像データに本発明を適用するのに際して特段の困難はない。濃度値の分布が相異なる複数の画像が得られさえすれば、上述した本発明の手順を撮影手法に関係なく略同様に適用することが可能である。

なお、画像間で濃度値分布を変えるため、上記実施形態では造影剤を用いているが、造影剤を用いることは本発明にとって必須の要件ではない。造影剤を用いずＸ線の照射強度や磁場強度などの撮影条件を変えることによって、画像間の濃度値分布を変えるようにしてもよい。また、位置合せされる画像間の撮影時刻の隔たりは、造影剤注前後といった短い間隔に限られない。撮影された日時が、数日・数ヶ月・数年といった単位で異なる画像間の位置合せに対しても本発明は適用可能である。

標的組織に関しては、肝臓以外に膵臓・腎臓・脾臓・心臓等の臓器や肺・脳・血管あるいは腫瘍といった他の種々の生体組織を、その対象とすることが可能である。本発明を肺野に適用した場合、位置合せにより得られた変形ベクトル場から、肺の呼吸による体積変化を推定することが可能となる。この点に関する技術は、本発明者等により特許庁に対して既に開示がなされている（特願２００３−１７２５５０号明細書）。他の組織を標的組織とした場合に本発明を適用するのに際して特段の困難はなく、上述した本発明の手順を標的組織の種別に関係なく略同様に適用することが可能である。

また、上記実施形態では２つの画像系列間での位置合せを実施しているが、本発明は３つ以上の多数の画像間での位置合せにも適用することが可能である。３つ以上の画像間の位置合せに本発明を適用する場合、位置合せの基準となる画像を選択しておき、他の画像を１つずつ基準の画像に対して位置合せする方法をとることができる。その場合、上述した本発明の手順を繰り返し行なえばよい。３つ以上の画像間の位置合せを一度に実施する場合は、上述した本発明の手順において、２次元の同時分布を３次元以上の多次元同時分布に置き換えたり、２変量正規分布を３変量以上の多変量正規分布に置き換えたりするなどの次元を拡張する処理を行なえばよい。

また、上記実施形態では自由曲面を表現するモデリング手法としてＢ−スプラインを用いているが、Ｂ−スプラインの手法に代えて、自由曲面を表現することが可能なベツィエの手法等その他のモデリング手法を用いることも可能である。

また、上記実施形態ではレジストレーション評価関数を最小化するためのアルゴリズムとして最急降下法を用いているが、最急降下法に代えて共役勾配法やニュートン・ラフソン法、準ニュートン法、レーベンベルグ・マルカート法等の種々の最適化法を用いることも可能である。

一対の画像において各生体組織に対応する同時分布の一例を示す図排他的条件を示す図局所的排他的条件を示す図肝臓領域の同時分布の推定手順を示す図８つのＣＴ画像データに対する位置合せの評価結果を示す図推定された変形ベクトル場を示す図腹部内の各生体組織の変位を示す図

Claims

画像間で濃度値の分布が互いに異なるように所定の標的組織を撮像した複数の画像間において、前記標的組織の非剛体としての変形を考慮して該標的組織の位置合せを行なう複数画像間の非剛体レジストレーション方法において、
前記複数の画像の各々において前記標的組織に対応する各画像点が持つ濃度値を各標本値とする確率変数を前記画像別に設定するとともに、これら画像別の各確率変数を組としてなる多次元確率ベクトルの同時分布を推定し、
所定の局所領域内において前記各画像点を前記画像間で対応づける変換を設定し、
該変換によって前記各画像点が前記画像間で適正に対応づけられたとすれば、前記各画像点が前記各画像において持つ各濃度値を組としてなる多次元ベクトルは、前記同時分布に従うと仮定して、
前記多次元ベクトルが観測される尤度を示す、前記変換のパラメータを母数とする尤度関数を求め、
求められた前記尤度関数を前記各画像間の類似性の尺度として用いることを特徴とする複数画像間の非剛体レジストレーション方法。
画像間で濃度値の分布が互いに異なるように所定の標的組織を他の組織と共に撮像した複数の画像間において、前記標的組織の非剛体としての変形を考慮して該標的組織の位置合せを行なう複数画像間の非剛体レジストレーション方法において、
前記複数の画像の各々において各画像点が持つ濃度値を各標本値とする確率変数を前記画像別に設定するとともに、これら画像別の各確率変数を組としてなる多次元確率ベクトルの同時分布を求め、
この同時分布に基づき、前記標的組織に略対応する対象分布範囲と他の非対象分布範囲とを選別し、
前記非対象分布範囲に対応する前記同時分布を均一化するとともに、
前記対象分布範囲に対応する前記同時分布を推定し、
所定の局所領域内において前記各画像点を前記画像間で対応づける変換を設定し、
該変換によって前記各画像点が前記画像間で適正に対応づけられたとすれば、前記各画像点が前記各画像において持つ各濃度値を組としてなる多次元ベクトルは、前記対象分布範囲における前記同時分布に従うと仮定して、
前記多次元ベクトルが観測される尤度を示す、前記変換のパラメータを母数とする尤度関数を求め、
求められた前記尤度関数を前記各画像間の類似性の尺度として用いることを特徴とする複数画像間の非剛体レジストレーション方法。
前記複数の画像間で前記標的組織の位置合せを行なうことにより、該標的組織に対応する空間領域の変形ベクトル場を推定することを特徴とする請求項１または２記載の複数画像間の非剛体レジストレーション方法。
前記変換の前記パラメータは、Ｂ−スプラインの制御点の位置を示すものであることを特徴とする請求項１〜３までのいずれか１項記載の複数画像間の非剛体レジストレーション方法。
前記複数の画像は、Ｘ線ＣＴ撮影により得られたものであることを特徴とする請求項１〜４までのいずれか１項複数画像間の非剛体レジストレーション方法。