JP5782314B2 - 生体計測装置および画像作成方法 - Google Patents

生体計測装置および画像作成方法 Download PDF

Info

Publication number
JP5782314B2
JP5782314B2 JP2011151086A JP2011151086A JP5782314B2 JP 5782314 B2 JP5782314 B2 JP 5782314B2 JP 2011151086 A JP2011151086 A JP 2011151086A JP 2011151086 A JP2011151086 A JP 2011151086A JP 5782314 B2 JP5782314 B2 JP 5782314B2
Authority
JP
Japan
Prior art keywords
pixel
reconstructed image
measurement site
image
light
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2011151086A
Other languages
English (en)
Other versions
JP2013019696A (ja
Inventor
博幸 工藤
博幸 工藤
直哉 齋藤
直哉 齋藤
之雄 上田
之雄 上田
健治 佳元
健治 佳元
山下　豊
豊 山下
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hamamatsu Photonics KK
Original Assignee
Hamamatsu Photonics KK
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hamamatsu Photonics KK filed Critical Hamamatsu Photonics KK
Priority to JP2011151086A priority Critical patent/JP5782314B2/ja
Priority to US14/130,690 priority patent/US9282932B2/en
Priority to PCT/JP2012/066552 priority patent/WO2013005635A1/ja
Priority to CN201280033086.7A priority patent/CN103635148B/zh
Priority to EP12807701.3A priority patent/EP2730231B1/en
Publication of JP2013019696A publication Critical patent/JP2013019696A/ja
Application granted granted Critical
Publication of JP5782314B2 publication Critical patent/JP5782314B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B10/00—Other methods or instruments for diagnosis, e.g. instruments for taking a cell sample, for biopsy, for vaccination diagnosis; Sex determination; Ovulation-period determination; Throat striking implements
    • A61B10/0041—Detection of breast cancer
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/0059—Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence
    • A61B5/0073—Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence by tomography, i.e. reconstruction of 3D images from 2D projections
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/0059—Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence
    • A61B5/0082—Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence adapted for particular medical purposes
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5211—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data
    • A61B6/5229—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data combining image data of a patient, e.g. combining a functional image with an anatomical image
    • A61B6/5247—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data combining image data of a patient, e.g. combining a functional image with an anatomical image combining images from an ionising-radiation diagnostic technique and a non-ionising radiation diagnostic technique, e.g. X-ray and ultrasound
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/13—Tomography
    • A61B8/15—Transmission-tomography
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/5207—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of raw data to produce diagnostic data, e.g. for generating an image
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/47—Scattering, i.e. diffuse reflection
    • G01N21/4795—Scattering, i.e. diffuse reflection spatially resolved investigating of object in scattering medium

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Health & Medical Sciences (AREA)
  • Surgery (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Pathology (AREA)
  • Molecular Biology (AREA)
  • Veterinary Medicine (AREA)
  • Animal Behavior & Ethology (AREA)
  • Medical Informatics (AREA)
  • Public Health (AREA)
  • Biophysics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Signal Processing (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Radiology & Medical Imaging (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Psychiatry (AREA)
  • Physiology (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Algebra (AREA)
  • Oncology (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)
  • Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)

Description

本発明は、生体計測装置および画像作成方法に関するものである。
頭部や乳房といった生体の内部情報を非侵襲的に計測する装置として、生体の光吸収特性を利用して内部情報を得る、いわゆる拡散光トモグラフィ（ＤＯＴ；Diffuse Optical Tomography）を用いたものが提案されている。このような計測装置においては、計測対象となる生体の部位に対して所定の照射位置から光を照射し、当該部位の内部を散乱されつつ伝播された光を所定の検出位置で検出し、その強度や時間波形などの測定結果から、当該部位の内部情報、すなわち当該部位の内部に存在する腫瘍などの光吸収体に関する情報を得ることができる。なお、特許文献１には、拡散光トモグラフィを用いた生体計測方法が記載されている。また、非特許文献１及び２には、拡散光トモグラフィにおける逐次近似画像再構成法が記載されている。
特開２００１−２６４２４５号公報
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) 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)
拡散光トモグラフィによる生体計測では、被計測部位内部における位置によって空間分解能や雑音特性が異なるため、不均一な画像が生成される。図１７は、そのような現象を説明するための模式図であり、被計測部位１００と、被計測部位１００の表面に設けられた光照射部１０１及び光検出部１０２が示されている。被計測部位１００の内部では、光照射位置１０１から出射されて光検出位置１０２に達する光子の飛行時間が短いほど、飛行距離が短く飛行経路が限定される。逆に、光子の飛行時間が長いほど、飛行距離が長くなり飛行経路は限定されない。そして、光子の飛行時間が短いデータには、図１７に示される領域Ａ１、すなわち被計測部位の表面に近い領域を通過する飛行経路Ｒ１が多く含まれる。また、光子の飛行時間が長いデータには、図１７に示される領域Ａ２、すなわち被計測部位の表面から遠い領域を通過する飛行経路Ｒ２が多く含まれる。したがって、被計測部位の表面から遠い領域の情報量は、被計測部位の表面に近い領域の情報量より少なくなってしまう。これにより、被計測部位の表面に遠い領域の空間分解能や雑音が、被計測部位の表面から近い領域の空間分解能や雑音よりも大きくなってしまう。
本発明は、このような問題点に鑑みてなされたものであり、被計測部位内部の位置による空間分解能や雑音特性の差を抑制し、より均一に近い画像を作成することができる生体計測装置および画像作成方法を提供することを目的とする。
上述した課題を解決するために、本発明に係る第１の生体計測装置は、被検者の被計測部位に光を照射する光照射部と、被計測部位からの拡散光を検出する光検出部と、光検出部からの出力信号に基づいて、被計測部位内の光吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する演算部とを備え、演算部は、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは再構成画像の画素数である。）を算出し、次の反復式

（但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘj (k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄj (k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより再構成画像を作成することを特徴とする。また、本発明に係る第２の生体計測装置は、被検者の被計測部位に放射線若しくは音波を照射する照射部と、被計測部位からの拡散した放射線若しくは音波を検出する検出部と、検出部からの出力信号に基づいて、被計測部位内における放射線若しくは音波の吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する演算部とを備え、演算部は、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは再構成画像の画素数である。）を算出し、上記反復式（１）を用いた逐次近似演算を行うことにより再構成画像を作成することを特徴とする。なお、本発明において、放射線とは例えばＸ線、γ線、若しくはマイクロ波といった短波長の電磁波をいい、音波とは例えば超音波のような波動をいう。
これらの生体計測装置では、再構成画像の各画素毎に設定されるＪ個の係数ｗ1〜ｗJを用いて、画像再構成のための逐次近似演算を行っている。例えば、これらの係数ｗ1〜ｗJを、Ｎ回の反復演算において収束速度が最も遅い領域に当該画素の収束速度を一致させるように設定することで、収束速度が均一化され、被計測部位内部の位置による空間分解能や雑音特性の差を抑制してより均一に近い画像を作成することが可能となる。
また、生体計測装置は、演算部が、再構成画像が分割されて成り複数の画素をそれぞれ含む複数の部分領域のうち、Ｎ回の反復演算において収束速度が最も遅い部分領域の収束速度（以下、最低収束速度という）ＣNを求め、０＜ｖm＜１を満たすＭ個の数値ｖm（但し、ｍは１からＭまでの整数）を用意し、１からＭまでの各ｍについて、次の反復式

を用いてＮ回の反復演算を行うことによりＪ個の各画素の画素値ｘ1 (N)〜ｘJ (N)を算出し、該画素値ｘ1 (N)〜ｘJ (N)から得られる各部分領域の収束速度が最低収束速度ＣNと略一致したときの数値ｖmを、該部分領域に含まれる複数の画素の係数ｗjとすることを特徴としてもよい。これにより、上述した効果をより好適に得ることができる。
また、生体計測装置は、演算部が、Ｎ回の反復演算において収束速度が最も遅い部分領域の係数ｗjを１とすることを特徴としてもよい。これにより、上述した効果をより好適に得ることができる。
また、本発明に係る第１の画像作成方法は、被検者の被計測部位に光を照射し、被計測部位からの拡散光を検出し、該検出信号に基づいて被計測部位内の光吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する方法であって、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは再構成画像の画素数である。）を算出し、次の反復式

（但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘj (k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄj (k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより再構成画像を作成することを特徴としてもよい。また、本発明に係る第２の画像作成方法は、被検者の被計測部位に放射線若しくは音波を照射し、被計測部位からの拡散した放射線若しくは音波を検出し、該検出信号に基づいて被計測部位内における放射線若しくは音波の吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する方法であって、再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは再構成画像の画素数である。）を算出し、上記反復式（３）を用いた逐次近似演算を行うことにより再構成画像を作成することを特徴とする。なお、本発明において、放射線とは例えばＸ線、γ線、若しくはマイクロ波といった短波長の電磁波をいい、音波とは例えば超音波のような波動をいう。
これらの画像作成方法では、再構成画像の各画素毎に設定されるＪ個の係数ｗ1〜ｗJを用いて、画像再構成のための逐次近似演算を行っている。例えば、これらの係数ｗ1〜ｗJを、Ｎ回の反復演算において収束速度が最も遅い領域に当該画素の収束速度を一致させるように設定することで、収束速度が均一化され、被計測部位内部の位置による空間分解能や雑音特性の差を抑制してより均一に近い画像を作成することが可能となる。
また、画像作成方法は、再構成画像が分割されて成り複数の画素をそれぞれ含む複数の部分領域のうち、Ｎ回の反復演算において収束速度が最も遅い部分領域の収束速度（以下、最低収束速度という）ＣNを求め、０＜ｖm＜１を満たすＭ個の数値ｖm（但し、ｍは１からＭまでの整数）を用意し、１からＭまでの各ｍについて、次の反復式

を用いてＮ回の反復演算を行うことによりＪ個の各画素の画素値ｘ1 (N)〜ｘJ (N)を算出し、該画素値ｘ1 (N)〜ｘJ (N)から得られる各部分領域の収束速度が最低収束速度ＣNと略一致したときの数値ｖmを、該部分領域に含まれる複数の画素の係数ｗjとすることを特徴としてもよい。これにより、上述した効果をより好適に得ることができる。
また、画像作成方法は、Ｎ回の反復演算において収束速度が最も遅い部分領域の係数ｗjを１とすることを特徴としてもよい。これにより、上述した効果をより好適に得ることができる。
本発明による生体計測装置および画像作成方法によれば、被計測部位内部の位置による空間分解能や雑音特性の差を抑制し、より均一に近い画像を作成することができる。
本発明の一実施形態に係る生体計測装置の構成を示す図である。 係数（ステップサイズ）の具体的な決定方法を示すフローチャートである。 （ａ），（ｂ）シミュレーションにおいて再構成の対象となる２種類の被測定部位を示す図である。 被計測部位の内部における光子の飛行をシミュレーションする順問題解析における条件を示す図表である。 検出された光子ヒストグラムから画像を再構成する逆問題解析における条件を示す図表である。 シミュレーションにおいてステップサイズを決定するために使用された４つの画像を示す図である。 （ａ）図３（ａ）に示された画像に対応する、シミュレーションによる再構成後の画像である。（ｂ）図３（ａ）に示された画像に対応する、ステップサイズを用いない従来の方法による再構成後の画像である。 （ａ）〜（ｃ）図９に示される３本のライン上の画素値の変化を示すグラフである。 画像上に想定された３本のラインを示す図である。 （ａ）図３（ｂ）に示された画像に対応する、シミュレーションによる再構成後の画像である。（ｂ）図３（ｂ）に示された画像に対応する、ステップサイズを用いない従来の方法による再構成後の画像である。 （ａ）,（ｂ）図１２に示される２本のライン上の画素値の変化を示すグラフである。 画像上に想定された２本のラインを示す図である。 （ａ）図３（ａ）に示された画像に対応する、統計雑音を付加した場合の本実施形態による再構成後の画像である。（ｂ）図３（ａ）に示された画像に対応する、統計雑音を付加した場合の従来の方法による再構成後の画像である。 （ａ）図３（ｂ）に示された画像に対応する、統計雑音を付加した場合の本実施形態による再構成後の画像である。（ｂ）図３（ｂ）に示された画像に対応する、統計雑音を付加した場合の従来の方法による再構成後の画像である。 部分領域毎のＮＳＤを示すグラフである。 標準偏差の算出結果を示す図表である。 拡散光トモグラフィにおいて不均一な画像が生成される現象を説明するための模式図である。
以下、添付図面を参照しながら本発明による生体計測装置および画像作成方法の実施の形態を詳細に説明する。なお、図面の説明において同一の要素には同一の符号を付し、重複する説明を省略する。
図１は、本発明の一実施形態に係る生体計測装置１０の構成を示す図である。本実施形態の生体計測装置１０は、計測対象である被検者の被計測部位Ｂに光を照射し、拡散光（戻り光）を検出し、その検出位置と計測された光量データ（例えば時間分解光子ヒストグラム）とに基づいて、光子の平均飛行経路と平均光路長を推定し、画像再構成問題として体内の情報を画像化する装置である。この装置によって得られる画像は、例えば腫瘍の位置や酸素化ヘモグロビン及び脱酸素化ヘモグロビンの分布を可視化したものであり、体組織の機能画像である。なお、被計測部位Ｂとしては、例えば頭部や女性の乳房等が想定される。
生体計測装置１０は、計測光を被計測部位Ｂの内部に照射する光照射部と、光照射部からの光の照射により被計測部位Ｂから生じた拡散光を検出する光検出部と、光検出部からの出力信号に基づいて被計測部位Ｂの吸収係数の空間的分布を計算し、被計測部位Ｂの再構成画像を作成する演算部１４とを備えている。
本実施形態の光照射部は、被計測部位Ｂに取り付けられるｎ個の光出射／検出端１６それぞれが有する光出射端、光源２２、および光スイッチ２４によって構成されている。光源２２としては、例えばレーザダイオードを使用することができる。計測光の波長としては、生体の透過率と定量すべき吸収体の分吸収係数との関係等から、７００ｎｍ〜９００ｎｍ程度の近赤外線領域の波長が好ましい。
計測光は、例えば連続光として光源２２から出射される。光源２２から出射された計測光は、光出射／検出端１６から被計測部位Ｂへ照射される。光スイッチ２４は、１入力ｎ出力の光スイッチであり、光源２２から光源用光ファイバ２６を介して光を入力し、この光を上記ｎ個の光出射／検出端１６それぞれに対して順に供給する。すなわち、光スイッチ２４は、各光出射／検出端１６に接続されたｎ本の出射用光ファイバ２８を１本ずつ順に選択し、当該出射用光ファイバ２８と光源２２とを光学的に接続する。
本実施形態の光検出部は、前述したｎ個の光出射／検出端１６それぞれが有する光検出端と、ｎ個の光出射／検出端１６それぞれに対応するｎ個の光検出器３０と、各光検出器の入力部前段に配置されたｎ個のシャッター３２とによって構成されている。ｎ個の光検出器３０それぞれには、各光出射／検出端１６の光検出端に入射した被計測部位Ｂからの拡散光が、検出用光ファイバ３４を介して入力される。光検出器３０は、対応する光出射／検出端１６に到達した拡散光の光強度に応じてアナログ信号を生成する。光検出器３０としては、光電子増倍管（ＰＭＴ：PhotomultiplierTube）の他、フォトダイオード、アバランシェフォトダイオード、ＰＩＮフォトダイオード等、様々なものを使用することができる。被計測部位Ｂからの拡散光が微弱であるときは、高感度あるいは高利得の光検出器を使用することが好ましい。光検出器３０の信号出力端には信号処理回路３６が接続されており、信号処理回路３６は、光検出器３０から出力されたアナログ信号をＡ／Ｄ変換して拡散光の光強度に応じたディジタル信号を生成し、該ディジタル信号を演算部１４へ提供する。
演算部１４は、信号処理回路３６から提供されたディジタル信号に基づいて、被計測部位Ｂ内の光吸収係数分布を演算し、被計測部位Ｂの内部に関する再構成画像を作成する。演算部１４は、例えば、ＣＰＵ（Central Processing Unit）といった演算手段及びメモリなどの記憶手段を有するコンピュータによって実現される。なお、演算部１４は、光源２２の発光、光スイッチ２４の動作及びシャッター３２の開閉を制御する機能を更に有するとよい。また、演算部１４には記録／表示部３８が接続されており、演算部１４における演算結果、すなわち被計測部位Ｂの再構成画像を可視化することが可能となっている。
被計測部位Ｂの内部情報の算出すなわち内部情報計測は、例えば次のようにして行われる。ｎ個の光出射／検出端１６のそれぞれから被計測部位Ｂの内部へ計測光を順に照射し、被計測部位Ｂを通過して拡散した光を、ｎ個の光出射／検出端１６を介してｎ個の光検出器３０により検出する。この検出結果に基づいて、被計測部位Ｂの内部における吸収係数の空間的分布を演算し、腫瘍などの吸収体の位置や形状に関する情報（内部情報）を含む再構成画像を作成する。
なお、演算部１４における吸収係数分布の算出には、例えば特許文献１に詳しく説明されているような周知の方法を用いるとよい。
続いて、光吸収係数の空間的分布に基づいて、再構成画像を作成する方法について説明する。なお、以下に説明する演算は、演算部１４において行われる。ここで、拡散光トモグラフィにおける画像再構成問題を定式化するために、未知の光吸収係数分布に基づく再構成画像を構成する各画素の値を、次のＪ次元列ベクトルｘによって表す。
ｘ＝（ｘ1，ｘ2，…，ｘJ）T
また、光検出部において検出された測定データである光子ヒストグラムを、次のＩ次元列ベクトルＴによって表す。
Ｔ＝（Ｔ1，Ｔ2，…，ＴI）T
また、これらｘとＴとを相互に関連づけるＩ×Ｊ型のシステム行列Ｌを
Ｌ＝｛ｌij｝
と定義する。ｌiをＬのｉ行目の要素ベクトルとする。更に、光吸収係数分布が一様で既知の画像を、Ｊ次元列ベクトルｘrefとし、そのｘrefに対応する測定データである光子ヒストグラムを、次のＩ次元列ベクトルＢによって表す。
Ｂ＝（Ｂ1，Ｂ2，…，ＢI）T
測定データである光子ヒストグラムに統計雑音が含まれていない場合、次式（５）が成り立つ。

しかし、統計雑音が混入した場合、上式（５）が成り立たない。したがって、統計雑音が混入している状態での最適なｘを求める必要がある。そこで、本実施形態ではいわゆる最尤推定法を用いてこのようなｘを求める。最尤推定法では、光検出部における光子の検出確率から尤度関数を立式し、その尤度関数を目的関数として最適化問題を解くことによって最適なｘを求めることができる。
拡散光トモグラフィにおける光子の検出確率はポアソン分布に従い、その統計雑音もポアソン分布に従う。したがって、拡散光トモグラフィの最適化問題は、次式（６）によって表される。

また、式（６）における目的関数F（ｘ）は、次式（７）に示される対数尤度関数によって表される。
本実施形態では、ＯＳ−Ｃｏｎｖｅｘ法を用いて上式（６）の最大化問題を勾配法で解くことにより、画像再構成を行う。ＯＳ−Ｃｏｎｖｅｘ法では、測定データをサブセットと呼ばれる部分データ集合に分割し、各サブセットに対応する評価関数が増加するようにサブセット毎に解を更新する。これにより、１回の反復演算においてサブセットの数だけ解を更新できるので、収束速度が向上する。
ｋを１からＮまでの整数（Ｎは反復演算の回数）とし、サブセット数をＱとし、ｑを１からＱまでの整数とし、第ｑ番目のサブセットのデータ集合をＳｑと定義する。具体的なＯＳ−Ｃｏｎｖｅｘ法の反復式は、次式（８）及び（９）によって表される。


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

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

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

例えば、上式（１３）に示される反復演算式を前述したＯＳ−Ｃｏｎｖｅｘ法に適用すると、式（９）に示された反復演算式は次式（１４）のようになる。
続いて、ステップサイズｗjを決定するための方法について説明する。
＜収束速度の評価＞
ステップサイズｗjを決定するためには、測定データに基づいて画素毎の収束速度を求めたのち、収束速度が最も遅い画素を判定する必要がある。画素毎の収束速度は、コントラスト回復率（ＣＲＣ）によって評価することができる。
拡散光トモグラフィの性格として、相互に隣接する画素の収束速度は互いに近いと考えられる。したがって本実施形態では、複数の画素をそれぞれ含む複数の部分領域に画像を分割し、スポットとなる領域を各部分領域内に配置し、該スポットにおけるＣＲＣをその領域の収束速度に相当する値として用いる。なお、ＣＲＣは、次式（１５）によって定義される。

但し、ＳＰはスポット領域における画素値を表し、ＢＧはスポット領域外（背景領域）の画素値を表す。また、添字ｍは領域内の平均値であることを示し、添字Ｒは画像再構成後の画素値であることを示し、添字Ｔｒは測定データに基づく画素値であることを示す。
＜ステップサイズｗjの決定＞
図２は、上述した収束速度の評価方法に基づく、ステップサイズｗjの具体的な決定方法を示すフローチャートである。この方法では、まず、画像をＥ個の部分領域に分割する（ステップＳ１）。なお、以下の説明において、或る任意の部分領域ｅの画素値の集合をＲｅとする。
次に、任意の反復回数Ｎを設定し、Ｎ回の反復演算においてＣＲＣが最も小さくなる（すなわち、収束速度が最も遅い）部分領域を判定し、そのＣＲＣを最低収束速度ＣNと定義する（ステップＳ２）。そして、０＜ｖm＜１を満たすＭ個の任意の数値ｖm（但し、ｍは１からＭまでの整数であり、ｖ1〜ｖMは互いに異なる値である）を用意する（ステップＳ３）。なお、このステップＳ３は、ステップＳ２やステップＳ１より前に行ってもよい。
続いて、１からＭまでの各ｍについて、次の反復式（１６）を用いてＮ回の反復演算を行うことにより、Ｊ個の各画素の画素値ｘ1 (N)〜ｘJ (N)を算出し、再構成画像を作成する（ステップＳ４）。

そして、該画素値ｘ1 (N)〜ｘJ (N)から得られる或る部分領域の収束速度が最低収束速度ＣNと略一致した場合に、そのときの数値ｖmを該部分領域に含まれる複数の画素のステップサイズｗjとする（ステップＳ５）。以降、前述したステップＳ４及びＳ５をＭ回繰り返すか、或いは、収束速度が最も遅い部分領域を除く全ての部分領域についてステップサイズｗjが決定するまで、ステップＳ４及びＳ５を繰り返す。
最後に、Ｎ回の反復演算において収束速度が最も遅い部分領域のステップサイズｗjを１とする（ステップＳ６）。こうして、全ての部分領域の全ての画素について、ステップサイズｗjが決定される。
以上に説明した本実施形態に係る生体計測装置１０及び画像作成方法によって得られる効果について、従来方法の課題とともに説明する。
拡散光トモグラフィは、生体内での近赤外光の拡散的な飛行経路を用いた画像再構成法である。拡散光トモグラフィの画像再構成では、解析的手法ではなく逐次近似画像再構成法が用いられるが、共役勾配法といった通常の逐次近似画像再構成法では、再構成画像内での空間分解能や雑音特性のばらつきが大きい不均一な画像が生成される問題がある。
画質の不均一性の原因は、拡散光トモグラフィにおいては測定データに含まれる情報量が被測定部位の内部の位置によって大きく異なるためである。すなわち、時間分解計測法を用いた拡散光トモグラフィーでは、被計測部位表面の光入射端から入射された光子は被計測部位内で散乱を繰り返しながら飛行し、被計測部位表面の光検出端に到達して検出される。この入射から検出までの光子の飛行時間が短いほど、飛行距離が短くなり飛行経路も限られる。逆に、飛行時間が長いほど、飛行距離が長くなり飛行経路は限定されない。このため、測定データに含まれる情報量は、光子の飛行時間に応じて変化する。典型的には、光子の飛行時間が短い場合には、被計測部位表面に近い飛行経路の割合が多くなり、また、飛行時間が長い場合には、被計測部位表面から遠い（被計測部位の中央付近の）飛行経路の割合が多くなる。そのため、被計測部位の中央付近の情報量は表面付近と比べて少なくなってしまう。これにより、再構成画像の周辺部における空間分解能及び雑音が再構成画像の中央付近における空間分解能及び雑音と比べて大きくなるといった、再構成画像内での空間分解能及び雑音のばらつきが生じてしまう。
本実施形態の生体計測装置１０及び画像作成方法では、再構成画像の各画素毎に設定されるＪ個の係数ｗ1〜ｗJを用いて、画像再構成のための逐次近似演算を行う。そして、本実施形態のように、これらの係数ｗ1〜ｗJを、Ｎ回の反復演算において収束速度が最も遅い部分領域に当該画素の収束速度を一致させるように設定することで、収束速度が均一化され、被計測部位内部の位置による空間分解能や雑音特性の差を抑制して、より均一に近い画像を作成することが可能となる。
ここで、本実施形態に係る画像作成方法による上記効果を確認するために行われたシミュレーションの結果について説明する。図３（ａ）及び図３（ｂ）は、このシミュレーションにおいて再構成の対象となる２種類の被測定部位を示す図であり、各図には、例えば腫瘍といった複数の光吸収体Ｄが示されている。なお、図３（ａ）に示される３つの光吸収体Ｄの中心座標（ｘ，ｙ）は、それぞれ（２１，２０）、（６４，１０６）、及び（７５，５１）である。また、図３（ｂ）に示される４つの光吸収体Ｄの中心座標（ｘ，ｙ）は、それぞれ（５６，７９）、（６６，２４）、（１０９，４９）、及び（１０９，１１０）である。
このシミュレーションでは、次のステップ（１）〜（３）を行った。
（１）各画素の収束速度を評価する為に、画像内において均等に配置されたスポットにおける画像再構成を行い、ステップサイズｗ1〜ｗJを決定した。
（２）得られたステップサイズｗ1〜ｗJを用いて画像の再構成を行い、再構成画像の全体で空間分解能が略均一となっていることを確認した。
（３）背景画像の再構成による領域間での雑音伝播を評価した。
なお、上記ステップ（１）〜（３）に共通する条件は、次のとおりである。まず、被計測部位の内部における光子の飛行をシミュレーションする順問題解析（Forward problem）と、検出された光子ヒストグラムから画像を再構成する逆問題解析（Inverse problem）という２種類の解析における条件を、それぞれ図４及び図５に示す。図４及び図５に示されるように、順問題解析と逆問題解析との間でグリッドサイズや画像サイズが異なるのは、連続系を離散系に変換したことによる。光検出部において検出される光子は、必ず離散的なデータとして取得される。一方、数値シミュレーションの値は連続関数に近い形で計算される。したがって、本シミュレーションでは、ダウンサンプリングすることにより連続系を離散系に変換している。
（１）ステップサイズｗ1〜ｗJの決定
測定データに基づく１３２行１３２列の画素を含む画像を１２行１２列の部分領域に分割し、各部分領域の中央にホットスポットを配置した。本シミュレーションでは、図６（ａ）〜図６（ｄ）に示される４つの画像を再構成することにより、各部分領域のＣＲＣを求め、前述した方法によって全画素のステップサイズｗ1〜ｗJを決定した。
（２）画像の再構成
次に、ステップ（１）により決定されたステップサイズｗ1〜ｗJを用い、式（１４）に示された反復演算式を演算することにより、画像の再構成を行った。なお、このステップでは、各部分領域のＣＲＣに差があると画像中にエッジとなって現れてしまうので、ステップサイズｗ1〜ｗJを各画素値とする画像（ステップサイズ画像）を平均値フィルターを用いてぼかし、滑らかに変化するステップサイズ画像を得た。本シミュレーションでは、カーネルサイズ９×９の平均値フィルターを用いた。
図７（ａ）は、図３（ａ）に示された画像に対応する、本シミュレーションによる再構成後の画像である。図７（ｂ）は、図３（ａ）に示された画像に対応する、ステップサイズｗ1〜ｗJを用いない従来の方法による再構成後の画像である。図８（ａ）〜図８（ｃ）は、図９に示されるように光吸収体Ｄを通る３本のラインＬ１〜Ｌ３上の画素値の変化を示すグラフであって、グラフＧ１１は本シミュレーションによる画素値の変化を示し、グラフＧ１２は従来の方法による画素値の変化を示している。また、図１０（ａ）は、図３（ｂ）に示された画像に対応する、本シミュレーションによる再構成後の画像である。図１０（ｂ）は、図３（ｂ）に示された画像に対応する、ステップサイズｗ1〜ｗJを用いない従来の方法による再構成後の画像である。図１１（ａ）及び図１１（ｂ）は、図１２に示されるように光吸収体Ｄを通る２本のラインＬ４及びＬ５上の画素値の変化を示すグラフであって、グラフＧ２１は本シミュレーションによる画素値の変化を示し、グラフＧ２２は従来の方法による画素値の変化を示している。
図７（ｂ）を参照すると、従来の再構成画像では、左上といった画像周辺部の光吸収体Ｄの画素値が比較的高くなっており、画像中央付近の光吸収体Ｄの画素値と大きく異なっていることがわかる。これに対し、ステップサイズｗ1〜ｗJを用いた本実施形態による再構成画像では、図７（ａ）に示されるように、各光吸収体Ｄの画素値が、画像の周辺部及び中央付近においてほぼ均一な値となっている。また、図８（ｄ）及び図８（ｅ）を参照すると、従来の再構成画像では、画像周辺部における画素値の変化が大きい。これに対し、図８（ａ）〜図８（ｃ）に示されるように、本実施形態による再構成画像では画像周辺部における画素値の変化が小さくなっている。これは、本実施形態において画像中央付近ではステップサイズｗjが大きくなって反復演算における更新量が大きくなり、画像周辺部ではステップサイズｗjが小さくなって更新量が小さくなったことによる。以上の結果から、本実施形態の生体計測装置１０及び画像作成方法によれば、収束速度が均一化されたことにより、被計測部位内部の位置による空間分解能の差を抑制して、より均一に近い画像を作成できることがわかる。
次に、統計雑音を付加した場合のシミュレーション結果を示す。このシミュレーションでは、各検出光子ヒストグラムの最大値が一定値（例えば５０）となるように光子ヒストグラムをフィッティングした後、この光子ヒストグラムにポアソン雑音を付加した。図１３（ａ）は、図３（ａ）に示された画像に対応する、統計雑音を付加した場合の本実施形態による再構成後の画像である。図１３（ｂ）は、図３（ａ）に示された画像に対応する、統計雑音を付加した場合の従来の方法による再構成後の画像である。また、図１４（ａ）は、図３（ｂ）に示された画像に対応する、統計雑音を付加した場合の本実施形態による再構成後の画像である。図１４（ｂ）は、図３（ｂ）に示された画像に対応する、統計雑音を付加した場合の従来の方法による再構成後の画像である。
図１３（ｂ）及び図１４（ｂ）を参照すると、統計雑音が存在する場合、ステップサイズｗ1〜ｗJを用いない従来の方法では、画像周辺部の空間分解能が高くなって雑音の影響が大きく出ているが、画像中央付近では空間分解能が高くなく雑音の影響が小さいことがわかる。これに対し、図１３（ａ）及び図１４（ａ）を参照すると、ステップサイズｗ1〜ｗJを用いた本実施形態による方法では、反復演算における画像周辺部の更新量が抑えられているので、雑音の影響が顕著に抑制され、左上に存在する光吸収体Ｄの形を明確に視認できる。一方、画像中央付近では、雑音の影響が従来方法と同程度に抑えられている。なお、画像中央付近の統計雑音が画像周辺部の統計雑音よりも大きいのは、統計雑音の有無によって分解能の向上速度に変化が生じるためであり、統計雑音がない場合を仮定して決定したステップサイズを用いたことに起因する現象である。
このように、計測データに統計雑音が含まれる場合、本実施形態の生体測定方法では従来方法と比較して画素値に変化はあるものの、画像全体の統計雑音の形態的特徴に大きな変化がないことがわかる。すなわち、本実施形態の生体測定方法によれば、画像周辺部において選択的に、分解能の向上を遅くすることができていることがわかる。
続いて、光吸収体Ｄを配置しない背景画像に関する雑音評価について説明する。ここでは、光吸収体Ｄが存在しない計測データにポアソン雑音を付加した。このときの再構成画像は、雑音によって歪んだ画像となる。この再構成画像を６行６列の３６個の部分領域に分割し、左上から右下へ順に番号を付与した。その後、それぞれの部分領域を雑音指標（ＮＳＤ）によって評価した。なお、ＮＳＤの計算式は次式（１７）である。また、図１５は、部分領域毎のＮＳＤを示すグラフであって、横軸は部分領域の番号を示し、縦軸はＮＳＤの値を示している。なお、添字ＳＤは標準偏差であることを示す。

また、部分領域毎のＮＳＤを母集団とする標準偏差を算出し、この標準偏差を部分領域間の雑音の不均一性を評価する値として用いる。図１６は、標準偏差の算出結果を示す図表である。
図１５を参照すると、従来方法では部分領域間のＮＳＤのばらつきが大きい。これに対し、本実施形態の方法では、ＮＳＤのばらつきが全体的に抑えられていることがわかる。このことから、本実施形態の方法では、空間分解能を均一に近づけるために用いられるステップサイズｗjが、雑音伝播の不均一をも緩和していることがわかる。更に、図１６に示されるように、雑音の不均一性についても、本実施形態の方法では従来方法よりも小さくなる。
本発明による生体計測装置および画像作成方法は、上述した実施形態に限られるものではなく、他に様々な変形が可能である。例えば、上記実施形態では被検者の被計測部位に光を照射し、被計測部位からの拡散光を検出し、該検出信号に基づいて被計測部位内における光吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成している。しかしながら、本発明は、光を用いた生体計測装置や画像作成方法に限らず、放射線や音波を用いた生体計測装置や画像作成方法にも適用可能である。すなわち、本発明に係る生体計測装置は、被検者の被計測部位に放射線若しくは音波を照射する照射部と、被計測部位からの拡散した放射線若しくは音波を検出する検出部と、検出部からの出力信号に基づいて、被計測部位内における放射線若しくは音波の吸収係数分布を演算し、被計測部位の内部に関する再構成画像を作成する演算部とを備えてもよい。この場合、演算部は、上述した実施形態と同様の手法により、再構成画像を作成するとよい。また、本発明に係る画像作成方法は、被検者の被計測部位に放射線若しくは音波を照射し、被計測部位からの拡散した放射線若しくは音波を検出し、該検出信号に基づいて被計測部位内における放射線若しくは音波の吸収係数分布を演算し、被計測部位の内部に関する再構成画像を上述した実施形態と同様の手法により作成する方法であってもよい。ここで、放射線とは例えばＸ線、γ線、若しくはマイクロ波といった短波長の電磁波をいい、音波とは例えば超音波のような波動をいう。
１０…生体計測装置、１４…演算部、１６…光出射／検出端、２２…光源、２４…光スイッチ、２６…光源用光ファイバ、２８…出射用光ファイバ、３０…光検出器、３２…シャッター、３４…検出用光ファイバ、３６…信号処理回路、３８…表示部。

Claims (8)

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

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

    （但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘj (k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄj (k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより前記再構成画像を作成することを特徴とする、生体計測装置。
  3. 前記演算部は、
    前記再構成画像が分割されて成り複数の画素をそれぞれ含む複数の部分領域のうち、Ｎ回の反復演算において収束速度が最も遅い部分領域の前記収束速度（以下、最低収束速度という）ＣNを求め、
    ０＜ｖm＜１を満たすＭ個の数値ｖm（但し、ｍは１からＭまでの整数）を用意し、
    １からＭまでの各ｍについて、次の反復式

    を用いてＮ回の反復演算を行うことによりＪ個の各画素の画素値ｘ1 (N)〜ｘJ (N)を算出し、
    該画素値ｘ1 (N)〜ｘJ (N)から得られる各部分領域の前記収束速度が前記最低収束速度ＣNと略一致したときの数値ｖmを、該部分領域に含まれる複数の画素の前記係数ｗjとすることを特徴とする、請求項１または２に記載の生体計測装置。
  4. 前記演算部は、Ｎ回の反復演算において収束速度が最も遅い部分領域の前記係数ｗjを１とすることを特徴とする、請求項３に記載の生体計測装置。
  5. 被検者の被計測部位に光を照射し、前記被計測部位からの拡散光を検出し、該検出信号に基づいて前記被計測部位内の光吸収係数分布を演算し、前記被計測部位の内部に関する再構成画像を作成する方法であって、
    前記再構成画像の各画素毎に設定される、０より大きく且つ１以下であるＪ個の係数ｗj（但し、添字ｊは１からＪまでの整数であり、Ｊは前記再構成画像の画素数である。）を算出し、
    次の反復式

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

    （但し、ｋは１からＮまでの整数であり、Ｎは反復演算の回数である。ｘj (k)は第ｊ番目の画素のｋ回目の反復演算時における画素値であり、ｄj (k)は第ｊ番目の画素のｋ回目の反復演算時における更新量である。）を用いた逐次近似演算を行うことにより前記再構成画像を作成することを特徴とする、画像作成方法。
  7. 前記再構成画像が分割されて成り複数の画素をそれぞれ含む複数の部分領域のうち、Ｎ回の反復演算において収束速度が最も遅い部分領域の前記収束速度（以下、最低収束速度という）ＣNを求め、
    ０＜ｖm＜１を満たすＭ個の数値ｖm（但し、ｍは１からＭまでの整数）を用意し、
    １からＭまでの各ｍについて、次の反復式

    を用いてＮ回の反復演算を行うことによりＪ個の各画素の画素値ｘ1 (N)〜ｘJ (N)を算出し、
    該画素値ｘ1 (N)〜ｘJ (N)から得られる各部分領域の前記収束速度が前記最低収束速度ＣNと略一致したときの数値ｖmを、該部分領域に含まれる複数の画素の前記係数ｗjとすることを特徴とする、請求項５または６に記載の画像作成方法。
  8. 前記演算ステップにおいて、Ｎ回の反復演算において収束速度が最も遅い部分領域の前記係数ｗjを１とすることを特徴とする、請求項７に記載の画像作成方法。
JP2011151086A 2011-07-07 2011-07-07 生体計測装置および画像作成方法 Active JP5782314B2 (ja)

Priority Applications (5)

Application Number Priority Date Filing Date Title
JP2011151086A JP5782314B2 (ja) 2011-07-07 2011-07-07 生体計測装置および画像作成方法
US14/130,690 US9282932B2 (en) 2011-07-07 2012-06-28 Biometric apparatus and image-generating method
PCT/JP2012/066552 WO2013005635A1 (ja) 2011-07-07 2012-06-28 生体計測装置および画像作成方法
CN201280033086.7A CN103635148B (zh) 2011-07-07 2012-06-28 生物体测量装置以及图像制作方法
EP12807701.3A EP2730231B1 (en) 2011-07-07 2012-06-28 Biometric apparatus and image-generating method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2011151086A JP5782314B2 (ja) 2011-07-07 2011-07-07 生体計測装置および画像作成方法

Publications (2)

Publication Number Publication Date
JP2013019696A JP2013019696A (ja) 2013-01-31
JP5782314B2 true JP5782314B2 (ja) 2015-09-24

Family

ID=47436990

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2011151086A Active JP5782314B2 (ja) 2011-07-07 2011-07-07 生体計測装置および画像作成方法

Country Status (5)

Country Link
US (1) US9282932B2 (ja)
EP (1) EP2730231B1 (ja)
JP (1) JP5782314B2 (ja)
CN (1) CN103635148B (ja)
WO (1) WO2013005635A1 (ja)

Families Citing this family (20)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6595211B2 (ja) * 2015-05-11 2019-10-23 キヤノンメディカルシステムズ株式会社 核医学診断装置及び核医学画像処理装置
CN104997533B (zh) * 2015-06-23 2017-11-14 武汉超信电子工程有限公司 超声探头几何参数自动校正方法和装置
US9730649B1 (en) 2016-09-13 2017-08-15 Open Water Internet Inc. Optical imaging of diffuse medium
WO2018189815A1 (ja) * 2017-04-11 2018-10-18 オリンパス株式会社 光学物性値分布の推定方法、光学物性値分布の推定プログラム、及び光学物性値分布推定装置
US10778911B2 (en) 2018-03-31 2020-09-15 Open Water Internet Inc. Optical transformation device for imaging
US10506181B2 (en) 2018-03-31 2019-12-10 Open Water Internet Inc. Device for optical imaging
US10778912B2 (en) 2018-03-31 2020-09-15 Open Water Internet Inc. System and device for optical transformation
CN108398245A (zh) * 2018-05-28 2018-08-14 南昌大学 一种同心发散光源的检测方法与装置
US10966612B2 (en) 2018-06-14 2021-04-06 Open Water Internet Inc. Expanding beam optical element
US10962929B2 (en) 2018-09-14 2021-03-30 Open Water Internet Inc. Interference optics for optical imaging device
US10874370B2 (en) 2019-01-28 2020-12-29 Open Water Internet Inc. Pulse measurement in optical imaging
US10955406B2 (en) 2019-02-05 2021-03-23 Open Water Internet Inc. Diffuse optical imaging with multiple beams
US11320370B2 (en) 2019-06-26 2022-05-03 Open Water Internet Inc. Apparatus for directing optical and acoustic signals
US11581696B2 (en) 2019-08-14 2023-02-14 Open Water Internet Inc. Multi-channel laser
US11622686B2 (en) 2019-11-22 2023-04-11 Open Water Internet, Inc. Optical imaging with unshifted reference beam
CN111332287B (zh) * 2020-03-11 2021-01-22 合肥鑫晟光电科技有限公司 一种驾驶辅助方法及系统
US20210330266A1 (en) * 2020-04-24 2021-10-28 Hi Llc Systems and Methods for Noise Removal in an Optical Measurement System
US11819318B2 (en) 2020-04-27 2023-11-21 Open Water Internet Inc. Optical imaging from light coherence
US11559208B2 (en) 2020-05-19 2023-01-24 Open Water Internet Inc. Imaging with scattering layer
US11259706B2 (en) 2020-05-19 2022-03-01 Open Water Internet Inc. Dual wavelength imaging and out of sample optical imaging

Family Cites Families (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4662222A (en) * 1984-12-21 1987-05-05 Johnson Steven A Apparatus and method for acoustic imaging using inverse scattering techniques
US5694938A (en) 1995-06-07 1997-12-09 The Regents Of The University Of California Methodology and apparatus for diffuse photon mimaging
US5905261A (en) * 1995-12-01 1999-05-18 Schotland; John Carl Imaging system and method using direct reconstruction of scattered radiation
JP4364995B2 (ja) 2000-03-21 2009-11-18 浜松ホトニクス株式会社 散乱吸収体内部の光路分布計算方法
CA2581592C (en) 2004-09-24 2014-04-15 New Art Advanced Research Technologies Inc. Method for fluorescence tomographic imaging
CN101312682B (zh) 2005-11-23 2011-05-25 皇家飞利浦电子股份有限公司 用于对混浊介质进行光学成像的方法和装置
US7983740B2 (en) * 2006-12-22 2011-07-19 Washington University High performance imaging system for diffuse optical tomography and associated method of use
WO2009077947A2 (en) 2007-12-17 2009-06-25 Koninklijke Philips Electronics N.V. Method for detecting the presence of inhomogeneities in an interior of a turbid medium and device for imaging the interior of turbid media
EP2239560B1 (en) 2007-12-27 2020-04-22 Omron Corporation X-ray examining apparatus and x-ray examining method
TWI394490B (zh) 2008-09-10 2013-04-21 Omron Tateisi Electronics Co Ｘ射線檢查裝置及ｘ射線檢查方法
JP2010094500A (ja) 2008-09-19 2010-04-30 Canon Inc 測定装置及び測定方法
CN101509871B (zh) 2009-03-25 2011-09-28 华中科技大学 适于小动物的荧光分子层析成像方法
JPWO2010143421A1 (ja) 2009-06-10 2012-11-22 パナソニック株式会社 光融合型イメージング方法、光融合型イメージング装置、プログラムおよび集積回路
JP5600946B2 (ja) 2010-01-28 2014-10-08 株式会社島津製作所 断層撮影装置

Also Published As

Publication number Publication date
EP2730231B1 (en) 2017-11-29
WO2013005635A1 (ja) 2013-01-10
CN103635148A (zh) 2014-03-12
US20140135620A1 (en) 2014-05-15
EP2730231A4 (en) 2015-03-25
EP2730231A1 (en) 2014-05-14
JP2013019696A (ja) 2013-01-31
CN103635148B (zh) 2015-11-25
US9282932B2 (en) 2016-03-15

Similar Documents

Publication Publication Date Title
JP5782314B2 (ja) 生体計測装置および画像作成方法
JP5975437B2 (ja) 生体計測装置の計測データ選択方法、生体計測装置の光出射位置決定方法、および生体計測装置
Yao et al. Direct approach to compute Jacobians for diffuse optical tomography using perturbation Monte Carlo-based photon “replay”
US8817257B2 (en) Method for reconstructing the optical properties of a medium using a combination of a plurality of mellin-laplace transforms of a magnitude comprising a time distribution of a received signal, and associated reconstruction system
JP5570674B2 (ja) 物体観測装置、物体観測方法、および記録媒体
US20070244395A1 (en) Systems and methods for multi-spectral bioluminescence tomography
WO2012153769A1 (ja) 光トモグラフィ装置
US20130197342A1 (en) Object information acquiring apparatus
Zuo et al. Spectral crosstalk in photoacoustic computed tomography
CN103815929B (zh) 被检体信息获取装置
JP5924658B2 (ja) 濃度定量装置、光吸収係数算出方法、等価散乱係数算出方法、濃度定量方法、光吸収係数の算出を行うプログラム及び濃度の算出を行うプログラム
JP5658979B2 (ja) 生体光計測装置及び生体における吸収係数の変動を推定する方法
JP5420163B2 (ja) 生体計測装置
CN106455994B (zh) 光声装置
Aspri et al. Mathematical and numerical challenges in diffuse optical tomography inverse problems
Kalyanov et al. Time-multiplexing approach for fast time-domain near-infrared optical tomography combined with neural-network-enhanced image reconstruction
Wong et al. Objective assessment and design improvement of a staring, sparse transducer array by the spatial crosstalk matrix for 3d photoacoustic tomography
Okawa et al. Reduction of Poisson noise in measured time-resolved data for time-domain diffuse optical tomography
Żołek et al. Analysis of estimation of optical properties of sub superficial structures in multi layered tissue model using distribution function method
Soonthornsaratoon Gradient-based methods for quantitative photoacoustic tomography
Lin et al. Anatomical Modeling and Optimization of Speckle Contrast Optical Tomography
Munoz-Barrutia et al. Sparse algebraic reconstruction for fluorescence mediated tomography
JP5834704B2 (ja) 濃度定量装置、光吸収係数算出方法、濃度定量方法、光吸収係数の算出を行うプログラム及び濃度の算出を行うプログラム
WO2018189815A1 (ja) 光学物性値分布の推定方法、光学物性値分布の推定プログラム、及び光学物性値分布推定装置
Paulsen et al. Near infrared spectroscopic imaging: theory

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20140702

A711 Notification of change in applicant

Free format text: JAPANESE INTERMEDIATE CODE: A711

Effective date: 20150108

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A821

Effective date: 20150108

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20150714

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20150717

R150 Certificate of patent or registration of utility model

Ref document number: 5782314

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150