特許 5782314 生体計測装置および画像作成方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】5782314

(24)【登録日】2015年7月24日

(45)【発行日】2015年9月24日

(54)【発明の名称】生体計測装置および画像作成方法

(51)【国際特許分類】

   A61B 10/00 20060101AFI20150907BHJP

   G01N 21/17 20060101ALI20150907BHJP

   G01N 22/02 20060101ALI20150907BHJP

   A61B 6/03 20060101ALI20150907BHJP

   A61B 8/00 20060101ALI20150907BHJP

【ＦＩ】

   A61B10/00 EZDM

   G01N21/17 620

   G01N22/02

   A61B6/03 350U

   A61B8/00

【請求項の数】8

【全頁数】25

(21)【出願番号】特願2011-151086(P2011-151086)

(22)【出願日】2011年7月7日

(65)【公開番号】特開2013-19696(P2013-19696A)

(43)【公開日】2013年1月31日

【審査請求日】2014年7月2日

【国等の委託研究の成果に係る記載事項】（出願人による申告）独立行政法人科学技術振興機構、先端計測分析技術・機器開発事業、産業技術力強化法第１９条の適用を受ける特許出願

(73)【特許権者】

【識別番号】000236436

【氏名又は名称】浜松ホトニクス株式会社

(74)【代理人】

【識別番号】100088155

【弁理士】

【氏名又は名称】長谷川 芳樹

(74)【代理人】

【識別番号】100113435

【弁理士】

【氏名又は名称】黒木 義樹

(74)【代理人】

【識別番号】100124291

【弁理士】

【氏名又は名称】石田 悟

(74)【代理人】

【識別番号】100174399

【弁理士】

【氏名又は名称】寺澤 正太郎

(72)【発明者】

【氏名】工藤 博幸

(72)【発明者】

【氏名】齋藤 直哉

(72)【発明者】

【氏名】上田 之雄

(72)【発明者】

【氏名】佳元 健治

(72)【発明者】

【氏名】山下 豊

【審査官】 門田 宏

(56)【参考文献】

【文献】 国際公開第２０１０／１４３４２１（ＷＯ，Ａ１）

【文献】 特開２０１０−９４５００（ＪＰ，Ａ）

【文献】 特表２０１１−５０５９５４（ＪＰ，Ａ）

【文献】 特表２０１０−５１４４７８（ＪＰ，Ａ）

【文献】 特表２００９−５１６５６１（ＪＰ，Ａ）

【文献】 特開２０１１−１５３９７６（ＪＰ，Ａ）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ａ６１Ｂ １０／００

Ａ６１Ｂ ６／０３

Ａ６１Ｂ ８／００

Ｇ０１Ｎ ２１／１７

Ｇ０１Ｎ ２２／０２

ＪＳＴＰｌｕｓ／ＪＭＥＤＰｌｕｓ／ＪＳＴ７５８０（ＪＤｒｅａｍＩＩＩ）

(57)【特許請求の範囲】

【請求項1】

被検者の被計測部位に光を照射する光照射部と、
前記被計測部位からの拡散光を検出する光検出部と、
前記光検出部からの出力信号に基づいて、前記被計測部位内の光吸収係数分布を演算し、前記被計測部位の内部に関する再構成画像を作成する演算部と
を備え、
前記演算部は、前記再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗ_j（但し、添字ｊは１からＪまでの整数であり、Ｊは前記再構成画像の画素数である。）を算出し、次の反復式

【数1】

（但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘ_j^(k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄ_j^(k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより前記再構成画像を作成することを特徴とする、生体計測装置。

【請求項2】

被検者の被計測部位に放射線若しくは音波を照射する照射部と、
前記被計測部位からの拡散した前記放射線若しくは前記音波を検出する検出部と、
前記検出部からの出力信号に基づいて、前記被計測部位内における前記放射線若しくは前記音波の吸収係数分布を演算し、前記被計測部位の内部に関する再構成画像を作成する演算部と
を備え、
前記演算部は、前記再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗ_j（但し、添字ｊは１からＪまでの整数であり、Ｊは前記再構成画像の画素数である。）を算出し、次の反復式

【数2】

（但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘ_j^(k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄ_j^(k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより前記再構成画像を作成することを特徴とする、生体計測装置。

【請求項3】

前記演算部は、
前記再構成画像が分割されて成り複数の画素をそれぞれ含む複数の部分領域のうち、Ｎ回の反復演算において収束速度が最も遅い部分領域の前記収束速度（以下、最低収束速度という）Ｃ_Nを求め、
０＜ｖ_m＜１を満たすＭ個の数値ｖ_m（但し、ｍは１からＭまでの整数）を用意し、
１からＭまでの各ｍについて、次の反復式

【数3】

を用いてＮ回の反復演算を行うことによりＪ個の各画素の画素値ｘ₁^(N)〜ｘ_J^(N)を算出し、
該画素値ｘ₁^(N)〜ｘ_J^(N)から得られる各部分領域の前記収束速度が前記最低収束速度Ｃ_Nと略一致したときの数値ｖ_mを、該部分領域に含まれる複数の画素の前記係数ｗ_jとすることを特徴とする、請求項１または２に記載の生体計測装置。

【請求項4】

前記演算部は、Ｎ回の反復演算において収束速度が最も遅い部分領域の前記係数ｗ_jを１とすることを特徴とする、請求項３に記載の生体計測装置。

【請求項5】

被検者の被計測部位に光を照射し、前記被計測部位からの拡散光を検出し、該検出信号に基づいて前記被計測部位内の光吸収係数分布を演算し、前記被計測部位の内部に関する再構成画像を作成する方法であって、
前記再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗ_j（但し、添字ｊは１からＪまでの整数であり、Ｊは前記再構成画像の画素数である。）を算出し、
次の反復式

【数4】

（但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘ_j^(k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄ_j^(k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより前記再構成画像を作成することを特徴とする、画像作成方法。

【請求項6】

被検者の被計測部位に放射線若しくは音波を照射し、前記被計測部位からの拡散した前記放射線若しくは前記音波を検出し、該検出信号に基づいて前記被計測部位内における前記放射線若しくは前記音波の吸収係数分布を演算し、前記被計測部位の内部に関する再構成画像を作成する方法であって、
前記再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗ_j（但し、添字ｊは１からＪまでの整数であり、Ｊは前記再構成画像の画素数である。）を算出し、
次の反復式

【数5】

（但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘ_j^(k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄ_j^(k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより前記再構成画像を作成することを特徴とする、画像作成方法。

【請求項7】

前記再構成画像が分割されて成り複数の画素をそれぞれ含む複数の部分領域のうち、Ｎ回の反復演算において収束速度が最も遅い部分領域の前記収束速度（以下、最低収束速度という）Ｃ_Nを求め、
０＜ｖ_m＜１を満たすＭ個の数値ｖ_m（但し、ｍは１からＭまでの整数）を用意し、
１からＭまでの各ｍについて、次の反復式

【数6】

を用いてＮ回の反復演算を行うことによりＪ個の各画素の画素値ｘ₁^(N)〜ｘ_J^(N)を算出し、
該画素値ｘ₁^(N)〜ｘ_J^(N)から得られる各部分領域の前記収束速度が前記最低収束速度Ｃ_Nと略一致したときの数値ｖ_mを、該部分領域に含まれる複数の画素の前記係数ｗ_jとすることを特徴とする、請求項５または６に記載の画像作成方法。

【請求項8】

前記演算ステップにおいて、Ｎ回の反復演算において収束速度が最も遅い部分領域の前記係数ｗ_jを１とすることを特徴とする、請求項７に記載の画像作成方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、生体計測装置および画像作成方法に関するものである。

【背景技術】

【0002】

頭部や乳房といった生体の内部情報を非侵襲的に計測する装置として、生体の光吸収特性を利用して内部情報を得る、いわゆる拡散光トモグラフィ（ＤＯＴ；Diffuse Optical Tomography）を用いたものが提案されている。このような計測装置においては、計測対象となる生体の部位に対して所定の照射位置から光を照射し、当該部位の内部を散乱されつつ伝播された光を所定の検出位置で検出し、その強度や時間波形などの測定結果から、当該部位の内部情報、すなわち当該部位の内部に存在する腫瘍などの光吸収体に関する情報を得ることができる。なお、特許文献１には、拡散光トモグラフィを用いた生体計測方法が記載されている。また、非特許文献１及び２には、拡散光トモグラフィにおける逐次近似画像再構成法が記載されている。

【先行技術文献】

【特許文献】

【0003】

【特許文献1】特開２００１−２６４２４５号公報

【非特許文献】

【0004】

【非特許文献1】Y. Ueda, K. Ohta, M. Oda, M. Miwa, Y. Tsuchiya, and Y. Yamashita,"Three-dimensional imaging of a tissuelike phantom by diffusion opticaltomography", Applied Optics, 40, 34, 6349-6355 (2001)

【非特許文献2】Y. Ueda, T. Yamanaka, D. Yamashita, T. Suzuki, E. Ohmae, M. Oda andY. Yamashita, "Reflectance Diffuse Optical Tomography: an application forhuman brain mapping" Japanese Journal of Applied Physic, Part 2, 44, 38,2717-2723 (2005)

【発明の概要】

【発明が解決しようとする課題】

【0005】

拡散光トモグラフィによる生体計測では、被計測部位内部における位置によって空間分解能や雑音特性が異なるため、不均一な画像が生成される。図１７は、そのような現象を説明するための模式図であり、被計測部位１００と、被計測部位１００の表面に設けられた光照射部１０１及び光検出部１０２が示されている。被計測部位１００の内部では、光照射位置１０１から出射されて光検出位置１０２に達する光子の飛行時間が短いほど、飛行距離が短く飛行経路が限定される。逆に、光子の飛行時間が長いほど、飛行距離が長くなり飛行経路は限定されない。そして、光子の飛行時間が短いデータには、図１７に示される領域Ａ１、すなわち被計測部位の表面に近い領域を通過する飛行経路Ｒ１が多く含まれる。また、光子の飛行時間が長いデータには、図１７に示される領域Ａ２、すなわち被計測部位の表面から遠い領域を通過する飛行経路Ｒ２が多く含まれる。したがって、被計測部位の表面から遠い領域の情報量は、被計測部位の表面に近い領域の情報量より少なくなってしまう。これにより、被計測部位の表面に遠い領域の空間分解能や雑音が、被計測部位の表面から近い領域の空間分解能や雑音よりも大きくなってしまう。

【0006】

本発明は、このような問題点に鑑みてなされたものであり、被計測部位内部の位置による空間分解能や雑音特性の差を抑制し、より均一に近い画像を作成することができる生体計測装置および画像作成方法を提供することを目的とする。

【課題を解決するための手段】

【0007】

上述した課題を解決するために、本発明に係る第１の生体計測装置は、被検者の被計測部位に光を照射する光照射部と、被計測部位からの拡散光を検出する光検出部と、光検出部からの出力信号に基づいて、被計測部位内の光吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する演算部とを備え、演算部は、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗ_j（但し、添字ｊは１からＪまでの整数であり、Ｊは再構成画像の画素数である。）を算出し、次の反復式

【数1】

（但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘ_j^(k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄ_j^(k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより再構成画像を作成することを特徴とする。また、本発明に係る第２の生体計測装置は、被検者の被計測部位に放射線若しくは音波を照射する照射部と、被計測部位からの拡散した放射線若しくは音波を検出する検出部と、検出部からの出力信号に基づいて、被計測部位内における放射線若しくは音波の吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する演算部とを備え、演算部は、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗ_j（但し、添字ｊは１からＪまでの整数であり、Ｊは再構成画像の画素数である。）を算出し、上記反復式（１）を用いた逐次近似演算を行うことにより再構成画像を作成することを特徴とする。なお、本発明において、放射線とは例えばＸ線、γ線、若しくはマイクロ波といった短波長の電磁波をいい、音波とは例えば超音波のような波動をいう。

【0008】

これらの生体計測装置では、再構成画像の各画素毎に設定されるＪ個の係数ｗ₁〜ｗ_Jを用いて、画像再構成のための逐次近似演算を行っている。例えば、これらの係数ｗ₁〜ｗ_Jを、Ｎ回の反復演算において収束速度が最も遅い領域に当該画素の収束速度を一致させるように設定することで、収束速度が均一化され、被計測部位内部の位置による空間分解能や雑音特性の差を抑制してより均一に近い画像を作成することが可能となる。

【0009】

また、生体計測装置は、演算部が、再構成画像が分割されて成り複数の画素をそれぞれ含む複数の部分領域のうち、Ｎ回の反復演算において収束速度が最も遅い部分領域の収束速度（以下、最低収束速度という）Ｃ_Nを求め、０＜ｖ_m＜１を満たすＭ個の数値ｖ_m（但し、ｍは１からＭまでの整数）を用意し、１からＭまでの各ｍについて、次の反復式

【数2】

を用いてＮ回の反復演算を行うことによりＪ個の各画素の画素値ｘ₁^(N)〜ｘ_J^(N)を算出し、該画素値ｘ₁^(N)〜ｘ_J^(N)から得られる各部分領域の収束速度が最低収束速度Ｃ_Nと略一致したときの数値ｖ_mを、該部分領域に含まれる複数の画素の係数ｗ_jとすることを特徴としてもよい。これにより、上述した効果をより好適に得ることができる。

【0010】

また、生体計測装置は、演算部が、Ｎ回の反復演算において収束速度が最も遅い部分領域の係数ｗ_jを１とすることを特徴としてもよい。これにより、上述した効果をより好適に得ることができる。

【0011】

また、本発明に係る第１の画像作成方法は、被検者の被計測部位に光を照射し、被計測部位からの拡散光を検出し、該検出信号に基づいて被計測部位内の光吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する方法であって、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗ_j（但し、添字ｊは１からＪまでの整数であり、Ｊは再構成画像の画素数である。）を算出し、次の反復式

【数3】

（但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘ_j^(k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄ_j^(k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより再構成画像を作成することを特徴としてもよい。また、本発明に係る第２の画像作成方法は、被検者の被計測部位に放射線若しくは音波を照射し、被計測部位からの拡散した放射線若しくは音波を検出し、該検出信号に基づいて被計測部位内における放射線若しくは音波の吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する方法であって、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗ_j（但し、添字ｊは１からＪまでの整数であり、Ｊは再構成画像の画素数である。）を算出し、上記反復式（３）を用いた逐次近似演算を行うことにより再構成画像を作成することを特徴とする。なお、本発明において、放射線とは例えばＸ線、γ線、若しくはマイクロ波といった短波長の電磁波をいい、音波とは例えば超音波のような波動をいう。

【0012】

これらの画像作成方法では、再構成画像の各画素毎に設定されるＪ個の係数ｗ₁〜ｗ_Jを用いて、画像再構成のための逐次近似演算を行っている。例えば、これらの係数ｗ₁〜ｗ_Jを、Ｎ回の反復演算において収束速度が最も遅い領域に当該画素の収束速度を一致させるように設定することで、収束速度が均一化され、被計測部位内部の位置による空間分解能や雑音特性の差を抑制してより均一に近い画像を作成することが可能となる。

【0013】

また、画像作成方法は、再構成画像が分割されて成り複数の画素をそれぞれ含む複数の部分領域のうち、Ｎ回の反復演算において収束速度が最も遅い部分領域の収束速度（以下、最低収束速度という）Ｃ_Nを求め、０＜ｖ_m＜１を満たすＭ個の数値ｖ_m（但し、ｍは１からＭまでの整数）を用意し、１からＭまでの各ｍについて、次の反復式

【数4】

を用いてＮ回の反復演算を行うことによりＪ個の各画素の画素値ｘ₁^(N)〜ｘ_J^(N)を算出し、該画素値ｘ₁^(N)〜ｘ_J^(N)から得られる各部分領域の収束速度が最低収束速度Ｃ_Nと略一致したときの数値ｖ_mを、該部分領域に含まれる複数の画素の係数ｗ_jとすることを特徴としてもよい。これにより、上述した効果をより好適に得ることができる。

【0014】

また、画像作成方法は、Ｎ回の反復演算において収束速度が最も遅い部分領域の係数ｗ_jを１とすることを特徴としてもよい。これにより、上述した効果をより好適に得ることができる。

【発明の効果】

【0015】

本発明による生体計測装置および画像作成方法によれば、被計測部位内部の位置による空間分解能や雑音特性の差を抑制し、より均一に近い画像を作成することができる。

【図面の簡単な説明】

【0016】

【図1】本発明の一実施形態に係る生体計測装置の構成を示す図である。

【図2】係数（ステップサイズ）の具体的な決定方法を示すフローチャートである。

【図3】（ａ），（ｂ）シミュレーションにおいて再構成の対象となる２種類の被測定部位を示す図である。

【図4】被計測部位の内部における光子の飛行をシミュレーションする順問題解析における条件を示す図表である。

【図5】検出された光子ヒストグラムから画像を再構成する逆問題解析における条件を示す図表である。

【図6】シミュレーションにおいてステップサイズを決定するために使用された４つの画像を示す図である。

【図7】（ａ）図３（ａ）に示された画像に対応する、シミュレーションによる再構成後の画像である。（ｂ）図３（ａ）に示された画像に対応する、ステップサイズを用いない従来の方法による再構成後の画像である。

【図8】（ａ）〜（ｃ）図９に示される３本のライン上の画素値の変化を示すグラフである。

【図9】画像上に想定された３本のラインを示す図である。

【図10】（ａ）図３（ｂ）に示された画像に対応する、シミュレーションによる再構成後の画像である。（ｂ）図３（ｂ）に示された画像に対応する、ステップサイズを用いない従来の方法による再構成後の画像である。

【図11】（ａ）,（ｂ）図１２に示される２本のライン上の画素値の変化を示すグラフである。

【図12】画像上に想定された２本のラインを示す図である。

【図13】（ａ）図３（ａ）に示された画像に対応する、統計雑音を付加した場合の本実施形態による再構成後の画像である。（ｂ）図３（ａ）に示された画像に対応する、統計雑音を付加した場合の従来の方法による再構成後の画像である。

【図14】（ａ）図３（ｂ）に示された画像に対応する、統計雑音を付加した場合の本実施形態による再構成後の画像である。（ｂ）図３（ｂ）に示された画像に対応する、統計雑音を付加した場合の従来の方法による再構成後の画像である。

【図15】部分領域毎のＮＳＤを示すグラフである。

【図16】標準偏差の算出結果を示す図表である。

【図17】拡散光トモグラフィにおいて不均一な画像が生成される現象を説明するための模式図である。

【発明を実施するための形態】

【0017】

以下、添付図面を参照しながら本発明による生体計測装置および画像作成方法の実施の形態を詳細に説明する。なお、図面の説明において同一の要素には同一の符号を付し、重複する説明を省略する。

【0018】

図１は、本発明の一実施形態に係る生体計測装置１０の構成を示す図である。本実施形態の生体計測装置１０は、計測対象である被検者の被計測部位Ｂに光を照射し、拡散光（戻り光）を検出し、その検出位置と計測された光量データ（例えば時間分解光子ヒストグラム）とに基づいて、光子の平均飛行経路と平均光路長を推定し、画像再構成問題として体内の情報を画像化する装置である。この装置によって得られる画像は、例えば腫瘍の位置や酸素化ヘモグロビン及び脱酸素化ヘモグロビンの分布を可視化したものであり、体組織の機能画像である。なお、被計測部位Ｂとしては、例えば頭部や女性の乳房等が想定される。

【0019】

生体計測装置１０は、計測光を被計測部位Ｂの内部に照射する光照射部と、光照射部からの光の照射により被計測部位Ｂから生じた拡散光を検出する光検出部と、光検出部からの出力信号に基づいて被計測部位Ｂの吸収係数の空間的分布を計算し、被計測部位Ｂの再構成画像を作成する演算部１４とを備えている。

【0020】

本実施形態の光照射部は、被計測部位Ｂに取り付けられるｎ個の光出射／検出端１６それぞれが有する光出射端、光源２２、および光スイッチ２４によって構成されている。光源２２としては、例えばレーザダイオードを使用することができる。計測光の波長としては、生体の透過率と定量すべき吸収体の分吸収係数との関係等から、７００ｎｍ〜９００ｎｍ程度の近赤外線領域の波長が好ましい。

【0021】

計測光は、例えば連続光として光源２２から出射される。光源２２から出射された計測光は、光出射／検出端１６から被計測部位Ｂへ照射される。光スイッチ２４は、１入力ｎ出力の光スイッチであり、光源２２から光源用光ファイバ２６を介して光を入力し、この光を上記ｎ個の光出射／検出端１６それぞれに対して順に供給する。すなわち、光スイッチ２４は、各光出射／検出端１６に接続されたｎ本の出射用光ファイバ２８を１本ずつ順に選択し、当該出射用光ファイバ２８と光源２２とを光学的に接続する。

【0022】

本実施形態の光検出部は、前述したｎ個の光出射／検出端１６それぞれが有する光検出端と、ｎ個の光出射／検出端１６それぞれに対応するｎ個の光検出器３０と、各光検出器の入力部前段に配置されたｎ個のシャッター３２とによって構成されている。ｎ個の光検出器３０それぞれには、各光出射／検出端１６の光検出端に入射した被計測部位Ｂからの拡散光が、検出用光ファイバ３４を介して入力される。光検出器３０は、対応する光出射／検出端１６に到達した拡散光の光強度に応じてアナログ信号を生成する。光検出器３０としては、光電子増倍管（ＰＭＴ：PhotomultiplierTube）の他、フォトダイオード、アバランシェフォトダイオード、ＰＩＮフォトダイオード等、様々なものを使用することができる。被計測部位Ｂからの拡散光が微弱であるときは、高感度あるいは高利得の光検出器を使用することが好ましい。光検出器３０の信号出力端には信号処理回路３６が接続されており、信号処理回路３６は、光検出器３０から出力されたアナログ信号をＡ／Ｄ変換して拡散光の光強度に応じたディジタル信号を生成し、該ディジタル信号を演算部１４へ提供する。

【0023】

演算部１４は、信号処理回路３６から提供されたディジタル信号に基づいて、被計測部位Ｂ内の光吸収係数分布を演算し、被計測部位Ｂの内部に関する再構成画像を作成する。演算部１４は、例えば、ＣＰＵ（Central Processing Unit）といった演算手段及びメモリなどの記憶手段を有するコンピュータによって実現される。なお、演算部１４は、光源２２の発光、光スイッチ２４の動作及びシャッター３２の開閉を制御する機能を更に有するとよい。また、演算部１４には記録／表示部３８が接続されており、演算部１４における演算結果、すなわち被計測部位Ｂの再構成画像を可視化することが可能となっている。

【0024】

被計測部位Ｂの内部情報の算出すなわち内部情報計測は、例えば次のようにして行われる。ｎ個の光出射／検出端１６のそれぞれから被計測部位Ｂの内部へ計測光を順に照射し、被計測部位Ｂを通過して拡散した光を、ｎ個の光出射／検出端１６を介してｎ個の光検出器３０により検出する。この検出結果に基づいて、被計測部位Ｂの内部における吸収係数の空間的分布を演算し、腫瘍などの吸収体の位置や形状に関する情報（内部情報）を含む再構成画像を作成する。

【0025】

なお、演算部１４における吸収係数分布の算出には、例えば特許文献１に詳しく説明されているような周知の方法を用いるとよい。

【0026】

続いて、光吸収係数の空間的分布に基づいて、再構成画像を作成する方法について説明する。なお、以下に説明する演算は、演算部１４において行われる。ここで、拡散光トモグラフィにおける画像再構成問題を定式化するために、未知の光吸収係数分布に基づく再構成画像を構成する各画素の値を、次のＪ次元列ベクトルｘによって表す。
ｘ＝（ｘ₁，ｘ₂，…，ｘ_J）^T
また、光検出部において検出された測定データである光子ヒストグラムを、次のＩ次元列ベクトルＴによって表す。
Ｔ＝（Ｔ₁，Ｔ₂，…，Ｔ_I）^T
また、これらｘとＴとを相互に関連づけるＩ×Ｊ型のシステム行列Ｌを
Ｌ＝｛ｌ_ij｝
と定義する。ｌ_iをＬのｉ行目の要素ベクトルとする。更に、光吸収係数分布が一様で既知の画像を、Ｊ次元列ベクトルｘ_refとし、そのｘ_refに対応する測定データである光子ヒストグラムを、次のＩ次元列ベクトルＢによって表す。
Ｂ＝（Ｂ₁，Ｂ₂，…，Ｂ_I）^T

【0027】

測定データである光子ヒストグラムに統計雑音が含まれていない場合、次式（５）が成り立つ。

【数5】

しかし、統計雑音が混入した場合、上式（５）が成り立たない。したがって、統計雑音が混入している状態での最適なｘを求める必要がある。そこで、本実施形態ではいわゆる最尤推定法を用いてこのようなｘを求める。最尤推定法では、光検出部における光子の検出確率から尤度関数を立式し、その尤度関数を目的関数として最適化問題を解くことによって最適なｘを求めることができる。

【0028】

拡散光トモグラフィにおける光子の検出確率はポアソン分布に従い、その統計雑音もポアソン分布に従う。したがって、拡散光トモグラフィの最適化問題は、次式（６）によって表される。

【数6】

また、式（６）における目的関数F（ｘ）は、次式（７）に示される対数尤度関数によって表される。

【数7】

【0029】

本実施形態では、ＯＳ−Ｃｏｎｖｅｘ法を用いて上式（６）の最大化問題を勾配法で解くことにより、画像再構成を行う。ＯＳ−Ｃｏｎｖｅｘ法では、測定データをサブセットと呼ばれる部分データ集合に分割し、各サブセットに対応する評価関数が増加するようにサブセット毎に解を更新する。これにより、１回の反復演算においてサブセットの数だけ解を更新できるので、収束速度が向上する。

【0030】

ｋを１からＮまでの整数（Ｎは反復演算の回数）とし、サブセット数をＱとし、ｑを１からＱまでの整数とし、第ｑ番目のサブセットのデータ集合をＳｑと定義する。具体的なＯＳ−Ｃｏｎｖｅｘ法の反復式は、次式（８）及び（９）によって表される。

【数8】

【数9】

ここで、一般的な逐次近似解法の反復式は、次式（１０）に示される。この反復式によれば、第（ｋ＋１）番目の値ｘ_j^(k+1)は、第ｋ番目の値ｘ_j^(k)と更新量ｄ_j^(k)との和によって表される。

【数10】

なお、上述したＯＳ−Ｃｏｎｖｅｘ法では、更新量ｄ_j^(k)は次式（１１）のように表される。

【数11】

【0031】

解の収束（収束速度）が早い画素では、Ｎ回の反復演算それぞれにおける更新量ｄ_j⁽¹⁾〜ｄ_j^(N)の値が大きい。一方、収束速度が遅い画素では、Ｎ回の反復演算それぞれにおける更新量ｄ_j⁽¹⁾〜ｄ_j^(N)の値が小さい。したがって、本来は同値となるべき２つの画素の収束速度が互いに異なると、一定回数の反復演算後におけるそれぞれの値が互いに異なる結果となる。つまり、画素毎の収束速度が大きく異なると、再構成画像における空間分解能の均一性が損なわれてしまうこととなる。

【0032】

すなわち、Ｊ個の画素における値ｘ₁^(N)〜ｘ_J^(N)の収束速度を互いに近づけることにより、Ｎ回の反復演算後における再構成画像の空間分解能（画質）を均一化することができる。その為には、値ｘ₁^(N)〜ｘ_J^(N)の収束速度が互いに近づく（好ましくは、略等しくなる）ように、上式（１０）における更新量ｄ_j^(k)に対し画素毎に異なる任意の係数を乗ずるとよい。但し、更新量ｄ_j^(k)は、当該第ｊ番目の画素における評価関数を増加させる最大の更新量であるので、収束速度が遅い画素の更新量ｄ_j^(k)に対して１より大きい係数を乗ずると、評価関数は増加せず減少してしまうので好ましくない。そこで、収束速度が早い画素の更新量ｄ_j^(k)に対し、１より小さい係数を乗算することによって、評価関数を増加させつつ、収束速度が遅い画素に合わせて収束速度を揃えることが可能となる。換言すれば、次式（１２）を満たす係数ｗ_j（以下、ステップサイズという）を各画素に設定し、このステップサイズｗ_jを更新量ｄ_j^(k)に乗ずることによって、各画素値ｘ₁^(N)〜ｘ_J^(N)の収束速度を任意に制御することが可能となる。

【数12】

したがって、上述した反復演算式（１０）は、次式（１３）のように書き換えられる。

【数13】

例えば、上式（１３）に示される反復演算式を前述したＯＳ−Ｃｏｎｖｅｘ法に適用すると、式（９）に示された反復演算式は次式（１４）のようになる。

【数14】

【0033】

続いて、ステップサイズｗ_jを決定するための方法について説明する。

【0034】

＜収束速度の評価＞
ステップサイズｗ_jを決定するためには、測定データに基づいて画素毎の収束速度を求めたのち、収束速度が最も遅い画素を判定する必要がある。画素毎の収束速度は、コントラスト回復率（ＣＲＣ）によって評価することができる。

【0035】

拡散光トモグラフィの性格として、相互に隣接する画素の収束速度は互いに近いと考えられる。したがって本実施形態では、複数の画素をそれぞれ含む複数の部分領域に画像を分割し、スポットとなる領域を各部分領域内に配置し、該スポットにおけるＣＲＣをその領域の収束速度に相当する値として用いる。なお、ＣＲＣは、次式（１５）によって定義される。

【数15】

但し、ＳＰはスポット領域における画素値を表し、ＢＧはスポット領域外（背景領域）の画素値を表す。また、添字ｍは領域内の平均値であることを示し、添字Ｒは画像再構成後の画素値であることを示し、添字Ｔｒは測定データに基づく画素値であることを示す。

【0036】

＜ステップサイズｗ_jの決定＞
図２は、上述した収束速度の評価方法に基づく、ステップサイズｗ_jの具体的な決定方法を示すフローチャートである。この方法では、まず、画像をＥ個の部分領域に分割する（ステップＳ１）。なお、以下の説明において、或る任意の部分領域ｅの画素値の集合をＲ_ｅとする。

【0037】

次に、任意の反復回数Ｎを設定し、Ｎ回の反復演算においてＣＲＣが最も小さくなる（すなわち、収束速度が最も遅い）部分領域を判定し、そのＣＲＣを最低収束速度Ｃ_Nと定義する（ステップＳ２）。そして、０＜ｖ_m＜１を満たすＭ個の任意の数値ｖ_m（但し、ｍは１からＭまでの整数であり、ｖ₁〜ｖ_Mは互いに異なる値である）を用意する（ステップＳ３）。なお、このステップＳ３は、ステップＳ２やステップＳ１より前に行ってもよい。

【0038】

続いて、１からＭまでの各ｍについて、次の反復式（１６）を用いてＮ回の反復演算を行うことにより、Ｊ個の各画素の画素値ｘ₁^(N)〜ｘ_J^(N)を算出し、再構成画像を作成する（ステップＳ４）。

【数16】

そして、該画素値ｘ₁^(N)〜ｘ_J^(N)から得られる或る部分領域の収束速度が最低収束速度Ｃ_Nと略一致した場合に、そのときの数値ｖ_mを該部分領域に含まれる複数の画素のステップサイズｗ_jとする（ステップＳ５）。以降、前述したステップＳ４及びＳ５をＭ回繰り返すか、或いは、収束速度が最も遅い部分領域を除く全ての部分領域についてステップサイズｗ_jが決定するまで、ステップＳ４及びＳ５を繰り返す。

【0039】

最後に、Ｎ回の反復演算において収束速度が最も遅い部分領域のステップサイズｗ_jを１とする（ステップＳ６）。こうして、全ての部分領域の全ての画素について、ステップサイズｗ_jが決定される。

【0040】

以上に説明した本実施形態に係る生体計測装置１０及び画像作成方法によって得られる効果について、従来方法の課題とともに説明する。

【0041】

拡散光トモグラフィは、生体内での近赤外光の拡散的な飛行経路を用いた画像再構成法である。拡散光トモグラフィの画像再構成では、解析的手法ではなく逐次近似画像再構成法が用いられるが、共役勾配法といった通常の逐次近似画像再構成法では、再構成画像内での空間分解能や雑音特性のばらつきが大きい不均一な画像が生成される問題がある。

【0042】

画質の不均一性の原因は、拡散光トモグラフィにおいては測定データに含まれる情報量が被測定部位の内部の位置によって大きく異なるためである。すなわち、時間分解計測法を用いた拡散光トモグラフィーでは、被計測部位表面の光入射端から入射された光子は被計測部位内で散乱を繰り返しながら飛行し、被計測部位表面の光検出端に到達して検出される。この入射から検出までの光子の飛行時間が短いほど、飛行距離が短くなり飛行経路も限られる。逆に、飛行時間が長いほど、飛行距離が長くなり飛行経路は限定されない。このため、測定データに含まれる情報量は、光子の飛行時間に応じて変化する。典型的には、光子の飛行時間が短い場合には、被計測部位表面に近い飛行経路の割合が多くなり、また、飛行時間が長い場合には、被計測部位表面から遠い（被計測部位の中央付近の）飛行経路の割合が多くなる。そのため、被計測部位の中央付近の情報量は表面付近と比べて少なくなってしまう。これにより、再構成画像の周辺部における空間分解能及び雑音が再構成画像の中央付近における空間分解能及び雑音と比べて大きくなるといった、再構成画像内での空間分解能及び雑音のばらつきが生じてしまう。

【0043】

本実施形態の生体計測装置１０及び画像作成方法では、再構成画像の各画素毎に設定されるＪ個の係数ｗ₁〜ｗ_Jを用いて、画像再構成のための逐次近似演算を行う。そして、本実施形態のように、これらの係数ｗ₁〜ｗ_Jを、Ｎ回の反復演算において収束速度が最も遅い部分領域に当該画素の収束速度を一致させるように設定することで、収束速度が均一化され、被計測部位内部の位置による空間分解能や雑音特性の差を抑制して、より均一に近い画像を作成することが可能となる。

【0044】

ここで、本実施形態に係る画像作成方法による上記効果を確認するために行われたシミュレーションの結果について説明する。図３（ａ）及び図３（ｂ）は、このシミュレーションにおいて再構成の対象となる２種類の被測定部位を示す図であり、各図には、例えば腫瘍といった複数の光吸収体Ｄが示されている。なお、図３（ａ）に示される３つの光吸収体Ｄの中心座標（ｘ，ｙ）は、それぞれ（２１，２０）、（６４，１０６）、及び（７５，５１）である。また、図３（ｂ）に示される４つの光吸収体Ｄの中心座標（ｘ，ｙ）は、それぞれ（５６，７９）、（６６，２４）、（１０９，４９）、及び（１０９，１１０）である。

【0045】

このシミュレーションでは、次のステップ（１）〜（３）を行った。
（１）各画素の収束速度を評価する為に、画像内において均等に配置されたスポットにおける画像再構成を行い、ステップサイズｗ₁〜ｗ_Jを決定した。
（２）得られたステップサイズｗ₁〜ｗ_Jを用いて画像の再構成を行い、再構成画像の全体で空間分解能が略均一となっていることを確認した。
（３）背景画像の再構成による領域間での雑音伝播を評価した。

【0046】

なお、上記ステップ（１）〜（３）に共通する条件は、次のとおりである。まず、被計測部位の内部における光子の飛行をシミュレーションする順問題解析（Forward problem）と、検出された光子ヒストグラムから画像を再構成する逆問題解析（Inverse problem）という２種類の解析における条件を、それぞれ図４及び図５に示す。図４及び図５に示されるように、順問題解析と逆問題解析との間でグリッドサイズや画像サイズが異なるのは、連続系を離散系に変換したことによる。光検出部において検出される光子は、必ず離散的なデータとして取得される。一方、数値シミュレーションの値は連続関数に近い形で計算される。したがって、本シミュレーションでは、ダウンサンプリングすることにより連続系を離散系に変換している。

【0047】

（１）ステップサイズｗ₁〜ｗ_Jの決定
測定データに基づく１３２行１３２列の画素を含む画像を１２行１２列の部分領域に分割し、各部分領域の中央にホットスポットを配置した。本シミュレーションでは、図６（ａ）〜図６（ｄ）に示される４つの画像を再構成することにより、各部分領域のＣＲＣを求め、前述した方法によって全画素のステップサイズｗ₁〜ｗ_Jを決定した。

【0048】

（２）画像の再構成
次に、ステップ（１）により決定されたステップサイズｗ₁〜ｗ_Jを用い、式（１４）に示された反復演算式を演算することにより、画像の再構成を行った。なお、このステップでは、各部分領域のＣＲＣに差があると画像中にエッジとなって現れてしまうので、ステップサイズｗ₁〜ｗ_Jを各画素値とする画像（ステップサイズ画像）を平均値フィルターを用いてぼかし、滑らかに変化するステップサイズ画像を得た。本シミュレーションでは、カーネルサイズ９×９の平均値フィルターを用いた。

【0049】

図７（ａ）は、図３（ａ）に示された画像に対応する、本シミュレーションによる再構成後の画像である。図７（ｂ）は、図３（ａ）に示された画像に対応する、ステップサイズｗ₁〜ｗ_Jを用いない従来の方法による再構成後の画像である。図８（ａ）〜図８（ｃ）は、図９に示されるように光吸収体Ｄを通る３本のラインＬ１〜Ｌ３上の画素値の変化を示すグラフであって、グラフＧ１１は本シミュレーションによる画素値の変化を示し、グラフＧ１２は従来の方法による画素値の変化を示している。また、図１０（ａ）は、図３（ｂ）に示された画像に対応する、本シミュレーションによる再構成後の画像である。図１０（ｂ）は、図３（ｂ）に示された画像に対応する、ステップサイズｗ₁〜ｗ_Jを用いない従来の方法による再構成後の画像である。図１１（ａ）及び図１１（ｂ）は、図１２に示されるように光吸収体Ｄを通る２本のラインＬ４及びＬ５上の画素値の変化を示すグラフであって、グラフＧ２１は本シミュレーションによる画素値の変化を示し、グラフＧ２２は従来の方法による画素値の変化を示している。

【0050】

図７（ｂ）を参照すると、従来の再構成画像では、左上といった画像周辺部の光吸収体Ｄの画素値が比較的高くなっており、画像中央付近の光吸収体Ｄの画素値と大きく異なっていることがわかる。これに対し、ステップサイズｗ₁〜ｗ_Jを用いた本実施形態による再構成画像では、図７（ａ）に示されるように、各光吸収体Ｄの画素値が、画像の周辺部及び中央付近においてほぼ均一な値となっている。また、図８（ｄ）及び図８（ｅ）を参照すると、従来の再構成画像では、画像周辺部における画素値の変化が大きい。これに対し、図８（ａ）〜図８（ｃ）に示されるように、本実施形態による再構成画像では画像周辺部における画素値の変化が小さくなっている。これは、本実施形態において画像中央付近ではステップサイズｗ_jが大きくなって反復演算における更新量が大きくなり、画像周辺部ではステップサイズｗ_jが小さくなって更新量が小さくなったことによる。以上の結果から、本実施形態の生体計測装置１０及び画像作成方法によれば、収束速度が均一化されたことにより、被計測部位内部の位置による空間分解能の差を抑制して、より均一に近い画像を作成できることがわかる。

【0051】

次に、統計雑音を付加した場合のシミュレーション結果を示す。このシミュレーションでは、各検出光子ヒストグラムの最大値が一定値（例えば５０）となるように光子ヒストグラムをフィッティングした後、この光子ヒストグラムにポアソン雑音を付加した。図１３（ａ）は、図３（ａ）に示された画像に対応する、統計雑音を付加した場合の本実施形態による再構成後の画像である。図１３（ｂ）は、図３（ａ）に示された画像に対応する、統計雑音を付加した場合の従来の方法による再構成後の画像である。また、図１４（ａ）は、図３（ｂ）に示された画像に対応する、統計雑音を付加した場合の本実施形態による再構成後の画像である。図１４（ｂ）は、図３（ｂ）に示された画像に対応する、統計雑音を付加した場合の従来の方法による再構成後の画像である。

【0052】

図１３（ｂ）及び図１４（ｂ）を参照すると、統計雑音が存在する場合、ステップサイズｗ₁〜ｗ_Jを用いない従来の方法では、画像周辺部の空間分解能が高くなって雑音の影響が大きく出ているが、画像中央付近では空間分解能が高くなく雑音の影響が小さいことがわかる。これに対し、図１３（ａ）及び図１４（ａ）を参照すると、ステップサイズｗ₁〜ｗ_Jを用いた本実施形態による方法では、反復演算における画像周辺部の更新量が抑えられているので、雑音の影響が顕著に抑制され、左上に存在する光吸収体Ｄの形を明確に視認できる。一方、画像中央付近では、雑音の影響が従来方法と同程度に抑えられている。なお、画像中央付近の統計雑音が画像周辺部の統計雑音よりも大きいのは、統計雑音の有無によって分解能の向上速度に変化が生じるためであり、統計雑音がない場合を仮定して決定したステップサイズを用いたことに起因する現象である。

【0053】

このように、計測データに統計雑音が含まれる場合、本実施形態の生体測定方法では従来方法と比較して画素値に変化はあるものの、画像全体の統計雑音の形態的特徴に大きな変化がないことがわかる。すなわち、本実施形態の生体測定方法によれば、画像周辺部において選択的に、分解能の向上を遅くすることができていることがわかる。

【0054】

続いて、光吸収体Ｄを配置しない背景画像に関する雑音評価について説明する。ここでは、光吸収体Ｄが存在しない計測データにポアソン雑音を付加した。このときの再構成画像は、雑音によって歪んだ画像となる。この再構成画像を６行６列の３６個の部分領域に分割し、左上から右下へ順に番号を付与した。その後、それぞれの部分領域を雑音指標（ＮＳＤ）によって評価した。なお、ＮＳＤの計算式は次式（１７）である。また、図１５は、部分領域毎のＮＳＤを示すグラフであって、横軸は部分領域の番号を示し、縦軸はＮＳＤの値を示している。なお、添字ＳＤは標準偏差であることを示す。

【数17】

また、部分領域毎のＮＳＤを母集団とする標準偏差を算出し、この標準偏差を部分領域間の雑音の不均一性を評価する値として用いる。図１６は、標準偏差の算出結果を示す図表である。

【0055】

図１５を参照すると、従来方法では部分領域間のＮＳＤのばらつきが大きい。これに対し、本実施形態の方法では、ＮＳＤのばらつきが全体的に抑えられていることがわかる。このことから、本実施形態の方法では、空間分解能を均一に近づけるために用いられるステップサイズｗ_jが、雑音伝播の不均一をも緩和していることがわかる。更に、図１６に示されるように、雑音の不均一性についても、本実施形態の方法では従来方法よりも小さくなる。

【0056】

本発明による生体計測装置および画像作成方法は、上述した実施形態に限られるものではなく、他に様々な変形が可能である。例えば、上記実施形態では被検者の被計測部位に光を照射し、被計測部位からの拡散光を検出し、該検出信号に基づいて被計測部位内における光吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成している。しかしながら、本発明は、光を用いた生体計測装置や画像作成方法に限らず、放射線や音波を用いた生体計測装置や画像作成方法にも適用可能である。すなわち、本発明に係る生体計測装置は、被検者の被計測部位に放射線若しくは音波を照射する照射部と、被計測部位からの拡散した放射線若しくは音波を検出する検出部と、検出部からの出力信号に基づいて、被計測部位内における放射線若しくは音波の吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する演算部とを備えてもよい。この場合、演算部は、上述した実施形態と同様の手法により、再構成画像を作成するとよい。また、本発明に係る画像作成方法は、被検者の被計測部位に放射線若しくは音波を照射し、被計測部位からの拡散した放射線若しくは音波を検出し、該検出信号に基づいて被計測部位内における放射線若しくは音波の吸収係数分布を演算し、被計測部位の内部に関する再構成画像を上述した実施形態と同様の手法により作成する方法であってもよい。ここで、放射線とは例えばＸ線、γ線、若しくはマイクロ波といった短波長の電磁波をいい、音波とは例えば超音波のような波動をいう。

【符号の説明】

【0057】

１０…生体計測装置、１４…演算部、１６…光出射／検出端、２２…光源、２４…光スイッチ、２６…光源用光ファイバ、２８…出射用光ファイバ、３０…光検出器、３２…シャッター、３４…検出用光ファイバ、３６…信号処理回路、３８…表示部。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図9】

【図12】

【図16】

【図17】

【図7】

【図8】

【図10】

【図11】

【図13】

【図14】

【図15】