JP4940340B2 - 椎骨セグメンテーション装置、方法及びプログラム - Google Patents

Description

この発明は、椎骨が含まれる三次元医用画像から各椎骨を分離・認識する椎骨セグメンテーション装置、方法及びプログラムに関する。

従来から、医療分野において、被写体を胸部から脚部にかけて複数の切断位置でスキャンして得られた断層画像に基づいて、複数の椎骨をそれぞれ分離・認識する画像処理（以下、この画像処理を「椎骨セグメンテーション」という。）のアルゴリズムが種々提案されている。

非特許文献１には、脊椎中心線に沿った平均濃度プロファイルに基づいて各椎間板領域の基準位置を推定した後、該基準位置近傍における断層画像に対して適切なモーフォロジ演算を行うことにより、椎体終板と椎間とを抽出し、各椎骨を分離・認識する方法が開示されている。

非特許文献２には、脊椎の中心線に沿うサジタル・コロナルのＣＲＰ（Ｃｕｒｖｅｄ Ｐｌａｎａｒ Ｒｅｆｏｒｍａｔｉｏｎ）画像を生成し、該中心線を法線とする断層画像の濃度値積分プロファイルに基づいて椎間板を抽出し、各椎骨を分離・認識する方法が開示されている。

しかしながら、スキャン方向（体軸方向）のスライス厚が大きい場合、パーシャルボリューム効果の影響から体軸方向における空間分解能が不足し、画像表現能が低下することが一般的に知られている。各椎骨の長さが比較的短い頚椎から胸椎上部までの範囲において、椎間板と椎体終板とを分離して断層画像上で表現することは特に困難である。その結果、椎間板の検出能（すなわち、各椎骨を分離・認識する位置精度）が低下する問題があった。

そこで、本願発明者は、特許文献１に記載されているように、椎骨の中心線を検出した後、該中心線近傍における体軸方向の濃度勾配を表す特徴量を算出し、該特徴量に基づいて椎体の中心位置を特定した後、椎間板の位置を推定する方法に関する発明を行った。椎間板よりも検出性が高い海綿体領域（椎骨の中空領域）に基づいて椎間板の位置を推定するので、各椎骨を分離・認識する位置精度をさらに向上できる。

また、本願発明者は、特許文献２に記載されているように、複数種類（例えば、サジタル断面及びコロナル断面の組合せ）の断面画像を並べて表示する場合、所定のマーク画像を各断面画像上の各対応箇所にそれぞれ付加する発明を行った。これにより、病巣等の位置を正確に把握できる。

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

"マルチスライスＣＴ画像を用いた骨粗しょう症診断支援アルゴリズムの構築"、電信情報通信学会技術報告ＭＩ１０６、ｐｐ．２５−２８、２００７年１月 "Ａｕｔｏｍａｔｅｄ Ｓｐｉｎｅ Ｃｏｌｕｍｎ Ｅｘｔｒａｃｔｉｏｎ ａｎｄ Ｐａｒｔｉｔｉｏｎ"、Ｉｎ ３ｒｄ ＩＥＥＥ Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｓｙｍｐｏｓｉｕｍ ｏｎ Ｂｉｏｍｅｄｉｃａｌ Ｉｍａｇｉｎｇ、ｐｐ．３９０−３９３、２００６

本発明は前記特許文献１及び２に開示されている方法に関連してなされたものであって、脊椎の形状が全体的又は局所的に変形している被写体の三次元医用画像に対しても、各椎骨を適切に分離・認識できる椎骨セグメンテーション装置、方法及びプログラムを提供する。

本発明に係る椎骨セグメンテーション装置は、各椎骨の中心軸に沿った断層画像を生成する画像生成部と、前記各椎骨についての前記中心軸方向に交叉する断面形状の鮮明度を表す第１特徴量を、前記画像生成部により生成された前記断層画像に基づいて該中心軸上の点毎に算出する第１特徴量算出部と、前記各椎骨についての前記中心軸方向に平行する断面形状の鮮明度を表す第２特徴量を、前記画像生成部により生成された前記断層画像に基づいて該中心軸上の点毎に算出する第２特徴量算出部と、前記各椎骨の配列の規則性を表す第３特徴量を、前記第１及び第２特徴量算出部により算出された前記第１及び第２特徴量に基づいて前記中心軸上の点毎に算出する第３特徴量算出部と、前記第３特徴量算出部により前記中心軸上の点毎に算出された前記第３特徴量に基づいて、該中心軸上における前記各椎骨の位置を推定する椎骨位置推定部とを有することを特徴とする。

このように、各椎骨の中心軸方向に交叉及び平行する断面の形状的特徴を考慮した第３特徴量を算出するようにしたので、異なる２方向からの断面形状に基づいた定量的な濃度解析が可能であり、脊椎の形状が全体的又は局所的に変形している被写体の三次元医用画像に対しても、各椎骨を適切に分離・認識できる。脊椎湾曲症等により全体的に脊椎の形状が変形している場合、特に効果がある。

また、前記第１特徴量は、前記中心軸上の各点を中心とする円環状の模様を抽出する特徴量であることが好ましい。これにより、略円筒形状の皮質を有する各椎骨の中心軸方向に交差する断面形状（すなわち、円環状の模様）を適切に定量化することができる。

さらに、前記第２特徴量は、前記中心軸方向に延在する管状の模様を抽出する特徴量であることが好ましい。これにより、略円筒形状の皮質を有する各椎骨の中心軸方向に平行する断面形状（すなわち、管状の模様）を適切に定量化することができる。

さらに、前記中心軸上の各点毎に算出された前記第３特徴量を所定の周期関数又は準周期関数に適合させる関数適合部をさらに有し、前記椎骨位置推定部は、前記関数適合部により適合させた前記所定の周期関数又は準周期関数に基づいて前記中心軸上における前記各椎骨の位置を推定することが好ましい。これにより、各椎骨の配列の規則性を表す第３特徴量を所定の規則的な関数形に適合させるので、前記第３特徴量が備える高空間周波数のノイズを除去可能であり、各椎骨の位置を高精度に推定できる。

さらに、前記所定の周期関数又は準周期関数には位相が変数として含まれ、動的計画法を用いて前記中心軸上の各点毎に前記位相を順次決定する位相決定部をさらに有することが好ましい。これにより、関数への適合誤差を低減可能であり、各椎骨の位置をさらに高精度に推定できる。骨折等により局所的に脊椎の形状が変形している場合、特に効果がある。

さらに、前記椎骨位置推定部により推定された前記中心軸上の位置に応じた前記各椎骨に対し、予め定義された識別情報をそれぞれ付与する識別情報付与部をさらに有することが好ましい。

さらに、前記識別情報付与部は、第１の時点で生成された前記断層画像に基づいて付与された前記識別情報と同一の識別情報を、第２の時点で生成された前記断層画像に基づいて推定された前記中心軸上の位置に応じた前記各椎骨に対してそれぞれ付与することが好ましい。

さらに、前記第１の時点で生成された前記断層画像に基づいて算出された前記第３特徴量の特性が、前記第２の時点で生成された前記断層画像に基づいて算出された前記第３特徴量の特性と合致するか否かを判別する特徴量判別部をさらに有し、前記特徴量判別部により前記第３特徴量の特性が合致すると判別された場合、識別情報付与部は前記同一の識別情報を付与することが好ましい。

さらに、前記識別情報は、解剖学上の椎骨の種別であることが好ましい。

本発明に係る椎骨セグメンテーション方法は、各椎骨の中心軸に沿った断層画像を生成する生成ステップと、前記各椎骨についての前記中心軸方向に交叉する断面形状の鮮明度を表す第１特徴量を、生成された前記断層画像に基づいて該中心軸上の点毎に算出する第１特徴量算出ステップと、前記各椎骨についての前記中心軸方向に平行する断面形状の鮮明度を表す第２特徴量を、生成された前記断層画像に基づいて該中心軸上の点毎に算出する第２特徴量算出ステップと、前記各椎骨の配列の規則性を表す第３特徴量を、算出された前記第１及び第２特徴量に基づいて前記中心軸上の点毎に算出する第３特徴量算出ステップと、前記中心軸上の点毎に算出された前記第３特徴量に基づいて前記中心軸上における前記各椎骨の位置を推定する椎骨位置推定ステップとを備えることを特徴とする。

本発明に係るプログラムは、コンピュータを、各椎骨の中心軸に沿った断層画像を生成する画像生成部、前記各椎骨についての前記中心軸方向に交叉する断面形状の鮮明度を表す第１特徴量を、生成された前記断層画像に基づいて該中心軸上の点毎に算出する第１特徴量算出部、前記各椎骨についての前記中心軸方向に平行する断面形状の鮮明度を表す第２特徴量を、生成された前記断層画像に基づいて該中心軸上の点毎に算出する第２特徴量算出部、前記各椎骨の配列の規則性を表す第３特徴量を、算出された前記第１及び第２特徴量に基づいて前記中心軸上の点毎に算出する第３特徴量算出部、前記中心軸上の点毎に算出された前記第３特徴量に基づいて前記中心軸上における前記各椎骨の位置を推定する椎骨位置推定部として機能させることを特徴とする。

本発明に係る椎骨セグメンテーション装置、方法及びプログラムによれば、各椎骨の中心軸に沿った断層画像を生成し、前記各椎骨についての前記中心軸方向に交叉する断面形状の鮮明度を表す第１特徴量を、生成された前記断層画像に基づいて該中心軸上の点毎に算出し、前記各椎骨についての前記中心軸方向に平行する断面形状の鮮明度を表す第２特徴量を、生成された前記断層画像に基づいて該中心軸上の点毎に算出し、前記各椎骨の配列の規則性を表す第３特徴量を、算出された前記第１及び第２特徴量に基づいて前記中心軸上の点毎に算出し、前記中心軸上の点毎に算出された前記第３特徴量に基づいて前記中心軸上における前記各椎骨の位置を推定するようにしたので、異なる２方向からの断面形状に基づいた定量的な濃度解析が可能であり、脊椎の形状が全体的又は局所的に変形している被写体の三次元医用画像に対しても、各椎骨を適切に分離・認識できる。脊椎湾曲症等により全体的に脊椎の形状が変形している場合、特に効果がある。

本実施の形態に係る椎骨セグメンテーション装置の電気的な概略ブロック図である。 椎骨中心線検出部により各椎骨の中心線が検出された椎体の模式図である。 椎骨セグメンテーション装置の動作を示すフローチャート図である。 図４Ａ及び図４Ｂは、各断面位置における第１特徴量（被積分関数）の二次元分布を表す模式図である。 椎体の三次元形状と、その内部点におけるヘッセ行列の固有ベクトルとの対応関係を示す模式図である。 図６Ａは、Ｚ’軸方向における第１及び第２特徴量のプロファイルである。図６Ｂは、第２特徴量にＬＯＧフィルタを作用させた後のプロファイルである。図６Ｃは、Ｚ’軸方向における第３特徴量のプロファイルである。 第３特徴量を所定の関数に適合させる方法を示すフローチャート図である。 第３特徴量の適合結果を表すグラフである。 図９Ａ及び図９Ｂは、位相の決定結果を示すグラフである。 図８のグラフ上に、胸椎、腰椎及び仙骨の対応位置を表記した概略説明図である。 腰椎周辺のサジタル断面画像の一例を表す図である。 第１の時点で算出されたプロファイルと、第２の時点で算出されたプロファイルとを照合する具体的方法を示す概略説明図である。

本発明に係る椎骨セグメンテーション方法について、それを実施する椎骨セグメンテーション装置との関係において好適な実施の形態を掲げ、添付の図面を参照しながら以下、詳細に説明する。

図１は、本実施の形態に係る椎骨セグメンテーション装置１０の電気的な概略ブロック図である。椎骨セグメンテーション装置１０は、入力部１２、制御部１４（ＣＰＵ等の情報処理装置）、出力部１６、及びメモリ１８（記録媒体）を備える。メモリ１８には、本実施の形態の椎骨セグメンテーション装置１０として機能させるためのプログラムが格納されている。入力部１２は、外部からの信号を受ける。入力部１２は、画像取得部２０を備える。画像取得部２０は、椎体を含む三次元医用画像を取得する。画像取得部２０は、ＣＴ画像、ＭＲＩ画像、ＲＩ画像、ＰＥＴ画像、及びＸ線画像等の医用画像のみならず、後述する三次元の生成画像又は人工的に作成された三次元医用画像等、二次元又は三次元医用画像等を取得してよい。

制御部１４は、画像生成部２２、椎骨中心線検出部２４、第１特徴量算出部２６、第２特徴量算出部２８、第３特徴量算出部３０、関数適合部３２、位相決定部３４、椎骨位置推定部３６、識別情報付与部３７、及び特徴量判別部３８を備える。

画像生成部２２は、画像取得部２０が取得した椎体を含む三次元医用画像に基づいて、所定の軸方向（例えば、体軸方向）に沿った複数の第１断層画像（例えば、アキシャル断層画像）を生成する。画像生成部２２は、生成した複数の第１断層画像をメモリ１８に記録する。なお、画像取得部２０が、椎体を含む三次元医用画像のアキシャル断面である複数の第１断層画像を直接取得した場合は、画像生成部２２は、そのまま、複数の第１断層画像を取得して、メモリ１８に記録する。

また、画像生成部２２は、画像取得部２０が取得した三次元医用画像に基づいて、後述する各椎骨の中心軸に沿った複数の第２断層画像を生成する。さらに、画像生成部２２は、生成する画像空間領域の範囲を適宜変更できる。

椎骨中心線検出部２４は、画像生成部２２が生成した複数の第１断層画像のそれぞれに含まれる各椎骨の中心線を検出する。

図２は、各椎骨の配列を表すサジタル画像の模式図である。椎体群４０は、略円筒状の皮質を有する１５個の椎骨４２から構成されている。各椎骨４２は、体軸（Ｚ軸）方向に対しＳ字状に湾曲するように配置されている。すなわち、各椎骨４２は、そのＳ字状の中心線Ａに沿って配置されているともいえる（以下、中心線Ａに沿った座標軸を「Ｚ’軸」という。）。

また、隣接する椎骨４２の間には椎間板４４が介在する。説明の便宜上、図２において椎間板４４の線図による表示を省略し、椎骨４２間の隙間として表現する。

さらに、椎骨４２、椎間板４４の区別をするため、これらの参照符号４２、４４の後にアルファベット文字をそれぞれ付加する。つまり、上方から順番に、椎骨４２ａ〜４２ｏ、椎間板４４ａ〜４４ｏとする。なお、最下部の椎間板４４ｏの下面側と骨盤４６とは接続されている。

複数の第１断層画像を用いて各椎骨４２の中心線Ａを直接的に検出するには、きわめて高度な画像処理技術を要する。そこで、椎体の構造的特徴に着目した種々の検出手法を採り得る。例えば、画像処理による検出が比較的容易である図示しない脊髄の中心線を予め求めておき、該脊髄と椎骨４２との相対的位置関係に基づいて、各椎骨４２の中心線Ａを正確に検出することができる（詳細は、特許文献１参照）。

なお、この脊髄の中心線の検出方法として、テンプレートマッチング手法や、画定手法を用いてもよいし、統合学習機械を作る手法であるＡｄａｂｏｏｓｔに基づいたラーニング手法を用いてもよい。

第１特徴量算出部２６は、椎骨４２（図２参照）についてのＺ’軸方向に交叉する断面形状の鮮明度を表す第１特徴量を算出する。第１特徴量として、例えば、Ｚ’軸上の所定の点Ｐを中心とする円環状の模様を抽出する特徴量を用いる。この特徴量の算出には、後述するヘッセ行列の固有値解析法を用いる。第１特徴量算出部２６は、算出した第１特徴量をメモリ１８に記録する。

第２特徴量算出部２８は、椎骨４２（図２参照）についてのＺ’軸方向に平行する断面形状の鮮明度を表す第２特徴量を算出する。第２特徴量として、例えば、Ｚ’軸上の所定の点Ｐを中心として該Ｚ’軸方向に延在する管状の模様を抽出する特徴量を用いる。この特徴量の算出には、後述するヘッセ行列の固有値解析法を用いる。第２特徴量算出部２８は、算出した第２特徴量をメモリ１８に記録する。

第３特徴量算出部３０は、各椎骨４２の配列の規則性を表す第３特徴量を算出する。第３特徴量として、例えば、第１及び第２特徴量の加重和を用いる。第３特徴量算出部３０は、算出した第３特徴量をメモリ１８に記録する。

関数適合部３２は、Ｚ’軸上の各点毎に算出された第３特徴量を所定の適合関数（周期関数又は準周期関数）に適合させる。ここで、「周期関数」とは、Ｚ’軸座標の位置に対し周期性を有し、該位置によらず周期が一定である関数をいう。また、「準周期関数」とは、Ｚ’軸座標の位置に対し周期性を有するが、該位置によって周期が変調する関数をいう。関数適合部３２は、適合させた適合関数の関数形を特定するための各種定数をメモリ１８に記録する。

位相決定部３４は、関数適合部３２が適合させた適合関数の中に変数として含まれている位相の値を決定する。位相決定部３４は、例えば、動的計画法を用いてＺ’軸上の各点毎に位相の値を決定する。位相決定部３４は決定した位相をメモリ１８に記録する。

椎骨位置推定部３６は、第３特徴量に基づいてＸＹＺ座標上での各椎骨４２の位置を推定する。椎骨位置推定部３６は、関数適合部３２及び位相決定部３４によって決定した適合関数から各椎骨４２のＺ’軸上の位置を推定し、該位置をＸＹＺ座標上の位置に変換する。椎骨位置推定部３６は推定した各椎骨４２の位置をメモリ１８に記録する。

識別情報付与部３７は、ＸＹＺ座標上で（又はＺ’軸上で）位置が推定された各椎骨４２に対し、予め定義された識別情報を付与する。ここで、識別情報は、各断層画像間での相対的位置関係を表す各種情報であっても、被写体の構造における絶対的位置関係を表す各種情報であってもよい。なお、絶対的位置関係を表す各種情報には、解剖学上の椎骨４２の種別（詳細は後述する。）が含まれてもよい。

特徴量判別部３８は、第１の時点ｔ１及び第２の時点ｔ２に、第３特徴量算出部３０によってそれぞれ算出した第３特徴量の特性（Ｚ’軸方向に対するプロファイル）を照合し、これらの特性が合致するか否かを判別する。特徴量判別部３８は、種々のアルゴリズムを用いて、両者のプロファイルのずれの程度を検出するとともに、そのシフト量Δｚ’を算出する。

出力部１６は、椎骨位置推定部３６が推定した各椎骨４２の位置を外部に出力する。出力部１６は、複数箇所の基準位置（椎骨４２又は椎間板４４の端点や中心位置等）を予め定めておき、基準位置をそれぞれ出力してもよい。また、出力部１６は、各椎骨４２の位置情報とともに、識別情報付与部３７により紐付けられた識別情報を併せて出力してもよい。

次に、椎骨セグメンテーション装置１０の動作を図３のフローチャートにしたがって説明する。

先ず、画像生成部２２は、画像取得部２０が取得した椎体を含む三次元医用画像に基づいて、体軸方向（Ｚ軸）に沿ったＸ−Ｙ断面画像、すなわち、複数の第１断層画像を生成する（ステップＳ１）。

なお、画像取得部２０がアキシャル断面である複数の第１断層画像を取得した場合は、画像生成部２２は、画像取得部２０が取得した複数の第１断層画像を取得する。

次いで、椎骨中心線検出部２４は、ステップＳ１で取得した複数の第１断層画像から、それぞれに含まれる椎骨４２の中心線Ａを検出する（ステップＳ２）。そして、椎骨中心線検出部２４は、上述した検出方法を用いて図２に示す椎骨４２の中心線Ａを検出する。その後、中心線Ａは、新たな座標軸（Ｚ’軸）として設定される。

次いで、画像生成部２２は、ステップＳ１で取得した三次元医用画像に基づいて、ステップＳ２で検出した中心線Ａ（Ｚ’軸）に沿ったＸ’−Ｙ’断面画像、すなわち、複数の第２断層画像を生成する（ステップＳ３）。なお、画像生成部２２は、取得した三次元医用画像のうちすべての画像情報を用いることなく、椎体群４０の存在領域を網羅するＺ’軸周辺の一部領域において再構成可能な画像情報のみを用いて、複数の第２断層画像を生成すればよい。そうすれば、メモリ１８の使用量や、制御部１４による演算時間を低減することができる。

次いで、第１特徴量算出部２６は、ステップＳ３で取得した複数の第２断層画像から第１特徴量を算出する（ステップＳ４）。

第１特徴量は、次の（１）式で示される。

なお、（ｘ’，ｙ’，ｚ’）は、Ｘ’Ｙ’Ｚ’座標系での位置を表す。また、積分領域Ｒは第２断層画像上の点であって、椎骨４２の皮質領域を表す。

さらに、３×３ヘッセ行列における固有値及び固有ベクトルは、それぞれ（λ₁，λ₂，λ₃）、及び（Ｅ₁，Ｅ₂，Ｅ₃）とする。さらにまた、λ₁≦λ₂≦λ₃である。ベクトルｄは、Ｚ’軸上の点（０，０，ｚ’）から（ｘ’，ｙ’，ｚ’）に向かうベクトルであり、ｄ＝（ｘ’，ｙ’，０）である。なお、εは、ゼロ割防止のために設けられた、微小な正の整数である。

図４Ａ及び図４Ｂは、各断面位置における第１特徴量（被積分関数）の二次元分布を表す模式図である。つまり、第２断層画像上での第１特徴量（被積分関数）の数値を画像の濃淡として可視化している。第１特徴量が大きい位置での濃度は濃く、第１特徴量が小さい位置での濃度は薄くなるように表示されている。

図４Ａは、断面位置Ｔ１（図２参照）における第１特徴量の算出値を表す可視画像５０の模式図である。断面位置Ｔ１は椎骨４２の中間位置であるから、その第２断層画像には椎骨４２の皮質の断面形状が含まれている。よって、可視画像５０上には、円環状の黒い模様５２が鮮明に現れている。その結果、積分領域Ｒ内での被積分関数の総和である第１特徴量は大きくなる。

一方、図４Ｂは、断面位置Ｔ２（図２参照）における第１特徴量の算出値を表す可視画像５４の模式図である。断面位置Ｔ２は椎間板４４ｍの位置であるから、その第２断層画像には椎骨４２の皮質の断面形状が含まれず、代わりに椎間板４４の辺縁部が含まれている。よって、可視画像５４上には、円環状の薄い模様５６が僅かに現れている。その結果、積分領域Ｒ内での被積分関数の総和である第１特徴量は小さくなる。

第１特徴量の別の算出例として、例えば、エッジの強度を検出するラプラシアンフィルタ（２次微分フィルタ）を用いてもよい。また、局所領域の画素にガウス分布の重み付けして平滑化した後、ラプラシアンを作用してそのゼロクロスをエッジとして検出するラプラシアンオブガウシアンフィルタ（Ｌａｐｌａｃｉａｎ Ｏｆ Ｇａｕｓｓｉａｎ Ｆｉｌｔｅｒ；以下、「ＬＯＧフィルタ」という。）を用いてもよい。

次いで、第２特徴量算出部２８は、ステップＳ３で取得した複数の第２断層画像から第２特徴量を算出する（ステップＳ５）。

第２特徴量は、次の（２）式で示される。

なお、（ｘ’，ｙ’，ｚ’）は、Ｘ’Ｙ’Ｚ’座標系での位置を表す。また、積分領域Ｃは中心線Ａ（Ｚ’軸）の周辺領域を表す。

（１）式と同様に、３×３ヘッセ行列における固有値及び固有ベクトルは、それぞれ（λ₁，λ₂，λ₃）、及び（Ｅ₁，Ｅ₂，Ｅ₃）とする。また、λ₁≦λ₂≦λ₃である。さらに、ｅｚは、Ｚ’軸における単位ベクトル（０，０，１）である。なお、εは、ゼロ割防止のために設けられた微小な正の整数である。

図５は、椎骨４２ｍの三次元形状と、その内部点Ｐにおけるヘッセ行列の固有ベクトルＥ₁〜Ｅ₃との対応関係を示す模式図である。ヘッセ行列の最小固有値λ₁に対応する固有ベクトルＥ₁は、椎骨４２ｍの延在方向（Ｚ’軸方向）を指向する。したがって、（２）式の被積分関数の形状から諒解されるように、第２特徴量は、点Ｐを基準として濃度勾配がない平坦な形状、特にＺ’軸方向のベクトル成分のみを抽出する。すなわち、第２特徴量は、Ｚ’軸方向に延在する管状の模様を抽出する特徴量であるといえる。

次いで、第３特徴量算出部３０は、ステップＳ４で算出した第１特徴量と、ステップＳ５で算出した第２特徴量とに基づいて第３特徴量を算出する（ステップＳ６）。

第３特徴量は、次の（３）式で示される。

なお、αは任意の重み付け係数であり、Ｇ（ｚ’，σ）は標準偏差をσとするガウス関数である。

図６Ａ〜図６Ｃは、Ｚ’軸方向における第１〜第３特徴量のプロファイルである。

図６Ａは、Ｚ’軸方向における第１及び第２特徴量のプロファイルである。第１特徴量（実線）及び第２特徴量（破線）は、椎骨４２の中央位置（図２に示す断面位置Ｔ１）で極大値をとり、椎間板４４の位置（図２に示す断面位置Ｔ２）で極小値をとる。

図６Ｂは、第２特徴量にＬＯＧフィルタを作用させた後のプロファイルである。第２特徴量のプロファイル（図６Ａ参照）のうち２次微分が０になっている箇所が抽出されるので、図６Ｂに示すプロファイルは周期的な関数形状を有している。このプロファイルは、椎骨４２の中央位置（図２に示す断面位置Ｔ１）で極大値をとり、椎間板４４の位置（図２に示す断面位置Ｔ２）で極小値をとる。

ところで、ガウス関数Ｇ（ｚ’，σ）におけるσの値に応じて、図６Ｂに示すプロファイルの形状が変化する。本実施の形態では、所定の範囲σ₀〜σ₁において、そのプロファイルの値が最大となるσを選択する。

図６Ｃは、Ｚ’軸方向における第３特徴量のプロファイルである。第３特徴量は、第１特徴量ｆ₁（ｚ’）と、第２特徴量ｆ₂（ｚ’）にＬＯＧフィルタを作用した特徴量とを加算した値である。

これにより、各椎骨４２の配列の規則性を表す第３特徴量ｆ₃（ｚ’）が生成される。この第３特徴量は、椎骨４２の中央位置（図２に示す断面位置Ｔ１）で極大値をとり、椎間板４４の位置（図２に示す断面位置Ｔ２）で極小値をとる。

次いで、関数適合部３２は、ステップＳ６で算出した第３特徴量を所定の関数に適合させる（ステップＳ７）。この適合方法について、図７のフローチャートを参照しながら詳細に説明する。

先ず、関数適合部３２は、適合関数の関数形を設定する（ステップＳ７１）。関数形が複数存在する場合は、メモリ１８からその情報を取得する。

図２に示すように、椎体群４０は、椎骨４２と椎間板４４が交互に配置されているという構造的特徴を有している。そこで、一例として三角関数のような周期関数を適合関数として用いることができる。

さらに、椎体群４０は、椎骨４２の軸方向の高さが頚椎（椎骨４２ａ）から腰椎（椎骨４２ｏ）にかけて次第に大きくなる、という構造的特徴も有している。そこで、次の（４）式に示す準周期関数ｇ（ｚ’）を適合関数として用いてもよい。

なお、ｃ₁〜ｃ₃は、ｇ（ｚ’）の形状を決定するための任意の定数である。また、ｃ₂＝０のとき、ｇ（ｚ’）は周期関数である。

このｇ（ｚ’）の値の定義は、上述した第３特徴量ｆ₃（ｚ’）の値の定義と一致する。すなわち、椎骨４２の中央位置（図２に示す断面位置Ｔ１）で極大値をとり、椎間板４４の位置（図２に示す断面位置Ｔ２）で極小値をとる。

次いで、関数適合部３２は、第３特徴量をグローバルに適合させる（ステップＳ７２）。ここで、「グローバルに適合させる」とは、Ｚ’軸上の位置ｚ’の採り得る範囲全体で適合させることをいう。

関数適合部３２は、最小２乗法等の多変数解析を行うことにより、最適な定数ｃ₁〜ｃ₃を決定することができる。例えば、評価値Ｖ₁は、次の（５）式で表される。

この評価値Ｖ₁が最大となるような定数ｃ₁〜ｃ₃を選択する。この場合、定数ｃ₁〜ｃ₃の採り得る範囲を予め定めておき、その範囲内にあるすべてのｃ₁〜ｃ₃の組み合わせで探索する。

図８は、第３特徴量の適合結果を表すグラフである。横軸はＺ’軸上の位置ｚ’であり、縦軸は第３特徴量である。このように、第３特徴量算出部３０が第２断層画像から算出した第３特徴量を、（４）式に示す適合関数に適合させることができる。

なお、第３特徴量ｆ₃（ｚ’）の平均値が０になるように（３）式に示すαの値を調整することが好ましい。かかる場合、（５）式で算出される評価値Ｖ₁は、相関係数の値に対応するからである。

また、Ｚ’軸上の位置ｚ’に応じて第３特徴量及び適合関数の振幅を略一致させることが好ましい。かかる場合、適合精度がさらに向上するからである。

次いで、位相決定部３４は、ステップＳ７２で得られた３つの定数ｃ₁〜ｃ₃を用いて、ローカルに適合させる（ステップＳ７３）。ここで、「ローカルに適合させる」とは、「グローバルに適合」させた結果を踏まえ、Ｚ’軸上の位置ｚ’毎にさらに適合させることをいう。具体的には、（４）式の適合関数に対し、位置ｚ’毎に位相δ（ｚ’）の項を導入して最適化を行うようにする。

位相決定部３４は、最小２乗法等の多変数解析を行うことにより、最適な定数δ（ｚ’）を決定することができる。例えば、評価値Ｖ₂は、次の（６）式で表される。

この評価値Ｖ₂が最大となるような位相δ（ｚ’）を決定する。ここでは、動的計画法を用いて、δ₁、δ₂、‥‥と位相を離散値として順次決定していき、最終的には全範囲（δ₁〜δ_N）での位相δ（ｚ’）を決定する。

図９Ａ及び図９Ｂは、位相の決定結果を示すグラフである。横軸はＺ’軸上の位置ｚ’であり、縦軸は位相δ（ｚ’）である。

図９Ａは、典型的な椎体群４０における位相δ（ｚ’）の決定結果を示す。椎骨４２はそれぞれ規則的に配置されているので、グローバルな適合（ステップＳ７２）のみであっても十分な適合精度を有する。したがって、位相δ（ｚ’）は位置ｚ’によらず一定である。

図９Ｂは、ある椎骨４２が一箇所だけ破損した椎体群４０における位相δ（ｚ’）の決定結果を示す。この場合は、破損位置を境界点として、左側（上方側）と右側（下方側）の椎骨４２の規則性が分断される。したがって、図９Ｂに示すように、破損位置を境界点として位相δ（ｚ’）が不連続となる。ただし、破損位置から左側（又は右側）領域に関しては椎骨４２の規則性は失われないので、位相δ（ｚ’）は位置ｚ’によらず一定である。

このように、局所的な椎骨４２の位置ずれをも考慮した関数の適合が可能であるので、各椎骨４２の位置の推定精度がさらに向上する。

次いで、椎骨位置推定部３６は、Ｚ’軸上における各椎骨４２の位置を推定する（ステップＳ８）。

椎骨位置推定部３６は、適合させた適合関数に基づいて、各椎骨４２の位置を推定する。例えば、式（４）に示す適合関数の例において、ｎ番目の椎間板４４の位置は次の（７）式のように求めることができる。

上記の（７）式はグローバルに適合した結果のみを考慮した位置である。ローカルに適合した結果をさらに考慮すると、δ（ｚ’）を一次微小量として、次の（８）式のように近似的に求めることができる。

次いで、椎骨位置推定部３６は、Ｚ軸上における各椎骨４２の位置を変換する（ステップＳ９）。椎骨位置推定部３６は、ステップＳ８で推定した各椎骨４２のＺ’軸上の位置をＸＹＺ座標の位置に変換する。なお、Ｚ’軸上の点ＰとＸＹＺ座標との位置関係は既知であるので、このような座標変換は容易である。

最後に、識別情報付与部３７は、分離・認識された各椎骨４２に対して識別情報を付与する（ステップＳ１０；第１付与ステップ）。

本実施の形態では、解剖学上の椎骨４２の種別を識別情報として用いる。より詳細には、脊椎（椎体群４０）を構成する各椎骨４２を、第１〜第７頚椎と、第１〜第１２胸椎と、第１〜第５腰椎と、第１〜第５仙椎（仙骨）と、数個の尾椎（尾骨）とにそれぞれ分類して識別する。

図１０は、図８のグラフ上に、胸椎、腰椎及び仙椎の対応位置を表記した概略説明図である。横軸はＺ’軸上の位置ｚ’であり、縦軸は第３特徴量（図８参照）である。

例えば、識別情報付与部３７は、外部からの入力指示に応じて、所定のピークＰｋに応じた椎骨４２に対して「第３腰椎」の識別情報を付与する。この際、識別情報付与部３７は、各椎骨４２の配列順番と種別とを対応付けたテーブルデータをメモリ１８から取得する。なお、前記入力指示は、医師又は技師（以下、ユーザという。）による図示しない入力装置を介しての入力指示であっても、椎骨セグメンテーション装置１０上で実現される画像認識処理の結果に基づいた自動的な入力指示であってもよい。

そして、識別情報付与部３７は、所定のピークＰｋを基準とした前後２個以内である５個のピークに応じた椎骨４２を、５個の腰椎として特定する。すなわち、識別情報付与部３７は、前記テーブルデータを参照することで、特定された各ピークに応じた椎骨４２に対し、位置ｚ’の小さい側から順に、第１〜第５腰椎（既に付与された第３腰椎を除く。）の識別情報をそれぞれ付与する。

また、識別情報付与部３７は、第１腰椎に応じたピークよりも左方（Ｚ’軸の負方向）に位置する１２個のピークに応じた椎骨４２を、１２個の胸椎として特定する。すなわち、識別情報付与部３７は、前記テーブルデータを参照することで、特定された各ピークに応じた椎骨４２に対し、位置ｚ’の小さい側から順に、第１〜第１２胸椎の識別情報をそれぞれ付与する。

一方、識別情報付与部３７は、第５腰椎に応じたピークよりも右方（Ｚ’軸の正方向）に位置する５個のピークに応じた椎骨４２を、５個の仙椎として特定する。すなわち、識別情報付与部３７は、前記テーブルデータを参照することで、特定された各ピークに応じた椎骨４２に対し、位置ｚ’の小さい側から順に、第１〜第５仙椎の識別情報をそれぞれ付与する。

その後、識別情報付与部３７は、図示しない頚椎又は尾椎に応じた各ピークを特定するとともに、特定された各ピークに応じた椎骨４２に対し、前記テーブルデータの記述に従って識別情報をそれぞれ付与する。

このようにして、識別情報付与部３７は、各椎骨４２に対して識別情報を付与する（ステップＳ１０）。なお、識別情報付与部３７は、各椎骨４２間に介在する椎間板４４に対しても、椎骨４２の場合と同様の方法を用いて所定の識別番号を付与してもよい。

また、第７頚椎と、その別名である「隆椎」とを併せて対応付けてもよく、どの用語を使用するかをユーザが選択できるようにしてもよい。

さらに、第１頚椎を基準として、最上部から最下部への順に、固有の椎骨番号を付与してもよい。あるいは、識別情報を付与する起点は、最上部の椎体（第１頚椎）又は最下部の椎体（第５仙骨）のいずれでもよい。

その後、図示しない表示装置に医用画像を表示する際、該医用画像と併せて、出力部１６側に出力した識別情報を可視化して表示してもよい。ここで、識別情報を可視化するために、公知の画像表示技術が種々用いられる。

図１１は、腰椎周辺のサジタル断面画像の一例を表す図である。このように、識別情報（第１〜第５腰椎）がそれぞれ付与された各椎骨４２上の位置に、第１〜第５腰椎なる文字列を断面画像１００に重畳して表示させてもよい。これにより、ユーザは、椎骨４２の解剖学上の種別（名称）を一見して把握可能であり、診断上便宜である。

以上のように、各椎骨４２の中心線Ａ（Ｚ’軸）に沿った断層画像を生成し、各椎骨４２についてのＺ’軸方向に交叉及び平行する断面形状の鮮明度を表す第１及び第２特徴量を、生成された前記断層画像に基づいてＺ’軸上の点毎に算出し、各椎骨４２の配列の規則性を表す第３特徴量を、算出された前記第１及び第２特徴量に基づいてＺ’軸上の点毎に算出し、Ｚ’軸上の点毎に算出された前記第３特徴量に基づいてＺ’軸上における各椎骨４２の位置を推定するようにしたので、異なる２方向からの断面形状に基づいた定量的な濃度解析が可能であり、脊椎の形状が全体的又は局所的に変形している被写体の三次元医用画像に対しても、各椎骨を適切に分離・認識できる。脊椎湾曲症等により全体的に脊椎の形状が変形している場合、特に効果がある。

また、所定の周期関数又は準周期関数には位相が変数として含まれ、動的計画法を用いて前記中心軸上の点毎に前記位相を順次決定することが好ましい。これにより、関数への適合誤差を低減可能であり、各椎骨４２の位置をさらに高精度に推定できる。骨折等により局所的に脊椎の形状が変形している場合、特に効果がある。

ところで、本実施の形態に係る椎骨セグメンテーション方法を、臨床に応用することができる。同一の被写体において過去に取得した椎骨４２の分離結果と、現在取得した椎骨４２の分離結果とを比較することで、前記被写体の脊椎の形状変化を経時で追跡することができる。すなわち、本方法は、脊椎湾曲症を含む脊椎疾患の診断（予防、検診、診察、経過観察等）を支援するツールとして有用である。

以下、異なる時点での椎骨４２の分離結果を比較する際における、椎骨セグメンテーション装置１０の動作について、図１、図３及び図１２を参照しながら詳細に説明する。

先ず、第１の時点ｔ１（例えば、過去）において、図３に示すフローチャートに従って、椎骨４２の分離・認識（ステップＳ１〜Ｓ９）、及び識別情報の付与（ステップＳ１０）に関する処理が実行される。その際、メモリ１８は、各椎骨４２の位置情報と、該位置情報に紐付けられた識別情報とをそれぞれ記録する。あるいは、図示しない外部装置（データベースサーバ等）にこれらを記録させてもよい。

次いで、第２の時点ｔ２（例えば、現在）において、図３に示すフローチャートに従って、椎骨４２の分離・認識（ステップＳ１〜Ｓ９）に関する処理が実行される。なお、第２の時点ｔ２での各処理においては、第１の時点ｔ１での各処理と同一の処理条件を適用することが好ましい。特に、三次元医用画像の取得前に画素値分布の正規化を予め行っておくことが好ましい。なぜならば、制御部１４による椎骨４２の分離・認識結果の再現性を高めるためである。

次いで、特徴量判別部３８は、第１の時点ｔ１で算出された第３特徴量と、第２の時点ｔ２で算出された第３特徴量とを照合する。第１の時点ｔ１と第２の時点ｔ２との間で被写体の脊椎（椎体群４０）の形状は殆ど変化しなかった場合の例を以下に示す。

図１２の上側のグラフは、第１の時点ｔ１での取得画像から算出された第３特徴量のプロファイル｛以下、プロファイルｆ_A（ｚ’）という。｝を表す。図１２の下側のグラフは、第２の時点ｔ２での取得画像から算出された第３特徴量のプロファイル｛以下、プロファイルｆ_B（ｚ’）という。｝を表す。図１２から諒解されるように、プロファイルｆ_A（ｚ’）及びプロファイルｆ_B（ｚ’）は略同じ傾向を表している。ただし、撮像開始位置のずれが原因で、プロファイルｆ_B（ｚ’）は、プロファイルｆ_A（ｚ’）と比べて、Ｚ’軸の正方向にΔｚ’だけ全体的にシフトしたとする。

特徴量判別部３８は、プロファイルｆ_A（ｚ’）及びｆ_B（ｚ’）のマッチング処理を行う。特徴量判別部３８は、例えば、プロファイルｆ_A（ｚ’）及びｆ_B（ｚ’−Δｚ’）のＳＳＤ（Ｓｕｍ ｏｆ Ｓｑｕａｒｅ Ｄｉｆｆｅｒｅｎｃｅ）が最小となるように、シフト量Δｚ’を算出する。なお、マッチング処理の際の評価関数として、ＳＳＤの他、ＮＣＣ（Ｎｏｒｍａｌｉｚｅｄ Ｃｒｏｓｓ Ｃｏｒｒｅｌａｔｉｏｎ）やＳＡＤ（Ｓｕｍ ｏｆ Ａｂｓｏｌｕｔｅ Ｄｉｆｆｅｒｅｎｃｅ）等公知の手法を用いてもよい。また、第３特徴量のプロファイルが有する各ピーク値自体を比較してもよい。

そして、特徴量判別部３８は、プロファイルｆ _Ｂ （ｚ’）の特性がプロファイルｆ_Ａ（ｚ’）の特性と合致するか否かを判別する。特徴量判別部３８は、プロファイルｆ_Ａ（ｚ’）及びｆ_Ｂ（ｚ’−Δｚ’）が十分近似しているか否か、例えばこれらの相関係数が所定の閾値を超えるか否かで判別してもよい。合致すると判別された場合、特徴量判別部３８は、シフト量Δｚ’をメモリ１８に供給する。

次いで、識別情報付与部３７は、第１の時点ｔ１で付与された識別情報と同一の識別情報を付与する（ステップＳ１０；第２付与ステップ）。識別情報付与部３７は、特徴量判別部３８による照合結果に基づいて、第１の時点ｔ１での各ピークに相当する第２の時点ｔ２での各ピークを特定する。そして、識別情報付与部３７は、特定された各ピークに応じた椎骨４２に対し、第１の時点ｔ１での識別情報と同一の識別情報を付与する。識別情報の再付与に際し、メモリ１８から第１の時点ｔ１での椎骨４２の位置情報、第２の時点ｔ２での椎骨４２の位置情報、シフト量Δｚ’及び上述したテーブルデータが参照される。

なお、図１２では、第１の時点ｔ１が過去であり、第２の時点ｔ２が現在である場合（ｔ１＜ｔ２）を例示したが、その逆（ｔ１＞ｔ２）であってもよい。すなわち、第１の時点ｔ１が現在であり、第２の時点ｔ２が過去であっても本発明を適用できる。

以上のように、推定されたＺ’軸上の位置に応じた各椎骨４２に対し、予め定義された識別情報（例えば、解剖学上の椎骨４２の種別）をそれぞれ付与することで、医療の専門家としてのユーザにとって椎骨４２の識別が容易になる。

また、第１の時点ｔ１で付与された識別情報を、第２の時点ｔ２で推定されたＺ’軸上の位置に応じた各椎骨４２に対してそれぞれ付与することで、被写体の脊椎の形状変化を経時で容易に追跡できる。

なお、この発明は、上述した実施形態に限定されるものではなく、この発明の主旨を逸脱しない範囲で自由に変更できることは勿論である。被写体は人間に限られず、脊椎動物であれば本発明を適用し得ることは言うまでもない。

１０…椎骨セグメンテーション装置 １４…制御部
２２…画像生成部 ２６…第１特徴量算出部
２８…第２特徴量算出部 ３０…第３特徴量算出部
３２…関数適合部 ３４…位相決定部
３６…椎骨位置推定部 ３７…椎骨情報付与部
３８…特徴量判別部 ４０…椎体群
４２、４２ａ〜４２ｏ…椎骨 ４４、４４ａ〜４４ｏ…椎間板
５０、５４…可視画像 ５２、５６…模様

Claims (11)

1. 各椎骨の中心軸に沿った断層画像を生成する画像生成部と、
   前記各椎骨についての前記中心軸方向に交叉する断面形状の鮮明度を表す第１特徴量を、前記画像生成部により生成された前記断層画像に基づいて該中心軸上の点毎に算出する第１特徴量算出部と、
   前記各椎骨についての前記中心軸方向に平行する断面形状の鮮明度を表す第２特徴量を、前記画像生成部により生成された前記断層画像に基づいて該中心軸上の点毎に算出する第２特徴量算出部と、
   前記各椎骨の配列の規則性を表す第３特徴量を、前記第１及び第２特徴量算出部により算出された前記第１及び第２特徴量に基づいて前記中心軸上の点毎に算出する第３特徴量算出部と、
   前記第３特徴量算出部により前記中心軸上の点毎に算出された前記第３特徴量に基づいて、該中心軸上における前記各椎骨の位置を推定する椎骨位置推定部と
   を有することを特徴とする椎骨セグメンテーション装置。
2. 請求項１記載の椎骨セグメンテーション装置において、
   前記第１特徴量は、前記中心軸上の各点を中心とする円環状の模様を抽出する特徴量であることを特徴とする椎骨セグメンテーション装置。
3. 請求項１又は２に記載の椎骨セグメンテーション装置において、
   前記第２特徴量は、前記中心軸方向に延在する管状の模様を抽出する特徴量であることを特徴とする椎骨セグメンテーション装置。
4. 請求項１〜３のいずれか１項に記載の椎骨セグメンテーション装置において、
   前記中心軸上の各点毎に算出された前記第３特徴量を所定の周期関数又は準周期関数に適合させる関数適合部をさらに有し、
   前記椎骨位置推定部は、前記関数適合部により適合させた前記所定の周期関数又は準周期関数に基づいて前記中心軸上における前記各椎骨の位置を推定する
   ことを特徴とする椎骨セグメンテーション装置。
5. 請求項４記載の椎骨セグメンテーション装置において、
   前記所定の周期関数又は準周期関数には位相が変数として含まれ、
   動的計画法を用いて前記中心軸上の各点毎に前記位相を順次決定する位相決定部をさらに有する
   ことを特徴とする椎骨セグメンテーション装置。
6. 請求項１〜５のいずれか１項に記載の椎骨セグメンテーション装置において、
   前記椎骨位置推定部により推定された前記中心軸上の位置に応じた前記各椎骨に対し、予め定義された識別情報をそれぞれ付与する識別情報付与部をさらに有することを特徴とする椎骨セグメンテーション装置。
7. 請求項６記載の椎骨セグメンテーション装置において、
   前記識別情報付与部は、第１の時点で生成された前記断層画像に基づいて付与された前記識別情報と同一の識別情報を、第２の時点で生成された前記断層画像に基づいて推定された前記中心軸上の位置に応じた前記各椎骨に対してそれぞれ付与することを特徴とする椎骨セグメンテーション装置。
8. 請求項７記載の椎骨セグメンテーション装置において、
   前記第１の時点で生成された前記断層画像に基づいて算出された前記第３特徴量の特性が、前記第２の時点で生成された前記断層画像に基づいて算出された前記第３特徴量の特性と合致するか否かを判別する特徴量判別部をさらに有し、
   前記特徴量判別部により前記第３特徴量の特性が合致すると判別された場合、識別情報付与部は前記同一の識別情報を付与する
   ことを特徴とする椎骨セグメンテーション装置。
9. 請求項６〜８のいずれか１項に記載の椎骨セグメンテーション装置において、
   前記識別情報は、解剖学上の椎骨の種別であることを特徴とする椎骨セグメンテーション装置。
10. 各椎骨の中心軸に沿った断層画像を生成する生成ステップと、
    前記各椎骨についての前記中心軸方向に交叉する断面形状の鮮明度を表す第１特徴量を、生成された前記断層画像に基づいて該中心軸上の点毎に算出する第１特徴量算出ステップと、
    前記各椎骨についての前記中心軸方向に平行する断面形状の鮮明度を表す第２特徴量を、生成された前記断層画像に基づいて該中心軸上の点毎に算出する第２特徴量算出ステップと、
    前記各椎骨の配列の規則性を表す第３特徴量を、算出された前記第１及び第２特徴量に基づいて前記中心軸上の点毎に算出する第３特徴量算出ステップと、
    前記中心軸上の点毎に算出された前記第３特徴量に基づいて前記中心軸上における前記各椎骨の位置を推定する椎骨位置推定ステップと
    を備えることを特徴とする椎骨セグメンテーション方法。
11. コンピュータを、
    各椎骨の中心軸に沿った断層画像を生成する画像生成部、
    前記各椎骨についての前記中心軸方向に交叉する断面形状の鮮明度を表す第１特徴量を、生成された前記断層画像に基づいて該中心軸上の点毎に算出する第１特徴量算出部、
    前記各椎骨についての前記中心軸方向に平行する断面形状の鮮明度を表す第２特徴量を、生成された前記断層画像に基づいて該中心軸上の点毎に算出する第２特徴量算出部、
    前記各椎骨の配列の規則性を表す第３特徴量を、算出された前記第１及び第２特徴量に基づいて前記中心軸上の点毎に算出する第３特徴量算出部、
    前記中心軸上の点毎に算出された前記第３特徴量に基づいて前記中心軸上における前記各椎骨の位置を推定する椎骨位置推定部
    として機能させることを特徴とするプログラム。

Images