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

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

Info

Publication number
JP2015146989A
JP2015146989A JP2014022779A JP2014022779A JP2015146989A JP 2015146989 A JP2015146989 A JP 2015146989A JP 2014022779 A JP2014022779 A JP 2014022779A JP 2014022779 A JP2014022779 A JP 2014022779A JP 2015146989 A JP2015146989 A JP 2015146989A
Authority
JP
Japan
Prior art keywords
image
alignment
imaging
subject
mri
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
JP2014022779A
Other languages
English (en)
Other versions
JP2015146989A5 (ja
JP6289142B2 (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 JP2014022779A priority Critical patent/JP6289142B2/ja
Priority to US14/613,464 priority patent/US9691150B2/en
Publication of JP2015146989A publication Critical patent/JP2015146989A/ja
Publication of JP2015146989A5 publication Critical patent/JP2015146989A5/ja
Application granted granted Critical
Publication of JP6289142B2 publication Critical patent/JP6289142B2/ja
Active 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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/60—Editing figures and text; Combining figures or text
    • 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
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Processing (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

【課題】 対象物である被検体のＭＲＩ画像やＣＴ画像などの３次元画像を、ＰＡＴ画像に対して高精度に位置合わせすること。
【解決手段】画像処理装置は、第一の撮像装置により対象物を撮像した第一の画像と、第二の撮像装置により対象物を撮像した第二の画像と、第二の画像撮像装置に対して位置が対応付けられた撮像部により対象物を撮像した第三の画像と、を取得する画像取得部と、第一の画像における対象物の観測情報と第三の画像における対象物の観測情報とが整合するように第一の画像における対象物と第二の画像における対象物との位置合わせをする位置合わせ部と、を備える。
【選択図】 図１

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年
しかし、医用画像収集装置（例えば、ＰＡＴ装置）により取得される画像と他の医用画像収集装置（例えば、ＭＲＩ）により取得される画像とでは画像化する特徴が異なるため、第一の画像（ＭＲＩ画像）に写る構造物と第二の画像（ＰＡＴ画像）上の構造物すべてが一致するわけではない。そのため、非特許文献１の技術を第二の画像（ＰＡＴ画像）と第一の画像（ＭＲＩ画像）の位置合わせに転用したとしても、画像マッチングのみで高精度な位置合わせを行うのは困難であった。一方、第二の画像（ＰＡＴ画像）には被検体の外形情報が十分に含まれていないため、非特許文献２の技術をそのまま用いるのは困難であった。
本発明は、対象物である被検体のＭＲＩ画像やＣＴ画像などの第一の画像（三次元画像（三次元断層画像））を、第二の画像（ＰＡＴ画像）に対して高精度に位置合わせする技術の提供を目的とする。
本発明の一様態に係る画像処理装置は、第一の撮像装置により対象物を撮像した第一の画像と、第二の撮像装置により前記対象物を撮像した第二の画像と、前記第二の画像撮像装置に対して位置が対応付けられた撮像手段により前記対象物を撮像した第三の画像と、を取得する画像取得手段と、
前記第一の画像における対象物の観測情報と前記第三の画像における対象物の観測情報とが整合するように前記第一の画像における対象物と前記第二の画像における対象物との位置合わせをする位置合わせ手段と、を備えることを特徴とする。
本発明によれば、対象物である被検体のＭＲＩ画像やＣＴ画像などの第一の画像（三次元画像）を、第二の画像（ＰＡＴ画像）に対して高精度に位置合わせすることが可能になる。
第１の実施形態の画像処理装置を含むモダリティシステムの構成例を示すブロック図。 医用画像ＤＢが保持する被検体のＭＲＩ画像を説明する図。 ＰＡＴによる被検者の撮像を説明する図。 ＰＡＴが撮像したＰＡＴ画像の例を示す図。 非保持状態における前面赤外線カメラの撮像画像ICAM1の一例を示す図。 第１の実施形態の画像処理装置の処理を説明する図。 表面形状の取得処理を説明する図。 変形ＭＲＩ画像とＰＡＴ画像の表示例を示す図。 非保持状態における位置合わせの詳細を説明する図。 仮想投影像生成部が部分表面領域を求める処理を説明する図。 ＭＲＩ画像内の被検体の体表近傍情報を利用したＭＩＰ像を示す図。 圧迫変形の推定の詳細を説明する図。 メッシュＭの生成方法を示す模式図。 保持板による圧迫変形シミュレーションを説明する図。 変形ＭＲＩ画像ID_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≦NS)を取得し、それら点群の位置をＭＲＩ画像座標系ＣMRIにおける三次元の位置座標ベクトルＶSk_MRIとして記録する。
Ｓ６０４で、ＰＡＴ画像取得部１０３は、ＰＡＴ１２が被検体を撮像したＰＡＴ画像を取得し、ＰＡＴ画像を変形画像評価部１１１、画像表示部１１２に出力する。また、ＰＡＴ画像のヘッダ部に含まれる付加的情報、例えば、座標変換行列ＴPtoD、ＴC1toD、ＴC2toD、ＴC3toDなどもＰＡＴ画像取得部１０３により変形画像評価部１１１に出力される。なお、ＰＡＴ画像取得部１０３が取得するＰＡＴ画像は、所定の波長に対して被検体内の光エネルギ堆積量分布を画像化した三次元画像とする。
Ｓ６０５で、カメラ画像取得部１０４は、ＰＡＴ１２の前面赤外線カメラ３０５、後面赤外線カメラ３０６および側面赤外線カメラ３０７が撮像した、非保持状態における被検体の赤外線カメラ画像と保持状態における被検体の赤外線カメラ画像とを取得する。そして、カメラ画像取得部１０４は、取得した赤外線カメラ画像を二次元形状取得部１０５、仮想投影像評価部１０８に出力する。ここで取得した赤外線カメラ画像は、保持状態における撮像画像ＩCAMi、および非保持状態における撮像画像I'CAMi(i=1, 2, 3)である。
ＰＡＴ画像取得部１０３およびカメラ画像取得部１０４は、ＰＡＴ１２の撮像に同期してＰＡＴ１２から直接画像を取得してもよいし、過去に撮像され記録された画像を図示しない医用画像記録装置から取得してもよい。ＰＡＴ画像取得部１０３、カメラ画像取得部１０４および先に説明した医用画像取得部１０１は画像を取得する画像取得部として機能する。
Ｓ６０６で、二次元形状取得部１０５は、入力された各赤外線カメラ画像に画像処理を施して、被検体の外形形状情報（表面形状情報）として、乳房輪郭形状（図５の５０１）と乳頭画像（図５の５０２）とを取得する。そして、二次元形状取得部１０５は、乳房輪郭形状および乳頭位置を剛体変換部１０６、変形画像評価部１１１に出力する。例えば、乳房輪郭形状の取得（検出）には一般的なエッジ検出手法を用いることができる。また、乳頭位置は、乳房領域の境界を表す曲線の曲率に基づき取得（検出）することができる。なお、乳房輪郭形状と乳頭位置の取得（検出）は、上記の方法に限られるわけではない。
Ｓ６０７で、剛体変換部１０６、仮想投影像生成部１０７、および仮想投影像評価部１０８は、ＭＲＩ画像と、非保持状態における赤外線カメラ画像に写った表在血管およびその周辺領域を含む線状領域の情報とを用いて、非保持状態における被検体とＭＲＩ画像との位置合わせを行う。ここで、非保持状態とは第二の状態（ＰＡＴ撮像時の被検体の位置や形状）の初期段階である。具体的には、剛体変換部１０６、仮想投影像生成部１０７、および仮想投影像評価部１０８は、仮定した位置合わせパラメータの候補値のそれぞれに基づき、赤外線カメラでＭＲＩ画像を仮想的に観測した仮想画像を生成し、これを赤外線カメラ画像と比較することで位置合わせパラメータを推定する。この位置合わせ（非保持状態における位置合わせ）の詳細は後述する。本ステップの処理により、ＭＲＩ画像座標ＣMRIから前面カメラ座標系ＣCAM1への剛体変換を表す座標変換行列ＴMtoC1が、位置合わせパラメータとして取得される。
Ｓ６０８で、剛体変換部１０６は、座標変換行列ＴMtoC1に基づき、ＭＲＩ画像座標系ＣMRIからＰＡＴ装置座標系ＣDEVへの剛体変換を表す座標変換行列ＴMtoDを算出する。剛体変換部１０６は、座標変換行列ＴMtoC1に、画像処理装置１０が保持する前面カメラ座標系ＣCAM1からＰＡＴ装置座標系CDEVへの座標変換行列ＴC1toDを掛けることで、座標変換行列ＴMtoDを算出する。
Ｓ６０９で、変形推定部１０９、変形画像生成部１１０、変形画像評価部１１１は、非保持状態における位置合わせ結果に基づき、保持状態における被検体とＭＲＩ画像の位置合わせ（以下、圧迫変形の推定）を行う。詳細は後述するが、物理変形シミュレーションを用いてＭＲＩ画像の圧迫変形を推定する。つまり、変形パラメータを様々に変化させて物理変形シミュレーションを行い、ＰＡＴ画像との比較によって変形の適切性を表す所定の評価値を得て、評価値を最小とする変形パラメータを位置合わせパラメータとして推定する。そして、推定した変形パラメータを用いて、変形ＭＲＩ画像ID_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にθ1を設定する。
次に、Ｓ９０４で、剛体変換部１０６は、平行移動後のＭＲＩ画像を、乳頭位置を基準に回転パラメータθiだけ回転移動したＭＲＩ画像IMRIonC1iと位置座標ベクトルvSk_MRIonC1iとを仮想投影像生成部１０７に出力する。つまり、剛体変換部１０６は、前面カメラ座標系ＣCAM1において回転パラメータθiだけ回転移動させる座標変換行列Ｔ２ｉを算出する。続いて、剛体変換部１０６は、ステップＳ９０１で導出した平行移動の座標変換行列Ｔ１MtoC1に座標変換行列Ｔ２ｉを掛けた剛体変換を表す座標変換行列ＴMtoC1iを導出する。そして、剛体変換部１０６は、ＭＲＩ画像および被検体の表面形状の位置座標ベクトルvSk_MRIを、それぞれ座標変換行列ＴMtoC1iを用いて座標変換したＭＲＩ画像IMRIonC1iと位置座標ベクトルvSk_MRIonC1iとを生成する。
次に、Ｓ９０５で、仮想投影像生成部１０７は、剛体変換されたＭＲＩ画像を前面赤外線カメラ３０５の視点から観測した場合に視野に入るであろう被検体の表面形状を部分表面領域として求める。言い換えれば、仮想投影像生成部１０７は、表面位置を示す情報である表面形状と、位置合わせパラメータの候補値である剛体変換パラメータ（剛体変換を表す座標変換行列ＴMtoC1i）に基づき、剛体変換されたＭＲＩ画像を前面赤外線カメラ３０５の視点から見た仮想画像を生成する。
図１０は仮想投影像生成部１０７が部分表面領域を求める処理(Ｓ９０５)を例示的に説明する図である。図１０は前面赤外線カメラ３０５の視点Ｐから観測されるＭＲＩ画像上の被検体の様子を示し、ＭＲＩ乳房領域１０００は、乳頭位置１００１が赤外線カメラにおける乳頭位置に一致するように剛体変換されている。ステップＳ９０５において、仮想投影像生成部１０７は、視点Ｐの位置・姿勢を基準に透視投影を行い、視点Ｐから観測される部分表面領域１００３を求める。仮想投影像生成部１０７は、視点Ｐからの投影線１００４を延長し、投影線１００４がＭＲＩ乳房領域１０００の体表形状を表す位置座標ベクトルvSk_MRIonC1iと交点をもつ観測範囲１００２を特定する。そして、仮想投影像生成部１０７は、観測範囲１００２において、位置座標ベクトルvSk_MRIonC1iの中から、投影線１００４が視点Ｐを起点にして最初に交わる体表点１００５を求める。そして、仮想投影像生成部１０７は、観測範囲１００２に含まれるすべての投影線１００４について体表点１００５を求めて部分表面領域１００３を決定する。
次に、仮想投影像生成部１０７は、剛体変換されたＭＲＩ画像における部分表面領域１００３の近傍情報を利用したＭＩＰ像を生成し、生成したＭＩＰ像を仮想投影像評価部１０８に出力する。ここで、ＭＩＰは「maximum intensity projection」の略であり、以下では、ステップＳ９０６で生成されるＭＩＰ像を「体表近傍ＭＩＰ像IMIPonC1i」と呼ぶ。
仮想投影像生成部１０７は、投影線１００４が最初に交わる体表点１００５を起点にして、視点Ｐから遠ざかる方向に所定の距離（例えば、５ｍｍ）をもつ体表近傍区間１００６を設定する。そして、仮想投影像生成部１０７は、観測範囲１００２に含まれるすべての投影線１００４について体表近傍区間１００６を設定することで体表近傍領域１００７を規定する。続いて、仮想投影像生成部１０７は、視点Ｐの位置・姿勢を基準に透視投影を行い、ＭＲＩ画像IMRIonC1iの中で体表近傍領域１００７に含まれる領域に限定したＭＩＰ像である体表近傍ＭＩＰ像IMIPonC1iを生成する。これにより、ＭＲＩ画像中の血管領域およびその周辺領域を含む線状領域（以下、血管領域１００８）の中で前面赤外線カメラ３０５側の体表近傍に存在する表在血管およびその周辺領域を含む線状領域（以下、表在血管１００９）の情報のみを可視化したＭＩＰ像が生成される。なお、体表近傍領域１００７の生成時に、被検体の皮膚が当該領域に含まれないようにしてもよい。具体的には、部分表面領域１００３に皮膚厚相当の所定の厚みを与えた領域を皮膚領域として導出し、上記の処理で得た体表近傍領域１００７からこれを除外してもよい。被検体の皮膚はＭＲＩ画像中において高い輝度値を有しているので、この除外処理により、生成するＭＩＰ像中において表在血管をより明瞭に描出できるようになる。
このように、体表近傍、かつ、前面赤外線カメラ３０５の視点Ｐから観測されるだろう領域の情報のみを可視化することで、乳房内部の血管領域や、前面赤外線カメラ３０５からみて反対側の体表近傍に存在する血管領域の可視化を防ぐことができる。つまり、実際の赤外線カメラ画像により近いＭＩＰ像である体表近傍ＭＩＰ像IMIPonC1iを生成することができる。
なお、体表近傍ＭＩＰ像IMIPonC1iの生成方法は、体表近傍の領域のみを、あるいは、体表近傍の領域を強調して可視化することができれば、上記の方法に限られるものではない。他の方法として、投影線１００４上において、体表点１００５からの距離が、視点Ｐから遠ざかる方向に大きくなるほど、輝度値に対する重みを小さくしてＭＩＰ像を生成する方法が挙げられる。この方法によれば、乳房内部の領域で、体表に近い領域ほど輝度値が強調されたＭＩＰ像が生成されるため、表在血管が強調されたＭＩＰ像が生成される。もちろん、本処理によるＭＩＰ像の生成を行う場合においても、皮膚領域を描画対象から除外することで、表在血管をより明瞭に描出できる。
図１１は、ＭＲＩ画像内の被検体の体表近傍情報を利用したＭＩＰ像１１０３を例示的に示す図である。図１１において、体表近傍ＭＩＰ像ＩMIPonC1iは、前面カメラ画像座標系ＣIMG1上の二次元画像として表される。これは、視点Ｐに基づく透視投影によって二次元のＭＩＰ像を生成する際の座標変換(Ｓ９０６)が、三次元の前面カメラ座標系ＣCAM1から二次元の前面カメラ画像座標系ＣIMG1への座標変換（カメラの撮像プロセス）と幾何的に等しいためである。
また、ＭＩＰ像の生成の際、体表の外側領域は処理対象に含まれないため、乳房輪郭形状１１００の外側には有意の輝度値は存在しない。さらに、ＭＲＩ画像における乳頭位置１００１が、前面赤外線カメラ３０５の撮像画像ＩCAM1における乳頭５０２の位置に一致するように剛体変換されているため、ＭＩＰ像の乳頭１１０１は、撮像画像ＩCAM1（図５）の乳頭５０２に一致する。そして、ＭＩＰ像においては、表在血管１１０２が特に高輝度な領域として可視化される。
次に、Ｓ９０７で、仮想投影像評価部１０８は、体表近傍ＭＩＰ像ＩMIPonC1iと、非保持状態の赤外線カメラ画像Ｉ'CAM1の両方において可視化された表在血管の輝度情報に基づき、両画像の類似度Ｓiを算出する。仮想投影像評価部１０８は、体表近傍ＭＩＰ像ＩMIPonC1i、非保持状態の赤外線カメラ画像Ｉ'CAM1ともに、類似度の算出領域から乳房の外側領域を除外し、類似度の算出領域を乳房の内部領域に限定する。
体表近傍ＭＩＰ像IMIPonC1iにおける表在血管１１０２（図１１）は、周囲の乳房領域に比べて高輝度に可視化されている。他方、赤外線カメラ画像Ｉ'CAM1における表在血管５０３（図５）は、周囲の乳房領域に比べて低輝度に可視化されている。そこで、仮想投影像評価部１０８は、それら画像の間の輝度情報を直接比較できるようにするため、体表近傍ＭＩＰ像ＩMIPonC1iの輝度値を反転する。そして、仮想投影像評価部１０８は、輝度値を反転した体表近傍ＭＩＰ像と赤外線カメラ画像Ｉ'CAM1の間で類似度Ｓｉ（0≦Si≦1）を算出する。体表近傍ＭＩＰ像ＩMIPonC1iの表在血管１１０２（図１１）と、赤外線カメラ画像Ｉ'CAM1の表在血管５０３（図５）とが類似すれば、類似度Ｓｉの値が高くなる（１に近付く）ものとする。
なお、本実施形態においては、類似度Ｓｉの評価尺度として画像間の相互情報量を適用するが、評価尺度はこれに限られず、相互相関係数やＳＳＤ (Sum of Squared Difference）など、既知の手法を用いてもよい。また、輝度値に直接的に基づく評価尺度である必要はなく、例えば両画像からエッジなどの画像特徴を検出して、それらの類似度や一致度を算出するような尺度でもよい。
Ｓ９０８で、仮想投影像評価部１０８は、類似度Ｓｉと類似度の最大値ＳMAXとを比較する。そして、類似度Ｓｉが最大値ＳMAXを超える（Si＞SMAX）場合（Ｓ９０８−Ｙｅｓ）、処理をステップＳ９０９に進める。
Ｓ９０９で、仮想投影像評価部１０８は、ＳMAXを更新し（SMAX=Si）、ＳMAXに対応する角度θMAXを更新する（θMAX=θi）。一方、ステップＳ９０８の判定で、類似度Ｓｉが最大値SMAX以下（Si≦SMAX）の場合（Ｓ９０８−Ｎｏ）、仮想投影像評価部１０８はステップＳ９０９の更新を行わずに処理をステップＳ９１０に進める。
Ｓ９１０で、仮想投影像評価部１０８は、ループ変数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に対して、ステップＳ９０５、Ｓ９０６で透視投影を行う視点Ｐを、後面赤外線カメラ３０６のカメラ視点または側面赤外線カメラ３０７のカメラ視点にそれぞれ置き換えればよい。
ここで、ＭＲＩ画像IMRIonC1iおよび位置座標ベクトルvSk_MRIonC1iはステップＳ９０４で前面カメラ画像座標系CIMG1に変換されている。従って、この透視投影に用いるカメラ視点の位置・姿勢としては、前面カメラ座標系ＣCAM1上で表現される位置姿勢を用いればよい。
上述したように、各赤外線カメラ（前面、側面、後面）の位置関係は、ＰＡＴ装置座標系ＣDEVを基準に対応付けられている。そのため、前面カメラ座標系ＣCAM1から後面カメラ座標系ＣCAM2または側面カメラ座標系ＣCAM3への座標変換行列を導出することができる。つまり、前面カメラ座標系ＣCAM1における後面または側面赤外線カメラのカメラ視点の位置姿勢を導出することができる。
前面、後面、側面の赤外線カメラそれぞれに基づく類似度を総合した総合類似度を算出して、類似度を評価してもよい。総合類似度としては、これら三種類の類似度の、例えば、重み付け平均値、最大値、最小値、中央値が挙げられる。また、上述の処理では、前面カメラ座標系ＣCAM1を基準にして、回転パラメータθiの設定、及びＭＲＩ画像を各カメラ座標系に透視投影するための座標変換を行った。しかし、ＰＡＴ装置座標系ＣDEVを基準にして、ＣMRIとCDEVとの間の剛体変換ＴMtoDを位置合わせパラメータとして、その候補値を設定するようにしてもよい。この場合、ＣMRIからＣDEVへの剛体変換に加えて各赤外線カメラへのビューイング変換を行った後に、ステップＳ９０６と同様の処理で体表近傍ＭＩＰ像を生成すればよい。そして、各赤外線画像との類似度評価に基づいて位置合わせパラメータを推定すればよい。
（圧迫変形の推定）
図１２のフローチャートにより圧迫変形の推定(Ｓ６０９)の詳細を説明する。
Ｓ１２０１で、変形推定部１０９は、ステップＳ６０２で取得された被検体の表面形状、および、ステップＳ６０８で取得された座標変換行列ＴMtoDを用いて、被検体の形状を表す三次元のメッシュ（以下、メッシュＭ）を生成する。変形推定部１０９は、座標変換行列ＴMtoDによる座標変換を、被検体の表面形状VSk_MRIに施して、ＰＡＴ装置座標系CDEVにおける被検体表面点群の位置座標ベクトルVSi_MRIonD（1≦i≦Ns）を算出する。そして、変形推定部１０９は、位置座標ベクトルVSi_MRIonDが表す表面形状に基づき、被検体の内部領域を判別し、内部領域にメッシュMを配置する。
図１３はメッシュＭの生成方法を例示的に説明する図である。図１３(ａ)は被検体の処理対象領域１３００のサジタル断面を示し、サジタル断面における被検体の表面位置１３０１と内部領域１３０２とを示す。図１３(ｂ)に示すように、メッシュＭは、被検体の内部領域１３０２に六面体や四面体などの立体構造をもつ要素１３０３を配置して生成される。メッシュＭは、これらの要素の頂点（ノード）１３０４の位置と連結情報によって記述される。
以下では、ステップＳ１２０１において配置されるメッシュＭのノード数をＮｍ、各ノードの位置をｓL（１≦Ｌ≦Ｎｍ）と表記する。各ノードの変位によって要素内の変位場を表現することができ、これに基づき、被検体内の任意の点の変位を求めることができる。
次に、Ｓ１２０２で、変形推定部１０９は、変形パラメータの各成分（被検体のヤング率やポアソン比など）が取り得る値を組み合わせた複数（Ｎｐ組）の変形パラメータの仮説ｐk（１≦ｋ≦Ｎｐ）を生成する。変形推定部１０９は、例えば、各成分が取り得る値の範囲を適切な間隔で分割し、それらすべてを組み合わせることで変形パラメータpkを生成する。例えば、ヤング率の比pyとポアソン比ppを変形パラメータpkの成分とし、ヤング率の比pyとポアソン比ppが取り得る値を次のような値とする。
py = {1, 2, 3, 4, 5}
pp = {0.0, 0.2, 0.4, 0.45, 0.499}
そして、ヤング率の比ｐyとポアソン比ｐpを組み合わせた変形パラメータを生成する。上記の例ではＮｐ＝２５（＝５×５）（組）になる。なお、ヤング率の比ｐyは、乳房の硬さの非等方性に対応するためのパラメータであり、人体のコロナル面（x-z平面）におけるヤング率に対する、人体の前後方向（y軸方向）のヤング率の割合を表す。
Ｓ１２０３で、変形推定部１０９は初期化処理を行う。変形推定部１０９は、ループ変数ｋに１を設定し、後述する評価値の最大値ＥMAXに０を設定し、後述する変形パラメータｐMAXにp1を設定する。
Ｓ１２０４で、変形推定部１０９は、変形パラメータｐkを用いて、有限要素法に基づく物理変形シミュレーションをメッシュＭに施し、変形後のメッシュである変形メッシュＤＭｋを生成する。このときの変形関数Ｆｋ(x, y, z)を、メッシュＭから変形メッシュＤＭｋへの各ノードの変位を与える変位ベクトルｄｋL（１≦Ｌ≦Ｎｍ）と定義する。
図１４は、ステップＳ１２０４における物理変形シミュレーションとして、保持板による圧迫変形シミュレーションを例示的に説明する図である。保持板による圧迫変形において、二枚の保持板を被検体の中心方向に移動すると、移動後の保持板に接する被検体の表面領域は保持板に押圧された（張り付いた）状態になる。
図１４（ａ）に示すように、二枚の保持板をそれぞれ、距離Δｄ１、Δｄ２分、矢印方向に移動すると仮定する。変形推定部１０９は、メッシュＭのノードのうち、体表面を表す表面ノードから、保持面ＰUd1、ＰLd2に接する外側表面ノード１４００と１４０１とを抽出し、外側表面ノード１４００および１４０１のそれぞれが保持面に張り付くための変位量を求める。そして、変形推定部１０９は、変位量を変形シミュレーションの境界条件Ｃとして有限要素法の計算を実行し、二枚の保持板がそれぞれ、距離Δｄ１、Δｄ２分、移動した場合の変形メッシュを生成する。
本実施形態において、変形推定部１０９は、図１４(ｂ)に示す二枚の保持板の最終的な保持位置ＰU、ＰLまで移動するまでの間、変形過程で生じる境界条件の変化に対応した複数回（Ｎ回）の変形シミュレーションを分割して行う。図１４（ｂ）はＮ回の変形シミュレーションを繰り返した結果の変形メッシュＤＭｋを表す。物理変形シミュレーションにより、図１４（ｂ）に示すように、被検体の変形メッシュが保持位置ＰU、ＰLの間でｚ軸方向に圧縮され、ｙ軸方向に伸長されることが、図１４（ａ）との比較により分かる。
Ｓ１２０５で、変形画像生成部１１０は、変形パラメータｐkに対応する変形をＭＲＩ画像に施して、変形ＭＲＩ画像を生成し、変形ＭＲＩ画像を変形画像評価部１１１に出力する。変形画像生成部１１０は、ＭＲＩ画像を、座標変換行列ＴMtoDを用いてＰＡＴ装置座標系ＣDEVに座標変換し、ステップＳ１２０４で算出した変形関数Ｆｋ(x, y, z)を用いて変形処理する。そして、変形画像生成部１１０は、座標変換行列ＴPtoDの逆行列を用いた座標変換により、ＰＡＴ画像座標系CPATにおける変形ＭＲＩ画像ID_MRIonPkを生成する。
図１５は、変形ＭＲＩ画像ＩD_MRIonPを例示的に示す図であり、図１５には変形ＭＲＩ画像８００、変形後の乳房領域１５０１、変形前の乳房形状１５０２が示されている。図１５（ａ）は変形ＭＲＩ画像ＩD_MRIonPをアキシャル断面でスライスした二次元画像であり、図１５（ｂ）は変形ＭＲＩ画像ＩD_MRIonPをサジタル断面でスライスした二次元画像である。変形後の乳房領域１５０１と変形前の乳房形状１５０２とを比較すると、ＰＡＴ画像座標系ＣPATのｚ軸方向への圧迫により、ｘ-ｙ平面上において乳房領域が伸長し、ｚ軸方向に圧縮されていることがわかる。
Ｓ１２０６で、変形画像評価部１１１は、ステップＳ６０４で取得されたＰＡＴ画像と、ステップＳ６０５で取得された保持状態の赤外線カメラ画像と、を用いて、変形ＭＲＩ画像ＩD_MRIonPkの適切性の評価値Ｅｋを算出する。そして、変形画像評価部１１１は、算出した評価値Ｅｋを変形推定部１０９に出力する。変形画像評価部１１１は、変形ＭＲＩ画像とＰＡＴ画像との間の類似度ＳMRIk（0≦SMRI、k≦1）、および、変形ＭＲＩ画像と保持状態の赤外線カメラ画像における乳房形状の間の残差Ｒｋに基づき、評価値Ｅｋを算出する。
なお、評価値Ｅｋの値が高いほど変形が適切であるものとする。また、類似度ＳMRIkの評価尺度として、ステップＳ９０７と同様に画像間の相互情報量を適用する。なお、評価尺度はこれに限らず、相互相関係数やＳＳＤ、血管分岐部等の特徴点の位置の一致度など、既知の何れの手法を用いてもよい。
図８は変形ＭＲＩ画像（ＩD_MRIonPk）８００とＰＡＴ画像４００の同一のアキシャル断面を示し、破線の矩形で示す対応領域８５０は、ＰＡＴ画像４００に対応する変形ＭＲＩ画像８００上の領域である。類似度ＳMRIkは、ＰＡＴ画像４００と対応領域８５０との間で計算される。ＰＡＴ画像４００と対応領域８５０との間で、可視化された、変形ＭＲＩ画像８００の血管領域８５１とＰＡＴ画像４００の血管領域８５２が類似すれば、類似度ＳMRIkの値が高くなる。
また、残差Ｒｋは、赤外線カメラ画像に写った被検体の輪郭（シルエット）の形状と、赤外線カメラ画像に投影した変形メッシュＤＭｋの外郭の形状との間の差として算出される。例えば、変形メッシュＤＭｋを保持状態の赤外線カメラ画像ＩCAM1へ投影する場合、変形メッシュＤＭｋは、座標変換行列ＴC1toDの逆行列を用いて前面カメラ座標系ＣCAM1に座標変換され、前面カメラ画像座標系ＣIMG1に投影変換される。また、同様の方法により、変形メッシュＤＭｋを後面カメラ画像ＩCAM2または側面カメラ画像ICAM3へ投影することができる。
残差Ｒｋは、例えば、変形メッシュＤＭｋを三つの赤外線カメラ画像それぞれに投影したメッシュの外郭と、三つの赤外線カメラ画像の間の各残差を総合した総合残差（例えば、三つの残差の重み付け平均値）として算出される。ただし、総合残差は、重み付け平均値に限らず、三つの残差の最大値、最小値、中央値などを用いてもよい。
また、残差Ｒｋは、例えば、赤外線カメラ画像に写った被検体の乳頭位置と、赤外線カメラ画像に投影した変形メッシュＤＭｋの乳頭位置との間の残差として算出してもよい。勿論、乳房形状の残差と、乳頭位置の残差との両方を統合した値（例えば、重み付け和）を残差Ｒｋとしてもよい。なお、乳房形状や乳頭位置に基づく残差は、赤外線カメラ画像から取得する例を説明したが、乳房の外観を撮像した一般的なカメラ画像から取得してもよい。
評価値Ｅｋは、例えば、類似度SMRIkと残差Ｒｋに基づく重み付け和として、下式で表される。
Ｅｋ = ａＳMRIk ＋ ｂ{1/(1+Rk)} ・・・(１)
ここで、ａ、ｂは重み係数（ａ＋ｂ=1）。
式(１)の第二項において(１＋Ｒｋ)の逆数を用いるのは次の理由からである。
・残差Ｒｋは、評価値Ｅｋとは逆に、変形が適切であるほど値が小さくなる指標である
・値の取り得る範囲を、類似度ＳMRIkと同様に、０から１の間にする
Ｓ１２０７で、変形推定部１０９は、入力される評価値Ｅｋと評価値の最大値ＥMAXとを比較する。評価値Ｅｋが評価値の最大値ＥMAXを超える（Ek＞EMAX）場合（Ｓ１２０７−Ｙｅｓ）、変形推定部１０９は、処理をステップＳ１２０８に進める。
Ｓ１２０８で、変形推定部１０９は、ＥMAXを（ＥMAX=Ｅｋ）と更新する。また、変形推定部１０９は、ＥMAXに対応する変形パラメータｐMAXを（ｐMAX＝ｐk）と更新する。
一方、ステップＳ１２０７の比較で、評価値Ｅｋが評価値の最大値ＥMAX以下（Ｅｋ≦ＥMAX）の場合（Ｓ１２０７−Ｎｏ）、変形推定部１０９は、ステップＳ１２０８の更新を行わずに処理をステップＳ１２０９に進める。
次に、Ｓ１２０９で、変形推定部１０９は、ループ変数ｋをインクリメントする。そして、Ｓ１２１０で、変形推定部１０９は、ループ変数ｋと仮説の総数Ｎｐとを比較する。ループ変数ｋが仮説の総数Ｎｐ以下（ｋ≦Ｎｐ）の場合（Ｓ１２１０−Ｎｏ）、変形推定部１０９は、処理をステップＳ１２０４に戻し、ループ変数ｋが仮説の総数Ｎｐを超える場合（ｋ＞Ｎｐ）の場合（Ｓ１２１０−Ｙｅｓ）、変形推定部１０９は、処理をステップＳ１２１１に進める。つまり、仮説の総数Ｎｐ分、ステップＳ１２０４からＳ１２１０の処理が繰り返される。
仮説の総数Ｎｐ分の処理が終了すると、Ｓ１２１１で、変形推定部１０９は、変形パラメータｐMAXを変形画像生成部１１０に出力する。言い換えれば、複数の変形パラメータから、評価値の最大値ＥMAXに対応する変形パラメータｐMAXが選択される。
Ｓ１２１２で、変形画像生成部１１０は、変形パラメータｐMAXに対応する変形ＭＲＩ画像ＩD_MRIonPMAXを圧迫変形の推定結果（ID_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)を最小化する（つまりSMIP(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（ΔxMAX）を求める。このΔxMAXは、f(x(k))の変化率が最大になるパラメータ空間の方向ベクトルであり、grad f(x(k))に等しい。
また、最適化アルゴリズムとして、ニュートン法など既知の手法の何れを用いてもよい。これにより、より少ない繰り返し回数で、最適な評価値を与える回転角度を推定することができ、処理の高速化を図ることができる。
（変形例２）
上記の剛体変換の推定処理（Ｓ６０８）では、ＭＲＩ画像の剛体変換の妥当性を評価する評価値として、赤外線カメラに写った内部情報との類似度を用いた。しかし、別の方法で評価値を求めてもよい。
例えば、ＭＲＩ画像から生成した体表近傍ＭＩＰ像に写る被検体の輪郭（シルエット）形状と赤外線カメラに写った被検体の輪郭形状に対する形状残差を算出し、これを評価値としても良い。
つまり、ステップＳ９０７において、体表近傍ＭＩＰ像ID_MIPonC1k及び赤外線カメラ画像Ｉ'CAM1の夫々から被検体である２次元の乳房形状を抽出し、それらの残差を評価値Ｅｋとする。このとき、体表近傍ＭＩＰ像ID_MIPonC1kは、生成元のＭＲＩ画像と同様、乳房外部は低輝度であるのに対し、乳房内部はそれよりも高輝度であり、両者の境界は明瞭である。従って、体表近傍ＭＩＰ像ID_MIPonC1kからの乳房形状の抽出は、ステップＳ６０２におけるＭＲＩ画像からの表面形状取得方法と同様の方法で行うことができる。また。赤外線カメラ画像Ｉ'CAM1からの乳房形状の抽出は、ステップＳ６０６における乳房輪郭形状の抽出方法と同様の方法で行うことができる。さらに、評価値Ｅｋの算出方法は、ステップＳ１２０６における変形ＭＲＩ画像と保持状態の赤外線カメラ画像における乳房形状の残差Ｒｋと同様の方法で算出することができる。
これにより、赤外線カメラに写った内部情報を使わずに、もしくは被検体の内部情報が写らない一般的なカメラを用いた場合でも、ＭＲＩ画像の剛体変換の処理を行うことができる。さらに、ステップＳ１２０６で算出する残差Ｒｋも、赤外線カメラではなく一般的なカメラを用いて算出可能であるため、処理全体を通して赤外線カメラを一般的なカメラに置き換えることができる。また、表在血管を有さない乳房以外の部位を対象物体とすることができる。
なお、赤外線カメラに写った内部情報の類似度と、輪郭形状の誤差を組み合わせた値を評価値としても良い。これは、ステップＳ１２０６における類似度ＳMRIkを赤外線カメラに写った内傍情報との類似度、残差Ｒｋを輪郭形状の形状残差に置き換えることで、ステップＳ１２０６と同様に算出することができる。
（変形例３）
上記の圧迫変形の推定処理（Ｓ６０９）では、ＭＲＩ画像の圧迫変形の妥当性を評価する評価値として、ＰＡＴ画像との類似度と、赤外線カメラに写った被検体の輪郭（シルエット）形状に対する形状誤差を用いた。しかし、別の方法で評価値を求めてもよい。
例えば、変形メッシュに基づき変形ＭＲＩ画像を生成し、変形ＭＲＩ画像を赤外線カメラ視点から投影した変形後の体表近傍ＭＩＰ像を生成し、体表近傍ＭＩＰ像と赤外線カメラ画像の類似度を算出して、この類似度を評価値に加えてもよい。
つまり、ステップＳ１２０５において、メッシュＭから変形メッシュＤＭｋへの各ノードの変位を表す変形関数Ｆｋを生成する。そして、ＭＲＩ画像に対して、ＭＲＩ画像座標系ＣMRIから前面カメラ装置座標系CCAM1への座標変換行列ＴMtoC1による剛体変換を施し、変形関数Ｆｋによる変形を施して、変形ＭＲＩ画像ID_MRIonC1kを生成する。
次に、ステップＳ９０５、Ｓ９０６と同様の方法により、変形ＭＲＩ画像ＩD_MRIonC1kから表面部分領域の近傍情報を利用した体表近傍変形ＭＩＰ像ID_MIPonC1kを生成する。そして、体表近傍ＭＩＰ像ＩD_MIPonC1kと保持状態の赤外線カメラ画像Ｉ'CAM1との間の類似度ＳMIPkを算出する。最後に、評価値ＥｋをステップＳ１２０６における類似度ＳMRIk、残差Ｒｋ、および、算出した類似度ＳMIPkに基づく重み付け和を下式により算出する。
Ｅｋ = ａSMRIk+ｂ{1/(1+Rk)}+ｃSMIPk …(3)
ここで、ａ、ｂ、ｃは重み係数（ａ＋ｂ＋ｃ＝１）。
これにより、体表近傍の表在血管の情報をさらに利用して変形パラメータを求めることができ、高精度な位置合わせが可能になる。
（変形例４）
上記の説明では、表在血管の可視化に、近赤外線による赤外線カメラ画像を利用する例を説明したが、例えば、体内における内部反射によって得られる偏光成分により表在血管を可視化した画像を利用することもできる。例えば、ハロゲンライトなどの光を体表に照射し、体表で反射する表面反射成分と、一旦、体内に侵入して体内で吸収と散乱を経た後に再び体表に出てくる内部反射成分とを区別する。こうすれば、体内の情報を可視化した画像が得られる。当該画像においては、赤外線カメラ画像と同様、内部反射成分のうち、血管内部のヘモグロビンの吸収箇所が周囲よりも暗く描出される。
このように、第一の撮像装置によって対象物を撮像した三次元画像と、第二の撮像装置によって対象物の表層部を撮像した二次元画像とが取得され、三次元画像から対象物の表面位置を示す情報が取得される。そして、表面位置を示す情報に基づいて、三次元画像を第二の撮像装置の視点から見た投影像が生成され、投影像および二次元画像を用いて、対象物に関して三次元画像と二次元画像との位置合わせが行われる。
つまり、本実施形態の画像処理装置は、第一の撮像装置（例えば、ＭＲＩやＣＴ）により対象物を撮像した第一の画像と、第一の撮像装置とは異なる第二の撮像装置１２０により対象物を撮像した第二の画像と、第二の画像撮像装置に対する位置関係が対応付けられている撮像部（３０５、３０６、３０７）により対象物を撮像した第三の画像と、を取得する画像取得部（１０１、１０３、１０４）と、第一の画像と第三の画像との比較結果を用いて、第一の画像と前記第二の画像との位置合わせをする位置合わせ部１１３と、を備える。これにより、対象物である被検体のＭＲＩ画像やＣＴ画像などの三次元画像を、ＰＡＴ画像に対して高精度に位置合わせすることが可能になる。
（第２の実施形態）
以下、本発明にかかる第２の実施形態の画像処理を説明する。なお、第２の実施形態において、第１の実施形態と同様の構成については、同一符号を付して、その詳細説明を省略する。
第１の実施形態では、非保持状態の被検体の画像とＭＲＩ画像との間で剛体変換パラメータを求め（非保持状態における位置合わせ）、その後、保持状態の被検体の画像とＭＲＩ画像との間で変形パラメータを推定（圧迫変形の推定）する方法を説明した。この方法は、非保持状態の被検体とＭＲＩ画像の被検体の形状とが異なる場合に生じる誤差や、非保持状態から保持状態へ移行する際の被検体の微動によって生じる誤差を含む可能性がある。
第２の実施形態では、非保持状態の被検体に関する計測情報を用いずに、保持状態の被検体（第二の状態の被検体）の画像とＭＲＩ画像（第一の状態の被検体）との間で、剛体変換パラメータと変形パラメータを同時に推定する方法を説明する。
図１６は、第２の実施形態の画像処理装置１０を含むモダリティシステムの構成例を示す図である。図１に示す第１の実施形態の構成と異なるのは、第１の実施形態の仮想投影像評価部１０８と変形画像評価部１１１との代わりに、評価部１１４を備える点である。また、第２の実施形態においては、位置合わせ部１１３内における情報の流れが異なるが、その点は、各部の動作と処理の説明において詳述する。
図１７のフローチャートにより第２の実施形態の画像処理装置１０の各部の動作と処理を説明する。なお、ステップＳ６０１からＳ６０６、Ｓ６１０における動作と処理は第１の実施形態と同様であり、その詳細説明を省略する。
ステップＳ１７０１において、赤外線カメラ画像、ＰＡＴ画像、ＭＲＩ画像の血管情報に基づき、ＭＲＩ画像における被検体の位置・姿勢と圧迫変形が推定される。この処理は、剛体変換部１０６、仮想投影像生成部１０７、変形推定部１０９、変形画像生成部１１０、評価部１１４によって行われる。
第２の実施形態では、赤外線カメラに対するＭＲＩ画像上の被検体の位置・姿勢と、被検体の圧迫変形を表す変形パラメータとを位置合わせパラメータとする。変形画像生成部１１０は、仮定した位置合わせパラメータの候補値に基づきＭＲＩ画像を赤外線カメラ座標系（前面カメラ座標系）に座標変換した後に圧迫変形を施すことで、当該位置合わせパラメータの候補値に基づく変形ＭＲＩ画像を生成する。
次に、仮想投影像生成部１０７は、生成された変形ＭＲＩ画像に対して、赤外線カメラの視点を基準に透視投影を行ったＭＩＰ像を生成する。その際、仮想投影像生成部１０７は、赤外線カメラの視点からの各投影線上において、乳房の三次元表面の近傍領域の輝度情報のみを可視化、あるいは、強調して可視化する画像処理を施して、乳房の表在血管を強調したＭＩＰ像を生成する。
次に、剛体変換部１０６は、変形ＭＲＩ画像を、ＰＡＴ１２に対して予め校正された赤外線カメラの位置・姿勢に基づき、赤外線カメラ座標系（前面カメラ座標系）からＰＡＴ画像座標系に座標変換する。
次に、評価部１１４は、血管情報が可視化された赤外線カメラ画像とＭＩＰ像の間、および、ＰＡＴ画像と座標変換後の変形ＭＲＩ画像との間の画像の類似度を算出し、それら類似度を総合した評価値を算出する。そして、仮定した位置合わせパラメータを変動させて、評価値が最大になる位置合わせパラメータを選択する。すなわち、当該位置合わせパラメータによって、ＭＲＩ画像とＰＡＴ画像の間で圧迫変形を含む位置合わせを行う。
図１８、図１９のフローチャートにより位置・姿勢と圧迫変形の推定(Ｓ１７０１)の詳細を説明する。
Ｓ１８０１で、変形推定部１０９は、ステップＳ６０２で取得された被検体の表面形状に基づき、被検体の形状を表すメッシュＭを生成する。この処理は、第１の実施形態のステップＳ１２０１と略同様であり、詳細説明を省略する。
Ｓ１８０２で、剛体変換部１０６は、ＭＲＩ画像を赤外線カメラ座標系に平行移動する。ステップＳ１８０２の処理は、ステップＳ６０３で取得されたＭＲＩ画像における乳頭位置と、ステップＳ６０６で取得された保持状態における赤外線カメラ画像上の乳頭位置に基づき行われる。この処理は、第１の実施形態のステップＳ９０１の処理と略同様であり、詳細説明を省略する。ただし、ステップＳ９０１の処理は非保持状態における赤外線カメラ画像に基づき行われるが、ステップＳ１８０２の処理は保持状態における赤外線カメラ画像に基づき行われる点で異なる。
次に、Ｓ１８０３で、剛体変換部１０６と変形推定部１０９は、剛体変換パラメータが取り得る値と、変形パラメータの各成分が取り得る値を組み合わせて複数（Ｎｔ組）の位置合わせパラメータ（変換パラメータ）の仮説ti（1≦i≦Nt）を設定する。
例えば、剛体変換部１０６は、剛体変換パラメータが取り得る値として、第１の実施形態のステップＳ９０２と同様に、複数（Ｎθ組）の回転パラメータθj（1≦j≦Nθ）を設定する。変形推定部１０９は、変形パラメータが取り得る値として、第１の実施形態のステップＳ１２０２と同様に、複数（Ｎｐ組）の変形パラメータpk（1≦k≦Np）を設定する。そして、回転パラメータθjと変形パラメータpkを組み合わせて複数（Ｎｔ=Ｎθ×Np組）の変換パラメータｔi（１≦i≦Ｎｔ）を設定する。なお、ＰＡＴ画像座標系と前面カメラ座標系との関係は既知であるので、これは、ＭＲＩ画像（第一の状態の被検体）からＰＡＴ画像（第二の状態の被検体）への位置合わせパラメータの候補値を設定することと等価である。また、変換パラメータｔiは、剛体変換部１０６と変形推定部１０９との間で共有されるものとする。
Ｓ１８０４で、評価部１１４は初期化処理を行う。つまり、評価部１１４は、ループ変数iに1を設定し、後述する評価値の最大値ＥMAXに０を設定し、後述する変換パラメータｔMAXにt1を設定する。
Ｓ１８０５で、剛体変換部１０６は、平行移動後のＭＲＩ画像を、乳頭位置を基準に変換パラメータｔi（つまりθj）だけ回転移動したＭＲＩ画像ＩMRIonC1iを生成する。そして、剛体変換部１０６は、ＭＲＩ画像IMRIonC1iと座標変換行列ＴMtoC1iを変形推定部１０９と変形画像生成部１１０に出力する。この処理は、第１の実施形態のステップＳ９０４と略同様であり、詳細説明を省略する。
Ｓ１８０６で、変形推定部１０９は、変換パラメータｔi（つまりθjとpk）を用いて、有限要素法に基づく物理変形シミュレーションをメッシュＭに施した変形メッシュＤＭｉを生成し、変形関数Ｆｊｋを変形画像生成部110に出力する。つまり、変形推定部１０９は、ステップＳ１８０５で導出された座標変換行列ＴMtoC1iを用いて、回転パラメータθjに対応する剛体変換をメッシュｍに施したメッシュＭｉを生成する。そして、変形推定部１０９は、メッシュＭｉに物理変形シミュレーションを施すことで、変形後のメッシュである変形メッシュＤＭｉを生成する。このときの変形関数Ｆｉ(x, y, z)を、メッシュＭから変形メッシュＤＭｉへの各ノードの変位を与える変位ベクトルｄｉL（１≦Ｌ≦Ｎｍ）と定義する。なお、この処理は、第１の実施形態のステップＳ１２０４と同様である。
Ｓ１８０７で、変形画像生成部１１０は、変換パラメータｔi（つまりθjとpk）に対応する変換をＭＲＩ画像に施した変形ＭＲＩ画像を生成し、変形ＭＲＩ画像を仮想投影像生成部１０７に出力する。つまり、変形画像生成部１１０は、ＭＲＩ画像に座標変換行列ＴMtoC1iを用いる剛体変換を施して、ＭＲＩ画像を回転パラメータθjに対応する前面カメラ座標系ＣCAM1に座標変換する。変形画像生成部１１０は、座標変換後のＭＲＩ画像に変形関数Ｆｊｋを用いる変形処理を施して、変形ＭＲＩ画像ＩD_MRIonC1iを生成する。
Ｓ１８０８で、仮想投影像生成部１０７は、前面カメラ座標系ＣCAM1において、変形ＭＲＩ画像ＩD_MRIonC1iを前面赤外線カメラ３０５の視点から観測した場合に視野に入るだろう被検体の表面形状（部分表面領域）を求める。この処理は、第１の実施形態のステップＳ９０５における剛体変換されたＭＲＩ画像を、変形ＭＲＩ画像ＩD_MRIonC1iに置き換えたものであり、詳細説明を省略する。
Ｓ１８０９で、仮想投影像生成部１０７は、変形ＭＲＩ画像ＩD_MRIonC1iにおける表面部分領域の近傍情報を利用した体表近傍変形ＭＩＰ像ID_MIPonC1iを生成し、体表近傍変形ＭＩＰ像を評価部１１４に出力する。この処理は、第１の実施形態のステップＳ９０６における剛体変換されたＭＲＩ画像を、変形ＭＲＩ画像ＩD_MRIonC1iに置き換えたものであり、詳細説明を省略する。
Ｓ１８１０で、変形画像生成部１１０は、ステップＳ１８０７で生成した変形ＭＲＩ画像ＩD_MRIonC1iを前面カメラ座標系からＰＡＴ画像座標系に座標変換した変形ＭＲＩ画像ID_MRIonPiを生成し、変形ＭＲＩ画像を評価部１１４に出力する。つまり、変形画像生成部１１０は、変形ＭＲＩ画像ＩD_MRIonC1iを、座標変換行列ＴC1toDを用いてＰＡＴ装置座標系CDEVに座標変換する。更に、変形画像生成部１１０は、座標変換行列ＴPtoDの逆行列を用いる座標変換により、ＰＡＴ画像座標系ＣPATにおける変形ＭＲＩ画像ID_MRIonPiを生成する。
Ｓ１８１１で、評価部１１４は、体表近傍変形ＭＩＰ像ＩD_MIPonC1iと保持状態の赤外線カメラ画像の類似度ＳMIPi（0≦SMIPi≦1）、および、変形ＭＲＩ画像ID_MRIonPiとＰＡＴ画像の類似度ＳMRIi（0≦SMRIi≦1）を算出する。さらに、評価部１１４は、変形ＭＲＩ画像ＩD_MRIonPiと保持状態の赤外線カメラ画像における乳房形状の残差Ｒｉを算出する。そして、評価部１１４は、それら類似度と残差を総合した評価値Ｅｉを算出する。
類似度ＳMIPiの算出方法は、第１の実施形態のステップＳ９０７における体表近傍ＭＩＰ像ＩMIPonC1jと非保持状態の赤外線カメラ画像を、体表近傍変形ＭＩＰ像ID_MIPonC1iと保持状態の赤外線カメラ画像に置き換えたものである。また、類似度ＳMRIiの算出方法は、第１の実施形態のステップＳ１２０６における変形ＭＲＩ画像ＩD_MRIkを、変形ＭＲＩ画像ID_MRIonPiに置き換えたものである。また、評価部１１４は、赤外線カメラ画像に写った被検体の輪郭（シルエット）の形状と、赤外線カメラ画像に投影した変形メッシュＤＭｊｋの外郭の形状との間の差として残差Ｒｉを算出する。残差Ｒｉの求め方はステップＳ１２０６と同様である。
評価値Ｅｋは、例えば、類似度ＳMIPi、ＳMRIiと残差Ｒｉに基づく重み付け和として、下式で表される。
Ｅｉ=ａSMIPi+ｂSMRIi＋ｃ{1/(1+Ri)} …(４)
ここで、ａ、ｂ、ｃは重み係数（ａ+ｂ+ｃ＝１）。
また、式(４)の第三項において（１＋Ｒｋ）の逆数を用いるのはステップＳ１２０６と同様の理由からである。
次に、Ｓ１８１２で、評価部１１４は、評価値Ｅｉと評価値の最大値ＥMAXを比較する。
そして、評価値Ｅｉが評価値の最大値ＥMAXを超える（Ei＞EMAX）場合（Ｓ１８１２−Ｙｅｓ）、評価部１１４は、処理をステップＳ１８１３に進める。
Ｓ１８１３で、評価部１１４は、ＥMAXを（ＥMAX＝Ｅｉ）と更新する。また、評価部１１４は、ＥMAXに対応する変換パラメータｔMAXを（ｔMAX＝ｔi）と更新する。
一方、ステップＳ１８１２の比較で、評価値Ｅｉが評価値の最大値ＥMAX以下（Ｅｉ≦ＥMAX）の場合（Ｓ１８１２−Ｎｏ）、評価部１１４は、Ｓ１８１３の更新を行わずに処理をステップＳ１８１４に進める。
Ｓ１８１４で、評価部１１４は、ループ変数iをインクリメントする。そして、Ｓ１８１５で、評価部１１４は、ループ変数ｉと仮説の総数Ｎｔとを比較する。ループ変数iが仮説の総数Ｎｔ以下（ｉ≦Ｎｔ）の場合（Ｓ１８１５−Ｎｏ）、評価部１１４は処理をステップＳ１８０５に戻す。一方、ステップＳ１８１５の比較で、ループ変数iが仮説の総数Ntを超える（ｉ＞Ｎｔ）の場合（Ｓ１８１５−Ｙｅｓ）、評価部１１４は処理をステップＳ１８１６に進める。つまり、仮説の総数Ｎｔ分、ステップＳ１８０５からＳ１８１５の処理が繰り返される。
仮説の総数Ｎｔ分の処理が終了すると、Ｓ１８１６で、評価部１１４は、評価値の最大値ＥMAXに対応する変換パラメータｔMAX（つまりθMAXとｐMAX）を変形画像生成部１１０に出力する。言い換えれば、複数の変換パラメータから、評価値の最大値ＥMAXに対応する変換パラメータｔMAXが選択される。
Ｓ１８１７で、変形画像生成部１１０は、変換パラメータｔMAXに対応する変形ＭＲＩ画像ＩD_MRIonPを画像表示部１１２に出力する。
以上で、剛体変換部１０６、仮想投影像生成部１０７、変形推定部１０９、変形画像生成部１１０、評価部１１４による位置・姿勢と圧迫変形の推定(Ｓ１７０１)が終了する。この処理によれば、様々な回転と変形を表す変換パラメータｔiを仮定して圧迫変形を含む変換を行った結果の中で変形の適切性の評価値Ｅｉが最大になる変換パラメータｔMAXによって、変形ＭＲＩ画像ID_MRIonPが生成される。
このように、ＭＲＩ画像を圧迫変形して生成した表在血管を強調したＭＩＰ像、保持状態の赤外線カメラ画像、ＰＡＴ画像の血管情報を比較して、ＭＲＩ画像上の被検体の位置・姿勢と圧迫変形を推定する。これにより、ＰＡＴ画像の撮像時、理想的な非保持状態の赤外線カメラ画像が撮像されていない場合や、非保持状態から保持状態に移る期間に乳房の位置・姿勢が変化したりする場合も、高精度な変形位置合わせを実現することができる。
（変形例）
上記の処理では、ＰＡＴ画像と赤外線カメラ画像の血管情報を同じ評価値Ｅｉに組み込んで利用する例を説明した。しかし、これら情報を同時に利用する必要はなく、例えば、ＰＡＴ画像の血管情報を先に利用し、その後、赤外線カメラ画像の血管情報を利用して、変形位置合わせを実施してもよい。
例えば、変換パラメータｔiを変動させ、変形ＭＲＩ画像とＰＡＴ画像との類似度を評価し、評価値が最大になる変換パラメータｔMAXを求める。次に、変換パラメータｔMAX付近に限定して変換パラメータtiを変動させて、体表近傍変形ＭＩＰ像ＩMIPonC1iと赤外線カメラ画像の類似度を評価し、評価値が最大になる変換パラメータtMAX2を求める。そして、変換パラメータtMAX2に対応する変形MRI画像を圧迫変形位置合わせ結果とする。これにより、ＭＲＩ画像に対して、形態画像として大まかな血管情報が描出されたＰＡＴ画像を利用して粗い変形位置合わせを実施した後、機能画像としてより詳細な血管情報が描出された赤外線カメラ画像を用いて精密な変形位置合わせを実施することができる。
また、第１の実施形態と同様に、非保持状態における剛体変換パラメータの推定を行った後に、θjの仮説の設定範囲を剛体変換パラメータの近傍に限定して、保持状態における剛体変換パラメータと変形パラメータを推定する構成でもよい。これにより、仮説の総数が減少し、処理の高速化を図ることができる。勿論、第１の実施形態の変形例１と同様に、仮説の総当りによる評価ではなく、一般的な最適化アルゴリズムを用いて変換パラメータを求める構成でもよい。この場合、非保持状態で推定した剛体変換パラメータをθjの初期値に用いることができる。
（第３の実施形態）
以下、本発明にかかる第３の実施形態の画像処理を説明する。なお、第３の実施形態において、第１、第２の実施形態と略同様の構成については、同一符号を付して、その詳細説明を省略する。
第１、第２の実施形態においては、ＰＡＴ１２が搭載する赤外線カメラによる乳房の撮像画像、および、ＭＲＩ画像から生成した乳房の表在血管の強調ＭＩＰ像に基づき、ＭＲＩ画像からＰＡＴ画像への位置合わせを行う例を説明した。しかし、ＭＲＩ画像などの予め撮像された被検体の三次元画像と赤外線カメラに写った被検体の表在血管の情報を用いて、当該三次元画像と位置合わせを行う対象画像はＰＡＴ画像に限られない。
例えば、診断支援を目的として超音波画像の撮像断面に対応する断面（以下、対応断面）の画像を、三次元画像であるＣＴ画像やＭＲＩ画像から生成して（切り出して）提示する場合が挙げられる。
第３の実施形態では、超音波画像、予め位置関係が対応付けられた赤外線カメラによる被検体の撮像画像、ＭＲＩ画像から生成した被検体の表在血管の強調ＭＩＰ像に基づき、ＭＲＩ画像撮像時の被検体から超音波画像撮像時の被検体への位置合わせ例を説明する。なお、以下では、乳房を被検体とし、ＭＲＩ画像を撮像した際の被検体を「第一の状態における被検体」、超音波画像を撮像した際の被検体を「第二の状態における被検体」と呼ぶ場合がある。
また、第３の実施形態において予め撮像されるＭＲＩ画像は、第一の状態における被検体として伏臥状態の乳房を撮像した画像であるが、超音波画像は、第二の状態における被検体として仰臥状態の乳房を撮像したものである。従って、ＭＲＩ画像と超音波画像の間での位置合わせを行うために、両画像間の剛体変換だけでなく、重力変形も推定する必要がある。以下、第３の実施形態の処理の流れを説明する。
超音波撮像装置に搭載された超音波探触子を仰臥状態の被検体に接触させて、超音波画像が撮影される。超音波探触子には、その位置と姿勢を計測するセンサ（磁気式、光学式など）が装着され、超音波撮影が行われている最中の超音波探触子の位置・姿勢が計測される。つまり、センサが基準とする座標系（以下、センサ座標系）における超音波画像の撮影領域が計測される。
また、被検体を撮影するように設置された赤外線カメラは、センサ座標系において、位置・姿勢の校正が行われている。従って、センサ座標系を介して、超音波画像と赤外線カメラ画像の位置・姿勢の対応付けが可能になる。そして、赤外線カメラによって撮像された、仰臥状態の乳房の表在血管が写った赤外線カメラ画像を取得する。
次に、第１、第２の実施形態と同様に、被検体を撮像したＭＲＩ画像から乳房の三次元表面形状を抽出する。続いて、赤外線カメラに対してＭＲＩ画像上の被検体の様々な位置・姿勢と重力変形を表す位置合わせパラメータの候補値を仮定し、位置合わせパラメータの候補値ごとに重力変形位置合わせを施した変形ＭＲＩ画像を生成する。
次に、変形ＭＲＩ画像ごとに画像内の体表近傍の情報を利用して赤外線カメラ側の表在血管を強調したＭＩＰ像を生成する。続いて、各ＭＩＰ像と赤外線カメラ画像の類似性を最大にする、赤外線カメラに対するＭＲＩ画像上の被検体の位置・姿勢と重力変形を推定する。
次に、予め対応付けられた超音波画像に対する赤外線カメラ画像の位置・姿勢に基づき、ＭＲＩ画像を超音波画像の座標系へ変形位置合わせした変形ＭＲＩ画像を生成し、超音波画像との間の変形位置合わせを行う。
このように、ＭＲＩ画像を重力変形させて赤外線カメラ側の表在血管を強調したＭＩＰ像と仰臥状態の赤外線カメラ画像を比較することで、赤外線カメラに対するMRI画像上の被検体の位置・姿勢と重力変形を推定することができる。従って、赤外線カメラとＭＲＩ画像の高精度な変形位置合わせを実現でき、超音波画像と赤外線カメラの位置関係を利用して、超音波画像とＭＲＩ画像の高精度な変形位置合わせを行うことができる。
（変形例１）
第３の実施形態では、伏臥状態で撮像されたＭＲＩ画像と、仰臥状態で撮像された超音波画像の間の乳房の変形状態を補正するために、赤外線カメラに対するMRI画像上の被検体の位置・姿勢に加え、重力変形も推定した。
しかし、例えば、超音波画像と位置合わせする対象であるＣＴ画像やＭＲＩ画像などの三次元画像が、超音波画像と同じ仰臥状態で撮像されていれば、両画像間において重力変形による補正を施す必要はない。この場合の処理を、第３の実施形態の処理との差分（赤外線カメラにおける第二の状態の被検体に対するＭＲＩ画像上の第一の状態の被検体の位置合わせ）についてのみ説明する。
まず、赤外線カメラに対するＭＲＩ画像上の被検体に対して位置・姿勢のみ様々に変動させた位置合わせパラメータの候補値を仮定し、仮定した位置・姿勢ごとにMRI画像内の体表近傍の情報を利用することで表在血管を強調したＭＩＰ像を生成する。
次に、生成した各ＭＩＰ像と赤外線カメラ画像の類似性を評価して、評価値が最大になる赤外線カメラに対するＭＲＩ画像上の被検体の位置・姿勢を推定する。これにより、被検体が剛体と見做せる場合は、超音波画像と位置合わせする対象である三次元画像の間で、高精度に剛体位置合わせを行うことができる。
（変形例２）
第３の実施形態では、赤外線カメラ画像を利用することで、ＭＲＩ画像などの被検体の三次元画像と、位置・姿勢センサで計測された超音波探触子によって取得された二次元の超音波画像との位置合わせを行う例を説明した。
しかし、被検体の三次元画像と位置合わせを行う対象は、例えば、超音波画像と同様に、位置・姿勢センサで計測された超音波探触子で取得された二次元のPAT画像でもよい。二次元のＰＡＴ画像は、超音波探触子に近赤外光の光源を備えることにより、人体に対して超音波ではなく近赤外光を照射するという点を除いて、超音波画像と同様の構成で画像を取得する。従って、第３の実施形態と同様に、赤外線カメラに対するＭＲＩ画像上の被検体の位置・姿勢および重力変形を推定する。そして、位置・姿勢センサにより分かっている二次元のＰＡＴ画像と赤外線カメラの位置関係を利用して、二次元のＰＡＴ画像とＭＲＩ画像の高精度な変形位置合わせを行うことができる。
（第４の実施形態）
以下、本発明にかかる第４の実施形態の画像処理を説明する。なお、第４の実施形態において、第１乃至第３の実施形態と略同様の構成については、同一符号を付して、その詳細説明を省略する。
第１乃至第３の実施形態においては、第一の処理として、赤外線カメラに対する、予め撮像されたＭＲＩ画像などの被検体の三次元画像の位置・姿勢および変形状態を推定した。そして、第二の処理として、目的とするモダリティと赤外線カメラとの位置関係を利用して、三次元画像とモダリティの位置合わせを行った。
しかし、実施する処理は第一の処理のみでもよい。つまり、赤外線カメラ画像と予め撮像された第二の状態の被検体の三次元画像の両画像に写った表在血管の情報を利用して、赤外線カメラに写った被検体に対する、予め撮像された第一の状態の被検体の三次元画像の位置・姿勢と変形を推定する処理のみでもよい。
上記の位置・姿勢および変形の推定により、例えば、被検体の手術時において、被検体を写す赤外線カメラに対して、予め撮像した被検体のＭＲＩ画像やＣＴ画像などの三次元画像を位置合わせすることができる。これにより、三次元画像に写った病変部位などを、手術中の被検体と対応付けて参照することが可能になる。
このように、赤外線カメラで撮像された二次元画像と、ＭＲＩ画像やＣＴ画像などの三次元画像から取得される被検体の表面近傍の内部構造の情報を利用して、三次元画像を被検体に対して高精度かつ自動的に位置合わせする仕組みを提供することができる。
（第５の実施形態）
以下、本発明にかかる第５の実施形態の画像処理を説明する。なお、第５の実施形態において、第１乃至第３の実施形態と略同様の構成については、同一符号を付して、その詳細説明を省略する。
第１乃至第３の実施形態においては、目的とするモダリティに位置が対応付けられた赤外線カメラ（一般的なカメラでも良い）に対する、予め撮像されたＭＲＩ画像などの被検体の三次元画像の位置・姿勢や変形状態を推定することで、三次元画像と目的とするモダリティの位置合わせを行った。
しかし、目的とするモダリティ（本実施形態ではＰＡＴ画像とする）と位置が対応付けられた撮像部は、被検体の距離画像を取得する距離画像カメラであっても良い。これにより、第１乃至第３の実施形態ではカメラ画像から取得できるのが被検体の２次元形状であったのに対し、被検体の距離画像を取得することによって、被検体の３次元形状を取得することができる。そして、取得した被検体の３次元形状と予め撮像された三次元画像（本実施形態ではＭＲＩ画像とする）に写る被検体の３次元形状を比較することで、距離画像カメラに対するＭＲＩ画像の位置・姿勢や変形状態を推定することができる。最後に、距離画像カメラとＰＡＴ装置の位置の対応関係を利用して、ＭＲＩ画像とＰＡＴ画像の位置合わせを行うことができる。以下、第５の実施形態の機能構成及び処理の流れを説明する。
ここで、第５の実施形態における機能構成及び処理の流れは、第１の実施形態における図１の機能構成及び図６の処理の流れの赤外線カメラに関する部分を、距離画像カメラに置き換えたものとなる。従って、共通する機能構成及び処理の流れは説明を省略し、差分のみを説明する。まず、機能構成について説明する。
距離画像カメラは、該カメラが撮像する画像に撮像されるシーン中の各点までの距離を画素毎に（各画素に写っているシーン中の部分毎に）計測し、該距離を表す情報を各画素の画素値として有する距離画像を生成するカメラである。距離画像カメラによって対象物体を計測することで、対象物体の表面の三次元形状を距離画像として撮像できる。なお、本実施形態で用いる距離画像カメラは、タイムオブフライト方式やパターン投光方式など、公知のいずれの方式に基づくものであっても良い。
本実施形態では、カメラ画像取得部１０４を、距離画像を取得する機能を備えた構成要素（距離画像取得部と定義する）に置き換える。
本実施形態では、二次元形状取得部１０５を、距離画像から被検体の三次元形状を取得する機能を備えた構成要素（形状取得部と定義する）に置き換える。
本実施形態では、仮想投影像生成部１０７と仮想投影像評価部１０８を、ＭＲＩ画像から取得した被検体の三次元形状と、非保持状態の被検体を撮像した距離画像から取得した被検体の三次元形状とを比較する機能を備えた構成要素（第一の評価部と定義する）に置き換える。
本実施形態では、変形画像評価部１１１を、変形ＭＲＩ画像から取得した被検体の三次元形状と、保持状態の被検体を撮像した距離画像から取得した被検体の三次元形状とを比較する機能（第二の評価部）に置き換える。
次に、処理の流れについて説明する。
本実施形態では、ステップＳ６０５を、距離画像取得部が非保持状態及び保持状態の被検体を撮像した距離画像を取得する処理に置き換える。
本実施形態では、ステップＳ６０６を、形状取得部が、ステップＳ６０５で取得した非保持状態の乳房の距離画像から被検体の乳房の三次元形状と乳頭位置を取得する処理に置き換える。被検体の三次元形状の取得については、例えば、被検体の写った距離画像と被検体の写っていない距離画像の差分画像を生成して、差分画像中の値が所定値より小さい領域は、被検体以外の領域（背景領域）であると判定して削除することで、被検体の領域における距離データ、つまり被検体の三次元形状を取得することができる。なお、乳頭位置は、ステップＳ６０２と同様、三次元形状の曲率を計算し、特徴点を検出することで取得することができる。
本実施形態では、ステップＳ６０７を、剛体変換部１０６及び第一の評価部が、ステップＳ６０２及びＳ６０３で取得したＭＲＩ画像に写る乳房の三次元形状及び乳頭位置と、ステップＳ６０６で取得した距離画像に写る非保持状態の乳房の三次元形状及び乳頭位置とを比較し、両者の三次元形状が最も整合するような剛体変換を推定する処理に置き換える。より具体的には、剛体変換部１０６は、ステップＳ６０３で取得したＭＲＩ画像上の三次元の乳頭位置がステップＳ６０６で取得した距離画像上の三次元の乳頭位置に一致するように、平行移動を表す座標変換行列を算出する。そして、第一の評価部は、ＭＲＩ画像上の三次元形状に対して平行移動を表す変換行列を適用することで、乳頭位置を一致させた上で、ＭＲＩ画像と距離画像の三次元形状の残差を算出する。そして、剛体変換部１０６は、残差が最小になるような回転移動の座標変換行列を探索して算出し、平行移動と回転移動を合わせた剛体変換の座標変換行列を算出する。
本実施形態では、ステップＳ６０９を、変形推定部１０９、変形画像生成部１１０、及び第二の評価部が、ＭＲＩ画像を変形させた変形ＭＲＩ画像から取得した乳房の三次元形状と、ステップＳ６０６で取得した距離画像に写る保持状態の乳房の三次元形状とを比較し、両者の三次元形状が最も整合するような変形を推定する処理に置き換える。より具体的には、変形推定部１０９は、第１の実施形態におけるステップＳ６０９と同様の方法で、ＭＲＩ画像を変形させた変形ＭＲＩ画像を生成する。そして、第二の評価部は、変形ＭＲＩ画像からステップＳ６０２と同様の方法で取得した乳房の三次元形状と、距離画像の保持状態の乳房の三次元形状と残差を算出する。そして、変形推定部１０９は、残差が最小になるような変形パラメータを探索して算出する。なお、ここでは、三次元形状の残差に基づいて変形パラメータを推定したが、無論、第１の実施形態と同様、残差に加えて、変形ＭＲＩ画像とＰＡＴ画像との間の類似度を加味した評価値に基づいて変形パラメータを推定しても良い。
このように、距離画像から取得した乳房の３次元形状を利用して、距離画像に対するＭＲＩ画像の位置・姿勢や変形状態を推定することで、カメラ画像から得られる２次元形状を利用するよりも、高精度な推定を行うことができる。
（その他の実施形態）
本発明は、以下の処理を実行することによっても実現される。即ち、上述した実施形態の機能を実現するソフトウェア（プログラム）を、ネットワーク又は各種記憶媒体を介してシステム或いは装置に供給し、そのシステム或いは装置のコンピュータ（またはＣＰＵやＭＰＵ等）がプログラムを読み出して実行する処理である。

Claims (16)

  1. 第一の撮像装置により対象物を撮像した第一の画像と、第二の撮像装置により前記対象物を撮像した第二の画像と、前記第二の撮像装置に対して位置が対応付けられた撮像手段により前記対象物を撮像した第三の画像と、を取得する画像取得手段と、
    前記第一の画像における対象物の観測情報と前記第三の画像における対象物の観測情報とが整合するように前記第一の画像における対象物と前記第二の画像における対象物との位置合わせをする位置合わせ手段と、
    を備えることを特徴とする画像処理装置。
  2. 前記位置合わせ手段は、前記第一の画像における対象物の観測情報と前記第三の画像における対象物の観測情報とが整合するように位置合わせを行った後に、該位置合わせ結果を利用して、前記第一の画像における対象物の観測情報と前記第二の画像における対象物の観測情報との比較を行って、前記第一の画像における対象物と前記第二の画像における対象物との位置合わせを行うことを特徴とする請求項１に記載の画像処理装置。
  3. 前記画像取得手段は、
    第一の状態の前記対象物を前記第一の撮像装置により撮像した前記第一の画像を取得する第一の画像取得手段と、
    前記第一の状態とは異なる第二の状態の前記対象物を前記第二の撮像装置により撮像した前記第二の画像を取得する第二の画像取得手段と、
    前記第二の状態の初期段階の前記対象物を前記撮像手段により撮像した前記第三の画像を取得する第三の画像取得手段と、
    を備え、
    前記位置合わせ手段は、前記第一の画像と前記第三の画像との比較により、前記第一の画像と前記第三の画像との位置合わせを行った後に、該位置合わせ結果を利用して、前記第一の画像と前記第二の画像との比較により前記第一の画像と前記第二の画像との位置合わせを行うことを特徴とする請求項２に記載の画像処理装置。
  4. 前記第一の撮像装置及び前記第二の撮像装置は、前記対象物の三次元画像を前記第一の画像及び前記第二の画像として撮像する撮像装置であり、
    前記撮像手段はカメラであり、
    前記位置合わせ手段は、
    前記第一の画像と前記第三の画像との間における前記対象物の表面形状情報の比較と、前記第一の画像と前記第三の画像との間における前記対象物の内部情報の比較とにより、前記第一の画像と前記第二の画像との位置合わせを行うことを特徴とする請求項２または３に記載の画像処理装置。
  5. 前記撮像手段は、前記対象物の表面近傍の血管領域を可視化する赤外線カメラ画像を撮像する赤外線カメラであり、
    前記位置合わせ手段は、前記第一の画像と前記赤外線カメラ画像との間で、前記対象物の表面近傍の血管領域を比較して前記位置合わせを行うことを特徴とする請求項１乃至４のいずれか１項に記載の画像処理装置。
  6. 前記撮像手段は、前記対象物の表面の三次元形状を距離画像として撮像する距離画像カメラであり、
    前記位置合わせ手段は、前記第一の画像と前記距離画像との間で、前記対象物の三次元の表面形状を比較して前記位置合わせを行うことを特徴とする請求項１乃至４のいずれか１項に記載の画像処理装置。
  7. 前記第一の撮像装置は、前記対象物の三次元画像を前記第一の画像として撮像する撮像装置であり、
    前記三次元画像から前記対象物の表面位置を示す情報を取得する情報取得手段と、
    前記表面位置を示す情報に基づいて、前記撮像手段の視点から見た前記三次元画像の投影像を生成する投影像生成手段と、
    を更に備え、
    前記位置合わせ手段は、前記投影像と前記第三の画像における対象物の観測情報とが整合するように前記第一の画像における対象物と前記第二の画像における対象物との位置合わせを行う
    ことを特徴とする請求項１に記載の画像処理装置。
  8. 前記第一の画像と前記第三の画像との位置合わせの結果を用いて、前記第一の状態で撮像された前記第一の画像を前記第二の状態の第一の画像に変形する変形画像生成手段を更に備え、
    前記変形画像生成手段は、前記第一の画像と前記第二の画像との位置合わせを行うための位置合わせパラメータを用いて、前記第一の状態から前記第二の状態に変形した第一の画像を生成し、
    前記投影像生成手段は、前記変形した第一の画像を用いて、前記撮像手段の視点から見た投影像を生成する
    ことを特徴とする請求項７に記載の画像処理装置。
  9. 前記位置合わせパラメータには、前記対象物の位置、姿勢、および前記対象物の変形を表すパラメータが含まれることを特徴とする請求項８に記載の画像処理装置。
  10. 前記変形画像生成手段は、前記位置合わせパラメータの複数の候補値に基づき、前記第二の状態に変形した複数の第一の画像を生成することを特徴とする請求項８に記載の画像処理装置。
  11. 前記第二の状態に変形した複数の第一の画像と、前記第二の画像との類似度を取得する評価手段を更に備え、
    前記位置合わせ手段は、前記類似度を用いた前記第一の画像と前記第二の画像との比較により、前記位置合わせを行うことを特徴とする請求項１０に記載の画像処理装置。
  12. 前記評価手段は、前記類似度の評価により前記複数の候補値から位置合わせパラメータを一つ選択することを特徴とする請求項１１に記載の画像処理装置。
  13. 前記第一の画像と前記第二の画像とを表示部に並べて表示する画像表示手段を更に備えることをと特徴とする請求項１乃至１２のいずれか１項に記載の画像処理装置。
  14. 画像処理装置の画像処理方法であって、
    画像取得手段が、第一の撮像装置により対象物を撮像した第一の画像と、第二の撮像装置により前記対象物を撮像した第二の画像と、前記第二の撮像装置に対して位置が対応付けられた撮像手段により前記対象物を撮像した第三の画像と、を取得する画像取得工程と、
    位置合わせ手段が、前記第一の画像における対象物の観測情報と前記第三の画像における対象物の観測情報とが整合するように前記第一の画像における対象物と前記第二の画像における対象物との位置合わせをする位置合わせ工程と、
    を有することを特徴とする画像処理方法。
  15. コンピュータを、請求項１乃至１３のいずれか１項に記載の画像処理装置の各手段として機能させるためのプログラム。
  16. 請求項１５に記載されたプログラムが記憶されたコンピュータにより読み取り可能な記憶媒体。
JP2014022779A 2014-02-07 2014-02-07 画像処理装置、画像処理方法、プログラムおよび記憶媒体 Active JP6289142B2 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
JP2014022779A JP6289142B2 (ja) 2014-02-07 2014-02-07 画像処理装置、画像処理方法、プログラムおよび記憶媒体
US14/613,464 US9691150B2 (en) 2014-02-07 2015-02-04 Image processing apparatus, image processing method and storage medium

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2014022779A JP6289142B2 (ja) 2014-02-07 2014-02-07 画像処理装置、画像処理方法、プログラムおよび記憶媒体

Publications (3)

Publication Number Publication Date
JP2015146989A true JP2015146989A (ja) 2015-08-20
JP2015146989A5 JP2015146989A5 (ja) 2017-03-16
JP6289142B2 JP6289142B2 (ja) 2018-03-07

Family

ID=53775367

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2014022779A Active JP6289142B2 (ja) 2014-02-07 2014-02-07 画像処理装置、画像処理方法、プログラムおよび記憶媒体

Country Status (2)

Country Link
US (1) US9691150B2 (ja)
JP (1) JP6289142B2 (ja)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2017047103A (ja) * 2015-09-04 2017-03-09 キヤノン株式会社 制御装置、乳房撮影システム、乳房撮影装置の制御方法
JP2017047079A (ja) * 2015-09-04 2017-03-09 東芝メディカルシステムズ株式会社 医用画像表示装置及びマンモグラフィ装置
JP2018011637A (ja) * 2016-07-19 2018-01-25 キヤノン株式会社 画像処理装置および画像処理方法
JP2018015262A (ja) * 2016-07-28 2018-02-01 国立大学法人京都大学 光音響イメージング装置
JP2018117956A (ja) * 2017-01-27 2018-08-02 キヤノン株式会社 音響波測定装置およびその制御方法
WO2018230409A1 (ja) * 2017-06-15 2018-12-20 キヤノン株式会社 情報処理装置、情報処理方法、及びプログラム
JP2019506949A (ja) * 2016-02-29 2019-03-14 コーニンクレッカ フィリップス エヌ ヴェＫｏｎｉｎｋｌｉｊｋｅ Ｐｈｉｌｉｐｓ Ｎ．Ｖ． 医療胸部画像の修正ための装置、撮像システム及び方法。
JP2019042304A (ja) * 2017-09-05 2019-03-22 株式会社ニデック 眼科用画像処理プログラム
JP2019165338A (ja) * 2018-03-19 2019-09-26 株式会社リコー 画像処理装置及び投影システム

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2676798C2 (ru) 2013-09-18 2019-01-11 Конинклейке Филипс Н.В. Протравленные лазером сцинтиллирующие кристаллы для повышенной эффективности
US10078889B2 (en) * 2015-08-25 2018-09-18 Shanghai United Imaging Healthcare Co., Ltd. System and method for image calibration
US10492694B2 (en) * 2015-08-27 2019-12-03 Canon Kabushiki Kaisha Object information acquisition apparatus
US10357216B2 (en) * 2015-12-21 2019-07-23 Toshiba Medical Systems Corporation Mammography apparatus
US10274935B2 (en) * 2016-01-15 2019-04-30 Honeywell Federal Manufacturing & Technologies, Llc System, method, and computer program for creating geometry-compliant lattice structures
JP2017164312A (ja) * 2016-03-16 2017-09-21 キヤノン株式会社 放射線撮影装置、挿入状態判定方法、及びプログラム
JP2018126389A (ja) * 2017-02-09 2018-08-16 キヤノン株式会社 情報処理装置、情報処理方法、およびプログラム
CN111339822B (zh) * 2017-07-17 2023-06-30 Oppo广东移动通信有限公司 活体检测方法及相关产品
JP7022584B2 (ja) * 2017-12-27 2022-02-18 キヤノン株式会社 放射線撮影装置、画像処理装置及び画像判定方法
JP2019191989A (ja) * 2018-04-26 2019-10-31 キヤノン株式会社 仮想視点画像を生成するシステム、方法及びプログラム
WO2020101019A1 (en) * 2018-11-16 2020-05-22 Ricoh Company, Ltd. Information processing system, information processing apparatus, and recording medium
CN110335271B (zh) * 2019-07-10 2021-05-25 浙江铁素体智能科技有限公司 一种电气部件故障的红外检测方法及装置

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH095440A (ja) * 1995-06-09 1997-01-10 Adac Lab Inc 非一様減衰補正を含む二重ｓｐｅｃｔ・ｐｅｔ撮影用のマルチヘッド核医学カメラ
JP2006320387A (ja) * 2005-05-17 2006-11-30 Univ Of Tsukuba 計算機支援診断装置および方法
US20100239150A1 (en) * 2008-12-05 2010-09-23 Canon Kabushiki Kaisha Information processing apparatus for registrating medical images, information processing method and program
US20100245766A1 (en) * 2009-03-17 2010-09-30 Zhang Hao F Systems and methods for photoacoustic opthalmoscopy
US20140350358A1 (en) * 2011-12-26 2014-11-27 Canon Kabushiki Kaisha Subject information accumulating apparatus

Family Cites Families (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7778686B2 (en) * 2002-06-04 2010-08-17 General Electric Company Method and apparatus for medical intervention procedure planning and location and navigation of an intervention tool
US7689003B2 (en) * 2006-03-20 2010-03-30 Siemens Energy, Inc. Combined 2D and 3D nondestructive examination
US7996060B2 (en) * 2006-10-09 2011-08-09 Biosense Webster, Inc. Apparatus, method, and computer software product for registration of images of an organ using anatomical features outside the organ
US8194936B2 (en) * 2008-04-25 2012-06-05 University Of Iowa Research Foundation Optimal registration of multiple deformed images using a physical model of the imaging distortion
JP2010088627A (ja) 2008-10-07 2010-04-22 Canon Inc 生体情報処理装置および生体情報処理方法
JP5147656B2 (ja) * 2008-11-20 2013-02-20 キヤノン株式会社 画像処理装置、画像処理方法、プログラム、及び記憶媒体
US8675996B2 (en) * 2009-07-29 2014-03-18 Siemens Aktiengesellschaft Catheter RF ablation using segmentation-based 2D-3D registration
JP5737858B2 (ja) * 2010-04-21 2015-06-17 キヤノン株式会社 画像処理装置、画像処理方法、及びプログラム
JP5685133B2 (ja) * 2011-04-13 2015-03-18 キヤノン株式会社 画像処理装置、画像処理装置の制御方法、およびプログラム
JP5822554B2 (ja) * 2011-06-17 2015-11-24 キヤノン株式会社 画像処理装置、画像処理方法、撮影システム及びプログラム
US8805038B2 (en) * 2011-06-30 2014-08-12 National Taiwan University Longitudinal image registration algorithm for infrared images for chemotherapy response monitoring and early detection of breast cancers
KR101572487B1 (ko) * 2013-08-13 2015-12-02 한국과학기술연구원 환자와 3차원 의료영상의 비침습 정합 시스템 및 방법
JP6182045B2 (ja) * 2013-10-11 2017-08-16 キヤノン株式会社 画像処理装置およびその方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH095440A (ja) * 1995-06-09 1997-01-10 Adac Lab Inc 非一様減衰補正を含む二重ｓｐｅｃｔ・ｐｅｔ撮影用のマルチヘッド核医学カメラ
JP2006320387A (ja) * 2005-05-17 2006-11-30 Univ Of Tsukuba 計算機支援診断装置および方法
US20100239150A1 (en) * 2008-12-05 2010-09-23 Canon Kabushiki Kaisha Information processing apparatus for registrating medical images, information processing method and program
US20100245766A1 (en) * 2009-03-17 2010-09-30 Zhang Hao F Systems and methods for photoacoustic opthalmoscopy
US20140350358A1 (en) * 2011-12-26 2014-11-27 Canon Kabushiki Kaisha Subject information accumulating apparatus

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2017047103A (ja) * 2015-09-04 2017-03-09 キヤノン株式会社 制御装置、乳房撮影システム、乳房撮影装置の制御方法
JP2017047079A (ja) * 2015-09-04 2017-03-09 東芝メディカルシステムズ株式会社 医用画像表示装置及びマンモグラフィ装置
JP2019506949A (ja) * 2016-02-29 2019-03-14 コーニンクレッカ フィリップス エヌ ヴェＫｏｎｉｎｋｌｉｊｋｅ Ｐｈｉｌｉｐｓ Ｎ．Ｖ． 医療胸部画像の修正ための装置、撮像システム及び方法。
JP2018011637A (ja) * 2016-07-19 2018-01-25 キヤノン株式会社 画像処理装置および画像処理方法
JP2018015262A (ja) * 2016-07-28 2018-02-01 国立大学法人京都大学 光音響イメージング装置
JP2018117956A (ja) * 2017-01-27 2018-08-02 キヤノン株式会社 音響波測定装置およびその制御方法
WO2018230409A1 (ja) * 2017-06-15 2018-12-20 キヤノン株式会社 情報処理装置、情報処理方法、及びプログラム
JP2019042304A (ja) * 2017-09-05 2019-03-22 株式会社ニデック 眼科用画像処理プログラム
JP2019165338A (ja) * 2018-03-19 2019-09-26 株式会社リコー 画像処理装置及び投影システム
JP7139637B2 (ja) 2018-03-19 2022-09-21 株式会社リコー 画像処理装置及び投影システム

Also Published As

Publication number Publication date
US20150228093A1 (en) 2015-08-13
JP6289142B2 (ja) 2018-03-07
US9691150B2 (en) 2017-06-27

Similar Documents

Publication Publication Date Title
JP6289142B2 (ja) 画像処理装置、画像処理方法、プログラムおよび記憶媒体
JP6182045B2 (ja) 画像処理装置およびその方法
CN111161326B (zh) 用于可变形图像配准的无监督深度学习的系统和方法
US9396576B2 (en) Method and apparatus for estimating the three-dimensional shape of an object
CN105025803B (zh) 从多个三维视图对大对象的分割
KR102114415B1 (ko) 의료 영상 정합 방법 및 장치
KR101267759B1 (ko) 화상 처리 장치, 화상 처리 방법 및 저장 매체
JP5486182B2 (ja) 情報処理装置および情報処理方法
CN103402453A (zh) 用于导航系统的自动初始化和配准的系统和方法
KR102273020B1 (ko) 의료 영상 정합 방법 및 그 장치
CN105899144B (zh) 图像处理装置、图像诊断系统、图像处理方法和存储介质
KR102233966B1 (ko) 의료 영상 정합 방법 및 그 장치
EP2854650A1 (en) Image generation apparatus
KR102250086B1 (ko) 의료 영상 정합 방법, 이를 포함하는 장치 및 컴퓨터 기록 매체
Ben-Hamadou et al. Construction of extended 3D field of views of the internal bladder wall surface: A proof of concept
JP6461257B2 (ja) 画像処理装置およびその方法
KR20160057024A (ko) 마커리스 ３차원 객체추적 장치 및 그 방법
JP2014150855A (ja) 乳房診断支援システム、および乳房データ処理方法
Nitkunanantharajah et al. Trackerless panoramic optoacoustic imaging: a first feasibility evaluation
KR102007316B1 (ko) 외장 의안의 3차원 영상 생성 장치 및 방법
Sun Ultrasound probe localization by tracking skin features
JP5808446B2 (ja) 情報処理装置および情報処理方法
KR102336446B1 (ko) 의료 영상 정합 방법 및 그 장치
Chen QUiLT (Quantitative Ultrasound in Longitudinal Tissue Tracking): Stitching 2D images into 3D Volumes for Organ Health Monitoring
Carminati et al. Automated registration of 3D TEE datasets of the descending aorta for improved examination and quantification of atheromas burden

Legal Events

Date Code Title Description
A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20170206

A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20170206

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20171020

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20171018

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20171218

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: 20180109

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20180206

R151 Written notification of patent or utility model registration

Ref document number: 6289142

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R151