(57)【要約】【課題】より精密に、対象物の内部の特性を解析する。【解決手段】解析装置は、対象物の形状情報を取得する形状情報取得部と、表面の少なくとも一部の領域に光が照射された対象物からの透過光または反射光を測定して得た測定情報を取得する測定情報取得部と、形状情報および測定情報に基づいて、対象物の内部の光学特性値を算出する解析部と、を備える。【選択図】図１
2019074366.tif 000002
【特許請求の範囲】
【請求項１】
対象物の形状情報を取得する形状情報取得部と、
表面の少なくとも一部の領域に光が照射された前記対象物からの透過光または反射光を測定して得た測定情報を取得する測定情報取得部と、
前記形状情報および前記測定情報に基づいて、前記対象物の内部の光学特性値を算出する解析部と、を備える解析装置。
【請求項２】
請求項１に記載の解析装置において、
前記解析部は、前記形状情報および前記測定情報に基づいて、前記対象物の少なくとも一部を体積要素に分割して数値計算を行うことにより物質中の光伝搬に関する方程式を解いた結果に基づき、前記対象物の内部の光学特性値を算出する解析装置。
【請求項３】
請求項１または２に記載の解析装置において、
前記解析部は、前記対象物の少なくとも一部と同一の形状を備え、一様な前記光学特性値を有する物体に、前記光を照射した場合または前記光を照射したと仮定した場合に得られる情報に基づいて、前記対象物の内部の光学特性値を算出する解析装置。
【請求項４】
請求項１から３までのいずれか一項に記載の解析装置において、
前記解析部は、前記形状情報に基づいて前記対象物の表面における前記光の照射位置を算出する解析装置。
【請求項５】
請求項１から４までのいずれか一項に記載の解析装置と、
前記光を照射する照射部と、
前記透過光または前記反射光を撮像する撮像部と、
を備える撮像システム。
【請求項６】
請求項５に記載の撮像システムにおいて、
前記照射部は、波長７００ｎｍ以上１０００ｎｍ以下の近赤外光を照射する撮像システム。
【請求項７】
請求項５または６に記載の撮像システムにおいて、
前記照射部は、光の進行方向に垂直な面での強度分布が縞状である前記光を照射する撮像システム。
【請求項８】
請求項５または６に記載の撮像システムにおいて、
前記照射部は、ライン状または点形状の前記光を前記対象物に走査して照射する撮像システム。
【請求項９】
請求項８に記載の撮像システムにおいて、
前記照射部は、前記光の強度を時間的に変調するか、前記光を走査する速度を調節することにより、前記対象物の所望の位置に所望の量の光を照射する撮像システム。
【請求項１０】
請求項８または９に記載の撮像システムにおいて、
ライン状または点形状の前記光を前記対象物に走査しながら複数枚取得した画像から、任意の照明パターンを合成する撮像システム。
【請求項１１】
請求項８から１０までのいずれか一項に記載の撮像システムにおいて、
前記照射部は、パルス光を照射し、
前記撮像部の露光時間は、前記パルス光が照射された前記対象物からの透過光または反射光の立ち上がりからピークまでの時間の２倍よりも短い撮像システム。
【請求項１２】
対象物の形状情報を取得することと、
前記対象物の表面の少なくとも一部の領域に光を照射することと、
前記光が照射されている前記対象物からの透過光または反射光を撮像することと、
前記形状情報および、撮像した前記透過光または前記反射光の強度情報に基づいて前記対象物の内部の光学特性の情報を算出することと、
を含む解析方法。
【請求項１３】
請求項１２に記載の解析方法において、
前記対象物は、前記対象物と比較して散乱性または吸収性の低い媒質の中に置かれており、
前記媒質は、気体または前記対象物の浸透圧以下の浸透圧を有する溶液である解析方法。
【請求項１４】
請求項１３に記載の解析方法において、
前記媒質は、空気、水、または生理食塩水である解析方法。
【請求項１５】
請求項１３または１４に記載の解析方法において、
前記媒質は、前記光に対して前記対象物の屈折率との差が所定内である液体であり、
前記光は、前記液体の内部で照射部から照射される解析方法。
【請求項１６】
請求項１２から１５までのいずれか一項に記載の解析方法において、
前記対象物の少なくとも一部と同一の形状を有する物体の表面の少なくとも一部の領域に光を照射し、前記光が照射された前記物体からの透過光または反射光を測定して得た測定情報を取得することを備える解析方法。
【請求項１７】
請求項１２から１６までのいずれか一項に記載の解析方法において、
前記対象物の表面の少なくとも一部の領域に対して、異なる複数の方向から前記光を照射する解析方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、解析装置、撮像システムおよび解析方法に関する。
【背景技術】
【０００２】
対象物の内部を含む構造の分析において、対象物に光を照射し、透過または反射した光を解析することが行われている（特許文献１参照）。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】国際公開第２０１２／１５３７６９号
【発明の概要】
【０００４】
本発明の第１の態様によると、解析装置は、対象物の形状情報を取得する形状情報取得部と、表面の少なくとも一部の領域に光が照射された前記対象物からの透過光または反射光を測定して得た測定情報を取得する測定情報取得部と、前記形状情報および前記測定情報に基づいて、前記対象物の内部の光学特性値を算出する解析部と、を備える。
本発明の第２の態様によると、撮像システムは、第１の態様の解析装置と、前記光を照射する照射部と、前記透過光または前記反射光を撮像する撮像部と、を備える。
本発明の第３の態様によると、解析方法は、対象物の形状情報を取得することと、前記対象物の表面の少なくとも一部の領域に光を照射することと、前記光が照射されている前記対象物からの透過光または反射光を撮像することと、前記形状情報および、撮像した前記透過光または前記反射光の強度情報に基づいて前記対象物の内部の光学特性の情報を算出することと、を含む。
【図面の簡単な説明】
【０００５】
【図１】一実施形態の解析装置および撮像システムの構成を示す概念図である。
【図２】図２（ａ）は、対象物に近赤外光を照射して測定をする場合の概念図であり、図２（ｂ）は、対象物のファントムに近赤外光を照射して測定をする場合の概念図である。
【図３】解析方法の流れを示すフローチャートである。
【図４】対象物を媒質中で測定する点を模式的に示す図である。
【図５】図５（ａ）は、対象物における照射光の受光位置と強度とを模式的に示す図であり、図５（ｂ）は、照射光の走査の態様を模式的に示す概念図である。
【図６】図６（ａ）は、対象物における照射光の受光位置と強度とを模式的に示す図であり、図６（ｂ）は、照射光の走査の態様を模式的に示す概念図である。
【図７】画像の取得方法の一例を説明するための概念図である。
【図８】対象物へ照射された入射パルスと出射パルスを示す概念図である。
【図９】対象物への近赤外光の照射の態様を示す概念図である。
【図１０】解析方法の流れを示すフローチャートである。
【図１１】反射光を撮像する場合の測定装置の構成を示す概念図である。
【発明を実施するための形態】
【０００６】
例えば特許文献１では、光ファイバーを対象物周囲に配置する等の課題がある。非接触にて、任意形状の対象物の内部を含む構造の分析を行うことができれば、様々な用途が広がる可能性がある。
以下では、適宜図面を参照しながら、一実施形態の解析装置、撮像システムおよび解析方法等について説明する。
【０００７】
図１は、本実施形態の撮像システム１の機能を模式的に示した図である。撮像システム１は、測定装置１０と、解析装置４０とを備える。測定装置１０は、対象物Ｓに光を照射する照射部２０と、対象物Ｓからの光を撮像する撮像部３０とを備える。照射部２０は、光源２１と、コンデンサレンズ２２と、照射側偏光フィルタ２３と、光変調器２４と、照射側レンズ２５と、走査ミラー２６とを備える。撮像部３０は、撮像側偏光フィルタ３１と、バンドパスフィルタ３２と、撮像側レンズ３３と、受光部３４とを備える。解析装置４０は、表示部４１と、通信部４２と、記憶部４３と、入力部４４と、制御部５０とを備える。制御部５０は、形状情報取得部５１と、測定情報取得部５２と、解析部５３と、画像作成部５４とを備える。
なお、図１では照射部２０が照射光Ｌ１を対象物Ｓに照射し、対象物Ｓからの透過光Ｌ２を撮像部３０により撮像する構成となっているが、後に詳述するような対象物Ｓからの反射光Ｌ３を撮像する構成（図１１参照）としてもよい。
【０００８】
解析装置４０は、測定装置１０の受光部３４で受光した光の情報に基づいて、対象物Ｓの内部の光学特性に関する情報を算出する。対象物Ｓは特に限定されないが、生体組織等の散乱性を有する対象が特に好ましい。測定装置１０において対象物Ｓは例えば空気等の任意の媒質中に配置して測定を行うことができるが、この媒質は、対象物Ｓより散乱性または吸収性の低い、気体、液体等が好ましい。解析装置４０が算出する上記光学特性に関する情報は、対象物Ｓの内部の光伝搬に関わるパラメータが好ましく、例えば吸収係数μａ、拡散係数Ｄ、あるいは等価散乱係数μ’s等である。
なお、対象物Ｓが生体の場合、複数波長λにおける吸収係数μa(λ)をそれぞれ算出することで、以下の式（１）を連立させて酸素化ヘモグロビン濃度 cHbO2、脱酸素化ヘモグロビン濃度 cHb を求めることができる。ここで、酸素化ヘモグロビンのモル吸光度εHbO2(λ)、脱酸素化ヘモグロビンのモル吸光度εHb(λ)の値は広く知られている。詳細は、ジェイクスらの文献（Steven L Jacques (2013) Phys. Med. Biol. P37-61）を参考にされたい。そして、これらの酸素化ヘモグロビン濃度cHbO2、脱酸素化ヘモグロビン濃度cHb を用いて、酸素飽和度 SO2 を求めることができる。
【数１】
2019074366.tif 000003
…（１）
なお、受光部３４で受光した光の情報を解析装置４０で解析できれば、測定装置１０と解析装置４０との配置の態様は特に限定されず、有線または無線の通信を介して情報をやりとりし遠隔に配置されてもよい。
【０００９】
光源２１は、発光ダイオード（以下、ＬＥＤと呼ぶ）等の光を出射する装置を含んで構成される。光源２１が出射する光の波長は、対象物Ｓに応じて適宜設定される。散乱媒質での測定を容易にする観点から近赤外光が好ましく、波長５００ｎｍ以上１６００ｎｍ以下の光がより好ましく、波長７８０ｎｍ以上９００ｎｍ以下の光がさらに好ましい。例えば、トータルヘモグロビンや酸素飽和度を計測する場合は、波長６００〜１０００ｎｍの可視領域から近赤外領域の光を用いることが好ましく、血糖値（グルコース濃度）を計測する場合は、波長７８０〜１０００ｎｍの近赤外領域の光や波長８〜１２μｍの中赤外領域の光を用いることが好ましい。また、樹脂で形成された物体の内部の欠陥を検出する場合は、波長５００〜１０００ｎｍの可視領域から近赤外領域の光を用いることが好ましい。光源２１が出射する光は、受光部３４で検出し、解析装置４０で解析することが可能であればその態様は特に限定されず、例えば可干渉性の光でも非干渉性の光でもよい。光源２１から出射された光はコンデンサレンズ２２で光源２１の像を変調器２４に結像する。このときコンデンサレンズ２２から出射した光は照射側偏光フィルタ２３を通過する。図１では、光源２１から出射される光の進路を一点鎖線で模式的に示した。
【００１０】
照射側偏光フィルタ２３は、光源２１からの光の特定の偏光成分のみを取り出して光変調器２４に出射する。照射側偏光フィルタ２３が透過させる偏光成分と撮像側偏光フィルタ３１が透過させる偏光成分とは、照射光Ｌ１の対象物Ｓによる透過光Ｌ２または反射光を、ノイズを少なく取り出すために所定の関係にある。例えば、照射側偏光フィルタ２３と撮像側偏光フィルタ３１とは、クロスニコル配置されている。これにより、信号対雑音比を高くすることができる。
【００１１】
光変調器２４は、照射側偏光フィルタ２３から出射された光を、時間的および／または空間的に所定のパターンの強度の光に変調して照射側レンズ２５に出射する。本実施形態では、光変調器２４は、空間的な強度分布が光の進行方向から見て縞状のパターンとなるような光を出射する。すなわち、光変調器２４は、対象物Ｓに照射される照射光Ｌ１の光軸に平行にＺ軸をとり、Ｚ軸に垂直にＸ軸およびＹ軸をとった場合、ＸＹ平面上において強度分布が縞状パターンの光が対象物Ｓに照射されるように入射した光を変調する。このように、光の進行方向に垂直な面での強度分布が縞状である照射光を用いることで本実施形態の解析方法を用いて効率的に対象物内部の光学特性を解析できる。本実施形態では、照射光Ｌ１の縞状パターンは、たとえばＸ軸に垂直（Ｙ軸に平行）な方向とし、Ｘ軸方向に沿って照射光Ｌ１の強度が正弦的に変化するものとする。すなわち、照射光Ｌ１の照射パターンは、対象物Ｓの表面∂Ｓを表す３次元ベクトルR＝（ｘR，ｙR，ｚR）上での照射光Ｌ１の光強度をＩ（R）として以下の式（２）で表現される。
【数２】
2019074366.tif 000004
…（２）
ここで、Ｉ０（R）は照射光Ｌ１の強度の最大値の半分の値を示す。ｆは正弦的に変化する際の空間周波数であり、δは任意の位相である。ｆは縞の一周期の幅の逆数に対応し、測定結果を解析して得られる光学特性値の分布に求める空間的解像度等に応じて適宜調節される。光変調器２４を出射した光は、照射側レンズ２５で屈折されて走査ミラー２６に入射する。
なお、照射光Ｌ１の縞状パターンの向き等の態様は特に限定されない。また、上記でいう、「縞状」とは、本実施形態の数値計算の方法が利用できればその態様は特に限定されないが、周期的、特に正弦的なパターンが好ましい。
【００１２】
走査ミラー２６は、ガルバノミラーまたはポリゴンミラー等を含んで構成され、光変調器２４で生成された所定のパターンの強度を有する光の出射方向を適宜調節して対象物Ｓに出射する。照射部２０は、対象物Ｓの表面の少なくとも一部の領域に照射光Ｌ１を照射する。本実施形態では、照射部２０は、照射光Ｌ１の一つの空間周波数ｆの縞について、少なくとも３つの異なる位相δ（式（２）参照）に対応する縞状パターンの照射光Ｌ１を対象物Ｓに照射する。
なお、一つの空間周波数ｆについてとる異なる位相の数は、特に限定されない。
【００１３】
照射部２０が、位相δ（式（２）参照）を１２０度ずつ変えて光の強度が空間周波数ｆで正弦的に変動する縞状パターンを対象物Ｓに照射した場合、対象物Ｓからの透過または反射した光の縞状パターンは、測定される対象物Ｓの表面∂Ｓを表す３次元ベクトルR＝（ｘR，ｙR，ｚR）からの光の強度をΦｍ（R,ｆ）として、例えば以下の式（３）のように表される。この式（３）は受光部３４から対象物Ｓを見たときの、対象物Ｓの表面∂Ｓ上での光強度分布を意味する。
【数３】
2019074366.tif 000005
…（３）
ここで、ΦＤＣ（R）およびΦＡＣ（R,f）はそれぞれ光の強度の平均値および正弦的に変動する変調成分を示し、電流の直流成分および交流成分になぞらえてＤＣおよびＡＣの添え字を付している。θAC（R,f）は、照射光Ｌ１からの光強度の縞状パターンに関する位相を示す。本実施形態のように空間周波数ｆに対し、１２０度ずつ異なる３つの位相δの照射光Ｌ１を対象物Ｓに照射する場合、ｍは０、１、２の３種類の値をとる。
【００１４】
上記式（３）からｍ＝０、１、２の３つの場合の式を連立させて、Φ０（R,ｆ）、Φ１（R,ｆ）、Φ２（R,ｆ）から、受光部３４から対象物Ｓの表面∂Ｓを見たときの光の強度のＡＣ変調成分ΦＡＣ（R,ｆ）および位相θAC（R,f）が以下の式（４）のように算出される。
【数４】
2019074366.tif 000006
…（４）
撮像部３０では、照射光Ｌ１が、対象物Ｓを透過した光Ｌ２または反射した光を撮像し、解析装置４０では、撮像部３０から送られた光の強度の情報を基に、各測定点での透過光Ｌ２または反射光のＡＣ変調成分ΦＡＣ（R,ｆ）と位相θAC（R,ｆ）を算出し、さらに対象物Ｓの内部の光学特性値を算出する。位相のθAC（R,ｆ）の測定は、周期的に強度が変動する縞の歪みを測定することに相当する。以下では、適宜、測定装置１０による透過光Ｌ２または反射光のこのような測定を位相シフト測定とよぶ。
なお、３つ以上の異なる位相δに対する位相シフト測定における計算方法の詳細は、マラカラらの文献（Daniel Malacara (2007)“Optical shop testing” Wiley-Interscience社）を参考にされたい。３つ以上の異なる位相δに対する位相シフト測定は、計測ノイズが多いときに有用であることが一般に知られている。
【００１５】
対象物Ｓからの透過光Ｌ２または反射光は、撮像側偏光フィルタ３１を透過した後、バンドパスフィルタ３２に入射する。バンドパスフィルタ３２は、透過光Ｌ２または反射光で想定される波長域に合わせ近赤外光の一部または全部の領域の光を透過させる。バンドパスフィルタ３２を透過した光は、撮像側レンズ３３で屈折されて受光部３４に集光する。
【００１６】
受光部３４は、ＣＣＤ等の撮像素子を含んで構成され、近赤外光を含む波長域に感度を有し、透過光Ｌ２または反射光を受光する。受光部３４は、異なる複数の空間周波数ｆについて、それぞれ３つの異なる位相の縞状パターンの強度を有する照射光Ｌ１を対象物Ｓに照射したときの、対象物Ｓからの透過光Ｌ２または反射光を受光する。対象物Ｓの各測定点に対応する光の強度のデータ、すなわち強度情報は、後述の解析装置４０の通信部４２若しくは入力部４４、または不図示の有線または無線通信等を介し、解析装置４０の制御部５０に送られる。
【００１７】
測定装置１０は、対象物Ｓと表面の形状が同じで、内部が既知の一様な光学特性を有するファントムに対し、上述の対象物Ｓに対して行ったと同様の測定を行う。ファントムは、少なくとも対象物Ｓの一部と同一の形状を備え、一様な光学特性値を有するものが好ましい。対象物Ｓが生体の一部である場合には、対応するファントムは、ポリオキシメチレン等の樹脂を含んで構成することができる。
なお、ファントムの光学特性については、解析可能であれば、特にその分布の態様は限定されない。
【００１８】
図２（ａ）は、測定装置１０において対象物Ｓを測定する際の模式図であり、図２（ｂ）は、測定装置１０においてファントムＰを測定する際の模式図である。照射部２０は、対象物Ｓと同様、ファントムＰに対しても、対象物Ｓへの照射と同様に、異なる複数の空間周波数ｆのそれぞれについて異なる複数の位相の縞状パターンの強度を有する照射光Ｌ１を照射し、ファントムＰからの透過光Ｌ２または反射光を撮像する。ファントムＰについての測定情報も対象物Ｓについての測定情報と同様、解析装置４０に送られる。後に詳述するが、解析装置４０は、ファントムＰの測定によって得られたパラメータを利用して対象物Ｓの内部の光学特性値を算出する。
【００１９】
解析装置４０の表示部４１（図１）は、液晶モニタ等の表示モニタを含んで構成され、測定装置１０の測定結果や、算出した対象物Ｓの内部における光学特性値の分布を示す三次元画像等の解析装置４０の解析結果等を表示する。
【００２０】
解析装置４０の通信部４２は、インターネット等の無線や有線接続により通信可能な通信装置を含んで構成され、解析装置４０の解析結果を送信したり、適宜必要なデータを送受信する。
【００２１】
解析装置４０の記憶部４３は、メモリやハードディスク等の不揮発性の記憶媒体で構成され、制御部５０に処理を実行させるプログラム、解析装置４０の解析結果等を記憶する。
【００２２】
解析装置４０の入力部４４は、各種ボタンおよび／またはタッチパネル等を含んで構成され、照射光Ｌ１が照射する縞状パターンの情報や、受光部３４の撮像条件その他の解析装置５０の動作に必要な情報を受け付ける。
【００２３】
解析装置４０の制御部５０は、ＣＰＵ等の処理装置を含んで構成され、解析装置４０を制御する動作の主体として機能する。制御部５０は、記憶部４３等に記憶されたプログラムを実行することにより、解析装置４０における対象物Ｓの光学特性値の算出を含む各種処理を行う。測定装置１０の照射部２０、撮像部３０の制御は、測定装置１０の不図示の制御部により行われる。
なお、解析装置４０の制御部５０が測定装置１０の制御を行う構成にしてもよい。
【００２４】
形状情報取得部５１は、通信部４２または入力部４４等を介して対象物Ｓの表面の形状情報を取得する。対象物Ｓの表面の形状は、例えば株式会社ニコン製ポータブル三次元形状計測装置Ｐ３Ｄ  ＮＣ−２３２３Ｓ等を用いて様々な状況で簡便に測定することができる。対象物Ｓの表面の形状情報の取得方法は特に限定されない。形状情報取得部５１は、得られた対象物Ｓの三次元の形状情報を保持したり、記憶部４３に記憶させたりする。
なお、測定された対象物Ｓの三次元形状情報に基づいて上述のファントムＰを製造してもよい。
【００２５】
測定情報取得部５２は、測定装置１０の受光部３４が受光した対象物Ｓからの透過光または反射光についての情報を、通信部４２若しくは入力部４４、または不図示の有線または無線通信等を介して取得する。測定情報取得部５２は、得られた対象物Ｓの測定情報を保持したり、記憶部４３に記憶させたりする。
【００２６】
解析部５３は、形状情報取得部５１が取得した対象物Ｓの形状情報と、測定情報取得部５２が取得した対象物Ｓの測定情報に基づいて、対象物Ｓの内部の光学特性値の分布を算出する。
【００２７】
解析部５３は、対象物Ｓの測定情報から、位相シフト解析により透過光Ｌ２または反射光のＡＣ変調成分ΦＡＣ（R,ｆ）を示す画像と位相θAC（R,ｆ）を示す画像とを抽出する。その後、解析部５３は、物質中の光伝搬を記述する方程式を解くことで、対象物Ｓの内部の光学特性値の分布を算出する。本実施形態では、解析部５３は、形状情報で示された対象物Ｓの表面上の照射光Ｌの照射位置を算出し、受光部３４が受光した透過光Ｌ２または反射光の情報に基づいて、対象物Ｓの表面∂Ｓ（対象物Ｓの表面∂Ｓ上を表す３次元ベクトルR＝（ｘR，ｙR，ｚR））および対象物Ｓの内部（対象物Ｓの内部を表す３次元ベクトルr＝（ｘ，ｙ，ｚ））について、以下の式（５）で示される拡散方程式を解くことで、光学特性値の分布を算出する。
【数５】
2019074366.tif 000007
…（５）
ここで、ｔは時間、Φ（ｒ,t）は光の強度、ベクトルsは光強度の流れの単位方向ベクトル、３次元ベクトル Ｊ （ｒ,ｓ,ｔ）は拡散放射流速、μａ（ｒ）は吸収係数、μ’ｓ（ｒ）は等価散乱係数、Ｄ（ｒ）は拡散係数、ｎは対象物Ｓの屈折率、ｃは光速であり、Ｉ（ｒ）は、式（２）に示した対象物Ｓの表面∂Ｓに照射されている照射光Ｌ１の光量強度分布である。光源は時刻ｔ＝０から照明されることを表すために、ステップ関数Θ（t）を用いた。本実施形態では、対象物Ｓの内部の光学特性値として、吸収係数μａ（ｒ）と拡散係数Ｄ（ｒ）、あるいは透過散乱係数μ’s（ｒ）を求める。受光部３４を用いて対象物Ｓを測定する場合、光の強度Φ（ｒ,t）に比例した画像が測定される。
なお、吸収係数μa（r）と等価散乱係数μ’ｓ（ｒ）から、式（５）の上から３番目の式を用いて拡散係数Ｄ（ｒ）が容易に導き出せるが、逆にこの式（５）の上から３番目の式を用いて拡散係数Ｄ（ｒ）と吸収係数μa（ｒ）から、等価散乱係数μ’ｓ（ｒ）を算出することもできる。
【００２８】
定常状態における式（５）を変形して、グリーン関数Ｇを用いて記述し、以下の式（６）を得る。この変形に関し、詳細はマーケルらの文献（Markel VA and Schotland JC (2001) J. Opt. Soc. Am. A, p1336-1347）を参考にされたい。
【数６】
2019074366.tif 000008
…（６）
ここで、δμａ（ｒ）および、δＤ（ｒ）はそれぞれ吸収係数μａ（ｒ）および拡散係数Ｄ（ｒ）の平均値μａ０、Ｄ０からのずれとしている。すなわち、μａ（ｒ）＝μａ０＋δμａ（ｒ）、Ｄ（ｒ）＝Ｄ０＋δＤ（ｒ）である。関数Ｇ（r1,r2）は、任意位置r2に点光源を置いたときの任意位置r1でのグリーン関数Ｇのことである。関数Ｇ０（r1,r2）は、δμａおよび、δＤが０とした場合の、任意位置r2に点光源を置いたときの任意位置r1でのグリーン関数Ｇのことである。３次元位置ベクトルRs＝（ｘR_s，ｙR_s，ｚR_s）は対象物Ｓの表面∂Ｓでの点光源の位置、すなわち対象物Ｓの表面∂Ｓで照射光Ｌ１が照射されている位置を示し、３次元ベクトルRｄ＝（ｘR_ｄ，ｙR_ｄ，ｚR_ｄ）は観測する対象物Ｓの表面∂Ｓの測定点の位置を示す。３次元ベクトルｒは、対象物Ｓの内部を示す。式（６）は、ダイソン方程式をボルン近似した方程式に相当する。
【００２９】
式（６）に以下の定義式（７）を用いて、以下の式（８）を得る。
【数７】
2019074366.tif 000009
…（７）
ここで式（２）のＩ０（R）を用いた。
【数８】
2019074366.tif 000010
…（８）
ここでΦＡＣ(Rd,f)、ΦDＣ(Rd)、θAC (Rd,f)は、対象物Ｓの表面∂Ｓ上での光のＡＣ変調成分、ＤＣ強度、ＡＣ成分の位相を示し、それぞれ式（３）に対応する。またΦＡＣ＿０(Rd,f) 、ΦDＣ＿０(Rd)、θ０(Rd,f)はファントムＰの表面∂Ｐ上での光のＡＣ変調成分、ＤＣ強度、ＡＣ成分の位相を示す。この式（８）の空間積分の範囲は、対象物Ｓの少なくとも一部の表面および内部の領域であり、解析部５３は当該領域を複数の体積要素に分割して数値計算する。
【００３０】
式（８）の左辺における、ΦＡＣ(Rd,f)およびθＡＣ(Rd,f)は観測する対象物Ｓからの透過光Ｌ２または反射光を測定することにより、ΦＡＣ＿０(Rd,f)およびθ０(Rd,f)は観測するファントムＰからの透過光Ｌ２または反射光を測定することにより実験的に求めることができる。式（８）の右辺のグリーン関数Ｇ０(Rd,r)は、位置 Rd に点光源を置いたときの対象物Ｓの内部位置rにおける光の強度Φ(r)をＦＤＴＤ（Ｆｉｎｉｔｅ−ｄｉｆｆｅｒｅｎｃｅ  ｔｉｍｅ−ｄｏｍａｉｎ）法や有限要素法（ＦＥＭ）等を用いて数値計算により算出することによって得ることができる。
式（８）の右辺におけるΦＡＣ＿０(r,f)およびθ０(r,f)は、ファントムＰを対象物とみなし、この対象物に縞状パターンの照明を行ったときのファントムＰ内部位置rの光量計算をＦＤＴＤ法や有限要素法等を用いて式（５）を解くことで数値計算により算出することができる。
【００３１】
式（８）の空間積分は数値計算に伴い離散化することで、以下の式（９）で表されるように行列計算の形にして線形方程式として吸収係数μａ（ｒ）および拡散係数Ｄ（ｒ）の値を求めることができる。このとき、各変数の空間的な解像度は解が導き出せるように適宜設定することができる。
【数９】
2019074366.tif 000011
…（９）
ここで各行列要素は、以下の式（１０）〜（１３）に示すように定義される。
【数１０】
2019074366.tif 000012
…（１０）
【数１１】
2019074366.tif 000013
…（１１）
【数１２】
2019074366.tif 000014
…（１２）
【数１３】
2019074366.tif 000015
…（１３）
ここで３次元ベクトルRｄi（i= 1〜m）は、受光部３４が測定する対象物Ｓの表面∂Ｓの測定点の位置をｍ箇所に離散化して指定したことを意味する。また３次元ベクトルｒi （i= 1〜ｋ）は、対象物Ｓの内部の位置をk箇所に離散化して指定したことを意味する。そして、ｆi（i= 1〜s）は用いるｓ種類の縞の空間周波数を表す。
なお、ベクトルΔΦＡＣは対象物ＳとファントムＰとの透過光Ｌ２または反射光の画像のＡＣ変調成分の差であり実測にて求められる。行列Ｆの要素Ｇ１、Ｇ２は、式（８）の積分を離散化して行列としたものに対応するから、上述の方法（式（１１）、式（１２）参照）によって、事前に数値計算しておくことができる。行列Ｇ１、Ｇ２が得られれば、双共役傾斜安定化法等を用いてこの方程式を解くことで対象物Ｓの内部の吸収係数μａ（ｒ）および拡散係数Ｄ（ｒ）、ならびにその他の光学特性値の分布を算出することができる。
【００３２】
ＦＤＴＤ法による数値計算では、数値が発散しないように安定条件が満たされる必要がある。本発明者らの研究の結果、拡散方程式をＦＤＴＤ法で解く場合の安定条件が見いだされた。式（５）の拡散方程式を、空間と時間に関してYee格子で離散化した式は、以下の式（１４）、（１５）のようになる。Yee格子に関しては、タフロフらの文献（Ａｌｌｅｎ  Ｔａｆｌｏｖｅ（２００５）“Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｅｌｅｃｔｒｏｄｙｎａｍｉｃｓ”Ａｒｔｅｃｈ  Ｈｏｕｓｅ社）を参考にされたい。
【数１４】
2019074366.tif 000016
…（１４）
【数１５】
2019074366.tif 000017
…（１５）
ここで拡散放射流速ベクトルＪ＝（ＪＸ，ＪＹ，ＪＺ）と光の強度Φに関して、離散化した時間変化を表すために上付き添え字mを用いて表した。またここで用いている係数は、以下の式（１６）のように定義し、またμｔｒ（ｉ，ｊ，ｋ）＝１／(３Ｄ)，ν＝ｃ／ｎとした。ここでYee格子は、対象物Ｓの内部の各座標軸の方向についてΔｘ、Δｙ、Δｚの離散間隔でサンプリングした。また時間発展は、時間間隔Δtでサンプリングして計算される。
【数１６】
2019074366.tif 000018
…（１６）
ここでΦm（ｉ，ｊ，ｋ）と（ＪmＸ（ｉ，ｊ，ｋ），ＪmＹ（ｉ，ｊ，ｋ），ＪmＺ（ｉ，ｊ，ｋ））は時間ステップm、格子点(i,ｊ,ｋ)での光の強度と拡散放射流速ベクトルである。またμa（ｉ，ｊ，ｋ）は格子点(i,ｊ,ｋ)でサンプリングされた吸収係数を表す。
そして、以下の式（１７）のように、拡散放射流速ベクトルＪと光の強度Φを平面波展開する。
【数１７】
2019074366.tif 000019
…（１７）
ここで平面波の任意の波数ベクトルを、ベクトルｕ＝（ｕｘ、ｕｙ、ｕｚ）とした。そして式（１７）を式（１４）、（１５）に代入し、任意の平面波で成り立つためにはΦm＋１＝ＺΦmが成り立たたなければならないことが分かる。ここで、以下の式（１８）のように表した。
【数１８】
2019074366.tif 000020
…（１８）
行列Ｚは、増幅行列と呼ばれる。この行列Ｚの固有値の絶対値が１以下である条件（フォン=ノイマンの安定性）から、数値計算が安定する以下の条件式（１９）が得られる。なお、フォン＝ノイマンの安定性の一般論については、ヒルシュらの文献（Charles Hirsch (1989) Numerical Computation of Internal and External Flows）を参考にされたい。
【００３３】
拡散方程式をＦＤＴＤ法により、対象の内部の各座標軸の方向についてΔｘ、Δｙ、Δｚの離散間隔でサンプリングし、時間発展についてΔｔのサンプリング間隔で数値計算する場合、以下の式（１９）の条件が安定に計算結果を得るために必要である。
【数１９】
2019074366.tif 000021
…（１９）
ここで、ｎは対象物Ｓの屈折率、ｃは光速である。解析部５３は、空間および時間に関し式（１９）の条件を常に満たすようにサンプリング間隔を設定してＦＤＴＤ法による数値計算を行う。
【００３４】
画像作成部５４（図１）は、得られた対象物Ｓの内部の吸収係数μａ（ｒ）や拡散係数Ｄ（ｒ）等の光学特性値の分布を示す画像を作成する。画像作成部５４は、形状情報に基づいて対象物Ｓの形状を示す画像を作成し、光学特性値の大きさを、色の違いや等高線等により区別して表示することができる。画像作成部５４は、形状情報に基づいて対象物Ｓの形状を示す三次元画像を作成し、光学特性値の大きさを、色の違いや光学特性値の値が等しい面等により区別して表示することができる。
【００３５】
図３は、本実施形態の解析装置４０を用いた対象物Ｓの内部の解析方法の流れを示すフローチャートである。ステップＳ１００１において、形状情報取得部５１は、対象物Ｓの三次元形状情報を取得する。ステップＳ１００１が終了したら、ステップＳ１００３に進む。
【００３６】
ステップＳ１００３において、測定装置１０の照射部２０は、対象物Ｓに近赤外光Ｌ１を照射し、透過光Ｌ２または反射光を撮像部３０が撮像し、解析装置４０の測定情報取得部５２は、近赤外光Ｌ１が照射された対象物Ｓからの透過光Ｌ２または反射光を測定して得た測定情報を取得する。ステップＳ１００３が終了したらステップＳ１００５に進む。
【００３７】
ステップＳ１００５において、解析部５３は、対象物Ｓの形状情報および測定情報、ならびにファントムＰに、近赤外光Ｌ１を照射した場合に得られる情報に基づいて、対象物Ｓの少なくとも一部を体積要素に分割して数値計算を行うことにより拡散方程式（５）を解く。すなわちＦＤＴＤ法や有限要素法等を用いて、式（８）で計算したときと同様にしてＧ0（Rd,r）、ΦAC_0（r,f）、θ0（r,f）を数値的に計算し、式（１１）、（１２）で用いている行列Ｇ1, Ｇ2を計算する。ステップＳ１００５が終了したら、ステップＳ１００７に進む。
【００３８】
ステップＳ１００７において、解析部５３は、ステップＳ１００５において拡散方程式を解いた結果に基づいて、対象物Ｓの内部の光学特性値を算出する。ステップＳ１００７が終了したら、ステップＳ１００９に進む。ステップＳ１００９において、表示部４１は、ステップＳ１００７で算出した光学特性値の分布を示す画像を表示する。ステップＳ１００９が終了したら、処理を終了する。
【００３９】
上述の実施の形態によれば、次の作用効果が得られる。
（１）本実施形態の解析装置４０は、対象物Ｓの形状情報を取得する形状情報取得部５１と、対象物Ｓの表面の少なくとも一部の領域に近赤外光Ｌ１を照射し、近赤外光Ｌ１が照射された対象物Ｓからの透過光Ｌ２または反射光を測定して得た測定情報を取得する測定情報取得部５２と、上記形状情報および上記測定情報に基づいて、対象物Ｓの内部の吸収係数μａ（ｒ）や拡散係数Ｄ（ｒ）等の光学特性値を算出する解析部５３と、を備える。これにより、様々な形状の対象物Ｓの内部の光学特性値を精密に算出することができる。
【００４０】
（２）本実施形態の解析装置４０において、解析部５３は、上記形状情報および上記測定情報に基づいて、対象物Ｓの少なくとも一部を体積要素に分割して数値計算を行うことにより光伝搬に関する方程式を解いた結果に基づき、対象物Ｓの内部の吸収係数μａ（ｒ）や拡散係数Ｄ（ｒ）等の光学特性値を算出する。これにより、対象物Ｓが直方体、円柱、球以外の任意の形状を有していても、本発明の拡散光トモグラフィ技術により精密に対象物Ｓの内部の光学特性値を測定することができる。
【００４１】
（３）本実施形態の解析装置４０において、解析部５３は、対象物Ｓの少なくとも一部と同一の形状を備え、一様な光学特性値を有するファントムＰに、近赤外光Ｌ１を照射した場合に得られる情報に基づいて、対象物Ｓの内部の光学特性値を算出する。これにより、任意の形状の対象物Ｓに対し、光伝搬に関する方程式に基づいて精密に対象物Ｓの内部の光学特性値を測定することができる。
【００４２】
（４）本実施形態の解析装置４０において、解析部５３は、対象物Ｓの形状情報に基づいて対象物Ｓの表面における近赤外光Ｌ１の照射位置を算出する。これにより、対象物Ｓを拡散方程式等を解くために理想的な形状の媒質で覆わなくても精密に対象物Ｓの内部の光学特性値を測定することができ、特に対象物Ｓを覆う媒質に白濁溶液を用いた際の光散乱による光量低下を防ぐことができる。
【００４３】
（５）本実施形態の撮像システム１は、解析装置４０を備え、照射部２０は、波長７００ｎｍ以上１０００ｎｍ以下の近赤外光Ｌ１を照射する。これにより、生体の一部等を対象としても、透過性が高いため深く、厚い対象物Ｓを解析することができる。
【００４４】
（６）本実施形態の解析方法は、対象物Ｓの形状情報を取得することと、対象物Ｓと比較して散乱性または吸収性の低い媒質の中に置かれた対象物Ｓの表面の少なくとも一部の領域に近赤外光Ｌ１を照射することと、近赤外光Ｌ１が照射されている対象物Ｓからの透過光Ｌ２または反射光を撮像することと、上記形状情報および、撮像した透過光Ｌ２または反射光の強度情報に基づいて対象物Ｓの内部の光学特性の情報を算出する。これにより、様々な形状の対象物Ｓの内部の光学特性値を精密に算出することができる。
【００４５】
次のような変形も本発明の範囲内であり、上述の実施形態と組み合わせることが可能である。以下の変形例において、上述の実施形態と同様の構造、機能を示す部位に関しては、同一の符号で参照し、適宜説明を省略する。
（変形例１）
上述の実施形態では、ファントムＰに対しても実際に測定装置１０により位相θ０等の測定を行っているが、ファントムＰに対して当該測定を行わず数値計算で代用してもよい。すなわち、解析部５３は、ファントムＰに照射光Ｌ１を照射した場合の透過光Ｌ２または反射光を撮像したファントムＰの表面∂Ｐからの変調成分ΦＡＣ＿０（R,f）および位相θ０（R,f）を、ＦＤＴＤ法や有限要素法等を用いて式（８）を計算したときと同様にして数値計算にて解くことで算出してもよい。解析部５３は、照射光Ｌ１の縞状パターンとファントムＰの形状、光学特性から、式（８）を計算したときと同様にして、ファントムＰ内部の光伝搬をシミュレートする等して、画像の変調成分ΦＡＣ＿０（r,f）および位相θ０（r,f）を求めることができる。
【００４６】
本変形例の解析装置１において、解析部５３は、対象物Ｓの少なくとも一部と同一の形状を備え、一様な光学特性値を有するファントムＰに、近赤外光Ｌ１を照射した場合または近赤外光Ｌ１を照射したと仮定した場合に得られる情報に基づいて、対象物Ｓの内部の光学特性値を算出する。これにより、対象物ＳのファントムＰが用意できなかったり、製造が難しい場合でも、精密に対象物Ｓの内部の光学特性値を算出することができる。
【００４７】
（変形例２）
上述の実施形態では、測定装置１０において、対象物Ｓを測定用の液体等の媒質中に配置せずに、直接対象物Ｓに照射光を照射した、すなわち照射光Ｌ１の経路が空気中であったが、対象物Ｓを、対象物Ｓと屈折率が等しいか、ほぼ等しい所定の範囲の屈折率を有する媒質中で測定してもよい。上記所定の範囲とは、対象物Ｓの屈折率または対象物Ｓの表面を構成する物質の屈折率との差が０．１以内、好ましくは０．０５以内の屈折率を有する媒質である。対象物Ｓがヒト等の生体の一部である場合、対象物Ｓは水や生理食塩水もしくは大豆油等の媒質で囲むことが好ましい。
【００４８】
図４は、対象物Ｓを所定の媒質で囲む場合を説明するための図である。対象物Ｓは、使用する近赤外光Ｌ１に対する対象物Ｓの屈折率とほぼ等しい屈折率を有する液体９０中に配置される。液体９０は、近赤外光Ｌ１を透過する不図示の容器により保持される。照射部２０から出射した近赤外光Ｌ１は、外気を通過した後に液体９０中を進み、対象物Ｓに入射する。対象物Ｓからの透過光Ｌ２または反射光は、液体９０、外気を通過し、撮像部３０へ入射する。
【００４９】
本変形例の撮像システム１では、対象物Ｓを囲む媒質は、近赤外光Ｌ１に対して対象物Ｓの屈折率との差が０．１以内、好ましくは０．０５以内である。これにより、対象物Ｓの表面での反射や散乱を抑え、効率よく位相シフト測定を行うことができる。
【００５０】
測定装置１０の照射光Ｌ１の照射口は、対象物Ｓを囲む液体媒質の内部に位置してもよい。この場合、近赤外光Ｌ１は、対象物Ｓを囲む液体の内部で照射部２０から照射される。これにより、上記液体と外気との界面による照射光Ｌ１の反射や散乱を防ぐことができ、効率よく位相シフト測定を行うことができる。
【００５１】
対象物Ｓの浸透圧と、対象物Ｓを囲む液体９０との浸透圧が等しいと、対象物Ｓと液体９０との水分等の移動が抑えられ、より安定的に長時間測定を行うことができる。対象物Ｓを囲む媒質９０は、対象物Ｓの浸透圧以下の浸透圧を有することが好ましく、対象物Ｓの浸透圧と媒質９０の浸透圧との差は対象物Ｓの浸透圧の１０％以内であることがより好ましい。
【００５２】
（変形例３）
上述の実施形態では、強度が空間的に正弦的に変化する縞状パターンの照射光Ｌ１を対象物Ｓに照射して位相シフト測定を行ったが、照射部２０は、ライン状または点形状の近赤外光Ｌ１を対象物Ｓに走査して照射する構成にしてもよい。これにより、近赤外光Ｌ１の走査パターンを時間的、空間的に変えることで様々な照射パターンを実現することができる。
【００５３】
図５は、ライン状または点状の光を対象物Ｓ上に走査させる例を説明するための図である。図５（ａ）は、対象物Ｓ上の走査方向に座標軸をとり、その座標軸上の位置を横軸に、照射光Ｌ１の強度を縦軸にグラフで模式的に示したものである。時刻ｔ１において照射光Ｌ１の照射が開始されると、対象物Ｓ上を照射光Ｌ１が走査する。このとき、対象物Ｓ上の各位置に所望の強度で照射光Ｌ１が照射されるように、照射光Ｌ１の強度が変調される。例えば、時刻ｔ２では照射光Ｌ１の強度が照射開始時よりも小さくなっている。時刻ｔ３で走査が終了し、照射光Ｌ１の照射が終了する。すなわち、この場合、時刻ｔ１に受光部３４のシャッターを解放し、露光を継続しながら照射光Ｌ１を走査し、時刻ｔ３に当該シャッターを閉じる。
【００５４】
図５（ｂ）は、照射光Ｌ１を走査している間の走査ミラー２５の動きを模式的に示す図である。走査ミラー２５は対象物Ｓ上を一定の速度で照射光Ｌ１の照射位置が移動するように照射し、対象物Ｓ上の走査している位置に所望の強度で光が当たるよう、照射光Ｌ１は強度変調される。
【００５５】
図６は、ライン状または点形状の光を対象物Ｓ上に走査させる他の例を説明するための図である。図６（ａ）は、対象物Ｓ上の走査方向に座標軸をとり、その座標軸上の位置を横軸に、照射光Ｌ１の強度を縦軸にグラフで模式的に示したものである。時刻ｔ１において照射光Ｌ１の照射が開始されると、対象物Ｓ上を照射光Ｌ１が走査する。このとき、対象物Ｓ上の各位置に所望の強度で照射光Ｌ１が照射されるように、照射光Ｌ１の走査速度が変調される。例えば、時刻ｔ２では照射光Ｌ１の走査速度Ｖ２が照射開始時の走査速度Ｖ１よりも大きくなっており、これにより時刻ｔ２で走査される対象物Ｓ上の位置は照射時間が短いため照射される照射光Ｌ１の累積的な強度は小さくなる。時刻ｔ３で走査が終了し、照射光Ｌ１の照射が終了する。すなわち、この場合も、時刻ｔ１に受光部３４のシャッターを解放し、露光を継続しながら照射光Ｌ１を走査し、時刻ｔ３に当該シャッターを閉じる。
【００５６】
図６（ｂ）は、照射光Ｌ１を走査している間の走査ミラー２５の動きを模式的に示す図である。走査ミラー２５は対象物Ｓ上を制御された速度で照射光Ｌ１の照射位置が移動するように照射し、対象物Ｓ上の走査している位置に所望の強度で光が当たるよう、照射光Ｌ１の走査速度は変調される。
なお、照射光Ｌ１の強度および走査速度の両方を変調して走査してもよい。
【００５７】
本変形例の撮像システム１において、照射部２０は、近赤外光Ｌ１の強度を時間的に変調するか、近赤外光Ｌ１を走査する速度を調節することにより、対象物Ｓの所望の位置に所望の量の光を照射する。これにより、ムラが無く精密に、様々な縞状パターンの照射光Ｌ１を高速に照射することができる。
【００５８】
（変形例４）
上述の変形例３では、受光部３４を露光したまま、照射光Ｌ１を走査させたが、照射光Ｌ１の走査中、複数回に分けて露光し、各回の露光で得られた画像を、後から画像処理で足し合わせて構成してもよい。画像処理では、適宜、画像を足し合わせる際に取得したそれぞれの画像の重み付けを変えてもよい。
【００５９】
図７は、離散的に画像を取得し、後で重みをつけてパターン化する点を説明するための図である。時刻ｔ１と、時刻ｔ２と、時刻ｔ３において、同じ強度で対象物Ｓを点状またはライン状に照射している。このとき、Ｘ軸に沿ってＩＤＣ＋ＩＡＣ・ｃｏｓ（２πfｘ+δ）の強度で照射した結果となる画像を合成したい場合には、ＩＤＣ＋ＩＡＣ・ｃｏｓ（２πfｘ+δ）をａ１・Ｉ１＋ａ２・Ｉ２＋ａ３・Ｉ３＋…と分解した場合の係数ａg（g＝１，２，３…）を、式（２）のＩに対応する照射光Ｌ１で得られた画像に重み付けして足し合わせることで実現できる。
【００６０】
本変形例の撮像システム１では、ライン状または点形状の近赤外光Ｌ１を対象物Ｓに走査しながら複数枚取得した画像から、任意の照明パターンを合成する。これにより、物理的に照明パターンを実現しなくても画像処理を利用して測定結果を実質的に得ることができ、また強度ムラのない照射を行うことができる。
【００６１】
（変形例５）
上述の実施形態において、照射光Ｌ１をピコ秒またはフェムト秒パルス光とし、受光部３４が透過光Ｌ２または反射光を受光する際、パルス光の立ち上がりから所定の時間内のみ露光し、スネーク光を取得する構成にしてもよい。これにより、透過光Ｌ２または反射光のパルス光のうちで、散乱の影響が小さい部分を取得することができ、高コントラストな画像を取得することができる。
【００６２】
図８は、スネーク光を説明するための図である。図８左のグラフは、ピコ秒またはフェムト秒パルスによる照射光Ｌ１の波形を、時間を横軸、振幅を縦軸にして示している。光源２１は、ピコ秒またはフェムト秒パルスレーザーを用いる。
【００６３】
図８右のグラフは、ピコ秒またはフェムト秒パルスによる照射光Ｌ１を対象物Ｓに照射したときの、対象物Ｓからの透過光Ｌ２または反射光の波形を、時間を横軸、振幅を縦軸にして示している。透過光Ｌ２または反射光では、照射光Ｌ１に比べて、散乱の影響を受けて遅れて到達した光により波形が時間方向に伸びていることがわかる。撮像部３０の露光時間ｗは、透過光Ｌ２または反射光のパルス光の立ち上がりからピークまでの時間をｔｒとすると、ｔｒの２倍より短いことが好ましく、ｔｒより短いことがさらに好ましい。露光時間ｗを制御して受光するため、受光部３４は、ＩＣＣＤ（Ｉｎｔｅｎｓｉｆｉｅｄ  ＣＣＤ）、もしくはＴｏＦ（Ｔｉｍｅ  ｏｆ  Ｆｌｉｇｈｔ）カメラ、またはＥＭＣＣＤ（Ｅｌｅｃｔｒｏｎ  Ｍｕｌｔｉｐｌｙｉｎｇ  ＣＣＤ）を含んで構成され、露光時間は照射時間と連動して制御される。パルス光による位相シフト測定の解析においては、上述または後述の実施形態または変形例で示された解析方法、例えば式（９）、（２０）、（２８）、（２９）、（３０）または（３１）に基づく数値計算を用いる方法に加え、ラドン変換による吸収像のイメージング等をすることができる。
【００６４】
本変形例の撮像システム１において、照射部２０は、パルス光を照射し、撮像部３０の露光時間は、当該パルス光が照射された対象物Ｓからの透過光Ｌ２または反射光の立ち上がりからピークまでの時間の２倍よりも短い。これにより、撮像部３０は、散乱の影響の少ない高コントラストの画像を取得することができる。
【００６５】
（変形例６）
上述の実施形態では、一方向から照射光Ｌ１を照射して位相シフト測定を行う構成にしたが、複数の方向から複数の照射部２０を用いて照射光Ｌ１を照射し、対象物Ｓからの透過光Ｌ２または反射光を測定して位相シフト測定を行う構成にしてもよい。
【００６６】
図９は、２方向から２つの照射部２０により照射光Ｌ１を照射する例を模式的に示す図である。第１照射部２０−１は対象物Ｓに第１照射光Ｌ１−１を照射し、第１撮像部３０−１がその透過光Ｌ２−１を撮像する。第２照射部２０−２は第１照射部２０−１とは異なる方向から対象物Ｓに第２照射光Ｌ１−２を照射し、第２撮像部３０−２はその透過光Ｌ２−２を撮像する。
【００６７】
それぞれ異なる２方向から照射光Ｌ１を照射して位相シフト測定を２回行った場合、上述の式（９）と同様の方法にて、以下の式（２０）を用いて吸収係数μaおよび拡散係数Ｄの値を求めることができる。
【数２０】
2019074366.tif 000022
…（２０）
ここで各行列要素は、以下の式（２１）〜（２７）式のように定義される。
【数２１】
2019074366.tif 000023
…（２１）
【数２２】
2019074366.tif 000024
…（２２）
【数２３】
2019074366.tif 000025
…（２３）
【数２４】
2019074366.tif 000026
…（２４）
【数２５】
2019074366.tif 000027
…（２５）
【数２６】
2019074366.tif 000028
…（２６）
【数２７】
2019074366.tif 000029
…（２７）
ここで、３次元ベクトルRdi_1（i= 1〜q）は第１方向についての測定に関する対象物Ｓの表面∂Ｓの測定点の位置をq箇所に離散化して指定したものであり、３次元ベクトルｒdi_2（i= １〜q）は第２方向についての測定に関する対象物Ｓの表面の測定点の位置をq箇所に離散化して指定したものである。また３次元ベクトルｒi （i = 1〜ｋ）は、対象物Ｓ内部の位置をｋ箇所に離散化して指定したことを意味する。そして、ｆi（i =1〜s）は用いるｓ種類の縞の空間周波数を表す。
また、ΔΦＡＣ＿１は、第１撮像部３０−１が撮像した透過光Ｌ２−１の、対象物ＳとファントムＰとの画像のＡＣ変調成分の差であり、ΔΦＡＣ＿２は、第２撮像部３０−２が撮像した透過光Ｌ２−２の、対象物ＳとファントムＰとの画像のＡＣ変調成分の差であり、これらは２回の測定によってそれぞれ求まる。Ｆ１およびＦ２は、それぞれ第１および第２撮像部について、式（９）におけるＦと同様にして、式（５）の拡散方程式を数値計算することによって求めることができる。
なお、対象物Ｓを支持する載置台を回転させることで、対象物Ｓに対し複数の方向から照射光Ｌ１を照射する構成としてもよい。
【００６８】
本変形例の解析方法では、対象物Ｓの表面の少なくとも一部の領域に対して、異なる複数の方向から近赤外光Ｌ１を照射する。これにより、一方向から位相シフト測定を行うのに比べ、より広い範囲から情報を取得してより精密に光学特性値を算出することができる。
【００６９】
（変形例７）
上述の実施形態では、グリーン関数を用いて拡散方程式（式（５））を変形し、縞状パターンの空間周波数ｆや位相δを変えた複数の測定結果をそれぞれあてはめ、得られた式を連立させて数値計算により光学特性値を算出した。しかしながら、解析部５３は、得られた実験結果に最も適合する光学特性値を反復的なフォワード計算により算出してもよい。
【００７０】
拡散方程式（式（５））は、吸収係数μａ（ｒ）、等価散乱係数μ’ｓ（ｒ）（または拡散係数Ｄ（ｒ））、および照射光Ｌ１の縞状パターン、強度が求まれば、対象物Ｓの表面∂Ｓからの光の強度のＡＣ変調成分ΦＡＣ（Rd,f）および位相θAC（Rd,f）を算出することができる。この算出をここでは、フォワード計算という。そして、解析部５３は、測定された対象物Ｓの表面∂Ｓからの光の強度のＡＣ変調成分ΦＡＣ_exp（Rd,f）および位相θAC_exp（rｄ,f）に基づいた以下の式（２８）で表される評価関数Ｈが最小値をとるように、吸収係数μａ（ｒ）、等価散乱係数μ’ｓ（ｒ）を最適化する。
【数２８】
2019074366.tif 000030
…（２８）
ここで、３次元ベクトルRｄは、測定に関する対象物Ｓの表面∂Ｓ上の位置である。
実際の測定では、検出器の量子効率や露出時間を使って、検出器で得られた強度画像を光の強度ΦAC_exp（Rｄ,f）に直す必要がある。このとき、検出器の量子効率が撮像素子内で空間的にばらついていることが測定誤差となってしまう恐れがある。これを避けるために光の強度を規格化して扱う、以下の式（２９）のような評価関数Ｈが最小値をとるように、吸収係数μａ（ｒ）、等価散乱係数μ’ｓ（ｒ）を最適化してもよい。
【数２９】
2019074366.tif 000031
…（２９）
ここで、ΦＡＣ_exp_0（Rｄ,f）は、ファントムＰに縞状パターンを照射したときの表面∂Ｐからの透過光または反射光を測定することによって得られる光強度のＡＣ変調成分である。ΦＡＣ_0（Rｄ,f）はフォワード計算によって得られるファントムＰの表面からの透過光または反射光のＡＣ変調成分である。このように光の強度画像を規格化することによって、光の強度を直接使わずに最適化計算ができるようになる。
なお、フォワード計算とは対象物Ｓの内部の吸収係数や等価散乱係数の空間分布が与えられたときに、照明光Ｌ１の縞パターンを照明したときに対象物Ｓの表面∂Ｓからの光の強度分布Φ(r)を拡散方程式（式（５））を解くことによって計算することを言う。また、拡散方程式（式（５））でフォワード計算が可能であれば、最適化する光学特性値は特に限定されない。
【００７１】
解析部５３は、公知のＬｅｖｅｎｂｅｒｇ−Ｍａｒｑｕａｒｄｔ法等の最適化計算のためのアルゴリズムを用いて、効率よく最適化を行うことができる。
なお、上述の変形例のように対象物Ｓにそれぞれ異なる２方向から照射光Ｌ１を照射し２回測定を行う場合には、以下の式（３０）で表される評価関数Ｈが最小値をとるように、光学特性値を最適化する。
【数３０】
2019074366.tif 000032
…（３０）
ここで、３次元ベクトルRｄ１は第１方向についての測定に関する対象物Ｓの表面∂Ｓ上の位置であり、３次元ベクトルRｄ２は第２方向についての測定に関する対象物Ｓの表面∂Ｓ上の位置である。
式（２９）と同様に実際の測定では、検出器の量子効率や露出時間を使って、検出器で得られた強度画像を光の強度ΦAC_exp（Rｄ,f）に直す必要がある。このとき、検出器の量子効率が撮像素子内で空間的にばらついていることが測定誤差となってしまう恐れがある。これを避けるために光の強度を規格化して扱う、以下の式（３１）ような評価関数Ｈが最小値をとるように、光学特性値を最適化してもよい。
【数３１】
2019074366.tif 000033
…（３１）
ΦＡＣ_exp_0（Rｄ,f）は、ファントムＰの表面∂Ｐからの透過光または反射光を測定することによって得られる光強度のＡＣ変調成分である。ΦＡＣ_0（Rｄ,f）はフォワード計算によって得られる、ファントムＰの表面∂Ｐからの透過光または反射光のＡＣ変調成分である。このように光の強度画像を規格化することによって、光の強度を直接使わずに最適化計算ができるようになる。
なお、拡散方程式（式（５））でフォワード計算が可能であれば、最適化する光学特性値は特に限定されない。
【００７２】
図１０は、本変形例の解析装置４０を用いた対象物Ｓの内部の解析方法の流れを示すフローチャートである。ステップＳ２００１からステップＳ２００３までは図３のフローチャートのステップＳ１００１からステップＳ１００３までと同じであるため説明を省略する。ステップＳ２００３が終了したら、ステップＳ２００５に進む。
【００７３】
ステップＳ２００５において、解析部５３は、対象物Ｓの形状情報および測定情報、ならびに拡散方程式から算出される評価関数Ｊに基づいて光学特性値を最適化する。ステップＳ２００５が終了したら、ステップＳ２００７に進む。ステップＳ２００７において、表示部４１は、ステップＳ２００５において最適化した光学特性値の分布を示す画像を表示する。ステップＳ２００７が終了したら、処理を終了する。
【００７４】
（反射光を撮像する場合の構成）
上述の実施形態において、反射光を撮像する構成の場合には、反射光を直接撮像したり、ハーフビームスプリッター等を用いて反射光を分離して位相シフト測定を行う。
図１１は、反射光を測定する場合の測定装置１０の構成を示す図である。走査ミラー２６から出射された入射光Ｌ１はハーフビームスプリッター２７を透過して対象物Ｓに入射する。入射光Ｌ１を対象物Ｓに照射したときの、対象物Ｓからの反射光Ｌ３はハーフビームスプリッターで反射され、撮像側偏光フィルタ３１に入射する。図１１において、透過光を測定する場合（図１）と同様の部分については、同一の符号で参照して説明を省略する。
【００７５】
本発明は上記実施形態の内容に限定されるものではない。本発明の技術的思想の範囲内で考えられるその他の態様も本発明の範囲内に含まれる。
【符号の説明】
【００７６】
１…撮像システム、１０…測定装置、２０…照射部、３０…撮像部、３４…受光部、４０…解析装置、５０…制御部、５１…形状情報取得部、５２…測定情報取得部、５３…解析部、５４…画像作成部、Ｓ…対象物、Ｐ…ファントム、Ｌ１…照射光、Ｌ２…透過光。
【図１】
2019074366.tif 000034
【図２】
2019074366.tif 000035
【図３】
2019074366.tif 000036
【図４】
2019074366.tif 000037
【図５】
2019074366.tif 000038
【図６】
2019074366.tif 000039
【図７】
2019074366.tif 000040
【図８】
2019074366.tif 000041
【図９】
2019074366.tif 000042
【図１０】
2019074366.tif 000043
【図１１】
2019074366.tif 000044