JP6212312B2 - 物体内部推定装置およびその方法 - Google Patents

物体内部推定装置およびその方法 Download PDF

Info

Publication number
JP6212312B2
JP6212312B2 JP2013146966A JP2013146966A JP6212312B2 JP 6212312 B2 JP6212312 B2 JP 6212312B2 JP 2013146966 A JP2013146966 A JP 2013146966A JP 2013146966 A JP2013146966 A JP 2013146966A JP 6212312 B2 JP6212312 B2 JP 6212312B2
Authority
JP
Japan
Prior art keywords
voting
light
unit
value
estimation
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
JP2013146966A
Other languages
English (en)
Other versions
JP2014055939A (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.)
Panasonic Corp
Panasonic Holdings Corp
Original Assignee
Panasonic Corp
Matsushita Electric Industrial Co Ltd
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 Panasonic Corp, Matsushita Electric Industrial Co Ltd filed Critical Panasonic Corp
Priority to JP2013146966A priority Critical patent/JP6212312B2/ja
Publication of JP2014055939A publication Critical patent/JP2014055939A/ja
Application granted granted Critical
Publication of JP6212312B2 publication Critical patent/JP6212312B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B11/00—Measuring arrangements characterised by the use of optical techniques
    • G01B11/24—Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)
  • Investigating Materials By The Use Of Optical Means Adapted For Particular Applications (AREA)

Description

本発明は、物体の内部を推定する物体内部推定装置およびその方法に関する。
物体の内部を推定する技術は重要である。この技術は、例えば、食品に混入された釘などの異物を検知すること、あるいは、皮膚に刺さった針などの異物の形状を推定することに応用される。高精度に物体の内部を推定できれば、食の安全性の向上や医療の高度化の実現をサポートすることができる。
物体の内部を推定する方法として、Ｘ線ＣＴ（Ｘ−ｒａｙ Ｃｏｍｐｕｔｅｄ Ｔｏｍｏｇｒａｐｈｙ）に基づく方法が幅広く利用されている。Ｘ線ＣＴを用いる方法では、基本的に、異物が混入された物体にＸ線を入射して、物体からのＸ線の出射量を観測する。そして、Ｘ線の入射量に対する出射量の割合（すなわち光の吸収量）に基づいて、物体内における光（Ｘ線）の吸収率の分布を推定する。実際には、３６０度方向から物体にＸ線を入射したときの、物体から出射されるＸ線の減衰量が観測され、逆投影法、あるいは、フィルタ逆投影法と呼ばれる手法を用いて物体の内部の吸収率分布が推定される。
このように線源から光線を入射して物体内を透過した光線の出射量を観測し、観測結果を解析することで、物体内に含まれる吸収物体の位置および形状を推定することができる。
特開２０００−１０５１９１号公報 特開平８−２９３２９号公報 特開２００２−３３３４００号公報
「Ｃ言語による画像再構成の基礎」、橋下、篠原、医療科学社、Ｐ２４５ 「ＭＣＭＬ−Ｍｏｎｔｅ Ｃａｒｌｏ ｍｏｄｅｌｉｎｇ ｏｆ ｌｉｇｈｔ ｔｒａｎｓｐｏｒｔ ｉｎ ｍｕｌｔｉ−ｌａｙｅｒｄ ｔｉｓｓｕｅｓ」、Ｗａｎｇ ｅｔ ａｌ．、Ｃｏｍｐｕｔｅｒ Ｍｅｔｈｏｄｓ ａｎｄ Ｐｒｏｇｒａｍｓ ｉｎ Ｂｉｏｍｅｄｉｃｉｎｅｓ ４７ （１９９５） ｐｐ．１３１−１４６． 「Ａｃｑｕｉｒｉｎｇ Ｓｃａｔｔｅｒｉｎｇ Ｐｒｏｐｅｒｔｉｅｓ ｏｆ Ｐａｒｔｉｃｉｐａｔｉｎｇ Ｍｅｄｉａ ｂｙ Ｄｉｌｕｔｉｏｎ」、Ｎａｒａｓｉｍｈａｍ ｅｔ ａｌ．、ＡＣＭ ＳＩＧＧＲＡＰＨ、２００６ 「Ｎｅａｒ−ｉｎｆｒａｒｅｄ ｏｐｔｉｃａｌ ｐｒｏｐｅｒｔｉｅｓ ｏｆ ｅｘ ｖｉｖｏ ｈｕｍａｎ ｓｋｉｎ ａｎｄｓｕｂｃｕｔａｎｅｏｕｓ ｔｉｓｓｕｅｓ ｍｅａｓｕｒｅｄ ｕｓｉｎｇ ｔｈｅ ＭｏｎｔｅＣａｒｌｏ ｉｎｖｅｒｓｉｏｎ ｔｅｃｈｎｉｑｕｅ」、Ｓｉｍｐｓｏｎ ｅｔ ａｌ．、Ｐｈｙｓ． Ｍｅｄ． Ｂｉｏｌ．４３、１９９８、ｐｐ．２４６５−２４７８
しかしながら、非特許文献１に示される方法では、物体内部を推定するための機械装置が大きくなり、出射線の計測に特殊な装置が必要であるという課題がある。また、Ｘ線を利用して人体の内部を推定する場合、いわゆる「医療被曝」のリスクがある。したがって、長時間のＸ線の利用や妊娠中の人に対してのＸ線の利用は避けるべきである。
近年、人体に安全な光である近赤外光を利用して物体内部の光吸収物体の位置を推定する技術が開発されている。近赤外光を用いる場合、安全性、装置の手軽さの面ではＸ線を用いる場合よりも優位である。しかしながら、近赤外光は、物体内で散乱するため、Ｘ線に比べて、出射光の解析が難しいという問題がある。
そこで、本発明は、人体に安全な光を用いて物体の内部を推定する場合に推定精度を向上させることができる物体内部推定装置およびその方法を提供する。
本発明の一態様に係る物体内部推定装置は、第１の物体の内部を推定する物体内部推定装置であって、前記第１の物体とは異なる第２の物体が前記第１の物体内に含まれていないと仮定した場合における、（i）光が前記第１の物体に入射してから出射するまでの間にとる複数の経路を表す複数の光路モデルと、（ii）前記複数の光路モデルに沿って前記第１の物体内を伝搬して前記第１の物体から出射した仮想的な光の所定面上における強度を示す第１の光強度分布とを取得する第１の取得部と、前記第１の物体に入射してから出射した光の強度分布であって前記所定面上において観測された光の強度分布である第２の光強度分布を取得する第２の取得部と、前記所定面上の位置ごとに、前記第１の光強度分布および前記第２の光強度分布間における光の強度の比較結果に基づいて、前記第１の物体内に前記第２の物体が存在しない非存在確率を示す投票値を算出する投票値算出部と、前記第１の物体の内部を複数の投票領域に仮想的に分割し、前記所定面上の位置ごとに、当該位置に到達する光の光路モデル上の投票領域に対して当該位置の前記投票値を投票する投票部と、前記複数の投票領域の各々に対して、前記投票部で得られた投票結果が閾値より小さい場合に、前記第２の物体が存在すると推定する推定部とを備える。
なお、これらの包括的または具体的な態様は、システム、方法、集積回路、コンピュータプログラムまたはコンピュータ読み取り可能なＣＤ−ＲＯＭなどの記録媒体で実現されてもよく、システム、方法、集積回路、コンピュータプログラムおよび記録媒体の任意な組み合わせで実現されてもよい。
本発明の一態様に係る物体内部推定装置によれば、人体に安全な光を用いて物体の内部を推定する場合に推定精度を向上させることができる。
実施の形態１に係る物体内部推定装置の機能構成を示すブロック図 実施の形態１に係る物体内部推定装置の処理動作を示すフローチャート 実施の形態１における光路モデルの導出処理の詳細を示すフローチャート 実施の形態１における光強度の減衰量を説明するための図 実施の形態１における光源部、散乱体、および撮像部の配置の一例を示す模式図 実施の形態１における観測画像の一例を示す図 実施の形態１における推定対象断面に対応する画素列の値の一例を示す図 実施の形態１における光路モデルの例を示す図 実施の形態１における光透過率の分布（第１の光強度分布）と観測値の分布（第２の光強度分布）の一例を示す図 実施の形態１における投票処理を説明するための図 実施の形態１における投票処理を説明するための図 実施の形態２に係る物体内部推定装置の機能構成を示すブロック図 実施の形態２に係る物体内部推定装置の処理動作を示すフローチャート 実施の形態２における複数の光路モデルを示す模式図 実施の形態２における補正処理を説明するための図 従来の散乱体内の吸収率分布推定方法を示すフローチャート
（本発明の基礎となった知見）
本発明者は、先行技術文献に記載の物体内部推定方法に関し、以下の問題が生じることを見いだした。なお、以下説明において、内部で散乱が生じる物体を散乱体と呼ぶ。また、散乱体内で光を吸収する物体を光吸収物体と呼ぶ。
図１６は、従来の散乱体内の吸収率分布推定方法を示すフローチャートである。図１６に示すように、特許文献１〜３では、複数の波長の入射光および出射光の強度を計測し、入射光と出射光の計測された強度の組合せから行列Ｙを生成する。その後、散乱体を複数の空間（以下、ボクセルと呼ぶ）に分割し、モンテカルロ・レイトレーシングで光の経路をシミュレーションする。シミュレーション結果が示す光の経路（光路モデル）から平均光路長を求める。複数の波長の光を利用することで未知数の削減を行い、行列Ｙと平均光路長とから散乱体中の吸収率分布を求める。
しかしながら、これらの方法では、散乱体内の吸収率分布を線形演算で算出するために、以下に示す２つの近似が行われている。第一に、散乱や散乱の広がりによって生じる光量の減衰が、複数の波長で等しいという近似である。第二に、異なる波長間で平均光路長が等しいという近似である。実際には、散乱と吸収係数の分布は波長依存性があるため、平均光路長は波長ごとに異なる。このような近似により、上記従来の方法では、散乱体内の正確な吸収率分布を求めることが難しい。
そこで、本発明の一態様に係る物体内部推定装置は、第１の物体の内部を推定する物体内部推定装置であって、前記第１の物体とは異なる第２の物体が前記第１の物体内に含まれていないと仮定した場合における、（i）光が前記第１の物体に入射してから出射するまでの間にとる複数の経路を表す複数の光路モデルと、（ii）前記複数の光路モデルに沿って前記第１の物体内を伝搬して前記第１の物体から出射した仮想的な光の所定面上における強度を示す第１の光強度分布とを取得する第１の取得部と、前記第１の物体に入射してから出射した光の強度分布であって前記所定面上において観測された光の強度分布である第２の光強度分布を取得する第２の取得部と、前記所定面上の位置ごとに、前記第１の光強度分布および前記第２の光強度分布間における光の強度の比較結果に基づいて、前記第１の物体内に前記第２の物体が存在する確率を示す投票値を算出する投票値算出部と、前記第１の物体の内部を複数の投票領域に仮想的に分割し、前記所定面上の位置ごとに、当該位置に到達する光の光路モデル上の投票領域に対して当該位置の前記投票値を投票する投票部と、投票結果に基づいて、前記複数の投票領域の各々に前記第２の物体が存在するか否かを推定する推定部とを備える。
この構成によれば、光が第１の物体に入射してから出射するまでの間にとる複数の経路を表す複数の光路モデルを用いて、第１の物体内の各領域に第２の物体が存在する否かを推定することができる。したがって、第１の物体内で光が散乱する場合でも、第１の物体内の第２の物体の位置および形状を推定することができる。つまり、近赤外光などＸ線以外の人体に安全な光を用いて物体の内部を推定する場合に推定精度を向上させることができる。
さらに、この構成によれば、第１の光強度分布と第２の光強度分布間における光の強度の比較結果に基づいて、第１の物体内に第２の物体が存在する確率を示す投票値を算出することができる。また、光路モデル上の投票領域に対する投票値の投票結果に基づいて、複数の投票領域の各々に第２の物体が存在するか否かを推定することができる。つまり、第１の物体の内部を推定するために、複数の波長の光を用いる必要がなく、線形演算の必要もない。したがって、上述したような吸収率分布の線形演算のための近似も不要となるので、物体の内部の推定精度を向上させることができる。
また例えば、前記物体内部推定装置は、さらに、前記複数の投票領域に対する投票回数のばらつきに基づいて、前記投票結果を補正する補正部を備え、前記推定部は、補正された前記投票結果に基づいて、前記複数の投票領域の各々に前記第２の物体が存在するか否かを推定してもよい。
この構成によれば、複数の投票領域に対する投票回数のばらつきに基づいて投票結果を補正することができる。したがって、例えば、第１の物体内の光が通過しにくい領域（投票回数が少ない領域）に第２の物体が存在すると誤って推定されることを減少させることができる。
また例えば、前記複数の投票領域は、前記第１の物体の断面上の複数の領域であり、前記推定部は、複数の断面の各々について、前記複数の投票領域の各々に前記第２の物体が存在するか否かを推定することにより、複数の２次元断層画像を生成し、前記物体内部推定装置は、さらに、前記複数の２次元断層画像を統合することにより、前記第１の物体内の３次元ボリュームデータを生成する統合部を備えてもよい。
この構成によれば、複数の２次元断層画像を統合することにより、３次元ボリュームデータを生成することができる。したがって、より正確に第１の物体の内部を推定することができる。
また例えば、前記第２の取得部は、可視光または近赤外光を前記第１の物体に照射する光源部と、前記第１の物体に入射してから出射した光を撮像することにより、前記第２の光強度分布を示す観測画像を取得する撮像部とを備えてもよい。
この構成によれば、第２の取得部は、光源部と撮像部とを用いて、第２の光強度分布を示す観測画像を取得することができる。
また例えば、前記第１の取得部は、前記複数の光路モデルを、モンテカルロ・レイトレーシング・シミュレーションによって導出してもよい。
この構成によれば、モンテカルロ・レイトレーシング・シミュレーションによって、複数の光路モデルを比較的容易に導出することができる。
また例えば、前記第１の取得部は、シミュレーション結果と観測結果との誤差を粗密探索によって削減することにより、前記モンテカルロ・レイトレーシング・シミュレーションで用いるパラメータを導出してもよい。
この構成によれば、モンテカルロ・レイトレーシング・シミュレーションで用いるパラメータを粗密探索によって導出することができ、パラメータの精度を向上させることができる。
また例えば、前記第１の取得部は、前記第１の物体に互いに異なる方向から入射する複数の光の各々について、前記複数の光路モデルと前記第１の光強度分布とを取得し、前記第２の取得部は、前記複数の光の各々について前記第２の光強度分布を取得し、前記投票値算出部は、前記複数の光の各々について前記投票値を算出し、前記投票部は、前記複数の光の各々について前記投票値を投票してもよい。
この構成によれば、第１の物体に複数の方向から入射する複数の光の各々について、投票結果を得ることができる。したがって、１つの方向から入射する光のみを用いる場合よりも、奥行き方向の推定精度を向上させることができる。
また例えば、前記物体内部推定装置は、さらに、前記複数の投票領域の各々の位置に基づいて前記投票結果を補正する補正部を備え、前記推定部は、補正された前記投票結果に基づいて、前記複数の投票領域の各々に前記第２の物体が存在するか否かを推定してもよい。
この構成によれば、各投票領域の位置に基づいて投票結果を補正することができる。したがって、投票領域の位置に依存して第２の物体が存在するか否かの推定精度が変化する場合に、推定精度の低下を抑制することができる。
また例えば、前記複数の投票領域は、第１の投票領域と、前記第１の投票領域よりも前記第１の物体の中心に近い第２の投票領域とを含み、前記補正部は、前記第１の投票領域における前記第２の物体の存在確率が前記第２の投票領域よりも相対的に低下するように前記投票結果を補正してもよい。
この構成によれば、第１の投票領域における第２の物体の存在確率が第２の投票領域よりも相対的に低下するように投票結果を補正することができる。ここで、第２の投票領域は、第１の投票領域よりも第１の物体の中心に近い。したがって、第２の物体を透過する光の光路モデル上の投票領域のうち、第２の物体が存在しない投票領域に第２の物体が存在すると推定されることを減少させることができる。
なお、これらの包括的または具体的な態様は、システム、方法、集積回路、コンピュータプログラムまたはコンピュータ読み取り可能なＣＤ−ＲＯＭなどの記録媒体で実現されてもよく、システム、方法、集積回路、コンピュータプログラムおよび記録媒体の任意な組み合わせで実現されてもよい。
以下、実施の形態について、図面を用いて詳細に説明する。
なお、以下で説明する実施の形態は、いずれも包括的または具体的な例を示す。以下の実施の形態で示される数値、形状、材料、構成要素、構成要素の配置および接続形態、ステップ、ステップの順序などは、一例であり、請求の範囲を限定する主旨ではない。また、以下の実施の形態における構成要素のうち、最上位概念を示す独立請求項に記載されていない構成要素は、任意の構成要素として説明される。
（実施の形態１）
実施の形態１に係る物体内部推定装置は、物体の内部を推定する。特に、本実施の形態においては、物体内部推定装置は、散乱体内に含まれる光吸収物体の位置および形状を推定する。
ここで、散乱体とは、その内部で散乱が生じる物体のことである。また、光吸収物体とは、光を吸収するあらゆる物体を含み、その吸収率等は問わない。散乱体は、固体、液体および気体のいずれであってもよい。例えば、散乱体が液体の水である場合、物体内部推定装置は、水中の吸収物体を推定することができる。また例えば、散乱体が気体の空気と液体の水滴との混合物である場合、物体内部推定装置は、霧の内部を推定することができる。なお、散乱体は、第１の物体の一例である。また、光吸収物体は、第２の物体の一例である。
＜機能構成＞
図１は、実施の形態１に係る物体内部推定装置１０の機能構成を示すブロック図である。図１に示すように、物体内部推定装置１０は、第１の取得部１１と、第２の取得部１２と、投票値算出部１３と、投票部１４と、推定部１５と、統合部１６とを含む。なお、物体内部推定装置１０は、統合部１６を含まなくてもよい。例えば、物体内部推定装置１０の外部の装置が統合部１６を備えてもよい。
（第１の取得部１１）
第１の取得部１１は、光吸収物体が散乱体内に含まれていないと仮定した場合における、複数の光路モデルと第１の光強度分布とを取得する。複数の光路モデルとは、光が散乱体に入射してから出射するまでの間にとる複数の経路を表す。第１の光強度分布とは、複数の光路モデルに沿って散乱体内を伝搬して散乱体から出射した仮想的な光の所定面上における強度を示す。所定面とは、散乱体外部の平面あるいは曲面であり、散乱体から出射した光が到達する面である。つまり、第１の取得部１１は、散乱体内における仮想的な複数の光路と、散乱体外における仮想的な光強度分布とを取得する。
本実施の形態では、第１の取得部１１は、散乱シミュレーションモデル（以下、単に散乱モデルと呼ぶ）を用いて、散乱体内の光が通過する分割領域（例えばボクセル）の集合を光路モデルとして求める。また、第１の取得部１１は、各光路モデルに沿って散乱体内を伝搬した光の透過率を求める。光の透過率とは、入射光の強度に対する出射光の強度の比率である。つまり、複数の光の透過率は、第１の光強度分布に相当する。
（第２の取得部１２）
第２の取得部１２は、散乱体に入射してから出射した光の強度分布であって所定面上において観測された光の強度分布である第２の光強度分布を取得する。上述したように、所定面は、散乱体に入射してから出射した光が到達する面である。例えば、所定面は、散乱体への光の入射方向に垂直な面である。
本実施の形態では、第２の取得部１２は、散乱体を透過した光（散乱光および直進光）を撮影することにより、第２の強度分布を示す観測画像を取得する。具体的には、第２の取得部１２は、可視光または近赤外光を散乱体に照射する光源部１２ａと、散乱体に入射してから出射した光を撮像することにより、第２の強度分布を示す観測画像を撮影する撮像部１２ｂとを含む。
例えば、光源部１２ａは、可視光線または近赤外光線を出力するレーザである。また例えば、撮像部１２ｂは、レンズユニットおよび撮像素子等を備えるカメラである。
（投票値算出部１３）
投票値算出部１３は、第２の光強度分布が観測された所定面上の位置ごとに、第１の光強度分布および第２の光強度分布間における光の強度の比較結果に基づいて投票値を算出する。この投票値は、散乱体内に光吸収物体が存在する確率を示す。逆に言えば、投票値は、散乱体内に光吸収物体が存在しない確率を示す。つまり、投票値は、所定面上の各位置に到達する光の光路上に光吸収物体が存在する確率および存在しない確率を示す。
なお、以下において、値が大きいほど光吸収物体が存在する可能性が高い確率を存在確率と呼ぶ。逆に、値が大きいほど光吸収物体が存在しない可能性が高い確率を非存在確率と呼ぶ。
投票値算出部１３は、例えば、第２の光強度分布における光の強度値（以下、第２の強度値と呼ぶ）を第１の光強度分布における光の強度値（以下、第１の強度値と呼ぶ）で除算した結果に比例する値を、非存在確率を示す投票値として算出してもよい。また例えば、投票値算出部１３は、第１の強度値と第２の強度値との差分値に比例する値を、存在確率を示す投票値として算出してもよい。また例えば、投票値算出部１３は、第１の強度値と第２の強度値との差分値を第１の強度値で除算した結果に比例する値を、存在確率を示す投票値として算出してもよい。
本実施の形態では、投票値算出部１３は、撮像部１２ｂによって撮影された観測画像の各画素値（観測値）と、散乱モデルにより得られる、各画素に対応する光透過率との比率を、各画素に到達する光の光路上に光吸収物体が存在しない確率（非存在確率）として求める。
散乱モデルにより得られる光透過率は、散乱体内に光吸収物体が存在しないと仮定して求められている。そのため、ある画素において光透過率と観測値とが等しい場合、その画素に到達する光の光路上に光吸収物体が存在しないと推定される。逆に、ある画素において観測値が光透過率よりも小さいほど、その画素に到達する光の光路上に光吸収物体が存在する確率が高いと推定される。
（投票部１４）
投票部１４は、散乱体の内部を複数の投票領域に仮想的に分割する。そして、投票部１４は、第２の光強度分布が観測された所定面上の位置ごとに、当該位置に到達する光の光路モデル上の投票領域に対して当該位置の投票値を投票する。具体的には、投票部１４は、例えば、各位置に到達する光の光路モデル上の各投票領域が有する値に各位置の投票値を加算する。
本実施の形態では、投票部１４は、投票値算出部１３で求めた、各画素における光吸収物体の非存在確率を示す投票値を、各画素に対応する光路モデル上の投票領域に投票する。
（推定部１５）
推定部１５は、投票結果の値の大きさに基づいて、散乱体内の光吸収物体の位置を推定する。つまり、推定部１５は、投票結果に基づいて、複数の投票領域の各々に光吸収物体が存在するか否かを推定する。具体的には、推定部１５は、例えば、各投票領域が有する値（投票結果）を閾値と比較する。そして、推定部１５は、比較結果に基づいて、当該投票領域に光吸収物体が存在するか否かを推定する。
本実施の形態では、複数の投票領域は、散乱体の断面上の複数の領域である。そして、推定部１５は、複数の断面の各々について、複数の投票領域の各々に光吸収物体が存在するか否かを推定することにより、複数の２次元断層画像を生成する。
（統合部１６）
統合部１６は、複数の２次元断層画像を統合することにより、散乱体内の３次元ボリュームデータを生成する。
＜動作＞
次に、以上のように構成された物体内部推定装置１０の処理動作を説明する。
図２は、実施の形態１に係る物体内部推定装置１０の処理動作を示すフローチャートである。
（ステップＳ１０１）
まず、第１の取得部１１によるモデル生成処理（Ｓ１０１）について説明する。モデル生成処理では、第１の取得部１１は、複数の光路モデルと光透過率とを算出する。散乱をモデル化する方法として、本実施の形態では、モンテカルロ・レイトレーシングを用いる。つまり、第１の取得部１１は、複数の光路モデルを、モンテカルロ・レイトレーシング・シミュレーションによって導出する。
モンテカルロ・レイトレーシングは、ＣＧ分野において、散乱や屈折などを伴う複雑な光線の経路をモデル化する手法として利用されている。具体的なモンテカルロ・レイトレーシングによる光路のモデル化の詳細については、例えば非特許文献２に記載されている。
図３は、モンテカルロ・レイトレーシング・シミュレーションによる光路モデルの導出処理の詳細を示すフローチャートである。第１の取得部１１は、この処理フローに従って光路モデルをシミュレートする。そして、第１の取得部１１は、後述する散乱モデルに従って光の透過率を算出する。
光路は光子の移動によって形成される。ここでは、３次元空間（ｘ、ｙ、ｚ）での光子の移動についてシミュレーションする方法について述べる。ある２次元断面（ｘ、ｙ）上を、光子がｚ＝０から散乱体の周縁まで散乱しながら移動するときの光路を求める。
（ステップＳ３０１）
光子生成処理（Ｓ３０１）では、光子の初期位置と初期移動方向とを設定する。例えば、光子の初期位置を（ｘ、ｙ、ｚ）＝（０、０、０）と設定し、移動方向に関するベクトルは（μｘ、μｙ、μｚ）＝（０、０、１）と設定する。
（ステップＳ３０２）
光子移動距離設定処理（Ｓ３０２）では、光子の移動ステップｓを求める。ｓは、Ｌａｍｂｅｒｔ−Ｂｅｅｒ則に従って、式（１）のように算出される。
ここで、μtは、散乱係数μsと吸収係数μaの和である消滅係数である。光子の現在位置における消滅係数と、乱数ξ（０＜ξ≦１）とを用いて、光子の移動ステップｓが定義される。移動ステップｓが決まれば、光子の移動先の位置が式（２）のように設定される。
（ステップＳ３０３）
吸収重み算出処理（Ｓ３０３）では、吸収係数と消滅係数との比から重みの減少量ΔＷを求める。この重みの初期値Ｗは対象物に応じて適切に決められる。重みが０になれば、光子が物体内部で吸収され、物体外に出射しないことを意味する。そのため、Ｗが大きければ、光子が物体内部で吸収されにくく、Ｗが小さければ、光子が物体内部で吸収されやすい。Ｗは、対象物体の媒質に応じて設定される。
（ステップＳ３０４）
散乱方向設定処理（Ｓ３０４）では、位相関数のパラメータｇと乱数ξ（０＜ξ≦１）とを用いて、式（４）のように、光子の移動角度θが決定される。位相関数とは、散乱の偏りを近似的に表現した関数である。この位相関数のパラメータｇは散乱の等方性を意味している。すなわち、ｇが正の値を有する場合は前方散乱となり、ｇが０の値を有する場合は等方散乱となり、ｇが負の値を有する場合は後方散乱となる。
ここで、ψ＝２πξとすると、光子の移動方向θは、式（５）、（６）、（７）で示される。
ただし、μz＝１の場合は、式（８）、（９）、（１０）が用いられ、μz＝−１の場合は、式（１１）、（１２）、（１３）が用いられる。
（ステップＳ３０５）
重み打ち切り処理（Ｓ３０５）は、ロシアンルーレット方式と呼ばれる光子の打ち切りを行う。具体的には、乱数ξ（０≦ξ≦１）と定数ｍとを用いて、式（１４）に示すように、乱数ξが１／ｍ以下であればＷをｍＷに更新し、乱数ξが１／ｍより大きければ、光子を消滅させる。
（ステップＳ３０６）
以上までの処理の後、生き残った光子（消滅しなかった光子）に対して、引き続き、光子移動距離設定処理が実行される。具体的には、更新されたμ’x、μ’y、μ’zと現在位置の消滅係数とを用いて、光子の移動ステップｓが求められる。このようにモンテカルロ・レイトレーシング・シミュレーションでは、移動距離と移動方向とが乱数によって確率的に摂動された光路モデルが求められる。
次に、第１の取得部１１は、ここまでの処理で求めた、散乱体内の複数の光路モデルに対する光の透過率を求める。ここでは、入射点から散乱点を経由して出射点につながる光路モデルが、モンテカルロ・レイトレーシング・シミュレーションで得られている。
図４は、散乱体内の３点ｐ１、ｐ２、ｐ３が与えられたときの、光強度の減衰量を説明するための図である。点ｐ１からの光が距離ｄ1だけ離れた散乱点ｐ２にて、角度θの方向に散乱し、距離ｄ2だけ離れた点ｐ３まで到達する光路において、光伝搬の割合Ｅp1,p2,p3は、式（１５）で表される。
ここで、ｐ（θ、ｇ）は、位相関数と呼ばれ、散乱の偏りを近似的に表現する関数であり、式（１６）で定義される。
モンテカルロ・レイトレーシング・シミュレーションの説明でも述べたように、位相関数のパラメータｇは散乱の等方性を意味している。すなわち、ｇが正の値を有する場合は前方散乱となり、ｇが０の値を有する場合は等方散乱となり、ｇが負の値を有する場合は後方散乱となる。
モンテカルロ・レイトレーシング・シミュレーションで得られる、散乱点ｒ1、ｒ2、…、ｒm、入射点ｒ0、および出射点ｒm+1をつなぐ光路Ｒｊにおける光の透過率について説明する。ある散乱点ｒkにおいて、角度θk-1,k,k+1の方向に光が散乱するとき、光の透過率は、式（１５）の積で表現されることから、光路Ｒｊにおける光透過率ＥＲｊは、式（１７）で定義される。
なお、距離ｄk,k+1と角度θk-1,k,k+1は式（１８）で定義される。
以上までの処理によって、第１の取得部１１は、光路Ｒｊと、その光路Ｒｊにおける光透過率ＥＲｊを求めることができる。ここで、ｊはモンテカルロ・レイトレーシング・シミュレーションにおけるサンプリング回数を示す。サンプリング回数が多いほど、高精度なモデル化が可能である。
なお、モンテカルロ・レイトレーシング・シミュレーションにて用いるパラメータ（散乱係数、吸収係数、および位相関数のパラメータｇ）は、内部推定の対象となる散乱体の実際の値に近いことが望ましい。一般に、これらのパラメータを実物体から推定することは困難である。
そこで、本実施の形態では、事前に光吸収物体を含まない散乱体から出射した光の強度分布を観測しておく。そして、第１の取得部１１は、モンテカルロ・レイトレーシング・シミュレーションの結果として得られる複数の光路モデルから導出される光透過率が示す光の強度分布が、実際の光の強度分布に近くなるようにパラメータ探索を行う。
例えば、非特許文献１では、散乱が一度だけ生じる単一散乱現象のときの、各種係数を推定した例が開示されている。また、非特許文献４では、皮膚などの生体に関する各種係数の例も開示されている文献がある。これらの文献に記載された各種係数を参考にして、第１の取得部１１は、実際に内部を推定する散乱体に近い媒質を有する物体のパラメータ（散乱係数、吸収係数、および位相関数のパラメータｇ）の初期値を設定する。そして、第１の取得部１１は、パラメータを微小に変化させながら、シミュレーション結果が示す光透過率によって得られる光強度分布と観測結果が示す光強度分布との誤差が最小となるパラメータを求める。つまり、第１の取得部１１は、シミュレーション結果と観測結果との誤差を粗密探索によって削減することにより、モンテカルロ・レイトレーシング・シミュレーションで用いるパラメータを導出する。これにより、第１の取得部１１は、パラメータの精度を向上させることができる。
（ステップＳ１０２）
次に、第２の取得部１２による撮影処理（Ｓ１０２）について説明する。撮影処理では、第２の取得部１２は、実際の散乱体の撮影を行う。まず、単方向からの散乱体の撮影方法について述べた後、信頼性が高い２次元断層画像を得るために多方向から散乱体を撮影する方法について述べる。
図５は、光源部１２ａ、散乱体、および撮像部１２ｂの配置の一例を示す模式図である。ここでは、光源にレーザを利用する場合を例にして説明する。レーザの種類は、可視光レーザでも近赤外光レーザでも構わない。
例えば、薄めた牛乳など可視光でもある程度透過するような媒質を有する散乱体の内部を推定する場合は、可視光レーザを利用できる。また例えば、食品や生体など可視光は透過しないが、近赤外光は透過するような媒質を有する散乱体の場合は、近赤外光レーザやＬＥＤ光源を用いることが望ましい。
可視光と近赤外光とでは、波長が異なるので、散乱係数、吸収係数、および位相関数のパラメータも異なる。そのため、モンテカルロ・レイトレーシング・シミュレーションに利用するパラメータを、光源と散乱体の媒質とに合わせて調整する必要がある。なお、可視光を用いる場合と近赤外光を用いる場合とで光吸収物体の推定処理に区別は無い。以下、散乱体に入射する光を「入射光」、散乱体から出射する光を「出射光」と呼ぶ。
図５では、光源部１２ａは、光を散乱体に照射する。また、撮像部１２ｂは、出射光を観測する。例えば、散乱体の水平断面上の光吸収物体の位置を推定する場合、図６および図７に示すように、観測画像から、推定対象断面に対応する水平方向（ｘ方向）の画素列が取得される。図７は、取得された画素列の各画素の位置（ｘ座標）を横軸に、画素値を縦軸に表したグラフである。
以上の処理で単方向から入射した光を観測できる。しかし、単方向のみの情報では、奥行き方向の情報が不十分であるため、光吸収物体の位置推定精度が低い。そこで、散乱体を回転台の中心に載せて散乱体を回転させながら撮影してもよい。つまり、撮像部１２ｂは、複数方向から光が散乱体に入射したときの出射光を撮像してもよい。
すなわち、散乱体に互いに異なる方向から入射する複数の光の各々について、ステップＳ１０１〜ステップＳ１０４が実行されてもよい。具体的には、第１の取得部１１は、複数の光の各々について、複数の光路モデルと第１の光強度分布とを取得してもよい。第２の取得部１２は、複数の光の各々について第２の光強度分布を取得してもよい。投票値算出部１３は、複数の光の各々について投票値を算出してもよい。投票部１４は、複数の光の各々について投票値を投票してもよい。
これにより、奥行き方向への推定精度を向上させることができる。なお、本実施の形態では、回転台を用いて散乱体を回転させる方法について述べたが、散乱体の周りを光源部１２ａと撮像部１２ｂとが回転しても構わない。
（ステップＳ１０３）
次に、投票値算出部１３による投票値算出処理（Ｓ１０３）について述べる。投票値算出処理（Ｓ１０３）と投票処理（Ｓ１０４）は、散乱体の回転角度に依存しないので、単方向から光が入射する場合を例に説明する。
図８は、モンテカルロ・レイトレーシング・シミュレーションによって導出された、画素Ａ、画素Ｂ、および画素Ｃにそれぞれ出射する３つの光路モデルの例を示している。本説明では、左の光路モデルを「画素Ａに対応する光路モデル」、中央の光路モデルを「画素Ｂに対応する光路モデル」、右の光路モデルを「画素Ｃに対応する光路モデル」と呼ぶ。
なお、１つの方向から散乱体に入射した光には、図８に示すような複数の光路モデルが導出される。つまり、１つの入射光は、複数の経路を伝搬して複数の位置から出射する。
図９は、散乱モデルによる光透過率の分布（第１の光強度分布）と観測値の分布（第２の光強度分布）の一例である。図９において、横軸は推定対象断面上の位置を示す。画素Ａ、画素Ｂ、および画素Ｃにそれぞれ対応する位置をＡ、Ｂ、Ｃと表示している。また、図９において、縦軸は画素値（光透過率および観測値）を示す。
光透過率と観測値とはスケールが異なるため、光透過率の最大値と観測値の最大値とが一致するように、光透過率の分布と観測値の分布とを正規化する。なお、図９に示す分布は、図６に示す円柱形の散乱体内に２つの円柱形の光吸収物体が配置されたときの例である。
式（１７）で得られる個々の光路モデルの光透過率を合算した後、正規化した正規化光透過率が求められる。ここで、画素ｉの正規化光透過率の値をＮＥＲｉとし、撮像部１２ｂで観測された画素ｉの画素値をＴｉとする。
モンテカルロ・レイトレーシング・シミュレーションで得られるＮ個の光路モデルＲｊ（ｊは１〜Ｎ）上に光吸収物体が存在しない確率Ｖ（Ｒｊ）は、式（１９）で定義される。
ここで、画素ｉの位置に出射する光線の光路における確率をＶ（Ｒｉｊ）と定義する。画素ｉの位置に出射する光線の正規化光透過率ＮＥＲｉおよび観測画素値Ｔｉと、光吸収物体が存在しない確率Ｖ（Ｒｉｊ）との関係は、式（２０）で表現される。
光吸収物体が存在しない確率Ｖ（Ｒｉｊ）は、画素ｉの画素値Ｔｉを用いて式（２１）のように表される。
光吸収物体の非存在確率は、観測した画素値と散乱モデルによる正規化光透過率との比で求められる。以上により、各画素に到達する光の経路上に光吸収物体が存在しない確率が求められる。投票値算出処理では、このように求められる非存在確率が投票値として画素ごとに出力される。
（ステップＳ１０４）
投票部１４による投票処理（Ｓ１０４）について説明する。投票処理では、投票部１４は、求められた非存在確率を、モンテカルロ・レイトレーシング・シミュレーションで求めた光路モデル上の投票領域に投票することで、散乱体内での光吸収物体の非存在確率の分布を求める。
投票処理について、図１０および図１１を用いて説明する。
図１０のように、投票値算出処理にて、画素Ａ、画素Ｂ、画素Ｃにおける光吸収物体の非存在確率（投票値）が０．２、０．９、０．１と算出されている。ここで、図１１のように、画素Ａ、画素Ｂ、画素Ｃに到達する光の光路モデル上の投票領域の値に、算出された投票値を加算していく。ここでは、投票領域は、散乱体を仮想的に分割して得られるボクセルである。図１１の（ａ）のように、画素Ａに到達する光の光路モデル上の各ボクセルには、投票値０．２が投票される。
画素Ａおよび画素Ｃに到達する光の光路モデル上の各投票領域に投票される投票値は、それぞれ０．２および０．１と比較的小さな値である（図１１の（ａ）および（ｃ））。逆に、画素Ｂに到達する光の光路モデル上の各投票領域に投票される投票値は、比較的大きな値である（図１１の（ｂ））。
観測値は、複数の出射光の和に対応するため、光吸収物体が存在していても完全に０になることは稀である。しかしながら、小さな投票値を複数回投票することにより、光吸収物体の非存在確率の差が投票領域ごとに明確になる。
（ステップＳ１０５）
推定部１５による推定処理（Ｓ１０５）について説明する。推定処理では、散乱体内の光吸収物体の推定を行う。具体的には、推定部１５は、投票領域ごとに、投票結果が閾値より小さければ、当該投票領域に光吸収物体が存在すると判定する。
以上の処理によって、ある断面における光吸収物体の推定結果が得られる。そして、推定部１５は、位置（例えば高さ）が異なる複数の断面の各々について、光吸収物体の推定を行うことなり、複数の２次元断層画像を生成する。統合部１６は、生成された複数の２次元断層画像を統合することで３次元ボリュームデータ（３次元画像）を生成する。
図６の例においては、画素列のｙ座標を変化させて、各断面における推定処理を行っていく。統合部１６は、単に、２次元断層画像を並べるだけでなく、用途に応じた信号処理（例えば、断層画像間の連続性を高めるための処理、あるいはノイズ除去のための平滑化処理など）を２次元断層画像に適用することが可能である。これらの信号処理としては、既存のボリュームレンダリングの高画質化処理が用いられればよい。
＜効果＞
以上のように、本実施の形態に係る物体内部推定装置１０によれば、光が散乱体に入射してから出射するまでの間にとる複数の経路を表す複数の光路モデルを用いて、散乱体内の各領域に光吸収物体が存在する否かを推定することができる。したがって、散乱体内で光が散乱する場合でも、散乱体内の光吸収物体の位置および形状を推定することができる。つまり、近赤外光などＸ線以外の人体に安全な光を用いて物体の内部を推定する場合に推定精度を向上させることができる。
さらに、本実施の形態に係る物体内部推定装置１０によれば、第１の光強度分布と第２の光強度分布間における光の強度の比較結果に基づいて、散乱体内に光吸収物体が存在する確率を示す投票値を算出することができる。また、光路モデル上の投票領域に対する投票値の投票結果に基づいて、複数の投票領域の各々に光吸収物体が存在するか否かを推定することができる。つまり、散乱体の内部を推定するために、複数の波長の光を用いる必要がなく、線形演算の必要もない。したがって、吸収率分布の線形演算のための近似も不要となるので、物体の内部の推定精度を向上させることができる。
また、本実施の形態に係る物体内部推定装置１０によれば、複数の２次元断層画像を統合することにより、３次元ボリュームデータを生成することができる。したがって、より正確に散乱体の内部を推定することができる。
（実施の形態２）
次に、実施の形態２に係る物体内部推定装置について説明する。実施の形態１と同様に、実施の形態２に係る物体内部推定装置は、物体の内部を推定する。ただし、本実施の形態では、投票結果を補正する点が実施の形態１と異なる。
以下、本実施の形態に係る物体内部推定装置について図面を参照しながら具体的に説明する。なお、本実施の形態において、実施の形態１と同様の機能あるいは処理については、説明が不必要に冗長になるのを避けるために説明を省略する場合がある。
＜機能構成＞
図１２は、実施の形態２に係る物体内部推定装置２０の機能構成を示すブロック図である。
補正部２１は、散乱体内の投票領域を通過する光路の数あるいは入射光の入射位置によって生じる投票値の偏りを補正する。具体的には、補正部２１は、例えば、複数の投票領域に対する投票回数のばらつきに基づいて投票結果を補正してもよい。また例えば、補正部２１は、複数の投票領域の各々の位置に基づいて投票結果を補正してもよい。
推定部１５は、補正された投票結果に基づいて、複数の投票領域の各々に光吸収物体が存在するか否かを推定する。
＜動作＞
次に、以上のように構成された物体内部推定装置２０の処理動作について説明する。図１３は、実施の形態２に係る物体内部推定装置２０の処理動作を示すフローチャートである。図１３では、図２のフローチャートに対して、補正処理（Ｓ２０１）が追加されている。
投票処理では、投票回数のばらつき、あるいは光の入射位置設定によって投票結果に偏りが生じる。そのため、補正処理（Ｓ２０１）では、以下のような処理を行うことにより、投票結果の偏りを削減する。
投票回数のばらつきに関しては、各投票領域の投票回数を記録しておき、投票回数を用いて投票結果の正規化を行う。例えば、補正部２１は、投票領域ごとに、当該領域に投票された投票値の和を投票回数で除算することにより、正規化された投票結果を得る。ここで、投票回数が０回となる投票領域は、不定領域と定義され、投票結果の計算対象から除外される。
例えば、図１４は、入射光が散乱体内で散乱したときに得られる複数の光路モデルを示す模式図である。図１４には、光の入射角が０度、４５度、および９０度であるときの複数の光路モデルの例が示されている。断面が円形である散乱体内に光吸収物体が２つ存在する。レーザのような直進性の高い光が散乱体に入射する場合、図１４に実線矢印で示すように、複数の光路モデルは扇状に形成される。このとき、入射位置の左右の領域は、光路が存在しないので不定領域となる。
第２の取得部１２は、図１４に示すように、散乱体を回転して撮影する。そして、投票部１４は、複数の入射角度に対して独立に投票処理を行う。そして、補正部２１は、不定領域となる投票領域の投票結果を除いて、複数の入射角度に対して得られた投票結果の平均値を算出する。つまり、いずれかの入射角度で不定領域とならなければ、投票領域は、投票結果を得られる。このように、補正部２１は、不定領域が投票結果に影響を及ぼさないように投票結果を算出する。
次に、光の入射位置設定による投票結果の偏りの補正について述べる。そもそも、散乱体に光を入射し撮影する場合、光の入射方向と入射位置における散乱体表面の法線方向とが対向していない場合、入射光が散乱体表面で屈折する。入射光が散乱体表面で屈折すると、物体の形状やフレネル係数などを考慮した散乱モデルが必要となるため、散乱モデル生成の難易度が高くなる。
そこで、光源部１２ａは、例えば、断面が円形の散乱体において、光の入射方向と入射位置における散乱体表面の法線方向とが対向する位置から光を散乱体に照射する。これにより、散乱体の表面における屈折の問題を回避する。しかしながら、このように光の入射位置に制限が加えられることによって、投票回数が十分な投票領域であっても、正しい投票結果が得られない場合がある。
例えば、図１４に示す「光線Ａ」の光路モデル上には、光吸収物体が存在する。そのため、「光線Ａ」の光路モデル上の各投票領域には、光吸収物体の存在確率が高いことを示す投票値が投票される。しかしながら、「光線Ａ」の光路モデルの下流側（出射位置付近）には、光吸収物体は存在しない。ここで、光の入射位置が限定される場合には、光吸収物体を横切らずに光線Ａの出射位置付近の投票領域を通過する光路モデルが存在しない場合がある。このような場合には、投票部１４は、実際には光吸収物体が存在しない投票領域に、光吸収物体が存在しないことを示す投票値を投票できない。
そこで、補正部２１は、複数の投票領域の各々の位置に基づいて投票結果を補正する。具体的には、補正部２１は、以下に示す散乱体内の位置に依存する信頼度に基づいて投票結果を補正する。例えば、断面が円形の散乱体に、入射位置における表面の法線と対向する位置から光が入射する場合、散乱体断面の中心部の投票値の信頼度が高く、周縁部（表面）に近いほど投票値の信頼度が低い。
そこで、補正部２１は、第１の投票領域における光吸収物体の存在確率が第２の投票領域よりも相対的に低下するように投票結果を補正する。ここでは、第２の投票領域は、第１の投票領域よりも散乱体の中心に近い。逆に言えば、第１の投票領域は、第２の投票領域よりも散乱体の表面に近い。具体的には、例えば投票値が非存在確率を示す場合、補正部２１は、第１の投票領域の投票結果の値を第１の値だけ増加させ、第２の投票領域の投票結果の値を第１の値より小さい第２の値だけ増加させる。
本実施の形態では、補正部２１は、図１５に示すような補正係数を用いて投票結果を補正する。この補正係数は、散乱体の中心から所定の距離ｒ以下の位置においては「１」の値を有し、距離ｒより離れた位置においてはｅｘｐ（α｜ｘ｜）の値を有する。補正部２１は、この補正係数を投票結果の値に掛けることにより、非存在確率を示す投票値による投票結果を補正する。これにより、補正部２１は、光路モデルが扇状に広がることにより、散乱体の周縁部における投票値の信頼度が低下する場合にも、適切な投票結果を得ることができる。なお、断面が円形以外の散乱体に対しても、補正部２１は、光路モデルの分布に基づく信頼度に応じた補正係数を用いることで、投票結果を適切に補正することが可能となる。
このように補正された投票結果に基づいて、推定処理（Ｓ１０５）において推定部１５は、各投票領域における光吸収物体の推定を行う。具体的には、推定部１５は、投票領域において投票値が閾値より小さければ、当該投票領域に光吸収物体が存在すると判定する。
＜効果＞
以上のように、本実施の形態に係る物体内部推定装置２０によれば、光の入射位置が限定される場合や、散乱度合いによって光が届かない場所が散乱体の表面近くに存在する場合でも、散乱体の内部を高精度に推定することが可能になる。
つまり、本実施の形態に係る物体内部推定装置２０によれば、複数の投票領域に対する投票回数のばらつきに基づいて投票結果を補正することができる。したがって、例えば、散乱体内の光が通過しにくい領域（投票回数が少ない領域）に光吸収物体が存在すると誤って推定されることを減少させることができる。
また、本実施の形態に係る物体内部推定装置２０によれば、各投票領域の位置に基づいて投票結果を補正することができる。したがって、投票領域の位置に依存して光吸収物体が存在するか否かの推定精度が変化する場合に、推定精度の低下を抑制することができる。
特に、本実施の形態に係る物体内部推定装置２０によれば、第１の投票領域における光吸収物体の存在確率が第２の投票領域よりも相対的に低下するように投票結果を補正することができる。ここで、第２の投票領域は、第１の投票領域よりも散乱体の中心に近い。したがって、光吸収物体を透過する光の光路モデル上の投票領域のうち、光吸収物体が存在しない投票領域に光吸収物体が存在すると推定されることを減少させることができる。
（他の実施の形態）
以上、１つまたは複数の態様に係る物体内部推定装置について、実施の形態に基づいて説明したが、本発明は、この実施の形態に限定されるものではない。本発明の趣旨を逸脱しない限り、当業者が思いつく各種変形を本実施の形態に施したものや、異なる実施の形態における構成要素を組み合わせて構築される形態も、１つまたは複数の態様の範囲内に含まれてもよい。
例えば、上記各実施の形態において、第１の取得部１１は、推定のたびに複数の光路モデルと光の透過率とを求めていたが、必ずしも推定のたびに求める必要はない。推定対象となる散乱体が事前にわかっていれば、その散乱体のために前に導出された複数の光路モデルと光の透過率とを記憶している記憶部（メモリ）から、複数の光路モデルと光の透過率とを読み込んでもよい。これにより、第１の取得部１１の処理負荷を軽減することができる。
また例えば、上記各実施の形態では、第１の取得部１１は、複数の光路モデルおよび光の透過率の導出処理を行っていたが、必ずしも導出処理を行わなくてもよい。例えば、第１の取得部１１は、推定対象となる散乱体のために既に求められた複数の光路モデルおよび光の透過率を外部の装置から取得してもよい。
また例えば、上記各実施の形態において、複数の光路モデルは、モンテカルロ・レイトレーシングによって求められていたが、他の散乱モデル化手法によって求められてもよい。
また例えば、上記各実施の形態では、第１の光強度分布が光の透過率である場合について説明したが、第１の光強度分布は光の透過率でなくてもよい。つまり、第１の光強度分布は、複数の光路モデルに沿って散乱体内を伝搬して散乱体から出射した仮想的な光の所定面上の強度を示すデータであれば、どのようなデータであってもよい。
また例えば、上記各実施の形態において、第２の取得部１２は、光源部１２ａおよび撮像部１２ｂを備えていたが、光源部１２ａおよび撮像部１２ｂを備えなくてもよい。例えば、物体内部推定装置の外部に設けられた光源が散乱体に光を照射してもよい。また例えば、第２の取得部１２は、物体内部推定装置の外部に設けられたカメラで撮影された観測画像をネットワークあるいは記録媒体等を介して取得してもよい。
また例えば、上記各実施の形態において、観測画像は、散乱体内を透過する光を撮像したものであったが、散乱体内で反射する光を撮像したものであっても構わない。
また例えば、上記各実施の形態では、第２の光強度分布が画像である場合について説明したが、第２の光強度分布は画像でなくてもよい。つまり、第２の光強度分布は、所定面上で観測された、散乱体から出射した光の強度分布を示すデータであれば、どのようなデータであってもよい。
また例えば、上記各実施の形態では、各投票領域に投票された投票値の和を投票結果として利用していたが、例えば、投票値の統計的な代表値（平均値、最大値あるいは最小値など）を投票結果として利用してもよい。適切な代表値を用いることで、投票回数のばらつきによる推定精度の低下を抑制することができる。
また例えば、上記各実施の形態において、散乱体の断面ごとに光吸収物体の推定処理が行われていたが、必ずしも断面毎に推定処理を行う必要はない。例えば、推定処理は、３次元で行われてもよい。この場合、複数の２次元断層画像を統合しなくても、散乱体の内部の３次元画像を得ることができる。つまり、物体内部推定装置は、統合部を備えなくてもよい。
また、上記各実施の形態における物体内部推定装置が備える構成要素の一部または全部は、１個のシステムＬＳＩ（Ｌａｒｇｅ Ｓｃａｌｅ Ｉｎｔｅｇｒａｔｉｏｎ：大規模集積回路）から構成されているとしてもよい。例えば、物体内部推定装置は、第１の取得部１１と、第２の取得部１２と、投票値算出部１３と、投票部１４と、推定部１５とを有するシステムＬＳＩから構成されてもよい。
システムＬＳＩは、複数の構成部を１個のチップ上に集積して製造された超多機能ＬＳＩであり、具体的には、マイクロプロセッサ、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）などを含んで構成されるコンピュータシステムである。前記ＲＯＭには、コンピュータプログラムが記憶されている。前記マイクロプロセッサが、前記コンピュータプログラムに従って動作することにより、システムＬＳＩは、その機能を達成する。
なお、ここでは、システムＬＳＩとしたが、集積度の違いにより、ＩＣ、ＬＳＩ、スーパーＬＳＩ、ウルトラＬＳＩと呼称されることもある。また、集積回路化の手法はＬＳＩに限るものではなく、専用回路または汎用プロセッサで実現してもよい。ＬＳＩ製造後に、プログラムすることが可能なＦＰＧＡ（Ｆｉｅｌｄ Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）、あるいはＬＳＩ内部の回路セルの接続や設定を再構成可能なリコンフィギュラブル・プロセッサを利用してもよい。
さらには、半導体技術の進歩または派生する別技術によりＬＳＩに置き換わる集積回路化の技術が登場すれば、当然、その技術を用いて機能ブロックの集積化を行ってもよい。バイオ技術の適用等が可能性としてありえる。
また、本発明は、このような物体内部推定装置として実現することができるだけでなく、物体内部推定装置に含まれる特徴的な構成要素をステップとする物体内部推定方法として実現することもできる。また、物体内部推定方法に含まれる特徴的な各ステップをコンピュータに実行させるコンピュータプログラムとして実現することもできる。そして、そのようなコンピュータプログラムを、コンピュータ読み取り可能な非一時的な記録媒体あるいはインターネット等の通信ネットワークを介して流通させることができるのは、言うまでもない。
また、上記各実施の形態における一部の処理をサーバ（クラウド）が行い、残りの処理を物体内部推定装置が行うという実施の形態でもよい。例えば、物体内部推定装置は、第１の取得部１１と、第２の取得部１２と、推定部１５のみを備え、投票値算出部１３と、投票部１４と、統合部１６は、サーバが備える形態等も考えられる。これに限らず、物体内部推定装置が処理量の小さい構成要素のみを備えるようにすれば、処理の効率化が図れる。
なお、上記各実施の形態において、各構成要素は、専用のハードウェアで構成されるか、各構成要素に適したソフトウェアプログラムを実行することによって実現されてもよい。各構成要素は、ＣＰＵまたはプロセッサなどのプログラム実行部が、ハードディスクまたは半導体メモリなどの記録媒体に記録されたソフトウェアプログラムを読み出して実行することによって実現されてもよい。ここで、上記各実施の形態の物体内部推定装置などを実現するソフトウェアは、次のようなプログラムである。
すなわち、このプログラムは、コンピュータに、第１の物体の内部を推定する物体内部推定方法であって、前記第１の物体とは異なる第２の物体が前記第１の物体内に含まれていないと仮定した場合における、（i）光が前記第１の物体に入射してから出射するまでの間にとる複数の経路を表す複数の光路モデルと、（ii）前記複数の光路モデルに沿って前記第１の物体内を伝搬して前記第１の物体から出射した仮想的な光の所定面上における強度を示す第１の光強度分布とを取得する第１の取得ステップと、前記第１の物体に入射してから出射した光の強度分布であって前記所定面上において観測された光の強度分布である第２の光強度分布を取得する第２の取得ステップと、前記所定面上の位置ごとに、前記第１の光強度分布および前記第２の光強度分布間における光の強度の比較結果に基づいて、前記第１の物体内に前記第２の物体が存在する確率を示す投票値を算出する投票値算出ステップと、前記第１の物体の内部を複数の投票領域に仮想的に分割し、前記所定面上の位置ごとに、当該位置に到達する光の光路モデル上の投票領域に対して当該位置の前記投票値を投票する投票ステップと、投票結果に基づいて、前記複数の投票領域の各々に前記第２の物体が存在するか否かを推定する推定ステップとを含む物体内部推定方法を実行させる。
本発明は、第１の物体内に含まれる第２の物体を正確に推定することができる物体内部推定装置、および、その装置を備える計測機器、監視カメラ、デジタルカメラなどの撮像装置として有用である。
１０、２０ 物体内部推定装置
１１ 第１の取得部
１２ 第２の取得部
１２ａ 光源部
１２ｂ 撮像部
１３ 投票値算出部
１４ 投票部
１５ 推定部
１６ 統合部
２１ 補正部

Claims (11)

  1. 第１の物体の内部を推定する物体内部推定装置であって、
    前記第１の物体とは異なる第２の物体が前記第１の物体内に含まれていないと仮定した場合における、（i）光が前記第１の物体に入射してから出射するまでの間にとる複数の経路を表す複数の光路モデルと、（ii）前記複数の光路モデルに沿って前記第１の物体内を伝搬して前記第１の物体から出射した仮想的な光の所定面上における強度を示す第１の光強度分布とを取得する第１の取得部と、
    前記第１の物体に入射してから出射した光の強度分布であって前記所定面上において観測された光の強度分布である第２の光強度分布を取得する第２の取得部と、
    前記所定面上の位置ごとに、前記第１の光強度分布および前記第２の光強度分布間における光の強度の比較結果に基づいて、前記第１の物体内に前記第２の物体が存在しない非存在確率を示す投票値を算出する投票値算出部と、
    前記第１の物体の内部を複数の投票領域に仮想的に分割し、前記所定面上の位置ごとに、当該位置に到達する光の光路モデル上の投票領域に対して当該位置の前記投票値を投票する投票部と、
    前記複数の投票領域の各々に対して、前記投票部で得られた投票結果が閾値より小さい場合に、前記第２の物体が存在すると推定する推定部とを備える
    物体内部推定装置。
  2. 前記物体内部推定装置は、さらに、
    前記複数の投票領域に対する投票回数のばらつきに基づいて、前記投票結果を補正する補正部を備え、
    前記推定部は、前記複数の投票領域の各々に対して、前記補正部で得られた補正された前記投票結果が閾値より小さい場合に、前記第２の物体が存在すると推定する
    請求項１に記載の物体内部推定装置。
  3. 前記複数の投票領域は、前記第１の物体の断面上の複数の領域であり、
    前記推定部は、複数の断面の各々について、前記複数の投票領域の各々に対して、前記投票結果が閾値より小さい場合に、前記第２の物体が存在すると推定することにより、複数の２次元断層画像を生成し、
    前記物体内部推定装置は、さらに、
    前記複数の２次元断層画像を統合することにより、前記第１の物体内の３次元ボリュームデータを生成する統合部を備える
    請求項１または２に記載の物体内部推定装置。
  4. 前記第２の取得部は、
    可視光または近赤外光を前記第１の物体に照射する光源部と、
    前記第１の物体に入射してから出射した光を撮像することにより、前記第２の光強度分布を示す観測画像を取得する撮像部とを備える
    請求項１〜３のいずれか１項に記載の物体内部推定装置。
  5. 前記第１の取得部は、前記複数の光路モデルを、モンテカルロ・レイトレーシング・シミュレーションによって導出する
    請求項１〜４のいずれか１項に記載の物体内部推定装置。
  6. 前記第１の取得部は、シミュレーション結果と観測結果との誤差を粗密探索によって削減することにより、前記モンテカルロ・レイトレーシング・シミュレーションで用いるパラメータを導出する
    請求項５に記載の物体内部推定装置。
  7. 前記第１の取得部は、前記第１の物体に互いに異なる方向から入射する複数の光の各々について、前記複数の光路モデルと前記第１の光強度分布とを取得し、
    前記第２の取得部は、前記複数の光の各々について前記第２の光強度分布を取得し、
    前記投票値算出部は、前記複数の光の各々について前記投票値を算出し、
    前記投票部は、前記複数の光の各々について前記投票値を投票する
    請求項１〜６のいずれか１項に記載の物体内部推定装置。
  8. 前記物体内部推定装置は、さらに、
    前記複数の投票領域の各々の位置に基づいて前記投票結果を補正する補正部を備え、
    前記推定部は、補正された前記投票結果に基づいて、前記複数の投票領域の各々に前記第２の物体が存在するか否かを推定する
    請求項１〜７のいずれか１項に記載の物体内部推定装置。
  9. 前記複数の投票領域は、第１の投票領域と、前記第１の投票領域よりも前記第１の物体の中心に近い第２の投票領域とを含み、
    前記補正部は、前記第１の投票領域における前記第２の物体の非存在確率が前記第２の投票領域よりも相対的に低下するように前記投票結果を補正する
    請求項８に記載の物体内部推定装置。
  10. 第１の物体の内部を推定する物体内部推定方法であって、
    前記第１の物体とは異なる第２の物体が前記第１の物体内に含まれていないと仮定した場合における、（i）光が前記第１の物体に入射してから出射するまでの間にとる複数の経路を表す複数の光路モデルと、（ii）前記複数の光路モデルに沿って前記第１の物体内を伝搬して前記第１の物体から出射した仮想的な光の所定面上における強度を示す第１の光強度分布とを取得する第１の取得ステップと、
    前記第１の物体に入射してから出射した光の強度分布であって前記所定面上において観測された光の強度分布である第２の光強度分布を取得する第２の取得ステップと、
    前記所定面上の位置ごとに、前記第１の光強度分布および前記第２の光強度分布間における光の強度の比較結果に基づいて、前記第１の物体内に前記第２の物体が存在しない非存在確率を示す投票値を算出する投票値算出ステップと、
    前記第１の物体の内部を複数の投票領域に仮想的に分割し、前記所定面上の位置ごとに、当該位置に到達する光の光路モデル上の投票領域に対して当該位置の前記投票値を投票する投票ステップと、
    前記複数の投票領域の各々に対して、前記投票ステップで得られた投票結果が閾値より小さい場合に、前記第２の物体が存在すると推定する推定ステップとを含む
    物体内部推定方法。
  11. 請求項１０に記載の物体内部推定方法をコンピュータに実行させるためのプログラム。
JP2013146966A 2012-08-13 2013-07-12 物体内部推定装置およびその方法 Active JP6212312B2 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2013146966A JP6212312B2 (ja) 2012-08-13 2013-07-12 物体内部推定装置およびその方法

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
JP2012179159 2012-08-13
JP2012179159 2012-08-13
JP2013146966A JP6212312B2 (ja) 2012-08-13 2013-07-12 物体内部推定装置およびその方法

Publications (2)

Publication Number Publication Date
JP2014055939A JP2014055939A (ja) 2014-03-27
JP6212312B2 true JP6212312B2 (ja) 2017-10-11

Family

ID=50065975

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2013146966A Active JP6212312B2 (ja) 2012-08-13 2013-07-12 物体内部推定装置およびその方法

Country Status (2)

Country Link
US (1) US9234747B2 (ja)
JP (1) JP6212312B2 (ja)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6126437B2 (ja) 2013-03-29 2017-05-10 キヤノン株式会社 画像処理装置および画像処理方法
CN107430263B (zh) 2014-12-23 2020-12-01 苹果公司 具有平均照明路径和平均收集路径的共焦检查系统
WO2016106350A1 (en) 2014-12-23 2016-06-30 Bribbla Dynamics Llc Confocal inspection system having non-overlapping annular illumination and collection regions
CN107209116B (zh) * 2014-12-23 2020-08-07 苹果公司 包括考虑样本内的光学路径长度的变化的光学检查系统和方法
JP6723510B2 (ja) * 2015-06-15 2020-07-15 株式会社リコー 光学モデルの選定方法、検出光量分布補正方法、光学モデルの補正方法及び光学検査装置
US10292590B2 (en) 2015-06-15 2019-05-21 Ricoh Company, Ltd. Optical examination method and optical examination device
WO2017040431A1 (en) 2015-09-01 2017-03-09 Bribbla Dynamics Llc Reference switch architectures for noncontact sensing of substances
AU2017253712B8 (en) 2016-04-21 2019-11-21 Apple Inc. Optical system for reference switching
CN116893160A (zh) 2017-09-29 2023-10-17 苹果公司 路径解析的光学采样架构
EP3752873A1 (en) 2018-02-13 2020-12-23 Apple Inc. Integrated photonics device having integrated edge outcouplers
EP4176304A1 (en) 2020-09-09 2023-05-10 Apple Inc. Optical system for noise mitigation

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5070455A (en) * 1989-11-22 1991-12-03 Singer Imaging, Inc. Imaging system and method using scattered and diffused radiation
JPH07200512A (ja) * 1993-09-13 1995-08-04 Ezel Inc 最適化問題解決装置
JP3433508B2 (ja) * 1993-12-01 2003-08-04 浜松ホトニクス株式会社 散乱吸収体計測方法及び散乱吸収体計測装置
JP3310782B2 (ja) 1994-07-14 2002-08-05 株式会社日立製作所 吸収物質濃度の空間分布画像化装置
JP3310784B2 (ja) 1994-08-01 2002-08-05 三菱重工業株式会社 ごみ焼却炉の２次空気吹込み方法
JP3276947B2 (ja) 1999-10-27 2002-04-22 株式会社日立製作所 吸収物質濃度空間分布画像装置
JP3730927B2 (ja) 2002-03-20 2006-01-05 株式会社日立製作所 吸収物質濃度空間分布画像化装置
JP3902999B2 (ja) * 2002-10-16 2007-04-11 日機装株式会社 光学的散乱特性推定装置およびその作動方法

Also Published As

Publication number Publication date
US9234747B2 (en) 2016-01-12
US20140043620A1 (en) 2014-02-13
JP2014055939A (ja) 2014-03-27

Similar Documents

Publication Publication Date Title
JP6212312B2 (ja) 物体内部推定装置およびその方法
CN108877892B (zh) 自动估算和减少计算机断层成像扫描中的散射
US10645259B2 (en) Machine learning based real-time radiation dose assessment
Thompson et al. Tissue perfusion measurements: multiple-exposure laser speckle analysis generates laser Doppler–like spectra
EP3180603B1 (en) Quantitative dark field imaging in tomography
Van de Sompel et al. Task-based performance analysis of FBP, SART and ML for digital breast tomosynthesis using signal CNR and Channelised Hotelling Observers
WO2014185078A1 (ja) Ｘ線ｃｔ画像処理方法，ｘ線ｃｔ画像処理プログラム及びｘ線ｃｔ画像装置
JP2013519452A (ja) 陽子コンピュータ断層撮影のためのシステム及び方法
US8290116B2 (en) Imaging apparatus including correction unit for scattered radiation
JP2016147061A (ja) マルチスペクトル光音響イメージングのためのデバイスおよび方法
US20160213345A1 (en) Methods and systems for estimating scatter
Leino et al. Perturbation Monte Carlo method for quantitative photoacoustic tomography
US20110268339A1 (en) Systems and methods for determining a location of a lesion in a breast
Niinimaki et al. Multiresolution parameter choice method for total variation regularized tomography
JP5924232B2 (ja) 散乱線補正方法、画像処理装置および断層撮影装置
Cierniak An analytical iterative statistical algorithm for image reconstruction from projections
WO2018082088A1 (zh) 一种用于锥束ct图像散射修正的阻挡光栅优化方法及装置
EP3803779B1 (en) Scatter correction for x-ray imaging
Premru et al. Monte Carlo simulation of radiation transfer in human skin with geometrically correct treatment of boundaries between different tissues
Jha Retrieving information from scattered photons in medical imaging
Funamizu et al. Simulation of spectral reflectances in human skin tissue using ray tracing and GPU-based Monte Carlo method
Rapsomanikis et al. A Monte Carlo simulator dedicated to a time-resolved optical tomographic modality based on the Henyey-Greenstein phase function
CN109561868B (zh) 用于采集暗场图像的体模设备、暗场成像系统和方法
Wong et al. Objective assessment and design improvement of a staring, sparse transducer array by the spatial crosstalk matrix for 3d photoacoustic tomography
Soloviev et al. Angularly selective mesoscopic tomography

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20160105

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20160930

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20161004

A521 Written amendment

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20161114

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20170411

A521 Written amendment

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20170501

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: 20170912

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20170915

R151 Written notification of patent or utility model registration

Ref document number: 6212312

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R151