JP6436442B2 - 光音響装置および画像処理方法 - Google Patents

Description

本発明は、光音響装置および画像処理方法に関する。

被検体の内部を可視化するイメージング技術の一つに、光音響イメージングがある。光音響イメージングでは、まず、光源から発生したパルス光が被検体に照射される。照射光が被検体内で伝播・拡散し、血液中のヘモグロビンなどの光吸収体に吸収されると、光音響効果により音響波（光音響波）が発生する。被検体の周囲に配置された複数の受信素子がこの光音響波を受信して受信信号（光音響信号）を出力する。プロセッサーが光音響信号を解析して画像再構成することで、被検体内部の光学特性値を示す画像データ（光音響画像）が得られる。

光音響信号には、様々な要因により、ＳＮ比低下の原因となるノイズが含まれていることがある。ＳＮ比の低い信号を用いて画像再構成を行うと光学特性情報の定量性が低くなる。特許文献１には、複数回に渡ってパルス光の照射と光音響信号の計測を行い、取得した複数の光音響信号を加算平均することでノイズを低減する方法が開示されている。さらに、特許文献１には、光音響画像と位置が対応付けられた超音波画像を取得し、その超音波画像をトラッキングして被検体の体動を推定して補正する方法が開示されている。さらに体動補正後の光音響信号を加算平均して、ノイズ低減効果を大きくしている。

また、被検体にパルス光を照射したときに、光音響信号から高い品質で画像化できる範囲には限りがある。この範囲はＦＯＶ（Ｆｉｅｌｄ ｏｆ Ｖｉｅｗ）と呼ばれる。従って、被検体の撮像対象領域がＦＯＶよりも広い場合には、パルス光の照射位置を移動させてＦＯＶを少しずつずらしながら光音響信号を取得する。そして移動量を考慮して夫々の光音響信号を加算平均することで、広い撮像対象領域を高精度に画像化できる。また、このように照射位置を移動させながら撮影したデータに対して、特許文献１のような体動推定を行う場合は、被検体の隣接ＦＯＶごとに、時系列順に位置ズレ補正を行う。

米国特許出願公開第２０１０／００４９０４４号明細書

しかしながら、上記の方法では、局所的な補正の誤差の積み重ねに起因して、加算平均後の光音響画像が本来の被検体の解剖構造に対して歪む可能性がある。

本発明は上記課題に鑑みてなされたものである。本発明の目的は、光音響イメージングにおいて被検体が動いた場合などにおいても、品質の高い光音響画像を生成する方法を提供することにある。

本発明は、以下の構成を採用する。すなわち、
光源と、
前記光源から光を照射された被検体から発生する音響波を受信する受信素子と、
前記音響波を用いて、前記被検体内部の特性情報を示す画像データを生成する処理手段と、
前記被検体における前記光の照射位置を変化させる変化手段と、
前記被検体の広域画像を取得する広域画像取得手段と、
を有し、
前記処理手段は、
前記変化手段により変化する前記照射位置ごとに、当該照射位置に対応する前記被検体の局所領域における前記画像データである局所画像を生成し、
前記照射位置ごとに得られた複数の前記局所画像の間の比較である第１の比較と、前記複数の局所画像のそれぞれと前記広域画像との比較である第２の比較と、に基づいて前記複数の局所画像を統合する
ことを特徴とする光音響装置である。

本発明はまた、以下の構成を採用する。すなわち、
照射位置を変化させながら光を照射された被検体から発生する音響波を用いて、前記被検体内部の特性情報を示す画像データを生成する画像処理方法であって、
前記被検体の広域画像を取得するステップと、
変化する前記照射位置ごとに、当該照射位置に対応する前記被検体の局所領域における前記画像データである局所画像を生成するステップと、
前記照射位置ごとに得られた複数の前記局所画像の間の比較と、前記複数の局所画像のそれぞれと前記広域画像との比較と、に基づいて前記複数の局所画像を統合するステップと、
を有することを特徴とする画像処理方法である。

本発明によれば、光音響イメージングにおいて被検体が動いた場合などにおいても、品質の高い光音響画像を生成する方法を提供できる。

第１の実施形態に係る画像処理装置の機能構成を示す図。 光音響信号計測装置が被検体を撮像可能な範囲を示す図。 赤外線カメラによって被検体が撮像される様子を示す図。 赤外線カメラによって撮影される赤外線カメラ画像の例を示す図。 第１の実施形態における全体の処理手順を示すフロー図。 比較ペア情報の生成方法を示す図。 赤外線カメラと投影画像とパルスボリュームの空間的関係を示す図。 パルス位置推定部によって評価される画像の関係の例を示す図。 第２の実施形態に係る画像処理装置の機能構成を示す図。 第２の実施形態における全体の処理手順を示すフロー図。 第３の実施形態に係る画像処理装置の機能構成を示す図。 第３の実施形態における全体の処理手順を示すフロー図。 第４の実施形態に係る画像処理装置の機能構成を示す図。 第４の実施形態における全体の処理手順を示すフロー図。 光音響信号計測装置の構成例を示す図。

以下に図面を参照しつつ、本発明の好適な実施の形態について説明する。ただし、以下に記載されている構成部品の寸法、材質、形状およびそれらの相対配置などは、発明が適用される装置の構成や各種条件により適宜変更されるべきものである。よって、この発明の範囲を以下の記載に限定する趣旨のものではない。

本発明は、被検体から伝播する音響波を検出し、被検体内部の特性情報を生成し、取得
する技術に関する。よって本発明は、被検体情報取得装置またはその制御方法、あるいは被検体情報取得方法や信号処理方法として捉えられる。本発明はまた、これらの方法をＣＰＵ等のハードウェア資源を備える情報処理装置に実行させるプログラムや、そのプログラムを格納した記憶媒体としても捉えられる。本発明はまた、被検体内部の特性情報を補正する装置や、その制御方法、特性情報の補正方法、補正プログラムとも捉えられる。本発明はまた、被検体内部の特性情報を示す画像データを処理する機能を有することから、画像処理装置、画像処理方法および画像処理プログラムとも捉えられる。

本発明の被検体情報取得装置は、被検体に光（電磁波）を照射し、光音響効果に従って被検体内または被検体表面の特定位置で発生して伝播した音響波を受信（検出）する、光音響トモグラフィ技術を利用した装置を含む。このような被検体情報取得装置は、光音響測定に基づき被検体内部の特性情報を画像データ等の形式で得ることから、光音響装置、または、光音響イメージング装置とも呼べる。

光音響装置における特性情報とは、光照射によって生じた音響波の発生源分布、被検体内の初期音圧分布、あるいは初期音圧分布から導かれる光エネルギー吸収密度分布や吸収係数分布、組織を構成する物質の濃度分布を示す。具体的には、酸化・還元ヘモグロビン濃度分布や、それらから求められる酸素飽和度分布などの血液成分分布、あるいは脂肪、コラーゲン、水分の分布などである。また、特性情報は、数値データとしてではなく、被検体内の各位置の分布情報として求めてもよい。すなわち、吸収係数分布や酸素飽和度分布などの分布情報を被検体情報としてもよい。これらの分布情報を画像化したものは、光音響画像とも呼べる。光音響画像は、３次元空間中の分布情報として、典型的には、３次元に配列したボクセルの夫々に当該位置における値を収めたボリュームデータの形式で取得される。ボリュームデータは、３次元ボリュームや３次元画像、３次元断層像とも呼べる。

本発明でいう音響波とは、典型的には超音波であり、音波、音響波と呼ばれる弾性波を含む。光音響効果により発生した音響波のことを、光音響波または光超音波と呼ぶ。探触子等により音響波から変換された電気信号を音響信号とも呼ぶ。光音響波に由来する音響信号を、光音響信号とも呼ぶ。

＜第１の実施形態＞
本実施形態に係る光音響装置は、複数回に渡って計測した一連の光音響信号を統合した画像を生成する。光音響装置は、計測中に生じる被検体の体動等の動きを推定して、その動きの影響を低減するように信号を補正する。より具体的には、光音響装置は、ＦＯＶが異なる複数の光音響画像（局所画像）同士を比較すると同時に、各光音響画像（局所画像）を、被検体全体を撮像した赤外線カメラ画像（広域画像）と比較する。これらの比較結果を用いた動き推定により、本来の被検体の解剖構造に対して全体として歪みが抑制された音響画像を取得できる。

（装置構成）
図１は、本実施形態に係る光音響装置の構成を示す。同図に示すように、本実施形態における光音響装置は、画像処理装置１０００、光音響信号計測装置１１０、赤外線カメラ１２０、および表示装置１３０を有する。画像処理装置１０００は、光音響信号計測装置１１０、赤外線カメラ１２０、および表示装置１３０と接続されている。

光音響信号計測装置１１０は、パルス光を被検体に照射し、光音響効果により発生した光音響波を被検体周辺の複数の受信素子で受信して電気信号（光音響信号）に変換する。光音響信号計測装置１１０は、パルス光照射と光音響信号の計測を複数回に渡って繰り返し行う。ここで、１回のパルス光照射（または１つのＦＯＶ）に対応する光音響信号を１
組とすると、パルス光の照射回数分だけの光音響信号の組が取得できる。計測した光音響信号は、画像処理装置１０００へと入力される。

（ＦＯＶ）
図２は、光音響信号計測装置が被検体を撮像可能な範囲（ＦＯＶ）を示す。ＦＯＶ２０２は、被検体２００に対して照射方向２０１の方向にパルス光を照射した時に撮像可能な範囲である。本実施形態の照射方向２０１は、光音響信号計測装置１１０の装置座標系C_DEVのＺ軸に平行である。図２のように、１つのＦＯＶでは被検体２００全体をカバーできない場合、光照射ごとにＦＯＶを少しずつずらしながら被検体全体から光音響信号を取得する必要がある。そのためには、光の射出端を可動式にする、光照射方向を変化可能にするなどの方法がある。このような構成によって、関心領域が１つのＦＯＶよりも広い場合に対応できる。なお、後の画像統合で利用するために、パルス照射ごとにＦＯＶの位置を記録することが好ましい。位置記録方法として例えば、光照射位置などを示す基準座標を記録する方法がある。

（光音響信号計測装置）
図１５に、光音響信号計測装置１１０の構成例を示す。測定対象となる被検体１５００は、例えば生体の乳房である。光音響信号計測装置１１０は、画像処理装置１０００と通信して、制御情報の受信や取得信号の送信を行う。光音響信号計測装置は、本発明の光音響装置に相当すると考えることもできるし、光音響信号計測装置と画像処理装置を合わせて光音響装置と考えてもよい。

図１５（ａ）は、被検体１５００を２枚の板状の保持部材１５０６で挟持するタイプの装置である。光源１５０２から発生したパルス光１５０５は、光学系１５０３により伝送され、射出端１５０４から被検体に照射される。探触子１５０１に配置された複数の受信素子１５０９（ａ〜ｎ）は、被検体から発生した光音響波を検出してアナログ電気信号である光音響信号に変換し、画像処理装置１０００に送信する。受信素子は、光音響信号を計測する信号計測部とも呼べる。

一方、図１５（ｂ）は、被検体１５００を筐体状の保持部材１５０６の開口部から下に垂らすタイプの装置である。本図では、お椀状の探触子１５０１に配置された複数の受信素子が光音響波を計測する。この形態では、受信素子の高感度方向が集中する高分解能領域が形成される。この場合は探触子および受信素子が信号計測部に当たる。被検体を支持して形状を安定させるために、光および音響波を透過するカップ状の支持部材を設けることも好ましい。

図１では画像処理装置１０００を機能ブロックの観点から説明したが、ここで図１５を用いて物理的構成の観点から説明する。画像処理装置１０００は例えば、ＣＰＵ１５０７ａ、ＧＰＵ１５０７ｂ、メモリ１５０７ｃ、ＦＰＧＡ１５０７ｄを備えて構成される。これらの各ブロックが情報処理や信号処理を行うことで、画像処理方法が実現される。また、オンライン上の演算資源を利用して情報処理を行っても良い。

光源１５０２としては、パルスレーザ光を発生可能なレーザ光源が好ましい。またフラッシュランプやＬＥＤ光源も利用できる。光学系１５０３や射出端１５０４としては、バンドルファイバ、プリズム、ミラー、レンズ等の光学部材が好適である。探触子１５０１が備える受信素子としては例えば、圧電素子、静電容量型素子、ファブリペロー型素子などを利用できる。探触子と画像処理装置の間や、画像処理装置内部に、光音響信号に対して増幅処理やデジタル変換処理を施す信号処理回路を設けることが好ましい。

射出端１５０４を移動させて光の照射位置を変化させるための変化手段１５０８を設け
、被検体１５００と射出端との相対位置を変化させることにより、ＦＯＶの位置が少しずつ移動し、被検体全域が計測される。変化手段は、探触子を射出端と共に移動させても良い。変化手段として例えば、モータや制御ステージなどの駆動部品が好適である。変化手段は、被検体の方を移動させても良い。あるいは変化手段は、光の照射方向を変化させることで光の照射位置、ひいてはＦＯＶを変化させても良い。

（赤外線カメラ）
赤外線カメラ１２０は、被検体の外観及び被検体の体表近傍の血管の様子を静止画像または動画像として撮影し、画像処理装置１０００に入力する。赤外線カメラ１２０は、被検体全体の外観を撮影可能な位置に設置される。本実施形態では３台の赤外線カメラを利用する。３台を区別する必要がある場合、それぞれ符号３０１，３０２，３０３と表す。

図３は、赤外線カメラによる撮像の様子を示す。赤外線カメラ３０１の撮像方向３０４は、装置座標系C_DEVのＺ方向に一致する（すなわち、カメラ座標系C_CAM1は、Ｚ軸が装置
座標系C_DEVの−Ｚ方向を向いている）。赤外線カメラ３０２の撮像方向３０５は、装置座標系C_DEVの−Ｘ方向に一致する（すなわち、カメラ座標系C_CAM2は、Ｚ軸が装置座標系C_DEVのＸ方向を向いている）。赤外線カメラ３０３の撮像方向３０６は、装置座標系C_DEVの
Ｙ方向に一致する（すなわち、カメラ座標系C_CAM3は、Ｚ軸が装置座標系C_DEVの−Ｙ方向
を向いている）。ここで、カメラ座標系C_CAM1、C_CAM2、C_CAM3から装置座標系C_DEVへの座
標変換を、それぞれ、T_C1toD、T_C2toD、T_C3toDと定義する。各カメラは装置座標系C_DEVにおいて校正済みであり、上記の座標変換の情報や夫々のカメラの内部パラメータは、予め画像処理装置１０００が保持している。

図４に、血管４０１が写し出された赤外線カメラ画像４００（赤外線画像とも呼ぶ）の例を示す。この赤外線カメラは、近赤外線の強度情報を画像化できる。ここで近赤外線には、皮膚をある程度透過する性質と、ヘモグロビンが含まれる血管部分に吸収される性質がある。そのため赤外線カメラを用いると、皮膚直下の静脈血管（表在血管）を可視化できる。また、血管部分は周囲よりも暗く写し出される。よって赤外線カメラ画像は、皮膚下の表在血管の形状を明瞭に描出する形態画像として扱える。

以下、赤外線カメラ３０１、３０２、３０３による撮影画像をそれぞれ、I_CAM1、I_CAM2、I_CAM3と表記する。カメラ画像I_CAM1、I_CAM2、I_CAM3の座標系（２次元）はそれぞれ、C_IMG1、C_IMG2、C_IMG3とする。例えば赤外線カメラ３０１に関して、カメラ画像座標系C_IMG1（２次元）上の座標は、カメラ座標系C_CAM1（３次元）におけるカメラの視線（原点と、
３次元空間中の投影面上の点と、を通る直線）と一対一の関係を持つ。このようなカメラ画像座標系とカメラ座標系の間の変換に関しては、一般的な座標変換方法を利用できる。カメラ座標系C_CAM1、C_CAM2、C_CAM3からカメラ画像座標系C_IMG1、C_IMG2、C_IMG3への変換を、それぞれ、T_C1toI1、T_C2toI2、T_C3toI3と定義する。

（表示装置）
表示装置１２０は、液晶装置やＣＲＴなどのディスプレイである。表示装置１２０は、画像処理装置１０００が出力する画像データにもとづき被検体の画像を表示する。表示装置はシステムの一部であってもよく、システム外部に存在してもよい。

（画像処理装置）
画像処理装置１０００は、信号取得部１０１０、光音響画像取得部１０２０、広域画像取得部１０３０、計測位置取得部１０４０、比較ペア生成部１０５０、投影画像生成部１０６０を備える。装置はさらに、比較部１０７０、パルス位置推定部１０８０、統合ボリューム生成部１０９０、表示制御部１１００を備える。画像処理装置は、ＣＰＵやメモリ等の演算資源を備え、プログラムの指示に従って所定の情報処理を行う情報処理装置（Ｐ
Ｃやワークステーション等）により構成される。画像処理装置の各ブロックは、それぞれ独立した回路として構成されても良く、プログラムモジュールとして仮想的に実現されても良い。

信号取得部１０１０は、光音響信号計測装置１１０が計測した光音響信号を取得し、光音響画像取得部１０２０に出力する。
光音響画像取得部１０２０は、１回のパルス光照射で取得した光音響信号に基づく画像再構成処理を行い、対応するＦＯＶの光音響画像をボリュームデータの形式で生成する。なお、パルス光照射毎に生成されるボリュームデータを以下ではパルスボリュームと呼ぶ。光音響画像取得部１０２０は、生成したパルスボリュームを投影画像生成部１０６０および統合ボリューム生成部１０９０に局所画像として出力する。

広域画像取得部１０３０は、赤外線カメラ１２０が撮像した赤外線カメラ画像（２次元画像）を広域画像として取得し、比較部１０７０に出力する。
計測位置取得部１０４０は、光音響信号計測装置１１０の信号計測部の位置情報（信号計測位置）を取得し、比較ペア生成部１０５０に出力する。

比較ペア生成部１０５０は、パルスボリュームのペアの情報（比較ペア情報）を生成し、比較部１０７０に出力する。
投影画像生成部１０６０は、ボリュームデータから投影画像を生成し、比較部１０７０に出力する。

比較部１０７０は、比較ペア情報に基づき、ペアとなったパルスボリュームを比較する処理（第１の比較）として、両者の投影画像の一致度に関する情報を算出する。この情報は、投影画像の元となった比較ペア間の一致度に関する情報であり、比較ペア間の位置ズレに関する情報とも呼べる。比較部１０７０はまた、夫々のパルスボリュームと赤外線カメラ画像とを比較する処理（第２の比較）として、夫々のパルスボリュームの投影画像と赤外線カメラ画像との一致度に関する情報を算出する。この情報は、投影画像の元となった夫々のパルスボリュームと赤外カメラ画像との一致度に関する情報であり、夫々のパルスボリュームと赤外カメラ画像の間の位置ズレに関する情報とも呼べる。算出された各情報は、パルス位置推定部１０８０に出力される。

パルス位置推定部１０８０は、比較部１０７０が算出したペア間の一致度に関する情報と、夫々のパルスボリュームと赤外線カメラ画像の一致度に関する情報と、の両方に基づいて、夫々のパルスボリュームの位置の推定値（パルス位置推定量）を算出する。算出した推定値は、統合ボリューム生成部１０９０に出力される。なお、パルス位置推定量を算出することは、信号計測位置に対するパルスボリュームの位置の補正量を算出することと等価である。信号計測位置の計測誤差が無視できる場合には、この補正量は、夫々のパルス光照射時の被検体の体動を表す。一方、被検体の体動が無視できる場合には、この補正量は、信号計測位置の計測誤差を表す値となる。

統合ボリューム生成部１０９０は、パルス位置推定量に基づいてパルスボリュームを統合した統合ボリュームを生成し、表示制御部１１００に出力する。
表示制御部１１００は、統合ボリューム生成部１０９０が生成した統合ボリュームを表示装置１３０に表示するための表示制御を行う。

（処理フロー）
図５は、画像処理装置１０００が行う全体の処理手順を示すフローチャートである。
(ステップＳ５０１０）信号データ取得
信号取得部１０１０は、光音響信号計測装置１１０が被検体を複数回撮像して得た光音
響信号データ（以下、信号データと称する）を取得する。本実施形態における撮像回数は、N_pulse回であり、信号データの個数は、N_pulse個である。１回の（i回目の）パルス
照射によって受信した信号データをP_i（1≦i≦N_pulse）と表記する。信号データP_iは
、探触子がある所定の位置に配置されたときに、複数の受信素子が受信した時系列の音響信号である。また、P_iの添え字i（1≦i≦N_pulse、iは正の整数）を「パルスインデックス」と呼ぶ。各信号データは、多次元のデータ列等の一般的な形態で表現される。

(ステップＳ５０２０）パルスボリューム再構成
光音響画像取得部１０２０は、信号データP_i（1≦i≦N_pulse）を用いて再構成処理を行い、パルスボリュームVp_i（1≦i≦N_pulse）を生成し、画像処理装置１０００に記録
する。再構成には一般的な方法、例えばトモグラフィ技術で普通に用いられるタイムドメインあるいはフーリエドメインでの逆投影法などを利用できる。この時、各受信素子の位置が考慮される。生成されるパルスボリュームVp_iは３次元ボリュームの画像データであり、Vp_i(x，y，z)のような関数表現によってデータの値（ボクセル値）を表記できる。
本実施形態でのパルスボリュームは、被検体を撮像可能な範囲（ＦＯＶ）を含む局所領域の３次元ボリュームである。

(ステップＳ５０３０）赤外線画像取得
広域画像取得部１０３０は、赤外線カメラ３０１、３０２、３０３が撮影した被検体の赤外線カメラ画像（I_CAM1，I_CAM2，I_CAM3）を取得する。ここでは、赤外線カメラ画像は
ある瞬間の静止画像とする。赤外線静止画像取得のタイミングは、被検体の状態（位置・姿勢など）が、光音響計測の時点と変わらない時間内とする。ただし、赤外線カメラが動画像を撮影し、広域画像取得部１０３０が動画像のフレーム中から好適な静止画像を選択しても良い。その場合、各フレーム間で画素値の差分をとるなどの手法で画像の変化が小さいフレームを検出することが好ましい。

(ステップＳ５０４０）計測位置取得
計測位置取得部１０４０は、光音響信号計測装置１１０の信号計測位置を取得する。信号計測位置とは、ステップＳ５０１０で取得したN_pulse個の信号データそれぞれの計測
時点における信号計測部の位置であり、装置座標系C_DEVにおける３次元の位置情報である。ステップＳ５０１０で取得した信号データや、それに基づいて生成されたパルスボリュームは、この信号計測位置を基準とした所定範囲（ＦＯＶ）における計測結果である。本実施形態では、信号計測位置をPosS_i（1≦i≦N_Pulse）と表記する。信号計測部が位置
変化手段によって移動しながら光音響計測を行う場合はパルスごとに信号計測位置が異なる。

(ステップＳ５０５０）比較ペア情報生成
比較ペア生成部１０５０は、信号計測位置PosS_iを取得し、後続の比較処理を行う比較ペアを決定して比較ペア情報を生成する。ペアの決定基準として例えば、互いにオーバーラップするパルスボリュームをペアとする方法がある。

図６は比較ペア情報の生成方法を説明する図である。紙面による説明の都合上、２次元に簡略化して説明する。図中の符号６０１から６０５は夫々、信号計測位置PosS_i（i=1
〜5）に応じたパルスボリュームV1からV5の領域を表す。比較ペア生成部１０５０は、こ
れらのパルスボリュームについて、オーバーラップ領域を有するパルスボリュームのペアを検出する。ここでは、「V1とV2」、「V2とV3」、「V3とV4」、「V3とV5」、「V4とV5」の領域がオーバーラップしている。そこで、比較ペア生成部１０５０は、これらのペアを示す情報を、パルスインデックスの組として生成する。図６であれば、R_1={1、2}、R_2={2、3}、R_3={3、4}、R_4={3、5}、R_5={4、5}となる。

ただし比較ペア情報の生成方法はこれに限られない。他に例えば、オーバーラップ領域の体積が所定の値以上となるパルスボリューム同士をペアとする方法がある。また、オーバーラップ領域の体積の割合が所定の値以上となるパルスボリューム同士を比較ペアとする方法がある。これらの方法によれば、オーバーラップ領域の量や割合が少ない場合はペアにならないため、パルスボリューム間の位置の推定処理が効率化・安定化する。

(ステップＳ５０６０）赤外線カメラ方向への投影画像生成
ステップＳ５０６０において、投影画像生成部１０６０は、パルスボリュームVi（1≦i≦N_pulse）の夫々について、ボリューム値を投影して投影画像を生成する。例えば、最
大輝度値投影（ＭＩＰ：Maximum Intensity Projection）画像を生成する。これには透視投影法など、一般的な手法を利用できる。投影画像の投影方向は、各赤外線カメラの撮像方向（３０４、３０５、３０６）と略同一とする（好ましくは一致させる）。つまり、装置座標系C_DEVのＺ軸方向、Ｘ軸方向、Ｙ軸方向夫々において投影画像を生成する。これは、後述するステップＳ５０８０において、投影画像（局所領域画像）と赤外線カメラ画像（被検体広域画像）とを同一方向から比較するためである。

透視投影法の場合についてより詳しく述べる。上述の通り、赤外線カメラ座標系から装置座標系への座標変換は、T_C1toD、T_C2toD、T_C3toDと定義されている。また、各カメラ座標系（３次元）からカメラ画像座標系（２次元）への変換は夫々、T_C1toI1、T_C2toI2、T_C3toI3と定義されている。従って、パルスボリュームの各ボクセルの座標に信号計測位置PosS_iだけの並進を与えることで装置座標系に変換し、夫々T_C1toD、T_C2toD、T_C3toDの逆
変換を掛けて赤外カメラ座標系に変換する。そして、その結果にさらに夫々T_C1toI1、T_C2toI2、T_C3toI3の変換を掛けることで、赤外線カメラ画像座標系に変換できる。すなわち
、パルスボリュームの各ボクセルが赤外カメラ画像上のどの位置に透視投影されるかを算出できる。こうして生成された透視投影による投影画像は、２次元座標への投影方法が赤外線カメラ画像と概ね一致する。

図７は、各赤外線カメラと、パルスボリュームと、生成された投影画像との空間的関係を示す。パルスボリューム７００は、装置座標系C_DEVのＸ、Ｙ、Ｚ軸で表される３次元空間中のボリュームデータである。投影画像７０１は、ボリュームデータ７００を装置座標系C_DEVのＺ軸方向、つまり赤外線カメラ３０１の撮像方向３０４に最大値投影した画像であり、装置座標系C_DEVのＸ軸とＹ軸によって表される２次元空間中の画像データである。同様に、投影画像７０２は、Ｙ軸とＺ軸によって表される２次元画像データ、投影画像７０３は、Ｚ軸とＸ軸によって表される２次元画像データである。

本実施形態では、ボリュームデータを装置座標系C_DEVのＸ、Ｙ、Ｚ方向の夫々に透視投影して生成した投影画像７０２、７０３、７０１を夫々、Ipx_i、Ipy_i、Ipz_i（1≦i≦N_pulse）と表記する。このとき、投影画像Ipx_iは赤外線カメラ画像I_CAM2、投影画像Ipy_iは赤外線カメラ画像I_CAM3、投影画像Ipz_iは赤外線カメラ画像I_CAM1、の撮像方向と同一方向に投影された画像である。また、本実施形態において画像データの値（画素値）に言及する場合、投影画像Ipx_i，Ipy_i，Ipz_iを、Ipx_i(y，z)，Ipy_i(x，z)，Ipz_i(x，y)のように関数として表記する。

本ステップでは、パルスボリュームをＸ、Ｙ、Ｚ軸の夫々の方向（本実施形態では各赤外線カメラの撮像方向と略一致する方向）に正射影によって投影した投影画像も生成する。これは、パルスボリュームに写る血管の陰影を強調して、後述するステップＳ５０７０でのパルスボリューム間の比較を容易にするためである。これら正射影の投影画像を夫々、Iox_i、Ioy_i、Ioz_i（1≦i≦N_pulse）と表記する。透視投影による投影画像（Ipx_i
、Ipy_i、Ipz_i）および正射影による投影画像（Iox_i、Ioy_i、Ioz_i）は、比較部１０
７０に出力される。

(ステップＳ５０７０）投影画像間の類似度のピーク位置取得
比較部１０７０は、比較ペア情報R_j（1≦j≦N_pair）が表す全てのペア（パルスイン
デックスの対）について、比較ペアとされたパルスボリュームを比較し、ペア間の一致度に関する情報を取得する。具体的には、両者の投影画像の類似度関数と、その類似度関数のピーク位置を取得する。そしてこの情報を、投影画像間の一致度に関する情報、すなわち、投影画像の元となったパルスボリュームのペア間の一致度に関する情報とする。これは、被検体の局所領域が写った画像同士の比較処理に当たる。

類似度関数としては、比較ペアとされたパルスボリュームの正射影による投影画像（Iox_i、Ioy_i、Ioz_i）同士を処理対象として、これらの位置をずらしながら重ねたときの
類似度を値とした画像類似度関数を用いる。具体的には、比較ペア情報R_jが表すパルス
インデックスの対（R_j，1およびR_j，2と表記する）の夫々の投影画像について、以下の３個の類似度関数を取得する。

ここで、f_simil(I1，I2，x，y)は、画像I1に対して画像I2の相対位置を(x，y)だけ並
進させた場合の画像間の類似度を算出する関数である。例えば、ＳＳＤ（Ｓｕｍ ｏｆ Ｓｑｕａｒｅｄ Ｄｉｆｆｅｒｅｎｃｅ）やＳＡＤ（Ｓｕｍ ｏｆ Ａｂｓｏｌｕｔｅ
Ｄｉｆｆｅｒｅｎｃｅ）、相互情報量、相互相関など、任意の類似度尺度が適用できる。ここで関数f_similは、画像間の類似度が高い場合には関数値として高い値を返すものと
する。

類似度関数(FLX_j，FLY_j，FLZ_j)の取得とは、具体的には、各関数の引数である並進
量(x，y，z)を所定の範囲内で離散的に変化させた場合の関数値の取得を意味する。比較
ペア（R_j，1およびR_j，2）の夫々のパルスボリュームの信号計測位置はPosS_R_j，1お
よびPosS_R_j，2であり、互いに異なっている。

本ステップでは、この信号計測位置の相対位置（R_j，1から見たR_j，2の位置、ΔPosS_j = PosS_R_j，2 - PosS_R_j，1）を並進量の基準値として、その前後の限定的な範囲内で投影画像を並進させて類似度を計算する。基準値は中央値などである。例えば、FLX_j(y，z)の場合には、yとzの値を、前記相対位置のy，z成分を中心として夫々を「-Kから+K
」までの整数値だけ前後に変化させた場合にFLX_jが返す、「(2K+1)×(2K+1)個」の値の
集合を取得する。この値の集合は、離散画像のようなメモリアレイ上に展開されるデータである。より発展的には、「(2K+1)×(2K+1)個」の値の集合に対してバイリニア法やバイキュービック法などを適用して、類似度関数FLX_j(y，z)を連続関数（またはそれに近い
関数）として導出してもよい。本実施形態では、比較ペアごとに、各類似度関数(FLX_j，FLY_j，FLZ_j：1≦j≦N_pair)を２次元の連続関数として取得する。

図８に、比較部１０７０によって比較される画像の関係を示す。符号８０１、８０２は夫々、比較ペアであるパルスボリュームの夫々を装置座標系C_DEVのＺ方向に投影した投影画像であり、正射影による投影画像と透視投影による投影画像の両方を含む。矢印８０３
は、投影画像８０１（正射影）と投影画像８０２（正射影）とを比較することを表す。このように本ステップでは、比較ペアとなっている投影画像の間の血管走行の類似性を比較することで、類似度関数が算出される。

次に比較部１０７０は、比較ペアの夫々に関する類似度関数が最大となる位置（類似度ピーク位置）を算出する。具体的には、上記のように算出した各比較ペアに関する類似度関数(FLX_j，FLY_j，FLZ_j：1≦j≦N_pair)の夫々について、次式に示すように、関数値
が最大となる位置PosLX_j，PosLY_j，PosLZ_jを算出する。

次に比較部は、これらの値を統合して、処理対象の比較ペアに関する類似度ピーク位置PosL_jを、パルスボリューム間の並進量を表す３次元ベクトルとして算出する。この統合は次式のように、座標軸毎に得られた最大値の平均を算出することで実施できる。

以上の処理により、比較部１０７０は、全ての比較ペアに関する類似度ピーク位置PosL_j（1≦j≦N_pair)を算出する。PosL_jの値は、局所画像に基づいて求めた、比較ペア間
の相対位置の推定値（個別最適値）を示している。比較ペア間で体動が生じていない場合には、類似度ピーク位置PosL_jの値は、信号計測位置のペア間の差異ΔPosS_jと一致することが期待される。言い換えると、PosL_jとΔPosS_jの差異（ΔPosS_j‐PosL_j）が、局所画像に基づいて取得した当該ペア間での（相対的な）体動の観測情報となる。

(ステップＳ５０８０）投影画像と赤外線画像間の類似度のピーク位置算出
比較部１０７０は、パルスボリュームの夫々と赤外カメラ画像を比較し、パルスボリュームの夫々と赤外カメラ画像との一致度に関する情報を取得する。つまり、各パルスボリュームの透視投影による投影画像Ipx_i、Ipy_i、Ipz_i（1≦i≦N_pulse）と、投影方向が夫々対応する赤外線カメラ画像I_CAM2、I_CAM3、I_CAM1との類似度関数と、その類似度関数
のピーク位置を得る。そしてこの情報を、投影画像の夫々と赤外線カメラ画像の一致度に関する情報、すなわち、投影画像の元となったパルスボリュームの夫々と赤外カメラ画像との一致度に関する情報とする。ここで、投影画像にも赤外線カメラ画像にも被検体の血管が描出されているため、血管の輝度情報を用いて両画像の類似性を評価できる。ただし本実施形態での血管の輝度は、投影画像中では周囲より明るく、赤外線カメラ画像中では周囲より暗い。そこで、赤外線カメラ画像と投影画像のいずれかの輝度値の明暗を反転させて比較を行うと良い。

パルスボリュームの夫々は被検体の局所領域が写った画像（局所画像）であり、赤外線カメラ画像は被検体の全体領域が写った画像（広域画像）である。本ステップの目的は、赤外線カメラ画像に対して各パルスボリュームの投影画像の位置をずらしながら重ねたと
きの類似度を計算することにより、被検体の広域画像と局所画像を比較することである。まず、各パルスボリュームに対応する夫々の投影画像について、以下の３個の類似度関数を取得する。

ここで、f_simil(I1，I2，x，y)は、ステップＳ５０７０と同様の関数である。本ステ
ップでも、各類似度関数(FGX_i，FGY_i，FGZ_i)の取得とは、各関数の引数である並進量(x，y，z)を所定の範囲内で離散的に変化させた場合の関数値の取得を意味する。上式より、赤外線カメラ画像に対して投影画像の位置を並進させた場合の、画像間の類似度を算出する関数が得られる。ここで、ステップＳ５０６０の処理において、信号計測位置を考慮して赤外線カメラの撮像方向への透視投影によって生成された投影画像は、既に赤外線カメラ画像と同一座標系で表されている。従って本ステップの類似度は、パルスボリュームの絶対位置ではなく、信号計測位置からの相対位置をパラメータとして算出する。すなわち、位置ズレのない状態を表す(x，y，z)＝(0，0，0)を基準として、限定的な範囲内で投影画像を並進させて類似度を計算する。本実施形態では、N_pulse個のパルスボリューム
に対応する夫々の投影画像について、各類似度関数(FGX_i，FGY_i，FGZ_i：1≦i≦N_pulse)を、２次元の連続関数として取得する。

図８において、赤外線カメラ画像４００ならびに投影画像８０１および８０２には、被検体の血管が描出されている。矢印８０４、８０５は夫々、赤外線カメラ画像４００に対する投影画像８０１（透視投影）の比較と、赤外線カメラ画像４００に対する投影画像８０２（透視投影）の比較を表す。このように本ステップでは、赤外線カメラ画像と投影画像の間の血管走行の類似性を比較して、類似度関数を算出する。

次に、比較部１０７０は、パルスボリュームの夫々に関する類似度関数が最大となる位置（類似度ピーク位置）を算出する。具体的には、算出した各パルスボリュームに関する類似度関数(FGX_i，FGY_i，FGZ_i：1≦i≦N_pulse)の夫々について、次式のように、関数値が最大となる位置PosGX_i，PosGY_i，PosGZ_iを算出する。

次に、これらの値を統合して、処理対象のパルスボリュームに関する類似度ピーク位置PosG_iを、パルスボリュームの信号計測位置からの並進量を表す３次元ベクトルとして算出する。この統合は、ステップＳ５０７０と同様、次式のように、座標軸毎に得られた最大値の平均を求めることにより実施できる。

以上の処理により、比較部１０７０は、赤外線カメラ画像に対する全ての投影画像の類似度ピーク位置PosG_i（1≦i≦N_pulse)を算出する。なお、当該パルスボリュームで体動が生じていない場合には、類似度ピーク位置PosG_iの値は (0，0，0)となることが期待される。言い換えると、信号計測位置PosS_i とPosG_iの和が、広域画像に基づく当該パル
スボリュームの推定位置（個別最適値）を示している。また、PosG_iの値そのものが、広域画像に基づいて取得した当該パルスボリュームにおける（絶対的な）体動の観測情報となる。

(ステップＳ５０９０）２種類のピーク位置に基づきパルス毎の体動推定
パルス位置推定部１０８０は、Ｓ５０７０の類似度ピーク位置PosL_j（1≦j≦N_pair）と、Ｓ５０８０の類似度ピーク位置PosG_i（1≦i≦N_pulse）から、パルス位置推定量Pos_i（1≦i≦N_pulse）を算出する。すなわち、被検体の局所画像（パルスボリューム）間
の相対的な位置に関する個別最適値と、被検体の広域画像（赤外線カメラ画像）に対する各局所画像の位置に関する個別最適値と、の両方をなるべく保ちつつ、全てのパルスボリュームの位置を全体最適化する。このパルス位置推定量Ｐｏｓ＿ｉは、パルスボリューム間の比較で最も類似する位置関係（矢印８０３：第１の比較）と、赤外線カメラ画像とパルスボリュームの比較で最も類似する位置関係（矢印８０４、８０５：第２の比較）の両方をなるべく満たす値である。

具体的には、コスト関数Eを次式のように定義し、このコスト関数Eを小さくするように各パルスボリュームの位置を最適化する。

ここで、αは第１項（第１の比較）と第２項（第２の比較）をバランスさせるための重み係数であり（0≦α≦1）、本実施形態では、α＝０．５とする。関数EL_jは、被検体の局所画像間の位置ズレによるコストを算出する関数である。この値は、類似度ピーク位置PosL_jと、当該比較ペアである２個のパルスボリュームの位置推定量の相対値（すなわち、Pos_R_j，2 - Pos_R_j，1）との相違等によって算出される。つまり、被検体の局所画
像間の相対的な位置の個別最適値からの、推定値のズレ量をコストとする関数である。

関数EL_jは、具体的には次式によって計算する。

ここで、a_jはj番目の比較ペアの1番目のパルスインデックス値R_j，1であり、同様にb_jはj番目の比較ペアの2番目のパルスインデックス値R_j，2である。

また、関数EG_iは、被検体の広域画像に対する局所画像の位置ズレによるコストを算出する関数である。この値は、類似度ピーク位置PosG_iのみに基づくパルスボリュームの推定位置（すなわち、PosS_i＋PosG_i）と、当該パルスボリュームの位置推定量Pos_iとの
相違等によって算出される。つまり、被検体の広域画像に対する局所画像の位置の個別最適値からの、推定値のズレ量をコストとする関数である。関数EG_iは、具体的には次式によって計算する。

なおコスト関数Eの最適化には、既知の様々な手法、例えば、一般的な最急降下法やニ
ュートン法などの繰り返し計算による非線形最適化法や、線形最適化法を利用できる。また、コスト関数Eの定義は上記に限らない。例えば、式１６や式１７はズレ量の評価尺度
としてＬ２ノルムを用いているが、Ｌ１ノルム等の他の距離尺度を用いても良い。また、例えば式１５に加え、各パルスボリュームの位置の変動（動き）に対する正則化を掛けても良い。一般的に撮像中の体動は極端に大きくはなく、連続的で滑らかだと想定される。よってパルスボリュームの位置の最適化の際には、上記の想定から逸脱するような動きが算出されることを防ぐために、正則化処理を行う。正則化手法として例えば、位置の補正量（信号計測位置PosS_iとパルス位置推定量Pos_iの差）の総和に所定の重み係数をかけ
て式１５に加算する方法がある。また、補正量の微分（加速度）の総和を用いる方法、補正量の周波数成分値に基づいて算出した値を用いる方法などがある。また、被検体の典型的な変動（体動）の仕方をモデルとして用意し、そのモデルと補正量の相違をコストとして式１５に加算する方法もある。

以上に説明した方法により、パルス位置推定部１０８０は、各パルスボリュームに関するパルス位置推定量Pos_i（1≦i≦N_pulse）を算出する。

(ステップＳ５１００）統合ボリューム生成 （最終出力ボリューム）
統合ボリューム生成部１０９０は、パルスボリュームVp_i（1≦i≦N_pulse）を統合し
て統合ボリュームVoを生成する。このとき統合ボリューム生成部１０９０は、各パルスボリュームVp_iの夫々を、ステップＳ５０９０で最後に得たパルス位置推定量Pos_iだけ並
進して共通の座標系に変換した後に統合する。統合ボリューム生成部１０９０は、まず、統合ボリュームの範囲として、並進後の各パルスボリュームの領域を包含する包含領域の範囲を算出する。次に、この包含領域において、同一座標に変換された各パルスボリュームのボクセル値を平均化して統合ボリュームを生成する。統合ボリュームは、表示制御部１１００に送信されると共に、不図示の記憶媒体に保存される。

(ステップＳ５１１０）表示
表示制御部１１００は、出力画像Voの情報を表示装置１３０に出力して表示させる。表示方法としては、統合ボリュームを所定の方向に最大値投影処理して得られる２次元画像（投影画像）を表示する方法、統合ボリュームをボリュームレンダリングする方法、統合ボリュームを任意の断面で切断した断層画像を表示する方法などがある。

以上、本実施形態の画像処理装置によれば、撮像中に生じる体動の補正処理において、被検体の局所領域の画像の間を比較すると同時に、被検体全体を撮像した画像と局所領域の画像を比較する。これにより、本来の被検体の解剖構造により合致する全体として歪みが抑制された画像を取得できる。

（変形例１−１）：並進量を変化させた投影画像を夫々生成してもよい
ステップＳ５０８０では、ステップＳ５０６０で生成された同じ投影画像Ipx_i，Ipy_i，Ipz_iの位置を並進させながら類似度関数を算出していた。しかし、装置座標系C_DEVに
おいてパルスボリュームの位置が変わると、透視投影によって生成される投影画像も変化する。そこで、ステップＳ５０６０において、パルスボリュームVi（1≦i≦N_pulse）の
夫々について、装置座標系C_DEVにおいて、信号計測位置からの並進量(x，y，z)を所定の
範囲内で離散的に変化させて、夫々の方向の投影画像を生成する。

そしてＳ５０８０において、各類似度関数を取得する際の信号計測位置からの並進量(x、y、z)を装置座標系C_DEV上で定義する。そして、個別の並進量に対応する類似度算出に
おいて、その並進量に対応する投影画像を選択した上で類似度を算出する。例えばステップＳ５０８０で類似度関数FGX_j(y_k，z_k)を算出する場合、赤外線カメラ画像I_CAM2との類似度を計算する投影画像Ipx_iとして、ステップＳ５０６０で生成した並進量(0,y_k,z_k)に対応する投影画像を選択する。他の方向についても同様に類似度関数を算出できる。本変形例の方法であれば、パルスボリュームの位置をより反映させた類似度が取得できる。

（変形例１−２）：赤外線カメラは装置座標系に直交していなくても良い
上記第１の実施形態では、３台の赤外線カメラの撮像方向は、装置座標系の座標軸と一致していた。しかし、赤外線カメラ画像内に被検体の外形全体が含まれるのであれば、撮像方向はこれに限定されない。例えば装置の構造によっては、赤外線カメラ３０１のＺ軸方向（装置座標系C_DEV）に被検体を配置できない場合がある。この場合、赤外線カメラ３０１の向く方向を、Ｙ軸を中心としてＺ軸を４５°回転させた方向とする。この方向を表すベクトルを、dと定義する。

ステップＳ５０６０で生成される投影画像の投影方向は、赤外線カメラの撮像方向と一致する。そこで、パルスボリュームを上述のベクトルdの方向に透視投影で投影すること
で、Ｚ方向への投影画像であるIz_iの代わりに、投影画像Id_iを生成する。次に、ステップＳ５０７０とＳ５０８０の類似度関数の算出では、Ｚ方向の投影画像Iz_iに対応する類似度関数FLZ_j(x,y)及びFGZ_i(x,y)の代わりに、下式のFLD_j(x,y,z)、FGD_i(x,y,z)を取得する。

ここで、f_simil(I1，I2，x，y，z)は、画像I1に対して画像I2の相対位置を(x，y，z)
だけ並進させた場合の画像間の類似度を算出する関数である。これは、ステップＳ５０７０で説明したf_simil(I1，I2，x，y)と比べ、並進軸が一つ多い。このときの並進量(x，y，z)は装置座標系C_DEV上で定義しておく。ここで、関数FLD_j(x，y，z)、及びFGD_i(x，y，z)を取得するために、並進量(x，y，z)の値を所定の範囲で離散的に変化させて類似度
を計算する。ただし、ベクトルdに平行な方向での並進では投影画像の相対位置は変わら
ない。そこで、並進量(x，y，z)を、ベクトルdに直交する平面上で、並進量が０である位
置を基準にした所定の領域内（例えば矩形領域）で変化させる。これにより、類似度関数FLD_j(x，y，z)及びFGD_i(x，y，z)が取得できる。他の赤外線カメラ３０２、３０３についても同様に処理することで、装置座標系と赤外線カメラの撮像方向が一致しない場合でも、赤外線カメラ画像とパルスボリュームを比較できる。

（変形例１−３）：類似度関数は輝度値を使った画像類似度で無くても良い
上記第１の実施形態では、ステップＳ５０６０の処理として、パルスボリュームの投影画像としてＭＩＰ画像を生成した。しかし、例えばＭＩＰ画像に代えて、ＭｉｎＰ画像（最小値投影画像）を生成しても良い。また、ＭＩＰ画像とＭｉｎＰ画像を両方生成し、以降の処理では両方または好ましい片方を用いても良い。

また、上記第１の実施形態のステップＳ５０７０やＳ５０８０では、画素値の類似度に基づいて画像間の類似度関数を求めた。しかし例えば、被検体の血管等の解剖学的な特徴を画像処理によって抽出して、その分布や位置の一致度に基づいて類似度関数を取得しても良い。このとき、エッジ検出やコーナー検出など、既知の画像処理技術によって特徴を抽出しても良い。またコンピュータビジョン等の技術分野で一般的な、SIFT特徴やSURF特徴などの高次な局所画像特徴を用いても良い。これらの方法によれば、画像間の輝度分布の相違やノイズの混入などに対して、より頑健な類似度関数を取得できる。さらにＳ５０７０では、投影画像ではなく、元のパルスボリューム同士をそのまま比較してもよい。

（変形例１−４）：並進だけでなく、回転も含んで良い
上記第１の実施形態では、ステップＳ５０７０、Ｓ５０８０で、パルスボリューム間の並進量と、赤外線カメラ画像に対する投影画像の並進量と、に関する類似度関数を取得した。しかし、位置関係として、並進に加えて回転を考慮しても良い。この場合、ステップＳ５０７０及びＳ５０８０で取得する類似度関数は、各投影画像間の、投影面上での並進量および投影面内の回転量に関する類似度関数となる。この場合、ピーク位置は、３次元ベクトルとしてではなく、例えば、投影画像間の類似度関数を最大とする剛体変換行列、及び赤外線カメラ画像と投影画像間の類似度関数を最大とする剛体変換行列として算出する。

またステップＳ５０９０では、パルス位置推定量の代わりに、各パルスボリュームの並進と回転（すなわち剛体変換）を算出する。この算出処理は、以下の２つのコストを有するコスト関数の最適化によって実施できる。１つ目は、投影画像間の剛体変換の差異を表す行列と、ステップＳ５０７０で求めた投影画像間の剛体変換行列との整合性を表すコストである。２つ目は、赤外線カメラ画像と投影画像の剛体変換の差異を表す行列と、ステップＳ５０８０で求めた赤外線カメラ画像と投影画像間の剛体変換行列との整合性を表すコストである。

さらに、ステップＳ５１００では、前記剛体変換に基づいてパルスボリュームを統合する。ただし、ステップＳ５０７０およびステップＳ５０８０において回転を含む類似度関数を取得する場合には、各投影面での投影画像の回転が互いに独立ではないことを考慮することが望ましい。すなわち、ある投影面における投影画像の面内回転は、当該パルスボリュームにおける前記投影面の法線方向を軸とした回転を意味する。この回転は、前記投影面とは異なる投影面における投影画像に対して、必ずしも並進や面内回転として表すことができない。そのため、他の投影面における投影画像は、当該回転を考慮して改めて生成し直すことが望ましい。

以上の方法によれば、撮像時の被検体の動きが並進だけではなく、回転も含むような場合に、より画質の高い統合ボリュームを生成できる。

（変形例１−５）１軸方向の投影画像を省略してもよい
上記第１の実施形態では、ステップＳ５０６０の処理において、直交三方位に投影した投影画像を生成していた。しかし例えば、Ｚ軸方向とＹ軸方向のみ投影画像を生成し、Ｘ軸方向への投影を省略しても良い。この場合でも、Ｚ軸方向の投影画像に基づいてＸＹ平面上の並進の情報が取得でき、またＹ軸方向の投影画像に基づいてＸＺ平面上の並進の情報が取得できる。そのため、全ての軸方向の動きの情報を取得できる。この方法により第１の実施形態に比べて計算コストを低減できる。この場合、赤外線カメラは、Ｚ軸方向とＹ軸方向の２台のみでよい。

また、被検体の体動に関して予め制約を仮定できる場合には、さらに計算コストを低減できる。例えば、撮影時の被検体の姿勢や、撮影装置自体による物理的な拘束等により被検体の体動がＺ軸方向には生じない（または極めて小さい）場合には、Ｚ軸方向の投影画像だけを生成し、ＸＹ方向の体動だけを対象として処理を実行する。これにより、体動が生じない方向の計算を省略できる。この場合、赤外線カメラは、Ｚ軸方向の１台のみでよい。また、詳細な体動を求める処理としての局所画像同士の（比較ペア同士の）比較は複数方向で行い、全体の整合性をとるための処理としての広域画像と局所画像の比較は１方向（例えばＸＹ方向）のみで行うという構成であってもよい。これによると、体動補正を３次元で行いつつ、最小の装置構成（１台の赤外カメラ）で発明の効果を得ることができる。

また、被検体の体動が１軸方向でしか生じないことが仮定できる場合には、少なくとも当該軸の動きが観察できる方位に投影した投影画像だけで処理が可能である。例えば、被検体がＸ軸方向にしか動かない場合、Ｙ軸方向またはＺ軸方向の投影画像だけでＸ軸方向の体動を算出できるので、計算コストを低減できる。

＜第２の実施形態＞
本発明の第２の実施形態について説明する。第１の実施形態では、計測された夫々の光音響画像の間の比較と、夫々の光音響画像と赤外線カメラ画像との比較と、を共に行って（これらの関係を同時に考慮するコスト関数を用いて）体動を推定し補正した。第２の実施形態では、最初に、計測された夫々の光音響信号の間の比較により第１の推定を行い、次に、第１の推定による補正を行った結果を赤外線カメラ画像と比較することで第２の推定を行う。したがって、最終結果に直結する第２の推定では、赤外線カメラ画像のみを目標とした推定が行われる。そのため、より赤外線カメラ画像に合致した光音響画像を取得できる。

（装置構成）
以下、図９を用いて本実施形態に係る光音響装置の構成について説明する。ただし、第１の実施形態と同様の構成については、同一の番号を付し、説明は省略する。
画像処理装置９０００は、信号取得部１０１０、光音響画像取得部１０２０、広域画像取得部１０３０、計測位置取得部１０４０、比較ペア生成部１０５０、投影画像生成部９０６０を備える。さらに、第１の比較部９０７０、パルス位置推定部９０８０、統合ボリューム生成部９０９０、変形仮説生成部９１１０、変形ボリューム生成部９１２０、第２の比較部９１３０、表示制御部１１００を備える。

投影画像生成部９０６０は、光音響画像取得部１０２０から取得した複数のパルスボリュームから投影画像（ＭＩＰ画像）を夫々生成し（以下、パルス投影画像と呼ぶ）、生成した複数のパルス投影画像を第１の比較部９０７０に出力する。また、変形ボリューム生成部９１２０から取得した複数の変形ボリュームから投影画像を夫々生成し（以下、変形投影画像と呼ぶ）、生成した複数の変形投影画像を第２の比較部９１３０に出力する。

第１の比較部９０７０は、比較ペア生成部１０５０から取得した比較ペア情報に基づいて、ペアとなったパルスボリュームを比較する処理（第１の比較）として、両者のパルス投影画像の画素値等を比較し、ペアの間の一致度に関する情報を算出する。そして、算出した情報をパルス位置推定部９０８０に出力する。
パルス位置推定部９０８０は、第１の比較部９０７０が算出したペア間の一致度に関する情報に基づいて、夫々のパルスボリュームの位置の推定値（パルス位置推定量）を算出して、統合ボリューム生成部９０９０に出力する。

統合ボリューム生成部９０９０は、パルス位置推定量に基づいてパルスボリュームを統合した統合ボリュームを生成して、変形仮説生成部９１１０に出力するとともに、不図示の記憶媒体に保存する。
変形仮説生成部９１１０は、後述する第２の比較部９１３０で算出された推定変形パラメータを基準として、統合ボリューム生成部９０９０から取得した統合ボリュームを変形させるための複数の変形パラメータの仮説（以下、変形仮説と呼ぶ）を生成する。そして、生成した複数の変形仮説を変形ボリューム生成部９１２０と第２の比較部９１３０に出力する。

変形ボリューム生成部９１２０は、変形仮説生成部９１１０から取得した複数の変形仮説ごとに、記憶媒体に保存されている統合ボリュームを変形させて複数の変形ボリュームを生成して、投影画像生成部９０６０に出力する。また、後述する第２の比較部９１３０から変形を完了させる命令情報（変形完了命令）と最終的な変形パラメータ（最終変形パラメータ）を取得した場合、最終変形パラメータに基づいて統合ボリュームを変形させた出力ボリュームを生成する。この出力ボリュームは、表示制御部１１００に出力される。

第２の比較部９１３０は、光音響画像取得部１０２０が生成した複数のパルスボリュームの夫々と、広域画像取得部１０３０から取得した赤外カメラ画像との比較処理（第２の比較）として、パルスボリューム群の統合画像と、赤外カメラ画像とを比較する。より具体的には、赤外線カメラ画像と、投影画像生成部９０６０から取得した複数の変形投影画像（パルスボリューム群の統合画像に由来する）の画素値等を比較し、両者の一致度に関する情報を算出する。さらに、算出した情報に基づくコスト値と、変形仮説生成部９１１０から取得した変形仮説と、に基づいて推定変形パラメータを算出（更新）して、変形仮説生成部９１１０に出力する。また、コスト値の収束状態に基づいて、変形を完了させるか否かを判断するとともに最終変形パラメータを算出する。変形を完了させると判断した場合には、変形完了命令の情報と最終変形パラメータを変形ボリューム生成部９１２０に出力する。

（処理フロー）
図１０は、画像処理装置９０００が行う全体の処理手順を示すフローチャートである。(ステップＳ１００１０〜Ｓ１００７０）
これらのステップにおいては、第１の実施形態のステップＳ５０１０からステップＳ５０７０と同様の処理が実行される。

(ステップＳ１００８０）類似度のピーク位置に基づきパルス毎の体動推定
パルス位置推定部９０８０は、ステップＳ１００７０で算出した類似度ピーク位置PosL_j（1≦j≦N_pair）に基づき、パルス位置推定量Pos_i（1≦i≦N_pulse）を算出する。すなわち、被検体の局所画像（パルスボリューム）間の相対的な位置に関する個別最適値をなるべく保ちつつ、全てのパルスボリュームの位置を全体最適化する。具体的には、コスト関数Eを次式のように定義し、このコスト関数Eを小さくするように各パルスボリュームの位置を最適化する。

関数EL_jは、第１の実施形態のステップＳ５０９０における式１６で定義した関数EL_jと同一である。なお、本実施形態におけるコスト関数Eは、第１の実施形態のステップＳ
５０９０におけるコスト関数E（式１５）の第２項をなくしたものに相当する。コスト関
数Eの最適化の方法は、第１の実施形態のステップＳ５０９０と同様である。

(ステップＳ１００９０）統合ボリューム生成
統合ボリューム生成部１０９０は、パルスボリュームVp_i（1≦i≦N_pulse）を統合し
て統合ボリュームVoを生成する。この生成方法は、第１の実施形態におけるステップＳ５１００と同様である。統合ボリュームは変形仮説生成部９１１０及び変形ボリューム生成部９１２０に出力されると共に、不図示の記憶媒体に保存される。

(ステップＳ１０１００）変形パラメータの仮説生成
変形仮説生成部９１１０は、推定変形パラメータの現在の値を基準として、統合ボリューム生成部９０９０から取得した統合ボリュームVoを変形させるための複数の変形パラメータの仮説（変形仮説）を生成する。ここで、推定変形パラメータの現在の値とは、本ステップを最初に処理する場合には、初期値、すなわち「変形なし」を表す値を取る。一方、本ステップを２回目以降に処理する場合には、後述するステップＳ１０１４０で算出される推定変形パラメータp_eが、推定変形パラメータの現在の値となる。

統合ボリュームの変形は、例えばFree Form Deformation（FFD）によって算出できる。この場合FFDの制御点を、統合ボリュームを包含する領域に対してグリッド状に配置する
。このとき変形パラメータは、FFDの各制御点の制御量で表される。具体的には、制御点
の(x,y,z)の各軸方向の制御量の集合を表す要素数N_grid（N_gridは制御点の数）のベク
トルを夫々、p_x、p_y、p_zとする。そして、それら全てを要素に持つベクトルp={p_x,p_y,p_z}^Tを変形パラメータと定義する。

次に、変形仮説の生成方法を説明する。ここで、現在の推定変形パラメータをp_eとし
たとき、p_eの各要素に対して、夫々微小な変化量を与えることを考える。但し、本ステ
ップを最初に処理する場合は変形仮説生成部９１１０から変形推定パラメータp_eが与え
られないため、全ての要素が０のベクトルp_0をp_eとする。このとき、各要素に与えられた微小な変化量からなるベクトルをΔpと定義する。

次に、変化量Δpを、ベクトルの方向がパラメータ空間中で均等に変化するように変動
させる。例えば制御点の数がN_gridのとき、パラメータ空間は、N_grid×３次元の空間になる。このとき、ある変化量Δpの組を、パラメータ空間中のある軸の要素のみを所定の
値dだけ移動させたベクトルの組として生成する。この処理を全てのパラメータ空間中の
軸（N_grid×３）に対して行うことにより、（N_grid×３）パターンの変化量Δp_i(1≦i≦N_grid×３)を取得できる。つまり、ベクトルの要素数分のΔp_iの組を取得する。なお、変化量Δpの組み合わせは、パラメータ空間中で均等に変化するように変動していれば
、この例に限定されない。

そして、現在の推定変形パラメータp_eに変化量Δp_iを加えた変形パラメータp_iを生
成する（本ステップを最初に処理する場合にはp_i=Δp_iとなる）。ここで、生成された
変形パラメータの数をN_hypoと定義する（本実施形態ではN_hypo=N_grid×３）。このと
き、生成された変形パラメータの集合{p_1，…，p_i，…，p_N_hypo}と、各制御点の位置情報を合わせたものを、変形仮説Hと定義する。そして、変形仮説部９１１０は、生成し
た変形仮説Hを変形ボリューム生成部９１２０に出力するとともに、不図示の記憶媒体に
保存する。なお、変形を表現するモデルはFFDに限られず、Thin Plate Splineなど、非線形な変形を表現する他の方法も利用できる。

(ステップＳ１０１１０）変形ボリュームの生成
変形ボリューム生成部９１２０は、変形仮説生成部９１１０から取得した変形仮説Hに
基づいて、不図示の記憶媒体に保存された統合ボリュームVoを変形させて変形ボリュームDV_iを生成する(1≦i≦N_hypo)。具体的には、FFDの変形アルゴリズムを用いて、変形仮
説Hに含まれる制御点の位置情報と夫々の変形パラメータΔp_iとに基づいて統合ボリュームVoを変形させる。変形ボリュームDV_iは投影画像生成部９０６０に出力される。

(ステップＳ１０１２０）赤外線カメラ方向への変形ボリュームの投影画像生成
投影画像生成部９０６０は、変形ボリューム生成部９１２０が生成した変形ボリュームDV_iから、赤外線カメラの撮像方向ごとに変形投影画像を生成する。これは、ステップＳ１００６０におけるパルスボリュームを、変形ボリュームDV_iに入れ替えただけの処理である。本ステップにおける投影方法は、第１の実施形態の「変形例１−１」のような透視投影法とする。装置座標系C_DEVのＸ、Ｙ、Ｚ方向の夫々での投影により生成された変形投影画像を夫々、DIpx_i、DIpy_i、DIz_i（1≦i≦N_hypo）と表記する。

(ステップＳ１０１３０）変形投影画像と赤外線画像の類似度算出
第２の比較部９１３０は、変形投影画像DIpx_i、DIpy_i、DIz_i（1≦i≦N_hypo）と、
夫々が対応する投影方向の赤外線カメラ画像I_CAM2、I_CAM3、I_CAM1との類似度を算出する
。算出された類似度を夫々Sx_i、Sy_i、Sz_iと表記する。ここでは第１の実施形態のステップＳ５０７０と同様の類似度尺度を適用できる。

(ステップＳ１０１４０）類似度に基づき変形パラメータ推定
第２の比較部９１３０は、ステップＳ１０１３０で算出された類似度Sx_i、Sy_i、Sz_i（1≦i≦N_hypo）に基づいて、推定変形パラメータを算出（更新）する。具体的には、まず、類似度Sx_i、Sy_i、Sz_iの夫々に基づいて、変形パラメータのi番目の要素を微小変
動させたときの（すなわち、p_iを仮定したときの）変形投影画像と赤外線カメラ画像と
の間の一致度のコスト関数E_sの値を算出する。コスト関数E_sは、例えばSx_i、Sy_i、Sz_iの３つの類似度の平均値として定義できる。コスト関数の定義方法はこれに限らず、例えば３つの類似度の中の最大値を採用してもよい。なお、コスト関数E_sは、変形パラメ
ータpを与えたときの類似度を表すため、pの関数としてE_s（p）と表記できる。従って、変形パラメータp=p_iに対応する（類似度Sx_i、Sy_i、Sz_iに基づく）コスト関数の値はE_s（p_i）と表記される。

次に、従来の推定変形パラメータp_eと、従来の推定変形パラメータのコスト値E_s（p_e）と、上記で算出した夫々の仮説のコスト値E_s（p_i）に基づいて、新たな推定変形パ
ラメータp_e^newを推定する。具体的には、パラメータp_e^newを、以下の式によって算出する。

ここで、αは小さな正の定数であり、本実施例では例えばα＝０．１を採用する。また、grad E_sはE_sの勾配ベクトルを表す。ここで、式２１の右辺のベクトルのi番目の要素は、コスト関数E_s（p）をベクトルpの要素iで偏微分した値を表す。より具体的には、その偏微分した値を、仮説のコスト値E_s（p_i）と従来のパラメータのコスト値E_s（p_e）の間の差分で近似的に求めた値である。これにより、パラメータ空間において変形パラメータを微小量だけ変動させたときに、コスト関数E_sの傾きがが最も大きくなる方向に変
動させた変形パラメータp_e^newを取得できる。

次に、上記で求めた推定変形パラメータp_e^newに対してステップＳ１０１１０〜Ｓ１０１３０と同様の処理を施すことで、パラメータp_e^newに基づく変形ボリュームの投影画像と赤外線カメラ画像との類似度を得る。そして、該類似度に基づくコスト値E_s（p_e^new
）を算出して、不図示の記憶媒体に保存する。このとき、本ステップが実行されて新しいコスト値E_s（p_e^new）が算出される度に、過去のコスト値を保持したまま、該コスト値
が不図示の記憶媒体に追加格納されるものとする。なお、上述した式２１の演算で使用するコスト値E_s（p_e）は、１回前に本ステップの処理を実行した際のコスト値E_s（p_e^new）であるため、不図示の記憶媒体から取得できる。ただし、本ステップの処理を最初に
行う場合には、コスト値E_s（p_e）が保持されていないので、式２１の演算に先立ち初期値p_0に関するコスト値E_s（p_0）を算出する必要がある。

以上のように、本ステップでは、コスト関数E_sの勾配方向に移動する変形パラメータ
を見つける処理を行っている。これは、最急降下法を用いてコスト関数E_sを最適化する
際の１ステップに相当する。従って、本ステップが繰り返されることで、推定変形パラメータp_e^newがより好ましい値になる。なお、コスト関数E_sの最適化方法はこれに限られ
ず、例えばニュートン法などでも良い。

(ステップＳ１０１５０）コスト関数が収束したか判定
第２の比較部９１３０は、不図示の記憶媒体に保存された歴代のコスト値E_s（p_e）を参照し、値が収束しているか否かを判定する。収束していれば、変形完了命令の情報とともに推定変形パラメータp_eを最終変形パラメータp_finalとして変形ボリューム生成部９１２０に出力し、処理をステップＳ１０１６０に移行する。収束していなければ、推定変形パラメータp_eを変形仮説生成部９１１０に出力し、ステップＳ１０１００に戻る。

収束の判定方法の例を挙げる。まず、コスト値E_s（p_e）の最新の値と一つ前の値を比較し、差分が所定の値よりも小さくなった場合に収束したと判定する方法がある。また、過去の所定の回数分のコスト値を抽出し、隣り合う２回分の変化量を夫々計算し、変化量の平均値が所定の値よりも小さくなった場合に収束したと判定する方法がある。また簡易的には、ステップＳ１０１４０を繰り返した回数が所定の値に達した時に収束したと判定する方法もある。

(ステップＳ１０１６０）最終変形ボリューム生成
変形ボリューム生成部９１２０は、第２の比較部９１３０から変形完了命令の情報と、最終変形パラメータp_finalを取得し、出力ボリュームV_finalを生成する。
以上のように本実施形態では、ステップＳ１０１００からＳ１０１５０を、赤外線カメラ画像に合致させることを目標とするコスト関数のコスト値が収束するまで繰り返した後に、ステップＳ１０１６０で最終変形ボリュームを生成する。これにより、コスト関数を最適にする変形パラメータを取得できる。

(ステップＳ１０１７０）表示
表示制御部１１００は、出力ボリュームV_finalの情報を表示装置１３０に出力して表
示する制御を行う。本処理は、第１の実施形態のステップＳ５１１０と同様である。

以上、本実施形態の画像処理装置によれば、被検体の動き補正の最終結果に直結する推定において、赤外線カメラ画像のみを目標としてパルスボリュームの変形を推定することで、より赤外線カメラ画像に写る被検体に合致した光音響画像を取得できる。

＜第３の実施形態＞
本発明の第３の実施形態について説明する。第１および第２の実施形態では、赤外線カメラ画像を広域画像として利用した。第３の実施形態では、広域画像として被検体の距離画像を利用する。これにより、ボリュームデータである光音響画像の情報を縮退させた２次元空間上ではなく、３次元空間上での情報の比較が可能になる。

（装置構成）
以下、図１１を用いて本実施形態に係る光音響装置の構成について説明する。第１の実施形態と同様の構成については、同一の番号を付し、説明は省略する。
距離画像カメラ１４０は、レーザを被検体に投射し、投射したレーザが被検体から反射して往復してくる時間に基づいて被検体までの距離を計測する（Time Of Flight）カメラである。これにより被検体の３次元表面形状が距離画像の形式で得られる。距離画像カメラ１４０は、被検体全体を計測できる位置に設置される。得られた距離画像(I_RNGと表記)は画像処理装置１１０００に入力される。

なお、距離画像カメラ１４０は装置座標系C_DEVにおいて校正済みであり、座標変換の情報やカメラの内部パラメータは、既知の情報として画像処理装置１１０００が保持している。距離画像I_RNG上の各画素値は、夫々の画素を通る視線上に存在する空間中の点までの距離を表す。すなわち、各画素の座標値と画素値から、３次元の距離画像カメラ座標系C_RNGにおける各画素に対応する点の座標を求めることができる。また、距離画像カメラ座標系C_RNGから装置座標系C_DEVへの変換（以降、T_RtoDと呼ぶ）には、既知の手法を利用でき
る。

画像処理装置１１０００は、信号取得部１０１０、光音響画像取得部１０２０、広域表面形状取得１１０３０、計測位置取得部１０４０、比較ペア生成部１０５０、表面形状抽出部１１０６０を備える。さらに、比較部１１０７０、パルス位置推定部１１０８０、統合ボリューム生成部１０９０、表示制御部１１００を備える。

広域表面形状取得１１０３０は、距離画像カメラ１４０が撮像した距離画像を広域画像として取得し、該距離画像から被検体全体の広域の表面形状を示す３次元点群情報を生成し（以下、広域表面形状と呼ぶ）、比較部１１７０に出力する。
局所表面形状抽出部１１０６０は、被検体の局所領域であるＦＯＶを描出したパルスボリュームから被検体の表面形状を抽出し（以下、局所表面形状と呼ぶ）、比較部１１０７０に出力する。

比較部１１０７０は、比較ペア情報に基づき、ペアとなったパルスボリュームを比較する処理（第１の比較）として、ペア間のボクセル値等を比較する処理を行い、ペア間の一
致度に関する情報として、ペアの間の類似度等を算出する。また、夫々のパルスボリュームと距離画像とを比較する処理（第２の比較）として、表面形状抽出部１１０６０から取得した局所表面形状と広域表面形状取得部１１０３０から取得した広域表面形状を比較し、一致度に関する情報を算出する。さらに、算出した一致度に関する情報をパルス位置推定部１１０８０に出力する。

パルス位置推定部１１０８０は、比較部１１０７０から入力されたペア間の一致度に関する情報と、夫々のパルスボリュームと距離画像の一致度に関する情報と、の両方に基づいて、夫々のパルスボリュームの位置の推定値（パルス位置推定量）を算出する。算出された推定値は、統合ボリューム生成部１０９０に出力される。

（処理フロー）
図１２は、画像処理装置１１０００が行う全体の処理手順を示すフローチャートである。
(ステップＳ１２０１０からＳ１２０２０）
これらのステップにおいては、第１の実施形態におけるステップＳ５０１０からステップＳ５０２０と同様の処理が行われる。

(ステップＳ１２０３０）広域表面形状取得
広域表面形状取得１１０３０は、距離画像カメラ１４０が撮像した被検体の距離画像（I_RNG）を取得する。ここでは距離画像として、被検体のある瞬間を写した静止画像を取得する。被検体の位置・姿勢が、光音響信号計測装置１１０による計測時とほぼ同じなら、静止画像の取得時点は問わない。なお、距離画像として動画像を取得し、画像の変化が小さいフレームを抽出しても良い。次に、広域表面形状取得１１０３０は、任意の既知の手法によって、距離画像I_RNGから被検体全体の表面形状を計測した広域表面形状（以下、Surf_Gと呼ぶ）を生成する。

(ステップＳ１２０４０からＳ１２０５０）
これらのステップにおいては、第１の実施形態におけるステップＳ５０４０からステップＳ５０５０と同様の処理が行われる。

(ステップＳ１２０６０）局所表面形状取得
表面形状取得部１１０６０は、パルスボリュームVi（1≦i≦N_pulse）の夫々について
、被検体の局所表面形状を抽出する。

ここでは表面形状を装置座標系C_DEVにおける３次元点群の情報として取得し、局所表面形状SurfL_i（1≦i≦N_pulse）として保存する。なお、局所表面形状の抽出は、体表をセグメンテーションする既知の技術を利用できる。また、局所表面形状を手動で抽出してもよい。例えば、操作者が、表示装置１３０上に表示されたパルスボリュームの断面画像から位置を指定する方法が考えられる。この場合、操作者より指定された３次元点群に補間処理を施すことが好ましい。

(ステップＳ１２０７０）パルスボリューム間の類似度のピーク位置取得
比較部１１０７０は、比較ペア情報R_j（1≦j≦N_pair）が表す全てのペアに対して、
比較ペアとされたパルスボリュームを比較し、ペア間の一致度に関する情報を取得する。具体的には、両パルスボリュームの類似度関数と、類似度関数のピーク位置を取得する。すなわち、下式のような、比較するパルスボリューム間の相対位置を(x，y，z)だけずら
しながら比較した場合のボリューム間の類似度関数FL_j(x，y，z)を算出する。

f_simil(I1，I2，x，y，z)は、画像I1に対して画像I2の相対位置を(x，y，z)だけ並進
させた場合の画像間の類似度を算出する関数である。式２２は、第１の実施形態のステップＳ５０７０における式１に対して、並進量の次元を１つ増やしたものであるので、詳細の説明は省略する。本実施形態では、比較ペア情報の夫々について、各類似度関数 (FL_j：1≦j≦N_pair)を３次元の連続関数として取得する。

次に、比較部１１０７０は、比較ペアの夫々に関する類似度関数が最大となる位置（類似度ピーク位置）を算出する。具体的には、上で算出した各比較ペアに関する類似度関数FL_j：1≦j≦N_pair)について、その関数値が最大となる位置PosL_jを次式で算出する。

以上の処理により、比較部１１０７０は、全てのペアに関する類似度ピーク位置PosL_j（1≦j≦N_pair)を算出する。

(ステップＳ１２０８０）局所表面形状と広域表面形状間の一致度のピーク位置取得
比較部１１０７０は、パルスボリュームの夫々と距離画像を比較し、パルスボリュームの夫々と距離画像との一致度に関する情報を取得する。具体的には、Ｓ１２０６０で得た局所表面形状SurfL_i（1≦i≦N_pulse）の夫々について、Ｓ１２０３０で得た広域表面形状SurfGとの間の一致度関数と、一致度関数のピーク位置を取得する。本ステップの目的
は、広域表面形状と各パルスボリュームの局所表面形状の類似度を計算することで、被検体の広域画像と局所画像を比較することである。各パルスボリュームに対応する夫々の局所表面形状SurfL_iについて、以下の一致度関数を取得する。

ここで、f_dist(Surf1，Surf2，x，y，z)は、表面形状Surf1に対して表面形状Surf2の
相対位置を(x，y，z)だけずらしながら比較した場合の形状間の一致度を算出する関数で
ある。本実施形態では表面形状は点群で表現されるため、点群間で互いに最近傍となる点同士を対応付け、対応付けられた点群間の平均距離を適用すれば良い。但し、一致度の算出方法はこれに限られない。例えば、表面形状を表す点群から三角パッチなどによるメッシュをはり、メッシュを構成する面同士の距離を比較して一致度を算出する方法もある。本実施形態では、パルスボリュームの夫々について、一致度関数(FG_i：1≦i≦N_pulse)
を３次元の連続関数として取得する。

次に、比較部１１０７０は、パルスボリュームの夫々に関して、距離画像に対する一致度関数が最大となる位置（一致度ピーク位置）を算出する。具体的には次式により、各パルスボリュームに関する類似度関数(FG_i：1≦i≦N_pulse)の夫々について、関数値が最
大となる位置PosG_iを算出する。

以上の処理により、比較部１１０７０は、全てのパルスボリュームに関して、距離画像に対する一致度のピーク位置PosG_i（1≦i≦N_pulse）を算出する。

(ステップＳ１２０９０）２種類のピーク位置に基づきパルス毎の体動推定
比較部１１０７０は、Ｓ１２０７０で得た類似度ピーク位置PosL_j（1≦j≦N_pair）と、Ｓ１２０８０で得た一致度ピーク位置PosG_i（1≦i≦N_pulse）からパルス位置推定量Pos_i（1≦i≦N_pulse）を求める。すなわち、被検体の局所画像（パルスボリューム）間
の相対的な位置に関する個別最適値と、被検体の広域画像（距離画像）に対する各局所画像の位置に関する個別最適値と、の両方をなるべく保ちつつ、全てのパルスボリュームの位置を全体最適化する。これは第１の実施形態のステップＳ５０９０と同様の処理である。その際、Ｓ５０９０における投影画像間の類似度ピーク位置をパルスボリューム間の類似度ピーク位置に、赤外線カメラ画像に対する投影画像の類似度ピーク位置を広域表面形状に対する局所表面形状の一致度ピーク位置に、夫々置き換える。

(ステップＳ１２１００からＳ１２１１０）
これらのステップにおいては、第１の実施形態のステップＳ５１００からステップＳ５１１０と同様の処理が行われる。

以上、本実施形態の画像処理装置によれば、被検体の距離画像から被検体の３次元表面形状を取得することで、３次元空間上での情報のみの比較に基づいて光音響画像を補正できる。これにより、比較する情報を低次元空間に縮退することで情報を落とすことなく、被検体の動きを補正できる。

（変形例３−１）広域表面形状は超音波画像から取得してもよい
上記の第３の実施形態では、ステップＳ１１０３０の処理として、距離画像から被検体の広域表面形状の抽出を行った。しかし、広域表面形状の取得方法はこれに限られない。例えば、被検体を超音波測定して得られた超音波画像から抽出してもよい。その場合、広域表面形状取得部１１０３０は、不図示の超音波撮像装置を用いて被検体全体の超音波ボリュームを取得し、広域表面形状を抽出する。この場合例えば、被検体の内部と外部での音響インピーダンスの違いに由来する境界を閾値処理によって抽出することで、被検体の広域表面形状を取得できる。これにより、例えば被検体が水に浸かっている場合は、距離画像カメラでは撮像条件が悪いのに対し、水と被検体（人体など）の間では音響インピーダンスの違いが大きいため、精度良く広域表面形状を取得できる。

（変形例３−２）パルスボリューム間の比較は局所表面形状の一致度を算出してもよい
上記の第３の実施形態では、ステップＳ１２０７０の処理として、パルスボリューム間の画像類似度を一致度として、そのピーク位置を算出していた。しかし、パルスボリューム間の一致度を算出する方法は他の何れの方法であってもよい。例えば、画像処理によって夫々のパルスボリュームから血管分岐等の解剖学的な特徴を抽出し、その位置や分布の一致度をパルスボリューム間の一致度として用いてもよい。あるいは、パルスボリューム間の局所表面形状の一致度を算出して、これをパルスボリューム間の一致度として用いてもよい。これは、ステップＳ１２０８０の広域表面形状と局所表面形状の一致度関数を算出する処理において、広域表面形状ではなく、異なる２つのパルスボリュームの局所表面形状を採用することで算出できるため、説明を省略する。これにより、３次元空間中に敷き詰められたボクセル単位での比較を要する画像類似度の計算よりも、３次元空間の限られた領域の比較しか必要としないため、計算コストを削減できる。

＜第４の実施形態＞
第４の実施形態における画像処理装置は、被検体の広域画像として、赤外線カメラ画像でなく、被検体の外観のみが描出されたカメラ画像を利用する。以下、上記各実施形態と
異なる点を中心に説明する。

（装置構成）
以下、図１３を用いて本実施形態に係る光音響装置の構成について説明する。ただし、上記実施形態と同様の構成については、同一の番号を付し、説明は省略する。
カメラ１５０は、被検体の外観を（可視光領域の光を）撮影するカメラであり、被検体全体の外観を撮影可能な位置に設置されている。カメラ１５０は第１の実施形態における赤外線カメラ１２０と違ってカラー画像を撮影する。それ以外の点では第１の実施形態における赤外線カメラ１２０の説明における「赤外線カメラ」の部分を「カメラ」に置き換えたものと考えて良い。

画像処理装置１０００は、信号取得部１０１０、光音響画像取得部１０２０、カメラ画像取得部１３０３０、計測位置取得部１０４０、比較ペア生成部１０５０、投影画像生成部１０６０、局所外形形状抽出部１３１２０を備える。さらに、比較部１３０７０、パルス位置推定部１３０８０、統合ボリューム生成部１０９０、表示制御部１１００を備える。

カメラ画像取得部１３０３０は、カメラ１５０が撮像したカメラ画像（２次元画像）を広域画像として取得し、広域外形形状抽出部１３１１０に出力する。
広域外形形状抽出部１３０３５は、カメラ画像から被検体の広域の２次元外形形状を抽出し、抽出した外形形状（以下、広域外形形状と呼ぶ）を比較部１３０７０に出力する。

局所外形形状抽出部１３０６５は、被検体の局所領域の投影画像から被検体の局所的な２次元外形形状を抽出し、抽出した外形形状（以下、局所外形形状と呼ぶ）を比較部１３０７０に出力する。
比較部１３０７０は、比較ペア情報に基づき、ペアとなったパルスボリュームを比較する処理（第１の比較）として、第１の実施形態の比較部１０７０と同様の処理を行う。また、夫々のパルスボリュームとカメラ画像とを比較する処理（第２の比較）を行う。具体的には、局所外形形状抽出部１３０６５から取得した局所外形形状と広域外形形状取得部１３０３５から取得した広域外形形状を比較し、局所外形形状と広域外形形状の間の一致度に関する情報を算出し、パルス位置推定部１３０８０に出力する。

パルス位置推定部１３０８０は、比較部が算出したペア間の一致度に関する情報と、夫々のパルスボリュームとカメラ画像の一致度に関する情報と、の両方に基づいて、夫々のパルスボリュームの位置の推定値（パルス位置推定量）を算出する。算出された推定値は、統合ボリューム生成部１０９０に出力される。

（処理フロー）
図１４は、画像処理装置１３０００が行う全体の処理手順を示すフローチャートである。
（ステップＳ１４０１０からＳ１４０２０）
これらのステップにおいては、第１の実施形態におけるステップＳ５０１０からステップＳ５０２０と同様の処理が実行される。

（ステップＳ１４０３０）
カメラ画像取得部１３０３０は、カメラ３０１、３０２、３０３が撮影したカメラ画像（I_CAM1，I_CAM2，及びI_CAM3）を取得する。この処理は、第１の実施形態におけるステッ
プＳ５０３０の「赤外線カメラ」を「カメラ」に置き換えたものと同様である。

（ステップＳ１４０３５）広域外形形状抽出
広域外形形状抽出部１３０３５は、カメラ画像（I_CAM1，I_CAM2，及びI_CAM3）の夫々か
ら被検体の外形形状を抽出する。外形形状の抽出には例えば、一般的なエッジ検出手法を使用できる。抽出した広域外形形状を夫々、Surf_G1、Surf_G2、Surf_G3と定義する。

（ステップＳ１４０４０からＳ１４０５０）
これらのステップにおいては、第１の実施形態におけるステップＳ５０４０からステップＳ５０５０と同様の処理が行われる。

（ステップＳ１４０６５）局所外形形状抽出
局所外形形状抽出部１３０６５は、投影画像の夫々から被検体の局所外形形状を抽出する。投影画像はパルスボリュームの投影画像であるため、写っている構造情報はパルスボリュームと変わらない。従って本ステップは、第３の実施形態のステップＳ１２０６０の処理を２次元領域において行う場合と同様に実行できる。抽出した局所外形形状を夫々、SurfLx_i、SurfLy_i、SurfLz_i（1≦i≦N_pulse）と表記する。
（ステップＳ１４０７０）
本ステップでは、第１の実施形態におけるステップＳ５０７０と同様の処理が行われる。

（ステップＳ１４０８０）局所外形形状と広域外形形状の一致度のピーク位置取得
比較部１３０７０は、パルスボリュームの夫々とカメラ画像を比較し、パルスボリュームの夫々とカメラ画像との一致度に関する情報を取得する。具体的には、局所外形形状SurfLx_i、SurfLy_i、SurfLz_iの夫々について、ステップＳ１４０３５で抽出された対応する広域外形形状Surf_G2、Surf_G3、Surf_G1との間で、形状の一致度関数を取得する。そ
して、取得した一致度関数のピーク位置を取得する。本ステップの目的は、カメラ画像から生成した広域外形形状に対する各パルスボリュームの局所外形形状の類似度を計算することで、被検体の広域画像と局所画像を比較することである。

一致度関数取得の際には、各方向の局所外形形状について、第３の実施形態のステップＳ１２０８０の式２４と同様に、関数FGx_i、FGy_i、FGz_iを取得する。但し、式２４が
３次元座標上で局所表面形状を並進させるのに対し、本実施形態では２次元座標上で局所外形形状を並進させる点で異なる。具体的な２次元座標は、FGx_iは（y,z）、FGy_iは（z,x）、FGz_iは（x,y）である。また、一致度関数のピーク位置PosG_i（1≦i≦N_pulse）
を算出する際は、第１の実施形態のステップＳ５０８０における各軸方向の類似度関数FGX_i，FGY_i，FGZ_iを一致度関数に置き換えれば良い。

（ステップＳ１４０９０）２種類のピーク位置に基づきパルス毎の体動推定
比較部１３０７０は、Ｓ１４０７０で得た類似度ピーク位置PosL_j（1≦j≦N_pair）とＳ１４０８０で得た一致度ピーク位置PosG_i（1≦i≦N_pulse）から、パルス位置推定量Pos_i（1≦i≦N_pulse）を算出する。すなわち、被検体の局所画像（投影画像）間の相対
的な位置に関する個別最適値と、被検体の広域画像（カメラ画像）に対する各局所画像の位置に関する個別最適値と、の両方をなるべく保ちつつ、全てのパルスボリュームの位置の全体最適化を行う。これによりパルス位置推定量Pos_i（1≦i≦N_pulse）を算出する。本ステップの処理は、第１の実施形態のステップＳ５０９０における赤外線カメラ画像に対する投影画像の類似度ピーク位置を、広域外形形状に対する局所外形形状の一致度ピーク位置に置き換えたものである。

以上、本実施形態によれば、被検体のカメラ画像から被検体の外形形状を取得し、光音響画像から抽出した形状と比較することで、光音響画像を補正できる。これにより、入手可能性やコスト面で有利な汎用的なカメラを利用して被検体の動きを補正できる。

１５０２：光源，１５０９：受信素子，１０００：画像処理装置，１５０８：変化手段，１０３０：広域画像取得部

Claims (14)

1. 光源と、
   前記光源から光を照射された被検体から発生する音響波を受信する受信素子と、
   前記音響波を用いて、前記被検体内部の特性情報を示す画像データを生成する処理手段と、
   前記被検体における前記光の照射位置を変化させる変化手段と、
   前記被検体の広域画像を取得する広域画像取得手段と、
   を有し、
   前記処理手段は、
   前記変化手段により変化する前記照射位置ごとに、当該照射位置に対応する前記被検体の局所領域における前記画像データである局所画像を生成し、
   前記照射位置ごとに得られた複数の前記局所画像の間の比較である第１の比較と、前記複数の局所画像のそれぞれと前記広域画像との比較である第２の比較と、に基づいて前記複数の局所画像を統合する
   ことを特徴とする光音響装置。
2. 前記処理手段は、前記局所領域ごとに、前記音響波を用いて、前記局所画像としてボリュームデータを生成する
   ことを特徴とする請求項１に記載の光音響装置。
3. 前記処理手段は、前記局所画像の投影画像を生成し、前記投影画像の間の一致度に関する情報を取得することにより、前記第１の比較を行う
   ことを特徴とする請求項２に記載の光音響装置。
4. 前記処理手段は、互いにオーバーラップする部分を有する前記局所画像をペアとし、前記第１の比較として、当該ペアとなった局所画像を比較する
   ことを特徴とする請求項１に記載の光音響装置。
5. 前記処理手段は、前記広域画像と前記局所画像の間の一致度が最大になる前記局所画像の位置を取得することにより、前記第２の比較を行う
   ことを特徴とする請求項１ないし４のいずれか１項に記載の光音響装置。
6. 前記処理手段は、前記第１の比較による前記複数の局所画像の間の一致度に関する情報と、前記第２の比較による前記複数の局所画像のそれぞれと前記広域画像との一致度に関する情報と、を用いて前記統合を行う
   ことを特徴とする請求項１ないし５のいずれか１項に記載の光音響装置。
7. 前記処理手段は、前記第１の比較により前記複数の局所画像の統合を行ったのち、前記第２の比較により前記統合した画像を変形させる
   ことを特徴とする請求項１ないし５のいずれか１項に記載の光音響装置。
8. 前記広域画像取得手段は、前記被検体の赤外線画像を取得する赤外線カメラである
   ことを特徴とする請求項１ないし７のいずれか１項に記載の光音響装置。
9. 前記広域画像取得手段は、前記広域画像として、前記被検体の距離画像に基づく３次元表面形状を取得する
   ことを特徴とする請求項１ないし７のいずれか１項に記載の光音響装置。
10. 前記広域画像取得手段は、前記被検体の外観が描出された画像を取得するカメラである
    ことを特徴とする請求項１ないし７のいずれか１項に記載の光音響装置。
11. 前記処理手段は、前記複数の局所画像のそれぞれの投影画像を生成し、前記投影画像と前記広域画像の間の一致度に関する情報を取得することにより、前記第２の比較を行う
    ことを特徴とする請求項２または請求項３に記載の光音響装置。
12. 前記広域画像取得手段はカメラであり、
    前記処理手段は、前記カメラの撮像方向と略同一の方向に前記局所画像を投影することで前記投影画像を生成する
    ことを特徴とする請求項１１に記載の光音響装置。
13. 前記処理手段は、前記複数の局所画像を統合する際に、前記第１の比較および前記第２の比較の結果に基づいて、前記複数の局所画像を生成するための音響波が取得される間に前記被検体の体動があったことによる前記複数の局所画像それぞれの位置ズレを補正することを特徴とする請求項１ないし５のいずれか１項に記載の光音響装置。
14. 照射位置を変化させながら光を照射された被検体から発生する音響波を用いて、前記被検体内部の特性情報を示す画像データを生成する画像処理方法であって、
    前記被検体の広域画像を取得するステップと、
    変化する前記照射位置ごとに、当該照射位置に対応する前記被検体の局所領域における前記画像データである局所画像を生成するステップと、
    前記照射位置ごとに得られた複数の前記局所画像の間の比較と、前記複数の局所画像のそれぞれと前記広域画像との比較と、に基づいて前記複数の局所画像を統合するステップと、
    を有することを特徴とする画像処理方法。

Images