JP6760603B2 - 流脈線可視化装置、流脈線可視化方法、および流脈線可視化プログラム - Google Patents
流脈線可視化装置、流脈線可視化方法、および流脈線可視化プログラム Download PDFInfo
- Publication number
- JP6760603B2 JP6760603B2 JP2017075671A JP2017075671A JP6760603B2 JP 6760603 B2 JP6760603 B2 JP 6760603B2 JP 2017075671 A JP2017075671 A JP 2017075671A JP 2017075671 A JP2017075671 A JP 2017075671A JP 6760603 B2 JP6760603 B2 JP 6760603B2
- Authority
- JP
- Japan
- Prior art keywords
- time
- fluid
- vector
- velocity
- analysis
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 210000003462 vein Anatomy 0.000 title claims description 165
- 238000012800 visualization Methods 0.000 title claims description 64
- 238000007794 visualization technique Methods 0.000 title claims description 3
- 238000004364 calculation method Methods 0.000 claims description 362
- 239000012530 fluid Substances 0.000 claims description 169
- 238000004458 analytical method Methods 0.000 claims description 137
- 238000012545 processing Methods 0.000 claims description 95
- 239000002245 particle Substances 0.000 claims description 63
- 238000004088 simulation Methods 0.000 claims description 63
- 230000001133 acceleration Effects 0.000 claims description 26
- 230000008859 change Effects 0.000 claims description 19
- 238000012937 correction Methods 0.000 claims description 11
- 239000013598 vector Substances 0.000 description 280
- 238000000034 method Methods 0.000 description 164
- 230000008569 process Effects 0.000 description 90
- 210000004165 myocardium Anatomy 0.000 description 64
- 230000017531 blood circulation Effects 0.000 description 36
- 230000006870 function Effects 0.000 description 27
- 238000010586 diagram Methods 0.000 description 26
- 230000002107 myocardial effect Effects 0.000 description 24
- 230000000694 effects Effects 0.000 description 8
- 230000006399 behavior Effects 0.000 description 7
- 238000002474 experimental method Methods 0.000 description 7
- 210000000709 aorta Anatomy 0.000 description 6
- 230000003287 optical effect Effects 0.000 description 6
- 230000008602 contraction Effects 0.000 description 5
- 210000004204 blood vessel Anatomy 0.000 description 4
- 238000011156 evaluation Methods 0.000 description 4
- 230000006872 improvement Effects 0.000 description 4
- 101100285518 Drosophila melanogaster how gene Proteins 0.000 description 3
- 230000001174 ascending effect Effects 0.000 description 3
- 238000004891 communication Methods 0.000 description 3
- 230000010354 integration Effects 0.000 description 3
- 230000002093 peripheral effect Effects 0.000 description 3
- 210000005245 right atrium Anatomy 0.000 description 3
- 238000007619 statistical method Methods 0.000 description 3
- 238000001356 surgical procedure Methods 0.000 description 3
- 210000001367 artery Anatomy 0.000 description 2
- 210000004369 blood Anatomy 0.000 description 2
- 239000008280 blood Substances 0.000 description 2
- 230000004087 circulation Effects 0.000 description 2
- 210000005246 left atrium Anatomy 0.000 description 2
- 210000005240 left ventricle Anatomy 0.000 description 2
- 230000000737 periodic effect Effects 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 239000000779 smoke Substances 0.000 description 2
- 230000007704 transition Effects 0.000 description 2
- 101100171060 Caenorhabditis elegans div-1 gene Proteins 0.000 description 1
- 208000002330 Congenital Heart Defects Diseases 0.000 description 1
- 230000005856 abnormality Effects 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 238000010009 beating Methods 0.000 description 1
- 238000005452 bending Methods 0.000 description 1
- 230000000747 cardiac effect Effects 0.000 description 1
- 210000004413 cardiac myocyte Anatomy 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000005094 computer simulation Methods 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 238000002001 electrophysiology Methods 0.000 description 1
- 230000007831 electrophysiology Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000014509 gene expression Effects 0.000 description 1
- 210000002837 heart atrium Anatomy 0.000 description 1
- 208000019622 heart disease Diseases 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 230000010349 pulsation Effects 0.000 description 1
- 238000005086 pumping Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 210000005241 right ventricle Anatomy 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000008685 targeting Effects 0.000 description 1
- 230000001225 therapeutic effect Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Description
記憶部は、流体シミュレーションの解析時間の進行に伴い解析空間内を加速度運動する複数のグリッド点の第1時間間隔の複数の第1時刻の各々における位置を示す位置情報と、複数のグリッド点上の流体の複数の第1時刻の各々における速度を示す速度情報とを含む流体情報を記憶する。処理部は、まず流体情報に基づいて、位置情報に表される複数のグリッド点の加速度運動に起因する誤差を補正する補正値を含む計算式により、複数のグリッド点上の流体の複数の第1時刻の各々における速度の時間微分値をそれぞれ算出する。次に処理部は、複数のグリッド点上の流体の複数の第1時刻の各々における速度と速度の時間微分値とに基づいて、解析時間の進行に伴い粒子発生源から順番に出力される複数の粒子の、第1時間間隔より短い第2時間間隔の複数の第2時刻の各々における位置をそれぞれ算出する。そして処理部は、複数の粒子の連なりを示す流脈線の表示情報を生成する。
〔第1の実施の形態〕
まず、第1の実施の形態について説明する。第1の実施の形態は、構造体が変形しても流脈線の追跡が可能な流脈線表示装置である。
〔第2の実施の形態〕
次に第2の実施の形態について説明する。第2の実施の形態は、心臓の動きと共に、心臓内の血流の流脈線を可視化する可視化装置である。
図4は、流脈線の計算例を示す図である。可視化装置100は、解析対象の空間中に粒子発生源30を定義する。解析が始まると、粒子群33が粒子発生源30から連続的に射出される。粒子群33は流れ場が時間に対して不変であれば、一定の曲線(流線)になる。流れ場が時間依存して変化すると、粒子群33が作り出す曲線も時々刻々と変化していく。この曲線が流脈線31,32である。図4では、時刻t0における粒子群33の連なりを示す流脈線31と時刻t1における粒子群33の連なりを示す流脈線32とが示されている。流脈線31,32は、障害物35の影響で複雑な曲線となる。
1.心筋(弾性体)が大きく変形するとき、心筋付近で流跡線、流脈線の挙動を正確に追跡することが困難なことである。
3.有限要素モデルの一部の品質の悪いメッシュが、全体の計算量を増加させてしまうことである。
1−1:可視化装置100は、流脈線の各点が可動領域外の心筋に入り込んだり、シミュレーション対象の系の外に出したりしたかどうかを正確に判定する。
1−3:可視化装置100は、任意の時刻における場の情報を推定するため、補間法を用いて場を補間する。
2−1:1−3の機能を適用すると計算量が増加するため、可視化装置100は、流脈線の点が動く最大距離を計算し、その距離を半径とする予測球内部の場の情報のみを計算することでデータ容量に関係なく計算量を一定にする。
3−1:大部分の計算は品質のよいメッシュに対して行われるのを利用し、可視化装置100は、品質のよいメッシュと決めて投機的に計算を実行し、失敗した場合には高精度計算を行うことで、計算量を削減する。投機的な計算とは、流脈線を構成する点の移動先が予測球の外部になる可能性があることを許容して、予測球の半径を小さくすることである。計算に失敗する場合とは、流脈線を構成する点の移動先が予測球内に存在しない場合である。
4−1:見かけの力を反映させた流脈線の計算を効率的に実施するため、見かけの力に応じた係数を事前に計算し、時間ステップごとに、その計算結果を参照して流脈線を計算する。
1.心筋(弾性体)が大きく変形をする場合にも心筋(弾性体)の運動を考慮して流脈線を計算できる。
3.確率モデルを用い、計算コストを最小化する予測球半径の決定ができる。
図5は、可視化装置の機能を示すブロック図である。可視化装置100は、情報の記憶機能として、シミュレーション結果記憶部110と、事前解析結果記憶部120とを有する。シミュレーション結果記憶部110は、心臓シミュレータ200から取得したシミュレーション結果を記憶する。例えば数値流体力学シミュレーションが心臓シミュレータ200で実施され、動的に変化する弾性体と流体の場のシミュレーション結果がL個(Lは1以上の整数)の時刻t0,t1,・・・,tLでファイルに保存される。例えば心筋の情報と血流の情報が別々のファイルとして、シミュレーション結果記憶部110に保存される。図5の例では、時刻ごとの心筋の情報が弾性体情報ファイル群111として保存され、時刻ごとの血流の情報が流体情報ファイル群112として保存されている。
次に、シミュレーション結果として得られる情報について詳細に説明する。
心筋側補間多項式係数テーブル121bには、グリッドのID(GRID ID)と時刻のインデックス(Time Index)と軸のインデックス(Direction Index)の組に対応付けて、軸方向の位置(Position)、心筋側補間多項式に用いる各係数(係数a、係数b、係数c、係数d)の値が設定されている。
可視化装置100は、図6と図7に示したシミュレーション結果と、図8に示した事前解析結果Ωとに基づいて、流脈線を計算し可視化する。以下、流脈線の可視化処理について詳細に説明する。
[ステップS101]流脈線計算部150は、事前解析結果記憶部120に事前解析ファイル121が格納されているか否かを判断する。流脈線計算部150は、事前解析ファイル121が格納されている場合、処理をステップS102に進める。また流脈線計算部150は、事前解析ファイル121が格納されていない場合、処理をステップS103に進める。
[ステップS104]事前解析部140は、事前解析処理によって算出された係数の値を含む事前解析ファイル121を生成し、事前解析ファイル121を事前解析結果記憶部120に格納する。
[ステップS105]流脈線計算部150は、計算する流脈線の本数M(Mは1以上の整数)を決定する。例えば流脈線計算部150は、ユーザから入力された値を、流脈線の本数Mとする。
j番目の流脈線ljは時間発展回数Nと同じ分だけの離散点で構成される。そこで流脈線計算部150は、流脈線ljに含まれる各離散点を示す点Pij(i=0,1,2,・・・,N)を生成する。流脈線計算部150は、各離散点の初期値の座標を、粒子発生源の座標とする。
[ステップS121]事前解析部140は、インデックスi=1から順番に、インデックスi=1,2,・・・,N−1について、ステップS122のループ処理を繰り返し実行する。
そして事前解析部140は、係数al,i,bl,i,cl,i,dl,iの配列を事前解析結果Ωに含めて保存する。Y,Z成分についても同様の手続きで補間関数を求めることができる。事前解析部140は、このような係数の算出を、すべての心筋のグリッド点と血流のグリッド点について行う。その結果、ベクトルdrl/dtの値が計算されたこととなる。
図11は、ALE制御されたグリッド点自身の運動による見かけの力の原理を示す図である。ALE制御されたグリッド点は、心臓の収縮に合わせて1拍動すると元の位置に戻り、心臓の拍動に合わせて運動する。このとき、グリッド点は速度が時々刻々と変化するので、加速度が生じる。ところがグリッド点自身に観測装置を設置すれば、観測者は自らを静止したと解釈してしまう。自分自身の加速度運動による効果は見かけの力としてあたかも流体に働くように見えてしまう。そのため、自分自身の加速度運動による補正をしなければならない。この見かけの力を補正するための手続きが式(5)、式(6)、式(7)には含まれている。このようにして、各グリッド点における時刻t=tiでの速度場の時間微分を計算することができる。事前解析部140は、このような速度場の時間微分をすべてのグリッド点での血流の速度場について求め、各グリッド点での速度場の時間微分を含むデータセットを、事前解析結果Ωとして事前解析結果記憶部120に保存する。
図12は、時間発展計算処理の手順を示すフローチャートの前半部分である。以下、図12に示す処理をステップ番号に沿って説明する。
[ステップS133]流脈線計算部150は、計算開始時刻t=tiおよび計算終了時刻t=ti+1を、メモリ102から読み込む。
[ステップS141]流脈線計算部150は、中間時刻(t=tk(k=1,2,・・・,Ndiv))の、k=1からk=Ndivまで、ステップS142〜S150のループ処理を繰り返し実行することで、時間発展計算を行う。これにより、各中間時刻t=tkにおける座標ベクトルrkが得られ、時刻t=ti+1における点Pi+1jの座標ベクトルrNdiv=ベクトルri+1が得られる。
ステップS142〜S144における各中間時刻での時間発展は以下のようにして行うことができる。点Pijが時刻t=tkで、ベクトルrk(tk)にあったとすると、流脈線の方程式は式(8)で与えられる。
図14は、ルンゲ・クッタ係数表の一例を示す図である。図14に示すルンゲ・クッタ係数表には、I(Iは1〜4の整数)の値ごとの、係数αIと係数βIの値が示されている。図14に示すようなルンゲ・クッタ係数表は、例えばメモリ102に予め記憶されている。
[ステップS145]流脈線計算部150は、ベクトルrk+1の心臓内部での位置判定を行う。この処理を行うのは、求めたベクトルrk+1は有限の時間幅で誤差を含んでおり、心筋内部に入り込んでしまう場合があるためである。サブルーチンとして位置判定を行った結果、位置判定の結果を示すステータスが返される。例えば心筋を跨いだり、心筋内部に入り込んだりしていなければ正常終了ステータスである「0」または「2」が返される。ベクトルrk+1が解析対象の流体の要素内の位置(例えば心房・心室中)を示している場合、ステータス「0」が返される。ベクトルrk+1が解析対象の流体の要素外の位置(例えば動脈中)を示している場合、ステータス「2」が返される。この処理の詳細は後述する(図23参照)。
[ステップS151]流脈線計算部150は、ステップS142〜S150のループ処理が1回終わるごとに、kに1を加算し、ループ処理を繰り返す。流脈線計算部150は、k=Ndivの処理が終了すると、処理をステップS152に進める。
図12、図13に示した処理により、流脈線上の点の座標が更新されていく。
[ステップS161]流脈線計算部150は、時刻tkにおける場の情報をメモリに設定する。設定される情報は、ベクトルB(ベクトルr,ti)、ベクトルv(ベクトルr,ti)、ベクトルB(ベクトルr,ti+1)、およびベクトルv(ベクトルr,ti+1)である。
[ステップS163]流脈線計算部150は、予測球半径R内部にある流体のグリッド点を検索し、予測球内部に存在するグリッド点の数をNB,elemに設定する。
[ステップS165]流脈線計算部150は、グリッド点iが境界面上か否かを判断する。境界面上であれば、流脈線計算部150は、処理をステップS167に進める。境界面上でなければ、流脈線計算部150は、処理をステップS166に進める。
グリッド点iが境界面S(t)上の点である場合、その速度場から構成される位相空間上の点を(ベクトルrID(t),ベクトルvID(t))とする。すると、時刻tiおよびti+1では値が出力されているため2点(ベクトルrID(ti),ベクトルvID(ti))、(ベクトルrID(ti+1),ベクトルvID(ti+1))を通る曲線として求めればよい。グリッド位置ベクトルrID(t)と速度場ベクトルvID(t)について、以下の式(11)の関係が成立する。
図18は、時刻t=tkにおける弾性体場の情報を求める処理の手順の一例を示すフローチャートである。以下、図18に示す処理をステップ番号に沿って説明する。
[ステップS174]流脈線計算部150は、グリッド点iが境界面S(t)の点であるか否かを判断する。境界面S(t)上であれば、流脈線計算部150は、処理をステップS176に進める。境界面S(t)上でなければ、流脈線計算部150は、処理をステップS175に進める。
なお流体の処理を先に行った場合、流体側に位置が同じになる対応するグリッド点iが存在する。そのため流脈線計算部150は、ステップS176,S177の処理に代えて、ステップS167における計算結果から、グリッド点iの座標を取得し、保存する処理を行ってもよい。
この速度場ベクトルvl(ベクトルrl(t),t)は、流脈線上の点の時間発展をルンゲ・クッタ法で行うための式(9)、式(10)の計算に用いられる。出力データが与えられている時刻をtiとすると、ti≦t≦ti+1を満たす整数iを決定することができる。事前解析結果Ωには、各頂点の時刻tiにおけるグリッド座標ベクトルrl(ti)の値、その点での速度場ベクトルvl(ベクトルrl(ti),ti)の値、および速度場の時間微分ベクトルdv(ベクトルrl(ti),ti)/dtが保存されている。従って、速度場は時刻tiとti+1で速度場とその微分が与えられていることになる。そのため、微分まで含めて連続な曲線を補間曲線として仮定することができる。与えられた2つの時刻での速度場とその時間微分が与えられているため、式(20)のような未知変数を4つとする3次の多項式を補間曲線として用いることができる。
次に、任意の時刻tにおける血流中の任意の座標ベクトルrでの速度場ベクトルv(ベクトルr,t)の計算方法について説明する。
図19は、四面体要素の一例を示す図である。四面体要素Ilの各頂点の座標が、それぞれP1(x1,y1,z1)、P2(x2,y2,z2)、P3(x3,y3,z3)、P4(x4,y4,z4)で与えられるものとする。このとき、構造関数は式(28)で与えられる。
図13のステップS142,S143で任意の時刻における心筋と血流のグリッド点iの位置を求めた後、ステップS144で流脈線上の点ベクトルr(tk)の時間発展がルンゲ・クッタ法で行われる。
[ステップS181]流脈線計算部150は、インデックスi=1から順番に、インデックスi=1,2,3,4について、ステップS182〜S183のループ処理を繰り返し実行する。
[ステップS183]流脈線計算部150は、ルンゲ・クッタ法で用いる中間の速度ベクトルvi(i=1,2,3,4)を、式(9)、式(10)で計算する。ここでベクトルviの計算には、データが与えられていない任意の時刻での流体中の任意の座標における速度場の計算も行われる。例えば、ベクトルv1は時刻tkにおける座標ベクトルr(tk)における速度場の値であるが、時刻tkは事前解析によりデータが与えられている時刻とは限らない。また、座標ベクトルrkもデータが与えられている流体中のグリッド点iの座標とは限らない。従って、任意の時刻tにおける、任意の座標ベクトルrでの速度場の計算を行うこととなる。
次に、心臓内部での位置判定について詳細に説明する。
2つ目は移動元の点Pと移動後の点Qを結んだ直線PQが心筋または心筋表面と交わるかどうかを判定する直線判定である。直線PQが心筋(表面)と交わればTを返し、交わらない場合はFを返す。また直線判定では、同時に交点の個数も返す。
点Qが流体中に存在せず、直線PQが心筋(表面)と交わる場合、点Qは(1)心筋を突き破って系外に出た場合と(2)心筋内部にめり込んだ場合との2通りが存在する。どちらも再計算をすることになるため、流脈線計算部150は、状態変数として「1」を設定する。
このような状態判定の処理手順を以下に示す。
[ステップS191]流脈線計算部150は、時間発展後の座標ベクトルrk+1が流体中に存在するかどうか判定する。この処理の詳細は後述する(図24参照)。
次に、時間発展後の位置が流体か否かの判定処理(ステップS191)について詳細に説明する。
[ステップS202]流脈線計算部150は、時刻tk+1での座標ベクトルrk+1を取得する。
[ステップS204]流脈線計算部150は、座標ベクトルrkを中心とした半径Rの予測球内部の要素のリスト(流体要素リストLf)を取得する。
[ステップS206]流脈線計算部150は、結果(result)の初期値として「F(False)」を設定する。
[ステップS208]流脈線計算部150は、要素リスト内の要素Liに対して、要素Liが座標ベクトルrk+1を含むかどうかを判定する。座標ベクトルrk+1を含む場合、流脈線計算部150は、処理をステップS210に進める。座標ベクトルrk+1を含まない場合、流脈線計算部150は、処理をステップS209に進める。
[ステップS211]流脈線計算部150は、resultをTrueに変更する。
図25は、移動ベクトルが作る線分と心筋の交点数の計数処理の手順の一例を示すフローチャートである。以下、図25に示す処理をステップ番号に沿って説明する。
[ステップS222]流脈線計算部150は、時刻tk+1での座標ベクトルrk+1を取得する。
[ステップS225]流脈線計算部150は、予測球半径内部にある弾性体の要素リスト(弾性体要素リストLe)を作成する。
[ステップS227]流脈線計算部150は、結果(result)の初期値として「F(False)」に設定する。
[ステップS229]流脈線計算部150は、i番目の要素Liが移動ベクトルdrと交わるか否かを判断する。要素Liは多面体であるため、流脈線計算部150は、要素Liに含まれるすべての面に対して移動ベクトルdrと交点を求める。要素内のすべての面に対して、交点が1点も存在しなければ、流脈線計算部150は、その要素とは交わらないと判定する。交わらない場合、流脈線計算部150は、処理をステップS230に進める。交わる場合、流脈線計算部150は、処理をステップS231に進める。
[ステップS232]流脈線計算部150は、線分を示すベクトルdrと心筋表面との交点数Z(Zは1以上の整数)を求め、メモリ102に保存する。
このようにして、交点が1点も存在しない場合、交わらないとしてFが返され、交点が1点以上存在する場合、交わるとしてTが返される。交わる場合は交点を持つ要素の番号と、心筋表面との交点数Zが保存される。
[ステップS241]流脈線計算部150は、要素の辺の長さと速度場の統計解析から予測球半径の最小値Rminを決定する。
[ステップS243]流脈線計算部150は、時刻tiでの速度の最大値|ベクトルvi,max|を取得する。
[ステップS245]流脈線計算部150は、|ベクトルvi,max|が|ベクトルvi+1,max|以上かどうかを判断する。流脈線計算部150は、|ベクトルvi,max|が|ベクトルvi+1,max|以上であれば、処理をステップS246に進める。流脈線計算部150は、|ベクトルvi,max|が|ベクトルvi+1,max|未満であれば、処理をステップS247に進める。
[ステップS247]流脈線計算部150は、|ベクトルvi+1,max|を|ベクトルvmax|に設定する。流脈線計算部150は、その後、処理をステップS248に進める。
[ステップS249]流脈線計算部150は、|ベクトルvmax|dtを計算し、計算結果を予測球半径Rに設定する。
以上のような処理により、適切な予測球半径Rを決定することができる。
[分割法による計算の高精度化と高速化]
図12のステップS137では、出力ファイルが与えられた時間ti≦t≦ti+1の時間分割をさらにNdiv分割している。この分割の仕方について詳細に説明する。
ここで時間ti≦t≦ti+1をNdiv分割することを考える。Ndiv分割すると、ルンゲ・クッタ法による時間発展回数はNdiv回となる。一方で、予測球半径は式(30)から決定される。そのため、時間刻み幅がNdiv分の1になると、予測球半径もNdiv分の1になる。1回の計算量は、予測球半径R3に比例するため、Ndiv -3になる。この計算をNdiv回繰り返すため、合計の計算量はNdiv -2になる。この計算の概要を模式化したのが図27である。
流体シミュレーションは有限個の離散点情報で行われる。従って、予測球半径を小さくしすぎると予測球内部に離散点情報が含まれない状況が生じうる。これは予測球半径に下限値が存在することを意味する。同時に、予測球半径は式(30)から決定されるため、時間刻み幅にも下限値が存在する。予測球を用いた計算では、予測球半径の下限値Rminを設定し、その予測球半径が下限値以下の場合になるときは、Rminを予測球半径として設定する。なお、予測球の下限値以下になった時に、それより大きな値であるRminを使用するので、流脈線上の点が予測球外部に出てしまうことはない。このように、計算を安定に進めるためにRminの決定は重要である。また、Rminの値は時間分割数Ndivの決定にも関与する。Rminの値が定まれば、R0=|Δtベクトルvmax|が最大移動距離であるため、天井関数を使って以下の式で分割数Ndivを決定できる。
しかし、下限値Rminの決め方には工夫がいる。具体的には、確率モデルを導入し、流脈線計算部150は、計算実行時に失敗することを許容した投機的な計算を行う。計算に失敗したときには計算が確実に成功するパラメータで計算を実行し、その計算にかかる時間がペナルティとしてかかると考える。そして、ペナルティまで含めて統計的に計算時間を最小にするように下限値Rminを決定する。詳細を以下に記す。
図29は、要素の辺の長さの分布の一例を示す図である。図29のグラフは、横軸が要素の辺の長さであり、縦軸が、辺の長さごとの該当する辺の存在確率を示す確率密度である。グラフ中には、2種類の心臓シミュレーション#1,#2における確率密度関数61,62を破線で示している。
図30は、予測球半径に応じた計算コストの変化を示す図である。図30の縦軸が計算コスト(T’/T)であり、横軸が予測球半径Rである。図30に示すように、R=0.003[m]付近でペナルティまで含めたコストが最悪値半径Rmax=0.008[m]と比べて約10%まで落ちることがわかる。そして、ペナルティまで含めて統計的に最小の計算量で済む半径が存在することがわかる。実際には、有限回の計算であるため、すべての要素を等確率に解析対象とするという状態にはならない。そのため、これよりも半径の値は小さくなる。
以下、心臓シミュレータによる心筋と冠循環の解析を用いて、流脈線を実際に計算した場合の計算速度について具体的に説明する。
流脈線の計算では、計算量の削減のため、および計算精度の向上のため、予測球が使用されている。予測球を使った計算では、時間刻み幅Δt=0.01[sec]をさらにNdiv分割して細分化した。Ndivの値は、統計的に計算量が最小になるように自動的に最適値が決定される。このシミュレーション結果の場合、Ndivの値はNdiv=3〜7程度に設定された。
図35は、ALEグリッド点の運動補正を加えた速度場の例を示す図である。図35には、速度場の評価を線形補間で行った場合の速度変化(点線)と、式(4)でALEグリッド座標自身の運動により生じた見かけの力の補正を入れた場合の速度変化(実線)との計算結果を示している。
〔その他の実施の形態〕
流脈線上の各点はお互いに相互作用せず、独立である。そのため、並列化に適している。そこでMPI(Message Passing Interface)やOpenMP(Multi-Processing)を用いて流脈線の計算を並列化させてもよい。これにより計算速度を向上させることができる。
2 構造体
3a,3b 流体領域
4 粒子発生源
5a,5b 流脈線
10 流脈線可視化装置
11 記憶部
12 処理部
13 表示部
Claims (7)
- 流体シミュレーションの解析時間の進行に伴い解析空間内を加速度運動する複数のグリッド点の第1時間間隔の複数の第1時刻の各々における位置を示す位置情報と、前記複数のグリッド点上の流体の前記複数の第1時刻の各々における速度を示す速度情報とを含む流体情報を記憶する記憶部と、
前記流体情報に基づいて、前記位置情報に表される前記複数のグリッド点の加速度運動に起因する誤差を補正する補正値を含む計算式により、前記複数のグリッド点上の前記流体の前記複数の第1時刻の各々における前記速度の時間微分値をそれぞれ算出し、前記複数のグリッド点上の前記流体の前記複数の第1時刻の各々における前記速度と前記速度の時間微分値とに基づいて、前記解析時間の進行に伴い粒子発生源から順番に出力される複数の粒子の前記第1時間間隔より短い第2時間間隔の複数の第2時刻の各々における位置をそれぞれ算出し、前記複数の粒子の連なりを示す流脈線の表示情報を生成する処理部と、
を有する流脈線可視化装置。 - 前記処理部は、前記位置情報に基づいて、前記複数のグリッド点の各々に対し、前記複数の第1時刻の各々におけるグリッド点移動速度とグリッド点移動加速度とをそれぞれ算出し、前記グリッド点移動速度と前記グリッド点移動加速度とを変数として含む前記計算式により、前記複数のグリッド点上の前記流体の前記複数の第1時刻の各々における前記速度の時間微分値をそれぞれ算出する、
請求項1記載の流脈線可視化装置。 - 前記処理部は、前記複数のグリッド点の各々について、グリッド点の前記複数の第1時刻の各々における前記流体の前記速度を示す点をそれぞれ接続し、前記点を通るとき、対応する第1時刻での前記流体の前記速度の時間微分値に応じた傾きとなる、前記流体の前記速度の時間変化を表す補間曲線を生成し、前記複数のグリッド点の各々の前記補間曲線に基づいて、前記複数の粒子の前記複数の第2時刻の各々における位置をそれぞれ算出する、
請求項1または2記載の流脈線可視化装置。 - 前記記憶部は、さらに、前記解析空間内の構造体の形状の時間変化を示す構造体情報を記憶し、
前記処理部は、第1解析時刻における第1流脈線に基づいて第2解析時刻における第2流脈線を計算する場合、
前記第1流脈線上の第1の位置にある離散点を含む、前記解析空間内の部分領域を、前記離散点についての解析対象領域に決定し、
前記流体情報に示される前記解析対象領域内の前記流体の前記速度に基づいて、前記第2解析時刻における前記離散点上の粒子の移動先を示す第2の位置を計算し、
前記構造体情報に示される前記解析対象領域内に存在する前記構造体の情報に基づいて、前記第2解析時刻における前記解析対象領域内の前記構造体の占有領域を特定し、
前記第1の位置と前記第2の位置とに基づいて、前記第2流脈線の前記占有領域への侵入の有無を判定し、
前記第2流脈線が前記占有領域に侵入しないと判定したとき、前記第2の位置を通過する前記第2流脈線の表示情報を生成する、
請求項1乃至3のいずれかに記載の流脈線可視化装置。 - さらに、前記表示情報に基づいて、前記流脈線を表示する表示部を有する、
請求項1ないし4のいずれかに記載の流脈線可視化装置。 - コンピュータが、
流体シミュレーションの解析時間の進行に伴い解析空間内を加速度運動する複数のグリッド点の第1時間間隔の複数の第1時刻の各々における位置を示す位置情報と、前記複数のグリッド点上の流体の前記複数の第1時刻の各々における速度を示す速度情報とを含む流体情報に基づいて、前記位置情報に表される前記複数のグリッド点の加速度運動に起因する誤差を補正する補正値を含む計算式により、前記複数のグリッド点上の前記流体の前記複数の第1時刻の各々における前記速度の時間微分値をそれぞれ算出し、
前記複数のグリッド点上の前記流体の前記複数の第1時刻の各々における前記速度と前記速度の時間微分値とに基づいて、前記解析時間の進行に伴い粒子発生源から順番に出力される複数の粒子の前記第1時間間隔より短い第2時間間隔の複数の第2時刻の各々における位置をそれぞれ算出し、
前記複数の粒子の連なりを示す流脈線の表示情報を生成する、
流脈線可視化方法。 - コンピュータに、
流体シミュレーション上の解析時間の進行に伴い解析空間内を加速度運動する複数のグリッド点の第1時間間隔の複数の第1時刻の各々における位置を示す位置情報と、前記複数のグリッド点上の流体の前記複数の第1時刻の各々における速度を示す速度情報とを含む流体情報に基づいて、前記位置情報に表される前記複数のグリッド点の加速度運動に起因する誤差を補正する補正値を含む計算式により、前記複数のグリッド点上の前記流体の前記複数の第1時刻の各々における前記速度の時間微分値をそれぞれ算出し、
前記複数のグリッド点上の前記流体の前記複数の第1時刻の各々における前記速度と前記速度の時間微分値とに基づいて、前記解析時間の進行に伴い粒子発生源から順番に出力される複数の粒子の前記第1時間間隔より短い第2時間間隔の複数の第2時刻の各々における位置をそれぞれ算出し、
前記複数の粒子の連なりを示す流脈線の表示情報を生成する、
処理を実行させる流脈線可視化プログラム。
Priority Applications (4)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2017075671A JP6760603B2 (ja) | 2017-04-06 | 2017-04-06 | 流脈線可視化装置、流脈線可視化方法、および流脈線可視化プログラム |
EP17203528.9A EP3333738A1 (en) | 2016-12-06 | 2017-11-24 | Streakline visualization apparatus, method, and program |
US15/827,579 US10799191B2 (en) | 2016-12-06 | 2017-11-30 | Streakline visualization apparatus and method |
CN201711261432.8A CN108170895B (zh) | 2016-12-06 | 2017-12-04 | 纹线可视化设备和方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2017075671A JP6760603B2 (ja) | 2017-04-06 | 2017-04-06 | 流脈線可視化装置、流脈線可視化方法、および流脈線可視化プログラム |
Publications (2)
Publication Number | Publication Date |
---|---|
JP2018180707A JP2018180707A (ja) | 2018-11-15 |
JP6760603B2 true JP6760603B2 (ja) | 2020-09-23 |
Family
ID=64276780
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2017075671A Active JP6760603B2 (ja) | 2016-12-06 | 2017-04-06 | 流脈線可視化装置、流脈線可視化方法、および流脈線可視化プログラム |
Country Status (1)
Country | Link |
---|---|
JP (1) | JP6760603B2 (ja) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113077435A (zh) * | 2021-03-30 | 2021-07-06 | 昆明同心医联科技有限公司 | 基于血流动力学的动脉瘤分析方法及装置 |
WO2023132392A1 (ko) * | 2022-01-07 | 2023-07-13 | 이에이트 주식회사 | 입자 기반의 시뮬레이션을 통해 경동맥에서의 혈류 특성을 분석하는 방법 및 시스템 |
CN116167221B (zh) * | 2023-02-21 | 2023-09-01 | 哈尔滨工程大学 | 基于完全信息熵的自适应步长流线生成方法、计算机设备和存储介质 |
CN116719982B (zh) * | 2023-08-08 | 2023-10-13 | 东莘电磁科技(成都)有限公司 | 一种基于电磁领域的线网格流场可视化方法 |
Family Cites Families (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP3174670B2 (ja) * | 1993-08-26 | 2001-06-11 | 三菱重工業株式会社 | 3次元流動解析結果表示方法及び装置 |
JP2625644B2 (ja) * | 1994-06-28 | 1997-07-02 | インターナショナル・ビジネス・マシーンズ・コーポレイション | 流線表示方法及びコンピュータ・システム |
JP2007135894A (ja) * | 2005-11-18 | 2007-06-07 | R Tech:Kk | ヒト血流データをもとにした血流解析装置及びシミュレーション方法 |
US8761474B2 (en) * | 2011-07-25 | 2014-06-24 | Siemens Aktiengesellschaft | Method for vascular flow pattern analysis |
JP6300244B2 (ja) * | 2014-07-03 | 2018-03-28 | 富士通株式会社 | 生体シミュレーション装置、生体シミュレーション装置の制御方法、および生体シミュレーション装置の制御プログラム |
-
2017
- 2017-04-06 JP JP2017075671A patent/JP6760603B2/ja active Active
Also Published As
Publication number | Publication date |
---|---|
JP2018180707A (ja) | 2018-11-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
García-Villalba et al. | Demonstration of patient-specific simulations to assess left atrial appendage thrombogenesis risk | |
JP6760603B2 (ja) | 流脈線可視化装置、流脈線可視化方法、および流脈線可視化プログラム | |
US10595790B2 (en) | Method and system for personalized non-invasive hemodynamic assessment of renal artery stenosis from medical images | |
US10366490B2 (en) | Highly integrated annotation and segmentation system for medical imaging | |
US10799191B2 (en) | Streakline visualization apparatus and method | |
EP3166032B1 (en) | Biometric simulation device, method for controlling biometric simulation device, and program for controlling biometric simulation device | |
EP2662794B1 (en) | Simulation method, simulator apparatus, and simulation program for hemodynamics in vessels | |
Georgiou et al. | Learning fluid flows | |
Oeltze‐Jafra et al. | Generation and visual exploration of medical flow data: Survey, research trends and future challenges | |
CN112150454B (zh) | 一种主动脉夹层评估方法、装置、设备及存储介质 | |
JP6741252B2 (ja) | 流脈線可視化装置、流脈線可視化方法、および流脈線可視化プログラム | |
Sankaran et al. | Physics driven real-time blood flow simulations | |
Pajaziti et al. | Shape-driven deep neural networks for fast acquisition of aortic 3D pressure and velocity flow fields | |
JP6253053B2 (ja) | データ探索装置、データ探索装置の制御方法およびデータ探索装置の制御プログラム | |
CN112381824B (zh) | 一种对图像的几何特征进行提取的方法和相关产品 | |
Mufti et al. | Shock wave prediction in transonic flow fields using domain-informed probabilistic deep learning | |
CN112381822B (zh) | 一种用于处理肺部病灶区图像的方法和相关产品 | |
JP6748368B2 (ja) | モデリング装置、モデリング方法、およびモデリングプログラム | |
Wang et al. | Fully parallelized Lattice Boltzmann scheme for fast extraction of biomedical geometry | |
Hong et al. | An implicit skeleton-based method for the geometry reconstruction of vasculatures | |
US20210319312A1 (en) | Deep learning acceleration of physics-based modeling | |
Tanade et al. | Cloud computing to enable wearable-driven longitudinal hemodynamic maps | |
Chen et al. | Fast coherent particle advection through time-varying unstructured flow datasets | |
KR102516382B1 (ko) | Lbm 기반의 혈류 시뮬레이션 장치, 방법 및 컴퓨터 프로그램 | |
Kretschmer et al. | Reliable Adaptive Modelling of Vascular Structures with Non‐Circular Cross‐Sections |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
RD01 | Notification of change of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7426 Effective date: 20170501 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A821 Effective date: 20170501 Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20170501 |
|
A80 | Written request to apply exceptions to lack of novelty of invention |
Free format text: JAPANESE INTERMEDIATE CODE: A80 Effective date: 20170501 |
|
RD02 | Notification of acceptance of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7422 Effective date: 20190927 |
|
RD04 | Notification of resignation of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7424 Effective date: 20190927 |
|
A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20200106 |
|
TRDD | Decision of grant or rejection written | ||
A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20200811 |
|
A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20200827 |
|
R150 | Certificate of patent or registration of utility model |
Ref document number: 6760603 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |