JP5451414B2 - 被検体情報処理装置および被検体情報処理方法 - Google Patents

被検体情報処理装置および被検体情報処理方法 Download PDF

Info

Publication number
JP5451414B2
JP5451414B2 JP2010008366A JP2010008366A JP5451414B2 JP 5451414 B2 JP5451414 B2 JP 5451414B2 JP 2010008366 A JP2010008366 A JP 2010008366A JP 2010008366 A JP2010008366 A JP 2010008366A JP 5451414 B2 JP5451414 B2 JP 5451414B2
Authority
JP
Japan
Prior art keywords
acoustic wave
receiving elements
subject
signal
wave source
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.)
Expired - Fee Related
Application number
JP2010008366A
Other languages
English (en)
Other versions
JP2011143175A (ja
JP2011143175A5 (ja
Inventor
晃一 丹治
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 JP2010008366A priority Critical patent/JP5451414B2/ja
Priority to US12/975,423 priority patent/US9289126B2/en
Priority to EP10197290A priority patent/EP2345364A1/en
Publication of JP2011143175A publication Critical patent/JP2011143175A/ja
Publication of JP2011143175A5 publication Critical patent/JP2011143175A5/ja
Application granted granted Critical
Publication of JP5451414B2 publication Critical patent/JP5451414B2/ja
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00Measuring for diagnostic purposes; Identification of persons
    • A61B5/0059Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00Measuring for diagnostic purposes; Identification of persons
    • A61B5/0093Detecting, measuring or recording by applying one single type of energy and measuring its conversion into another type of energy
    • A61B5/0095Detecting, 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
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/13Tomography
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • G06T11/006Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Veterinary Medicine (AREA)
  • Pathology (AREA)
  • Public Health (AREA)
  • Biophysics (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Radiology & Medical Imaging (AREA)
  • Algebra (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)

Description

本発明は、被検体情報処理装置および被検体情報処理方法に関する。
従来、乳癌などの診断では、エックス線(マンモグラフィー)、音響波、MRI(核磁気共鳴法)を用いた生体情報処理装置が多く使われている。近年、パルス光を生体内に伝播させ、その伝播光の吸収によって発生した音響波(典型的には超音波)を生体表面で検知し、生体内の初期発生圧力分布あるいは光吸収係数値分布をイメージングする生体情報イメージング装置が注目を浴びている。
生体内の光吸収体の画像を取得するためには、受信した音響波の信号から、画像を再構成する必要がある。画像再構成の手法としては、上記信号の逆投影によって画像を構成する方法(非特許文献1;逆投影法)や、反復的に音響波源を推定し、上記信号との最小二乗解を求めて画像を構成する方法(非特許文献2;反復法)等が提案されている。
Photoacoustic imaging in biomedicine, REVIEW OF SCIENTIFIC INSTRUMENTS Vol.77, 041101, (2006) Comparison of iterative reconstruction approaches for photoacoustic tomography, Proc. of SPIE Vol. 6437 64370(2007)
非特許文献1によると、通常の光音響波イメージング装置においては、閉球体表面上、無限平面上、または、無限長を有する円柱表面上の全ての位置(方向)で音響波を検出することで、完全に生体情報を画像化できる。
しかしながら、生体を診断する場合には、上記のような方法で音響波を検出することはできない(生体周りの全方向で音響波を検出することはできない)。そのため、例えば、生体の一部に音響波検出器を配置して、音響波を検出することになる。この場合、光吸収体と音響波検出器との距離が大きければ大きい程、音響波検出器と光吸収体とがなす立体角は小さくなり、取得する情報(生体の内部情報;音響波源の分布の情報)の不確実性が増してしまう。
このように、取得する情報の不確実性が増加すると、例えば、非特許文献1で示された逆投影法により画像再構成を行った場合に、実際には存在しない虚像(アーティファクト)が発生してしまう。さらに、解像度(横解像度)が低下し、診断精度が低下してしまう。
一方、反復法によれば、生体周りの全方向で音響波を検出できない場合においても、アーティファクトの発生や横解像度の低下を抑制できることが報告されている。しかしながら、非特許文献2に開示の方法は、反復計算の度に、解析解を用いた音響波の算出が行われるため、処理負荷が大きい。また、そのような方法において、反復計算の度に、受信素子のサイズや受信することのできる音響波の周波数帯域など、実際の測定条件を取り込むことは、きわめて非効率的である。実際の測定条件を考慮したシミュレーションを実施することにより音響波を算出する場合においても、計算時間やメモリが膨大になってしまう(三次元の画像を再構成することが困難となる)。
本発明は、生体周りの全方向で音響波を検出できない場合において、簡易な方法で正確な音響波源の分布を推定することのできる技術を提供することを目的とする。
本発明の被検体情報処理装置は、被検体に光を照射する光源と、光吸収により発生する音響波を受信し、信号を出力する複数の受信素子と、予め定められた複数の位置のそれぞれについて、その位置に音響波源が存在すると仮定したときに前記複数の受信素子から出力される信号を表す情報を記憶する記憶手段と、前記複数の受信素子から出力された信号と、前記記憶手段に記憶されている情報とから、前記被検体内の音響波源の分布を推定する推定手段と、を有し、前記記憶手段に記憶されている前記情報は、第1の位置に音響波源が存在すると仮定したときに前記複数の受信素子から出力される第1の信号と、第2の位置に音響波源が存在すると仮定したときに前記複数の受信素子から出力される第2の信号とを、前記複数の受信素子の位置と前記複数の受信素子が音響波を受信する時刻とから構成されている信号空間上にて平行移動または回転移動することによって一致する信号を共通化した情報であることを特徴とする。
本発明の被検体情報処理方法は、予め定められた複数の位置のそれぞれについて、その位置に音響波源が存在すると仮定したときに複数の受信素子から得られる信号を表す情報を記憶する記憶手段を有する被検体情報処理装置における被検体情報処理方法であって、被検体に光を照射することにより被検体内で発生する音響波を前記複数の受信素子で受信して、信号を出力するステップと、前記複数の受信素子から出力された信号と、前記記憶部に記憶されている情報とから、前記被検体内の音響波源の分布を推定するステップと、を有し、前記記憶手段に記憶されている前記情報は、第1の位置に音響波源が存在すると仮定したときに前記複数の受信素子から出力される第1の信号と、第2の位置に音響波源が存在すると仮定したときに前記複数の受信素子から出力される第2の信号とを、前記複数の受信素子の位置と前記複数の受信素子が音響波を受信する時刻とから構成されている信号空間上にて平行移動または回転移動することによって一致する信号を共通化した情報であることを特徴とする。
本発明によれば、生体周りの全方向で音響波を検出できない場合において、簡易な方法で正確な音響波源の分布を推定することができる。
本実施形態に係る生体情報イメージング装置の構成の一例を示す図。 音響波信号の一例を示す図。 再構成画像を表示する際の処理の流れの一例を示すフローチャート。 信号処理部による処理の流れの一例を示すフローチャート。 LUTデータの一例を示す図。 音響波源分布の推定に用いる方程式の一例をマトリクス形式で表した図。 省メモリ化を実現するための方法の一例を示す図。 小ヤコビアン行列の一例を示す図。 対称性を利用して用意された小ヤコビアン行列群の一例を示す図。 対称性を利用して用意された小ヤコビアン行列群の一例を示す図。 処理の高速化を実現するための方法の一例を示す図。 コンパクト化による処理の高速化の効果を示す図。 受信素子の配列パターンと得られる再構成画像の一例を示す図。 シミュレーションによる音響波源分布の推定の結果の一例を示す図。
以下、本実施形態にかかる生体情報処理装置および生体情報処理方法について、図面を参照しながら具体的に説明する。なお、以下の実施形態は一例であり、本発明はこれに限定されるものではない。また、本発明において、音響波とは、音波、超音波、光音響波と呼ばれるものを含み、被検体に近赤外線等の光(電磁波)を照射して被検体内部で発生する弾性波のことを示す。
(装置構成)
まず、本実施形態に係る生体情報処理装置の概要について説明する。特に、生体情報処理装置の例として、生体内の音響波源の分布を推定し、イメージングする生体情報イメージング装置について説明する。図1に、本実施形態に係る生体情報イメージング装置の構成の一例を示す。
本実施形態に係る生体情報イメージング装置は、パルス光2を被検体である生体1に照射する光源4を備える。通常、光源4から発せられたパルス光2は光ファイバや液体ライトガイドなどの光伝播装置3を通って生体1の表面に照射される。
また、本実施形態に係る生体情報イメージング装置は、光吸収により発生する音響波を受信し、信号を出力する複数の受信素子(音響波検出器9)を更に備える。具体的には、複数の受信素子は、生体1内の光吸収体8が光のエネルギーの一部を吸収して発生した音響波7を検出し電気信号に変換する。なお、複数の受信素子は、二次元アレイ状に配列していることが好ましい。音響波検出器9の前部(生体1側)には、被検体を確実に固定するための圧迫板10が設けられている。
また、本実施形態に係る生体情報イメージング装置は、音響波検出器9で得られた電気信号を解析する信号処理部6、及び、その処理信号に基づいた画像(再構成画像)を表示する画像表示部5を更に備える。
具体的には、信号処理部6は、情報生成部、記憶部、および、推定部を有する。情報生成部は、記憶部で記憶する情報を数値解析により生成する。記憶部は、予め定められた複数の位置のそれぞれについて、その位置に音響波源が存在すると仮定したときに複数の受信素子から出力される信号を表す情報を記憶する。推定部は、実際の測定において複数の受信素子から出力される信号と、記憶部に記憶されている情報とから、生体内の音響波源の分布(音響波源分布)を推定する。
画像表示部5は、液晶ディスプレイ、プラズマディスプレイ、有機ELディスプレイ、電子放出素子を有するディスプレイなどの画像表示装置である。
次に、本実施形態に係る生体情報イメージング装置が再構成画像を出力するまでの過程を、図3のフローチャートを用いて説明する。
S301では、光源4からパルス光2が生体1の表面に照射される。S302では、音響波検出器9が、生体1内の光吸収体から発生する音響波を検出する。通常、球状の光吸収体から発生する音響波は、横軸を時間、縦軸を音響波の強度とすると、図2に示すようなN字型の波形を示す。S303では、信号処理部6が、音響波検出器9で得られた信号(音響波信号;測定信号)から、反復法により光吸収体の位置(音響波源の分布)を推定し、画像を生成(再構成)する。S304では、画像表示部5が、信号処理部6で生成さ
れた画像を出力(表示)する。
(音響波源の分布の推定方法)
次に、信号処理部6による音響波源の分布の推定方法について説明する。なお、以下では、説明を容易にするために二次元の音響波源の分布を推定する場合(二次元画像を再構成する場合)について説明するが、三次元の音響波源の分布を推定する場合(三次元画像を再構成する場合)にも同様の考え方を適用することができる。また、一般的に、信号処理部6の処理は、コンピュータの演算処理装置がソフトウェア(プログラム)を実行することで実現される。
本実施形態では、画素サイズの各音響波源から発生する音響波が生体に固有であること、及び、任意の音響波源分布のとき発生する音響波は、各音響波源が発生する音響波の重ね合わせで表現できる(音響波の線形性)ことを利用する。それにより、測定したデータ(音響波信号;測定信号)から、音響波源分布を逆問題として推定することができる。具体的に、本実施形態では、推定部が、音響波源の分布を仮定し、該仮定した分布における各音響波源の位置に対応する信号を記憶部に記憶されている情報から取得し、取得した信号を重ね合わせる。そして、重ね合わせた信号と実際の測定により得られた信号との一致度が最も高くなるような音響波源の分布を、生体内の音響波源の分布とする。なお、本実施形態では、2次元の再構成画像の1要素である“ピクセル”と、3次元の再構成画像の1要素である“ボクセル”の両方を“画素”と呼ぶ。
図4は、本実施形態に係る信号処理部6による処理の流れの一例を示すフローチャートである。
S401では、信号処理部6(情報生成部)が、画素位置毎に、その位置に音響波源が存在すると仮定し、各受信素子で測定されうる音響波信号を順問題解析により算出する。S402では、信号処理部6(記憶部)が、S401での算出結果を配列データ(LUT(Look Up Table)データ;小ヤコビアン行列)として記憶する。
順問題解析では、式1の波動方程式の解を計算することにより、各受信素子で測定されうる音響波信号を算出する。式1において、p(x,y,t)は、位置(x,y)、時刻tにおいて観測される圧力(音響波)を意味している。cは音速である。なお、音響波信号は、解析的な表式を用いて算出されてもよいし、FDM(Finite Difference Method)やFEM(Finite Element Method)を用いたシミュレーション解析によって算出されて
もよい。
Figure 0005451414
なお、この順問題解析は、受信素子のサイズ(エレメントサイズ)、受信素子で受信される音響波の周波数帯域、受信素子の指向性、音響波源から受信素子までの音速分布等を考慮して行ってもよい。例えば、数値解析において、実際の測定に用いる受信素子の受信面内(音響波を受信する面内)の複数点に到達する音響波の強度を計算し、それら強度の平均値を当該受信素子に到達する音響波の強度とすればよい。数値解析において複数の受信素子で受信される音響波の周波数帯域を、実際の測定に用いる受信素子で受信可能な音響波の周波数帯域に制限すればよい。具体的には、算出した音響波信号をフーリエ変換し、フィルター処理により、周波数空間で実際の測定に用いる受信素子で受信可能な周波数帯域に制限した後、逆フーリエ変換を行えばよい。また、数値解析において、実際の測定における音響波源から受信素子までの音速分布(例えば、生体と圧迫板を考慮した音速分布)を想定したモデルを用いればよい。
図5に、4画素×4画素の大きさの画像を再構成する場合のLUTデータの一例を示す
。図5(A)は、生体内を示す図(実空間(x,y))であり、図5(B)は、各受信素子で得られる音響波信号の一例を示す図(信号空間(x,t);受信素子の位置と受信時刻とから構成される空間)である。図中、実空間において塗りつぶされている位置は音響波源が存在する位置であり、信号空間において塗りつぶされている位置(小ヤコビアン行列の要素)は値を有する位置(要素)である。
S401では、例えば、実空間座標(1,1)に音響波が存在する場合の音響波信号として、図5(B)に示す様な音響波信号を算出する。なお、音響波は球面波として時間とともに広がるため、信号空間内では音響波信号の軌跡は二次曲線を描くことになる。また、4画素×4画素の大きさの画像を再構成する場合には、実空間座標(画素位置)と音響波信号の対応関係が4×4=16個得られることになる。
S403では、信号処理部6(推定部)が、S402で記憶されたLUTデータから、画像再構成(音響波源分布の推定)に用いる方程式の構成要素を生成する。
図6は、音響波源分布の推定に用いる方程式を、マトリクス形式で表したものである。図6の(a)は、S401で算出された画素位置毎の音響波信号を並べたものである(以下、マトリクスJ(ヤコビアン行列)と記す)。図6の(b)は、推定する音響波源分布を示している(以下、音響波源分布ΔPと記す)。図6の(c)は、音響波信号の実測値を示している(以下、音響波信号ΔPと記す)。音響波の線形性から、音響波源分布ΔPのときに測定される音響波信号は、画素位置毎の音響波信号(マトリクスJ)と音響波源分布ΔPとの積JΔPとなる。従って、正しい音響波源分布ΔPを推定した場合には、JΔPとΔPは完全に一致、または、近い値になる。これを式で表現すると式2の様になる。
Figure 0005451414
ここで、マトリクスJは一般には正方行列ではないが、両辺にJ(Jの転置行列)をかけると式3の様になる。これを正規方程式という。
Figure 0005451414
S404では、信号処理部6(推定部)が、正規方程式(式3)の解を求める。
正規方程式の解法としては、連立一次方程式の解法である従来手法(直説解法や反復解法)を用いることができる。再構成画像として三次元の画像を生成する場合の様に、音響波源分布を推定する領域のサイズが大きい場合には、例えば、共役勾配法等の反復解法を用いることが好ましい(本実施例では反復解法を用いるものとする)。正規方程式(式3)の解は、
Figure 0005451414
を最小にする最小2乗解であることが保証される。すなわち、正規方程式(式3)の解は、最小2乗解として(式2)を最もよく満たす最適解である。それにより、実際に測定された音響波信号との一致度が最も高くなるような音響波信号を得るための音響波源分布ΔPが推定される。
S405では、信号処理部6(推定部)が、式4で得られた値を評価する(反復解法の収束判定)。具体的には、式4で得られた値が所定の閾値(例えばε=1.0e−6)以
下か否かを判定する。そして、所定の閾値以下となるまで、S404の処理を繰り返す。
なお、得られた音響波信号にノイズが存在する場合には、式3の左辺の対角成分に一定の値αを加えた
Figure 0005451414
を解けばよい。この場合の解は、
Figure 0005451414
を最小にする最小2乗解であることが保証される。
以上の処理により、予め用意されたマトリクスJと実測値(音響波信号ΔP)から、音響波源分布ΔPを推定し、画像再構成を行うことができる。
なお、再構成画像のサイズが大きい場合(推定する音響波源分布の領域が大きい場合;例えば、サイズの大きな二次元や三次元の画像を再構成する場合)に、図6(a)のようにマトリクス全体をメモリに記憶することは、現実的なメモリ容量を鑑みると困難である。そこで、本願発明者らは、後述する対称性を利用して大幅な省メモリ化を実現する方法を考案した。更に、本願発明者らは、マトリクス要素の局在性を利用して計算の高速化をする方法を考案した。これにより、汎用PCなどであっても(メモリ容量や演算能力が限られていても)、サイズの大きい再構成画像を生成することが可能となる。
(省メモリ化)
まず、省メモリ化(記憶部に記憶する情報の圧縮方法)について説明する。
予め定められた複数の位置の内、第1の位置に音響波源が存在すると仮定したときに複数の受信素子から出力される信号を第1の信号とし、第2の位置に音響波源が存在すると仮定したときに複数の受信素子から出力される信号を第2の信号とする。この場合、第1の信号と、第2の信号とは、信号空間上で平行移動または回転移動することにより、互いに一致する部分を有することがある(このような性質を対称性と呼ぶ)。本実施形態では、第1の信号と第2の信号の間の一致する部分を共通化することにより、記憶部に記憶する情報のデータ量が圧縮された情報を用いる。
図7は、省メモリ化を実現するための方法の一例を示す図である。図7の方法では、上述した画素位置毎のLUTデータよりも大きいLUTデータ(図7の例では、図5(B)の二倍のサイズのLUTデータ)を用いる。そして、上記対称性を利用して、受信素子の配列方向(図中x方向)に並ぶ画素位置毎にLUTデータ内の参照する領域を変化させる。その結果、LUTデータの数を減らすことができ、省メモリ化が可能になる。
但し、図8に示すように、受信素子までの距離が異なる画素位置を音響波源とした場合の音響波信号は、それぞれ、信号空間内において音響波の軌跡の曲率が異なる。そのため、複数の受信素子までの距離が互いに異なる複数の小ヤコビアン行列(図8中y方向(深さ方向)の位置毎の小ヤコビアン行列)が必要となる。
図9は、4画素×4画素の二次元の画像を再構成する場合に用いる小ヤコビアン行列群(上記対称性を利用して用意された小ヤコビアン行列群)の一例を示す図である。なお、図9では、小ヤコビアン行列内の各要素の値の有無について省略している。
図6の(a)の場合には、小ヤコビアン行列群の総要素数は、4(受信素子の数)×6(サンプリング数)×16(再構成画像の画素数)=384個となる。一方、図9の例では、2×(受信素子の数)×6(サンプリング数)×4(深さ方向の画素数)=192個
の要素で済む。したがって、図9の場合には、図6の(a)の場合に比べ、必要なメモリサイズを2分の1にすることができる。この省メモリ化の効果は、音響波源分布を推定する領域のサイズが大きい場合に顕著になる。例えば、64画素×64画素の再構成画像を生成する場合には、対称性を利用することにより、画素位置毎の小ヤコビアン行列を用いる場合よりも、必要なメモリサイズを32分の1にすることができる。
図10は、4画素×4画素×4画素の三次元の画像を再構成する場合に用いる小ヤコビアン行列群(上記対称性を利用して用意された小ヤコビアン行列群)の一例を示す図である。なお、図10は、三次元の画像の1つの面(4画素×4画素)に対応する16個の受信素子を用い、サンプリング数を6とした場合の例である。また、図10では、図9と同様に、小ヤコビアン行列内の各要素の値の有無について省略している。
図10から、上記対称性を利用した小ヤコビアン行列を用いる場合の方が、画素位置毎の小ヤコビアン行列を用いる場合よりも小ヤコビアン行列群の総要素数(必要なメモリサイズ)が少ないことがわかる。例えば、64画素×64画素×64画素の再構成画像を生成する場合には、上記対称性を利用することにより、画素位置毎の小ヤコビアン行列を用いる場合よりも、必要なメモリサイズを1024分の1にすることができる。それにより、従来は扱えなかったサイズの領域について音響波源分布を推定することが可能になる。
この様に、省メモリ化の効果は、音響波源分布を推定する領域のサイズが大きいほど顕著になり、特に三次元の音響波源分布を推定する際により顕著になる。
(処理の高速化)
次に、処理の高速化について説明する。
図11は、計算(処理)の高速化を実現するための方法の一例を示す図である(図中、数値0〜5はtの値を意味する)。受信素子で検出される音響波はパルス波であるため、小ヤコビアン行列の要素の大部分は、実質的に値を有していない。そこで、本実施形態では、小ヤコビアン行列において、値を有していない領域を圧縮する(省く)ことにより、計算を高速化する。具体的には、記憶部に記憶する情報として、予め定められた複数の位置のそれぞれについて、その位置に音響波源が存在すると仮定したときの音響波を受信する受信素子とその受信時刻との組み合わせ(小ヤコビアン行列内での位置)を表す情報を用いる。以下、このような小ヤコビアン行列の圧縮を“コンパクト化”と呼ぶ。
図12は、小ヤコビアン行列のコンパクト化による処理の高速化の効果を模式的に示す。図12の例では、コンパクト化前の小ヤコビアン行列は8×6の要素を有する。そのため、そのようなヤコビアン行列を用いる場合には、参照する要素の数は48個となる。一方、コンパクト化後の小ヤコビアン行列では、値のある要素のみが参照されるため、参照する要素の数が9個となる。それにより、小ヤコビアン行列の要素を参照する際の演算量を48分の9に減少することができ、処理を高速化することができる。
なお、小ヤコビアン行列の要素を値の大きい順に並べ替え、復元の際に所定値以上の大きさの値のみを値の大きい順に参照すれば、精度は低下するが演算量をより低減させることができる(処理を高速化することができる)。
図13は、受信素子の配列パターンと得られる再構成画像の一例を示す図である。図13(A)は、二次元的に配列された複数の受信素子で測定したデータから、三次元の画像を再構成する(三次元の音響波源分布を推定する)例を示している。図13(B)は、二次元的に配列された複数の受信素子で測定したデータから、二次元の画像を再構成する(二次元の音響波源分布を推定する)例を示している。図13(C)は、一次元的に配列された複数の受信素子で測定したデータから、二次元の画像を再構成する(二次元の音響波源分布を推定する)例を示している。
上述した対称性を利用した方法を用いれば、図13(A)〜(C)のいずれの場合においても省メモリ化の効果を得ることができる。特に、図13(A)に示す場合において高
い省メモリ化の効果を得ることができる。なお、図13の例では、立方体の1つの面に複数の受信素子が配置されているが、複数の受信素子は立方体の複数面に配置されていてもよい。
なお、本実施形態では、記憶部に記憶する情報を数値解析により生成するものとしたが、そのような情報は、生体を模擬したサンプルを用いて実際に測定された音響波の測定値を用いることにより生成してもよい。具体的には、生体内の音速と圧迫板の音速を模擬するサンプルを作製し、その中に単位画素の大きさに相当する音響波源を配置すればよい。そして、そのようなサンプルを用いて、音響波(配置した音響波源から発生する音響波)を測定することにより、小ヤコビアン行列を生成してもよい。
また、本実施形態では、信号処理部6が、記憶部に記憶する情報を生成する情報生成部を有する構成としたが、情報は別の装置などで生成され、記憶部に予め記憶されていてもよい。
また、本実施形態では、対称性が並進対称性である場合(受信素子の配列方向が直線的である場合)について説明したが、対称性は軸対称性や球対称性であっても、それらの対称性を利用することによって省メモリ化を実現できる。ただし、対称性が軸対称性や球対称性の場合には、受信素子の位置や形状と、測定される信号との対応関係の情報が必要となる。
<実施例1>
以下、シミュレーションによる画像の再構成(音響波源分布の推定)の結果について図14を用いて説明する。
なお、本実施例では64画素×64画素の再構成画像を生成した。図14において、各画像の左上の座標を(0,0)、右下の座標を(63,63)とする。
図14(A)は、実測値(音響波信号ΔP)のデータを生成するための音響波源分布を示している。本実施例では、座標(15,15),(31,31),(47,47)に音響波源を設定した。また、受信素子は4辺の内の1辺(下辺)にのみ配置した。受信素子の数は64個とした。そして、サンプリング数を256として実測値(音響波信号ΔP)のデータを生成した。
図14(B)は、従来の逆投影法により、この実測値(音響波信号ΔP)のデータから得られた再構成画像を示す。図14(B)に示すように、従来の逆投影法では、はっきりとアーティファクトが観測され、元の音響波源分布を正しく推定することができていないことがわかる。また、像のコントラストも低く不鮮明である。これらは、信号を取得する受信素子が被検体周りの全方向に存在しないために生じる(非特許文献1)。
図14(C)は、本発明の方法(順問題解析結果の利用、対称性を利用した省メモリ化、及び、小ヤコビアン行列のコンパクト化)によって得られた再構成画像を示す。図14(C)に示すように、本実施例では、予め用意された順問題解析結果を用いることにより、従来の逆投影法を用いた場合(図14(B))と比較して、受信素子が被検体周りの全方向に存在しない場合でも正確な音響波源分布を推定することができた。
また、従来の反復法によって、図14(C)のように正確な音響波源分布の推定が可能な領域のサイズは、メモリ容量や計算速度等の都合上かなり限定されていた。
例えば、汎用PCを用いて、従来の反復法により音響波源分布を計算すると、10分の計算時間を要した。それに対し、本実施例では、小ヤコビアン行列をコンパクト化しているため、同等のPCを用いた場合に5秒程度(従来の120分の1程度)の計算時間で済んだ。
また、本実施例では、対称性を利用した小ヤコビアン行列を用いているため、画素位置
毎の小ヤコビアン行列を用いる場合に比べ、使用するメモリ容量を32分の1に低減することができた。
以上の結果から、本発明が診断装置の実用化に大きく寄与していることがわかる。
<実施例2>
64画素×64画素×64画素の再構成画像を生成することにより、省メモリ化の効果を検証した。なお、本実施例では、三次元の画像の1つの面(64画素×64画素)に対応する4096個の受信素子を用い、サンプリング数を512とした。
この場合、対称性を利用しない場合に必要となるメモリの容量は、小ヤコビアン行列群の総要素数である、64(受信素子数)×512(サンプリング数)×64(画素数)に比例することとなる。一方、対称性を利用する場合に必要となるメモリの容量は、4×64(受信素子の数)×512(サンプリング数)×64(深さ方向の画素数)に比例することとなる。すなわち、本実施例によれば、対称性を利用することにより、利用しない場合のメモリ容量の1024分の1のメモリ容量で同じ計算を行うことができる。
<実施例3>
64画素×64画素×64画素の再構成画像を生成することにより、小ヤコビアン行列のコンパクト化による処理の高速化の効果を検証した。なお、本実施例では、三次元の画像の1つの面(64画素×64画素)に対応する4096個の受信素子を用い、サンプリング数を512とした。
この場合、小ヤコビアン行列のコンパクト化により、小ヤコビアン行列の参照される要素の数は、4×(64×64×512)から、4×(64×64×30)に減少する。すなわち、参照される要素の数は、コンパクト化前に比べ512分の30となる。その結果、小ヤコビアン行列をコンパクト化するための処理を増やしたとしても、トータルで15倍程度、処理を高速化することができる。
以上述べたように、本実施形態に係る生体情報処理装置および生体情報処理方法によれば、各音響波源から発生する音響波が生体に固有であること、及び、音響波の線形性を利用して、簡易な方法で正確な音響波源の分布を推定することができる。具体的には、音響波源の分布を仮定し、その分布における各音響波源の位置に対応する信号を重ね合わせる。そして、重ね合わせた信号と実際の測定により得られた信号との一致度が最も高くなるような音響波源の分布が生体内の音響波源の分布とする。それにより、正確な音響波源の分布を推定することができる。また、予め定められた複数の位置に対応する信号が記憶されており、実際の測定で得られた信号はそれらの重ね合わせで表現される。そのため、従来のように、反復計算の度に、解析解を用いた音響波の算出を行う必要はなく、簡易に音響波源の分布を推定することができる。
4・・・光源 6・・・信号処理部 9・・・音響波検出器

Claims (10)

  1. 被検体に光を照射する光源と、
    光吸収により発生する音響波を受信し、信号を出力する複数の受信素子と、
    予め定められた複数の位置のそれぞれについて、その位置に音響波源が存在すると仮定したときに前記複数の受信素子から出力される信号を表す情報を記憶する記憶手段と、
    記複数の受信素子から出力された信号と、前記記憶手段に記憶されている情報とから、前記被検体内の音響波源の分布を推定する推定手段と、
    を有し、
    前記記憶手段に記憶されている前記情報は、第1の位置に音響波源が存在すると仮定したときに前記複数の受信素子から出力される第1の信号と、第2の位置に音響波源が存在すると仮定したときに前記複数の受信素子から出力される第2の信号とを、前記複数の受信素子の位置と前記複数の受信素子が音響波を受信する時刻とから構成されている信号空間上にて平行移動または回転移動することによって一致する信号を共通化した情報であることを特徴とする被検体情報処理装置。
  2. 前記推定手段は、音響波源の分布を仮定し、該仮定した分布における各音響波源の位置に対応する信号を前記記憶手段に記憶されている情報から取得し、取得した信号を重ね合わせ、重ね合わせた信号と実際の測定により得られた信号との一致度が最も高くなるような音響波源の分布を、前記被検体内の音響波源の分布とする
    ことを特徴とする請求項1に記載の被検体情報処理装置。
  3. 前記記憶手段に記憶されている前記情報は、前記予め定められた複数の位置のそれぞれについて、その位置に音響波源が存在すると仮定したときの音響波を受信する前記複数の受信素子の位置前記複数の受信素子のそれぞれが音響波を受信する時刻との組み合わせを表す情報である
    ことを特徴とする請求項1または2に記載の被検体情報処理装置。
  4. 前記記憶手段で記憶する前記情報を数値解析により生成する生成手段をさらに有する
    ことを特徴とする請求項1〜3のいずれか1項に記載の被検体情報処理装置。
  5. 前記生成手段は、前記数値解析において、前記複数の受信素子のそれぞれについて、
    信素子の受信面内の複数点に到達する音響波の強度を計算し、該音響波の強度の平均値を該受信素子に到達する音響波の強度とする
    ことを特徴とする請求項4に記載の被検体情報処理装置。
  6. 前記生成手段は、前記数値解析において前記複数の受信素子で受信される音響波の周波数帯域を、前記複数の受信素子受信可能な音響波の周波数帯域に制限することを特徴とする請求項4または5に記載の被検体情報処理装置。
  7. 前記生成手段は、前記数値解析において、前記音響波源から前記複数の受信素子までの音速分布を想定したモデルを用いる
    ことを特徴とする請求項4〜6のいずれか1項に記載の被検体情報処理装置。
  8. 前記記憶手段で記憶する前記情報を、前記被検体を模擬したサンプルを用いて実際に測定された音響波の測定値を用いることにより生成する生成手段をさらに有する
    ことを特徴とする請求項1〜3のいずれか1項に記載の被検体情報処理装置。
  9. 予め定められた複数の位置のそれぞれについて、その位置に音響波源が存在すると仮定したときに複数の受信素子から得られる信号を表す情報を記憶する記憶手段を有する被検体情報処理装置における被検体情報処理方法であって、
    被検体に光を照射することにより被検体内で発生する音響波を前記複数の受信素子で受信して、信号を出力するステップと、
    記複数の受信素子から出力された信号と、前記記憶手段に記憶されている情報とから、前記被検体内の音響波源の分布を推定するステップと、
    を有し、
    前記記憶手段に記憶されている前記情報は、第1の位置に音響波源が存在すると仮定したときに前記複数の受信素子から出力される第1の信号と、第2の位置に音響波源が存在すると仮定したときに前記複数の受信素子から出力される第2の信号とを、前記複数の受信素子の位置と前記複数の受信素子が音響波を受信する時刻とから構成されている信号空間上にて平行移動または回転移動することによって一致する信号を共通化した情報であることを特徴とする被検体情報処理方法。
  10. コンピュータに、請求項9に記載の被検体情報処理方法を実行させるためのプログラム。
JP2010008366A 2010-01-18 2010-01-18 被検体情報処理装置および被検体情報処理方法 Expired - Fee Related JP5451414B2 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
JP2010008366A JP5451414B2 (ja) 2010-01-18 2010-01-18 被検体情報処理装置および被検体情報処理方法
US12/975,423 US9289126B2 (en) 2010-01-18 2010-12-22 Subject information obtaining apparatus, subject information obtaining method, and program
EP10197290A EP2345364A1 (en) 2010-01-18 2010-12-29 Subject information obtaining apparatus, subject information obtaining method, and program

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2010008366A JP5451414B2 (ja) 2010-01-18 2010-01-18 被検体情報処理装置および被検体情報処理方法

Publications (3)

Publication Number Publication Date
JP2011143175A JP2011143175A (ja) 2011-07-28
JP2011143175A5 JP2011143175A5 (ja) 2013-02-28
JP5451414B2 true JP5451414B2 (ja) 2014-03-26

Family

ID=43904005

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2010008366A Expired - Fee Related JP5451414B2 (ja) 2010-01-18 2010-01-18 被検体情報処理装置および被検体情報処理方法

Country Status (3)

Country Link
US (1) US9289126B2 (ja)
EP (1) EP2345364A1 (ja)
JP (1) JP5451414B2 (ja)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2013116809A1 (en) * 2012-02-03 2013-08-08 Los Alamos National Security, Llc Ultrasound waveform tomography with tv regularization
JP5762995B2 (ja) * 2012-02-28 2015-08-12 富士フイルム株式会社 光音響画像生成装置及び方法
JP6168802B2 (ja) * 2012-04-12 2017-07-26 キヤノン株式会社 処理装置、処理方法、及びプログラム
EP2732756B1 (en) 2012-11-15 2019-09-11 Canon Kabushiki Kaisha Object information acquisition apparatus
US20140182383A1 (en) 2012-12-28 2014-07-03 Canon Kabushiki Kaisha Object information obtaining device, display method, and non-transitory computer-readable storage medium
CN108261209B (zh) * 2018-01-23 2021-07-23 湖南大学 改进型的高分辨声聚焦光声内窥成像反投影成像的方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP4739363B2 (ja) * 2007-05-15 2011-08-03 キヤノン株式会社 生体情報イメージング装置、生体情報の解析方法、及び生体情報のイメージング方法
WO2009073979A1 (en) 2007-12-12 2009-06-18 Carson Jeffrey J L Three-dimensional photoacoustic imager and methods for calibrating an imager
JP5284129B2 (ja) * 2008-02-06 2013-09-11 キヤノン株式会社 イメージング装置、及び解析方法

Also Published As

Publication number Publication date
JP2011143175A (ja) 2011-07-28
US20110178739A1 (en) 2011-07-21
EP2345364A1 (en) 2011-07-20
US9289126B2 (en) 2016-03-22

Similar Documents

Publication Publication Date Title
Waibel et al. Reconstruction of initial pressure from limited view photoacoustic images using deep learning
US20190133451A1 (en) Apparatus for acquiring biofunctional information, method for acquiring biofunctional information, and program therefor
JP5451414B2 (ja) 被検体情報処理装置および被検体情報処理方法
Li et al. Refraction corrected transmission ultrasound computed tomography for application in breast imaging
Huang et al. Bezier interpolation for 3-D freehand ultrasound
JP6504826B2 (ja) 情報処理装置および情報処理方法
CN104510495B (zh) 被检体信息获取装置及其控制方法
Wen et al. An accurate and effective FMM-based approach for freehand 3D ultrasound reconstruction
CN103025248B (zh) 图像信息获取装置和图像信息获取方法
Sandhu et al. 3D frequency-domain ultrasound waveform tomography breast imaging
JP6929365B2 (ja) 波形インバージョンを用いる非侵襲性医療撮像のための方法および装置
US20170032702A1 (en) Method and Apparatus For Generating an Ultrasound Scatterer Representation
Poudel et al. Mitigation of artifacts due to isolated acoustic heterogeneities in photoacoustic computed tomography using a variable data truncation-based reconstruction method
Perez-Liva et al. Speed of sound ultrasound transmission tomography image reconstruction based on Bézier curves
Kretzek et al. GPU-based 3D SAFT reconstruction including attenuation correction
US11940531B2 (en) Crosstalk-free source encoding for ultrasound tomography
US20170238898A1 (en) Device, system, and method for hemispheric breast imaging
Fincke et al. Towards ultrasound travel time tomography for quantifying human limb geometry and material properties
KR20160054360A (ko) 영상 장치 및 영상화 방법
US20190313984A1 (en) Elastography based on x-ray computed tomography and sound wave integration
Zhou et al. Frequency-domain full-waveform inversion-based musculoskeletal ultrasound computed tomography
Zhu et al. Mitigating the limited view problem in photoacoustic tomography for a planar detection geometry by regularised iterative reconstruction
Peters et al. Digital image–based elasto–tomography: Nonlinear mechanical property reconstruction of homogeneous gelatine phantoms
Xu et al. Compact reverse time migration: A real-time approach for full waveform ultrasound imaging for breast
Matthews et al. Synergistic image reconstruction for hybrid ultrasound and photoacoustic computed tomography

Legal Events

Date Code Title Description
A521 Written amendment

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20130109

A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20130109

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20131028

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

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20131226

R151 Written notification of patent or utility model registration

Ref document number: 5451414

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R151

LAPS Cancellation because of no payment of annual fees