(57)【要約】【課題】  本発明は、被検体とプローブとの相対位置の変動量を精度よく推定することのできる装置を提供することを目的とする。【解決手段】  本発明の一実施形態に係る光音響装置は、第１の波長の光及び第２の波長の光のそれぞれの光照射により発生する光音響波を受信することにより、信号群を出力する受信手段と、信号群に基づいて、第１の波長に対応する第１の画像データを生成し、信号群に基づいて、第２の波長に対応する第２の画像データを生成し、第１の画像データと第２の画像データとの画像特性の差が小さくなるように、第１の画像データ及び第２の画像データの少なくとも一方を正規化処理し、正規化処理された第１及び第２の画像データに基づいて、第１及び第２の波長の光の照射タイミングにおける位置ずれ情報を取得する処理手段と、を有する。【選択図】  図１
2018050776.tif 000002
【特許請求の範囲】
【請求項１】
第１の波長の光、及び、前記第１の波長とは異なる第２の波長の光のそれぞれを被検体に照射する光照射手段と、
前記第１の波長の光及び前記第２の波長の光のそれぞれの光照射により発生する光音響波を受信することにより、信号群を出力する受信手段と、
前記信号群に基づいて、前記第１の波長に対応する第１の画像データを生成し、
前記信号群に基づいて、前記第２の波長に対応する第２の画像データを生成し、
前記第１の画像データと前記第２の画像データとの画像特性の差が小さくなるように、前記第１の画像データ及び前記第２の画像データの少なくとも一方を正規化処理し、
正規化処理された前記第１及び第２の画像データに基づいて、前記第１及び第２の波長の光の照射タイミングにおける位置ずれ情報を取得する処理手段と、
を有する
ことを特徴とする光音響装置。
【請求項２】
前記処理手段は、
前記位置ずれ情報及び前記信号群に基づいて、合成画像データを生成する
ことを特徴とする請求項１に記載の光音響装置。
【請求項３】
前記処理手段は、
前記位置ずれ情報に基づいて、前記第１の画像データを位置ずれ補正し、
前記位置ずれ情報に基づいて、前記第２の画像データを位置ずれ補正し、
位置ずれ補正された前記第１及び第２の画像データに基づいて、合成画像データを生成する
ことを特徴とする請求項１に記載の光音響装置。
【請求項４】
前記処理手段は、
前記第１の波長の光の照射タイミングにおける前記受信手段の第１の位置情報を取得し、
前記第１の位置情報及び前記信号群を用いて、前記第１の画像データを生成する
ことを特徴とする請求項１から３のいずれか１項に記載の光音響装置。
【請求項５】
前記処理手段は、
前記第１の波長の光の照射タイミングにおける前記受信手段の第１の位置情報を取得し、
前記位置ずれ情報に基づいて、前記第１の位置情報を位置ずれ補正し、
前記第２の波長の光の照射タイミングにおける前記受信手段の第２の位置情報を取得し、
前記位置ずれ情報に基づいて、前記第２の位置情報を位置ずれ補正し、
位置ずれ補正された前記第１及び第２の位置情報及び前記信号群を用いて、合成画像データを取得する
ことを特徴とする請求項１に記載の光音響装置。
【請求項６】
前記処理手段は、
酸素飽和度の空間分布を示すデータを、前記合成画像データとして生成する
ことを特徴とする請求項２から５のいずれか１項に記載の光音響装置。
【請求項７】
前記処理手段は、
前記第１の画像データと前記第２の画像データとの画像強度の差が小さくなるように、前記第１の画像データ及び前記第２の画像データの少なくとも一方の画像強度を変更する
ことを特徴とする請求項１から６のいずれか１項に記載の光音響装置。
【請求項８】
前記処理手段は、
前記第１の画像データ及び前記第２の画像データの画像強度の最大値を決定し、前記最大値が所定値となるように、前記第１の画像データ及び前記第２の画像データの少なくとも一方の画像強度を変更する
ことを特徴とする請求項７に記載の光音響装置。
【請求項９】
前記処理手段は、
前記第１の画像データ及び前記第２の画像データの画像強度の数値範囲を所定の数値範囲とするように、前記第１の画像データ及び前記第２の画像データの少なくとも一方の画像強度を変更する
ことを特徴とする請求項７に記載の光音響装置。
【請求項１０】
前記処理手段は、
前記第１の画像データと前記第２の画像データとの解像度の差が小さくなるように、前記第１の画像データ及び前記第２の画像データの少なくとも一方の解像度を変更する
ことを特徴とする請求項１から９のいずれか１項に記載の光音響装置。
【請求項１１】
第１の波長の光、及び、前記第１の波長とは異なる第２の波長の光のそれぞれを複数回被検体に照射する光照射手段と、
前記第１の波長の光及び前記第２の波長の光のそれぞれの複数回の光照射により発生する光音響波を受信することにより、信号群を出力する受信手段と、
前記信号群に基づいて得られる前記第１の波長に対応する第１の画像データと、前記信号群に基づいて得られる前記第２の波長に対応する第２の画像データとの画像特性の差が小さくなるように、前記信号群を正規化処理し、
正規化処理された前記信号群に基づいて、前記第１の画像データ及び前記第２の画像データを生成し、
前記第１の画像データ及び前記第２の画像データに基づいて、前記第１及び第２の波長の光の照射タイミングにおける位置ずれ情報を取得する処理手段と、
を有する
ことを特徴とする光音響装置。
【請求項１２】
第１の波長の光が被検体に照射されることにより発生する音響波に基づいて生成された、前記第１の波長に対応する第１の画像データを取得し、
前記第１の波長とは異なる第２の波長の光が前記被検体に照射されることにより発生する音響波に基づいて生成された、前記第２の波長に対応する第２の複数の画像データを取得し、
前記第１の画像データと前記第２の画像データとの画像特性の差が小さくなるように、前記第１の画像データ及び前記第２の画像データの少なくとも一方を正規化処理し、
正規化処理された前記第１及び第２の画像データに基づいて、前記第１及び第２の波長の光の照射タイミングにおける位置ずれ情報を取得する
ことを特徴とする情報処理方法。
【請求項１３】
第１の波長の光、及び、前記第１の波長とは異なる第２の波長の光のそれぞれが被検体に照射されることにより発生する音響波に起因する信号群を取得し、
前記信号群に基づいて得られる前記第１の波長に対応する第１の画像データと、前記信号群に基づいて得られる前記第２の波長に対応する第２の画像データとの画像特性の差が小さくなるように、前記信号群を正規化処理し、
正規化処理された前記信号群に基づいて、前記第１の画像データ及び前記第２の画像データを生成し、
前記第１の画像データ及び前記第２の画像データに基づいて、前記第１及び第２の波長の光の照射タイミングにおける位置ずれ情報を取得する処理手段と、
を有する
ことを特徴とする情報処理方法。
【請求項１４】
請求項１２または１３に記載の情報処理方法をコンピュータに実行させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、光音響装置に関する。
【背景技術】
【０００２】
医療分野において生体の生理的情報である機能情報のイメージングの研究が近年行われている。機能情報のイメージング技術の一つとして光音響イメージング（Ｐｈｏｔｏａｃｏｕｓｔｉｃ  Ｉｍａｇｉｎｇ：ＰＡＩ）がある。
【０００３】
光音響イメージングでは、光が被検体に照射される。照射光は被検体を伝播・拡散し、そのエネルギーは被検体内で吸収される。その結果、光音響効果により音響波（以降、光音響波と呼ぶ）が発生する。この光音響波の受信信号が解析処理されることにより、被検体内部の光学特性値の空間分布が画像データとして取得される。
【０００４】
複数回光を被検体に照射し、複数回の光照射により発生する光音響波に基づいて画像データを生成する場合、複数回の光照射の間に被検体とプローブとの相対位置が変動すると、画像データの解像度が低下する。
【０００５】
特許文献１には、複数回の光照射によって得られる複数の画像データ間の位置ずれ量を算出することにより、被検体とプローブとの相対位置の変動量を推定する技術が開示されている。また、特許文献１には、推定された相対位置の変動量に基づいて、複数の画像データを位置合わせする技術が開示されている。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２０１４−１４０７１６号公報
【発明の概要】
【発明が解決しようとする課題】
【０００７】
ところで、光音響イメージングにおいては、互いに異なる複数の波長の光を被検体に照射することにより機能情報を取得することができる。
【０００８】
しかしながら、複数の波長のそれぞれの光の照射により発生する光音響波に基づいて複数の画像データを生成する場合、複数の画像データ間の位置ずれ量の推定精度が低下してしまう場合がある。
【０００９】
そこで、本発明は、複数の波長の光照射により発生する光音響波に基づいて複数の画像データを生成する場合にも、複数の画像データ間の位置ずれ情報を精度よく推定することのできる装置を提供することを目的とする。すなわち、本発明は、複数の波長の光照射により発生する光音響波に基づいて複数の画像データを生成する場合にも、被検体とプローブとの相対位置の変動量を精度よく推定することのできる装置を提供することを目的とする。
【課題を解決するための手段】
【００１０】
発明の一実施形態に係る光音響装置は、第１の波長の光、及び、第１の波長とは異なる第２の波長の光のそれぞれを被検体に照射する光照射手段と、第１の波長の光及び第２の波長の光のそれぞれの光照射により発生する光音響波を受信することにより、信号群を出力する受信手段と、信号群に基づいて、第１の波長に対応する第１の画像データを生成し、信号群に基づいて、第２の波長に対応する第２の画像データを生成し、第１の画像データと第２の画像データとの画像特性の差が小さくなるように、第１の画像データ及び第２の画像データの少なくとも一方を正規化処理し、正規化処理された第１及び第２の画像データに基づいて、第１及び第２の波長の光の照射タイミングにおける位置ずれ情報を取得する処理手段と、を有する。
【発明の効果】
【００１１】
本発明によれば、被検体とプローブとの相対位置の変動量を精度よく推定することができる。
【図面の簡単な説明】
【００１２】
【図１】第１の実施形態に係る光音響装置の構成を示す模式図
【図２】第１の実施形態に係るプローブの模式図
【図３】第１の実施形態に係るコンピュータ及びその周辺機器の構成図
【図４】第１の実施形態に係る光音響装置の作動のフローチャート
【図５】第１の実施形態に係る光音響装置の測定位置を示す模式図
【図６】第１の実施形態に係るパルスボリュームデータを示す模式図
【図７】第１の実施形態に係る画像間の位置ずれ量を示す模式図
【図８】第１の実施形態に係る位置ずれ量の時間的な補間を説明するためのシーケンス図
【図９】第１の実施形態に係る光音響装置の測定位置を示す別の模式図
【図１０】第１の実施形態に係る光音響装置の測定位置を示す別の模式図
【図１１】第１の実施形態に係る位置ずれ補正を説明するための模式図
【図１２】第２の実施形態に係る光音響装置の作動のフローチャート
【図１３】第２の実施形態に係るＧＵＩを示す模式図
【図１４】第３の実施形態に係る光音響装置の作動のフローチャート
【図１５】第４の実施形態に係るＧＵＩを示す模式図
【発明を実施するための形態】
【００１３】
以下、本発明の実施の形態について、図面を用いて説明する。同一の構成要素には原則として同一の符号を付して、説明を省略する。
【００１４】
［第１の実施形態］
本明細書においては、撮像対象とプローブとの相対位置が変動する事象を「位置ずれ」と呼ぶ。また、撮像対象とプローブとの相対位置の変動量を「位置ずれ量」と呼ぶ。また、位置ずれに起因する、複数の画像データに現れる撮像対象の位置の変化量を、「複数の画像データ間の位置ずれ量」と呼ぶ。並進量、回転量、または変形量などの位置ずれを表現するパラメータを「位置ずれ情報」と呼ぶ。
【００１５】
位置ずれは、複数回の光照射の間に被検体が動くことやプローブが動くことにより発生する。例えば、ユーザがハンドヘルド型のプローブを把持して走査することや走査手段がプローブを機械的に走査することなどによりプローブは動く。これらの位置ずれに起因して、複数回の光照射により発生する光音響波に基づいて画像データを生成するときに解像度の低下が生じてしまう。
【００１６】
ところで、光照射により発生する光音響波の強度は、撮像対象の吸収係数に比例して決定される。また、吸収係数は波長依存性を有しているため、波長毎に光フルエンスが同じであっても、波長によって発生する光音響波の強度は変化する。
【００１７】
一方、特許文献１に記載されたように複数の画像データ間の位置ずれを推定する場合、複数の画像データ間で同じ撮像対象の画像強度が異なると、位置ずれの推定精度が低下してしまう場合がある。そのため、互いに異なる複数の波長のそれぞれの光照射によって得られた複数の画像データを用いて位置ずれ量を推定しようとすると、波長間で画像強度や解像度が異なる等の理由から、位置ずれの推定精度が低下してしまう場合がある。
【００１８】
そこで、本発明の一実施形態に係る光音響装置は、互いに異なる複数の波長の一部の波長に対応する複数の画像データ間の位置ずれ情報を取得し、当該位置ずれ情報をその他の波長に対応する位置ずれ情報の取得に用いる。これにより、特定の波長に対応する位置ずれ情報を基準に複数の波長に対応する位置ずれ情報を推定するため、波長間の画像強度の違いによる位置ずれ情報の取得精度に与える影響を抑制することができる。また、特定の波長に対応する位置ずれ情報を取得し、その他の波長に対応する位置ずれ情報を補間により取得することにより、その他の波長に対応する位置ずれ情報の取得に要する処理量を低減することができる。
【００１９】
すなわち、本発明の一実施形態に係る光音響装置においては、互いに異なる第１の波長及び第２の波長のそれぞれの光が複数回被検体に照射される。そして、第１の波長に対応する第１の複数の画像データが生成され、第１の複数の画像データ間の位置ずれ情報が取得される。ここで取得された位置ずれ情報が、第１の波長の光の照射タイミングにおける被検体とプローブとの相対位置の変動量（位置ずれ量）に相当する。そして、第１の位置ずれ情報に基づいて、第２の波長の光照射タイミングにおける被検体とプローブとの変動量（位置ずれ量）が推定される。
【００２０】
また、本発明の一実施形態に係る光音響装置は、上記のように推定された複数波長の光照射に対応する位置ずれ情報と、複数波長の光照射により発生した光音響波に起因する受信データ（信号群または画像データ）とに基づいて、機能情報が取得される。機能情報としては、オキシヘモグロビン濃度、デオキシヘモグロビン濃度、総ヘモグロビン量、または酸素飽和度等の物質の濃度を示す情報が挙げられる。総ヘモグロビン量とは、オキシヘモグロビンおよびデオキシヘモグロビンの和の総量のことである。酸素飽和度とは、全ヘモグロビンに対するオキシヘモグロビンの割合のことである。
【００２１】
なお、機能情報は、空間分布を表す画像データに限らず、数値や文字を表す情報であってもよい。例えば、機能情報は、被検体を構成する物質の平均的な濃度、物質濃度の空間分布の特定位置の値、物質濃度の空間分布の画素値の統計値（平均値や中央値など）などの情報を含む概念である。例えば、機能情報を示す画像として、被検体を構成する物質の平均的な濃度の数値が、表示部１６０に表示されてもよい。
【００２２】
本発明の一実施形態においては、オキシヘモグロビンのモル吸収係数とデオキシヘモグロビンのモル吸収係数が等しい波長の光を利用して得られる複数の画像データに基づいて、位置ずれ情報を取得してもよい。血管を測定対象とする場合、このような波長の光を用いて血管から発生する光音響波の強度は、酸素飽和度に依存しないため、血管毎の画像強度のばらつきが小さくなる。そのため、このような波長の光を用いることにより、位置ずれ情報の推定精度が高くなる傾向がある。なお、オキシヘモグロビンのモル吸収係数とデオキシヘモグロビンのモル吸収係数が完全に等しい波長（等モル吸収係数の波長）に限らず、略等しい波長であっても位置ずれ情報の推定精度は高くなる傾向がある。例えば、等モル吸収係数の波長の±１０ｎｍ以内の波長を、略等しい波長として採用してもよい。また、例えば、オキシヘモグロビンのモル吸収係数とデオキシヘモグロビンのモル吸収係数が、等モル吸収係数から±１０％以内に相当する波長を略等しい波長として採用してもよい。すなわち、位置ずれ情報の取得に好適な波長の光を利用して得られた位置ずれ情報に基づいて、その他の波長に対応する位置ずれ情報を取得することが好ましい。
【００２３】
以下、第１の実施形態の光音響装置の構成の構成及び処理について説明する。
【００２４】
本実施形態において、光音響装置を用いた例を説明する。図１を用いて本実施形態に係る光音響装置の構成を説明する。図１は、光音響装置全体の概略ブロック図である。本実施形態に係る光音響装置は、駆動部１３０、信号収集部１４０、コンピュータ１５０、表示部１６０、入力部１７０、及びプローブ１８０を有する。
【００２５】
図２は、本実施形態に係るプローブ１８０の模式図である。プローブ１８０は、光照射部１１０、及び受信部１２０を有する。測定対象は、被検体１００である。
【００２６】
駆動部１３０は、光照射部１１０と受信部１２０を駆動し、機械的な走査を行う。光照射部１１０がパルス光１１３を被検体１００に照射し、被検体１００内で音響波が発生する。光に起因して光音響効果により発生する音響波を光音響波とも呼ぶ。受信部１２０は、光音響波を受信することによりアナログ信号としての電気信号（光音響信号）を出力する。
【００２７】
信号収集部１４０は、受信部１２０から出力されたアナログ信号をデジタル信号に変換し、コンピュータ１５０に出力する。コンピュータ１５０は、信号収集部１４０から出力されたデジタル信号を、超音波または光音響波に由来する信号データとして記憶する。
【００２８】
コンピュータ１５０は、記憶されたデジタル信号に対して信号処理を行うことにより、光音響画像を表す画像データを生成する。また、コンピュータ１５０は、得られた画像データに対して画像処理を施した後に、画像データを表示部１６０に出力する。表示部１６０は、光音響画像を表示する。操作者としての医師や技師等は、表示部１６０に表示された光音響画像を確認することにより、診断を行うことができる。表示画像は、操作者やコンピュータ１５０からの保存指示に基づいて、コンピュータ１５０内のメモリや、モダリティとネットワークで接続されたデータ管理システムなどに保存される。
【００２９】
また、コンピュータ１５０は、光音響装置に含まれる構成の駆動制御も行う。また、表示部１６０は、コンピュータ１５０で生成された画像の他にＧＵＩなどを表示してもよい。入力部１７０は、操作者が情報を入力できるように構成されている。操作者は、入力部１７０を用いて測定開始や終了、作成画像の保存指示などの操作を行うことができる。
【００３０】
本実施形態に係る光音響装置により得られる光音響画像は、光照射により発生した光音響波に由来するあらゆる画像を含む概念である。光音響画像は、光音響波の発生音圧（初期音圧）、光吸収エネルギー密度、及び光吸収係数、などの少なくとも１つの空間分布を表す画像データである。
【００３１】
以下、本実施形態に係る光音響装置の各構成の詳細を説明する。
【００３２】
（光照射部１１０）
光照射部１１０は、パルス光１１３を発する光源と、光源から射出されたパルス光１１３を被検体１００へ導く光学系とを含む。なお、パルス光は、いわゆる矩形波、三角波などの光を含む。
【００３３】
光源が発する光のパルス幅としては、１ｎｓ以上、１００ｎｓ以下のパルス幅であってもよい。また、光の波長として４００ｎｍから１６００ｎｍ程度の範囲の波長であってもよい。血管を高解像度でイメージングする場合は、血管での吸収が大きい波長（４００ｎｍ以上、７００ｎｍ以下）を用いてもよい。生体の深部をイメージングする場合には、生体の背景組織（水や脂肪など）において典型的に吸収が少ない波長（７００ｎｍ以上、１１００ｎｍ以下）の光を用いてもよい。
【００３４】
光源としては、レーザーや発光ダイオードを用いることができる。また、複数波長の光を用いて測定する際には、波長の変換が可能な光源であってもよい。なお、複数波長を被検体に照射する場合、互いに異なる波長の光を発生する複数台の光源を用意し、それぞれの光源から交互に照射することも可能である。複数台の光源を用いた場合もそれらをまとめて光源として表現する。レーザーとしては、固体レーザー、ガスレーザー、色素レーザー、半導体レーザーなど様々なレーザーを使用することができる。例えば、Ｎｄ：ＹＡＧレーザーやアレキサンドライトレーザーなどのパルスレーザーを光源として用いてもよい。また、Ｎｄ：ＹＡＧレーザー光を励起光とするＴｉ：ｓａレーザーやＯＰＯ（Ｏｐｔｉｃａｌ  Ｐａｒａｍｅｔｒｉｃ  Ｏｓｃｉｌｌａｔｏｒｓ）レーザーを光源として用いてもよい。また、光源としてマイクロウェーブ源を用いてもよい。
【００３５】
光学系には、レンズ、ミラー、光ファイバ等の光学素子を用いることができる。乳房等を被検体１００とする場合、パルス光のビーム径を広げて照射するために、光学系の光出射部は光を拡散させる拡散板等で構成されていてもよい。一方、光音響顕微鏡においては、解像度を上げるために、光学系の光出射部はレンズ等で構成し、ビームをフォーカスして照射してもよい。
【００３６】
なお、光照射部１１０が光学系を備えずに、光源から直接被検体１００にパルス光１１３を照射してもよい。
【００３７】
（受信部１２０）
受信部１２０は、音響波を受信することにより電気信号を出力するトランスデューサ１２１と、トランスデューサ１２１を支持する支持体１２２とを含む。また、トランスデューサ１２１は、音響波を送信することもできる。便宜上、図２においてはトランスデューサ１２１を１つしか図示していないが、受信部１２０は複数のトランスデューサを含んでいてもよい。
【００３８】
トランスデューサ１２１を構成する部材としては、ＰＺＴ（チタン酸ジルコン酸鉛）に代表される圧電セラミック材料や、ＰＶＤＦ（ポリフッ化ビニリデン）に代表される高分子圧電膜材料などを用いることができる。また、圧電素子以外の素子を用いてもよい。例えば、静電容量型トランスデューサ（ＣＭＵＴ：Ｃａｐａｃｉｔｉｖｅ  Ｍｉｃｒｏ−ｍａｃｈｉｎｅｄ  Ｕｌｔｒａｓｏｎｉｃ  Ｔｒａｎｓｄｕｃｅｒｓ）、ファブリペロー干渉計を用いたトランスデューサなどを用いることができる。なお、音響波を受信することにより電気信号を出力できる限り、いかなるトランスデューサを採用してもよい。また、トランスデューサにより得られる信号は時間分解信号である。つまり、トランスデューサにより得られる信号の振幅は、各時刻にトランスデューサで受信される音圧に基づく値（例えば、音圧に比例した値）を表したものである。
【００３９】
光音響波を構成する周波数成分は、典型的には１００ＫＨｚから１００ＭＨｚであり、トランスデューサ１２１として、これらの周波数を検出することのできるものを採用することができる。
【００４０】
支持体１２２は、機械的強度が高い金属材料などから構成されていてもよい。照射光を被検体に多く入射させるために、支持体１２２の被検体１００側の表面に鏡面もしくは光散乱させる加工が行われていてもよい。本実施形態において支持体１２２は半球殻形状であり、半球殻上に複数のトランスデューサ１２１を支持できるように構成されている。この場合、支持体１２２に配置されたトランスデューサ１２１の指向軸は半球の曲率中心付近に集まる。そして、複数のトランスデューサ１２１から出力された信号を用いて画像化したときに曲率中心付近の画質が高くなる。なお、支持体１２２はトランスデューサ１２１を支持できる限り、いかなる構成であってもよい。支持体１２２は、１Ｄアレイ、１．５Ｄアレイ、１．７５Ｄアレイ、２Ｄアレイと呼ばれるような平面又は曲面内に、複数のトランスデューサを並べて配置してもよい。
【００４１】
また、支持体１２２は音響マッチング材を貯留する容器として機能してもよい。すなわち、支持体１２２をトランスデューサ１２１と被検体１００との間に音響マッチング材を配置するための容器としてもよい。
【００４２】
また、受信部１２０が、トランスデューサ１２１から出力される時系列のアナログ信号を増幅する増幅器を備えてもよい。また、受信部１２０が、トランスデューサ１２１から出力される時系列のアナログ信号を時系列のデジタル信号に変換するＡ／Ｄ変換器を備えてもよい。すなわち、受信部１２０が後述する信号収集部１４０を備えてもよい。
【００４３】
なお、音響波を様々な角度で検出できるようにするために、理想的には被検体１００を全周囲から囲むようにトランスデューサ１２１を配置してもよい。ただし、被検体１００が大きく全周囲を囲むようにトランスデューサを配置できない場合は、図２に示したように半球状の支持体１２２上にトランスデューサを配置して全周囲を囲む状態に近づけてもよい。図２においては、トランスデューサ１２１を１つしか図示していないが、複数のトランスデューサを半球上の支持体１２２上に配置してもよい。
【００４４】
なお、トランスデューサの配置や数及び支持体の形状は被検体に応じて最適化すればよく、本発明に関してはあらゆる受信部１２０を採用することができる。
【００４５】
受信部１２０と被検体１００との間の空間は、光音響波が伝播することができる媒質で満たす。この媒質には、音響波が伝搬でき、被検体１００やトランスデューサ１２１との界面において音響特性が整合し、できるだけ光音響波の透過率が高い材料を採用する。例えば、この媒質には、水、超音波ジェルなどを採用することができる。
【００４６】
（駆動部１３０）
駆動部１３０は、被検体１００と受信部１２０との相対位置を変更する部分である。本実施形態では、駆動部１３０は、支持体１２２をＸＹ方向に移動させる装置であり、ステッピングモータを搭載した電動のＸＹステージある。駆動部１３０は、駆動力を発生させるステッピングモータなどのモータと、駆動力を伝達させる駆動機構と、受信部１２０の位置情報を検出する位置センサとを含む。駆動機構としては、リードスクリュー機構、リンク機構、ギア機構、油圧機構、などを用いることができる。また、位置センサとしては、エンコーダー、可変抵抗器、などを用いたポテンショメータなどを用いることができる。
【００４７】
なお、駆動部１３０は被検体１００と受信部１２０との相対位置をＸＹ方向（二次元）に変更させるものに限らず、一次元または三次元に変更させてもよい。
【００４８】
なお、駆動部１３０は、被検体１００と受信部１２０との相対的な位置を変更できれば、受信部１２０を固定し、被検体１００を移動させてもよい。被検体１００を移動させる場合は、被検体１００を指示する被検体支持部（不図示）を動かすことで被検体１００を移動させる構成などが考えられる。また、被検体１００と受信部１２０の両方を移動させてもよい。
【００４９】
また駆動部１３０は、相対位置を連続的に移動させてもよいし、ステップアンドリピートによって移動させてもよい。駆動部１３０は、電動ステージであってもよいし、手動ステージであってもよい。
【００５０】
また、本実施形態では、駆動部１３０は光照射部１１０と受信部１２０を同時に駆動して走査を行っているが、光照射部１１０だけを駆動したり、受信部１２０だけを駆動したりしてもよい。
【００５１】
（信号収集部１４０）
信号収集部１４０は、トランスデューサ１２１から出力されたアナログ信号である電気信号を増幅するアンプと、アンプから出力されたアナログ信号をデジタル信号に変換するＡ／Ｄ変換器とを含む。信号収集部１４０は、ＦＰＧＡ（Ｆｉｅｌｄ  Ｐｒｏｇｒａｍｍａｂｌｅ  Ｇａｔｅ  Ａｒｒａｙ）チップなどで構成されてもよい。信号収集部１４０から出力されるデジタル信号は、コンピュータ１５０内の記憶部１５２に記憶される。信号収集部１４０は、Ｄａｔａ  Ａｃｑｕｉｓｉｔｉｏｎ  Ｓｙｓｔｅｍ（ＤＡＳ）とも呼ばれる。本明細書において電気信号は、アナログ信号もデジタル信号も含む概念である。なお、信号収集部１４０は、光照射部１１０の光射出部に取り付けられた光検出センサと接続されており、パルス光１１３が光照射部１１０から射出されたことをトリガーに、同期して処理を開始してもよい。また、信号収集部１４０は、フリーズボタンなどを用いてなされる指示をトリガーに同期して、当該処理を開始してもよい。
【００５２】
（コンピュータ１５０）
コンピュータ１５０は、演算部１５１、記憶部１５２、制御部１５３を含む。各構成の機能については処理フローの説明の際に説明する。
【００５３】
演算部１５１としての演算機能を担うユニットは、ＣＰＵやＧＰＵ（Ｇｒａｐｈｉｃｓ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）等のプロセッサ、ＦＰＧＡ（Ｆｉｅｌｄ  Ｐｒｏｇｒａｍｍａｂｌｅ  Ｇａｔｅ  Ａｒｒａｙ）チップ等の演算回路で構成されることができる。これらのユニットは、単一のプロセッサや演算回路から構成されるだけでなく、複数のプロセッサや演算回路から構成されていてもよい。演算部１５１は、入力部１７０から、被検体音速や保持カップの構成などの各種パラメータを受けて、受信信号を処理してもよい。
【００５４】
記憶部１５２は、ＲＯＭ（Ｒｅａｄ  ｏｎｌｙ  ｍｅｍｏｒｙ）、磁気ディスクやフラッシュメモリなどの非一時記憶媒体で構成することができる。また、記憶部１５２は、ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）などの揮発性の媒体であってもよい。なお、プログラムが格納される記憶媒体は、非一時記憶媒体である。なお、記憶部１５２は、１つの記憶媒体から構成されるだけでなく、複数の記憶媒体から構成されていてもよい。
【００５５】
記憶部１５２は、後述する方法で演算部１５１により生成される光音響画像を示す画像データを保存することができる。
【００５６】
制御部１５３は、ＣＰＵなどの演算素子で構成される。制御部１５３は、光音響装置の各構成の動作を制御する。制御部１５３は、入力部１７０からの測定開始などの各種操作による指示信号を受けて、光音響装置の各構成を制御してもよい。また、制御部１５３は、記憶部１５２に格納されたプログラムコードを読み出し、光音響装置の各構成の作動を制御する。
【００５７】
コンピュータ１５０は専用に設計されたワークステーションであってもよい。また、コンピュータ１５０の各構成は異なるハードウェアによって構成されてもよい。また、コンピュータ１５０の少なくとも一部の構成は単一のハードウェアで構成されてもよい。
【００５８】
図３は、本実施形態に係るコンピュータ１５０の具体的な構成例を示す。本実施形態に係るコンピュータ１５０は、ＣＰＵ１５４、ＧＰＵ１５５、ＲＡＭ１５６、ＲＯＭ１５７、外部記憶装置１５８から構成される。また、コンピュータ１５０には、表示部１６０としての液晶ディスプレイ１６１、入力部１７０としてのマウス１７１、キーボード１７２が接続されている。
【００５９】
また、コンピュータ１５０および複数のトランスデューサ１２１は、共通の筺体に収められた構成で提供されてもよい。ただし、筺体に収められたコンピュータで一部の信号処理を行い、残りの信号処理を筺体の外部に設けられたコンピュータで行ってもよい。この場合、筺体の内部および外部に設けられたコンピュータを総称して、本実施形態に係るコンピュータとすることができる。すなわち、コンピュータを構成するハードウェアが一つの筺体に収められていなくてもよい。
【００６０】
（表示部１６０）
表示部１６０は、液晶ディスプレイや有機ＥＬ（Ｅｌｅｃｔｒｏ  Ｌｕｍｉｎｅｓｃｅｎｃｅ）などのディスプレイである。コンピュータ１５０により得られた被検体情報等に基づく画像や特定位置の数値等を表示する装置である。表示部１６０は、画像や装置を操作するためのＧＵＩを表示してもよい。なお、被検体情報の表示にあたっては、表示部１６０またはコンピュータ１５０において画像処理（輝度値の調整等）を行った上で表示することもできる。
【００６１】
（入力部１７０）
入力部１７０としては、ユーザーが操作可能な、マウスやキーボードなどで構成される操作コンソールを採用することができる。入力できる内容としては、画像再構成の条件の選択であったり、位置ずれ補正の方法の選択であったり、補間方法の選択であったりしてもよい。また、スライダーバーのような形で、合成された画像を確認しながら、位置ずれ量の足し合わせの重みなどを変更できるようにしてもよい。また、表示部１６０をタッチパネルで構成し、表示部１６０を入力部１７０として利用してもよい。
【００６２】
なお、光音響装置の各構成はそれぞれ別の装置として構成されてもよいし、一体となった１つの装置として構成されてもよい。また、光音響装置の少なくとも一部の構成が一体となった１つの装置として構成されてもよい。
【００６３】
（被検体１００）
被検体１００は光音響装置を構成するものではないが、以下に説明する。本実施形態に係る光音響装置は、人や動物の悪性腫瘍や血管疾患などの診断や化学治療の経過観察などを目的として使用できる。よって、被検体１００としては、生体、具体的には人体や動物の乳房や各臓器、血管網、頭部、頸部、腹部、手指および足指を含む四肢などの診断の対象部位が想定される。例えば、人体が測定対象であれば、オキシヘモグロビンあるいはデオキシヘモグロビンやそれらを含む多く含む血管あるいは腫瘍の近傍に形成される新生血管などを光吸収体の対象としてもよい。また、頸動脈壁のプラークなどを光吸収体の対象としてもよい。また、メチレンブルー（ＭＢ）、インドシニアングリーン（ＩＣＧ）などの色素、金微粒子、またはそれらを集積あるいは化学的に修飾した外部から導入した物質を光吸収体としてもよい。
【００６４】
以下、図４に示すフローチャートに沿って、本実施形態に係る情報処理を含む光音響装置の作動を説明する。
【００６５】
（Ｓ１００：波長λ１の光及び波長λ２の光をそれぞれ複数回照射し、光音響波を受信する工程）
光照射部１１０は、互いに異なる波長λ１及び波長λ２の光のそれぞれを複数回照射し、受信部１２０は、それらの光照射により発生する光音響波を受信する。
【００６６】
制御部１５３は、走査情報と光照射を示す情報（制御信号）をプローブ１８０に送信する。駆動部１３０が受信部１２０を移動させながら、光照射部１１０は複数回にわたって複数の波長のパルス光を被検体１００に照射する。すなわち、複数回の光照射が行われる期間に、駆動部１３０は受信部１２０を移動させる。その結果、駆動部１３０は、受信部１２０が各光照射時に互いに異なる位置に位置するように受信部１２０を移動させることができる。トランスデューサ１２１は、光照射部１１０が複数回パルス光を照射することにより発生した光音響波を受信することにより、光照射回数分だけ信号を出力する。以下、複数波長による複数回の光照射回数分出力された信号を総称して、複数波長に対応する信号群と呼ぶ。
【００６７】
以下、Ｎ回の光照射を行う場合を説明する。なお、波長λ１で、ｉ回目の光照射により得られた信号を
【００６８】
【数１】
2018050776.tif 000003
【００６９】
と表記する。添字ｉが付いたアイテムは、ｉ回目の光照射に対応するアイテムであることを表す。ｉは正の整数であり、パルスインデックスとも呼ぶ。
【００７０】
例えば、図５に示すように、プローブ１８０を移動させながら光照射及び光音響波の受信を行うことができる。符号５０１で示された丸は、波長λ１の光を照射したときのプローブ１８０の位置（測定位置）を表す。符号５０２で示された三角は、波長λ２の光を照射したときのプローブ１８０の位置（測定位置）を表す。符号５０３で示された実線は、プローブ１８０の軌跡を表す。図５に示されたように、本実施形態では、光照射部１１０が複数の波長の光を交互に照射しながら、駆動部１３０がプローブ１８０を動かしてもよい。図５に示す例では最外周の軌跡から内側に向かって渦巻状にプローブ１８０を走査している。図５では便宜上ＸＹ平面で描写しているが、平面にかぎらず立体的にプローブを走査してもよい。
【００７１】
信号収集部１４０は、トランスデューサ１２１から出力されたアナログ信号群である複数波長に対応する信号群に対してＡＤ変換処理等を行い、処理後の光音響信号をコンピュータ１５０に送信する。デジタル信号群としての光音響信号は、記憶部１５２に記憶される。
【００７２】
（Ｓ２００：波長λ１に対応する複数の画像データを取得する工程）
演算部１５１は、Ｓ１００で取得された信号群に基づいて、第１の波長の光照射に対応する複数の画像データを取得する。演算部１５１は、各光照射の信号のそれぞれから画像データを生成し、その中から波長λ１の光照射により得られた画像データを取得してもよい。また、演算部１５１は、信号群の中から波長λ１の光照射に対応する信号を選択的に用いて、波長λ１の光照射により得られた画像データを取得してもよい。
【００７３】
演算部１５１は、光音響信号に対してＵｎｉｖｅｒｓａｌ  Ｂａｃｋ−Ｐｒｏｊｅｃｔｉｏｎ（ＵＢＰ）等の再構成処理を行うことにより、光音響画像を生成する。なお、光音響画像を生成したところで、記憶部１５２に保存された光音響信号を削除してもよい。１回のパルス光の照射によって得られる画像データをパルスボリュームデータとも呼ぶ。パルスボリュームデータは、２次元または３次元に配列したボクセル（２次元の場合はピクセルとも呼ぶ）のそれぞれに当該位置における値を収めたボリュームデータの形式で取得される。ボリュームデータは、２次元または３次元ボリューム、２次元または３次元画像、２次元または３次元断層像とも呼べる。
【００７４】
再構成手法については、タイムドメイン再構成手法、フーリエドメイン再構成手法、モデルベース再構成手法（繰り返し再構成手法）などの高知の再構成手法を採用することができる。例えば、ＰＨＹＳＩＣＡＬ  ＲＥＶＩＥＷ  Ｅ７１，０１６７０６（２００５）に記載されたようなＵｎｉｖｅｒｓａｌ  Ｂａｃｋ−Ｐｒｏｊｅｃｔｉｏｎ（ＵＢＰ）と呼ばれるタイムドメイン再構成手法を採用することができる。
【００７５】
演算部１５１は、画像データとして初期音圧分布データ
【００７６】
【数２】
2018050776.tif 000004
【００７７】
を取得してもよい。ａは、波長λａの光照射に対応するアイテムであることを示す波長インデックスである。演算部１５１は、信号群に加え、光照射時のトランスデューサの位置情報に基づいて、初期音圧分布データを取得してもよい。演算部１５１は、記憶部１５２に予め記憶された各光照射時のトランスデューサの位置情報を読み出すことにより位置情報を取得することができる。また、演算部１５１は、光照射をトリガーとして、駆動部１３０に備えられた位置センサからの受信部１２０の位置情報を受け取ることにより、トランスデューサの位置情報を取得してもよい。
【００７８】
図６は、本実施形態に係る波長λ１でのパルスボリュームデータの一部（Ｐλ１，１〜Ｐλ１，１０）を示す。本実施形態におけるパルスボリュームデータは３次元空間中のボリュームデータであるが、紙面の説明の都合からパルスボリュームデータをＸＹ面で表す。本実施形態においては、時間的に隣接する初期音圧分布データの少なくとも一部の領域が重畳するように、再構成領域が設定される。本実施形態では、半球状の支持体１２２の曲率中心を中心とした６０ｍｍ角の立方体領域を、１回の光照射、すなわち１つの電気信号群に基づいて再構成される再構成領域とする。この場合、１回の光照射による再構成領域の大きさ（６０ｍｍ）は、光照射間の受信部１２０の移動量と比べて大きい。そのため、図６に示すように、時間的に一続きの光照射に対応する２つ以上のパルスボリュームデータが重畳することとなる。再構成領域の大きさや形状は、予め設定されていてもよい。また、ユーザが入力部１７０を用いて再構成領域の大きさや形状を指定してもよい。基準位置に対する各パルスボリュームデータの中心位置を、各パルスボリュームデータの位置とする。図６では、一例としてパルスボリュームデータＰλ１，１０の位置ＰｏｓＰλ１，１０を示した。この図の例では光照射ごとに受信部１２０の位置が異なるため、図６に示された本実施形態で得られる各パルスボリュームデータは、基準位置に対して互いに異なる位置に位置する。
【００７９】
なお、本工程において、演算部１５１は、被検体内での光フルエンス分布データΦ［Ｐａ・ｍ３／Ｊ］を取得してもよい。そして、演算部１５１は、初期音圧分布データを、光フルエンス分布データとグリュナイゼン係数分布データとで除算することにより、被検体内の光吸収係数分布データμａ［１／ｍ］を取得してもよい。この場合、光吸収係数分布データをパルスボリュームデータとしてもよい。
【００８０】
例えば、演算部１５１は、Ｐｒｏｃ．ｏｆ  ＳＰＩＥ  Ｖｏｌ．７５６１  ７５６１１７−１に記載されたように、光拡散方程式を解くことにより光フルエンス分布データを取得してもよい。
【００８１】
また、例えば、グリュナイゼン係数は被検体の種類が決定するとほぼ一意に値が決定されることが知られているため、被検体に対応するグリュナイゼン係数分布データΓを予め記憶部１５２に記憶しておくことができる。そして、演算部１５１は、予め記憶部１５２に記憶されたグリュナイゼン係数分布データΓを読み出すことにより、取得してもよい。
【００８２】
なお、ユーザが、把持部を有するプローブ１８０を把持し、プローブ１８０を移動させてもよい。また、複数回の光照射を行っている期間にプローブ１８０を移動させなくてもよい。また、演算部１５１は、１回の光照射により得られる電気信号群に基づいて全画像化領域の画像データを取得し、それを複数回の光照射に繰り返してもよい。
【００８３】
（Ｓ３００：波長λ１の光の照射タイミングにおける位置ずれ情報を取得する工程）
演算部１５１は、Ｓ２００で取得された波長λ１に対応する複数の画像データに基づいて、波長λ１の光の照射タイミングにおける位置ずれ情報を取得する。例えば、以下説明するように、パルスボリュームデータ間の位置ずれ量を算出することや、パルスボリュームデータが合成されたボリュームデータの位置ずれ量を算出することにより、当該位置ずれ情報を取得することができる。
【００８４】
まず、演算部１５１は、波長λ１の光照射間における、対象物と受信部との相対的な位置関係の変動による、各パルスボリュームデータの位置ずれ量を推定する。その際、得られたパルスボリュームデータから、任意のパルスボリュームデータのペアを選択する。ｋ個目のペアをＲ＿ｋと表記する。また、ペアＲ＿ｋを構成する波長λ１のパルスボリュームデータの一方をＰλ１，ｋ１、もう一方をＰλ１，ｋ２と表記する。以下、本実施形態では、Ｋ個のペアが選択された場合について説明する。なお、オーバーラップ領域を有する２つのパルスボリュームデータをペアとすることが好ましい。これにより、共通の特徴を有さないパルスボリュームデータ同士を比較することを避けることができるため、冗長な計算を減らすことができる。さらに、オーバーラップ領域の大きいパルスボリュームデータ同士をペアとすることが好ましい。そこで、例えば、演算部１５１は、パルスボリュームデータ間のオーバーラップ領域の体積の割合が所定の値以上のペアを選択してもよい。また、後述するように、位置ずれ補正に合成ボリュームデータを用いる場合には、パルスボリュームデータの重畳数が多い領域同士が重なるようなペアを選択してもよい。
【００８５】
また、あるパルスボリュームデータに対して、当該パルスボリュームデータのインデックスから、インデックスが所定の範囲に含まれるパルスボリュームデータをペアの対象として選択してもよい。また、インデックスが連続する、すなわち時間的に連続するパルスボリュームデータをペアの対象として選択してもよい。例えば、本実施形態では、演算部１５１は、オーバーラップ領域が５０％以上となるパルスボリュームデータを対象としてペアを選択する。
【００８６】
以下、各パルスボリュームデータの位置ずれ量を推定する方法の例を説明する。
【００８７】
演算部１５１は、式１に示すように、Ｐλ１，ｋ１とＰλ１，ｋ２との間の類似度関数Ｆ＿ｋを取得する。
Ｆ＿ｋ（ｘ，ｙ，ｚ）＝ｆｓｉｍｉｌ（Ｐ＿ｋ，ｘ，ｙ，ｚ）・・・式１
【００８８】
ここで、類似度関数Ｆ＿ｋは、ペアＲ＿ｋを構成する片方のパルスボリュームデータＰλ１，ｋ１に対する、もう一方のパルスボリュームデータＰλ１，ｋ２の相対位置を（ｘ，ｙ，ｚ）だけ並進させた場合の類似度を算出する関数である。ここで関数ｆｓｉｍｉｌは、画像間の類似度が高い場合には関数値として高い値を返すものとする。類似度関数Ｆの取得とは、各関数の引数である並進量（ｘ，ｙ，ｚ）、すなわち画像データ間の相対位置を所定の範囲内で離散的に変化させた場合の関数値の取得を意味する。例えば、ｘ，ｙ，ｚのそれぞれの値を−Ｌから＋Ｌまでの整数値として変化させた場合のそれぞれについてＦが返す（２Ｌ＋１）×（２Ｌ＋１）×（２Ｌ＋１）個の値の集合の取得を意味する。より発展的には、（２Ｌ＋１）×（２Ｌ＋１）×（２Ｌ＋１）個の値の集合をさらにバイリニア法やバイキュービック法などを用いて、より連続関数に近い情報として類似度関数Ｆを導出し、これを取得するものとしてよい。
【００８９】
なお、Ｐλ１，ｋ１に対するＰλ１，ｋ２の相対位置（光照射間の受信部１２０の移動量）だけ並進させた位置を基準として、Ｐλ１，ｋ２の位置を所定の範囲内で離散的に変化させた場合の関数値を取得してもよい。
【００９０】
例えば、類似度を算出する関数としては、ＳＳＤ（Ｓｕｍ  ｏｆ  Ｓｑｕａｒｅｄ  Ｄｉｆｆｅｒｅｎｃｅ）やＳＡＤ（Ｓｕｍ  ｏｆ  Ａｂｓｏｌｕｔｅ  Ｄｉｆｆｅｒｅｎｃｅ）、相互情報量、相互相関など、任意の類似度尺度が適用できる。また、例えば、パルスボリュームデータから特徴的な形態を抽出し、それらの位置との一致度を測ることによって類似度関数を取得してもよい。抽出する特徴としては、血管などの解剖学的な特徴、エッジ検出やコーナー検出などの画像処理分野で一般的に用いられる公知の技術によって抽出された特徴を用いてもよい。また、抽出する特徴としては、コンピュータビジョン等の技術分野で一般的に使われるＳＩＦＴ特徴やＳＵＲＦ特徴などのより高次の局所画像特徴等を用いてもよい。これらの方法によれば、パルスボリュームデータ間の輝度分布の相違やノイズの混入によって、より頑健な類似度関数を取得することができると考えられる。
【００９１】
なお、演算部１５１は、類似度の計算の結果に、重みを乗ずる処理を施すことにより類似度関数を取得してもよい。
【００９２】
また、類似度算出の対象となるパルスボリュームデータ間で正しく類似度を算出することができなかった場合、その結果を以降の処理には使用しなくてもよい。正しく類似度を算出することができない場合としては、いずれに並進させても類似度が十分小さい、または変わらない場合などが考えられる。この処理によれば、同一の特徴が十分に現れているパルスボリュームデータ同士の比較結果（類似度関数）を選択的に、以降の処理に使用することができる。
【００９３】
続いて、演算部１５１は、式２に示すように、類似度関数Ｆ＿ｋの関数値が最大となる並進量Ｍ＿ｋの取得を行う。
Ｍ＿ｋ＝ａｒｇ  ｍａｘ（Ｆ＿ｋ（ｘ，ｙ，ｚ））・・・式２
【００９４】
演算部１５１は、各ペアについて類似度関数Ｆ＿ｋの関数値が最大となる並進量Ｍ＿ｋの取得を行う。
【００９５】
パルスボリュームデータの位置を推定する場合に、ペアＲ＿ｋに対する個別最適値である並進量Ｍ＿ｋをなるべく保つような評価関数を定義する。すなわち、Ｐλ１，ｋ１に対するＰλ１，ｋ２の位置が並進量Ｍ＿ｋから離れるにつれて、値が低下する評価関数を定義する。式３は、この場合の評価関数Ｅ＿ｋの例を表す。
Ｅ＿ｋ＝（Ｍ＿ｋ−（ＰｏｓＰλ１，ｋ１−ＰｏｓＰλ１，ｋ２）２）
＝（Ｍ＿ｋ（ｘ）−（ＰｏｓＰλ１，ｋ１（ｘ）−ＰｏｓＰλ１，ｋ２（ｘ））２）
＋（Ｍ＿ｋ（ｙ）−（ＰｏｓＰλ１，ｋ１（ｙ）−ＰｏｓＰλ１，ｋ２（ｙ））２）
＋（Ｍ＿ｋ（ｚ）−（ＰｏｓＰλ１，ｋ１（ｚ）−ＰｏｓＰλ１，ｋ２（ｚ））２）
・・・式３
【００９６】
ＰｏｓＰλ１，ｋ１は、基準位置に対するＰλ１，ｋ１の位置を表す。ＰｏｓＰλ１，ｋ２は、基準位置に対するＰλ１，ｋ２の位置を表す。なお、評価関数を定義する際に、類似度関数Ｆ＿ｋを、当該類似度関数Ｆ＿ｋにフィットするような二次関数に近似してもよい。また、類似度関数Ｆ＿ｋが、並進量Ｍ＿ｋの周辺において、二次関数に従って低下すると禁じできる場合には、式３はＰλ１，ｋ１とＰλ１，ｋ２の位置関係から類似度関数Ｆ＿ｋの値を並進量Ｍ＿ｋの周辺で近似する関数となる。
【００９７】
続いて、演算部１５１は、式４のように定義されたコスト関数Ｅが最小化した時の、基準位置に対する全てのパルスボリュームデータの位置ＰｏｓＰ′λ１，ｊを取得する。ここで、ｊはあるパルスに対するパルスインデックスである。
【００９８】
【数３】
2018050776.tif 000005
【００９９】
コスト関数が最小化したときの、基準位置に対するパルスボリュームデータの位置は、被検体１００と受信部１２０との相対的な位置関係の変動による位置ずれ後のパルスボリュームデータの位置情報を表す。
【０１００】
例えば、演算部１５１は、式４に示すコスト関数Ｅを最小化する（０に最近接する）解を線形最小二乗法を解くことにより得る。これにより、一意に各パルスボリュームデータの位置ＰｏｓＰ′λ１，ｊを求めることができるため、計算コストが小さい。
【０１０１】
なお、上記で説明した線形最適化によるコスト関数の最適化に限らず、コスト関数の最適化は公知のいかなる方法を用いてもよい。例えば、最急降下法、ニュートン法のような繰り返し計算による非線形最適化の方法などによって最適化してもよい。すなわち、演算部１５１は、コスト関数が最小化するような各パルスボリュームデータの位置を探索することにより、基準位置に対するパルスボリュームデータの位置ずれ後の位置情報を取得する。
【０１０２】
なお、コスト関数は、想定される各パルスボリュームデータの位置の光照射間の変動（動き）に対して正則化をかけるように定義してもよい。被検体として乳房を考えた場合、呼吸による動きが支配的であると想定される。この場合、被検体の動きは最大で数ｍｍ程度の動きであり、その動きは時間的に連続で滑らかなものであることが想定される。また、その動きは周期的な動きとなることが想定される。上記のように想定される被検体の動きから逸脱するような動きが算出されることに対して抑制を働かせるような正則化をかけることができる。
【０１０３】
正則化のかけ方はいかなる方法であってもよい。例えば、導出過程の被検体の変動量（移動距離）の総和に所定の重み係数をかけてコスト関数に加算することで、正則化することができる。また、被検体の変動の周波数成分値に基づいて算出された値をコスト関数に加算してもよい。また、被検体の典型的な変動のしかたをモデルとして用意し、そのモデルにおける変動との相違をコストとしてコスト関数に加算するようにしてもよい。
【０１０４】
また、「コスト関数を最小化させる」とは、コスト関数が厳密に最小となる場合だけではなく、解の候補を変化させたときにコスト関数の値が所定の値以下となる場合やコスト関数の変化量が所定の値以下となる場合も含む。すなわち、演算部１５１は、コスト関数が所定の条件を満たすことをもって、コスト関数が最小化したと判断してもよい。また、ユーザが入力部１７０を用いて、コスト関数が最小化したことを指示してもよい。この場合、演算部１５１は、入力部１７０からの指示を受けてコスト関数が最小化したと判断する。
【０１０５】
続いて、演算部１５１は、各パルスボリュームデータに対する、コスト関数が最小化したときの位置ずれ量Ｍｏｐｔλ１，ｊを取得する。この位置ずれ量Ｍｏｐｔλ１，ｊは、被検体１００と受信部１２０との相対的な位置関係の変動による、各パルスボリュームデータの位置ずれ量を表す。
【０１０６】
図７は、パルスボリュームデータＰλ１，２の位置ＰｏｓＰλ１，２と、前述の方法でコスト関数が最小化したときのパルスボリュームデータＰ′λ１，２の位置ＰｏｓＰ′λ１，２とを示す。図７において、パルスボリュームデータＰλ１，２を実線で表し、コスト関数が最小化したときのパルスボリュームデータＰ′λ１，２を破線で表す。
【０１０７】
なお、本工程においては、被検体１００と受信部１２０との相対位置の変動による、パルスボリュームデータの位置ずれ情報を取得することができる限り、いかなる手法を用いてもよい。
【０１０８】
なお、演算部１５１は、パルスボリュームデータを２つ以上合成することで得られる波長λ１でのｋ番目の合成ボリュームデータＧλ１，ｋを用いて位置ずれ情報を取得してもよい。その際、合成したパルスボリュームデータをすべて包含し、かつ、最小となる矩形の領域を合成ボリュームデータとしたり、少なくとも２つ以上のパルスボリュームデータが重畳した領域を包含する任意の領域を合成ボリュームデータとしてもよい。すなわち、合成ボリュームデータの領域は、合成対象のパルスボリュームデータを全て包含しなくてもよい。
【０１０９】
合成ボリュームデータは、選択されたパルスボリュームデータを位置に従って加算することで取得する。その他、選択されたパルスボリュームデータを加算し、重なるパルスボリュームの数で除算することで平均化してもよいし、被検体の特徴をより正確に再現したボリュームデータを取得できるかぎり、あらゆる手法を用いることができる。ただし、光照射間の被検体１００と受信部１２０との相対位置の変動を補正する処理（例えば、パルスボリュームの位置を変更する処理）については、本明細書に係る「合成」に含まれない。
【０１１０】
例えば、演算部１５１は、合成対象のパルスボリュームデータのそれぞれを重み付けした後に加算することにより合成してもよい。また、演算部１５１は、外れ値除去法などによってノイズを多く含む値を除外したパルスボリュームデータに対して加算値や平均値を算出するようにしてもよい。
【０１１１】
これらの合成処理により、各パルスボリュームデータに含まれるノイズが低減され、被検体の特徴をより正確に再現した合成ボリュームデータを取得することができる。さらに、パルスボリュームデータの数よりも合成ボリュームデータの数のほうが少ない場合は、パルスボリュームデータ同士を比較することにより全パルスボリュームデータの位置ずれを推定する手法に比べて、計算量と計算コストを少なくすることもできる。
【０１１２】
ただし、このように合成して得られる合成ボリュームデータには、複数回の光照射間の被検体と光音響波の受信部との相対的な位置関係の変動による影響が含まれている。そのため、合成ボリュームデータには当該変動による品質の低下が生じている可能性がある。以下、この品質の低下を抑制するために、合成ボリュームデータの推定位置からパルスボリュームデータの位置を推定する処理を説明する。
【０１１３】
演算部１５１は、波長λ１での各合成ボリュームデータの位置ずれ量ＭＧｏｐｔλ１，ｊに基づいて、各パルスボリュームデータの位置ずれ量Ｍｏｐｔλ１，ｉを推定する。
【０１１４】
演算部１５１は、合成ボリュームデータに対応付けられたパルスボリュームデータの位置ずれ量については、推定された合成ボリュームデータの位置ずれ量割り当てることができる。その他のパルスボリュームデータの位置ずれ量については、演算部１５１が、割り当てられたパルスボリュームデータの位置ずれ量に対して補間処理を行うことにより推定することができる。補間処理の手法について、線形補間やスプライン補間などの公知の手法を採用することができる。また、想定される被検体の動きから逸脱するような位置を算出しないような制約をかけて補間処理を行ってもよい。
【０１１５】
合成の対象となったパルスボリュームデータのうち、任意のパルスボリュームデータを、合成ボリュームデータに対応付けられたパルスボリュームデータとしてもよい。例えば、合成対象のパルスボリュームデータが奇数個である場合、時間的に中心に位置するパルスボリュームデータを合成ボリュームデータに対応付けてもよい。
【０１１６】
また、例えば、本実施形態のように合成対象のパルスボリュームデータが偶数個である場合、時間的に中心付近に位置するいずれかのパルスボリュームデータを合成ボリュームデータに対応付けてもよい。例えば、本実施形態のように１０個のパルスボリュームデータを合成対象とする場合、波長λ１での合成ボリュームデータＧλ１，ｊの位置ずれ量ＭＧｏｐｔλ１，ｊを、パルスボリュームデータＰλ１，５ｊの位置ずれ量ＭＧｏｐｔλ１，５ｊとして割り当ててもよい。
【０１１７】
また、合成対象のパルスボリュームデータが偶数個である場合、時間的に中心に位置する仮想のパルスボリュームデータを合成ボリュームデータに対応付けてもよい。例えば、本実施形態のように１０個のパルスボリュームデータを合成対象とする場合、波長λ１での合成ボリュームデータＧλ１，ｊの位置ずれ量を、パルスインデックスが５．５ｊの仮想パルスボリュームデータの位置ずれ量に割り当ててもよい。
【０１１８】
また、重み付けを伴って合成される場合、合成対象のパルスボリュームデータのうち、最も高い重み係数で重み付けられたパルスボリュームデータを、合成ボリュームデータに対応付けてもよい。また、合成対象のパルスボリュームデータのうち、重み係数が中央値となるパルスボリュームデータを、合成ボリュームデータに対応付けてもよい。
【０１１９】
以上の処理により、波長λ１に対応する合成ボリュームデータの位置ずれ情報に基づいて、波長λ１に対応するパルスボリュームデータの位置ずれ情報を取得することができる。このようにして、波長λ１の光の照射タイミングにおける位置ずれ情報を取得してもよい。
【０１２０】
なお、パルスボリュームデータや合成ボリュームデータの最大値投影画像（Ｍａｘｉｍｕｍ  Ｉｎｔｅｎｓｉｔｙ  Ｐｒｏｊｅｃｔｉｏｎ：ＭＩＰ）等の２次元の投影データを用いて上記と同様に位置ずれ情報を取得してもよい。以下、その処理の一例を説明する。
【０１２１】
演算部１５１は、波長λ１での合成された初期音圧分布データＧλ１，ｊについて、Ｘ方向、Ｙ方向、Ｚ方向のそれぞれの方向に投影した投影データとしてＭＩＰデータを取得する。Ｘ方向に投影したＭＩＰデータはＹ軸とＺ軸によって表される２次元の空間情報であり、ＩＧλ１，ｊ（ｙ，ｚ）と表記される。Ｙ方向に投影したＭＩＰデータはＺ軸とＸ軸によって表される２次元の空間分布情報であり、ＩＧλ１，ｊ（ｚ，ｘ）と表記される。Ｚ方向に投影したＭＩＰデータはＸ軸とＹ軸によって表される２次元の空間分布情報であり、ＩＧλ１，ｊ（ｘ，ｙ）と表記する。
【０１２２】
なお、３次元の画像データを２次元の画像データへと変換できるかぎり、ＭＩＰ画像以外の投影手法を採用してもよい。例えば、ＭＩＰ画像にかえて最小値投影（Ｍｉｎｉｍｕｍ  Ｉｎｔｅｎｓｉｔｙ  Ｐｒｏｊｅｃｔｉｏｎ：ＭＩＮＰ）画像を生成して用いるようにしてもよい。また、投影方向の複数のスライドを加算することにより、投影データを取得してもよい。
【０１２３】
続いて、ＸＹ面、ＹＺ面、ＺＸ面のそれぞれについて、ペアを構成する片方のＭＩＰデータに対するもう一方のＭＩＰデータの相対位置を並進させ、類似度を算出する。類似度算出の手法については、前述した手法を用いることができる。そして、ＸＹ面、ＹＺ面、ＺＸ面のそれぞれについて、Ｐλ１，ｋ１に対するＰλ１，ｋ２の類似度が最大となる並進量ＭＸ＿ｋ、ＭＹ＿ｋ、ＭＺ＿ｋと、それぞれの並進量の各座標軸の成分の平均値を、類似度が最大となる三次元の並進量Ｍの各成分値として算出する。
【０１２４】
【数４】
2018050776.tif 000006
【０１２５】
続いて、演算部１５１は、式５に示す並進量Ｍ＿ｋを用いて、式４に示すコスト関数が最小化したときの各合成ボリュームデータの位置を推定することができる。
【０１２６】
以上の処理により、３次元の画像データから変換された２次元の画像データに基づいて、基準位置に対する合成ボリュームデータの位置を取得することができる。３次元画像データから２次元画像データに変換することにより、３次元画像データのまま処理を行う場合と比べて少ない計算コストで、位置ずれ後のボリュームデータの位置を取得することができる。
【０１２７】
なお、パルスボリュームとしての３次元画像データを２次元の画像データに変換することにより、上記の方法で波長λ１に対応するパルスボリュームデータの位置ずれ情報を取得してもよい。このようにして、波長λ１の光の照射タイミングにおける位置ずれ情報を取得してもよい。
【０１２８】
ここまでは、被検体と受信部との相対的な位置関係の変動として並進が生じる場合の例を説明した。ただし、当該変動として回転や変形が生じる場合についても同様にそれらによる位置ずれ量を推定することができる。
【０１２９】
例えば、回転を考慮する場合、演算部１５１は、並進量に加えて回転量を引数として、各パルスボリュームデータの位置および回転量（位置ずれ量）を推定することができる。続いて、演算部１５１は、推定された位置および回転量に基づいて各パルスボリュームデータを剛体変換処理（位置ずれ補正処理）した後に合成することにより、合成ボリュームデータを取得することができる。なお、回転量のみを位置ずれ量としてもよいし、算出した２次元もしくは３次元の並進・回転行列等の変換行列や変換のための各種パラメータを位置ずれ量としてもよい。
【０１３０】
また、例えば、変形を考慮する際、演算部１５１は、パルスボリュームデータに設定された各点での変位量（並進および回転量の少なくとも一つ）を引数として変位量を推定することができる。続いて、演算部１５１は、推定された変位量に基いて各パルスボリュームデータを変形処理（位置ずれ補正処理）した後に合成することにより、合成ボリュームデータを取得することができる。例えば、Ｆｒｅｅ  Ｆｏｒｍ  Ｄｅｆｏｒｍａｔｉｏｎ（ＦＦＤ）やＴｈｉｎ  Ｐｌａｔｅ  Ｓｐｌｉｎｅなどの変形を表現する手法によって、パルスボリュームデータ間の変位量を算出することができる。これらの処理により、変形を含む高次な変動を考慮して、品質の高い合成ボリュームデータを取得することができる。
【０１３１】
本実施形態では全光照射による光音響波の測定が完了した後にパルスボリュームデータの取得を開始する例を説明したが、光照射のたびに逐次パルスボリュームデータを取得してもよい。後者の場合、取得されたパルスボリュームデータを逐次表示部１６０に表示させてもよい。これにより、ユーザは全測定が完了する前に取得済のパルスボリュームデータを確認することができる。このとき、パルスボリュームデータが重畳した領域については、前述の方法で合成してもよい。
【０１３２】
また、本実施形態では全てのパルスボリュームデータを取得した後に位置ずれ量を算出したが、光照射のたびに得られたパルスボリュームデータを用いて逐次位置ずれ量を算出してもよい。また、逐次得られるパルスボリュームデータから一定の数のパルスボリュームデータを合成し合成ボリュームデータを作成して逐次位置ずれ量を算出してもよい。
【０１３３】
また、位置ずれ量を算出するためのボリュームデータは、負の値を除去したり、画像強度を正規化したり等の事前処理を行ってから使用してもよい。
（Ｓ４００：波長λ１の光の照射タイミングにおける位置ずれ情報に基づいて、波長λ２の光の照射タイミングにおける位置ずれ情報を取得する工程）
演算部１５１は、Ｓ３００で取得された波長λ１の光の照射タイミングにおける位置ずれ情報に基づいて、波長λ２の光の照射タイミングにおける位置ずれ情報を取得する。すなわち、Ｓ３００に記載されたような方法で算出された特定の波長に対応する位置ずれ情報から、その他の波長に対応する位置ずれ情報を算出する。例えば、光の照射タイミングが既知である場合、λ１の光の照射タイミングにおける位置ずれ情報を時間的に補間してもよい。また、光の照射タイミングにおけるプローブ１８０の位置、すなわち受信部１２０の位置が既知である場合、波長λ１の光の照射タイミングにおける位置ずれ情報を空間的に補間してもよい。演算部１５１は、これらの補間により、波長λ２の光の照射タイミングにおける位置ずれ情報を取得してもよい。
【０１３４】
時間的に補間して算出する場合に関して、図８を用いて説明する。図８は、時間ｔにおける複数の波長のパルス光の照射タイミングを表すシーケンス図を含む。また、図８は、波長λ１での各パルス光に対応する位置ずれ量（並進量）を表すグラフを含む。符号８０１の逆三角形は、パルス光の照射タイミングを模式的に表したものである。波長λ１のパルス光は、Ｔ１１、Ｔ２１、Ｔ３１、Ｔ４１、及びＴ５１の時間にそれぞれ照射されていることが理解される。また、波長λ２のパルス光は、Ｔ１２、Ｔ２２、Ｔ３２、Ｔ４２、及びＴ５２の時間にそれぞれ照射されることが理解される。また、Ｓ３００において波長λ１のパルス光の照射タイミングにおける位置ずれ量が算出されているため、その位置ずれ量が図８に丸印８０２でプロットされている。そこで、演算部１５１は、波長λ２に対応する位置ずれ量は波長λ１での位置ずれ量から、時間的に補間して算出する。図８においてバツ印８０３は、時間的に補間して算出された波長λ１のパルス光の照射タイミングにおける位置ずれ量を示す。補間方法はランツォシュ補間を用いるのが理想的に正しいが、その他線形補間、キュービック補間、スプライン補間等、どのような補間方法でもよい。ここでは例として、線形補間、ランツォシュ補間を用いた場合で説明する。例えば、波長λ２でのＴ１２のタイミングでの位置ずれ量は、観測している被検体の動きが波長λ１のパルス照射時間（サンプリング時間）よりも遅い場合、波長λ１に対応するＴ１１とＴ２１での位置ずれ量から算出できる。Ｔ１１、Ｔ２１での位置ずれ量をそれぞれＭ１１、Ｍ２１とした場合、Ｔ１２での位置ずれ量Ｍ１２は、下記式６で算出してもよい。
【０１３５】
【数５】
2018050776.tif 000007
【０１３６】
ここで波長λ２の時刻Ｔ１２のパルスの位置ずれ量Ｍ１２を算出するのに、Ｍ１１とＭ２１を用いて算出を行った。しかし、補間方法によっては、さらに時間的に離れたパルスの位置ずれ量を用いてもよい。例えば、ランツォシュ補間は多変数補間の一種で、デジタル信号のサンプリングレートを向上するのに用いられ、最もよく補間を行えることが知られている。１次元の信号に対するランツォシュ補間は、式２のランツォシュカーネルＬ（ｔ）を用い、式７によって行う。
【０１３７】
【数６】
2018050776.tif 000008
【０１３８】
【数７】
2018050776.tif 000009
【０１３９】
ここでｔは時刻、ａはカーネルサイズを決定する正の整数、ｓｉは正の整数ｉに対する１次元信号のサンプル、Ｓ（ｔ）は補間された値、［  ］は床関数である。線形補間の場合と異なり、最近傍の２つのパルスの位置ずれ量のみならず、近傍にある複数のパルスの位置ずれ量を用いて補間を行っていることがわかる。
【０１４０】
このようにすることで、１つの波長の時間的に近傍にあるパルスの位置ずれ量から、他の波長のパルスの位置ずれ量を算出することができる。また、その過程で、複数の波長の位置ずれ補正を行う必要がなくなるため、計算時間も短縮することができる。
【０１４１】
時間的に補間を行う方法の例として、２つの２次元のＡｆｆｉｎｅ変換行列から補間して１つの２次元のＡｆｆｉｎｅ変換行列を算出する方法を説明する。
【０１４２】
取得した複数波長のパルス光による信号群のうち、波長λ１のパルス光を照射して取得した信号から、対応するパルスボリュームデータを取得する。得られた波長λ１のパルスボリュームデータから、Ａｆｆｉｎｅ変換をベースにした位置ずれ補正を行い、波長λ２での位置ずれ量を算出する。波長λ２での位置ずれ量は、パルスボリュームデータに対してＡｆｆｉｎｅ変換行列として算出される。
【０１４３】
まず、２つの２次元のＡｆｆｉｎｅ変換行列Ａ１とＡ２があるとする。ここで、Ａ１は、ｔ１の時刻に波長λ１のパルス光を照射して得られるパルスボリュームデータＰλ１，０と、ｔ１＋Δの時刻に波長λ１のパルス光を照射して得られるパルスボリュームデータＰλ１，１間での位置ずれ量（Ａｆｆｉｎｅ変換行列）である。また、Ａ２は、ｔ１＋Δの時刻に波長λ２のパルス光を照射して得られるパルスボリュームデータＰλ２，１と、ｔ１＋２Δの時刻に波長λ２のパルス光を照射して得られるパルスボリュームデータＰλ２，２間での位置ずれ量である。
【０１４４】
このとき、求めたいＡｆｆｉｎｅ変換行列Ｂは、ｔ１＋Δ／２の時刻に波長λ１のパルスを照射して得られるパルスボリュームデータＰλ１，１／２と、ｔ１＋３Δ／２の時刻に波長のパルスを照射したパルスボリュームデータＰλ１，３／２間での位置ずれ量であるとする。
【０１４５】
Ａ１はＡ１Ｒという回転成分とＡ１Ｓという拡大縮小成分に分解することができ、Ａ２も同様にＡ２Ｒという回転成分とＡ２Ｓという拡大縮小成分に分解できる。よって、Ｂの回転成分ＢＲはＡ１ＲとＡ２Ｒの行列の各要素を時間的に補間することで、Ｂの拡大縮小成分ＢＳはＡ１ＳとＡ２Ｓの行列の各要素を時間的に補間することで得られる。このようにして得られたＢＲとＢＳを積算することで、位置ずれ量Ｂが算出できる。
【０１４６】
また、空間的に補間を行い、空間的に近傍にあるパルスの位置ずれ量を算出してもよい。図９は、図５と同様に空間的に走査しながら複数の波長のパルスを照射するときのプローブ１８０の位置（測定位置）を表す模式図である。波長λ１に対応する位置ずれ情報を空間的に補間することにより、波長λ２の光照射のタイミングにおける測定位置５０２Ａに対応する位置ずれ情報を取得する場合を説明する。演算部１５１は、測定位置５０２Ａの近傍の測定位置に対応する波長λ１の測定位置を決定する。例えば、波長λ２の測定位置５０２Ａから所定の距離に含まれる波長λ１の測定位置を決定する。符号９０１の一点鎖線は、測定位置５０２Ａから所定の距離の範囲を示す。そして、波長λ１の測定位置５０１Ａ、５０１Ｂ、５０１Ｃ、５０１Ｄ、及び５０１Ｅは、波長λ２の測定位置５０２Ａから所定の距離に含まれる波長λ１の測定位置を示す。よって、演算部１５１は、波長λ１の測定位置５０１Ａ、５０１Ｂ、５０１Ｃ、５０１Ｄ、及び５０１Ｅに対応する位置ずれ情報を空間的に補間することにより、波長λ２の測定位置５０２Ａに対応する位置ずれ情報を取得する。なお、空間的に補間する方法としては、すでに算出されている位置ずれ情報に対して関数のフィッティングを行い、求めたい位置での位置ずれ情報を関数から求めたてもよい。また、空間的に補間する方法としては、距離に応じた重みをつけ、求めたい位置での位置ずれ情報を算出してもよい。以下では、位置ずれ情報を取得したい測定位置から所定の距離内にあるパルスボリュームデータの位置ずれ情報に対して距離に応じた重みをつけることで、算出したいパルスボリュームデータの位置ずれ情報を算出する方法について説明する。
【０１４７】
ここではパルスを照射した位置をパルスボリュームデータの中心の座標とする。また、位置ずれ情報を算出したい波長λ２のパルスボリュームデータＰｒｅｓｕｌｔの中心座標ＰＯＳｒｅｓｕｌｔから、半径Ｒｍｍの範囲に含まれる波長λ１のパルスボリュームデータの位置ずれ情報のパルス数をＮとする。また、位置ずれ情報としての位置ずれ量（並進量）を
【０１４８】
【数８】
2018050776.tif 000010
【０１４９】
、中心座標を
【０１５０】
【数９】
2018050776.tif 000011
【０１５１】
とすると、Ｐｒｅｓｕｌｔの位置ずれ量Ｍｒｅｓｕｌｔは、以下の式で算出できる。
【０１５２】
【数１０】
2018050776.tif 000012
【０１５３】
例えば、図９においては、符号５０２Ａが、位置ずれ情報を算出したい波長λ２のパルスボリュームデータＰｒｅｓｕｌｔの中心座標ＰＯＳｒｅｓｕｌｔに対応する測定位置を表す。また、符号９０１が、測定位置５０２Ａから半径Ｒｍｍの範囲を表す。また、符号５０１Ａ〜Ｅが、範囲９０１に含まれるパルスパルスボリュームデータの中心座標に対応する測定位置を表す。
【０１５４】
ここでωｉは、距離から計算される重みであり、ＰＯＳｉとＰＯＳｒｅｓｕｌｔの間の空間的な距離
【０１５５】
【数１１】
2018050776.tif 000013
【０１５６】
を用いて次の式で計算できる。
【０１５７】
【数１２】
2018050776.tif 000014
【０１５８】
なお、空間的に補間する場合、上記方法に限らず、あらゆるパルスボリュームを用いて空間的に補間してもよい。例えば、再近傍のパルスを用いたり、渦巻状の走査が複数周回している場合は各周回で近傍２つのパルスを用いたり、全パルスを用いたりしてもよい。
【０１５９】
このようにすることで、特定の波長の測定位置と空間的に近傍にある測定位置に対応する位置ずれ情報から、他の波長のパルスの位置ずれ情報を算出することができる。
【０１６０】
３つ以上の波長の光を照射する場合に本実施形態に係る位置ずれ情報の取得方法を適用してもよい。この場合、特定の波長で得られた位置ずれ情報を、その他の複数波長の位置ずれ情報の取得に用いてもよい。また、位置ずれ情報の取得に好適な波長を含む複数の波長で得られた位置ずれ情報を用いて、その他の波長の位置ずれ情報を取得してもよい。すなわち、少なくとも一つ以上の波長により得られた位置ずれ情報を用いて、その他の波長の位置ずれ情報を求めてもよい。
【０１６１】
また、ＦＦＤ等の変形位置合わせ技術を用いて、変形を加味した位置ずれ情報を取得する場合、その位置ずれ情報は変形場でもよい。変形場とは、変形用画像から基準画像への変形の場である。
【０１６２】
位置ずれ情報が変形場の場合の例を、図１０を用いて説明する。ここでは、簡便のために３次元画像データ領域を、ＸＹの２次元平面にして説明する。符号１０１１、符号１０３１は波長λ１でのパルス照射位置、符号１０１２、符号１０３２は前述のパルスの画像データ領域、符号１０５０の点線の領域は符号１０１２と符号１０３２の画像の重なり領域である。また、符号１０２１は波長λ２でのパルス照射位置、符号１０２２は符号１０２１のパルスに対応する画像データ領域である。ここで、符号１０１１、符号１０３１、符号１０２１はそれぞれ、時刻ｔ１、ｔ２、（ｔ１＋ｔ２）／２の時刻に照射されたとする。
【０１６３】
ここで、符号１０１２と符号１０３２の重なり領域１０５０において、符号１０３２の画像を基準として符号１０１２の画像を変形位置合わせを行うことを考える。すると、重なり領域符号１０５０のサイズの、３次元的な変形場を算出することができる。
【０１６４】
そして、符号１０２２の画像データ領域のうち、符号１０５０の領域の変形場を、符号１０３２の画像を基準とした符号１０１２の画像の変形場から算出することができる。算出方法としては、符号１０３２の画像を基準として符号１０１２の画像の変形が、ｔ２−ｔ１の時間で時間の経過に対して線形に発生したものと想定して、変形場を時間で補間することで得られる。より具体的には、（ｔ２−（ｔ１＋ｔ２）／２）／（ｔ２−ｔ１）を積算すればよい。
【０１６５】
この例では、２画像間の変形場から、その中間の時刻での画像の片方の画像からの変形場を算出した。これによると、複数波長で交互に照射されたパルスの画像データそれぞれに対して、１つの波長の中で重なり領域がある２つの画像の撮像時刻の間に撮像した、その他の波長の画像の変形場を算出できる。これを複数波長の全撮像パルスに適用することで、全ての波長でいずれかのパルスの画像データを基準とした各パルスへの変形場を算出できる。
【０１６６】
その他の例としては、複数のパルスの画像データ間で複数の変形場が求められている際に、それら複数の変形場から時間的に補間することで、その近傍の時刻で照射したパルスの変形場を算出することができる。変形場の算出のためには、画像データ中のあるボクセルにおける変形場での変位量に着目し、離散的な時刻で撮像した画像データ間の複数の変形場において、変形場中の各ボクセルがどのように変位していくかを追跡する。そしてその経時的な変位量から、とある時刻での変形場を算出してもよい。
【０１６７】
位置ずれ情報を取得するための波長λ１としての基準は、位置ずれ情報の取得精度の高いものを選択するのが好ましい。そこで、例えば、演算部１５１は、複数波長のそれぞれを用いて生成されたパルスボリュームデータや合成ボリュームデータ、ＭＩＰ画像等の解像度やＳＮ（Ｓｉｇｎａｌ  ｔｏ  Ｎｏｉｓｅ  Ｒａｔｉｏ）などの画質を示す評価指標を算出する。そして、演算部１５１が、算出された画質を示す評価指標が良好であった波長を波長λ１として、位置ずれ情報を取得してもよい。
【０１６８】
また、時間的に補間して位置ずれ情報を取得する場合、演算部１５１は、所定の期間内に多くの回数照射された光の波長を判定し、その波長を波長λ１として設定してもよい。例えば、演算部１５１は、所定の期間内に最も多くの回数が照射された光の波長を、当該所定の期間内においては波長λ１として設定してもよい。また、演算部１５１は、所定の期間内に所定の回数より多く照射された光の波長を、当該所定の期間内においては波長λ１として設定してもよい。なお、波長λ１として複数の波長が設定されてもよい。
【０１６９】
また、波長間の照射タイミング実質的に同じと見なせる場合には、波長λ１の光の照射タイミングにおける位置ずれ情報を、波長λ２の光の照射タイミングにおける位置ずれ情報として取得してもよい。例えば、Ｓ３００で得られた位置ずれ量をその他の波長のパルスにそのまま用いることができる。
【０１７０】
例えば、複数波長間でのパルス光の照射間隔が被検体の動きに対して無視できるほど短い場合などには、波長間の照射タイミングが実質的に同じとみなせる。例えば、生体の呼吸による体動を考えた場合に、波長間の照射タイミングを実質的に同じとみなせる照射間隔を説明する。周期が３秒、最大変位量が３ｍｍの呼吸により生体が線形に変位すると仮定し、音響画像の解像度の許容誤差を０．２５ｍｍと設定した場合、照射間隔が１２５ｍｓ以内であれば、波長間の照射タイミングが実質的に同じであるとみなしてもよい。なお、このとき、生体の変位量が光音響画像の解像度の誤差にそのまま反映されると仮定している。また、光音響画像の解像度の許容誤差を０．１ｍｍと設定すると、照射間隔が５０ｍｓ以内であれば、波長間の照射タイミングが実質的に同じであるとみなしてもよい。
【０１７１】
また、演算部１５１は、Ｓ３００またはＳ４００で得られた位置ずれ情報を、時間的または空間的に補間することにより、その他のモダリティで得られた画像の位置ずれ情報として使用してもよい。
【０１７２】
（Ｓ５００：波長λ１及び波長λ２の光照射タイミングにおける位置ずれ情報に基づいて、位置ずれ補正を行う工程）
演算部１５１は、Ｓ３００で得られた波長λ１に対応する位置ずれ情報、及び、Ｓ４００で得られた波長λ２に対応する位置ずれ情報に基づいて位置ずれ補正を行う。
【０１７３】
例えば、演算部１５１は、Ｓ３００で得られた波長λ１に対応する位置ずれ情報に基づいて、Ｓ２００で取得された波長λ１のパルスボリュームデータの位置情報を位置ずれ補正する処理を行ってもよい。また、演算部１５１は、Ｓ１００で得られた信号群に基づいて、波長λ２に対応するパルスボリュームデータを生成してもよい。そして、演算部１５１は、Ｓ４００で得られた波長λ２に対応する位置ずれ情報に基づいて、波長λ２のパルスボリュームデータの位置を位置ずれ量だけ補正する処理を行ってもよい。さらに、演算部１５１は、位置ずれが補正された波長λ１のパルスボリュームデータ及び波長λ２のパルスボリュームデータを合成することにより位置合わせし、合成ボリュームデータを生成してもよい。
【０１７４】
以下、波長λ１のパルスボリュームの位置ずれを補正する例を説明する。図１１は、本工程における位置ずれ補正処理（並列処理）の一例を示す。
【０１７５】
図１１Ａは、本実施形態における並列処理前のパルスボリュームデータの一部（Ｐλ１，２３〜Ｐλ１，２５）を示す。破線は、Ｐλ１，２５の外周およびＰλ１，２５内の特徴６０１を表す。実線は、Ｐλ１，２４の外周および内部の特徴６０２表す。点線は、Ｐλ１，２３の外周および内部の特徴６０３を表す。なお、特徴６０１、特徴６０２、特徴６０３はいずれも同一の特徴を表している。図１１Ａの状態では、各パルスボリュームデータ内の特徴は異なる位置に位置している。
【０１７６】
図１１Ｂは、並進前のパルスボリュームデータＰλ１，２５を、前述の方法で推定された位置ずれ量Ｍｏｐｔλ１，２５だけ並進させた後のパルスボリュームデータＰ′λ１，２５を示す。図１１Ｃは、並進前のパルスボリュームデータＰλ１，２４を、前述の方法で推定された位置ずれ量Ｍｏｐｔλ１，２４だけ並進させた後のパルスボリュームデータＰ′λ１，２４を示す。
【０１７７】
図１１Ｄは、並進前のパルスボリュームデータＰλ１，２３を、前述の方法で推定された位置ずれ量Ｍｏｐｔλ１，２３だけ並進させた後のパルスボリュームデータＰ′λ１，２３を示す。
【０１７８】
図１１Ｅは、並進後のパルスボリュームデータＰ′λ１，２３、Ｐλ１，２４、およびＰ′λ１，２５を重ねあわせた様子を示す。図１１Ｅにおいては、各パルスボリュームデータ内の特徴６０１、６０２、および６０３はほぼ同じ位置で重なっている。演算部１５１は、図１１Ｅに示すように並進処理された各パルスボリュームデータを合成することにより、位置合わせされたボリュームデータを取得することができる。位置合わせされたボリュームデータは、合成されたボリュームデータに相当する。「位置合わせ」とは、位置ずれ補正処理と合成処理の両処理を行うことを指す。
【０１７９】
また、例えば、演算部１５１は、再構成処理において、Ｓ３００で得られた波長λ１に対応する位置ずれ情報に基づいて、波長λ１の照射タイミングにおける受信部１２０の位置情報を位置ずれ量だけ補正する処理を行ってもよい。また、演算部１５１は、再構成処理において、Ｓ４００で得られた波長λ２に対応する位置ずれ情報に基づいて、波長λ２の照射タイミングにおける受信部１２０の位置情報を位置ずれ量だけ補正する処理を行ってもよい。そして、演算部１５１は、Ｓ１００で得られた信号群、及び、波長λ１及び波長λ２の光の照射タイミングにおける位置ずれが補正された受信部１２０の位置情報に基づいて再構成処理を行い、合成ボリュームデータを生成してもよい。
【０１８０】
ここで、受信部１２０の位置情報の位置ずれ補正を行う方法の例を説明する。まず、演算部１５１は、位置ずれを考慮しないときの受信部１２０の位置情報を取得する。例えば、演算部１５１は、記憶部１５２に予め記憶された光照射時の受信部１２０の位置情報を読み出すことにより、位置ずれを考慮しないときの位置情報を取得してもよい。また、演算部１５１は、光照射をトリガーとして、駆動部１３０に備えられた位置センサから受信部１２０の位置情報を受け取ることにより、位置ずれを考慮しないときの受信部１２０の位置情報を取得してもよい。
【０１８１】
続いて、演算部１５１は、光照射時における位置ずれを考慮しないときの受信部１２０の位置情報を、Ｓ３００またはＳ４００で得られた位置ずれ情報が示す位置ずれ量だけ補正（例えば、並進処理）する。これにより、演算部１５１は、光照射時における位置ずれ補正された受信部１２０の位置情報を取得することができる。すなわち、演算部１５１は、Ｓ３００またはＳ４００で得られた位置ずれ情報に基づいて、位置ずれ補正された受信部１２０の位置情報を取得することができる。
【０１８２】
演算部１５１は、Ｓ１００で得られた信号群と、位置ずれ補正された受信部１２０の位置情報とに基づいて、合成画像データを取得する。本工程では、演算部１５１が、１回の光照射に対応する信号からは全画像化領域よりも小さい領域を再構成し、これを複数の光照射に繰り返すことにより１つのボリュームデータを生成してもよい。この場合、本工程において、演算部１５１は、複数の光照射に対応する複数のパルスボリュームデータを取得し、当該複数の画像データを合成する。また、演算部１５１は、複数の光照射に対応する信号群から全画像化領域を再構成することにより、１つのボリュームデータを生成してもよい。
【０１８３】
なお、複数の波長に対応するデータの合成処理を行う場合にも、上記と同様の処理を行うことができる。また、複数波長に対応するデータの合成処理には、複数波長に対応するデータの加算処理や平均化処理の他に、複数波長に対応するデータを比較演算することにより酸素飽和度などの機能情報を取得する処理も含まれる。以下、複数波長に対応するデータを用いて機能情報としての酸素飽和度を取得する方法を説明する。
【０１８４】
波長λ１と波長λ２では、ヘモグロビン以外の光吸収が無視できるほど低いと仮定すると、波長λ１と波長λ２の吸収係数はそれぞれ、オキシヘモグロビンのモル吸光係数とデオキシヘモグロビンのモル吸収係数を用いて式１１、式１２のように表わされる。
μａ（λ１）＝εｏｘ（λ１）Ｃｏｘ＋εｄｅ（λ１）Ｃｄｅ・・・式１１
μａ（λ２）＝εｏｘ（λ２）Ｃｏｘ＋εｄｅ（λ２）Ｃｄｅ・・・式１２
【０１８５】
ここで、μａ（λ１）は、位置（ｉ、ｊ、ｋ）における波長λ１の光の吸収係数、μａ（λ２）は位置（ｉ、ｊ、ｋ）における波長λ２の光の吸収係数を示し、単位は［ｍｍ−１］で示すことができる。Ｃｏｘはオキシヘモグロビンの量［ｍｏｌ］、Ｃｄｅはデオキシヘモグロビンの量［ｍｏｌ］である。いずれも位置（ｉ、ｊ、ｋ）における値を示すものとする。
【０１８６】
εｏｘ（λ１）とεｄｅ（λ１）はそれぞれ波長λ１におけるオキシヘモグロビン、デオキシヘモグロビンのモル吸収係数［ｍｍ−１ｍｏｌ−１］を示す。εｏｘ（λ２）とεｄｅ（λ２）はそれぞれ波長λ２におけるオキシヘモグロビン、デオキシヘモグロビンのモル吸収係数［ｍｍ−１ｍｏｌ−１］を示す。εｏｘ（λ１）、εｄｅ（λ１）、εｏｘ（λ２）、εｄｅ（λ２）はあらかじめ測定や文献値によって得ることができる。
【０１８７】
よって、ＣｏｘとＣｄｅはそれぞれ、モル吸光係数と、μａ（λ１）及びμａ（λ２）と、を用いて式１１、１２の連立方程式を解くことにより求められる。用いる波長の数が多い場合は、最小二乗法を用いるとよい。また、酸素飽和度ＳＯ２は式１３に示すように、全ヘモグロビン中のオキシヘモグロビンの割合で定義される。よって、酸素飽和度ＳＯ２は、式１１、１２、１３に基づき、式１４で示すことができる。
【０１８８】
よって、演算部１５１は、式１４にしたがって、モル吸光係数と、μａ（λ１）及びμａ（λ２）とに基づき、位置（ｉ、ｊ、ｋ）における酸素飽和度ＳＯ２を得ることができる。
【０１８９】
【数１３】
2018050776.tif 000015
【０１９０】
【数１４】
2018050776.tif 000016
【０１９１】
このような処理を複数位置に対して行うことにより、複数位置での酸素飽和度が求められ、酸素飽和度分布を取得できる。酸素飽和度分布は、吸収係数分布の比較演算（例えば、非を算出する処理）で得られるものであり、複数の波長での吸収係数の値が相対的に正しければ、酸素飽和度分布は適切に求められる。よって、吸収係数分布が絶対値として正確に求まっていなくてもよい。
【０１９２】
ここでは、位置ずれの影響が低減した初期音圧画像を算出し、そこから吸収係数画像、酸素飽和度画像を算出した。なお、その方法に限らず、複数の波長でのパルスボリュームデータを吸収係数分布画像として求め、そのうちの一つの波長での吸収係数分布画像で位置ずれ量を算出してもよい。また、そのようにして算出された位置ずれ量をその他の波長に適用することで、位置ずれの影響が低減した吸収係数画像を算出し、それらの吸収係数画像から酸素飽和度画像を算出してもよい。
【０１９３】
このようにすることで、複数の波長での画像間で、位置ずれの影響が低減した吸収係数画像を取得し、さらに位置ずれの影響が低減された吸収係数画像を合成することにより、酸素飽和度画像を算出することが可能となる。
【０１９４】
本工程では、光照射間の被検体と光音響波の受信部との相対的な位置関係の変動の影響が抑制された画像データ（位置合わせされたボリュームデータ）を取得することができる。
【０１９５】
なお、光照射間で受信部１２０が移動しない場合であっても上記処理を適用することができる。すなわち、光音響装置が駆動部１３０を備えない場合であっても上記処理を適用することができる。この場合も、複数回の光照射間の被検体と光音響波の受信部との相対的な位置関係の変動の影響が抑制された画像データを取得することができる。
【０１９６】
また、光音響装置とは異なるモダリティ（例えば、超音波診断装置など）で、光音響装置を用いた撮像と併せて撮像を行う場合、光音響装置で得られた位置ずれ情報を他のモダリティにおける位置ずれ情報の取得に用いてもよい。例えば、演算部１５１は、Ｓ３００、Ｓ４００、またはＳ５００で得られた位置ずれ情報を時間的または空間的に補間することにより、他のモダリティで得られた画像データの位置ずれ情報を取得してもよい。また、光音響装置と他のモダリティとが実質的に同じ時間に実質的に同じ位置の画像データを生成する場合、光音響装置で得られた位置ずれ情報を他のモダリティでの位置ずれ情報としてもよい。
【０１９７】
例えば、他のモダリティとして超音波診断装置を想定する場合、光音響装置を含む検査システムは、超音波の送受信部を備える。送受信部は、被検体に対して超音波の送信し、送信された超音波のエコー波を受信することにより超音波信号を出力する。送受信部は、音響波を受信することにより電気信号を出力するトランスデューサを含む。送受信部は複数のトランスデューサを含んでいてもよい。なお、送受信部は、超音波を送信するトランスデューサと、音響波を受信するためのトランスデューサとを別に用意してもよい。また、超音波を送信するトランスデューサと、音響波を受信するためのトランスデューサとが、同じトランスデューサで構成されていてもよい。また、超音波を送受信するためのトランスデューサと、光音響波を受信するためのトランスデューサとを別に用意してもよい。また、超音波を送受信するトランスデューサと光音響波を受信するトランスデューサとが、同じトランスデューサで構成されていてもよい。
【０１９８】
［第２の実施形態］
以下、第２の実施形態の光音響装置の構成の構成及び処理について説明する。第２の実施形態では、第１の実施形態の光音響装置と同様の装置を用いる。第２の実施形態においては、第１の実施形態の光音響装置と同様の構成については同じ符号を付し、詳細な説明を省略する。
【０１９９】
前述したように、複数の波長の光を用いる場合、複数の画像データ間の位置ずれの推定精度が低い波長が含まれている可能性がある。この場合、その波長により得られた複数の画像データのみを用いると、その波長の光の照射タイミングにおける位置ずれ情報を精度よく取得することができない。すなわち、複数の波長の光を用いる場合、波長間での位置ずれの推定精度にバラつきが生じてしまう。
【０２００】
そこで、本発明の一実施形態に係る光音響装置は、互いに異なる複数の波長のそれぞれについて、複数の画像データ間の位置ずれ情報を取得する。そして、複数の波長のそれぞれで得られた位置ずれ情報を組み合わせることにより、先に得られた位置ずれ情報を更新する。これにより、推定精度が低くなってしまう波長が含まれている場合であっても、位置ずれの推定精度が比較的高い波長により得られた位置ずれ情報も用いて、得られた位置ずれ情報を更新するため、位置ずれ情報を精度良く取得することができる。
【０２０１】
すなわち、本発明の一実施形態に係る光音響装置においては、互いに異なる第１の波長及び第２の波長のそれぞれの光が複数回被検体に照射される。そして、第１の波長に対応する第１の複数の画像データが生成され、第１の複数の画像データ間の位置ずれ情報が取得される。ここで取得された位置ずれ情報が、第１の波長の光の照射タイミングにおける被検体とプローブとの相対位置の変動量（位置ずれ量）に相当する。また、第２の波長に対応する第２の複数の画像データが生成され、第２の複数の画像データ間の位置ずれ情報が取得される。そして、第１の位置ずれ情報と第２の位置ずれ情報とに基づいて、第１及び第２の波長の光照射タイミングにおける位置ずれ情報が更新される。これにより、位置ずれの推定精度の低い波長の位置ずれ情報の取得の際に、位置ずれの推定精度の高い波長により得られた位置ずれ情報も用いるため、波長間での位置ずれの推定精度のバラつきが小さくなる。
【０２０２】
さらに、それぞれの波長の位置ずれ情報を組み合わせるときの位置ずれ情報に対する重みを変更することにより、推定精度の高い波長の推定結果を優先させることもできる。
【０２０３】
なお、信号群の取得や各波長の光照射タイミングにおける位置ずれ情報の取得等の方法については、第１の実施形態に記載の方法を適用することができる。
【０２０４】
以下、図１２に示すフローチャートに沿って、本実施形態に係る情報処理を含む光音響装置の作動を説明する。図４に示すステップと同様のステップには同一の符号を付し、詳細な説明を省略する。
【０２０５】
（Ｓ２１０：波長λ２に対応する複数の画像データを生成する工程）
演算部１５１は、Ｓ２００と同様に、Ｓ１００で取得された信号群に基づいて、第２の波長の光照射に対応する複数の画像データを取得する。
【０２０６】
（Ｓ３１０：波長λ２の光の照射タイミングにおける位置ずれ情報を取得する工程）
演算部１５１は、Ｓ３００と同様に、Ｓ２１０で取得された波長λ２に対応する複数の画像データに基づいて、波長λ２の光の照射タイミングにおける位置ずれ情報を取得する。
（Ｓ６００：波長λ１及び波長λ２の光の照射タイミングにおける位置ずれ情報に基づいて、位置ずれ情報を更新する）
演算部１５１は、Ｓ３００で得られた波長λ１の光の照射タイミングにおける位置ずれ情報と、波長λ２の光の照射タイミングにおける位置ずれ情報とに基づいて、それぞれの位置ずれ情報を更新する。
【０２０７】
例えば、演算部１５１が、Ｓ３００において波長λ１の光の照射タイミングにおける位置ずれ情報として並進量Ｍｏｐｔ＿ｓｔλ１を算出し、Ｓ３１０において波長λ２の光の照射タイミングにおける位置ずれ情報として並進量Ｍｏｐｔ＿ｓｔλ２を算出した場合を考える。このとき、本工程において、演算部１５１は、式１５に示すように位置ずれ情報を平均化することにより、各波長の位置ずれ情報Ｍｏｐｔ＿ｓｔλ１及びＭｏｐｔ＿ｓｔλ２を更新してもよい。
Ｍｏｐｔ＿ｓｔλ１＝（Ｍｏｐｔλ１＋Ｍｏｐｔλ２）／２
Ｍｏｐｔ＿ｓｔλ２＝（Ｍｏｐｔλ１＋Ｍｏｐｔλ２）／２・・・式１５
【０２０８】
その他、演算部１５１は、各波長の位置ずれ情報を用いて、第１の実施形態で説明したように、時間的または空間的に補間することにより、各波長の位置ずれ情報を更新してもよい。
【０２０９】
また、演算部１５１は、各波長の位置ずれ情報に重みづけを行うことにより、各波長の位置ずれ情報を更新してもよい。演算部１５１は、所定の重みを用いて各波長の位置ずれ情報を重みづけすることにより、各波長の位置ずれ情報を更新してもよい。また、演算部１５１は、ユーザが入力部１７０を用いて行った指示によって決定された重みを用いて各波長の位置ずれ情報を重みづけしてもよい。
【０２１０】
（Ｓ７００：更新された位置ずれ情報に基づいて、位置ずれ補正を行う工程）
演算部１５１は、Ｓ６００で更新された、波長λ１及び波長λ２の光の照射タイミングにおける位置ずれ情報に基づいて、Ｓ５００で説明した同様の方法で、位置ずれ補正を行う。
【０２１１】
本実施形態では、位置ずれの推定精度の低い波長により得られた位置ずれ情報だけでなく、位置ずれの推定精度の高い波長により得られた位置ずれ情報も用いて組み合わせられた位置ずれ情報を利用するため、波長間での位置ずれの推定精度のバラつきが小さくなる。
【０２１２】
図１３は、各波長の位置ずれ情報に対する重みを決定するためのスライダーを有するＧＵＩの模式図を示す。例えば、表示部１６０に表示される図１３に示すようなＧＵＩを用いて重みを決定することができる。
【０２１３】
ＧＵＩ１００１には、光音響画像の表示領域１００２、波長λ１の位置ずれ量のグラフ１００３、及び波長λ２の位置ずれ量のグラフ１００４が表示されている。グラフ１００３及び１００４の縦軸は位置ずれ量、横軸はパルスインデックスを示している。すなわち、グラフ１００３及び１００４は、波長λ１及び波長λ２でのパルスボリュームデータに対する位置ずれ量をプロットしたものである。スライダーバー１００７は、複数波長の位置ずれ量の重みを決定するためのアイテムである。スライダーバー１００７を左に操作すると波長λ１の位置ずれ量が優位に、右に操作すると波長λ２の位置ずれ量が優位に足しあわされる。すなわち、スライダーバー１００７を左に操作すると波長λ１の重みが大きくなり、右に操作すると波長λ２の重みが大きくなり、重みづけされた最終的な位置ずれ量が決定される。そして、重みづけされた最終的な位置ずれ量を用いて得られる光音響画像が表示領域１００２に表示される。さらに、ユーザがスライダーバー１００７を操作すると、最終的な位置ずれ量が再計算され、更新される。そして、更新された位置ずれ量を用いて再度光音響画像が取得され、表示領域１００２に表示される光音響画像が更新される。図１３においては、表示領域１００２に足の光音響画像１００８が表示されており、光音響画像１００８には血管像１００９が含まれている。
【０２１４】
このようなＧＵＩを用いることにより、複数の波長での位置ずれ量に対する重みをスライダーバーで変化させて決定することで、表示領域に表示される光音響画像の画質（解像度等）の変化を確認しながら最終的な位置ずれ量を決定することができる。
【０２１５】
［第３の実施形態］
以下、第３の実施形態の光音響装置の構成の構成及び処理について説明する。第３の実施形態では、第１また第２の実施形態の光音響装置と同様の装置を用いる。第３の実施形態においては、第１または第２の実施形態の光音響装置と同様の構成については同じ符号を付し、詳細な説明を省略する。
【０２１６】
前述したように、互いに異なる複数の波長の光を用いる場合、波長間で画像強度等の画像特性が異なる場合がある。その場合、波長間で画像強度が異なるために位置ずれの推定精度が低下してしまう場合がある。
【０２１７】
そこで、本発明の一実施形態に係る光音響装置は、互いに異なる複数の波長のそれぞれの光照射により発生する光音響波に基づいて、各波長の光による画像データを生成する。そして、波長間の画像特性の差が小さくなるような処理を行った後に、複数波長の複数の画像データを用いて位置ずれ情報を取得する。これにより、波長間の画像特性の違いにより生じる位置ずれの推定精度の低下を抑制することができる。
【０２１８】
すなわち、本発明の一実施形態に係る光音響装置においては、互いに異なる第１の波長及び第２の波長のそれぞれの光が複数回被検体に照射される。そして、第１の波長に対応する第１の複数の画像データが生成される。また、第２の波長に対応する第２の複数の画像データが生成される。そして、第１の複数の画像データと第２の複数の画像データとの画像特性の差が小さくなるように、第１の複数の画像データ及び第２の複数の画像データの少なくとも一方に画像処理が行われる。そして、画像特性の差が小さくなる処理が行われた後の、第１の複数の画像データ及び第２の複数の画像データを用いて、第１及び第２の波長の光照射タイミングにおける位置ずれ情報が取得される。
【０２１９】
なお、上記では、画像データに対する画像処理により、波長間の画像特性の差を小さくする例を説明したが、画像データとする前の信号に対する信号処理により結果的に波長間の画像特性の差が小さくしてもよい。
【０２２０】
以下、図１４に示すフローチャートに沿って、本実施形態に係る情報処理を含む光音響装置の作動を説明する。図４に示すステップと同様のステップには同一の符号を付し、詳細な説明を省略する。
（Ｓ８００：波長λ１に対応する複数の画像データと波長λ２に対応する複数の画像データとの画像特性の差を小さくする処理を行う）
演算部１５１は、Ｓ２００で得られた波長λ１の複数の画像データと、Ｓ２１０で得られた波長λ２の複数の画像データとの画像特性の差が小さくなるように、波長λ１の複数の画像データ及び波長λ２の複数の画像データに対して画像処理を行う。例えば、波長間の画像特性の差を小さくする処理としては、画像強度の最大値または最小値を揃える処理や複数波長の画像強度の画像平均や分散が同程度となるような処理を行うことが挙げられる。なお、本明細書においては、画像データの画像特性の差を小さくする処理のことを正規化処理と呼ぶ。
【０２２１】
ここでは、画像データの画像強度を正規化する例を説明する。まず、Ｓ２００及びＳ３００で得られた全パルスボリュームデータの画像強度の最大値ｖａｌｍａｘ（全パルスに対して１つの値）を求める。そして、パルスボリュームデータの画像強度のうち０以下の値を０に丸め、その後、パルスボリュームデータのボクセルの画像強度Ｐｘ，ｙ，ｚを、求められた画像強度の最大値ｖａｌｍａｘが所定値ｖａｌになるように正規化する。すなわち、演算部１５１は、式１６に示すようにパルスボリュームデータを正規化する。
Ｐ′ｘ，ｙ，ｚ＝Ｐｘ，ｙ，ｚ＊ｖａｌ／ｖａｌｍａｘ・・・式１６
【０２２２】
ここでＰ′ｘ，ｙ，ｚは、各ボクセルの正規化後の画像強度の値である。ここでは複数波長のパルスボリュームデータの最大値を所定値とするように正規化したが、各波長について最大値ｖａｌｍａｘを求め、最大値ｖａｌｍａｘが所定値ｖａｌになるように正規化してもよい。また、０以下の値を０に丸めず、パルスボリュームデータの最小値を取得し、最小値と最大値の間の強度を０からｖａｌに変換したり、最小値と最大値の間の強度をｖａｌ’からｖａｌに変換したりしてもよい。すなわち、画像強度が所望の数値範囲に含まれるように正規化してもよい。なお、波長間の画像強度の差が小さくなる限り、どのような方法で正規化してもよい。
【０２２３】
また、その他の手法としては、画像の強度の平均が０、分散が１などの、特定の値になるように画像を正規化してもよいし、それぞれの値を０、１ではなく、特定の値になるように正規化してもよい。
【０２２４】
（Ｓ９００：画像特性の差を小さくする処理後の画像データに基づいて位置ずれ情報を取得する工程）
演算部１５１は、Ｓ８００で波長間の画像特性の差を小さくする処理を行った後の、波長λ１及び波長λ２の複数の画像データに基づいて、第１及び第２の波長の光照射タイミングにおける位置ずれ情報を取得する。画像データを用いた位置ずれ情報の取得方法については、Ｓ３００で説明した方法と同様の方法を採用することができる。
【０２２５】
このように、波長間の画像特性の差が低減された複数の画像データを用いることにより、波長間の画像特性の差に起因する位置ずれの推定精度の低下を抑制することができる。そして、このようにして得られた位置ずれ情報を用いて、精度良く位置ずれ補正を行うことができる（Ｓ５００）。
【０２２６】
なお、本実施形態では、各波長において複数の画像データが生成される例を説明したが、１つの波長では１つの画像データが生成されない場合であっても、複数の画像データが生成される限り、本実施形態で説明した位置ずれ情報の取得方法を適用することができる。
【０２２７】
［第４の実施形態］
以下、第４の実施形態の光音響装置の構成の構成及び処理について説明する。第４の実施形態では、第１、第２、または第３の実施形態の光音響装置と同様の装置を用いる。第４の実施形態においては、第１、２、または第３の実施形態の光音響装置と同様の構成については同じ符号を付し、詳細な説明を省略する。
【０２２８】
本実施形態では、入力部１７０を用いたユーザの指示に基づいて、第１、第２、及び第３で説明した位置ずれ情報の取得方法の少なくとも一つの方法を実行する例を説明する。
【０２２９】
図１５は、本実施形態の表示部１６０に表示されるＧＵＩを示す。ＧＵＩ１５０１には、光音響画像の表示領域１５０２、波長λ１の位置ずれ量のグラフ１５０３、及び波長λ２の位置ずれ量のグラフ１５０４が表示されている。グラフ１５０３及び１５０４の縦軸は位置ずれ量、横軸はパルスインデックスを示している。すなわち、グラフ１５０３及び１５０４は、波長λ１及び波長λ２でのパルスボリュームデータに対する位置ずれ量をプロットしたものである。図１５においては、表示領域１５０２に足の光音響画像１５０８が表示されており、光音響画像１５０８には血管像１５０９が含まれている。
【０２３０】
ユーザは、位置ずれ情報の取得手法の選択ボタン１５１１、１５１２、及び１５１３の中から入力部１７０を用いて所望の取得方法を選択する。
【０２３１】
選択ボタン１５１１は、第１の実施形態で説明した特定の波長に対応する位置ずれ情報を、その他に対応する位置ずれ情報の取得に用いる方法に対応するラジオボタンである。本実施形態においては、選択ボタン１５１１が選択されると、波長λ１に対応する位置ずれ情報を補間することにより、波長λ２に対応する位置ずれ情報を取得するモードとなる。さらに、補間方法についても、ユーザが選択ボタン１５１４を選択することにより、演算部１５１が時間的に補間するのか、空間的に補間するのかを決定することができる。
【０２３２】
また、選択ボタン１５１２は、第２の実施形態で説明した複数波長のそれぞれに対応する位置ずれ情報を合成し、各位置ずれ情報を更新する方法に対応するラジオボタンである。本実施形態においては、選択ボタン１５１２が選択されると、波長λ１に対応する位置ずれ情報と波長λ２による位置ずれ情報とを重みづけした後に合成するモードとなる。さらに、ユーザがスライダーバー１５１５を操作することにより、各位置ずれ情報の重みを変更することができる。なお、スライダーバー１５１５の機能は、スライダーバー１００７と同様である。
【０２３３】
また、選択ボタン１５１３は、第３の実施形態で説明した波長間の画像特性の差が小さくなるように正規化処理された後に位置ずれ情報を取得する方法に対応するラジオボタンである。本実施形態においては、選択ボタン１５１３が選択されると、波長λ１の光に基づいて得られたデータ及び波長λ２の光に基づいて得られたデータに対して正規化処理を行うモードとなる。さらに、正規化処理の方法についても、ユーザが選択ボタン１５１６を選択することにより、演算部１５１が画像強度の差を小さくする正規化処理を行うのか、解像度の差を小さくする正規化処理を行うのかを決定することができる。
【０２３４】
このように、ユーザが所望の位置ずれ情報の取得方法を選択し、選択された位置ずれ情報が適用された光音響画像を確認することにより、画像データの特性に適した位置ずれ情報の取得方法で得られた合成画像データを確認することができる。
【０２３５】
［その他の実施形態］
また、本発明は、以下の処理を実行することによっても実現される。即ち、上述した実施形態の機能を実現するソフトウェア（プログラム）を、ネットワーク又は各種記憶媒体を介してシステム或いは装置に供給し、そのシステム或いは装置のコンピュータ（またはＣＰＵやＭＰＵ等）がプログラムを読み出して実行する処理である。
【符号の説明】
【０２３６】
１１０  光照射部
１２０  受信部
１５０  コンピュータ
【図１】
2018050776.tif 000017
【図２】
2018050776.tif 000018
【図３】
2018050776.tif 000019
【図４】
2018050776.tif 000020
【図５】
2018050776.tif 000021
【図６】
2018050776.tif 000022
【図７】
2018050776.tif 000023
【図８】
2018050776.tif 000024
【図９】
2018050776.tif 000025
【図１０】
2018050776.tif 000026
【図１１】
2018050776.tif 000027
【図１２】
2018050776.tif 000028
【図１３】
2018050776.tif 000029
【図１４】
2018050776.tif 000030
【図１５】
2018050776.tif 000031