(58)【調査した分野】(Int.Cl.,DB名)
【背景技術】
【0003】
従来、物体の内部情報、つまり後方散乱、反射率分布及び屈折率分布の微分構造を非破壊、高分解能で捉えるために、OCTを用いることが行われている。
【0004】
医療分野等で用いられる非破壊断層計測技術の1つとして、光断層画像化法「光コヒーレンストモグラフィー」(OCT)がある(特許文献1参照)。OCTは、光を計測プローブとして用いるため、被計測物体の反射率分布、屈折率分布、分光情報、偏光情報(複屈折率分布)等が計測できるという利点がある。
【0005】
基本的なOCT43は、マイケルソン干渉計を基本としており、その原理を
図8で説明する。光源44から射出された光は、コリメートレンズ45で平行化された後に、ビームスプリッター46により参照光と物体光に分割される。物体光は、物体アーム内の対物レンズ47によって被計測物体48に集光され、そこで散乱・反射された後に再び対物レンズ47、ビームスプリッター46に戻る。
【0006】
一方、参照光は参照アーム内の対物レンズ49を通過した後に参照鏡50によって反射され、再び対物レンズ49を通してビームスプリッター46に戻る。このようにビームスプリッター46に戻った物体光と参照光は、物体光とともに集光レンズ51に入射し光検出器52(フォトダイオード等)に集光される。
【0007】
OCTの光源44は、時間的に低コヒーレンスな光(異なった時刻に光源から出た光同士は極めて干渉しにくい光)の光源を利用する。時間的低コヒーレンス光を光源としたマイケルソン型の干渉計では、参照アームと物体アームの距離がほぼ等しいときにのみ干渉信号が現れる。この結果、参照アームと物体アームの光路長差(τ)を変化させながら、光検出器52で干渉信号の強度を計測すると、光路長差に対する干渉信号(インターフェログラム)が得られる。
【0008】
そのインターフェログラムの形状が、被計測物体48の奥行き方向の反射率分布を示しており、1次元の軸方向走査により被計測物体48の奥行き方向の構造を得ることができる。このように、OCT43では、光路長走査により、被計測物体48の奥行き方向の構造を計測できる。
【0009】
このような軸方向(A方向)の走査のほかに、横方向(B方向)の機械的走査(Bスキャン)を加え、2次元の走査を行うことで被計測物体の2次元断面画像が得られる。この横方向の走査を行う走査装置としては、被計測物体を直接移動させる構成、物体は固定したままで対物レンズをシフトさせる構成、被計測物体も対物レンズも固定したままで、対物レンズの瞳面付近においたガルバノミラーの角度を回転させる構成等が用いられている。
【0010】
以上の基本的なOCTが発展したものとして、光源の波長を走査してスペクトル干渉信号を得る波長走査型OCT(Swept Source OCT、略して「SS−OCT」という。)と、分光器を用いてスペクトル信号を得るスペクトルドメインOCTがある。後者としてフーリエドメインOCT(Fourier Domain OCT、略して「FD−OCT」という。特許文献2参照)、及びPS−OCT(特許文献3参照)がある。
【0011】
波長走査型OCTは、高速波長スキャニングレーザーにより光源の波長を変え、スペクトル信号と同期取得された光源走査信号を用いて干渉信号を最配列し、信号処理を加えることで3次元光断層画像を得るものである。なお、光源の波長を変える手段として、モノクロメーターを利用したものでも、波長走査型OCTとして利用可能である。
【0012】
フーリエドメインOCTは、被計測物体からの反射光の波長スペクトルを、スペクトロメーター(スペクトル分光器)で取得し、このスペクトル強度分布に対してフーリエ変換することで、実空間(OCT信号空間)上での信号を取り出すことを特徴とするものであり、このフーリエドメインOCTは、奥行き方向の走査を行う必要がなく、x軸方向の走査を行うことで被計測物体の断面構造を計測可能である。
【0013】
PS−OCTは、B−スキャンと同時に直線偏光したビームの偏光状態を連続変調し、試料(被検物体)のもつ偏光情報を捉え、試料のより微細な構造および屈折率の異方性を計測可能とする光コヒーレンストモグラフィー装置である。
【0014】
さらに詳しくは、PS−OCTは、フーリエドメインOCTと同様に、被計測物体からの反射光の波長スペクトルをスペクトル分光器で取得するものであるが、入射光及び参照光をそれぞれ1/2波長板、1/4波長板等を通して水平直線偏光、垂直直線偏光、45°直線偏光、円偏光として、被計測物体からの反射光と参照光を重ねて1/2波長板、1/4波長板等を通して、例えば水平偏光成分だけをスペクトル分光器に入射させて干渉させ、物体光の特定偏光状態をもつ成分だけを取り出してフーリエ変換するものである。このPS−OCTも、奥行き方向の走査を行う必要がない。
【発明の概要】
【発明が解決しようとする課題】
【0016】
一般的に、OCTを生体計測へ応用することを考える場合、物体の屈折率分布の微分構造は、非破壊、高分解能で捉えることはできるものの、繊維状の構造(繊維の伸長方向等)や歯のエナメル質の相違に起因する複屈折による偏光依存性を持つ生物試料の測定においては、解像度の低下とともに、構造を捉えられない。
【0017】
この点、PS−OCTは、ある特定部分からの散乱光成分とある偏光状態の参照光とを干渉させて、その干渉成分には偏光特性が強く反映され、位相遅延情報(リターデーション量)も得ることができ、その結果、奥行き方向の断面のある特定部分の偏光情報、複屈折情報を捉えることが可能となる。従って、PS−OCTは、眼科、循環器、呼吸器、消化器など広く非侵襲な組織判別を必要とする臨床分野における検査装置として有用である。
【0018】
このように、PS−OCTは、組織の偏光特性を非侵襲に可視化する光トモグラフィー装置として、組織コントラストという定性観察には有用である。しかしながら、本発明者らは、PS−OCTの高精度化を目的に鋭意研究開発する過程において、PS−OCTは、病変部位を観察し、病期の正確な分類などの定量観察においては精度が必ずしも十分とは言えず、さらなる改善すべき点があるということが判った。
【0019】
本発明は、上記従来のPS−OCTについて得られた改善すべき点を解決することを目的とするものであり、PS−OCTにおいて計測されたデータを非線形に補正し、PS−OCTの定量解析能力を高め、病期の正確な定量診断等を可能とするプログラム及び該プログラムを搭載したPS−OCTシステムを実現することを課題とする。
【0020】
より具体的には、PS−OCTで計測されるリターデーション(位相遅延量、位相差の量)は誤差を含んでいて、真値の周りに分布(これをリターデーションの分布、或いは位相差の分布と称す)している。この分布は一般的に対称ではなく、必ずしも真値が中心にあるとは限らない非対称である場合も含まれる。従って、標準的な平均値をとる位相解析では正しいリターデーションを得ることができない。
【0021】
本発明は、このような従来のPS−OCTの問題を解決しようとするものであり、リターデーションが誤差を含みノイズとなり、その分布が正規分布でもなく、真の値の周りに対象に分布していない場合であっても、モンテカルロシミュレーションでノイズの特性を解析することによって得られた分布変換関数を用いて、計測データを変換することにより、系統的な誤差を除去し、ノイズに埋もれた真の値を推測し、PS−OCTの画像をより鮮明に補正するプログラムを実現するとともに、このプログラムを搭載したPS−OCTシステムを実現することを課題とする。
【課題を解決するための手段】
【0022】
本発明は上記課題を解決するために、PS−OCTで得られた計測データを処理するコンピュータに搭載されるプログラムであって、コンピュータを、PS−OCTで計測したリターデーションについての位相差の分布のデータを、非線形の変換関数によって、対称な位相差の分布のデータに変換し、平均することにより、真の位相量を求める手段として機能させることを特徴とするプログラムを提供する。
【0023】
本発明は上記課題を解決するために、PS−OCTで得られた計測データを処理するコンピュータに搭載されるプログラムであって、コンピュータを、予め、0〜πまでの複数の位相量それぞれについての分布から成る位相差の分布のデータのセットを、異なるESNR毎に複数のセット作成し記憶する手段、PS−OCTで計測した位相差の分布のデータを対称な位相差の分布のデータに変換する非線形の変換関数の係数を決める手段、PS−OCTで計測した位相差の分布のデータからESNRを見積もり、該見積もられたESNRに基づき、前記異なるESNR毎に作成した位相差の分布の複数のセットのいずれかを選択する手段、PS−OCTで計測した位相差の分布のデータを、前記選択された位相差の分布のセットに応じて決められた係数の下で、前記変換関数によって対称な位相差の分布のデータに変換する手段、及び前記対称な位相差の分布のデータを平均することにより、真の位相量を算出する手段、として機能させることを特徴とするプログラム。
【0024】
本発明は上記課題を解決するために、PS−OCTと、該PS−OCTで得られた計測データを処理するコンピュータとを備えたPS−OCTシステムであって、前記コンピュータは、入力装置、出力装置、CPU及び記憶装置を備え、PS−OCTで計測した位相差の分布のデータを、非線形の変換関数によって、対称な位相差の分布のデータに変換し、平均することにより、真の位相量を求める手段として機能することを特徴とするPS−OCTシステムを提供する。
【0025】
本発明は上記課題を解決するために、PS−OCTと、該PS−OCTで得られた計測データを処理するコンピュータとを備えたPS−OCTシステムであって、前記コンピュータは、入力装置、出力装置、CPU及び記憶装置を備えており、予め、0〜πまでの複数の位相量それぞれについての分布から成る複数の位相差の分布のセットを、異なるESNR毎に複数のセット作成し記憶装置に記憶させる手段、PS−OCTで計測した位相差の分布のデータを対称な位相差の分布のデータに変換する非線形の変換関数の係数を決める手段、PS−OCTで計測した位相差の分布のデータからESNRを見積もり、該見積もられたESNRに基づき、前記異なるESNR毎に作成した位相差の分布の複数のセットのいずれかを選択する手段、PS−OCTで計測した位相差の分布のデータを対称な位相差の分布を、前記選択された位相差の分布のセットに応じて決められた係数の下で、前記変換関数によって対称な位相差の分布のデータに変換する手段、前記対称な位相差の分布のデータを平均することにより、真の位相量を算出する手段、として機能することを特徴とするPS−OCTシステムを提供する。
【0026】
前記非線形の変換関数は、後記する式(2)である。
【0027】
前記変換関数の係数を決める手段は、前記非線形の変換関数について、0〜πまで、π/mきざみで変え、次のm+1の連立方程式である後記する式(3)を作成し、該連立方程式を解いて、変換関数の係数biのテーブルを作る構成である。但し、mは、πをきざむ回数である。
【0028】
前記連立方程式を解いて前記変換関数の係数を決める手段においては、後記する式(4)の行列DMの疑似逆行列DM
+を数値的に求めて、後記する式(5)で計算し、最適の変換関数の係数biを一義的に決める構成としてもよい。
【0029】
前記変換関数の係数を決める手段は、変分法を用い、後記する式(6)で示される二乗誤差が最小となるように最適の変換関数の係数biを一義的に決める構成としてもよい。
【発明の効果】
【0030】
本発明によれば、PS−OCTにおいて計測されたデータを非線形に補正し、PS−OCTの定量解析能力を高めることができるので、従来のPS−OCTは、病変部位の形態観察等の用途に限られていたが、本発明のPS−OCTシステムは、病変部位の病期を含め、正確な定量診断を可能とし、コンピュータ診断のための有用な手段として利用することが可能となる。
【発明を実施するための形態】
【0032】
本発明に係るPS−OCTの計測データを補正するプログラム及び該プログラムを搭載したPS−OCTシステムを実施するための形態を図面を参照して、以下説明する。
【0033】
本発明に係るPS−OCTシステムは、PS−OCTと、PS−OCTで得られた画像データを処理する画像処理装置とを備えている。画像処理装置は、通常のコンピュータが使用され、本発明に係るプログラムは、PS−OCTで得られた画像データを補正する手段として機能させるものである。
【0034】
(PS−OCT)
PS−OCTは、すでに特許第4344829号公報等で知られているが、本発明の前提となる技術であるから、その概要を説明する。
【0035】
本発明に係る偏光感受光画像計測装置は、Bスキャンと同時に(同期して)光源からの偏光ビーム(偏光子により直線的に偏光されたビーム)をEO変調器(偏光変調器、電気光学変調器)によって連続的に変調し、この連続的に偏光を変調した偏光ビームを分けて、一方を入射ビームとして走査して試料に照射し、その反射光(物体光)を得ると共に、他方を参照光として、両者のスペクトル干渉によりOCT計測を行うものである。
【0036】
そして、このスペクトル干渉成分のうち、垂直偏光成分(H)と水平偏光成分(V)を同時に2つの光検出器で測定することにより、試料の偏光特性を表すジョーンズベクトルを得る(H画像とV画像)構成を特徴とするものである。
【0037】
図1は、本発明に係る偏光感受光画像計測装置の光学系の全体構成を示す図である。
図1に示す偏光感受光画像計測装置1は、光源2、偏光子3、EO変調器4、ファイバーカプラー(光カプラー)5、参照アーム6、試料アーム7、分光器8等の光学要素を備えている。この偏光感受光画像計測装置1の光学系は、光学要素が互いにファイバー9で結合されているが、ファイバーで結合されていないタイプの構造(フリースペース型)であってもよい。
【0038】
光源2は、広帯域スペクトルを有するスーパールミネッセントダイオード(SLD:Super Luminessent Diode)を使用する。なお、光源2は、パルスレーザでもよい。光源2には、コリメートレンズ11、光源2からの光を直線偏光にする偏光子3、進相軸を45°の方向にセットされたEO変調器(偏光変調器、電気光学変調器)4、集光レンズ13及びファイバーカプラー5が、順次、接続されている。
【0039】
EO変調器4は、進相軸を45°の方向に固定して、該EO変調器4にかける電圧を正弦的に変調することで、進相軸とそれに直交する遅相軸との間の位相差(リターデーション)を連続的に変えるもので、これにより、光源2から出て偏光子3で(縦)直線偏光となった光がEO変調器4に入射すると、上記変調の周期で、直線偏光→楕円偏光→直線偏光………などのように変調される。EO変調器4は、市販されているEO変調器を使用すればよい。
【0040】
ファイバーカプラー5には分岐するファイバー9を介して、参照アーム6と試料アーム7が接続されている。参照アーム6には、偏波コントローラ(polarization controller)10、コリメートレンズ11、偏光子12、集光レンズ13及び参照鏡(固定鏡)14が、順次、設けられている。参照アーム6の偏光子12は、上記のとおり偏光状態を変調しても参照アーム6から戻ってくる光の強度が変化しないような方向を選択するために用いている。この偏光子12の方向(直線偏光の偏光方向)の調整は偏波コントローラ10とセットで行う。
【0041】
試料アーム7では、偏波コントローラ15、コリメートレンズ11、固定鏡24、ガルバノ鏡16、集光レンズ13が、順次、設けられ、ファイバーカプラー5からの入射ビームが2軸のガルバノ鏡16により走査されて試料17に照射される。試料17からの反射光は物体光として再びファイバーカプラー5に戻り、参照光と重畳されて干渉ビームとして分光器8に送られる。
【0042】
分光器8は、順次接続される偏波コントローラ18、コリメートレンズ11、(偏光感受型体積位相ホログラフィック)回折格子19、フーリエ変換レンズ20、偏光ビームスプリッター21及び2つの光検出器22、23を備えている。この実施の形態では、光検出器22、23として、ラインCCDカメラ(1次元CCDカメラ)を利用する。ファイバーカプラー5から送られてくる干渉ビームは、コリメートレンズ11でコリメートされ、回折格子19によって干渉スペクトルに分光される。
【0043】
回折格子19で分光された干渉スペクトルビームは、フーリエ変換レンズ20でフーリエ変換され偏光ビームスプリッター21で水平及び垂直成分に分けられ、それぞれ2つラインCCDカメラ(光検出器)22、23で検出される。この2つラインCCDカメラ22、23は、水平および垂直偏光信号両方の位相情報を検知するために使われるので、2つのラインCCDカメラ22、23は同一の分光器の形成に寄与するものでなくてはならない。
【0044】
なお、光源2、参照アーム6、試料アーム7及び分光器8には、それぞれ偏波コントローラ10、15、18が設けられているが、これらは、光源2から参照アーム6、試料アーム7、分光器8に送られるそれぞれのビームの初期偏光状態を調整して、EO変調器4で連続的に変調された偏光状態が、参照光と物体光においても互いに一定の振幅と一定の相対偏光状態の関係が維持され、さらにファイバーカプラー5に接続された分光器8において一定の振幅と一定の相対偏光状態を保たれるようにコントロールする。
【0045】
また、2つラインCCDカメラ22、23を含む分光器8を校正するときはEO変調器4は止める。参照光をブロックし、スライドガラスと反射鏡を試料アーム7におく。この配置は水平および垂直偏光成分のピークの位置が同じであることを保証する。そして、スライドガラスの後ろの面と反射鏡からのOCT信号は2つの分光器8で検知される。OCT信号のピークの位相差はモニターされる。
【0046】
この位相差はすべての光軸方向の深さでゼロであるべきである。次に、信号は2つラインCCDカメラ22、23を含む分光器8で複素スペクトルを得るために、ウィンドウされ逆フーリエ変換される。この位相差はすべての周波数でゼロであるべきなので、これらの値をモニターすることによって2つラインCCDカメラ22、23の物理的な位置は位相差が最小になるようにアライメントされる。
【0047】
PS−OCTの原理:
本発明の特徴は、次のとおりである。光源2からの光を直線偏光し、この直線偏光されたビームをEO変調器4により連続的に偏光状態の変調を行う。即ち、EO変調器4は、進相軸を45°の方向に固定して、該EO変調器4にかける電圧を正弦的に変調することで、進相軸とそれに直交する遅相軸との間の位相差(偏光角:リターデーション)を連続的に変えるもので、これにより、光源2から出て直線偏光子で(縦)直線偏光となった光がEO変調器4に入射すると、上記変調の周期で、直線偏光→楕円偏光→直線偏光………などのように変調される。
【0048】
そして、直線偏光された偏光ビームをEO変調器4により連続的に偏光状態の変調を行うと同時に、Bスキャンを同期して行う。即ち、1回のBスキャンの間に、EO変調器4による偏光の連続的な変調を複数周期行う。ここで、1周期とは、偏光角(リターデーション)φが0〜2πと変化する期間である。要するに、この1周期の間に、偏光子からの光の偏光が、直線偏光(垂直偏光)→楕円偏光→直線偏光(水平偏光)………などのように連続的に変調する。
【0049】
このように偏光ビームの偏光を連続的に変調しながら、試料アーム7では、入射ビームをガルバノ鏡16により試料17に走査してBスキャンを行い、分光器8において、その反射光である物体光と参照光の干渉スペクトルについて、その水平偏光成分および垂直偏光成分を2つのラインCCDカメラ22、23で検出する。これにより、1回のBスキャンによって、それぞれ水平偏光成分及び垂直偏光成分に対応する2枚のA−Bスキャン画像が得られる。
【0050】
上記のとおり、1回のBスキャンの間に、偏光ビームの偏光の連続的な変調を複数周期行うが、各周期(1周期)の連続的な変調の間に2つのラインCCDカメラ22、23で検出した水平偏光成分および垂直偏光成分それぞれの偏光情報が1画素分の偏光情報となる。1周期の連続的な変調の間に2つのラインCCDカメラ22、23で偏光情報を検出タイミング信号に同期して行い、1周期に検出回数(取込回数)を、4回、8回等、適宜決めればよい。
【0051】
このようにして1回のBスキャンの間に得た2枚のA−Bスキャン画像のデータを、Bスキャン方向に1次元フーリエ変換を行う。すると、0次、1次、−1次のピークが出る。ここで、0次のピークをそれぞれ抽出し、そのデータのみを用いて逆フーリエ変換すると、H0、V0画像が得られる。同様に、1次のピークをそれぞれ抽出し、そのデータのみを用いて逆フーリエ変換すると、H1、V1画像が得られる。
【0052】
H0、H1画像から、試料17の偏光特性を示すジョーンズマトリックスの成分のうち、J(1,1)およびJ(1,2)を求める事ができる。そして、V0、V1画像から、試料17の偏光特性を示すジョーンズマトリックスの成分のうち、J(2,1)およびJ(2,2)を求める事ができる。
【0053】
なお、ここでは、ジョーンズマトリックスの成分J(1,1)、J(1,2)、J(2,1)、J(2,2)を特許文献1の記載に基づいて記載したが、それぞれ、本発明の後記する式(1)のジョーンズマトリックスの成分I
0H(z)、I
1H(z)、I
0V(z)、I
0V(z)に相当する。
【0054】
このようにして、1回のBスキャンにおいて4つの偏光特性を含む情報が得られる。そして、この4つの情報をそれぞれ、通常のFD−OCTと同様にAスキャン方向にフーリエ変換すると、1次のピークが試料17の深さ方向の情報を有し、しかもそれぞれ偏光特性に応じた4枚のA−B画像の計測データが得られる。
【0055】
(本発明のプログラム及びPS−OCTシステム)
上記構成のPS−OCTで得られた計測データは、画像処理装置として使用するコンピュータ30に入力される。このコンピュータ30は通常のコンピュータであり、
図2に示すように、入力部31、出力部32、CPU33、記憶装置34及びデータバス35を備えている。
【0056】
本発明に係るPS−OCTシステムに搭載されるプログラムは、コンピュータ30の記憶装置34に記憶されるプログラムであって、コンピュータ30に入力されたPS−OCTで得られた画像の計測データを補正し、PS−OCTの画像をより鮮明に補正する手段として、コンピュータ30を機能させるプログラムである。このプログラムを搭載することにより、本発明に係るPS−OCTシステムは、画像をより鮮明に補正する手段を備えたシステムとなる。
【0057】
即ち、従来はPS−OCTで計測された生データ (生の偏光位相遅延データ)から得られた画像を、直接、診断その他の各種の用途に用いていた。本発明のプログラム及びPS−OCTシステムでは、PS−OCTで計測された生データを、非線形関数で補正して高精度な生体複屈折計測を可能とするものであり、この非線形関数は、モンテカルロ法(Monte-Carlo analysis)で、予め本発明のプログラムを用い、コンピュータ30においてシミュレーションで設計される。
【0058】
以下、本発明のプログラム及び該プログラム搭載したPS−OCTシステムを詳細に説明する。
【0059】
図3は、リターデーション(位相差、位相遅延量、)の分布図(位相差の分布図)である。横軸は位相であり、縦軸はリターデーションの生じる頻度である。
図3の右欄は、PS−OCTで計測したリターデーションに基づく位相差の分布であり、
図3の左欄は、シミュレーションで作成した位相差の分布である。
【0060】
位相差の分布は、横軸をπ/mずつ刻んだm+1個のデータ(m+1個の位相差の分布)を1セットとして、ESNRを変えて複数のセット作成する。例えば、m=60として、横軸をπ/60ずつ刻んだ61個のデータを1セットとして、ESNRを変えて複数のセット作成する。
【0061】
図3では、2つの異なるESNRのセットについて、それぞれπ/2及びπについての位相差の分布を示す図である。
【0062】
ところで、PS−OCTによるリターデーションのノイズモデルは、ジョーンズマトリックスのそれぞれの要素に、加算的なノイズを設定して、次の式(1)で示される。ここで、Sは真のリターデーションの値であり、Nはノイズであり、Iは測定値である。添え字の0、1は、計測チャンネル(
図1の光検出器22、23参照)である。前記段落0052に説明したジョーンズマトリックスの成分は、真の値SにノイズNが加わった測定値Iに対応する。ここでzは深さ方向の座標で、ジョーンズマトリックスのそれぞれの要素は表面(z=0)から深さzまでのリターデーション等の積分量に相当する。
【0064】
計測データのノイズを表す評価値として、通常は、S/N比(SN比、SNR等とも呼ばれている)を使用するが、本発明では、より効果的にノイズを表す評価値としてESNRを使用する。本明細書では、ESNRはγとして表し、次の式で定義される。
1/γ=1/4(1/SNR
S,0+1/SNR
S,1+1/SNR
r,0+1/SNR
r,1 )
なお、この式中、sはPS−OCTのサンプルアームであり、rは参照アームである。0と1は、回数である。
【0065】
複数の独立したSNRをもつものを、単一の値で示すためにESNRを導入した。添え字の0、1は、計測チャンネル(
図1の光検出器22、23参照)である。SNR
S,0、 SNR
S,1、SNR
r,0、SNR
r,1 は、それぞれ添え字の組み合わせによって、サンプルアームから0チャンネルへ寄与するSNR、サンプルアームから1チャンネルへ寄与するSNR、参照アームから0チャンネルへ寄与するSNR、参照アームから1チャンネルへ寄与するSNRを表す。
【0066】
図3の右欄の上から2番目の図に示すように、PS−OCTで計測されるリターデーション(位相遅延量、位相差)は誤差を含んでいて、真値(π/2)の周りに分布している。そして、ESNRが小さくなると、
図3の右欄の上から1番目の図に示すように、誤差の幅を大きくなる。
【0067】
また、この位相差の分布は、
図3の右欄の上から3、4番目の図に示すように、必ずしも、対称ではなく、真値が中心にあるとは限らず、非対称である場合も含まれる。
【0068】
要するに、リターデーションが誤差を含みノイズとなり、位相差の分布が正規分布でもなく、真の位相量の周りに対象に分布していない場合も含まれる。このような場合は、標準的な平均値をとる位相解析では正しいリターデーションを得ることができない。
【0069】
本発明に係るプログラムは、モンテカルロシミュレーションでノイズの特性を解析することによって得られた分布変換関数を用いて、PS−OCTで計測された画像の計測データを変換することにより、系統的な誤差を除去し、ノイズに埋もれた真の位相量を推測し、PS−OCTの画像をより鮮明に補正する手段としてコンピュータ30を機能させるプログラムであり、本発明のシステムは、このプログラムを搭載し上記補正を行うPS−OCTシステムである。
【0070】
原理:
本発明に係るプログラムは、PS−OCTで計測された画像の計測データによる実測のリターデーションの分布(位相差の分布)を補正し、真の位相量δ
T(真のリターデーションの値)を求める手段とし、コンピュータ30を機能させるものであるが、まず、その原理(アルゴリズム)を、以下に説明する。
【0071】
本発明では、PS−OCTで計測された画像の計測データによる実測の位相量δ
Mから、非線形関数fを用いて対称な位相差の分布δ
E に変換する。要するに、対称な位相差の分布δ
E は、実測の位相量δ
M の冪関数であるδ
E=f(δ
M)として表される。この関数を本明細書では、変換関数と言い、次の式(2)で表される。これは、計測値の非対称な分布δ
Mを、正しい値を重みの中心とする分布δ
Eに変換するために行われる非線形な写像を意味している。
【0073】
この式(2)中、biは変換関数の係数である。最適な変換関数の係数bi(iは、0、1、2、…………n)は、ESNRγの関数であるが、次のように、連立方程式を解いて決めるか、或いは最小二乗法を用いて決め、bi(γ)のテーブルを作っておく。
【0074】
そのために、位相値を、0からπまでπ/mきざみで変え、モンテカルロ法によってδ
Mを求める。係数biは、式(3)に示す連立方程式を解いて決める。式(3)中、〈δ
M,π/m〉は、π/mきざみ、例えば、π/60きざみでシミュレートした、0からπにおけるリターデーションのシミュレーション値δ
Mを示す。ここで、mは、例えば、1、2、……60である。なお、式(3)中の〈 〉の意味は、多数のモンテカルロシミュレーションによって得られた値の平均値を示す。具体的には、式(3)を、式(4)に示すように行列で書き下し、係数bi(γ)を求める。
【0077】
そして、係数を決めるには、式(4)のシミュレーションで得られたδ
M,π/mnのn×mの行列DMの疑似逆行列DM
+を数値的に求めて、次の式(5)で計算し係数を求める。なお、式(4)中の〈 〉の意味は、多数のモンテカルロシミュレーションによって得られた値の平均値を示す。
【0079】
最小二乗法を用いて係数を決める手段は、変分法を用い、次の式(6)で示す、二乗誤差R
2が最小になるように、モンテカルロ法によって、最適の変換関数の係数bi(γ)を一義的に決めることができる。ここでも、位相値を、0からπまでπ/mきざみで、例えばm=60として、π/60きざみで変え、モンテカルロ法によってシミュレーション値δ
Mを求める。kは、位相をこきざみに計測した計測回を示し、0、1、2……mであり、例えばm=60等である。
【0081】
最適の変換関数の係数biが決まれば、対称な位相差の分布δ
Eを表す変換関数δ
E=f(δ
M)が決まり、PS−OCTによる実測の位相量δ
Mを代入すれば、δ
E が求まる。
【0082】
そして、δ
Eを次の式(7)に示すように平均する(式(7)中の〈 〉は平均の意味)ことにより、真の位相量δ
Tを求めることができる。平均は、一つの位相量を求めるために行った実測回数でおこなう。
【0084】
構成:
本発明の原理は以上のとおりであるが、この原理に基づき、本発明のPS−OCTのシステムにおいては、本発明に係るプログラムは、画像処理装置として使用するコンピュータ30を次の手段として機能させる。
【0085】
(1)位相差の分布の作成手段
位相差の分布の作成手段は、予め、
図3の左欄に示すようなリターデーションについての位相差の分布をシミュレーションで複数作成する手段である。
【0086】
具体的には、作成する位相差の分布は、π/m、例えばm=60として、横軸をπ/60ずつ刻んだ、61個のデータ(61個の位相差の分布)を1セットとして、ESNRを変えて位相差の分布のセットを複数セット作成する。この複数セットのデータは、作成後、記憶装置34に記憶される。
【0087】
(2)変換関数の係数の決定手段
変換関数の係数の決定手段は、複数の位相差の分布のセットについて、それぞれ変換関数の係数biを決める手段である。
【0088】
この変換関数の係数biは、前記したとおり、式(3)の連立方程式を解いて決めるか、或いは式(6)において最小二乗法を用いて決める。biが決められたら、そのテーブルを作成する。この変換関数の係数biのテーブルは、記憶装置34に記憶される。
【0089】
(3)位相差の分布の選択手段
位相差の分布の選択手段は、PS−OCTで計測された画像の実測の位相差の分布、位相量の補正に先立ち、その計測データ及びシステム設計等からESNRの値を見積もり、上記(1)に示す複数の位相差の分布(
図3の左欄のような位相差の分布)のセットについて、どのセットを用いるかを決定する。
【0090】
(4)変換関数による位相差の分布の変換手段
変換関数の算出手段は、前記記憶装置34に記憶された変換関数の係数bi又はそのテーブルのうち、選択手段で選択された位相差の分布のセットに対応する係数bi又はそのテーブルを用いて、式(2)のbiに代入し、PS−OCTで計測した実測のリターデーションの位相量δ
Mから、対称な位相差の分布を表す変換関数δ
E=f(δ
M)を求める手段である。
【0091】
(5)真の位相量算出手段
真の位相量算出手段は、上記のとおり求めた対称な位相差の分布を表すδ
E=f(δ
M)から、前記式(7)に基づき平均することにより、真の位相量δ
Tを算出する手段である。
【0092】
(試験例1)
本発明者らは、本発明のプログラム及び該プログラム搭載したPS−OCTシステムを評価するために、シミュレーション値と実測値を比較する比較試験を行った。この比較試験では、試料(計測対象物)は、ガラス板、1/8波長板及び1/4波長板であり、それぞれに生じるリターデーションの真値は、0、π/2及びπである。
【0093】
なお、シミュレーション値は、前記したとおり、予めbiの係数をモンテカルロシミュレーションで決めておき、これによって得られた係数を用いてこのシミュレーションデータを変換したものである。実測データを上記係数で変換したものが実測値である。モンテカルロシミュレーションにはランダムノイズが入っていることと、biを有限の項数で打ち切っているので、シミュレーションといえども完全には一致はしない。
【0094】
この試験では、PS−OCTによる計測のために、1.3μmのプローブを使用した。このPS−OCTは、空気中での深度解析は8.3μmであり、30,000ライン/秒の測定速度を有する。計測は、試料の前方に光強度を変えることの可能な非偏光フィルタを置くことによって、ESNRを変えて行った。
【0095】
図4は、シミュレーションと試験を比較した比較試験の結果を示す。点線は、真のリターデーションの値(0、π/2及びπ)を示す。実線(b)、(c)、(e)は、本発明による補正を行わずに、リターデーションの実測値を単に平均のみをとったものである。実線(a)、(d)、(f)は、本発明によるシミュレーションの結果であり、係数を作るのに用いたシミュレーション数は65536個であり、これは、係数を求めるためのデータでもあり、結果を得るために用いた位相データでもある。
【0096】
実線(b)、(c)、(e)及び(a)、(d)、(f)の近辺に位置する●■+のマークは、それぞれ試験によって得られたリターデーションの実測値及び該実測値を本発明により補正した試験補正値あって、ここでは64個の測定データを用いている。シミュレーションと試験補正値は、絶対値及びESNR依存性等において、ほぼ一致している。
【0097】
(試験例2)
本発明者らは、本発明の評価をするために、PS−OCTの評価に通常使用されている鶏の胸肉(ex−vivo)を試料として計測する比較試験を行った。
【0098】
この試験では、試料の10mmの部分について、512本のAラインから成るBスキャンを行った。要するに、A方向(深さ方向)の計測で得られるAラインを、Bスキャン(A方向に直交するB方向のスキャン)によって、512本得る。しかも、この計測は、試料について同じ箇所を128回行った。8dbより小さいESNRは不安定であるので除去した。
【0099】
図5は、試験例2の試験結果を示す画像図である。
図5(a)はOCT強度像を示し、
図5(b)は段落0064のESNRの値を画像濃度で表した像を示し、
図5(c)は単純平均位相像を示し、
図5(d)は本発明により補正され最適化された画像を示す。
【0100】
試料の深い部分では、
図5(c)の単純平均位相像ではコントラストが低下するが、
図5(d)はそれほどコントラストは低下しないという結果が得られた。
【0101】
図6は、
図5の(c)、(d)における点線で示す断面像であり、細線はESNRの値を示しており、太い点線はリターデーションを示している。
図5(c)、(d)において、ESNRの値(細線で示す)は、試料の奥(深い部分)にいくほど低くなっていく、要するに信号強度が小さくなりノイズが大きくなることが示されている。
【0102】
リターデーションは、単純平均による
図6(c)に示すものより、本発明を用いた
図6(d)に示すものの方が、より明確に示されている。特に、
図6(d)に示すものは、0及びπ付近に近い値が再生されている事がわかる。これによりPS−OCTの画像の解析において、本発明が有用であるということが評価された。
【0103】
また、変換を行わない処理(
図5(c)参照)では0およびπ付近の位相データが欠落し位相変化の連続性が失われているが、本発明によって処理した結果(
図5(d)参照)では、0およびπ付近の位相も再現されていて深さ方向のリターデーションの絶対値が正しく計測されていることがわかる。
【0104】
(試験例3)
本発明者らは、本発明の評価をするために、PS−OCTにより、人の眼底(in−vivo)について計測を行った。この計測は、1.0μmのプローブビームを用いて、11.0μmの眼底深さについて、30,000ライン/秒の測定速度を有するPS−OCTによって行った。
【0105】
図7は、この試験によって得られた画像を示す。
図7(a)はこれまでの単純平均処理画像であり、
図7(b)は本発明による補正処理をおこなったものである。
図7(b)に示すように、本発明によれば、位相量が0の表面付近が正確に再生されていて、より真実に近い位相量、位相差の分布が得られているものと考えられる。
【0106】
以上、本発明に係るプログラム及び該プログラムを搭載したPS−OCTシステムを実施するための形態を説明したが、本発明はこのような実施の形態に限定されるものではなく、特許請求の範囲に記載された技術的事項の範囲内でいろいろな実施の形態があることは言うまでもない。