JP6480104B2 - 光コヒーレンストモグラフィー装置及び光コヒーレンストモグラフィーによる変位測定方法 - Google Patents

Description

本発明は、光コヒーレンストモグラフィーを用いて被測定物体の変位を測定する変位測定装置及び測定方法に関する。

光凝固治療は、いくつかの眼疾患（例えば、糖尿病、網膜血管閉塞、等）のために好適に利用される効果的なレーザ治療である。しかしながら、光凝固治療は、医原性の合併症が生じる可能性がある。

大半の重大な障害は、治療のしすぎによる過度のレーザダメージである。そのようなダメージを防止するため、レーザ照射パラメータは、高い再現性で、ダメージのある閾値より少ない条件にて設定されることが要求されている。実際、この要求は厳しい。これは、レーザ治療は、再生産性が乏しいからである。

大半のケースの場合、レーザ照射は、白化する病変の視認性に基づく術者の主観的な評価に依存する。加えて、組織の散乱・吸収特性の変化は、評価をより難しくさせる。色素（pigments）を含む眼内媒体の光学的な特性は、空間的、時間的に変化するし、また、個人差によっても変化する。レーザ照射中における組織変化において多くの不確実性があるため、なんらかのモニタリング技術が要求される。

近年、リアルタイムで非侵襲のモニタリング方法は、再生産性を改善する潜在的な能力を示している。例えば、光音響信号は、レーザによって誘発された温度上昇を測定するために用いられる。この技術は、温度コントロールされた光凝固を可能とする。

もう一つの例としては、光干渉断層計（光コヒーレンストモグラフィー：optical coherence tomography。本明細書では、略して「ＯＣＴ」とも言う。）である。ＯＣＴは、散乱の変化を測定するだけでなく、温度膨張による組織変位を、位相感度測定を用いて測定することも可能である。特に、ＯＣＴは、光音響測定よりもいくつかの利点を持つ。

第１に、ＯＣＴは、非接触での測定が可能である。非接触での測定は、ルーティンでの臨床実務において重要である。第２として、ＯＣＴは、局所的な組織変化を観察するのに必要な、空間的に分解された構造変化を可視化できる。第３として、ＯＣＴは、機械的な変位を含む、組織の温度変化を直接的にモニタリングできる。

そのような変位は、位相感度測定（例えば、非特許文献１参照）、スペックル追跡、相関係数（例えば、非特許文献２、３参照）のいずれかを用いて測定できる。つまり、ＯＣＴ測定は、温度膨張測定において重要な役割を持つ。

B.J.Vakoc,G.J.Tearney,and B.E.Bouma,"Real-time microscopic visualization of tissue response to laser thermal therapy,"Journalof Biomedical Optics 12,020501-020501(2007). A.Ahmad,S.G.Adie,E.J.Chaney,U.Sharma,and S.A.Boppart,"Cross-correlation-based image acquisition technique for manually-scanned optical coherence tomography,"Optics Express 17,8125-8136(2009). X.Liu,Y.Huang,and J.U.Kang,"Distortion-free freehand-scanning OCT implemented with real-time scanning speed variance correction,"Optics Express 20,16567-16583(2012).

しかしながら、ＯＣＴを用いた従来の測定方法では、Ａスキャン信号の軸方向に関する変位のみを測定する方法（従来の位相感度測定）、Ａスキャン信号の横方向に関する変位のみを測定する方法（従来の相関係数を用いた変位測定）しか存在しなかった。つまり、従来の方法は、断層画像中のある微小領域の変位を二次元的に測定することは困難であった。

本発明は、上記問題点を解決することを目的とし、ＯＣＴを用いた新規な変位測定装置を提供することを課題とする。

本発明は上記課題を解決するために、被測定物体の奥行き方向の軸に平行な２次元断層画像を同一位置に関して複数取得し、取得された複数の断層画像に基づいて被測定物体の変位の測定が可能な光コヒーレンストモグラフィー装置であって、前記複数の断層画像のうちの１つの断層画像を基準画像とし、前記基準画像の微小領域に対して他の断層画像の対応する微小領域を前記軸方向及び前記軸方向に直交する断層画像面内の方向にシフトすることによって、相関係数を計算する相関係数算出部と、前記相関係数算出部からの出力に基づいて、前記軸方向の変位と前記軸方向に直交する断層画像面内の方向の変位とを相関係数に関係づける相関係数の数学モデルのパラメータを算出する相関係数モデル算出部と、変位測定対象の断層画像と前記基準画像との位相シフトを求める位相シフト算出部と、前記位相シフト算出部からの出力に基づいて前記軸方向の変位を求める軸方向変位算出部と、前記軸方向変位算出部からの出力及び前記相関係数モデル算出部からの出力に基づいて前記軸方向に直交する断層画像面内の方向の変位を求める横方向変位算出部と、を有する光コヒーレンストモグラフィー装置を提供する。

前記相関係数算出部は、前記基準画像の信号対ノイズ比と前記他の断層画像の信号対ノイズ比とをそれぞれ算出し、後述する（２）式に基づいて、前記信号対ノイズ比の影響なしで信号相関係数を計算し、前記信号相関係数を前記相関係数に置き換えることが好ましい。

また、本発明は上記課題を解決するために、被測定物体の奥行き方向の軸に平行な２次元断層画像を同一位置に関して複数取得し、取得された複数の断層画像に基づいて被測定物体の変位を測定する光コヒーレンストモグラフィーによる変位測定方法であって、前記複数の断層画像のうちの１つの断層画像を基準画像とし、前記基準画像の微小領域に対して他の断層画像の対応する微小領域を前記軸方向及び前記軸方向に直交する断層画像面内の方向にシフトすることによって、相関係数を計算する相関係数算出工程と、前記相関係数算出工程で算出された前記相関係数に基づいて、前記軸方向の変位と前記軸方向に直交する断層画像面内の方向の変位とを相関係数に関係づける相関係数の数学モデルのパラメータを算出する相関係数モデル算出工程と、変位測定対象の断層画像と前記基準画像との位相シフトを求める位相シフト算出工程と、前記位相シフト算出工程で算出された前記位相シフトに基づいて前記軸方向の変位を求める軸方向変位算出工程と、前記軸方向変位算出工程で算出された前記軸方向の変位及び前記相関係数モデル算出工程で算出された前記パラメータに基づいて前記軸方向に直交する断層画像面内の方向の変位を求める横方向変位算出工程と、を有する光コヒーレンストモグラフィーによる変位測定方法を提供する。

前記相関係数算出工程において、前記基準画像の信号対ノイズ比と前記他の断層画像の信号対ノイズ比とをそれぞれ算出し、後述する（２）式に基づいて、前記信号対ノイズ比の影響なしで信号相関係数を計算し、前記信号相関係数を前記相関係数に置き換えることが好ましい。

被測定物体に向けてレーザを照射した場合の被測定物体の変位を測定する構成であることが好ましい。

前記微小領域を、前記断層画像中に二次元的に複数設定し、前記対応関係を用いて、各微小領域に関して測定された相関係数の減衰から各微小領域の変位量を求め、各微小領域での変位量の二次元的な分布を示す変位マップを取得する構成であることが好ましい。

前記複数の断層画像は、時間的に連続する複数の断層画像であって、測定された相関係数の減衰から単位時間当たりの前記微小領域の変位量を求める構成であることが好ましい。

従来の測定では、Ａスキャン信号の軸方向に関する変位のみを測定する方法（従来の位相感度測定）、Ａスキャン信号の横方向に関する変位のみを測定する方法（従来の相関係数を用いた変位測定）しか存在しなかったが、本発明に係る光コヒーレンストモグラフィーによる変位測定装置及び測定方法によれば、断層画像中のある微小領域の変位を二次元的に測定することが可能となった。

本発明に係る光コヒーレンストモグラフィーによる変位測定装置及び測定方法の処理ステップの一例を示すフローチャートである。 本発明に係る光コヒーレンストモグラフィーによる変位測定装置及び測定方法の実施例として使用するＦＤ−ＯＣＴの全体構成を示す図であり、（ａ）は実験確認のためのサンプルアームの例であり、（ｂ）は、（ａ）のレンズ、ガルバノミラー及び被測定物体の部分に関する他の態様を示し、レーザ照射光学系を設けた場合の例である。 （ａ）、（ｂ）は互いに同じ写真で、ともに組織ファントムのＢスキャン画像の代表図であり、（ａ）の丸で囲われた１、２に示す十字は、それぞれランダムに選択された位置を示し、（ｂ）の四角は関心領域を示す。 測定された相関係数及び計算上シミュレートされた相関係数が、横方向の変位の関数としてプロットされたグラフである。それらは、図３（ａ）の選択された位置で観察された。 本発明についての実験の結果を示す図であり、（ａ）は測定された横方向における変位をリニアにプロットしたものであり、上記関心領域での最頻値がプロットされている。（ｂ）は測定された横方向の変位のＲＭＳエラーがプロットされている。（ｃ）は測定された軸方向の変位のリニアなプロットである。（ｄ）は測定された軸方向の変位のＲＭＳエラーのプロットである。 （ａ）はＯＣＴのＢスキャンの強度画像、（ｂ）はＨＳＬでの変位マップであり、マップのＨｕｅと膨張は、変位ベクトルの方向と大きさによって与えられる。（ｃ）は測定された相関係数を示す画像であり、（ｄ）は測定されたドップラー位相シフトを示す画像であり、アンラッピング前のドップラー位相シフトを示す画像である。 レーザ照射による被測定物体の変位を示す図であり、（ａ）は、測定された軸方向の変位であり、（ｂ）は、横方向の変位である。 レーザ照射されたＲＰＥ複合体での横方向の平均と、軸方向における内網膜の最大変位が、ＢスキャンＩＤの関数としてプロットされている。

本発明に係る光コヒーレンストモグラフィーによる変位測定装置及び測定方法（以下、略して「本発明の変位測定装置及び測定方法」とも言う）を実施するための形態を実施例に基づいて図面を参照して、以下に説明する。

本発明者らは、本発明において、ＯＣＴを用いた二次元的な組織変位測定装置及び測定方法を提案し、レーザ照射中での温度膨張測定も可能とした。本発明の変位測定装置及び測定方法における主な特徴は、本発明の変位測定装置及び測定方法は、ドップラ位相シフトと相関係数の両方を用いて、軸方向と横（側面）方向の変位を測定する構成にある。

以下の実施例の記載では、次の順で、本発明の実施例を明確にする。
（１）まず、本発明の変位測定装置を構成する、また測定方法に使用するための光コヒーレンストモグラフィーを含むシステムについて説明する。

（２）次に、本発明の原理及び処理動作について順次説明する。即ち、本発明の変位測定装置及び変位測定方法は、その特徴的な構成は、上記のとおり、ドップラ位相シフトと相関係数の両方を用いて、軸方向と横（側面）方向の変位を測定する構成にある。

この構成は、変位測定装置を構成するシステムに含まれるコンピュータに搭載される変位測定のためのプログラムが、コンピュータを、ドップラ位相シフトと相関係数の両方を用いて、軸方向と横（側面）方向の変位を算出するように処理動作（機能）させることで具現化される。従って、変位測定のためのプログラムによる動作について、その原理（本発明の原理）及び処理動作を説明する。

（３）次に、本発明の変位測定装置及び変位測定方法の効果を確認するために、本発明者らは実験を行った。以下、その実験例について、実験装置、実験結果等について説明する。

（変位測定装置のシステム）
本発明に係る光コヒーレンストモグラフィーによる変位測定装置及び測定方法では、一例として、図２に示すＦＤ−ＯＣＴ１（ＦＤ−ＯＣＴはフーリエドメインＯＣＴの略であり、特開２００７−２９８４６１号公報参照）を使用する。

ＦＤ−ＯＣＴ１は、図２（ａ）に示すように、広帯域光源２、干渉計３、及び分光器４（スペクトロメータ）とを備えている。このＦＤ−ＯＣＴ１は、低コヒーレンス干渉の原理を用いて奥行き方向の分解能を得ているため、光源として、ＳＬＤ（スーパールミネッセントダイオード）、超短パルスレーザ等の広帯域光源２が用いられる。

広帯域光源２から出た光は、まずビームスプリッター５で物体光と参照光に分割される。このうち物体光は、レンズ６を通してガルバノミラー７で反射され被測定物体８（眼底などの生体試料）を照射し、そこで反射、散乱された後に分光器４に導かれる。

一方、参照光はレンズ９を通して参照鏡１０（平面鏡）で反射された後に物体光と並行に分光器４に導かれる。これらの二つの光は分光器４の回折格子によって同時に分光され、スペクトル領域で干渉し、結果、スペクトル干渉縞がＣＣＤによって計測される。

このスペクトル干渉縞に対して適当な信号処理を行うことで、被測定物体８のある点における深さ方向１次元の屈折率分布の微分、つまり、反射率分布を得ることが可能となる。さらに、ガルバノミラー７を駆動することによって被測定物体８上の計測点を１次元走査することにより２次元断層画像（ＦＤ−ＯＣＴ画像）を得ることができる。

通常のＯＣＴでは、２次元断層画像を得るために、深さ（光軸）方向の走査（この走査を「Ａ−スキャン」と言い、この方向を「Ａ−方向」、「Ａスキャン方向」とも言う。）と、横方向の操作（この走査を「Ｂ−スキャン」と言い、この方向を「Ｂ−方向」、「Ｂスキャン方向」とも言う。）の２次元の機械的走査が必要なのに対して、ＦＤ−ＯＣＴ１では、Ａ−スキャンは不要で一回の測定で深さ方向の後方散乱データを取得することができるから、Ｂ−スキャンの１次元の機械的走査しか必要とされない。

なお、物体光を走査するための光スキャナとしては、ガルバノミラーに限定されず、ポリゴンミラー、レゾナントスキャナ、音響光学素子（ＡＯＭ）であってもよい。

また、物体光と合成させる参照光を生成するための参照光学系としては、マイケルソンタイプであってもよいし、マッハツェンダタイプであっても良い。

なお、上記例では、参照ミラーを用いた反射光学系を用いたが、透過光学系（例えば、光ファイバー）によって構成されてもよい。透過光学系は、例えば、カップラーからの光を戻さず透過させることにより検出器へと導く。

ＦＤ−ＯＣＴ１には、追加的に、図２（ｂ）に示すように、被測定物体８に向けてレーザを照射するためのレーザ照射光学系１００が設けられていてもよい。レーザ照射光学系１００は、例えば、治療用レーザ光を照射し、被測定物体８を凝固するために用いられる。なお、レーザ照射光学系１００は、被測定物体を低侵襲にて治療可能な構成であってもよい（例えば、マイクロパルスレーザー装置、低出力パルスレーザ装置等）。

レーザ照射光学系１００には、治療用レーザ光を被測定物体８に対して走査させるための治療レーザ用光スキャナが設けられてもよい。なお、治療レーザ用光スキャナと、ＦＤ−ＯＣＴ１の光スキャナは、同期して制御されてもよい。また、同じ光スキャナが用いられても良い。

なお、レーザ照射光学系を備えるＯＣＴの詳細な構成については、例えば、特開２０１２−２１３６３４号公報、特開２０１２−１３５５５０号公報等に記載されている。

なお、レーザ照射光学系１００の光軸は、光路結合部材（例えば、ダイクロイックミラー、ハーフミラー等）２００によってＦＤ−ＯＣＴ１の光軸と同軸に配置される。

以下の説明では、被測定物体１が眼の網膜であって、レーザ照射中の網膜の変位を測定する場合を例として説明する。被測定物体１は、眼に限定されない生体であってもよいし、生体でなくてもよい。

ＦＤ−ＯＣＴ１は、装置全体を制御するコンピュータ３００のプロセッサ（ＣＰＵ）によって制御される。コンピュータ３００は、取得された画像を処理する画像処理部として用いられる。また、ＦＤ−ＯＣＴ１によって取得される画像は、コンピュータ３００の記憶部に記憶されるように構成する。

この記憶部には、ＦＤ−ＯＣＴ１を動作させるためのプログラムが格納されており、また後述する変位測定のための処理動作を行うためのプログラムが格納されている。このように、コンピュータに搭載される変位測定のためのプログラムが、コンピュータを、ドップラ位相シフトと相関係数の両方を用いて、軸方向と横（側面）方向の変位を算出する処理動作（機能）させるようにする。

（本発明の原理）
特定の画素でのＯＣＴ信号sは、下記の数１に示す式（１）のように、ランダムに分布した散乱体の複素反射率ａ（ｒ）と、点像分布関数（ＰＳＦ：point spread function）ｈ（ｒ）と、のコンボリューション（＊）によって表現される。

ここで、ｒ＝（ｘ、ｙ、ｚ）は、デカルト（Ｃａｒｔｅｓｉａｎ）座標における位置ベクトルである。（ｘ、ｙ）は、横方向の位置を示す。ｚは深さ位置を示す。

仮に、散乱現象の数が十分大きい場合、確率変数Ｓは、複素の円周標準正規分布に従う確率変数として考慮される。加えて、測定における加算的ノイズＮは、複素の円周標準正規分布に従い、ＯＣＴ信号Ｓから独立する。

それらの加算、すなわちＧ＝Ｓ＋Ｎによって与えられる測定信号Gは、また、複素の円周標準正規分布に従う確率変数である。ここで、２つの測定された信号Ｇ１とＧ２との間の母集団の複素相関係数は、次の数2に示す数式で定義される。

ここで、振幅ρは、母集団の複素相関係数を表し、位相差Φは、母集団の位相シフトを示し、添字（subcript）＊は、複素共役を示す。

測定信号Ｇ_ｉ（ｉ＝１、２）を、ＯＣＴ信号Ｓ_ｉと追加的ノイズＮ_ｉに置換することによって、ＯＣＴ信号Ｓ_１とＳ_２との間の母集団の信号相関係数は、次の数３に示す数式で定義され、下記の数４に示す式（２）のように書き換えられる。

ここで、ＳＮＲ _１ とＳＮＲ _２ は、各測定における信号対ノイズ比の期待値であり、次の数５に示す数式として定義される。

この数式によって、信号対ノイズ比の影響無しで、信号相関係数が求められる。この数式に関するさらなる詳細は、17. S. Makita, F. Jaillon, I. Jahan, and Y. Yasuno, “Noise statistics of phase-resolved optical coherence tomography imaging: single-and dual-beam-scan Doppler optical coherence tomography,”Opt. Express 22(4), 4830-4848 (2014)に記載されている。

二次元変位測定において、変位前後の信号相関係数は、次の数６に示す数式として書き換えられる。Ｃは、散乱過程と、走査システムの再現性の影響を示した定数である。ρ_ｈは、ＰＳＦ_ｓの相関係数であり、Δｘは横方向の変位であり、Δｚは軸方向の変位である。

ここで、画像面内（in-plane）と剛体（rigid-body）の変位のみが考慮される。言い換えれば、他のファクターによって引き起こされる相関係数の減少が考慮されない。他のファクターとしては、画像面外（out-of-plane）の変位Δｙ、非剛体の変位、すなわち、散乱体分布の歪み（deformation of scatterer distribution）が挙げられる。

サンプル上でのビームプロファイルをガウス関数、光源のスペクラム形状をガウス関数と仮定すると、ＰＳＦは、次の数７に示す数式として表される。

ここで、ｗは、１／ｅ^２でのガウシアンビームスポット半径であり、Δｋは、光源のガウシアンスペクトラムの１／ｅ^２での最大幅であり、ｚ_０は、干渉計におけるゼロディレイでの深さ位置である。

信号相関係数は、以前の研究（例えば、X. Liu, Y. Huang, and J. U. Kang, “Distortion-free freehand-scanning OCT implemented with real-time scanning speed variance correction,” Optics Express 20, 16567-16583 (2012)）に類似した手法で、次の数８に示す式（３）として書き換えられる。

信号相関係数は、変位の関数として減少し、それは、ガウシアン関数に従う。つまり、相関係数は、変位量に対してガウス関数的に減衰する。式（３）を横方向の変位Δｘについて解くと、次の数９に示す式（４）として表される。

よって、式（３）からＯＣＴ画像の平面内での移動量は、相関係数の減衰から理論的に推定できる。

軸方向の変位は、複素相関係数の位相ΔΦを得ることによって取得される。このドップラシフトは、アンラップされ、以下の数１０の式（５）を用いた軸方向の変位に変換される。ここで、ｎはサンプルの屈折率、λは光源の中心波長、Φunwrappedは、アンラップされた位相シフトである。

（処理動作）
＜処理ステップ＞
ＯＣＴ Ａラインは、スペクトル干渉縞のフーリエ変換を用いることによって取得される。スペクトルインターフェログラムは、波数空間においてリサンプルされ、再形成（reshaped）されてガウシアンプロファイルを持つ。

その分散は、計算上補正される。複数のＯＣＴ Ｂスキャンは、下記数１１の式で表されるように、サンプルの同一位置で繰り返しプローブビームを走査することによって取得される。ここで、Ｎは、Ｂスキャンの数である。

本発明の変位測定装置及び測定方法によって、ＯＣＴＢスキャン間において、空間的に分解された二次元変位を測定する。変位測定のための処理ステップは、図１のフローチャートにおいて要約される。

前処理において、分解能と、式（３）内の定数パラメータ(w;Δk;C)と、が求められる。分解能パラメータは、Ｂスキャンにおいて不変にセットされうる。しかし実際には、分解能パラメータは、収差による空間的な位置、散乱過程、ノイズの寄与、によって変化する。そこで、変位測定前の全画素に関するパラメータを求めてもよい。

後処理では、相関係数、ドップラー位相シフト、求められたパラメータを用いて、全画素に関する横方向と軸方向の変位が計算されてもよい。各画素での変位は、Ｋｘ×Ｋｚの画素ブロック（微小領域）の中で近接する画素を用いて計算される。ＫｘとＫｚは、それぞれ横方向と軸方向におけるＫｅｒｎｅｌサイズである。

＜前処理＞
前処理において、以下のステップをとることによってデータ処理される。最初のステップでは、相関係数ρsimulationは、参照Ｂスキャンと、デジタル処理によってシフトされたＢスキャンと、の間で計算される。

この目的のため、変位がない２つの連続するＢスキャンｇ_ｒ(ｘ，ｙ)とｇ_ｒ+１(ｘ，ｙ)が選択される。一方は、参照Ｂスキャンとして用いられ、他方は、対象Ｂスキャンとして用いられる（図１（ア）参照）。

対象Ｂスキャンは、デジタル処理によってシフトされたＢスキャンを生成するため、ゼロパッディング（zero-padding）を用いてリサンプルされる。ここで、対象Ｂスキャンは、横方向にシフトされる。シフトされたＢスキャンは、次の数１２の数式として示される。ｄは、一定の変位量である。

ｇ_ｒ(ｘ，ｙ)とｇ_ｒ+１(ｘ＋Δｘ_ｌ，ｙ)との間の相関は、次の数１３に示す式（６）を用いて計算される（図１（イ）参照）。ρ_ｃは、計算された相関係数であり、Ｒは参照Ｂスキャンの識別情報、Ｔは対象Ｂスキャンの識別情報である。

ここで、例えば、７×７ピクセルのＫｅｒｎｅｌサイズが用いられる。ＳＮＲ_ｎ(ｘ，ｙ)は、次の数１４に示す式（７）によって与えられる。Ｎは、ノイズフロアである。

この計算は、式（２）に基づくので、ノイズに対してよりロバストである。実際、ノイズフロアは、測定前に記録される。しかし、ノイズ成分は、確率論的である。計算されたＳＮＲは、真のＳＮＲといつも対応しない。同様なことは、計算された相関係数でも起こる。例えば、いくつかのケースにおいて相関係数が１より大きくなることがある。

第２のステップでは、最小二乗フィッティングは、式（３）における横方向の分解能パラメータを見つけるために、全ての画素に関して行われる。

第１のステップにおいて取得され、計算上シミュレートされた相関係数ρ_{ｓｉｍｕｌａｔｉｏｎ}(Δｘ_ｌ、０)は、全ての画素に関して、数学モデルρ_ｓ(Δｘ_ｌ，０;Ｃ，ｗ)にフィットされる。次の数１５に示されるフィットされたパラメータは、下記の数１６に示される等式（８）によって与えられる。

同時に、上記処理ステップ１と２は、全ての画素に関して、軸方向の分解能パラメータΔｋを見つけるために行われる。

最後のステップでは、定数パラメータＣ´は、選択された連続するＢスキャンｇ_ｒ(ｘ，ｙ)とｇ_ｒ+１(ｘ，ｙ)との間の相関によって決定される。結果として求められた次の数１７で示されるパラメータは、ロバスト変位測定を提供する（図１（ウ）参照）。

＜後処理＞

後処理において、以下のようにデータ処理される。第１に、変位後に取得される対象Ｂスキャンｇ_ｔ（ｘ，ｙ）が選択され（図１（エ）参照）、相関係数ρ_{ｍｅａｓｕｒｅｍｅｎｔ}(ｘ，ｙ)（図１（オ）参照）とドップラー位相シフトΔφ(ｘ，ｙ)（図１（カ）参照）は、参照Ｂスキャンｇ_ｒ（ｘ，ｙ）と対象Ｂスキャンｇ_ｔ（ｘ、ｙ）との間で計算される。

相関係数ρ_{ｍｅａｓｕｒｅｍｅｎｔ}(ｘ，ｙ)は、式（６）を用いて取得され、測定されたドップラー位相シフトΔφ(ｘ，ｙ)は、次の数１８に示される式（９）によって与えられる（図１（カ）参照）。

測定された位相シフトは、優先順位付き経路追従法を用いてアンラップされる（[20] D. C. Ghiglia and M. D. Pritt, Two-dimensional phase unwrapping: theory, algorithms, and software, (Wiley, New York, 1998). 参照）。

軸方向の変位Δｚ（ｘ，ｙ）は、式（５）を用いてアンラップされた位相シフトから変換される（図１（キ）参照）。横方向の変位Δｘ（ｘ，ｙ）は、測定された相関係数ρ_{ｍｅａｓｕｒｅｍｅｎｔ}(ｘ，ｙ)、測定された軸方向の変位Δｚ（ｘ，ｙ）、求められた次の数１９に示される定数及び分解能パラメータを式（４）に代入することによって取得される（図１（ク）参照）。

（実験例）
本発明に係る光コヒーレンストモグラフィーによる変位測定装置及び測定方法の効果を確認するために、本発明者らは実験を行った。以下、その実験例について説明する。

実験装置として利用した光コヒーレンストモグラフィーは、図２に示したＦＤ−ＯＣＴ１である。その具体的な仕様等について以下説明するが、もちろん、以下の装置は、実験装置を具体的に示したに過ぎず、本実施形態の装置は、下記の実験装置に限定されない。

実験装置のＦＤ−ＯＣＴ１としては、広帯域光源２として１μｍＳＬＤ光源（スーパー・ルミネッセント・ダイオード）を使用したＳＤ−ＯＣＴ（スペクトラルドメインＯＣＴ）が用いられる。

使用した１μｍＳＬＤ光源は、１．０２μｍの中心波長を持ち、１００ｎｍのスペクトル帯域を持つ。理論的な軸方向の分解能は、組織の中で３．４μｍであり、測定された軸方向の解像度は、組織の中で５．９μｍであった。

光源からの光は、ビームスプリッター５（例えば、５０／５０ファイバーカップラ）によって分割され、サンプル（物体）アームと参照アームに送られる。サンプルアームにおいて、光は、ＯＣＴのＢスキャンを取得するために、ガルバノミラー７によって横方向に走査される。

画像サイズは、１０２４画素（水平）×５１２画素（垂直）である（１．５ｍｍ×１．３ｍｍ）。スペクトラル干渉縞は、分光器（スペクトロメータ）４によって検出される。スペクトロメータには、高速InGaAsラインセンサが用いられた。

実験装置は、９１，９１１ Ａ−ｓｃａｎｓ/ｓのラインレートで駆動される。プローブビームパワーは、１．８６ｍＷ、測定された感度は、２．２ｄｂの再結合損失を含む９３ｄｂであった。

被計測物体８として組織ファントム（phantom）（疑似組織）を使用し、これは、図２（ａ）に示すように、サンプルアームで測定される。ファントムは、ゼラチン粉１０ｇ、intralipid（脂肪乳剤イントラリピッド）、温水５０ｍｌによって作成された。

組織ファントムは、移動ステージによって変位される。最大移動範囲は２０μｍであり、その際の軸上の絶対精度が１μｍである。対物レンズの焦点距離は、１６ｍｍであり、ガウシアンビームウェストサイズは、１０μｍであった。

また、被計測物体８として豚眼の網膜を使用し、これも、図２（ｂ）に示すように、サンプルアームで測定された。豚眼は、死後１２時間内に摘出され測定された。５３２ｎｍのレーザ照射光学系１００からの光凝固レーザーは、豚眼を照射するために用いられる。レーザ出力は、０．２ｍＷ、照射時間は、０．０２秒であった。

光凝固レーザーと広帯域光源２からのＯＣＴ照射ビームとの中心軸は、図２（ｂ）に示すように、同軸に位置合わせされる。レーザが照射された組織は、照射スポットの中心から膨張あるいは収縮し、横方向の変位は、照射ビームの軸を中心に軸対称となると推測される。

相関係数を用いた変位測定は、実験的に確認された。移動ステージに配置された組織ファントムは、複数のＢスキャンを取得中において、段階的に、横方向又は軸方向に変位される。

サンプルの位置は、０．２μｍの大きさで増加され、各位置で毎回１６Ｂスキャンが取得される。全体として、２５６Ｂスキャンが取得され、変位の合計は、３μｍであった。変位は、参照Ｂスキャンと他のＢスキャンとの間で測定された。

相関係数は、式（６）を用いて、参照Ｂスキャンと他のＢスキャンとの間で計算される。計算上シミュレートされた相関係数は、上記のように、参照Ｂスキャンと、デジタル処理によってシフトされたＢスキャンと、の間で計算される。それらは、図３（ａ）に示されるように、ランダムに選択された位置で観察される。

測定された相関係数（図４において矩形の棒状部に示される）と、計算上シミュレートされた相関係数（線と棒状部の交点で示される）は、図４に示すように、変位の関数としてプロットされる。なお、図４（ａ）は組織ファントムについての結果であり、図４（ｂ）は豚眼の網膜に関する結果である。

計算上シミュレートされた相関係数は、単調に減少し、測定された相関係数と一致する。それらの結果は、相関係数の低下が、サンプルの位置に依存していることを示している。これは、最小二乗フィッティングが、全ての画素に対して要求されることを示している。

しかしながら、もし、ある閾値を変位が超えると、非相関の影響が重大となる。相関係数は、振動を開始し、もはや信頼性がなくなる。そこで、フィッティングは、信頼性がある範囲内でのみ行われる。この制限は、本手段の測定可能範囲を示している。実験装置のシステムでは、閾値は、大体、１０μｍより少なかった。フィッティングの範囲は、３μｍより少なかった。

なお、レーザ照射された被計測物体である組織ファントム、豚眼の網膜等の変位は、数ミクロンよりも小さいことが予想されるので、測定範囲の制限は、大きな問題ではない。つまり、本発明の変位測定装置及び測定方法は、温度膨張のモニタリングに適している。

測定された横方向と軸方向の変位は、移動ステージの位置検知出力と比較され、全てのＢスキャンに関して記憶される。横方向と軸方向の変位測定の能力は、図３（ｂ）に示されるように、選択された関心領域において評価される。強度が１０ｄｂより少ない画素と、移動ステージの過渡的な反応（transient response）を受けた画素は、無視される。

以下のパラメータは、変位を取得するために用いられる。

上記数２０で示されるパラメータは、計算上シミュレートされた相関係数を、数学的なモデル関数にフィッティングすることによって与えられる。

上記数２１で示されるパラメータは、測定された相関係数を、数学的なモデル関数にフィッティングすることによって与えられる。

上記数２２で示されるパラメータは、従来（conventionally）で使用されるように一定に設定される。変位測定の能力は、図５に要約される。図５（ａ）は測定された横方向における変位をリニアにプロットしたものであり、上記関心領域での最頻値がプロットされている。

図５（ｂ）は測定された横方向の変位のＲＭＳエラーがプロットされている。図５（ｃ）は測定された軸方向の変位のリニアなプロットである。図５（ｄ）は測定された軸方向の変位のＲＭＳエラーのプロットである。

上記数２３に示されるフィッティングされたパラメータは、最小のＲＭＳエラーを提供する。

上記数２４に示されるフィッティングされたパラメータは、最大のＲＭＳエラーを示した。

言い換えれば、固定されたパラメータを用いて測定された変位は、より大きな変化があった。それは、パラメータの算出がロバストな変位測定に導くことを示している。

一方、フィッティングされた次の数２５に示されるパラメータは、固定されたパラメータよりも良い能力を示した。それは、本発明の変位測定装置及び測定方法が、ロバストな変位測定が可能であることを示している。ＲＭＳエラーは、２．８μｍの変位で、０．６μｍより少ない。

位置センサーからの大きなずれは、０．４μｍより小さな変位で見つかった。これは、１より大きな相関係数の変位は、０に強制的に変換されるためである。

また、フィッティングされた次の数２６に示されるパラメータの場合、位置センサーからのずれは、変位と共に増加する。これは、主に、計算上のシミュレーションが、追加的なノイズの寄与を考慮しないからである。

レーザ照射は、Ｂスキャンの連続が取得される間において、豚眼の網膜に行われた。レーザ照射による温度に伴う組織の変化は、ＯＣＴ強度、ドップラー位相シフト、相関係数、及び変位マップによってモニタリングされる。色分けされた変位マップは、色（Ｈｕｅ）が、ずれの方向を示し、その彩度が、変位の大きさを示しており、明度は、ＯＣＴ強度によって与えられる。

レーザ照射中の代表的な画像は、図６に示される。図６（ａ）はＯＣＴのＢスキャンの強度画像であり、図６（ｂ）はＨＳＬ（色相（Hue）、彩度（Saturation）、輝度（Lightness/Luminance または Intensity））での変位マップであり、マップのＨｕｅと膨張は、変位ベクトルの方向と大きさによって与えられる。

図６（ｃ）は測定された相関係数を示す画像であり、図６（ｄ）は測定されたドップラー位相シフトを示す画像であり、アンラッピング前のドップラー位相シフトを示す画像である。

これらの画像からみて、また、組織ファントムについて動的な温度変化が、見られた。また、内網膜での大きな位相シフトと、網膜色素上皮複合体（ＲＰＥ complex）での重大な相関係数の低下と、が見られた。ＯＣＴ強度は、十分な変化を示さなかった。測定された変位は、明らかに変位マップによって視覚化された。

豚の網膜等についての横方向と縦方向の変位は、図７（ａ）と図７（ｂ）にそれぞれ示された。内網膜層の照射領域は、図７（ａ）に示すように軸方向に変位される。網膜色素上皮複合体は、図７（ｂ）に示すように、非常に大きく、横方向に変位された。

時間依存性を理解するため、内網膜層の軸方向での変位と、網膜色素上皮複合体の横方向での変位は、図８に示すように、ＢスキャンＩＤの関数としてプロットされた。変位の増加は、レーザ照射中において明らかである。冷却状態において、軸方向における変位の緩やかな減少が見受けられ、漸近的に一定値に向かっている。一方、横方向の変位は増加した状態のままであり、軸方向の変位よりも十分に大きい。

＜レーザ照射中の温度変化＞
図８に示すように、豚の網膜内層の軸方向における変位の緩やかな減少は、レーザ照射後だけに見られた。横方向の変位は、増加した状態のままであり、軸方向の変位より十分に大きかった。

それは、凝固レーザは、ＲＰＥ複合体の散乱・吸収分布を大きく変形させ、非弾性変形を生じさせることを示している。変形は、横方向の変位から区別することは不可能であるがこの結果は、本発明に係る光コヒーレンストモグラフィーによる変位測定装置及び測定方法がレーザ照射された組織の直接的な温度変化を検出するのに有用である事を示している。

（変形例）
なお、上記説明では、ＦＤ−ＯＣＴ１として、広帯域光源２及び分光器４を用いたＳＤ−ＯＣＴ（スペクトラルドメインＯＣＴ）を例にして説明したが、これに限定されない。例えば、ＦＤ−ＯＣＴ１として、波長走査型光源及びポイントセンサを用いたＳＳ−ＯＣＴ（スウィプトソースＯＣＴ）であっても、本実施形態の本発明に係る光コヒーレンストモグラフィーによる変位測定装置及び測定方法の適用は可能である。

以上、本発明に係る光コヒーレンストモグラフィーによる変位測定装置及び測定方法を実施するための形態を実施例に基づいて説明したが、本発明はこのような実施例に限定されることなく、特許請求の範囲記載の技術的事項の範囲内で、いろいろな実施例があることは言うまでもない。

上記構成から成る本発明に係る光コヒーレンストモグラフィーによる変位測定装置及び測定方法は、眼科治療、動植物の研究、半導体研究、工業材料の研究等の分野での利用が最適である。

１ ＦＤ−ＯＣＴ
２ 広帯域光源
３ 干渉計
４ 分光器
５ ビームスプリッタ
６、９ レンズ
７ ガルバノミラー
８ 被測定物体（（疑似）組織等）
１０ 参照鏡
１００ レーザ照射光学系
３００ コンピュータ

Claims (4)

1. 被測定物体の奥行き方向の軸に平行な２次元断層画像を同一位置に関して複数取得し、取得された複数の断層画像に基づいて被測定物体の変位の測定が可能な光コヒーレンストモグラフィー装置であって、
   前記複数の断層画像のうちの１つの断層画像を基準画像とし、前記基準画像の微小領域に対して他の断層画像の対応する微小領域を軸方向及び前記軸方向に直交する断層画像面内の方向にシフトすることによって、相関係数を計算する相関係数算出部と、
   前記相関係数算出部からの出力に基づいて、前記軸方向の変位と前記軸方向に直交する断層画像面内の方向の変位とを相関係数に関係づける相関係数の数学モデルのパラメータを算出する相関係数モデル算出部と、
   変位測定対象の断層画像と前記基準画像との位相シフトを求める位相シフト算出部と、
   前記位相シフト算出部からの出力に基づいて前記軸方向の変位を求める軸方向変位算出部と、
   前記軸方向変位算出部からの出力及び前記相関係数モデル算出部からの出力に基づいて前記軸方向に直交する断層画像面内の方向の変位を求める横方向変位算出部と、
   を有する光コヒーレンストモグラフィー装置。
2. 前記相関係数算出部は、前記基準画像の信号対ノイズ比と前記他の断層画像の信号対ノイズ比とをそれぞれ算出し、次に示す（２）式に基づいて、前記信号対ノイズ比の影響なしで信号相関係数を計算し、前記信号相関係数を前記相関係数に置き換える請求項１に記載の光コヒーレンストモグラフィー装置。
   なお、上記の（２）式において、ρ_ｓ は前記信号相関係数を表し、ρは前記複数の断層画像の複素相関係数を表し、ＳＮＲ_１は前記基準画像の信号対ノイズ比の期待値を表し、ＳＮＲ_２は前記他の断層画像の信号対ノイズ比の期待値を表す。
3. 被測定物体の奥行き方向の軸に平行な２次元断層画像を同一位置に関して複数取得し、取得された複数の断層画像に基づいて被測定物体の変位を測定する光コヒーレンストモグラフィーによる変位測定方法であって、
   前記複数の断層画像のうちの１つの断層画像を基準画像とし、前記基準画像の微小領域に対して他の断層画像の対応する微小領域を軸方向及び前記軸方向に直交する断層画像面内の方向にシフトすることによって、相関係数を計算する相関係数算出工程と、
   前記相関係数算出工程で算出された前記相関係数に基づいて、前記軸方向の変位と前記軸方向に直交する断層画像面内の方向の変位とを相関係数に関係づける相関係数の数学モデルのパラメータを算出する相関係数モデル算出工程と、
   変位測定対象の断層画像と前記基準画像との位相シフトを求める位相シフト算出工程と、
   前記位相シフト算出工程で算出された前記位相シフトに基づいて前記軸方向の変位を求める軸方向変位算出工程と、
   前記軸方向変位算出工程で算出された前記軸方向の変位及び前記相関係数モデル算出工程で算出された前記パラメータに基づいて前記軸方向に直交する断層画像面内の方向の変位を求める横方向変位算出工程と、
   を有する光コヒーレンストモグラフィーによる変位測定方法。
4. 前記相関係数算出工程において、前記基準画像の信号対ノイズ比と前記他の断層画像の信号対ノイズ比とをそれぞれ算出し、次に示す（２）式に基づいて、前記信号対ノイズ比の影響なしで信号相関係数を計算し、前記信号相関係数を前記相関係数に置き換える請求項３に記載の光コヒーレンストモグラフィーによる変位測定方法。
   なお、上記の（２）式において、ρ_ｓ は前記信号相関係数を表し、ρは前記複数の断層画像の複素相関係数を表し、ＳＮＲ_１は前記基準画像の信号対ノイズ比の期待値を表し、ＳＮＲ_２は前記他の断層画像の信号対ノイズ比の期待値を表す。