JP2015073855A - 画像処理装置およびその方法 - Google Patents

画像処理装置およびその方法 Download PDF

Info

Publication number
JP2015073855A
JP2015073855A JP2013214140A JP2013214140A JP2015073855A JP 2015073855 A JP2015073855 A JP 2015073855A JP 2013214140 A JP2013214140 A JP 2013214140A JP 2013214140 A JP2013214140 A JP 2013214140A JP 2015073855 A JP2015073855 A JP 2015073855A
Authority
JP
Japan
Prior art keywords
image
dimensional image
mri
deformation
subject
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
JP2013214140A
Other languages
English (en)
Other versions
JP6182045B2 (ja
Inventor
和大 宮狭
Kazuhiro Miyasa
和大 宮狭
卓郎 宮里
Takuro Miyazato
卓郎 宮里
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Canon Inc
Original Assignee
Canon Inc
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Canon Inc filed Critical Canon Inc
Priority to JP2013214140A priority Critical patent/JP6182045B2/ja
Priority to US14/501,729 priority patent/US9311709B2/en
Publication of JP2015073855A publication Critical patent/JP2015073855A/ja
Application granted granted Critical
Publication of JP6182045B2 publication Critical patent/JP6182045B2/ja
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G06T7/344—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods involving models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/14—Transformations for image registration, e.g. adjusting or mapping for alignment of images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10048—Infrared image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10088—Magnetic resonance imaging [MRI]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30068—Mammography; Breast
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30108—Industrial image inspection
    • G06T2207/30148—Semiconductor; IC; Wafer

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)
  • Image Processing (AREA)
  • Processing Or Creating Images (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)
  • Image Analysis (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Quality & Reliability (AREA)

Abstract

【課題】 対象物の三次元画像を、当該対象物自体、あるいは、当該対象物を撮像した他の画像と高精度に位置合わせする。【解決手段】 医用画像取得部101は、第一の撮像装置によって対象物を撮像した三次元画像を取得する。カメラ画像取得部104は、第二の撮像装置によって対象物の表層部を撮像した二次元画像を取得する。三次元形状取得部102は、三次元画像から対象物の表面位置を示す情報を取得する。仮想投影像生成部107は、表面位置を示す情報に基づいて、三次元画像を第二の撮像装置の視点から見た投影像を生成する。位置合わせ部113は、投影像および二次元画像を用いて、対象物に関して三次元画像と二次元画像との位置合わせを行う。【選択図】 図1

Description

本発明は、種々の医用画像収集装置（モダリティ）で撮像した医用画像の画像処理に関する。
光音響断層撮像装置(photoacoustc tomograph: PAT)は、被検体に光パルスを照射して被検体内の吸収物質を励起し、吸収物質の熱弾性膨張により生じる光音響信号を検出して、被検体の光吸収に関わる性質を画像化する。つまり、PATは、照射光に対する被検体内の光エネルギ堆積量分布（光エネルギ吸収密度分布）を画像化する。また、これに基づき、照射波長に関する被検体の光吸収係数分布を画像化する。さらに、複数の波長に関する光吸収係数分布に基づき、被検体を構成する物質の状態（例えば、ヘモグロビンの酸素飽和度など）を画像化することも可能である。
これら画像が、癌などの悪性腫瘍の内外に生じる新生血管に関する情報を可視化すると期待されている。以下では、これら画像を総称して「光音響断層画像（PAT画像）」と呼ぶ。
PATは、エネルギが小さい近赤外光パルスを照射するため、X線などと比べて、人体の深部の画像化が難しい。特許文献1は、乳房を測定対象としたPATの一形態として、乳房を二枚の平板（以下、保持板）で保持して、乳房の厚さを薄くした状態で撮像を行うことを記載する。そのため、PATと、磁気共鳴撮像装置(MRI)などの他のモダリティの併用診断を行う場合、保持による圧迫変形を考慮した変形位置合わせを行う（一方の画像を、他方の画像に合わせるように、変形する）ことで、医師による効率的な診断が可能になる。
PAT画像とMRI画像の位置合わせ方法として、画像マッチングによる方法が挙げられる。例えば、非特許文献1は、平板圧迫された乳房を撮像したX線マンモグラフィ(MMG)と、乳房のMRI画像の間の位置合わせ技術を記載する。具体的には、MRI画像に対して平板圧迫による物理変形シミュレーションを施した変形MRI画像を生成し、変形MRI画像から疑似MMG画像を生成し、疑似MMG画像と実際に撮像されたMMGの画像マッチングにより位置合わせを行う。
また、非特許文献2は、MRI画像に対して平板圧迫による物理変形シミュレーションを施した結果から得られる変形後の乳房の形状を、MMG画像から抽出した乳房の二次元形状に基づいて評価する技術を開示する。
また、上記と異なる複数モダリティの位置合わせによる診断支援として、超音波画像の撮像断面に対応する断面（以下、対応断面）の画像を、三次元画像データであるCT画像やMRI画像などの参照画像から生成して（切り出して）提示する試みがある。例えば、特許文献2は、参照画像であるCT画像やMRI画像を被検体と位置合わせしておき、被検体に対する超音波探触子の位置と姿勢を計測して、超音波画像と参照画像の位置合わせを行う技術を開示する。
特開2010-088627号公報 特許第03871747号公報
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年
しかし、PATが画像化する特徴とMRIが画像化する特徴は異なるため、MRI画像に写る構造物と、PAT画像上の構造物すべてが一致するわけではない。よって、画像マッチングのみで高精度な位置合わせを行うのは困難である。そのため、両方の座標系の間で対応する複数の点を手入力する手間が必要になる。
本発明は、対象物の三次元画像を、当該対象物自体、あるいは、当該対象物を撮像した他の画像と高精度に位置合わせすることを目的とする。
本発明は、前記の目的を達成する一手段として、以下の構成を備える。
本発明にかかる画像処理は、第一の撮像装置によって対象物を撮像した三次元画像と、第二の撮像装置によって前記対象物の表層部を撮像した二次元画像とを取得し、前記三次元画像から前記対象物の表面位置を示す情報を取得し、前記表面位置を示す情報に基づいて、前記三次元画像を前記第二の撮像装置の視点から見た投影像を生成し、前記投影像および前記二次元画像を用いて、前記対象物に関して前記三次元画像と前記二次元画像との位置合わせを行う。
本発明によれば、対象物の三次元画像を、当該対象物自体、あるいは、当該対象物を撮像した他の画像と高精度に位置合わせすることができる。
実施例1の画像処理装置を含むモダリティシステムの構成例を示すブロック図。 医用画像DBが保持する被検体のMRI画像を説明する図。 PATによる被検者の撮像を説明する図。 PATが撮像したPAT画像の例を示す図。 非保持状態における前面赤外線カメラの撮像画像ICAM1の一例を示す図。 実施例1の画像処理装置の各部の動作と処理を説明するフローチャート。 表面形状の取得処理を説明する図。 変形MRI画像とPAT画像の表示例を示す図。 非保持状態における位置合わせの詳細を説明するフローチャート。 仮想投影像生成部が部分表面領域を求める処理を説明する図。 MRI画像内の被検体の体表近傍情報を利用したMIP像を示す図。 圧迫変形の推定の詳細を説明するフローチャート。 メッシュMの生成方法を示す模式図。 保持板による圧迫変形シミュレーションを説明する図。 変形MRI画像ID_MRIonPを示す模式図。 実施例2の画像処理装置を含むモダリティシステムの構成例を示すブロック図。 実施例2の画像処理装置の各部の動作と処理を説明するフローチャート。 位置・姿勢と圧迫変形の推定の詳細を説明するフローチャート。 位置・姿勢と圧迫変形の推定の詳細を説明するフローチャート。
以下、本発明にかかる実施例の画像処理を図面を参照して詳細に説明する。ただし、発明の範囲は図示する例に限定されない。
実施例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の初期値に用いることができる。
以下、本発明にかかる実施例3の画像処理を説明する。なお、実施例3において、実施例1、2と略同様の構成については、同一符号を付して、その詳細説明を省略する。
実施例1、2においては、PAT12が搭載する赤外線カメラによる乳房の撮像画像、および、MRI画像から生成した乳房の表在血管の強調MIP像に基づき、MRI画像からPAT画像への位置合わせを行う例を説明した。しかし、MRI画像などの予め撮像された被検体の三次元画像と赤外線カメラに写った被検体の表在血管の情報を用いて、当該三次元画像と位置合わせを行う対象画像はPAT画像に限られない。
例えば、診断支援を目的として超音波画像の撮像断面に対応する断面（以下、対応断面）の画像を、三次元画像であるCT画像やMRI画像から生成して（切り出して）提示する場合が挙げられる。
実施例3では、超音波画像、予め位置関係が対応付けられた赤外線カメラによる被検体の撮像画像、MRI画像から生成した被検体の表在血管の強調MIP像に基づき、MRI画像撮像時の被検体から超音波画像撮像時の被検体への位置合わせ例を説明する。なお、以下では、乳房を被検体とし、MRI画像を撮像した際の被検体を「第一の状態における被検体」、超音波画像を撮像した際の被検体を「第二の状態における被検体」と呼ぶ場合がある。
また、実施例3において予め撮像されるMRI画像は、第一の状態における被検体として伏臥状態の乳房を撮像した画像であるが、超音波画像は、第二の状態における被検体として仰臥状態の乳房を撮像したものである。従って、MRI画像と超音波画像の間での位置合わせを行うために、両画像間の剛体変換だけでなく、重力変形も推定する必要がある。以下、実施例3の処理の流れを説明する。
超音波撮像装置に搭載された超音波探触子を仰臥状態の被検体に接触させて、超音波画像が撮影される。超音波探触子には、その位置と姿勢を計測するセンサ（磁気式、光学式など）が装着され、超音波撮影が行われている最中の超音波探触子の位置・姿勢が計測される。つまり、センサが基準とする座標系（以下、センサ座標系）における超音波画像の撮影領域が計測される。
また、被検体を撮影するように設置された赤外線カメラは、センサ座標系において、位置・姿勢の校正が行われている。従って、センサ座標系を介して、超音波画像と赤外線カメラ画像の位置・姿勢の対応付けが可能になる。そして、赤外線カメラによって撮像された、仰臥状態の乳房の表在血管が写った赤外線カメラ画像を取得する。
次に、実施例1、2と同様に、被検体を撮像したMRI画像から乳房の三次元表面形状を抽出する。続いて、赤外線カメラに対してMRI画像上の被検体の様々な位置・姿勢と重力変形を表す位置合わせパラメータの候補値を仮定し、位置合わせパラメータの候補値ごとに重力変形位置合わせを施した変形MRI画像を生成する。
次に、各変形MRI画像ごとに画像内の体表近傍の情報を利用して赤外線カメラ側の表在血管を強調したMIP像を生成する。続いて、各MIP像と赤外線カメラ画像の類似性を最大にする、赤外線カメラに対するMRI画像上の被検体の位置・姿勢と重力変形を推定する。
次に、予め対応付けられた超音波画像に対する赤外線カメラ画像の位置・姿勢に基づき、MRI画像を超音波画像の座標系へ変形位置合わせした変形MRI画像を生成し、超音波画像との間の変形位置合わせを行う。
このように、MRI画像を重力変形させて赤外線カメラ側の表在血管を強調したMIP像と仰臥状態の赤外線カメラ画像を比較することで、赤外線カメラに対するMRI画像上の被検体の位置・姿勢と重力変形を推定することができる。従って、赤外線カメラとMRI画像の高精度な変形位置合わせを実現でき、超音波画像と赤外線カメラの位置関係を利用して、超音波画像とMRI画像の高精度な変形位置合わせを行うことができる。
［変形例1］
実施例3では、伏臥状態で撮像されたMRI画像と、仰臥状態で撮像された超音波画像の間の乳房の変形状態を補正するために、赤外線カメラに対するMRI画像上の被検体の位置・姿勢に加え、重力変形も推定した。
しかし、例えば、超音波画像と位置合わせする対象であるCT画像やMRI画像などの三次元画像が、超音波画像と同じ仰臥状態で撮像されていれば、両画像間において重力変形による補正を施す必要はない。この場合の処理を、実施例3の処理との差分（赤外線カメラにおける第二の状態の被検体に対するMRI画像上の第一の状態の被検体の位置合わせ）についてのみ説明する。
まず、赤外線カメラに対するMRI画像上の被検体に対して位置・姿勢のみ様々に変動させた位置合わせパラメータの候補値を仮定し、仮定した位置・姿勢ごとにMRI画像内の体表近傍の情報を利用することで表在血管を強調したMIP像を生成する。
次に、生成した各MIP像と赤外線カメラ画像の類似性を評価して、評価値が最大になる赤外線カメラに対するMRI画像上の被検体の位置・姿勢を推定する。これにより、被検体が剛体と見做せる場合は、超音波画像と位置合わせする対象である三次元画像の間で、高精度に剛体位置合わせを行うことができる。
［変形例2］
実施例3では、赤外線カメラ画像を利用することで、MRI画像などの被検体の三次元画像と、位置・姿勢センサで計測された超音波探触子によって取得された二次元の超音波画像との位置合わせを行う例を説明した。
しかし、被検体の三次元画像と位置合わせを行う対象は、例えば、超音波画像と同様に、位置・姿勢センサで計測された超音波探触子で取得された二次元のPAT画像でもよい。二次元のPAT画像は、超音波探触子に近赤外光の光源を備えることにより、人体に対して超音波ではなく近赤外光を照射するという点を除いて、超音波画像と同様の構成で画像を取得する。従って、実施例3と同様に、赤外線カメラに対するMRI画像上の被検体の位置・姿勢および重力変形を推定する。そして、位置・姿勢センサにより分かっている二次元のPAT画像と赤外線カメラの位置関係を利用して、二次元のPAT画像とMRI画像の高精度な変形位置合わせを行うことができる。
以下、本発明にかかる実施例4の画像処理を説明する。なお、実施例4において、実施例1-3と略同様の構成については、同一符号を付して、その詳細説明を省略する。
実施例1から3においては、第一の処理として、赤外線カメラに対する、予め撮像されたMRI画像などの被検体の三次元画像の位置・姿勢および変形状態を推定した。そして、第二の処理として、目的とするモダリティと赤外線カメラとの位置関係を利用して、三次元画像とモダリティの位置合わせを行った。
しかし、実施する処理は第一の処理のみでもよい。つまり、赤外線カメラ画像と予め撮像された第二の状態の被検体の三次元画像の両画像に写った表在血管の情報を利用して、赤外線カメラに写った被検体に対する、予め撮像された第一の状態の被検体の三次元画像の位置・姿勢と変形を推定する処理のみでもよい。
上記位置・姿勢および変形の推定により、例えば、被検体の手術時において、被検体を写す赤外線カメラに対して、予め撮像した被検体のMRI画像やCT画像などの三次元画像を位置合わせすることができる。これにより、三次元画像に写った病変部位などを、手術中の被検体と対応付けて参照することが可能になる。
このように、赤外線カメラで撮像された二次元画像と、MRI画像やCT画像などの三次元画像から取得される被検体の表面近傍の内部構造の情報を利用して、三次元画像を被検体に対して高精度かつ自動的に位置合わせする仕組みを提供することができる。従って、PAT画像と、MRI画像やCT画像などの三次元画像の間の高精度な位置合わせが可能になる。
［その他の実施例］
また、本発明は、以下の処理を実行することによっても実現される。即ち、上述した実施形態の機能を実現するソフトウェア（プログラム）を、ネットワーク又は各種記録媒体を介してシステム或いは装置に供給し、そのシステムあるいは装置のコンピュータ（又はCPUやMPU等）がプログラムを読み出して実行する処理である。

Claims (15)

  1. 第一の撮像装置によって対象物を撮像した三次元画像と、第二の撮像装置によって前記対象物の表層部を撮像した二次元画像とを取得する画像取得手段と、
    前記三次元画像から前記対象物の表面位置を示す情報を取得する情報取得手段と、
    前記表面位置を示す情報に基づいて、前記三次元画像を前記第二の撮像装置の視点から見た投影像を生成する生成手段と、
    前記投影像および前記二次元画像を用いて、前記対象物に関して前記三次元画像と前記二次元画像との位置合わせを行う位置合わせ手段とを有する画像処理装置。
  2. 前記生成手段は、前記表面位置を示す情報と、前記三次元画像と前記二次元画像の位置合わせを表す位置合わせパラメータの複数の候補値に基づき、複数の前記投影像を生成する請求項1に記載された画像処理装置。
  3. 前記生成手段は、前記複数の候補値のそれぞれに基づき、前記第二の撮像装置によって観測される前記表面位置の部分領域を求め、前記部分領域の情報を用いて、前記複数の投影像を生成する請求項2に記載された画像処理装置。
  4. 前記生成手段は、前記第二の撮像装置の視点から見た前記部分領域のMIP像を前記投影像として生成する請求項3に記載された画像処理装置。
  5. 前記位置合わせ手段は、前記複数の投影像および前記二次元画像を用いて、前記三次元画像と前記二次元画像の位置合わせパラメータを前記複数の候補値から推定する請求項2から請求項4の何れか一項に記載された画像処理装置。
  6. さらに、前記複数の投影像のそれぞれと前記二次元画像の類似度を評価する第一の評価手段を有し、
    前記位置合わせ手段は、前記類似度に基づき、前記三次元画像と前記二次元画像の位置合わせパラメータとして前記複数の候補値の一つを選択する請求項2から請求項4の何れか一項に記載された画像処理装置。
  7. 前記第一の評価手段は、前記投影像および前記二次元画像において可視化された前記対象物の表在血管の輝度情報に基づき前記類似度を評価する請求項6に記載された画像処理装置。
  8. 前記位置合わせパラメータは、前記対象物の変形を表す変形パラメータを少なくとも含み、
    前記生成手段は、前記対象物に関する変形パラメータの複数の候補値に基づき前記三次元画像に変形を施した変形三次元画像を生成し、前記変形三次元画像を前記第二の撮像装置の視点から見た前記投影像として生成する請求項1から請求項7の何れか一項に記載された画像処理装置。
  9. 前記画像取得手段は、第一の状態の前記対象物を前記第一の撮像装置によって撮像した前記三次元画像、第二の状態の前記対象物を前記第二の撮像装置によって撮像した前記二次元画像、および、前記第二の状態の対象物を前記第二の撮像装置によって撮像した三次元画像を取得し、
    さらに、前記複数の候補値のそれぞれに基づき前記第一の状態の三次元画像に変換を施した三次元画像と、前記第二の状態の三次元画像との類似度を評価する第二の評価手段を有し、
    前記位置合わせ手段は、前記第一および第二の評価手段の評価の結果に基づき、前記第一の状態の三次元画像と前記二次元画像の位置合わせパラメータとして前記複数の候補値の一つを選択する請求項6または請求項7に記載された画像処理装置。
  10. さらに、前記第一の状態の三次元画像と前記二次元画像の位置合わせパラメータを用いて前記第一の状態の三次元画像に変換を施した変形三次元画像を生成する手段を有する請求項9に記載された画像処理装置。
  11. さらに、前記第二の状態の三次元画像と前記変形三次元画像を表示する手段を有する請求項10に記載された画像処理装置。
  12. 前記第二の撮像装置は赤外線カメラを有する請求項1から請求項11の何れか一項に記載された画像処理装置。
  13. 第一の撮像装置によって対象物を撮像した三次元画像と、第二の撮像装置によって前記対象物の表層部を撮像した二次元画像とを取得し、
    前記三次元画像から前記対象物の表面位置を示す情報を取得し、
    前記表面位置を示す情報に基づいて、前記三次元画像を前記第二の撮像装置の視点から見た投影像を生成し、
    前記投影像および前記二次元画像を用いて、前記対象物に関して前記三次元画像と前記二次元画像との位置合わせを行う画像処理方法。
  14. コンピュータを請求項1から請求項12の何れか一項に記載された画像処理装置の各手段として機能させるためのプログラム。
  15. 請求項14に記載されたプログラムが記録されたコンピュータが読み取り可能な記録媒体。
JP2013214140A 2013-10-11 2013-10-11 画像処理装置およびその方法 Expired - Fee Related JP6182045B2 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
JP2013214140A JP6182045B2 (ja) 2013-10-11 2013-10-11 画像処理装置およびその方法
US14/501,729 US9311709B2 (en) 2013-10-11 2014-09-30 Image processing apparatus and image processing method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2013214140A JP6182045B2 (ja) 2013-10-11 2013-10-11 画像処理装置およびその方法

Related Child Applications (1)

Application Number Title Priority Date Filing Date
JP2017142840A Division JP6461257B2 (ja) 2017-07-24 2017-07-24 画像処理装置およびその方法

Publications (2)

Publication Number Publication Date
JP2015073855A true JP2015073855A (ja) 2015-04-20
JP6182045B2 JP6182045B2 (ja) 2017-08-16

Family

ID=52809722

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2013214140A Expired - Fee Related JP6182045B2 (ja) 2013-10-11 2013-10-11 画像処理装置およびその方法

Country Status (2)

Country Link
US (1) US9311709B2 (ja)
JP (1) JP6182045B2 (ja)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2017082078A1 (ja) * 2015-11-11 2017-05-18 ソニー株式会社 画像処理装置および画像処理方法
JP2017127623A (ja) * 2016-01-15 2017-07-27 キヤノン株式会社 画像処理装置、画像処理方法、およびプログラム
CN112529960A (zh) * 2020-12-17 2021-03-19 珠海格力智能装备有限公司 目标对象的定位方法、装置、处理器和电子装置

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6289142B2 (ja) * 2014-02-07 2018-03-07 キヤノン株式会社 画像処理装置、画像処理方法、プログラムおよび記憶媒体
EP3109824B1 (en) * 2015-06-24 2019-03-20 RaySearch Laboratories AB System and method for handling image data
EP3206189B1 (en) * 2016-01-15 2018-11-14 Canon Kabushiki Kaisha Image processing apparatus, image processing method, and program
JP2017164198A (ja) * 2016-03-15 2017-09-21 キヤノン株式会社 情報処理システムおよび表示制御方法
JP6843632B2 (ja) * 2017-01-27 2021-03-17 キヤノン株式会社 音響波測定装置およびその制御方法
JP6934734B2 (ja) 2017-03-17 2021-09-15 キヤノン株式会社 画像処理装置、画像処理装置の制御方法およびプログラム
WO2019065466A1 (ja) * 2017-09-29 2019-04-04 キヤノン株式会社 画像処理装置、画像処理方法、およびプログラム
JP7262340B2 (ja) * 2019-07-31 2023-04-21 富士フイルムヘルスケア株式会社 超音波ｃｔ装置、画像処理装置、および、画像処理プログラム
CN118138889B (zh) * 2024-05-06 2024-06-28 南京维赛客网络科技有限公司 虚拟场景中摄像机的智能跟拍方法、系统及存储介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007159643A (ja) * 2005-12-09 2007-06-28 Canon Inc 画像処理装置および方法
JP2009157767A (ja) * 2007-12-27 2009-07-16 Nippon Telegr & Teleph Corp <Ntt> 顔画像認識装置、顔画像認識方法、顔画像認識プログラムおよびそのプログラムを記録した記録媒体
JP2010088627A (ja) * 2008-10-07 2010-04-22 Canon Inc 生体情報処理装置および生体情報処理方法
US20100266220A1 (en) * 2007-12-18 2010-10-21 Koninklijke Philips Electronics N.V. Features-based 2d-3d image registration

Family Cites Families (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP3871747B2 (ja) 1996-11-25 2007-01-24 株式会社日立メディコ 超音波診断装置
US7706633B2 (en) * 2004-04-21 2010-04-27 Siemens Corporation GPU-based image manipulation method for registration applications
US20100158332A1 (en) * 2008-12-22 2010-06-24 Dan Rico Method and system of automated detection of lesions in medical images
US20140082542A1 (en) * 2010-07-19 2014-03-20 Qview, Inc. Viewing and correlating between breast ultrasound and mammogram or breast tomosynthesis images
EP2579777B1 (en) * 2010-06-11 2019-07-31 The Florida International University Board of Trustees Second generation hand-held optical imager
US20130094742A1 (en) * 2010-07-14 2013-04-18 Thomas Feilkas Method and system for determining an imaging direction and calibration of an imaging apparatus
US9014456B2 (en) * 2011-02-08 2015-04-21 University Of Louisville Research Foundation, Inc. Computer aided diagnostic system incorporating appearance analysis for diagnosing malignant lung nodules
JP5950619B2 (ja) * 2011-04-06 2016-07-13 キヤノン株式会社 情報処理装置
US8948487B2 (en) * 2011-09-28 2015-02-03 Siemens Aktiengesellschaft Non-rigid 2D/3D registration of coronary artery models with live fluoroscopy images
US11109835B2 (en) * 2011-12-18 2021-09-07 Metritrack Llc Three dimensional mapping display system for diagnostic ultrasound machines
JP6304970B2 (ja) * 2013-08-09 2018-04-04 キヤノン株式会社 画像処理装置、画像処理方法
WO2015103388A1 (en) * 2014-01-02 2015-07-09 Metritrack, Inc. System and method for tracking completeness of co-registered medical image data

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007159643A (ja) * 2005-12-09 2007-06-28 Canon Inc 画像処理装置および方法
US20100266220A1 (en) * 2007-12-18 2010-10-21 Koninklijke Philips Electronics N.V. Features-based 2d-3d image registration
JP2011508620A (ja) * 2007-12-18 2011-03-17 コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ 特徴に基づいた２次元／３次元画像のレジストレーション
JP2009157767A (ja) * 2007-12-27 2009-07-16 Nippon Telegr & Teleph Corp <Ntt> 顔画像認識装置、顔画像認識方法、顔画像認識プログラムおよびそのプログラムを記録した記録媒体
JP2010088627A (ja) * 2008-10-07 2010-04-22 Canon Inc 生体情報処理装置および生体情報処理方法

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2017082078A1 (ja) * 2015-11-11 2017-05-18 ソニー株式会社 画像処理装置および画像処理方法
US11290698B2 (en) 2015-11-11 2022-03-29 Sony Corporation Image processing apparatus and image processing method
JP2017127623A (ja) * 2016-01-15 2017-07-27 キヤノン株式会社 画像処理装置、画像処理方法、およびプログラム
CN112529960A (zh) * 2020-12-17 2021-03-19 珠海格力智能装备有限公司 目标对象的定位方法、装置、处理器和电子装置

Also Published As

Publication number Publication date
JP6182045B2 (ja) 2017-08-16
US9311709B2 (en) 2016-04-12
US20150104091A1 (en) 2015-04-16

Similar Documents

Publication Publication Date Title
JP6289142B2 (ja) 画像処理装置、画像処理方法、プログラムおよび記憶媒体
JP6182045B2 (ja) 画像処理装置およびその方法
ES2706542T3 (es) Sistema de diagnóstico para mamografía multimodal
CN105025803B (zh) 从多个三维视图对大对象的分割
WO2010064348A1 (ja) 医用画像の位置あわせのための情報処理装置、情報処理方法、プログラム
US9396576B2 (en) Method and apparatus for estimating the three-dimensional shape of an object
JP5685133B2 (ja) 画像処理装置、画像処理装置の制御方法、およびプログラム
KR102273020B1 (ko) 의료 영상 정합 방법 및 그 장치
KR20140096919A (ko) 의료 영상 정합 방법 및 장치
JP2009273597A (ja) 位置合わせ処理装置、位置合わせ方法、プログラム、及び記憶媒体
KR20140032810A (ko) 의료 영상들의 정합 방법 및 장치
KR20150118484A (ko) 의료 영상 정합 방법 및 장치
KR102233966B1 (ko) 의료 영상 정합 방법 및 그 장치
JP2014530348A (ja) 元の放射線画像を更新するための放射線画像システム及び方法
KR102285007B1 (ko) 초음파 스캐너의 탐촉자의 위치 및 자세 추적을 이용한 초음파 영상 제공 장치 및 방법
US9990725B2 (en) Medical image processing apparatus and medical image registration method using virtual reference point for registering images
JP7023196B2 (ja) 検査支援装置、方法およびプログラム
JP6145870B2 (ja) 画像表示装置および方法、並びにプログラム
JP2020014711A (ja) 検査支援装置、方法およびプログラム
JP6461257B2 (ja) 画像処理装置およびその方法
KR20160057024A (ko) 마커리스 ３차원 객체추적 장치 및 그 방법
KR101282008B1 (ko) 초음파 영상을 이용한 운동상태의 장기 및 병변 위치추정시스템 및 위치추정방법과, 그 방법을 수행하는 명령어를 포함하는 컴퓨터 판독가능 기록매체
Punithakumar et al. Cardiac ultrasound multiview fusion using a multicamera tracking system
JP2008259706A (ja) 投影画像生成装置、方法およびそのプログラム
Zenbutsu et al. 3D ultrasound assisted laparoscopic liver surgery by visualization of blood vessels

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20161007

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20170526

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20170623

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20170721

R151 Written notification of patent or utility model registration

Ref document number: 6182045

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R151

LAPS Cancellation because of no payment of annual fees