以下に添付図面を参照して、処理装置、血管模型、画像処理方法、プログラム、および、造形装置の実施の形態を詳細に説明する。
なお、実施の形態では、一例として、処理装置を、医用画像診断装置で生成された画像を構造流体解析するためのコンピュータ装置に適用した場合を説明する。なお、本実施の形態では、一例として、医用画像診断装置は、医用画像を解析する場合を説明する。このコンピュータ装置は、医用画像診断装置に組み込まれていても良いし、医用画像診断装置とは別体のワークステーション等であっても良い。以下、実施の形態の処理装置の適用されたコンピュータ装置を組み込んだ医用画像診断装置を、図面を参照して詳細に説明する。なお、処理装置は、医用画像診断装置とは別体として構成されていてもよい。
実施の形態の医用画像診断装置は、被検体をスキャンするための撮像機構を装備する如何なる種類の画像診断装置にも適用可能である。実施の形態の医用画像診断装置としては、例えばX線コンピュータ断層撮影装置(X線CT装置)、磁気共鳴診断装置、超音波診断装置、SPECT(Single Photon Emission CT)装置、PET(Positron Emission Tomography)装置、および放射線治療装置等に適宜利用可能である。以下、説明を具体的に行うため実施の形態の医用画像診断装置は、X線コンピュータ断層撮影装置であるものとする。
図1は、実施の形態の医用画像診断装置(X線コンピュータ断層撮影装置)の概略的なハードウェア構成図である。図1に示すように、X線コンピュータ断層撮影装置は、CT架台10と、コンソール20と、造形部72と、を有する。CT架台10と、造形部72と、はコンソール20に信号授受可能に接続されている。
造形部72は、三次元造形物を製造する公知の装置である。造形部72は、三次元造形物を造形可能な装置であればよく、例えば、熱溶解積層方式、粉末固着方式の何れであってもよい。
なお、本実施の形態では、造形部72が造形に用いる材料は、人間の血管の機械的特性の取り得る範囲を満たす材料を用いることが好ましい。具体的には、造形部72が造形に用いる材料は、ヤング率0.01MPaから20MPa程度の材料であることが好ましい。例えば、造形部72は、シリコーンゴムなどの、可撓性、弾発性を有する模擬血管材料を用いて造形することが好ましい。
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は、架台制御部23、再構成装置25、血管解析装置50、処理装置70、入力部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に記憶される。
血管解析装置50は、システム制御部21と、解析装置27と、を含む。
システム制御部21は、架台制御部23、再構成装置25、解析装置27、入力部29、表示部31、および記憶部33に接続されており、これらを制御する。
解析装置27は、画像を用いて構造流体解析を実行する。本実施の形態では、一例として、解析装置27は、時系列のCT画像(医用画像)を用いて、構造流体解析を実行する場合を説明する。そして、解析装置27は、被検体の血管の所定の位置の材料特性(詳細後述)を算出する。具体的には、解析装置27は、被検体の医用画像から後述する処理によって最終的に構築した力学モデルを解析することで、材料特性を算出する。本実施の形態では、一例として、解析装置27は、被検体の血管の各位置の材料特性を算出する場合を説明する。解析装置27の詳細は後述する。
医用画像は、被検体の血管の形状を同定可能な画像である。医用画像は、例えば、時系列のCT画像や、MRI画像や、超音波エコー画像などである。
本実施の形態では、医用画像として、時系列のCT画像を用いる場合を一例として説明する。
血管解析装置50は、医用画像診断装置(X線コンピュータ断層撮影装置)に組み込まれていても良いし、医用画像診断装置とは別体のコンピュータ装置であっても良い。血管解析装置50が医用画像診断装置とは別体の場合、血管解析装置50は、医用画像診断装置やPACS(picturearchiving and communication systems)からネットワークを介して時系列のCT画像等の医用画像を収集すれば良い。
記憶部33は、ハードディスク装置等の種々の記憶媒体により構成される。記憶部33は、時系列の投影データや時系列のCT画像などの種々のデータを記憶する。例えば、記憶部33は、時系列のCT画像をDICOM(digital imaging and communications in medicine)規格に準拠した医用画像ファイル形式で記憶する。また、記憶部33は、外部機器により収集された医用画像を記憶しても良い。
入力部29は、ユーザからの各種指示や情報入力を受け付ける。入力部29は、例えば、キーボード、マウス、スイッチ等である。
表示部31は、各種画像を表示する。表示部31は、例えばCRTディスプレイや、液晶ディスプレイ、有機ELディスプレイ、プラズマディスプレイ等である。
なお、入力部29と表示部31とを一体的に構築したUI部73としてもよい。UI部73には、表示機能と入力機能の双方を備えたタッチパネルなどがある。
処理装置70は、解析装置27、再構成装置25、入力部29、表示部31、及び記憶部33に電気的に接続されている。また、処理装置70は、造形部72に電気的に接続されている。
処理装置70は、血管モデルを生成する装置である。
血管モデルは、被検体の血管の血管模型の生成に用いる画像データであると共に、被検体の血管の血管画像の生成に用いる画像データである。
図2は、処理装置70の機能ブロック図である。
処理装置70は、第1取得部70Aと、第2取得部70Bと、構築部70Cと、第1生成部70Dと、第1制御部70Eと、受付部70Fと、第2生成部70Gと、第2制御部70Hと、第3制御部70Iと、を含む。
第1取得部70A、第2取得部70B、構築部70C、第1生成部70D、第1制御部70E、受付部70F、第2生成部70G、第2制御部70H、及び第3制御部70Iの一部またはすべては、例えば、CPU(Central Processing Unit)などの処理装置にプログラムを実行させること、−すなわち、ソフトウェアにより実現してもよいし、IC(Integrated Circuit)などのハードウェアにより実現してもよいし、ソフトウェアおよびハードウェアを併用して実現してもよい。
第1取得部70Aは、被検体の血管に関する画像を取得する。第1取得部70Aは、再構成装置25(図1参照)から、被検体の血管に関する画像を取得する。なお、第1取得部70Aは、記憶部33から、被検体の血管に関する画像を取得してもよい。
第2取得部70Bは、上記医用画像の被検体の血管の所定の位置の材料特性を取得する。第2取得部70Bは、解析装置27から、被検体の血管の所定の位置の材料特性を取得する。本実施の形態では、第2取得部70Bは、解析装置27から、被検体の血管の各位置の材料特性を取得する場合を説明する。
材料特性は、被検体の血管の、圧力に対する血管変形率を示す。材料特性は、具体的には、被検体の血管の、芯線方向及び周方向に沿った各位置の、弾性率、可撓性、強度、硬度(硬さ)、靱性などを示す。
弾性率は、例えば、ヤング率やポアソン比である。強度は、例えば、剛性、引張強さ、圧縮強さ、せん断強さ、などを示す。硬さは、例えば、ビッカース硬さ、ブリネル硬さ、ロックウェル硬さ、ショア硬さ、ヌープ硬さ、モース硬さ、などを示す。靭性は、物質破壊に対する仕事量、粘り強さ、脆性(もろさ)などを示す。
なお、第2取得部70Bは、解析装置27から、力学モデルを取得してもよい。そして、力学モデルを解析することで、被検体の血管の材料特性を取得してもよい。
力学モデルは、時系列の医用画像から構築された、血管の三次元形状を示す形状モデルに、血管の形状履歴(時系列の血管形態指標)や、強制変位履歴(時系列の血管形状変形指標)を割当てることで、医用画像からは直接的には導き出すことの困難な、血管の変形挙動などを付加したデータである。
構築部70Cは、第1取得部70Aで取得した被検体の時系列の医用画像から、形状モデルを構築する。形状モデルは、被検体の血管の三次元形状を示すデータである。形状モデルは、具体的には、各時刻における血管の幾何学的構造を示すデータである。形状モデルの構築には、公知の方法を用いればよい。
図3は、形状モデル80の一例を示す図である。図3(A)に示すように、形状モデル80は、被検体の血管の三次元の形状そのものを示すデータであり、医用画像から構築される。このため、形状モデル80は、被検体の実際の血管の弾性率などの材料特性を含むものではない。
実際には、図3(B)に示すように、形状モデル80によって示される形状の血管に対応する、実際の被検体の血管の材料特性は、血管における各位置によって異なる。例えば、図3(B)に示すように、形状モデル80によって示される形状の血管の血管壁82の内、例えば、プラークの形成された位置A1の剛性は高く、プラークの形成されていない位置A3、プラークに連続する位置A2、の順に、剛性が低い。しかし、医用画像から構築された形状モデル80は、形状そのものを示すため、形状モデル80から実際の被検体の血管の材料特性を得ることは出来ない。
そこで、第1生成部70Dは、構築部70Cが構築した形状モデルによって示される、被検体の血管の所定の位置に、第2取得部70Bで取得した被検体の血管の材料特性に関する付加情報を付加した血管モデルを生成する。本実施の形態では、一例として、第1生成部70Dは、該形状モデルによって示される被検体の血管の各位置に、各位置に対応する材料特性に関する付加情報を付加した血管モデルを生成する場合を説明する。
血管モデルは、被検体の血管の三次元形状を示すデータである形状モデルに、該血管の各位置の材料特性に関する付加情報を付加したものである。
付加情報は、形状モデルによって示される血管の位置ごとの、材料特性に関する情報であればよい。付加情報は、材料特性そのものであってもよいし、材料特性を実現可能な手段や方法を示す情報であってもよいし、材料特性を異なる種類の特性に変換した情報であってもよい。
また、付加情報は、形状モデルを構成する各画素の三次元座標ごとに、被検体の血管における対応する位置の材料特性に関する情報を示したものであればよい。
材料特性に関する付加情報は、例えば、血管壁の内壁に対する外壁の相対位置を補正した情報、血管モデルに応じた血管模型の各位置の造形材料の種類を示す情報、及び、補強部材の付与位置を示す情報、の少なくとも1つである。これらの付加情報はいずれも、材料特性に応じて、血管の位置毎に生成される。
まず、付加情報が、血管壁の内壁に対する外壁の相対位置を補正した情報である場合を説明する。例えば、第1生成部70Dは、生成する血管モデルによって示される血管壁の厚みが、被検体の血管の各位置の材料特性に応じた厚みとなるように、血管モデルによって示される血管の内壁に対する外壁の相対位置を補正した付加情報を付加した、血管モデルを生成する。
図4は、血管モデル生成の説明図である。図4(A)は、形状モデル40の一例を示す図である。図4(B)は、図4(A)に示す形状モデル40に、付加情報48を付加した血管モデルM102の一例の説明図である。
例えば、第2取得部70Bで取得した被検体の血管の材料特性としての硬度が、図4(A)に示すように、血管における、位置A1>位置A3>位置A2の関係を示すとする。すなわち、位置A1の硬度が最も高く、位置A3、位置A3のこの順に硬度が低いとする。
この場合、第1生成部70Dは、図4(B)に示すように、生成する血管モデルM102によって示される血管壁42の厚みが、被検体の血管の各位置A1〜位置A3の材料特性を実現可能な厚みとなるように、血管モデルM102によって示される血管の内壁40Aに対する外壁40Bの相対位置を補正した付加情報(図4(B)中、付加情報48A〜48C参照)を付加した、血管モデルM102を生成する。
具体的には、第1生成部70Dは、材料特性が硬度を示す場合、硬度が高くなるほど、血管壁42の厚みが厚くなるように、血管の内壁40Aに対する外壁40Bの相対位置を補正した付加情報を付加する。同様に、第1生成部70Dは、材料特性が弾性率を示す場合、弾性率が低いほど、血管壁42の厚みが厚くなるように、血管の内壁40Aに対する外壁40Bの相対位置を補正した付加情報を付加する。同様に、第1生成部70Dは、材料特性が可撓性を示す場合、可撓性が低くなるほど、血管壁42の厚みが厚くなるように、血管の内壁40Aに対する外壁40Bの相対位置を補正した付加情報を付加する。また、第1生成部70Dは、材料特性が強度を示す場合、強度が高くなるほど、血管壁42の厚みが厚くなるように、血管の内壁40Aに対する外壁40Bの相対位置を補正した付加情報を付加する。
また、第1生成部70Dは、材料特性が、弾性率、可撓性、強度、硬度、および、靱性などの材料特性を示すパラメータの複数を組み合わせた特性を示す場合、これらの複数のパラメータの材料特性を実現可能な厚みとなるように、血管の内壁40Aに対する外壁40Bの相対位置を補正した付加情報を付加する。
すなわち、第1生成部70Dは、形状モデルによって示される被検体の血管の内径は変えず、血管の外径を変更するように(すなわち、外壁の位置を変更するように)、血管の内壁40Aに対する外壁40Bの相対位置を補正した付加情報48(付加情報48A〜48C参照)を付加する。
このため、血管モデルM102は、図4(B)に示すように、血管の内径については、形状モデル40によって示される内径と略同じである。一方、血管モデルM102によって示される血管の外径については、形状モデル40によって示される血管の外径とは異なる。すなわち、血管モデルM102は、対応する位置の材料特性を実現可能な厚みとなるように、内壁40Aに対する外壁40Bの相対位置が補正されたものとなる。
次に、付加情報が、血管モデルに応じた血管模型の各位置の造形材料の種類を示す情報である場合を説明する。
例えば、第1生成部70Dは、血管モデルに応じた血管模型の各位置の造形材料が、対応する各位置の材料特性を満たすように、血管模型の各位置の造形材料の種類を示す付加情報を付加した血管モデルを生成する。
図5は、血管モデル生成の説明図である。図5(A)は、形状モデル40の一例を示す図であり、図4(A)と同様である。図4(A)と同様に、例えば、第2取得部70Bで取得した被検体の血管の材料特性としての硬度が、図5(A)に示すように、位置A1>位置A3>位置A2の関係を示すとする。すなわち、位置A1の硬度が最も高く、位置A3、位置A3のこの順に硬度が低いとする。
図5(B)は、図5(A)に示す形状モデル40に、付加情報46を付加した血管モデルM100の一例の説明図である。
この場合、第1生成部70Dは、図5(B)に示すように、生成する血管モデルM100に応じた血管模型D100によって示される血管壁42の造形材料が、被検体の血管の各位置A1〜位置A3の材料特性を満たすように、血管模型D100の各位置の造形材料の種類を示す付加情報46(付加情報46A〜46C)を付加した、血管モデルM100を生成する。
具体的には、第1生成部70Dは、材料特性が硬度を示す場合、材料特性によって示される各位置の硬度を実現可能な材料であって、且つ、血管模型D100を造形する造形部72(図1参照)で使用可能な造形材料の種類を示す情報を付加する。同様に、第1生成部70Dは、材料特性が弾性率を示す場合、材料特性によって示される各位置の弾性率を実現可能な材料であって、且つ、血管模型D100を造形する造形部72(図1参照)で使用可能な造形材料の種類を示す付加情報を付加する。また、同様に、第1生成部70Dは、材料特性が可撓性を示す場合、材料特性によって示される各位置の可撓性を実現可能な材料であって、且つ、血管模型D100を造形する造形部72(図1参照)で使用可能な造形材料の種類を示す付加情報を付加する。また、第1生成部70Dは、材料特性が強度を示す場合、材料特性によって示される各位置の強度を実現可能な材料であって、且つ、血管模型D100を造形する造形部72(図1参照)で使用可能な造形材料の種類を示す情報を付加する。
また、第1生成部70Dは、材料特性が、弾性率、可撓性、強度、硬度、および、靱性などの材料特性を示すパラメータの複数を組み合わせた特性を示す場合、これらの複数のパラメータの材料特性に応じた造形材料の種類を示す付加情報を付加する。
すなわち、第1生成部70Dは、材料特性に応じた機械的強度を実現するために、形状モデルによって示される被検体の血管の形状は変えず、血管模型D100の造形に用いる造形材料の種類を位置ごとに定めた付加情報を付加する。
このため、血管モデルM100は、図5(B)に示すように、血管の内径、外径、及び形状については、形状モデル40によって示される形状及び大きさと略同じである。しかし、血管モデルM100には、形状モデル40に、対応する位置の材料特性を実現可能な造形材料の種類が付加情報として付与されている。
次に、付加情報が、補強部材の付与位置を示す情報である場合を説明する。
例えば、第1生成部70Dは、血管モデルに応じた血管模型の各位置の機械特性が、対応する各位置の材料特性を満たすように、補強部材の付与位置を示す付加情報を付加した血管モデルを生成する。
図6は、血管モデル生成の説明図である。図6(A)は、形状モデル40の一例を示す図であり、図4(A)と同様である。図4(A)と同様に、例えば、第2取得部70Bで取得した被検体の血管の材料特性としての硬度が、図6(A)に示すように、位置A1>位置A3>位置A2の関係を示すとする。すなわち、位置A1の硬度が最も高く、位置A3、位置A3のこの順に硬度が低いとする。
図6(B)は、図6(A)に示す形状モデル40に、付加情報49を付加した血管モデルM104の一例の説明図である。
この場合、第1生成部70Dは、図6(B)に示すように、生成する血管モデルM104に応じた血管模型D104の各位置の機械強度が、被検体の血管の各位置A1〜位置A3の材料特性を満たすように、血管模型D104によって示される血管を外側から補強する補強部材Pの種類、および補強部材Pの付与位置を示す、付加情報49(付加情報49A〜49C)を付加した、血管モデルM104を生成する。
更に詳細には、第1生成部70Dは、図6(B)に示すように、生成する血管モデルM104に応じた血管模型D104における、血管壁42の内壁側から測定したときの強度が、対応する各位置の材料特性を満たすように、補強部材Pの種類、および、血管壁42の外側への補強部材Pの付与位置、を示す付加情報49(付加情報49A〜49C)を付加した、血管モデルM104を生成する。
具体的には、第1生成部70Dは、材料特性が硬度を示す場合、硬度が高い位置ほど、より厚みの厚い補強部材Pを設置するように、補強部材Pの種類(例えば、厚み)と、血管壁42の外側への補強部材Pの付与位置と、を示す付加情報を付加した、血管モデルM104を生成する。同様に、第1生成部70Dは、材料特性が弾性率を示す場合、弾性率が低い位置ほど、より厚みの厚い補強部材Pを設置するように、補強部材Pの種類と、血管壁42の外側への補強部材Pの付与位置と、を示す付加情報を付加した、血管モデルM104を生成する。同様に、第1生成部70Dは、材料特性が可撓性を示す場合、可撓性が低い位置ほど、より厚みの厚い補強部材Pを設置するように、補強部材Pの種類と、血管壁42の外側への補強部材Pの付与位置と、を示す付加情報を付加した、血管モデルM104を生成する。同様に、第1生成部70Dは、材料特性が強度を示す場合、強度が高い位置ほど、より厚みの厚い補強部材Pを設置するように、補強部材Pの種類と、血管壁42の外側への補強部材Pの付与位置と、を示す付加情報を付加した、血管モデルM104を生成する。
なお、付加情報49は、少なくとも補強部材Pの付与位置を示す情報であればよく、補強部材Pの付与位置と、補強部材Pの種類と、の双方を含む形態に限定されない。また、補強部材Pの種類は、例えば、厚み、材質、形状、などを示す。
また、第1生成部70Dは、材料特性が、弾性率、可撓性、強度、硬度、および、靱性、などの材料特性を示すパラメータの複数を組み合わせた特性を示す場合、これらの複数のパラメータの材料特性に応じた位置に、材料特性に応じた厚みまたは材質の補強部材Pを付与するように、補強部材Pの種類と、補強部材Pの付与位置と、を示す付加情報を付加した、血管モデルM104を生成する。
図6(B)に示す例では、第1生成部70Dは、血管模型D104によって示される血管の、被検体の血管の位置A3に対応する位置には、補強部材P2及び補強部材P3を設置することを示す、付加情報を付加する。なお、補強部材P2と補強部材P3の材質は同じであり、厚みも同じであるとする。また、第1生成部70Dは、血管模型D104によって示される血管の、被検体の血管の位置A1に対応する位置には、補強部材P2より厚みの厚い補強部材P1を設置することを示す付加情報を付加する。そして、第1生成部70Dは、血管模型D104によって示される血管の、被検体の血管の位置A2に対応する位置には、何も補強部材Pを設置しないことを示す付加情報を付加する。
なお、各位置に設置する補強部材Pは、同じ材料で構成され、厚みを調整することによって材料特性を実現可能としてもよい。また、各位置に設置する補強部材Pは、厚みは略同じとし、位置毎に、材料特性を実現可能な異なる材料で構成してもよく、位置ごとに、異なる材料で構成されていてもよい。
図2に戻り、第1制御部70Eは、構築部70Cで構築された形状モデルによって示される血管画像を表示部31に表示する。図7は、血管画像90の一例を示す模式図である。
図2に戻り、受付部70Fは、ユーザによる指示を受け付ける。具体的には、ユーザは、表示部31に表示された血管画像90を参照しながら、入力部29を操作することで、血管画像90における特定の位置を指示する。すると、受付部70Fは、血管画像90における、ユーザによって指示された位置を示す位置情報を、ユーザによる指示として受け付ける。
第2生成部70Gは、表示された血管画像90における、ユーザによって指示された位置に対応する付加情報によって示される材料特性に応じた強度の出力信号を生成する。
例えば、第2生成部70Gは、表示された血管画像90における、ユーザによって指示された位置に対応する付加情報によって示される材料特性を、第1生成部70Dによって生成された血管モデルから読取る。そして、第2生成部70Gは、読取った材料特性に応じた強度の出力信号を生成する。
例えば、第2生成部70Gは、読取った材料特性が弾性率を示す場合、弾性率が低いほど、強度の大きい出力信号(例えば、電圧値の大きい出力信号)を生成する。また、例えば、第2生成部70Gは、読取った材料特性が硬さを示す場合、硬度が高いほど、強度の大きい出力信号(例えば、電圧値の大きい出力信号)を生成する。また、例えば、第2生成部70Gは、読取った機械的特性が靱性を示す場合、靱性が高いほど、強度の大きい出力信号を生成する。
そして、第2制御部70Hは、第2生成部70Gで生成した出力信号の強度に応じた応力を生じさせるように、入力部29を制御する。
本実施の形態では、応力は、入力部29における、ユーザによって操作される予め定められた領域がユーザの押圧による荷重をうけたときに、該領域内に生じる単位面積当たりの内力を示す。例えば、第2制御部70Hは、第2生成部70Gで生成した出力信号の強度が高いほど、応力が大きくなるように、入力部29を制御する。
例えば、図7に示す血管画像90上に表示されたポインタP1の位置を、ユーザによる入力部29の操作指示によって操作することで、血管画像90上のある位置を指示したとする。このときに、入力部29を操作するユーザの手などは、該位置に対応する材料特性が高いほど、強い応力が伝えられる。具体的には、入力部29におけるユーザによって把握される部分(例えば、マウスのボタン部分)に、血管画像90上のポインタP1の位置の材料特性に応じた応力が伝えらえる。
このため、ユーザが、入力部29を操作することによって、血管画像90における、材料特性によって示される機械強度の高い位置を指示するほど、入力部29における応力が大きくなる。このため、ユーザは、血管画像90における、材料特性によって示される機械強度の高い位置を指示するほど、強い応力を感じることとなる。このため、ユーザは、画面上の血管画像90を指示することで、実際の被検体の血管の材料特性に応じた血管の状態を把握することができる。
図2に戻り、第3制御部70Iは、血管モデルに応じた血管模型を造形するように、造形部72を制御する。
例えば、血管モデルが、図4(B)に示す血管モデルM102であったとする。
この場合、第3制御部70Iは、該血管モデルM102によって示される形状の血管模型D102を造形するように、造形部72を制御する。造形部72は、血管モデルM102に応じた血管模型D102を造形する。このため、造形部72は、図4(B)に示すように、血管の内径については、形状モデル40によって示される内径と略同じ血管模型D102を造形する。しかし、造形部72は、血管模型D102の血管壁の厚みについては、対応する位置の材料特性を実現可能な厚みに補正された、血管模型D102を造形する。
すなわち、血管模型D102は、被検体の血管を模した血管模型D102であって、被検体の血管に応じた内径であり、且つ、血管の内壁に対する外壁の相対位置が被検体の血管の材料特性を実現可能な位置とされた、血管模型である。
また、例えば、血管モデルが、図5(B)に示す血管モデルM100であったとする。
この場合、第3制御部70Iは、血管モデルM100によって示される形状で、且つ、血管モデルM100の付加情報46によって示される種類の造形材料で、血管模型D100を造形するように、造形部72を制御する。このため、造形部72は、図5(B)に示すように、血管の内径、外径、及び形状については、形状モデル40によって示される形状及び大きさと略同じであるが、位置ごとに材料特性を実現可能な種類の造形材料で造形した血管模型D100を造形する。
このとき、造形部72は、同一の材料によって血管模型を造形した後に、付加情報46によって示される位置に、該付加情報46によって示される種類の造形材料の特性を発現するように、光や熱などの刺激を付与することで材料を変性させ、血管模型D100としてもよい。
また、例えば、血管モデルが、図6(B)に示す血管モデルM104であったとする。
この場合、第3制御部70Iは、形状モデル40によって示される血管の模型における、血管の外側の、付加情報によって示される付与位置に、付加情報によって示される補強部材Pを付与した血管模型D104を造形するように、造形部72を制御する。
なお、この場合、第3制御部70Iは、形状モデル40によって示される形状部分については、上述した、模擬血管材料を用いるように造形部72を制御することが好ましい。一方、補強部材Pの部分については、上述した模擬血管材料の特性を満たす材料以外で造形してもよい。また、第3制御部70Iは、形状モデル40によって示される血管の模型を造形するように造形部72を制御し、補強部材Pの種類、および補強部材Pの付与位置を示す付与情報を、表示部31に表示してもよい。この場合、ユーザは、造形部72によって造形された模型に、表示部31に表示された付与情報によって示される補強部材Pを設置することで、血管模型D104を造形すればよい。
次に、処理装置70が実行する画像処理の手順を説明する。
図8は、処理装置70が実行する、血管モデル生成処理の手順を示すフローチャートである。
まず、第1取得部70Aが、被検体の血管に関する医用画像を取得する(ステップS100)。次に、構築部70Cが、ステップS100で取得した医用画像から、血管の三次元形状を示す形状モデルを構築する(ステップS102)。
次に、第2取得部70Bが、血管の各位置の材料特性を取得する(ステップS104)。次に、第1生成部70Dが、ステップS102で構築された形状モデルによって示される血管の各位置に、ステップS104で取得した材料特性に関する付加情報を付加した血管モデルを生成する(ステップS106)。
次に、第1生成部70Dは、ステップS106で生成した血管モデルと、ステップS102で構築した形状モデルと、被検体の識別情報と、を対応づけて、記憶部33(図1参照)に記憶する(ステップS108)。そして、本ルーチンを終了する。
図9は、血管模型の造形処理の手順を示すフローチャートである。
まず、受付部70Fは、造形指示を入力部29から受け付けたか否かを判別する(ステップS200)、ステップS200で否定判断すると(ステップS200:No)、本ルーチンを終了する。一方、ステップS200で肯定判断すると(ステップS200:Yes)、ステップS202へ進む。
ステップS202では、第3制御部70Iが、第1生成部70Dで生成した血管モデルに応じた血管模型を造形するように、造形部72を制御する(ステップS202)。そして、本ルーチンを終了する。
図10は、処理装置70が実行する表示制御の手順を示すフローチャートである。
まず、受付部70Fが、入力部29から血管画像の表示指示を受け付けたか否かを判断する(ステップS300)。例えば、第1制御部70Eは、生成された血管モデルと、血管モデルに対応する被検体の識別情報と、を対応づけて、これらの一覧を表示部31に表示する。ユーザは、入力部29を操作することで、表示部31に表示された、血管モデル及び識別情報の一覧の中から、血管画像の表示対象を選択する。すると、受付部70Fは、ユーザによって選択された識別情報を含む表示指示を、入力部29から受け付ける。
ステップS300で否定判断すると(ステップS300:No)、本ルーチンを終了する。一方、ステップS300で肯定判断すると(ステップS300:Yes)、ステップS302へ進む。
ステップS302では、第1制御部70Eが、ステップS300で受け付けた表示指示に含まれる識別情報に対応する、血管モデルおよび形状モデルを記憶部33から読取り、該形状モデルによって示される血管画像を表示部31に表示する制御を行う(ステップS302)。
次に、受付部70Fは、表示部31に表示された血管画像における何れかの位置の指示を受け付けたか否かを判断する(ステップS304)。ステップS304で否定判断すると(ステップS304:No)、本ルーチンを終了する。ステップS304で肯定判断すると(ステップS304:Yes)、ステップS306へ進む。
ステップS306では、第2生成部70Gが、ステップS304で受け付けた、血管画像における指示された位置に対応する付加情報によって示される、材料特性に応じた強度の出力信号を生成する(ステップS306)。
次に、第2制御部70Hは、ステップS306で生成された出力信号の強度に応じた応力を生じさせるように、入力部29を制御する(ステップS308)。そして、本ルーチンを終了する。
以上説明したように、本実施の形態の処理装置70は、第1取得部70Aと、第2取得部70Bと、構築部70Cと、第1生成部70Dと、を備える。第1取得部70Aは、被検体の血管に関する画像を取得する。第2取得部70Bは、血管の各位置の材料特性を取得する。構築部70Cは、医用画像から、血管の三次元形状を示す形状モデルを構築する。第1生成部70Dは、形状モデルによって示される血管の各位置に、材料特性に関する付加情報を付加した血管モデルを生成する。
このように、本実施の形態の処理装置70は、被検体の血管の材料特性に応じた血管モデルを生成する。
従って、本実施の形態の処理装置70は、各々の被検体の血管の材料特性に応じた血管モデルを提供することができる。
具体的には、血管内にカテーテルやステントなどを挿入する手術練習などを実行する場合、特に、有用である。すなわち、この血管モデルを、バーチャルシミュレータや血管模型の造形に用いることで、実際の検証対象の被検体の血管と同等の材料特性(弾性率、可撓性、強度、硬度、および、靱性、の少なくとも1つ)の血管模型やバーチャルシミュレータを用いて、手術練習などを実行することができる。
また、本実施の形態の処理装置70の第1生成部70Dは、血管モデルによって示される血管壁の厚みが、血管の各位置の材料特性に応じた厚みとなるように、血管壁の内壁に対する外壁の相対位置を補正した付加情報を付加した血管モデルを生成する。
このため、血管モデルは、血管の内径については、被検体の医用画像から構築した形状モデルによって示される血管の内径と略同じであるが、血管の外径については、対応する位置の材料特性を実現可能な厚みとなるように補正されたものとなる。
このため、本実施の形態の処理装置70は、上記効果に加えて、更に、各々の被検体の材料特性に応じた血管モデルを提供することができる。
また、第3制御部70Iが、血管モデルに応じた血管模型を造形するように、造形部72を制御することによって、血管モデルに応じた血管模型を造形することができる。
このため、造形部72は、血管の内径については、被検体の医用画像から構築した形状モデルによって示される血管の内径と略同じであるが、血管の外径については、対応する位置の材料特性を実現可能な厚みとなるように補正された血管模型を造形することができる。
すなわち、本実施の形態で造形された血管模型は、被検体の血管の材料特性と同等の、材料特性(弾性率、可撓性、強度、硬度、および、靱性、の少なくとも1つ)を有することとなる。
このため、本実施の形態の処理装置70は、上記効果に加えて、実際の検証対象の被検体の血管そのものを用いた状態に近い環境で、手術練習などを実行可能な、血管模型を提供することができる。
また、造形部72によって造形された血管模型を用いることで、被検体の血管のプラークなどの存在部位にカテーテルなどの医療器具を挿入した手術などの手術練習を行うことができる。また、本実施の形態において造形した血管模型は、被検体の血管と同等の材料特性を有していることから、血管内の病変部位に対して、使用するステント等の医療器具がどの程度血管拡張に有効であるか、どの位置まで血管拡張を行うことが有効であるか、などの被検体の血管ごとの検討を行うことができる。このため、手術を実施前に、治療効果の推定が可能となる。
また、第1生成部70Dは、血管モデルに応じた血管模型の各位置の造形材料が、対応する各位置の材料特性を満たす材料となるように、血管模型の各位置の造形材料の種類を示す付加情報を付加した、血管モデルを生成する。
このため、本実施の形態の処理装置70は、上記効果に加えて、更に、各々の被検体の材料特性に応じた血管モデルを提供することができる。
また、第3制御部70Iが、血管モデルに応じた血管模型を造形するように、三次元造形物を造形する造形部72を制御することによって、血管モデルに応じた血管モデルを造形することができる。
このため、造形部72は、血管の内径、外径、及び形状については、形状モデルによって示される形状及び大きさと略同じであるが、被検体の血管の材料特性を実現可能な種類の造形材料を位置ごとに用いた、血管模型を造形することができる。
また、第1生成部70Dは、血管モデルに応じた血管模型の各位置の強度が、対応する各位置の材料特性を満たすように、血管模型によって示される血管を外側から補強する補強部材の付与位置を示す付加情報を付加した血管モデルを生成する。
このため、血管モデルは、血管の内側から血管を操作したときの実際の材料特性を満たす機械特性を備えたものとなる。
このため、本実施の形態の処理装置70は、上記効果に加えて、実際の検証対象の被検体の血管そのものを用いた状態に近い環境で、手術練習などを実行可能な、血管模型を提供することができる。
また、第1制御部70Eは、形状モデルによって示される血管画像を表示部31に表示する。そして、第2生成部70Gは、表示された血管画像における、ユーザによって指示された位置に対応する付加情報によって示される材料特性に応じた強度の出力信号を生成する。第2制御部70Hは、出力信号の強度に応じた応力を生じさせるように、ユーザによって操作される入力部29を制御する。
このため、ユーザは、血管画像における、材料特性の高い位置を指示するほど、強い応力を感じることとなる。
従って、表示部31に表示された血管画像を介して、被検体の実際の血管の材料特性を触感的にユーザに提供することができる。
次に、実施の形態の処理装置70及び詳細を後述する血管解析装置50のハードウェア構成について説明する。図11は、処理装置70及び血管解析装置50のハードウェア構成例を示すブロック図である。
処理装置70及び血管解析装置50は、CPU800、ROM(Read Only Memory)820、RAM(Random Access Memory)840、HDD(Hard Disk Drive)(図示省略)、及び通信I/F(Interface)860を有する。CPU800、ROM820、RAM840、HDD(図示省略)、及び通信I/F860は、バスにより相互に接続されており、通常のコンピュータを利用したハードウェア構成となっている。
実施の形態の処理装置70及び血管解析装置50で実行される各種処理を実行するためのプログラムは、ROM820等に予め組み込んで提供される。
なお、実施の形態の処理装置70及び血管解析装置50で実行される上記各種処理を実行するためのプログラムは、これらの装置にインストール可能な形式又は実行可能な形式のファイルでCD−ROM、フロッピー(登録商標)ディスク(FD)、CD−R、DVD(Digital Versatile Disk)等のコンピュータで読み取り可能な記録媒体に記録されて提供するように構成してもよい。
また、実施の形態の処理装置70及び血管解析装置50で実行される上記各種処理を実行するためのプログラムを、インターネット等のネットワークに接続されたコンピュータ上に格納し、ネットワーク経由でダウンロードさせることにより提供するように構成してもよい。また、実施の形態の処理装置70及び血管解析装置50で実行される上記各種処理を実行するためのプログラムを、インターネット等のネットワーク経由で提供または配布するように構成してもよい。
実施の形態の処理装置70及び血管解析装置50で実行される上記各種処理を実行するためのプログラムは、上述した各機能部を含むモジュール構成となっている。実際のハードウェアとしてはCPU800がROM820等の記憶媒体から各プログラムを読み出して実行することにより上記各機能部が主記憶装置上にロードされ、主記憶装置上に生成されるようになっている。
(血管解析装置)
次に、血管解析装置50について詳細に説明する。
なお、実施の形態の医用画像診断装置は、心臓血管、頸動脈、または脳動脈等の人体のあらゆる部位の血管を解析対象とすることができる。以下、一例として、心臓の血管が解析対象であることとして説明を進める。
心臓の血管としては、例えば冠動脈と大動脈とが挙げられる。冠動脈は、大動脈の冠動脈起始部から始まり心筋表面を走行し、心外膜側から内膜側に入り込む。冠動脈は、心筋の内膜において無数の毛細管に分岐する。分岐後、無数の毛細管は、再び統合して大心静脈を形成し、冠静脈洞に接続する。冠血管系は、他の臓器と異なり、心筋の収縮および弛緩という力学的変化のなかで、灌流が保障されなければならないという点で特徴的である。
冠血流の特徴としては、心筋収縮による機械的血流阻害作用で冠動脈起始部の内圧が高くなる収縮期よりも、左心室拡張期に灌流圧が低下したときに多く流れることである。そのため、正常の冠動脈血流速波形は収縮期と拡張期の二峰性であり、拡張期血流が優位である。肥大型心筋症や大動脈弁狭窄症では、収縮期に逆行性波が認められ、大動脈逆流症では、収縮期順行波が大きくなる等、疾患によって特異な血流波形となることが知られている。また、拡張期の順行性波形は左室拡張機能、特に左室弛緩と密接な関係がある。左室弛緩遅延例では拡張期波形のピークが後ろにずれ、また減速脚がゆるやかになる傾向がある。また、このような症例では、頻拍時には拡張期の冠血流は十分に増大できず、心筋虚血を助長すると考えられている。
解剖学的に大動脈起始部から分岐する左右冠動脈に、大動脈圧に等しい冠灌流圧(すなわち、冠動脈が分枝する大動脈起始部の圧力)がかかることにより、冠血流が生じる。冠血流を決定するのは、大動脈圧である駆動圧と共に冠血管抵抗が重要である。140〜180μm以上の太い冠血管には、冠血管低抗の20%程度が存在するのに対し、100〜150μm以下の微小血管には、抵抗成分の残りの多くが存在するといわれる。従って、いわゆる冠狭窄等の無い場合には、抵抗値は冠微小血管の緊張性(トーヌス)に左右される。
血管抵抗因子としては、血管特性、動脈硬化、管狭窄、血液粘性、機械的因子があげられる。冠微小血管のトーヌスは、血管特性、心筋代謝(心筋酸素消費)、神経体液性因子、機械的因子、体液因子としての各種の血管作動性物質、血液粘性に規定され、さらに、心肥大、冠動脈硬化等を含めた様々な病変によっても影響され冠循環障害を起こす。
冠動脈血流拍動は、冠動脈血流の拍動パターン、心筋収縮による心筋内血流の制御、機械的刺激に対する心筋内血管の反応に影響される。心筋収縮が血流を阻害する機序としては、心筋内圧の上昇、心筋内血管容量の変化、心筋内血管の圧迫が挙げられる。心筋拡張期の血流規定因子には、拡張期の冠動脈圧、拡張期の血管外力、心拍数、心周期に占める拡張期の割合、心筋弛緩が存在する。
血管解析装置50は、時系列のCT画像に基づいて、力学モデルを構築し、力学モデルを利用して心臓の血管についての構造流体解析を実行し、血管内の力学的指標や血管流量指標を高精度に算出する。精度の高い力学的指標や血管流量指標を算出するためには、力学モデルに精度の高い潜在変数を割り当てる必要がある。血管解析装置50は、力学モデルを構築する際、初期的な力学モデルに逆解析を施して潜在変数を統計的に同定する。これにより血管解析装置50は、高精度の潜在変数を決定することができる。
力学的指標は、血管壁や血液に関する力学的な指標を意味する。血管壁に関する力学的指標としては、例えば血管壁の変位に関する指標、血管壁に生じる応力やひずみに関する指標、血管内腔に負荷される内圧分布に関する指標、血管の硬さ等をあらわす材料特性に関する指標等に分類される。血管の硬さ等をあらわす材料特性に関する指標は、血管組織の応力とひずみの関係を示す曲線の平均的な傾き等が挙げられる。
血液に関する力学的な指標における血管流量指標は、血管を流れる血液に関する血行動態の指標を意味する。血管流量指標としては、例えば血液の流量、血液の流速、血液の粘性等が挙げられる。なお、潜在変数は、例えば、血管の材料構成式、または血液の材料構成式といった材料モデルのパラメータ(例えばヤング率やポアソン比等)、血管内腔に負荷される内圧分布等の負荷条件パラメータ、構造解析や流体解析の境界条件パラメータ、時系列の形態指標や形状変形指標の不確定性に関連するばらつき分布パラメータの少なくとも一つを含む。
なお、材料モデルのパラメータ、材料パラメータは、上述した材料特性に相当する。
ここで、時系列の形態指標や形状変形指標の不確定性に関連するばらつき分布パラメータとは、医用画像データには、各CT値のノイズに起因したばらつき分布や、生体組織の境界閾値の曖昧性に起因する確率分布等が存在することを考慮して、血管組織や血液の境界座標および特徴点(血管分岐部や造影剤分散配置等)の空間座標における不確定性、または幾何学的構造パラメータ(芯線に垂直な断面の内腔半径等)の不確定性、または医用画像データ自体(CT値や境界閾値等)の不確定性を、確率分布として表現したものである。
力学モデルは、血管や血液の挙動を表現するための数値モデルである。力学モデルは、構造流体解析の手法に応じて異なるタイプを有している。例えば、力学モデルは、連続体力学モデルと簡易的力学モデルとに分類される。連続体力学モデルは、例えば、有限要素法(FEM:finite element method)や境界要素法に用いられる。簡易的力学モデルは、例えば材料力学に基づく材料力学モデルと流れ学に基づく流体力学モデルとに分類される。
なお、以下の説明において特に言及しない場合、力学モデルのタイプについては特に限定しないものとする。初期的な力学モデルは、潜在変数の確率分布や変数範囲から得られる潜在変数のパラメータに関するサンプリング集合(各パラメータの組み合わせの集合)が割り当てられた力学モデルを意味するものとする。
図12は、構造流体解析の対象領域(以下、解析対象領域と呼ぶ)に関する力学モデルM1の一例を示す図である。図12に示すように、力学モデルM1は、大動脈領域R1と右冠動脈領域R2と左冠動脈領域R3とを有している。血液は、大動脈から右冠動脈または左冠動脈へ流れる。
図12に示すように、力学モデルM1の大動脈起始部側の末端が血流の入口に設定され、右冠動脈領域の末端と左冠動脈領域の末端とが血流の出口に設定される。入口と出口との各々に境界条件が設定される。入口に関する境界条件は、例えば入口における血流の流速または血流による圧力、またはそれらの変化率を含む。出口に関する境界条件は、例えば、出口における血流の流速または血流による圧力、またはそれらの変化率を含む。
大動脈、右冠動脈、および左冠動脈の変形は、血流に起因する血管壁への力学的作用、心臓の拍動による血管壁への力学的作用(外力)、血管断面境界の負荷条件、血管壁の材料モデル、血管の無応力状態、および血管壁の幾何学的形状等の様々な要因に依存する。
ここで、血流に起因する血管壁への力学的作用は、例えば、血流に起因する内圧と血流に起因するせん断力とを含む。血流に起因する内圧により、血管半径方向あるは血管内腔面の垂直方向に変形する。心臓の拍動による血管壁への力学的作用と、血流に起因するせん断力は、血管芯線方向に関する伸縮、ねじり、曲げといった血管壁への力学的作用による変形を生じさせる。
血管芯線方向に関する伸縮、ねじり、曲げといった血管の全体的および局所的な変形挙動は、大動脈領域R1、右冠動脈領域R2、および左冠動脈領域R3、またはその一部、またはその周辺組織に、強制変位(移動ベクトルや回転変位)または荷重ベクトルの時間的変化を与えることで、負荷条件を設定する。また、血流に起因する内圧に基づく血管半径方向または内腔面に垂直な方向への変形は、血管内腔への圧力分布の時間的変化として与えることで、負荷条件を設定する。
大動脈領域R1、右冠動脈領域R2、および左冠動脈領域R3には、構造流体解析において、強制変位による変位拘束条件が割り当てられる。これにより、構造流体解析における血管壁の変形自由度を縮小することができ、計算収束性の安定化や、解析時間の短縮を実現できる。
また、例えば、血管の形状の変形度合は、血管壁の材料に依存する。そのため、大動脈領域R1、右冠動脈領域R2、および左冠動脈領域R3に材料モデルが割り当てられる。また、例えば、血管の形状の変形度合は、血管の無応力状態に依存する。負荷条件の初期値として、血管の残留応力分布を負荷条件として与えても良い。
ここで、血管解析対象の数値計算用力学モデルの空間を離散化した節点集合および節点から構成される要素集合において、材料モデルや境界条件や負荷条件等の解析条件を同定しない領域と、同定する領域にわけ、解析条件を同定しない領域内の節点には強制変位履歴の変位拘束条件を与え、材料モデルを同定する領域では、血管壁表面(外表面)の節点のみ強制変位履歴の変位拘束条件を与え、血管壁内部は変位自由度を確保し変位拘束を与えない。これにより、構造流体解析の変形自由度をおさえることができ、安定的かつ効率的に解析を行うことができる。
ただし、血管壁表面上に緩衝用のダミーの要素集合を設け、その表面の節点に強制変位を与えてもよい。これにより、血管内腔にプラーク等による突起がある場合や、血管分岐部等、内圧による負荷が、芯線方向の断面外の変形に影響を与える場合には、血管内腔だけではなく血管壁の形態指標も参照して血管に負荷される荷重ベクトルと内圧を分離して同定することができる。また、生理学的には壁表面の脂肪層を模擬しており、一方、数値計算上は、血管壁表面に強制変位を与えることで血管壁内部に局所的に実際とは異なる高い応力が発生することを避ける効果もある。
これら材料モデル、境界条件、および負荷条件等の潜在変数に関するパラメータは、後述する力学モデルに基づく逆解析(統計的同定処理)により同定される。逆解析により同定された精確な潜在変数は、力学モデルに割り当てられる。精確な潜在変数が割り当てられた力学モデルにより、解析対象血管領域外の血管や心臓等の外部要因による当該解析対象血管領域への影響を加味した構造流体解析または流体解析または構造解析または画像解析に基づく血行動態解析を実行することができる。
血管解析装置50は、力学モデルの構築に関し、逆解析による潜在変数の同定により、次の4点の困難を解決することができる。困難1.冠動脈の材料モデルの同定方法。困難2.心臓の形状の変形の冠動脈への影響の組み込み。困難3.冠動脈の境界条件の同定方法。困難4.医用画像データの不確定性に起因したばらつきを有する血管形状による画像解析や構造流体解析。この4点の困難の克服により、血管解析装置50は、逆解析による潜在変数の同定を行わない従来の血管構造流体解析に比して、解析精度の向上を実現する。
次に、実施の形態の医用画像診断装置における構造流体解析処理の詳細について説明する。図13は、システム制御部21の制御のもとに行われる構造流体解析処理の典型的な流れを示す図である。図14は、解析装置27の機能ブロック図である。
図13に示すように、構造流体解析処理においては、まず、システム制御部21により記憶部33から処理対象の医用画像ファイルが読み出され、解析装置27に供給される。医用画像ファイルは、時系列のCT画像の他に、当該被検体に関する血管内腔に関する圧力値のデータ、血液流量指標の観測値のデータ、およびプラーク指標を含んでいる。CT画像以外ではMRI画像や超音波エコー画像であってもよい。時系列のCT画像は、時系列のCT値の3次元空間分布を表現するデータである。時系列のCT画像は、例えば1心拍で約20枚、すなわち、約20心位相分のCT画像を含んでいる。
図13に示すように、システム制御部21は、解析装置27に領域設定処理を行わせる(ステップS1)。ステップS1において解析装置27の領域設定部51は、時系列のCT画像に含まれる血管領域に構造流体解析の解析対象領域を設定する。解析対象領域は、冠動脈に関する血管領域の任意の一部分に設定される。例えば、領域設定部51は、ユーザによる入力部29を介した指示、または、画像処理により血管領域に解析対象領域と同定対象領域を設定する。
ここで、図15を参照しながら、血管の構造について説明する。図15は、血管の芯線に略直交する断面(以下、血管断面と呼ぶ)を模式的に示す図である。図15に示すように、血管は、管状の血管壁を有している。血管壁の中心軸は芯線と呼ばれている。血管壁の内側は内腔と呼ばれている。内腔に血液が流れる。内腔と血管壁との境は血管内壁と呼ばれている。血管壁の外側には心筋等の血管周辺組織が分布している。血管壁と血管周辺組織との境は血管外壁と呼ばれている。血管壁内部にプラークが発生することがある。
図15に示すように、プラークは、例えば、石灰化した石灰化プラーク、粥状プラーク等に分類される。粥状プラークは、やわらかく、血管内壁が破裂して血栓として血管内部に染み出す危険性があり、不安定プラークと呼ばれることもある。従って、プラークの性状を把握することは臨床的に有用である。プラークの性状や存在領域は、医用画像ファイルに含まれるプラーク指標により特定可能である。プラーク指標は、例えば骨の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は、力学モデルの血管壁表面または血管壁内部または血管内腔における節点の変位の時間的変化から、血管の芯線方向に関する伸縮やねじりや曲げに関する変形を抽出し、血管芯線と芯線に垂直な断面における節点の強制変位として与えることで表現しても良い。このように、血管形状変形指標としては、力学モデルにおける各時刻の節点の強制変位データ(強制変位履歴)を特定する。
以下、図16および図17を参照しながら、画像解析・追尾処理を説明する。図16は、血管芯線の形態の時間的変化を示す図である。図16に示すように、例えば時系列の医用画像は、1心拍につき20枚のCT画像を含んでいるものとする。すなわち、心位相0%から95%まで5%おきにCT画像が得られることとする。画像解析・追尾処理部53により血管領域の芯線が抽出される。図16に示すように、芯線の形態は、心位相の経過に従って変化している。
図17は、時刻tと時刻t+Δtとの間における追尾処理の一例を示す図である。図17に示すように、血管芯線上にP1からP10の力学モデルにおける節点が設定されており、各断面上の血管の力学モデルの節点と力学的につながっている。ただし、血液の力学モデルの節点とは独立である。血管の追跡点の変位データをもとに、血管芯線上のP1からP20の節点の変位データを補間等の処理により算出し、各節点に強制変位が設定されているものとする。血管形状変形指標と血管形態指標について説明するため、節点P13と節点P14とにより規定される局所血管領域RAを考える。
時刻tにおいて、芯線方向に関する節点P13と節点P14との間の距離がLであり、血管領域の半径がrであるとする。画像解析・追尾処理部53から、節点P13と節点P14の血管芯線方向の伸縮、ねじり、または曲げ等の強制変位を抽出し、節点P13の強制変位(3次元空間における移動変位と、芯線方向の回転変位)と節点P14の強制変位(3次元空間における移動変位と、芯線方向の回転変位)を算出する。
図18は、血管芯線の曲げ変形および回転変位の算出例を示す図である。図18に示すように、例えばねじり角は、線aおよび線bから構成される面の法線方向ベクトルcの変化から算出しても良い。
図17および図18に示すように、画像解析・追尾処理部53は、追跡点の座標と移動ベクトルとに基づいて、芯線上の各節点の強制変位(3次元空間における移動変位と、芯線方向の回転変位)を算出し、血管形状変形指標を算出する。例えば、画像解析・追尾処理部53は、隣り合う2つの節点の座標差の時間変化を芯線方向に関する伸縮距離ΔLとして算出する。
また、画像解析・追尾処理部53は、各芯線上の節点について、当該節点と当該節点を含む血管領域断面上の他の節点(血管内腔または血管壁またはプラーク領域における節点)との間の距離の時間変化を半径方向に関する伸縮距離Δrとして算出する。また、画像解析・追尾処理部53は、各追跡点について、当該追跡点の近傍の複数の追跡点の座標と移動ベクトルとに基づいて、芯線上の当該節点の芯線方向のねじれ角度Δθを算出する。
また、画像解析・追尾処理部53は、血液領域の造影剤やプロトンの画像追尾により、血液流量指標として、流速または芯線方向断面の平均流速または平均流量を算出してもよい。
血管形状変形指標は、力学モデルにおける強制変位として利用される。以下、時系列の血管形態指標を「形状履歴」と呼び、時系列の血管形状変形指標を「強制変位履歴」と呼ぶことにする。
ステップS2が行われるとシステム制御部21は、解析装置27に構築処理を行わせる(ステップS3)。ステップS3において解析装置27の力学モデル構築部55は、形状履歴(時系列の血管形態指標)と、強制変位履歴(時系列の血管形状変形指標)と、時系列の医用画像(CT画像,MRI画像,超音波エコー画像等のDICOMデータ)とに基づいて、解析対象領域に関する力学モデルを暫定的に構築する。力学モデルは、構造流体解析を行うための解析対象領域に関する数値モデルである。
以下、ステップS3について詳細に説明する。力学モデル構築部55は、まず、医用画像と形状履歴とに基づいて、力学モデル(数理モデル)を解くための形状モデルを構築する。形状モデルは、各時刻における血管領域の幾何学的構造を模式的に表現したものである。形状モデルは、例えば複数の離散化領域に区分されている。各離散化領域の頂点は、節点と呼ばれる。
力学モデル構築部55は、時刻毎の医用画像に含まれる血管領域と血管形態指標とに基づいて、時刻毎の形状モデルを構築しても良いし、特定の時相の医用画像に含まれる血管領域と血管形態指標とに基づいて、時刻毎の形状モデルを構築しても良い。また、例えば、初期の負荷状態として、解析対象領域に対応する血管に残留応力が存在しないと仮定する場合、無応力状態の時相として、解析対象領域に対応する血管が最も収縮した時相を無応力状態であると仮定する。
図19は、形状モデルの芯線に直交する断面を示す図である。図19に示すように、形状モデルは、芯線から外側に向けて血管内腔領域、血管壁領域を有している。プラークが存在する場合、血管壁領域にプラーク領域を設けても良い。また、血管周辺組織による血管への影響を考慮する場合、血管周辺組織のダミー要素を血管壁領域の外側に設けても良い。
形状モデルが構築されると力学モデル構築部55は、各潜在変数の確率分布や変数範囲から得られる潜在変数のパラメータに関するサンプリング値(例えばマルコフ連鎖モンテカルロ法等による、各パラメータの組み合わせの集合からのサンプリング)を力学モデルに設定する。
例えば、力学モデル構築部55は、図12に示すように、大動脈領域R1の大動脈起始部側の末端に入口に関する境界条件の同定対象の領域(以下、境界条件同定領域)を設定し、右冠動脈領域R2の末端と左冠動脈領域R3の末端とに出口に関する境界条件同定領域を設定する。力学モデル構築部55は、各境界条件同定領域に境界条件の確率分布または変数範囲から得られる境界条件のパラメータに関するサンプリング値を割り当てる。
また、力学モデル構築部55は、大動脈領域R1、右冠動脈領域R2、および左冠動脈領域R3に材料モデルの同定対象の領域(以下、材料モデル同定領域と呼ぶ)および負荷条件の同定対象の領域(以下、負荷条件同定領域と呼ぶ)を設定する。力学モデル構築部55は、各材料モデル同定領域に材料モデルの確率分布または変数範囲から得られる材料モデルのパラメータに関するサンプリング値を割り当て、各負荷条件同定領域に負荷条件の確率分布や変数範囲から得られる負荷条件のパラメータに関するサンプリング値を割り当てる。
血管は、流量が0でも残留応力が存在すると言われている。例えば、力学モデル構築部55は、流量が0の場合の残留応力を負荷条件の初期値として解析対象領域に割り当てても良い。また、力学モデル構築部55は、幾何学的構造に不確定性がある部位に、幾何学的構造の同定対象の領域(以下、幾何学的構造同定領域と呼ぶ)を設定しても良い。なお、幾何学的構造のパラメータは、幾何学的構造の不確定性に関連するばらつき分布パラメータ、または医用画像データに内在するばらつき分布パラメータであり、各CT値のばらつき分布や生体組織の境界閾値のばらつき分布等であってよい。詳細は、後述するが、力学モデル構築部55は、プラーク領域に材料モデルを設定しても良い。材料モデルの詳細については後述する。
形状モデルが構築されると力学モデル構築部55は、ステップS2において算出された時系列の血管形状変形指標、すなわち、強制変位履歴を形状モデルに割り当てる。潜在変数および強制変位履歴が割り当てられた形状モデルを「力学モデル」と呼ぶことにする。
図20は、形状モデルM2は血管や血液の力学モデルの一部を示しており、力学モデルの節点への強制変位履歴の割り当てを説明するための図である。図20は、形状モデルM2の一部分を示している。ただし、図20では芯線がM2内に位置する場合を示しているが、芯線がM2外に位置する場合でも良い。図20に示すように、形状モデル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に時系列の血管形状変形指標を強制変位履歴として割り当てられる。これにより、血管の全体および局所に関する伸縮変形やねじれ変形や曲げ変形が表現される。
なお、強制変位履歴の割当対象は、芯線上の節点およびビーム要素に限定されない。図21および図22は、形状モデルへの強制変位履歴の他の割り当て方法を示す図である。図21は、血管形状変形指標がねじれの場合における割り当て例を示している。また、図22は、血管形状変形指標が曲げの場合における割り当て例を示している。図21および図22に示すように、力学モデル構築部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において算出された血液流量指標の予測値とが、事前に収集された血管形態指標の観測値と血液流量指標の観測値とに整合するように力学的モデルの潜在変数のパラメータを統計的に同定する。
図14に示すように、統計的同定部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は、各潜在変数について、事後分布の最頻値や平均値等の統計値から各潜在変数のパラメータを同定する。例えば、上述の例の場合、血管内腔圧力値に関する事後分布が算出され、この事後分布から血管内腔圧力値の同定値が算出される。
図23は、階層ベイズモデルおよびマルコフ連鎖モンテカルロ法による負荷条件(血管内の平均圧力)に関する事後分布算出と平均内圧の同定とを説明するための図である。図23に示すように、血管起始部から延びる血管に石灰化プラーク領域と粥状プラーク領域とが設定されているものとする。石灰化プラーク領域は材料モデル同定領域に設定され、粥状プラーク領域に材料モデル同定領域に設定される。血管起始部から血管芯線方向に沿って進むにつれて血管内圧が降下する。血管芯線に沿って複数の節点が設定される。各節点を含む直交断面(節点断面)において内腔内圧の事後分布が算出され、事後分布の最頻値が同定される。
なお、血管形態指標の観測値としては、例えば、ステップS2において算出された血管形態指標が用いられる。
第2統計的同定部61−2による処理は、データ分布の算出に用いる指標が異なるだけで第1統計的同定部61−1による処理と同様である。すなわち、第2統計的同定部61−2は、まず、ステップS5において算出された血液流量指標の予測値と観測値とに基づくデータ分布を設定する。
次に、第2統計的同定部61−2は、力学モデルの潜在変数に事前分布を割り当てる。例えば、血管に関する材料モデルのパラメータや血液に関する材料モデルのパラメータ、プラークに関する材料モデルのパラメータに関する事前分布が割り当てられる。これら材料モデルのパラメータとしては、例えば、弾性率等の材料モデルパラメータや、血液の構成式における粘性に関するパラメータが挙げられる。材料モデルのパラメータの想定範囲や確率分布は、経験的に予め設定することができる。
第2統計的同定部61−2は、各離散化領域について材料モデルのパラメータの確率分布、すなわち、事前分布を設定し、これら想定範囲内に限定して、想定した確率分布に従って、材料モデルのパラメータに関するモンテカルロシミュレーションを実行することができ、力学モデルへ設定するための材料モデルパラメータ(潜在変数)のサンプリング値を得ることができる。
次に、第2統計的同定部61−2は、各潜在変数について、事前分布とデータ分布とに統計的同定処理を施すことにより事後分布を算出し、算出された事後分布の統計値から各潜在変数のパラメータを同定する。例えば、上述の例の場合、材料モデルのパラメータに関する事後分布が算出され、この事後分布から材料モデルのパラメータの同定値が算出される。
図24は、階層ベイズモデルおよびマルコフ連鎖モンテカルロ法による材料モデルパラメータに関する事後分布算出と材料モデルのパラメータ(血管壁の等価弾性率)の同定とを説明するための図である。図24に示すように、血管モデルは、図23と同様であるとする。材料モデル同定領域に限定して、血管壁の材料モデルのパラメータ(例えば、等価弾性率)の事後分布が算出され、事後分布の最頻値が同定される。
なお、血液流量指標の観測値は、例えば、大動脈に送り出される血流量変化であると仮定し、血管形態指標の観測値を、時系列の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は、予測値に応じたカラー値をカラーテーブルを利用して特定し、特定されたカラー値に応じた色で当該予測値に対応する離散化領域を表示する。
図25は、力学的指標の一つである内圧の空間分布の表示例を示す図である。図25に示すように、表示部31は、力学モデルを構成する各離散化領域を当該離散化領域に関する内圧値に応じた色で動画的に表示する。ユーザは、力学モデルを観察することにより、経時的且つ空間的に変化する力学的指標を色で把握することができる。
図26は、血液流量指標の一つである流速値の空間分布の表示例を示す図である。図26に示すように、表示部31は、力学モデルを構成する各離散化領域を当該離散化領域に関する流速値に応じた色で動画的に表示する。ユーザは、力学モデルを観察することにより、経時的且つ空間的に変化する血液流量指標を色で把握することができる。
例えば、血管内が完全に狭窄している場合、狭窄部位の内圧は、非狭窄部位の内圧よりも小さい。力学的指標として内圧が指定された場合、ユーザは、力学モデル上での局所的な色の違いにより狭窄の有無を判断することができる。また、狭窄部位の流量は、非狭窄部位の流用よりも小さい。血液流量指標として流量が指定された場合、ユーザは、力学モデル上での局所的な色の違いにより狭窄の有無を判断することができる。
また、血管応力解析部57は、プラーク領域の材料モデルのパラメータの同定結果に基づいて、力学的指標として硬さ値の空間分布を算出しても良い。この場合、表示部31は、プラーク領域に関する硬さ値の空間分布を力学モデル上において表示しても良い。また、表示部31は、プラーク領域周辺の内圧分布や応力分布、ひずみ分布を表示して良い。ユーザは、これらの表示をプラークの性状と破綻しやすさとを推定することに活用することができる。
力学的指標および血液流量指標の予測値は、力学モデルの離散化領域の色で表現する方法のみに限定されない。例えば、図27、図28、および図29に示すように、グラフで表示しても良い。なお、図27は、左冠動脈起始部の血圧に関するグラフである。図27のグラフの縦軸は正規化した血圧に規定され、横軸は心位相[%]に規定される。図28は、LCXとLDAとの分岐点付近の血圧に関するグラフである。図28のグラフの縦軸は正規化した血圧に規定され、横軸は心位相[%]に規定される。図29は、芯線方向に関する血圧変化に関するグラフである。図29のグラフの縦軸は血圧比に規定され、横軸は大動脈からの距離[mm]に規定される。表示部31は、力学的指標および血液流量指標の予測値をグラフで表示することにより、これら値をユーザに簡便に把握させることができる。
ステップS11が行われると構造流体解析処理が終了する。
なお、図20において、強制変位履歴は、形状モデルの芯線部と外壁部とに設定されるとしたが、強制変位履歴の設定箇所は、これに限定されない。例えば、強制変位履歴は、芯線部と外壁部との間の血管壁領域に設定されても良い。
また、強制変位履歴の拘束条件の割り当て対象は、境界条件および材料モデルを同定するか否かに応じて切り分けられても良い。図30は、強制変位履歴の他の割り当て例を示す図であり、形状モデルの断面を示している。例えば、図30(a)に示すように、境界条件および材料モデルを同定する場合、形状モデルの外壁部OW上の節点PN2のみに強制変位履歴を割り当て、血管壁領域RVの節点PN3には強制変位履歴を割り当てなければよい。
また、図30(b)に示すように、境界条件および材料モデルを同定しない場合、形状モデルの外壁部の節点PN2と血管壁領域RVの節点PN3との両方に強制変位履歴を割り当てればよい。この場合、芯線上の節点PN1に強制変位履歴が割り当てられる。また、外壁部OWの節点PN2と節点PN1とをビーム要素EBで結び、ビーム要素EB上の節点PN2およびPN3にも強制変位履歴を割り当ててもよい。このとき、半径方向に関する収縮および膨張は、ビーム要素EBの伸縮変位で表現する。なお、血管内腔領域RIには、強制変位履歴が割り当てられなくて良い。
図31は、強制変位履歴の割り当ての他の例を示す図であり、血管周辺組織のダミー要素RDを含む形状モデルの断面を示している。図31に示すように、ダミー要素RDは、血管壁領域RNの外側に設定される。形状モデルがダミー要素RDを含む場合、血管壁領域RNに加え、ダミー要素RDにも節点PN4が設定される。節点PN4にも強制変位履歴が割り当てられる。
力学モデル構築部55は、境界条件および材料モデルを同定する場合、血管壁領域RVに含まれる節点PN3に強制変位履歴を割り当て、境界条件および材料モデルを同定しない場合、強制変位履歴を割り当てなくても良い。節点PN3に強制変位履歴を割り当てる場合、内腔領域RIに関する形状指標以外にも血管壁領域RVに関する形状指標も参照して材料モデルの同定が行われる。
図32は、強制変位履歴の割り当ての他の例を示す図であり、プラーク領域RPを含む形状モデルの断面を示している。図32に示すように、プラーク領域RPは、血管壁領域RVに含まれる。プラーク領域RPは、材料モデル同定領域に設定される。プラーク領域RPについては、内腔形状指標、血管壁形状指標、およびプラーク指標を考慮して材料モデルが同定される。
既述のように、プラーク指標は、例えば、超音波診断装置による組織性状診断により得られたプラークの性状に関するデータである。力学モデル構築部55は、性状に応じてプラーク領域を複数の部分領域に区分し、複数の部分領域に個別に材料モデル同定領域を設定する。各部分領域には、予め、当該部分領域の性状に応じたパラメータ範囲を設定することが好ましい。既述の統計的同定処理により、各部分領域についての材料モデルパラメータが統計的同定部61により同定される。そして、ステップS10において表示部31が力学的指標として血管の材料特性に関する指標を表示することにより、ユーザは、プラークの性状を正確且つ容易に把握することができる。
次に、潜在変数の一つである材料モデルについて詳細に説明する。血管の材料モデルとしては、弾性モデル、超弾性モデル、異方性超弾性モデル、粘性特性を考慮した超弾性モデル等が適用可能である。異方性超弾性モデルとしては、例えば、Y.C.Funにより提案された数理モデルや、Holzapfel−Gasser構成式と呼ばれる数理モデルを適用できる。
単位参照体積あたりのひずみエネルギーは、以下の(1)式で表わされる。(1)式の第1項は、コラーゲンを含まない等方性基礎材料のせん断変形に関するエネルギーを表している。また、(1)式の第2項は、コラーゲンを含まない等方性基礎材料の体積変形に関するエネルギーを表している。また、(1)式の第3項は、コラーゲン繊維各グループの寄与(繊維方向の分散を考慮)を表している。
図33は、繊維グループの変形を示す図である。図33に示すように、円筒形上の外膜を仮定する。芯線方向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」等の境界面補足型の手法でもよい。
次に、簡易的な力学モデルの例として、内圧と外圧とを受ける厚肉円筒の材料力学の式と、ハーゲン・ポアズイユ流れおよび修正ベルヌーイの式とについて詳細に説明する。
まず、図34と図35とを参照しながら、厚肉円筒の材料力学の式について説明する。図34は、肉厚円筒の力学モデルの直交断面を示す図である。図35は、図34の微小扇形要素の拡大図である。内半径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画像に描出されない血管や心臓等の外部要因による影響を加味した構造流体解析を実行することができる。このように、実施の形態の医用画像診断装置によれば、血管の構造流体解析の精度の向上を図ることができる。
そして、血管解析装置50は、同定した潜在変数を形状モデルに割当てることで、最終的な力学モデルを構築する。さらに、血管解析装置50は、構築した力学モデルを解析することで、被検体の血管の各位置の材料特性(材料モデル、材料パラメータ)を同定する。そして、同定した材料特性を、処理装置70へ出力する。
以上、実施の形態を説明したが、実施の形態および変形例は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施の形態および変形例は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。実施の形態およびその変形は、発明の範囲や要旨に含まれると共に、特許請求の範囲に記載された発明とその均等の範囲に含まれる。