(57)【要約】【課題】心臓外科医に対する心臓触診の訓練を支援するシステムを提供する。【解決手段】この心筋触診訓練システムは、リアルタイム処理を行うため、仮想心臓表示と心臓弾性計算処理とに分けられ、２台のＰＣによる並列処理を行なう。このシステムは、仮想心臓処理用ＰＣ１０、力覚伝達機構を構成するＤＣサーボモータとポテンショメータをリアルタイム制御するためのＰＣ１２、仮想心臓の弾性力を訓練者１４の手に伝達するための力覚伝達機構１６、力覚伝達機構のモータを駆動するモータドライバ１８によって構成されている．ＰＣ１０とＰＣ１２とは、データ通信ボード８によって接続される。【選択図】    図１
2004181176.tif 000002
【特許請求の範囲】
【請求項１】
横方向の心筋要素モデルと縦方向の心筋要素モデルと質量モデルを横方向の心筋要素モデルは碁盤目状の交点相互を結ぶように、縦方向の心筋要素モデルは碁盤目状の交点と地盤を結ぶように、質量モデルは碁盤目状の交点に、平面状に配設した心筋モデルを作成し、かつ健常者や心筋梗塞の患者の心筋の個人差に合わせて変更が出来るパラメータを設定することによって、健常者や心筋梗塞の患者の心臓の筋肉組織を模擬させるプログラムをＰＣ内に入れておいて、
アクチュエーターのアタッチメントを指で触ったことにより生じた変位を計測し、この変位量をＰＣに入力して、変位量の数値に対応して健常者や心筋梗塞の患者の筋肉組織を指で触った時に感じる反作用と近似した反作用を計算させ、この反作用の計算値をアクチュエーターへ送信し、アクチュエーターのアタッチメントを駆動させることによって、心筋触診と近似の感触を得るようにしたことを特徴とする心筋触診訓練システム。
【請求項２】
心室表面に複数個の座標を取得し、座標間をワイヤーフレームで結び、一心拍複数フレームを連続描画させることで仮想心臓の拍動運動の動画とするシステムにおいて、各座標の間を等分した座標を想定して、スプライン補間を行い、ワイヤーフレームの形状を滑らかにするとともに、
横方向の心筋要素モデルと縦方向の心筋要素モデルと質量モデルを横方向の心筋要素モデルは碁盤目状の交点相互を結ぶように、縦方向の心筋要素モデルは碁盤目状の交点と地盤を結ぶように、質量モデルは碁盤目状の交点に、平面状に配設した心筋モデルを作成し、かつ健常者や心筋梗塞の患者の心筋の個人差に合わせて変更が出来るパラメータを設定することによって、健常者や心筋梗塞の患者の心臓の筋肉組織を模擬させるプログラムをＰＣ内に入れておいて、
アクチュエーターのアタッチメントを指で触ったことにより生じた変位を計測し、この変位量をＰＣに入力して、変位量の数値に対応して健常者や心筋梗塞の患者の筋肉組織を指で触った時に凹むのと近似した変形量を計算させ、この変形量の計算値をディスプレーへ送信し、ディスプレー側では、ワイヤーフレームの形状を滑らかにするためにスプライン補間したことによって出来た座標を用い、指の中心点を仮想心臓の適切な座標に決め、その中心点を中央とするパッチ状のワイヤーフレーム画像を形成させ、送信されてきた変形量の計算値によって、拍動運動の動きに加えて心筋触診時と近似の変形した動画が描画されるようにしたことを特徴とする心筋触診訓練システム。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、心臓外科医に対する心臓触診の訓練を支援するシステムに関する。
【０００２】
【従来の技術】
心臓外科手術、特に左室形成術において心筋触診は手術プランニング上重要な役割を果たす。ＭＲＩ、ＰＥＴ、ＣＴ、超音波やＸ線などの医用画像診断データと医師の経験に基づいて、術前の心臓疾患の定性的推定はある程度可能であるが、最終的な手術プランニングの決め手は実際の心筋触診であり、心臓外科医は指先の感覚を頼りに心筋梗塞で薄くなった心筋部を認識しながら切開縫合を行なう。この左室形成術は胸部切開後、人工心肺下において行われるため、心筋触診の機会は術中の短時間のみである。経験の浅い医師にとって、触診による手術プランニング学習とその訓練の機会はこの短時間の一部のさらに短時間のみである。熟練した医師にとっても、触診による手術プランニングと手術との間に時間的な余裕がない。高齢化が進むわが国において、心臓外科手術の安全性向上は心臓外科にとって重要な課題である。このようなことから、心臓外科医より術前の触診が可能となる模擬触診システムの開発が待望されている。
１９９０年代後半から、ロボット技術を用いた手術支援システムおよび手術訓練システム技術が急速に世界各地で開発されてきている。その手術対象の多くは頭部、腹部、胸部、目などにある静止した臓器である。ロボット技術を用いている手術訓練システムは、実臓器の代わりに画面上のモデル臓器を用い、訓練者はアクチュエーター付きの手袋であるデータグローブを通じて画面の臓器に触れるようになっているが、その時、微小反作用が訓練者に提示されないだけでなく、標準的な臓器モデルを用いて個別患者の臓器モデルを用いないため、術前の訓練システムとしては余り役に立っていない。
またこれまで、教育用にせよ臨床用にせよ動く臓器である心臓の触診訓練システムの報告はない。
【０００３】
【発明が解決しようとする課題】
以上述べたように心臓外科医より、術前の心筋触診訓練が可能となる心筋触診訓練システムの開発が待望されているにもかかわらず、いまだに開発が実現していないのは技術的に難しい問題点があるからである。
【０００４】
心筋触診訓練システムで必要なことは、視覚と力覚とによって、健常な人の心臓、心筋梗塞患者などの人の心臓の心筋触診の臨場感を実現することである。
視覚については、ディスプレー上に心臓の動画像を表示させ、指でアクチュエーターのアタッチメントを触った時、その変位を計測し変位信号をＰＣに送り、ＰＣによって心臓筋肉表面の張り具合や緩み具合の状態によって生じる指の周辺の凹み具合を計算して、その結果を時間遅れなくディスプレーに送って、心臓の動画像の指で触った部分の凹む映像が、力覚の動きと整合した動画像で表示されることが出来れば臨場感を実現できる。
力覚については、指でアクチュエーターのアタッチメントを触った時、その変位を計測し変位信号をＰＣに送り、ＰＣによって心臓筋肉表面の張り具合や緩み具合の状態を感じ取れる力覚を発生するように計算して、その結果を時間遅れなくアクチュエーターに返して、力に変換してアタッチメントから指に心臓筋肉表面の力覚を伝えるこが出来れば臨場感を実現できる。
【０００５】
心臓の動きを画面表示することは写真技術を使用すれば簡単に出来るが、各種の触診の状況を再現出来ない上、指の力覚の再現も出来ない。
ＰＣの制御が可能な形で心臓の動きを画面表示するには、先ず三次元座標を取得しなければならない。静止した臓器であれば簡単であるが心臓は動く臓器であり、困難である。また、後述のＰＣ内での心筋の動きと力覚を得るための計算を行なった結果に対応した情報を三次元座標に送って心筋表面の変形を再現するには、送り側と受け入れ側のマッチングを図らなければならないのであるが、相手が筋肉から得られる情報であり、このマッチングが困難である。さらに、画像の動きは、滑らかなものしなければならず、かつリアルタイムに表示させなければならないので、情報処理が困難である。
【０００６】
心臓筋肉表面の張り具合や緩み具合に応じた画像の動きと、その状態を感じ取れる力覚とを発生するには、ＰＣ内で処理して、計算結果を出力する必要がある。この計算が出来るようにするためには、筋肉を模擬出来る力学モデルを作らなければならない。筋肉というものは、多数の細い筋肉繊維が多層構造で構成されている。従来の技術では、このような複雑な筋肉組織を模擬するためには有限要素法という技術を使うのが一般的である。これは模擬するものの全体を有限個の要素に分割する計算技術である。例えば５［ｃｍ３］の肉片であったとすれば、それを１０万個の要素に分割して、各要素の挙動を計算する手法である。
しかし、有限要素法では要素の数が多いため計算量が必然的に増加する。計算量が大きいと時間遅れが発生する。触診システムでは、手の接触に対する変形量や反作用を時間遅れなく計算する必要があるのに対して、有限要素法では、時間遅れがあって対応できないのである。
【０００７】
本発明は、このような心筋触診訓練システムを開発しようとする時に障害となる問題点を解決して、触診力覚構築部は実際の心臓触診に近い感触が得られるように、仮想心臓構築部は実際の心臓の拍動に近い画面表示が出来るように、アクチュエーターのアタッチメントに指が触れて指が心臓を押したことによる心臓の表面の凹みが実際の心臓触診時の心臓の表面の凹みに近い姿でディスプレー上に描画されるようにして、このシステムを使用する人が心臓触診の臨場感を体感出来るように、そして手入力によって、健常者や心筋梗塞の患者の心筋の個人差に合わせて、その出力の変更が出来るようにした心筋触診訓練システムを提供することを目的とするものである。
【０００８】
【課題を解決するための手段】
そして、本発明は上記目的を達成するために、上記課題に対して次のような手段を講じる。
【０００９】
心室表面に複数個の座標を取得し、座標間を直線で結び、仮想心臓はワイヤーフレームで構築される．一心拍複数フレームを連続描画させることで仮想心臓の拍動運動の動画とするシステムにおいて、各座標の間を等分した座標を想定して、スプライン補間を行ない、ワイヤーフレームの形状を滑らかにするとともに、横方向の心筋要素モデルと縦方向の心筋要素モデルと質量モデルを横方向の心筋要素モデルは碁盤目状の交点相互を結ぶように、縦方向の心筋要素モデルは碁盤目状の交点と地盤を結ぶように、質量モデルは碁盤目状の交点に、平面状に配設した心筋モデルを作成し、かつ健常者や心筋梗塞の患者の心筋の個人差に合わせて変更が出来るパラメータを設定することによって、健常者や心筋梗塞の患者の心臓の筋肉組織を模擬させるプログラムをＰＣ１０内に入れておいて、アクチュエーター１６のアタッチメントを指で触ったことにより生じた変位をＰＣ１２で計測し、この変位量をＰＣ１０にデータ送信して、変位量の数値に対応して健常者や心筋梗塞の患者の筋肉組織を指で触った時の凹みに近似した変形量を計算させ、この変形量の計算値をディスプレーへ送信し、ディスプレー側では、ワイヤーフレームの形状を滑らかにするためにスプライン補間したことによって出来た座標を用い、指の中心点を仮想心臓の適切な座標に決め、その中心点を中央とするパッチ状のワイヤーフレーム画像を形成させ、送信されてきた変形量の計算値によって、拍動運動の動きに加えて心筋触診時と近似の変形した動画が描画されるようにする。
心臓は左心室が拡張と収縮を繰り返すことで全身に血液を送っているため、心臓形状はその動きを繰り返して常に動いている。そこで拍動１周期を等時間間隔に分け、それぞれの時の形状ごとに三次元座標を取得する。三次元座標取得の方法として、高速シネＣＴ画像より取得したヒト拍動心臓の断面動画を利用する。リアルタイム制御に支障をきたすという問題を克服する心筋の三次元力学高速モデルとして、心筋要素モデルを設ける。心筋要素モデルは質量、バネ、ダンパーによって構成する。心筋を解剖学知見で観察すると、横方向のつながりを形成する筋繊維による層と縦方向のつながりを形成する層が重なって出来る心筋壁が形成されて心筋となっていることに気付く。このことを足がかりに心筋を横方向と縦方向で異なった心筋要素モデルとして模擬した。
力学モデルに用いられる機械的要素はバネとダンパーである。心筋の横方向モデルは、基本的にはバネとダンパーによる直列モデルとバネによる並列モデルで構成される。前述のように力学モデルは縦方向と横方向とで異なったモデルを用いた。そして、心筋収縮期に心筋が硬くなる生理的性質をモデルとして組み込むため、横方向の心筋モデルは拍動周期中にモデルを変化させる。
各力学モデルのパラメータの同定は、横方向パラメータは拍動心臓の血流データから行ない、縦方向パラメータは摘出心臓を切り開き、力センサを用いて弾性データを取得した。これを初期データとして動きや力などの関係式を作成し、パラメータに数値を入れ、結果を実証するという、式の作成と実証試験との繰り返しによって、人の心臓の拍動運動の姿の再現、力覚の再現が可能となったもので、本力学モデルは、有限要素法から得られるものと同等の再現性があり、かつ計算速度は実用性のある速さである。
【００１１】
仮想心臓の表面は点を格子状に配置することで表現している。そこに力学モデルパッチを仮想心臓の表面を構成するうちの例えば９×９の点を指の接触点を中心として当てはめる。このパッチの範囲が力学モデルの対応する範囲としてみなす。次に力学モデルは質量をバネとダンパーで構成される心筋要素を結んだモデルとなっており、質量の位置変化により変形を知ることができる構成となっているが、このままでは力学モデルパッチの点の数と力学モデルの点の数が合わないため、力学モデルの質量の間を前述の表面の補間と同じやり方でスプライン補間法を用いて９×９の変化量を求める。９×９の変化量が求まれば、力学モデルパッチの対応する表面部分において、指の接触している表面の構成点から表面座標を求める際に決めたｚ軸への方向を力学モデルの表面に対して垂直な方向として、パッチの対応する表面の９×９の点をそれぞれ対応する変化量分移動させて仮想心臓を描画することで表面の変形を行なう。
このようにして、人の心臓のＭＲＩからとったデータを元にして作ってきた仮想心臓構築側と動物からとったデータを元にして作ってきた力学モデルとのマッチングの問題は解決される。
また画像の表面を滑らかにすることは、スプライン補間法を用いることで解決される。
【００１２】
【発明の実施の形態】
以下、本発明の実施の一形態について、図面を参照しながら説明する。
この心筋触診訓練システムは、技術面で大きく分けると（Ａ）仮想心臓構築と（Ｂ）弾性伝達機構構築の二つの技術要素とに分別できる。本システムは、リアルタイム処理を行なうため、仮想心臓表示と心臓弾性計算処理とに分けられ、２台のＰＣ（パーソナルコンピュータ）による並列処理を行なう。そして本システムは、図１に示すように、仮想心臓処理用ＰＣ１０、力覚伝達機構を構成するＤＣサーボモータとポテンショメータをリアルタイム制御するためのＰＣ１２、仮想心臓の弾性力を訓練者１４の手に伝達するための力覚伝達機器１６、力覚伝達機構のモータを駆動するモータドライバ１８によって構成されている。ＰＣ１０とＰＣ１２とは、データ通信ボード８によって接続される。
【００１３】
利用の例として、熟練医師がグラフィカルユーザインターフェイスを利用して、過去の患者のデータを入力する事項について述べる。グラフィカルユーザインターフェイスは、以下ＧＵＩと記述する。ここで入力されるデータは、患者の体重、心拍数（６０〜１８０［ｂｐｍ］）、左心房圧（０〜１５［ｍｍＨｇ］）、病変部設定入力（力学パラメータの補正係数、心筋部位による心筋の弾性、及び収縮力を補正する。仮想心臓上の色で判断できる。）である。熟練医師によって設定された仮想心臓はＧＵＩ上で拍動運動を続け、その様子はＰＶ線図と体循環パラメータと共に観察することができる。また仮想心臓は訓練者が希望する視点や大きさで自由に観察できる。訓練者は、人差し指と親指に力覚伝達マニピュレータを装備し、腕と手首、指の動きで仮想心臓に触ることができる。力覚伝達マニピュレータを通して仮想心臓に触れた指には仮想心臓の拍動に沿った弾性がモニターの視覚情報に時間遅れなく伝達される。ここで視覚情報とは、仮想心臓の指による変形である。
【００１４】
仮想心臓処理用ＰＣ１０はＲｅｄｈａｔ−Ｌｉｎｕｘ  ＯＳを利用している。ＰＣのスペックはＣＰＵはＤｕａｌＣＰＵであり、Ａｔｈｒｏｎ１８００＋を用いている。メモリサイズは５１２ＭＢである。また、力覚伝達機器制御用ＰＣ１２はＲＴ−ＬｉｎｕｘＯＳを利用している。ＰＣのスペックは、ＣＰＵにＰｅｎｔｉｕｍＩＩＩ  １ＧＨｚ、メモリサイズは３８４ＭＢである。システム実行方法について述べる。システムを利用するためのフローチャートを図２に示す。図２のフローチャートを参照しながら、制御開始のために仮想心臓処理用ＰＣ１０と力覚伝達機構制御用ＰＣ１２の電源をそれぞれＯＮにした後の仮想心臓処理用ＰＣ１０のプロセス、力覚伝達機器制御用ＰＣ１２のプロセスについて述べる。
【００１５】
仮想心臓処理用ＰＣ１０のプロセスについて述べる。まず電源ＯＮとして、端末処理ＷｉｎｄｏｗであるＫ−ｔｅｒｍを立ち上げる。プログラムが保存されたフォルダに移動する。フォルダ内には、メインプログラムが書かれたｍａｉｎ．ｃ、仮想心臓描画プログラムが書かれたｆｒａｍｅ．ｃ、血液体循環プログラムが書かれたｂｏｄｙ．ｃが格納されており、Ｋ−ｔｅｒｍ上で．／ｇｃｃと入力することでプログラムがコンパイルされ、実行ファイルが生成される。実行ファイル名はｒｕｎであり、．／ｒｕｎと打ち込むことで実行される。実行されると心筋触診訓練システムのＧＵＩが画面に表示される。画面の様子は、図３に模式図として示している。
【００１６】
力覚伝達機器制御用ＰＣ１２のプロセスについて述べる。まず電源ＯＮとして、管理者（ｒｏｏｔ）でログインする、これはリアルタイム制御のためにＰＣのモジュール操作を行う必要があるため、一般ユーザに制限を加えてからである。端末処理ＷｉｎｄｏｗであるＫ−ｔｅｒｍを立ち上げ、プログラムが保存されたフォルダに移動する。メインプログラムが書かれたｅｘｅ．ｃ、制御プログラムが書かれたｍｏｄ．ｃがある。Ｋ−ｔｅｒｍ上でｍａｋｅと打ち込み、ｍｏｄ．ｏというモジュールができる。ｉｎｓｍｏｄ．ｏと打ち込むことでＰＣにリアルタイム制御モジュールが組み込まれる。次に．／ｇｃｃ  ｅｘｅ．ｃと打ち込むとａ．ｏｕｔという実行ファイルが形成される。この実行ファイルを．／ａ．ｏｕｔと入力して実行することでシステムは制御開始となる。
【００１７】
図３の心筋触診訓練システムのＧＵＩ画面の模式図と図４のフローチャートを参照しながらＧＵＩの構築方法について述べる。ＧＵＩはＬｉｎｕｘのグラフィカルライブラリＧＴＫ＋を利用して構築している。本ＧＵＩは（ａ）仮想心臓、（ｂ）メニューバー、（ｃ）ＰＶ線図、（ｄ）病変部設定Ｗｉｎｄｏｗ、（ｅ）血液循環パラメータ表示枠、（ｆ）訓練者の指と指へ伝達される弾性情報の表示枠、（ｇ）仮想心臓視点設定パネル、（ｈ）スタートストップボタンパネルが搭載されている。これらの表示部分はＷｉｄｇｅｔ（ウィジェット）と呼ばれ、ＧＴＫ＋ライブラリのコンポーネントである。仮想心臓表示はＦｒｅｅの三次元描画ソフトＯｐｅｎＧＬグラフィックスライブラリを利用している。各ウィジェット、ボタンやフレームの名称はすべて１つ１つそれぞれのｌａｂｅｌｗｉｄｇｅｔとして定義する必要がある。プログラム記述はｍａｉｎ．ｃのＧＴＫ＋コード記述部分において、始めにＧＵＩ全体のＷｉｎｄｏｗサイズを指定する。続いてＷｉｎｄｏｗに割り当てるＴａｂｌｅを縦と横の枠数を指定する。使用するウィジェットをそれぞれ指定する。ウィジェットをＴａｂｌｅに配置する。Ｔａｂｌｅに配置後、ＴａｂｌｅをＷｉｎｄｏｗに乗せる。
【００１８】
大きく分別できる二つの技術要素の中の一つ目である（Ａ）仮想心臓の構築について説明する。本システムの仮想心臓の特徴は、心臓という動く臓器をモデリングすることである。ＧＵＩに示される仮想心臓は、三次元描画グラフィックスライブラリＯｐｅｎＧＬを使ったワイヤーフレーム形式にて描画している。ＯｐｅｎＧＬは既存の技術であり、Ｆｒｅｅのソフトウェアである。ＯｐｅｎＧＬでワイヤーフレームを描くためには、大きく分けると、次のような手順が必要である。
第一に、描くもの、本システムでは仮想心臓の三次元座標を取得する必要がある。
第二に、三次元座標を取得してＯｐｅｎＧＬによって仮想空間に座標を映し出し、それらの点間を直線で結んだだけでは仮想心臓は角ばったものになるので、それらの点間をスプライン補間という方法で、補間を行なうことで表面を滑らかに描画する。
第三に、仮想心臓は、仮想指の位置（力覚伝達機構からデータ転送される指の位置データ）によって、仮想指が仮想心臓表面に接触すると、仮想心臓は指の位置と、接触部周辺の心臓の筋肉の硬さの具合によって変形する。そのため点間を五つの要素に分けて質点を配置し、点間で構成される１格子を２５個の質点要素で構成する力学モデルパッチを割り当てる。
心臓の弾性及び形状処理をするためには、筋肉を模擬出来るモデルを作らなければならない。筋肉というものは、多数の細い筋肉繊維が多層構造で構成されている。従来の技術では、このような複雑な筋肉組織を模擬するためには有限要素法という技術を使うのが一般的であるが、これは模擬するものの全体を有限個の要素に分割する計算技術である。例えば５［ｃｍ３］の肉片であったとすれば、それを１０万個の要素に分割して、各要素の挙動を計算する手法である。
しかし、有限要素法では要素の数が多いため計算量が必然的に増加する。計算量が大きいと時間遅れが発生する。触診システムでは、手の接触に対して時間遅れなく力を計算する必要があるため、有限要素法では、時間遅れがあってリアルタイム弾性計算に対応できない。
そのため、本発明では、心筋弾性計算処理に有限要素法は用いず、質量、バネ、ダンパーの機械的要素を用いて低次元力学モデルを設計した。
【００１９】
仮想心臓の構築の手順を大きく分けたときの手順第一である仮想心臓の三次元座標取得について述べる。
心臓は左心室が拡張と収縮を繰り返すことで全身に血液を送っているため、心臓形状はその動きを繰り返して常に動いている。そこで拍動１周期を等時間間隔に分け、それぞれの時の形状ごとに三次元座標を取得する。三次元座標取得の方法として、高速シネＣＴ画像より取得したヒト拍動心臓の断面動画を利用する。断面動画について分かり易くするため図５の模式図により説明する。模式図の（イ）は心臓を上から見た図、（ロ）は正面から見た図である。
図５の（イ）のＡ方向の断面動画より、１心拍の動画より等時間間隔に拡張期５枚、収縮期９枚の画像を抜き出す。画像取得のタイミングを図６示す。この１４枚の断面図の垂直方向座標を、図５の（ロ）のＢ方向の拡張末期と収縮末期の断面図より補間し、心臓表面の三次元座標を取得する。
座標取得の方法につき順を追って説明する。まず高速シネＣＴ動画より取得した１４枚の心臓の断面図より、心室と心房の境界に、心房が入らないよう図７の（イ）のように線を引く。境界線の位置を決定すると、１４枚の断面図すべてに画像の同じ位置に線を引く。次に境界線の引いた後の収縮末期の断面画像より、心室の最下部を通過する境界線との垂直線を図７の（ロ）のように引く。この垂直線の位置を決定すると、境界線と同じ様に１４枚の断面図すべてに画像の同じ位置に線を引く。このように画像の同じ位置に線を引くと、高速シネＣＴの動画の心臓は動いているが、画像を取得する空間の位置は変化していないため、測定する際の空間の基準位置を決定することができる。この２本の線の交点を原点とし、心室と心房の境界線をｘ軸、心室の最下部への垂直線をｚ軸、原点を通りこの２軸に直交する線をｙ軸とする。ここでｙ軸に相当する線は、図を見た時、前後方向になるため図には表示していない。このようにすると、高速シネＣＴ画像による断面図はｘ−ｚ平面による断面図、図５の（ロ）のＢ方向の拡張末期と収縮末期の断面図はｘ−ｙ平面による断面図となる。
このようにして座標軸を決定した後、高速シネＣＴ画像による収縮末期の断面図より、図７の（ハ）のように原点と心室の最下部を通過するｘ軸に平行な線とｚ軸との交点を結び、その線を８等分する点を取得する。これにより求まる９個の点を通過するｘ軸に平行な線を引き、それぞれの平行線と心臓表面との交点より表面のｘ−ｚ平面における位置座標が２個ずつ求まり、計１８個の座標が決定される。この作業を残りの１３枚の断面図においても同様に行う。
次に拡張末期と収縮末期の２枚のｘ−ｙ平面による断面図より、ｘ−ｙ平面の座標を求める。拡張末期と収縮末期の２枚の断面図はｘ−ｙ平面の断面図に対応し、収縮末期の断面図において、高速シネＣＴの切断面との交点と等しくなるようにｘ軸を決定し、ｘ軸と垂直になるように断面図の中心に引いた線をｙ軸として決定する。このように決定したｘ軸とｙ軸を拡張末期の断面図においても同様に決定する。それぞれの断面図における軸が決定されると、図８のように原点からｘ軸の正の方向に線を引き、その線と心臓の表面との接点のｘ、ｙ座標を測定する。このように線と心臓表面との交点のｘ、ｙ座標を測定すると、次は原点を中心に線を２２．５°回転させ、同様に心臓表面との交点の座標を取得する。この作業を繰り返すことにより、１つの断面図につき１６個の心臓表面の座標を取得する。
これにより拡張末期と収縮末期におけるｘ−ｙ平面断面図の表面座標が取得できたが、このままでは収縮末期〜拡張末期、拡張末期〜収縮末期の間の断面図がないため、その間の断面座標が必要となる。そこで高速シネＣＴ断面画像におけるｘ軸との交点２点間の座標の長さは収縮末期〜拡張末期の間にかけて長くなり、拡張末期〜収縮末期の間で短くなっていくため、ｘ−ｙ平面における断面の表面座標は中心を決定し、そこからの長さを測ることで決定していることから、収縮末期と拡張末期のｘ軸との交点の２点の長さをそれぞれ半分にし、収縮末期の長さを０％、拡張末期の長さを１００％として、高速シネＣＴによる１４枚の断面図すべてにおいて式（１−１）により変化の割合を求め、拡張末期と収縮末期の両方の断面図で、中心点から同じ角度で引いた線との交点より求めた座標において、その求めた割合分収縮末期の位置から拡張末期の位置まで表面が移動したとして、式（１−２）のようにして中心点から表面までの長さを求め収縮末期と拡張末期以外のときのｘ−ｙ平面断面の表面座標１６点を求める。これにより高速シネＣＴより取得した１４枚の断面図すべてにｘ−ｙ平面断面の表面座標１６点が求まる。ｘ−ｙ平面断面の表面座標１６点を求めると、中心から０°と１８０°回転させた線との交点より求めた表面座標は、ｘ軸との交点と一致するため、その２点の長さと、高速シネＣＴ断面図より測定したｘ軸との交点２点間の長さが等しくなるように、取得した中心点からの表面までの長さを同じ割合だけ変化させる。これにより表面の形状を変えることなく、高速シネＣＴより取得した座標にあわせることができる。前に高速シネＣＴ画像１４枚においてそれぞれｘ軸から心室最下部までを８等分するようにｘ軸に平行に引いた線と心臓表面との交点を１６点求めておいたので、それを用いてそれぞれの平行線との交点２点の間隔に合わせて同様にｘ−ｙ平面の表面座標１６点のサイズを変更し、それぞれ０°と１８０°のときの表面座標を高速シネＣＴの座標とあわせることで、９つの等間隔に配置したｘ−ｙ平面それぞれに対し心臓表面座標１６点を配置できるため、高速シネＣＴ画像１枚につき１６×９の座標が求まる。これを１４枚作成するので、１４パターンの座標が取得できる
【数１】
2004181176.tif 000003
ここでＲｎ：パターンｎにおける変化率、ｌｃｏ：収縮末期における高速シネＣＴ画像のｘ軸との交点の２点間隔、ｌｅｘ：拡張末期における高速シネＣＴ画像のｘ軸との交点の２点間隔、ｌｎ：パターンｎにおける高速シネＣＴ画像のｘ軸との交点の２点間隔、Ｐｎ：パターンｎの断面図における中心から表面までの長さ、Ｐｃｏ：収縮末期の断面図における中心から表面までの長さ、Ｐｅｘ：拡張末期の断面図における中心から表面までの長さとする。
【００２０】
仮想心臓の構築の手順を大きく分けたときの手順第二であるスプライン補間について述べる。
スプライン補間とは、複数の点を順に繋いでいく際、各々の点の区間に前後の多項式と連続となるような多項式を当てはめることにより、複数の点を滑らかに繋ぐ線の多項式を求める方法である。ここでは各点間を３次の多項式で補間する、最も一般的な３次のスプラインを用いる。部分区間Ｉｉ＝［ｘｉ−１，ｘｉ］について考えると、スプラインＳ（ｘ）は（ｘｉ−１，ｙｉ−１）、（ｘｉ，ｙｉ）を通る３次関数となるので、
【外１】
2004181176.tif 000004
とすると、
【数２】
2004181176.tif 000005
と表せる。式（２−２）を２回積分すると、
【数３】
2004181176.tif 000006
となり、Ｓ（ｘｉ−１）＝ｙｉ−１かつＳ（ｘｉ）＝ｙｉよりｃｉとｄｉが求まり、
【数４】
2004181176.tif 000007
となる。また、式（２−４）を積分すると区間Ｉｉで、
【数５】
2004181176.tif 000008
が得られ、ｘｉについての左方微分係数と右方微分係数は次のようになる。
【数６】
2004181176.tif 000009
【外２】
2004181176.tif 000010
は連続なので、
【数７】
2004181176.tif 000011
とならねばならない。このことから、
【数８】
2004181176.tif 000012
となり、Ｍ１，．．．，Ｍｎに関するｎ−２個の一次方程式が求まる。しかしＭ１，．．．，Ｍｎを決定するにはさらに二つの条件が必要となるため、ここではＭ１＝０、Ｍｎ＝０としてＭ１，．．．，Ｍｎを求める。Ｍ１，．．．，Ｍｎが求まれば、式（２−４）よりスプラインは求まる。
今回は計算の簡略化のため連続する４点よりスプラインを求めることで、間の２点間の補間を行うことを考える。ｘ−ｙ平面上の点は中心座標より同じ角度ずつ回転させたときの座標なので、それぞれの点の中心点との角度の間隔は等しい。そこで４点の座標を図９のように各点の間隔を角度として取り、その角度間隔を１、各々の点の中心座標からの長さをｙｎとして、それぞれ（−１，ｙｉ−１）、（０，ｙｉ）、（１，ｙｉ＋１）、（２，ｙｉ＋２）とすると、部分区間［０，１］でのスプラインは、
【数９】
2004181176.tif 000013
と表される
次にｚ軸方向の点の繋がりを考えると、座標を取得しているｘ−ｙ平面は、ｚ軸方向の間隔と等間隔にとってあるので、連続する点の間隔をｚ軸の値として取ると前回と同様間隔は等しくなり、同様にして各点の間隔を１とする。式（２−９）より曲面の場合を考えると、曲面のスプラインは、４×４の計１６点より、
【数１０】
2004181176.tif 000014
と求めることができる。曲面を形成する４×４の計１６点のｘ、ｙ、ｚ座標ごとに式（２−１０）を当てはめることで、点ｙｉ，ｊ、ｙｉ，ｊ＋１、ｙｉ＋１，ｊ、ｙｉ＋１，ｊ＋１で囲まれた範囲の曲面の方程式を求めることができ、図１０のようにスプラインより求めた滑らかな曲面を通過する点を求めて、点の数を縦横５倍に増やすことで仮想心臓の表面を滑らかできる。またこの場合、ｚ軸方向の最上部の点を含む面を求めるとき、さらにｚ軸の上方向を構成する点が存在しないため、点ｙｉ−１ｊ−１、ｙｉ−１ｊ、ｙｉ−１ｊ＋１、ｙｉ−１ｊ＋２が存在せず、スプラインを求めるのに必要な点を確保することができない。そこでその４点を、点ｙｉｊ−１、ｙｉｊ、ｙｉｊ＋１、ｙｉｊ＋２と大きさが等しいとして曲面を求める。これは最下部についても同様に曲面を求めるための点が足りないという問題が起こるため、同様に足りない点を連続する点と同じ大きさとして求める。
【００２１】
仮想心臓の構築の手順を大きく分けたときの手順第三である力学モデルパッチについて述べる。
力学モデルパッチによる仮想心臓表面の変形描画について述べる。力学モデルについては後に詳しく示すが、仮想心臓の表面は点を格子状に配置することで表現している。そこに力学モデルパッチを図１１のように仮想心臓の表面を構成するうちの９×９の点を指の接触点を中心として当てはめる。このパッチの範囲を力学モデルの対応する範囲としてみなす。次に力学モデルは質量をバネとダンパーで構成される心筋要素５×５を結んだモデルとなっており、質量の位置変化により変形を知ることができる構成となっているが、このままでは力学モデルパッチの点の数と力学モデルの点の数が合わないため、図１２のように力学モデルの質量の間を前述の表面の補間と同じやり方でスプライン補間を用いて９×９の変化量を求める。９×９の変化量が求まれば、力学モデルパッチの対応する表面部分において、指の接触している表面の構成点から表面座標を求める際に決めたｚ軸への方向を力学モデルの表面に対して垂直な方向として、パッチの対応する表面の９×９の点をそれぞれ対応する変化量分移動させて仮想心臓を描画することで表面の変形を行なう。
【００２２】
大きく分別できる二つの技術要素の中の二つ目である（Ｂ）弾性伝達機構構築に係る仮想心臓力学モデルについて述べる。
実際の心臓は心筋によって構成されるが、心筋自体は多数の筋繊維による複雑な多層構造である。従来の技術を使うのならば、有限要素法によって指で押すことによる心臓形状の変化や弾性力計算を行なうことが可能である。しかし、有限要素法では計算時間が課題になりリアルタイム制御に支障をきたすという問題がある。
この問題を克服する心筋の三次元力学高速モデルとして、心筋要素モデルを設ける。心筋要素モデルは質量、バネ、ダンパーによって構成される。心筋を解剖学知見で観察すると、横方向のつながりを形成する筋繊維による層と縦方向のつながりを形成する層が重なって出来る心筋壁が形成されて心筋となっていることに気付く。本発明では、このことを足がかりに心筋を横方向と縦方向で異なった心筋要素モデルとして模擬した。
図１３は心筋モデルの姿を表す図であり、（イ）は横方向の心筋要素モデルを表す図の内で心筋弛緩期のモデルを表す図、（ロ）は横方向の心筋要素モデルを表す図の内で心筋収縮期のモデルを表す図、（ハ）は縦方向の心筋要素モデルを表す図、（ニ）は横方向の心筋要素モデルと縦方向の心筋要素モデルと質量モデルが配設されている心筋モデルの姿を表す図、（ホ）は仮想的な指の力と左心室圧に対して心筋モデルが変形して凹む横からの姿を表す図である。
【００２３】
図１３に示すように力学モデルに用いられる機械的要素はバネとダンパーである。心筋の横方向モデルは、基本的にはバネＫＨ１とダンパーＣＨによる直列モデルとバネＫＨ２による並列モデルで構成される。前述のように力学モデルは縦方向と横方向とで異なったモデルを用いた。そして、心筋収縮期に心筋が硬くなる生理的性質をモデルとして組み込むため、横方向の心筋モデルは拍動周期中にモデルを変化させる。
【００２４】
心筋の拍動周期について図１４のＰＶ線図を用いて説明する。心臓拍動における左心室容積と左心室圧について、横軸に左心室容積、縦軸に左心室圧をプロットして線で結ぶと独特の周回が描かれる。これをＰＶ線図という。仮に拍動の始点をＰＶ線図中のＡ点から開始したとすると、Ａ−Ｂ−Ｃ区間は心筋弛緩期と呼ばれ、Ｃ−Ｄ−Ａ区間は心筋興奮期と呼ばれる。Ａ−Ｂ−Ｃ区間のうちＡ−Ｂ区間は等容弛緩期と呼ばれる。ここでは心筋が血液を駆出後もっとも収縮したところで、Ａ−Ｂの区間では容積は変わらず心室の圧力だけ減少する。Ｂ−Ｃでは血液が流入する区間であり、拡張期と呼ばれる。これは左心房圧のほうが左心室圧よりも大となり左心室への血液流入が起こる期間である。この血圧の上下関係は後に詳細に述べる血液循環のＷｉｎｄｋｅｓｓｅｌモデルによって計算される。Ｃ−Ｄは左心室に溜まった血液を拍出するために徐々に心筋が収縮を始める期間であり、このとき容積変化はない。これは左心室圧よりも大動脈圧のほうが大であるからである。心筋が徐々に収縮することで血液が入った左心室の圧力が上昇を始める。左心室圧が大動脈圧を超えると血液は拍出され、容積が減少する。これがＤ−Ａの区間である。心筋の収縮は、横方向に時間、縦方向に心筋収縮力（心筋エラスタンス）をとった場合、収縮力の時間変化を示すカーブはエラスタンスカーブと呼ばれ、図１５に示される。エラスタンスカーブは心臓力学的に式（３−１）で示される。
【数１１】
2004181176.tif 000015
ここでＥｌｔ（ｔ）は心筋エラスタンスの時間変化［ｍｍＨｇ／ｍｌ］、ＰＬＶ（ｔ）は左心室圧の時間変化［ｍｍＨｇ］、ＶＬＶ（ｔ）は左心室容積の時間変化［ｍｌ］、Ｖ０は自然長左室容積［ｍｌ］（左心室に血液が入っていないときの容積）である。図１５のエラスタンスカーブは犬の拍動心臓から得た実データである。
【００２５】
パラメータの同定手段について説明する。横方向は前述したように心筋の層であるため、解体された心臓ではその特性を得ることは不可能である。そのため横方向パラメータは拍動心臓の血流データから同定を行なう。縦方向パラメータは摘出心臓を切り開き、力センサを用いて弾性データの取得を行なう。
説明内容を分かり易くするため、測定方法について少し説明する。犬の拍動心臓から、心臓に拘束条件を加えた上で血液データを測定し、血液データを反映する心臓モデルをパラメータ解析用に設けて同定を行なった。拘束条件は、図１６に示される心筋変位測定器を利用する。測定器はアルミフレーム２０とポテンショメータ（角度センサ）２２、力センサ２４、スライド式固定用板２６で構成される。アルミフレーム２０は、角度固定用ネジ２８でスライド式固定用板２６に固定される。ポテンショメータ２２で心筋への押しこみ変位を測定しつつ、力センサ２４で弾性を得ることは、心臓表面と力センサ表面が滑って弾性データを得ることが不可能であるので、血液データの解析による方法のみに焦点を絞った。縦方向の弾性データは豚の摘出心臓から力センサ２４を用いて測定した。横方向は犬、縦方向は豚としたのは、心臓外科医による動物実験の途中に行うなどの試験の都合によるもので、データに一貫性はないが、心筋触診訓練システムの初期値としては十分である。
横方向力学的パラメータ、バネＫＨ１とダンパーＣＨによる直列モデルとバネＫＨ２の同定のための血圧データの取得について説明する。測定の方法は、心筋変位測定器で心筋の基部、中央部、心尖部をそれぞれ固定する。それぞれの固定場所において、変位測定器の変位を狭めることで心筋を摘むように拘束する。拘束後、レールに取り付けられたネジ２８によって固定する。固定された部位は、心筋として稼動しないと考えると、拘束部位と拘束変位によって血流データに現れる影響から、拘束部位による心筋の弾性特性を割り出せるという考えである。
【００２６】
横方向パラメータの同定手段について説明する。図１７は心臓を円筒と機械的モデルを使って模擬したもので、円筒モデルと呼ぶことにする。右側の円筒３０が左心室容積を示し、左側の機械的機構３２が心筋モデルを示す。この心筋モデルは、図１３に示す力学モデルの横方向モデルの弛緩期版と同形である。ここで始めに、円筒の底面積をＳ０とする。図１４のＰＶ線図は犬の実験から得たデータによって書かれたものである。本発明では拍動周期を４つの区間に分け、等容弛緩期、拡張期、等容収縮期、駆出期としている。単純に弾性のみを解析できる区間は、心筋が弛緩状態にある等容性弛緩期と拡張期であると考えることができる。これは収縮期に入ると、心筋は自発的に硬くなり、血流データでは弾性パラメータを導出できない。
等容弛緩期では、容積は一定であるので、ｘ１は不動である。この期間において、左心室の圧力は低下する。圧力は左心室壁に内圧として働き、これは左心室圧と円筒の底面積との積で求めることができる。これは式（３−２）となる。
【数１２】
2004181176.tif 000016
但し、ｒは容積の円筒モデル上底面積であるため、半径ｒ＝０．０４［ｍ］とした。
この力は心筋の力、すなわち心筋弾性でつりあいが取れているといえるので、式（３−４）が成り立つ。ｘ１は不動であるので、ここではＫ２は関与しない。
【数１３】
2004181176.tif 000017
また、心筋モデルのほうでは、式（３−５）の釣り合い式が成立する。
【数１４】
2004181176.tif 000018
式（３−５）の微分方程式を解くと式（３−６）が導かれる。
【数１５】
2004181176.tif 000019
式（３−４）、（３−６）から式（３−７）によって左心室圧の変化を数式で表すことができる。
【数１６】
2004181176.tif 000020
式（３−６）を用いて、等容弛緩期の左心室圧データを近似すると、バネＫ１とダンパーＣの比を得ることができる。心臓の基部を固定した室圧データでは、図１８のように近似できる。基部固定ではバネＫ１とダンパーＣの比を得ることが出来た。
【００２７】
等容弛緩期においてバネＫ１とダンパーＣの比を求めたので次に、拡張期のデータ解析を行なう。拡張期では、左心室圧と左心室容積は共に動的である。このときの心筋モデルの弾性力は式（３−８）に示される。
【数１７】
2004181176.tif 000021
ここで実データからｘ１の変動を求める必要がある。円筒モデルにおいて、容積の変化は式（３−９）で示される。
【数１８】
2004181176.tif 000022
ここでｘ１を２次関数で定義するために
【数１９】
2004181176.tif 000023
とすると、ｘ１（ｔ）は図１９に示されるように近似される。図１９は左心室容積変化の実データを２次関数で近似したグラフである。
心筋モデルの直列要素の釣り合い式は
【数２０】
2004181176.tif 000024
式（３−１１）をラプラス変換すると
【数２１】
2004181176.tif 000025
Ｘ２について整理する。
【数２２】
2004181176.tif 000026
また、式（３−１３）をラプラス変換して整理する。
【数２３】
2004181176.tif 000027
Ｘ２を消去するために式（３−１３）を式（３−１４）へ代入する
【数２４】
2004181176.tif 000028
式（３−１６）を逆ラプラス変換すると
【数２５】
2004181176.tif 000029
実験データの力を
【外３】
2004181176.tif 000030
とおき、ｔ、Ｋ１、Ｋ２を入力変数として得られる力を
【外４】
2004181176.tif 000031
とし式（３−１６）との最小二乗誤差をとる。ここで、
【数２６】
2004181176.tif 000032
（５−１８）を最小にするＫ１、Ｋ２を求める。
【数２７】
2004181176.tif 000033
これらより
【数２８】
2004181176.tif 000034
とＫ１、Ｋ２に関する連立方程式を得る。上記した式において、積分時間域ｔＢは拡張期開始時刻、ｔＣは拡張期終了時刻を示す。ここまでで、基部を拘束した実験データで得られたＫ１、Ｋ２の値を求めた。式（３−７）のＫ１とＣとの比が分っているので、Ｃの値を求めることができる。
ここまでの流れで求めたパラメータは基部を抑えることで、中央部のバネと心尖部のバネが伸びて得られた結果であることが重要である。基部、中央部、心尖部を抑えることで得られたバネ定数をそれぞれＫＢ、ＫＭ、ＫＡとすると、バネ定数の逆数で表すことができるコンプライアンスの概念に基づいて、各部位のコンプライアンスをＣＢ、ＣＭ、ＣＡとすると、それぞれ
【数２９】
2004181176.tif 000035
基部計測で得られた値の逆数はＣ１、中央部計測で得られた値の逆数をＣ２、心尖部計測データで得られた値の逆数をＣ３とすると式（３−２６）の関係となる。
【数３０】
2004181176.tif 000036
心筋の各部位による硬さの順序は、解剖学的にも、また医師の経験からも基部、中央部、心尖部であるため、ＣＡ＞ＣＭ＞ＣＢ、Ｃ１＞Ｃ２＞Ｃ３の大小関係とする。以上より各パラメータを同定することができる。
【００２８】
縦方向パラメータ同定手段について説明する。
心筋の縦方向モデルはバネＫＶとダンパ−ＣＶとの並列要素である。ここでＶはＶｅｒｔｉｃａｌの意を示す。これら二つの機械的パラメータを同定する方法について述べる。
パラメータの同定には、まず弾性係数ＫＶは心筋表面に対して図１６の心筋変位測定器の力センサを押し付けて得られる弾性データとポテンショメータで得られる心筋変位データとの関係から、変位・弾性力すなわち弾性係数バネＫＶ［Ｎ／ｍ］を求めることができる。この実験を基部、基部と中央部の中間点、中央部、中央部と心尖部の中間点、心尖部において行った。測定点は５点であるが、仮想心臓は３５点の力点がある。測定した５点間を滑らかに補間するように残りの３１点を決定した。図２０は仮想心臓の力学モデルに与えた縦方向弾性係数を示すグラフであり、このグラフは心筋変位測定器の力センサを押し付けて得られる弾性データとポテンショメータで得られる心筋変位データとの関係から、変位・弾性力すなわち弾性係数を求め作成した。横軸に心筋部位（左から基部、中央部、心尖部）、縦軸に弾性係数を示す。
粘性係数Ｃｖ［Ｎ／ｓ／ｍ］はスティックを使って、心筋表面を押し込み、瞬間的にスティックを離し、心筋表面の復元時間を測定する。この測定方法をＲａｐｉｄＲｅｌｅａｓｅ方式とよぶ。弾性係数は既に求められた、縦方向モデルの運動方程式は式（３−２７）で表すことができる。
【数３１】
2004181176.tif 000037
質量ｍを１０［ｇ］と仮定して、式（３−２７）をシミュレーション解析し、復元時間が一致するようにＣｖを調整することで同定する。図２１は仮想心臓の力学モデルに与えた縦方向粘性係数を示しており、横軸に心筋部位（左から基部、中央部、心尖部）、縦軸に粘性係数を示している。
【００２９】
収縮要素について説明する。
ここまでで、拍動周期中の弛緩期における力学モデルのパラメータの同定について述べた。縦方向は拍動周期中すべての時期において同一パラメータを使用してもよいが、心筋収縮期に心筋が硬くなる生理的性質の影響が横方向は大きいため、横方向の心筋モデルは拍動周期中、前記のように弛緩期と収縮期で違うパラメータを使用する必要がある。
収縮期の横方向において、Ｋ２とＣに置き換わる収縮要素について説明する。収縮要素は、心筋が収縮時に硬くなることを模擬するために用いる時間変数式バネである。心筋が収縮する強さは一般に心筋エラスタンスと呼ばれる。また、収縮期における収縮強さの時間変化をエラスタンスカーブとよび、既に数式は式（３−１）、グラフは図１５に示した。収縮要素は、エラスタンス依存型バネということで式（３−２８）に示される。
【数３２】
2004181176.tif 000038
Ｓｎは指が接触した面積であり、Ｋｉｎｉｔは場所による値である。Ｋ２の部分に置き換わった収縮要素であれば、Ｋｉｎｉｔ＝Ｋ２であり、Ｃの変わりに置き換わった収縮要素であれば、ここではテスト的に１としている。それはＳｎは触れていないときはＫＥｌｔとなるので、モデルとして成立しない。このための補正値である。補正値＝１は仮の値であり、熟練した医師により調整されるべきである。
【００３０】
病変部と正常部との違い、訓練者の手入力による差別化について述べる。例えば正常心筋部は拍動周期において、上述した力学モデルに同定したパラメータを利用しているとする。すべての係数に対してｋαを乗じることにする。ここでｋαは補正係数である。ｋα＝１であれば通常の同定したままのパラメータである。つまり正常心筋弾性力および収縮強さを模擬していることになる。もしもｋα＜１であれば、何らかの原因で収縮力が弱まった疾患心筋を模擬したことになる。訓練者の手入力による病変部と正常部との差別化はこのｋαの入力になる。現在のプログラムでは、縦方向及び横方向のすべてのパラメータにｋαを乗じており、初期値はｋα＝１である。現在このｋαは０．０１刻みで変化させることができる。補正係数ｋαの影響、訓練者の手入力のｋαを視覚的に捕らえることができる。これは本発明のインターフェイス上で捉えることができる。仮想心臓に色をつけて識別する方法である。色のつけ方は実際の病気に倣っている。心筋梗塞が発生直後は、心筋は内出血状態であり青紫色になる。このときｋα＝０．８、０．９程度、その後心筋は徐々に繊維化を初め、白色に近い状態に陥っていく。
【００３１】
血液循環モデルについて述べる。
本発明に用いる仮想心臓には血液循環モデルを組み込む。このモデルはＷｉｎｄｋｅｓｓｅｌモデルといい、人体の血液循環を電気回路で模擬したものである。血流を電流、血管抵抗を電気抵抗、血管柔らかさをコンデンサ、心筋収縮力を可変コンデンサでそれぞれ模擬される。計算されるパラメータは、左心房圧Ｐｌａ、左心室圧Ｐｌｖ、左心室容積Ｖｌｖ、左心室流入量Ｑｌ、左心室流出量Ｑａｏ、体循環血流量Ｑｓ、気管支動脈血流量Ｑｂ、冠動脈流量Ｑｃｏ、大動脈圧Ｐａｏ、右心房圧Ｐｒａ、右心室圧Ｐｒｖ、肺動脈圧Ｐｐａ、右室流入量Ｑｒ、右室流出量Ｑｐａ、肺循環血流量Ｑｐである。
【００３２】
Ｗｉｎｄｋｅｓｓｅｌモデルに使用される血管抵抗値は、過去に本願発明者らによってシミュレーション推定された値を用いた。例えば、体循環を並列回路とみなし、全体の抵抗値を１として，脳であれば全体血流量の数％程度などとして，血液が流れる割合から血管抵抗値を決定する方法である。
【００３３】
心筋エラスタンスについて述べる。
Ｗｉｎｄｋｅｓｓｅｌモデルは、心臓の収縮強さである心筋エラスタンスを与えることでシミュレーションは実行される。無論、実際の心臓においても心筋が収縮拡張を繰り返すことで血液は循環する。シミュレーションの起動ソースとなる心筋エラスタンスは心臓力学において式（３−１）において示される。
【数３３】
2004181176.tif 000039
ここでＥｌｔ（ｔ）は心筋エラスタンスの時間変化、Ｐｌｖ（ｔ）は左心室圧の時間変化、Ｖｌｖ（ｔ）は左心室容積の時間変化、ＶＬ０は自然無血状態の左心室容積である。これらＰｌｖ（ｔ）、Ｖｌｖ（ｔ）の値には犬の心臓から得た左心室圧と左心室容積のデータを入力し、犬の心筋エラスタンスの時間変化を得ることができた。
【００３４】
以下に、シミュレーションするためのＷｉｎｄｋｅｓｓｅｌモデルの計算式を示す。
左心室容積の計算式、左心室内での血流量に関する式：［ｍｌ／ｓｅｃ］。
【数３４】
2004181176.tif 000040
大動脈圧の計算式、大動脈での血流量に関する式：［ｍｌ／ｓｅｃ］
【数３５】
2004181176.tif 000041
中心静脈での血流量に関する式：［ｍｌ／ｓｅｃ］
【数３６】
2004181176.tif 000042
左心室圧：［ｍｍＨｇ］
【数３７】
2004181176.tif 000043
ＶＬ０：無血状態自然長容積、左心室容積の初期値
左室流入弁抵抗での流量：
【３８】
2004181176.tif 000044
（左心室流入弁の開閉条件を満たす時（Ｐｌａ＞Ｐｌｖ）に、Ｑｌは流れる。）
左室流出弁での流量（大動脈流量）：［ｍｌ／ｓｅｃ］
【数３９】
2004181176.tif 000045
（左心室流出弁の開閉条件を満たす時（Ｐｌｖ＞Ｐａｏ）に、Ｑａｏは流れる。）
体循環血流量：［ｍｌ／ｓｅｃ］
【数４０】
2004181176.tif 000046
気管支動脈流量：［ｍｌ／ｓｅｃ］
【数４１】
2004181176.tif 000047
冠動脈流量：［ｍｌ／ｓｅｃ］
【数４２】
2004181176.tif 000048
右心室内での血流量に関する式：［ｍｌ／ｓｅｃ］
【数４３】
2004181176.tif 000049
肺動脈での血流量に関する式：［ｍｌ／ｓｅｃ］
【数４４】
2004181176.tif 000050
肺静脈での血流量に関する式：［ｍｌ／ｓｅｃ］
【数４５】
2004181176.tif 000051
右心室圧：［ｍｍＨｇ］
【数４６】
2004181176.tif 000052
ＶＲ０：右心室容積の初期値
右室流入弁抵抗での流量：
【数４７】
2004181176.tif 000053
（右心室流入弁の開閉条件を満たす時（Ｐｒａ＞Ｐｒｖ）に、Ｑｒは流れる。）
右室流出弁での流量：［ｍｌ／ｓｅｃ］
【数４８】
2004181176.tif 000054
（左心室流出弁の開閉条件を満たす時（Ｐｒｖ＞Ｐｐａ）に、Ｑｐａは流れる。）
肺循環血流量：［ｍｌ／ｓｅｃ］
【数４９】
2004181176.tif 000055
但し、ここで右心室に関する測定データはなく、右心室エラスタンスは左心室エラスタンスの２割とした。
【００３５】
シミュレーションにおけるエラスタンスについて述べる。
心筋エラスタンスを横軸に１拍動周期時間、縦軸にエラスタンスをとるとカーブが示される。これをエラスタンスカーブと呼ぶ。本発明では、心筋収縮期のみ心筋エラスタンスの実データから与えている。心筋弛緩期のエラスタンスはバネとダンパーの特性からの計算結果から導出する形としている。このことの利点は、収縮期においては本発明では、収縮期において心筋が硬くなる性質を模擬するためにエラスタンス依存型の時間変化的バネ（収縮要素）として与えており、これは心筋弛緩期には現れない。心筋弛緩期、等容性弛緩期では、ただ心筋の弾性に従って左心室圧が下がり、拡張期には、左心室圧と容積が拡張する。これらの特性は、仮想心臓力学モデルで述べて機械的パラメータに依存する。
心筋収縮期における左心室圧の計算式は式（３−５）であるが、等容弛緩期において、左心室圧の計算式は、
【数５０】
2004181176.tif 000056
但し、ここで使用するＫとＣの値は、基部、中央部、心尖部の同定された弾性係数Ｋ１の平均値であり、Ｃは基部、中央部、心尖部の同定された粘性係数Ｃである。またＰＬＶ（０）は等容性弛緩期に入った直後の左心室圧である。
拡張期における左心室圧の変化は、式（３−１９）で示される。この式は左心室圧と円筒の底面積Ｓ０との積で得られる力の変化から得られる。ＶＬＶ（０）は拡張期に入った直後の左心室容積である。
【数５１】
2004181176.tif 000057
【００３６】
前述の説明に関連する事項として、左心房圧Ｐｌａについて、ＰＶ線図について、心拍数ＨＲについて、シミュレーション刻み幅について説明する。
【００３７】
心筋触診を要する手術は、心筋のみの弾性力を触診する必要があるため、左心室圧を０［ｍｍＨｇ］にするために人工心肺を使って血液を外部循環させる。このことを模擬するために、計算中では式（３−１３）を省いて、Ｐｌａを固定値にする。現システムにおけるＰｌａの変動範囲は０〜１５［ｍｍＨｇ］である。Ｐｌａは触診システムのグラフィカルユーザインターフェイス上から変化させることが可能で、左心室圧を０［ｍｍＨｇ］に定着させることが出来ている。
【００３８】
心臓拍動における左心室容積と左心室圧について、横軸に左心室容積、縦軸に左心室圧をプロットして線で結ぶと独特の周回が描かれる。これをＰＶ線図という。触診システムでは、上記のＷｉｎｄｋｅｓｓｅｌモデルで計算された左心室容積と左心室圧のデータをグラフにプロットしたＰＶ線図を実行中に確認することができる。但し、Ｐｌａを０［ｍｍＨｇ］とした場合、ＰＶ線図には点が描かれる。
【００３９】
触診システムのインターフェイスでは、左心房圧の他に心拍数ＨＲを設定することができる。心拍数を変更すると、仮想心臓の拍動周期が短くなり、エラスタンスカーブの間隔が短くなる。仮想心臓の１心拍は１４枚の画像で構成される。また拍動周期の基準は実験で用いた犬の拍動心臓とし、犬の自然心拍数１４０［ｂｐｍ］であったことを用いると、１回心拍時間は０．４２秒である。よって一枚の画像を切り替えるタイミングは０．４２／１４である。これにＨＲ／１４０を乗じて得られる式（３−２０）がＨＲによる画像切り替え時間である。
【数５２】
2004181176.tif 000058
【００４０】
触診システムにおける仮想心臓処理用ＰＣ１０は、Ｒｅｄｈａｔ−Ｌｉｎｕｘによるタイマー機能により、１秒間に１００回のループを繰り返す。さらに１回のループの中でＷｉｎｄｋｅｓｓｅｌモデルの式（３−２９）〜（３−４７）を１０回繰り返す。式中に用いられる刻み幅ｄｔはｄｔ＝０．００１である。
【００４１】
心筋触診訓練システムにおいて、訓練者に力覚情報を与える部分である力覚伝達機器について説明する。
図２２に力覚伝達機器を示す。図２２において、（イ）は力覚伝達機器の模式図、（ロ）は関係を説明するための仮想心臓の模式図である。
駆動部分はＤＣサーボモータを利用している。力覚伝達機器を構成する各機器はＤＣサーボモータ４０、角度検出にはポテンショメータ４２、駆動力伝達はプーリー４４とタイミングベルト４６を用いた。図中、４８はアタッチメント、４９はギアホイール、５０はアルミフレームを示している。ＤＣサーボモータ制御にはハイブリッド制御を用いる。ハイブリッド制御は角度とトルクを同時制御するために用いられ、角度θｒは仮想空間上で計算される仮想指の変位、トルクはトルク制御のための電流制御であり、目標電流値ｉは臓器モデルの反力に相応した値である。ＤＣモータの運動方程式を式（３−２７）、ハイブリッド制御の制御則を式（３−４９）、（３−５０）に示す。指令値は電流である。
【数５３】
2004181176.tif 000059
ここで、ｚはポテンショメータ検出角度のノイズ補償のためのソフトウェアフィルタ式である。ｚはフィルタリング後の角度となり、制御則計算にはｚを用いる。Ｋθ、Ｋω、Ｋｆはフィードバックゲインであり、Ｔｐはソフトウェアフィルタの時定数である。書くフィードバックゲインは式（３−４９）、（３−５０）、（３−５１）から得られる微分方程式をラプラス変換し、極配置法を用いて決定する。フィードバックゲインは、
【数５４】
2004181176.tif 000060
極配置はｓ１，ｓ２，ｓ３＝−５５、ｓ４＝−１００としてＴｐ＝６．０×０．００１、Ｋθ＝０．３１８、Ｋω＝０．０１８、Ｋｆ＝１００とした。ここでＴｐ＝６．０×０．００１は、制御周波数が１０００［Ｈｚ］であるためソフトウェアフィルタ効率は６割である。この数値は実験による調整で得た適正値であり、およそ６割の効率とすればスムーズな応答が得られている。
力覚伝達機器の操作により仮想指の位置を決定するが、三次元空間上の仮想指の移動方向と、力覚伝達機器操作の対応を図２２に示す。図２２中、力覚伝達機器図の矢印と仮想心臓図の矢印は同意である。矢印▲１▼は指５２の仮想心臓５４表面を押す動きであり、（イ）の矢印▲１▼▲２▼をそれぞれ中央方向に指を移動することで仮想指は（ロ）の矢印▲１▼▲２▼の仮想心臓中心方向に移動する。この方向において仮想心臓の弾力を取得する。但し、指は心臓によって押し戻されるので双方向に移動する。▲２▼は指を心臓表面に対して上下方向の移動を行なう。▲３▼は手首の動きに合わせて２本指をｘｚ平面に対して回転させる。▲４▼は手首の動きに合わせて２本指を同ｘｙ平面に対して方向回転させる。▲５▼は手首をｙ方向、奥行き方向に移動させる。
なお、図２２の（ロ）において、５６は力学モデルパッチを示している。
力覚伝達機器の機構から、弾性伝達駆動はＤＣサーボモータで行なっており、ＤＣサーボモータは指の動き（矢印▲１▼）にのみ対応しているため、現在のところ仮想心臓弾性は指先でのみ得ることができる。
力覚伝達機器のリアルタイム制御はμ秒単位のリアルタイム処理が補償されているＲＴ−ＬｉｎｕｘＯＳを使用している。ＲＴとはＲｅａｌＴｉｍｅを意味する。ここで仮想心臓処理用ＰＣ１０とのデータ通信の折り合いも含めてリアルタイム処理プログラムのフローチャートを図２３に示す。ＲＴＬｉｎｕｘは大きくユーザ空間とモジュール空間とで制御を行なう。ユーザ空間では、画面表示のｐｒｉｎｔｆ関数などの使用は許可されているが、力覚伝達機器の制御であるＡＤ／ＤＡボードの制御はモジュール空間で行なわれる。モジュール空間では基本的に整数計算のみ許されるが、それはリアルタイムタスク内の規則である。モジュール空間は、Ｌｉｎｕｘ処理内でモジュールプログラムとして組み込まれる。Ｋ−ｔｅｒｍ上において、ｍａｋｅ  ｍｏｄｕｌｅ．ｃでｍｏｄｕｌｅ．ｏモジュールプログラムができる。これをｉｎｓｍｏｄ  ｍｏｄｕｌｅ．ｏでＰＣに組み込む。あとはユーザ空間での操作となる。
ユーザ空間では、プログラム構築中に制御時間、制御サンプリング周波数を設定する。プログラムを実行すると始めにデータ通信ボード使用のため、ボードがオープン状態となる。このとき、モジュール空間のリアルタイムタスク中に実数計算が許されないので、リアルタイムタスクに入る前にフィードバックゲインの実数計算を行ない、リアルタイムタスク中に利用できる整数型に変換される（スケーリング）。その後すぐにモジュール空間にタスク指令が出て、リアルタイム制御に入る。このときからモジュール空間はリアルタイムタスクを繰り返す。モジュール空間とユーザ空間とのやり取りは可能である。データ通信ボードはユーザ空間で取り扱う。ユーザ空間でスタートタスク（リアルタイム開始指令）をモジュール空間へ発し、その後、処理ループに入り、１．仮想心臓処理ＰＣから仮想指の位置をデータ通信する。２．仮想指の情報をモジュール空間へ送信する。３．モジュール空間から送信される訓練者の指の位置を仮想心臓処理ＰＣにデータ通信する。４．ｐｒｉｎｔｆ関数により画面表示を行なう。
一方モジュール空間では、ユーザ空間からスタートタスクを受け取るとリアルタイムタスクの繰り返しループに入る。１．ユーザ空間からの仮想指の位置と弾性情報を受信する。２．仮想指の位置と仮想心臓弾性データを目標値として制御則の計算をする。３．ＤＡボードからモータ駆動電流に相当するデジタル信号（電圧−１０〜＋１０Ｖの範囲、５Ｖで１０２４）を出力する。但し、ＤＡ出力されたアナログ信号（電圧−１０〜＋１０Ｖ）はモータドライバを介してＤＣモータを制御する。
制御終了処理は、Ｋ−ｔｅｒｍ上でＣｔｒｌ＋Ｃで強引に終了するか、当初設定した制御終了カウントと現在の制御タイムカウントとを比較して終了するかを決定する。
【００４２】
【発明の効果】
以上説明したように、本心筋触診訓練システムは、高速シネＣＴ動画像からのデータによって構築した仮想心臓によって、実際の心臓の拍動に近い画面表示が出来るとともに、筋肉を模擬する力学モデルを作ることによって、有限要素法を用いることなく心筋触診における指の押えに対する凹みの変動量および反作用の数値計算が可能になって、アクチュエーターのアタッチメントに指が触れて指が心臓を押したことによって、リアルタイムに、心臓の表面の凹みが実際の心臓触診時の心臓表面の凹みに近い姿で仮想心臓画像の上に加わって描画され、また実際の心臓触診に近い感触がその反作用として指に戻されることで心臓触診の臨場感を体感できる。
また変動量および反作用の数値計算にはパラメータを同定し、手入力によって数値変更により、健常者や心筋梗塞患者の心筋の個人差に合わせて、その出力の変更が出来るようにしている。
上記のように、手入力によって、健常者や心筋梗塞患者の心筋の個人差に合わせて、その出力の変更が出来るようにした有用で心臓外科医待望の心筋触診訓練システムを提供することができる。
また、将来、手入力によるとはいえ、実際の心筋梗塞患者を開胸することのない、近似触診の実現につながるものであり、ＭＲＩなどの医療機器の高速化が実現して、手によって入力する部分のデータ入力がＭＲＩなどから直接に入力可能になれば、心筋梗塞患者などのリアルタイムの近似触診の実現につながるものである。
【００４３】
【図面の簡単な説明】
【図１】本発明の実施形態の心筋触診訓練システムの構成を示す模式図。
【図２】本発明の実施形態の心筋触診訓練システムを利用するためのフローチャート。
【図３】本発明の実施形態の心筋触診訓練システムを利用するとき、画面に表示されるＧＵＩの様子を示す模式図。
【図４】本発明の実施形態の心筋触診訓練システムにおけるＧＵＩ構築のフローチャート。
【図５】本発明の実施形態の仮想心臓構築における三次元座標の取得の説明を分かり易くするために作った心臓の模式図であり、（イ）は心臓を上から見た模式図、（ロ）は心臓を正面から見た模式図。
【図６】本発明の実施形態の仮想心臓構築における三次元座標の取得の説明を分かり易くするために作った高速シネＣＴによる断面図の画像取得のタイミングを示す図
【図７】本発明の実施形態の仮想心臓構築における三次元座標の取得の説明を分かり易くするために作った心臓を正面から見た模式図であり、（イ）は心室と心房の境界に、心房が入らないよう境界線を引くところの模式図、（ロ）は心室の最下部を通過する境界線との垂直線を引くところの模式図、（ハ）は原点と心室の最下部を通過するｘ軸に平行な線とｚ軸との交点を結び、その線を８等分する点を取得し、これにより求まる９個の点を通過するｘ軸に平行な線を引き、それぞれの平行線と心臓表面との交点より表面のｘ−ｚ平面における位置座標が２個ずつ求まり、計１８の座標が決定されるところの模式図。
【図８】本発明の実施形態の仮想心臓構築における三次元座標の取得の説明を分かり易くするために作った模式図。
【図９】本発明の実施形態の仮想心臓構築におけるスプライン補間についての説明を分かり易くするために作った模式図。
【図１０】本発明の実施形態の仮想心臓構築におけるスプライン補間についての説明を分かり易くするために作った模式図。
【図１１】本発明の実施形態の仮想心臓構築における仮想心臓表面の力学モデルパッチの配置を示す模式図。
【図１２】本発明の実施形態の仮想心臓構築におけるスプラインによる力学モデルの変位の補間を示す模式図。
【図１３】本発明の実施形態の心筋モデルの姿を表す図を示し、（イ）は横方向の心筋要素モデルを表す図の内で心筋弛緩期のモデルを表す図、（ロ）は横方向の心筋要素モデルを表す図の内で心筋収縮期のモデルを表す図、（ハ）は縦方向の心筋要素モデルを表す図、（ニ）は横方向の心筋要素モデルと縦方向の心筋要素モデルと質量モデルが配設されている心筋モデルの姿を表す図、（ホ）は仮想的な指の力と左心室圧に対して心筋モデルが変形して凹む横からの姿を表す図。
【図１４】本発明の実施形態の心筋触診訓練システムの発明に利用した犬の実験から得たデータによって書かれたＰＶ線図。
【図１５】本発明の実施形態の心筋触診訓練システムの発明に利用した犬の拍動心臓から得た実データから作られたエラスタンスカーブのグラフ。
【図１６】本発明の実施形態の心筋触診訓練システムの発明に利用したもので、犬の拍動心臓に拘束条件を加えた上で血液データを測定するときに利用した心筋変位測定器の模式図。
【図１７】本発明の実施形態の心臓を円筒と機械的モデルを使って模擬した円筒モデルの姿を表す図。
【図１８】本発明の実施形態の心筋触診訓練システムの発明に利用したもので、犬の拍動心臓に心筋変位測定器の力センサを押し付けて得られる弾性データから作られたグラフ。
【図１９】本発明の実施形態の心筋触診訓練システムの発明に利用したもので、左心室容積変化の実データを２次関数で近似したグラフ。
【図２０】本発明の実施形態の心筋触診訓練システムの発明に利用したもので、心筋変位測定器の力センサを押し付けて得られる弾性データとポテンショメータで得られる心筋変位データとの関係から、変位・弾性力すなわち弾性係数を求め作成したグラフ。
【図２１】本発明の実施形態の心筋触診訓練システムの発明に利用したもので、仮想心臓の力学モデルに与えた縦方向粘性係数のグラフ。
【図２２】本発明の実施形態の心筋触診訓練システムにおける訓練者に力覚情報を与える部分である力覚伝達機器における機器の操作と仮想心臓との関係を説明するための模式図であり、（イ）は力覚伝達機器の模式図、（ロ）は関係を説明するための仮想心臓の模式図。
【図２３】本発明の実施形態の心筋触診訓練システムにおける力覚伝達機器のリアルタイム処理プログラムのフローチャート。
【符号の説明】
８    データ通信ボード
１０  仮想心臓処理用ＰＣ
１２  力覚伝達機器制御用ＰＣ
１４  訓練者
１６  力覚伝達機器
１８  モータドライバ
２０  アルミフレーム
２２  ポテンショメータ（角度センサ）
２４  力センサ
２６  スライド式固定用板
２８  角度固定用ネジ
４０  ＤＣサーボモータ
４２  ポテンショメータ
４４  プーリー
４６  タイミングベルト
４８  アタッチメント
４９  ギアホイール
５０  アルミフレーム
【図１】
2004181176.tif 000061
【図２】
2004181176.tif 000062
【図３】
2004181176.tif 000063
【図４】
2004181176.tif 000064
【図５】
2004181176.tif 000065
【図６】
2004181176.tif 000066
【図７】
2004181176.tif 000067
【図８】
2004181176.tif 000068
【図９】
2004181176.tif 000069
【図１０】
2004181176.tif 000070
【図１１】
2004181176.tif 000071
【図１２】
2004181176.tif 000072
【図１３】
2004181176.tif 000073
【図１４】
2004181176.tif 000074
【図１５】
2004181176.tif 000075
【図１６】
2004181176.tif 000076
【図１７】
2004181176.tif 000077
【図１８】
2004181176.tif 000078
【図１９】
2004181176.tif 000079
【図２０】
2004181176.tif 000080
【図２１】
2004181176.tif 000081
【図２２】
2004181176.tif 000082
【図２３】
2004181176.tif 000083