JP5683213B2 - 画像形成装置及び画像形成方法 - Google Patents

Description

本発明は、被検体から発生する音響波を受信して、各部位の組織や特性分布を画像データ化するための画像形成装置及び画像形成方法に関するものである。

音響波を利用して被検体内部の情報を画像データ化する技術として、超音波エコー法、光音響トモグラフィー（ＰＡＴ：Photoacoustic Tomography）などが知られている。超音波エコー法は、被検体内部に向けて超音波を送信し、被検体内部の組織で反射された超音波を受信して画像データ化する方法である。光音響トモグラフィーは、光音響効果によって被検体内部から放射される音響波を受信して画像データ化する方法である。光音響効果（Photoacoustic effect）とは、可視光、近赤外光、マイクロ波等の電磁波パルスを被検体に照射した際に、被検体が電磁波を吸収して発熱し熱膨張する事によって音響波（典型的には超音波）が発生する現象である。

超音波エコー法やＰＡＴの画像形成におけるビームフォーミング法として、タイムドメイン法とフーリエドメイン法が知られている。フーリエドメイン法は、時系列信号をフーリエ変換して、変数を時間から周波数に空間座標を波数に変換した後に処理を行う方式である（非特許文献１）。フーリエドメイン法は高速フーリエ変換を利用する事によって高速処理ができる特徴がある。またＣＰＵの高速化が進んでソフトウェア上での数値処理も実用的になり、最近では並列処理の進歩も著しいので、将来はソフトウェアによる画像処理、即ちフーリエドメイン法が有望と見られている。特にＰＡＴの画像形成装置ではフーリエドメイン法が好ましく利用される。

Minghua Xu and Lihong V. Wang：REVIEW OF SCIENTIFIC INSTRUMENTS 77巻 041101 2006年4月17日

従来、超音波エコーやＰＡＴでは主にリニア型やセクタ型等、１次元に音響変換器が配列され電子的に走査を行うタイプのプローブが用いられる事が多かった。このタイプのプローブを用いると被検体断面の２次元面像が形成できる。しかしながら、患部の見落としを少なくするためには、３次元画像の取得が効果的である。また３次元画像であれば、診断後の画像の見直しや第三者による確認も容易である。３次元画像を取得するには、１次元状に音響変換器が配列されたプローブを電子走査方向と垂直な方向に機械走査する、音響変換器が２次元状に配列されたタイプのプローブを使用する、等の方法を採用すれば良い。

しかし３次元画像データの形成に伴って処理すべきデータ量が莫大になる。特に高解像度画像が求められると実時間処理は不可能で、従来の装置を使い慣れた操作者には大きな違和感を与える。画像処理をソフトウェアで行う事が多いＰＡＴではさらに事情が深刻で、高速なフーリエドメイン法でも、現状では３次元画像データの形成には数十秒以上を要し、実時間動作とはほど遠い状態であった。

本発明は、光音響トモグラフィーや超音波エコーにおいて、フーリエドメイン法による画像データの形成処理を高速化し、高解像度の３次元画像の実時間処理を可能にするため
の技術を提供することを目的とする。

本発明の第１態様は、被検体内部より放射される音響波に基づき被検体内部の情報を画像データ化する画像形成装置であって、音響波を受信し、信号に変換する複数の音響変換器と、前記複数の音響変換器から出力された受信信号を、フーリエドメイン法で処理することによって、被検体内部の複数の部位それぞれから放射された音響波の強度を算出する信号処理部と、を備え、前記信号処理部は、前記複数の音響変換器の受信信号からある部位が放射した音響波の強度を求める演算式の中で、音響変換器の位置と、部位の位置と、音響波が前記部位から到達するのに要する時間と、波数の次元を有するフーリエ変換の変数と、で定まる項の値を予め計算し、係数として記憶している係数メモリと、前記音響変換器の受信信号に対して、対応する係数を前記係数メモリから取得し、乗算する乗算手段と、前記乗算手段の乗算結果を部位ごとに累積するボクセルメモリと、を備える画像形成装置である。

本発明の第２態様は、被検体内部の複数の部位より放射された音響波を複数の音響変換器により所定の期間に亘って受信して受信信号に変換し、該受信信号をフーリエドメイン法で処理することによって各部位が放射した音響波の強度を画像データ化する画像形成方法であって、複数の音響変換器の受信信号からある部位が放射した音響波の強度を求める演算式の中で、音響変換器の位置と、部位の位置と、音響波が前記部位から到達するのに要する時間と、波数の次元を有するフーリエ変換の変数と、で定まる項の値を予め計算し、係数として係数メモリに格納するステップと、前記音響変換器の受信信号に対して、対応する係数を前記係数メモリから取得し、乗算するステップと、前記受信信号と前記係数の乗算結果を部位ごとにボクセルメモリに累積するステップと、を含む画像形成方法である。

本発明の第３態様は、上記画像形成方法の各ステップをコンピュータに実行させるためのプログラムである。

本発明によれば、光音響トモグラフィーや超音波エコーにおいて、フーリエドメイン法による画像データ形成処理を高速化し、高解像度の３次元画像データの実時間処理が可能となる。

本発明の計算アルゴリズム１を説明する図。 本発明の計算アルゴリズム２を説明する図。 本発明の計算アルゴリズム３を説明する図。 本発明の第１の実施形態のシステム構成を示す図。 本発明の第２の実施形態のシステム構成を示す図。 本発明の第３の実施形態のシステム構成を示す図。 フーリエドメイン法の原理説明における座標系を示す図。 フーリエドメイン法の従来の計算アルゴリズムを説明する図。

本発明の方法を説明する前に、現状のフーリエドメイン法によるＰＡＴの画像データ形成のアルゴリズムとその問題点について、Wangの解説（非特許文献１）に従って説明する。Wangは被検体を取り囲む音響変換器が配列される面として、球面、平面、円筒面の３種類をあげて説明しているが、ここでは説明が簡単な平面の場合を取り上げる。なお、説明の都合上非特許文献１とは異なる記号を用いる場合がある。

図７の１０１は音響変換器が配列される平面である。この平面１０１に含まれるｘ軸と
ｙ軸は直交しており、ｚ軸は平面１０１に垂直であり、この平面１０１よりｚが正の領域に被検体が存在する。平面１０１上にＭ個の音響変換器が配列され、ｚ＞０の領域内のＮ個の部位について測定を行うものとする。１０２はｉ番目の音響変換器を示している（i=1,2,・・・,M）。１０３は被検体内部のｊ番目の部位を示している（j=1,2,・・・,N）。被検体に各種のパルスレーザーやマグネトロン管からの電磁パルスが照射されて部位ｊに吸収されると、部位ｊは発熱し急激に膨張して音響波１０４を放射する。音響波の初期圧力をp₀(x_j,y_j,z_j)とする。(x_j,y_j,z_j)はｊ番目の部位の座標である。部位の大きさにも依存するが、電磁パルスの照射で発生した音響波（光音響波と呼ぶ。）は一般に周波数が数百ｋＨｚ〜数ＭＨｚ程度の超音波である。光音響波１０４は被検体内部を伝播しｉ番目の音響変換器１０２により受信され、受信信号p_d(x_i,y_i,t^*)として出力される。(x_i,y_i)は
この音響変換器１０２の座標であり、t^*（=ct）は被検体内部の音速cで規格化した電磁パルス照射後の経過時間である。なお、被検体内部の特定部位だけに限定して電磁パルスを照射できないＰＡＴでは、p_d(x_i,y_i,t^*)は複数の部位からの光音響波が足し合わされた値となる。

以下、図８を使って従来のフーリエドメイン法のアルゴリズムをステップ毎に説明する。式の番号は非特許文献１に従う。また、変数t,u,v,w,jなどはステップ２０１の前に初
期化されているものとする。

ステップ２０１：時刻tにおいてＭ個の音響変換器にて光音響波を受信し、受信信号p_d(x_i,y_i,t^*)を出力する。
ステップ２０２：受信信号を各音響変換器に対応付けられたラインメモリに蓄積する。ここでラインメモリと称するのは、通常FIFO（First In First Out）と呼ばれる、時系列情報を配列の順に入力し、入力された順に出力するメモリの事である。
ステップ２０３：所定の期間に亘って受信とメモリへの蓄積を繰返す。

ステップ２０４：波数の次元を持つ３個のフーリエ変換の変数(u,v,w)を決定する。
ステップ２０５：所定の(u,v,w)についてk=(u²+v²+w²)^1/2で決まる変数kによって時間
に関するフーリエ変換

を行う。(14)式は無限領域の積分として表されているが、実際には有限な範囲で離散的な数値積分が行われる。以下説明で積分が現れる場合は同様に解釈する。続いて、
ステップ２０６：長さに関するフーリエ変換

を行う。(22’)は非特許文献１においては逆フーリエ変換(22)として記載されているが、ここでは説明の都合上フーリエ変換として表現した。(x_i,y_i)は音響変換器ｉの座標であ
る。(22’)式は無限領域の積分として表されているが、実際にはＭ個の音響変換器ｉに関する和として計算される。
ステップ２０７，２０８：ステップ２０５と２０６を予め設定した全ての(u,v,w)の組
に対して実行し、

の組が得られる。ここで、sgn(w)は、w>0のとき1、w=0のとき0、w<0のとき-1をとる関数
である。

ステップ２０９：得られたＰ₀(u,v,w)の組を用いて(u,v,w)に関する逆フーリエ変換

を行うと、部位ｊに対する初期圧力p₀(x_j,y_j,z_j)が求められる。
ステップ２１０：さらに被検体内のＮ個の部位に対して逆フーリエ変換(21)を繰り返すと被検体全体が画像データ化できる。

次に上記全工程で必要となる掛け算の回数を概算する。音響変換器をＭ個、見ようとする被検体内の部位を１次元あたりＬ個としてＬ^３（＝Ｎ）個、音響波の受信回数をＴ回、フーリエ変換で使用する３個の変数(u,v,w)の組をＦ^３個（１次元あたりＦ個）とする。
ある(u,v,w)に対しステップ２０５でＴ回、ステップ２０６でＭ回、両者を組み合わせる
とＴ×Ｍ回、さらにこの処理を(u,v,w)の全ての組に対して行うとＴ×Ｍ×Ｆ^３回の掛け
算が必要となる。またある(u,v,w)に対してステップ２０９がＦ^３回であるから、これを
全部位に対して行うとＬ^３×Ｆ^３回となり、合わせて（Ｔ×Ｍ×Ｆ^３＋Ｌ^３×Ｆ^３）＝（Ｔ×Ｍ＋Ｌ^３）Ｆ^３回となる。精度を保つ上で、Ｆ≒Ｌ、Ｔ≒Ｌ、Ｍ≒Ｌ^２とするのが効率的と考えるので、掛け算回数は約２Ｌ^６回、１回の掛け算を１ステップで行うとすれば演算ステップ数は２Ｌ^６となる。ここで(x,y)に関しては高速フーリエ変換が適用できる
ので、このステップ数はＬ^２［Ｌ×log（Ｌ）］^２＝Ｌ^４log^２（Ｌ）程度に節約できるが、高解像度を得るためにＬを大きくすると、これでも膨大な回数となり実時間での画像再構成が困難になる。

本発明者は鋭意検討の結果、上記アルゴリズムを変形すると処理の一部が並列化でき、全体としての画像形成時間を大きく短縮出来る事を見出した。以下、本発明の画像形成装置及び画像形成方法で採用したアルゴリズムの詳細を説明する。

(21)式のＰ₀(u,v,w)に（23）式を代入すると、

となる（なお、100番台は本明細書で独自に使用する番号である）。

(101)式のＱに(22’)式を代入すると、

となる。

(102)式のＰ_dに(14)式を代入すると、

である。

(103)式は、複数の音響変換器の受信信号から部位ｊが放射した音響波の強度p₀を求め
るための演算式である。この(103)式の中の項Ａ(x_j,y_j,z_j:x_i,y_i,t^*)には、(104)式に示
すように、p_d(x_i,y_i,t^*)が含まれておらず、また(u,v,w)に関する積分は全てこの式の中
で完結している。従って、音響変換器の位置(x_i,y_i)と、被検体の中で画像データ化した
い部位の位置(x_j,y_j,z_j)（又は音響変換器と部位の相対位置）と、全部位からの音響波の到来時刻（受信時刻）が分かっていれば、(104)式により係数Ａを予め計算する事が出来
る。

（アルゴリズム１）
得られた係数Ａを用いて全部位(x_j,y_j,z_j)(j=0,1,・・・,N)に対する初期圧力p₀(x_j,y_j,z_j)を求めるためのアルゴリズムの一例を図１を用いて示す。なお以下の説明において（
）内の番号は後述の図４における装置の構成要素の番号である。係数メモリ（４０８）には、予め計算されたi,j,t^*の全ての組み合わせに対応する係数Ａが格納されている。また、ボクセルメモリ（４１０）、変数t,iなどはステップ７０１の前に初期化されている
ものとする。

ステップ７０１：選択手段（４０７）で１つの音響変換器（４０５）ｉを選択し、音響波を受信し、受信信号p_d(x_i,y_i,t^*)を出力する。
ステップ７０２：係数メモリ（４０８）から「音響変換器ｉ，時刻t^*」における各部位ｊ（j=0,1,・・・,N）に対応するＮ個の係数Ａ(x_j,y_j,z_j:x_i,y_i,t^*)を読み込み、Ｎ個の
乗算手段（４０９）で各係数Ａをそれぞれp_d(x_i,y_i,t^*)に掛ける。このＮ個の処理は同時に（並列に）実行される。
ステップ７０３：Ｎ個の乗算結果Ａ(x_j,y_j,z_j:x_i,y_i,t^*)・p_d(x_i,y_i,t^*)をＮ個のボク
セルメモリ（４１０）にそれぞれ累積（加算）する。このＮ個の処理も同時に（並列に）実行される。
ステップ７０４：Ｍ個の音響変換器ｉを順次選択し、ステップ７０１〜７０３の処理を繰り返す。その結果、各ボクセルメモリ（４１０）には、初期圧力p₀の時刻t^*における瞬時値∬Ａ(x_j,y_j,z_j:x_i,y_i,t^*)・p_d(x_i,y_i,t^*)dx_idy_iが蓄積される。
ステップ７０５：各部位からの音響波が受信される全期間に亘ってステップ７０１〜７０４の処理を繰り返すことにより、(103)式で示される各部位ｊの初期圧力p₀(x_j,y_j,z_j)
の値が各ボクセルメモリ（４１０）に累積される。これにより被検体内部の３次元画像データが形成される。

（アルゴリズム２）
図２は、画像データを形成するためのアルゴリズムの別の例を示している。図１のアルゴリズム１は、１つの音響変換器の受信信号からＮ個の部位それぞれの値を求める処理を並列に行うものであるのに対し、図２のアルゴリズム２は、Ｍ個の音響変換器のそれぞれの受信信号から１つの部位の値を求める処理を並列に行うものである。以下の説明において（ ）内の番号は後述の図６における装置の構成要素の番号である。係数メモリ（６０９）には、予め計算されたi,j,t^*の全ての組み合わせに対応する係数Ａが格納されている。また、ボクセルメモリ（６１２）、変数t,jなどはステップ８０１の前に初期化されて
いるものとする。

ステップ８０１：Ｍ個の音響変換器（６０３）にてそれぞれ音響波を受信し、受信信号p_d(x_i,y_i,t^*)を出力する。
ステップ８０２：部位ｊを選択し、係数メモリ（６０９）から「部位ｊ，時刻t^*」における各音響変換器ｉ（i=0,1,・・・,M）に対応するＭ個の係数Ａ(x_j,y_j,z_j:x_i,y_i,t^*)を
読み込む。Ｍ個の乗算手段（６０８）にて係数Ａを対応する受信信号p_d(x_i,y_i,t^*)にそれぞれ掛ける。このＭ個の処理は同時に（並列に）実行される。
ステップ８０３：Ｍ個の乗算結果Ａ(x_j,y_j,z_j:x_i,y_i,t^*)・p_d(x_i,y_i,t^*)を加算手段（
６１０）にて加算し、結果を選択手段（６１１）にて部位ｊに対応するボクセルメモリ（６１２）に累積（加算）する。
ステップ８０４：Ｎ個の部位ｊを順次選択し、ステップ８０２、８０３の処理を繰り返す。その結果、各ボクセルメモリ（６１２）には、初期圧力p₀の時刻t^*における瞬時値∬Ａ(x_j,y_j,z_j:x_i,y_i,t^*)・p_d(x_i,y_i,t^*)dx_idy_iが蓄積される。
ステップ８０５：各部位からの音響波が受信される全期間に亘ってステップ８０１〜８０４の処理を繰り返すことで、(103)式で示される各部位ｊの初期圧力p₀(x_j,y_j,z_j)の値
が、各ボクセルメモリ（６１２）に累積される。これにより被検体内部の３次元画像データが形成される。

（アルゴリズム３）
図３は、画像データを形成するためのアルゴリズムの別の例を示している。前述のアルゴリズム１、２では、時間ステップΔtごとに全部位の圧力p₀の瞬時値を算出したのに対
し、図３のアルゴリズム３では、全期間の受信信号から初期圧力p₀の値を求める処理を部位ごとに順番に行う。このアルゴリズムは送信ビームフォーカスを行う超音波エコー装置に好ましく適用できる。なお以下の説明において（ ）内の番号は後述の図６における装置の構成要素の番号である。係数メモリ（６０９）には、予め計算されたi,j,t^*の全ての組み合わせに対応する係数Ａが格納されている。また、ボクセルメモリ（６１２）、変数t,jなどはステップ９０１の前に初期化されているものとする。

ステップ９０１：Ｍ個の音響変換器（６０３）にてそれぞれ音響波を受信し、受信信号p_d(x_i,y_i,t^*)を出力する。
ステップ９０２：部位ｊを選択し、係数メモリ（６０９）から「部位ｊ，時刻t^*」における各音響変換器ｉ（i=0,1,・・・,M）に対応するＭ個の係数Ａ(x_j,y_j,z_j:x_i,y_i,t^*)を
読み込む。Ｍ個の乗算手段（６０８）にて係数Ａを対応する信号p_d(x_i,y_i,t^*)に掛ける。ここで、送信ビームフォーカスを行っている場合には、送信された音響波の焦点に位置する部位ｊが選択される。このＭ個の処理は同時に（並列に）実行される。
ステップ９０３：Ｍ個の乗算結果Ａ(x_j,y_j,z_j:x_i,y_i,t^*)・p_d(x_i,y_i,t^*)を加算手段（
６１０）にて加算し、結果を選択手段（６１１）にて部位ｊに対応するボクセルメモリｊ（６１２）に累積（加算）する。
ステップ９０４：部位ｊからの信号が受信される全期間に亘ってステップ９０１〜９０３の処理を繰り返す。その結果、部位ｊに対応するボクセルメモリｊには、(103)式で示
される部位ｊの初期圧力p₀(x_j,y_j,z_j)の値が累積される。
ステップ９０５：すべての部位ｊに対して（つまり、各部位ｊに焦点を形成する様に音響波を送信しつつ）ステップ９０１〜９０４の処理を繰り返すことで、ボクセルメモリ（６１２）内に被検体内部の３次元画像データが形成される。

一般に、超音波エコー装置では、信号強度を高めるために、各部位ｊに順次焦点を形成する様に超音波を送信し、部位ｊから反射してくる反射波を順次処理する事が多い。図３に示したアルゴリズム３によれば図６の装置でこのシーケンスを実行できる。従って、図６に示す一つの装置構成を、超音波エコー画像の生成（図３のアルゴリズム３）とＰＡＴ画像の生成（図２のアルゴリズム２）の両方に利用することができる。

（本アルゴリズムの利点）
上述したアルゴリズム１〜３に必要な掛け算の回数を概算する。ただし、係数Ａは予め計算され、係数メモリに格納されているものとする。アルゴリズム１では、ステップ７０２のＮ回の掛け算をステップ７０４のループでＭ回繰り返し、さらにステップ７０５のループでＴ回繰り返す。よって、掛け算の回数はＮ×Ｍ×Ｔ回となり、前述した様にＮ＝Ｌ^３、Ｔ×Ｍ≒Ｌ^３とおくと、約Ｌ^６回となる。同様にアルゴリズム２では、ステップ８０２のＭ回の掛け算をステップ８０４のループでＮ回繰り返し、さらにステップ８０５のループでＴ回繰り返すので、掛け算の回数はＭ×Ｎ×Ｔ回≒Ｌ^６回である。また、アルゴリズム３では、ステップ９０２のＭ回の掛け算をステップ９０４のループでＴ回、ステップ９０５のループでＮ回繰り返し、掛け算の回数はＭ×Ｔ×Ｎ回≒Ｌ^６回である。

掛け算の回数（Ｌ^６）は従来のフーリエドメイン法の回数（２Ｌ^６）に比べ約半分に削減されている。さらに、本アルゴリズムの注目すべき点は、ステップ７０２、８０２、９０２の処理が並列実行し易い形になっていることにある。例えば、アルゴリズム１においてステップ７０２のＮ個の掛け算を並列処理すると、演算ステップ数はＬ^３となる。またアルゴリズム２、３においてステップ８０２、９０２のＭ個の掛け算を並列処理すると、演算ステップ数はＬ^４となる。これらのステップ数は、従来のフーリエドメイン法で高速フーリエ変換を利用した場合のステップ数Ｌ^４log^２（Ｌ）に比べてもはるかに少ない。
従って、本アルゴリズムによれば、フーリエドメイン法による画像形成処理を大幅に高速化することができ、高解像度の３次元画像データの実時間処理が可能となる。

なお、音響変換器の位置(x_i,y_i)はプローブの構造に応じて決まっており、測定部位の
座標(x_j,y_j,z_j)や測定の時間ステップt^*についても複数のパターンの中から用途や要求解像度に応じたパターンを選択すれば良いので事前に決定可能である。従って、係数Ａは予め計算しておくことが可能である。

以上の説明では、音響変換器に到着した光音響波の受信信号p_d(x_i,y_i,t^*)を完全に実時間処理することを前提にしている。しかし全部位から音響波が到着し終わるまでの100μ
秒程度の時間は人間の感覚からすれば極端に短いので、各音響変換器ｉの出力を一旦バッファメモリに蓄積し、受信の際より長い時間間隔で、受信した順番に出力して上記の処理を行う事もできる。この場合、例えば時間間隔を１０００倍に長くしても１００ｍ秒で処理が行えるので、人間には殆ど実時間動作と感じられる。その一方で、信号処理（例えばステップ７０１〜７０５）に対して１０００倍の時間的余裕を与えることができる。なお、この場合バッファメモリの内容は一度読み出すだけで良く基本的には繰り返し使用する必要はない。

上記アルゴリズムを実施する際に、係数Ａを蓄積するためのメモリ（係数メモリ）の容量が膨大になることが考えられる。係数ＡはＴ×Ｍ×Ｎ（≒Ｌ^６）個あるから、例えば、１次元について１０００点のデータを取得する場合は（Ｌ＝１０００）、１つの係数Ａを４バイトとすると、４×１０^１８バイトものメモリ容量が必要になる。しかし、以下の手段で係数メモリの容量を減らす事ができる。

手段１）係数Ａの（104）式は(u,v)について対称である。従って(u,v)が張る平面の半
分について積分を行えば、（104）式はその２倍として求められる。
手段２）被検体の内部で画像データ化しようとする部位のx,y配列と音響変換器のx,y配列を一致させると、これらの変数は(104)式の中で(x_i-x_j)や(y_i-y_j)という相対値として
現れてくるだけなので、Ｌ^２個の係数Ａが同じ値になる。言い換えると、１つの係数Ａを、(x_i-x_j)及び(y_i-y_j)が等しいＬ^２個の計算に利用するようにすることで、メモリ容量を１／Ｌ^２に削減できる。

なお以上の説明は、簡単のため音響変換器が平面状に配列された場合について行った。しかし本発明の趣旨は、音響変換器が平面以外の曲面に沿って配列される場合にも適用できる。

次に、本発明の効果を示すために、発明の実施に好適な画像形成装置のシステム構成とその動作について説明する。以下の画像形成装置は、被検体内部より放射される音響波に基づき被検体内部の情報を画像データ化する装置であり、光音響トモグラフィー（ＰＡＴ）や超音波エコー診断装置などに利用されるものである。画像形成装置により生成された画像データは不図示の表示部（モニタ）により３次元画像として表示される。

＜第１の実施形態＞
図４は、第１の実施形態の画像形成装置のシステム構成を示している。本実施形態は、上述したアルゴリズム１を光音響トモグラフィー（ＰＡＴ）の画像データ形成に適用した例である。本システムは、概略、電磁波パルス源４０３と、Ｍ個の音響変換器４０５をもつ音響波プローブと、音響変換器４０５で所定の期間に亘って受信された光音響波の受信信号から被検体内部の各部位から放射された光音響波の強度を求める信号処理部とを備える。本実施形態の信号処理部は、選択手段４０７、係数メモリ４０８、Ｎ個の乗算手段４０９、Ｎ個のボクセルメモリ４１０を備える。このシステムにより、被検体４０１の内部のＮ個の部位ｊ４０２の初期圧力を求めて、３次元画像データを生成する。

電磁波パルス源４０３は、パルス幅数ｎ秒〜数十ｎ秒程度の、可視光や波長700〜1160
ｎｍ程度の近赤外光や周波数2.45ＧＨｚのマイクロ波のパルスを放射する電磁波パルス照射手段である。このパルスが被検体４０１に照射され、そのエネルギーの一部が部位ｊに吸収されると、光音響効果により部位ｊから光音響波４０４が放射される。光音響波４０４はＭ個の音響変換器４０５によって電気信号（受信信号）に変換される。図４では音響変換器４０５は１次元に配列されている様に描かれているが、実際は音響変換器４０５は２次元に配列されている。図示されていないが、各音響変換器にはプレアンプが接続されていてもよい。また音響変換器にＡ／Ｄ変換手段が接続されていても良い。さらに音響変換器にサンプルホールド手段が接続されていても良い。以下これらの手段をも含んで、単に音響変換器と呼ぶ。

各音響変換器４０５の出力は選択手段４０７に接続される。ある時刻t^*において、選択手段４０７はＭ個の音響変換器４０５の中からｉ番目を選択してその出力である受信信号p_d(x_i,y_i,t^*)を、Ｎ個の部位４０２に対応付けられたＮ個の乗算手段４０９の全てに同時に出力する。これに同期して、ｊ番目の乗算手段４０９には、係数メモリ４０８から係数
Ａ(x_j,y_j,z_j:x_i,y_i,t^*)が読み出される。乗算手段４０９は音響変換器４０５の出力p_dと
係数Ａの積Ａ(x_j,y_j,z_j:x_i,y_i,t^*)・p_d(x_i,y_i,t^*)を計算する。この結果は各々の乗算手
段に対応するボクセルメモリ４１０に蓄積される。Ｍ個の音響変換器４０５が順に処理され、時刻t^*における∬Ａ(x_j,y_j,z_j:x_i,y_i,t^*)・p_d(x_i,y_i,t^*)dx_idy_iが計算される。さら
に次の時刻t^*+Δt^*の音響変換器の信号p_d(x_i,y_i,t^*+Δt^*)を同様に処理する。所定の期間に亘り以上の処理を繰り返すと、全ボクセルに対して(103)式で示される初期音圧p₀(x_j,y_j,z_j)が得られる。

本実施形態によれば、Ｎ個の乗算を並列処理できるので、演算ステップ数はＭ×Ｔ（＝Ｌ^３）回となる。従来のフーリエドメイン法で同じ演算を行うにはＬ^６ステップが必要となることと比較すると、本実施形態の構成により大幅に計算時間が短縮されることがわかる。

＜第２の実施形態＞
図５は、第２の実施形態のシステム構成を示している。本実施形態のシステムは、音響変換器４０５と選択手段４０７の間に、音響変換器４０５の出力信号をバッファするラインメモリ４０６を有している。その他の構成要素は第１の実施形態のものと共通であるため、図４と同じ符号で図示している。

ラインメモリ４０６はＦＩＦＯメモリであり、各音響変換器４０５から出力される受信信号を順次蓄積し、蓄積された受信信号を受信（入力）された順番に選択手段４０７へ出力する。選択手段４０７以降の機能は第１の実施形態と同じである。ラインメモリ４０６が信号を出力する間隔は適宜設定できる。よって、ラインメモリ４０６に信号を蓄積する第１の時間間隔よりもラインメモリ４０６から信号を出力する第２の時間間隔を長く設定することで、信号処理部（乗算手段４０９等）の動作クロックを第１の実施形態よりも遅くできる。

一般に音響波が被検体内の部位から到着するのに要する時間は１００μ秒程度である。この１００μ秒の期間の受信信号を１μ秒ごとに取得すると１００個の信号がラインメモリ４０６に蓄積される。これを１０００倍の時間間隔で読み出したとしても、ラインメモリ４０６から全信号を出力するのに要する時間は１００ｍ秒程度に過ぎない。この程度の時間であれば、人間には実時間動作に感じられる。一方、このように信号の出力間隔を広げると、選択手段４０７による走査、係数の読み出し、乗算、ボクセルメモリへの書き込みなどの処理を余裕を持って行える様になる。これは回路設計上、有利である。

＜第３の実施形態＞
図６は、第３の実施形態のシステム構成を示している。本実施形態は、上述したアルゴリズム３を超音波エコーの画像データ形成に適用した例である。本システムは、送信回路６０４、Ｍ個の音響変換器６０３を備える音響波プローブ、Ｍ個の乗算手段６０８、係数メモリ６０９、加算手段６１０、選択手段６１１、Ｎ個のボクセルメモリ６１２を備える。なお、破線で図示したように、必要に応じて音響変換器６０３と乗算手段６０８の間にラインメモリ６０７を挿入し、第２の実施形態と同様の処理を行ってもよい。

被検体６０１内の部位ｊ６０２に音響変換器６０３から音響波（超音波）を送信する。圧電素子や容量結合を用いた音響変換器６０３は、音響波の受信だけでなく、送信にも使用可能である。しかし、図６の構成に限られず、音響変換器６０３とは別に超音波送信手段を設けても良い。個々の音響変換器６０３（または超音波送信手段）には、送信回路６０４より適切な遅延がかかった送信信号が送られ、その結果、部位ｊに収束する超音波が送信される。部位ｊが周囲の組織と異なる音響特性を持つと、その差に応じた強度の超音波が反射され、音響変換器６０３により受信される。

時刻t^*における音響変換器６０３またはラインメモリ６０７の出力p_d(x_i,y_i,t^*)は、各音響変換器６０３に応じた乗算手段６０８に入力される。各乗算手段６０８には、係数メモリ６０９から部位ｊに対応した係数Ａが読み出される。乗算手段６０８は音響変換器６０３の出力p_dと係数Ａの積Ａ(x_j,y_j,z_j:x_i,y_i,t^*)・p_d(x_i,y_i,t^*)を計算する。Ｍ個の乗
算手段６０８の出力が加算手段６１０で合計され、∬Ａ(x_j,y_j,z_j:x_i,y_i,t^*)・p_d(x_i,y_i,t^*)dx_idy_iが求められる。これが選択手段６１１を介して部位ｊに対応したボクセルメモ
リ６１２に加算される。この処理を所定の期間に亘って繰り返すと、(103)式で示される
部位ｊに対する反射強度p₀(x_j,y_j,z_j)が得られる。次に他の部位ｊ’に収束する様に超音波を送信し、上記の処理を行うと部位ｊ’の反射強度が得られる。これを全部位について行うと、被検体内部の超音波反射強度の分布が画像データ化できる。

本実施形態によれば、Ｍ個の乗算を並列処理できるので、演算ステップ数はＮ×Ｔ（＝Ｌ^４）回となる。従来のフーリエドメイン法で同じ演算を行うにはＬ^６ステップが必要となることと比較すると、本実施形態の構成により大幅に計算時間が短縮されることがわかる。

本実施形態のシステム構成はＰＡＴに適用する事も可能である。その場合には、上述したアルゴリズム２を適用し、時刻t^*ごとに部位ｊに関する走査を行う事も好ましい。この場合も従来のフーリエドメイン法に比べてはるかに少ない演算ステップ数（Ｌ^４）で画像データの形成が可能である。

さらに本実施形態は、同じシステム構成でＰＡＴと超音波エコーの画像データを形成することが可能であるという利点もある。例えば、ＰＡＴの画像と超音波エコーの画像を重ね合わせて表示すれば、被検体内部の様子に関してより詳細な情報を得る事ができる。

なお、上記実施形態では、乗算手段、加算手段、選択手段、係数メモリ、ボクセルメモリなどを専用回路で構成する例を示したが、同様の機能及び処理をコンピュータのソフトウェア（プログラム）で実現することもできる。即ち、係数メモリ及びボクセルメモリの機能はコンピュータの記憶装置（メモリ）が担い、音響変換器から取り込んだデジタル信号の乗算、加算、選択等の処理をプログラムが担うのである。具体的には、上述した実施形態の機能を実現するソフトウェア（プログラム）を、ネットワーク又は各種記憶媒体を介してシステム或いは装置に供給し、そのシステム或いは装置のコンピュータ（またはＣＰＵやＭＰＵ）がプログラムを読み込み実行することで、上述した機能及び処理が実現される。このとき、並列処理可能なＣＰＵを用いたり、複数のＣＰＵで並列処理を行うことで、計算時間の大幅な短縮を実現できる。

４０１，６０１：被検体、４０２，６０２：部位、４０５，６０３：音響変換器、４０８，６０９：係数メモリ、４０９，６０８：乗算手段、４１０，６１２：ボクセルメモリ

Claims (10)

1. 被検体内部より放射される音響波に基づき被検体内部の情報を画像データ化する画像形成装置であって、
   音響波を受信し、信号に変換する複数の音響変換器と、
   前記複数の音響変換器から出力された受信信号を、フーリエドメイン法で処理することによって、被検体内部の複数の部位それぞれから放射された音響波の強度を算出する信号処理部と、を備え、
   前記信号処理部は、
   前記複数の音響変換器の受信信号からある部位が放射した音響波の強度を求める演算式の中で、音響変換器の位置と、部位の位置と、音響波が前記部位から到達するのに要する時間と、波数の次元を有するフーリエ変換の変数と、で定まる項の値を予め計算し、係数として記憶している係数メモリと、
   前記音響変換器の受信信号に対して、対応する係数を前記係数メモリから取得し、乗算する乗算手段と、
   前記乗算手段の乗算結果を部位ごとに累積するボクセルメモリと、
   を備えることを特徴とする画像形成装置。
2. 前記乗算手段は、前記複数の部位のそれぞれに対応した複数の乗算手段からなり、
   前記複数の乗算手段は、１つの音響変換器から得られた受信信号に対して、前記複数の部位のそれぞれに対応した係数を乗算する処理を並列に実行する
   ことを特徴とする請求項１に記載の画像形成装置。
3. 前記乗算手段は、前記複数の音響変換器のそれぞれに対応した複数の乗算手段からなり、
   前記複数の乗算手段は、複数の音響変換器のそれぞれから得られた受信信号に対して、１つの部位に対応した係数を乗算する処理を並列に実行する
   ことを特徴とする請求項１に記載の画像形成装置。
4. 前記音響変換器と前記信号処理部の間に、各音響変換器の受信信号を順次蓄積し、蓄積された受信信号を受信された順番に前記信号処理部へ出力するバッファメモリをさらに備
   える
   ことを特徴とする請求項１乃至３のうちいずれか１項に記載の画像形成装置。
5. 前記バッファメモリは、第１の時間間隔で前記音響変換器の受信信号を蓄積し、蓄積された受信信号を前記第１の時間間隔よりも長い第２の時間間隔で前記信号処理部へ出力する
   ことを特徴とする請求項４に記載の画像形成装置。
6. 被検体に電磁波パルスを照射するための電磁波パルス照射手段をさらに備え、
   前記音響変換器は、照射された電磁波パルスのエネルギーを吸収した部位から光音響効果によって放射される光音響波を受信する
   ことを特徴とする請求項１乃至５のうちいずれか１項に記載の画像形成装置。
7. 前記音響変換器は、被検体に超音波を送信して、送信された超音波の反射波を受信する
   ことを特徴とする請求項１乃至５のうちいずれか１項に記載の画像形成装置。
8. 音響変換器ｉの位置を(xi,yi)、被検体内の部位ｊの位置を(xj,yj,zj)、音響波が前記部位から到達するのに要する時間をt*としたときに、前記係数Ａは下記式で与えられる
   ことを特徴とする請求項１乃至７のうちいずれか１項に記載の画像形成装置。
   

   

   
   但し、u,v,wは、フーリエ変換の変数であり、sgn(w)は、w>0のとき1、w=0のとき0、w<0のとき-1をとる関数である。
9. 被検体内部の複数の部位より放射された音響波を複数の音響変換器により所定の期間に亘って受信して受信信号に変換し、該受信信号をフーリエドメイン法で処理することによって各部位が放射した音響波の強度を画像データ化する画像形成方法であって、
   複数の音響変換器の受信信号からある部位が放射した音響波の強度を求める演算式の中で、音響変換器の位置と、部位の位置と、音響波が前記部位から到達するのに要する時間と、波数の次元を有するフーリエ変換の変数と、で定まる項の値を予め計算し、係数として係数メモリに格納するステップと、
   前記音響変換器の受信信号に対して、対応する係数を前記係数メモリから取得し、乗算するステップと、
   前記受信信号と前記係数の乗算結果を部位ごとにボクセルメモリに累積するステップと、
   を含むことを特徴とする画像形成方法。
10. 請求項９に記載の画像形成方法の各ステップをコンピュータに実行させるためのプログラム。

Images