実施の形態の血管解析装置および血管解析方法は、医用画像診断装置で生成された医用画像に含まれる血管領域を構造流体解析するためのコンピュータ装置に適用できる。このコンピュータ装置は、医用画像診断装置に組み込まれていても良いし、医用画像診断装置とは別体のワークステーション等であっても良い。以下、実施の形態の血管解析装置および血管解析方法が適用されたコンピュータ装置が組み込まれた医用画像診断装置を、図面を参照して詳細に説明する。
(概要)
冠動脈等の動脈の流体解析では、血管またはプラークの硬さの影響、心臓拍動が動脈流へ及ぼす影響、動脈入口および動脈出口の境界条件の影響、を適切に考慮する必要がある。
このため、実施の形態の医用画像診断装置は、以下の構造流体解析を行う。すなわち、医用画像診断装置は、医用3次元画像に時間軸を加えた医用4次元画像を画像解析して得られる3次元変形履歴情報、および画像追尾から得られる変位履歴情報から、動脈の血管壁の形状変形指標と変位拘束条件を算出する。また、医用画像診断装置は、硬さ算定領域では、血管壁表面のみ、血管壁の力学モデルの節点に、強制変位履歴による変位自由度拘束条件を与える。
また、医用画像診断装置は、硬さ算定領域以外の領域では、血管壁の力学モデルにおける表面と内部の節点に強制変位履歴を与え、変位拘束条件を与えた構造流体解析を行う。これにより、硬さ算定領域の材料モデルおよび入口、出口の境界条件を同定する構造流体解析を行うことができ、高精度な血管および血行動態の解析が可能となる。そして、医用画像診断装置は、このような高精度な構造流体解析で得られた正常時および狭窄時のノミナルモデルを作成し、診断対象データとのモデル同定から狭窄度を診断する。
具体的には、実施の形態の医用画像診断装置は、CT(Computed Tomography)診断装置から血管の4次元画像を取得し、血管の画像追尾から得られた血管部の各点の移動ベクトルを強制変位履歴として算出する。また、医用画像診断装置は、画像追尾情報から血管内腔の容積変化、直径変化、形状変化、または位相変化等の物理量を算出する。また、医用画像診断装置は、算出した物理量を指標として、指定領域部の材料モデル、または境界条件を同定する。
医用画像診断装置は、このような画像追尾と構造流体解析と逆解析(統計的同定法)に基づく血管解析により、血行動態における血圧、血流、応力、またはひずみ等の狭窄等の血管病変の予防、診断に必要な指標を高精度に解析する。そして、医用画像診断装置は、狭窄の有無が判明している典型例を対象に解析を行い、低次の標準モデル(ノミナルモデル)を予め作成しておき、診断対象データとのモデル同定から狭窄度を診断する。
このような実施の形態の医用画像診断装置は、被検体をスキャンするための撮像機構を装備する如何なる種類の画像診断装置にも適用可能である。実施の形態の医用画像診断装置としては、例えばX線コンピュータ断層撮影装置(X線CT装置)、磁気共鳴診断装置、超音波診断装置、SPECT(Single Photon Emission CT)装置、PET(Positron Emission Tomography)装置、および放射線治療装置等に適宜利用可能である。以下、説明を具体的に行うため実施の形態の医用画像診断装置は、X線コンピュータ断層撮影装置であるものとする。
図1は、実施の形態の医用画像診断装置(X線コンピュータ断層撮影装置)の概略的なハードウェア構成図である。図1に示すように、X線コンピュータ断層撮影装置は、CT架台10とコンソール20とを有する。CT架台10は、コンソール20の架台制御部23の制御に従って、X線で被検体の撮像部位をスキャンする。撮像部位は、例えば心臓である。
CT架台10は、X線管11、X線検出器13、およびデータ収集装置15を有している。X線管11とX線検出器13とは、回転軸Z回りに回転可能にCT架台10に装備されている。X線管11は、造影剤が注入された被検体にX線を照射する。X線検出器13は、X線管11から発生され被検体を透過したX線を検出し、検出されたX線の強度に応じた電気信号を発生する。
データ収集装置15は、X線検出器13から電気信号を読み出してデジタルデータに変換する。1ビュー毎のデジタルデータのセットは、生データセットと呼ばれている。複数のスキャン時刻に関する時系列の生データセットは、非接触データ伝送装置(図示しない)によりコンソール20に伝送される。
コンソール20は、システム制御部21を中枢として、架台制御部23、再構成装置25、画像処理装置27、入力機器29、表示機器31、および記憶装置33を有している。
架台制御部23は、ユーザにより入力機器29を介して設定されたスキャン条件に応じてコンソール20内の各装置を制御する。
再構成装置25は、生データセットに基づいて、被検体に関するCT画像のデータを発生する。具体的には、まず、再構成装置25は、生データセットに前処理を施して投影データセットを発生する。前処理としては、対数変換や不均一補正、キャリブレーション補正等が含まれる。次に、再構成装置25は、投影データセットに画像再構成処理を施してCT画像のデータを発生する。画像再構成アルゴリズムとしては、FBP(filtered backprojection)法等の解析学的画像再構成法や、ML−EM(maximum likelihood expectation maximization)法やOS−EM(ordered subset expectation maximization)法等の逐次近似画像再構成等の既存のアルゴリズムが採用可能である。
実施の形態において再構成装置25は、時系列の投影データセットに基づいて、時系列のCT画像のデータを発生する。CT画像は、造影剤により造影された血管に関する画素領域(以下、血管領域と呼ぶことにする。)を含んでいる。なお、CT画像は、CT値の2次元空間分布を表現するスライスデータであっても良いし、CT値の3次元空間分布を表現するボリュームデータであっても良い。以下、CT画像はボリュームデータであるとする。時系列のCT画像のデータは、記憶装置33に記憶される。
画像処理装置27は、時系列のCT画像に基づいて、力学モデルを構築して構造流体解析を実行する。画像処理装置27の処理の詳細については後述する。
入力機器29は、ユーザからの各種指令や情報入力を受け付ける。入力機器29としては、キーボードやマウス、スイッチ等が利用可能である。
表示機器31は、CT画像や構造流体解析結果等の種々の情報を表示する。表示機器31としては、例えばCRTディスプレイや、液晶ディスプレイ、有機ELディスプレイ、プラズマディスプレイ等が適宜利用可能である。
記憶装置33は、ハードディスク装置等の種々の記憶媒体により構成される。記憶装置33は、時系列の投影データや時系列のCT画像データ等の種々のデータを記憶する。例えば、記憶装置33は、時系列のCT画像データをDICOM(digital imaging and communications in medicine)規格に準拠した医用画像ファイル形式で記憶する。また、記憶装置33は、外部機器により収集された医用データを時系列のCT画像データに医用画像ファイル内において関連付けて記憶しても良い。
システム制御部21は、中央演算処理装置(CPU:central processing unit)や読み出し専用メモリ(ROM:read only memory)、ランダムアクセスメモリ(RAM:random access memory)を有する。システム制御部21は、X線コンピュータ断層撮影装置の中枢として機能する。システム制御部21は、ROMやRAMに記憶されている血管解析プログラムを実行して実施の形態の血管構造解析処理を実行する。
なお、システム制御部21、画像処理装置27、入力機器29、表示機器31、および記憶装置33は、血管解析装置50を構成する。実施の形態のように血管解析装置50は、医用画像診断装置(X線コンピュータ断層撮影装置)に組み込まれていても良いし、医用画像診断装置とは別体のコンピュータ装置であっても良い。血管解析装置50が医用画像診断装置とは別体の場合、血管解析装置50は、医用画像診断装置やPACS(picturearchiving and communication systems)からネットワークを介して時系列のCT画像等の医用データを収集すれば良い。
次に、実施の形態の動作例を詳細に説明する。なお、実施の形態の医用画像診断装置は、心臓血管、頸動脈、または脳動脈等の人体のあらゆる部位の血管を解析対象とすることができる。以下、一例として、心臓の血管が解析対象であることとして説明を進める。
心臓の血管としては、例えば冠動脈と大動脈とが挙げられる。冠動脈は、大動脈の冠動脈起始部から始まり心筋表面を走行し、心外膜側から内膜側に入り込む。冠動脈は、心筋の内膜において無数の毛細管に分岐する。分岐後、無数の毛細管は、再び統合して大心静脈を形成し、冠静脈洞に接続する。冠血管系は、他の臓器と異なり、心筋の収縮および弛緩という力学的変化のなかで、灌流が保障されなければならないという点で特徴的である。
冠血流の特徴としては、心筋収縮による機械的血流阻害作用で冠動脈起始部の内圧が高くなる収縮期よりも、左心室拡張期に灌流圧が低下したときに多く流れることである。そのため、正常の冠動脈血流速波形は収縮期と拡張期の二峰性であり、拡張期血流が優位である。肥大型心筋症や大動脈弁狭窄症では、収縮期に逆行性波が認められ、大動脈逆流症では、収縮期順行波が大きくなる等、疾患によって特異な血流波形となることが知られている。また、拡張期の順行性波形は左室拡張機能、特に左室弛緩と密接な関係がある。左室弛緩遅延例では拡張期波形のピークが後ろにずれ、また減速脚がゆるやかになる傾向がある。また、このような症例では、頻拍時には拡張期の冠血流は十分に増大できず、心筋虚血を助長すると考えられている。
解剖学的に大動脈起始部から分岐する左右冠動脈に、大動脈圧に等しい冠灌流圧(すなわち、冠動脈が分枝する大動脈起始部の圧力)がかかることにより、冠血流が生じる。冠血流を決定するのは、大動脈圧である駆動圧と共に冠血管抵抗が重要である。140〜180μm以上の太い冠血管には、冠血管低抗の20%程度が存在するのに対し、100〜150μm以下の微小血管には、抵抗成分の残りの多くが存在するといわれる。従って、いわゆる冠狭窄等の無い場合には、抵抗値は冠微小血管の緊張性(トーヌス)に左右される。
血管抵抗因子としては、血管特性、動脈硬化、管狭窄、血液粘性、機械的因子があげられる。冠微小血管のトーヌスは、血管特性、心筋代謝(心筋酸素消費)、神経体液性因子、機械的因子、体液因子としての各種の血管作動性物質、血液粘性に規定され、さらに、心肥大、冠動脈硬化等を含めた様々な病変によっても影響され冠循環障害を起こす。
冠動脈血流拍動は、冠動脈血流の拍動パターン、心筋収縮による心筋内血流の制御、機械的刺激に対する心筋内血管の反応に影響される。心筋収縮が血流を阻害する機序としては、心筋内圧の上昇、心筋内血管容量の変化、心筋内血管の圧迫が挙げられる。心筋拡張期の血流規定因子には、拡張期の冠動脈圧、拡張期の血管外力、心拍数、心周期に占める拡張期の割合、心筋弛緩が存在する。
血管解析装置50は、時系列のCT画像に基づいて、力学モデルを構築し、力学モデルを利用して心臓の血管についての構造流体解析を実行し、血管内の力学的指標や血管流量指標を高精度に算出する。精度の高い力学的指標や血管流量指標を算出するためには、力学モデルに精度の高い潜在変数を割り当てる必要がある。血管解析装置50は、力学モデルを構築する際、初期的な力学モデルに逆解析を施して潜在変数を統計的に同定する。これにより血管解析装置50は、高精度の潜在変数を決定することができる。
力学的指標は、血管壁や血液に関する力学的な指標を意味する。血管壁に関する力学的指標としては、例えば血管壁の変位に関する指標、血管壁に生じる応力やひずみに関する指標、血管内腔に負荷される内圧分布に関する指標、血管の硬さ等をあらわす材料特性に関する指標等に分類される。血管の硬さ等をあらわす材料特性に関する指標は、血管組織の応力とひずみの関係を示す曲線の平均的な傾き等が挙げられる。
血液に関する力学的な指標における血管流量指標は、血管を流れる血液に関する血行動態の指標を意味する。血管流量指標としては、例えば血液の流量、血液の流速、血液の粘性等が挙げられる。なお、潜在変数は、例えば、血管の材料構成式、または血液の材料構成式といった材料モデルのパラメータ(例えばヤング率やポアソン比等)、血管内腔に負荷される内圧分布等の負荷条件パラメータ、構造解析や流体解析の境界条件パラメータ、時系列の形態指標や形状変形指標の不確定性に関連するばらつき分布パラメータの少なくとも一つを含む。
ここで、時系列の形態指標や形状変形指標の不確定性に関連するばらつき分布パラメータとは、医用画像データには、各CT値のノイズに起因したばらつき分布や、生体組織の境界閾値の曖昧性に起因する確率分布等が存在することを考慮して、血管組織や血液の境界座標および特徴点(血管分岐部や造影剤分散配置等)の空間座標における不確定性、または幾何学的構造パラメータ(芯線に垂直な断面の内腔半径等)の不確定性、または医用画像データ自体(CT値や境界閾値等)の不確定性を、確率分布として表現したものである。
力学モデルは、血管や血液の挙動を表現するための数値モデルである。力学モデルは、構造流体解析の手法に応じて異なるタイプを有している。例えば、力学モデルは、連続体力学モデルと簡易的力学モデルとに分類される。連続体力学モデルは、例えば、有限要素法(FEM:finite element method)や境界要素法に用いられる。簡易的力学モデルは、例えば材料力学に基づく材料力学モデルと流れ学に基づく流体力学モデルとに分類される。
なお、以下の説明において特に言及しない場合、力学モデルのタイプについては特に限定しないものとする。初期的な力学モデルは、潜在変数の確率分布や変数範囲から得られる潜在変数のパラメータに関するサンプリング集合(各パラメータの組み合わせの集合)が割り当てられた力学モデルを意味するものとする。
図2は、構造流体解析の対象領域(以下、解析対象領域と呼ぶ)に関する力学モデルM1の一例を示す図である。図2に示すように、力学モデルM1は、大動脈領域R1と右冠動脈領域R2と左冠動脈領域R3とを有している。血液は、大動脈から右冠動脈または左冠動脈へ流れる。
図2に示すように、力学モデルM1の大動脈起始部側の末端が血流の入口に設定され、右冠動脈領域の末端と左冠動脈領域の末端とが血流の出口に設定される。入口と出口との各々に境界条件が設定される。入口に関する境界条件は、例えば入口における血流の流速または血流による圧力、またはそれらの変化率を含む。出口に関する境界条件は、例えば、出口における血流の流速または血流による圧力、またはそれらの変化率を含む。
大動脈、右冠動脈、および左冠動脈の変形は、血流に起因する血管壁への力学的作用、心臓の拍動による血管壁への力学的作用(外力)、血管断面境界の負荷条件、血管壁の材料モデル、血管の無応力状態、および血管壁の幾何学的形状等の様々な要因に依存する。
ここで、血流に起因する血管壁への力学的作用は、例えば、血流に起因する内圧と血流に起因するせん断力とを含む。血流に起因する内圧により、血管半径方向あるは血管内腔面の垂直方向に変形する。心臓の拍動による血管壁への力学的作用と、血流に起因するせん断力は、血管芯線方向に関する伸縮、ねじり、曲げといった血管壁への力学的作用による変形を生じさせる。
血管芯線方向に関する伸縮、ねじり、曲げといった血管の全体的および局所的な変形挙動は、大動脈領域R1、右冠動脈領域R2、および左冠動脈領域R3、またはその一部、またはその周辺組織に、強制変位(移動ベクトルや回転変位)または荷重ベクトルの時間的変化を与えることで、負荷条件を設定する。また、血流に起因する内圧に基づく血管半径方向または内腔面に垂直な方向への変形は、血管内腔への圧力分布の時間的変化として与えることで、負荷条件を設定する。
大動脈領域R1、右冠動脈領域R2、および左冠動脈領域R3には、構造流体解析において、強制変位による変位拘束条件が割り当てられる。これにより、構造流体解析における血管壁の変形自由度を縮小することができ、計算収束性の安定化や、解析時間の短縮を実現できる。
また、例えば、血管の形状の変形度合は、血管壁の材料に依存する。そのため、大動脈領域R1、右冠動脈領域R2、および左冠動脈領域R3に材料モデルが割り当てられる。また、例えば、血管の形状の変形度合は、血管の無応力状態に依存する。負荷条件の初期値として、血管の残留応力分布を負荷条件として与えても良い。
ここで、血管解析対象の数値計算用力学モデルの空間を離散化した節点集合および節点から構成される要素集合において、材料モデルや境界条件や負荷条件等の解析条件を同定しない領域と、同定する領域にわけ、解析条件を同定しない領域内の節点には強制変位履歴の変位拘束条件を与え、材料モデルを同定する領域では、血管壁表面(外表面)の節点のみ強制変位履歴の変位拘束条件を与え、血管壁内部は変位自由度を確保し変位拘束を与えない。これにより、構造流体解析の変形自由度をおさえることができ、安定的かつ効率的に解析を行うことができる。
ただし、血管壁表面上に緩衝用のダミーの要素集合を設け、その表面の節点に強制変位を与えてもよい。これにより、血管内腔にプラーク等による突起がある場合や、血管分岐部等、内圧による負荷が、芯線方向の断面外の変形に影響を与える場合には、血管内腔だけではなく血管壁の形態指標も参照して血管に負荷される荷重ベクトルと内圧を分離して同定することができる。また、生理学的には壁表面の脂肪層を模擬しており、一方、数値計算上は、血管壁表面に強制変位を与えることで血管壁内部に局所的に実際とは異なる高い応力が発生することを避ける効果もある。
これら材料モデル、境界条件、および負荷条件等の潜在変数に関するパラメータは、後述する力学モデルに基づく逆解析(統計的同定処理)により同定される。逆解析により同定された精確な潜在変数は、力学モデルに割り当てられる。精確な潜在変数が割り当てられた力学モデルにより、解析対象血管領域外の血管や心臓等の外部要因による当該解析対象血管領域への影響を加味した構造流体解析または流体解析または構造解析または画像解析に基づく血行動態解析を実行することができる。
血管解析装置50は、力学モデルの構築に関し、逆解析による潜在変数の同定により、次の4点の困難を解決することができる。困難1.冠動脈の材料モデルの同定方法。困難2.心臓の形状の変形の冠動脈への影響の組み込み。困難3.冠動脈の境界条件の同定方法。困難4.医用画像データの不確定性に起因したばらつきを有する血管形状による画像解析や構造流体解析。この4点の困難の克服により、血管解析装置50は、逆解析による潜在変数の同定を行わない従来の血管構造流体解析に比して、解析精度の向上を実現する。
次に、実施の形態の医用画像診断装置における構造流体解析処理の詳細について説明する。図3は、システム制御部21の制御のもとに行われる構造流体解析処理の典型的な流れを示す図である。図4は、画像処理装置27の機能ブロック図である。
図3に示すように、構造流体解析処理においては、まず、システム制御部21により記憶装置33から処理対象の医用画像ファイルが読み出され、画像処理装置27に供給される。医用画像ファイルは、時系列のCT画像のデータの他に、当該被検体に関する血管内腔に関する圧力値のデータ、血液流量指標の観測値のデータ、およびプラーク指標を含んでいる。CT画像以外ではMRI画像や超音波エコー画像であってもよい。時系列のCT画像のデータは、時系列のCT値の3次元空間分布を表現するデータである。時系列のCT画像は、例えば1心拍で約20枚、すなわち、約20心位相分のCT画像を含んでいる。
図3に示すように、システム制御部21は、画像処理装置27に領域設定処理を行わせる(ステップS1)。ステップS1において画像処理装置27の領域設定部51は、時系列のCT画像に含まれる血管領域に構造流体解析の解析対象領域を設定する。解析対象領域は、冠動脈に関する血管領域の任意の一部分に設定される。例えば、領域設定部51は、ユーザによる入力機器29を介した指示、または、画像処理により血管領域に解析対象領域と同定対象領域を設定する。
ここで、図5を参照しながら、血管の構造について説明する。図5は、血管の芯線に直交する断面(以下、血管断面と呼ぶ)を模式的に示す図である。図5に示すように、血管は、管状の血管壁を有している。血管壁の中心軸は芯線と呼ばれている。血管壁の内側は内腔と呼ばれている。内腔に血液が流れる。内腔と血管壁との境は血管内壁と呼ばれている。血管壁の外側には心筋等の血管周辺組織が分布している。血管壁と血管周辺組織との境は血管外壁と呼ばれている。血管壁内部にプラークが発生することがある。
図5に示すように、プラークは、例えば、石灰化した石灰化プラーク、粥状プラーク等に分類される。粥状プラークは、やわらかく、血管内壁が破裂して血栓として血管内部に染み出す危険性があり、不安定プラークと呼ばれることもある。従って、プラークの性状を把握することは臨床的に有用である。プラークの性状や存在領域は、医用画像ファイルに含まれるプラーク指標により特定可能である。プラーク指標は、例えば骨のCT値を基準に正規化したCT値の大きさ等により相対的に判別することができる。しかし、血管内部のプラークの変形特性や硬さを解析するのは容易ではない。
ステップS1が行われるとシステム制御部21は、画像処理装置27に画像解析・追尾処理を行わせる(ステップS2)。ステップS2において画像処理装置27の画像解析・追尾処理部53は、時系列のCT画像に画像処理を施して時系列の血管形態指標と時系列の血管形状変形指標とを算出する。具体的には、画像解析・追尾処理部53は、時系列のCT画像に画像解析処理を施すことにより時系列の血管形態指標を算出し、時系列のCT画像に追尾処理を施すことにより時系列の血管形状変形指標を算出する。
より詳細には、画像解析・追尾処理部53は、画像解析処理において、各CT画像から血管領域を抽出し、血管の内腔に関する画素領域(以下、血管内腔領域と呼ぶ)と血管壁に関する画素領域(以下、血管壁領域と呼ぶ)とを特定する。画像解析・追尾処理部53は、血管形態指標として、血管の芯線に垂直な断面、または血管内腔面に垂直な面が、血管内腔、血管壁、プラーク領域に交わる領域上の複数の画素の3次元座標を特定する。
なお、血管形態指標は、3次元座標だけでなく、芯線に垂直な断面における一定角度ごとの血管内腔の半径や直径および0°の方向ベクトル、または断面における全角度に対する平均面積や平均半径、または、芯線方向に垂直な複数の断面で囲まれた血管内腔容積、または内腔面に垂直な複数断面で囲まれた血管壁容積やプラーク容積等の幾何学的指標でも良い。
追尾処理において、画像解析・追尾処理部53は、ユーザからの入力機器29を介した示または画像処理により、血管領域や血液や造影剤やプロトンにおける特徴点や特徴形状、代表点、画素等の複数の追跡点を設定する。例えば、画像解析・追尾処理部53は、血管分岐部や表面の特徴形状等の追跡点集合を設定する。各時刻(各心位相)における画像解析・追尾処理部53の追尾処理により得られた追跡点集合の変位データから、力学モデルの血管壁表面または血管壁内部または血管内腔における節点の変位の時間的変化を補間処理等により算出し、強制変位として与える。
また、例えば画像解析・追尾処理部53は、力学モデルに血管芯線上の節点を定義する。画像解析・追尾処理部53は、力学モデルの血管壁表面または血管壁内部または血管内腔における節点の変位の時間的変化から、血管の芯線方向に関する伸縮やねじりや曲げに関する変形を抽出し、血管芯線と芯線に垂直な断面における節点の強制変位として与えることで表現しても良い。このように、血管形状変形指標としては、力学モデルにおける各時刻の節点の強制変位データ(強制変位履歴)を特定する。
以下、図6および図7を参照しながら、画像解析・追尾処理を説明する。図6は、血管芯線の形態の時間的変化を示す図である。図6に示すように、例えば時系列の医用画像は、1心拍につき20枚のCT画像を含んでいるものとする。すなわち、心位相0%から95%まで5%おきにCT画像が得られることとする。画像解析・追尾処理部53により血管領域の芯線が抽出される。図6に示すように、芯線の形態は、心位相の経過に従って変化している。
図7は、時刻tと時刻t+Δtとの間における追尾処理の一例を示す図である。図7に示すように、血管芯線上にP1からP10の力学モデルにおける節点が設定されており、各断面上の血管の力学モデルの節点と力学的につながっている。ただし、血液の力学モデルの節点とは独立である。血管の追跡点の変位データをもとに、血管芯線上のP1からP20の節点の変位データを補間等の処理により算出し、各節点に強制変位が設定されているものとする。血管形状変形指標と血管形態指標について説明するため、節点P13と節点P14とにより規定される局所血管領域RAを考える。
時刻tにおいて、芯線方向に関する節点P13と節点P14との間の距離がLであり、血管領域の半径がrであるとする。画像解析・追尾処理部53から、節点P13と節点P14の血管芯線方向の伸縮、ねじり、または曲げ等の強制変位を抽出し、節点P13の強制変位(3次元空間における移動変位と、芯線方向の回転変位)と節点P14の強制変位(3次元空間における移動変位と、芯線方向の回転変位)を算出する。
図8は、血管芯線の曲げ変形および回転変位の算出例を示す図である。図8に示すように、例えばねじり角は、線aおよび線bから構成される面の法線方向ベクトルcの変化から算出しても良い。
図7および図8に示すように、画像解析・追尾処理部53は、追跡点の座標と移動ベクトルとに基づいて、芯線上の各節点の強制変位(3次元空間における移動変位と、芯線方向の回転変位)を算出し、血管形状変形指標を算出する。例えば、画像解析・追尾処理部53は、隣り合う2つの節点の座標差の時間変化を芯線方向に関する伸縮距離ΔLとして算出する。
また、画像解析・追尾処理部53は、各芯線上の節点について、当該節点と当該節点を含む血管領域断面上の他の節点(血管内腔または血管壁またはプラーク領域における節点)との間の距離の時間変化を半径方向に関する伸縮距離Δrとして算出する。また、画像解析・追尾処理部53は、各追跡点について、当該追跡点の近傍の複数の追跡点の座標と移動ベクトルとに基づいて、芯線上の当該節点の芯線方向のねじれ角度Δθを算出する。
また、画像解析・追尾処理部53は、血液領域の造影剤やプロトンの画像追尾により、血液流量指標として、流速または芯線方向断面の平均流速または平均流量を算出してもよい。
血管形状変形指標は、力学モデルにおける強制変位として利用される。以下、時系列の血管形態指標を「形状履歴」と呼び、時系列の血管形状変形指標を「強制変位履歴」と呼ぶことにする。
ステップS2が行われるとシステム制御部21は、画像処理装置27に構築処理を行わせる(ステップS3)。ステップS3において画像処理装置27の力学モデル構築部55は、形状履歴(時系列の血管形態指標)と、強制変位履歴(時系列の血管形状変形指標)と、時系列の医用画像(CT画像,MRI画像,超音波エコー画像等のDICOMデータ)とに基づいて、解析対象領域に関する力学モデルを暫定的に構築する。力学モデルは、構造流体解析を行うための解析対象領域に関する数値モデルである。
以下、ステップS3について詳細に説明する。力学モデル構築部55は、まず、医用画像と形状履歴とに基づいて、力学モデル(数理モデル)を解くための形状モデルを構築する。形状モデルは、各時刻における血管領域の幾何学的構造を模式的に表現したものである。形状モデルは、例えば複数の離散化領域に区分されている。各離散化領域の頂点は、節点と呼ばれる。
力学モデル構築部55は、時刻毎の医用画像に含まれる血管領域と血管形態指標とに基づいて、時刻毎の形状モデルを構築しても良いし、特定の時相の医用画像に含まれる血管領域と血管形態指標とに基づいて、時刻毎の形状モデルを構築しても良い。また、例えば、初期の負荷状態として、解析対象領域に対応する血管に残留応力が存在しないと仮定する場合、無応力状態の時相として、解析対象領域に対応する血管が最も収縮した時相を無応力状態であると仮定する。
図9は、形状モデルの芯線に直交する断面を示す図である。図9に示すように、形状モデルは、芯線から外側に向けて血管内腔領域、血管壁領域を有している。プラークが存在する場合、血管壁領域にプラーク領域を設けても良い。また、血管周辺組織による血管への影響を考慮する場合、血管周辺組織のダミー要素を血管壁領域の外側に設けても良い。
形状モデルが構築されると力学モデル構築部55は、各潜在変数の確率分布や変数範囲から得られる潜在変数のパラメータに関するサンプリング値(例えばマルコフ連鎖モンテカルロ法等による、各パラメータの組み合わせの集合からのサンプリング)を力学モデルに設定する。
例えば、力学モデル構築部55は、図2に示すように、大動脈領域R1の大動脈起始部側の末端に入口に関する境界条件の同定対象の領域(以下、境界条件同定領域)を設定し、右冠動脈領域R2の末端と左冠動脈領域R3の末端とに出口に関する境界条件同定領域を設定する。力学モデル構築部55は、各境界条件同定領域に境界条件の確率分布または変数範囲から得られる境界条件のパラメータに関するサンプリング値を割り当てる。
また、力学モデル構築部55は、大動脈領域R1、右冠動脈領域R2、および左冠動脈領域R3に材料モデルの同定対象の領域(以下、材料モデル同定領域と呼ぶ)および負荷条件の同定対象の領域(以下、負荷条件同定領域と呼ぶ)を設定する。力学モデル構築部55は、各材料モデル同定領域に材料モデルの確率分布または変数範囲から得られる材料モデルのパラメータに関するサンプリング値を割り当て、各負荷条件同定領域に負荷条件の確率分布や変数範囲から得られる負荷条件のパラメータに関するサンプリング値を割り当てる。
血管は、流量が0でも残留応力が存在すると言われている。例えば、力学モデル構築部55は、流量が0の場合の残留応力を負荷条件の初期値として解析対象領域に割り当てても良い。また、力学モデル構築部55は、幾何学的構造に不確定性がある部位に、幾何学的構造の同定対象の領域(以下、幾何学的構造同定領域と呼ぶ)を設定しても良い。なお、幾何学的構造のパラメータは、幾何学的構造の不確定性に関連するばらつき分布パラメータ、または医用画像データに内在するばらつき分布パラメータであり、各CT値のばらつき分布や生体組織の境界閾値のばらつき分布等であってよい。詳細は、後述するが、力学モデル構築部55は、プラーク領域に材料モデルを設定しても良い。材料モデルの詳細については後述する。
形状モデルが構築されると力学モデル構築部55は、ステップS2において算出された時系列の血管形状変形指標、すなわち、強制変位履歴を形状モデルに割り当てる。潜在変数および強制変位履歴が割り当てられた形状モデルを「力学モデル」と呼ぶことにする。
図10は、形状モデルM2は血管や血液の力学モデルの一部を示しており、力学モデルの節点への強制変位履歴の割り当てを説明するための図である。図10は、形状モデルM2の一部分を示している。ただし、図10では芯線がM2内に位置する場合を示しているが、芯線がM2外に位置する場合でも良い。図10に示すように、形状モデルM2には複数の節点PN(PN1,PN2)が設定されている。芯線上の節点をPN1と称し、血管や血液の力学モデルにおける節点をPN2と称することにする。形状モデルM2は、ダミー要素表面、血管外壁、血管内壁、プラーク領域表面、プラーク領域内部、または血液部、に設定される。力学モデル構築部55は、形状モデルM2の各節点PN1に強制変位、すなわち、血管形状変化指標を時刻毎に割り当てる。
具体的には、力学モデル構築部55は、芯線上の隣り合う節点PN1と節点PN1とをビーム要素(またはリジッド要素)EBで結ぶ。また、力学モデル構築部55は、節点PN1と当該節点PN1を通る直交断面に含まれる他の節点PN2とをビーム要素EBで結ぶ。力学モデル構築部55は、節点PN1およびビーム要素EBに各血管形状変形指標の形状変位方向に関する拘束条件を割り当てる。
材料モデルまたは血管内腔の内圧を同定する領域では、強制変位としては、芯線方向に関する血管壁(またはダミー要素)表面の伸縮、血管壁(またはダミー要素)表面のねじれ、血管壁(またはダミー要素)表面の曲げ変形が挙げられる。例えば、材料モデルや血管内腔の内圧を同定しない領域では、芯線方向に関する強制変位だけではなく、半径方向に関する血管壁の時系列の伸縮(変位)も強制変位履歴として割り当てられる。
また、血管内腔に突起がある場合や血管分岐部等のように、内圧が、芯線方向断面外への変形に寄与する場合は、その領域には強制変位を与えずに、その領域の周辺部(例えば、ダミー要素の表面節点)のみに強制変位履歴を割り当てる。また、力学モデル構築部55、節点PN1およびビーム要素EBに時系列の血管形状変形指標を強制変位履歴として割り当てられる。これにより、血管の全体および局所に関する伸縮変形やねじれ変形や曲げ変形が表現される。
なお、強制変位履歴の割当対象は、芯線上の節点およびビーム要素に限定されない。図11および図12は、形状モデルへの強制変位履歴の他の割り当て方法を示す図である。図11は、血管形状変形指標がねじれの場合における割り当て例を示している。また、図12は、血管形状変形指標が曲げの場合における割り当て例を示している。図11および図12に示すように、力学モデル構築部55は、形状モデルの表面または内部の節点に直接的に強制変位履歴を割り当てても良い。例えば、ステップS2において画像解析・追尾処理部53は、特徴点についての伸縮量やねじり量や曲げ量等の血管形状変形指標を算出する。力学モデル構築部55は、算出された血管形状変形指標を特徴点の周囲の節点へ補間(内挿や外挿)することにより直接的に割り当てる。
実施の形態の画像処理装置27は、ステップS3において暫定的に構築された力学モデルを用いて逆解析を施し、力学モデルに設定される潜在変数を統計的に同定する。統計的同定処理は、後述のステップS6において行われる。ステップS4およびS5は、それぞれ統計的同定処理に用いる血管形態指標および血液流量指標を算出するために設けられている。
ステップ3が行われるとシステム制御部21は、画像処理装置に27に血管応力解析処理を行わせる(ステップS4)。ステップS4において画像処理装置27の血管応力解析部57は、現段階の力学モデルに血管応力解析を施して時系列の血管形態指標の予測値を算出する。血管形態指標は、既述の血管形態指の何れであっても良いが、例えば血管芯線方向に関する内腔領域の断面形状指標や血管壁の断面形状指標が用いられると良い。
具体的には、内腔領域の断面形状指標は、内腔領域の注目画素の座標値、内腔領域の幾何学的構造パラメータ(内腔領域の半径、内腔領域の直径等)の少なくとも一つである。また、血管壁領域の断面形状指標は、具体的には、血管壁領域の注目画素の座標値、血管壁領域の幾何学的構造パラメータ(血管壁領域の半径、壁領域の直径等)の少なくとも一つである。なお、予測値は、力学モデルに血管応力解析を施して算出された血管形態指標の算出値を意味する。
また、ステップ3が行われるとシステム制御部21は、画像処理装置27に血液流体解析処理を行わせる(ステップS5)。ステップS5において画像処理装置27の血液流体解析部59は、暫定的に構築された力学モデルに血液流体解析を施して時系列の血液流量指標の予測値を算出する。血液流量指標は、血流量または流速、またはその空間的・時間的な平均値の少なくとも一つである。なお、予測値は、力学モデルに血液流体解析を施して算出された血液流体指標の算出値を意味する。
ステップS4およびS5が行われるとシステム制御部21は、画像処理装置27に同定処理を行わせる(ステップS6)。ステップS6において画像処理装置27の統計的同定部61は、ステップS4において算出された血管形態指標の予測値とステップS5において算出された血液流量指標の予測値とが、事前に収集された血管形態指標の観測値と血液流量指標の観測値とに整合するように力学的モデルの潜在変数のパラメータを統計的に同定する。
図4に示すように、統計的同定部61は、第1統計的同定部61−1と第2統計的同定部61−2とを有している。第1統計的同定部61−1は、血管形態指標の予測値が血管形態指標の観測値に整合するように潜在変数のパラメータを統計的に同定する。第2統計的同定部61−2は、血液流量指標の予測値が血液流量指標の観測値に整合するように潜在変数のパラメータを統計的に同定する。以下、第1統計的同定部61−1と第2統計的同定部61−2とを順番に説明する。
具体的には、ステップS6において第1統計的同定部61−1は、ステップS4において算出された血管形態指標の予測値と観測値とに基づくデータ分布を設定する。データ分布は、例えば、血管形態指標の予測値と観測値との誤差に関する多変量正規分布関数を示す。
力学モデルにおける各節点または各要素での予測値と観測値の誤差に関する正規分布関数値を算出し、それらの積をデータ分布として統計的同定処理部において設定する。データ分布は、時刻毎に個別に設定されても良いし、複数時刻まとめて設定されても良い。
次に、第1統計的同定部61−1は、力学モデルの潜在変数に事前分布(事前確率分布)を割り当てる。具体的には、材料モデル、境界条件、負荷条件、および時系列の形態指標や形状変形指標の不確定性に関連するパラメータの各々に事前分布が割り当てられる。例えば、負荷条件のパラメータの一つである血管内腔に関する圧力値に関する事前分布が割り当てられる。圧力値の取り得る値の範囲(想定範囲)は、経験的に予め限定することができる。第1統計的同定部61−1は、これら想定範囲内に限定して内圧値に関するモンテカルロシミュレーションを実行することにより、各離散化領域について内圧値の確率分布、すなわち、事前分布を算出する。
また、力学モデル構築部55は、事前分布として、芯線方向の圧力分布は滑らかであること、また、時間経過に伴う圧力変化が滑らかであること、血流の逆流がないことが観測されている場合には芯線方向の平均的な圧力変化の傾きが負であることを、例えば多変量正規分布関数により数学的に表現した確率分布を事前分布として設定しても良い。
これら想定範囲内に限定して、想定した確率分布に従って、負荷条件のパラメータに関するモンテカルロシミュレーションを実行することができ、力学モデルへ設定するための負荷条件(潜在変数)のサンプリング値を得ることができる。
次に、第1統計的同定部61−1は、各潜在変数について、事前分布とデータ分布とに統計的同定処理を施すことにより事後分布(事後確率分布)を算出する。統計的同定処理は、例えば、階層ベイズモデルやマルコフ連鎖モデルが挙げられる。そして、第1統計的同定部61−1は、各潜在変数について、事後分布の最頻値や平均値等の統計値から各潜在変数のパラメータを同定する。例えば、上述の例の場合、血管内腔圧力値に関する事後分布が算出され、この事後分布から血管内腔圧力値の同定値が算出される。
図13は、階層ベイズモデルおよびマルコフ連鎖モンテカルロ法による負荷条件(血管内の平均圧力)に関する事後分布算出と平均内圧の同定とを説明するための図である。図13に示すように、血管起始部から延びる血管に石灰化プラーク領域と粥状プラーク領域とが設定されているものとする。石灰化プラーク領域は材料モデル同定領域に設定され、粥状プラーク領域に材料モデル同定領域に設定される。血管起始部から血管芯線方向に沿って進むにつれて血管内圧が降下する。血管芯線に沿って複数の節点が設定される。各節点を含む直交断面(節点断面)において内腔内圧の事後分布が算出され、事後分布の最頻値が同定される。
なお、血管形態指標の観測値としては、例えば、ステップS2において算出された血管形態指標が用いられる。
第2統計的同定部61−2による処理は、データ分布の算出に用いる指標が異なるだけで第1統計的同定部61−1による処理と同様である。すなわち、第2統計的同定部61−2は、まず、ステップS5において算出された血液流量指標の予測値と観測値とに基づくデータ分布を設定する。
次に、第2統計的同定部61−2は、力学モデルの潜在変数に事前分布を割り当てる。例えば、血管に関する材料モデルのパラメータや血液に関する材料モデルのパラメータ、プラークに関する材料モデルのパラメータに関する事前分布が割り当てられる。これら材料モデルのパラメータとしては、例えば、弾性率等の材料モデルパラメータや、血液の構成式における粘性に関するパラメータが挙げられる。材料モデルのパラメータの想定範囲や確率分布は、経験的に予め設定することができる。
第2統計的同定部61−2は、各離散化領域について材料モデルのパラメータの確率分布、すなわち、事前分布を設定し、これら想定範囲内に限定して、想定した確率分布に従って、材料モデルのパラメータに関するモンテカルロシミュレーションを実行することができ、力学モデルへ設定するための材料モデルパラメータ(潜在変数)のサンプリング値を得ることができる。
次に、第2統計的同定部61−2は、各潜在変数について、事前分布とデータ分布とに統計的同定処理を施すことにより事後分布を算出し、算出された事後分布の統計値から各潜在変数のパラメータを同定する。例えば、上述の例の場合、材料モデルのパラメータに関する事後分布が算出され、この事後分布から材料モデルのパラメータの同定値が算出される。
図14は、階層ベイズモデルおよびマルコフ連鎖モンテカルロ法による材料モデルパラメータに関する事後分布算出と材料モデルパラメータ(血管壁の等価弾性率)の同定とを説明するための図である。図14に示すように、血管モデルは、図13と同様であるとする。材料モデル同定領域に限定して、血管壁の材料モデルのパラメータ(例えば、等価弾性率)の事後分布が算出され、事後分布の最頻値が同定される。
なお、血液流量指標の観測値は、例えば、大動脈に送り出される血流量変化であると仮定し、血管形態指標の観測値を、時系列のCT画像から画像処理により計測される左心室の容積変化値(CFA)として用いることができる。造影剤の冠動脈内注入後の造影剤の画像追尾により特徴点の移動量の時間的変化を算出することにより、流速や流量を算出してもよい。また、造影剤の血管芯線方向または時間的な特定領域の濃度変化量を取得し、その濃度変化を各領域の芯線方向距離間隔で除した値や、濃度変化の時間的変化率から、流速や流量を算出してもよい。MRIの場合はプロトンの画像追尾を用い、超音波エコーの場合には、コントラストエコー図法等により流量を算出する。
また、解析対象領域の各画素の座標値が確定値であることを前提としない場合、すなわち、解析対象領域の幾何学的構造に不確定性があると仮定する場合、潜在変数に幾何学的構造を含めても良い。この場合、統計的同定部61は、各節点の座標値の芯線方向に関する所定範囲内の変動、または、解析対象領域の径の所定範囲内の変動を表現した正規分布等の確率分布を事前分布として設定すると良い。この場合、解析対象領域の形状が滑らかであること、また、芯線における節点の順番が不変であるという制約も事前分布として設定しても良い。これら想定範囲内に限定して、想定した確率分布に従って、幾何学的構造のパラメータに関するモンテカルロシミュレーションを実行することができ、力学モデルへ設定するための幾何学的構造の不確定性パラメータ(潜在変数)のサンプリング値を得ることができる。
なお、各ステップS6において第1統計的同定部61−1による統計的同定処理と第2統計的同定部61−2による統計的同定処理との両方が行われなくても良い。すなわち、各ステップS6においては、第1統計的同定部61−1による統計的同定処理と第2統計的同定部61−2による統計的同定処理との何れか一方が行われても良い。
また、上記の例においては、第1統計的同定部61−1は、血管形態指標の予測値が血管形態指標の観測値に整合するように潜在変数のパラメータを統計的に同定し、第2統計的同定部61−2は、血液流量指標の予測値が血液流量指標の観測値に整合するように潜在変数のパラメータを統計的に同定するとした。しかし、統計的同定部61は、構造−流体連成解析に基づいて、血管形態指標の予測値と血液流量指標の予測値とが血管形態指標の観測値と血液流量指標の観測値とに整合するように潜在変数のパラメータを統計的に同定しても良い。統計的同定部61による統計的同定処理のさらなる詳細については後述する。
ステップS6が行われるとシステム制御部21は、画像処理装置27に設定処理を行わせる(ステップS7)。ステップS7において画像処理装置27の力学モデル構築部55は、ステップS6において算出された潜在変数のパラメータを力学モデルに設定する。
ステップS7が行われるとシステム制御部21は、同定終了条件が満たされたか否かを判定する(ステップS8)。ステップS8において同定終了条件が満たされていないと判定した場合(ステップS8:NO)、システム制御部21は、ステップS4、S5、S6、S7、およびS8を繰り返す。
ここで、同定終了条件は、同定終了を判定するための指標(以下、同定終了指標と呼ぶ)が規定値に達するか否かにより表現される。同定終了指標としては、例えば、血管形態指標の予測値と観測値との差分値が挙げられる。この場合、システム制御部21は、この差分値が既定値よりも大きい場合、同定終了条件が満たされていないと判定し、差分値が既定値よりも小さい場合、同定終了条件が満たされたと判定する。
また、同定終了指標は、例えば、モンテカルロ法のサンプリング点の数でも良い。この場合、システム制御部21は、このサンプリング点の数が既定値よりも小さい場合、同定終了条件が満たされていないと判定し、サンプリング点の数が既定値よりも大きい場合、同定終了条件が満たされたと判定する。同定終了条件が満たされた場合、力学モデル構築部55は、その時点の最新の力学モデルを最終的な力学モデルに設定する。
最終的な力学モデルが構築されると、力学モデル構築部55は、血管形状変形指標の観測値、最終的な力学モデルにおける負荷条件のパラメータ、および材料モデルのパラメータを関連付けたモデル(以下、関連モデルと呼ぶ)を算出する。関連モデルは、記憶装置33に記憶される。関連モデルは、検索の容易性等のため患者情報や検査情報等に関連付けて記憶されると良い。なお、血管形態指標や血液流量指標の観測値、最終的な力学モデルにおける負荷条件のパラメータ、および材料モデルのパラメータは、必ずしもモデルの形態で関連付けられる必要はなく、例えば、テーブルまたはデータベースであっても良い。
上記のステップS4、S5、S6、S7、およびS8は、同一の同定法で反復しても良いし、異なる同定法で反復しても良い。異なる同定法で反復する場合、例えば、まず、簡易的力学モデルを利用して潜在変数を暫定的に同定し、次に、連続体力学モデルを利用して潜在変数を正確に同定しても良い。このように統計的同定処理を異なる手法で2段階に分けて行うことにより、潜在変数のパラメータを短時間で収束させることができる。簡易的力学モデルを利用する方法としては、内圧および外圧を厚肉円筒の材料力学の式と、ハーゲン・ポアズイユ流れおよび修正ベルヌーイの式とを用いる方法が挙げられる。連続体力学モデルを利用する方法としては、FEM構造流体解析が挙げられる。簡易的力学モデルを利用する同定法と連続体力学モデルを利用する同定法との詳細については後述する。
ステップS8において同定終了受件が満たされたと判定した場合(ステップS8:YES)、システム制御部21は、画像解析・追尾処理部53に修正処理を行わせても良い(ステップS9)。ステップS9において画像解析・追尾処理部53は、統計的同定法による逆解析で得られた潜在変数のもとで実施した構造流体解析結果(力学的指標の予測値および血液流体指標の予測値)が観測値(力学的指標の観測値および血液流体指標の観測値)に整合するように、時系列の医用画像に含まれる血管領域の形状を修正しても良い。
表示機器31は、修正後の時系列の医用画像に基づく診断結果を表示する。これにより、血管解析装置50は、最終的な力学モデルを考慮した診断結果を表示することができる。または、表示機器31は、逆解析による同定とその構造流体解析により観測結果とが整合しない血管箇所・領域を画面に表示しても良い。
例えば、血管の挙動の動きの速い心位相の画像はぼやけることが多く、医用画像をもとにした画像解析により観測した血管形状には誤差が大きい箇所や領域が存在する。比較的、血管の挙動が安定した心位相のデータでは、ノイズが少ない画像から得られる。そのような誤差分布の小さい血管形状データに基づき、誤差が大きな心位相での血管形状を、力学モデルを用いることで、正しく内挿することができ、そのような誤差が大きい血管箇所や領域について、正しく内挿した形状とともに、元データからのばらつき分布とともに、表示することができる。これにより、血管形状表示の安定性を確保できるとともに、形状の不確定性をユーザが認識できる。
ステップS9が行われるとシステム制御部21は、画像処理装置27に血管応力解析処理を行わせる(ステップS10)。ステップS10において画像処理装置27の血管応力解析部57は、最終的な力学モデルに血管応力解析を施して時系列の力学的指標の予測値の空間分布を算出する。具体的には、離散化領域毎に力学的指標の予測値が算出される。
また、ステップS9が行われるとシステム制御部21は、画像処理装置27に血液流体解析処理を行わせる(ステップS11)。ステップS11において画像処理装置27の血液流体解析部59は、暫定的に構築された力学モデルに血液流体解析を施して時系列の血液流量指標の予測値の空間分布を算出する。具体的には、離散化領域毎に血液流量指標の予測値が算出される。
なお、力学的指標または血液流量指標として、FFRが算出されても良い。
ステップS9およびS10が行われるとシステム制御部21は、表示機器31に表示処理を行わせる(ステップS11)。ステップS11において表示機器31は、ステップS9において算出された時系列の力学的指標の予測値とステップS10において算出された時系列の血液流量指標の予測値とを表示する。
例えば、表示機器31は、時系列の力学的指標または時系列の血管流量指標を、時系列の力学モデルを当該予測値に応じた色で動画的に表示する。このため、表示機器31は、各種の予測値とカラー値(例えば、RGB)との関係を示すカラーテーブルを保持している。表示機器31は、予測値に応じたカラー値をカラーテーブルを利用して特定し、特定されたカラー値に応じた色で当該予測値に対応する離散化領域を表示する。
図15は、力学的指標の一つである内圧の空間分布の表示例を示す図である。図15に示すように、表示機器31は、力学モデルを構成する各離散化領域を当該離散化領域に関する内圧値に応じた色で動画的に表示する。ユーザは、力学モデルを観察することにより、経時的且つ空間的に変化する力学的指標を色で把握することができる。
図16は、血液流量指標の一つである流速値の空間分布の表示例を示す図である。図16に示すように、表示機器31は、力学モデルを構成する各離散化領域を当該離散化領域に関する流速値に応じた色で動画的に表示する。ユーザは、力学モデルを観察することにより、経時的且つ空間的に変化する血液流量指標を色で把握することができる。
例えば、血管内が完全に狭窄している場合、狭窄部位の内圧は、非狭窄部位の内圧よりも小さい。力学的指標として内圧が指定された場合、ユーザは、力学モデル上での局所的な色の違いにより狭窄の有無を判断することができる。また、狭窄部位の流量は、非狭窄部位の流用よりも小さい。血液流量指標として流量が指定された場合、ユーザは、力学モデル上での局所的な色の違いにより狭窄の有無を判断することができる。
また、血管応力解析部57は、プラーク領域の材料モデルパラメータの同定結果に基づいて、力学的指標として硬さ値の空間分布を算出しても良い。この場合、表示機器31は、プラーク領域に関する硬さ値の空間分布を力学モデル上において表示しても良い。また、表示機器31は、プラーク領域周辺の内圧分布や応力分布、ひずみ分布を表示して良い。ユーザは、これらの表示をプラークの性状と破綻しやすさとを推定することに活用することができる。
力学的指標および血液流量指標の予測値は、力学モデルの離散化領域の色で表現する方法のみに限定されない。例えば、図17、図18、および図19に示すように、グラフで表示しても良い。なお、図17は、左冠動脈起始部の血圧に関するグラフである。図17のグラフの縦軸は正規化した血圧に規定され、横軸は心位相[%]に規定される。図18は、LCXとLDAとの分岐点付近の血圧に関するグラフである。図18のグラフの縦軸は正規化した血圧に規定され、横軸は心位相[%]に規定される。図19は、芯線方向に関する血圧変化に関するグラフである。図19のグラフの縦軸は血圧比に規定され、横軸は大動脈からの距離[mm]に規定される。表示機器31は、力学的指標および血液流量指標の予測値をグラフで表示することにより、これら値をユーザに簡便に把握させることができる。
ステップS11が行われると構造流体解析処理が終了する。
なお、図10において、強制変位履歴は、形状モデルの芯線部と外壁部とに設定されるとしたが、強制変位履歴の設定箇所は、これに限定されない。例えば、強制変位履歴は、芯線部と外壁部との間の血管壁領域に設定されても良い。
また、強制変位履歴の拘束条件の割り当て対象は、境界条件および材料モデルを同定するか否かに応じて切り分けられても良い。図20は、強制変位履歴の他の割り当て例を示す図であり、形状モデルの断面を示している。例えば、図20の(a)に示すように、境界条件および材料モデルを同定する場合、形状モデルの外壁部OW上の節点PN2のみに強制変位履歴を割り当て、血管壁領域RVの節点PN3には強制変位履歴を割り当てなければよい。
また、図20の(b)に示すように、境界条件および材料モデルを同定しない場合、形状モデルの外壁部の節点PN2と血管壁領域RVの節点PN3との両方に強制変位履歴を割り当てればよい。この場合、芯線上の節点PN1に強制変位履歴が割り当てられる。また、外壁部OWの節点PN2と節点PN1とをビーム要素EBで結び、ビーム要素EB上の節点PN2およびPN3にも強制変位履歴を割り当ててもよい。このとき、半径方向に関する収縮および膨張は、ビーム要素EBの伸縮変位で表現する。なお、血管内腔領域RIには、強制変位履歴が割り当てられなくて良い。
図21は、強制変位履歴の割り当ての他の例を示す図であり、血管周辺組織のダミー要素RDを含む形状モデルの断面を示している。図21に示すように、ダミー要素RDは、血管壁領域RNの外側に設定される。形状モデルがダミー要素RDを含む場合、血管壁領域RNに加え、ダミー要素RDにも節点PN4が設定される。節点PN4にも強制変位履歴が割り当てられる。
力学モデル構築部55は、境界条件および材料モデルを同定する場合、血管壁領域RVに含まれる節点PN3に強制変位履歴を割り当て、境界条件および材料モデルを同定しない場合、強制変位履歴を割り当てなくても良い。節点PN3に強制変位履歴を割り当てる場合、内腔領域RIに関する形状指標以外にも血管壁領域RVに関する形状指標も参照して材料モデルの同定が行われる。
図22は、強制変位履歴の割り当ての他の例を示す図であり、プラーク領域RPを含む形状モデルの断面を示している。図22に示すように、プラーク領域RPは、血管壁領域RVに含まれる。プラーク領域RPは、材料モデル同定領域に設定される。プラーク領域RPについては、内腔形状指標、血管壁形状指標、およびプラーク指標を考慮して材料モデルが同定される。
既述のように、プラーク指標は、例えば、超音波診断装置による組織性状診断により得られたプラークの性状に関するデータである。力学モデル構築部55は、性状に応じてプラーク領域を複数の部分領域に区分し、複数の部分領域に個別に材料モデル同定領域を設定する。各部分領域には、予め、当該部分領域の性状に応じたパラメータ範囲を設定することが好ましい。既述の統計的同定処理により、各部分領域についての材料モデルパラメータが統計的同定部61により同定される。そして、ステップS10において表示機器31が力学的指標として血管の材料特性に関する指標を表示することにより、ユーザは、プラークの性状を正確且つ容易に把握することができる。
次に、潜在変数の一つである材料モデルについて詳細に説明する。血管の材料モデルとしては、弾性モデル、超弾性モデル、異方性超弾性モデル、粘性特性を考慮した超弾性モデル等が適用可能である。異方性超弾性モデルとしては、例えば、Y.C.Funにより提案された数理モデルや、Holzapfel−Gasser構成式と呼ばれる数理モデルを適用できる。
単位参照体積あたりのひずみエネルギーは、以下の(1)式で表わされる。(1)式の第1項は、コラーゲンを含まない等方性基礎材料のせん断変形に関するエネルギーを表している。また、(1)式の第2項は、コラーゲンを含まない等方性基礎材料の体積変形に関するエネルギーを表している。また、(1)式の第3項は、コラーゲン繊維各グループの寄与(繊維方向の分散を考慮)を表している。
図23は、繊維グループの変形を示す図である。図23に示すように、円筒形上の外膜を仮定する。芯線方向zと周方向θとにより規定される面における平均方向Aでの繊維グループの変形は、下記の(2)式で表される。
(1)式および(2)式中の材料モデルに関するパラメータとしては、以下の表1に示すように、材料パラメータと繊維分散パラメータとがある。材料パラメータとしては、C10、D,K1、K2等が用いられ、繊維分散パラメータとしては、Kappaやγ等が用いられる。各パラメータのデフォルト値と制約条件とは表1に示す通りである。
血液の材料モデルは、以下の(3)式のようなCasson構成式や(4)式のようなHB構成式が好適である。
これら材料モデルのパラメータは、上述のステップS6において統計的同定部61により、血管形態指標および血液流体指標を用いた統計的同定処理により同定される。
次に、統計的同定部61により行われる統計的同定処理の詳細について説明する。
時系列の医用画像から計測される血管形態指標および血液流量指標のような観測変数は、不確定性を有している。統計的同定部61は、このような不確定性が存在する状況下における潜在変数の統計的同定法として、階層ベイズモデルとマルコフ連鎖モンテカルロ法とに基づく統計的手法を活用する。
既述のように、ステップS6において統計的同定部61は、ステップS4において算出された血管形態指標または血液流量指標の予測値と観測値とに基づくデータ分布を設定する。データ分布は、例えば血管形態指標または血液流量指標の予測値と観測値との誤差に関する多変量正規分布関数を示す。データ分布は、時刻毎に個別に設定されても良いし、複数時刻まとめて設定されても良い。
次に、統計的同定部61は、形状モデルの強制変位と潜在変数とに事前分布を割り当てる。事前分布は、取り得る値の確率分布を示す。次に、統計的同定部61は、潜在変数に関する数値シミュレーションのパラメータサーベイを実行し、潜在変数と血管形態指標または血液流量指標との関係を表現するモデルを構築する。例えば、材料モデルパラメータと内圧分布パラメータと血管形態指標または血液流量指標との関係がモデルに規定される。なお、血管形態指標または血液流量指標と潜在変数との関係は、モデルという形態ではなく、データベースまたはテーブルにより規定されても良い。これらモデル、データベース、または、テーブルは、記憶装置33に記憶される。
統計的同定部61は、モデル、データベース、または、テーブルを利用して事前分布から血管形態指標または血液流量指標の確率分布を算出する。統計的同定部61は、階層ベイズモデルとマルコフ連鎖モンテカルロ法とにより得られる事後分布から潜在変数を統計的に同定する。
具体的には、この同定問題は、次の3条件を満足しない不良設定問題となる。3条件は、(1)解の存在が保証される、(2)解が唯一に定まる、(3)解がデータに対して連続的に変化し測定誤差に対して解が安定している、である。不良設定問題は、正規化理論とその拡張とからなる枠組み内で捉えると扱いやすい。不良設定問題を解くことは、標準正規化理論では不十分である。不良設定問題を解くためには、内部状態の不連続点を検出し、検出された不連続点を内部状態の推定に役立てる方法論が必要となる。この点で、本実施形態の統計的同定処理においてマルコフ確率場理論が有効となる。
血管形態指標と血液流量指標とに不確定性が存在する環境下において統計的同定部61は、適切な制約条件の下において潜在変数の確率分布パラメータを同定する。適切な制約条件を決定するためには、解の性質を事前に知っている必要がある。統計的同定部61は、解空間の制約条件に関するデータベースをシミュレーションと観測値とに基づいて、発生する。
統計的同定部61は、発生されたデータベースを利用して、超多自由度大規模問題に対してマルコフ確率場理論と階層ベイズモデルとに基づく統計的同定処理を実行する。制約条件となる事前分布の設定では、多くの数値実験結果に基づいて、これらの要因に関するパラメータの確率分布が並列的に個別に構成される。統計的同定部61は、複数の確率分布を統合しデータの欠損を補間することで潜在変数のパラメータを同定する。この処理のために、統計的同定部61は、マルコフ確率場理論を用いたモデルに基づく階層ベイズ法による推定を行う。解析対象とする構造の変形状態の実測結果から、同定した中間変数をもとに、ある負荷条件と境界条件における圧力や流量分布が推定できるという仕組みである。
冠動脈の構造流体解析における材料モデル、境界条件、および負荷条件の同定問題は、非線形逆解析と位置付けられ、解の一意性および安定性が保証されない場合が多い。生体組織の材料特性や血圧の現実的に取り得る範囲は先験情報として想定できるため、これらを事前分布の確率分布として設定できる。また、圧力や変位は空間的時間的に滑らかであることも想定できるためこの情報も先験情報として事前分布の確率分布として設定できる。
または、血液の流れに逆流が生じていない事実を考慮できる場合は血管芯線方向の全体的な圧力分布の傾きは負(圧力降下が存在)であることも制約条件として用いて良い。負荷条件(内圧分布等)、境界条件、および材料モデルに対して、時系列のCT画像に基づく血管形状変形指標の観測値と力学モデルに基づく血管形状変形指標の予測値との2乗誤差分布をデータ分布として設定できる。
観測可能な平均流量に関する2乗誤差分布もデータ分布として追加してもよい。これらの事前分布とデータ分布とに基づいて、階層ベイズモデルとモンテカルロ法とを利用して事後分布を算出することができる。事後分布の発生確率や分散により、潜在変数のパラメータの同定値を得ることができる。発生確率が高く、分散が小さいほうが確信度合の高い同定値といえる。
事後分布が多峰性分布となる場合でも、複数の同定値のうちの分散が小さい同定値を選択すれば良い。または、複数の同定値が存在し得る場合、それぞれの同定条件で、構造流体解析を実施し、それぞれの可能性を認識して、同定値や解析結果を、診断や予防の指針情報として活用することができる。時系列のCT画像にも誤差が含まれていることから、力学モデルの各節点に関する血管形態指標にも誤差が含まれる。このため、各血管形態指標を、例えば、時系列のCT画像から計測された血管形態指標の予測値を平均値とした正規分布の確率変数として扱い、位置の空間的順序を保つという制約を含めた上で、事前分布を設定してもよい。
また、潜在変数のパラメータの同定において、一意性がなく、複数の候補が考えられる場合がある。この場合、時系列のCT画像から計測された血管形態指標の不確定性に従う乱数のサンプリング点に対する、潜在変数の同定値のサンプル集合の変動幅をチェックすることで、各同定値の候補のロバスト性(安定性)を判定する。各同定値の候補のロバスト性に基づいて、最終的な同定値を決定しても良い。
次に、力学モデルの詳細について説明する。既述のように、力学モデル構築部55は、力学モデルの種類に応じて異なるタイプの力学モデルを構築することができる。連続体力学に基づくFEMを用いる場合、力学モデル構築部55は、血管壁の応力解析用のための形状モデル(FEMモデル)と血液の流体解析用のための形状モデル(FEMモデル)との両方を構築する。
材料力学に基づく簡易的な同定法を用いる場合、材料力学における内圧を受ける厚肉円筒の式から近似的に圧力と弾性率と変位との関係を求める。この場合、芯線方向に配列された複数の離散化領域の各々について厚肉円筒近似が形状モデルとして用いられる。具体的には、力学モデル構築部55は、芯線上に離散的に配列された節点を通る断面上の血管内腔形状と血管壁表面形状と断面中心とを特定する。
次に力学モデル構築部55は、血管内腔形状と血管壁表面形状とに基づいて、平均面積、内腔の平均半径、および平均壁厚を算出する。そして力学モデル構築部55は、平均面積、内腔の平均半径、および平均壁厚に基づいて、各離散化領域の血管領域に厚肉円筒近似を施して形状モデルを構築する。
流れ力学に基づく簡易的な同定法を用いる場合、血流の平均圧力と平均流量とを近似的に求めるため、流体力学における修正ベルヌーイの式、または、ハーゲン・ポアズイユ流れ(Hagen−Poiseuille Flow)の式を用いる。この場合、複数の離散化領域の各々の血圧差と流量との関係を近似的に求めるための形状モデルが構築される。ここで、血圧差は、入口の血圧と出口の圧力との圧力差であり、流量は、単位時間あたり入口流量(または流速)と出口流量(または流速)とを意味する。ただし、力学モデル構築部55は、ステップS2において算出された芯線方向に関する伸縮距離とねじれ角とに基づいて、各節点に関する移動ベクトルと当該節点に隣接する節点に関する回転変位とを各時刻で対応させると良い。
次に、連続体力学モデルを用いる構造流体解析について説明する。
血管および血液(血管と血液を構成する物質)についての運動学は、その運動を引き起こしている力とは無関係である。血管および血液についての運動学の基本概念は、位置、時間、物体、運動、および変形し得る物質の集合についての直感的概念を数学的用語へ抽象化したものである。血管および血液に関する変形および運動の局所的解析を支配する基礎的運動学テンソルは、変形勾配テンソルFと速度勾配テンソルLとである。
変形勾配テンソルFは、運動する血管および血液の物質要素に生じる大きさおよび形の変化を規定する。変形勾配テンソルFは、回転テンソル(正格直交テンソル)Rとストレッチングテンソル(正値対称テンソル)U,Vとの積で表わされる。ストレッチングテンソルU,Vは、まず、基本形態での正規直交ベクトルRによって定められる方向に関するストレッチを課し、次に、正規直交ベクトルRによって与えられる剛体回転を課すことによりもたらされる。なお、ストレッチを課す順番と剛体回転を課す順番とは逆でも良い。
速度勾配テンソルLは、基準形態に依存せず現在形態にのみ依存する。速度勾配テンソルLは、運動する血管および血液の物質要素に生じる大きさおよび形の変化が生じる速度を規定する。速度勾配テンソルLは、ひずみ速度テンソルD(対称テンソル)とスピンテンソル(反対称テンソル)とに分離できる。ひずみ速度テンソルDは、物体がその現在形態をちょうど通過するときのストレッチの変化率を表す。スピンテンソルは、物体がその現在形態をちょうど通過するときの回転の変化率を表す。
血管の力学モデルの外表面節点(または積分点)の一部に関する変形勾配テンソルと速度勾配テンソルとに変位拘束条件(時間的変化を含む)を割り当てる。力学モデルの内腔節点(または積分点)に関する予測値(変形勾配テンソル、速度勾配テンソル、またはそれらの関数値(例えば、変位や面積でも良い))と、観測値(観測データから得られた変形勾配テンソル、速度勾配テンソル、またはそれらの関数値)とが整合するように、材料モデルのパラメータ、負荷条件(力学モデルの内腔における表面力ベクトル)、および境界条件(血管境界における力ベクトル)を同定する。ここで、血管内部の応力の初期状態は、予め仮定しておいてもよいし、同定してもよい。
連続体力学に基づく力学モデルは、運動中の血管および血液における質量、運動学、運動量、角運動量、およびエネルギーの平衡を表す方程式を基礎としている。質量、力、熱、および内部エネルギーという概念が基本である。平衡則とは、運動量の時間全微分が物体力と接触力との和に等しくなり、また、角運動量の時間全微分が物体トルクと接触トルクとの和に等しくなり、また、運動エネルギーと内部エネルギーとの時間的変化が、仕事率(力学的エネルギー)と単位時間あたりの熱供給と熱流束との和に等しくなる、ということを意味している。
平衡則、構成式、および跳躍条件から、変形勾配テンソル、速度勾配テンソル、および応力テンソルといった場の方程式を導き、血管および血液の力学モデルを記述できる。また、血管および血液におけるひずみ場は、適合条件を満たす。
ここで、構成式とは、密度、内部エネルギー、速度ベクトル、応力テンソル、熱流束ベクトル、および温度からなる10個のスカラ方程式の組の関係を与えるものである。場の釣り合い式の中の17個のスカラ場、つまり、密度、内部エネルギー、速度ベクトル、応力テンソル、および熱流束ベクトルのうち、場の釣り合いにより8個のスカラ関係が与えられ、温度を加えた残りの未知量の関係を与える。ただし、物体力bと熱源rとは既知としている。これらのスカラ場を与える方程式のパラメータが材料モデルパラメータである。
連続体力学に基づく力学モデルは、有限要素法又は境界要素法による数値解析法により変位ベクトル、応力テンソル、ひずみテンソル、および速度ベクトル等の場の近似解を、与えられた境界条件、負荷条件、および材料モデルのもとに算出することができる。
構造−流体連成解析において、構造および流体の方程式を解く方法は、一体型解法(monolithic method)と分離型解法(partitioned method)とのいずれでも良い。また、構造と流体との間の境界面での連成は、弱連成でも強連成でもよい。また、流体解析においては、ALE法に代表されるような境界面追跡型の手法で移動境界を扱ってもよいし、「Immersed Boundary method」、「Immersed Finite Element Method」、または「Fictitious Domain Method」等の境界面補足型の手法でもよい。
次に、簡易的な力学モデルの例として、内圧と外圧とを受ける厚肉円筒の材料力学の式と、ハーゲン・ポアズイユ流れおよび修正ベルヌーイの式とについて詳細に説明する。
まず、図24と図25とを参照しながら、厚肉円筒の材料力学の式について説明する。図24は、肉厚円筒の力学モデルの直交断面を示す図である。図25は、図24の微小扇形要素の拡大図である。内半径ra、外半径rbの厚肉円筒に内圧paと外圧pbとが作用する場合の応力やひずみ、変位等の式を説明する。Eおよびνは材料モデルパラメータである。Eは弾性率、νはポアソン比を表している。厚肉円筒では,半径応力σrも考慮し、円周応力σθ の半径方向分布も考慮する必要がある。以降では、軸方向のひずみεzは、断面の位置および向きに関して一様とする。
また円筒断面は、軸対称である場合について説明するが任意形状であっても良い。円筒断面が軸対称である場合、平衡条件は、任意の断面の半径方向についてのみ考えれば良い。任意の断面上で半径rとr+drの同心円筒と中心角dθで切り取られた単位厚さ1の微小扇形要素について、半径方向に関する力の平衡を考える。変形も軸対称であるため、ab面およびbc面には、せん断応力が生じないので、垂直応力のみ作用する。従って、半径方向に関する力の平衡は、以下の(5)式のように表現することができる。
σr rdθ+2σθr sin(dθ/2)-(σr+(dσr/dr)dr)(r+dr)dθ=0 …(5)
ここで、drはrより小さく、dσr はσrより小さいので、(5)式に含まれる高次の微小項を省略し、sin(dθ/2)≒dθ/2とすれば、(5)式は、以下の(6)式のように表現することができる。
rdσr/dr+σr-σθ=0 …(6)
半径rにおける半径方向の変位をuとすれば、u+drでの同方向の変位はu+(du/dr)drとなるので、半径rにおける半径方向のひずみεrはεr=du/drとなる。また、半径方向の変位uによって、半径rの円は半径r+uの円になる。従って、円周ひずみεθは、以下の(7)式のように表現することができる。
εθ =(2π(r+u)-2πr)/2πr=u/r …(7)
また、応力とひずみの関係式から、以下の(8)式または(9)式が得られる。
d2u/dr2+(1/r)(du/dr)-u/r2=0 …(8)
d2u/dr2+d(u/r)/dr=0 …(9)
(8)式または(9)式を積分すると以下の(10)式を得ることができる。
u=c1r+c2/r …(10)
これにより以下の(11)式、(12)式、(13)式が得られる。
σr=(E/((1+ν)(1-2ν)))(c1-(1-2ν)(c2/r2)+νεz) …(11)
σθ=(E/((1+ν)(1-2ν)))(c1+(1-2ν)(c2/r2)+νεz) …(12)
σz=(Eν/((1+ν)(1-2ν)))(2c1+((1-ν)/ν)εz) …(13)
(11)式、(12)式、(13)式における定数c1,c2は、周辺条件、すなわち、円筒の内周r=raでσr=−pa、外周r=rbでσr=−pbから定めることができる。この周辺条件により、(11)式、(12)式、および(13)式から以下の(14)式、(15)式、(16)式をそれぞれ得ることができる。また、変位uは、以下の(17)式のように表現することができる。
σr=(1/(rb2-ra2))(ra2(1-rb2/r2)pa-rb2(1-ra2/r2)pb) …(14)
σθ=(1/(rb2-ra2))(ra2(1+rb2/r2)pa-rb2(1+ra2/r2)pb) …(15)
σz=2ν(ra2pa-rb2pb)/(rb2-ra2)+Eεz=ν(σr+σθ)+Eεz …(16)
u=((1+ν)(1-2ν)/E)((ra2pa-rb2pb)/(rb2-ra2))r+((1+ν)/E)((ra2rb2)/((rb2-ra2)r))(pa-pb)-νεzr …(17)
(16)式および(17)式は、εzの項が含んでいる。従って、(16)式は、対象とする円筒の境界における拘束条件に応じて異なる。例えば、対象とする円筒の両端が拘束される場合、εz =0となり、両端が開放され、σz=0となる。
rθ面内のせん断応力τrは、以下の(18)式のように表現することができる。
τr=(1/2)(σr-σθ)=((ra2rb2)/((rb2-ra2)r2))pb …(18)
θz面内のせん断応力τr’は、σz=0のとき最大になるのでσz=0とすると、以下の(19)式のように表現することができる。
τr’=(1/2)|σθ|=(1/2)((rb2)/(rb2-ra2))(1+ra2/r2)pb …(19)
次に、ハーゲン・ポアズイユ流れおよび修正ベルヌーイの式について説明する。流体が円管内へ流入すると、下流に進むにつれて圧力は降下するとともに、流れの速度分布も徐々に変化する。流れが管内へ流入すると管壁から境界層が発達し、下流へ進むにつれて境界層は厚さを増し、ついには管内の流れは境界層に覆われる。このため、速度分布は管入口のほぼ平らな分布から下流の放物形分布へと変化し、それ以降の速度分布は変化しなくなる。この状態を完全に発達した流れと呼び、管摩擦損失による圧力降下も一定の割合になる。流れが管入口から発達した流れに達する区間を助走区間または入口区間といい、その区間の長さを助走距離または入口長さと呼ぶ。完全に発達した流れの速度分布は下流方向へ変化しないから、以下の(20)式に示すように、管摩擦損失によって生じる圧力損失ΔPの作用力と流体の粘性によって生じるせん断応力τの摩擦力は釣り合うことになる。円管内流れの場合、レイノルズ数Reがおよそ2300以下のとき層流となるといわれている。
速度分布uは層流の場合、以下の(21)式のように軸対称な回転放物面で表わされる。
流量Qは、速度分布uを管断面全体にわたって積分することにより、以下の(22)式のように表現することができる。
断面平均流速vは、以下の(23)式のように表現することができる。
圧力勾配は、流れ方向へ一定で、圧力は減少する。圧力勾配は、管長lの圧力降下Δpを用いると、以下の(24)式のように表現できる。(24)式に(22)式を代入すると(25)式が得られる。(25)式は、流量Qが圧力損失Δpに比例することを意味する。この関係を満たす流れをハーゲン・ポアズイユ流れという。
損失がある場合の修正ベルヌーイの式は、以下の(26)式、(27)式のように表現することができる。
ここで、lは管の長さ、dは管内径、vは管内平均速度、λは管摩擦係数である。管摩擦係数λは、流れが層流の場合にはレイノルズ数Reによって、乱流の場合にはレイノルズ数Reと表面粗さとによって定まる値になる。流体の粘性によって摩擦抵抗が必ず作用する。この摩擦抵抗は、流れを駆動する動力またはエネルギーを消費することになるのでエネルギー損失になる。
ハーゲンポワズイユ流れの式をダルシー・ワイズバッハの式のように変形すると、圧力差と流速および管内径とに関する以下の(28)式が得られる。
以上でハーゲン・ポアズイユ流れおよび修正ベルヌーイの式の説明を終了する。
上述のように、画像処理装置27は、この材料力学の式とハーゲン・ポアズイユ流れおよび修正ベルヌーイの式とを利用して潜在変数を同定することができる。例えば、力学モデルとして、血管の変形を厚肉管の材料力学の式を用い、管径変化を内圧変化と弾性率とにより表現する場合について考える。
無応力状態を初期形状(例えば、血管が最も収縮する状態)と仮定した場合、血管壁およびプラークの弾性率をある値に設定すると、血管内腔の平均半径等の血管形状変形指標の観測値の時間的変化量と内圧の変化量との関係式が得られる。血管形状変形指標の観測値は、時系列のCT画像から計測される。この血管形状変形指標の観測値の時間的変化量に合致するように血管の内圧分布の時間的変化が決定される。この内圧分布の下に血液の流体解析を行うことで血管流量指標の予測値が計測される。この血管流量指標の予測値が観測値に一致しない場合、画像処理装置27は、最初に決めた血管壁またはプラークの弾性率を変更して、さらに同様の解析を行う。
これを繰り返すことにより、画像処理装置27は、血管形状変形指標の観測値と血液流量指標の観測値とに整合する血管壁およびプラークの弾性率、内圧分布、流体解析の圧力境界条件等の潜在変数を決定することができる。この決定方法をより効率的かつ安定的に行うために、階層ベイズモデルとマルコフ連鎖モンテカルロ法とによる統計的同定手法を用いてもよい。
上記のように、血管解析装置50は、記憶装置33、領域設定部51、画像解析・追尾処理部53、力学モデル構築部55、および統計的同定部61を有している。記憶装置33は、被検体の血管に関する時系列の医用画像のデータを記憶する。領域設定部51は、時系列の医用画像に含まれる血管領域に解析対象領域を設定する。画像解析・追尾処理部53は、時系列の医用画像を画像処理して解析対象領域の時系列の形態指標と時系列の形状変形指標とを算出する。
力学モデル構築部55は、時系列の形態指標と時系列の形状変形指標と時系列の医用画像とに基づいて、解析対象領域の構造流体解析に関する力学モデルを暫定的に構築する。統計的同定部61は、暫定的に構築された力学モデルに基づく血管形態指標および血液流量指標の予測値が予め計測された血管形態指標および血液流量指標の観測値に整合するように解析対象領域に関する力学モデルの潜在変数を同定する。
上記構成により、実施の形態の血管解析装置50は、材料モデル、境界条件、負荷条件、および幾何学的構造等の潜在変数を血管形状変形指標と血液流量指標とを用いた逆解析により同定することができる。
血管解析装置50は、潜在変数を変更しながら逆解析を反復して行うことにより、上述の4点の困難、すなわち、1.冠動脈の材料モデルの同定方法、2.心臓の形状の変形の冠動脈への影響の組み込み、3.冠動脈の境界条件の同定方法、4.不確定性を有する血管形状を利用した材料モデルや負荷条件、境界条件の同定方法、を全て加味した潜在変数を同定することができる。
従って、血管解析装置50は、CT画像に描出されない血管や心臓等の外部要因による影響を加味した構造流体解析を実行することができる。このように、実施の形態の医用画像診断装置によれば、血管の構造流体解析の精度の向上を図ることができる。
次に、このような血管の力学モデルによる構造解析、流体解析、または構造流体連成解析結果から得られた様々な指標(血管の断面積、応力の変化、血流の流速・流量・圧力の変化等)を用いてノミナルモデルを作成し、診断対象のデータとのモデル同定により血管の異常を診断する方法について説明する。なお、以下、血管の異常診断として、血管の狭窄度の診断を行う例を説明する。
ある2組の指標を選択し、一方を入力、他方を出力としたとき、入出力関係に因果関係があるならば入出力モデルを仮定し、モデルを同定することが可能となる。すなわち、図26に示すように、狭窄の有無に依存しない指標を入力u1とし、狭窄の有無で影響を受ける指標を出力とする。また、狭窄なしの出力をy1、狭窄ありの出力をy2とすると、狭窄なしの場合は入力u1、出力y1としたモデルG1が仮定される。同様に、狭窄ありの場合は、入力u1、出力y2としたモデルG2が仮定される。
このとき、モデルG1およびモデルG2が有限個のパラメータで表現され、モデル同士を比較できるならば、血管の狭窄度を評価できる。モデルの同定精度が高ければ、推定したモデルG1hatに、入力u1を入力した時に出力される出力y1hatは、観測したy1との誤差が小さい。なお、「hat」は、推定値を意味する。
一般的に推定精度を上げるためには、モデルの次数を上げる必要がある。高次のモデル同士の比較は、評価対象となるパラメータが多くなるため、結果が不安定なものになりやすい。そこで、予め同定したモデルを「ノミナルモデル」とし、その出力y1hatを計算しておく。y1hatを入力とし、比較対象となる指標y2を出力とする新たな入出力モデルG12を仮定し、同定する。このモデルは狭窄の有無を表現したものであるため、以下、「狭窄モデル」と呼ぶこととする。y1hatには、高次のモデルの要素が含まれている。このため、狭窄モデルとして、次数を減らした低次モデルを用いることができる。
一例を図27に示す。図27のグラフは、横軸が時間、縦軸がCT値(造影剤の濃度)を示している。この図27のグラフにおいて、○で示すグラフは、狭窄ありの出力y2のグラフである。実線のグラフは、ノミナルモデルG2を同定し、共通の入力u1を入力したときの出力y2hatのグラフである。これに対し、点線のグラフは、狭窄なしの出力y1から同定した別のノミナルモデルG1で推定した出力y1hatを入力とし、y2を出力とした狭窄モデルG12を同定し、y1hatを入力した時の推定値のグラフである。ノミナルモデルG1は3/6次、狭窄モデルG12は2/2次であり、少ないパラメータでありながら推定精度が高いのが分かる。
図28に、実施の形態の医用画像診断装置の、血管の狭窄度の診断を行う場合における、血管解析装置50の機能ブロック図を示す。この場合、血管解析装置50は、構造流体解析を行う上述の画像処理装置27、およびノミナルモデル(後述するARXモデルの各種パラメータ)が記憶される記憶装置33と共に、造影剤注入プロファイル取得部71を有する。また、血管解析装置50は、第1入出力モデル算出部72、時系列CT画像取得部73、血管芯線抽出部74、時系列CT値抽出部75、第2入出力モデル算出部76、狭窄モデル算出部77、および狭窄度診断部78を有している。
造影剤注入プロファイル取得部71は、注入プロファイル取得部の一例である。時系列CT画像取得部73は、時系列CT値取得部の一例である。画像処理装置27は、構造流体解析部の一例である。第1入出力モデル算出部72は、第1算出部の一例である。造影剤注入プロファイル取得部71、第1入出力モデル算出部72、時系列CT画像取得部73、画像処理装置27、および記憶装置33は、時系列モデル生成部の一例である。記憶装置33は、記憶装置の一例である。時系列CT値抽出部75は、時系列CT値抽出部の一例である。狭窄モデル算出部77は、狭窄モデル算出部の一例である。狭窄度診断部78は、診断部の一例である。第2入出力モデル算出部76は、第2算出部の一例である。
造影剤注入プロファイル取得部71〜狭窄度診断部78は、全てをハードウェアで構成してもよいし、全てをソフトウェアで構成してもよい。または、一部をハードウェアで構成し、残りをソフトウェアで構成してもよい。
まず、血管の狭窄度の診断を行う場合、上述のノミナルモデルを事前に用意する。順を追って説明すると、CT画像の撮像には、血管造影のために造影剤を用いることが多い。造影剤が、例えば上腕部の静脈より一定量注入されると、造影剤が通過した領域のCT値(HU値)が上昇し、血管以外の組織との識別が容易になる。CT架台10は、時系列でCT画像を撮影するように、システム制御部21および架台制御部23により制御される。これにより、CT値が変化する速度を観測することができる。
一例として、図29に造影剤の注入パターンを示す。図29の実線のグラフ、点線のグラフ、および一点鎖線のグラフは、それぞれ造影剤の注入パターンを示している。図29の実線のグラフは、造影剤を約7秒間かけてゆっくり注入した注入パターンを示している。図29の点線のグラフは、造影剤を約4.5秒間かけて高速に注入した注入パターンを示している。図29の一点鎖線のグラフは、造影剤を約6秒間かけて中速で注入した注入パターンを示している。なお、造影剤の注入は、担当医師等が手動で行ってもよいし、担当医師監督のもと、注入機を用いて行ってもよい。
図30は、造影剤注入後の注目領域におけるCT値の変化のグラフである。図29の実線のグラフ、点線のグラフ、および一点鎖線のグラフに示した注入パターンのうち、どの注入パターンで造影剤を注入しても、大体、図30と同様のグラフを得られる。この図30において、○のグラフは、狭窄なしの場合のCT値の変化を示している。また、図30において、●のグラフは、狭窄ありの場合のCT値の変化を示している。この図30は、注目領域に狭窄が無い場合、CT値の変化は、○のグラフに示すように立ち上がりが早く、また、ピークも高くなる。これに対して、注目領域の上流側に狭窄部が存在する場合、CT値の変化は、●のグラフに示すように立ち上がりが遅く、また、ピークも低くなる。実施の形態の医用画像診断装置では、造影剤の注入パターンを入力、注目領域でのCT値を出力とした入出力モデルを用いている。
図28の説明に戻り、ノミナルモデルを生成する場合、造影剤注入プロファイル取得部71は、図29を用いて説明した造影剤の注入プロファイルを取得する。取得された造影剤の注入プロファイルは、画像処理装置27、第1入出力モデル算出部72、および第2入出力モデル算出部76に供給される。
上述の画像処理装置27は、時系列CT画像取得部73で取得された時系列のCT画像(MRI画像等でもよい)を、構造流体解析する。なお、時系列CT画像取得部73は、記憶装置33に記憶されている時系列のCT画像を取得してもよいし、CT架台10でリアルタイムかつ時系列で撮像されたCT画像を取得してもよい。画像処理装置27の解析により、流体の一つである造影剤の注入プロファイルに対応する出力(造影剤の濃度変化)をシミュレーションすることが可能となる。構造流体解析部は、図29を用いて説明した注入パターンで注入された造影剤の濃度変化をシミュレーションする。そして、このシミュレーション結果となる図30に示す造影剤の濃度変化を示す情報を、第1入出力モデル算出部72に供給する。
第1入出力モデル算出部72(および第2入出力モデル算出部76)は、有限個のパラメータからなる線形差分方程式で表現されるモデルを算出する。すなわち、第1入出力モデル算出部72(および第2入出力モデル算出部76)は、有限個のパラメータで表現する入出力モデルの一例として、いわゆるARXモデルを算出する。ARXは、「Auto Regressive with eXogenous inputs)、または(Auto Regressive with eXternal inputs)の略記である。ARXモデルは、入力をu(t)、出力をy(t)、式誤差ε(t)とすると、以下の(29)式で表される。
A(q)y(t) = B(q)u(t) + ε(t) …(29)
ただし、A(q)は、以下の(30)式、B(q)は、以下の(31)式のとおりである。
A(q) = 1 + a1q-1 + … + amq-m …(30)
B(q) = b1q-1 + … + bnq-n …(31)
(29)式における「t」は、各時間でサンプリングされた離散的なサンプリングデータ(図30の○または●のCT値(サンプリングデータ)を参照)を意味している。また、(30)式および(31)式における「q」が、時間軸に沿ってシフトするシフトパラメータを意味している。例えば、「q-1」は、時間的に一つ前の値を用いることを意味している。また、この例のモデルは、合計m+n個のパラメータで表現されている。このモデルのパラメータA、Bは、式誤差の二乗和を最小化するように最小二乗法を用いて決定される。
このように算出されたARXモデルの各パラメータは、ノミナルモデルとして記憶装置33に記憶される。一例として、記憶装置33に記憶されるノミナルモデルとしては、狭窄なしの場合のノミナルモデル、狭窄度30パーセントの場合のノミナルモデル、狭窄度50パーセントの場合のノミナルモデル等の複数のノミナルモデルが記憶される。なお、この例では、入出力モデルとしてARXモデルを用いたが、パラメータで表現可能なモデルであれば、どのようなモデルを用いてもよい。
次に、このようなノミナルモデルが記憶装置33に記憶されると、血管の狭窄度診断が可能となる。血管の狭窄度診断時となると、造影剤が注入された被検体の時系列CT画像を、時系列CT画像取得部73が取得する。また血管芯線抽出部74は、取得された時系列CT画像から、上述のように血管の芯線を抽出する。時系列CT値抽出部75は、芯線が抽出された血管の注目領域のCT値を時系列で抽出する。この時系列で抽出された各CT値が、図27のグラフの○印で示すCT値に相当する。
次に、第2入出力モデル算出部76は、第1入出力モデル算出部72と同様に、有限個のパラメータで表現する入出力モデルの一例として、いわゆるARXモデルを算出する。第2入出力モデル算出部76は、各CT値を取得した際に被検体に注入した造影剤の注入プロファイルを、造影剤注入プロファイル取得部71を介して取得する。そして、第2入出力モデル算出部76は、造影剤の注入プロファイルおよび時系列で抽出された各CT値から、記憶装置33に記憶されているノミナルモデルと同様のパラメータセットとなるARXモデルを、第2入出力モデルとして算出する。
最後に、狭窄モデル算出部77が、記憶装置33に記憶されているノミナルモデル(図27の実線のグラフに相当するパラメータ)と第2入出力モデルとの各パラメータの差分となる狭窄モデルを算出し、狭窄度診断部78は、算出された狭窄モデルのパラメータを、血管の狭窄度として診断する。
すなわち、造影剤注入プロファイルをノミナルモデルに入力した時の出力結果を入力とし、同じく造影剤注入プロファイルを第2入出力モデルに入力したときの出力結果を出力としてARXモデルで計算を行う。仮に、第2入出力モデルとノミナルモデルとの間に差が無ければ、得られた狭窄モデルの伝達関数は「1」となり、第2入出力モデルとノミナルモデルとが同一であることを示すようになる。
狭窄度診断部78は、狭窄モデルをそのまま評価(伝達関数1からどれくらい離れているか)してもよいし、例えば30パーセントの狭窄、10パーセントの狭窄等のように、血管の狭窄度を表示してもよい。狭窄度の異なるデータを使って作成した複数のノミナルモデルで評価を行い、「どのノミナルモデルに近いか」ということを表示してもよい。さらには、狭窄度診断部78は、狭窄が無しの正常時の値を「0」として、スケールで狭窄の度合いを表示してもよい。すなわち、狭窄の度合いが分かれば、どのような表示形態で表示してもよい。
以上の説明から明らかなように、実施の形態の医用画像診断装置は、画像追尾、構造流体解析および逆解析(統計的同定法)に基づく血管解析により、血行動態における血圧、血流、応力、またはひずみ等の狭窄等の血管病変の予防、診断に必要な指標を高精度に解析する。そして、実施の形態の医用画像診断装置は、狭窄の有無が判明している典型例を対象に解析を行い、低次の標準モデル(ノミナルモデル)を予め作成しておき、診断対象データとのモデル同定から狭窄度を診断する。
これにより、画像処理装置で計算された解析結果を有効に利用して、血管の異常診断を可能とすることができる。また、低次の標準モデル(ノミナルモデル)と診断対象データとのモデル同定をから狭窄度を診断することで、少ない計算リソースおよび極めて短時間で血管の異常診断を行うことができる。
また、時系列CT値抽出部75で抽出された各CT値、および造影剤注入プロファイルから第2入出力モデル算出部76でARXモデルを算出し、算出したARXモデルと、記憶装置33に記憶されているノミナルモデルとを、狭窄モデル算出部77で比較している。このため、同じARXモデル同士を比較することができるため、精度よく、狭窄度を診断することができる。
(変形例1)
上述の実施の形態の説明では、同じARXモデル同士を比較するために、第2入出力モデル算出部76でARXモデルを算出したが、時系列CT値抽出部75で抽出された各CT値と、記憶装置33に記憶されているノミナルモデルとを狭窄モデル算出部77で比較して、狭窄度を診断してもよい。この場合、第2入出力モデル算出部76を省略することができるため、実施の形態の医用画像診断装置の構成の簡略化および低コスト化を図ることができる。
(変形例2)
時系列のCT値として、血管の狭窄前後を含む芯線方向にCT値を取得し、血管の走行方向における勾配(造影剤濃度勾配)を算出した勾配の時系列データを用いてもよい。図31は、造影剤濃度勾配のグラフの一例である。図31中、実線で示すグラフは、狭窄なしの場合の造影剤濃度勾配のグラフである。造影剤濃度勾配は、負の勾配となり、狭窄なしの場合、早い時間で立下り、早い時間で0レベルに復帰する勾配のグラフとなる。これに対して、図31中、点線で示すグラフは、狭窄ありの場合の造影剤濃度勾配のグラフである。狭窄ありの場合、緩やかに時間をかけて立下り、また、緩やかに時間をかけて0レベルに復帰する勾配のグラフとなる。このような勾配の時系列データを用いても、上述と同様の効果を得ることができる。
以上、実施の形態および変形例を説明したが、実施の形態および変形例は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施の形態および変形例は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。実施の形態およびその変形は、発明の範囲や要旨に含まれると共に、特許請求の範囲に記載された発明とその均等の範囲に含まれる。