(57)【要約】【課題】二光束干渉計において二光束の間の光路長差が時間的に変動している状態であっても対象物の位相画像を取得することができる装置および方法を提供する。【解決手段】干渉画像取得部２は、二光束干渉計２０において二光束の間の光路長差が時間的に変動している状態で対象物についてＭ個の干渉画像を取得する。演算部３は、これらＭ個の干渉画像それぞれについて画像全体のオフセット位相シフト量を推定し、Ｍ個の干渉画像のうちからオフセット位相シフト量が互いに異なるＮ個の干渉画像を含む組を設定し、その組に含まれるＮ個の干渉画像に基づいて位相画像を生成する。【選択図】図１
2019124679.tif 000002
【特許請求の範囲】
【請求項１】
光源と、前記光源から出力された光を二光束に分岐した後に該二光束を合波して干渉光を出力する二光束干渉計と、前記二光束干渉計から出力された干渉光を受光する撮像器とを含み、前記二光束干渉計における何れかの光束の光路上に配置された対象物の干渉画像を前記撮像器により撮像して取得する干渉画像取得部と、
前記干渉画像取得部により取得された干渉画像に基づいて前記対象物の位相画像を生成する演算部と、
を備え、
前記干渉画像取得部は、前記二光束干渉計において前記二光束の間の光路長差が時間的に変動している状態で前記対象物についてＭ個の干渉画像を取得し、
前記演算部は、
前記Ｍ個の干渉画像それぞれについて画像全体のオフセット位相シフト量を推定し、
前記Ｍ個の干渉画像のうちから前記オフセット位相シフト量が互いに異なるＮ個の干渉画像を含む組を設定し、
前記組に含まれる前記Ｎ個の干渉画像に基づいて位相画像を生成する、
位相画像取得装置（ただし、３≦Ｎ≦Ｍ）。
【請求項２】
前記演算部は、
前記組を複数設定して、これら複数の組それぞれについて各組に含まれる前記Ｎ個の干渉画像に基づいて位相画像を生成し、
前記複数の組それぞれで生成された位相画像についてオフセット位相シフト量を補正した後に積算する、
請求項１に記載の位相画像取得装置。
【請求項３】
前記演算部は、前記Ｍ個の干渉画像をオフセット位相シフト量について昇順または降順に並び替えた後に前記複数の組に順に振り分ける、
請求項２に記載の位相画像取得装置。
【請求項４】
前記演算部は、前記位相画像を生成するために干渉画像に対して乗算される係数を反復的な最適化計算によって調整する、
請求項１〜３の何れか１項に記載の位相画像取得装置。
【請求項５】
光源と、前記光源から出力された光を二光束に分岐した後に該二光束を合波して干渉光を出力する二光束干渉計と、前記二光束干渉計から出力された干渉光を受光する撮像器とを含み、前記二光束干渉計における何れかの光束の光路上に配置された対象物の干渉画像を前記撮像器により撮像して取得する干渉画像取得部と、
前記干渉画像取得部により取得された干渉画像に基づいて前記対象物の位相画像を生成する演算部と、
を備え、
前記干渉画像取得部は、前記二光束干渉計において前記二光束の間の光路長差が時間的に変動している状態で前記対象物についてＭ個の干渉画像を取得し、
前記演算部は、
前記Ｍ個の干渉画像それぞれについてオフセット位相シフト量の分布を推定し、
各画素位置について、前記Ｍ個の干渉画像それぞれの当該画素のうちから前記オフセット位相シフト量が互いに異なるＮ個の画素を含む組を設定し、
各画素位置について、前記組に含まれる前記Ｎ個の画素に基づいて位相値を求め、
各画素位置の位相値に基づいて位相画像を生成する、
位相画像取得装置（ただし、３≦Ｎ≦Ｍ）。
【請求項６】
前記演算部は、
前記組を複数設定して、これら複数の組それぞれについて各組に含まれる前記Ｎ個の画素に基づいて位相値を求め、
前記複数の組それぞれで求めた位相値についてオフセット位相シフト量を補正した後に積算し、
各画素位置の位相値の積算値に基づいて位相画像を生成する、
請求項５に記載の位相画像取得装置。
【請求項７】
前記演算部は、各画素位置について、前記Ｍ個の干渉画像それぞれの当該画素をオフセット位相シフト量について昇順または降順に並び替えた後に前記複数の組に順に振り分ける、
請求項６に記載の位相画像取得装置。
【請求項８】
前記二光束干渉計は、前記二光束の双方または何れか一方の光路上に配置された液体の液面の揺らぎにより、前記二光束の間の光路長差の時間的な変動を生じさせる、
請求項１〜７の何れか１項に記載の位相画像取得装置。
【請求項９】
前記光源はインコヒーレント光を出力する、
請求項１〜８の何れか１項に記載の位相画像取得装置。
【請求項１０】
光源と、前記光源から出力された光を二光束に分岐した後に該二光束を合波して干渉光を出力する二光束干渉計と、前記二光束干渉計から出力された干渉光を受光する撮像器とを含む干渉画像取得部を用い、前記二光束干渉計において前記二光束の間の光路長差が時間的に変動している状態で、前記二光束干渉計における何れかの光束の光路上に配置された対象物についてＭ個の干渉画像を前記撮像器により撮像して取得する干渉画像取得ステップと、
前記Ｍ個の干渉画像それぞれについて画像全体のオフセット位相シフト量を推定するオフセット位相シフト量推定ステップと、
前記Ｍ個の干渉画像のうちから前記オフセット位相シフト量が互いに異なるＮ個の干渉画像を含む組を設定する組設定ステップと、
前記組に含まれる前記Ｎ個の干渉画像に基づいて位相画像を生成する位相画像生成ステップと、
を備える位相画像取得方法（ただし、３≦Ｎ≦Ｍ）。
【請求項１１】
前記組設定ステップにおいて前記組を複数設定し、
前記位相画像生成ステップにおいて前記複数の組それぞれについて各組に含まれる前記Ｎ個の干渉画像に基づいて位相画像を生成し、
前記複数の組それぞれで生成された位相画像についてオフセット位相シフト量を補正した後に積算する積算ステップを更に備える、
請求項１０に記載の位相画像取得方法。
【請求項１２】
前記組設定ステップにおいて、前記Ｍ個の干渉画像をオフセット位相シフト量について昇順または降順に並び替えた後に前記複数の組に順に振り分ける、
請求項１１に記載の位相画像取得方法。
【請求項１３】
前記位相画像生成ステップにおいて、前記位相画像を生成するために干渉画像に対して乗算される係数を反復的な最適化計算によって調整する、
請求項１０〜１２の何れか１項に記載の位相画像取得方法。
【請求項１４】
光源と、前記光源から出力された光を二光束に分岐した後に該二光束を合波して干渉光を出力する二光束干渉計と、前記二光束干渉計から出力された干渉光を受光する撮像器とを含む干渉画像取得部を用い、前記二光束干渉計において前記二光束の間の光路長差が時間的に変動している状態で、前記二光束干渉計における何れかの光束の光路上に配置された対象物についてＭ個の干渉画像を前記撮像器により撮像して取得する干渉画像取得ステップと、
前記Ｍ個の干渉画像それぞれについてオフセット位相シフト量の分布を推定するオフセット位相シフト量推定ステップと、
各画素位置について、前記Ｍ個の干渉画像それぞれの当該画素のうちから前記オフセット位相シフト量が互いに異なるＮ個の画素を含む組を設定する組設定ステップと、
各画素位置について、前記組に含まれる前記Ｎ個の画素に基づいて位相値を求める位相値算出ステップと、
各画素位置の位相値に基づいて位相画像を生成する位相画像生成ステップと、
を備える位相画像取得方法（ただし、３≦Ｎ≦Ｍ）。
【請求項１５】
前記組設定ステップにおいて前記組を複数設定し、
前記位相値算出ステップにおいて前記複数の組それぞれについて各組に含まれる前記Ｎ個の画素に基づいて位相値を求め、
前記複数の組それぞれで求めた位相値についてオフセット位相シフト量を補正した後に積算する積算ステップを更に備え、
前記位相画像生成ステップにおいて各画素位置の位相値の積算値に基づいて位相画像を生成する、
請求項１４に記載の位相画像取得方法。
【請求項１６】
前記組設定ステップにおいて、各画素位置について、前記Ｍ個の干渉画像それぞれの当該画素をオフセット位相シフト量について昇順または降順に並び替えた後に前記複数の組に順に振り分ける、
請求項１５に記載の位相画像取得方法。
【請求項１７】
前記干渉画像取得ステップにおいて、前記二光束干渉計における前記二光束の双方または何れか一方の光路上に配置された液体の液面の揺らぎにより、前記二光束の間の光路長差の時間的な変動を生じさせる、
請求項１０〜１６の何れか１項に記載の位相画像取得方法。
【請求項１８】
前記光源はインコヒーレント光を出力する、
請求項１０〜１７の何れか１項に記載の位相画像取得方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、対象物の位相画像を取得する装置および方法に関するものである。
【背景技術】
【０００２】
特許文献１に、対象物の干渉画像および位相画像を取得することができる装置および方法の発明が開示されている。この文献に記載された装置は、光源と、この光源から出力された光を二光束に分岐した後に該二光束を合波して干渉光を出力する二光束干渉計と、この二光束干渉計から出力された干渉光を受光する撮像器とを備える。この装置は、二光束干渉計における二光束の間の光路長差を互いに異なる複数の設定値それぞれで安定化した状態とし、各状態において対象物の干渉画像を撮像器により取得し、これら取得した複数の干渉画像に基づいて位相画像を求める。
【０００３】
従来の位相画像取得装置は、二光束干渉計における二光束の間の光路長差を安定化する為に、この光路長差を検出するとともに、この光路長差が設定値となるように、二光束のうちの何れかの光路上に配置したミラーの位置をフィードバック制御する。
【０００４】
しかし、例えば、培養液に入れられている対象物としての細胞の位相画像を求める場合、培養液の液面の揺らぎに因り、その培養液を透過する光の光路長は時間的に変動する。その液面の変動の幅は１波長分以上であり、その変動の周波数は数十Ｈｚである。この液面の変動により、二光束干渉計における二光束の間の光路長差を安定化することが困難であった。
【０００５】
従来では、このような問題に対処する為、培養液に接するように蓋またはカバーガラスを被せたり或いは水没型水浸対物レンズを用いたりすることで、培養液と気体（一般には空気）との界面を光が通過することがないよう工夫をしていた（非特許文献１参照）。この工夫により、光路上においては液面の変動が生じないようにすることができ、二光束干渉計における二光束の間の光路長差を安定化することができる。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】国際公開第２０１６／１２１２４８号
【非特許文献】
【０００７】
【非特許文献１】山内豊彦、他、「定量位相顕微鏡の小型化と各種応用計測 ＝空気定盤不要な卓上二光束干渉顕微鏡＝」、光アライアンス、２０１７年２月、ｐｐ.２６〜２９.
【非特許文献２】Daniel Malacara, et al, "Interferogram Analysis for Optical Testing," SecondEdition, Taylor & Francis, (2005).
【非特許文献３】武田光夫、他、「高速フーリエ変換法による立体計測と干渉計測」、光学、第10巻第6号、pp.476-480、(1981).
【発明の概要】
【発明が解決しようとする課題】
【０００８】
従来では、液面の変動を抑制するために設けた蓋、カバーガラスまたは水没型水浸対物レンズが培養液に接することから、これらに因り細胞が汚染されるという問題が生じる場合があった。
【０００９】
本発明は、上記問題点を解消する為になされたものであり、二光束干渉計において二光束の間の光路長差が時間的に変動している状態であっても対象物の位相画像を取得することができる装置および方法を提供することを目的とする。
【課題を解決するための手段】
【００１０】
本発明の位相画像取得装置は、(1) 光源と、光源から出力された光を二光束に分岐した後に該二光束を合波して干渉光を出力する二光束干渉計と、二光束干渉計から出力された干渉光を受光する撮像器とを含み、二光束干渉計における何れかの光束の光路上に配置された対象物の干渉画像を撮像器により撮像して取得する干渉画像取得部と、(2) 干渉画像取得部により取得された干渉画像に基づいて対象物の位相画像を生成する演算部と、を備える。干渉画像取得部は、二光束干渉計において二光束の間の光路長差が時間的に変動している状態で対象物についてＭ個の干渉画像を取得する。演算部は、Ｍ個の干渉画像それぞれについて画像全体のオフセット位相シフト量を推定し、Ｍ個の干渉画像のうちからオフセット位相シフト量が互いに異なるＮ個の干渉画像を含む組を設定し、組に含まれるＮ個の干渉画像に基づいて位相画像を生成する。ただし、３≦Ｎ≦Ｍ。
【００１１】
本発明の一側面による位相画像取得装置では、演算部は、組を複数設定して、これら複数の組それぞれについて各組に含まれるＮ個の干渉画像に基づいて位相画像を生成し、複数の組それぞれで生成された位相画像についてオフセット位相シフト量を補正した後に積算する。演算部は、Ｍ個の干渉画像をオフセット位相シフト量について昇順または降順に並び替えた後に複数の組に順に振り分けてもよい。或いは、演算部は、位相画像を生成するために干渉画像に対して乗算される係数を反復的な最適化計算によって調整してもよい。
【００１２】
また、本発明の位相画像取得装置は、上記と同様の干渉画像取得部と、次のような演算部と、を備えてもよい。すなわち、演算部は、Ｍ個の干渉画像それぞれについてオフセット位相シフト量の分布を推定し、各画素位置について、Ｍ個の干渉画像それぞれの当該画素のうちからオフセット位相シフト量が互いに異なるＮ個の画素を含む組を設定し、各画素位置について、組に含まれるＮ個の画素に基づいて位相値を求め、各画素位置の位相値に基づいて位相画像を生成してもよい。ただし、３≦Ｎ≦Ｍ。このとき、一側面において、演算部は、組を複数設定して、これら複数の組それぞれについて各組に含まれるＮ個の画素に基づいて位相値を求め、複数の組それぞれで求めた位相値についてオフセット位相シフト量を補正した後に積算し、各画素位置の位相値の積算値に基づいて位相画像を生成する。また、演算部は、各画素位置について、Ｍ個の干渉画像それぞれの当該画素をオフセット位相シフト量について昇順または降順に並び替えた後に複数の組に順に振り分けてもよい。
【００１３】
本発明の他の一側面による位相画像取得装置では、二光束干渉計は、二光束の双方または何れか一方の光路上に配置された液体の液面の揺らぎにより、二光束の間の光路長差の時間的な変動を生じさせる。
【００１４】
本発明の更に他の一側面による位相画像取得装置では、光源はインコヒーレント光を出力する。
【００１５】
本発明の位相画像取得方法は、(1) 光源と、光源から出力された光を二光束に分岐した後に該二光束を合波して干渉光を出力する二光束干渉計と、二光束干渉計から出力された干渉光を受光する撮像器とを含む干渉画像取得部を用い、二光束干渉計において二光束の間の光路長差が時間的に変動している状態で、二光束干渉計における何れかの光束の光路上に配置された対象物についてＭ個の干渉画像を撮像器により撮像して取得する干渉画像取得ステップと、(2) Ｍ個の干渉画像それぞれについて画像全体のオフセット位相シフト量を推定するオフセット位相シフト量推定ステップと、(3) Ｍ個の干渉画像のうちからオフセット位相シフト量が互いに異なるＮ個の干渉画像を含む組を設定する組設定ステップと、(4) 組に含まれるＮ個の干渉画像に基づいて位相画像を生成する位相画像生成ステップと、を備える。ただし、３≦Ｎ≦Ｍ。
【００１６】
本発明の一側面による位相画像取得方法では、組設定ステップにおいて組を複数設定し、位相画像生成ステップにおいて複数の組それぞれについて各組に含まれるＮ個の干渉画像に基づいて位相画像を生成し、複数の組それぞれで生成された位相画像についてオフセット位相シフト量を補正した後に積算する積算ステップを更に備える。組設定ステップにおいて、Ｍ個の干渉画像をオフセット位相シフト量について昇順または降順に並び替えた後に複数の組に順に振り分けてもよい。或いは、位相画像生成ステップにおいて、位相画像を生成するために干渉画像に対して乗算される係数を反復的な最適化計算によって調整してもよい。
【００１７】
また、本発明の位相画像取得方法は、(1) 上記と同様の干渉画像取得ステップと、(2) Ｍ個の干渉画像それぞれについてオフセット位相シフト量の分布を推定するオフセット位相シフト量推定ステップと、(3) 各画素位置について、Ｍ個の干渉画像それぞれの当該画素のうちからオフセット位相シフト量が互いに異なるＮ個の画素を含む組を設定する組設定ステップと、(4) 各画素位置について、組に含まれるＮ個の画素に基づいて位相値を求める位相値算出ステップと、(5) 各画素位置の位相値に基づいて位相画像を生成する位相画像生成ステップと、を備えてもよい。ただし、３≦Ｎ≦Ｍ。このとき、一側面において、本発明の位相画像取得方法では、組設定ステップにおいて組を複数設定し、位相値算出ステップにおいて複数の組それぞれについて各組に含まれるＮ個の画素に基づいて位相値を求め、複数の組それぞれで求めた位相値についてオフセット位相シフト量を補正した後に積算する積算ステップを更に備え、位相画像生成ステップにおいて各画素位置の位相値の積算値に基づいて位相画像を生成する。組設定ステップにおいて、各画素位置について、Ｍ個の干渉画像それぞれの当該画素をオフセット位相シフト量について昇順または降順に並び替えた後に複数の組に順に振り分けてもよい。
【００１８】
本発明の他の一側面による位相画像取得方法では、干渉画像取得ステップにおいて、二光束干渉計における二光束の双方または何れか一方の光路上に配置された液体の液面の揺らぎにより、二光束の間の光路長差の時間的な変動を生じさせる。
【００１９】
本発明の更に他の一側面による位相画像取得方法では、光源はインコヒーレント光を出力する。
【発明の効果】
【００２０】
本発明によれば、二光束干渉計において二光束の間の光路長差が時間的に変動している状態であっても対象物の位相画像を取得することができる。
【図面の簡単な説明】
【００２１】
【図１】図１は、位相画像取得装置１の構成例を示す図である。
【図２】図２は、サンプルの構成例を示す図である。
【図３】図３は、実施例で取得された１００個の干渉画像のうちの３個の干渉画像を示す図である。
【図４】図４は、実施例で取得された１００個の干渉画像それぞれのオフセット位相シフト量を示すグラフである。
【図５】図５は、実施例で取得された位相アンラピング処理前の位相画像である。
【図６】図６は、実施例で取得された位相アンラピング処理および背景歪み補正の後の位相画像である。
【図７】図７は、Ｍ個の干渉画像それぞれのオフセット位相シフト量αを単位円上にフェーザ表示した例を示す図である。
【図８】図８は、Ｍ個の干渉画像それぞれのオフセット位相シフト量αを単位円上にフェーザ表示した図において、Ｍ個の干渉画像のうちから３個の干渉画像を選択する方法を説明する図である。
【図９】図９は、α２＝−１２０°としα３＝＋１２０°としてノイズを付与しなかった場合に得られる位相画像（正解画像）である。
【図１０】図１０は、α２＝−１２０°としα３＝＋１２０°としてノイズを付与した場合の干渉画像Ｉ１〜Ｉ３である。
【図１１】図１１は、図１０の干渉画像Ｉ１〜Ｉ３に基づいて得られた位相画像である。
【図１２】図１２は、α２＝−２０°としα３＝＋３０°としてノイズを付与した場合の干渉画像Ｉ１〜Ｉ３である。
【図１３】図１３は、図１２の干渉画像Ｉ１〜Ｉ３に基づいて得られた位相画像である。
【図１４】図１４は、図９の正解画像に対する図１３の位相画像の誤差を表す画像である。
【図１５】図１５は、３個の干渉画像それぞれのオフセット位相シフト量αを単位円上にフェーザ表示した図である。
【図１６】図１６は、オフセット位相シフト量α２，α３と標準偏差stdとの間の関係を濃淡で示す図である。
【図１７】図１７は、オフセット位相シフト量α２，α３と標準偏差stdとの間の関係を等高線で示す図である。
【図１８】図１８は、オフセット位相シフト量α２，α３を様々な値に設定して得られた三角形の面積areaおよび標準偏差stdをプロットしたグラフである。
【図１９】図１９は、三角形の面積areaに対する標準偏差stdの最小値の関係を示すグラフである。
【図２０】図２０は、三角形の面積areaに対する標準偏差stdの最小値の関係を示す表である。
【図２１】図２１は、シミュレーションで用意した１２個の干渉画像Ｉ１〜Ｉ12のオフセット位相シフト量、組分け及び複素係数列を纏めた表である。
【図２２】図２２は、第１の方法で得られた位相画像である。
【図２３】図２３は、図９の正解画像に対する図２２の位相画像の誤差を表す画像である。
【図２４】図２４は、第２の方法で得られた位相画像である。
【図２５】図２５は、図９の正解画像に対する図２４の位相画像の誤差を表す画像である。
【図２６】図２６は、第３の方法で得られた位相画像である。
【図２７】図２７は、図９の正解画像に対する図２６の位相画像の誤差を表す画像である。
【図２８】図２８は、或る期間において各干渉画像のオフセット位相シフト量が時間的に概ね単調に増加している場合を示す図である。
【図２９】図２９は、第４の方法で得られた位相画像である。
【図３０】図３０は、図９の正解画像に対する図２９の位相画像の誤差を表す画像である。
【図３１】図３１は、第２〜第４の各方法によるシミュレーションで用意した１２個の干渉画像Ｉ１〜Ｉ12のオフセット位相シフト量、組分け及び複素係数列を纏めた表である。
【図３２】図３２は、第１〜第４の各方法によるシミュレーションで得られた誤差の標準偏差を比較して示すグラフである。
【図３３】図３３は、第２の方法による複素係数列および第４の方法による複素係数列を複素平面上にプロットした図である。
【図３４】図３４は１２個の干渉画像それぞれのオフセット位相シフト量αを単位円上にフェーザ表示した図において、第２の方法に基づいて１２個の干渉画像のうちから４組の３個の干渉画像を選択した方法を説明する図である。
【図３５】図３５は１２個の干渉画像それぞれのオフセット位相シフト量αを単位円上にフェーザ表示した図において、第３の方法に基づいて１２個の干渉画像のうちから４組の３個の干渉画像を選択した方法を説明する図である。
【図３６】図３６は、シミュレーションにおいて用いた位相画像を示す図である。
【図３７】図３７は、シミュレーションにおいて用いた位相画像を示す図である。
【図３８】図３８（ａ）は、位置(x,y)に依らず位相値が１である場合の位相分布関数を示す図である。図３８（ｂ）は、ｘ座標値に比例する位相値を有する場合の位相分関数布を示す図である。図３８（ｃ）は、ｙ座標値に比例する位相値を有する場合の位相分布関数を示す図である。図３８（ｄ）は、ｘ２＋ｙ２−１の値に比例する位相値を有する場合の位相分布関数を示す図である。
【図３９】図３９（ａ）〜（ｄ）は、オフセット位相シフト量分布Δφｎ(x,y)の例を示す図である。
【図４０】図４０は、図３９（ａ）〜（ｄ）のオフセット位相シフト量分布Δφｎ(x,y)の各例において設定した係数p1(n)，p2(n)，p3(n)，p4(n)の各値を纏めた表である。
【図４１】図４１は、シミュレーションで用いた３６個の干渉画像のうちの３個の干渉画像を示す図である。
【図４２】図４２は、或る干渉画像Ｉｎ(x,y)を二次元フーリエ変換することで得られたフーリエ変換像を示す図である。
【図４３】図４３は、図４２のフーリエ変換像の原点付近を拡大して示す図である。
【図４４】図４４は、図４２のフーリエ変換像における一方のスペクトルピーク位置を含む一定領域の空間周波数情報のみを用いて二次元逆フーリエ変換することで得られた二次元の位相分布画像を示す図である。
【図４５】図４５は、従来のフーリエ変換法により得られた位相画像の中央領域を示す図である。
【図４６】図４６は、基準となる位相分布画像φLPF0(x,y)の一例を示す図である。
【図４７】図４７は、基準となる位相分布画像φLPF0(x,y)に対する相対的な位相分布画像ΔφLPFｎ(x,y)の一例を示す図である。
【図４８】図４８（ａ）は、図４７の相対的な位相分布画像ΔφLPFｎ(x,y)の中央領域における相対的な位相分布を示す図である。図４８（ｂ）は、図４８（ａ）の相対的な位相分布を近似して得られた相対的な位相分布を示す図である。
【図４９】図４９は、シミュレーションにおいて生成された位相画像を示す図である。
【図５０】図５０は、シミュレーションにおいて生成された位相画像を示す図である。
【図５１】図５１は、シミュレーションにおいて生成された位相画像と正解画像（図３６）との差分の画像を示す図である。
【発明を実施するための形態】
【００２２】
以下、添付図面を参照して、本発明を実施するための形態を詳細に説明する。なお、図面の説明において同一の要素には同一の符号を付し、重複する説明を省略する。本発明は、これらの例示に限定されるものではなく、特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれることが意図される。
【００２３】
図１は、位相画像取得装置１の構成例を示す図である。位相画像取得装置の光学系は様々な構成を有するものであってもよく、この図に示した光学系の構成は一例にすぎない。この図に示される位相画像取得装置１は、干渉画像取得部２および演算部３を備える。干渉画像取得部２は、光源１０、レンズ１１、レンズ１２、分岐用ビームスプリッタ２１、合波用ビームスプリッタ２２、ビームスプリッタ３１、第１反射ミラー３２、ビームスプリッタ４１、第２反射ミラー４２、ピエゾ素子４３、ステージ４４、レンズ５１、レンズ５２、チューブレンズ５３、撮像器６１および制御部７２を含む。分岐用ビームスプリッタ２１から合波用ビームスプリッタ２２に到る迄の二光束干渉計２０は、マッハツェンダ干渉計を構成している。
【００２４】
この位相画像取得装置１は、容器９０に入れられた対象物の透過光に基づいて干渉画像を取得する。対象物は、特定の細胞や生体サンプルに限定されるものではない。例えば、対象物として、培養細胞、不死化細胞、初代培養細胞、がん細胞、脂肪細胞、肝臓細胞、心筋細胞、神経細胞、グリア細胞、体性幹細胞、胚性幹細胞、多能性幹細胞、ｉＰＳ細胞、およびこれらの細胞のうち少なくとも１つの細胞をもとに作られた細胞塊（スフェロイド）などが挙げられる。また、対象物として、生体に限らず、透過型の構成で計測可能な工業サンプル、たとえばガラス内部、半導体素子の内部、樹脂素材、液晶、高分子化合物、光学素子なども挙げられる。
【００２５】
光源１０は光を出力する。スペックルノイズを低減するには、光源１０はインコヒーレントな光を出力するのが好適である。光源１０は、例えばハロゲンランプなどのランプ系光源、ＬＥＤ（Light emitting diode）光源、ＳＬＤ（Superluminescent diode）光源、ＡＳＥ（Amplified spontaneous emission）光源等である。レンズ１１，１２は、光源１０から出力された光を容器９０内の対象物に集光する。
【００２６】
分岐用ビームスプリッタ２１は、光源１０と光学的に結合され、光源１０から出力されてレンズ１１，１２を経た光を入力し、その光を２分岐して第１光束および第２光束とする。分岐用ビームスプリッタ２１は例えばハーフミラーであってもよい。分岐用ビームスプリッタ２１は、第１光束を測定側光学系のビームスプリッタ３１へ出力し、第２光束を参照側光学系のビームスプリッタ４１へ出力する。
【００２７】
測定側光学系には、ビームスプリッタ３１および第１反射ミラー３２が設けられている。ビームスプリッタ３１は、分岐用ビームスプリッタ２１から出力された第１光束を入力して第１反射ミラー３２へ反射させ、また、第１反射ミラー３２で反射された第１光束を入力して合波用ビームスプリッタ２２へ透過させる。ビームスプリッタ３１は例えばハーフミラーであってもよい。
【００２８】
参照側光学系には、ビームスプリッタ４１、第２反射ミラー４２、ピエゾ素子４３およびステージ４４が設けられている。ビームスプリッタ４１は、分岐用ビームスプリッタ２１から出力された第２光束を入力して第２反射ミラー４２へ反射させ、また、第２反射ミラー４２で反射された第２光束を入力して合波用ビームスプリッタ２２へ透過させる。ビームスプリッタ４１は例えばハーフミラーであってもよい。
【００２９】
ピエゾ素子４３は、第２反射ミラー４２の反射面に垂直な方向に該第２反射ミラー４２を移動させることができる。ステージ４４は、第２反射ミラー４２の反射面に垂直な方向に該第２反射ミラー４２およびピエゾ素子４３を移動させることができる。ピエゾ素子４３およびステージ４４は、参照側光学系の光路長を調整することができ、第１光束と第２光束との間の光路長差を調整することができる。ステージ４４は光路長差を粗調整することができ、ピエゾ素子４３は光路長差を微調整することができる。
【００３０】
合波用ビームスプリッタ２２は、ビームスプリッタ３１から出力され容器９０内の対象物およびレンズ５１を経た第１光束を入力するとともに、ビームスプリッタ４１から出力されレンズ５２を経た第２光束を入力して、これら第１光束と第２光束とを合波して合波光（干渉光）を出力する。合波用ビームスプリッタ２２は例えばハーフミラーであってもよい。
【００３１】
チューブレンズ５３は、合波用ビームスプリッタ２２から出力された干渉光を撮像器６１へ導き、その干渉光を撮像器６１の撮像面上に結像させる。撮像器６１は、その干渉光を受光して検出信号を出力するものであり、特に撮像面上の干渉光の強度分布を表す検出信号を出力する。撮像器６１は、例えばＣＣＤエリアイメージセンサやＣＭＯＳエリアイメージセンサなどのイメージセンサである。
【００３２】
演算部３は、撮像器６１から出力された検出信号を入力し、その検出信号に基づいて容器９０内の対象物の干渉画像を取得し、更に位相画像を取得する。演算部３は、ＦＰＧＡ（field-programmable gate array）やＧＰＵ（GraphicsProcessing Unit）などの画像処理プロセッサを含んで構成され、或いは、パーソナルコンピュータやタブレット端末などのコンピュータであってもよい。また、演算部３は、干渉画像等を表示する表示部を備えていてもよい。
【００３３】
制御部（コントローラ）７２は、ピエゾ素子４３およびステージ４４の双方または何れか一方を駆動することで第２反射ミラー４２を移動させて、参照側光学系の光路長を調整する。これにより、制御部７２は、合波用ビームスプリッタ２２による合波の際の第１光束と第２光束との間の位相差を調整することができる。
【００３４】
なお、演算部３および制御部７２は、プロセッサ及びメモリ等を含むコンピュータである。また、演算部３および制御部７２は、個別のコンピュータであってもよいし、１つのコンピュータであってもよい。コンピュータは、例えば、パーソナルコンピュータやタブレット端末などのスマートデバイスであってもよい。また、演算部３または制御部７２は、利用者からの入力を受け付ける入力部（キーボード、マウス、タブレット端末等）や、干渉強度等を表示する表示部（ディスプレイ、タブレット端末、スピーカー、バイブレータ）を備えていてもよい。なお、表示部がディスプレイやタブレット端末などのように、画面表示できる装置であれば、干渉強度と合わせて干渉画像などを表示してもよい。
【００３５】
位相画像取得装置１を用いて以下のようにして容器９０内の対象物の干渉画像を取得することができる。光源１０から出力された光は、レンズ１１，１２を経て、分岐用ビームスプリッタ２１により２分岐されて第１光束および第２光束とされる。分岐用ビームスプリッタ２１から出力された第１光束は、ビームスプリッタ３１を透過して、第１反射ミラー３２で反射される。第１反射ミラー３２で反射された第１光束は、ビームスプリッタ３１で反射されて、容器９０内の対象物に集光され、対象物を透過する。対象物を透過した第１光束は、レンズ５１を経て合波用ビームスプリッタ２２に入力される。この第１光束は、対象物における透過の際に光学的な遅れを有する。分岐用ビームスプリッタ２１から出力された第２光束は、ビームスプリッタ４１で反射されて、第２反射ミラー４２で反射される。第２反射ミラー４２で反射された第２光束は、ビームスプリッタ４１を透過して、レンズ５２を経て合波用ビームスプリッタ２２に入力される。
【００３６】
レンズ５１から合波用ビームスプリッタ２２に入力された第１光束、および、レンズ５２から合波用ビームスプリッタ２２に入力された第２光束は、合波用ビームスプリッタ２２により合波される。その合波による干渉光は、チューブレンズ５３を経て、撮像器６１により受光される。演算部３により、干渉光を受光した撮像器６１から出力される検出信号に基づいて干渉画像が取得され、更に、複数の干渉画像に基づいて位相画像が取得される。また、制御部７２により駆動されたピエゾ素子４３またはステージ４４により第２反射ミラー４２の位置が制御されることにより、第１光束と第２光束との間の光路長差が調整され、合波用ビームスプリッタ２２による合波の際の第１光束と第２光束との間の位相差が調整される。
【００３７】
従来技術では、干渉画像を取得する際に、二光束干渉計２０における二光束の間の光路長差を互いに異なる複数の設定値それぞれで安定化した状態とし、各状態において対象物の干渉画像を撮像器６１により取得し、これら取得した複数の干渉画像に基づいて位相画像を求める。また、二光束干渉計２０における二光束の間の光路長差を安定化する為に、この光路長差を検出するとともに、この光路長差が設定値となるように、二光束のうちの何れかの光路上に配置したミラーの位置をフィードバック制御する。しかし、このようなフィードバック制御をしても、例えば対象物としての細胞が培養液に入れられていて該培養液の液面が揺らいでいる場合には、その液面の変動による光路長の変動量が、フィードバック制御により安定化しうる光路長の変動量よりも大きいことにより、二光束干渉計２０における二光束の間の光路長差を安定化することが困難であった。この問題に対処するため、従来技術では、培養液の液面の変動を抑制するために、培養液に接するように蓋またはカバーガラスをかぶせたり、水没型水浸対物レンズを用いたりすることで、培養液と気体（例えば、空気）との界面を光が通過しないよう工夫していたが、その場合には細胞が汚染されるという問題が生じる場合があった。
【００３８】
これに対して、以下に説明する実施形態では、例えば対象物としての細胞が培養液に入れられていて該培養液の液面が揺らいでいる場合等のように、二光束干渉計において二光束の間の光路長差が時間的に変動している状態であっても、対象物の位相画像を取得することができる。なお、光路長差の時間的な変動は、培養液とは別に第１光束または第２光束の光路上に配置された液体の液面の揺らぎに因るものであってもよいし、ピエゾ素子４３により駆動された第２反射ミラー４２の振動に因るものであってもよいし、また、外部から装置に加えられる振動に因るものであってもよい。光路長差の時間的な変動は、非周期的でランダムなものであってもよいし、周期的なものであってもよいし、また、一方向に漸増または漸減するものであってもよい。十分な大きさの位相差を得る為に外部から装置に振動を加えてもよいし、十分な大きさの位相差が得られるまで干渉画像の取得を続けてもよい。
【００３９】
図２は、サンプルの構成例を示す図である。対象物は、容器９０に入れられた培養液９２中の細胞９３である。容器９０の上部は蓋９１が被せられている。培養液９２と蓋９１との間に空間があり、培養液９２の液面は気体（一般には空気）と接している。したがって、培養液９２の液面の揺らぎに因り、培養液９２を透過する光の光路長は時間的に変動し、二光束干渉計２０において二光束の間の光路長差は時間的に変動する。
【００４０】
本実施形態では、干渉画像取得部２は、二光束干渉計２０において二光束の間の光路長差が時間的に変動している状態で対象物についてＭ個の干渉画像を取得する（干渉画像取得ステップ）。Ｍ個の干渉画像それぞれの取得タイミングは互いに異なる。そして、演算部３は、これらＭ個の干渉画像それぞれについて画像全体のオフセット位相シフト量を推定し（オフセット位相シフト量推定ステップ）、Ｍ個の干渉画像のうちからオフセット位相シフト量が互いに異なるＮ個の干渉画像を含む組を設定し（組設定ステップ）、その組に含まれるＮ個の干渉画像に基づいて位相画像を生成する（位相画像生成ステップ）。
【００４１】
好適には、演算部３は、上記のような組を複数設定して（組設定ステップ）、これら複数の組それぞれについて各組に含まれるＮ個の干渉画像に基づいて位相画像を生成し（位相画像生成ステップ）、これら複数の組それぞれで生成された位相画像についてオフセット位相シフト量を補正した後に積算する（積算ステップ）。なお、積算ステップにおいて、積算値を組の数で除算することで平均値を求めてもよく、積算値と平均値との間には実質的な差異はない。Ｍ個の干渉画像のうち、何れの組にも含まれない干渉画像があってもよい。また、Ｍ個の干渉画像のうち、２以上の組に含まれる干渉画像があってもよい。
【００４２】
Ｎは３以上の整数である。Ｎ個の干渉画像を含む組を１つのみ設定するのであれば、ＭはＮ以上の整数である。Ｎ個の干渉画像を含む組を複数設定するのであれば、Ｍは(Ｎ＋１)以上の整数である。なお、ＳＮが良好な位相画像を求めるには、組の数は多いほど好ましく、Ｍは大きいほど好ましい。
【００４３】
各干渉画像について画像全体のオフセット位相シフト量の推定（オフセット位相シフト量推定ステップ）は例えば次のようにして行われる。干渉画像において干渉縞が直線状であって等間隔に現れている領域に注目し、この領域の干渉強度分布を空間的に複素フーリエ変換する。この複素フーリエ変換で得られた空間周波数分布においてピークとなる複素空間周波数の位相値を求め、この位相値を干渉画像のオフセット位相シフト量として推定する。複素フーリエ変換の対象となる領域は、二次元的に干渉画像を切り出したものであってもよいし、一次元的に干渉画像の一部のラインプロファイルを抽出したものであってもよい。
【００４４】
各組に含まれるＮ個の干渉画像に基づく位相画像の生成（位相画像生成ステップ）は次のようにして行われる。以下では、Ｎ＝３の場合について説明する。或る組に含まれる３個の干渉画像Ｉ１〜Ｉ３は下記（１）式で表される。Ｉ１〜Ｉ３、ａ、ｂ、φは、位置（ｘ，ｙ）の関数である。φは、位相画像の画素値として求めるべき位相値である。α１〜α３は、各干渉画像の取得タイミングでの二光束干渉計における光路長差に起因するオフセット位相シフト量である。この（１）式から下記（２）式が得られる（非特許文献２参照）。この（２）式の右辺は、オフセット位相シフト量α１〜α３および位相値φで表され、ａ、ｂに依存しない。各干渉画像のオフセット位相シフト量α１〜α３が既に推定されているので、この（２）式から位相値φを求めることができる。すなわち、位相画像を求めることができる。
【００４５】
【数１】
2019124679.tif 000003
【００４６】
【数２】
2019124679.tif 000004
【００４７】
複数の組それぞれで生成された位相画像の積算（積算ステップ）において、積算により画素値が減衰することがないように、各位相画像のオフセット位相シフト量を補正した後に、各位相画像を複素振幅画像として積算する。このように積算することで、位相画像の各画素値の誤差を小さくすることができる。
【００４８】
次に、実施例について説明する。容器として、３５ｍｍφのガラスボトムディッシュを用いた。容器の底部内側に対象物としてイマージョンオイルを付着させ、これに深さ７ｍｍ程度に水を入れて、蓋を被せた。光源として中心波長６５０ｎｍのレーザダイオードを用いた。このレーザダイオードへの注入電流を周波数１００ｋＨｚで変調することで、レーザダイオードから出力される光の線幅を拡げた。このように線幅を拡げることで出力光をインコヒーレントとして、スペックルノイズおよびレーザ光の不安定性を低減した。各干渉画像を取得する際の露光時間は３ミリ秒であり、撮像フレームレートは３０ｆｐｓであった。１００個の干渉画像を取得した。
【００４９】
図３は、実施例で取得された１００個の干渉画像のうちの３個の干渉画像を示す図である。二光束干渉計において二光束の間の光路長差が時間的に変動していることにより、オフセット位相シフト量も時間的に変動していることから、これら３個の干渉画像は互いに異なっている。
【００５０】
図４は、実施例で取得された１００個の干渉画像それぞれのオフセット位相シフト量を示すグラフである。なお、この図では時間方向の位相アンラッピング処理を行っていない。横軸は干渉画像の番号である。縦軸は、オフセット位相シフト量推定ステップにおいて各干渉画像について推定されたオフセット位相シフト量（ラジアン）である。このように、オフセット位相シフト量が時間的に変動していることが分かる。
【００５１】
組設定ステップにおいて、各組に含まれる干渉画像の個数Ｎを３として、２０組を設定した。２０組のうちの第ｋの組は、１００個の干渉画像のうちの第(５ｋ−４)の干渉画像と、この第(５ｋ−４)の干渉画像に対してオフセット位相シフト量の差が＋１２０°に最も近い干渉画像と、この第(５ｋ−４)の干渉画像に対してオフセット位相シフト量の差が−１２０°に最も近い干渉画像と、を含むように設定した。ｋは１〜２０の各整数である。
【００５２】
位相画像生成ステップにおいて、２０組それぞれについて各組に含まれる３個の干渉画像に基づいて位相画像を生成した。そして、積算ステップにおいて、２０組それぞれで生成された位相画像についてオフセット位相シフト量を補正した後に積算して、最終的な位相画像を得た。
【００５３】
図５および図６は、実施例で取得された位相画像を示す図である。図５は、位相アンラピング処理前の位相画像である。図６は、位相アンラピング処理および背景歪み補正の後の位相画像である。この位相画像は、容器の底部内側に付着させた対象物としてのイマージョンオイルの形状をよく表している。なお、対象物の厚さの最大値は４０ラジアンの位相差に相当するが、このような大きな位相値であっても、位相アンラッピング処理に問題は生じない。
【００５４】
このように、例えば対象物としての細胞が培養液に入れられていて該培養液の液面が揺らいでいる場合等のように、二光束干渉計において二光束の間の光路長差が時間的に変動している状態であっても、対象物の位相画像を取得することができる。従来のように培養液の液面の変動を抑制するために培養液に接するようにカバーガラス等を設ける必要がないので、カバーガラス等による細胞の汚染の問題が生じない。また、機械的な可動部が無くても位相画像を取得することができるので、機械的部品の劣化の問題を抑制することができ、簡易な構成とすることができる。
【００５５】
次に、Ｎ個の干渉画像を含む組の設定（組設定ステップ）について詳細に説明する。Ｎ個の干渉画像を含む組を設定するに際して、上記（２）式に基づく位相値φの算出が精度よく行われるように、Ｎ個の干渉画像を選択するのが好ましい。その為には、次のようにしてＮ個の干渉画像を選択することができる。まず、Ｍ個の干渉画像それぞれのオフセット位相シフト量αについて、横軸をcosαとし縦軸をsinαとする直交座標系において原点を中心とする半径１の単位円の周上にフェーザ表示する。図７は、Ｍ個の干渉画像それぞれのオフセット位相シフト量αを単位円上にフェーザ表示した例を示す図である。
【００５６】
この図において、単位円上の黒丸が、フェーザ表示したオフセット位相シフト量αの位置を表している。単位円上にフェーザ表示されたＭ個のオフセット位相シフト量αのうちから、Ｎ個のオフセット位相シフト量αの点を順に線分で結ぶことで形成されるＮ角形の面積が出来る限り大きくなるように（最も好適には正Ｎ角形に近い形状となるように）、Ｎ個のオフセット位相シフト量αを選択し、これらに対応するＮ個の干渉画像を選択して組とするのが好ましい。
【００５７】
図８は、Ｍ個の干渉画像それぞれのオフセット位相シフト量αを単位円上にフェーザ表示した図において、Ｍ個の干渉画像のうちから３個の干渉画像を選択する方法を説明する図である。Ｎ＝３の場合、単位円上にフェーザ表示されたＭ個のオフセット位相シフト量αのうちから、３個のオフセット位相シフト量αの点を順に線分で結ぶことで形成される三角形の面積が出来る限り大きくなるように（最も好適には正三角形に近い形状となるように）、３個のオフセット位相シフト量αを選択し、これらに対応する３個の干渉画像を選択して組とする。この図では、各組に含まれる３個の干渉画像のオフセット位相シフト量αの点が、実線、破線または一点鎖線で結ばれている。３個の干渉画像を含む組が３つ設定されている。
【００５８】
Ｎ＝３の場合に、単位円上にフェーザ表示された３個のオフセット位相シフト量αの点を順に線分で結ぶことで形成される三角形の面積と、これらに対応する３個の干渉画像に基づいて生成される位相画像における各画素値の誤差との間の関係について、シミュレーションを行った。３個の干渉画像のうちの第ｎの干渉画像Ｉｎを次式で表した。ａ＝１.５とし、ｂ＝１.０とした。干渉画像Ｉｎの各画素値にノイズを付与した。付与したノイズは、ｎｆ＝０.１５sqrt(Ｉｎ) として、−ｎｆ／２〜＋ｎｆ／２ の範囲の一様乱数とした。すなわち、光子数の平方根に比例するショットノイズを模擬した。そして、３個の干渉画像のうちの第１の干渉画像Ｉ１のオフセット位相シフト量α１を０とし、第２の干渉画像Ｉ２のオフセット位相シフト量α２および第３の干渉画像Ｉ３のオフセット位相シフト量α３を様々な値に設定した。
【００５９】
【数３】
2019124679.tif 000005
【００６０】
図９は、α２＝−１２０°としα３＝＋１２０°としてノイズを付与しなかった場合に得られる位相画像（正解画像）である。図１０は、α２＝−１２０°としα３＝＋１２０°としてノイズを付与した場合の干渉画像Ｉ１〜Ｉ３である。図１１は、図１０の干渉画像Ｉ１〜Ｉ３に基づいて得られた位相画像である。図１２は、α２＝−２０°としα３＝＋３０°としてノイズを付与した場合の干渉画像Ｉ１〜Ｉ３である。図１３は、図１２の干渉画像Ｉ１〜Ｉ３に基づいて得られた位相画像である。図１４は、図９の正解画像に対する図１３の位相画像の誤差を表す画像である。これらの図から、干渉画像にノイズが含まれている場合に位相画像の誤差を小さくするには、３個のオフセット位相シフト量αの点を順に線分で結ぶことで形成される三角形の面積が大きいのが好ましいことが分かる。
【００６１】
図１５は、３個の干渉画像それぞれのオフセット位相シフト量αを単位円上にフェーザ表示した図である。３個の干渉画像のうちの第１の干渉画像Ｉ１のオフセット位相シフト量α１を０とし、第２の干渉画像Ｉ２のオフセット位相シフト量を＋α２とし、第３の干渉画像Ｉ３のオフセット位相シフト量を−α３としている。α２の各値およびα３の各値に対して、３個のオフセット位相シフト量の点を順に線分で結ぶことで形成される三角形の面積areaを求めることができ、また、これらに対応する３個の干渉画像に基づいて生成される位相画像における各画素値の誤差の標準偏差stdを求めることができる。
【００６２】
図１６は、オフセット位相シフト量α２，α３と標準偏差stdとの間の関係を濃淡で示す図である。図１７は、オフセット位相シフト量α２，α３と標準偏差stdとの間の関係を等高線で示す図である。図１８は、オフセット位相シフト量α２，α３を様々な値に設定して得られた三角形の面積areaおよび標準偏差stdをプロットしたグラフである。図１９は、三角形の面積areaに対する標準偏差stdの最小値の関係を示すグラフである。図２０は、三角形の面積areaに対する標準偏差stdの最小値の関係を示す表である。これらの図から分かるように、α２＝α３＝１２０°のとき、すなわち、３点を結ぶ三角形が正三角形であるときに、誤差が最も小さくなる。三角形の面積areaが大きいほど誤差が小さい。
【００６３】
次に、Ｎ個の干渉画像を含む組の設定（組設定ステップ）の例について説明する。先ず、Ｍ個の干渉画像Ｉｍ（ｍ＝１〜Ｍ）が取得されており、各干渉画像Ｉｍのオフセット位相シフト量αｍが既知であるとする。前述したとおり、Ｍ個の干渉画像Ｉ１〜ＩＭのうちから任意の３個の干渉画像を選んで位相画像を生成することできる。ここで、上記（２）式を変形することで、生成される位相画像に対応する偏角を持つ複素振幅画像Comp(x,y) を定式化することができる。式変形の過程を以下に記す。
【００６４】
（２）式の右辺に現れる係数について複素数Ａ，Ｂを下記（４）式で表す。（２）式の右辺の分母および分子にcosφを乗じることにより、下記（５）式が得られる。未知数はcosφおよびsinφである。
【００６５】
【数４】
2019124679.tif 000006
【００６６】
【数５】
2019124679.tif 000007
【００６７】
下記（６）式の行列Ｍを用いると、（５）式は下記（７）式で表される。行列Ｍが正則行列であれば、（７）式は下記（８）式で表される。
【００６８】
【数６】
2019124679.tif 000008
【００６９】
【数７】
2019124679.tif 000009
【００７０】
【数８】
2019124679.tif 000010
【００７１】
複素振幅画像Comp(x,y) は、下記（９）式で表すと、下記（１０）式の形になる。すなわち、複素振幅画像Compは、干渉画像と複素数列との内積で表されることとなる。
【００７２】
【数９】
2019124679.tif 000011
【００７３】
【数１０】
2019124679.tif 000012
【００７４】
複数の組それぞれにおいて３個の干渉画像から複素振幅画像を生成して、これら複数個の複素振幅画像を積算した上で偏角をとるという演算は、下記（１１）式で複素振幅画像を生成して、この複素振幅画像の偏角をとって位相画像を求めることに等しい。
【００７５】
【数１１】
2019124679.tif 000013
【００７６】
したがって、複数の組それぞれにおいて３個の干渉画像から位相画像を生成することは、各組で選ばれた３個の干渉画像をＩm1，Ｉm2，Ｉm3としたときに、複素係数coef_m1，coef_m2，coef_m3を順次決定して行くという操作に相当する。これらより、本実施形態において複数の干渉画像から位相画像を生成すること、上記（１１）式の複素係数ベクトル（coef_1  …  coef_N）を決定するという操作に一般化することができる。
【００７７】
以下では、この係数を決定する方法として４つの方法を挙げて、それぞれのシミュレーション結果を示す。第１の方法では、３個の干渉画像を完全にランダムに選ぶ。第２の方法では、３個の干渉画像を適切に選ぶ。第３の方法では、３個の干渉画像の選び方が不適切である。第４の方法では、１２個の干渉画像の全てを用いた場合の係数ベクトルの各要素について反復的な最適化を行う。
【００７８】
このシミュレーションでは１２個の干渉画像Ｉ１〜Ｉ12を用意した。各干渉画像Ｉｍは上記（３）式で表されるものとして、ａ＝１.５とし、ｂ＝１.０とした。各干渉画像Ｉｍのオフセット位相シフト量αｍを、−π〜＋πの範囲で一様分布する乱数列とした。干渉画像Ｉｍの各画素値にノイズを付与した。付与したノイズは、ｎｆ＝０.１５sqrt(Ｉｍ) として、−ｎｆ／２〜＋ｎｆ／２ の範囲の一様乱数とした。すなわち、光子数の平方根に比例するショットノイズを模擬した。また、ノイズがない場合の位相画像（正解画像）は図９と同じである。
【００７９】
第１の方法によるシミュレーションの内容および結果は以下のとおりである。それぞれにランダムなオフセット位相シフト量αｍが与えられた１２個の干渉画像Ｉ１〜Ｉ12のうちから任意に３個の干渉画像を選ぶことで位相画像を生成することができるが、最も単純な方法として、干渉画像Ｉ１〜Ｉ３を第１の組とし、干渉画像Ｉ４〜Ｉ６を第２の組とし、干渉画像Ｉ７〜Ｉ９を第３の組とし、干渉画像Ｉ10〜Ｉ12を第４の組とした。その結果、図２１に記したような複素係数列が得られた。図２１は、シミュレーションで用意した１２個の干渉画像Ｉ１〜Ｉ12のオフセット位相シフト量、組分け及び複素係数列を纏めた表である。これらの複素係数列に基づいて上記（１１）式によって複素振幅画像を生成し、さらに位相画像を生成した。図２２は、第１の方法で得られた位相画像である。図２３は、図９の正解画像に対する図２２の位相画像の誤差を表す画像である。なお、これらの画像における単位はラジアンである。誤差の標準偏差は３３.４ミリラジアンであった。
【００８０】
第２の方法によるシミュレーションの内容および結果は以下のとおりである。１２個の干渉画像Ｉ１〜Ｉ12をオフセット位相シフト量について昇順または降順に並び替え、その並び替えた後の１２個の干渉画像を順にＪ１〜Ｊ12とした。そして、１２個の干渉画像Ｊ１〜Ｊ12を４個の組に順に振り分けた。すなわち、干渉画像Ｊ１，Ｊ５，Ｊ９を第１の組とし、干渉画像Ｊ２，Ｊ５，Ｊ10を第２の組とし、干渉画像Ｊ３，Ｊ６，Ｊ11を第３の組とし、干渉画像Ｊ４，Ｊ７，Ｊ12を第４の組とした（図３４参照）。図２４は、第２の方法で得られた位相画像である。図２５は、図９の正解画像に対する図２４の位相画像の誤差を表す画像である。誤差の標準偏差は２３.６ミリラジアンであった。
【００８１】
第２の方法についてより一般化して説明する。Ｍ個の干渉画像Ｊ１〜ＪＭに対して、Ｍ／３を超えない最大の整数をＭＧとしたときに、これらの干渉画像をＭＧ個の組に振り分けることができる。振り分けに際しては、あらかじめ干渉画像Ｊ１〜ＪＭをオフセット位相シフト量について昇順または降順に並び替えておく。そのうえで、干渉画像Ｊ１，Ｊ1+ＭＧ，Ｊ1+２×ＭＧを第１の組とし、干渉画像Ｊ２，Ｊ２＋ＭＧ，Ｊ２＋２×ＭＧを第２の組とし、干渉画像Ｊ３，Ｊ３＋ＭＧ，Ｊ３＋２×ＭＧを第３の組とし、以下第ＭＧの組まで振り分けることができる。すなわち、第ｋ番目の組は、干渉画像Ｊｋ，Ｊｋ＋ＭＧ，Ｊｋ＋２×ＭＧによって構成される。このように振り分けることで、各組に含まれる３個の干渉画像の相対的な位相差が、概ね均等に１２０°ずつとなるように組分けすることができる。
【００８２】
第３の方法によるシミュレーションの内容および結果は以下のとおりである。１２個の干渉画像Ｉ１〜Ｉ12をオフセット位相シフト量について昇順または降順に並び替え、その並び替えた後の１２個の干渉画像を順にＪ１〜Ｊ12とした。そして、干渉画像Ｊ１〜Ｊ３を第１の組とし、干渉画像Ｊ４〜Ｊ６を第２の組とし、干渉画像Ｊ７〜Ｊ９を第３の組とし、干渉画像Ｊ10〜Ｊ12を第４の組とした（図３５参照）。図２６は、第３の方法で得られた位相画像である。図２７は、図９の正解画像に対する図２６の位相画像の誤差を表す画像である。誤差の標準偏差は２３５ミリラジアンであった。なお、このような干渉画像の組み分けは、図２８に示されるように或る期間において各干渉画像のオフセット位相シフト量が時間的に単調に増加（または減少）している場合に、連続して取得された３個の干渉画像を各組に含めることに相当する。
【００８３】
第４の方法によるシミュレーションの内容および結果は以下のとおりである。第４の方法では、１２個の干渉画像に対して組み分けを行わず、前記各組に含まれる干渉画像の個数ＮをＮ＝Ｍ＝１２として画像再構成を行っている。第１から第３の方法によって示されたように、Ｍ個の干渉画像から位相分布を再構成するための複素係数列としては、様々な値を取ることができ、元となる干渉画像にノイズが無い場合は、位相分布を与える複素係数列は無数に存在するが、この複素係数列の取り方を工夫することで、よりノイズに対してロバストな複素係数列を導くことができる。
【００８４】
まず、Ｍ個の干渉画像Ｉ１〜ＩＭの複素係数列（coef_1  …  coef_Ｍ）に関して、初期条件としての係数列を設定する。たとえば、第２の方法で求められた複素係数列を初期条件とすることができる。その後、複素係数列のうち一部の要素を微小量だけ複素平面上で変化させることを考える。例えば、Ｌ番目の係数coef_Lについて、その実部を0.1だけ増加させたり、その虚部を0.1だけ減少させたりする。この際に新たに作られた複素係数列（coef_1  …  coef_Ｍ）では位相画像の画像再構成ができないが、Ｌ番目ではない別の一つ以上の係数について、その実部および虚部の双方または何れか一方を微小量だけ増加または減少させることで、複素係数列を微小に補正して第２の係数列を作成することができる。この計算は、コンピュータ上でオフセット位相量がα１〜αＭであり−π〜＋πにわたって均等に位相が分布するような干渉画像をノイズの無い条件で作成し、それらに対して位相画像の画像再構成ができるように複素係数列を調整することで実現してもよい。
【００８５】
次に、新たに得られた複素係数列のノイズに対するロバスト性を評価する。ノイズに対するロバスト性の評価は、たとえばコンピュータ上でオフセット位相量がα１〜αＭでありノイズを含み、−π〜＋πにわたって均等に位相が分布するような干渉画像を作成し、それらに対して新たに得られた複素係数列を適用して位相画像を求め、理想的な正解となる位相画像との間で二乗平均誤差を求めることで実現してもよい。ここで、初期条件としての係数列を用いた位相画像における二乗平均誤差と、初期条件から微小に変化させて新たに得られた条件での係数列を用いた位相画像における二乗平均誤差とを比較し、新たに得られた複素係数列を用いた方がノイズに対して良好な再構成結果をもたらす場合は、新たに得られた複素係数列を暫定的な最善の複素係数列として採用する。
【００８６】
逆に、新たに得られた複素係数列を用いた方がノイズに対して劣悪な再構成結果をもたらす場合は、前段階の複素係数列を暫定的な最善の複素係数列として採用する。複素係数列の中で何番目の要素を微小変化させるか、また微小変化させるのは実部か虚部か、また微小変化の方向は正の方向か負の方向かについては、幾通りもの試行があり得る。したがって、複素係数列の微小変化とノイズに対するロバスト性の評価を幾通りも試行し、それらの中で最も良い結果をもたらす複素係数列を暫定的な最善の複素係数列として漸次採用していくことができる。
【００８７】
例えばＭ個の要素を持つ複素係数列について、そのＬ番目の係数の実部または虚部を、正または負の方向に微小量変化させる場合、その取り方は４×Ｍ通りある。また、或るＬ番目の係数を変化させた場合に、Ｌ番目ではない別の係数について、その実部および虚部の双方または何れか一方を微小量だけ増加または減少させることで、複素係数列を微小に補正して新たな係数列を作成することができる。この際に、Ｌ番目の係数の微小変化分をどの係数で補償するかの選び方は、（Ｍ−１）通りある。したがって、ある係数の実部または虚部を正または負の方向に微小量変化させ、その変化分を補償するために別の係数を微小量変化させる場合、４×Ｍ×（Ｍ−１）通りの組み合わせが存在する。これらの組み合わせ全てについて、ノイズに対するロバスト性を評価し、組み合わせの中で最も良い結果をもたらす複素係数列を暫定的な最善の複素係数列として漸次採用してもよい。
【００８８】
暫定的な最善の複素係数列を求めることができた後には、その暫定的な最善の複素係数列を起点として、よりノイズに対してロバストな複素係数を同様に求めることができる。この試行をコンピュータ上で反復的に繰り返すことで、最も良い結果をもたらすと予想される複素係数列を求めることができる。
【００８９】
なお、この複素係数列を求めるにあたって必要な情報は、オフセット位相量α１〜αＭの数値であって、実際の干渉画像そのものではない。したがって、ノイズに対するロバスト性の評価においてコンピュータ上でシミュレーションされたノイズを含む干渉画像を用い、最終的な画像再構成においては実際に取得された干渉画像を用いても問題はない。
【００９０】
１２個の干渉画像Ｉ１〜Ｉ12の複素係数列（coef_1  …  coef_12）に関して、第２の方法で得られた複素係数列を初期値とし、よりノイズに対してロバストな複素係数列を反復的に求めていくことによって最適化計算によって各要素の複素数を決定した。図２９は、第４の方法で得られた位相画像である。図３０は、図９の正解画像に対する図２９の位相画像の誤差を表す画像である。誤差の標準偏差は２２.０ミリラジアンであった。
【００９１】
図３１は、第２〜第４の各方法によるシミュレーションで用意した１２個の干渉画像Ｉ１〜Ｉ12のオフセット位相シフト量、組分け及び複素係数列を纏めた表である。以上のシミュレーション結果から、次のことが分かる。
【００９２】
干渉画像にノイズが存在していない理想的な条件では、第１〜第４の方法の何れによっても上記（１１）式に基づいて位相画像を生成することができる。しかし、実際の干渉画像にはノイズが存在し、その中でも光のショットノイズは物理的に避けられない。ノイズを含んだ干渉画像を用いたシミュレーションの結果、生成された位相画像の精度は、第４の方法で最も良好であり、これに第２の方法および第１の方法が続いた。第３の方法では、位相画像の精度は非常に悪かった。図３２は、第１〜第４の各方法によるシミュレーションで得られた誤差の標準偏差を比較して示すグラフである。
【００９３】
図２８に示したようにオフセット位相シフト量の変化の速度が撮像器の撮像フレームレートと比べて遅い場合、取得順に並べられた干渉画像は、第３の方法のＪ１〜Ｊ12のような順番で並んでしまう。このような場合、単純に得られた画像の順で３個の干渉画像を含む組を作って位相画像を生成すると、好ましくない結果をもたらす。すなわち、３個の干渉画像から複素係数列を求める方法では、得られた干渉画像を適切に再配置することが好ましい。
【００９４】
コンピュータ上での反復計算を用いた最適化によっても、上記（１１）式の複素係数列を求めることができる。この場合、３個の干渉画像から複素係数列を求める場合の最良な結果と比べても、更にノイズの少ない位相画像を生成することできる。図３３は、第２の方法による複素係数列および第４の方法による複素係数列を複素平面上にプロットした図である。
【００９５】
これまでに説明した実施形態は、Ｍ個の干渉画像それぞれについて画像全体のオフセット位相シフト量を推定し（オフセット位相シフト量推定ステップ）、Ｍ個の干渉画像のうちからオフセット位相シフト量が互いに異なるＮ個の干渉画像を含む組を設定し（組設定ステップ）、その組に含まれるＮ個の干渉画像に基づいて位相画像を生成する（位相画像生成ステップ）、というものであった。また、好適には、上記のような組を複数設定して（組設定ステップ）、これら複数の組それぞれについて各組に含まれるＮ個の干渉画像に基づいて位相画像を生成し（位相画像生成ステップ）、これら複数の組それぞれで生成された位相画像についてオフセット位相シフト量を補正した後に積算する（積算ステップ）、というものであった。
【００９６】
処理対象の各干渉画像は、位相画像取得装置により取得された画像の全体であってもよいし、位相画像取得装置により取得された画像のうちの一部である部分干渉画像であってもよい。位相画像取得装置により取得された干渉画像を幾つかの部分干渉画像に分割して、各々の部分干渉画像について上記の各ステップを行って位相画像を生成してもよい。このとき、部分干渉画像によって組分けが異なっていてもよい。
【００９７】
これまで説明した実施形態は、処理対象の干渉画像（または部分干渉画像）においてオフセット位相シフト量が時間的に変化するものの画素位置に依らず一定であると見做すことができる場合に好適なものである。例えば、位相画像取得装置の被撮像視野の範囲が容器の面積と比べて十分に小さい場合、取得される干渉画像においてオフセット位相量シフト量は画素位置に依らず一定であると見做すことができる。
【００９８】
これに対して、処理対象の干渉画像（または部分干渉画像）においてオフセット位相シフト量が画素位置に依らず一定であると見做すことができず、オフセット位相シフト量の分布の不均一性が無視できない場合がある。例えば、位相画像取得装置の被撮像視野の範囲が容器の面積と同程度である場合、水面の揺らぎにより、取得される干渉画像においてオフセット位相量シフト量の分布の不均一性が無視できないことがある。このような場合には、以下に説明する実施形態のようにして干渉画像から位相画像を生成するのが好ましい。
【００９９】
以下では、干渉画像を画素単位に分割して各ステップを行って位相画像を生成する方法について、シミュレーション結果とともに説明する。図３６および図３７は、シミュレーションにおいて用いた位相画像を示す図である。図３６は、位相物体の位相画像を示す。図３７は、干渉画像取得時に傾斜が存在する場合の位相画像を示す。図３７に示される位相画像をφtrue(x,y)とする。本シミュレーションでは、このような位相画像で表される位相物体について位相画像取得装置により３６個の干渉画像Ｉ１〜Ｉ３６を取得するものとする。また、干渉画像取得時に二光束干渉計において二光束の間の光路長差が例えば水面の揺らぎにより時間的に変動しているとする。
【０１００】
３６個の干渉画像Ｉ１〜Ｉ３６のうち第ｎの干渉画像Ｉｎ(x,y)は下記（１２）式で表されるものとする。この（１２）式の干渉画像Ｉｎ(x,y)は、真の位相画像φtrue(x,y)にオフセット位相シフト量分布αｎ(x,y)が加えられた場合に得られる干渉画像に、ランダムノイズが加えられたものである。ランダムノイズとして、０〜０.０５の範囲の一様乱数を用いた。オフセット位相シフト量分布αｎ(x,y)は、液面の揺らぎを模擬して下記（１３）式で表されるものとする。この（１３）式中の係数p1(n)，p2(n)，p3(n)およびp4(n)それぞれとして乱数を用いた。
【０１０１】
【数１２】
2019124679.tif 000014
【０１０２】
【数１３】
2019124679.tif 000015
【０１０３】
（１３）式の右辺の第１項は干渉縞の平行移動をもたらし、第２項および第３項は干渉縞の傾斜の変化をもたらし、第４項は干渉縞の二次元的な歪みをもたらす。図３８（ａ）は、位置(x,y)に依らず位相値が１である場合の位相分布関数を示す図である。図３８（ｂ）は、ｘ座標値に比例する位相値を有する場合の位相分布関数を示す図である。図３８（ｃ）は、ｙ座標値に比例する位相値を有する場合の位相分布関数を示す図である。図３８（ｄ）は、ｘ２＋ｙ２−１の値に比例する位相値を有する場合の位相分布関数を示す図である。なお、ｘ，ｙそれぞれの範囲を−１〜＋１の範囲で規格化している。
【０１０４】
上記（１３）式は、これらの位相分布関数それぞれに係数p1(n)，p2(n)，p3(n)，p4(n)を乗じたものの総和に相当している。図３９（ａ）〜（ｄ）は、オフセット位相シフト量分布αｎ(x,y)の例を示す図である。図４０は、図３９（ａ）〜（ｄ）のオフセット位相シフト量分布αｎ(x,y)の各例において設定した係数p1(n)，p2(n)，p3(n)，p4(n)の各値を纏めた表である。
【０１０５】
なお、オフセット位相シフト量分布αｎ(x,y)は、上記（１３）式では定数項、ｘまたはｙに比例する項、ｘまたはｙの二乗に比例する項を含むものであるが、これらに加えて、更に高次の項、指数関数の項、三角関数（sin、cos）の項を含んでいてもよい。一般に、オフセット位相シフト量分布αｎ(x,y)は、有限個の空間的な分布関数の線形和で表すことができる。
【０１０６】
図４１は、シミュレーションで用いた３６個の干渉画像のうちの３個の干渉画像を示す図である。これらの干渉画像は、真の位相画像φtrue(x,y)を用いて上記（１２）式および（１３）式に基づいて作成したものである。
【０１０７】
演算部３は、次のようなオフセット位相シフト量推定ステップ、組設定ステップ、位相値算出ステップおよび位相画像生成ステップを行うことで、干渉画像から位相画像を生成することができる。
【０１０８】
オフセット位相シフト量推定ステップにおいて、各干渉画像についてオフセット位相シフト量の分布を推定する。オフセット位相シフト量推定ステップでは、例えば次のように、フーリエ変換による空間的縞解析法による位相分布の算出処理、基準となる位相分布との差分処理、および、平滑化処理を順に行うのが好ましい。
【０１０９】
フーリエ変換による空間的縞解析法（以下「フーリエ変換法」という。）は、非特許文献３に記載されているとおり、干渉画像を二次元フーリエ変換してフーリエ変換像を作成し、このフーリエ変換像におけるスペクトルピーク位置を含む一定領域の空間周波数情報のみを用いて二次元逆フーリエ変換することで、二次元の位相分布画像を得るものである。これにより得られる二次元の位相分布画像は、空間的周波数が低く、なだらかである。
【０１１０】
一般に、フーリエ変換法では、対象物の空間的周波数分布に比べて干渉縞の空間的周波数が十分に高い必要があり、レーザ光の使用が不可欠である。しかし、レーザ光を用いた場合、スペックルノイズの問題が生じる。干渉縞の空間的周波数を高くしてインコヒーレント光を用いる場合には、画像の一部の領域にしか干渉縞が得られない。本実施形態では、対象物の空間的周波数分布に比べて干渉縞の空間的周波数が低い場合に特に有効である。そのような制限が発生するのは、特に光源のスペクトル幅が広いインコヒーレント光による干渉計測の場合である。
【０１１１】
図４２は、或る干渉画像Ｉｎ(x,y)を二次元フーリエ変換することで得られたフーリエ変換像を示す図である。図４３は、図４２のフーリエ変換像の原点付近を拡大して示す図である。干渉画像における干渉縞の空間的周波数が低いので、フーリエ変換像においてスペクトルピークは原点の近くに存在する。また、フーリエ変換像において原点を挟んで互いに対称な位置に２つのスペクトルピークが存在する。
【０１１２】
図４４は、図４２のフーリエ変換像における一方のスペクトルピーク位置を含む一定領域（図４３中の破線の円で示した領域）の空間周波数情報のみを用いて二次元逆フーリエ変換することで得られた二次元の位相分布画像を示す図である。ここで得られる位相分布画像（図４４）は、干渉縞の空間的周波数より高い周波数成分を有する構造を示すものではなく、空間的周波数が低く、なだらかである。この位相分布画像（図４４）は、正解画像（図３７）の低周波数成分しか表していない。
【０１１３】
３６個の干渉画像Ｉ１〜Ｉ３６それぞれからフーリエ変換法により位相分布画像を求めた後、これら３６個の位相分布画像について二次元位相アンラッピング後に平均をとって傾き補正を行うと、これにより得られる画像は、従来法による位相分布の再生結果となる。図４５は、従来のフーリエ変換法により得られた位相画像の中央領域を示す図である。この位相画像（図４５）は、正解画像（図３６）と比べると非常に画質が劣化したものとなっており、位相分布を再現できていない。
【０１１４】
本実施形態では、各干渉画像Iｎ(x,y)からフーリエ変換法により得られた位相分布画像φLPFｎ(x,y)を、該干渉画像Iｎ(x,y)のオフセット位相シフト量分布の推定に用いる。３６個の位相分布画像φLPF１(x,y)〜φLPF３６(x,y)は一般に互いに異なっている。そこで、基準となる位相分布画像φLPF0(x,y)に対して相対的な位相分布画像ΔφLPF１(x,y)〜ΔφLPF３６(x,y)を求める。基準となる位相分布画像φLPF0(x,y)は、３６個の位相分布画像φLPF１(x,y)〜φLPF３６(x,y)のうちの何れかの画像であってもよいし、３６個の位相分布画像φLPF１(x,y)〜φLPF３６(x,y)を平均した画像であってもよい。相対的位相分布画像ΔφLPFｎ(x,y)は下記（１４）式で表される。ここで、argは複素数の偏角を取る演算であり、ｉは虚数単位であり、expは指数関数である。
【０１１５】
【数１４】
2019124679.tif 000016
【０１１６】
図４６は、基準となる位相分布画像φLPF0(x,y)の一例を示す図である。図４７は、基準となる位相分布画像φLPF0(x,y)に対する相対的な位相分布画像ΔφLPFｎ(x,y)の一例を示す図である。相対的な位相分布画像ΔφLPFｎ(x,y)は、上記（１３）式によって近似することができる。相対的な位相分布画像ΔφLPFｎ(x,y)を上記（１３）式で近似したときの係数p1，p2，p3，p4の各値は、最小二乗法により求めることができる。この近似により、各干渉画像Ｉｎ(x,y)についてオフセット位相シフト量分布αｎ(x,y)を推定することができる。
【０１１７】
図４８（ａ）は、図４７の相対的な位相分布画像ΔφLPFｎ(x,y)の中央領域における相対的な位相分布を示す図である。図４８（ｂ）は、図４８（ａ）の相対的な位相分布を近似して得られた相対的な位相分布を示す図である。この例では、p1=3.11、p2=-0.21、p3=-0.31、p4=-0.06であった。なお、相対的な位相分布画像ΔφLPFｎ(x,y)を上記（１３）式で近似するのは、相対的な位相分布画像ΔφLPFｎ(x,y)を平滑化するためである。他の手法による平滑化処理を行ってもよく、例えは平滑化フィルタを用いた処理を行ってもよい。
【０１１８】
このようにして、オフセット位相シフト量推定ステップにおいて、各干渉画像Ｉｎ(x,y)についてオフセット位相シフト量分布αｎ(x,y)を推定することができる。以降は、各画素位置について組設定ステップおよび位相値算出ステップを行った後、位相画像生成ステップを行う。
【０１１９】
３６個の干渉画像I１(x,y)〜I３６(x,y)それぞれにおいて、任意の位置（x1,y1）にある画素の輝度値Iｎ(x1,y1)およびオフセット位相シフト量αｎ(x1,y1)を互いに対応させることができる。したがって、画素位置（x1,y1）毎に、３６個の画素（輝度値Iｎ(x1,y1)、オフセット位相シフト量αｎ(x1,y1)）のうちからオフセット位相シフト量が互いに異なるＮ個の画素を含む組を設定し（組設定ステップ）、その組に含まれるＮ個の画素に基づいて上記（２）式により位相値を求めることができる（位相値算出ステップ）。なお、組の設定および位相値の算出は、前述した方法と同様である。そして、各画素位置の位相値に基づいて位相画像を生成することができる（位相画像生成ステップ）。
【０１２０】
また、画素位置（x1,y1）毎に、このような組を複数設定し（組設定ステップ）、複数の組それぞれについて各組に含まれるＮ個の画素に基づいて位相値を求め（位相値算出ステップ）、複数の組それぞれで求めた位相値についてオフセット位相シフト量を補正した後に積算して（積算ステップ）、各画素位置の位相値の積算値に基づいて位相画像を生成してもよい（位相画像生成ステップ）。さらに、画素位置（x1,y1）毎に、３６個の画素（輝度値Iｎ(x1,y1)、オフセット位相シフト量αｎ(x1,y1)）をオフセット位相シフト量について昇順または降順に並び替えた後に複数の組に順に振り分けてもよい（組設定ステップ）。
【０１２１】
以下では、画素位置（x1,y1）毎に３６個の画素（輝度値Iｎ(x1,y1)、オフセット位相シフト量αｎ(x1,y1)）をオフセット位相シフト量について昇順に並び替えた後に複数の組に順に振り分けた場合について、シミュレーション結果を説明する。各組に含まれる画素の数を３とし、組の数を１２とした。昇順に並び替えた後に、第１画素，第１３画素および第２５画素を第１組に振り分け、第２画素，第１４画素および第２６画素を第２組に振り分け、以下同様にして、第１２画素，第２４画素および第３６画素を第１２組に振り分けた（組設定ステップ）。
【０１２２】
そして、１２組それぞれについて各組に含まれる３個の画素に基づいて位相値を求め（位相値算出ステップ）、１２組それぞれで求めた位相値についてオフセット位相シフト量を補正した後に平均値を求めた（積算ステップ）。このようにして画素位置（x1,y1）毎に位相値を求め、各画素位置の位相値に基づいて位相画像を生成した（位相画像生成ステップ）。
【０１２３】
画素位置（x1,y1）によって３６個の画素の組分けが異なる場合がある。しかし、画素位置（x1,y1）毎に高精度の位相値が得られるように組分けをすることができるので、最終的に得られる位相画像も高精度のものとなる。
【０１２４】
図４９および図５０は、シミュレーションにおいて生成された位相画像を示す図である。図４９は、位相アンラピングおよび傾き補正をする前の位相画像を示す。図５０は、位相アンラピングおよび傾き補正をした後の位相画像の中央領域を示す。これらの図に示される位相画像は、図４５に示された従来のフーリエ変換法による位相画像と比べて、高周波成分まで正確に再生されていることが分かる。
【０１２５】
図５１は、シミュレーションにおいて生成された位相画像と正解画像（図３６）との差分の画像を示す図である。画像の周縁領域では主としてオフセット位相シフト量の推定誤差に因り精度が幾らか低いものの、画像の中央領域では高精度に位相分布が再生されている。位相画像（４００×４００画素）のうちの中央領域（２００×２００画素）における誤差の標準偏差は、５３ミリラジアンであった。例えば光源の波長が５５０ナノメートルである場合、この誤差量は４.６ナノメートルの光路長差に相当する。
【０１２６】
なお、各干渉画像Ｉｎ(v,y)について推定されたオフセット位相シフト量分布αｎ(x1,y1)において上記（１３）式の右辺第１項が支配的である場合（すなわち、オフセット位相シフト量の空間的変化が小さい場合）には、干渉画像を幾つかの部分干渉画像に分割して、部分干渉画像毎に近似的なオフセット位相シフト量を設定し、各々の部分干渉画から位相画像を生成してもよい。このとき、部分干渉画像によって組分けが異なっていてもよい。このようにすることで、処理の負荷を軽減することができる。
【符号の説明】
【０１２７】
１…位相画像取得装置、２…干渉画像取得部、３…演算部、１０…光源、１１，１２…レンズ、２０…二光束干渉計、２１…分岐用ビームスプリッタ、２２…合波用ビームスプリッタ、３１…ビームスプリッタ、３２…第１反射ミラー、４１…ビームスプリッタ、４２…第２反射ミラー、４３…ピエゾ素子、４４…ステージ、５１，５２…レンズ、５３…チューブレンズ、６１…撮像器、７２…制御部。
【図１】
2019124679.tif 000017
【図２】
2019124679.tif 000018
【図３】
2019124679.tif 000019
【図４】
2019124679.tif 000020
【図５】
2019124679.tif 000021
【図６】
2019124679.tif 000022
【図７】
2019124679.tif 000023
【図８】
2019124679.tif 000024
【図９】
2019124679.tif 000025
【図１０】
2019124679.tif 000026
【図１１】
2019124679.tif 000027
【図１２】
2019124679.tif 000028
【図１３】
2019124679.tif 000029
【図１４】
2019124679.tif 000030
【図１５】
2019124679.tif 000031
【図１６】
2019124679.tif 000032
【図１７】
2019124679.tif 000033
【図１８】
2019124679.tif 000034
【図１９】
2019124679.tif 000035
【図２０】
2019124679.tif 000036
【図２１】
2019124679.tif 000037
【図２２】
2019124679.tif 000038
【図２３】
2019124679.tif 000039
【図２４】
2019124679.tif 000040
【図２５】
2019124679.tif 000041
【図２６】
2019124679.tif 000042
【図２７】
2019124679.tif 000043
【図２８】
2019124679.tif 000044
【図２９】
2019124679.tif 000045
【図３０】
2019124679.tif 000046
【図３１】
2019124679.tif 000047
【図３２】
2019124679.tif 000048
【図３３】
2019124679.tif 000049
【図３４】
2019124679.tif 000050
【図３５】
2019124679.tif 000051
【図３６】
2019124679.tif 000052
【図３７】
2019124679.tif 000053
【図３８】
2019124679.tif 000054
【図３９】
2019124679.tif 000055
【図４０】
2019124679.tif 000056
【図４１】
2019124679.tif 000057
【図４２】
2019124679.tif 000058
【図４３】
2019124679.tif 000059
【図４４】
2019124679.tif 000060
【図４５】
2019124679.tif 000061
【図４６】
2019124679.tif 000062
【図４７】
2019124679.tif 000063
【図４８】
2019124679.tif 000064
【図４９】
2019124679.tif 000065
【図５０】
2019124679.tif 000066
【図５１】
2019124679.tif 000067