JP6436442B2 - 光音響装置および画像処理方法 - Google Patents
光音響装置および画像処理方法 Download PDFInfo
- Publication number
- JP6436442B2 JP6436442B2 JP2015080675A JP2015080675A JP6436442B2 JP 6436442 B2 JP6436442 B2 JP 6436442B2 JP 2015080675 A JP2015080675 A JP 2015080675A JP 2015080675 A JP2015080675 A JP 2015080675A JP 6436442 B2 JP6436442 B2 JP 6436442B2
- Authority
- JP
- Japan
- Prior art keywords
- image
- comparison
- subject
- local
- pulse
- 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.)
- Active
Links
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/0093—Detecting, measuring or recording by applying one single type of energy and measuring its conversion into another type of energy
- A61B5/0095—Detecting, measuring or recording by applying one single type of energy and measuring its conversion into another type of energy by applying light and detecting acoustic waves, i.e. photoacoustic measurements
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/0033—Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room
- A61B5/0035—Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room adapted for acquisition of images from more than one imaging mode, e.g. combining MRI and optical tomography
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
- A61B5/7207—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal of noise induced by motion artifacts
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/103—Detecting, measuring or recording devices for testing the shape, pattern, colour, size or movement of the body or parts thereof, for diagnostic purposes
- A61B5/11—Measuring movement of the entire body or parts thereof, e.g. head or hand tremor, mobility of a limb
- A61B5/1126—Measuring movement of the entire body or parts thereof, e.g. head or hand tremor, mobility of a limb using a particular sensing technique
- A61B5/1128—Measuring movement of the entire body or parts thereof, e.g. head or hand tremor, mobility of a limb using a particular sensing technique using image analysis
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/145—Measuring characteristics of blood in vivo, e.g. gas concentration, pH value; Measuring characteristics of body fluids or tissues, e.g. interstitial fluid, cerebral tissue
- A61B5/1455—Measuring characteristics of blood in vivo, e.g. gas concentration, pH value; Measuring characteristics of body fluids or tissues, e.g. interstitial fluid, cerebral tissue using optical sensors, e.g. spectral photometrical oximeters
- A61B5/14551—Measuring characteristics of blood in vivo, e.g. gas concentration, pH value; Measuring characteristics of body fluids or tissues, e.g. interstitial fluid, cerebral tissue using optical sensors, e.g. spectral photometrical oximeters for measuring blood gases
Description
光源と、
前記光源から光を照射された被検体から発生する音響波を受信する受信素子と、
前記音響波を用いて、前記被検体内部の特性情報を示す画像データを生成する処理手段と、
前記被検体における前記光の照射位置を変化させる変化手段と、
前記被検体の広域画像を取得する広域画像取得手段と、
を有し、
前記処理手段は、
前記変化手段により変化する前記照射位置ごとに、当該照射位置に対応する前記被検体の局所領域における前記画像データである局所画像を生成し、
前記照射位置ごとに得られた複数の前記局所画像の間の比較である第1の比較と、前記複数の局所画像のそれぞれと前記広域画像との比較である第2の比較と、に基づいて前記複数の局所画像を統合する
ことを特徴とする光音響装置である。
照射位置を変化させながら光を照射された被検体から発生する音響波を用いて、前記被検体内部の特性情報を示す画像データを生成する画像処理方法であって、
前記被検体の広域画像を取得するステップと、
変化する前記照射位置ごとに、当該照射位置に対応する前記被検体の局所領域における前記画像データである局所画像を生成するステップと、
前記照射位置ごとに得られた複数の前記局所画像の間の比較と、前記複数の局所画像のそれぞれと前記広域画像との比較と、に基づいて前記複数の局所画像を統合するステップと、
を有することを特徴とする画像処理方法である。
する技術に関する。よって本発明は、被検体情報取得装置またはその制御方法、あるいは被検体情報取得方法や信号処理方法として捉えられる。本発明はまた、これらの方法をCPU等のハードウェア資源を備える情報処理装置に実行させるプログラムや、そのプログラムを格納した記憶媒体としても捉えられる。本発明はまた、被検体内部の特性情報を補正する装置や、その制御方法、特性情報の補正方法、補正プログラムとも捉えられる。本発明はまた、被検体内部の特性情報を示す画像データを処理する機能を有することから、画像処理装置、画像処理方法および画像処理プログラムとも捉えられる。
本実施形態に係る光音響装置は、複数回に渡って計測した一連の光音響信号を統合した画像を生成する。光音響装置は、計測中に生じる被検体の体動等の動きを推定して、その動きの影響を低減するように信号を補正する。より具体的には、光音響装置は、FOVが異なる複数の光音響画像(局所画像)同士を比較すると同時に、各光音響画像(局所画像)を、被検体全体を撮像した赤外線カメラ画像(広域画像)と比較する。これらの比較結果を用いた動き推定により、本来の被検体の解剖構造に対して全体として歪みが抑制された音響画像を取得できる。
図1は、本実施形態に係る光音響装置の構成を示す。同図に示すように、本実施形態における光音響装置は、画像処理装置1000、光音響信号計測装置110、赤外線カメラ120、および表示装置130を有する。画像処理装置1000は、光音響信号計測装置110、赤外線カメラ120、および表示装置130と接続されている。
組とすると、パルス光の照射回数分だけの光音響信号の組が取得できる。計測した光音響信号は、画像処理装置1000へと入力される。
図2は、光音響信号計測装置が被検体を撮像可能な範囲(FOV)を示す。FOV202は、被検体200に対して照射方向201の方向にパルス光を照射した時に撮像可能な範囲である。本実施形態の照射方向201は、光音響信号計測装置110の装置座標系CDEVのZ軸に平行である。図2のように、1つのFOVでは被検体200全体をカバーできない場合、光照射ごとにFOVを少しずつずらしながら被検体全体から光音響信号を取得する必要がある。そのためには、光の射出端を可動式にする、光照射方向を変化可能にするなどの方法がある。このような構成によって、関心領域が1つのFOVよりも広い場合に対応できる。なお、後の画像統合で利用するために、パルス照射ごとにFOVの位置を記録することが好ましい。位置記録方法として例えば、光照射位置などを示す基準座標を記録する方法がある。
図15に、光音響信号計測装置110の構成例を示す。測定対象となる被検体1500は、例えば生体の乳房である。光音響信号計測装置110は、画像処理装置1000と通信して、制御情報の受信や取得信号の送信を行う。光音響信号計測装置は、本発明の光音響装置に相当すると考えることもできるし、光音響信号計測装置と画像処理装置を合わせて光音響装置と考えてもよい。
、被検体1500と射出端との相対位置を変化させることにより、FOVの位置が少しずつ移動し、被検体全域が計測される。変化手段は、探触子を射出端と共に移動させても良い。変化手段として例えば、モータや制御ステージなどの駆動部品が好適である。変化手段は、被検体の方を移動させても良い。あるいは変化手段は、光の照射方向を変化させることで光の照射位置、ひいてはFOVを変化させても良い。
赤外線カメラ120は、被検体の外観及び被検体の体表近傍の血管の様子を静止画像または動画像として撮影し、画像処理装置1000に入力する。赤外線カメラ120は、被検体全体の外観を撮影可能な位置に設置される。本実施形態では3台の赤外線カメラを利用する。3台を区別する必要がある場合、それぞれ符号301,302,303と表す。
座標系CDEVの−Z方向を向いている)。赤外線カメラ302の撮像方向305は、装置座標系CDEVの−X方向に一致する(すなわち、カメラ座標系CCAM2は、Z軸が装置座標系CDEVのX方向を向いている)。赤外線カメラ303の撮像方向306は、装置座標系CDEVの
Y方向に一致する(すなわち、カメラ座標系CCAM3は、Z軸が装置座標系CDEVの−Y方向
を向いている)。ここで、カメラ座標系CCAM1、CCAM2、CCAM3から装置座標系CDEVへの座
標変換を、それぞれ、TC1toD、TC2toD、TC3toDと定義する。各カメラは装置座標系CDEVにおいて校正済みであり、上記の座標変換の情報や夫々のカメラの内部パラメータは、予め画像処理装置1000が保持している。
3次元空間中の投影面上の点と、を通る直線)と一対一の関係を持つ。このようなカメラ画像座標系とカメラ座標系の間の変換に関しては、一般的な座標変換方法を利用できる。カメラ座標系CCAM1、CCAM2、CCAM3からカメラ画像座標系CIMG1、CIMG2、CIMG3への変換を、それぞれ、TC1toI1、TC2toI2、TC3toI3と定義する。
表示装置120は、液晶装置やCRTなどのディスプレイである。表示装置120は、画像処理装置1000が出力する画像データにもとづき被検体の画像を表示する。表示装置はシステムの一部であってもよく、システム外部に存在してもよい。
画像処理装置1000は、信号取得部1010、光音響画像取得部1020、広域画像取得部1030、計測位置取得部1040、比較ペア生成部1050、投影画像生成部1060を備える。装置はさらに、比較部1070、パルス位置推定部1080、統合ボリューム生成部1090、表示制御部1100を備える。画像処理装置は、CPUやメモリ等の演算資源を備え、プログラムの指示に従って所定の情報処理を行う情報処理装置(P
Cやワークステーション等)により構成される。画像処理装置の各ブロックは、それぞれ独立した回路として構成されても良く、プログラムモジュールとして仮想的に実現されても良い。
光音響画像取得部1020は、1回のパルス光照射で取得した光音響信号に基づく画像再構成処理を行い、対応するFOVの光音響画像をボリュームデータの形式で生成する。なお、パルス光照射毎に生成されるボリュームデータを以下ではパルスボリュームと呼ぶ。光音響画像取得部1020は、生成したパルスボリュームを投影画像生成部1060および統合ボリューム生成部1090に局所画像として出力する。
計測位置取得部1040は、光音響信号計測装置110の信号計測部の位置情報(信号計測位置)を取得し、比較ペア生成部1050に出力する。
投影画像生成部1060は、ボリュームデータから投影画像を生成し、比較部1070に出力する。
表示制御部1100は、統合ボリューム生成部1090が生成した統合ボリュームを表示装置130に表示するための表示制御を行う。
図5は、画像処理装置1000が行う全体の処理手順を示すフローチャートである。
(ステップS5010)信号データ取得
信号取得部1010は、光音響信号計測装置110が被検体を複数回撮像して得た光音
響信号データ(以下、信号データと称する)を取得する。本実施形態における撮像回数は、N_pulse回であり、信号データの個数は、N_pulse個である。1回の(i回目の)パルス
照射によって受信した信号データをP_i(1≦i≦N_pulse)と表記する。信号データP_iは
、探触子がある所定の位置に配置されたときに、複数の受信素子が受信した時系列の音響信号である。また、P_iの添え字i(1≦i≦N_pulse、iは正の整数)を「パルスインデックス」と呼ぶ。各信号データは、多次元のデータ列等の一般的な形態で表現される。
光音響画像取得部1020は、信号データP_i(1≦i≦N_pulse)を用いて再構成処理を行い、パルスボリュームVp_i(1≦i≦N_pulse)を生成し、画像処理装置1000に記録
する。再構成には一般的な方法、例えばトモグラフィ技術で普通に用いられるタイムドメインあるいはフーリエドメインでの逆投影法などを利用できる。この時、各受信素子の位置が考慮される。生成されるパルスボリュームVp_iは3次元ボリュームの画像データであり、Vp_i(x,y,z)のような関数表現によってデータの値(ボクセル値)を表記できる。
本実施形態でのパルスボリュームは、被検体を撮像可能な範囲(FOV)を含む局所領域の3次元ボリュームである。
広域画像取得部1030は、赤外線カメラ301、302、303が撮影した被検体の赤外線カメラ画像(ICAM1,ICAM2,ICAM3)を取得する。ここでは、赤外線カメラ画像は
ある瞬間の静止画像とする。赤外線静止画像取得のタイミングは、被検体の状態(位置・姿勢など)が、光音響計測の時点と変わらない時間内とする。ただし、赤外線カメラが動画像を撮影し、広域画像取得部1030が動画像のフレーム中から好適な静止画像を選択しても良い。その場合、各フレーム間で画素値の差分をとるなどの手法で画像の変化が小さいフレームを検出することが好ましい。
計測位置取得部1040は、光音響信号計測装置110の信号計測位置を取得する。信号計測位置とは、ステップS5010で取得したN_pulse個の信号データそれぞれの計測
時点における信号計測部の位置であり、装置座標系CDEVにおける3次元の位置情報である。ステップS5010で取得した信号データや、それに基づいて生成されたパルスボリュームは、この信号計測位置を基準とした所定範囲(FOV)における計測結果である。本実施形態では、信号計測位置をPosS_i(1≦i≦N_Pulse)と表記する。信号計測部が位置
変化手段によって移動しながら光音響計測を行う場合はパルスごとに信号計測位置が異なる。
比較ペア生成部1050は、信号計測位置PosS_iを取得し、後続の比較処理を行う比較ペアを決定して比較ペア情報を生成する。ペアの決定基準として例えば、互いにオーバーラップするパルスボリュームをペアとする方法がある。
〜5)に応じたパルスボリュームV1からV5の領域を表す。比較ペア生成部1050は、こ
れらのパルスボリュームについて、オーバーラップ領域を有するパルスボリュームのペアを検出する。ここでは、「V1とV2」、「V2とV3」、「V3とV4」、「V3とV5」、「V4とV5」の領域がオーバーラップしている。そこで、比較ペア生成部1050は、これらのペアを示す情報を、パルスインデックスの組として生成する。図6であれば、R_1={1、2}、R_2={2、3}、R_3={3、4}、R_4={3、5}、R_5={4、5}となる。
ステップS5060において、投影画像生成部1060は、パルスボリュームVi(1≦i≦N_pulse)の夫々について、ボリューム値を投影して投影画像を生成する。例えば、最
大輝度値投影(MIP:Maximum Intensity Projection)画像を生成する。これには透視投影法など、一般的な手法を利用できる。投影画像の投影方向は、各赤外線カメラの撮像方向(304、305、306)と略同一とする(好ましくは一致させる)。つまり、装置座標系CDEVのZ軸方向、X軸方向、Y軸方向夫々において投影画像を生成する。これは、後述するステップS5080において、投影画像(局所領域画像)と赤外線カメラ画像(被検体広域画像)とを同一方向から比較するためである。
変換を掛けて赤外カメラ座標系に変換する。そして、その結果にさらに夫々TC1toI1、TC2toI2、TC3toI3の変換を掛けることで、赤外線カメラ画像座標系に変換できる。すなわち
、パルスボリュームの各ボクセルが赤外カメラ画像上のどの位置に透視投影されるかを算出できる。こうして生成された透視投影による投影画像は、2次元座標への投影方法が赤外線カメラ画像と概ね一致する。
、Ipy_i、Ipz_i)および正射影による投影画像(Iox_i、Ioy_i、Ioz_i)は、比較部10
70に出力される。
比較部1070は、比較ペア情報R_j(1≦j≦N_pair)が表す全てのペア(パルスイン
デックスの対)について、比較ペアとされたパルスボリュームを比較し、ペア間の一致度に関する情報を取得する。具体的には、両者の投影画像の類似度関数と、その類似度関数のピーク位置を取得する。そしてこの情報を、投影画像間の一致度に関する情報、すなわち、投影画像の元となったパルスボリュームのペア間の一致度に関する情報とする。これは、被検体の局所領域が写った画像同士の比較処理に当たる。
類似度を値とした画像類似度関数を用いる。具体的には、比較ペア情報R_jが表すパルス
インデックスの対(R_j,1およびR_j,2と表記する)の夫々の投影画像について、以下の3個の類似度関数を取得する。
進させた場合の画像間の類似度を算出する関数である。例えば、SSD(Sum of Squared Difference)やSAD(Sum of Absolute
Difference)、相互情報量、相互相関など、任意の類似度尺度が適用できる。ここで関数f_similは、画像間の類似度が高い場合には関数値として高い値を返すものと
する。
量(x,y,z)を所定の範囲内で離散的に変化させた場合の関数値の取得を意味する。比較
ペア(R_j,1およびR_j,2)の夫々のパルスボリュームの信号計測位置はPosS_R_j,1お
よびPosS_R_j,2であり、互いに異なっている。
」までの整数値だけ前後に変化させた場合にFLX_jが返す、「(2K+1)×(2K+1)個」の値の
集合を取得する。この値の集合は、離散画像のようなメモリアレイ上に展開されるデータである。より発展的には、「(2K+1)×(2K+1)個」の値の集合に対してバイリニア法やバイキュービック法などを適用して、類似度関数FLX_j(y,z)を連続関数(またはそれに近い
関数)として導出してもよい。本実施形態では、比較ペアごとに、各類似度関数(FLX_j,FLY_j,FLZ_j:1≦j≦N_pair)を2次元の連続関数として取得する。
は、投影画像801(正射影)と投影画像802(正射影)とを比較することを表す。このように本ステップでは、比較ペアとなっている投影画像の間の血管走行の類似性を比較することで、類似度関数が算出される。
が最大となる位置PosLX_j,PosLY_j,PosLZ_jを算出する。
の相対位置の推定値(個別最適値)を示している。比較ペア間で体動が生じていない場合には、類似度ピーク位置PosL_jの値は、信号計測位置のペア間の差異ΔPosS_jと一致することが期待される。言い換えると、PosL_jとΔPosS_jの差異(ΔPosS_j‐PosL_j)が、局所画像に基づいて取得した当該ペア間での(相対的な)体動の観測情報となる。
比較部1070は、パルスボリュームの夫々と赤外カメラ画像を比較し、パルスボリュームの夫々と赤外カメラ画像との一致度に関する情報を取得する。つまり、各パルスボリュームの透視投影による投影画像Ipx_i、Ipy_i、Ipz_i(1≦i≦N_pulse)と、投影方向が夫々対応する赤外線カメラ画像ICAM2、ICAM3、ICAM1との類似度関数と、その類似度関数
のピーク位置を得る。そしてこの情報を、投影画像の夫々と赤外線カメラ画像の一致度に関する情報、すなわち、投影画像の元となったパルスボリュームの夫々と赤外カメラ画像との一致度に関する情報とする。ここで、投影画像にも赤外線カメラ画像にも被検体の血管が描出されているため、血管の輝度情報を用いて両画像の類似性を評価できる。ただし本実施形態での血管の輝度は、投影画像中では周囲より明るく、赤外線カメラ画像中では周囲より暗い。そこで、赤外線カメラ画像と投影画像のいずれかの輝度値の明暗を反転させて比較を行うと良い。
きの類似度を計算することにより、被検体の広域画像と局所画像を比較することである。まず、各パルスボリュームに対応する夫々の投影画像について、以下の3個の類似度関数を取得する。
ップでも、各類似度関数(FGX_i,FGY_i,FGZ_i)の取得とは、各関数の引数である並進量(x,y,z)を所定の範囲内で離散的に変化させた場合の関数値の取得を意味する。上式より、赤外線カメラ画像に対して投影画像の位置を並進させた場合の、画像間の類似度を算出する関数が得られる。ここで、ステップS5060の処理において、信号計測位置を考慮して赤外線カメラの撮像方向への透視投影によって生成された投影画像は、既に赤外線カメラ画像と同一座標系で表されている。従って本ステップの類似度は、パルスボリュームの絶対位置ではなく、信号計測位置からの相対位置をパラメータとして算出する。すなわち、位置ズレのない状態を表す(x,y,z)=(0,0,0)を基準として、限定的な範囲内で投影画像を並進させて類似度を計算する。本実施形態では、N_pulse個のパルスボリューム
に対応する夫々の投影画像について、各類似度関数(FGX_i,FGY_i,FGZ_i:1≦i≦N_pulse)を、2次元の連続関数として取得する。
スボリュームの推定位置(個別最適値)を示している。また、PosG_iの値そのものが、広域画像に基づいて取得した当該パルスボリュームにおける(絶対的な)体動の観測情報となる。
パルス位置推定部1080は、S5070の類似度ピーク位置PosL_j(1≦j≦N_pair)と、S5080の類似度ピーク位置PosG_i(1≦i≦N_pulse)から、パルス位置推定量Pos_i(1≦i≦N_pulse)を算出する。すなわち、被検体の局所画像(パルスボリューム)間
の相対的な位置に関する個別最適値と、被検体の広域画像(赤外線カメラ画像)に対する各局所画像の位置に関する個別最適値と、の両方をなるべく保ちつつ、全てのパルスボリュームの位置を全体最適化する。このパルス位置推定量Pos_iは、パルスボリューム間の比較で最も類似する位置関係(矢印803:第1の比較)と、赤外線カメラ画像とパルスボリュームの比較で最も類似する位置関係(矢印804、805:第2の比較)の両方をなるべく満たす値である。
像間の相対的な位置の個別最適値からの、推定値のズレ量をコストとする関数である。
相違等によって算出される。つまり、被検体の広域画像に対する局所画像の位置の個別最適値からの、推定値のズレ量をコストとする関数である。関数EG_iは、具体的には次式によって計算する。
ュートン法などの繰り返し計算による非線形最適化法や、線形最適化法を利用できる。また、コスト関数Eの定義は上記に限らない。例えば、式16や式17はズレ量の評価尺度
としてL2ノルムを用いているが、L1ノルム等の他の距離尺度を用いても良い。また、例えば式15に加え、各パルスボリュームの位置の変動(動き)に対する正則化を掛けても良い。一般的に撮像中の体動は極端に大きくはなく、連続的で滑らかだと想定される。よってパルスボリュームの位置の最適化の際には、上記の想定から逸脱するような動きが算出されることを防ぐために、正則化処理を行う。正則化手法として例えば、位置の補正量(信号計測位置PosS_iとパルス位置推定量Pos_iの差)の総和に所定の重み係数をかけ
て式15に加算する方法がある。また、補正量の微分(加速度)の総和を用いる方法、補正量の周波数成分値に基づいて算出した値を用いる方法などがある。また、被検体の典型的な変動(体動)の仕方をモデルとして用意し、そのモデルと補正量の相違をコストとして式15に加算する方法もある。
統合ボリューム生成部1090は、パルスボリュームVp_i(1≦i≦N_pulse)を統合し
て統合ボリュームVoを生成する。このとき統合ボリューム生成部1090は、各パルスボリュームVp_iの夫々を、ステップS5090で最後に得たパルス位置推定量Pos_iだけ並
進して共通の座標系に変換した後に統合する。統合ボリューム生成部1090は、まず、統合ボリュームの範囲として、並進後の各パルスボリュームの領域を包含する包含領域の範囲を算出する。次に、この包含領域において、同一座標に変換された各パルスボリュームのボクセル値を平均化して統合ボリュームを生成する。統合ボリュームは、表示制御部1100に送信されると共に、不図示の記憶媒体に保存される。
表示制御部1100は、出力画像Voの情報を表示装置130に出力して表示させる。表示方法としては、統合ボリュームを所定の方向に最大値投影処理して得られる2次元画像(投影画像)を表示する方法、統合ボリュームをボリュームレンダリングする方法、統合ボリュームを任意の断面で切断した断層画像を表示する方法などがある。
ステップS5080では、ステップS5060で生成された同じ投影画像Ipx_i,Ipy_i,Ipz_iの位置を並進させながら類似度関数を算出していた。しかし、装置座標系CDEVに
おいてパルスボリュームの位置が変わると、透視投影によって生成される投影画像も変化する。そこで、ステップS5060において、パルスボリュームVi(1≦i≦N_pulse)の
夫々について、装置座標系CDEVにおいて、信号計測位置からの並進量(x,y,z)を所定の
範囲内で離散的に変化させて、夫々の方向の投影画像を生成する。
おいて、その並進量に対応する投影画像を選択した上で類似度を算出する。例えばステップS5080で類似度関数FGX_j(y_k,z_k)を算出する場合、赤外線カメラ画像ICAM2との類似度を計算する投影画像Ipx_iとして、ステップS5060で生成した並進量(0,y_k,z_k)に対応する投影画像を選択する。他の方向についても同様に類似度関数を算出できる。本変形例の方法であれば、パルスボリュームの位置をより反映させた類似度が取得できる。
上記第1の実施形態では、3台の赤外線カメラの撮像方向は、装置座標系の座標軸と一致していた。しかし、赤外線カメラ画像内に被検体の外形全体が含まれるのであれば、撮像方向はこれに限定されない。例えば装置の構造によっては、赤外線カメラ301のZ軸方向(装置座標系CDEV)に被検体を配置できない場合がある。この場合、赤外線カメラ301の向く方向を、Y軸を中心としてZ軸を45°回転させた方向とする。この方向を表すベクトルを、dと定義する。
で、Z方向への投影画像であるIz_iの代わりに、投影画像Id_iを生成する。次に、ステップS5070とS5080の類似度関数の算出では、Z方向の投影画像Iz_iに対応する類似度関数FLZ_j(x,y)及びFGZ_i(x,y)の代わりに、下式のFLD_j(x,y,z)、FGD_i(x,y,z)を取得する。
だけ並進させた場合の画像間の類似度を算出する関数である。これは、ステップS5070で説明したf_simil(I1,I2,x,y)と比べ、並進軸が一つ多い。このときの並進量(x,y,z)は装置座標系CDEV上で定義しておく。ここで、関数FLD_j(x,y,z)、及びFGD_i(x,y,z)を取得するために、並進量(x,y,z)の値を所定の範囲で離散的に変化させて類似度
を計算する。ただし、ベクトルdに平行な方向での並進では投影画像の相対位置は変わら
ない。そこで、並進量(x,y,z)を、ベクトルdに直交する平面上で、並進量が0である位
置を基準にした所定の領域内(例えば矩形領域)で変化させる。これにより、類似度関数FLD_j(x,y,z)及びFGD_i(x,y,z)が取得できる。他の赤外線カメラ302、303についても同様に処理することで、装置座標系と赤外線カメラの撮像方向が一致しない場合でも、赤外線カメラ画像とパルスボリュームを比較できる。
上記第1の実施形態では、ステップS5060の処理として、パルスボリュームの投影画像としてMIP画像を生成した。しかし、例えばMIP画像に代えて、MinP画像(最小値投影画像)を生成しても良い。また、MIP画像とMinP画像を両方生成し、以降の処理では両方または好ましい片方を用いても良い。
上記第1の実施形態では、ステップS5070、S5080で、パルスボリューム間の並進量と、赤外線カメラ画像に対する投影画像の並進量と、に関する類似度関数を取得した。しかし、位置関係として、並進に加えて回転を考慮しても良い。この場合、ステップS5070及びS5080で取得する類似度関数は、各投影画像間の、投影面上での並進量および投影面内の回転量に関する類似度関数となる。この場合、ピーク位置は、3次元ベクトルとしてではなく、例えば、投影画像間の類似度関数を最大とする剛体変換行列、及び赤外線カメラ画像と投影画像間の類似度関数を最大とする剛体変換行列として算出する。
上記第1の実施形態では、ステップS5060の処理において、直交三方位に投影した投影画像を生成していた。しかし例えば、Z軸方向とY軸方向のみ投影画像を生成し、X軸方向への投影を省略しても良い。この場合でも、Z軸方向の投影画像に基づいてXY平面上の並進の情報が取得でき、またY軸方向の投影画像に基づいてXZ平面上の並進の情報が取得できる。そのため、全ての軸方向の動きの情報を取得できる。この方法により第1の実施形態に比べて計算コストを低減できる。この場合、赤外線カメラは、Z軸方向とY軸方向の2台のみでよい。
本発明の第2の実施形態について説明する。第1の実施形態では、計測された夫々の光音響画像の間の比較と、夫々の光音響画像と赤外線カメラ画像との比較と、を共に行って(これらの関係を同時に考慮するコスト関数を用いて)体動を推定し補正した。第2の実施形態では、最初に、計測された夫々の光音響信号の間の比較により第1の推定を行い、次に、第1の推定による補正を行った結果を赤外線カメラ画像と比較することで第2の推定を行う。したがって、最終結果に直結する第2の推定では、赤外線カメラ画像のみを目標とした推定が行われる。そのため、より赤外線カメラ画像に合致した光音響画像を取得できる。
以下、図9を用いて本実施形態に係る光音響装置の構成について説明する。ただし、第1の実施形態と同様の構成については、同一の番号を付し、説明は省略する。
画像処理装置9000は、信号取得部1010、光音響画像取得部1020、広域画像取得部1030、計測位置取得部1040、比較ペア生成部1050、投影画像生成部9060を備える。さらに、第1の比較部9070、パルス位置推定部9080、統合ボリューム生成部9090、変形仮説生成部9110、変形ボリューム生成部9120、第2の比較部9130、表示制御部1100を備える。
パルス位置推定部9080は、第1の比較部9070が算出したペア間の一致度に関する情報に基づいて、夫々のパルスボリュームの位置の推定値(パルス位置推定量)を算出して、統合ボリューム生成部9090に出力する。
変形仮説生成部9110は、後述する第2の比較部9130で算出された推定変形パラメータを基準として、統合ボリューム生成部9090から取得した統合ボリュームを変形させるための複数の変形パラメータの仮説(以下、変形仮説と呼ぶ)を生成する。そして、生成した複数の変形仮説を変形ボリューム生成部9120と第2の比較部9130に出力する。
図10は、画像処理装置9000が行う全体の処理手順を示すフローチャートである。(ステップS10010〜S10070)
これらのステップにおいては、第1の実施形態のステップS5010からステップS5070と同様の処理が実行される。
パルス位置推定部9080は、ステップS10070で算出した類似度ピーク位置PosL_j(1≦j≦N_pair)に基づき、パルス位置推定量Pos_i(1≦i≦N_pulse)を算出する。すなわち、被検体の局所画像(パルスボリューム)間の相対的な位置に関する個別最適値をなるべく保ちつつ、全てのパルスボリュームの位置を全体最適化する。具体的には、コスト関数Eを次式のように定義し、このコスト関数Eを小さくするように各パルスボリュームの位置を最適化する。
5090におけるコスト関数E(式15)の第2項をなくしたものに相当する。コスト関
数Eの最適化の方法は、第1の実施形態のステップS5090と同様である。
統合ボリューム生成部1090は、パルスボリュームVp_i(1≦i≦N_pulse)を統合し
て統合ボリュームVoを生成する。この生成方法は、第1の実施形態におけるステップS5100と同様である。統合ボリュームは変形仮説生成部9110及び変形ボリューム生成部9120に出力されると共に、不図示の記憶媒体に保存される。
変形仮説生成部9110は、推定変形パラメータの現在の値を基準として、統合ボリューム生成部9090から取得した統合ボリュームVoを変形させるための複数の変形パラメータの仮説(変形仮説)を生成する。ここで、推定変形パラメータの現在の値とは、本ステップを最初に処理する場合には、初期値、すなわち「変形なし」を表す値を取る。一方、本ステップを2回目以降に処理する場合には、後述するステップS10140で算出される推定変形パラメータp_eが、推定変形パラメータの現在の値となる。
。このとき変形パラメータは、FFDの各制御点の制御量で表される。具体的には、制御点
の(x,y,z)の各軸方向の制御量の集合を表す要素数N_grid(N_gridは制御点の数)のベク
トルを夫々、p_x、p_y、p_zとする。そして、それら全てを要素に持つベクトルp={p_x,p_y,p_z}Tを変形パラメータと定義する。
たとき、p_eの各要素に対して、夫々微小な変化量を与えることを考える。但し、本ステ
ップを最初に処理する場合は変形仮説生成部9110から変形推定パラメータp_eが与え
られないため、全ての要素が0のベクトルp_0をp_eとする。このとき、各要素に与えられた微小な変化量からなるベクトルをΔpと定義する。
させる。例えば制御点の数がN_gridのとき、パラメータ空間は、N_grid×3次元の空間になる。このとき、ある変化量Δpの組を、パラメータ空間中のある軸の要素のみを所定の
値dだけ移動させたベクトルの組として生成する。この処理を全てのパラメータ空間中の
軸(N_grid×3)に対して行うことにより、(N_grid×3)パターンの変化量Δp_i(1≦i≦N_grid×3)を取得できる。つまり、ベクトルの要素数分のΔp_iの組を取得する。なお、変化量Δpの組み合わせは、パラメータ空間中で均等に変化するように変動していれば
、この例に限定されない。
成する(本ステップを最初に処理する場合にはp_i=Δp_iとなる)。ここで、生成された
変形パラメータの数をN_hypoと定義する(本実施形態ではN_hypo=N_grid×3)。このと
き、生成された変形パラメータの集合{p_1,…,p_i,…,p_N_hypo}と、各制御点の位置情報を合わせたものを、変形仮説Hと定義する。そして、変形仮説部9110は、生成し
た変形仮説Hを変形ボリューム生成部9120に出力するとともに、不図示の記憶媒体に
保存する。なお、変形を表現するモデルはFFDに限られず、Thin Plate Splineなど、非線形な変形を表現する他の方法も利用できる。
変形ボリューム生成部9120は、変形仮説生成部9110から取得した変形仮説Hに
基づいて、不図示の記憶媒体に保存された統合ボリュームVoを変形させて変形ボリュームDV_iを生成する(1≦i≦N_hypo)。具体的には、FFDの変形アルゴリズムを用いて、変形仮
説Hに含まれる制御点の位置情報と夫々の変形パラメータΔp_iとに基づいて統合ボリュームVoを変形させる。変形ボリュームDV_iは投影画像生成部9060に出力される。
投影画像生成部9060は、変形ボリューム生成部9120が生成した変形ボリュームDV_iから、赤外線カメラの撮像方向ごとに変形投影画像を生成する。これは、ステップS10060におけるパルスボリュームを、変形ボリュームDV_iに入れ替えただけの処理である。本ステップにおける投影方法は、第1の実施形態の「変形例1−1」のような透視投影法とする。装置座標系CDEVのX、Y、Z方向の夫々での投影により生成された変形投影画像を夫々、DIpx_i、DIpy_i、DIz_i(1≦i≦N_hypo)と表記する。
第2の比較部9130は、変形投影画像DIpx_i、DIpy_i、DIz_i(1≦i≦N_hypo)と、
夫々が対応する投影方向の赤外線カメラ画像ICAM2、ICAM3、ICAM1との類似度を算出する
。算出された類似度を夫々Sx_i、Sy_i、Sz_iと表記する。ここでは第1の実施形態のステップS5070と同様の類似度尺度を適用できる。
第2の比較部9130は、ステップS10130で算出された類似度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の3つの類似度の平均値として定義できる。コスト関数の定義方法はこれに限らず、例えば3つの類似度の中の最大値を採用してもよい。なお、コスト関数E_sは、変形パラメ
ータpを与えたときの類似度を表すため、pの関数としてE_s(p)と表記できる。従って、変形パラメータp=p_iに対応する(類似度Sx_i、Sy_i、Sz_iに基づく)コスト関数の値はE_s(p_i)と表記される。
ラメータp_enewを推定する。具体的には、パラメータp_enewを、以下の式によって算出する。
動させた変形パラメータp_enewを取得できる。
)を算出して、不図示の記憶媒体に保存する。このとき、本ステップが実行されて新しいコスト値E_s(p_enew)が算出される度に、過去のコスト値を保持したまま、該コスト値
が不図示の記憶媒体に追加格納されるものとする。なお、上述した式21の演算で使用するコスト値E_s(p_e)は、1回前に本ステップの処理を実行した際のコスト値E_s(p_enew)であるため、不図示の記憶媒体から取得できる。ただし、本ステップの処理を最初に
行う場合には、コスト値E_s(p_e)が保持されていないので、式21の演算に先立ち初期値p_0に関するコスト値E_s(p_0)を算出する必要がある。
を見つける処理を行っている。これは、最急降下法を用いてコスト関数E_sを最適化する
際の1ステップに相当する。従って、本ステップが繰り返されることで、推定変形パラメータp_enewがより好ましい値になる。なお、コスト関数E_sの最適化方法はこれに限られ
ず、例えばニュートン法などでも良い。
第2の比較部9130は、不図示の記憶媒体に保存された歴代のコスト値E_s(p_e)を参照し、値が収束しているか否かを判定する。収束していれば、変形完了命令の情報とともに推定変形パラメータp_eを最終変形パラメータp_finalとして変形ボリューム生成部9120に出力し、処理をステップS10160に移行する。収束していなければ、推定変形パラメータp_eを変形仮説生成部9110に出力し、ステップS10100に戻る。
変形ボリューム生成部9120は、第2の比較部9130から変形完了命令の情報と、最終変形パラメータp_finalを取得し、出力ボリュームV_finalを生成する。
以上のように本実施形態では、ステップS10100からS10150を、赤外線カメラ画像に合致させることを目標とするコスト関数のコスト値が収束するまで繰り返した後に、ステップS10160で最終変形ボリュームを生成する。これにより、コスト関数を最適にする変形パラメータを取得できる。
表示制御部1100は、出力ボリュームV_finalの情報を表示装置130に出力して表
示する制御を行う。本処理は、第1の実施形態のステップS5110と同様である。
本発明の第3の実施形態について説明する。第1および第2の実施形態では、赤外線カメラ画像を広域画像として利用した。第3の実施形態では、広域画像として被検体の距離画像を利用する。これにより、ボリュームデータである光音響画像の情報を縮退させた2次元空間上ではなく、3次元空間上での情報の比較が可能になる。
以下、図11を用いて本実施形態に係る光音響装置の構成について説明する。第1の実施形態と同様の構成については、同一の番号を付し、説明は省略する。
距離画像カメラ140は、レーザを被検体に投射し、投射したレーザが被検体から反射して往復してくる時間に基づいて被検体までの距離を計測する(Time Of Flight)カメラである。これにより被検体の3次元表面形状が距離画像の形式で得られる。距離画像カメラ140は、被検体全体を計測できる位置に設置される。得られた距離画像(IRNGと表記)は画像処理装置11000に入力される。
る。
局所表面形状抽出部11060は、被検体の局所領域であるFOVを描出したパルスボリュームから被検体の表面形状を抽出し(以下、局所表面形状と呼ぶ)、比較部11070に出力する。
致度に関する情報として、ペアの間の類似度等を算出する。また、夫々のパルスボリュームと距離画像とを比較する処理(第2の比較)として、表面形状抽出部11060から取得した局所表面形状と広域表面形状取得部11030から取得した広域表面形状を比較し、一致度に関する情報を算出する。さらに、算出した一致度に関する情報をパルス位置推定部11080に出力する。
図12は、画像処理装置11000が行う全体の処理手順を示すフローチャートである。
(ステップS12010からS12020)
これらのステップにおいては、第1の実施形態におけるステップS5010からステップS5020と同様の処理が行われる。
広域表面形状取得11030は、距離画像カメラ140が撮像した被検体の距離画像(IRNG)を取得する。ここでは距離画像として、被検体のある瞬間を写した静止画像を取得する。被検体の位置・姿勢が、光音響信号計測装置110による計測時とほぼ同じなら、静止画像の取得時点は問わない。なお、距離画像として動画像を取得し、画像の変化が小さいフレームを抽出しても良い。次に、広域表面形状取得11030は、任意の既知の手法によって、距離画像IRNGから被検体全体の表面形状を計測した広域表面形状(以下、Surf_Gと呼ぶ)を生成する。
これらのステップにおいては、第1の実施形態におけるステップS5040からステップS5050と同様の処理が行われる。
表面形状取得部11060は、パルスボリュームVi(1≦i≦N_pulse)の夫々について
、被検体の局所表面形状を抽出する。
比較部11070は、比較ペア情報R_j(1≦j≦N_pair)が表す全てのペアに対して、
比較ペアとされたパルスボリュームを比較し、ペア間の一致度に関する情報を取得する。具体的には、両パルスボリュームの類似度関数と、類似度関数のピーク位置を取得する。すなわち、下式のような、比較するパルスボリューム間の相対位置を(x,y,z)だけずら
しながら比較した場合のボリューム間の類似度関数FL_j(x,y,z)を算出する。
させた場合の画像間の類似度を算出する関数である。式22は、第1の実施形態のステップS5070における式1に対して、並進量の次元を1つ増やしたものであるので、詳細の説明は省略する。本実施形態では、比較ペア情報の夫々について、各類似度関数 (FL_j:1≦j≦N_pair)を3次元の連続関数として取得する。
比較部11070は、パルスボリュームの夫々と距離画像を比較し、パルスボリュームの夫々と距離画像との一致度に関する情報を取得する。具体的には、S12060で得た局所表面形状SurfL_i(1≦i≦N_pulse)の夫々について、S12030で得た広域表面形状SurfGとの間の一致度関数と、一致度関数のピーク位置を取得する。本ステップの目的
は、広域表面形状と各パルスボリュームの局所表面形状の類似度を計算することで、被検体の広域画像と局所画像を比較することである。各パルスボリュームに対応する夫々の局所表面形状SurfL_iについて、以下の一致度関数を取得する。
相対位置を(x,y,z)だけずらしながら比較した場合の形状間の一致度を算出する関数で
ある。本実施形態では表面形状は点群で表現されるため、点群間で互いに最近傍となる点同士を対応付け、対応付けられた点群間の平均距離を適用すれば良い。但し、一致度の算出方法はこれに限られない。例えば、表面形状を表す点群から三角パッチなどによるメッシュをはり、メッシュを構成する面同士の距離を比較して一致度を算出する方法もある。本実施形態では、パルスボリュームの夫々について、一致度関数(FG_i:1≦i≦N_pulse)
を3次元の連続関数として取得する。
大となる位置PosG_iを算出する。
比較部11070は、S12070で得た類似度ピーク位置PosL_j(1≦j≦N_pair)と、S12080で得た一致度ピーク位置PosG_i(1≦i≦N_pulse)からパルス位置推定量Pos_i(1≦i≦N_pulse)を求める。すなわち、被検体の局所画像(パルスボリューム)間
の相対的な位置に関する個別最適値と、被検体の広域画像(距離画像)に対する各局所画像の位置に関する個別最適値と、の両方をなるべく保ちつつ、全てのパルスボリュームの位置を全体最適化する。これは第1の実施形態のステップS5090と同様の処理である。その際、S5090における投影画像間の類似度ピーク位置をパルスボリューム間の類似度ピーク位置に、赤外線カメラ画像に対する投影画像の類似度ピーク位置を広域表面形状に対する局所表面形状の一致度ピーク位置に、夫々置き換える。
これらのステップにおいては、第1の実施形態のステップS5100からステップS5110と同様の処理が行われる。
上記の第3の実施形態では、ステップS11030の処理として、距離画像から被検体の広域表面形状の抽出を行った。しかし、広域表面形状の取得方法はこれに限られない。例えば、被検体を超音波測定して得られた超音波画像から抽出してもよい。その場合、広域表面形状取得部11030は、不図示の超音波撮像装置を用いて被検体全体の超音波ボリュームを取得し、広域表面形状を抽出する。この場合例えば、被検体の内部と外部での音響インピーダンスの違いに由来する境界を閾値処理によって抽出することで、被検体の広域表面形状を取得できる。これにより、例えば被検体が水に浸かっている場合は、距離画像カメラでは撮像条件が悪いのに対し、水と被検体(人体など)の間では音響インピーダンスの違いが大きいため、精度良く広域表面形状を取得できる。
上記の第3の実施形態では、ステップS12070の処理として、パルスボリューム間の画像類似度を一致度として、そのピーク位置を算出していた。しかし、パルスボリューム間の一致度を算出する方法は他の何れの方法であってもよい。例えば、画像処理によって夫々のパルスボリュームから血管分岐等の解剖学的な特徴を抽出し、その位置や分布の一致度をパルスボリューム間の一致度として用いてもよい。あるいは、パルスボリューム間の局所表面形状の一致度を算出して、これをパルスボリューム間の一致度として用いてもよい。これは、ステップS12080の広域表面形状と局所表面形状の一致度関数を算出する処理において、広域表面形状ではなく、異なる2つのパルスボリュームの局所表面形状を採用することで算出できるため、説明を省略する。これにより、3次元空間中に敷き詰められたボクセル単位での比較を要する画像類似度の計算よりも、3次元空間の限られた領域の比較しか必要としないため、計算コストを削減できる。
第4の実施形態における画像処理装置は、被検体の広域画像として、赤外線カメラ画像でなく、被検体の外観のみが描出されたカメラ画像を利用する。以下、上記各実施形態と
異なる点を中心に説明する。
以下、図13を用いて本実施形態に係る光音響装置の構成について説明する。ただし、上記実施形態と同様の構成については、同一の番号を付し、説明は省略する。
カメラ150は、被検体の外観を(可視光領域の光を)撮影するカメラであり、被検体全体の外観を撮影可能な位置に設置されている。カメラ150は第1の実施形態における赤外線カメラ120と違ってカラー画像を撮影する。それ以外の点では第1の実施形態における赤外線カメラ120の説明における「赤外線カメラ」の部分を「カメラ」に置き換えたものと考えて良い。
広域外形形状抽出部13035は、カメラ画像から被検体の広域の2次元外形形状を抽出し、抽出した外形形状(以下、広域外形形状と呼ぶ)を比較部13070に出力する。
比較部13070は、比較ペア情報に基づき、ペアとなったパルスボリュームを比較する処理(第1の比較)として、第1の実施形態の比較部1070と同様の処理を行う。また、夫々のパルスボリュームとカメラ画像とを比較する処理(第2の比較)を行う。具体的には、局所外形形状抽出部13065から取得した局所外形形状と広域外形形状取得部13035から取得した広域外形形状を比較し、局所外形形状と広域外形形状の間の一致度に関する情報を算出し、パルス位置推定部13080に出力する。
図14は、画像処理装置13000が行う全体の処理手順を示すフローチャートである。
(ステップS14010からS14020)
これらのステップにおいては、第1の実施形態におけるステップS5010からステップS5020と同様の処理が実行される。
カメラ画像取得部13030は、カメラ301、302、303が撮影したカメラ画像(ICAM1,ICAM2,及びICAM3)を取得する。この処理は、第1の実施形態におけるステッ
プS5030の「赤外線カメラ」を「カメラ」に置き換えたものと同様である。
広域外形形状抽出部13035は、カメラ画像(ICAM1,ICAM2,及びICAM3)の夫々か
ら被検体の外形形状を抽出する。外形形状の抽出には例えば、一般的なエッジ検出手法を使用できる。抽出した広域外形形状を夫々、Surf_G1、Surf_G2、Surf_G3と定義する。
これらのステップにおいては、第1の実施形態におけるステップS5040からステップS5050と同様の処理が行われる。
局所外形形状抽出部13065は、投影画像の夫々から被検体の局所外形形状を抽出する。投影画像はパルスボリュームの投影画像であるため、写っている構造情報はパルスボリュームと変わらない。従って本ステップは、第3の実施形態のステップS12060の処理を2次元領域において行う場合と同様に実行できる。抽出した局所外形形状を夫々、SurfLx_i、SurfLy_i、SurfLz_i(1≦i≦N_pulse)と表記する。
(ステップS14070)
本ステップでは、第1の実施形態におけるステップS5070と同様の処理が行われる。
比較部13070は、パルスボリュームの夫々とカメラ画像を比較し、パルスボリュームの夫々とカメラ画像との一致度に関する情報を取得する。具体的には、局所外形形状SurfLx_i、SurfLy_i、SurfLz_iの夫々について、ステップS14035で抽出された対応する広域外形形状Surf_G2、Surf_G3、Surf_G1との間で、形状の一致度関数を取得する。そ
して、取得した一致度関数のピーク位置を取得する。本ステップの目的は、カメラ画像から生成した広域外形形状に対する各パルスボリュームの局所外形形状の類似度を計算することで、被検体の広域画像と局所画像を比較することである。
3次元座標上で局所表面形状を並進させるのに対し、本実施形態では2次元座標上で局所外形形状を並進させる点で異なる。具体的な2次元座標は、FGx_iは(y,z)、FGy_iは(z,x)、FGz_iは(x,y)である。また、一致度関数のピーク位置PosG_i(1≦i≦N_pulse)
を算出する際は、第1の実施形態のステップS5080における各軸方向の類似度関数FGX_i,FGY_i,FGZ_iを一致度関数に置き換えれば良い。
比較部13070は、S14070で得た類似度ピーク位置PosL_j(1≦j≦N_pair)とS14080で得た一致度ピーク位置PosG_i(1≦i≦N_pulse)から、パルス位置推定量Pos_i(1≦i≦N_pulse)を算出する。すなわち、被検体の局所画像(投影画像)間の相対
的な位置に関する個別最適値と、被検体の広域画像(カメラ画像)に対する各局所画像の位置に関する個別最適値と、の両方をなるべく保ちつつ、全てのパルスボリュームの位置の全体最適化を行う。これによりパルス位置推定量Pos_i(1≦i≦N_pulse)を算出する。本ステップの処理は、第1の実施形態のステップS5090における赤外線カメラ画像に対する投影画像の類似度ピーク位置を、広域外形形状に対する局所外形形状の一致度ピーク位置に置き換えたものである。
Claims (14)
- 光源と、
前記光源から光を照射された被検体から発生する音響波を受信する受信素子と、
前記音響波を用いて、前記被検体内部の特性情報を示す画像データを生成する処理手段と、
前記被検体における前記光の照射位置を変化させる変化手段と、
前記被検体の広域画像を取得する広域画像取得手段と、
を有し、
前記処理手段は、
前記変化手段により変化する前記照射位置ごとに、当該照射位置に対応する前記被検体の局所領域における前記画像データである局所画像を生成し、
前記照射位置ごとに得られた複数の前記局所画像の間の比較である第1の比較と、前記複数の局所画像のそれぞれと前記広域画像との比較である第2の比較と、に基づいて前記複数の局所画像を統合する
ことを特徴とする光音響装置。 - 前記処理手段は、前記局所領域ごとに、前記音響波を用いて、前記局所画像としてボリュームデータを生成する
ことを特徴とする請求項1に記載の光音響装置。 - 前記処理手段は、前記局所画像の投影画像を生成し、前記投影画像の間の一致度に関する情報を取得することにより、前記第1の比較を行う
ことを特徴とする請求項2に記載の光音響装置。 - 前記処理手段は、互いにオーバーラップする部分を有する前記局所画像をペアとし、前記第1の比較として、当該ペアとなった局所画像を比較する
ことを特徴とする請求項1に記載の光音響装置。 - 前記処理手段は、前記広域画像と前記局所画像の間の一致度が最大になる前記局所画像の位置を取得することにより、前記第2の比較を行う
ことを特徴とする請求項1ないし4のいずれか1項に記載の光音響装置。 - 前記処理手段は、前記第1の比較による前記複数の局所画像の間の一致度に関する情報と、前記第2の比較による前記複数の局所画像のそれぞれと前記広域画像との一致度に関する情報と、を用いて前記統合を行う
ことを特徴とする請求項1ないし5のいずれか1項に記載の光音響装置。 - 前記処理手段は、前記第1の比較により前記複数の局所画像の統合を行ったのち、前記第2の比較により前記統合した画像を変形させる
ことを特徴とする請求項1ないし5のいずれか1項に記載の光音響装置。 - 前記広域画像取得手段は、前記被検体の赤外線画像を取得する赤外線カメラである
ことを特徴とする請求項1ないし7のいずれか1項に記載の光音響装置。 - 前記広域画像取得手段は、前記広域画像として、前記被検体の距離画像に基づく3次元表面形状を取得する
ことを特徴とする請求項1ないし7のいずれか1項に記載の光音響装置。 - 前記広域画像取得手段は、前記被検体の外観が描出された画像を取得するカメラである
ことを特徴とする請求項1ないし7のいずれか1項に記載の光音響装置。 - 前記処理手段は、前記複数の局所画像のそれぞれの投影画像を生成し、前記投影画像と前記広域画像の間の一致度に関する情報を取得することにより、前記第2の比較を行う
ことを特徴とする請求項2または請求項3に記載の光音響装置。 - 前記広域画像取得手段はカメラであり、
前記処理手段は、前記カメラの撮像方向と略同一の方向に前記局所画像を投影することで前記投影画像を生成する
ことを特徴とする請求項11に記載の光音響装置。 - 前記処理手段は、前記複数の局所画像を統合する際に、前記第1の比較および前記第2の比較の結果に基づいて、前記複数の局所画像を生成するための音響波が取得される間に前記被検体の体動があったことによる前記複数の局所画像それぞれの位置ズレを補正することを特徴とする請求項1ないし5のいずれか1項に記載の光音響装置。
- 照射位置を変化させながら光を照射された被検体から発生する音響波を用いて、前記被検体内部の特性情報を示す画像データを生成する画像処理方法であって、
前記被検体の広域画像を取得するステップと、
変化する前記照射位置ごとに、当該照射位置に対応する前記被検体の局所領域における前記画像データである局所画像を生成するステップと、
前記照射位置ごとに得られた複数の前記局所画像の間の比較と、前記複数の局所画像のそれぞれと前記広域画像との比較と、に基づいて前記複数の局所画像を統合するステップと、
を有することを特徴とする画像処理方法。
Priority Applications (2)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2015080675A JP6436442B2 (ja) | 2015-04-10 | 2015-04-10 | 光音響装置および画像処理方法 |
US15/080,779 US10682060B2 (en) | 2015-04-10 | 2016-03-25 | Photoacoustic apparatus and image processing method |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2015080675A JP6436442B2 (ja) | 2015-04-10 | 2015-04-10 | 光音響装置および画像処理方法 |
Related Child Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2018209146A Division JP6598963B2 (ja) | 2018-11-06 | 2018-11-06 | 画像処理装置、画像処理方法およびプログラム |
Publications (2)
Publication Number | Publication Date |
---|---|
JP2016198297A JP2016198297A (ja) | 2016-12-01 |
JP6436442B2 true JP6436442B2 (ja) | 2018-12-12 |
Family
ID=57112209
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2015080675A Active JP6436442B2 (ja) | 2015-04-10 | 2015-04-10 | 光音響装置および画像処理方法 |
Country Status (2)
Country | Link |
---|---|
US (1) | US10682060B2 (ja) |
JP (1) | JP6436442B2 (ja) |
Families Citing this family (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP6714019B2 (ja) * | 2015-05-07 | 2020-06-24 | コーニンクレッカ フィリップス エヌ ヴェKoninklijke Philips N.V. | 医療手順における動き補償のためのシステム及び方法 |
GB2542114B (en) * | 2015-09-03 | 2018-06-27 | Heartfelt Tech Limited | Method and apparatus for determining volumetric data of a predetermined anatomical feature |
JP2018082763A (ja) * | 2016-11-21 | 2018-05-31 | キヤノン株式会社 | 画像生成装置およびその制御方法 |
JP6929048B2 (ja) * | 2016-11-30 | 2021-09-01 | キヤノン株式会社 | 表示制御装置、表示方法、及びプログラム |
JP7320352B2 (ja) * | 2016-12-28 | 2023-08-03 | パナソニック インテレクチュアル プロパティ コーポレーション オブ アメリカ | 三次元モデル送信方法、三次元モデル受信方法、三次元モデル送信装置及び三次元モデル受信装置 |
JP2018126389A (ja) * | 2017-02-09 | 2018-08-16 | キヤノン株式会社 | 情報処理装置、情報処理方法、およびプログラム |
US10062584B1 (en) * | 2017-06-05 | 2018-08-28 | United Microelectronics Corp. | Method for forming semiconductor structure |
US20200163612A1 (en) * | 2017-07-21 | 2020-05-28 | Helmholtz Zentrum Munchen Deutsches Forschungzentrum Fur Gesundheit Und Umwelt (Gmbh) | System for optoacoustic imaging, in particular for raster-scan optoacoustic mesoscopy, and method for optoacoustic imaging data processing |
Family Cites Families (13)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
RU2480147C2 (ru) | 2006-12-19 | 2013-04-27 | Конинклейке Филипс Электроникс Н.В. | Комбинированная система фотоакустического и ультразвукового формирования изображений |
JP5335280B2 (ja) | 2008-05-13 | 2013-11-06 | キヤノン株式会社 | 位置合わせ処理装置、位置合わせ方法、プログラム、及び記憶媒体 |
JP5737858B2 (ja) | 2010-04-21 | 2015-06-17 | キヤノン株式会社 | 画像処理装置、画像処理方法、及びプログラム |
EP2386998B1 (en) * | 2010-05-14 | 2018-07-11 | Honda Research Institute Europe GmbH | A Two-Stage Correlation Method for Correspondence Search |
JP5858636B2 (ja) | 2011-04-13 | 2016-02-10 | キヤノン株式会社 | 画像処理装置、その処理方法及びプログラム |
US20130030288A1 (en) * | 2011-07-28 | 2013-01-31 | Electronics And Telecommunications Research Institute | Image diagnosis apparatus including x-ray image tomosynthesis device and photoacoustic image device and image diagnosis method using the same |
JP5917037B2 (ja) * | 2011-07-29 | 2016-05-11 | キヤノン株式会社 | 被検体情報取得装置および被検体情報取得方法 |
JP5995449B2 (ja) | 2012-01-24 | 2016-09-21 | キヤノン株式会社 | 情報処理装置及びその制御方法 |
JP6146956B2 (ja) * | 2012-03-13 | 2017-06-14 | キヤノン株式会社 | 装置、表示制御方法、及びプログラム |
JP6004714B2 (ja) * | 2012-04-12 | 2016-10-12 | キヤノン株式会社 | 被検体情報取得装置およびその制御方法 |
JP6000705B2 (ja) | 2012-07-17 | 2016-10-05 | キヤノン株式会社 | データ処理装置及びデータ処理方法 |
JP6238549B2 (ja) * | 2013-04-16 | 2017-11-29 | キヤノン株式会社 | 被検体情報取得装置、被検体情報取得装置の制御方法 |
JP6362301B2 (ja) * | 2013-04-30 | 2018-07-25 | キヤノン株式会社 | 被検体情報取得装置、および、被検体情報取得装置の作動方法 |
-
2015
- 2015-04-10 JP JP2015080675A patent/JP6436442B2/ja active Active
-
2016
- 2016-03-25 US US15/080,779 patent/US10682060B2/en active Active
Also Published As
Publication number | Publication date |
---|---|
US10682060B2 (en) | 2020-06-16 |
JP2016198297A (ja) | 2016-12-01 |
US20160296120A1 (en) | 2016-10-13 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP6436442B2 (ja) | 光音響装置および画像処理方法 | |
US10278584B2 (en) | Method and system for three-dimensional imaging | |
US10004462B2 (en) | Systems, methods, and devices for removing prospective motion correction from medical imaging scans | |
Kovalski et al. | Three-dimensional automatic quantitative analysis of intravascular ultrasound images | |
US9691150B2 (en) | Image processing apparatus, image processing method and storage medium | |
US7869663B2 (en) | Methods, systems and computer program products for analyzing three dimensional data sets obtained from a sample | |
CN103502770B (zh) | 使用光学相干断层扫描进行实时跟踪的改进成像 | |
JP5334415B2 (ja) | 試料の機械的歪み及び弾性的性質を測定するプロセス、システム及びソフトウェア | |
US9396576B2 (en) | Method and apparatus for estimating the three-dimensional shape of an object | |
US11432875B2 (en) | Left atrial appendage closure guidance in medical imaging | |
RU2582055C2 (ru) | Совмещение данных изображения для динамической перфузионной компьютерной томографии | |
US9311709B2 (en) | Image processing apparatus and image processing method | |
US20080095417A1 (en) | Method for registering images of a sequence of images, particularly ultrasound diagnostic images | |
WO2011015952A1 (en) | Method and system for stabilizing a series of intravascular ultrasound images and extracting vessel lumen from the images | |
CN116456897A (zh) | 运动补偿激光散斑对比成像 | |
JP7023196B2 (ja) | 検査支援装置、方法およびプログラム | |
JP7071240B2 (ja) | 検査支援装置、方法およびプログラム | |
US10229494B2 (en) | Automated analysis of intravascular OCT image volumes | |
EP2498222A2 (en) | Method and system for regression-based 4D mitral valve segmentation from 2D+T magnetic resonance imaging slices | |
JP7267337B2 (ja) | 眼の画像データ処理 | |
Rosales et al. | Modelling of image-catheter motion for 3-D IVUS | |
US11937977B2 (en) | Compounding and non-rigid image registration for ultrasound speckle reduction | |
JP6598963B2 (ja) | 画像処理装置、画像処理方法およびプログラム | |
WO2018036893A1 (en) | Image processing apparatus and method for segmenting a region of interest | |
JP2023542754A (ja) | 変形を追跡するためのシステムおよび方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20171219 |
|
A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20180926 |
|
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: 20181009 |
|
RD02 | Notification of acceptance of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7422 Effective date: 20181102 |
|
A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20181106 |
|
R151 | Written notification of patent or utility model registration |
Ref document number: 6436442 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R151 |