実施例1の画像処理装置は、乳房を被検体(検査の対象物)として、光音響断層撮像装置(PAT)に搭載された赤外線カメラの画像およびPAT画像と、MRI画像とを比較して、PAT画像とMRI画像の変形位置合わせを行う。すなわち、MRI撮像時の被検体の位置や形状を「第一の状態」、PAT撮像時の被検体の位置や形状を「第二の状態」と呼ぶ場合に、第一の状態の被検体を表すMRI画像に変形が施され、第二の状態の被検体に対して位置合わせが行われる。
具体的な処理としては、最初に、PAT撮像前の非保持状態(以下、第二の状態の初期段階)における被検体をPAT搭載の赤外線カメラで撮像した二次元画像が取得され、MRI画像と該二次元画像との位置合わせが行われる。すなわち、第一の状態と第二の状態(正確には、第二の状態の初期段階)の被検体の位置合わせパラメータとして、両者の間の剛体変換が推定される。そして、この剛体変換を初期値として、第一の状態とPAT撮像時(すなわち、第二の状態)の被検体の位置合わせパラメータとして、圧迫変形の変形パラメータが推定される。実施例1では、これら二段階の処理により、PAT画像とMRI画像の間の位置合わせパラメータが導出される。
[装置の構成]
図1のブロック図により実施例1の画像処理装置10を含むモダリティシステムの構成例を示す。画像処理装置10は、医用画像データベース(DB)11と光音響断層撮像装置(PAT)12に接続する。医用画像DB11は、被検体である乳房をMRIで撮像した三次元画像データを保持する。PAT12は、PAT画像を撮像する装置であり、被検体のPAT画像と赤外線カメラ画像を保持する。
●医用画像DB
図2により医用画像DB11が保持する被検体のMRI画像を説明する。図2(a)に示す被検体のMRI画像300は、人体の頭尾方向に垂直な断面(アキシャル断面)でスライスした二次元画像(乳頭304を含む断面)の集合(三次元画像データ)である。MRI画像300を構成する画素は、MRI画像座標系CMRIにおける位置が定義されている。また、MRI画像300は、被検体の体外領域302と、被検体の体内領域303の撮像結果を含む。
図2(b)に示すMRI画像301は、人体の左右方向に垂直な断面(サジタル断面)でスライスした二次元画像の集合(三次元画像データ)である。MRI画像301は、MRI画像300と同様に、被検体の対外領域302と体内領域303に撮像結果を含む。なお、本実施例においては、患者の右手側から左手側をx軸の正方向、胸側から背側をy軸の正方向、足側から頭側をz軸の正方向とするMRI画像座標系CMRIを定義する。
●PAT
図3によりPAT12による被検者の撮像を説明する。被検者500は、PAT12の上面のベッドに伏臥位の体位をとる。そして、被検体である片方の乳房501をPAT12上面の開口部502に挿入する。このとき、照射光が乳房の内部まで届くように、乳房501は透明な二枚の保持板(足側の固定保持板503と頭側の可動保持板504)により圧迫された状態で保持され、乳房501の厚みが薄くなった状態で撮像される。保持は、可動保持板504を足方向(固定保持板503方向)に移動することで行われる。
固定保持板503と可動保持板504は何れも平板であり、乳房501との接触面(以下、保持面)は平面とする。また、保持時の固定保持板503と可動保持板504の間の距離(以下、保持厚)がPAT12によって計測され、保持厚は、PAT画像の付加的情報として当該画像のヘッダ部に保存される。
照射光である近赤外光パルスが、保持板の平面に直交する方向から図示しない光源によって照射される。そして、被検体内で発生した光音響信号が、保持板の平面に直交するように配置された図示しない超音波探触子によって受信される。
PAT12にはPAT装置座標系CDEVが定義されており、x-y面は、固定保持板503と可動保持板504の平面に平行な面であり、z軸は保持された乳房501の厚みの方向である。例えば、MRI画像座標系CMRIと同様に、被検者500の右手側から左手側をx軸の正方向、胸側(下)から背側(上)をy軸の正方向、足側から頭側をz軸の正方向と定義する。PAT装置座標系CDEVの原点は、例えば、固定保持板503上の右手側の下端位置に設定される。以降、PAT12において、上記座標系を基準にして、他の座標系との関係を扱うものとする。
図4によりPAT12が撮像したPAT画像の例を示す。PAT画像600は、図2(a)と同様にアキシャル断面の二次元画像の集合(三次元画像データ)である。本実施例では、MRI画像座標系CMRIと同様に、被検者500の右手側から左手側をx軸の正方向、胸側から背側をy軸の正方向、足側から頭側をz軸の正方向とするPAT画像座標系CPATを定義する。
PAT画像座標系CPATからPAT装置座標系CDEVへの変換を行う座標変換行列を「TPtoD」と定義する。なお、TPtoDを含め、以降に登場する座標変換行列は、すべて座標系の並進と回転を表す4×4行列とする。PAT画像座標系CPATは、PAT装置座標系CDEVに対して平行座標系であり、被検体501の撮像範囲に応じてCPATの原点位置が変化する。つまり、座標変換行列TPtoDには回転成分がなく、撮像範囲に基づき一意に算出することができる。座標変換行列TPtoDは、PAT画像の付加的情報として当該画像のヘッダ部分に保存される。
図3に示すように、PAT12には、被検体501の外観と体表近傍の血管の様子を撮像するための三台の赤外線カメラ(前面赤外線カメラ505、後面赤外線カメラ506、側面赤外線カメラ507)が搭載されている。前面赤外線カメラ505は、可動保持板504を通して頭側から被検体501の外観を撮像可能な位置に設置される。後方赤外線カメラ506は、固定保持板503を通して足側から被検体501の外観を撮像可能な位置に設置される。側面赤外線カメラ507は、側面から被検体501の外観を撮像可能な位置に設置される。
PAT12は、赤外線カメラ505から507によって撮像される、被検体501を保持していない状態(以下、非保持状態)と、被検体501を保持した状態(以下、保持状態)における被検体501の画像を保存する機能を有す。以下、保持状態において、前面赤外線カメラ505が撮像する画像をICAM1、後面赤外線カメラ506が撮像する画像をICAM2、側面赤外線カメラ507が撮像する画像をICAM3とする。また、非保持状態において、前面赤外線カメラ505が撮像する画像をI'CAM1、後面赤外線カメラ506が撮像する画像をI'CAM2、側面赤外線カメラ507が撮像する画像をI'CAM3とする。
前面赤外線カメラ505の座標系(前面カメラ座標系)CCAM1のz軸(視軸の負の方向を示す)は、PAT装置座標系CDEVのz軸と略同一方向を向く。同様に、後面赤外線カメラ506の座標系(後面カメラ座標系)CCAM2のz軸は、PAT装置座標系CDEVのz軸と略反対方向を向く。また、側面赤外線カメラ507の座標系(側面カメラ座標系)CCAM3のz軸は、PAT装置座標系CDEVの-x軸方向を向く。
カメラ座標系CCAM1、CCAM2、CCAM3からPAT装置座標系CDEVへの座標変換行列をそれぞれ、TC1toD、TC2toD、TC3toDと定義する。赤外線カメラ505から507は、PAT装置座標系CDEVにおいて校正済み(言い替えれば、PAT12との位置関係が既知)である。また、上記の座標変換行列や赤外線カメラ505から507の内部パラメータは、既知の情報として、画像処理装置10に保持されている。
図5により非保持状態における前面赤外線カメラ505の撮像画像ICAM1の一例を示す。赤外線カメラは、近赤外線の強度情報が可視化された画像を撮像する装置であり、近赤外線の以下の性質により、被検体の表層部である皮膚直下の静脈血管(表在血管)を可視化する性質がある。
・近赤外線によって皮膚がある程度透ける
・近赤外線がヘモグロビンを含む血管部分に吸収され、血管が周囲よりも暗く写る
赤外線カメラ画像700は、皮膚下の表在血管の形状を明瞭に描出するため、形態画像として扱うことができる。図5において、赤外線カメラ画像700上には、保持状態の乳房輪郭形状701が写り、さらに乳頭画像702、そして表在血管画像703が写っている。
前面赤外線カメラ505の撮像画像ICAM1の二次元座標系CIMG1上の座標は、三次元のカメラ座標系CCAM1においては、原点である焦点位置と、三次元空間中のカメラの投影面上の点を通る直線、つまり視線と一対一の関係を有す。前面赤外線カメラ505の撮像画像の座標系CIMG1とカメラ座標系CCAM1の間の変換については、一般的な撮像画像と三次元空間の間の座標変換方法を用いるため説明は省略する。また、後面赤外線カメラ506と側面赤外線カメラ507の撮像画像は、視点位置が異なる点を除き、前面赤外線カメラ505の撮像画像と同様であるから、詳細説明を省略する。
●画像処理装置
図6のフローチャートにより実施例1の画像処理装置10の各部の動作と処理を説明する。
医用画像取得部101は、医用画像DB11が保持する被検体のMRI画像を取得し、MRI画像を三次元形状取得部102、剛体変換部106、変形画像生成部111に出力する(S201)。
三次元形状取得部102は、入力されるMRI画像に画像処理を施して、被検体の表面に当る各画素の位置(表面位置)を検出し、被検体の表面形状を示す情報を取得する(S202)。さらに、検出した表面位置から得られる形状の三次元曲率に基づき、MRI画像中の特徴点の位置を取得する(S203)。被検体が乳房の場合、MRI画像中の特徴点は乳頭であり、以下では、三次元形状取得部102が乳頭位置を示す情報を取得する情報取得部として説明を続ける。三次元形状取得部102は、取得した表面形状と乳頭位置を剛体変換部106、変形推定部109に出力する。本実施例において、MRI画像から取得される被検体の表面形状は、非保持状態における被検体の形状モデルになる。
図7により表面形状の取得処理を説明する。図7(a)(b)は、図2(a)(b)に示すMRI画像300、301から被検体の体外領域302と体内領域303の境界304(表面位置)を検出した表面検出画像400、401を示す。表面検出画像400、401は、例えば、被検体の表面とそれ以外を区別可能な二値画像などでよい。
本実施例では、被検体の表面形状として、NS個の点群PSk(1≦k≦NS)を取得し、それら点群の位置をMRI画像座標系CMRIにおける三次元の位置座標ベクトルVSk_MRIとして記録する。
PAT画像取得部103は、PAT12が被検体を撮像したPAT画像を取得し、PAT画像を変形画像評価部111、画像表示部112に出力する(S204)。また、PAT画像のヘッダ部に含まれる付加的情報、例えば、座標変換行列TPtoD、TC1toD、TC2toD、TC3toDなども変形画像評価部111に出力される。なお、PAT画像取得部103が取得するPAT画像は、所定の波長に対して被検体内の光エネルギ堆積量分布を画像化した三次元画像とする。
カメラ画像取得部104は、PAT12の赤外線カメラ505から507が撮像した、非保持状態と保持状態における被検体の赤外線カメラ画像を取得し、赤外線カメラ画像を二次元形状取得部105、仮想投影像評価部108に出力する(S205)。ここで取得される赤外線カメラ画像はICAMi、I'CAMi(i=1, 2, 3)である。
PAT画像取得部103とカメラ画像取得部104は、PAT12の撮像に同期してPAT12から直接画像を取得してもよいし、過去に撮像され記録された画像を図示しない医用画像記録装置から取得してもよい。
二次元形状取得部105は、入力される各赤外線カメラ画像に画像処理を施して、乳房輪郭形状(図6の701)と乳頭画像(図6の702)の位置を取得し、乳房輪郭形状と乳頭位置を剛体変換部106、変形画像評価部111に出力する(S206)。例えば、乳房輪郭形状の検出には一般的なエッジ検出手法を用いることができる。また、乳頭位置は、乳房領域の境界を表す曲線の曲率に基づき検出することができる。なお、乳房輪郭形状と乳頭位置の検出は、上記の方法に限られるわけではない。
剛体変換部106、仮想投影像生成部107、仮想投影像評価部108は、MRI画像と、非保持状態における赤外線カメラ画像に写った表在血管の情報を用いて、非保持状態における被検体とMRI画像の位置合わせを行う(S207)。なお、非保持状態とは第二の状態の初期段階である。具体的には、仮定した位置合わせパラメータの候補値のそれぞれに基づき、該赤外線カメラでMRI画像を仮想的に観測した仮想画像を生成し、これを赤外線カメラ画像と比較することで該位置合わせパラメータを推定する。この位置合わせ(非保持状態における位置合わせ)の詳細は後述するが、本処理により、MRI画像座標CMRIから前面カメラ座標系CCAM1への剛体変換を表す変換行列TMtoC1を、位置合わせパラメータとして取得する。
剛体変換部106は、変換行列TMtoC1に基づき、MRI画像座標系CMRIからPAT装置座標系CDEVへの剛体変換を表す変換行列TMtoDを算出する(S208)。つまり、変換行列TMtoC1に、画像処理装置10が保持する前面カメラ座標系CCAM1からPAT装置座標系CDEVへの変換行列TC1toDを掛けることで、変換行列TMtoDが算出される。
変形推定部109、変形画像生成部110、変形画像評価部111は、非保持状態における位置合わせ結果に基づき、保持状態における被検体とMRI画像の位置合わせ(以下、圧迫変形の推定)を行う(S209)。詳細は後述するが、物理変形シミュレーションを用いてMRI画像の圧迫変形を推定する。つまり、変形パラメータを様々に変化させて物理変形シミュレーションを行い、PAT画像との比較によって変形の適切性を表す所定の評価値を得て、該評価値を最小とする変形パラメータを位置合わせパラメータとして推定する。そして、推定した変形パラメータを用いて、変形MRI画像ID_MRIonP(変形三次元画像)を推定結果として生成する。
画像表示部109は、生成された変形三次元画像(変形MRI画像)と、ステップS203で取得されたPAT画像を図示しないモニタに並べて表示する(S210)。図8により変形MRI画像とPAT画像の表示例を示す。図8の例は、同一のアキシャル断面の変形MRI画像1400とPAT画像600を上下に並置した例である。変形MRI画像1400に重畳された破線の矩形1500は、PAT画像600の表示領域に対応する領域をユーザに示す表示情報である。
なお、三次元形状取得部102、二次元形状取得部105、剛体変換部106、仮想投影像生成部107、仮想投影像評価部108、変形推定部109、変形画像生成部110、変形画像評価部111は、位置合わせ部113を構成する。
●非保持状態における位置合わせ
非保持状態における位置合わせにおいて、MRI画像座標系CMRIから前面カメラ座標系CCAM1への剛体変換を推定する。図9のフローチャートにより非保持状態における位置合わせ(S207)の詳細を説明する。
剛体変換部106は、MRI画像を赤外線カメラ座標系(前面カメラ座標系)に平行移動するためのパラメータを算出する(S801)。まず、非保持状態における各赤外線カメラ画像から得た二次元の乳頭位置から、三角測量の原理に基づき、前面カメラ座標系CCAM1における三次元の乳頭位置を算出する。そして、MRI画像中の乳頭位置が、赤外線カメラ画像から求めた非保持状態における三次元の乳頭位置に一致するように、MRI画像座標系CMRIから前面カメラ座標系CCAM1への平行移動を表す変換行列T1MtoC1を算出する。
次に、剛体変換部106は、MRI画像中の被検体を回転移動するための各成分(三軸回りの回転角度)が取り得る値を組み合わせた複数(Nθ組)の回転パラメータの候補値(仮説)θi={θx, θy, θz}(1≦i≦nθ)を設定する(S802)。これは、言い換えると、本処理における回転パラメータの候補値θiと、ステップS801で算出した平行移動のパラメータとを組み合わせた、剛体変換のパラメータの候補値を設定したことになる。また、PAT画像座標系と前面カメラ座標系との関係が既知であることを考えると、これは、MRI画像(第一の状態の被検体)からPAT画像(第二の状態の被検体)への剛体変換の候補値を設定することとも等価である。例えば、x軸周り回転角度をθx、z軸回りの回転角度θzとすると、-10度から+10度の範囲で五度刻みに下記の五つの角度を設定する。
θx = {-10, -5, 0, +5, +10}
θz = {-10, -5, 0, +5, +10}
また、y軸回りの回転角度をθyとすると、-180度から+180度の範囲で五度刻みに下記の72の角度を設定する。
θy = {-180, -175, …, -5, 0, +5, …, +175, +180}
このとき、回転パラメータθiの取り得る値の数(候補値(仮説)の総数)Nθは1800になる(つまり、1≦i≦1800)。
次に、剛体変換部106は初期化を行う(S803)。つまり、ループ変数iに1を、後述する類似度Siの最大値SMAXに0、後述する角度θMAXにθ1を設定する。
次に、剛体変換部106は、平行移動後のMRI画像を、乳頭位置を基準に回転パラメータθiだけ回転移動したMRI画像IMRIonC1iと位置座標ベクトルvSk_MRIonC1iを仮想投影像生成部107に出力する(S804)。つまり、前面カメラ座標系CCAM1において回転パラメータθiだけ回転移動させる変換行列T2iを算出する。続いて、ステップS801で導出した平行移動の変換行列T1MtoC1に変換行列T2iを掛けた剛体変換を表す変換行列TMtoC1iを導出する。続いて、MRI画像と被検体の表面形状の位置座標ベクトルvSk_MRIを、それぞれ変換行列TMtoC1iを用いて座標変換したMRI画像IMRIonC1iと位置座標ベクトルvSk_MRIonC1iを生成する。
仮想投影像生成部107は、剛体変換されたMRI画像を前面赤外線カメラ505の視点から観測した場合に視野に入るであろう被検体の表面形状を部分表面領域として求める(S805)。言い換えれば、表面位置を示す情報である表面形状と、位置合わせパラメータの候補値である剛体変換パラメータ(剛体変換を表す変換行列TMtoC1i)に基づき、剛体変換されたMRI画像を前面赤外線カメラ505の視点から見た仮想画像が生成される。
図10により仮想投影像生成部107が部分表面領域を求める処理(S805)を説明する。図10は前面赤外線カメラ505の視点Pから観測されるMRI画像上の被検体の様子を示し、MRI乳房領域900は、乳頭位置901が赤外線カメラにおける乳頭位置に一致するように剛体変換されている。ステップS805において、仮想投影像生成部107は、視点Pの位置・姿勢を基準に透視投影を行い、視点Pから観測される部分表面領域903を求める。つまり、まず視点Pからの投影線904を延長し、投影線904がMRI乳房領域900の体表形状を表す位置座標ベクトルvSk_MRIonC1iと交点をもつ観測範囲902を特定する。そして、観測範囲902において、位置座標ベクトルvSk_MRIonC1iの中から、投影線904が視点Pを起点にして最初に交わる体表点905を求める。そして、観測範囲902に含まれるすべての投影線904について体表点905を求めると部分表面領域903が決定する。
次に、仮想投影像生成部107は、剛体変換されたMRI画像における部分表面領域903の近傍情報を利用したMIP像を生成し、生成したMIP像を仮想投影像評価部108に出力する(S806)。なお、MIPは「maximum intensity projection」の略であり、以下では、ステップS806で生成されるMIP像を「体表近傍MIP像IMIPonC1i」と呼ぶ。
ステップS806において、仮想投影像生成部107は、投影線904が最初に交わる体表点905を起点にして、視点Pから遠ざかる方向に所定の距離(例えば5mm)をもつ体表近傍区間906を設定する。そして、観測範囲902に含まれるすべての投影線904について体表近傍区間906を設定することで体表近傍領域907を規定する。続いて、視点Pの位置・姿勢を基準に透視投影を行い、MRI画像IMRIonC1iの中で体表近傍領域907に含まれる領域に限定したMIP像である体表近傍MIP像IMIPonC1iを生成する。これにより、MRI画像中の血管領域908の中で前面赤外線カメラ505側の体表近傍に存在する表在血管909の情報のみを可視化したMIP像が生成される。なお、体表近傍領域907の生成時に、被検体の皮膚が当該領域に含まれないようにしてもよい。具体的には、部分表面領域903に皮膚厚相当の所定の厚みを与えた領域を皮膚領域として導出し、前記の処理で得た体表近傍領域907からこれを除外してもよい。被検体の皮膚はMRI画像中において高い輝度値を有しているので、この除外処理により、生成するMIP像中において表在血管をより明瞭に描出できるようになる。
このように、体表近傍、かつ、前面赤外線カメラ505の視点Pから観測されるだろう領域の情報のみを可視化することで、乳房内部の血管領域や、前面赤外線カメラ505からみて反対側の体表近傍に存在する血管領域の可視化を防ぐことができる。つまり、実際の赤外線カメラ画像により近いMIP像である体表近傍MIP像IMIPonC1iを生成することができる。
なお、体表近傍MIP像IMIPonC1iの生成方法は、体表近傍の領域のみを、あるいは、体表近傍の領域を強調して可視化することができれば、上記の方法に限られるものではない。他の方法として、投影線904上において、体表点905からの距離が、視点Pから遠ざかる方向に大きくなるほど、輝度値に対する重みを小さくしてMIP像を生成する方法が挙げられる。この方法によれば、乳房内部の領域で、体表に近い領域ほど輝度値が強調されたMIP像が生成されるため、表在血管が強調されたMIP像が生成される。もちろん、本処理によるMIP像の生成を行う場合においても、皮膚領域を描画対象から除外することで、表在血管をより明瞭に描出できる。
図11によりMRI画像内の被検体の体表近傍情報を利用したMIP像1003を示す。図11において、体表近傍MIP像IMIPonC1iは、前面カメラ画像座標系CIMG1上の二次元画像として表される。これは、視点Pに基づく透視投影によって二次元のMIP像を生成する際の座標変換(S806)が、三次元の前面カメラ座標系CCAM1から二次元の前面カメラ画像座標系CIMG1への座標変換(カメラの撮像プロセス)と幾何的に等しいためである。
また、MIP像の生成の際、体表の外側領域は処理対象に含まれないため、乳房輪郭形状1000の外側には有意の輝度値は存在しない。さらに、MRI画像における乳頭位置901が、前面赤外線カメラ505の撮影画像ICAM1における乳頭位置702に一致するように剛体変換されているため、MIP像の乳頭1001は、撮影画像ICAM1(図7)の乳頭703に一致する。そして、MIP像においては、表在血管1002がとくに高輝度な領域として可視化される。
次に、仮想投影像評価部108は、体表近傍MIP像IMIPonC1iと、非保持状態の赤外線カメラ画像I'CAM1の両方において可視化された表在血管の輝度情報に基づき、両画像の類似度Siを算出する(S807)。
ステップS807において、仮想投影像評価部108は、体表近傍MIP像IMIPonC1i、非保持状態の赤外線カメラ画像I'CAM1ともに、類似度の算出領域から乳房の外側領域を除外し、類似度の算出領域を乳房の内部領域に限定する。
体表近傍MIP像IMIPonC1iにおける表在血管1002は、周囲の乳房領域に比べて高輝度に可視化されている。他方、赤外線カメラ画像I'CAM1における表在血管703は、周囲の乳房領域に比べて低輝度に可視化されている。そこで、仮想投影像評価部108は、それら画像の間の輝度情報を直接比較できるようにするため、体表近傍MIP像IMIPonC1iの輝度値を反転する。そして、輝度値を反転した体表近傍MIP像と赤外線カメラ画像I'CAM1の間で類似度Si(0≦Si≦1)を算出する。体表近傍MIP像IMIPonC1iの表在血管1002(図11)と、赤外線カメラ画像I'CAM1の表在血管703(図5)が類似すれば、類似度Siの値が高くなる(1に近付く)ものとする。
なお、本実施例においては、類似度Siの評価尺度として画像間の相互情報量を適用するが、評価尺度はこれに限られず、相互相関係数やSSD (sum of squared difference)など、既知の手法を用いてもよい。また、輝度値に直接的に基づく評価尺度である必要はなく、例えば両画像からエッジなどの画像特徴を検出して、それらの類似度や一致度を算出するような尺度でもよい。
次に、仮想投影像評価部108は、類似度Siと類似度の最大値SMAXと比較する(S808)。そして、類似度Siが最大値SMAXを超える(Si>SMAX)場合は、SMAXを更新し(SMAX=Si)、SMAXに対応する角度θMAXを更新する(θMAX=θi)(S809)。また、類似度Siが最大値SMAX以下(Si≦SMAX)の場合は更新を行わない。
次に、仮想投影像評価部108は、ループ変数iをインクリメントし(S810)、ループ変数iと仮説の総数Nθを比較する(S811)。ループ変数iが仮説の総数Nθ以下(i≦Nθ)の場合は処理をステップS804に戻し、ループ変数iが仮説の総数Nθを超える(i>Nθ)場合は処理をステップS812に進める。つまり、仮説の総数Nθ分、ステップS804からS811の処理が繰り返される。
仮説の総数Nθ分の処理が終了すると、剛体変換部106は、角度θMAXにおける変換行列TMtoC1MAXを、MRI画像座標CMRIから前面カメラ座標系CCAM1への剛体変換を表す最終的な変換行列TMtoC1とする(S812)。言い換えれば、複数の回転パラメータから、類似度の最大値SMAXに対応する回転パラメータθMAXが選択される。
以上で、剛体変換部106、仮想投影像生成部107、仮想投影像評価部108が行う、非保持状態における位置合わせ(S207)の処理が終了する。この処理により、様々な剛体変換パラメータの仮説(つまり回転パラメータの仮説θi)に基づき生成した体表近傍MIP像IMIPonC1iが、赤外線カメラ画像I'CAM1に最も類似する角度θMAXに基づく変換行列TMtoC1が得られる。
上記では、前面カメラ画像座標系CIMG1に変換した体表近傍MIP像IMIPonC1iを生成し、前面赤外線カメラ505の画像I'CAM1との類似度Siを評価する例を説明した。しかし、類似度Siの評価対象は、前面赤外線カメラ505に限られるわけではない。
例えば、後面カメラ画像座標系CIMG2や側面カメラ画像座標系CIMG3に変換した体表近傍MIP像IMIPonC1iを生成する。そして、体表近傍MIP像IMIPonC1iと後面赤外線カメラ506の画像I'CAM2または側面赤外線カメラ507の画像I'CAM3との類似度を評価してもよい。その場合、ステップS804で得た変換行列の候補値TMtoC1iに対して、ステップS805、S806で透視投影を行う視点Pを、後面または側面赤外線カメラ506、507のカメラ視点にそれぞれ置き換えればよい。
ここで、MRI画像IMRIonC1iと位置座標ベクトルvSk_MRIonC1iはステップS804で前面カメラ画像座標系CIMG1に変換されている。従って、この透視投影に用いるカメラ視点の位置・姿勢としては、前面カメラ座標系CCAM1上で表現される位置姿勢を用いればよい。
なお、上述したように、各赤外線カメラ(前面、側面、後面)の位置関係は、PAT装置座標系CDEVを基準に対応付けられている。そのため、前面カメラ座標系CCAM1から後面カメラ座標系CCAM2または側面カメラ座標系CCAM3への変換行列を導出することができる。つまり、前面カメラ座標系CCAM1における後面または側面赤外線カメラのカメラ視点の位置姿勢を導出することができる。
なお、これらの前面、後面、側面の赤外線カメラそれぞれに基づく類似度を総合した総合類似度を算出して、類似度を評価してもよい。総合類似度としては、これら三種類の類似度の、例えば、重み付け平均値、最大値、最小値、中央値が挙げられる。また、上述では、前面カメラ座標系CCAM1を基準にして、回転パラメータθiの設定、及びMRI画像を各カメラ座標系に透視投影するための座標変換を行った。しかし、PAT装置座標系CDEVを基準にして、CMRIとCDEVとの間の剛体変換TMtoDを位置合わせパラメータとしてその候補値を設定するようにしてもよい。この場合、CMRIからCDEVへの剛体変換に加えて各赤外線カメラへのビューイング変換を行った後に、ステップS806と同様の処理で体表近傍MIP像を生成する。そして、各赤外線画像との類似度評価に基づいて位置合わせパラメータを推定すればよい。
●圧迫変形の推定
図12のフローチャートにより圧迫変形の推定(S209)の詳細を説明する。
変形推定部109は、ステップS202で取得された被検体の表面形状、および、ステップS208で取得された座標変換行列TMtoDを用いて、被検体の形状を表す三次元のメッシュ(以下、メッシュM)を生成する(S901)。つまり、座標変換行列TMtoDによる座標変換を、被検体の表面形状VSk_MRIに施して、PAT装置座標系CDEVにおける被検体表面点群の位置座標ベクトルVSi_MRIonD(1≦i≦Ns)を算出する。そして、位置座標ベクトルVSi_MRIonDが表す表面形状に基づき、被検体の内部領域を判別し、内部領域にメッシュMを配置する。
図13の模式図によりメッシュMの生成方法を示す。図13(a)は被検体の処理対象領域1200のサジタル断面を示し、サジタル断面における被検体の表面位置1201と、対応する内部領域1202を示す。図13(b)に示すように、メッシュMは、被検体の内部領域1202に六面体や四面体などの立体構造をもつ要素1203を配置して生成される。メッシュMは、これらの要素の頂点(ノード)1204の位置と連結情報によって記述される。
以下では、ステップS901において配置されるメッシュMのノード数をNm、各ノードの位置をsL(1≦L≦Nm)と表記する。各ノードの変位によって要素内の変位場を表現することができ、これに基づき、被検体内の任意の点の変位を求めることができる。
次に、変形推定部109は、変形パラメータの各成分(被検体のヤング率やポアソン比など)が取り得る値を組み合わせた複数(Np組)の変形パラメータの仮説pk(1≦k≦Np)を生成する(S902)。例えば、各成分が取り得る値の範囲を適切な間隔で分割し、それらすべてを組み合わせることで変形パラメータpkを生成する。例えば、ヤング率の比pyとポアソン比ppを変形パラメータpkの成分とし、ヤング率の比pyとポアソン比ppが取り得る値を次とする。
py = {1, 2, 3, 4, 5}
pp = {0.0, 0.2, 0.4, 0.45, 0.499}
そして、ヤング率の比pyとポアソン比ppを組み合わせた変形パラメータを生成する。上記の例ではNp=25になる。なお、ヤング率の比pyは、乳房の硬さの非等方性に対応するためのパラメータであり、人体のコロナル面(x-z平面)におけるヤング率に対する、人体の前後方向(y軸方向)のヤング率の割合を表す。
次に、変形推定部109は初期化を行う(S903)。つまり、ループ変数kに1を、後述する評価値の最大値EMAXに0、後述する変形パラメータpMAXにp1を設定する。
次に、変形推定部109は、変形パラメータpkを用いて、有限要素法に基づく物理変形シミュレーションをメッシュMに施し、変形後のメッシュである変形メッシュDMkを生成する(S904)。このときの変形関数Fk(x, y, z)を、メッシュMから変形メッシュDMkへの各ノードの変位を与える変位ベクトルdkL(1≦L≦Nm)と定義する。
図14によりステップS904における物理変形シミュレーションである、保持板による圧迫変形シミュレーションを説明する。保持板による圧迫変形において、二枚の保持板を被検体の中心方向に移動すると、移動後の保持板に接する被検体の表面領域は保持板に張り付いた状態になる。
図14(a)に示すように、二枚の保持板をそれぞれ、距離Δd1、Δd2分、移動すると仮定する。メッシュMのノードのうち、体表面を表す表面ノードから、保持面PUd1、PLd2に接する外側表面ノード1300と1301を抽出し、外側表面ノード1300と1301それぞれが保持面に張り付くための変位量を求める。そして、変位量を変形シミュレーションの境界条件Cとして有限要素法の計算を実行し、二枚の保持板がそれぞれ、距離Δd1、Δd2分、移動した場合の変形メッシュを生成する。
本実施例においては、図14(b)に示す二枚の保持板の最終的な保持位置PU、PLまで移動する間を複数回(N回)の変形シミュレーションに分割して、変形過程で生じる境界条件の変化に対応する。つまり、図14(b)はN回の変形シミュレーションを繰り返した結果の変形メッシュDMkを表す。物理変形シミュレーションにより、図14(b)に示すように、被検体の変形メッシュが保持位置PU、PLの間でz軸方向に圧縮され、y軸方向に伸長されることが分かる。
変形画像生成部110は、変形パラメータpkに対応する変形をMRI画像に施して、変形MRI画像を生成し、変形MRI画像を変形画像評価部111に出力する(S905)。つまり、MRI画像を、座標変換行列TMtoDを用いてPAT装置座標系CDEVに座標変換し、ステップS904で算出した変形関数Fkを用いて変形処理する。そして、座標変換行列TPtoDの逆行列を用いた座標変換により、PAT画像座標系CPATにおける変形MRI画像ID_MRIonPkを生成する。
図15の模式図により変形MRI画像ID_MRIonPを示す。図15は変形MRI画像1400、変形後の乳房領域1401、変形前の乳房形状1402を示す。また、図15(a)は変形MRI画像ID_MRIonPをアキシャル断面でスライスした二次元画像、図15(b)は変形MRI画像ID_MRIonPをサジタル断面でスライスした二次元画像である。変形後の乳房領域1401と変形前の乳房形状1402を比較すると、PAT画像座標系CPATのz軸方向への圧迫により、x-y平面上において乳房領域が伸長し、z軸方向に圧縮されていることがわかる。
変形画像評価部111は、ステップS203で取得されたPAT画像と、ステップS205で取得された保持状態の赤外線カメラ画像を用いて算出した変形MRI画像ID_MRIonPkの適切性の評価値Ekを変形推定部109に出力する(S906)。つまり、変形画像評価部111は、変形MRI画像とPAT画像の間の類似度SMRIk(0≦SMRI、k≦1)、および、変形MRI画像と保持状態の赤外線カメラ画像における乳房形状の間の残差Rkに基づき、評価値Ekを算出する。
なお、評価値Ekの値が高いほど変形が適切であるものとする。また、類似度SMRIkの評価尺度として、ステップS807と同様に画像間の相互情報量を適用する。なお、評価尺度はこれに限らず、相互相関係数やSSD、血管分岐部等の特徴点の位置の一致度など、既知の何れの手法を用いてもよい。
図8は変形MRI画像(ID_MRIonPk)1400とPAT画像600の同一のアキシャル断面を示し、破線の矩形で示す対応領域1500は、PAT画像600に対応する変形MRI画像1400上の領域である。類似度SMRIkは、PAT画像600と対応領域1500の間で計算される。PAT画像600と対応領域1500の間で、可視化された、変形MRI画像1400の血管領域1501とPAT画像600の血管領域1502が類似すれば、類似度SMRIkの値が高くなる。
また、残差Rkは、赤外線カメラ画像に写った被検体の輪郭(シルエット)の形状と、赤外線カメラ画像に投影した変形メッシュDMkの外郭の形状の間の差として算出される。例えば、変形メッシュDMkを保持状態の赤外線カメラ画像ICAM1へ投影する場合、変形メッシュDMkは、座標変換行列TC1toDの逆行列を用いて前面カメラ座標系CCAM1に座標変換され、前面カメラ画像座標系CIMG1に投影変換される。また、同様の方法により、変形メッシュDMkを後面カメラ画像ICAM2または側面カメラ画像ICAM3へ投影することができる。
残差Rkは、例えば、変形メッシュDMkを三つの赤外線カメラ画像それぞれに投影したメッシュの外郭と、三つの赤外線カメラ画像の間の各残差を総合した総合残差(例えば、三つの残差の重み付け平均値)として算出される。ただし、総合残差は、重み付け平均値に限らず、三つの残差の最大値、最小値、中央値などを用いてもよい。
また、残差Rkは、例えば、赤外線カメラ画像に写った被検体の乳頭位置と、赤外線カメラ画像に投影した変形メッシュDMkの乳頭位置の間の残差として算出してもよい。勿論、乳房形状の残差と、乳頭位置の残差の両方を統合した値(例えば、重み付け和)を残差Rkとしてもよい。なお、乳房形状や乳頭位置に基づく残差は、赤外線カメラ画像から取得する例を説明したが、乳房の外観を撮像した一般的なカメラ画像から取得してもよい。
評価値Ekは、例えば、類似度SMRIkと残差Rkに基づく重み付け和として、下式で表される。
Ek = aSMRIk + b{1/(1+Rk)} …(1)
ここで、a、bは重み係数(a+b=1)。
式(1)の第二項において(1+Rk)の逆数を用いるのは次の理由からである。
・残差Rkは、評価値Ekとは逆に、変形が適切であるほど値が小さくなる指標である
・値の取り得る範囲を、類似度SMRIkと同様に、0から1の間にする
変形推定部109は、入力される評価値Ekと評価値の最大値EMAXを比較する(S907)。そして、評価値Ekが評価値の最大値EMAXを超える(Ek>EMAX)場合は、EMAXを更新し(EMAX=Ek)、EMAXに対応する変形パラメータpMAXを更新する(pMAX=pk)(S808)。また、評価値Ekが評価値の最大値EMAX以下(Ek≦EMAX)の場合は更新を行わない。
次に、変形推定部109は、ループ変数kをインクリメントし(S909)、ループ変数kと仮説の総数Npを比較する(S910)。ループ変数kが仮説の総数Np以下(k≦Np)の場合は処理をステップS904に戻し、ループ変数kが仮説の総数Npを超える(k>Np)の場合は処理をステップS911に進める。つまり、仮説の総数Np分、ステップS904からS910の処理が繰り返される。
仮説の総数Np分の処理が終了すると、変形推定部109は、変形パラメータpMAXを変形画像生成部110に出力する(S911)。言い換えれば、複数の変形パラメータから、評価値の最大値EMAXに対応する変形パラメータpMAXが選択される。変形画像生成部110は、変形パラメータpMAXに対応する変形MRI画像ID_MRIonPMAXを圧迫変形の推定結果(ID_MRIonP)とし、変形MRI画像ID_MRIonPを画像表示部112に出力する(S912)。
以上で、変形推定部109、変形画像生成部110、変形画像評価部111による圧迫変形の推定(S209)が終了する。この処理によれば、様々な変形パラメータの仮説pkの下で変形シミュレーションを実行する。そして、それらシミュレーション結果の中で変形の適切性の評価値Ekが最大になる変形パラメータpMAXによって変形MRI画像ID_MRIonPが生成される。
このように、乳房を被検体として、PAT12に搭載された赤外線カメラの画像およびPAT画像と、MRI画像との比較により、PAT画像とMRI画像の位置合わせ精度の向上を図ることができる。つまり、MRI画像から生成した表在血管を強調した二次元のMIP像と、非保持状態の赤外線カメラ画像を比較して、赤外線カメラに対するMRI画像上の被検体の位置・姿勢を推定し、MRI画像と赤外線カメラ画像の高精度な剛体位置合わせを実現する。さらに、MRI画像を赤外線カメラの座標系からPAT12の座標系に座標変換することで、剛体位置合わせの結果を、MRI画像とPAT画像の間の変形位置合わせ処理の初期状態として活用する。つまり、MRI画像とPAT画像を比較して位置合わせを行う際は、非保持状態から保持状態への圧迫変形を推定するだけでよい。そして、赤外線カメラ画像に写った乳房形状と、PAT画像に写った内部構造の両方に合致するようにMRI画像を変形することで、高精度な圧迫変形の推定が可能になる。
なお、本実施例では、MRI画像を変形位置合わせする手法として、有限要素法に基づく物理変形シミュレーションを用いたが、必ずしもこの手法に限られるものではない。例えば、一般的な変形手法である、FFD (Free Form Deformation)などの手法を用いてもよい。FFDを用いた変形処理では、まず画像内の被検体を直方体で取り囲むように、格子状の制御点を配置する。そして、この制御点を動かすことにより、直方体内部に存在する画像領域を変形させることができる。ここで、各制御点の変位量のセットを、変形パラメータの候補値(仮説)pk(1≦k≦Np)と定義する。そして、変形パラメータpkの値を様々に変動させ、上述の変形画像評価部111における評価値Ekを最大化するような変形パラメータpkを算出することで、上述の目的に合致する変形位置合わせを実現することができる。
なお、本実施形態では、被検体として人体の乳房を適用したが、これに限られるものではなく、表在血管を有する生体の部位であれば被検体は何であってもよい。また、医用画像データベース110に登録されている画像としてMRI画像を適用したが、生体を撮影した三次元画像であれば、何れのモダリティによる画像であってもよい。
[変形例1]
上記では、赤外線カメラに対するMRI画像内の被検体の位置・姿勢を求めるための回転角度を網羅的に変動させて評価値を算出し、最適な評価値を与える回転角度を取得した(S207)。しかし、MRI画像内の被検体の位置・姿勢の取得に別の方法を用いることができる。つまり、一般的な最適化アルゴリズムを用いて、評価値が最適になる回転角度を推定してもよい。例えば、最適化アルゴリズムの一種である最急降下法を用いる方法を説明する。
三次元ベクトルをx、赤外線カメラに対するMRI画像内の被検体の回転角度を表すパラメータを(θx, θy, θz)とする。回転角度を表すベクトルxを与えたときに、ステップS807で生成される赤外線カメラ画像と体表近傍MIP像との類似度をSMIP(x)とする。最急降下法において、ベクトルxを与えたときに、最小化する関数f(x)を類似度SMIP(x)の逆数1/SMIP(x)とする。ここで、f(x)をSMIP(x)の逆数とする理由は、類似度S(x)が最大になる回転角度を表すパラメータを求めるためである。以上のように設定された各変数を、次式を用いて更新し、収束させることでf(x)を最小化する(つまりSMIP(x)を最大化する)パラメータxを算出する。
x(k+1) = x(k) - αgrad f(x(k)) …(2)
ここで、αは一回に更新する数値の割合を決めるパラメータ(通常、小さな正の定数)、
kは更新回数、
grad f(x(k))はk回目の更新における関数f(x)の勾配ベクトル(関数f(x)の変化率が最も大きい方向を向く)。
勾配ベクトルgrad f(x(k))は以下の方法で求める。k回目の更新におけるベクトルx(k)=(θx(k), θy(k), θz(k))Tとする。ベクトルx(k)の各要素に対して、それぞれ微小な変化量Δx=(Δθx, Δθy, Δθz)を与えたときの、関数f(x(k)+Δx)を算出する。
そのベクトルの方向がパラメータ空間において均等に変化するように変化量Δxを変動させた場合の、各Δxに対応するf(x(k)+Δx)を算出する。算出したf(x(k)+Δx)の集合から、f(x(k))-f(x(k)+Δx)が最大になるΔx(ΔxMAX)を求める。このΔxMAXは、f(x(k))の変化率が最大になるパラメータ空間の方向ベクトルであり、grad f(x(k))に等しい。
また、最適化アルゴリズムとして、ニュートン法など既知の手法の何れを用いてもよい。これにより、より少ない繰り返し回数で、最適な評価値を与える回転角度を推定することができ、処理の高速化を図ることができる。
[変形例2]
上記では、MRI画像の圧迫変形の妥当性を評価する評価値として、PAT画像との類似度と、赤外線カメラに写った被検体の輪郭(シルエット)形状に対する形状誤差を用いた(S209)。しかし、別の方法で評価値を求めてもよい。
例えば、変形メッシュに基づき変形MRI画像を生成し、変形MRI画像を赤外線カメラ視点から投影した変形後の体表近傍MIP像を生成し、体表近傍MIP像と赤外線カメラ画像の類似度を算出して、この類似度を評価値に加えてもよい。
つまり、ステップS905において、メッシュMから変形メッシュDMkへの各ノードの変位を表す変形関数Fkを生成する。そして、MRI画像に対して、MRI画像座標系CMRIから前面カメラ装置座標系CCAM1への変換行列TMtoC1による剛体変換を施し、変形関数Fkによる変形を施して、変形MRI画像ID_MRIonC1kを生成する。
次に、ステップS805、S806と同様の方法により、変形MRI画像ID_MRIonC1kから表面部分領域の近傍情報を利用した体表近傍変形MIP像ID_MIPonC1kを生成する。そして、体表近傍MIP像ID_MIPonC1kと保持状態の赤外線カメラ画像I'CAM1との間の類似度SMIPkを算出する。最後に、評価値EkをステップS906における類似度SMRIk、残差Rk、および、算出した類似度SMIPkに基づく重み付け和を下式により算出する。
Ek = aSMRIk + b{1/(1+Rk)} + cSMIPk …(3)
ここで、a、b、cは重み係数(a+b+c=1)。
これにより、体表近傍の表在血管の情報をさらに利用して変形パラメータを求めることができ、高精度な位置合わせが可能になる。
[変形例3]
上記では、表在血管の可視化に、近赤外線による赤外線カメラ画像を利用する例を説明したが、例えば、体内における内部反射によって得られる偏光成分により表在血管を可視化した画像を利用することもできる。例えばハロゲンライトなどの光を体表に照射し、体表で反射する表面反射成分と、一旦、体内に侵入して体内で吸収と散乱を経た後に再び体表に出てくる内部反射成分を区別する。こうすれば、体内の情報を可視化した画像が得られる。当該画像においては、赤外線カメラ画像と同様、内部反射成分のうち、血管内部のヘモグロビンの吸収箇所が周囲よりも暗く描出される。
このように、第一の撮像装置によって対象物を撮像した三次元画像と、第二の撮像装置によって対象物の表層部を撮像した二次元画像とが取得され、三次元画像から対象物の表面位置を示す情報が取得される。そして、表面位置を示す情報に基づいて、三次元画像を第二の撮像装置の視点から見た投影像が生成され、投影像および二次元画像を用いて、対象物に関して三次元画像と二次元画像との位置合わせが行われる。
以下、本発明にかかる実施例2の画像処理を説明する。なお、実施例2において、実施例1と略同様の構成については、同一符号を付して、その詳細説明を省略する。
実施例1では、非保持状態の被検体の画像とMRI画像の間で剛体変換パラメータを求め(非保持状態における位置合わせ)、その後、保持状態の被検体の画像とMRI画像の間で変形パラメータを推定(圧迫変形の推定)する方法を説明した。この方法は、非保持状態の被検体とMRI画像の被検体の形状が異なる場合に生じる誤差や、非保持状態から保持状態へ移行する際の被検体の微動によって生じる誤差を含む可能性がある。
実施例2では、非保持状態の被検体に関する計測情報を用いずに、保持状態の被検体(第二の状態の被検体)の画像とMRI画像(第一の状態の被検体)の間で、剛体変換パラメータと変形パラメータを同時に推定する方法を説明する。
図16のブロック図により実施例2の画像処理装置10を含むモダリティシステムの構成例を示す。図1に示す実施例1の構成と異なるのは、実施例1の仮想投影像評価部108と変形画像評価部111の代わりに、評価部114を備える点である。また、実施例2においては、位置合わせ部113内における情報の流れが異なるが、その点は、各部の動作と処理の説明において詳述する。
図17のフローチャートにより実施例2の画像処理装置10の各部の動作と処理を説明する。なお、ステップS201からS206、S210における動作と処理は実施例1と同様であり、その詳細説明を省略する。
ステップS220において、赤外線カメラ画像、PAT画像、MRI画像の血管情報に基づき、MRI画像における被検体の位置・姿勢と圧迫変形が推定される。この処理は、剛体変換部106、仮想投影像生成部107、変形推定部109、変形画像生成部110、評価部114によって行われる。
実施例2では、赤外線カメラに対するMRI画像上の被検体の位置・姿勢と、被検体の圧迫変形を表す変形パラメータを位置合わせパラメータとする。そして、仮定した位置合わせパラメータの候補値に基づきMRI画像を赤外線カメラ座標系(前面カメラ座標系)に座標変換した後に圧迫変形を施すことで、当該位置合わせパラメータの候補値に基づく変形MRI画像を生成する。
次に、生成した変形MRI画像に対して、赤外線カメラの視点を基準に透視投影を行ったMIP像を生成する。その際、赤外線カメラの視点からの各投影線上において、乳房の三次元表面の近傍領域の輝度情報のみを可視化、あるいは、強調して可視化する画像処理を施して、乳房の表在血管を強調したMIP像を生成する。
次に、変形MRI画像を、PAT12に対して予め校正された赤外線カメラの位置・姿勢に基づき、赤外線カメラ座標系(前面カメラ座標系)からPAT画像座標系に座標変換する。
次に、血管情報が可視化された赤外線カメラ画像とMIP像の間、および、PAT画像と座標変換後の変形MRI画像の間の画像の類似度を算出し、それら類似度を総合した評価値を算出する。そして、仮定した位置合わせパラメータを変動させて、評価値が最大になる位置合わせパラメータを選択する。すなわち、当該位置合わせパラメータによって、MRI画像とPAT画像の間で圧迫変形を含む位置合わせを行う。
図18、図19のフローチャートにより位置・姿勢と圧迫変形の推定(S220)の詳細を説明する。
変形推定部109は、ステップS202で取得された被検体の表面形状に基づき、被検体の形状を表すメッシュMを生成する(S1101)。この処理は、実施例1のステップS901と略同様であり、詳細説明を省略する。
剛体変換部106は、MRI画像を赤外線カメラ座標系に平行移動する(S1102)。ステップS1102の処理は、ステップS203で取得されたMRI画像における乳頭位置と、ステップS206で取得された保持状態における赤外線カメラ画像上の乳頭位置に基づき行われる。この処理は、実施例1のステップS801の処理と略同様であり、詳細説明を省略する。ただし、ステップS801の処理は非保持状態における赤外線カメラ画像に基づき行われるが、ステップS1102の処理は保持状態における赤外線カメラ画像に基づき行われる点で異なる。
次に、剛体変換部106と変形推定部109は、剛体変換パラメータが取り得る値と、変形パラメータの各成分が取り得る値を組み合わせて複数(Nt組)の位置合わせパラメータ(変換パラメータ)の仮説ti(1≦i≦Nt)を設定する(S1103)。
例えば、剛体変換部106は、剛体変換パラメータが取り得る値として、実施例1のステップS8020と同様に、複数(Nθ組)の回転パラメータθj(1≦j≦Nθ)を設定する。変形推定部109は、変形パラメータが取り得る値として、実施例1のステップS902と同様に、複数(Np組)の変形パラメータpk(1≦k≦Np)を設定する。そして、回転パラメータθjと変形パラメータpkを組み合わせて複数(Nt=Nθ×Np組)の変換パラメータti(1≦i≦Nt)を設定する。なお、PAT画像座標系と前面カメラ座標系との関係は既知であるので、これは、MRI画像(第一の状態の被検体)からPAT画像(第二の状態の被検体)への位置合わせパラメータの候補値を設定することと等価である。また、変換パラメータtiは、剛体変換部106と変形推定部109の間で共有されるものとする。
評価部114は初期化を行う(S1104)。つまり、ループ変数iに1を、後述する評価値の最大値EMAXに0、後述する変換パラメータtMAXにt1を設定する。
剛体変換部106は、平行移動後のMRI画像を、乳頭位置を基準に変換パラメータti(つまりθj)だけ回転移動したMRI画像IMRIonC1iを生成する。そして、MRI画像IMRIonC1iと座標変換行列TMtoC1iを変形推定部109と変形画像生成部110に出力する(S1105)。この処理は、実施例1のステップS804と略同様であり、詳細説明を省略する。
変形推定部109は、変換パラメータti(つまりθjとpk)を用いて、有限要素法に基づく物理変形シミュレーションをメッシュMに施した変形メッシュDMiを生成し、変形関数Fjkを変形画像生成部110に出力する(S1106)。つまり、ステップS1105で導出された座標変換行列TMtoC1iを用いて、回転パラメータθjに対応する剛体変換をメッシュMに施したメッシュMiを生成する。そして、メッシュMiに物理変形シミュレーションを施すことで、変形後のメッシュである変形メッシュDMiを生成する。このときの変形関数Fi(x, y, z)を、メッシュMから変形メッシュDMiへの各ノードの変位を与える変位ベクトルdiL(1≦L≦Nm)と定義する。なお、この処理は、実施例1のステップS904と同様である。
変形画像生成部110は、変換パラメータti(つまりθjとpk)に対応する変換をMRI画像に施した変形MRI画像を生成し、変形MRI画像を仮想投影像生成部107に出力する(S1107)。つまり、MRI画像に座標変換行列TMtoC1iを用いる剛体変換を施して、MRI画像を回転パラメータθjに対応する前面カメラ座標系CCAM1に座標変換する。座標変換後のMRI画像に変形関数Fjkを用いる変形処理を施して、変形MRI画像ID_MRIonC1iを生成する。
仮想投影像生成部107は、前面カメラ座標系CCAM1において、変形MRI画像ID_MRIonC1iを前面赤外線カメラ505の視点から観測した場合に視野に入るだろう被検体の表面形状(部分表面領域)を求める(S1108)。この処理は、実施例1のステップS805における剛体変換されたMRI画像を、変形MRI画像ID_MRIonC1iに置き換えたものであり、詳細説明を省略する。
次に、仮想投影像生成部107は、変形MRI画像ID_MRIonC1iにおける表面部分領域の近傍情報を利用した体表近傍変形MIP像ID_MIPonC1iを生成し、体表近傍変形MIP像を評価部114に出力する(S1109)。この処理は、実施例1のステップS806における剛体変換されたMRI画像を、変形MRI画像ID_MRIonC1iに置き換えたものであり、詳細説明を省略する。
変形画像生成部110は、ステップ1107で生成した変形MRI画像ID_MRIonC1iを前面カメラ座標系からPAT画像座標系に座標変換した変形MRI画像ID_MRIonPiを生成し、変形MRI画像を評価部114に出力する(S1110)。つまり、変形MRI画像ID_MRIonC1iを、座標変換行列TC1toDを用いてPAT装置座標系CDEVに座標変換する。さらに、座標変換行列TPtoDの逆行列を用いる座標変換により、PAT画像座標系CPATにおける変形MRI画像ID_MRIonPiを生成する。
評価部114は、体表近傍変形MIP像ID_MIPonC1iと保持状態の赤外線カメラ画像の類似度SMIPi(0≦SMIPi≦1)、および、変形MRI画像ID_MRIonPiとPAT画像の類似度SMRIi(0≦SMRIi≦1)を算出する。さらに、変形MRI画像ID_MRIonPiと保持状態の赤外線カメラ画像における乳房形状の残差Riを算出する。そして、それら類似度と残差を総合した評価値Eiを算出する(S1111)。
類似度SMIPiの算出方法は、実施例1のステップS807における体表近傍MIP像IMIPonC1jと非保持状態の赤外線カメラ画像を、体表近傍変形MIP像ID_MIPonC1iと保持状態の赤外線カメラ画像に置き換えたものである。また、類似度SMRIiの算出方法は、実施例1のステップS906における変形MRI画像ID_MRIkを、変形MRI画像ID_MRIonPiに置き換えたものである。また、残差Riの算出方法は、赤外線カメラ画像に写った被検体の輪郭(シルエット)の形状と、赤外線カメラ画像に投影した変形メッシュDMjkの外郭の形状の間の差として算出する。残差Riの求め方はステップS906と同様である。
評価値Ekは、例えば、類似度SMIPi、SMRIiと残差Riに基づく重み付け和として、下式で表される。
Ei = aSMIPi + bSMRIi + c{1/(1+Ri)} …(4)
ここで、a、b、cは重み係数(a+b+c=1)。
また、式(4)の第三項において(1+Rk)の逆数を用いるのはステップS906と同様の理由からである。
次に、評価部114は、評価値Eiと評価値の最大値EMAXを比較する(S1112)。そして、評価値Eiが評価値の最大値EMAXを超える(Ei>EMAX)場合は、EMAXを更新し(EMAX=Ei)、EMAXに対応する変換パラメータtMAXを更新する(tMAX=ti)(S1113)。また、評価値Eiが評価値の最大値EMAX以下(Ei≦EMAX)の場合は更新を行わない。
次に、評価部114は、ループ変数iをインクリメントし(S1114)、ループ変数iと仮説の総数Ntを比較する(S1115)。ループ変数iが仮説の総数Nt以下(i≦Nt)の場合は処理をステップS1105に戻し、ループ変数iが仮説の総数Ntを超える(i>Nt)の場合は処理をステップS1116に進める。つまり、仮説の総数Nt分、ステップS1105からS1115の処理が繰り返される。
仮説の総数Nt分の処理が終了すると、評価部114は、評価値の最大値EMAXに対応する変換パラメータtMAX(つまりθMAXとpMAX)を変形画像生成部110に出力する(S1116)。言い換えれば、複数の変換パラメータから、評価値の最大値EMAXに対応する変換パラメータtMAXが選択される。変形画像生成部110は、変換パラメータtMAXに対応する変形MRI画像ID_MRIonPを画像表示部112に出力する(S1117)。
以上で、剛体変換部106、仮想投影像生成部107、変形推定部109、変形画像生成部110、評価部114による位置・姿勢と圧迫変形の推定(S220)が終了する。この処理によれば、様々な回転と変形を表す変換パラメータtiを仮定して圧迫変形を含む変換を行った結果の中で変形の適切性の評価値Eiが最大になる変換パラメータtMAXによって、変形MRI画像ID_MRIonPが生成される。
このように、MRI画像を圧迫変形して生成した表在血管を強調したMIP像、保持状態の赤外線カメラ画像、PAT画像の血管情報を比較して、MRI画像上の被検体の位置・姿勢と圧迫変形を推定する。これにより、PAT画像の撮像時、理想的な非保持状態の赤外線カメラ画像が撮像されていない場合や、非保持状態から保持状態に移る期間に乳房の位置・姿勢が変化したりする場合も、高精度な変形位置合わせを実現することができる。
[変形例]
上記では、PAT画像と赤外線カメラ画像の血管情報を同じ評価値Eiに組み込んで利用する例を説明した。しかし、これら情報を同時に利用する必要はなく、例えば、PAT画像の血管情報を先に利用し、その後、赤外線カメラ画像の血管情報を利用して、変形位置合わせを実施してもよい。
例えば、変換パラメータtiを変動させ、変形MRI画像とPAT画像との類似度を評価し、評価値が最大になる変換パラメータtMAXを求める。次に、変換パラメータtMAX付近に限定して変換パラメータtiを変動させて、体表近傍変形MIP像IMIPonC1iと赤外線カメラ画像の類似度を評価し、評価値が最大になる変換パラメータtMAX2を求める。そして、変換パラメータtMAX2に対応する変形MRI画像を圧迫変形位置合わせ結果とする。これにより、MRI画像に対して、形態画像として大まかな血管情報が描出されたPAT画像を利用して粗い変形位置合わせを実施した後、機能画像としてより詳細な血管情報が描出された赤外線カメラ画像を用いて精密な変形位置合わせを実施することができる。
また、実施例1と同様に、非保持状態における剛体変換パラメータの推定を行った後に、θjの仮説の設定範囲を剛体変換パラメータの近傍に限定して、保持状態における剛体変換パラメータと変形パラメータを推定する構成でもよい。これにより、仮説の総数が減少し、処理の高速化を図ることができる。勿論、実施例1の変形例1と同様に、仮説の総当りによる評価ではなく、一般的な最適化アルゴリズムを用いて変換パラメータを求める構成でもよい。この場合、非保持状態で推定した剛体変換パラメータをθjの初期値に用いることができる。