(57)【要約】      （修正有）【課題】従来に比較して高精度で対象物を画像化する技術を提供する。【解決手段】対象物に対して光を照射して対象物から発生する超音波を複数のセンサ２２で検出し、音圧空間分布を示す複数の検出信号からなる音圧データに基づいて対象物を画像化する画像化処理を行う光超音波画像化装置１０において、所定の初期音速データに基づいて、音速データのプロジェクション行列に初期の音圧ベクトルを乗算したときの乗算結果を音圧データとするときに、初期の音圧ベクトルに関する目的関数であって、対象物の画像データのスパース性、もしくは、画像データに対して所定の変換を行ったときの変換の基底のデータのスパース性を表す目的関数の値が所定のしきい値以下になるように、音速が異なる複数の領域の音速を含む音速ベクトルを変化させて、音圧ベクトルと音速ベクトルを同時に推定することで、音速ベクトルに基づく画像化された画像データを計算する。【選択図】図１
2018161318.tif 000002
【特許請求の範囲】
【請求項１】
光超音波画像化法により、対象物に対して光を照射して前記対象物から発生する超音波を複数のセンサで複数の検出信号として検出し、音圧空間分布を示す複数の検出信号からなる音圧データに基づいて前記対象物を画像化する画像化処理を行う制御手段を備えた光超音波画像化装置において、
前記制御手段は、所定の初期音速データに基づいて、音速データのプロジェクション行列に初期の音圧ベクトルを乗算したときの乗算結果を音圧データとするときに、前記初期の音圧ベクトルに関する目的関数であって、前記対象物の画像データのスパース性、もしくは、前記画像データに対して所定の変換を行ったときの変換の基底のデータのスパース性を表す目的関数の値が所定のしきい値以下になるように、音速が異なる複数の領域の音速を含む音速ベクトルを変化させて、音圧ベクトルと音速ベクトルを同時に推定することで、前記音速ベクトルに基づく画像化された画像データを計算することを特徴とする光超音波画像化装置。
【請求項２】
前記制御手段は、前記目的関数の値が実質的に最小になるように、前記音速ベクトルに基づく画像化された画像データを計算することを特徴とする請求項１記載の光超音波画像化装置。
【請求項３】
前記目的関数は、前記対象物の画像データ、もしくは前記画像データに対して所定の変換を行ったときの変換の基底のデータのＬ１ノルム、Ｌ２ノルム、又は全変動、又はこれらの組み合わせであることを特徴とする請求項１又は２記載の光超音波画像化装置。
【請求項４】
前記制御手段は、音圧ベクトルと音速ベクトルを同時に推定する途中の前記対象物の画像データ、もしくは、前記画像化後の画像データを表示部に出力することで可視化することを特徴とする請求項１〜３のうちのいずれか１つに記載の光超音波画像化装置。
【請求項５】
前記制御手段は、前記光超音波画像化法の空間解像度と時間解像度のうちの少なくとも一方を昇順で変化させて、前記光超音波画像化法により画像化処理をそれぞれ実行し、各画像化処理で得られた音速ベクトルを次の初期音速ベクトルとして用いて画像化処理を行うことを特徴とする請求項１〜４のうちのいずれか１つに記載の光超音波画像化装置。
【請求項６】
前記対象物は互いに音速が異なる第１及び第２の領域を含み、
前記制御手段は、
（１）所定の初期音速ベクトルを用いて、所定の画像化法を用いて前記対象物を画像化して画像データを計算し、
（２）前記計算された画像データに対して所定の画像処理を行うことで、前記第２の領域を推定し、
（３）所定の位置合わせ法により、前記推定した第２の領域を、前記対象物の各領域の音速初期値を地図で示す所定の音速初期値アトラスに位置合わせすることで初期音速マップを計算し、
（４）前記位置合わせされた初期音速マップを用いて、請求項１〜５のうちのいずれか１つに記載の画像化処理を行うことを特徴とする光超音波画像化装置。
【請求項７】
光超音波画像化法により、対象物に対して光を照射して前記対象物から発生する超音波を複数のセンサで複数の検出信号として検出し、音圧空間分布を示す複数の検出信号からなる音圧データに基づいて前記対象物を画像化する画像化処理を行う制御手段を備えた光超音波画像化装置のための光超音波画像化方法において、
前記制御手段が、所定の初期音速データに基づいて、音速データのプロジェクション行列に初期の音圧ベクトルを乗算したときの乗算結果を音圧データとするときに、前記初期の音圧ベクトルに関する目的関数であって、前記対象物の画像データのスパース性、もしくは、前記画像データに対して所定の変換を行ったときの変換の基底のデータのスパース性を表す目的関数の値が所定のしきい値以下になるように、音速が異なる複数の領域の音速を含む音速ベクトル又はその変化量変化させて、音圧ベクトルと音速ベクトルを同時に推定することで、前記音速ベクトルに基づく画像化された画像データを計算するステップを含むことを特徴とする光超音波画像化方法。
【請求項８】
前記制御手段が、前記目的関数の値が実質的に最小になるように、前記音速ベクトルに基づく画像化された画像データを計算するステップを含むことを特徴とする請求項７記載の光超音波画像化方法。
【請求項９】
前記目的関数は、前記対象物の画像データ、もしくは前記画像データに対して所定の変換を行ったときの変換の基底のデータのＬ１ノルム、Ｌ２ノルム、又は全変動、又はこれらの組み合わせであることを特徴とする請求項７又は８記載の光超音波画像化方法。
【請求項１０】
前記制御手段が、音圧ベクトルと音速ベクトルを同時に推定する途中の前記対象物の画像データ、もしくは、前記画像化後の画像データを表示部に出力することで可視化するステップを含むことを特徴とする請求項７〜９のうちのいずれか１つに記載の光超音波画像化方法。
【請求項１１】
前記制御手段が、前記光超音波画像化法の空間解像度と時間解像度のうちの少なくとも一方を昇順で変化させて、前記光超音波画像化法により画像化処理をそれぞれ実行し、各画像化処理で得られた音速ベクトルを次の初期音速ベクトルとして用いて画像化処理を行うことを特徴とする請求項７〜１０のうちのいずれか１つに記載の光超音波画像化方法。
【請求項１２】
前記対象物は互いに音速が異なる第１及び第２の領域を含み、
前記制御手段が、
（１）所定の初期音速ベクトルを用いて、所定の画像化法を用いて前記対象物を画像化して画像データを計算し、
（２）前記計算された画像データに対して所定の画像処理を行うことで、前記第２の領域を推定し、
（３）所定の位置合わせ法により、前記推定した第２の領域を、前記対象物の各領域の音速初期値を地図で示す所定の音速初期値アトラスに位置合わせすることで初期音速マップを計算し、
（４）前記位置合わせされた初期音速マップを用いて、請求項７〜１１のうちのいずれか１つに記載の画像化処理を行うステップを含むことを特徴とする光超音波画像化方法。
【請求項１３】
請求項７〜１２のうちのいずれか１つに記載の光超音波画像化方法の各ステップを含むことを特徴とする、コンピュータにより実行可能な光超音波画像化装置の制御プログラム。
【請求項１４】
請求項１３記載の光超音波画像化装置の制御プログラムを格納することを特徴とする、コンピュータにより読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、光超音波画像化法を用いて検査対象物を画像化する光超音波画像化装置及び方法、光超音波画像化装置の制御プログラム並びに記録媒体に関する。
【背景技術】
【０００２】
光音響イメージング（Photoacoustic Imaging: 以下、ＰＡイメージング又はＰＡ画像化ともいう）は、癌、腫瘍血管新生物、及びその他の多くの疾患の早期臨床診断のための有望な新技術である（例えば、非特許文献１参照）。ＰＡイメージングは、造影剤なしでも高空間分解能で、生体の血管を非侵襲的に視覚化することができ、多くの医療用途に有用である。ここで、ＰＡイメージングは光音響効果を利用する。すなわち、画像化された組織（すなわち、血管）は、短パルス近赤外線照射のレーザエネルギーを吸収して熱に変換し、熱弾性膨張と超音波の放出をもたらす。放射された光音響波を検知することにより物体の３次元構造を再構成することができる（例えば、非特許文献２参照）。
【０００３】
再構成問題には以下の２つの大きな困難性がある。第１の困難性として、ＰＡイメージングのデータ取得速度は、レーザの繰り返し速度と超音波センサの数によって制限される。サンプリングの回数はリアルタイム撮像のために制限され、システムコストが低減される。この不十分なサンプリングは、通常、ストリーキングアーチファクトを引き起こす。第２の困難性として、ほとんどのＰＡ再構成方法で音速が正しく与えられていると仮定し、再構成精度は、与えられた速度の精度に大きく依存する。音速はしばしば温度と組織の変化に伴って変化し、それに応じてそれを調整することが不可欠である。
【０００４】
例えば、特許文献１では、略球冠形状または略球帯形状に沿って多数の振動子を配置して構成される探触子を用いて、高精度な被検体情報を取得するための被検体情報取得装置が開示されている。この被検体情報取得装置は、測定対象からの音響波を受信して電気信号を生成する複数の振動子と、複数の振動子の少なくとも一部の指向軸が集まるように支持する支持体と、点音源と、測定対象と支持体の相対位置を制御する手段と、複数の振動子それぞれから測定対象の注目位置までの距離を求めて注目位置における特性情報を生成する処理手段と、点音源から複数の振動子への音響波の受信時刻と、音速を用いて、複数の振動子それぞれの位置情報の補正データを求める補正手段を有する被検体情報取得装置を用いることを特徴としている。
【０００５】
また、特許文献２では、超音波イメージングに用いる音速を最適化することで、従来に比して高分解能な超音波画像を取得することができる超音波イメージング装置を提供することが開示されている。当該超音波イメージング装置では、分解能最適化ユニットは、走査断面内の各位置毎の組織成分に応じた最適音速を判定し、この最適音速を用いて、走査断面内の各位置からの受信ビーム毎の受信遅延時間等を計算する。制御プロセッサは、最適音速を用いて計算される受信遅延時間を用いて、実際に診断に用いられる超音波画像を取得するためのスキャンにおける遅延加算処理を実行する。これにより、受信遅延時間の計算に用いる設定音速と実際の生体内音速とのずれを修正し、分解能が最適化された超音波画像を取得することができる。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２０１６−０５５１６２号公報
【特許文献２】特開２００８−２６４５３１号公報
【非特許文献】
【０００７】
【非特許文献１】Kitai, T. et al., "Photoacoustic mammography: initial clinicalresults," Breast Cancer, 21, pp.146-153 (2014).
【非特許文献２】Li C. et al. "Photoacoustic tomography and sensing in biomedicine," Physics in Medicine and Biology, 54 (19), pp.59-97 (2009).
【非特許文献３】Candes, E.J et al., "Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information," IEEE Transactions on Information Theory, pp. 489-509, (2006).
【非特許文献４】Zhang, Z. et al., "Total variation based gradient descent algorithm for sparse-view photoacoustic image reconstruction," Ultrasonics, 52, pp.1046-1055, (2012).
【非特許文献５】Xu M et al., "Time-domain reconstruction algorithms and numerical simulations for thermoacoustic tomography in various geometries, IEEE Trans on Biomedical Engineering, 50(9), pp.1086-1099, (2003).
【非特許文献６】Xu Y et al., "Exact frequency-domain reconstruction for thermoacoustic tomography I: Planar geometry," IEEE Medical Imaging, 21(7), pp.823-828, (2002).
【非特許文献７】Xu, M. et al., "Universal back-projection algorithm for photoacoustic computed tomography," Physical Review, 71, 016706, (2005).
【非特許文献８】Guo, Z. et al., "Compressed sensing in photoacoustic tomography in vivo," Journal of Biomedical Optics, 15(2), 021311, (2010).
【非特許文献９】Treeby, BE. et al., "Automatic sound speed selection in photoacoustic image reconstruction using an autofocus approach," Journal of Biomedical Optics, 16(9), 090501, (2011).
【非特許文献１０】Li, C. et al., "An Efficient Augmented Lagrangian Method with Applications to Total Variation Minimization," Computational Optimization and Applications, 56(3), pp.507-530, (2013).
【非特許文献１１】Treeby, BE. et al., "K-wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields," Journal of Biomedical Optics, 15(2), 021314, (2010).
【非特許文献１２】Kruger, RA. et al., "Dedicated 3D photoacoustic breast imaging," Medical Physics, 40(11), 113301, pp.1-8 (2013).
【非特許文献１３】Mast, TD. et al., "Simulation of ultrasonic pulse propagation through the abdominal wall," Journal of Acoustic Society of America, 102(2), pp.1177-1190, (1997).
【非特許文献１４】三村和史，「圧縮センシング−疎情報の再構成とそのアルゴリズム」，数理解析研究所講究録，第１８０３巻，２０１２年２６−５６，（２０１２）．
【非特許文献１５】R. Baraniuk et al., “Random projections of smooth manifolds,” 11 Foundations of Computational Mathematics, Vol. 9, No. 1, pp. 51-77, 2009.
【非特許文献１６】Y. C. Eldar et al., "Robust recovery of signals from a union of subspaces,” IEEE Transactions on Information Theory, Vol. 55, No. 11, pp. 5302-5316, 2009.
【非特許文献１７】R. Baraniuk, V. et al., "Model-based compressive sensing,” Preprint, 2008
【非特許文献１８】A. W. Fitzgibbon, "Robust Registration of 2D and 3D Point Sets," Proceedings of the British Machine Vision Conference, 2001.
【非特許文献１９】Rueckert, D. et al., "Nonrigid registration using free-form deformations: application to breast MR images," IEEE Transactions on Medical Imaging, 8, pp.712-721 (1999).
【非特許文献２０】R. Bise et al., " Vascular Registration in Photoacoustic Imaging by Low-Rank Alignment via Foreground, Background and Complement Decomposition," Proceedings of International Conference on Medical Image Computing and Computer-Assisted Intervention, pp.326-334, (2016).
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、特許文献１又は２の装置を用いて、検査対象物を画像化する場合において、音速は既知である必要があるが、音速は検査対象物の組織又は温度でゆらぎがあり、想定された音速とは異なり、画像が劣化し、画像の精度がいまだ低いという問題点があった。
【０００９】
本発明の目的は以上の問題点を解決し、従来技術に比較して高精度で対象物を画像化できる光超音波画像化装置及び方法、光超音波画像化装置の制御プログラム並びに記録媒体を提供することにある。
【課題を解決するための手段】
【００１０】
第１の発明にかかる光超音波画像化装置は、光超音波画像化法により、対象物に対して光を照射して前記対象物から発生する超音波を複数のセンサで複数の検出信号として検出し、音圧空間分布を示す複数の検出信号からなる音圧データに基づいて前記対象物を画像化する画像化処理を行う制御手段を備えた光超音波画像化装置において、
前記制御手段は、所定の初期音速データに基づいて、音速データのプロジェクション行列に初期の音圧ベクトルを乗算したときの乗算結果を音圧データとするときに、前記初期の音圧ベクトルに関する目的関数であって、前記対象物の画像データのスパース性、もしくは、前記画像データに対して所定の変換を行ったときの変換の基底のデータのスパース性を表す目的関数の値が所定のしきい値以下になるように、音速が異なる複数の領域の音速を含む音速ベクトルを変化させて、音圧ベクトルと音速ベクトルを同時に推定することで、前記音速ベクトルに基づく画像化された画像データを計算することを特徴とする。
【００１１】
前記光超音波画像化装置において、前記制御手段は、前記目的関数の値が実質的に最小になるように、前記音速ベクトルに基づく画像化された画像データを計算することを特徴とする。
【００１２】
また、前記光超音波画像化装置において、前記目的関数は、前記対象物の画像データ、もしくは前記画像データに対して所定の変換を行ったときの変換の基底のデータのＬ１ノルム、Ｌ２ノルム、又は全変動、もしくはこれらの組み合わせであることを特徴とする。
【００１３】
さらに、前記光超音波画像化装置において、前記制御手段は、音圧ベクトルと音速ベクトルを同時に推定する途中の前記対象物の画像データ、もしくは、前記画像化後の画像データを表示部に出力することで可視化することを特徴とする。
【００１４】
またさらに、前記光超音波画像化装置において、前記制御手段は、前記光超音波画像化法の空間解像度と時間解像度のうちの少なくとも一方を昇順で変化させて、前記光超音波画像化法により画像化処理をそれぞれ実行し、各画像化処理で得られた音速ベクトルを次の初期音速ベクトルとして用いて画像化処理を行うことを特徴とする。
【００１５】
またさらに、前記光超音波画像化装置において、
前記対象物は互いに音速が異なる第１及び第２の領域を含み、
前記制御手段は、
（１）所定の初期音速ベクトルを用いて、所定の画像化法を用いて前記対象物を画像化して画像データを計算し、
（２）前記計算された画像データに対して所定の画像処理を行うことで、前記第２の領域を推定し、
（３）所定の位置合わせ法により、前記推定した第２の領域を、前記対象物の各領域の音速初期値を地図で示す所定の音速初期値アトラスに位置合わせすることで初期音速マップを計算し、
（４）前記位置合わせされた初期音速マップを用いて、請求項１〜５のうちのいずれか１つに記載の画像化処理を行うことを特徴とする。
【００１６】
第２の発明に係る光超音波画像化方法は、光超音波画像化法により、対象物に対して光を照射して前記対象物から発生する超音波を複数のセンサで複数の検出信号として検出し、音圧空間分布を示す複数の検出信号からなる音圧データに基づいて前記対象物を画像化する画像化処理を行う制御手段を備えた光超音波画像化装置のための光超音波画像化方法において、
前記制御手段が、所定の初期音速データに基づいて、音速データのプロジェクション行列に初期の音圧ベクトルを乗算したときの乗算結果を音圧データとするときに、前記初期の音圧ベクトルに関する目的関数であって、前記対象物の画像データのスパース性、もしくは、前記画像データに対して所定の変換を行ったときの変換の基底のデータのスパース性を表す目的関数の値が所定のしきい値以下になるように、音速が異なる複数の領域の音速を含む音速ベクトルを変化させて、音圧ベクトルと音速ベクトルを同時に推定することで、前記音速ベクトルに基づく画像化された画像データを計算するステップを含むことを特徴とする。
【００１７】
前記光超音波画像化方法において、前記制御手段が、前記目的関数の値が実質的に最小になるように、前記音速ベクトルに基づく画像化された画像データを計算するステップを含むことを特徴とする。
【００１８】
また、前記光超音波画像化方法において、前記目的関数は、前記対象物の画像データ、もしくは前記画像データに対して所定の変換を行ったときの変換の基底のデータのＬ１ノルム、Ｌ２ノルム、又は全変動、もしくはこれらの組み合わせであることを特徴とする。
【００１９】
さらに、前記光超音波画像化方法において、前記制御手段が、音圧ベクトルと音速ベクトルを同時に推定する途中の前記対象物の画像データ、もしくは、前記画像化後の画像データを表示部に出力することで可視化するステップを含むことを特徴とする。
【００２０】
またさらに、前記光超音波画像化方法において、前記制御手段が、前記光超音波画像化法の空間解像度と時間解像度のうちの少なくとも一方を昇順で変化させて、前記光超音波画像化法により画像化処理をそれぞれ実行し、各画像化処理で得られた音速ベクトルを次の初期音速ベクトルとして用いて画像化処理を行うことを特徴とする。
【００２１】
またさらに、前記光超音波画像化方法において、
前記対象物は互いに音速が異なる第１及び第２の領域を含み、
前記制御手段が、
（１）所定の初期音速ベクトルを用いて、所定の画像化法を用いて前記対象物を画像化して画像データを計算し、
（２）前記計算された画像データに対して所定の画像処理を行うことで、前記第２の領域を推定し、
（３）所定の位置合わせ法により、前記推定した第２の領域を、前記対象物の各領域の音速初期値を地図で示す所定の音速初期値アトラスに位置合わせすることで初期音速マップを計算し、
（４）前記位置合わせされた初期音速マップを用いて、請求項１〜５のうちのいずれか１つに記載の画像化処理を行うステップを含むことを特徴とする。
【００２２】
第３の発明に係る、コンピュータにより実行可能な光超音波画像化装置の制御プログラムは、前記光超音波画像化方法の各ステップを含むことを特徴とする。
【００２３】
第４の発明に係る、コンピュータにより読み取り可能な記録媒体は、前記光超音波画像化装置の制御プログラムを格納することを特徴とする。
【発明の効果】
【００２４】
従って、本発明に係る光超音波画像化装置等によれば、従来技術に比較して高精度で検査対象物を画像化できる。
【図面の簡単な説明】
【００２５】
【図１】本発明の一実施形態に係る光超音波画像化システムの構成例を示す模式断面図及びブロック図である。
【図２】光超音波画像化処理の結果画像であって、（ａ）は真値画像、（ｂ）はバックプロジェクション法による画像、（ｃ）は不正確な音速を用いた圧縮センシング法（以下、ＣＳ法という）による画像、及び（ｄ）は正確な音速を用いたＣＳ法による画像を示す図である。
【図３】一定の音速ｃの場合における時刻ｔ１での超音波伝搬を示す平面図である。
【図４】実施形態１に係る光超音波画像化処理を示すフローチャートである。
【図５】真値画像と複数のセンサ２２の位置との関係を示す正面図である。
【図６】従来例であるバックプロジェクション法及びＣＳ法、並びに実施形態１に係る画像化法による音速に対するＰＳＮＲ（Peak Signal to Noise Ratio）特性を示すグラフである。
【図７】真値画像及び実施形態１に係る画像化法による画像におけるデータセットの音速に対する推定された音速を示すグラフである。
【図８】種々の音速における画像化画像であって、（ａ）はバックプロジェクション法による画像、（ｂ）はＣＳ法による画像、及び（ｃ）は実施形態１に係る画像化法による画像を示す図である。
【図９】（ａ）の音速の組み合わせを用いた画像化結果画像であって、（ｂ）はバックプロジェクション法による画像、（ｃ）はＣＳ法による画像、（ｄ）は実施形態１に係る画像化による画像、（ｅ）は（ｃ）の画像の拡大画像、及び（ｆ）は（ｄ）の画像の拡大画像を示す図である。
【図１０】本発明の実施形態２に係る高精度化光超音波画像化処理の一例を示すフローチャートである。
【図１１】本発明の実施形態３に係る高精度化光超音波画像化処理の一例を示すフローチャートである。
【図１２】本発明の実施形態４に係る音速初期値アトラスを用いた画像化処理を示すフローチャートである。
【図１３】図１２の画像化処理で用いる音速初期値アトラスを示す模式図である。
【図１４Ａ】図１２のステップＳ３１における処理を示す模式図である。
【図１４Ｂ】図１２のステップＳ３２における処理を示す模式図である。
【図１４Ｃ】図１２のステップＳ３３の位置合わせ前における処理を示す模式図である。
【図１４Ｄ】図１２のステップＳ３３の位置合わせ後における処理を示す模式図である。
【図１４Ｅ】図１２のステップＳ３４における処理を示す模式図である。
【図１４Ｆ】図１２のステップＳ３５における処理を示す模式図である。
【発明を実施するための形態】
【００２６】
以下、比較例及び本発明に係る実施形態について図面を参照して説明する。なお、以下の各実施形態において、同様の構成要素については同一の符号を付している。
【００２７】
本発明の実施形態の目的は、限定された数のサンプリングデータから画像化された検査対象物（オブジェクト）の３次元構造を回復するときに、同時に音速を推定しながら画像化するための新しい圧縮センシング法（以下、ＣＳ法という）による画像化法を提案することにある。ＣＳ法は、多くの信号が適切な基準又は辞書でわずか数の非ゼロ係数で表されるという基本的な事実に基づいており、限られた数の観測からでも再構成精度を向上させるのに有効である。ここで、非線形最適化は、そのような信号を少数の測定から回復することを可能にする（例えば、非特許文献３参照）。ＣＳ法による画像化は、様々な医用画像化方式に適用されてきたが、ＭＲＩ、ＣＴ、及びＰＡなどのような様々な方法がある（例えば、非特許文献４参照）。現在のＰＡイメージングのＣＳ法は、誤った音速の場合には有効に機能しないという問題点があった。
【００２８】
図２は光超音波画像化処理の結果画像であって、図２（ａ）は真値画像、図２（ｂ）はバックプロジェクション法による画像、図２（ｃ）は不正確な音速を用いたＣＳ法による画像、及び図２（ｄ）は正確な音速を用いたＣＳ法による画像を示す図である。
【００２９】
図２（ｄ）から明らかなように、正しい音速が与えられたとき、ＣＳ法は真値画像（図２（ａ））と同様の良好な結果を得ることができる。対照的に、与えられた音速が正しくない場合、バックプロジェクション法とＣＳ法の両方では、あまり良い画像を得ることができない。（図２（ｂ）及び図２（ｃ）参照）。
【００３０】
以上の問題点を解決するために、本実施形態では、スパース性を有する超音波信号から再構成問題を同時に解決し、音速を推定することができるＰＡイメージングのための新しいＣＳ法による画像化処理を提案する。様々なＣＳ法の最適化技術を当該実施形態内で使用することができる。特に、本発明者らは、拡張されたラグランジュ乗数法を用いて全変動（Total Variation）を最小化するＣＳ法を最適化する方法を使用した。発明者による実験の結果は、この提案された方法の有効性を実証することができ、音速を評価して画像化画質を大幅に改善することができたものであり、本実施形態は、音速を同時に推定しながら画像再構成の問題に取り組むための最初の試みである。なお、従来例では、均一な音速を想定してもこの問題は解決されていなかったものである。
【００３１】
図１は本発明の一実施形態に係る光超音波画像化システムの構成例を示す模式断面図及びブロック図である。
【００３２】
図１において、本実施形態に係る光超音波画像化システムは、光超音波光源及び検出装置２０と、光超音波画像化装置１０とを備えて構成される。
【００３３】
光超音波光源及び検出装置２０において、例えば半球面形状のセンサ支持部２１に互いに所定の間隔を有して、超音波を検出する複数のセンサ２２が設けられる。上記半球面の内部に水２４を浸して、検査対象物１を載置するメッシュ板２３を支持している。なお、検査対象物１が水２４内に沈められているものとする。センサ支持部２１の最底部には、ＰＡイメージングのための所定波長のレーザ光を検査対象物１に対して放射するためのレーザ光源２５が設けられる。ここで、レーザ光源２５は後述するコントローラ１１により制御され、ＰＡイメージングによりレーザ光に応答して検査対象物１から放射される超音波を複数のセンサ２２により検出する。複数のセンサ２２により検出された超音波は電気信号に変化された後、後述する信号受信回路１２に送られる。
【００３４】
検査対象物１は撮像対象である。具体例としては、乳房等の生体や、装置の調整などにおいては生体の音響特性と光学特性を模擬したファントムが挙げられる。音響特性とは具体的には音響波の伝播速度および減衰率であり、光学特性とは具体的には光の吸収係数および散乱係数である。なお、メッシュ板２３は検査対象物１を支持し、光、超音波に対する透過特性を持つことが好ましい。
【００３５】
光超音波画像化装置１０は、装置の動作を制御するコントローラ１１と、信号受信回路１２と、データメモリ１３と、操作部１４と、表示部１５と、画像処理部１６とを備えて構成される。信号受信回路１２は、複数のセンサ２２で変換された複数の電気信号を受信して所定のサンプリング周波数でＡ／Ｄ変換してデータメモリ１３に格納する。操作部１４は光超音波画像化装置１０の操作者がデータ又は指示を入力するためのマウス、キーボードを含み構成される。表示部１５は、入力されたデータ又は指示を表示するとともに、光超音波画像化処理の途中又は処理後の結果画像を可視化して表示する。画像処理部１６はコントローラ１１により制御され、詳細後述する光超音波画像化処理を実行する。
【００３６】
なお、図１の実施形態では、半球形状のセンサ支持部２１を用いているが、本発明はこれに限らず、複数のセンサ２２を所定の領域に指向軸が集まり、所定の検出領域を形成できる配置であればよい。
【００３７】
図１の光超音波画像化システムでは、光超音波画像化法を用いて、光音響効果を利用して例えば生体の血管を非侵襲的に視覚化することができ、癌及び他の多くの疾患の早期臨床診断に有用である。検査対象物１である画像化すべき組織はレーザエネルギーを吸収し、そのエネルギーを熱弾性膨張に変換し、超音波を放射する。光超音波を検出することによって検査対象物１の物体の３次元構造を再構成することができる。ほとんどのＰＡイメージングでは、音速は正しく与えられていると想定されているが、音速は伝播媒体及び組織の温度に依存して変化することが知られている。従って、生体でこの値を正確に知ることは困難であるが、本実施形態では、再構成問題を同時に解決し、音速を推定するＰＡのための新しい圧縮センシング法を提案する。
【００３８】
例えば、本実施形態に関連する文献において、逆球面ラドン変換に基づくＰＡ再構成法は、時間領域（例えば、非特許文献５参照）と周波数領域（例えば、非特許文献６参照）の両方で提案されている。非特許文献７において提案されたフィルタリングされたバックプロジェクション法は、その利便性のために最も普及している。これらの方法の１つの問題点は、サンプリング数が不十分であるときに生成される深刻なアーチファクトである。スパース性を有するサンプリングデータを用いて画質を改善するために、ＣＳ法による再構成にＣＳ理論を適用した。非特許文献８では、ＰＡ再構成においてＣＳ理論を用いる可能性を実証した。また、非特許文献４では、勾配降下による目的関数としての全変動（ＴＶ）を反復的に最小化する方法を提案した。これらの方法は、音速が正しく与えられていることを前提としている。音速を選択するために、非特許文献９は、様々な候補音速を用いて画像を繰り返し再構築する自動焦点アプローチを提案したが、音速が空間的に均一であると仮定される場合、
次に、鮮明度のメトリックを最大にする最良の音速を選択する。しかし、音速が異なる複数の領域がある場合には、組み合わせの数は指数関数的に大きくなるという問題点があった。本発明者らが知るうる限り、この挑戦的な問題に最初に取り組んでおり、本実施形態に係る方法は最適化において音速のセットを反復的に更新することを特徴としている。
【００３９】
まず、音速推定によるＰＡイメージングのための圧縮センシングについて以下に説明する。ここで、ＰＡ伝播の離散化された定式化について簡単に説明する初期の音圧がどのようにして超音波として検出器に伝播するかを示す。ＰＡイメージングでは、再構成される目標画像マップは、音響初期音圧の空間分布であると考えることができる。
【００４０】
図３は一定の音速ｃの場合における時刻ｔ１での超音波伝搬を示す平面図である。図３において、光吸収物体Ａ及びＣからの超音波信号は、時刻ｔ１で別のセンサ２２に到着したばかりであり、音速ｃと時刻ｔ１との積から距離Ｒを算出することができる。時刻ｔ１でのセンサ２２の観測信号は丁度到着した超音波信号の合計であると考えることができる。非特許文献４によれば、この前方伝搬問題は次式で表される。
【００４１】
ｂ＝Ａ（ｃ）ｕ      （１）
【００４２】
ここで、ＴＫ×１ベクトルｂは、複数のセンサ２２からなるセンサアレイによって測定された電気信号を示し、Ｋはセンサ２２の数であり、Ｔはタイムスロットの数である。ベクトルｕは初期音圧の空間分布を表し、ここで、ｕ（Ｎ＝Ｎｘ×Ｎｙ×Ｎｚ）の大きさは視野に依存する。Ｎｘはｘ方向のセンサ２２の数であり、Ｎｙはｙ方向のセンサ２２の数であり、Ｎｚはｚ方向のセンサ２２の数である。
【００４３】
ここで、ＴＫ×Ｎ行列Ａ（ｃ）をプロジェクション（射影）行列と定義し、音圧がどのように各位置と時間に伝播するかを記述する。Ａ（ｃ）の実施例は補足説明Ａ２に記載されている。
【００４４】
再構成問題は、ｂ＝Ａ（ｃ）ｕを満たす初期音圧ｕの分布を解くことを含む。サンプリングデータが不十分な場合は、その逆の問題は悪い考えである。これは、再構成された画像に縞模様のアーチファクトをもたらすからである。この問題に対処するために、ＣＳ法において、再構成問題の最適化における目標画像の特性を考慮に入れる。ここでは、例えば、フーリエドメイン（例えば、非特許文献３参照）、小さな全変動（例えば、非特許文献４参照）のような、より基本的なドメインにおけるスパース性を仮定して、ＣＳ法における再構成問題を以下のように定式化することができる。
【００４５】
【数１】
2018161318.tif 000003
（２）
【００４６】
ここで、「ｓ．ｔ．」は「ｓｏ  ｔｈａｔ」の略で制約条件を示す。ｆは目標画像の特性（スパース性などの特徴）を記述する目的関数である。典型的には、Ｌ１ノルム｜Ψｕ｜１又は全変動ＴＶ（ｕ）は目的関数ｆのために使用される。ここで、Ψはスパース変換行列である。音速マップｃが正しく与えられていれば、プロジェクション行列Ａを計算して固定することができる。上記式（２）では、ｂ＝Ａ（ｃ）ｕの条件のもとで、音圧の空間分布ｕと音速ｃとを変化させて目的関数ｆ（ｕ）の値が小さくなるような、音圧の空間分布ｕと音速ｃを求める。従って、ＣＳ法（例えば、非特許文献８参照）を用いて画像を再構成することができる。しかし、この問題を解くことは、未知数の数が多く、プロジェクション行列Ａ（ｃ）は通常の音速ｃでは非線形であるために、ｕとｃを同時に推定することは困難である。
【００４７】
次いで、音速ｃを推定するための反復線形化について以下に説明する。
【００４８】
上述のように、式（２）の最適化問題は直接扱いにくい。ここで、大きな困難性の１つは、２Ｎ個の多数の未知パラメータである。この問題に対処するために、発明者らは追加的だが現実的な仮定を導入することによってそれを緩和することができる。撮像対象では、本発明者らは異なる音速を持つ領域はほんのわずかであり、これらの領域の位置は大体知られている。すなわち、例えば検査対象物１の組織の位置は患者によって大きく異なるわけではないため、小さな位置の差異は再構成に深刻な影響を与えない。ここで、音速ベクトルｃを次式で表す。
【００４９】
ｃ＝［ｃ１，ｃ２，…，ｃＬ］      （３）
【００５０】
ここで、ｃｌはｌ番目の領域の音速であり、Ｌは速度の異なる領域の数である。以下、本発明者らは再構成問題におけるＬ個の速度（音速）を推定する。もう一つの難点は非線形性であり、音響伝搬関数Ａ（ｃ）は音速ｃに対して線形ではない。音速ｃの変化が小さい場合には、音速ｃの現在の推定値について線形化することによって、Δｃ＝［Δｃ１｜…｜ｃＬ］に対して、プロジェクション行列Ａ（ｃ）を近似することができる。ここで、Ａｉ，ｊ，ｋ（ｃ＋Δｃ）の要素は、次式のように近似することができる。
【００５１】
【数２】
2018161318.tif 000004
（４）
【００５２】
これは、未知数ｕ，Δｃにおける最適化問題につながり、目的関数は次式のように記述できる。
【００５３】
【数３】
2018161318.tif 000005
（５）
【００５４】
線形化は局所的にしか保持されないので、式（５）からの解（ｃ＋Δｃ）が式（１）を正確に解くことを期待すべきではない。式（１）の局所最小値を見つけるために、本発明者らは音速ベクトルｃの現在の推定値を繰り返し線形化し、式（５）の一連の最適化を解くことができる。後述する実験結果でわかるように、初期の音速の速度差が大きすぎない限り、この反復は、初期音圧ｕ及び音速ｃを推定することができる。完全な最適化手順は、図４の光超音波画像化処理として要約される。なお、式（５）では、目的関数の値が最小値であるときの結果データ（ｕ，Δｃ）を求めているが、本発明はこれに限らず、実質的に最小値、もしくは当該最小値に近傍である、所定のしきい値以下の値である上記結果データを求めてもよい。
【００５５】
図４は実施形態１に係る光超音波画像化処理を示すフローチャートである。当該光超音波画像化処理は、線形化された最適化のシーケンスを表す。
【００５６】
図４のステップＳ１において、まず、初期音速ベクトルをｃとする。次いで、ステップＳ２において、ベクトル化された検出信号データをｂとする。
【００５７】
【数４】
2018161318.tif 000006
【００５８】
ここで、ｂは、Ｔ個のタイムスロットと、Ｋ個のセンサ数とで形成される空間におけるベクトル化された検出信号データである。
【００５９】
ステップＳ３では、現在の音速ベクトルｃに基づいて、次式を用いてプロジェクション行列Ａ（ｃ）及びヤコビアン行列を計算する。
【００６０】
【数５】
2018161318.tif 000007
（６）
【００６１】
次いで、ステップＳ４では、次式を用いて、線形最適化の計算を行うことで結果データ（ｕ＊，Δｃ＊）を計算する。ここで、Δｃは音速変化量ベクトルであり、Δｃ＊は線形最適化後の音速変化量ベクトルである。
【００６２】
【数６】
2018161318.tif 000008
（７）
【００６３】
式（７）では、具体的には、補足説明Ａ１における式（１０）を用いる。
【００６４】
ステップＳ５では、ｃ＋Δｃ＊を計算して次回の音速ベクトルｃとおく。ステップＳ６では、現在の音圧データｕ＊を用いて信号データｂ＝Ａ（ｃ）ｕ＊の計算を行って信号データｂを計算して表示部１５に出力して中間画像化画像を表示部１５に表示する。次いで、ステップＳ７では、収束条件（Δｃ＊＜Δｃｔｈ）が満たすか否かを判断し、ＹＥＳのときはステップＳ８進む一方、ＮＯのときはステップＳ３に戻る。ここで、Δｃｔｈは、音速変化量ベクトルのしきい値である。さらに、ステップＳ８では、結果データ（ｕ＊，Δｃ＊）及び信号データｂを表示部１５に出力して画像化画像を表示して当該処理を終了する。
【００６５】
次いで、各繰り返しの最適化について以下に説明する。
【００６６】
提案された最適化方法では、様々なＣＳ法を使用することができる。本実施形態では、全変動（ＴＶ：Total Variation）の目的関数を用いる。全変動の目的関数を最小化するＴＶ最小化法は、現代の生物医学的画像再構成において広く使用されており、ＰＡイメージングについても適用されている（例えば、非特許文献４参照）。全変動の目的関数では、基礎となる信号の勾配が疎である（スパース性を有する）ことが前提である。これは、ｌ１の最小化と比較してより正確にエッジ又は境界を保存することによって、再構成された画像をより鮮明にするためである。全変動の目的関数を用いた最小化問題を解決するために、本発明者らは、拡張ラグランジュ乗数（ＡＬＭ）（例えば、非特許文献１０参照）を含む方法（拡張ラグランジュ乗数法）を用いる。具体的には、ＣＳ再構成のための離散ＴＶモデルは次式で表される。
【００６７】
【数７】
2018161318.tif 000009
（８）
【００６８】
ここで
【数８】
2018161318.tif 000010
の和は全変動の目的関数における期間を示し、Ｄｉｕはボクセルｉにおける音圧空間分布ｕの離散的勾配である。式（８）の最適化問題は、次式の拡張ラグランジュ関数を最小化することによって解決することができる。
【００６９】
【数９】
2018161318.tif 000011
（９）
【００７０】
ここで、ＤはＮ×Ｎ行列であり、ｉ行目がＤｉであり、Σ及びΔは各制約条件に対するラグランジュ乗数行列であり、βとμは正のスカラー値である。なお、最適化の詳細については、補足説明Ａ１を参照してください。
【００７１】
次いで、実験的評価について以下に説明する。
【００７２】
本発明者らは、ＰＡ再構成に対する本方法の有効性を検証するためにシミュレーションを行った。ＭＡＴＬＡＢのｋ−ｗａｖｅツールボックス（例えば、非特許文献１１参照）を使用して光音響伝播をシミュレートした。なお、このツールボックスは光音響シミュレーションに広く使用されている。ここで、シミュレーションは２次元で行った。
【００７３】
図５は真値画像と複数のセンサ２２の位置との関係を示す正面図である。図５において、右端のグレースケールは最初の音響音圧の強さを示す。黒い点はセンサ２２の位置を示す。ここで、センサ２２の数は７０である。再構成画像のサイズは２１６×２１６ピクセルであり、各画素のサイズは０．０４６３［ｍｍ］であった。
【００７４】
まず、音速が一定の場合の手法を評価した。１４６０〜１５５０［ｍ／ｓ］の種々の音速を使用してテストデータセットを生成し、音速の増分は１０［ｍ／ｓ］であり、生体組織の音速はこの範囲内にあることが多いからである。その結果を、時間逆転法（例えば、非特許文献１１参照）を用いたバックプロジェクション法、及びＡＬＭを用いた全変動目的関数を最小化するＣＳ法（例えば、非特許文献１０参照）で得られた結果と比較した。これらの方法は、各シミュレートされたデータセットに適用され、音速は１５００［ｍ／ｓ］とした。
初期の音速として１５００［ｍ／ｓ］を設定した再構成プロセスにおける速度を推定した。
【００７５】
図８は種々の音速における画像化画像であって、図８（ａ）はバックプロジェクション法による画像、図８（ｂ）はＣＳ法による画像、及び図８（ｃ）は実施形態１に係る画像化法による画像を示す図である。ここで、バックプロジェクション及びＣＳ法は、誤って与えられた音速のためにアーチファクトを被ることが分かる。対照的に、本実施形態に係る方法は１４８０〜１５４０［ｍ／ｓ］の音速範囲で画像を正常に再生できた。画質を定量的に評価するために、元のファントムを参照して再構成画像のピーク信号対雑音比（ＰＳＮＲ）を計算した。ここで、ＰＳＮＲの価値が大きいほど、画像の品質が良くなる。
【００７６】
図６は従来例であるバックプロジェクション法及びＣＳ法、並びに実施形態１に係る画像化法による音速に対するＰＳＮＲ（Peak Signal to Noise Ratio）特性を示すグラフである。図６から明らかなように、本実施形態に係る結果の画質は１４８０〜１５４０［ｍ／ｓ］の中で最高である。
【００７７】
図７は真値画像及び実施形態１に係る画像化法による画像におけるデータセットの音速に対する推定された音速を示すグラフである。すなわち、図７は本実施形態の音速推定を評価したものである。本実施形態に係る方法は１４８０から１５４０［ｍ／ｓ］までの正確な音速を正常に推定した。次に、速度の異なる複数の領域が存在する場合の方法を評価した。
【００７８】
図９（ａ）の音速の組み合わせを用いた画像化結果画像であって、図９（ｂ）はバックプロジェクション法による画像、図９（ｃ）はＣＳ法による画像、図９（ｄ）は実施形態１に係る画像化による画像、図９（ｅ）は図９（ｃ）の画像の拡大画像、及び図９（ｆ）は図９（ｄ）の画像の拡大画像を示す図である。ここで、図９（ａ）に示す音速の組み合わせを用いてテストデータセットを生成し、音速の範囲は１４７０〜１５４０［ｍ／ｓ］であった。図９から明らかなように、音速が１５４０〜１５１０［ｍ／ｓ］でバックプロジェクション法（図９（ｂ））及びＣＳ法（図９（ｃ））の結果においてアーチファクトが発生している。対照的に、本実施形態に係る方法は、画像を正常に回復し、すべての組み合わせについて音速を推定した。ここで、平均ＰＳＮＲは以下の通りであった。
（１）バックプロジェクション法：２１．７５；
（２）ＣＳ法：２３．４６；
（３）実施形態１に係る方法：２８．４５、
なお、平均音推定誤差は２．１７［ｍ／ｓ］（１．４％）であった。
【００７９】
以上説明したように本実施形態に係る方法によれば、ＰＡイメージングのための新しい再構成方法を提案し、スパース性を有する超音波信号から再構成問題を、ＣＳ法を使用して音速を同時に推定することで解くことができた。実験の結果は、実施形態に係る方法の有効性を実証し、音速と画質の大幅な改善に成功した。
【００８０】
（補足説明Ａ１）ＡＬＭの最適化．
式（９）の最適化問題は、次式（１０）のように反復的に最小化される。
【００８１】
【数１０】
2018161318.tif 000012
【数１１】
2018161318.tif 000013
【数１２】
2018161318.tif 000014
（１０）
【００８２】
ここで、「ｓｇｎ」は符号関数であり、ｋは反復のインデックスであり、Ａｍｉｊｏのような非モノトーンラインサーチを用いて決定される。
【００８３】
（補足説明Ａ２）プロジェクション行列．
非特許文献１３によれば、プロジェクション行列Ａ（ｃ）に次式の減衰正弦波を用いる。
【００８４】
【数１３】
2018161318.tif 000015
（１１）
【００８５】
ここで、ω０は超音波信号の周波数であり、ｉは位置指数であり、ｔｊはｊ番目のタイムスロットの時間、ｄｌは、ｋ番目のセンサ２２と音源との間の経路に音速ｃｌを有するエリア上の距離である。
【００８６】
なお、プロジェクション行列Ａ（ｃ）に関して、以下を利用してもよい。
平均を
【数１４】
2018161318.tif 000016
、分散σとするＧａｕｓｓ分布関数等のタイムスロットｔｊが
【数１５】
2018161318.tif 000017
のとき最大となり、これから離れるにつれて指数的に減衰する形状の関数。
【００８７】
（実施形態１の変形例）
以上の実施形態１において用いる目的関数としては以下の種々のもの、もしくは、前記目的関数又はこれらの組み合わせを用いることができる。
（１）再構成画像を表す音圧空間分布ｕを所定の変換関数により変換した後の基底ドメイン（ベースドメイン）でのスパース性を用いる場合：
Ｌ１ノルム｜Ψｕ｜１を利用する。ここで、例えば以下の変換関数の基底を用いることができる。
（ａ）フーリエ変換：Ψ  は逆フーリエ変換の基底である。
（ｂ）ウエーブレット（Ｗａｖｅｌｅｔ）変換：Ψは逆Ｗａｖｅｌｅｔ変換の基底である。
（ｃ）離散コサイン変換
（ｄ）Ｄｉｃｔｉｏｎａｒｙ  ｏｆ  ｏｒｔｈｏｇｏｎａｌ  ｂａｓｅｓ（例えば、非特許文献１６参照）
（２）再構成画像を表す音圧空間分布ｕ自体のスパース性を用いる場合：
Ｌ１ノルム｜ｕ｜１を利用する。ここで、例えば以下の関数を用いる。
（ａ）ＴＶ（ｔｏｔａｌ  ｖａｒｉａｔｉｏｎ）：全変動；
（ｂ）Ｌ１ノルム：｜ＴＶ（ｕ）｜１：局所的な輝度の変化がスパースである；
（ｃ）Ｌ２ノルム：｜ＴＶ（ｕ）｜２：局所的な輝度の変化量の総和が少ない；
（ｄ）ｌｏｗ  ｄｉｍｅｎｓｉｏｎａｌ  ｓｍｏｏｔｈ  ｍａｎｉｆｏｌｄｓ（例えば、非特許文献１７参照）；
（ｅ）ｕｎｉｏｎ  ｏｆ  ｓｕｂ−ｓｐａｃｅｓ（例えば、非特許文献１８参照）；
（ｆ）ｂｌｏｃｋ−ｓｐａｒｓｉｔｙ（例えば、非特許文献１９参照）。
（３）上記（１）及び（２）の組み合わせ（例えば線形結合など）を用いてもよい。
【００８８】
実施形態２．
図１０は本発明の実施形態２に係る高精度化光超音波画像化処理の一例を示すフローチャートである。実施形態２では、センサ２２の数Ｎを例えばＮ＝２０×２０×２０、Ｎ＝４０×４０×４０、Ｎ＝１００×１００×１００と少ない数から多い数に変化させることで空間解像度を昇順で変化させて、実施形態１に係る光超音波画像化処理（Ｓ１００）をそれぞれを実行して各処理で得られた音速ベクトルｃ＊を次の初期音速ベクトルｃとして用いて画像化処理を行うことを特徴としている。
【００８９】
図１０において、ステップＳ１１において、センサ数Ｎ＝２０×２０×２０として、ステップＳ１２でステップＳ１００の光超音波画像化処理を実行して、音速ベクトルｃ＊を求める。次いで、ステップＳ１３において、上記求めた音速ベクトルｃ＊を次の画像化処理のための初期音速ベクトルとして用いる。
【００９０】
次いで、ステップＳ１４において、センサ数Ｎ＝４０×４０×４０として、ステップＳ１５でステップＳ１００の光超音波画像化処理を実行して、音速ベクトルｃ＊を求める。次いで、ステップＳ１６において、上記求めた音速ベクトルｃ＊を次の画像化処理のための初期音速ベクトルとして用いる。
【００９１】
次いで、ステップＳ１５において、センサ数Ｎ＝１００×１００×１００として、ステップＳ１６でステップＳ１００の光超音波画像化処理を実行して、音速ベクトルｃ＊を含む結果データ（ｕ，ｃ＊）を求めて表示部１５に出力して可視化して表示する。
【００９２】
以上説明したように実施形態２によれば、空間解像度を昇順で変化させて、実施形態１に係る光超音波画像化処理（Ｓ１００）をそれぞれを実行して各処理で得られた音速ベクトルｃ＊を次の初期音速ベクトルｃとして用いて画像化処理を行うことで、画像化画像の精度を高めることができる。
【００９３】
実施形態３．
図１１は本発明の実施形態３に係る高精度化光超音波画像化処理の一例を示すフローチャートである。実施形態３では、所定の単位時間当たりのタイムスロットの数Ｔ（サンプリング周波数の逆数に対応する）を例えばＴ＝５０、Ｔ＝２００、Ｔ＝１０００と少ない数から多い数に変化させることで時間解像度を昇順で変化させて、実施形態１に係る光超音波画像化処理（Ｓ１００）をそれぞれを実行して各処理で得られた音速ベクトルｃ＊を次の初期音速ベクトルｃとして用いて画像化処理を行うことを特徴としている。
【００９４】
図１１において、ステップＳ２１において、タイムスロット数Ｔ＝５０として、ステップＳ２２でステップＳ１００の光超音波画像化処理を実行して、音速ベクトルｃ＊を求める。次いで、ステップＳ２３において、上記求めた音速ベクトルｃ＊を次の画像化処理のための初期音速ベクトルとして用いる。
【００９５】
次いで、ステップＳ２４において、タイムスロット数Ｔ＝２００として、ステップＳ２５でステップＳ１００の光超音波画像化処理を実行して、音速ベクトルｃ＊を求める。次いで、ステップＳ１６において、上記求めた音速ベクトルｃ＊を次の画像化処理のための初期音速ベクトルとして用いる。
【００９６】
次いで、ステップＳ２５において、タイムスロット数Ｔ＝１０００として、ステップＳ２６でステップＳ１００の光超音波画像化処理を実行して、音速ベクトルｃ＊を含む結果データ（ｕ，ｃ＊）を求めて表示部１５に出力して可視化して表示する。
【００９７】
以上説明したように実施形態３によれば、時間解像度を昇順で変化させて、実施形態１に係る光超音波画像化処理（Ｓ１００）をそれぞれを実行して各処理で得られた音速ベクトルｃ＊を次の初期音速ベクトルｃとして用いて画像化処理を行うことで、画像化画像の精度を高めることができる。
【００９８】
なお、実施形態２及び３を組み合わせてもよい。すなわち、空間解像度と時間解像度の少なくとも一方を昇順で変化させて、実施形態１に係る光超音波画像化処理（Ｓ１００）をそれぞれを実行して各処理で得られた音速ベクトルｃ＊を次の初期音速ベクトルｃとして用いて画像化処理を行ってもよい。
【００９９】
実施形態４．
図１２は本発明の実施形態４に係る音速初期値アトラスを用いた画像化処理を示すフローチャートである。実施形態４に係る画像化処理は、音速初期値アトラスを用いて画像化精度を高めることを特徴としている。ここで、音速初期値アトラスとは予め作成及び保持しており実施形態１において初期音速（標準値）を決めるために利用される２次元又は３次元の音速アトラスであって、標準音速エリアが示されている地図（マップ）をいう。
【０１００】
図１３は図１２の画像化処理で用いる音速初期値アトラスを示す模式図である。図１３から明らかなように、検査対象物１内の脂肪と筋肉の各領域で音速初期値が異なる。ここで、脂肪と筋肉は音速が互いに異なる一例の領域である。
【０１０１】
図１２のステップＳ３１において、まず、音速を例えば１５００ｍ／ｓで固定とし、公知のバックプロジェクション法（又はＣＳ法）を用いて画像データの再構成を行う（図１４Ａ）。ここで、音速が異なると画質には影響が出るが、大まかな組織の場所を知る上では、影響は少ないと考えられる。次いで、ステップＳ３２において、ステップＳ３１で得られた画像データに対して、公知のスムージング、膨張等の処理と、しきい値処理等の公知の画像処理を行って、血管分布から、血管を含む体組織領域を推定する（図１４Ｂ）。ステップＳ３３において、ステップＳ３２で、公知の位置合わせ法（例えば、非特許文献２０参照）を用いて推定した体組織領域と、音速初期値アトラスとの間で位置合わせを行う（図１４Ｃ）ことで初期音速マップ（図１４Ｄ）を計算する。具体的には、例えばステップＳ２２で求めた体組織領域（一点鎖線内）をアトラスの領域（点線内）に対して位置合わせを以下の方法で行う。
（１）例えばアフィン変換で、体組織領域を平行移動、スケール、回転等によりアトラス領域に位置合わせする（例えば、非特許文献１８参照）。
（２）例えば格子点をコントロール点とした非剛体位置合わせ（例えば、非特許文献１９参照）により、体組織領域をアトラス領域に位置合わせする。
【０１０２】
次いで、ステップＳ３４において、ステップＳ３３で求めた位置合わせ後の初期音速マップを用いて、異なる音速を有する領域及び初期音速マップを生成する。さらに、ステップＳ３５において、ステップＳ３４で生成した初期音速マップを初期音速として用いて光超音波画像化処理（Ｓ１００）を実行することで音速マップ及び再構築画像データを推定する。
【０１０３】
以上説明したように本実施形態によれば、音速初期値アトラスを用いて、実施形態１と組み合わせて位置合わせすることで、実施形態１に係る方法に比較して高精度で画像化できる。
【０１０４】
以上の各実施形態に係る画像化処理又は方法の各ステップをコンピュータで実行可能なプログラムとして構成することができる。当該プログラムは例えばＣＤ−ＲＯＭ又はＤＶＤ−ＲＯＭなどのコンピュータにより読取可能な記録媒体に格納することができる。
【産業上の利用可能性】
【０１０５】
以上詳述したように、本発明によれば、従来技術に比較して高精度で、例えば人体などの生体組織の対象物を画像化できる医用画像化装置等を提供することができる。
【符号の説明】
【０１０６】
１…検査対象物、
１０…光超音波画像化装置、
１１…コントローラ、
１２…信号受信回路、
１３…データメモリ、
１４…操作部、
１５…表示部、
１６…画像処理部、
２０…光超音波光源及び検出装置、
２１…センサ支持部、
２２…センサ、
２３…メッシュ板、
２４…水、
２５…レーザ光源。
【図１】
2018161318.tif 000018
【図２】
2018161318.tif 000019
【図３】
2018161318.tif 000020
【図４】
2018161318.tif 000021
【図５】
2018161318.tif 000022
【図６】
2018161318.tif 000023
【図７】
2018161318.tif 000024
【図８】
2018161318.tif 000025
【図９】
2018161318.tif 000026
【図１０】
2018161318.tif 000027
【図１１】
2018161318.tif 000028
【図１２】
2018161318.tif 000029
【図１３】
2018161318.tif 000030
【図１４Ａ】
2018161318.tif 000031
【図１４Ｂ】
2018161318.tif 000032
【図１４Ｃ】
2018161318.tif 000033
【図１４Ｄ】
2018161318.tif 000034
【図１４Ｅ】
2018161318.tif 000035
【図１４Ｆ】
2018161318.tif 000036