JP6289142B2

JP6289142B2 - 画像処理装置、画像処理方法、プログラムおよび記憶媒体

Info

Publication number: JP6289142B2
Application number: JP2014022779A
Authority: JP
Inventors: 和大宮狭; 卓郎宮里
Original assignee: Canon Inc
Current assignee: Canon Inc
Priority date: 2014-02-07
Filing date: 2014-02-07
Publication date: 2018-03-07
Anticipated expiration: 2034-02-07
Also published as: US20150228093A1; JP2015146989A; US9691150B2

Description

本発明は、画像処理装置、画像処理方法、プログラムおよび記憶媒体に関する。特に、本発明は、種々の医用画像収集装置（モダリティ）で撮像した医用画像の画像処理技術に関する。

医用画像収集装置の一種である光音響断層撮像装置（ＰＡＴ（Photoacoustic Tomography）装置；以下、ＰＡＴ装置ともいう。）が、特許文献１に記載されている。ＰＡＴ装置は、被検体に光パルスを照射することで被検体内の吸収物質を励起し、吸収物質の熱弾性膨張により生じる光音響信号を検出することで被検体の光吸収に関わる性質を三次元画像（三次元断層画像）として画像化する装置である。ＰＡＴ装置では、照射光に対する被検体内の光エネルギ堆積量分布（光エネルギ吸収密度分布）が画像化される。また、これに基づき、照射波長に関する被検体の光吸収係数分布が画像化される。さらに、複数の波長に関する光吸収係数分布に基づいて、被検体を構成する物質の状態（例えばヘモグロビンの酸素飽和度など）を画像化することも可能である。

これらの画像は、癌などの悪性腫瘍の内外に生じる新生血管に関する情報を可視化するものであると期待されている。以下では、これらの画像を総称して、光音響断層画像（ＰＡＴ画像）と呼ぶ。

ＰＡＴ装置は、エネルギが小さい近赤外光パルスを照射するため、Ｘ線などと比べて人体の深部の画像化が難しい。そこで、特許文献１のＰＡＴ装置では、乳房を測定対象としたＰＡＴ装置の一形態として、乳房を２枚の平板（以降、保持板と呼ぶ）で保持して乳房の厚さを薄くした状態で撮像を行っている。そのため、ＰＡＴ装置と、磁気共鳴撮像装置（ＭＲＩ）などの他の医用画像収集装置との併用診断を行う場合には、保持による圧迫変形を考慮した変形位置合わせを行う（一方の画像を他方に合うように変形させる）ことで、医師は診断を効率的に行うことが可能になる。

ＰＡＴ画像とＭＲＩ画像の位置合わせの方法として、画像マッチングによる方法が挙げられる。例えば、非特許文献１には、平板により圧迫された乳房を撮像したＸ線マンモグラフィ（ＭＭＧ）と、乳房のＭＲＩ画像との間の位置合わせ技術が記載されている。より具体的には、ＭＲＩ画像に対して平板圧迫による物理変形シミュレーションを施した変形ＭＲＩ画像を生成し、変形ＭＲＩ画像から疑似ＭＭＧ画像を生成し、疑似ＭＭＧ画像と実際に撮像されたＭＭＧとの画像マッチングによる位置合わせを行う技術が開示されている。

また、非特許文献２には、ＭＲＩ画像に対して平板圧迫による物理変形シミュレーションを施した結果から得られる変形後の乳房の形状を、ＭＭＧ画像から抽出した乳房の２次元形状に基づいて評価する技術が開示されている。

特開２０１０−８８６２７号公報

Angela Leeほか「Breast X-ray and MR image fusion using finite element modeling」Proc. Workshop on Breast Image Analysis in conjunction with MICCAI 2011、129-136頁、2011年 C. Tannerほか「Breast Shapes on Real and Simulated Mammograms」Proc. Int. Workshop on Digital Mammography 2010 (IWDM 2010), LNCS 6136、540-547頁、2010年

しかし、医用画像収集装置（例えば、ＰＡＴ装置）により取得される画像と他の医用画像収集装置（例えば、ＭＲＩ）により取得される画像とでは画像化する特徴が異なるため、第一の画像（ＭＲＩ画像）に写る構造物と第二の画像（ＰＡＴ画像）上の構造物すべてが一致するわけではない。そのため、非特許文献１の技術を第二の画像（ＰＡＴ画像）と第一の画像（ＭＲＩ画像）の位置合わせに転用したとしても、画像マッチングのみで高精度な位置合わせを行うのは困難であった。一方、第二の画像（ＰＡＴ画像）には被検体の外形情報が十分に含まれていないため、非特許文献２の技術をそのまま用いるのは困難であった。

本発明は、対象物である被検体のＭＲＩ画像やＣＴ画像などの第一の画像（三次元画像（三次元断層画像））を、第二の画像（ＰＡＴ画像）に対して高精度に位置合わせする技術の提供を目的とする。

本発明の一態様に係る画像処理装置は、第一の撮像手段により対象物を撮像した第一の画像と、第二の撮像手段により前記対象物を撮像した第二の画像と、前記第二の撮像手段に対して位置が対応付けられた第三の撮像手段により前記対象物を撮像した第三の画像と、を取得する画像取得手段と、
前記第一の画像と前記第三の画像との位置合わせを行い、前記第一の画像と前記第三の画像との位置合わせの結果に基づいて、前記第一の画像と前記第二の画像との位置合わせを行う位置合わせ手段と、を備えることを特徴とする。

本発明によれば、対象物である被検体のＭＲＩ画像やＣＴ画像などの第一の画像（三次元画像）を、第二の画像（ＰＡＴ画像）に対して高精度に位置合わせすることが可能になる。

第１の実施形態の画像処理装置を含むモダリティシステムの構成例を示すブロック図。医用画像ＤＢが保持する被検体のＭＲＩ画像を説明する図。ＰＡＴによる被検者の撮像を説明する図。ＰＡＴが撮像したＰＡＴ画像の例を示す図。非保持状態における前面赤外線カメラの撮像画像I_CAM1の一例を示す図。第１の実施形態の画像処理装置の処理を説明する図。表面形状の取得処理を説明する図。変形ＭＲＩ画像とＰＡＴ画像の表示例を示す図。非保持状態における位置合わせの詳細を説明する図。仮想投影像生成部が部分表面領域を求める処理を説明する図。ＭＲＩ画像内の被検体の体表近傍情報を利用したＭＩＰ像を示す図。圧迫変形の推定の詳細を説明する図。メッシュＭの生成方法を示す模式図。保持板による圧迫変形シミュレーションを説明する図。変形ＭＲＩ画像I_{D_MRIonP}を示す模式図。第２の実施形態の画像処理装置を含むモダリティの構成例を示すブロック図。第２の実施形態の画像処理装置の処理手順を説明する図。位置・姿勢と圧迫変形の推定の詳細を説明する図。位置・姿勢と圧迫変形の推定の詳細を説明する図。

以下、図面を参照して、本発明の実施形態を例示的に詳しく説明する。ただし、この実施形態に記載されている構成要素はあくまで例示であり、本発明の技術的範囲は、特許請求の範囲によって確定されるのであって、以下の個別の実施形態によって限定されるわけではない。

（第１の実施形態）
画像処理装置１０は、画像取得部（１０１、１０３、１０４）と、位置合わせ部１１３と、を有する。画像取得部（１０１、１０３、１０４）は、第一の撮像装置（例えば、ＭＲＩ）により対象物を撮像した第一の画像と、第二の撮像装置（例えば、ＰＡＴ装置）により対象物を撮像した第二の画像と、第二の画像撮像装置に対して位置が対応付けられた撮像部により対象物を撮像した第三の画像と、を取得する。

位置合わせ部１１３は、第一の画像における対象物の観測情報と第三の画像における対象物の観測情報とが整合するように第一の画像における対象物と第二の画像における対象物との位置合わせをする。

本実施形態に係る画像処理装置は、乳房を被検体として、光音響断層撮像装置(ＰＡＴ装置)に搭載された赤外線カメラの画像およびＰＡＴ画像と、ＭＲＩ画像とを比較して、ＰＡＴ画像とＭＲＩ画像の変形位置合わせを行う。ＭＲＩ撮像時の被検体の位置や形状を第一の状態、ＰＡＴ撮像時の被検体の位置や形状を第二の状態とした場合に、画像処理装置は、第一の状態の被検体を示すＭＲＩ画像に変形を施し、第二の状態の被検体の画像に対して位置合わせを行う。

具体的な処理として、画像処理装置は、最初に、ＰＡＴ撮像前の非保持状態（以下、第二の状態の初期段階）における被検体をＰＡＴ搭載の赤外線カメラで撮像した二次元画像を得て、ＭＲＩ画像と該二次元画像との位置合わせを行う。すなわち、画像処理装置は、第一の状態と第二の状態（正確には、第二の状態の初期段階）の被検体の位置合わせパラメータとして、両者の間の剛体変換を推定する。画像処理装置は、この剛体変換を初期値として、第一の状態とＰＡＴ撮像時（すなわち、第二の状態）の被検体の位置合わせパラメータとして、圧迫変形の変形パラメータを推定する。そして、画像処理装置は、これら２段階の処理により、ＰＡＴ画像とＭＲＩ画像の間の位置合わせパラメータを導出する。

図１は第１の実施形態に係る画像処理装置１０を含むモダリティシステムの構成例を示す図である。画像処理装置１０の画像取得部は、第一の状態の対象物を第一の撮像装置により撮像した第一の画像を取得する第一の画像取得部（例えば、医用画像取得部１０１）を備える。また、画像処理装置１０の画像取得部は、第一の状態とは異なる第二の状態の対象物を第二の撮像装置により撮像した第二の画像を取得する第二の画像取得部（例えば、ＰＡＴ画像取得部１０３）を備える。また、画像処理装置１０の画像取得部は、第二の状態の初期段階の対象物を撮像部（例えば、赤外線カメラ）により撮像した第三の画像（赤外線カメラ画像）を取得する第三の画像取得部（例えば、カメラ画像取得部１０４）を備える。

位置合わせ部１１３は、第一の画像における対象物の観測情報と第三の画像における対象物の観測情報とが整合するように位置合わせを行った後に、該位置合わせ結果を利用して、第一の画像における対象物の観測情報と第二の画像における対象物の観測情報との比較を行って、第一の画像における対象物と第二の画像における対象物との位置合わせを行う。すなわち、位置合わせ部１１３は、第一の画像と第三の画像との比較により、第一の画像と第三の画像との位置合わせを行った後に、該位置合わせ結果を利用して、第一の画像と第二の画像との比較により第一の画像と第二の画像との位置合わせを行う。

また、画像処理装置１０は、医用画像データベース(ＤＢ)１１と光音響断層撮像装置(ＰＡＴ装置)１２に接続する。医用画像ＤＢ１１は、被検体を事前にＭＲＩで撮像した三次元画像データを保持する。ＰＡＴ１２は、ＰＡＴ画像を撮像する装置であり、被検体のＰＡＴ画像と赤外線カメラ画像を保持する。なお、医用画像ＤＢ１１が保持する三次元画像データは、被検体をＭＲＩで撮像したものに限定されるものではない。ＰＡＴ画像との対比を行いたい三次元画像データあればモダリティを限定するものではなく、例えば、Ｘ線ＣＴ装置で被検体を撮像したＸ線ＣＴ画像であってもよい。また、同じＰＡＴ装置で過去に被検体を撮像したＰＡＴ画像であってもよい。

図２により医用画像ＤＢ１１が保持する被検体のＭＲＩ画像を説明する。図２（ａ）に示す被検体のＭＲＩ画像２００は、人体の頭尾方向に垂直な断面（アキシャル断面）でスライスした二次元画像（乳頭２０４を含む断面）の集合（三次元画像データ）である。ＭＲＩ画像２００を構成する画素は、ＭＲＩ画像座標系Ｃ_MRIにおける位置が定義されている。また、ＭＲＩ画像２００は、被検体の体外領域２０２および被検体の体内領域２０３の撮像結果を含む。

図２（ｂ）に示すＭＲＩ画像２００は、人体の左右方向に垂直な断面（サジタル断面）でスライスした二次元画像の集合（三次元画像データ）である。ＭＲＩ画像２００は、ＭＲＩ画像２００と同様に、被検体の体外領域２０２および体内領域２０３の撮像結果を含む。なお、本実施形態においては、患者の右手側から左手側をx軸の正方向、胸側から背側をy軸の正方向、足側から頭側をz軸の正方向とするＭＲＩ画像座標系Ｃ_MRIを定義する。

図３によりＰＡＴ１２による被検者の撮像を説明する。被検者３００は、ＰＡＴ１２の上面のベッドに伏臥位の体位をとる。そして、被検体である片方の乳房３０１をＰＡＴ１２上面の開口部３０２に挿入する。このとき、照射光が乳房の内部まで届くように、乳房３０１は透明な二枚の保持板（足側の固定保持板３０３と頭部側の可動保持板３０４）により圧迫された状態で保持され、乳房３０１の厚みが薄くなった状態で撮像される。保持は、可動保持板３０４を足方向（固定保持板３０３方向）に移動することで行われる。

固定保持板３０３と可動保持板３０４は何れも平板であり、乳房３０１との接触面（以下、保持面）は平面とする。また、乳房３０１を保持する際の固定保持板３０３と可動保持板３０４との間の距離（以下、保持厚）がＰＡＴ１２によって計測され、保持厚は、ＰＡＴ画像の付加的情報として当該画像のヘッダ部に保存される。

照射光である近赤外光パルスが、保持板の平面に対して直交する方向から図示しない光源によって照射される。そして、被検体内で発生した光音響信号が、保持板の平面に直交するように配置された図示しない超音波探触子によって受信される。

ＰＡＴ１２にはＰＡＴ装置座標系Ｃ_DEVが定義されている。x-y面は、固定保持板３０３および可動保持板３０４の平面に平行な面であり、z軸は、保持された乳房３０１の厚み（保持厚）の方向である。例えば、ＭＲＩ画像座標系Ｃ_MRIと同様に、被検者３００の右手側から左手側をx軸の正方向、胸側（下）から背側（上）をy軸の正方向、足側から頭側をz軸の正方向と定義する。ＰＡＴ装置座標系Ｃ_DEVの原点は、例えば、固定保持板３０３上の右手側の下端位置に設定される。以降、ＰＡＴ１２において、上記の座標系を基準にして、他の座標系との関係を扱うものとする。

図４によりＰＡＴ１２が撮像したＰＡＴ画像を例示的に説明する。ＰＡＴ画像４００は、図２（ａ）に示したＭＲＩ画像と同様にアキシャル断面の二次元画像の集合（三次元画像データ）である。本実施形態では、ＭＲＩ画像座標系Ｃ_MRIと同様に、被検者３００の右手側から左手側をx軸の正方向、胸側から背側をy軸の正方向、足側から頭側をz軸の正方向とするＰＡＴ画像座標系Ｃ_PATを定義する。

ＰＡＴ画像座標系Ｃ_PATからＰＡＴ装置座標系Ｃ_DEVへの変換を行う座標変換行列を「Ｔ_PtoD」と定義する。なお、Ｔ_PtoDを含め、以降に登場する座標変換行列は、すべて座標系の並進と回転を表す４×４行列とする。ＰＡＴ画像座標系Ｃ_PATは、ＰＡＴ装置座標系Ｃ_DEVに対して平行移動した平行座標系であり、被検体である乳房３０１の撮像範囲に応じてＣ_PATの原点位置が平行移動により変化する。つまり、座標変換行列Ｔ_PtoDには回転成分がなく、撮像範囲に基づき一意に算出できる。座標変換行列Ｔ_PtoDはＰＡＴ画像の付加的情報としてＰＡＴ画像のヘッダ部分に保存される。

図３に示すように、ＰＡＴ１２には、被検体である乳房３０１の外観と体表近傍の血管の様子を撮像するための三台の赤外線カメラ（前面赤外線カメラ３０５、後面赤外線カメラ３０６、側面赤外線カメラ３０７）が搭載されている。前面赤外線カメラ３０５は、可動保持板３０４を通して頭部側から乳房３０１の外観を撮像可能な位置に配置されている。後面赤外線カメラ３０６は、固定保持板３０３を通して足側から被検体である乳房３０１の外観を撮像可能な位置に配置されている。側面赤外線カメラ３０７は、側面から乳房３０１の外観を撮像可能な位置に配置されている。

ＰＡＴ１２は、前面赤外線カメラ３０５、後面赤外線カメラ３０６および側面赤外線カメラ３０７によって撮像される、被検体である乳房３０１を保持していない状態（以下、非保持状態）における乳房３０１の画像と、乳房３０１を保持した状態（以下、保持状態）における乳房３０１の画像と、を保存する機能を有する。

以下、保持状態において、前面赤外線カメラ３０５が撮像する画像をＩ_CAM1、後面赤外線カメラ３０６が撮像する画像をＩ_CAM2、側面赤外線カメラ３０７が撮像する画像をＩ_CAM3とする。また、非保持状態において、前面赤外線カメラ３０５が撮像する画像をＩ'_CAM1、後面赤外線カメラ３０６が撮像する画像をＩ'_CAM2、側面赤外線カメラ３０７が撮像する画像をＩ'_CAM3とする。

前面赤外線カメラ３０５の座標系（前面カメラ座標系）Ｃ_CAM1のz軸（視軸の負の方向を示す）は、ＰＡＴ装置座標系Ｃ_DEVのz軸と略同一方向に設定される。同様に、後面赤外線カメラ３０６の座標系（後面カメラ座標系）Ｃ_CAM2のz軸は、ＰＡＴ装置座標系Ｃ_DEVのz軸と略反対方向に設定される。また、側面赤外線カメラ３０７の座標系（側面カメラ座標系）Ｃ_CAM3のz軸は、ＰＡＴ装置座標系Ｃ_DEVの-x軸方向に設定される。

カメラ座標系Ｃ_CAM1、Ｃ_CAM2、Ｃ_CAM3からＰＡＴ装置座標系Ｃ_DEVへの座標変換行列をそれぞれ、Ｔ_C1toD、Ｔ_C2toD、Ｔ_C3toDと定義する。前面赤外線カメラ３０５から側面赤外線カメラ３０７は、ＰＡＴ装置座標系Ｃ_DEVにおいて校正済みである（言い替えれば、ＰＡＴ１２との位置関係が対応付けられている）。また、上記の座標変換行列や前面赤外線カメラ３０５から側面赤外線カメラ３０７の内部パラメータは、既知の情報として、画像処理装置１０に保持されている。

図５により非保持状態における前面赤外線カメラ３０５の撮像画像Ｉ_CAM1を例示的に説明する。赤外線カメラは、近赤外線の強度情報が可視化された画像を撮像する装置であり、近赤外線の以下の性質により、被検体の表層部である皮膚直下の静脈血管（表在血管）を可視化する性質がある。

・近赤外線によって皮膚がある程度透ける
・近赤外線がヘモグロビンを含む血管部分に吸収され、血管が周囲よりも暗く写る
赤外線カメラ画像５００は、皮膚下の表在血管の形状を明瞭に描出するため、形態画像として扱うことができる。図５において、赤外線カメラ画像５００上には、保持状態の乳房輪郭形状５０１が写り、さらに乳頭５０２、そして表在血管５０３が写っている。

前面赤外線カメラ３０５の撮像画像Ｉ_CAM1の二次元座標系Ｃ_IMG1上の座標は、三次元のカメラ座標系Ｃ_CAM1においては、原点である焦点位置と、三次元空間中のカメラの投影面上の点を通る直線、つまり視線と一対一の関係を有する。前面赤外線カメラ３０５の撮像画像の座標系Ｃ_IMG1とカメラ座標系Ｃ_CAM1の間の変換については、一般的な撮像画像と三次元空間の間の座標変換方法を用いるため説明は省略する。また、後面赤外線カメラ３０６と側面赤外線カメラ３０７の撮像画像は、視点位置が異なる点を除き、前面赤外線カメラ３０５の撮像画像と同様であるから説明を省略する。

画像処理装置の処理手順を図６のフローチャートにより説明する。まず、Ｓ６０１で、医用画像取得部１０１は、医用画像ＤＢ１１が保持する被検体のＭＲＩ画像を取得し、ＭＲＩ画像を三次元形状取得部１０２、剛体変換部１０６、変形画像生成部１１０に出力する。

Ｓ６０２で、三次元形状取得部１０２は、入力されたＭＲＩ画像に画像処理を施して、被検体の外形形状情報（表面形状情報）として以下の情報を取得する。すなわち、三次元形状取得部１０２は、被検体の表面に当る各画素の位置（表面位置）を検出し、被検体の表面形状を示す情報を取得する。さらに、Ｓ６０３で、三次元形状取得部１０２は、検出した表面位置から得られる形状の三次元曲率に基づき、ＭＲＩ画像中の特徴点の位置を取得する。被検体が乳房の場合、ＭＲＩ画像中の特徴点は乳頭であり、以下では、三次元形状取得部１０２が乳頭位置を示す情報を取得する情報取得部として説明を続ける。三次元形状取得部１０２は、取得した表面形状と乳頭位置とを剛体変換部１０６、変形推定部１０９に出力する。本実施形態において、ＭＲＩ画像から取得される被検体の表面形状は、非保持状態における被検体の形状モデルになる。

図７は三次元形状取得部１０２による表面形状の取得処理を説明する図である。図７（ａ）、（ｂ）は、図２（ａ）、（ｂ）に示したＭＲＩ画像２００、３０１から被検体の体外領域２０２と体内領域２０３との境界７０４（表面位置）を検出した表面検出画像７００、７０１を示す。表面検出画像７００、７０１は、例えば、被検体の表面とそれ以外とを区別可能な二値画像などでよい。

本実施形態では、被検体の表面形状として、Ｎ_S個の点群Ｐ_Sk(1≦k≦N_S)を取得し、それら点群の位置をＭＲＩ画像座標系Ｃ_MRIにおける三次元の位置座標ベクトルＶ_{Sk_MRI}として記録する。

Ｓ６０４で、ＰＡＴ画像取得部１０３は、ＰＡＴ１２が被検体を撮像したＰＡＴ画像を取得し、ＰＡＴ画像を変形画像評価部１１１、画像表示部１１２に出力する。また、ＰＡＴ画像のヘッダ部に含まれる付加的情報、例えば、座標変換行列Ｔ_PtoD、Ｔ_C1toD、Ｔ_C2toD、Ｔ_C3toDなどもＰＡＴ画像取得部１０３により変形画像評価部１１１に出力される。なお、ＰＡＴ画像取得部１０３が取得するＰＡＴ画像は、所定の波長に対して被検体内の光エネルギ堆積量分布を画像化した三次元画像とする。

Ｓ６０５で、カメラ画像取得部１０４は、ＰＡＴ１２の前面赤外線カメラ３０５、後面赤外線カメラ３０６および側面赤外線カメラ３０７が撮像した、非保持状態における被検体の赤外線カメラ画像と保持状態における被検体の赤外線カメラ画像とを取得する。そして、カメラ画像取得部１０４は、取得した赤外線カメラ画像を二次元形状取得部１０５、仮想投影像評価部１０８に出力する。ここで取得した赤外線カメラ画像は、保持状態における撮像画像Ｉ_CAMi、および非保持状態における撮像画像I'_CAMi(i=1, 2, 3)である。

ＰＡＴ画像取得部１０３およびカメラ画像取得部１０４は、ＰＡＴ１２の撮像に同期してＰＡＴ１２から直接画像を取得してもよいし、過去に撮像され記録された画像を図示しない医用画像記録装置から取得してもよい。ＰＡＴ画像取得部１０３、カメラ画像取得部１０４および先に説明した医用画像取得部１０１は画像を取得する画像取得部として機能する。

Ｓ６０６で、二次元形状取得部１０５は、入力された各赤外線カメラ画像に画像処理を施して、被検体の外形形状情報（表面形状情報）として、乳房輪郭形状（図５の５０１）と乳頭画像（図５の５０２）とを取得する。そして、二次元形状取得部１０５は、乳房輪郭形状および乳頭位置を剛体変換部１０６、変形画像評価部１１１に出力する。例えば、乳房輪郭形状の取得（検出）には一般的なエッジ検出手法を用いることができる。また、乳頭位置は、乳房領域の境界を表す曲線の曲率に基づき取得（検出）することができる。なお、乳房輪郭形状と乳頭位置の取得（検出）は、上記の方法に限られるわけではない。

Ｓ６０７で、剛体変換部１０６、仮想投影像生成部１０７、および仮想投影像評価部１０８は、ＭＲＩ画像と、非保持状態における赤外線カメラ画像に写った表在血管およびその周辺領域を含む線状領域の情報とを用いて、非保持状態における被検体とＭＲＩ画像との位置合わせを行う。ここで、非保持状態とは第二の状態（ＰＡＴ撮像時の被検体の位置や形状）の初期段階である。具体的には、剛体変換部１０６、仮想投影像生成部１０７、および仮想投影像評価部１０８は、仮定した位置合わせパラメータの候補値のそれぞれに基づき、赤外線カメラでＭＲＩ画像を仮想的に観測した仮想画像を生成し、これを赤外線カメラ画像と比較することで位置合わせパラメータを推定する。この位置合わせ（非保持状態における位置合わせ）の詳細は後述する。本ステップの処理により、ＭＲＩ画像座標Ｃ_MRIから前面カメラ座標系Ｃ_CAM1への剛体変換を表す座標変換行列Ｔ_MtoC1が、位置合わせパラメータとして取得される。

Ｓ６０８で、剛体変換部１０６は、座標変換行列Ｔ_MtoC1に基づき、ＭＲＩ画像座標系Ｃ_MRIからＰＡＴ装置座標系Ｃ_DEVへの剛体変換を表す座標変換行列Ｔ_MtoDを算出する。剛体変換部１０６は、座標変換行列Ｔ_MtoC1に、画像処理装置１０が保持する前面カメラ座標系Ｃ_CAM1からＰＡＴ装置座標系C_DEVへの座標変換行列Ｔ_C1toDを掛けることで、座標変換行列Ｔ_MtoDを算出する。

Ｓ６０９で、変形推定部１０９、変形画像生成部１１０、変形画像評価部１１１は、非保持状態における位置合わせ結果に基づき、保持状態における被検体とＭＲＩ画像の位置合わせ（以下、圧迫変形の推定）を行う。詳細は後述するが、物理変形シミュレーションを用いてＭＲＩ画像の圧迫変形を推定する。つまり、変形パラメータを様々に変化させて物理変形シミュレーションを行い、ＰＡＴ画像との比較によって変形の適切性を表す所定の評価値を得て、評価値を最小とする変形パラメータを位置合わせパラメータとして推定する。そして、推定した変形パラメータを用いて、変形ＭＲＩ画像I_{D_MRIonP}（変形三次元画像）を推定結果として生成する。

Ｓ６１０で、画像表示部１１２は、生成された変形三次元画像（変形ＭＲＩ画像）と、ステップＳ６０３で取得されたＰＡＴ画像とを図示しない表示部に並べて表示する。図８は、変形ＭＲＩ画像およびＰＡＴ画像の表示例を示す図である。図８の例は、同一のアキシャル断面の変形ＭＲＩ画像８００とＰＡＴ画像４００を上下に並置した例である。変形ＭＲＩ画像８００に重畳された破線の矩形で示す対応領域８５０は、ＰＡＴ画像４００の表示領域に対応する領域をユーザに示す表示情報である。

なお、三次元形状取得部１０２、二次元形状取得部１０５、剛体変換部１０６、仮想投影像生成部１０７、仮想投影像評価部１０８、変形推定部１０９、変形画像生成部１１０、変形画像評価部１１１は、位置合わせ部１１３を構成する。画像処理装置の画像取得部（１０１、１０３、１０４）は、第一の撮像装置（例えば、ＭＲＩやＣＴ）により対象物を撮像した第一の画像（三次元画像）と、第一の撮像装置とは異なる第二の撮像装置１２０（光音響断層撮像装置）により対象物を撮像した第二の画像（光音響断層画像）と、第二の画像撮像装置に対する位置関係が対応付けられている撮像部（３０５、３０６、３０７）により対象物を撮像した第三の画像（赤外線カメラ画像）と、を取得する。位置合わせ部１１３は、第一の画像と第三の画像との間における対象物の形状情報の比較と、形状情報の比較結果を用いた第一の画像と第二の画像（光音響断層画像）との間における対象物の内部情報の比較とにより、第一の画像と第二の画像との位置合わせを行う。

（非保持状態における位置合わせ処理）
次に、Ｓ６０７で説明した非保持状態における位置合わせ処理を具体的に説明する。非保持状態における位置合わせ処理において、ＭＲＩ画像座標系Ｃ_MRIから前面カメラ座標系Ｃ_CAM1への剛体変換を推定する。図９は、非保持状態における位置合わせ(Ｓ６０７)の具体的な処理の流れを説明する図である。

Ｓ９０１で、剛体変換部１０６は、ＭＲＩ画像を赤外線カメラ座標系（前面カメラ座標系）に平行移動するためのパラメータを算出する。まず、非保持状態における各赤外線カメラ画像から得た二次元の乳頭位置から、三角測量の原理に基づき、剛体変換部１０６は、前面カメラ座標系Ｃ_CAM1における三次元の乳頭位置を算出する。そして、剛体変換部１０６は、ＭＲＩ画像中の乳頭位置が、赤外線カメラ画像から求めた非保持状態における三次元の乳頭位置に一致するように、ＭＲＩ画像座標系Ｃ_MRIから前面カメラ座標系Ｃ_CAM1への平行移動を表す座標変換行列Ｔ１_MtoC1を算出する。

Ｓ９０２で、剛体変換部１０６は、ＭＲＩ画像中の被検体を回転移動するための各成分（三軸回りの回転角度）が取り得る値を組み合わせた複数（Ｎ_θ組）の回転パラメータの候補値（仮説）θi={θx, θy, θz}(1≦i≦n_θ)を設定する。これは、言い換えると、本処理における回転パラメータの候補値θiと、ステップＳ９０１で算出した平行移動のパラメータとを組み合わせた、剛体変換のパラメータの候補値を設定したことになる。また、ＰＡＴ画像座標系と前面カメラ座標系との関係が既知であることを考えると、これは、ＭＲＩ画像（第一の状態の被検体）からＰＡＴ画像（第二の状態の被検体）への剛体変換の候補値を設定することとも等価である。例えば、x軸周り回転角度をθx、z軸回りの回転角度θzとすると、-10度から+10度の範囲で五度刻みに下記の５つの角度を設定する。

θx = {-10, -5, 0, +5, +10}
θz = {-10, -5, 0, +5, +10}
また、y軸回りの回転角度をθyとすると、-180度から+180度の範囲で五度刻みに下記の７２の角度を設定する。

θy = {-180, -175, …, -5, 0, +5, …, +175, +180}
このとき、回転パラメータθiの取り得る値の数（候補値（仮説）の総数）N_θは1800（＝５×５×７２）になる（つまり、1≦i≦1800）。

次に、Ｓ９０３で、剛体変換部１０６は初期化を行う。つまり、剛体変換部１０６は、ループ変数iに１を、後述する類似度Ｓｉの最大値Ｓ_MAXに０、後述する角度θ_MAXにθ₁を設定する。

次に、Ｓ９０４で、剛体変換部１０６は、平行移動後のＭＲＩ画像を、乳頭位置を基準に回転パラメータθiだけ回転移動したＭＲＩ画像I_MRIonC1iと位置座標ベクトルv_{Sk_MRIonC1}iとを仮想投影像生成部１０７に出力する。つまり、剛体変換部１０６は、前面カメラ座標系Ｃ_CAM1において回転パラメータθiだけ回転移動させる座標変換行列Ｔ２ｉを算出する。続いて、剛体変換部１０６は、ステップＳ９０１で導出した平行移動の座標変換行列Ｔ１_MtoC1に座標変換行列Ｔ２ｉを掛けた剛体変換を表す座標変換行列Ｔ_MtoC1iを導出する。そして、剛体変換部１０６は、ＭＲＩ画像および被検体の表面形状の位置座標ベクトルv_{Sk_MRI}を、それぞれ座標変換行列Ｔ_MtoC1iを用いて座標変換したＭＲＩ画像I_MRIonC1iと位置座標ベクトルv_{Sk_MRIonC1}iとを生成する。

次に、Ｓ９０５で、仮想投影像生成部１０７は、剛体変換されたＭＲＩ画像を前面赤外線カメラ３０５の視点から観測した場合に視野に入るであろう被検体の表面形状を部分表面領域として求める。言い換えれば、仮想投影像生成部１０７は、表面位置を示す情報である表面形状と、位置合わせパラメータの候補値である剛体変換パラメータ（剛体変換を表す座標変換行列Ｔ_MtoC1i）に基づき、剛体変換されたＭＲＩ画像を前面赤外線カメラ３０５の視点から見た仮想画像を生成する。

図１０は仮想投影像生成部１０７が部分表面領域を求める処理(Ｓ９０５)を例示的に説明する図である。図１０は前面赤外線カメラ３０５の視点Ｐから観測されるＭＲＩ画像上の被検体の様子を示し、ＭＲＩ乳房領域１０００は、乳頭位置１００１が赤外線カメラにおける乳頭位置に一致するように剛体変換されている。ステップＳ９０５において、仮想投影像生成部１０７は、視点Ｐの位置・姿勢を基準に透視投影を行い、視点Ｐから観測される部分表面領域１００３を求める。仮想投影像生成部１０７は、視点Ｐからの投影線１００４を延長し、投影線１００４がＭＲＩ乳房領域１０００の体表形状を表す位置座標ベクトルv_{Sk_MRIonC1}iと交点をもつ観測範囲１００２を特定する。そして、仮想投影像生成部１０７は、観測範囲１００２において、位置座標ベクトルv_{Sk_MRIonC1}iの中から、投影線１００４が視点Ｐを起点にして最初に交わる体表点１００５を求める。そして、仮想投影像生成部１０７は、観測範囲１００２に含まれるすべての投影線１００４について体表点１００５を求めて部分表面領域１００３を決定する。

次に、仮想投影像生成部１０７は、剛体変換されたＭＲＩ画像における部分表面領域１００３の近傍情報を利用したＭＩＰ像を生成し、生成したＭＩＰ像を仮想投影像評価部１０８に出力する。ここで、ＭＩＰは「maximum intensity projection」の略であり、以下では、ステップＳ９０６で生成されるＭＩＰ像を「体表近傍ＭＩＰ像I_MIPonC1i」と呼ぶ。

仮想投影像生成部１０７は、投影線１００４が最初に交わる体表点１００５を起点にして、視点Ｐから遠ざかる方向に所定の距離（例えば、５ｍｍ）をもつ体表近傍区間１００６を設定する。そして、仮想投影像生成部１０７は、観測範囲１００２に含まれるすべての投影線１００４について体表近傍区間１００６を設定することで体表近傍領域１００７を規定する。続いて、仮想投影像生成部１０７は、視点Ｐの位置・姿勢を基準に透視投影を行い、ＭＲＩ画像I_MRIonC1iの中で体表近傍領域１００７に含まれる領域に限定したＭＩＰ像である体表近傍ＭＩＰ像I_MIPonC1iを生成する。これにより、ＭＲＩ画像中の血管領域およびその周辺領域を含む線状領域（以下、血管領域１００８）の中で前面赤外線カメラ３０５側の体表近傍に存在する表在血管およびその周辺領域を含む線状領域（以下、表在血管１００９）の情報のみを可視化したＭＩＰ像が生成される。なお、体表近傍領域１００７の生成時に、被検体の皮膚が当該領域に含まれないようにしてもよい。具体的には、部分表面領域１００３に皮膚厚相当の所定の厚みを与えた領域を皮膚領域として導出し、上記の処理で得た体表近傍領域１００７からこれを除外してもよい。被検体の皮膚はＭＲＩ画像中において高い輝度値を有しているので、この除外処理により、生成するＭＩＰ像中において表在血管をより明瞭に描出できるようになる。

このように、体表近傍、かつ、前面赤外線カメラ３０５の視点Ｐから観測されるだろう領域の情報のみを可視化することで、乳房内部の血管領域や、前面赤外線カメラ３０５からみて反対側の体表近傍に存在する血管領域の可視化を防ぐことができる。つまり、実際の赤外線カメラ画像により近いＭＩＰ像である体表近傍ＭＩＰ像I_MIPonC1iを生成することができる。

なお、体表近傍ＭＩＰ像I_MIPonC1iの生成方法は、体表近傍の領域のみを、あるいは、体表近傍の領域を強調して可視化することができれば、上記の方法に限られるものではない。他の方法として、投影線１００４上において、体表点１００５からの距離が、視点Ｐから遠ざかる方向に大きくなるほど、輝度値に対する重みを小さくしてＭＩＰ像を生成する方法が挙げられる。この方法によれば、乳房内部の領域で、体表に近い領域ほど輝度値が強調されたＭＩＰ像が生成されるため、表在血管が強調されたＭＩＰ像が生成される。もちろん、本処理によるＭＩＰ像の生成を行う場合においても、皮膚領域を描画対象から除外することで、表在血管をより明瞭に描出できる。

図１１は、ＭＲＩ画像内の被検体の体表近傍情報を利用したＭＩＰ像１１０３を例示的に示す図である。図１１において、体表近傍ＭＩＰ像Ｉ_MIPonC1iは、前面カメラ画像座標系Ｃ_IMG1上の二次元画像として表される。これは、視点Ｐに基づく透視投影によって二次元のＭＩＰ像を生成する際の座標変換(Ｓ９０６)が、三次元の前面カメラ座標系Ｃ_CAM1から二次元の前面カメラ画像座標系Ｃ_IMG1への座標変換（カメラの撮像プロセス）と幾何的に等しいためである。

また、ＭＩＰ像の生成の際、体表の外側領域は処理対象に含まれないため、乳房輪郭形状１１００の外側には有意の輝度値は存在しない。さらに、ＭＲＩ画像における乳頭位置１００１が、前面赤外線カメラ３０５の撮像画像Ｉ_CAM1における乳頭５０２の位置に一致するように剛体変換されているため、ＭＩＰ像の乳頭１１０１は、撮像画像Ｉ_CAM1（図５）の乳頭５０２に一致する。そして、ＭＩＰ像においては、表在血管１１０２が特に高輝度な領域として可視化される。

次に、Ｓ９０７で、仮想投影像評価部１０８は、体表近傍ＭＩＰ像Ｉ_MIPonC1iと、非保持状態の赤外線カメラ画像Ｉ'_CAM1の両方において可視化された表在血管の輝度情報に基づき、両画像の類似度Ｓiを算出する。仮想投影像評価部１０８は、体表近傍ＭＩＰ像Ｉ_MIPonC1i、非保持状態の赤外線カメラ画像Ｉ'_CAM1ともに、類似度の算出領域から乳房の外側領域を除外し、類似度の算出領域を乳房の内部領域に限定する。

体表近傍ＭＩＰ像I_MIPonC1iにおける表在血管１１０２（図１１）は、周囲の乳房領域に比べて高輝度に可視化されている。他方、赤外線カメラ画像Ｉ'_CAM1における表在血管５０３（図５）は、周囲の乳房領域に比べて低輝度に可視化されている。そこで、仮想投影像評価部１０８は、それら画像の間の輝度情報を直接比較できるようにするため、体表近傍ＭＩＰ像Ｉ_MIPonC1iの輝度値を反転する。そして、仮想投影像評価部１０８は、輝度値を反転した体表近傍ＭＩＰ像と赤外線カメラ画像Ｉ'_CAM1の間で類似度Ｓｉ（0≦Si≦1）を算出する。体表近傍ＭＩＰ像Ｉ_MIPonC1iの表在血管１１０２（図１１）と、赤外線カメラ画像Ｉ'_CAM1の表在血管５０３（図５）とが類似すれば、類似度Ｓｉの値が高くなる（１に近付く）ものとする。

なお、本実施形態においては、類似度Ｓｉの評価尺度として画像間の相互情報量を適用するが、評価尺度はこれに限られず、相互相関係数やＳＳＤ (Sum of Squared Difference）など、既知の手法を用いてもよい。また、輝度値に直接的に基づく評価尺度である必要はなく、例えば両画像からエッジなどの画像特徴を検出して、それらの類似度や一致度を算出するような尺度でもよい。

Ｓ９０８で、仮想投影像評価部１０８は、類似度Ｓｉと類似度の最大値Ｓ_MAXとを比較する。そして、類似度Ｓｉが最大値Ｓ_MAXを超える（Si＞S_MAX）場合（Ｓ９０８−Ｙｅｓ）、処理をステップＳ９０９に進める。

Ｓ９０９で、仮想投影像評価部１０８は、Ｓ_MAXを更新し（S_MAX=Si）、Ｓ_MAXに対応する角度θ_MAXを更新する（θ_MAX=θi）。一方、ステップＳ９０８の判定で、類似度Ｓｉが最大値S_MAX以下（Si≦S_MAX）の場合（Ｓ９０８−Ｎｏ）、仮想投影像評価部１０８はステップＳ９０９の更新を行わずに処理をステップＳ９１０に進める。

Ｓ９１０で、仮想投影像評価部１０８は、ループ変数iをインクリメントし、ループ変数iと仮説の総数Ｎ_θとを比較する(Ｓ９１１)。ループ変数iが仮説の総数N_θ以下（i≦N_θ）の場合（Ｓ９１１−Ｎｏ）、仮想投影像評価部１０８は、処理をステップＳ９０４に戻し、ループ変数iが仮説の総数N_θを超える（i＞N_θ）場合（Ｓ９１１−Ｙｅｓ）、仮想投影像評価部１０８は、処理をステップＳ９１２に進める。つまり、仮説の総数N_θ分、ステップＳ９０４からＳ９１１の処理が繰り返される。

仮説の総数Ｎ_θ分の処理が終了すると、Ｓ９１２で、剛体変換部１０６は、角度θ_MAXにおける座標変換行列Ｔ_MtoC1MAXを、ＭＲＩ画像座標Ｃ_MRIから前面カメラ座標系Ｃ_CAM1への剛体変換を表す最終的な座標変換行列Ｔ_MtoC1とする。言い換えれば、複数の回転パラメータから、類似度の最大値Ｓ_MAXに対応する回転パラメータθ_MAXが選択される。

以上で、剛体変換部１０６、仮想投影像生成部１０７、仮想投影像評価部１０８が行う、非保持状態における位置合わせ(Ｓ６０７)の処理が終了する。この処理により、様々な剛体変換パラメータの仮説（つまり回転パラメータの仮説θi）に基づき生成した体表近傍ＭＩＰ像Ｉ_MIPonC1iが、赤外線カメラ画像Ｉ'_CAM1に最も類似する角度θ_MAXに基づく座標変換行列Ｔ_MtoC1が得られる。

上記の処理では、前面カメラ画像座標系Ｃ_IMG1に変換した体表近傍ＭＩＰ像Ｉ_MIPonC1iを生成し、前面赤外線カメラ３０５の画像Ｉ'_CAM1との類似度Ｓｉを評価する例を説明した。しかし、類似度Ｓｉの評価対象は、前面赤外線カメラ３０５に限られるわけではない。

例えば、後面カメラ画像座標系Ｃ_IMG2や側面カメラ画像座標系Ｃ_IMG3に変換した体表近傍ＭＩＰ像Ｉ_MIPonC1iを生成する。そして、体表近傍ＭＩＰ像Ｉ_MIPonC1iと後面赤外線カメラ３０６の画像Ｉ'_CAM2、または側面赤外線カメラ３０７の画像Ｉ'_CAM3との類似度を評価してもよい。その場合、ステップＳ９０４で得た座標変換行列の候補値Ｔ_MtoC1iに対して、ステップＳ９０５、Ｓ９０６で透視投影を行う視点Ｐを、後面赤外線カメラ３０６のカメラ視点または側面赤外線カメラ３０７のカメラ視点にそれぞれ置き換えればよい。

ここで、ＭＲＩ画像I_MRIonC1iおよび位置座標ベクトルv_{Sk_MRIonC1}iはステップＳ９０４で前面カメラ画像座標系C_IMG1に変換されている。従って、この透視投影に用いるカメラ視点の位置・姿勢としては、前面カメラ座標系Ｃ_CAM1上で表現される位置姿勢を用いればよい。

上述したように、各赤外線カメラ（前面、側面、後面）の位置関係は、ＰＡＴ装置座標系Ｃ_DEVを基準に対応付けられている。そのため、前面カメラ座標系Ｃ_CAM1から後面カメラ座標系Ｃ_CAM2または側面カメラ座標系Ｃ_CAM3への座標変換行列を導出することができる。つまり、前面カメラ座標系Ｃ_CAM1における後面または側面赤外線カメラのカメラ視点の位置姿勢を導出することができる。

前面、後面、側面の赤外線カメラそれぞれに基づく類似度を総合した総合類似度を算出して、類似度を評価してもよい。総合類似度としては、これら三種類の類似度の、例えば、重み付け平均値、最大値、最小値、中央値が挙げられる。また、上述の処理では、前面カメラ座標系Ｃ_CAM1を基準にして、回転パラメータθiの設定、及びＭＲＩ画像を各カメラ座標系に透視投影するための座標変換を行った。しかし、ＰＡＴ装置座標系Ｃ_DEVを基準にして、Ｃ_MRIとC_DEVとの間の剛体変換Ｔ_MtoDを位置合わせパラメータとして、その候補値を設定するようにしてもよい。この場合、Ｃ_MRIからＣ_DEVへの剛体変換に加えて各赤外線カメラへのビューイング変換を行った後に、ステップＳ９０６と同様の処理で体表近傍ＭＩＰ像を生成すればよい。そして、各赤外線画像との類似度評価に基づいて位置合わせパラメータを推定すればよい。

（圧迫変形の推定）
図１２のフローチャートにより圧迫変形の推定(Ｓ６０９)の詳細を説明する。

Ｓ１２０１で、変形推定部１０９は、ステップＳ６０２で取得された被検体の表面形状、および、ステップＳ６０８で取得された座標変換行列Ｔ_MtoDを用いて、被検体の形状を表す三次元のメッシュ（以下、メッシュＭ）を生成する。変形推定部１０９は、座標変換行列Ｔ_MtoDによる座標変換を、被検体の表面形状V_{Sk_MRI}に施して、ＰＡＴ装置座標系C_DEVにおける被検体表面点群の位置座標ベクトルV_{Si_MRIonD}（1≦i≦Ns）を算出する。そして、変形推定部１０９は、位置座標ベクトルV_{Si_MRIonD}が表す表面形状に基づき、被検体の内部領域を判別し、内部領域にメッシュMを配置する。

図１３はメッシュＭの生成方法を例示的に説明する図である。図１３(ａ)は被検体の処理対象領域１３００のサジタル断面を示し、サジタル断面における被検体の表面位置１３０１と内部領域１３０２とを示す。図１３(ｂ)に示すように、メッシュＭは、被検体の内部領域１３０２に六面体や四面体などの立体構造をもつ要素１３０３を配置して生成される。メッシュＭは、これらの要素の頂点（ノード）１３０４の位置と連結情報によって記述される。

以下では、ステップＳ１２０１において配置されるメッシュＭのノード数をＮｍ、各ノードの位置をｓ_L（１≦Ｌ≦Ｎｍ）と表記する。各ノードの変位によって要素内の変位場を表現することができ、これに基づき、被検体内の任意の点の変位を求めることができる。

次に、Ｓ１２０２で、変形推定部１０９は、変形パラメータの各成分（被検体のヤング率やポアソン比など）が取り得る値を組み合わせた複数（Ｎｐ組）の変形パラメータの仮説ｐ_k（１≦ｋ≦Ｎｐ）を生成する。変形推定部１０９は、例えば、各成分が取り得る値の範囲を適切な間隔で分割し、それらすべてを組み合わせることで変形パラメータp_kを生成する。例えば、ヤング率の比p_yとポアソン比p_pを変形パラメータp_kの成分とし、ヤング率の比p_yとポアソン比p_pが取り得る値を次のような値とする。

p_y = {1, 2, 3, 4, 5}
p_p = {0.0, 0.2, 0.4, 0.45, 0.499}
そして、ヤング率の比ｐ_yとポアソン比ｐ_pを組み合わせた変形パラメータを生成する。上記の例ではＮｐ＝２５（＝５×５）（組）になる。なお、ヤング率の比ｐ_yは、乳房の硬さの非等方性に対応するためのパラメータであり、人体のコロナル面（x-z平面）におけるヤング率に対する、人体の前後方向（y軸方向）のヤング率の割合を表す。

Ｓ１２０３で、変形推定部１０９は初期化処理を行う。変形推定部１０９は、ループ変数ｋに１を設定し、後述する評価値の最大値Ｅ_MAXに０を設定し、後述する変形パラメータｐ_MAXにp₁を設定する。

Ｓ１２０４で、変形推定部１０９は、変形パラメータｐ_kを用いて、有限要素法に基づく物理変形シミュレーションをメッシュＭに施し、変形後のメッシュである変形メッシュＤＭｋを生成する。このときの変形関数Ｆｋ(x, y, z)を、メッシュＭから変形メッシュＤＭｋへの各ノードの変位を与える変位ベクトルｄｋ_L（１≦Ｌ≦Ｎｍ）と定義する。

図１４は、ステップＳ１２０４における物理変形シミュレーションとして、保持板による圧迫変形シミュレーションを例示的に説明する図である。保持板による圧迫変形において、二枚の保持板を被検体の中心方向に移動すると、移動後の保持板に接する被検体の表面領域は保持板に押圧された（張り付いた）状態になる。

図１４（ａ）に示すように、二枚の保持板をそれぞれ、距離Δｄ１、Δｄ２分、矢印方向に移動すると仮定する。変形推定部１０９は、メッシュＭのノードのうち、体表面を表す表面ノードから、保持面Ｐ_Ud1、Ｐ_Ld2に接する外側表面ノード１４００と１４０１とを抽出し、外側表面ノード１４００および１４０１のそれぞれが保持面に張り付くための変位量を求める。そして、変形推定部１０９は、変位量を変形シミュレーションの境界条件Ｃとして有限要素法の計算を実行し、二枚の保持板がそれぞれ、距離Δｄ１、Δｄ２分、移動した場合の変形メッシュを生成する。

本実施形態において、変形推定部１０９は、図１４(ｂ)に示す二枚の保持板の最終的な保持位置Ｐ_U、Ｐ_Lまで移動するまでの間、変形過程で生じる境界条件の変化に対応した複数回（Ｎ回）の変形シミュレーションを分割して行う。図１４（ｂ）はＮ回の変形シミュレーションを繰り返した結果の変形メッシュＤＭｋを表す。物理変形シミュレーションにより、図１４（ｂ）に示すように、被検体の変形メッシュが保持位置Ｐ_U、Ｐ_Lの間でｚ軸方向に圧縮され、ｙ軸方向に伸長されることが、図１４（ａ）との比較により分かる。

Ｓ１２０５で、変形画像生成部１１０は、変形パラメータｐ_kに対応する変形をＭＲＩ画像に施して、変形ＭＲＩ画像を生成し、変形ＭＲＩ画像を変形画像評価部１１１に出力する。変形画像生成部１１０は、ＭＲＩ画像を、座標変換行列Ｔ_MtoDを用いてＰＡＴ装置座標系Ｃ_DEVに座標変換し、ステップＳ１２０４で算出した変形関数Ｆｋ(x, y, z)を用いて変形処理する。そして、変形画像生成部１１０は、座標変換行列Ｔ_PtoDの逆行列を用いた座標変換により、ＰＡＴ画像座標系C_PATにおける変形ＭＲＩ画像I_{D_MRIonP}kを生成する。

図１５は、変形ＭＲＩ画像Ｉ_{D_MRIonP}を例示的に示す図であり、図１５には変形ＭＲＩ画像８００、変形後の乳房領域１５０１、変形前の乳房形状１５０２が示されている。図１５（ａ）は変形ＭＲＩ画像Ｉ_{D_MRIonP}をアキシャル断面でスライスした二次元画像であり、図１５（ｂ）は変形ＭＲＩ画像Ｉ_{D_MRIonP}をサジタル断面でスライスした二次元画像である。変形後の乳房領域１５０１と変形前の乳房形状１５０２とを比較すると、ＰＡＴ画像座標系Ｃ_PATのｚ軸方向への圧迫により、ｘ-ｙ平面上において乳房領域が伸長し、ｚ軸方向に圧縮されていることがわかる。

Ｓ１２０６で、変形画像評価部１１１は、ステップＳ６０４で取得されたＰＡＴ画像と、ステップＳ６０５で取得された保持状態の赤外線カメラ画像と、を用いて、変形ＭＲＩ画像Ｉ_{D_MRIonP}kの適切性の評価値Ｅｋを算出する。そして、変形画像評価部１１１は、算出した評価値Ｅｋを変形推定部１０９に出力する。変形画像評価部１１１は、変形ＭＲＩ画像とＰＡＴ画像との間の類似度Ｓ_MRIk（0≦S_MRI、k≦1）、および、変形ＭＲＩ画像と保持状態の赤外線カメラ画像における乳房形状の間の残差Ｒｋに基づき、評価値Ｅｋを算出する。

なお、評価値Ｅｋの値が高いほど変形が適切であるものとする。また、類似度Ｓ_MRIkの評価尺度として、ステップＳ９０７と同様に画像間の相互情報量を適用する。なお、評価尺度はこれに限らず、相互相関係数やＳＳＤ、血管分岐部等の特徴点の位置の一致度など、既知の何れの手法を用いてもよい。

図８は変形ＭＲＩ画像（Ｉ_{D_MRIonP}k）８００とＰＡＴ画像４００の同一のアキシャル断面を示し、破線の矩形で示す対応領域８５０は、ＰＡＴ画像４００に対応する変形ＭＲＩ画像８００上の領域である。類似度Ｓ_MRIkは、ＰＡＴ画像４００と対応領域８５０との間で計算される。ＰＡＴ画像４００と対応領域８５０との間で、可視化された、変形ＭＲＩ画像８００の血管領域８５１とＰＡＴ画像４００の血管領域８５２が類似すれば、類似度Ｓ_MRIkの値が高くなる。

また、残差Ｒｋは、赤外線カメラ画像に写った被検体の輪郭（シルエット）の形状と、赤外線カメラ画像に投影した変形メッシュＤＭｋの外郭の形状との間の差として算出される。例えば、変形メッシュＤＭｋを保持状態の赤外線カメラ画像Ｉ_CAM1へ投影する場合、変形メッシュＤＭｋは、座標変換行列Ｔ_C1toDの逆行列を用いて前面カメラ座標系Ｃ_CAM1に座標変換され、前面カメラ画像座標系Ｃ_IMG1に投影変換される。また、同様の方法により、変形メッシュＤＭｋを後面カメラ画像Ｉ_CAM2または側面カメラ画像I_CAM3へ投影することができる。

残差Ｒｋは、例えば、変形メッシュＤＭｋを三つの赤外線カメラ画像それぞれに投影したメッシュの外郭と、三つの赤外線カメラ画像の間の各残差を総合した総合残差（例えば、三つの残差の重み付け平均値）として算出される。ただし、総合残差は、重み付け平均値に限らず、三つの残差の最大値、最小値、中央値などを用いてもよい。

また、残差Ｒｋは、例えば、赤外線カメラ画像に写った被検体の乳頭位置と、赤外線カメラ画像に投影した変形メッシュＤＭｋの乳頭位置との間の残差として算出してもよい。勿論、乳房形状の残差と、乳頭位置の残差との両方を統合した値（例えば、重み付け和）を残差Ｒｋとしてもよい。なお、乳房形状や乳頭位置に基づく残差は、赤外線カメラ画像から取得する例を説明したが、乳房の外観を撮像した一般的なカメラ画像から取得してもよい。

評価値Ｅｋは、例えば、類似度S_MRIkと残差Ｒｋに基づく重み付け和として、下式で表される。

Ｅｋ = ａＳ_MRIk ＋ｂ{1/(1+Rk)} ・・・(１)
ここで、ａ、ｂは重み係数（ａ＋ｂ=1）。

式(１)の第二項において(１＋Ｒｋ)の逆数を用いるのは次の理由からである。

・残差Ｒｋは、評価値Ｅｋとは逆に、変形が適切であるほど値が小さくなる指標である
・値の取り得る範囲を、類似度Ｓ_MRIkと同様に、０から１の間にする
Ｓ１２０７で、変形推定部１０９は、入力される評価値Ｅｋと評価値の最大値Ｅ_MAXとを比較する。評価値Ｅｋが評価値の最大値Ｅ_MAXを超える（Ek＞E_MAX）場合（Ｓ１２０７−Ｙｅｓ）、変形推定部１０９は、処理をステップＳ１２０８に進める。

Ｓ１２０８で、変形推定部１０９は、Ｅ_MAXを（Ｅ_MAX=Ｅｋ）と更新する。また、変形推定部１０９は、Ｅ_MAXに対応する変形パラメータｐ_MAXを（ｐ_MAX＝ｐ_k）と更新する。

一方、ステップＳ１２０７の比較で、評価値Ｅｋが評価値の最大値Ｅ_MAX以下（Ｅｋ≦Ｅ_MAX）の場合（Ｓ１２０７−Ｎｏ）、変形推定部１０９は、ステップＳ１２０８の更新を行わずに処理をステップＳ１２０９に進める。

次に、Ｓ１２０９で、変形推定部１０９は、ループ変数ｋをインクリメントする。そして、Ｓ１２１０で、変形推定部１０９は、ループ変数ｋと仮説の総数Ｎｐとを比較する。ループ変数ｋが仮説の総数Ｎｐ以下（ｋ≦Ｎｐ）の場合（Ｓ１２１０−Ｎｏ）、変形推定部１０９は、処理をステップＳ１２０４に戻し、ループ変数ｋが仮説の総数Ｎｐを超える場合（ｋ＞Ｎｐ）の場合（Ｓ１２１０−Ｙｅｓ）、変形推定部１０９は、処理をステップＳ１２１１に進める。つまり、仮説の総数Ｎｐ分、ステップＳ１２０４からＳ１２１０の処理が繰り返される。

仮説の総数Ｎｐ分の処理が終了すると、Ｓ１２１１で、変形推定部１０９は、変形パラメータｐ_MAXを変形画像生成部１１０に出力する。言い換えれば、複数の変形パラメータから、評価値の最大値Ｅ_MAXに対応する変形パラメータｐ_MAXが選択される。

Ｓ１２１２で、変形画像生成部１１０は、変形パラメータｐ_MAXに対応する変形ＭＲＩ画像Ｉ_{D_MRIonP}MAXを圧迫変形の推定結果（I_{D_MRIonP}）とし、変形ＭＲＩ画像Ｉ_{D_MRIonP}を画像表示部１１２に出力する。

以上で、変形推定部１０９、変形画像生成部１１０、変形画像評価部１１１による圧迫変形の推定処理(Ｓ６０９)が終了する。この処理によれば、様々な変形パラメータの仮説ｐ_kの下で変形シミュレーションを実行する。そして、それらシミュレーション結果の中で変形の適切性の評価値Ｅｋが最大になる変形パラメータｐ_MAXによって変形ＭＲＩ画像Ｉ_{D_MRIonP}が生成される。

このように、乳房を被検体として、ＰＡＴ１２に搭載された赤外線カメラの画像およびＰＡＴ画像と、ＭＲＩ画像との比較により、ＰＡＴ画像とＭＲＩ画像との位置合わせ精度の向上を図ることができる。つまり、ＭＲＩ画像から生成した表在血管を強調した二次元のＭＩＰ像と、非保持状態の赤外線カメラ画像とを比較して、赤外線カメラに対するＭＲＩ画像上の被検体の位置・姿勢を推定し、ＭＲＩ画像と赤外線カメラ画像との高精度な剛体位置合わせを実現する。さらに、ＭＲＩ画像を赤外線カメラの座標系からＰＡＴ１２の座標系に座標変換することで、剛体位置合わせの結果を、ＭＲＩ画像とＰＡＴ画像との間の変形位置合わせ処理の初期状態として活用する。つまり、ＭＲＩ画像とＰＡＴ画像とを比較して位置合わせを行う際は、非保持状態から保持状態への圧迫変形を推定するだけでよい。そして、赤外線カメラ画像に写った乳房形状と、ＰＡＴ画像に写った内部構造の両方に合致するようにＭＲＩ画像を変形することで、高精度な圧迫変形の推定が可能になる。

なお、本実施形態では、ＭＲＩ画像を変形位置合わせする手法として、有限要素法に基づく物理変形シミュレーションを用いたが、必ずしもこの手法に限られるものではない。例えば、一般的な変形手法である、ＦＦＤ (Free Form Deformation)などの手法を用いてもよい。ＦＦＤを用いた変形処理では、まず画像内の被検体を直方体で取り囲むように、格子状の制御点を配置する。そして、この制御点を動かすことにより、直方体内部に存在する画像領域を変形させることができる。ここで、各制御点の変位量のセットを、変形パラメータの候補値（仮説）ｐ_k(１≦ｋ≦Ｎｐ)と定義する。そして、変形パラメータｐ_kの値を様々に変動させ、上述の変形画像評価部１１１における評価値Ｅｋを最大化するような変形パラメータｐ_kを算出することで、上述の目的に合致する変形位置合わせを実現することができる。

なお、本実施形態では、被検体として人体の乳房を適用したが、これに限られるものではなく、表在血管を有する生体の部位であれば被検体は何であってもよい。また、医用画像データベース１１に登録されている画像としてＭＲＩ画像を適用したが、生体を撮像した三次元画像であれば、何れのモダリティによる画像であってもよい。

（変形例１）
上記の非保持状態における位置合わせ処理（Ｓ６０７）では、赤外線カメラに対するMRI画像内の被検体の位置・姿勢を求めるための回転角度を網羅的に変動させて評価値を算出し、最適な評価値を与える回転角度を取得した。しかし、ＭＲＩ画像内の被検体の位置・姿勢の取得に別の方法を用いることができる。つまり、一般的な最適化アルゴリズムを用いて、評価値が最適になる回転角度を推定してもよい。例えば、最適化アルゴリズムの一種である最急降下法を用いる方法を説明する。

三次元ベクトルをｘ、赤外線カメラに対するＭＲＩ画像内の被検体の回転角度を表すパラメータを(θx, θy, θz)とする。回転角度を表すベクトルｘを与えたときに、ステップＳ９０７で生成される赤外線カメラ画像と体表近傍ＭＩＰ像との類似度をＳ_MIP(x)とする。最急降下法において、ベクトルｘを与えたときに、最小化する関数ｆ(x)を類似度Ｓ_MIP(x)の逆数１/Ｓ_MIP(x)とする。ここで、ｆ(ｘ)をＳ_MIP(x)の逆数とする理由は、類似度Ｓ(x)が最大になる回転角度を表すパラメータを求めるためである。以上のように設定された各変数を、次式を用いて更新し、収束させることでｆ(x)を最小化する（つまりS_MIP(x)を最大化する）パラメータxを算出する。

x^(k+1) = x^(k) - αgrad f(x^(k)) ・・・(２)
ここで、αは一回に更新する数値の割合を決めるパラメータ（通常、小さな正の定数）、
kは更新回数、
grad f(x^(k))はk回目の更新における関数f(x)の勾配ベクトル（関数f(x)の変化率が最も大きい方向を向く）。

勾配ベクトルgrad f(x^(k))は以下の方法で求める。ｋ回目の更新におけるベクトルｘ^(k)=(θx^(k), θy^(k), θz^(k))^Tとする。ベクトルx^(k)の各要素に対して、それぞれ微小な変化量Δx=(Δθx, Δθy, Δθz)を与えたときの、関数ｆ(x^(k)+Δx)を算出する。

そのベクトルの方向がパラメータ空間において均等に変化するように変化量Δxを変動させた場合の、各Δxに対応するf(x^(k)+Δx)を算出する。算出したf(x^(k)+Δx)の集合から、f(x^(k))-f(x^(k)+Δx)が最大になるΔx（Δx_MAX）を求める。このΔx_MAXは、f(x^(k))の変化率が最大になるパラメータ空間の方向ベクトルであり、grad f(x^(k))に等しい。

また、最適化アルゴリズムとして、ニュートン法など既知の手法の何れを用いてもよい。これにより、より少ない繰り返し回数で、最適な評価値を与える回転角度を推定することができ、処理の高速化を図ることができる。

（変形例２）
上記の剛体変換の推定処理（Ｓ６０８）では、ＭＲＩ画像の剛体変換の妥当性を評価する評価値として、赤外線カメラに写った内部情報との類似度を用いた。しかし、別の方法で評価値を求めてもよい。

例えば、ＭＲＩ画像から生成した体表近傍ＭＩＰ像に写る被検体の輪郭（シルエット）形状と赤外線カメラに写った被検体の輪郭形状に対する形状残差を算出し、これを評価値としても良い。

つまり、ステップＳ９０７において、体表近傍ＭＩＰ像I_{D_MIPonC1}k及び赤外線カメラ画像Ｉ'_CAM1の夫々から被検体である２次元の乳房形状を抽出し、それらの残差を評価値Ｅｋとする。このとき、体表近傍ＭＩＰ像I_{D_MIPonC1}kは、生成元のＭＲＩ画像と同様、乳房外部は低輝度であるのに対し、乳房内部はそれよりも高輝度であり、両者の境界は明瞭である。従って、体表近傍ＭＩＰ像I_{D_MIPonC1}kからの乳房形状の抽出は、ステップＳ６０２におけるＭＲＩ画像からの表面形状取得方法と同様の方法で行うことができる。また。赤外線カメラ画像Ｉ'_CAM1からの乳房形状の抽出は、ステップＳ６０６における乳房輪郭形状の抽出方法と同様の方法で行うことができる。さらに、評価値Ｅｋの算出方法は、ステップＳ１２０６における変形ＭＲＩ画像と保持状態の赤外線カメラ画像における乳房形状の残差Ｒｋと同様の方法で算出することができる。

これにより、赤外線カメラに写った内部情報を使わずに、もしくは被検体の内部情報が写らない一般的なカメラを用いた場合でも、ＭＲＩ画像の剛体変換の処理を行うことができる。さらに、ステップＳ１２０６で算出する残差Ｒｋも、赤外線カメラではなく一般的なカメラを用いて算出可能であるため、処理全体を通して赤外線カメラを一般的なカメラに置き換えることができる。また、表在血管を有さない乳房以外の部位を対象物体とすることができる。

なお、赤外線カメラに写った内部情報の類似度と、輪郭形状の誤差を組み合わせた値を評価値としても良い。これは、ステップＳ１２０６における類似度Ｓ_MRIkを赤外線カメラに写った内傍情報との類似度、残差Ｒｋを輪郭形状の形状残差に置き換えることで、ステップＳ１２０６と同様に算出することができる。

（変形例３）
上記の圧迫変形の推定処理（Ｓ６０９）では、ＭＲＩ画像の圧迫変形の妥当性を評価する評価値として、ＰＡＴ画像との類似度と、赤外線カメラに写った被検体の輪郭（シルエット）形状に対する形状誤差を用いた。しかし、別の方法で評価値を求めてもよい。

例えば、変形メッシュに基づき変形ＭＲＩ画像を生成し、変形ＭＲＩ画像を赤外線カメラ視点から投影した変形後の体表近傍ＭＩＰ像を生成し、体表近傍ＭＩＰ像と赤外線カメラ画像の類似度を算出して、この類似度を評価値に加えてもよい。

つまり、ステップＳ１２０５において、メッシュＭから変形メッシュＤＭｋへの各ノードの変位を表す変形関数Ｆｋを生成する。そして、ＭＲＩ画像に対して、ＭＲＩ画像座標系Ｃ_MRIから前面カメラ装置座標系C_CAM1への座標変換行列Ｔ_MtoC1による剛体変換を施し、変形関数Ｆｋによる変形を施して、変形ＭＲＩ画像I_{D_MRIonC1}kを生成する。

次に、ステップＳ９０５、Ｓ９０６と同様の方法により、変形ＭＲＩ画像Ｉ_{D_MRIonC1}kから表面部分領域の近傍情報を利用した体表近傍変形ＭＩＰ像I_{D_MIPonC1}kを生成する。そして、体表近傍ＭＩＰ像Ｉ_{D_MIPonC1}kと保持状態の赤外線カメラ画像Ｉ'_CAM1との間の類似度Ｓ_MIPkを算出する。最後に、評価値ＥｋをステップＳ１２０６における類似度Ｓ_MRIk、残差Ｒｋ、および、算出した類似度Ｓ_MIPkに基づく重み付け和を下式により算出する。

Ｅｋ = ａS_MRIk+ｂ{1/(1+Rk)}+ｃS_MIPk …(3)
ここで、ａ、ｂ、ｃは重み係数（ａ＋ｂ＋ｃ＝１）。

これにより、体表近傍の表在血管の情報をさらに利用して変形パラメータを求めることができ、高精度な位置合わせが可能になる。

（変形例４）
上記の説明では、表在血管の可視化に、近赤外線による赤外線カメラ画像を利用する例を説明したが、例えば、体内における内部反射によって得られる偏光成分により表在血管を可視化した画像を利用することもできる。例えば、ハロゲンライトなどの光を体表に照射し、体表で反射する表面反射成分と、一旦、体内に侵入して体内で吸収と散乱を経た後に再び体表に出てくる内部反射成分とを区別する。こうすれば、体内の情報を可視化した画像が得られる。当該画像においては、赤外線カメラ画像と同様、内部反射成分のうち、血管内部のヘモグロビンの吸収箇所が周囲よりも暗く描出される。

このように、第一の撮像装置によって対象物を撮像した三次元画像と、第二の撮像装置によって対象物の表層部を撮像した二次元画像とが取得され、三次元画像から対象物の表面位置を示す情報が取得される。そして、表面位置を示す情報に基づいて、三次元画像を第二の撮像装置の視点から見た投影像が生成され、投影像および二次元画像を用いて、対象物に関して三次元画像と二次元画像との位置合わせが行われる。

つまり、本実施形態の画像処理装置は、第一の撮像装置（例えば、ＭＲＩやＣＴ）により対象物を撮像した第一の画像と、第一の撮像装置とは異なる第二の撮像装置１２０により対象物を撮像した第二の画像と、第二の画像撮像装置に対する位置関係が対応付けられている撮像部（３０５、３０６、３０７）により対象物を撮像した第三の画像と、を取得する画像取得部（１０１、１０３、１０４）と、第一の画像と第三の画像との比較結果を用いて、第一の画像と前記第二の画像との位置合わせをする位置合わせ部１１３と、を備える。これにより、対象物である被検体のＭＲＩ画像やＣＴ画像などの三次元画像を、ＰＡＴ画像に対して高精度に位置合わせすることが可能になる。

（第２の実施形態）
以下、本発明にかかる第２の実施形態の画像処理を説明する。なお、第２の実施形態において、第１の実施形態と同様の構成については、同一符号を付して、その詳細説明を省略する。

第１の実施形態では、非保持状態の被検体の画像とＭＲＩ画像との間で剛体変換パラメータを求め（非保持状態における位置合わせ）、その後、保持状態の被検体の画像とＭＲＩ画像との間で変形パラメータを推定（圧迫変形の推定）する方法を説明した。この方法は、非保持状態の被検体とＭＲＩ画像の被検体の形状とが異なる場合に生じる誤差や、非保持状態から保持状態へ移行する際の被検体の微動によって生じる誤差を含む可能性がある。

第２の実施形態では、非保持状態の被検体に関する計測情報を用いずに、保持状態の被検体（第二の状態の被検体）の画像とＭＲＩ画像（第一の状態の被検体）との間で、剛体変換パラメータと変形パラメータを同時に推定する方法を説明する。

図１６は、第２の実施形態の画像処理装置１０を含むモダリティシステムの構成例を示す図である。図１に示す第１の実施形態の構成と異なるのは、第１の実施形態の仮想投影像評価部１０８と変形画像評価部１１１との代わりに、評価部１１４を備える点である。また、第２の実施形態においては、位置合わせ部１１３内における情報の流れが異なるが、その点は、各部の動作と処理の説明において詳述する。

図１７のフローチャートにより第２の実施形態の画像処理装置１０の各部の動作と処理を説明する。なお、ステップＳ６０１からＳ６０６、Ｓ６１０における動作と処理は第１の実施形態と同様であり、その詳細説明を省略する。

ステップＳ１７０１において、赤外線カメラ画像、ＰＡＴ画像、ＭＲＩ画像の血管情報に基づき、ＭＲＩ画像における被検体の位置・姿勢と圧迫変形が推定される。この処理は、剛体変換部１０６、仮想投影像生成部１０７、変形推定部１０９、変形画像生成部１１０、評価部１１４によって行われる。

第２の実施形態では、赤外線カメラに対するＭＲＩ画像上の被検体の位置・姿勢と、被検体の圧迫変形を表す変形パラメータとを位置合わせパラメータとする。変形画像生成部１１０は、仮定した位置合わせパラメータの候補値に基づきＭＲＩ画像を赤外線カメラ座標系（前面カメラ座標系）に座標変換した後に圧迫変形を施すことで、当該位置合わせパラメータの候補値に基づく変形ＭＲＩ画像を生成する。

次に、仮想投影像生成部１０７は、生成された変形ＭＲＩ画像に対して、赤外線カメラの視点を基準に透視投影を行ったＭＩＰ像を生成する。その際、仮想投影像生成部１０７は、赤外線カメラの視点からの各投影線上において、乳房の三次元表面の近傍領域の輝度情報のみを可視化、あるいは、強調して可視化する画像処理を施して、乳房の表在血管を強調したＭＩＰ像を生成する。

次に、剛体変換部１０６は、変形ＭＲＩ画像を、ＰＡＴ１２に対して予め校正された赤外線カメラの位置・姿勢に基づき、赤外線カメラ座標系（前面カメラ座標系）からＰＡＴ画像座標系に座標変換する。

次に、評価部１１４は、血管情報が可視化された赤外線カメラ画像とＭＩＰ像の間、および、ＰＡＴ画像と座標変換後の変形ＭＲＩ画像との間の画像の類似度を算出し、それら類似度を総合した評価値を算出する。そして、仮定した位置合わせパラメータを変動させて、評価値が最大になる位置合わせパラメータを選択する。すなわち、当該位置合わせパラメータによって、ＭＲＩ画像とＰＡＴ画像の間で圧迫変形を含む位置合わせを行う。

図１８、図１９のフローチャートにより位置・姿勢と圧迫変形の推定(Ｓ１７０１)の詳細を説明する。

Ｓ１８０１で、変形推定部１０９は、ステップＳ６０２で取得された被検体の表面形状に基づき、被検体の形状を表すメッシュＭを生成する。この処理は、第１の実施形態のステップＳ１２０１と略同様であり、詳細説明を省略する。

Ｓ１８０２で、剛体変換部１０６は、ＭＲＩ画像を赤外線カメラ座標系に平行移動する。ステップＳ１８０２の処理は、ステップＳ６０３で取得されたＭＲＩ画像における乳頭位置と、ステップＳ６０６で取得された保持状態における赤外線カメラ画像上の乳頭位置に基づき行われる。この処理は、第１の実施形態のステップＳ９０１の処理と略同様であり、詳細説明を省略する。ただし、ステップＳ９０１の処理は非保持状態における赤外線カメラ画像に基づき行われるが、ステップＳ１８０２の処理は保持状態における赤外線カメラ画像に基づき行われる点で異なる。

次に、Ｓ１８０３で、剛体変換部１０６と変形推定部１０９は、剛体変換パラメータが取り得る値と、変形パラメータの各成分が取り得る値を組み合わせて複数（Ｎｔ組）の位置合わせパラメータ（変換パラメータ）の仮説t_i（1≦i≦Nt）を設定する。

例えば、剛体変換部１０６は、剛体変換パラメータが取り得る値として、第１の実施形態のステップＳ９０２と同様に、複数（Ｎ_θ組）の回転パラメータθj（1≦j≦N_θ）を設定する。変形推定部１０９は、変形パラメータが取り得る値として、第１の実施形態のステップＳ１２０２と同様に、複数（Ｎｐ組）の変形パラメータp_k（1≦k≦Np）を設定する。そして、回転パラメータθjと変形パラメータp_kを組み合わせて複数（Ｎｔ=Ｎ_θ×Np組）の変換パラメータｔ_i（１≦i≦Ｎｔ）を設定する。なお、ＰＡＴ画像座標系と前面カメラ座標系との関係は既知であるので、これは、ＭＲＩ画像（第一の状態の被検体）からＰＡＴ画像（第二の状態の被検体）への位置合わせパラメータの候補値を設定することと等価である。また、変換パラメータｔ_iは、剛体変換部１０６と変形推定部１０９との間で共有されるものとする。

Ｓ１８０４で、評価部１１４は初期化処理を行う。つまり、評価部１１４は、ループ変数iに1を設定し、後述する評価値の最大値Ｅ_MAXに０を設定し、後述する変換パラメータｔ_MAXにt₁を設定する。

Ｓ１８０５で、剛体変換部１０６は、平行移動後のＭＲＩ画像を、乳頭位置を基準に変換パラメータｔ_i（つまりθj）だけ回転移動したＭＲＩ画像Ｉ_MRIonC1iを生成する。そして、剛体変換部１０６は、ＭＲＩ画像I_MRIonC1iと座標変換行列Ｔ_MtoC1iを変形推定部１０９と変形画像生成部１１０に出力する。この処理は、第１の実施形態のステップＳ９０４と略同様であり、詳細説明を省略する。

Ｓ１８０６で、変形推定部１０９は、変換パラメータｔ_i（つまりθjとp_k）を用いて、有限要素法に基づく物理変形シミュレーションをメッシュＭに施した変形メッシュＤＭｉを生成し、変形関数Ｆｊｋを変形画像生成部110に出力する。つまり、変形推定部１０９は、ステップＳ１８０５で導出された座標変換行列Ｔ_MtoC1iを用いて、回転パラメータθjに対応する剛体変換をメッシュｍに施したメッシュＭｉを生成する。そして、変形推定部１０９は、メッシュＭｉに物理変形シミュレーションを施すことで、変形後のメッシュである変形メッシュＤＭｉを生成する。このときの変形関数Ｆｉ(x, y, z)を、メッシュＭから変形メッシュＤＭｉへの各ノードの変位を与える変位ベクトルｄｉ_L（１≦Ｌ≦Ｎｍ）と定義する。なお、この処理は、第１の実施形態のステップＳ１２０４と同様である。

Ｓ１８０７で、変形画像生成部１１０は、変換パラメータｔ_i（つまりθjとp_k）に対応する変換をＭＲＩ画像に施した変形ＭＲＩ画像を生成し、変形ＭＲＩ画像を仮想投影像生成部１０７に出力する。つまり、変形画像生成部１１０は、ＭＲＩ画像に座標変換行列Ｔ_MtoC1iを用いる剛体変換を施して、ＭＲＩ画像を回転パラメータθjに対応する前面カメラ座標系Ｃ_CAM1に座標変換する。変形画像生成部１１０は、座標変換後のＭＲＩ画像に変形関数Ｆｊｋを用いる変形処理を施して、変形ＭＲＩ画像Ｉ_{D_MRIonC1}iを生成する。

Ｓ１８０８で、仮想投影像生成部１０７は、前面カメラ座標系Ｃ_CAM1において、変形ＭＲＩ画像Ｉ_{D_MRIonC1}iを前面赤外線カメラ３０５の視点から観測した場合に視野に入るだろう被検体の表面形状（部分表面領域）を求める。この処理は、第１の実施形態のステップＳ９０５における剛体変換されたＭＲＩ画像を、変形ＭＲＩ画像Ｉ_{D_MRIonC1}iに置き換えたものであり、詳細説明を省略する。

Ｓ１８０９で、仮想投影像生成部１０７は、変形ＭＲＩ画像Ｉ_{D_MRIonC1}iにおける表面部分領域の近傍情報を利用した体表近傍変形ＭＩＰ像I_{D_MIPonC1}iを生成し、体表近傍変形ＭＩＰ像を評価部１１４に出力する。この処理は、第１の実施形態のステップＳ９０６における剛体変換されたＭＲＩ画像を、変形ＭＲＩ画像Ｉ_{D_MRIonC1}iに置き換えたものであり、詳細説明を省略する。

Ｓ１８１０で、変形画像生成部１１０は、ステップＳ１８０７で生成した変形ＭＲＩ画像Ｉ_{D_MRIonC1}iを前面カメラ座標系からＰＡＴ画像座標系に座標変換した変形ＭＲＩ画像I_{D_MRIonP}iを生成し、変形ＭＲＩ画像を評価部１１４に出力する。つまり、変形画像生成部１１０は、変形ＭＲＩ画像Ｉ_{D_MRIonC1}iを、座標変換行列Ｔ_C1toDを用いてＰＡＴ装置座標系C_DEVに座標変換する。更に、変形画像生成部１１０は、座標変換行列Ｔ_PtoDの逆行列を用いる座標変換により、ＰＡＴ画像座標系Ｃ_PATにおける変形ＭＲＩ画像I_{D_MRIonP}iを生成する。

Ｓ１８１１で、評価部１１４は、体表近傍変形ＭＩＰ像Ｉ_{D_MIPonC1}iと保持状態の赤外線カメラ画像の類似度Ｓ_MIPi（0≦S_MIPi≦1）、および、変形ＭＲＩ画像I_{D_MRIonP}iとＰＡＴ画像の類似度Ｓ_MRIi（0≦S_MRIi≦1）を算出する。さらに、評価部１１４は、変形ＭＲＩ画像Ｉ_{D_MRIonP}iと保持状態の赤外線カメラ画像における乳房形状の残差Ｒｉを算出する。そして、評価部１１４は、それら類似度と残差を総合した評価値Ｅｉを算出する。

類似度Ｓ_MIPiの算出方法は、第１の実施形態のステップＳ９０７における体表近傍ＭＩＰ像Ｉ_MIPonC1jと非保持状態の赤外線カメラ画像を、体表近傍変形ＭＩＰ像I_{D_MIPonC1}iと保持状態の赤外線カメラ画像に置き換えたものである。また、類似度Ｓ_MRIiの算出方法は、第１の実施形態のステップＳ１２０６における変形ＭＲＩ画像Ｉ_{D_MRI}kを、変形ＭＲＩ画像I_{D_MRIonP}iに置き換えたものである。また、評価部１１４は、赤外線カメラ画像に写った被検体の輪郭（シルエット）の形状と、赤外線カメラ画像に投影した変形メッシュＤＭｊｋの外郭の形状との間の差として残差Ｒｉを算出する。残差Ｒｉの求め方はステップＳ１２０６と同様である。

評価値Ｅｋは、例えば、類似度Ｓ_MIPi、Ｓ_MRIiと残差Ｒｉに基づく重み付け和として、下式で表される。

Ｅｉ=ａS_MIPi+ｂS_MRIi＋ｃ{1/(1+Ri)} …(４)
ここで、ａ、ｂ、ｃは重み係数（ａ+ｂ+ｃ＝１）。

また、式(４)の第三項において（１＋Ｒｋ）の逆数を用いるのはステップＳ１２０６と同様の理由からである。

次に、Ｓ１８１２で、評価部１１４は、評価値Ｅｉと評価値の最大値Ｅ_MAXを比較する。

そして、評価値Ｅｉが評価値の最大値Ｅ_MAXを超える（Ei＞E_MAX）場合（Ｓ１８１２−Ｙｅｓ）、評価部１１４は、処理をステップＳ１８１３に進める。

Ｓ１８１３で、評価部１１４は、Ｅ_MAXを（Ｅ_MAX＝Ｅｉ）と更新する。また、評価部１１４は、Ｅ_MAXに対応する変換パラメータｔ_MAXを（ｔ_MAX＝ｔ_i）と更新する。

一方、ステップＳ１８１２の比較で、評価値Ｅｉが評価値の最大値Ｅ_MAX以下（Ｅｉ≦Ｅ_MAX）の場合（Ｓ１８１２−Ｎｏ）、評価部１１４は、Ｓ１８１３の更新を行わずに処理をステップＳ１８１４に進める。

Ｓ１８１４で、評価部１１４は、ループ変数iをインクリメントする。そして、Ｓ１８１５で、評価部１１４は、ループ変数ｉと仮説の総数Ｎｔとを比較する。ループ変数iが仮説の総数Ｎｔ以下（ｉ≦Ｎｔ）の場合（Ｓ１８１５−Ｎｏ）、評価部１１４は処理をステップＳ１８０５に戻す。一方、ステップＳ１８１５の比較で、ループ変数iが仮説の総数Ntを超える（ｉ＞Ｎｔ）の場合（Ｓ１８１５−Ｙｅｓ）、評価部１１４は処理をステップＳ１８１６に進める。つまり、仮説の総数Ｎｔ分、ステップＳ１８０５からＳ１８１５の処理が繰り返される。

仮説の総数Ｎｔ分の処理が終了すると、Ｓ１８１６で、評価部１１４は、評価値の最大値Ｅ_MAXに対応する変換パラメータｔ_MAX（つまりθ_MAXとｐ_MAX）を変形画像生成部１１０に出力する。言い換えれば、複数の変換パラメータから、評価値の最大値Ｅ_MAXに対応する変換パラメータｔ_MAXが選択される。

Ｓ１８１７で、変形画像生成部１１０は、変換パラメータｔ_MAXに対応する変形ＭＲＩ画像Ｉ_{D_MRIonP}を画像表示部１１２に出力する。

以上で、剛体変換部１０６、仮想投影像生成部１０７、変形推定部１０９、変形画像生成部１１０、評価部１１４による位置・姿勢と圧迫変形の推定(Ｓ１７０１)が終了する。この処理によれば、様々な回転と変形を表す変換パラメータｔ_iを仮定して圧迫変形を含む変換を行った結果の中で変形の適切性の評価値Ｅｉが最大になる変換パラメータｔ_MAXによって、変形ＭＲＩ画像I_{D_MRIonP}が生成される。

このように、ＭＲＩ画像を圧迫変形して生成した表在血管を強調したＭＩＰ像、保持状態の赤外線カメラ画像、ＰＡＴ画像の血管情報を比較して、ＭＲＩ画像上の被検体の位置・姿勢と圧迫変形を推定する。これにより、ＰＡＴ画像の撮像時、理想的な非保持状態の赤外線カメラ画像が撮像されていない場合や、非保持状態から保持状態に移る期間に乳房の位置・姿勢が変化したりする場合も、高精度な変形位置合わせを実現することができる。

（変形例）
上記の処理では、ＰＡＴ画像と赤外線カメラ画像の血管情報を同じ評価値Ｅｉに組み込んで利用する例を説明した。しかし、これら情報を同時に利用する必要はなく、例えば、ＰＡＴ画像の血管情報を先に利用し、その後、赤外線カメラ画像の血管情報を利用して、変形位置合わせを実施してもよい。

例えば、変換パラメータｔ_iを変動させ、変形ＭＲＩ画像とＰＡＴ画像との類似度を評価し、評価値が最大になる変換パラメータｔ_MAXを求める。次に、変換パラメータｔ_MAX付近に限定して変換パラメータt_iを変動させて、体表近傍変形ＭＩＰ像Ｉ_MIPonC1iと赤外線カメラ画像の類似度を評価し、評価値が最大になる変換パラメータt_MAX2を求める。そして、変換パラメータt_MAX2に対応する変形MRI画像を圧迫変形位置合わせ結果とする。これにより、ＭＲＩ画像に対して、形態画像として大まかな血管情報が描出されたＰＡＴ画像を利用して粗い変形位置合わせを実施した後、機能画像としてより詳細な血管情報が描出された赤外線カメラ画像を用いて精密な変形位置合わせを実施することができる。

また、第１の実施形態と同様に、非保持状態における剛体変換パラメータの推定を行った後に、θjの仮説の設定範囲を剛体変換パラメータの近傍に限定して、保持状態における剛体変換パラメータと変形パラメータを推定する構成でもよい。これにより、仮説の総数が減少し、処理の高速化を図ることができる。勿論、第１の実施形態の変形例１と同様に、仮説の総当りによる評価ではなく、一般的な最適化アルゴリズムを用いて変換パラメータを求める構成でもよい。この場合、非保持状態で推定した剛体変換パラメータをθjの初期値に用いることができる。

（第３の実施形態）
以下、本発明にかかる第３の実施形態の画像処理を説明する。なお、第３の実施形態において、第１、第２の実施形態と略同様の構成については、同一符号を付して、その詳細説明を省略する。

第１、第２の実施形態においては、ＰＡＴ１２が搭載する赤外線カメラによる乳房の撮像画像、および、ＭＲＩ画像から生成した乳房の表在血管の強調ＭＩＰ像に基づき、ＭＲＩ画像からＰＡＴ画像への位置合わせを行う例を説明した。しかし、ＭＲＩ画像などの予め撮像された被検体の三次元画像と赤外線カメラに写った被検体の表在血管の情報を用いて、当該三次元画像と位置合わせを行う対象画像はＰＡＴ画像に限られない。

例えば、診断支援を目的として超音波画像の撮像断面に対応する断面（以下、対応断面）の画像を、三次元画像であるＣＴ画像やＭＲＩ画像から生成して（切り出して）提示する場合が挙げられる。

第３の実施形態では、超音波画像、予め位置関係が対応付けられた赤外線カメラによる被検体の撮像画像、ＭＲＩ画像から生成した被検体の表在血管の強調ＭＩＰ像に基づき、ＭＲＩ画像撮像時の被検体から超音波画像撮像時の被検体への位置合わせ例を説明する。なお、以下では、乳房を被検体とし、ＭＲＩ画像を撮像した際の被検体を「第一の状態における被検体」、超音波画像を撮像した際の被検体を「第二の状態における被検体」と呼ぶ場合がある。

また、第３の実施形態において予め撮像されるＭＲＩ画像は、第一の状態における被検体として伏臥状態の乳房を撮像した画像であるが、超音波画像は、第二の状態における被検体として仰臥状態の乳房を撮像したものである。従って、ＭＲＩ画像と超音波画像の間での位置合わせを行うために、両画像間の剛体変換だけでなく、重力変形も推定する必要がある。以下、第３の実施形態の処理の流れを説明する。

超音波撮像装置に搭載された超音波探触子を仰臥状態の被検体に接触させて、超音波画像が撮影される。超音波探触子には、その位置と姿勢を計測するセンサ（磁気式、光学式など）が装着され、超音波撮影が行われている最中の超音波探触子の位置・姿勢が計測される。つまり、センサが基準とする座標系（以下、センサ座標系）における超音波画像の撮影領域が計測される。

また、被検体を撮影するように設置された赤外線カメラは、センサ座標系において、位置・姿勢の校正が行われている。従って、センサ座標系を介して、超音波画像と赤外線カメラ画像の位置・姿勢の対応付けが可能になる。そして、赤外線カメラによって撮像された、仰臥状態の乳房の表在血管が写った赤外線カメラ画像を取得する。

次に、第１、第２の実施形態と同様に、被検体を撮像したＭＲＩ画像から乳房の三次元表面形状を抽出する。続いて、赤外線カメラに対してＭＲＩ画像上の被検体の様々な位置・姿勢と重力変形を表す位置合わせパラメータの候補値を仮定し、位置合わせパラメータの候補値ごとに重力変形位置合わせを施した変形ＭＲＩ画像を生成する。

次に、変形ＭＲＩ画像ごとに画像内の体表近傍の情報を利用して赤外線カメラ側の表在血管を強調したＭＩＰ像を生成する。続いて、各ＭＩＰ像と赤外線カメラ画像の類似性を最大にする、赤外線カメラに対するＭＲＩ画像上の被検体の位置・姿勢と重力変形を推定する。

次に、予め対応付けられた超音波画像に対する赤外線カメラ画像の位置・姿勢に基づき、ＭＲＩ画像を超音波画像の座標系へ変形位置合わせした変形ＭＲＩ画像を生成し、超音波画像との間の変形位置合わせを行う。

このように、ＭＲＩ画像を重力変形させて赤外線カメラ側の表在血管を強調したＭＩＰ像と仰臥状態の赤外線カメラ画像を比較することで、赤外線カメラに対するMRI画像上の被検体の位置・姿勢と重力変形を推定することができる。従って、赤外線カメラとＭＲＩ画像の高精度な変形位置合わせを実現でき、超音波画像と赤外線カメラの位置関係を利用して、超音波画像とＭＲＩ画像の高精度な変形位置合わせを行うことができる。

（変形例１）
第３の実施形態では、伏臥状態で撮像されたＭＲＩ画像と、仰臥状態で撮像された超音波画像の間の乳房の変形状態を補正するために、赤外線カメラに対するMRI画像上の被検体の位置・姿勢に加え、重力変形も推定した。

しかし、例えば、超音波画像と位置合わせする対象であるＣＴ画像やＭＲＩ画像などの三次元画像が、超音波画像と同じ仰臥状態で撮像されていれば、両画像間において重力変形による補正を施す必要はない。この場合の処理を、第３の実施形態の処理との差分（赤外線カメラにおける第二の状態の被検体に対するＭＲＩ画像上の第一の状態の被検体の位置合わせ）についてのみ説明する。

まず、赤外線カメラに対するＭＲＩ画像上の被検体に対して位置・姿勢のみ様々に変動させた位置合わせパラメータの候補値を仮定し、仮定した位置・姿勢ごとにMRI画像内の体表近傍の情報を利用することで表在血管を強調したＭＩＰ像を生成する。

次に、生成した各ＭＩＰ像と赤外線カメラ画像の類似性を評価して、評価値が最大になる赤外線カメラに対するＭＲＩ画像上の被検体の位置・姿勢を推定する。これにより、被検体が剛体と見做せる場合は、超音波画像と位置合わせする対象である三次元画像の間で、高精度に剛体位置合わせを行うことができる。

（変形例２）
第３の実施形態では、赤外線カメラ画像を利用することで、ＭＲＩ画像などの被検体の三次元画像と、位置・姿勢センサで計測された超音波探触子によって取得された二次元の超音波画像との位置合わせを行う例を説明した。

しかし、被検体の三次元画像と位置合わせを行う対象は、例えば、超音波画像と同様に、位置・姿勢センサで計測された超音波探触子で取得された二次元のPAT画像でもよい。二次元のＰＡＴ画像は、超音波探触子に近赤外光の光源を備えることにより、人体に対して超音波ではなく近赤外光を照射するという点を除いて、超音波画像と同様の構成で画像を取得する。従って、第３の実施形態と同様に、赤外線カメラに対するＭＲＩ画像上の被検体の位置・姿勢および重力変形を推定する。そして、位置・姿勢センサにより分かっている二次元のＰＡＴ画像と赤外線カメラの位置関係を利用して、二次元のＰＡＴ画像とＭＲＩ画像の高精度な変形位置合わせを行うことができる。

（第４の実施形態）
以下、本発明にかかる第４の実施形態の画像処理を説明する。なお、第４の実施形態において、第１乃至第３の実施形態と略同様の構成については、同一符号を付して、その詳細説明を省略する。

第１乃至第３の実施形態においては、第一の処理として、赤外線カメラに対する、予め撮像されたＭＲＩ画像などの被検体の三次元画像の位置・姿勢および変形状態を推定した。そして、第二の処理として、目的とするモダリティと赤外線カメラとの位置関係を利用して、三次元画像とモダリティの位置合わせを行った。

しかし、実施する処理は第一の処理のみでもよい。つまり、赤外線カメラ画像と予め撮像された第二の状態の被検体の三次元画像の両画像に写った表在血管の情報を利用して、赤外線カメラに写った被検体に対する、予め撮像された第一の状態の被検体の三次元画像の位置・姿勢と変形を推定する処理のみでもよい。

上記の位置・姿勢および変形の推定により、例えば、被検体の手術時において、被検体を写す赤外線カメラに対して、予め撮像した被検体のＭＲＩ画像やＣＴ画像などの三次元画像を位置合わせすることができる。これにより、三次元画像に写った病変部位などを、手術中の被検体と対応付けて参照することが可能になる。

このように、赤外線カメラで撮像された二次元画像と、ＭＲＩ画像やＣＴ画像などの三次元画像から取得される被検体の表面近傍の内部構造の情報を利用して、三次元画像を被検体に対して高精度かつ自動的に位置合わせする仕組みを提供することができる。

（第５の実施形態）
以下、本発明にかかる第５の実施形態の画像処理を説明する。なお、第５の実施形態において、第１乃至第３の実施形態と略同様の構成については、同一符号を付して、その詳細説明を省略する。

第１乃至第３の実施形態においては、目的とするモダリティに位置が対応付けられた赤外線カメラ（一般的なカメラでも良い）に対する、予め撮像されたＭＲＩ画像などの被検体の三次元画像の位置・姿勢や変形状態を推定することで、三次元画像と目的とするモダリティの位置合わせを行った。

しかし、目的とするモダリティ（本実施形態ではＰＡＴ画像とする）と位置が対応付けられた撮像部は、被検体の距離画像を取得する距離画像カメラであっても良い。これにより、第１乃至第３の実施形態ではカメラ画像から取得できるのが被検体の２次元形状であったのに対し、被検体の距離画像を取得することによって、被検体の３次元形状を取得することができる。そして、取得した被検体の３次元形状と予め撮像された三次元画像（本実施形態ではＭＲＩ画像とする）に写る被検体の３次元形状を比較することで、距離画像カメラに対するＭＲＩ画像の位置・姿勢や変形状態を推定することができる。最後に、距離画像カメラとＰＡＴ装置の位置の対応関係を利用して、ＭＲＩ画像とＰＡＴ画像の位置合わせを行うことができる。以下、第５の実施形態の機能構成及び処理の流れを説明する。

ここで、第５の実施形態における機能構成及び処理の流れは、第１の実施形態における図１の機能構成及び図６の処理の流れの赤外線カメラに関する部分を、距離画像カメラに置き換えたものとなる。従って、共通する機能構成及び処理の流れは説明を省略し、差分のみを説明する。まず、機能構成について説明する。

距離画像カメラは、該カメラが撮像する画像に撮像されるシーン中の各点までの距離を画素毎に（各画素に写っているシーン中の部分毎に）計測し、該距離を表す情報を各画素の画素値として有する距離画像を生成するカメラである。距離画像カメラによって対象物体を計測することで、対象物体の表面の三次元形状を距離画像として撮像できる。なお、本実施形態で用いる距離画像カメラは、タイムオブフライト方式やパターン投光方式など、公知のいずれの方式に基づくものであっても良い。

本実施形態では、カメラ画像取得部１０４を、距離画像を取得する機能を備えた構成要素（距離画像取得部と定義する）に置き換える。

本実施形態では、二次元形状取得部１０５を、距離画像から被検体の三次元形状を取得する機能を備えた構成要素（形状取得部と定義する）に置き換える。

本実施形態では、仮想投影像生成部１０７と仮想投影像評価部１０８を、ＭＲＩ画像から取得した被検体の三次元形状と、非保持状態の被検体を撮像した距離画像から取得した被検体の三次元形状とを比較する機能を備えた構成要素（第一の評価部と定義する）に置き換える。

本実施形態では、変形画像評価部１１１を、変形ＭＲＩ画像から取得した被検体の三次元形状と、保持状態の被検体を撮像した距離画像から取得した被検体の三次元形状とを比較する機能（第二の評価部）に置き換える。

次に、処理の流れについて説明する。

本実施形態では、ステップＳ６０５を、距離画像取得部が非保持状態及び保持状態の被検体を撮像した距離画像を取得する処理に置き換える。

本実施形態では、ステップＳ６０６を、形状取得部が、ステップＳ６０５で取得した非保持状態の乳房の距離画像から被検体の乳房の三次元形状と乳頭位置を取得する処理に置き換える。被検体の三次元形状の取得については、例えば、被検体の写った距離画像と被検体の写っていない距離画像の差分画像を生成して、差分画像中の値が所定値より小さい領域は、被検体以外の領域（背景領域）であると判定して削除することで、被検体の領域における距離データ、つまり被検体の三次元形状を取得することができる。なお、乳頭位置は、ステップＳ６０２と同様、三次元形状の曲率を計算し、特徴点を検出することで取得することができる。

本実施形態では、ステップＳ６０７を、剛体変換部１０６及び第一の評価部が、ステップＳ６０２及びＳ６０３で取得したＭＲＩ画像に写る乳房の三次元形状及び乳頭位置と、ステップＳ６０６で取得した距離画像に写る非保持状態の乳房の三次元形状及び乳頭位置とを比較し、両者の三次元形状が最も整合するような剛体変換を推定する処理に置き換える。より具体的には、剛体変換部１０６は、ステップＳ６０３で取得したＭＲＩ画像上の三次元の乳頭位置がステップＳ６０６で取得した距離画像上の三次元の乳頭位置に一致するように、平行移動を表す座標変換行列を算出する。そして、第一の評価部は、ＭＲＩ画像上の三次元形状に対して平行移動を表す変換行列を適用することで、乳頭位置を一致させた上で、ＭＲＩ画像と距離画像の三次元形状の残差を算出する。そして、剛体変換部１０６は、残差が最小になるような回転移動の座標変換行列を探索して算出し、平行移動と回転移動を合わせた剛体変換の座標変換行列を算出する。

本実施形態では、ステップＳ６０９を、変形推定部１０９、変形画像生成部１１０、及び第二の評価部が、ＭＲＩ画像を変形させた変形ＭＲＩ画像から取得した乳房の三次元形状と、ステップＳ６０６で取得した距離画像に写る保持状態の乳房の三次元形状とを比較し、両者の三次元形状が最も整合するような変形を推定する処理に置き換える。より具体的には、変形推定部１０９は、第１の実施形態におけるステップＳ６０９と同様の方法で、ＭＲＩ画像を変形させた変形ＭＲＩ画像を生成する。そして、第二の評価部は、変形ＭＲＩ画像からステップＳ６０２と同様の方法で取得した乳房の三次元形状と、距離画像の保持状態の乳房の三次元形状と残差を算出する。そして、変形推定部１０９は、残差が最小になるような変形パラメータを探索して算出する。なお、ここでは、三次元形状の残差に基づいて変形パラメータを推定したが、無論、第１の実施形態と同様、残差に加えて、変形ＭＲＩ画像とＰＡＴ画像との間の類似度を加味した評価値に基づいて変形パラメータを推定しても良い。

このように、距離画像から取得した乳房の３次元形状を利用して、距離画像に対するＭＲＩ画像の位置・姿勢や変形状態を推定することで、カメラ画像から得られる２次元形状を利用するよりも、高精度な推定を行うことができる。

（その他の実施形態）
本発明は、以下の処理を実行することによっても実現される。即ち、上述した実施形態の機能を実現するソフトウェア（プログラム）を、ネットワーク又は各種記憶媒体を介してシステム或いは装置に供給し、そのシステム或いは装置のコンピュータ（またはＣＰＵやＭＰＵ等）がプログラムを読み出して実行する処理である。

Claims

第一の撮像手段により対象物を撮像した第一の画像と、第二の撮像手段により前記対象物を撮像した第二の画像と、前記第二の撮像手段に対して位置が対応付けられた第三の撮像手段により前記対象物を撮像した第三の画像と、を取得する画像取得手段と、
前記第一の画像と前記第三の画像との位置合わせを行い、前記第一の画像と前記第三の画像との位置合わせの結果に基づいて、前記第一の画像と前記第二の画像との位置合わせを行う位置合わせ手段と、
を備えることを特徴とする画像処理装置。
前記第三の撮像手段は、前記対象物の表面近傍の血管領域を可視化する赤外線カメラ画像を撮像する赤外線カメラであり、
前記位置合わせ手段は、前記第一の画像と前記第三の画像としての前記赤外線カメラ画像との間で、前記対象物の表面近傍の血管領域を比較して、前記第一の画像と前記第三の画像との位置合わせを行うことを特徴とする請求項１に記載の画像処理装置。
前記第三の撮像手段は、前記対象物の表面の三次元形状を距離画像として撮像する距離画像カメラであり、
前記位置合わせ手段は、前記第一の画像と前記第三の画像としての前記距離画像との間で、前記対象物の三次元の表面形状を比較して、前記第一の画像と前記第三の画像との位置合わせを行うことを特徴とする請求項１に記載の画像処理装置。
前記第一の撮像手段は、前記対象物の三次元画像を前記第一の画像として撮像する撮像手段であり、
前記第一の画像としての前記三次元画像から前記対象物の表面位置を示す情報を取得する情報取得手段と、
前記表面位置を示す情報に基づいて、前記第三の撮像手段の視点から見た前記三次元画像の投影画像を生成する投影画像生成手段と、
を更に備え、
前記位置合わせ手段は、前記投影画像と前記第三の画像との位置合わせの結果に基づいて、前記第一の画像と前記第二の画像との位置合わせを行うことを特徴とする請求項１に記載の画像処理装置。
前記第一の画像は、第一の状態の前記対象物を前記第一の撮像手段により撮像することにより生成された画像であり、
前記第二の画像は、前記第一の状態とは異なる第二の状態の前記対象物を前記第二の撮像手段により撮像することにより生成された画像であり、
前記第三の画像は、前記第二の状態の前記対象物を前記第三の撮像手段により撮像することにより生成された画像であり、
前記第一の画像と前記第三の画像との位置合わせの結果を用いて、前記第一の状態で撮像された前記第一の画像を前記第二の状態の第一の画像に変形する変形画像生成手段を更に備え、
前記変形画像生成手段は、前記第一の画像と前記第二の画像との位置合わせを行うための位置合わせパラメータを用いて、前記第二の状態の第一の画像を生成し、
前記投影画像生成手段は、前記第二の状態の第一の画像を用いて、前記第三の撮像手段の視点から見た投影画像を生成することを特徴とする請求項４に記載の画像処理装置。
前記位置合わせパラメータには、前記対象物の位置、姿勢、および前記対象物の変形を表すパラメータが含まれることを特徴とする請求項５に記載の画像処理装置。
前記変形画像生成手段は、前記位置合わせパラメータの複数の候補値に基づき、前記第二の状態の複数の第一の画像を生成することを特徴とする請求項５に記載の画像処理装置。
前記第二の状態の複数の第一の画像と、前記第二の画像との類似度を取得する評価手段を更に備え、
前記位置合わせ手段は、前記類似度を用いて、前記第一の画像と前記第二の画像との位置合わせを行うことを特徴とする請求項７に記載の画像処理装置。
前記評価手段は、前記類似度を用いて、前記複数の候補値から位置合わせパラメータを一つ選択することを特徴とする請求項８に記載の画像処理装置。
位置合わせされた前記第一の画像と前記第二の画像とを表示部に並べて表示させる表示制御手段を更に備えることをと特徴とする請求項１乃至９のいずれか１項に記載の画像処理装置。
第一の撮像手段により対象物を撮像した第一の画像と、第二の撮像手段により前記対象物を撮像した第二の画像と、前記第二の撮像手段に対して位置が対応付けられた第三の撮像手段により前記対象物を撮像した第三の画像と、を取得する画像取得手段と、
前記第一の画像と前記第三の画像との間の第一の類似度を算出し、前記第一の画像と前記第二の画像との間の第二の類似度を算出する評価手段と、
前記第一の類似度と前記第二の類似度とに基づいて、前記第一の画像と前記第二の画像とを位置合わせする位置合わせ手段と、
を備えることを特徴とする画像処理装置。
前記画像取得手段は、前記第一の画像として三次元画像を取得し、
前記第三の撮像手段の視点から見た前記第一の画像の投影画像を生成する投影画像生成手段を有し、
前記評価手段は、前記第一の画像の前記投影画像と前記第三の画像との間の前記第一の類似度を算出することを特徴とする請求項１１に記載の画像処理装置。
前記第一の撮像手段は、ＭＲＩ装置またはＣＴ装置であることを特徴とする請求項１乃至１２のいずれか１項に記載の画像処理装置。
前記第二の撮像手段は、前記対象物に光を照射し、当該照射した光により発生した光音響信号を検出し、当該光音響信号に基づいて前記第二の画像として光音響画像を生成する光音響撮像装置であることを特徴とする請求項１乃至１３のいずれか１項に記載の画像処理装置。
前記第三の撮像手段は、赤外線を前記対象物に照射することにより前記第三の画像として赤外線カメラ画像を生成する赤外線カメラであることを特徴とする請求項１乃至１４のいずれか１項に記載の画像処理装置。
第一の撮像手段により対象物を撮像した第一の画像を取得し、第二の撮像手段により前記対象物を撮像した第二の画像を取得し、前記第二の撮像手段に対して位置が対応付けられた第三の撮像手段により前記対象物を撮像した第三の画像を取得し、
前記第一の画像と前記第三の画像との位置合わせを行い、
前記第一の画像と前記第三の画像との位置合わせの結果に基づいて、前記第一の画像と前記第二の画像との位置合わせを行うことを特徴とする画像処理方法。
第一の撮像手段により対象物を撮像した第一の画像を取得し、
第二の撮像手段により前記対象物を撮像した第二の画像を取得し、
前記第二の撮像手段に対して位置が対応付けられた第三の撮像手段により前記対象物を撮像した第三の画像を取得し、
前記第一の画像と前記第三の画像との間の第一の類似度を算出し、
前記第一の画像と前記第二の画像との間の第二の類似度を算出し、
前記第一の類似度と前記第二の類似度とに基づいて、前記第一の画像と前記第二の画像とを位置合わせすることを特徴とする画像処理方法。
前記第一の画像として三次元画像を取得し、
前記第三の撮像手段の視点から見た前記第一の画像の投影画像を生成し、
前記第一の画像の前記投影画像と前記第三の画像との間の前記第一の類似度を算出することを特徴とする請求項１７に記載の画像処理方法。
前記第一の撮像手段は、ＭＲＩ装置またはＣＴ装置であることを特徴とする請求項１６乃至１８のいずれか１項に記載の画像処理方法。
前記第二の撮像手段は、前記対象物に光を照射し、当該照射した光により発生した光音響信号を検出し、当該光音響信号に基づいて前記第二の画像として光音響画像を生成する光音響撮像装置であることを特徴とする請求項１６乃至１９のいずれか１項に記載の画像処理方法。
前記第三の撮像手段は、赤外線を前記対象物に照射することにより前記第三の画像として赤外線カメラ画像を生成する赤外線カメラであることを特徴とする請求項１６乃至２０のいずれか１項に記載の画像処理方法。
コンピュータを、請求項１６乃至２１のいずれか１項に記載の画像処理方法をコンピュータに実行させるためのプログラム。
請求項２２に記載されたプログラムが記憶されたコンピュータにより読み取り可能な記憶媒体。