JPWO2016031697A1 - 軟骨診断装置および診断用プローブ - Google Patents

Abstract

軟骨診断装置１は、ＯＣＴの光学系を含む光学ユニット２と、光学ユニット２に接続される一方、先端部が関節腔に挿入可能に構成され、光学ユニット２からの光を軟骨に導いて走査させるための光学機構と、軟骨に対して所定の荷重を負荷可能な荷重機構とを含むプローブ４と、荷重機構および光学機構の駆動を制御し、それらの駆動に基づいて光学ユニット２から出力された光干渉信号を処理し、荷重負荷による軟骨内部の変形に伴う所定の力学特徴量の変化をその軟骨の断層位置に対応づけて演算し、その力学特徴量の変化に基づいて軟骨組織の損傷度合いを演算する制御演算部６と、軟骨組織の損傷度合いを断層可視化する態様で表示する表示装置８と、を備える。

Description

本発明は、軟骨の変性度を診断するための装置に関する。

軟骨は荷重衝撃の緩和や関節滑動性の向上等の重要な役割を担っているが、血液の循環がなく、自己治癒が困難な組織である。高齢者の多くは軟骨の磨耗による変形性膝関節症（Osteoarthritis：以下「ＯＡ」と表記する）を発症しており、特に高齢化社会を迎える国ではその診断治療法の確立が求められている。関節軟骨は８０％の水分と２０％のマトリクスを含み、そのマトリクスはコラーゲンとプロテオグリカンから構成される。特に、そのコラーゲン線維に拘束されるプロテオグリカンは軟骨内部の流動特性を決定し、軟骨の優れた粘弾性特性に大きく関与すると考えられている。ＯＡは、その軟骨の粘弾性特性の損失によって発症する。

近年の医療診断技術の発達に伴い、光コヒーレンス断層画像（Optical Coherence Tomography：以下「ＯＣＴ」という）が開発されている。このＯＣＴによれば、非侵襲、非接触にて生体組織内部をマイクロ断層可視化できる。また、２次元ＯＣＴ断層画像の取得レートはビデオレート以上であり、高時間分解能を有している。そこで、このＯＣＴを用いて軟骨力学特性を断層可視化する手法も提案されている（非特許文献１）。この手法によれば、軟骨に所定の応力を負荷してその応力緩和過程を断層可視化することで、正常軟骨と変性軟骨とを識別できる可能性がある。

「Dynamic Optical Coherence Straingraphyを用いた酵素処理による変性軟骨のマイクロ力学断層可視化による検討」岡本祐太，佐伯壮一，峯松孝幸、日本機械学会2013年度年次大会講演論文集(J024032)

非特許文献１に記載の技術は、軟骨に関するマイクロメカニクス可視化診断への有力な足掛かりとなるものの、診断効率や診断精度の面で改善の余地があった。また、生きた人体の軟骨を診断対象とする場合の具体的装置構成について開示するものではなかった。このため、未だ実用化には到っていない。

本発明はこのような課題に鑑みてなされたものであり、その目的の一つは、ＯＣＴを用いた軟骨診断をより実用に供し得るものとすることにある。

本発明のある態様は、関節軟骨を診断するための軟骨診断装置である。この軟骨診断装置は、光コヒーレンストモグラフィーを用いる光学系を含む光学ユニットと、光学ユニットに接続される一方、先端部が関節腔に挿入可能に構成され、光学ユニットからの光を軟骨に導いて走査させるための光学機構と、軟骨に対して所定の変形エネルギーを付与するための負荷装置とを含むプローブと、負荷装置および光学機構の駆動を制御し、それらの駆動に基づいて光学ユニットから出力された光干渉信号を処理し、変形エネルギーの付与による軟骨内部の変形に伴う所定の力学特徴量の変化をその軟骨の断層位置に対応づけて演算し、その力学特徴量の変化に基づいて軟骨組織の損傷度合いを演算する制御演算部と、軟骨組織の損傷度合いを断層可視化する態様で表示する表示装置と、を備える。

この態様によると、所定の変形エネルギーが付与されたときに軟骨に生じる力学特徴量の変化と、軟骨組織の損傷度合いとの間に対応関係があることを利用した演算処理により、診断対象である軟骨の損傷度合いが断層可視化される。このため、医師等がその断層可視化された表示を確認することにより、軟骨診断を容易に行えるようになる。すなわち、ＯＣＴを用いた軟骨診断をより実用に供し得るものとすることができる。

本発明の別の態様は、診断用プローブである。この診断用プローブは、光コヒーレンストモグラフィーを用いる光学系を含む光学ユニットに接続されることにより、対象部位の診断を可能とするプローブであって、対象部位に当接可能に構成された先端部と、光学ユニットからの光を軟骨に導いて走査させるための光学機構と、対象部位に対して所定の変形エネルギーを付与するための負荷装置と、を備える。

この態様の診断用プローブを対象部位としての軟骨に当接させ、光学ユニットを作動させることにより、上述した軟骨の損傷度合いを断層可視化することが可能となる。すなわち、ＯＣＴを用いた軟骨診断をより実用に供し得るものとすることができる。

本発明によれば、ＯＣＴを用いた軟骨診断をより実用に供し得るものとできる。

第１実施例に係る軟骨診断装置の構成を概略的に表す図である。 軟骨診断装置を構成する診断プローブの構成を概略的に表す図である。 再帰的相互相関法による処理手順を概略的に示す図である。 サブピクセル解析による処理手順を概略的に示す図である。 各処理により得られる結果を示す図である。 軟骨診断処理の処理手順を表す図である。 軟骨診断処理の処理手順を表す図である。 実験装置の構成を概略的に表す図である。 実験に用いた試験片を表す図である。 実験に際して試験片に与えた荷重負荷状態を示すグラフである。 応力緩和試験によるひずみ速度分布の時間変化を断層可視化した図である。 応力緩和試験によるひずみ速度分布の時間変化を断層可視化した図である。 正常軟骨と模擬変性軟骨の減衰係数の比較を表す図である。 制御演算部により実行される軟骨診断処理の流れを示すフローチャートである。 図１４におけるＳ３６の変性度合い演算提示処理を詳細に示すフローチャートである。 第２実施例に係る変性度合い演算提示処理を詳細に示すフローチャートである。

本発明の一実施形態は、関節軟骨を診断するための軟骨診断装置である。この軟骨診断装置は、軟骨に所定の変形エネルギー（荷重）を負荷する一方でＯＣＴを用いて軟骨の断層画像を撮影し、その荷重に対する軟骨組織の挙動からその損傷度合い（変性度）を演算するものである。その演算結果が軟骨の断層画像の態様で表示装置に可視化されるため、医師等がその断層画像を見ることで軟骨診断を行うことができる。

この軟骨診断装置は、ＯＣＴの光学系を含む光学ユニットと、その光学ユニットに接続されるプローブを備える。このプローブは、関節腔に挿入可能な先端部を有するとともに、ＯＣＴの演算処理に用いられる物理量を取得するための光学機構と荷重機構（「負荷装置」として機能する）を備える。プローブの先端部は、例えばシリンジ針に挿通されるなどして関節腔に挿入されるものでもよい。プローブの先端を軟骨の表面（又はその近傍）に当接させた状態で荷重機構を駆動することにより、軟骨に対して所定の荷重が負荷される。その荷重に対する軟骨組織の応答をＯＣＴにより断層計測することにより、軟骨の変性度を評価することができる。後述のように、変性度の進行度合いに応じてその応答が異なる形で表れるからである。光学機構は、この軟骨組織の挙動を断層画像として取得するために光学ユニットからの物体光を軟骨に照射して走査させ、その反射光を取得して光学ユニットに送る。

光学ユニットは、その反射光と参照光とを合波し、その光干渉信号を制御演算部へ送る。制御演算部は、荷重機構および光学機構の駆動を制御する一方でその光干渉信号を処理し、上述した荷重負荷による軟骨内部の変形に伴う力学特徴量の変化をその軟骨の断層位置に対応づけて演算する。ここでいう「力学特徴量」は、軟骨組織の変形ベクトルの空間分布に基づいて得られるものでよい。例えば、変形ベクトルを時間微分した変形速度ベクトルであってもよいし、変形速度ベクトルをさらに空間微分したひずみ速度テンソルであってもよい。

荷重機構による荷重負荷方法については、測定対象（軟骨）に対して一定のひずみを与えて応力の時間変化を測定する応力緩和法に基づくものでもよい。あるいは、測定対象に対して動的ひずみを与えて応力の最大値および位相差を測定する動的粘弾性法に基づくものでもよい。あるいは、測定対象に対して一定の大きさの応力を与えてひずみの時間変化を測定するクリープ法に基づくものでもよい。制御演算部は、上記力学特徴量の変化に基づいて軟骨組織の損傷度合いを演算し、それを表示装置の画面に断層可視化する。

制御演算部は、応力緩和法に基づいて以下の処理を実行してもよい。すなわち、荷重機構の駆動により軟骨に所定の荷重を負荷した後に応力緩和させ、その応力緩和中に光学ユニットにより撮影された前後の断層画像データに基づき、軟骨の断層位置に対応した変形速度ベクトルを力学特徴量として演算してもよい。そして、その変形速度ベクトルの変化から得られる減衰係数の断層分布を、軟骨組織の損傷度合いの分布と等価として断層可視化させてもよい。なお、ここでいう「前後の断層画像データ」は、連続撮影された２枚の断層画像によるものでもよいし、連続しないものの所定時間をあけて撮影された２枚の断層画像によるものでもよい。

あるいは、制御演算部は、力学特徴量として演算した変形速度ベクトルを空間微分することによりひずみ速度テンソルを演算し、そのひずみ速度の減衰係数の断層分布を軟骨組織の損傷度合いの分布と等価として断層可視化させてもよい。このようにひずみ速度テンソルを演算することで、減衰係数の断層分布に軟骨組織の圧縮、膨張、剪断などの物理的意味合いをもたせることができる。なお、この物理的な意味合いから、軟骨変性度を力学的変形能として直接的に医師等に伝えることもできる。

制御演算部は、順次取得する断層画像データの相互相関、輝度勾配の変位および輝度パターンの変形に基づき変形速度ベクトルをサブピクセル精度にて演算してもよい。ＯＣＴの解像度には限界があるため、このようなサブピクセル解析を導入することにより、より高精度な軟骨診断が実現可能となる。

制御演算部は、動的粘弾性法に基づいて以下の処理を実行してもよい。すなわち、荷重機構の駆動により軟骨に所定の動的荷重（動的力）を負荷し、その動的荷重による軟骨の振動中に光学ユニットにより撮影された前後の断層画像データに基づき、軟骨の各断層位置における変形速度ベクトルを力学特徴量として演算し、その変形速度ベクトルのうち加振周波数成分である特定変形速度ベクトルの振幅および位相差の少なくとも一方に基づき軟骨の粘弾性の断層分布を演算し、その粘弾性の断層分布を軟骨組織の損傷度合いの分布と等価として断層可視化させてもよい。軟骨の粘弾性がその軟骨の変性度に応じて変化することに着目したものである。なお、ここでいう「前後の断層画像データ」も、連続撮影された２枚の断層画像によるものでもよいし、連続しないものの所定時間をあけて撮影された２枚の断層画像によるものでもよい。

あるいは、軟骨の各断層位置におけるひずみ速度テンソルを力学特徴量として演算し、そのひずみ速度テンソルのうち加振周波数成分である特定ひずみ速度テンソルの振幅および位相差の少なくとも一方に基づき軟骨の粘弾性の断層分布を演算してもよい。

また、荷重機構が、軟骨に付与される荷重（力）を計測可能なセンサを含んでもよい。制御演算部は、演算された力学特徴量の加振周波数成分の振幅および位相差の少なくとも一方と、センサにより計測された荷重値に基づき軟骨の粘弾性の断層分布を演算し、その粘弾性の断層分布を軟骨組織の損傷度合いの分布と等価として断層可視化させてもよい。

なお、上記プローブは、例えば手術中の内臓を診断するなど、軟骨以外の対象部位の診断に適用することも可能である。すなわち、プローブが、ＯＣＴの光学系を含む光学ユニットに接続されることにより、所定の対象部位の診断を可能とする診断用プローブとして構成されてもよい。このプローブは、対象部位に当接可能に構成された先端部と、光学ユニットからの光を軟骨に導いて走査させるための光学機構と、対象部位に対して所定の変形エネルギー（応力）を負荷するための荷重機構（負荷装置）とを備える。

また、上記技術を利用した軟骨診断方法を構築してもよい。この方法は、軟骨に対して所定の変形エネルギー（荷重）を負荷するステップと、変形エネルギーの付与（荷重の負荷）に応じた軟骨の変形度合いを光コヒーレンストモグラフィーにより断層画像として表示させるステップと、断層画像に基づいて軟骨組織の損傷度合いを診断するステップとを備えるものでよい。

以下、図面を参照しつつ本発明を具体化した実施例について詳細に説明する。
［第１実施例］
図１は、第１実施例に係る軟骨診断装置の構成を概略的に表す図である。図２は、軟骨診断装置を構成する診断プローブの構成を概略的に表す図である。本実施例の軟骨診断装置は、診断対象である関節軟骨に所定の応力を負荷し、その応力に対する軟骨組織の変形度合いを断層可視化することにより、軟骨組織の損傷度合い（変性度）を診断可能とするものである。この断層可視化にＯＣＴを利用する。

図１に示すように、軟骨診断装置１は、ＯＣＴを用いる光学系を含む光学ユニット２、光学ユニット２に接続される診断用のプローブ４、ＯＣＴにより得られた光干渉データに基づいて軟骨組織の損傷度合いを演算する制御演算部６、および軟骨組織の損傷度合いを断層可視化する態様で表示する表示装置８を備える。図示の例では、光学ユニット２としてマッハツェンダー干渉計をベースとした光学系が示されているが、マイケルソン干渉計その他の光学系を採用することもできる。

本実施例では、ＯＣＴとして波長走査型レーザを光源とするＳＳ−ＯＣＴ（Swept Source OCT）を用いるが、ＴＤ−ＯＣＴ（Time Domain OCT）、ＳＤ−ＯＣＴ（Spectral Domain OCT）その他のＯＣＴを用いることもできる。ＳＳ−ＯＣＴは、ＴＤ−ＯＣＴとは異なり、参照鏡走査などの機械的な光遅延走査を必要としないため、高い時間分解能と高い位置検出精度が得られる点で好ましい。

光学ユニット２は、光源１０、オブジェクトアーム１２、リファレンスアーム１４、および光検出器１６を備える。各光学要素は、光ファイバーにて互いに接続されている。光源１０から出射された光は、カプラ１８（ビームスプリッタ）にて分けられ、その一方がオブジェクトアーム１２に導かれる物体光となり、他方がリファレンスアーム１４に導かれる参照光となる。オブジェクトアーム１２に導かれた物体光は、サーキュレータ２０を介してプローブ４に導かれ、診断対象である軟骨に照射される。この物体光は、軟骨の表面および断面にて後方散乱光として反射されてサーキュレータ２０に戻り、カプラ２２に導かれる。

一方、リファレンスアーム１４に導かれた参照光は、サーキュレータ２４を介して反射鏡２６に導かれる。このとき、参照光は、コリメートレンズ２８を経て集光レンズ３０にて反射鏡２６に集光される。この参照光は、反射鏡２６にて反射されてサーキュレータ２４に戻り、カプラ２２に導かれる。すなわち、物体光と参照光とがカプラ２２にて合波（重畳）され、その干渉光が光検出器１６により検出される。

プローブ４は、光学ユニット２のオブジェクトアーム１２を構成し、その本体３２にはサーキュレータ２０から延びる光ファイバー３４が挿通されている。本体３２の先端にはシリンジ針３６が取り付けられており、プローブ４の先端部は、そのシリンジ針３６を介して生体内（患者の膝Ｋ）に挿入可能とされている。本体３２には、光学ユニット２からの光を膝Ｋの軟骨に導いて走査させるための光学機構と、軟骨に対して所定の荷重（応力）を負荷するための荷重機構と、これらの機構を駆動するための駆動部３８（アクチュエータ）が設けられている。光学機構は、軟骨に向けて光を照射し、その反射光をオブジェクトアーム１２に導く。各機構の詳細については後述する。

カプラ２２にて合波された干渉光は、光検出器１６に入力される。光検出器１６は、これを光干渉信号（干渉光強度を示す信号）として検出する。この光干渉信号は、Ａ／Ｄ変換器４０を介して制御演算部６に入力される。Ａ／Ｄ変換器４０は、光検出器１６から出力されたアナログ信号をデジタル信号に変換して制御演算部６へ出力する。

制御演算部６は、ＣＰＵ、ＲＯＭ、ＲＡＭ、ハードディスクなどを有し、これらのハードウェア、およびソフトウェアによって、光学ユニット２の光学系全体の制御と、プローブ４の駆動制御と、ＯＣＴによる画像出力のための演算処理を行う。制御演算部６は、プローブ４の荷重機構および光学機構の駆動を制御し、それらの駆動に基づいて光学ユニット２から出力された光干渉信号を処理し、ＯＣＴによる軟骨の断層画像を取得する。そして、その断層画像データに基づき、後述の手法により軟骨組織の損傷度合いを演算する。

表示装置８は、例えば液晶ディスプレイからなり、制御演算部６にて演算された軟骨組織の損傷度合いを断層可視化する態様で画面に表示させる。

図２に示すように、プローブ４は、筒状の本体３２の先端に穿刺部としてのシリンジ針３６が取り付けられる。シリンジ針３６の内方には有底円筒状のシース５０が同軸状に設けられ、本体３２の内方にはシース５０を軸線方向に駆動可能な圧電素子５２（ピエゾ素子）が設けられている。圧電素子５２は、本体３２の内部に設けられた図示略の支持機構に支持されている。シース５０は、その先端が閉止されており、後端開口部が圧電素子５２の先端部に固定されている。シース５０の先端面は、患者の軟骨Ｊの表面に当接可能な当接面となっている。

シース５０の内方には光ファイバー３４の先端部が挿通され、その先端にＧＲＩＮレンズ５４が連設され、さらにその先端にプリズムミラー５６が連設されている。すなわち、光ファイバー３４、ＧＲＩＮレンズ５４およびプリズムミラー５６は、軸線方向に一体化された導光体を構成している。ＧＲＩＮレンズ５４は、レンズ内の屈折率を連続的に変えることで集光機能を発揮する。プリズムミラー５６は、その先端面が軸線に対して傾斜した反射面とされている。なお、本実施例では、照射ビームの径と焦点距離・ビームスポット径などを整える光学素子としてＧＲＩＮレンズ５４を採用しているが、同様の機能を発揮できる他の光学素子を採用してもよい。また、ビームを指定方向に反射する光学素子としてプリズムミラー５６を採用しているが、同様の機能を発揮できる他の光学素子を採用してもよい。

光ファイバー３４によりプローブ４に導かれた光は、ＧＲＩＮレンズ５４およびプリズムミラー５６を介してシース５０の先端から出射され、軟骨Ｊに照射される（点線矢印参照）。それにより軟骨Ｊの表面又は断面にて反射された光がＧＲＩＮレンズ５４に取り込まれ、光ファイバー３４を介して光学ユニット２のオブジェクトアーム１２に導かれる。

圧電素子５２には軸線に沿って貫通する挿通孔５８が設けられている。光ファイバー３４は、この挿通孔５８に挿通され、シース５０内に導入されている。このような構成により、光ファイバー３４に機械的負荷をかけることなく、圧電素子５２による荷重をシース５０に負荷することができる。圧電素子５２への通電によりその圧電素子５２が軸線方向に伸長し、シース５０を軸線方向先端側に駆動することができる。この駆動力は、シース５０の先端面による押圧力（圧縮応力）となって軟骨Ｊに負荷される。すなわち、シース５０および圧電素子５２が、「負荷装置（荷重機構）」として機能する。

なお、シース５０は、図示しない付勢部材により後方に付勢されており、圧電素子５２の非通電時にはシリンジ針３６の内方に収容されるように構成されている。圧電素子５２に通電されると、図示のようにシース５０の先端面がシリンジ針３６から突出できる構成とされている。

本体３２の内方における圧電素子５２の後方には、光ファイバー３４を軸線周りに回動させるための回動機構６０（光ロータリージョイント）が設けられている。回動機構６０は、光ファイバー３４と同軸状に固定された回転子を有する。回動機構６０への通電により、光ファイバー３４およびＧＲＩＮレンズ５４を自軸周りに回動させ、プリズムミラー５６の向きを変化させることができる。すなわち、光ファイバー３４の先端部、ＧＲＩＮレンズ５４、プリズムミラー５６および回動機構６０が、「光学機構」として機能する。回動機構６０の駆動により、物体光を軟骨Ｊに対して走査させることができる。圧電素子５２および回動機構６０への通電制御は、制御演算部６によってなされる。

以下、軟骨組織の損傷度合いを提示するまでの演算処理方法について説明する。
上述のように、オブジェクトアーム１２を経た物体光（軟骨Ｊからの反射光）と、リファレンスアーム１４を経た参照光とが合波され、光検出器１６により光干渉信号として検出される。制御演算部６は、この光干渉信号を干渉光強度に基づく軟骨Ｊの断層画像として取得することができる。

ところで、ＯＣＴの光軸方向（奥行き方向）の分解能であるコヒーレンス長ｌ_ｃは、光源の自己相関関数によって決定される。ここでは、コヒーレンス長ｌ_ｃを自己相関関数の包括線の半値半幅とし、下記式（１）にて表すことができる。

ここで、λ_ｃはレーザビームの中心波長であり、Δλはレーザビームの半値全幅である。

一方、光軸垂直方向（ビーム走査方向）の分解能は、集光レンズによる集光性能に基づき、ビームスポット径Ｄの１／２とされる。そのビームスポット径Ｄは、下記式（２）にて表すことができる。

ここで、ｄ_Ｂは集光レンズに入射するビーム径、ｆは集光レンズの焦点である。このようにＯＣＴによる分解能には限界があるところ、本実施例では後述するサブピクセル解析の導入などにより、軟骨組織の診断をマイクロスケールにて行うことを可能にしている。以下、その詳細について説明する。

まず、ＯＣＴを利用したマイクロスケールの力学特性検出法について説明する。この検出法は、計測対象の変形前後の２枚のＯＣＴ断層画像にデジタル相互相関法を適用することで変形ベクトル分布を算出し、マイクロスケールにて生体組織内部のひずみ速度テンソル分布を断層検出する手法である。

変形ベクトル分布を算出する際には、繰り返し相互相関処理を実施する再帰的相互相関法(Recursive Cross-correlation method)を適用する。これは、低解像度において算出された変形ベクトルを参照し、探査領域を限定するとともに階層的に検査領域を縮小して相互相関法を適用する手法である。これにより、高解像度な変形ベクトルを取得することができる。さらに、スペックル・ノイズ低減法として、隣接検査領域の相関値分布との乗算を行う隣接相互相関乗法(Adjacent Cross-correlation Multiplication)を用いる。そして、乗算されることによって高ＳＮ化した相関値分布から最大相関値を探索する。

また、マイクロスケールにおける微小変形解析では、変形ベクトルのサブピクセル精度が重要となる。このため、輝度勾配を利用する風上勾配法(Up-stream Gradientmethod)と、伸縮および剪断を考慮した画像変形法(Image Deformation method)の両サブピクセル解析法を併用し、変形ベクトルの高精度検出を実現する。なお、ここでいう「風上勾配法」は、勾配法（オプティカルフロー法）の一種である。

このようにして得られた変形ベクトルを時間微分することにより変形速度ベクトル分布を演算し、それをさらに空間微分することによりひずみ速度テンソル分布を演算する。後述のように、軟骨組織のひずみ速度テンソルの減衰係数はその軟骨の損傷度合い（変性度）と等価とみることができるため、これを軟骨の損傷度合いとして断層可視化する。

以下、各手法について詳細に説明する。図３は、再帰的相互相関法による処理手順を概略的に示す図である。図４は、サブピクセル解析による処理手順を概略的に示す図である。図５は、各処理により得られる結果を示す図である。

（再帰的相互相関法）
図３（Ａ）〜（Ｃ）は、再帰的相互相関法による処理過程を示している。各図にはＯＣＴにより連続的に撮影される前後の断層画像が示されている。左側には先の断層画像（Image1）が示され、右側には後の断層画像（Image2）が示されている。

相互相関法とは、局所的なスペックル・パターンの類似度を下記式（３）に基づく相関値Ｒ_ｉｊより評価する方法である。そのため、図３（Ａ）に示すように、連続的に撮影された前後のＯＣＴ画像について、先の断層画像（Image1）に類似度の検査対象となる検査領域Ｓ１が設定され、後の断層画像（Image2）に類似度の探査範囲となる探査領域Ｓ２が設定される。

ここで、空間座標として、光軸方向にＺ軸、光軸と垂直方向にＸ軸を設定している。ｆ(Ｘ_ｉ,Ｚ_ｊ)とｇ(Ｘ_ｉ,Ｚ_ｊ)は、変形前後のＯＣＴ画像において設置された中心位置(Ｘ_ｉ,Ｚ_ｊ)の検査領域Ｓ１(Ｎ_ｘ×Ｎ_ｚピクセル)のスペックル・パターンを表している。

探査領域Ｓ２(Ｍ_ｘ×Ｍ_ｚピクセル)内における相関値分布Ｒ_i,j(ΔＸ,ΔＺ)を算出し、図３（Ｂ）に示すようにパターンマッチングを行う。実際には、下記式（４）に示すように、最大相関値を与える移動量Ｕ_ｉ,jを変形前後の変形ベクトルとして決定する。

なお、ｆ⁻とｇ⁻はｆ(Ｘ_ｉ,Ｚ_ｊ)とｇ(Ｘ_ｉ,Ｚ_ｊ)の検査領域Ｓ１内の平均値を表している。

本手法では、検査領域Ｓ１を縮小しながら相互相関処理を繰り返して空間解像度を高める再帰的相互相関法を採用している。なお、本実施例では解像度を上げる際には空間解像度が倍になるようにしている。図３（Ｃ）に示すように、検査領域Ｓ１を１／４に分割し、前階層にて算出された変形ベクトルを参照し、探査領域Ｓ２を縮小している。ここで探査領域Ｓ２も１／４に分割している。再帰的相互相関法を用いることで、高解像度において多発する過誤ベクトルの抑制を可能にしている。このような再帰的相互相関処理を施すことにより、図５（Ａ）に示すように変形ベクトルの解像度を高めることができる。

また、これに加え、下記式（５）により、演算中の座標を中心とする周囲８つの座標を含む合計９つの変形ベクトルの平均偏差σを用いた閾値を設定し、過誤ベクトルの除去を行い再帰処理に伴う誤差伝播を抑制する。

ここで、Ｕmはベクトル量の中央値を表し、閾値となる係数κは任意に設定した。

（隣接相互相関乗法）
本実施例では、スペックルノイズの影響を受けたランダム性の強い相関値分布から正確な最大相関値を決定する方法として、隣接相互相関乗法を導入している。下記式（６）により、隣接相互相関乗法では検査領域Ｓ１における相関値分布Ｒi,j(Δx,Δz)と、その検査領域Ｓ１にオーバーラップする隣接検査領域に対するＲi＋Δi,j(Δx,Δz)とＲi,j＋Δj(Δx,Δz)の乗算を行い、新たな相関値分布Ｒ'i,j(Δx,Δz)を用いて最大相関値を検索する。

これにより、相関値同士の乗算によってランダム性を低減させることが可能になる。上述した検査領域Ｓ１の縮小と共に干渉強度分布の情報量も減少するため、スペックル・ノイズを原因とする複数相関ピークの出現が計測精度の悪化を招いていると考えられる。一方、隣接境界同士の移動量には相関があるため、最大相関値座標付近では強い相関値が残存する。この隣接相互相関乗法の導入によって最大相関値ピークが明瞭化され計測精度が向上し、正確な移動座標を抽出することが可能となる。また、この隣接相互相関乗法をＯＣＴの各ステージに導入することで、誤差伝播が抑制され、スペックル・ノイズに対するロバスト性が向上する。それにより、高空間解像度においても高精度な変形ベクトルとひずみ分布の算出が可能となる。

（風上勾配法）
図４（Ａ）〜（Ｃ）は、サブピクセル解析による処理過程を示している。各図にはＯＣＴにより連続的に撮影される前後の断層画像が示されている。左側には先の断層画像（Image1）が示され、右側には後の断層画像（Image2）が示されている。

本実施例では、サブピクセル解析のために風上勾配法と画像変形法を採用する。最終的な移動量の算出は後述の画像変形法によるが、計算の収束性の問題から、画像変形法より先に風上勾配法を適用する。検査領域サイズが小さく高空間解像度の条件において、サブピクセル移動量を高精度検出する画像変形法及び風上勾配法を適用している。画像変形法におけるサブピクセル移動量の検出が困難な場合において、風上勾配法によりサブピクセル移動量を算出する。

サブピクセル解析では、注目点における変形前後の輝度差が各成分の輝度勾配と移動量によって表される。このため、検査領域Ｓ１内の輝度勾配データより最小二乗法を用いてサブピクセル移動量を決定することができる。本実施例では、輝度勾配を求める際に、サブピクセル変形前の風上側の輝度勾配を与える風上差分法を採用している。

サブピクセル解析は計測誤差からだけではなく、組織の散乱効果や偏光特性に複雑に関係したスペックル・ノイズに強く依存している。さらに、ひずみテンソル分布の算出には移動量の空間微係数が必要となるため、サブピクセルエラーが微係数算出の数値不安定を増大させ、ひずみテンソル分布の不安定化を招いてしまう。サブピクセル解析は様々な手法が存在するが、本実施例では検査領域サイズが小さく高空間解像度の条件においても、サブピクセル移動量を高精度検出する勾配法を採用している。

風上勾配法は、検査領域Ｓ１内の注目点の移動を、図４（Ａ）に示すピクセル精度に留まらず、図４（Ｂ）に示すサブピクセル精度にて算出するものである。なお、図中の各格子は１ピクセルを表している。実際には図示の断層画像と比較して相当小さいが、説明の便宜上、大きく表記している。この風上勾配法は、微小変形前後における輝度分布の変化を輝度勾配と移動量によって定式化する手法であり、ｆを輝度とすると、微小変形ｆ(x＋Δx,z＋Δz)をテイラー展開する下記式（７）として表される。

上記式（７）は、注目点の変形前後の輝度差が変形前の輝度勾配と移動量によって表されることを示している。なお、移動量(Δx,Δz)については上記式（７）のみでは決定できないため、検査領域Ｓ１内で移動量が一定と考え、最小二乗法を適用して算出している。

上記式（７）を用いて移動量を算出する際には、右辺の各注目点における移動前後の輝度差は一意にしか求まらない。そのため輝度勾配をどれだけ正確に算出するかが移動量の精度に直結する。輝度勾配の差分化では、一次精度風上差分を用いている。差分化において高次差分を適用すると、多くのデータが必要になり、ノイズが含まれていた際に影響を大きく受けてしまうためである。また、検査領域Ｓ１内の各点を基準とした高次差分では、検査領域Ｓ１外のデータを多く使用することとなり、検査領域Ｓ１そのものの移動量ではなくなってしまうという問題点も存在するからである。

輝度勾配を求める際に変形前の風上側の輝度勾配が移動することによって注目点の輝度差が生まれると考えることができるので、変形前は風上側の差分を適用する。ここでいう風上は、実際の移動方向ではなく、ピクセル移動量に対するサブピクセル移動量の向きのことであり、最大相関値ピークに放物線近似を施すことによって風上側を決定する。逆に、変形後の風下側の輝度勾配が逆に移動することによって注目点の輝度差が生じると考えることができるため、変形後は風下側の差分を適用する。

変形前の風上差分と変形後の風下差分を用いて２通りの解を求め、それらの平均をとった。さらに、実際には移動量が軸方向に沿わない場合には、変形前や変形後の輝度勾配が注目点と同一軸上に無く、ずれた位置の勾配を求める必要がある。つまり、Ｘ方向の勾配を求めたい場合には、Ｚ方向の移動も考慮して勾配を求めなければならない。そのため、輝度の内挿による輝度勾配の推定をすることで、精度向上を図っている。基本的には変形前（もしくは変形後）の位置を予測し、その位置での勾配を内挿により求める。

変形前（後）の注目点の位置は放物線近似を施した際のサブピクセル移動量により求める。その注目点位置が囲まれる８つの座標を用い、それらの比によって輝度勾配を算出する。具体的には、下記式（８）を用いる。そのようにして算出された輝度勾配と、輝度変化を用いて最小二乗法を適用し移動量を決定した。

（画像変形法）
上述した風上勾配法までは検査領域Ｓ１の形状は変更せず、正方形を保ったまま変形ベクトルの算出を行っている。しかし、現実には計測対象の変形に合わせて検査領域Ｓ１も変形していると考えられるため、検査領域Ｓ１の微小変形を考慮したアルゴリズムを導入し、変形ベクトル算出を高精度にて算出する必要がある。このため、本実施例ではサブピクセル精度での変形ベクトルの算出に画像変形法を導入している。すなわち、組織変形前の検査領域Ｓ１と組織変形後の伸縮及びせん断変形を考慮した検査領域Ｓ１とで相互相関を実施し、相関値ベースの反復計算によってサブピクセル変形量を決定している。なお、検査領域Ｓ１の伸縮及びせん断変形は線形で近似している。

画像変形法は一般的に、材料の表面ひずみ計測法に用いられ、ランダムパターンを塗布した材料表面を高空間分解能なカメラで撮影した画像に対して適用される。一方、ＯＣＴ断層像はスペックルノイズを多く含むだけでなく、特に生体組織では組織内の基質や水分の流動も伴い屈折率が変化するため、スペックルパターンに対する変形が大きい。このため、複雑組織にて局所的な変形が発生し、検査領域Ｓ１に膨張、収縮、剪断等の変形が伴う場合には解の収束解が著しく低下する。本手法における検査領域Ｓ１の縮小は、局所的な組織力学特性を検出するために必要不可欠である。そこで画像変形法では、風上勾配法で得られた変形量を収束計算の初期値として採用し、さらに、輝度分布の双三次関数補間によって検査領域Ｓ１を縮小した際でも低ロバスト性を実現している。なお、変形例においては、双三次関数補間以外の補間関数を用いてもよい。

より詳細には、以下の手順にて演算を実行する。まず、組織変形前のＯＣＴ断層像の輝度分布に双三次関数補間法を適用し、輝度分布の連続化を実施する。双三次関数補間法とは、ｓｉｎｃ関数を区分的に三次関数近似した畳み込み関数を用い、輝度情報の空間連続性を再現する手法である。本来は連続的な輝度分布を画像計測する際には光学系に依存した点広がり関数が畳み込まれるため、ｓｉｎｃ関数を用いた逆畳み込みを行うことにより、本来の連続的な輝度分布が復元される。離散的な一軸信号ｆ(x)の補間をする場合、畳み込み関数ｈ(x)は下記式（９）にて表される。

なお、ＯＣＴ計測条件の違いによって輝度補間関数ｈ(x)の形状も変更する必要がある。そこで輝度補間関数ｈ(x)のｘ＝１での微係数ａを可変とし、ａの値を変更することで輝度補間関数ｈ(x)の形状を変更可能なアルゴリズムとした。本実施例では、擬似ＯＣＴ断層像を用いた数値実験による検証結果を元にし、ａの値を決定した。以上のように画像補間をすることで、伸縮及びせん断変形を考慮した検査領域Ｓ１の各点にて、ＯＣＴ輝度値を求めることが可能となる。

伸縮及びせん断変形を考慮して算出した検査領域Ｓ１は、図４（Ｃ）に示すように、移動とともに変形を伴う。組織変形前のＯＣＴ断層像におけるある検査領域Ｓ１内の整数ピクセル位置での座標(x,z)が組織変形後に座標(x^＊,z^＊)に移動すると考えると、x^＊,z^＊の値は下記式（１０）にて表される。

ここで、Δx，Δzはそれぞれ検査領域Ｓ１中心から座標x，zまでの距離、u，vはそれぞれx，z方向への変形量、∂u／∂x，∂v／∂zはそれぞれx，z方向における検査領域Ｓ１の垂直方向への変形量、∂u／∂z，∂v／∂xはそれぞれx，ｚ方向における検査領域Ｓ１のせん断方向への変形量である。数値解法にはNewton-Raphson法を用い、６変数(u，v，∂u／∂x，∂u／∂z，∂v／∂x，∂v／∂z)での相関値微係数が０となるように、すなわち最大相関値を得るように反復計算を行う。なお、反復計算の収束性を高めるため、x，z方向の移動量初期値には風上勾配法で得られたサブピクセル移動量を用いる。相関値Ｒに対するヘッセ行列をＨ、相関値対するヤコビベクトルを▽Ｒとすると、１回の反復で得られる更新量ΔＰiは下記式（１１）にて表される。

収束の判定には、反復計算で随時得られる漸近解が収束解の近傍で十分小さくなることを用いる。しかし、スペックルパターンの変化が激しい領域においては、線形変形では追従できないために正しい収束解が得られない場合がある。その場合、本実施例では風上勾配法によって求めたサブピクセル移動量を採用している。

以上のようにして得られるサブピクセル精度の変形ベクトルを時間微分することにより、図５（Ｂ）に示すような変形速度ベクトル分布を算出することができる。そして、その変形速度ベクトル分布を空間微分することにより、ひずみ速度テンソルを算出することが可能となる。

（時空間移動最小二乗法）
ひずみ速度テンソルの算出には移動最小二乗法（Moving Least Square Method:以下「ＭＬＳＭ」という）を用いる。ＭＬＳＭは、移動量分布を平滑化すると共に、微係数の安定算出を可能とする手法である。ＭＬＳＭにおいて用いる二乗誤差式は下記式（１２）にて表される。

上記式（１２）において、Ｓ(x,z,t)を最小とするａ〜ｋのパラメータを求める。すなわち、近似関数として水平方向ｘ，奥行き方向ｚ，時間方向ｔの３変数２次多項式として下記式（１３）を採用する。そして、最小二乗近似に基づき、下記式（１４）から最適な微係数を算出して平滑化する。

この微係数を用いることにより、下記式（１５）に示すひずみ速度テンソルを算出することができる。ｆx，ｆzは各軸のひずみ増分を示し、その時間変化量からひずみ速度を算出している。

次に、上述した各手法を用いて行われる軟骨診断処理について説明する。図６および図７は、軟骨診断処理の処理手順を表す図である。
本実施例では、上述した各手法を軟骨診断装置１に適用する。図６に示すように、制御演算部６は、上述のようにして演算された変形速度ベクトル分布を順次記憶するとともに（図６（Ａ））、各変形速度ベクトル分布を空間微分して得られたひずみ速度テンソル分布を順次記憶する（図６（Ｂ））。なお、ひずみ速度テンソル分布は、図７（Ａ）に示すように、所定範囲ごとに識別可能に断層画像の態様で分布表示することができる。図中においてひずみ速度がマイナスの領域は圧縮領域を示し、プラスの領域は膨張領域を示す。

そして、断層画像における個々の座標についてのひずみ速度の時間変化を演算すると、図７（Ｂ）に示すような演算結果が得られる。同図の横軸は時間の経過を示し、縦軸はひずみ速度を示している。このように演算されたひずみ速度の時間変化に対して平滑化処理を行い、各座標についてひずみ速度の減衰係数を演算すると、図７（Ｃ）に示すような減衰係数分布を得ることができる。この減衰係数の断層分布は、後述の実験結果からも分かるように、軟骨組織の損傷度合いの分布に対応するものと評価することができる。そこで、制御演算部６は、この減衰係数の断層分布を軟骨組織の損傷度合いの分布と等価として、表示装置８に断層可視化させる。

本実施例では、減衰係数の断層分布そのものを軟骨組織の損傷度合い（変性度）として提示する。このため、減衰係数が大きい領域は損傷度（変性度）が大きい領域として示され、減衰係数が小さい領域は損傷度（変性度）が小さい領域として示される。なお、図中のハッチング部はひずみ速度が小さく、正負の値に振れる領域（精度が低い領域）を示し、軟骨診断においては無視することができる。

次に、減衰係数の断層分布が軟骨組織の損傷度合いの分布に対応することを検証するために行った実験結果について説明する。図８は、実験装置の構成を概略的に表す図である。図９は、実験に用いた試験片を表す図である。図１０は、実験に際して試験片に与えた荷重負荷状態を示すグラフである。なお、図８において、図１に示すものと同様の構成については同一の符号を付している。

図８に示すように、本実験装置は、図１に示した光学ユニット２に対して圧縮試験機７０を接続して構成される。圧縮試験機７０は、図１に示したプローブ４に代わり、試験片Ｗに対して圧縮荷重を負荷するものである。試験片Ｗとしては、ブタの膝関節の軟骨を用いている。圧縮試験機７０は、試験片Ｗに接触するガラス窓７２、ガラス窓７２との間に試験片Ｗを挟む金属圧子７４、金属圧子７４を軸線方向に駆動するためのリニアアクチュエータ７６、試験片Ｗに負荷される荷重を検出するロードセル７８等を備えている。

本実験では、制御演算部６を２台のパーソナルコンピュータ８０，８２にて構成し、その一方により光学ユニット２の駆動およびＯＣＴの演算処理を実行し、他方により圧縮試験機７０を駆動した。オブジェクトアーム１２には、試験用プローブ８４が接続されている。試験用プローブ８４は、図１に示したプローブ４に代わり、物体光を試験片Ｗに向けて走査可能な光学機構を備えている。

図９（Ａ）は圧縮試験機７０への試験片Ｗの設置態様を示し、図９（Ｂ）は試験片Ｗの側面視およびその大きさを示し、図９（Ｃ）は試験片Ｗの正面視およびその大きさを示している。試験片Ｗとして、生後約６ヶ月のブタの膝関節大腿骨外側顆の軟骨組織サンプルを直径５．５ｍｍの円筒形状にて採取した。また、変性軟骨の力学特性を診断評価するため、酵素処理により模擬変性軟骨を作製し、正常軟骨との比較を行った。酵素処理には、関節軟骨の弾性を保っているコラーゲン線維を分解するcollagenase type2(C6885,Sigma Ardrich)を用いた。ｐＨ７．４のリン酸緩衝生理食塩水(10010,Invitrogen)にcollagenaseを溶解させて作製した３０[unit/ml]のcollagenase溶液に正常関節軟骨を入れ，３７［℃］の状態下で０．５時間，１時間，２時間の時間設定で処理を施し，コラーゲン線維変性をＯＡ変性度として模擬している。

図８に戻り、本実験では、試験片Ｗの下骨側を金属圧子７４に固定し、圧縮試験機７０内に設置したガラス窓７２に押し当て、軟骨組織に圧縮変形を発生させた。ガラス窓７２を介してＯＣＴビームを打ち込むことで応力緩和試験中のＯＣＴ断層画像を２７．０３［flame/sec］にて連続取得した。この圧縮実験では、圧縮量を１０［％］ひずみ、圧縮速度を０．１［％/sec］として固定した。圧縮荷重負荷時間は１００［sec］となり、１０［％］ひずみに達した後は応力緩和時の軟骨組織挙動の観察を行った。また、再帰的相互相関法の設定は検査領域を９［pixel］から５［pixel］に縮小する２階層の設定とし、６．５×１５．６２［μｍ］毎のベクトル密度にてサブピクセル精度の変形速度ベクトル分布を算出した。

図１０は、この応力緩和試験においてロードセル７８にて検出された応力時系列データを示す。同図の横軸は荷重負荷開始からの時間［sec］を示し、縦軸は圧縮応力[MPa]を示している。図示のように、圧縮中（０〜１００秒）では応力が増加し、圧縮終了（１００秒）以降は応力緩和が発生している。そして、１２００秒付近で平衡状態に達することが観察された。なお、このように応力が指数関数的に減衰していることから、軟骨はその緩和弾性率が時間変化する粘弾性体であることが分かる。

図１１および図１２は、応力緩和試験によるひずみ速度分布の時間変化を断層可視化した図である。図１１は正常軟骨についての演算結果を示し、（Ａ）は荷重負荷開始から０．３１[sec]、（Ｂ）は５０．３５[sec]、（Ｃ）は１００．２９[sec]、（Ｄ）は１０１．９６[sec]、（Ｅ）は１１０．９５[sec]、（Ｆ）は１２０．１２[sec]、（Ｇ）は１３０．５２[sec]、（Ｈ）は１４０．５１[sec]、（Ｉ）は１５０．４９[sec]の結果を示している。

一方、図１２は模擬変性軟骨についての演算結果を示し、（Ａ）は荷重負荷開始から１．０７[sec]、（Ｂ）は５０．０７[sec]、（Ｃ）は１００．９９[sec]、（Ｄ）は１０１．９９[sec]、（Ｅ）は１１０．９８[sec]、（Ｆ）は１２０．９７[sec]、（Ｇ）は１３０．９６[sec]、（Ｈ）は１４０．９５[sec]、（Ｉ）は１５０．９４[sec]の結果を示している。すなわち、各図の（Ａ）および（Ｂ）は圧縮中を示し、（Ｃ）〜（Ｉ）は応力緩和中を示している（図１０参照）。

図１１および図１２から、試験片Ｗの表層において大きな圧縮ひずみ速度が発生していることが確認できる。この領域は軟骨の表層領域（表面から１５％）に対応し、変形速度ベクトルが急激に変化する領域と対応している。一方、軟骨表面から１５〜７０％は中間層領域と呼ばれ、検出されたひずみ速度は小さく、変形速度ベクトルの急激な変化が見られない領域と対応している。コラーゲン線維の配向特性は、表層と中間層において大きく異なっていることが知られている。特に、表層では軟骨表面に水平配向しているため圧縮方向に対する剛性が低く、圧縮ひずみ速度の集中が発生したと考えられる。

また、図１１と図１２とを対比すると、正常軟骨では１５０．４９[sec]でもひずみ速度の残存が確認できるのに対し、模擬変性軟骨では１３０．９６[sec]でほとんどひずみ速度が見られないことがわかる．したがって、緩和時間に着目し比較することで関節軟骨の変性の度合いを評価することが可能であることが分かる。このことから、ひずみ速度の経時変化から算出される減衰係数の比較によって軟骨変性度を評価できると考えられる。本実施例ではこの点に着目し、図７に示したように、ひずみ速度テンソル分布から減衰係数の分布を算出している。すなわち、ひずみ速度の減衰係数の断層分布を、軟骨組織の損傷度合いの分布と等価として画面に断層可視化することで、軟骨診断が可能となるようにした。

図１３は、正常軟骨と模擬変性軟骨の減衰係数の比較を表す図である。図中の横軸は酵素処理時間（Ｈ：hour）を示し、縦軸が減衰係数を示す。酵素処理時間として０Ｈ，０．５Ｈ，１Ｈ，２Ｈの結果が示されている。なお、０Ｈは酵素処理時間がゼロ、つまり正常軟骨を示している。図示の結果は、表層領域部分における減衰係数分布を空間平均したものであり、応力緩和開始から５０秒までの緩和挙動に対して対数化処理が施されている。

図示の結果より、酵素処理時間が長くなるにつれて、減衰係数が大きくなっていることが分かる。コラゲナーゼ酵素処理では、弾性を支配するコラーゲン線維の分解に伴い、粘性を支配するプロテオグリカンの離脱も発生する。このため、コラーゲン線維分解による多孔度の増大に加え、プロテオグリカン離脱による保水特性の低下によって、透水特性が著しく上昇し粘性特性が損なわれると考えられる。透水特性の増大は粘性による荷重支持機構の破綻を意味することから、図示のように緩和弾性係数の減衰係数が低下すると考えられる。そして図示のように、正常軟骨に対していずれの酵素処理時間の結果においても有意差があり、変性の検出が可能であった。本実験で用いた模擬変性軟骨試料は正常軟骨と比べ、外見上は変化がなくレントゲンやＣＴでの目視診断では判断できない程度の変性である。したがって、従来診断できなかった初期の変形性膝関節症の診断が、本手法により可能となることが分かる。

次に、制御演算部６が実行する具体的処理の流れについて説明する。
図１４は、制御演算部６により実行される軟骨診断処理の流れを示すフローチャートである。なお、この軟骨診断処理は、医師等によりプローブ４の先端部が患者の膝関節に挿入された状態にて開始される。制御演算部６は、駆動部３８を駆動して図１０に示した応力負荷および応力緩和処理を実行する一方で、図１４に示す処理を実行する。

すなわち、制御演算部６は、連続撮影された時刻の異なる前後２枚のＯＣＴ断層画像Ｉ(x,z,t)とＩ(x,z,t+Δt)を読み込む（Ｓ１０）。続いて、再帰的相互相関法による処理を実行する。ここではまず、最小解像度（最大サイズの検査領域）での相互相関処理を実行し、相関係数分布を求める（Ｓ１２）。続いて、隣接相互相関乗法により、隣接する相関係数分布の積を演算する（Ｓ１４）。このとき、標準偏差フィルタ等の空間フィルタにより過誤ベクトルの除去をし（Ｓ１６）、最小二乗法等により除去ベクトルの内挿補間を実行する（Ｓ１８）。続いて、検査領域を小さくすることによって解像度を上げて相互相関処理を継続する（Ｓ２０）。すなわち、低解像度での参照ベクトルを基に相互相関処理を実行する。このときの解像度が予め定める最高解像度でなければ（Ｓ２２のＮ）、Ｓ１４に戻る。

そして、Ｓ１４〜Ｓ２０の処理を繰り返し、最高解像度での相互相関処理が完了すると（Ｓ２２のＹ）、サブピクセル解析を実行する。すなわち、最高解像度（最小サイズの検査領域）での変形ベクトルの分布に基づき、風上勾配法によるサブピクセル移動量を演算する（Ｓ２４）。そして、このとき算出されたサブピクセル移動量に基づき、画像変形法によるサブピクセル変形量を演算する（Ｓ２６）。続いて、最大相互相関値によるフィルタ処理により過誤ベクトルの除去をし（Ｓ２８）、最小二乗法等により除去ベクトルの内挿補間を実行する（Ｓ３０）。そして、このようにして得られた変形ベクトルの時間微分を行い、その断層画像について変形速度ベクトルの断層分布Ｕ(x,z,t)，Ｗ(x,z,t)を算出し、Ｓ１０に戻る。ここで、Ｕ(x,z,t)は変形速度ベクトル分布のｚ方向成分であり、Ｗ(x,z,t)は変形速度ベクトル分布のｘ方向成分である。以上の処理を後続の断層画像についても実行する（Ｓ３４のＮ）。そして、設定数の断層画像について変形ベクトルの断層分布が求まると（Ｓ３４のＹ）、軟骨診断のための変性度合い演算提示処理を実行する（Ｓ３６）。

図１５は、図１４におけるＳ３６の変性度合い演算提示処理を詳細に示すフローチャートである。変性度合い演算提示処理において、制御演算部６は、時空間移動最小二乗法により、Ｓ３４にて算出した変形速度ベクトル分布の平滑化を実行する（Ｓ４０）。そして、その平滑化された変形速度ベクトルに対して空間微分をすることにより、ひずみ速度テンソルを算出する（Ｓ４２）。そして、ひずみ速度の時系列データについて対数変換を行い（Ｓ４４）、最小二乗法によりひずみ速度の減衰係数を演算する（Ｓ４６）。このようにして得られた減衰係数の分布を軟骨組織の損傷度合いと等価とみなして表示装置８に断層可視化する（Ｓ４８）。すなわち、例えば図７（Ｃ）に示した態様にて軟骨の変性度を提示する。

以上に説明したように、本実施例によれば、応力緩和によるひずみ速度の変化（減衰係数）と、軟骨組織の損傷度合いとの間に対応関係があることを利用した演算処理により、診断対象である軟骨の損傷度合いが断層可視化される。このため、医師等がその断層可視化された画像を確認することにより、軟骨診断を容易に行うことが可能となる。すなわち、ＯＣＴを用いた軟骨診断をより実用に供し得るものとすることができる。

［第２実施例］
次に、本発明の第２実施例について説明する。本実施例は、応力負荷方法として応力緩和法ではなく、動的粘弾性法を適用した点が異なる以外は第１実施例とほぼ同様である。このため、上記第１実施例と同様の構成部分については必要に応じて同一の符号を付す等して適宜その説明を省略する。

本実施例では、図１および図２に示したプローブ４により所定の振動荷重（動的力、動的荷重）を発生させる。すなわち、プローブ４への通電により、圧電素子５２が所定の加振周波数にて振動し、軟骨Ｊに対して振動荷重を付与する。それにより、軟骨Ｊの動的粘弾性を評価し、その損傷度合い（変性度）を断層可視化する。

図１６は、第２実施例に係る変性度合い演算提示処理を詳細に示すフローチャートである。この処理は、第１実施例の図１５に示した処理に置き換えて実行される。制御演算部６は、駆動部３８を駆動して軟骨Ｊに動的荷重（動的力、振動荷重）を付与する一方で、図１４および図１６に示す処理を実行する。

なお、本実施例では、図１４に示すＳ２４の処理において、上記式（７）に代えて、下記式（１６）を用いた輝度分布の定式化を行う。
δf/δt＝δf/δx・Δx+δf/δy・Δy+δf/δz・Δz ・・・（１６）

そして図１６に示すように、変性度合い演算提示処理において、制御演算部６は、時空間移動最小二乗法により、Ｓ３４にて算出した変形速度ベクトル分布の平滑化を実行する（Ｓ２４０）。そして、その平滑化された変形速度ベクトルに対してフーリエ変換を行う（Ｓ２４２）。それにより、変形速度ベクトルにおいてプローブ４の加振周波数ωに同期する成分のみを抽出するバンドパス・フィルタ処理を実行する（Ｓ２４４）。その後、フーリエ逆変換により、加振周波数成分の変形速度ベクトルＵ_ω(x,z,t)，Ｗ_ω(x,z,t)の断層分布を演算する（Ｓ２４６）。ここで、Ｕ_ω(x,z,t)は変形速度ベクトルのｚ方向成分であり、Ｗ_ω(x,z,t)は変形速度ベクトルのｘ方向成分である。Ｗ_ω(x,z,t)は、例えば下記式（１７）の形で表すことができる。
Ｗ_ω(x,z,t)＝Ａ_ω(x,z)sin(ωｔ+Φ_ω(x,z)) ・・・（１７）

続いて、圧電素子５２の変位情報から、加振周波数成分に係る変形速度ベクトルの規格化振幅Ａ'_ω(x,z)と位相差Φ_ω(x,z)の断層分布を算出する（Ｓ２４８）。ここで、規格化振幅Ａ'_ω(x,z)は、振幅Ａ_ω(x,z)を圧電素子５２の変位情報Ｗ_ｐｚｔ(t)＝Ａ_ｐｚｔsin(ωｔ)の振幅Ａ_ｐｚｔにて除算することにより得られる。なお、変形例においては、規格化振幅Ａ'_ω(x,z)および位相差Φ_ω(x,z)を最小二乗法により直接算出してもよい。

続いて、規格化振幅Ａ'_ω(x,z)および位相差Φ_ω(x,z)の少なくとも一方の断層分布に基づいて軟骨の動的粘弾性を演算する（Ｓ２５０）。そして、演算された動的粘弾性に基づいて軟骨組織の損傷度合いを断層可視化する（Ｓ４８）。

なお、この動的粘弾性の評価法については種々考えられる。規格化振幅Ａ'_ω(x,z)の断層分布に基づいて粘弾性評価を行う場合、以下のように評価できる。すなわち、規格化振幅Ａ'_ω(x,z)が大きい箇所は動的粘弾性率が小さくなる。すなわち、規格化振幅Ａ'_ω(x,z)は、軟骨において奥行き方向の剛性が低い表層では相対的に大きくなり、中層では相対的に小さくなる。軟骨の変性が進行すると、表層の奥行き方向の剛性がさらに低くなるため、規格化振幅Ａ'_ω(x,z)はさらに大きくなると考えられる。

位相差Φ_ω(x,z)の断層分布に基づいて粘弾性評価を行う場合には、位相接続処理（アンラッピング処理）を施してもよい。この場合、以下のように評価できる。すなわち、奥行き方向の位相差Φ_ω(x,z)は、軟骨の表面ではゼロであり、表層ではその表層の粘弾性効果により深部に向かうにしたがって徐々に大きくなる。つまり、位相が遅れていく。変性軟骨では、表層の粘性特性が損なわれるため、奥行き方向の位相差Φ_ω(x,z)が小さくなる（つまり弾性挙動を示すようになる）。

位相差Φ_ω(x,z)の空間微分（∂Φ_ω／∂z）の大きさにより変性度を診断することもできる。∂Φ_ω／∂zが大きい箇所は変性度が低く（つまり正常）、∂Φ_ω／∂zが小さい箇所は変性度が高いといえる（つまり粘性効果が損なわれている）。このような知見のもと、軟骨組織の損傷度合いを断層可視化することができる。

以上、本発明の好適な実施例について説明したが、本発明はその特定の実施例に限定されるものではなく、本発明の技術思想の範囲内で種々の変形が可能であることはいうまでもない。

上記第１実施例では、ひずみ速度の減衰係数を演算し、軟骨の損傷度合いと等価として断層可視化する例を示した。変形例においては、変形速度の減衰係数を演算し、軟骨の損傷度合いと等価として断層可視化してもよい。

上記第２実施例では、規格化振幅Ａ'_ω(x,z)と位相差Φ_ω(x,z)から粘弾性評価を行う例を示したが、これら規格化振幅Ａ'_ω(x,z)および位相差Φ_ω(x,z)から粘弾性物理量を表現する値を新たに定義し、それに基づいて粘弾性評価を行ってもよい。

上記第２実施例では、図１６に示した変性度合い演算提示処理において、変形速度ベクトルの断層分布に基づいて軟骨の粘弾性を評価することとした。変形例においては、ひずみ速度テンソルの断層分布を演算し、それに基づいて同様に軟骨の粘弾性を評価し、断層可視化してもよい。このひずみ速度テンソルは、変形速度ベクトルを空間微分することにより算出することができる。

上記第２実施例では明示しなかったが、プローブ４にロードセルなどのセンサを取り付け、軟骨に付与する荷重を計測できる構成としてもよい。これにより、軟骨組織の複素弾性率（貯蔵弾性率、損失弾性率、複素粘性率等）を算出することができ、規格化振幅Ａ'_ω(x,z)と位相差Φ_ω(x,z)から粘弾性評価をするための物理量を定義し易くなる。

上記実施例では述べなかったが、「力学特徴量」として、ひずみ速度の振幅値やひずみ速度の時間遅れ（位相遅れ）を採用してもよい。すなわち、上述した動的粘弾性法を例に挙げると、応力の負荷および緩和を繰り返す過程でひずみ速度が正の値と負の値を行き来するように変動する。この負荷の変動に追従するように、ひずみ速度も変動するようになる。この点につき、その負荷変動に対する、ひずみ速度の変動の大きさ（振幅値）とひずみ速度変動の追従性（応答性）が、軟骨の損傷度合いに対応して変化する傾向がある。具体的には、軟骨の損傷度合いが進むほど、ひずみ速度の振幅値は大きくなり、時間遅れ（位相遅れ）は小さくなる傾向にある。そこで、そのひずみ速度の振幅値や時間遅れ（位相遅れ）について断層分布を演算してもよい。その振幅や時間遅れ（位相遅れ）の断層分布を、軟骨組織の損傷度合いの分布と等価として断層可視化させてもよい。

あるいは、「力学特徴量」として、ひずみ速度の周期的な変動の中央値を採用してもよい。すなわち、上述した動的粘弾性法を例に挙げると、応力の負荷および緩和を繰り返す過程でひずみ速度が正の値と負の値を行き来するように変動する。そのひずみ速度の変動の中心は、粘性力と弾性力との釣り合いから、ゼロからややずれる傾向がある。そして、そのずれ量が、軟骨の損傷度合いに対応して変化する傾向がある。そこで、ひずみ速度の変動中心（中央値）のゼロからのずれ量について断層分布を演算してもよい。そのずれ量の断層分布を、軟骨組織の損傷度合いの分布と等価として断層可視化させてもよい。

上記実施例では、軟骨に所定の変形エネルギー（負荷）を付与するための負荷装置として、圧電素子等の接触により応力を付与する荷重機構を例示した。変形例においては、超音波（音圧）、光音響波、電磁波等によって非接触にて軟骨に負荷（加振力）を付与する負荷装置を採用してもよい。

なお、本発明は上記実施例や変形例に限定されるものではなく、要旨を逸脱しない範囲で構成要素を変形して具体化することができる。上記実施例や変形例に開示されている複数の構成要素を適宜組み合わせることにより種々の発明を形成してもよい。また、上記実施例や変形例に示される全構成要素からいくつかの構成要素を削除してもよい。

【０００７】
Ｔの光学系を含む光学ユニットに接続されることにより、所定の対象部位の診断を可能とする診断用プローブとして構成されてもよい。このプローブは、対象部位に当接可能に構成された先端部と、光学ユニットからの光を対象部位に導いて走査させるための光学機構と、対象部位に対して所定の変形エネルギー（応力）を負荷するための荷重機構（負荷装置）とを備える。
［００２４］
また、上記技術を利用した軟骨診断方法を構築してもよい。この方法は、軟骨に対して所定の変形エネルギー（荷重）を負荷するステップと、変形エネルギーの付与（荷重の負荷）に応じた軟骨の変形度合いを光コヒーレンストモグラフィーにより断層画像として表示させるステップと、断層画像に基づいて軟骨組織の損傷度合いを診断するステップとを備えるものでよい。
［００２５］
以下、図面を参照しつつ本発明を具体化した実施例について詳細に説明する。
［第１実施例］
図１は、第１実施例に係る軟骨診断装置の構成を概略的に表す図である。図２は、軟骨診断装置を構成する診断プローブの構成を概略的に表す図である。本実施例の軟骨診断装置は、診断対象である関節軟骨に所定の応力を負荷し、その応力に対する軟骨組織の変形度合いを断層可視化することにより、軟骨組織の損傷度合い（変性度）を診断可能とするものである。この断層可視化にＯＣＴを利用する。
［００２６］
図１に示すように、軟骨診断装置１は、ＯＣＴを用いる光学系を含む光学ユニット２、光学ユニット２に接続される診断用のプローブ４、ＯＣＴにより得られた光干渉データに基づいて軟骨組織の損傷度合いを演算する制御演算部６、および軟骨組織の損傷度合いを断層可視化する態様で表示する表示装置８を備える。図示の例では、光学ユニット２としてマッハツェンダー干渉計をベースとした光学系が示されているが、マイケルソン干渉計その他の光学系を採用することもできる。
［００２７］
本実施例では、ＯＣＴとして波長走査型レーザを光源とするＳＳ−ＯＣＴ（Ｓｗｅｐｔ Ｓｏｕｒｃｅ ＯＣＴ）を用いるが、ＴＤ−ＯＣＴ（Ｔｉｍｅ Ｄｏｍａｉｎ ＯＣＴ）、ＳＤ

Claims (8)

1. 関節軟骨を診断するための軟骨診断装置であって、
   光コヒーレンストモグラフィーを用いる光学系を含む光学ユニットと、
   前記光学ユニットに接続される一方、先端部が関節腔に挿入可能に構成され、前記光学ユニットからの光を軟骨に導いて走査させるための光学機構と、軟骨に対して所定の変形エネルギーを付与するための負荷装置とを含むプローブと、
   前記負荷装置および前記光学機構の駆動を制御し、それらの駆動に基づいて前記光学ユニットから出力された光干渉信号を処理し、変形エネルギーの付与による軟骨内部の変形に伴う所定の力学特徴量の変化をその軟骨の断層位置に対応づけて演算し、その力学特徴量の変化に基づいて軟骨組織の損傷度合いを演算する制御演算部と、
   前記軟骨組織の損傷度合いを断層可視化する態様で表示する表示装置と、
   を備えることを特徴とする軟骨診断装置。
2. 前記制御演算部は、
   前記負荷装置の駆動により軟骨に所定の荷重を負荷した後に応力緩和させ、
   その応力緩和中に前記光学ユニットにより撮影された前後の断層画像データに基づき、軟骨の断層位置に対応した変形速度ベクトルを前記力学特徴量として演算し、その変形速度ベクトルの変化から得られる減衰係数の断層分布を、前記軟骨組織の損傷度合いの分布と等価として断層可視化させることを特徴とする請求項１に記載の軟骨診断装置。
3. 前記制御演算部は、前記力学特徴量として演算した変形速度ベクトルを空間微分することによりひずみ速度テンソルを演算し、そのひずみ速度の減衰係数の断層分布を前記軟骨組織の損傷度合いの分布と等価として断層可視化させることを特徴とする請求項２に記載の軟骨診断装置。
4. 前記制御演算部は、順次取得する断層画像データの相互相関、輝度勾配の変位および輝度パターンの変形に基づき前記変形速度ベクトルをサブピクセル精度にて演算することを特徴とする請求項２または３に記載の軟骨診断装置。
5. 前記制御演算部は、
   前記負荷装置の駆動により軟骨に所定の動的力を負荷し、
   その動的力による軟骨の振動中に前記光学ユニットにより撮影された前後の断層画像データに基づき、軟骨の各断層位置における変形速度ベクトルを前記力学特徴量として演算し、その変形速度ベクトルのうち加振周波数成分である特定変形速度ベクトルの振幅および位相差の少なくとも一方に基づき軟骨の粘弾性の断層分布を演算し、その粘弾性の断層分布を前記軟骨組織の損傷度合いの分布と等価として断層可視化させることを特徴とする請求項１に記載の軟骨診断装置。
6. 前記制御演算部は、
   前記負荷装置の駆動により軟骨に所定の動的力を負荷し、
   その動的力による軟骨の振動中に前記光学ユニットにより撮影された前後の断層画像データに基づき、軟骨の各断層位置におけるひずみ速度テンソルを前記力学特徴量として演算し、そのひずみ速度テンソルのうち加振周波数成分である特定ひずみ速度テンソルの振幅および位相差の少なくとも一方に基づき軟骨の粘弾性の断層分布を演算し、その粘弾性の断層分布を前記軟骨組織の損傷度合いの分布と等価として断層可視化させることを特徴とする請求項１に記載の軟骨診断装置。
7. 前記負荷装置が、軟骨に付与される力又は荷重を計測可能なセンサを含み、
   前記制御演算部は、
   演算された前記力学特徴量の加振周波数成分の振幅および位相差の少なくとも一方と、前記センサによる計測値に基づき軟骨の粘弾性の断層分布を演算し、その粘弾性の断層分布を前記軟骨組織の損傷度合いの分布と等価として断層可視化させることを特徴とする請求項５または６に記載の軟骨診断装置。
8. 光コヒーレンストモグラフィーを用いる光学系を含む光学ユニットに接続されることにより、対象部位の診断を可能とする診断用プローブであって、
   前記対象部位に当接可能に構成された先端部と、
   前記光学ユニットからの光を軟骨に導いて走査させるための光学機構と、
   前記対象部位に対して所定の変形エネルギーを付与するための負荷装置と、
   を備えることを特徴とする診断用プローブ。

Images