JP2013019773A - Ｐｓ−ｏｃｔの計測データを補正するプログラム及び該プログラムを搭載したｐｓ−ｏｃｔシステム - Google Patents

Abstract

【課題】ＰＳ−ＯＣＴにおいて計測されたデータを非線形に補正し、ＰＳ−ＯＣＴの定量解析能力を高め、病変部位の病期を含めた正確な定量診断を可能とし、コンピュータ診断のための有用な手段として利用を可能とする。
【解決手段】ＰＳ−ＯＣＴ１におけるリターデーションが誤差を含みノイズとなり、その分布が正規分布でもなく、真の値の周りに対象に分布していない場合であっても、モンテカルロシミュレーションでノイズの特性を解析することによって得られた分布変換関数を用いて計測データを変換することにより、系統的な誤差を除去し、ノイズに埋もれた真の値を推測し、ＰＳ−ＯＣＴ１の画像をより鮮明に補正する。
【選択図】図１

Description

本発明は、光コヒーレンストモグラフィー（ＯＣＴ：Optical coherence tomography）による計測データを補正して高精度化するプログラム及び該プログラムを搭載したＯＣＴシステムに関する。

特に、本発明は、偏光感受型光コヒーレンストモグラフィー装置（偏光感受光画像計測装置、偏光感受型ＯＣＴ又は偏光光干渉断層計とも呼ばれているが、通常、ＰＳ−ＯＣＴ（Polarization-Sensitive OCT）と呼ばれているので、以下、「ＰＳ−ＯＣＴ」と称する。）の計測データを補正して高精度化するプログラム及び該プログラムを搭載したＰＳ−ＯＣＴシステムに関する。

従来、物体の内部情報、つまり後方散乱、反射率分布及び屈折率分布の微分構造を非破壊、高分解能で捉えるために、ＯＣＴを用いることが行われている。

医療分野等で用いられる非破壊断層計測技術の１つとして、光断層画像化法「光コヒーレンストモグラフィー」（ＯＣＴ）がある（特許文献１参照）。ＯＣＴは、光を計測プローブとして用いるため、被計測物体の反射率分布、屈折率分布、分光情報、偏光情報（複屈折率分布）等が計測できるという利点がある。

基本的なＯＣＴ４３は、マイケルソン干渉計を基本としており、その原理を図８で説明する。光源４４から射出された光は、コリメートレンズ４５で平行化された後に、ビームスプリッター４６により参照光と物体光に分割される。物体光は、物体アーム内の対物レンズ４７によって被計測物体４８に集光され、そこで散乱・反射された後に再び対物レンズ４７、ビームスプリッター４６に戻る。

一方、参照光は参照アーム内の対物レンズ４９を通過した後に参照鏡５０によって反射され、再び対物レンズ４９を通してビームスプリッター４６に戻る。このようにビームスプリッター４６に戻った物体光と参照光は、物体光とともに集光レンズ５１に入射し光検出器５２（フォトダイオード等）に集光される。

ＯＣＴの光源４４は、時間的に低コヒーレンスな光（異なった時刻に光源から出た光同士は極めて干渉しにくい光）の光源を利用する。時間的低コヒーレンス光を光源としたマイケルソン型の干渉計では、参照アームと物体アームの距離がほぼ等しいときにのみ干渉信号が現れる。この結果、参照アームと物体アームの光路長差（τ）を変化させながら、光検出器５２で干渉信号の強度を計測すると、光路長差に対する干渉信号（インターフェログラム）が得られる。

そのインターフェログラムの形状が、被計測物体４８の奥行き方向の反射率分布を示しており、１次元の軸方向走査により被計測物体４８の奥行き方向の構造を得ることができる。このように、ＯＣＴ４３では、光路長走査により、被計測物体４８の奥行き方向の構造を計測できる。

このような軸方向（Ａ方向）の走査のほかに、横方向（Ｂ方向）の機械的走査（Ｂスキャン）を加え、２次元の走査を行うことで被計測物体の２次元断面画像が得られる。この横方向の走査を行う走査装置としては、被計測物体を直接移動させる構成、物体は固定したままで対物レンズをシフトさせる構成、被計測物体も対物レンズも固定したままで、対物レンズの瞳面付近においたガルバノミラーの角度を回転させる構成等が用いられている。

以上の基本的なＯＣＴが発展したものとして、光源の波長を走査してスペクトル干渉信号を得る波長走査型ＯＣＴ（Swept Source OCT、略して「ＳＳ−ＯＣＴ」という。）と、分光器を用いてスペクトル信号を得るスペクトルドメインＯＣＴがある。後者としてフーリエドメインＯＣＴ（Fourier Domain OCT、略して「ＦＤ−ＯＣＴ」という。特許文献２参照）、及びＰＳ−ＯＣＴ（特許文献３参照）がある。

波長走査型ＯＣＴは、高速波長スキャニングレーザーにより光源の波長を変え、スペクトル信号と同期取得された光源走査信号を用いて干渉信号を最配列し、信号処理を加えることで３次元光断層画像を得るものである。なお、光源の波長を変える手段として、モノクロメーターを利用したものでも、波長走査型ＯＣＴとして利用可能である。

フーリエドメインＯＣＴは、被計測物体からの反射光の波長スペクトルを、スペクトロメーター（スペクトル分光器）で取得し、このスペクトル強度分布に対してフーリエ変換することで、実空間（ＯＣＴ信号空間）上での信号を取り出すことを特徴とするものであり、このフーリエドメインＯＣＴは、奥行き方向の走査を行う必要がなく、ｘ軸方向の走査を行うことで被計測物体の断面構造を計測可能である。

ＰＳ−ＯＣＴは、Ｂ−スキャンと同時に直線偏光したビームの偏光状態を連続変調し、試料（被検物体）のもつ偏光情報を捉え、試料のより微細な構造および屈折率の異方性を計測可能とする光コヒーレンストモグラフィー装置である。

さらに詳しくは、ＰＳ−ＯＣＴは、フーリエドメインＯＣＴと同様に、被計測物体からの反射光の波長スペクトルをスペクトル分光器で取得するものであるが、入射光及び参照光をそれぞれ１／２波長板、１／４波長板等を通して水平直線偏光、垂直直線偏光、４５°直線偏光、円偏光として、被計測物体からの反射光と参照光を重ねて１／２波長板、１／４波長板等を通して、例えば水平偏光成分だけをスペクトル分光器に入射させて干渉させ、物体光の特定偏光状態をもつ成分だけを取り出してフーリエ変換するものである。このＰＳ−ＯＣＴも、奥行き方向の走査を行う必要がない。

特開２００２−３１０８９７号公報 特開平１１−３２５８４９号公報 特開２００４−０２８９７０号公報

一般的に、ＯＣＴを生体計測へ応用することを考える場合、物体の屈折率分布の微分構造は、非破壊、高分解能で捉えることはできるものの、繊維状の構造（繊維の伸長方向等）や歯のエナメル質の相違に起因する複屈折による偏光依存性を持つ生物試料の測定においては、解像度の低下とともに、構造を捉えられない。

この点、ＰＳ−ＯＣＴは、ある特定部分からの散乱光成分とある偏光状態の参照光とを干渉させて、その干渉成分には偏光特性が強く反映され、位相遅延情報（リターデーション量）も得ることができ、その結果、奥行き方向の断面のある特定部分の偏光情報、複屈折情報を捉えることが可能となる。従って、ＰＳ−ＯＣＴは、眼科、循環器、呼吸器、消化器など広く非侵襲な組織判別を必要とする臨床分野における検査装置として有用である。

このように、ＰＳ−ＯＣＴは、組織の偏光特性を非侵襲に可視化する光トモグラフィー装置として、組織コントラストという定性観察には有用である。しかしながら、本発明者らは、ＰＳ−ＯＣＴの高精度化を目的に鋭意研究開発する過程において、ＰＳ−ＯＣＴは、病変部位を観察し、病期の正確な分類などの定量観察においては精度が必ずしも十分とは言えず、さらなる改善すべき点があるということが判った。

本発明は、上記従来のＰＳ−ＯＣＴについて得られた改善すべき点を解決することを目的とするものであり、ＰＳ−ＯＣＴにおいて計測されたデータを非線形に補正し、ＰＳ−ＯＣＴの定量解析能力を高め、病期の正確な定量診断等を可能とするプログラム及び該プログラムを搭載したＰＳ−ＯＣＴシステムを実現することを課題とする。

より具体的には、ＰＳ−ＯＣＴで計測されるリターデーション（位相遅延量、位相差の量）は誤差を含んでいて、真値の周りに分布（これをリターデーションの分布、或いは位相差の分布と称す）している。この分布は一般的に対称ではなく、必ずしも真値が中心にあるとは限らない非対称である場合も含まれる。従って、標準的な平均値をとる位相解析では正しいリターデーションを得ることができない。

本発明は、このような従来のＰＳ−ＯＣＴの問題を解決しようとするものであり、リターデーションが誤差を含みノイズとなり、その分布が正規分布でもなく、真の値の周りに対象に分布していない場合であっても、モンテカルロシミュレーションでノイズの特性を解析することによって得られた分布変換関数を用いて、計測データを変換することにより、系統的な誤差を除去し、ノイズに埋もれた真の値を推測し、ＰＳ−ＯＣＴの画像をより鮮明に補正するプログラムを実現するとともに、このプログラムを搭載したＰＳ−ＯＣＴシステムを実現することを課題とする。

本発明は上記課題を解決するために、ＰＳ−ＯＣＴで得られた計測データを処理するコンピュータに搭載されるプログラムであって、コンピュータを、ＰＳ−ＯＣＴで計測したリターデーションについての位相差の分布のデータを、非線形の変換関数によって、対称な位相差の分布のデータに変換し、平均することにより、真の位相量を求める手段として機能させることを特徴とするプログラムを提供する。

本発明は上記課題を解決するために、ＰＳ−ＯＣＴで得られた計測データを処理するコンピュータに搭載されるプログラムであって、コンピュータを、予め、０〜πまでの複数の位相量それぞれについての分布から成る位相差の分布のデータのセットを、異なるＥＳＮＲ毎に複数のセット作成し記憶する手段、ＰＳ−ＯＣＴで計測した位相差の分布のデータを対称な位相差の分布のデータに変換する非線形の変換関数の係数を決める手段、ＰＳ−ＯＣＴで計測した位相差の分布のデータからＥＳＮＲを見積もり、該見積もられたＥＳＮＲに基づき、前記異なるＥＳＮＲ毎に作成した位相差の分布の複数のセットのいずれかを選択する手段、ＰＳ−ＯＣＴで計測した位相差の分布のデータを、前記選択された位相差の分布のセットに応じて決められた係数の下で、前記変換関数によって対称な位相差の分布のデータに変換する手段、及び前記対称な位相差の分布のデータを平均することにより、真の位相量を算出する手段、として機能させることを特徴とするプログラム。

本発明は上記課題を解決するために、ＰＳ−ＯＣＴと、該ＰＳ−ＯＣＴで得られた計測データを処理するコンピュータとを備えたＰＳ−ＯＣＴシステムであって、前記コンピュータは、入力装置、出力装置、ＣＰＵ及び記憶装置を備え、ＰＳ−ＯＣＴで計測した位相差の分布のデータを、非線形の変換関数によって、対称な位相差の分布のデータに変換し、平均することにより、真の位相量を求める手段として機能することを特徴とするＰＳ−ＯＣＴシステムを提供する。

本発明は上記課題を解決するために、ＰＳ−ＯＣＴと、該ＰＳ−ＯＣＴで得られた計測データを処理するコンピュータとを備えたＰＳ−ＯＣＴシステムであって、前記コンピュータは、入力装置、出力装置、ＣＰＵ及び記憶装置を備えており、予め、０〜πまでの複数の位相量それぞれについての分布から成る複数の位相差の分布のセットを、異なるＥＳＮＲ毎に複数のセット作成し記憶装置に記憶させる手段、ＰＳ−ＯＣＴで計測した位相差の分布のデータを対称な位相差の分布のデータに変換する非線形の変換関数の係数を決める手段、ＰＳ−ＯＣＴで計測した位相差の分布のデータからＥＳＮＲを見積もり、該見積もられたＥＳＮＲに基づき、前記異なるＥＳＮＲ毎に作成した位相差の分布の複数のセットのいずれかを選択する手段、ＰＳ−ＯＣＴで計測した位相差の分布のデータを対称な位相差の分布を、前記選択された位相差の分布のセットに応じて決められた係数の下で、前記変換関数によって対称な位相差の分布のデータに変換する手段、前記対称な位相差の分布のデータを平均することにより、真の位相量を算出する手段、として機能することを特徴とするＰＳ−ＯＣＴシステムを提供する。

前記非線形の変換関数は、後記する式（２）である。

前記変換関数の係数を決める手段は、前記非線形の変換関数について、０〜πまで、π／ｍきざみで変え、次のｍ＋１の連立方程式である後記する式（３）を作成し、該連立方程式を解いて、変換関数の係数ｂｉのテーブルを作る構成である。但し、ｍは、πをきざむ回数である。

前記連立方程式を解いて前記変換関数の係数を決める手段においては、後記する式（４）の行列DMの疑似逆行列DM⁺を数値的に求めて、後記する式（５）で計算し、最適の変換関数の係数ｂｉを一義的に決める構成としてもよい。

前記変換関数の係数を決める手段は、変分法を用い、後記する式（６）で示される二乗誤差が最小となるように最適の変換関数の係数ｂｉを一義的に決める構成としてもよい。

本発明によれば、ＰＳ−ＯＣＴにおいて計測されたデータを非線形に補正し、ＰＳ−ＯＣＴの定量解析能力を高めることができるので、従来のＰＳ−ＯＣＴは、病変部位の形態観察等の用途に限られていたが、本発明のＰＳ−ＯＣＴシステムは、病変部位の病期を含め、正確な定量診断を可能とし、コンピュータ診断のための有用な手段として利用することが可能となる。

本発明に係るＰＳ−ＯＣＴシステムの全体構成を示す図である。 本発明に係るＰＳ−ＯＣＴシステムの画像処理を行うコンピュータを示す図である。 リターデーション（位相差、位相遅延量、）の分布図（位相差の分布図）であり、右欄の４つの図は、ＰＳ−ＯＣＴで計測した位相差の分布図であり、左欄はシミュレーションで作成した位相差の分布図である。 本発明における試験例１におけるシミュレーションと計測値を比較した比較試験の結果を示すものである。 本発明における試験例２の試験結果を示す画像である。 図５の（ｃ）、（ｄ）における点線で示す断面像である。 本発明における試験例３の試験結果を示す画像である。 従来のＯＣＴを説明する図である。

本発明に係るＰＳ−ＯＣＴの計測データを補正するプログラム及び該プログラムを搭載したＰＳ−ＯＣＴシステムを実施するための形態を図面を参照して、以下説明する。

本発明に係るＰＳ−ＯＣＴシステムは、ＰＳ−ＯＣＴと、ＰＳ−ＯＣＴで得られた画像データを処理する画像処理装置とを備えている。画像処理装置は、通常のコンピュータが使用され、本発明に係るプログラムは、ＰＳ−ＯＣＴで得られた画像データを補正する手段として機能させるものである。

（ＰＳ−ＯＣＴ）
ＰＳ−ＯＣＴは、すでに特許第４３４４８２９号公報等で知られているが、本発明の前提となる技術であるから、その概要を説明する。

本発明に係る偏光感受光画像計測装置は、Ｂスキャンと同時に（同期して）光源からの偏光ビーム（偏光子により直線的に偏光されたビーム）をＥＯ変調器（偏光変調器、電気光学変調器）によって連続的に変調し、この連続的に偏光を変調した偏光ビームを分けて、一方を入射ビームとして走査して試料に照射し、その反射光（物体光）を得ると共に、他方を参照光として、両者のスペクトル干渉によりＯＣＴ計測を行うものである。

そして、このスペクトル干渉成分のうち、垂直偏光成分（Ｈ）と水平偏光成分（Ｖ）を同時に２つの光検出器で測定することにより、試料の偏光特性を表すジョーンズベクトルを得る（Ｈ画像とＶ画像）構成を特徴とするものである。

図１は、本発明に係る偏光感受光画像計測装置の光学系の全体構成を示す図である。図１に示す偏光感受光画像計測装置１は、光源２、偏光子３、ＥＯ変調器４、ファイバーカプラー（光カプラー）５、参照アーム６、試料アーム７、分光器８等の光学要素を備えている。この偏光感受光画像計測装置１の光学系は、光学要素が互いにファイバー９で結合されているが、ファイバーで結合されていないタイプの構造（フリースペース型）であってもよい。

光源２は、広帯域スペクトルを有するスーパールミネッセントダイオード（ＳＬＤ：Super Luminessent Diode）を使用する。なお、光源２は、パルスレーザでもよい。光源２には、コリメートレンズ１１、光源２からの光を直線偏光にする偏光子３、進相軸を４５°の方向にセットされたＥＯ変調器（偏光変調器、電気光学変調器）４、集光レンズ１３及びファイバーカプラー５が、順次、接続されている。

ＥＯ変調器４は、進相軸を４５°の方向に固定して、該ＥＯ変調器４にかける電圧を正弦的に変調することで、進相軸とそれに直交する遅相軸との間の位相差（リターデーション）を連続的に変えるもので、これにより、光源２から出て偏光子３で（縦）直線偏光となった光がＥＯ変調器４に入射すると、上記変調の周期で、直線偏光→楕円偏光→直線偏光………などのように変調される。ＥＯ変調器４は、市販されているＥＯ変調器を使用すればよい。

ファイバーカプラー５には分岐するファイバー９を介して、参照アーム６と試料アーム７が接続されている。参照アーム６には、偏波コントローラ（polarization controller）１０、コリメートレンズ１１、偏光子１２、集光レンズ１３及び参照鏡（固定鏡）１４が、順次、設けられている。参照アーム６の偏光子１２は、上記のとおり偏光状態を変調しても参照アーム６から戻ってくる光の強度が変化しないような方向を選択するために用いている。この偏光子１２の方向（直線偏光の偏光方向）の調整は偏波コントローラ１０とセットで行う。

試料アーム７では、偏波コントローラ１５、コリメートレンズ１１、固定鏡２４、ガルバノ鏡１６、集光レンズ１３が、順次、設けられ、ファイバーカプラー５からの入射ビームが２軸のガルバノ鏡１６により走査されて試料１７に照射される。試料１７からの反射光は物体光として再びファイバーカプラー５に戻り、参照光と重畳されて干渉ビームとして分光器８に送られる。

分光器８は、順次接続される偏波コントローラ１８、コリメートレンズ１１、（偏光感受型体積位相ホログラフィック）回折格子１９、フーリエ変換レンズ２０、偏光ビームスプリッター２１及び２つの光検出器２２、２３を備えている。この実施の形態では、光検出器２２、２３として、ラインＣＣＤカメラ（１次元ＣＣＤカメラ）を利用する。ファイバーカプラー５から送られてくる干渉ビームは、コリメートレンズ１１でコリメートされ、回折格子１９によって干渉スペクトルに分光される。

回折格子１９で分光された干渉スペクトルビームは、フーリエ変換レンズ２０でフーリエ変換され偏光ビームスプリッター２１で水平及び垂直成分に分けられ、それぞれ２つラインＣＣＤカメラ（光検出器）２２、２３で検出される。この２つラインＣＣＤカメラ２２、２３は、水平および垂直偏光信号両方の位相情報を検知するために使われるので、２つのラインＣＣＤカメラ２２、２３は同一の分光器の形成に寄与するものでなくてはならない。

なお、光源２、参照アーム６、試料アーム７及び分光器８には、それぞれ偏波コントローラ１０、１５、１８が設けられているが、これらは、光源２から参照アーム６、試料アーム７、分光器８に送られるそれぞれのビームの初期偏光状態を調整して、ＥＯ変調器４で連続的に変調された偏光状態が、参照光と物体光においても互いに一定の振幅と一定の相対偏光状態の関係が維持され、さらにファイバーカプラー５に接続された分光器８において一定の振幅と一定の相対偏光状態を保たれるようにコントロールする。

また、２つラインＣＣＤカメラ２２、２３を含む分光器８を校正するときはＥＯ変調器４は止める。参照光をブロックし、スライドガラスと反射鏡を試料アーム７におく。この配置は水平および垂直偏光成分のピークの位置が同じであることを保証する。そして、スライドガラスの後ろの面と反射鏡からのＯＣＴ信号は２つの分光器８で検知される。ＯＣＴ信号のピークの位相差はモニターされる。

この位相差はすべての光軸方向の深さでゼロであるべきである。次に、信号は２つラインＣＣＤカメラ２２、２３を含む分光器８で複素スペクトルを得るために、ウィンドウされ逆フーリエ変換される。この位相差はすべての周波数でゼロであるべきなので、これらの値をモニターすることによって２つラインＣＣＤカメラ２２、２３の物理的な位置は位相差が最小になるようにアライメントされる。

ＰＳ−ＯＣＴの原理：
本発明の特徴は、次のとおりである。光源２からの光を直線偏光し、この直線偏光されたビームをＥＯ変調器４により連続的に偏光状態の変調を行う。即ち、ＥＯ変調器４は、進相軸を４５°の方向に固定して、該ＥＯ変調器４にかける電圧を正弦的に変調することで、進相軸とそれに直交する遅相軸との間の位相差（偏光角：リターデーション）を連続的に変えるもので、これにより、光源２から出て直線偏光子で（縦）直線偏光となった光がＥＯ変調器４に入射すると、上記変調の周期で、直線偏光→楕円偏光→直線偏光………などのように変調される。

そして、直線偏光された偏光ビームをＥＯ変調器４により連続的に偏光状態の変調を行うと同時に、Ｂスキャンを同期して行う。即ち、１回のＢスキャンの間に、ＥＯ変調器４による偏光の連続的な変調を複数周期行う。ここで、１周期とは、偏光角（リターデーション）φが０〜２πと変化する期間である。要するに、この１周期の間に、偏光子からの光の偏光が、直線偏光（垂直偏光）→楕円偏光→直線偏光（水平偏光）………などのように連続的に変調する。

このように偏光ビームの偏光を連続的に変調しながら、試料アーム７では、入射ビームをガルバノ鏡１６により試料１７に走査してＢスキャンを行い、分光器８において、その反射光である物体光と参照光の干渉スペクトルについて、その水平偏光成分および垂直偏光成分を２つのラインＣＣＤカメラ２２、２３で検出する。これにより、１回のＢスキャンによって、それぞれ水平偏光成分及び垂直偏光成分に対応する２枚のＡ−Ｂスキャン画像が得られる。

上記のとおり、１回のＢスキャンの間に、偏光ビームの偏光の連続的な変調を複数周期行うが、各周期（１周期）の連続的な変調の間に２つのラインＣＣＤカメラ２２、２３で検出した水平偏光成分および垂直偏光成分それぞれの偏光情報が１画素分の偏光情報となる。１周期の連続的な変調の間に２つのラインＣＣＤカメラ２２、２３で偏光情報を検出タイミング信号に同期して行い、１周期に検出回数（取込回数）を、４回、８回等、適宜決めればよい。

このようにして１回のＢスキャンの間に得た２枚のＡ−Ｂスキャン画像のデータを、Ｂスキャン方向に１次元フーリエ変換を行う。すると、０次、１次、−１次のピークが出る。ここで、０次のピークをそれぞれ抽出し、そのデータのみを用いて逆フーリエ変換すると、Ｈ０、Ｖ０画像が得られる。同様に、１次のピークをそれぞれ抽出し、そのデータのみを用いて逆フーリエ変換すると、Ｈ１、Ｖ１画像が得られる。

Ｈ０、Ｈ１画像から、試料１７の偏光特性を示すジョーンズマトリックスの成分のうち、Ｊ（１，１）およびＪ（１，２）を求める事ができる。そして、Ｖ０、Ｖ１画像から、試料１７の偏光特性を示すジョーンズマトリックスの成分のうち、Ｊ（２，１）およびＪ（２，２）を求める事ができる。

なお、ここでは、ジョーンズマトリックスの成分Ｊ（１，１）、Ｊ（１，２）、Ｊ（２，１）、Ｊ（２，２）を特許文献１の記載に基づいて記載したが、それぞれ、本発明の後記する式（１）のジョーンズマトリックスの成分Ｉ_０Ｈ（ｚ）、Ｉ_１Ｈ（ｚ）、Ｉ_０Ｖ（ｚ）、Ｉ_０Ｖ（ｚ）に相当する。

このようにして、１回のＢスキャンにおいて４つの偏光特性を含む情報が得られる。そして、この４つの情報をそれぞれ、通常のＦＤ−ＯＣＴと同様にＡスキャン方向にフーリエ変換すると、１次のピークが試料１７の深さ方向の情報を有し、しかもそれぞれ偏光特性に応じた４枚のＡ−Ｂ画像の計測データが得られる。

（本発明のプログラム及びＰＳ−ＯＣＴシステム）
上記構成のＰＳ−ＯＣＴで得られた計測データは、画像処理装置として使用するコンピュータ３０に入力される。このコンピュータ３０は通常のコンピュータであり、図２に示すように、入力部３１、出力部３２、ＣＰＵ３３、記憶装置３４及びデータバス３５を備えている。

本発明に係るＰＳ−ＯＣＴシステムに搭載されるプログラムは、コンピュータ３０の記憶装置３４に記憶されるプログラムであって、コンピュータ３０に入力されたＰＳ−ＯＣＴで得られた画像の計測データを補正し、ＰＳ−ＯＣＴの画像をより鮮明に補正する手段として、コンピュータ３０を機能させるプログラムである。このプログラムを搭載することにより、本発明に係るＰＳ−ＯＣＴシステムは、画像をより鮮明に補正する手段を備えたシステムとなる。

即ち、従来はＰＳ−ＯＣＴで計測された生データ (生の偏光位相遅延データ)から得られた画像を、直接、診断その他の各種の用途に用いていた。本発明のプログラム及びＰＳ−ＯＣＴシステムでは、ＰＳ−ＯＣＴで計測された生データを、非線形関数で補正して高精度な生体複屈折計測を可能とするものであり、この非線形関数は、モンテカルロ法（Monte-Carlo analysis）で、予め本発明のプログラムを用い、コンピュータ３０においてシミュレーションで設計される。

以下、本発明のプログラム及び該プログラム搭載したＰＳ−ＯＣＴシステムを詳細に説明する。

図３は、リターデーション（位相差、位相遅延量、）の分布図（位相差の分布図）である。横軸は位相であり、縦軸はリターデーションの生じる頻度である。図３の右欄は、ＰＳ−ＯＣＴで計測したリターデーションに基づく位相差の分布であり、図３の左欄は、シミュレーションで作成した位相差の分布である。

位相差の分布は、横軸をπ／ｍずつ刻んだｍ＋１個のデータ（ｍ＋１個の位相差の分布）を１セットとして、ＥＳＮＲを変えて複数のセット作成する。例えば、ｍ＝６０として、横軸をπ／６０ずつ刻んだ６１個のデータを１セットとして、ＥＳＮＲを変えて複数のセット作成する。

図３では、２つの異なるＥＳＮＲのセットについて、それぞれπ／２及びπについての位相差の分布を示す図である。

ところで、ＰＳ−ＯＣＴによるリターデーションのノイズモデルは、ジョーンズマトリックスのそれぞれの要素に、加算的なノイズを設定して、次の式（１）で示される。ここで、Ｓは真のリターデーションの値であり、Ｎはノイズであり、Ｉは測定値である。添え字の０、１は、計測チャンネル（図１の光検出器２２、２３参照）である。前記段落００５２に説明したジョーンズマトリックスの成分は、真の値ＳにノイズＮが加わった測定値Ｉに対応する。ここでｚは深さ方向の座標で、ジョーンズマトリックスのそれぞれの要素は表面（ｚ＝０）から深さｚまでのリターデーション等の積分量に相当する。

計測データのノイズを表す評価値として、通常は、Ｓ／Ｎ比（ＳＮ比、ＳＮＲ等とも呼ばれている）を使用するが、本発明では、より効果的にノイズを表す評価値としてＥＳＮＲを使用する。本明細書では、ＥＳＮＲはγとして表し、次の式で定義される。

１／γ＝１／４（１／ＳＮＲ_Ｓ，０＋１／ＳＮＲ_Ｓ，１＋１／ＳＮＲ_ｒ，０＋１／ＳＮＲ_ｒ，１ ）

なお、この式中、ｓはＰＳ−ＯＣＴのサンプルアームであり、ｒは参照アームである。０と１は、回数である。

複数の独立したＳＮＲをもつものを、単一の値で示すためにＥＳＮＲを導入した。添え字の０、１は、計測チャンネル（図１の光検出器２２、２３参照）である。ＳＮＲ_Ｓ，０、 ＳＮＲ_Ｓ，１、ＳＮＲ_ｒ，０、ＳＮＲ_ｒ，１ は、それぞれ添え字の組み合わせによって、サンプルアームから０チャンネルへ寄与するＳＮＲ、サンプルアームから１チャンネルへ寄与するＳＮＲ、参照アームから０チャンネルへ寄与するＳＮＲ、参照アームから１チャンネルへ寄与するＳＮＲを表す。

図３の右欄の上から２番目の図に示すように、ＰＳ−ＯＣＴで計測されるリターデーション（位相遅延量、位相差）は誤差を含んでいて、真値（π／２）の周りに分布している。そして、ＥＳＮＲが小さくなると、図３の右欄の上から１番目の図に示すように、誤差の幅を大きくなる。

また、この位相差の分布は、図３の右欄の上から３、４番目の図に示すように、必ずしも、対称ではなく、真値が中心にあるとは限らず、非対称である場合も含まれる。

要するに、リターデーションが誤差を含みノイズとなり、位相差の分布が正規分布でもなく、真の位相量の周りに対象に分布していない場合も含まれる。このような場合は、標準的な平均値をとる位相解析では正しいリターデーションを得ることができない。

本発明に係るプログラムは、モンテカルロシミュレーションでノイズの特性を解析することによって得られた分布変換関数を用いて、ＰＳ−ＯＣＴで計測された画像の計測データを変換することにより、系統的な誤差を除去し、ノイズに埋もれた真の位相量を推測し、ＰＳ−ＯＣＴの画像をより鮮明に補正する手段としてコンピュータ３０を機能させるプログラムであり、本発明のシステムは、このプログラムを搭載し上記補正を行うＰＳ−ＯＣＴシステムである。

原理：
本発明に係るプログラムは、ＰＳ−ＯＣＴで計測された画像の計測データによる実測のリターデーションの分布（位相差の分布）を補正し、真の位相量δ_Ｔ（真のリターデーションの値）を求める手段とし、コンピュータ３０を機能させるものであるが、まず、その原理（アルゴリズム）を、以下に説明する。

本発明では、ＰＳ−ＯＣＴで計測された画像の計測データによる実測の位相量δ_Ｍから、非線形関数ｆを用いて対称な位相差の分布δ_Ｅ に変換する。要するに、対称な位相差の分布δ_Ｅ は、実測の位相量δ_Ｍ の冪関数であるδ_Ｅ＝ｆ（δ_Ｍ）として表される。この関数を本明細書では、変換関数と言い、次の式（２）で表される。これは、計測値の非対称な分布δ_Ｍを、正しい値を重みの中心とする分布δ_Ｅに変換するために行われる非線形な写像を意味している。

この式（２）中、ｂｉは変換関数の係数である。最適な変換関数の係数ｂｉ（ｉは、０、１、２、…………ｎ）は、ＥＳＮＲγの関数であるが、次のように、連立方程式を解いて決めるか、或いは最小二乗法を用いて決め、ｂｉ(γ)のテーブルを作っておく。

そのために、位相値を、０からπまでπ／ｍきざみで変え、モンテカルロ法によってδ_Ｍを求める。係数biは、式（３）に示す連立方程式を解いて決める。式（３）中、〈δ_M，π／ｍ〉は、π／ｍきざみ、例えば、π／６０きざみでシミュレートした、０からπにおけるリターデーションのシミュレーション値δ_Ｍを示す。ここで、ｍは、例えば、１、２、……６０である。なお、式（３）中の〈 〉の意味は、多数のモンテカルロシミュレーションによって得られた値の平均値を示す。具体的には、式（３）を、式（４）に示すように行列で書き下し、係数ｂｉ(γ)を求める。

そして、係数を決めるには、式（４）のシミュレーションで得られたδ_{Ｍ，π／ｍ} ^ｎのｎ×ｍの行列ＤＭの疑似逆行列ＤＭ^＋を数値的に求めて、次の式（５）で計算し係数を求める。なお、式（４）中の〈 〉の意味は、多数のモンテカルロシミュレーションによって得られた値の平均値を示す。

最小二乗法を用いて係数を決める手段は、変分法を用い、次の式（６）で示す、二乗誤差Ｒ^２が最小になるように、モンテカルロ法によって、最適の変換関数の係数ｂｉ(γ)を一義的に決めることができる。ここでも、位相値を、０からπまでπ／ｍきざみで、例えばｍ＝６０として、π／６０きざみで変え、モンテカルロ法によってシミュレーション値δ_Ｍを求める。ｋは、位相をこきざみに計測した計測回を示し、０、１、２……ｍであり、例えばｍ＝６０等である。

最適の変換関数の係数ｂｉが決まれば、対称な位相差の分布δ_Ｅを表す変換関数δ_Ｅ＝ｆ（δ_Ｍ）が決まり、ＰＳ−ＯＣＴによる実測の位相量δ_Ｍを代入すれば、δ_Ｅ が求まる。

そして、δ_Ｅを次の式（７）に示すように平均する（式（７）中の〈 〉は平均の意味）ことにより、真の位相量δ_Ｔを求めることができる。平均は、一つの位相量を求めるために行った実測回数でおこなう。

構成：
本発明の原理は以上のとおりであるが、この原理に基づき、本発明のＰＳ−ＯＣＴのシステムにおいては、本発明に係るプログラムは、画像処理装置として使用するコンピュータ３０を次の手段として機能させる。

（１）位相差の分布の作成手段
位相差の分布の作成手段は、予め、図３の左欄に示すようなリターデーションについての位相差の分布をシミュレーションで複数作成する手段である。

具体的には、作成する位相差の分布は、π／ｍ、例えばｍ＝６０として、横軸をπ／６０ずつ刻んだ、６１個のデータ（６１個の位相差の分布）を１セットとして、ＥＳＮＲを変えて位相差の分布のセットを複数セット作成する。この複数セットのデータは、作成後、記憶装置３４に記憶される。

（２）変換関数の係数の決定手段
変換関数の係数の決定手段は、複数の位相差の分布のセットについて、それぞれ変換関数の係数ｂｉを決める手段である。

この変換関数の係数ｂｉは、前記したとおり、式（３）の連立方程式を解いて決めるか、或いは式（６）において最小二乗法を用いて決める。ｂｉが決められたら、そのテーブルを作成する。この変換関数の係数ｂｉのテーブルは、記憶装置３４に記憶される。

（３）位相差の分布の選択手段
位相差の分布の選択手段は、ＰＳ−ＯＣＴで計測された画像の実測の位相差の分布、位相量の補正に先立ち、その計測データ及びシステム設計等からESNRの値を見積もり、上記（１）に示す複数の位相差の分布（図３の左欄のような位相差の分布）のセットについて、どのセットを用いるかを決定する。

（４）変換関数による位相差の分布の変換手段
変換関数の算出手段は、前記記憶装置３４に記憶された変換関数の係数ｂｉ又はそのテーブルのうち、選択手段で選択された位相差の分布のセットに対応する係数ｂｉ又はそのテーブルを用いて、式（２）のｂｉに代入し、ＰＳ−ＯＣＴで計測した実測のリターデーションの位相量δ_Ｍから、対称な位相差の分布を表す変換関数δ_Ｅ＝ｆ（δ_Ｍ）を求める手段である。

（５）真の位相量算出手段
真の位相量算出手段は、上記のとおり求めた対称な位相差の分布を表すδ_Ｅ＝ｆ（δ_Ｍ）から、前記式（７）に基づき平均することにより、真の位相量δ_Ｔを算出する手段である。

（試験例１）
本発明者らは、本発明のプログラム及び該プログラム搭載したＰＳ−ＯＣＴシステムを評価するために、シミュレーション値と実測値を比較する比較試験を行った。この比較試験では、試料（計測対象物）は、ガラス板、１／８波長板及び１／４波長板であり、それぞれに生じるリターデーションの真値は、０、π／２及びπである。

なお、シミュレーション値は、前記したとおり、予めｂｉの係数をモンテカルロシミュレーションで決めておき、これによって得られた係数を用いてこのシミュレーションデータを変換したものである。実測データを上記係数で変換したものが実測値である。モンテカルロシミュレーションにはランダムノイズが入っていることと、ｂｉを有限の項数で打ち切っているので、シミュレーションといえども完全には一致はしない。

この試験では、ＰＳ−ＯＣＴによる計測のために、１．３μｍのプローブを使用した。このＰＳ−ＯＣＴは、空気中での深度解析は８．３μｍであり、３０，０００ライン／秒の測定速度を有する。計測は、試料の前方に光強度を変えることの可能な非偏光フィルタを置くことによって、ＥＳＮＲを変えて行った。

図４は、シミュレーションと試験を比較した比較試験の結果を示す。点線は、真のリターデーションの値（０、π／２及びπ）を示す。実線（ｂ）、（ｃ）、（ｅ）は、本発明による補正を行わずに、リターデーションの実測値を単に平均のみをとったものである。実線（ａ）、（ｄ）、（ｆ）は、本発明によるシミュレーションの結果であり、係数を作るのに用いたシミュレーション数は６５５３６個であり、これは、係数を求めるためのデータでもあり、結果を得るために用いた位相データでもある。

実線（ｂ）、（ｃ）、（ｅ）及び（ａ）、（ｄ）、（ｆ）の近辺に位置する●■＋のマークは、それぞれ試験によって得られたリターデーションの実測値及び該実測値を本発明により補正した試験補正値あって、ここでは６４個の測定データを用いている。シミュレーションと試験補正値は、絶対値及びＥＳＮＲ依存性等において、ほぼ一致している。

（試験例２）
本発明者らは、本発明の評価をするために、ＰＳ−ＯＣＴの評価に通常使用されている鶏の胸肉（ｅｘ−ｖｉｖｏ）を試料として計測する比較試験を行った。

この試験では、試料の１０ｍｍの部分について、５１２本のＡラインから成るＢスキャンを行った。要するに、Ａ方向（深さ方向）の計測で得られるＡラインを、Ｂスキャン（Ａ方向に直交するＢ方向のスキャン）によって、５１２本得る。しかも、この計測は、試料について同じ箇所を１２８回行った。８ｄｂより小さいＥＳＮＲは不安定であるので除去した。

図５は、試験例２の試験結果を示す画像図である。図５（ａ）はＯＣＴ強度像を示し、図５（ｂ）は段落００６４のＥＳＮＲの値を画像濃度で表した像を示し、図５（ｃ）は単純平均位相像を示し、図５（ｄ）は本発明により補正され最適化された画像を示す。

試料の深い部分では、図５（ｃ）の単純平均位相像ではコントラストが低下するが、図５（ｄ）はそれほどコントラストは低下しないという結果が得られた。

図６は、図５の（ｃ）、（ｄ）における点線で示す断面像であり、細線はESNRの値を示しており、太い点線はリターデーションを示している。図５（ｃ）、（ｄ）において、ESNRの値（細線で示す）は、試料の奥（深い部分）にいくほど低くなっていく、要するに信号強度が小さくなりノイズが大きくなることが示されている。

リターデーションは、単純平均による図６（ｃ）に示すものより、本発明を用いた図６（ｄ）に示すものの方が、より明確に示されている。特に、図６（ｄ）に示すものは、０及びπ付近に近い値が再生されている事がわかる。これによりＰＳ−ＯＣＴの画像の解析において、本発明が有用であるということが評価された。

また、変換を行わない処理（図５（ｃ）参照）では０およびπ付近の位相データが欠落し位相変化の連続性が失われているが、本発明によって処理した結果（図５（ｄ）参照）では、０およびπ付近の位相も再現されていて深さ方向のリターデーションの絶対値が正しく計測されていることがわかる。

（試験例３）
本発明者らは、本発明の評価をするために、ＰＳ−ＯＣＴにより、人の眼底（ｉｎ−ｖｉｖｏ）について計測を行った。この計測は、１．０μｍのプローブビームを用いて、１１．０μｍの眼底深さについて、３０，０００ライン／秒の測定速度を有するＰＳ−ＯＣＴによって行った。

図７は、この試験によって得られた画像を示す。図７（ａ）はこれまでの単純平均処理画像であり、図７（ｂ）は本発明による補正処理をおこなったものである。図７（ｂ）に示すように、本発明によれば、位相量が０の表面付近が正確に再生されていて、より真実に近い位相量、位相差の分布が得られているものと考えられる。

以上、本発明に係るプログラム及び該プログラムを搭載したＰＳ−ＯＣＴシステムを実施するための形態を説明したが、本発明はこのような実施の形態に限定されるものではなく、特許請求の範囲に記載された技術的事項の範囲内でいろいろな実施の形態があることは言うまでもない。

従来は、ＰＳ−ＯＣＴは、病変部位の形態観察用途に限られていたが、本発明に係るプログラム及び該プログラムを搭載したＰＳ−ＯＣＴシステムは、上記のような構成であるから、眼科、循環器、呼吸器、消化器など広く被侵襲な組織判別を必要とする臨床分野において、病変部位の病期を含めた正確な定量診断を可能とし、コンピュータ診断のための有用な手段として利用することが可能である。

例えば、前眼部ＰＳ−ＯＣＴ装置に組み込んで利用することで、緑内障手術後の手術痕の非接触・非侵襲な観察が可能になり、結果、眼科治療の方針が劇的に効率化が可能となる。

１ ＰＳ−ＯＣＴ（偏光感受光画像計測装置）
２ 光源
３、１２ 偏光子
４ ＥＯ変調器（偏光変調器、電気光学変調器）
５ ファイバーカプラー（光カプラー）
６ 参照アーム
７ 試料アーム
８ 分光器
９ ファイバー
１０、１５、１８ 偏波コントローラ（polarization controller）
１１ コリメートレンズ
１３ 集光レンズ
１４ 参照鏡（固定鏡）
１６ ガルバノ鏡
１７ 試料
１９ 回折格子
２０ フーリエ変換レンズ
２１ 偏光ビームスプリッター
２２、２３ 光検出器（ラインＣＣＤカメラ）
２４ 固定鏡
３０ コンピュータ（画像処理装置）
３１ 入力部
３２ 出力部
３３ ＣＰＵ
３４ 記憶装置
３５ データバス

Claims (12)

1. ＰＳ−ＯＣＴで得られた計測データを処理するコンピュータに搭載されるプログラムであって、
   コンピュータを、
   ＰＳ−ＯＣＴで計測したリターデーションについての位相差の分布のデータを、非線形の変換関数によって、対称な位相差の分布のデータに変換し、平均することにより、真の位相量を求める手段として機能させることを特徴とするプログラム。
2. ＰＳ−ＯＣＴで得られた計測データを処理するコンピュータに搭載されるプログラムであって、
   コンピュータを、
   予め、０〜πまでの複数の位相量それぞれについての分布から成る位相差の分布のデータのセットを、異なるＥＳＮＲ毎に複数のセット作成し記憶する手段、
   ＰＳ−ＯＣＴで計測した位相差の分布のデータを対称な位相差の分布のデータに変換する非線形の変換関数の係数を決める手段、
   ＰＳ−ＯＣＴで計測した位相差の分布のデータからＥＳＮＲを見積もり、該見積もられたＥＳＮＲに基づき、前記異なるＥＳＮＲ毎に作成した位相差の分布の複数のセットのいずれかを選択する手段、
   ＰＳ−ＯＣＴで計測した位相差の分布のデータを、前記選択された位相差の分布のセットに応じて決められた係数の下で、前記変換関数によって対称な位相差の分布のデータに変換する手段、及び
   前記対称な位相差の分布のデータを平均することにより、真の位相量を算出する手段、
   として機能させることを特徴とするプログラム。
3. 前記非線形の変換関数は、次の式（２）であることを特徴とする請求項２記載のプログラム。
   但し、式（２）中、δＥは対称な位相差の分布であり、δＭはＰＳ−ＯＣＴで計測した位相量であり、ｂｉ（ｉ＝０、１、２、３、………ｎ）は変換関数の係数であり、γはＥＳＮＲである。
   
4. 前記変換関数の係数を決める手段は、前記非線形の変換関数について、０〜πまで、π／ｍきざみで変え、次のｍ＋１の連立方程式である式（３）を作成し、該連立方程式を解いて、変換関数の係数ｂｉのテーブルを作る構成であることを特徴とする請求項３記載のプログラム。
   但し、式（３）中のｍは、πをきざむ回数である。
   
5. 前記連立方程式を解いて前記変換関数の係数を決める手段においては、式（４）の行列DMの疑似逆行列DM⁺を数値的に求めて、式（５）で計算し、最適の変換関数の係数ｂｉを一義的に決める構成であることを特徴とする請求項４記載のプログラム。
   但し、式（４）中のｍは、πをきざむ回数である。
   
6. 前記変換関数の係数を決める手段は、変分法を用い、次の式（６）で示される二乗誤差が最小となるように最適の変換関数の係数ｂｉを一義的に決める構成であることを特徴とする請求項３記載のプログラム。
   但し、式（６）中のｍは、πをきざむ回数である。
   
7. ＰＳ−ＯＣＴと、該ＰＳ−ＯＣＴで得られた計測データを処理するコンピュータとを備えたＰＳ−ＯＣＴシステムであって、
   前記コンピュータは、入力装置、出力装置、ＣＰＵ及び記憶装置を備え、ＰＳ−ＯＣＴで計測した位相差の分布のデータを、非線形の変換関数によって、対称な位相差の分布のデータに変換し、平均することにより、真の位相量を求める手段として機能することを特徴とするＰＳ−ＯＣＴシステム。
8. ＰＳ−ＯＣＴと、該ＰＳ−ＯＣＴで得られた計測データを処理するコンピュータとを備えたＰＳ−ＯＣＴシステムであって、
   前記コンピュータは、
   入力装置、出力装置、ＣＰＵ及び記憶装置を備えており、
   予め、０〜πまでの複数の位相量それぞれについての分布から成る複数の位相差の分布のセットを、異なるＥＳＮＲ毎に複数のセット作成し記憶装置に記憶させる手段、
   ＰＳ−ＯＣＴで計測した位相差の分布のデータを対称な位相差の分布のデータに変換する非線形の変換関数の係数を決める手段、
   ＰＳ−ＯＣＴで計測した位相差の分布のデータからＥＳＮＲを見積もり、該見積もられたＥＳＮＲに基づき、前記異なるＥＳＮＲ毎に作成した位相差の分布の複数のセットのいずれかを選択する手段、
   ＰＳ−ＯＣＴで計測した位相差の分布のデータを対称な位相差の分布を、前記選択された位相差の分布のセットに応じて決められた係数の下で、前記変換関数によって対称な位相差の分布のデータに変換する手段、及び
   前記対称な位相差の分布のデータを平均することにより、真の位相量を算出する手段、
   として機能することを特徴とするＰＳ−ＯＣＴシステム。
9. 前記非線形の変換関数は、次の式（２）であることを特徴とする請求項８記載のＰＳ−ＯＣＴシステム。
   但し、式（２）中、δＥは対称な位相差の分布であり、δＭはＰＳ−ＯＣＴで計測した位相量であり、ｂｉ（ｉ＝０、２、３、………ｎ）は変換関数の係数であり、γはＥＳＮＲである。
   
10. 前記変換関数の係数を決める手段は、前記非線形の変換関数について、０〜πまで、π／ｍきざみで変え、次のｍ＋１の連立方程式である式（３）を作成し、該連立方程式を解いて、変換関数の係数ｂｉのテーブルを作る構成であることを特徴とする請求項９記載のＰＳ−ＯＣＴシステム。
    但し、式（３）中のｍは、πをきざむ回数である。
    
11. 前記連立方程式を解いて変換関数の係数を決める手段においては、式（４）の行列DMの疑似逆行列DM⁺を数値的に求めて、式（５）で計算し、最適の変換関数の係数ｂｉを一義的に決める構成であることを特徴とする請求項１０記載のＰＳ−ＯＣＴシステム。
    但し、式（４）中のｍは、πをきざむ回数である。
    
12. 前記変換関数の係数を決める手段は、変分法を用い、次の式（６）で示される二乗誤差が最小となるように最適の変換関数の係数ｂｉを一義的に決める構成であることを特徴とする請求項９記載のＰＳ−ＯＣＴシステム。
    但し、式（６）中のｍは、πをきざむ回数である。