JP2015156108A - 非線形変形計算システム及び手術シミュレータ - Google Patents

Abstract

【課題】木構造で表現される多重解像度の変形モデルの共回転系変形計算を行う場合、木構造の適切な階層において回転成分を抽出できるようにする。【解決手段】変形モデルを構成する複数のメッシュ要素の木構造を示す第１記憶部３５ａと、回転成分を抽出すべき複数の回転成分抽出対象メッシュ要素に対応する複数の対象ノードを示す第２記憶部３５ｂと、木構造のリーフノードそれぞれに対応する共回転系変形計算を、各リーフノードの先祖ノードである対象ノードに対応する変形計算部３２ａと、回転成分抽出対象メッシュ要素の回転抽出部３２ｂと、前記対象ノードを動的に変更する対象ノード変更部３２ｃとを備え、対象ノード変更部３２ｃは、回転成分抽出対象メッシュ要素から抽出される回転成分の大きさの評価を行い、その評価結果に基づいて、当該回転成分抽出対象メッシュ要素に対応する対象ノードの子ノード又は親ノードを新たな対象ノードに設定する。【選択図】図２

Description

本発明は、非線形変形計算システム及び手術シミュレータに関するものである。

手術シミュレータには、臓器の高速な変形シミュレーションが必要となる。臓器に触れたり臓器を掴んだりすることで臓器に生じる動的な変形に基づく手応え（力覚）をユーザに返すには、高い更新レートで変形計算及び反力計算を行う必要がある。
また、手術シミュレーションにおいては、臓器等の実物体と整合性のとれた変形シミュレーションが可能であること、特に、幾何学的な非線形性を考慮可能であること（回転を含む変形が発生するため）が望まれる。
したがって、幾何学的非線形性を考慮可能な、低計算量の変形モデルが必要とされている。

本発明者らは、幾何学的非線形特性の考慮が可能なモデルとして共回転系変形モデル（非特許文献１参照）を用いることを、非特許文献２において提案している。
共回転系変形モデルでは、変形計算に用いる変形前の四面体要素と変形後の四面体要素との間で形状のマッチングを行い、ローカルな座標を定義し直すことで、回転を含む運動成分を除去し、その後、線形有限要素法を用いて変形計算を行う。
しかし、共回転系変形モデルを用いる場合には、前述のマッチング（回転を含む運動成分の抽出）に多くの計算量を要する問題点がある。

ここで、本発明者らは、変形モデルとして適応的四面体を用いた多重解像度モデルを用いる場合に、適応的四面体メッシュの二分木階層構造を用いて、前述のマッチング計算の削減を行う手法を、非特許文献２にて提案している。ここで、多重解像度モデルでは、モデルの境界部分や変形の大きい部分など特徴が大きく変化する部分はメッシュ要素の分割数を大きくして、指定された近似精度を維持することができる変形モデルである。

適応的四面体メッシュにおいて構築される二分木構造が構築される過程では、一つの親四面体から２つの子四面体、４つの孫四面体等が作り出される。そのため、子、孫やひ孫等の子孫の四面体の回転成分（回転行列）は、これら先祖の四面体の回転成分（回転行列）に類似しているという考えの下、非特許文献２では、子孫の四面体の回転行列を先祖の四面体の回転行列で代用することでマッチング計算の回数を削減している。

M.Muller, J.Dorsey, L.McMillan, R.Jagnow and B.Cutler,"Stable real-time deformations", Proc.ACM SCA ’02, pp.49-54, 2002. 田川和義，山田隆洋，田中弘美，「オンラインリメッシュ型共回転系変形モデルを用いたインタラクティブな手術シミュレーション」，可視化情報学会，可視化情報シンポジウム，新宿区・工学院大学，2013年7月16日〜17日

しかし、非特許文献２は、二分木のリーフノードに相当する四面体の回転行列を、所定の分割レベルにおける先祖の四面体の回転行列で代用したものを開示しているにすぎず、木構造で表現される多重解像度モデルにおいて、具体的に、木構造のどの階層において回転成分を抽出すべきであるかという基準については開示されていない。

本発明は、木構造で表現される多重解像度の変形モデルの共回転系変形計算を行う場合において、木構造の適切な階層において回転成分を抽出できるようにすることを目的とする。

一の観点からみた本発明は、動的に変形が生じる多重解像度の変形モデルの変形計算を行う非線形変形計算システムであって、前記変形モデルを構成する複数のメッシュ要素の木構造を示す木構造情報を記憶する第１記憶部と、前記木構造において、回転成分を抽出すべき複数の回転成分抽出対象メッシュ要素に対応する複数の対象ノードを示す対象ノード情報を記憶する第２記憶部と、前記木構造のリーフノードそれぞれに対応する複数のリーフメッシュ要素それぞれについての共回転系変形計算を、各リーフノードの先祖ノードである対象ノードに対応する前記回転成分抽出対象メッシュ要素の回転成分を用いて行う変形計算部と、前記回転成分抽出対象メッシュ要素の前記回転成分を抽出する回転抽出部と、前記対象ノード情報が示す前記対象ノードを動的に変更する対象ノード変更部と、を備え、前記対象ノード変更部は、回転成分抽出対象メッシュ要素から抽出される回転成分の大きさの評価を行い、その評価結果に基づいて、当該回転成分抽出対象メッシュ要素に対応する対象ノードの子ノード又は親ノードを新たな対象ノードに設定する非線形変形計算システムである。
他の観点からみた本発明は、手術対象部位を示す変形モデルの変形計算を、前記非線形変形計算システムによって行う手術シミュレータである。

本発明によれば、木構造で表現される多重解像度の変形モデルの共回転系変形計算を行う場合において、木構造の適切な階層において回転成分を抽出することができる。

手術シミュレータの構成図である。 コンピュータの機能ブロック図である。 初期立方体とそれを分割したルート四面体を示す図である。 四面体の二分割二統合の説明図である。 変形モデルの二分木構造図である。 変形モデルの局所領域における変形を示す図である。 共回転系変形計算の説明図である。 リーフノードと回転成分抽出対象ノードの関係を示す説明図である。 対象ノード変更処理のフローチャートである。 第１隣接メッシュ情報と第２隣接メッシュ情報のデータ構造図である。

以下、本発明の好ましい実施形態について添付図面を参照しながら説明する。
［１．手術シミュレータの概要］
図１は、実施形態に係る手術シミュレータ（非線形変形計算システム）１を示している。この手術シミュレータ１は、内視鏡下の低侵襲手術のシミュレートを行うものである。手術シミュレータ１は、遠隔地にある複数の拠点それぞれに設置された複数（図１では２つ）の訓練システム１０Ａ，１０Ｂによって構成され、遠隔地にいる複数のユーザ（例えば、指導医と指導を受ける若手医師）それぞれが、仮想空間上の臓器モデル（手術対象部位）に対して、触覚を伴った協働作業行うことで、手術訓練を行うことができる。

各訓練システム１０Ａ，１０Ｂは、それぞれ、ネットワーク２を互いに介して接続されたコンピュータ３Ａ，３Ｂ、及び、コンピュータ３Ａ，３Ｂに接続された力覚入出力デバイス４Ａ，４Ｂを備えている。各拠点のコンピュータ３Ａ，３Ｂそれぞれは、同一の３次元手術対象部位モデル（３次元変形モデル；臓器モデル）を有し、各拠点で共有される当該モデルの変形のための演算処理を行う。

力覚入出力デバイス４Ａ，４Ｂは、手術において医師が使用する鉗子・電気メスなどの手術器具としての役割を手術シミュレータ１において果たすものであり、各拠点のユーザが力覚入出力デバイス４Ａ，４Ｂに設けられたスタイラス４Ａ−１，４Ｂ−１を操作することで、手術器具を介して現実の臓器に対する処置を行うのと同様の操作を、仮想空間上の臓器モデルに対して行うことができる。また、力覚入出力デバイス４Ａ，４Ｂは、各ユーザが力覚入出力デバイス４Ａ，４Ｂを介して仮想空間上の臓器モデルに操作（インタラクション）した際の力覚を各ユーザに提示することができる。

図２に示すように、各コンピュータ３Ａ，３Ｂは、通信処理部３１、物理計算処理部３２、力覚処理部３３、描画処理部３４、及び記憶部３５を備えている。
通信処理部３１は、手術シミュレータ１を構成する他のコンピュータの通信処理部３１と通信するためのものである。通信処理部３１は、各ユーザが力覚入出力デバイス４Ａ，４Ｂを操作することで得られた操作パラメータなどを送受信する。操作パラメータには、スタイラス４Ａ−１，４Ｂ−１の位置・姿勢情報のほか、当該位置・姿勢情報のタイムスタンプなどが含まれる。通信処理部３１によって遠隔地の力覚入出力デバイス４Ａ，４Ｂによって得られた操作パラメータも受信できるため、遠隔地の力覚入出力デバイス４Ａ，４Ｂによる操作も臓器モデルの変形に反映させることができる。

物理計算処理部３２は、４面体メッシュ要素からなる臓器（手術対象部位）のボリュームモデル（臓器モデル；変形モデル）の変形計算などの物理計算を行う。
物理計算処理部３２は、複数の力覚入出力デバイス４Ａ，４Ｂからの操作パラメータに基づいてスタイラス４Ａ−１，４Ｂ−１によって操作される仮想ツール（仮想空間における鉗子・電気メスなどの手術ツール）が、変形モデルの表面に接触しているか否かを判定し、もし接触しているならば、変形モデルの表面を共有している４面体メッシュ要素に対して、仮想ツールから変形モデルに与えている力の算出を行う。変形計算では、仮想ツールの衝突時に変形モデルにかかる力を基に、例えば質点バネダンパモデルを用いて、各メッシュ要素にかかる力の算出を行う。
物理計算処理部３２は、変形計算部３２ａ、回転抽出部３２ｂ、及び対象ノード変更部３２ｃを備えている。これらの各部３２ａ，３２ｂ，３２ｃについては後述する。

力覚処理部３３は、物理計算部３２によって計算された変形モデルにかかる力を、力覚入出力デバイス４Ａ，４Ｂに出力する。これにより、臓器モデルの変形に伴う力覚をユーザへフィードバックすることができる。
描画処理部３４は、物理計算部３２によって計算された変形モデルの変形状態を示す画像を生成し、モニタ３６へ出力する。これにより、ユーザは、仮想空間上の変形モデル（臓器モデル）を視覚的に認識することができる。なお、描画処理部３４は、変形モデルに加えて、スタイラス４Ａ−１，４Ｂ−１によって操作される仮想ツールを示す画像も、変形モデルの画像とともに生成する。すなわち、モニタ３６には、内視鏡手術において内視鏡から得られる映像と同様の仮想空間映像を表示することができる。

記憶部３５には、変形モデルを構成する四面体メッシュ要素の木構造を示す木構造情報を記憶する木構造情報記憶部（第１記憶部）３５ａのほか、対象ノード情報を記憶する対象ノード情報記憶部（第２記憶部）３５ｂ、第１隣接メッシュ情報を記憶する第１隣接メッシュ情報記憶部３５ｃ、第２隣接メッシュ情報を記憶する第２隣接メッシュ情報記憶部（第３記憶部）３５ｄを備えている。これらの各記憶部３５ａ，３５ｂ，３５ｃ，３５ｄについては後述する。
なお、記憶部３５には、コンピュータ３Ａ，３Ｂに、手術シミュレータ１としての処理を実行させるためのコンピュータプログラムも格納されている。このコンピュータプログラムがコンピュータ３Ａ，３Ｂにおいて実行されることで、コンピュータ３Ａ，３Ｂは、手術シミュレータにおける必要な機能を発揮する。

［２．変形モデル］
［２．１ 多重解像度モデルとオンラインリメッシュ］
力覚を伴う手術シミュレーションでは、臓器のような柔軟物体と手術ツールとのインタラクションを高精度かつリアルタイムに再現し表示する必要がある。柔軟物体の変形モデルとしては、質点バネダンパモデル（ＭＳＤＭ）と有限要素法（ＦＥＭ）とを用いたモデルが代表的である。
本発明者らは、高い精度と計算コスト低減とを両立するため、必要に応じて複数の精度のモデルを使い分ける適応的四面体を用いた多重解像度モデルとオンラインリメッシュを提案しており（非特許文献２参照）、本実施形態においても、適応的四面体を用いた多重解像度モデルとオンラインリメッシュが採用されている。

多重解像度モデルは以下の手順にて作成される。まず、臓器などの手術対象部位を示すボリュームデータが持つ空間全体を初期立方体（図３（ａ）参照）の格子で区切る。ボリュームデータは、ＣＴやＭＲＩで撮影された断層画像から得られる。それぞれの初期立方体（初期立方体メッシュ要素）は、図３（ｂ）に示すように、６等分され、６つ同形の四面体メッシュ要素（ルート四面体）から構成される。

四面体メッシュ要素（以下、単に「四面体」という）は、ボリュームデータの複雑さに応じて、再帰的に二分割される。四面体メッシュ要素の二分割は、四面体の最長辺（基辺）の中点において行われる。
四面体の再帰的二分割の過程では、図４に示すように、親四面体の分割レベルを３ｋとすると、分割レベルが３ｋ，３ｋ＋１，３ｋ＋２，３（ｋ＋１）と進むごとに、四面体の形状は、ＴＹＰＥ０，ＴＹＰＥ１，ＴＹＰＥ２、さらに、（分割レベル＝３ｋのＴＹＰＥ０の各辺を半分にした）ＴＹＰＥ０の順に循環する。なお、初期立方体を６等分した四面体は、ＴＹＰＥ０である。
二分割は、ボリュームデータの等値面の曲率などの局所的特徴の変化に応じて、指定された近似精度を満たすまでに再帰的に行われる。より具体的には、ある四面体を二分割するか否かは、その四面体とその近傍領域の局所的特徴の変化に基づいて決定され、もし、その四面体と近傍領域内の近似精度が、指定の近似精度に達していなければ、その四面体は二分割される。

このような再帰的二分割によって、局所的特徴の変化の大きい箇所（モデルの境界部分など）は、より小さな四面体によって構成されるとともに、局所的特徴の変化の小さい箇所は、より大きな四面体によって構成される。すなわち、メッシュ要素の解像度（分割レベル）がモデルの部分によって異なる多重解像度の変形モデルが得られる。

上記のようにして得られた多重解像度の変形モデルは、図５に示すように、二分木階層構造を用いて表現される。二分木構造の各ノードが四面体に対応する。図５に示す二分木構造おいて、分割レベルがＬｅｖｅ１のノード（ルートノード）が、初期立方体を構成するルート四面体（ＴＹＰＥ０）の一つであるとすると、Ｌｅｖｅｌ２の２つの子ノードは、ルート四面体を二分割した２つの四面体（ＴＹＰＥ１）を示している。Ｌｅｖｅｌ３の孫ノードは、Ｌｅｖｅｌ２の四面体（ＴＹＰＥ１）をさらに二分割した四面体（ＴＹＰＥ２）を示し、Ｌｅｖｅｌ４のひ孫ノードは、Ｌｅｖｅｌ３の四面体（ＴＹＰＥ２）をさらに二分割した四面体（ＴＹＰＥ０）を示し、Ｌｅｖｅｌ５（の四面体（ＴＹＰＥ０）の玄孫ノードは、Ｌｅｖｅｌ４の四面体（ＴＹＰＥ０）をさらに二分割した四面体（ＴＹＰＥ１）を示している。

図５の木構造において、黒丸はリーフノードであり、多重解像度の変形モデルは、各リーフノードに対応する四面体（リーフメッシュ要素；リーフ四面体）によって構成される。変形モデルを構成する四面体は、図５のリーフノードのレベルが異なるレベルに分散していることからも明らかなように、局所的に分割レベル（解像度）が異なっている。
なお、変形モデルは、多数のルート四面体を再帰的に二分割することによって構成されているため、一つの変形モデルは、複数の二分木構造の集合によって表現される。

変形モデルの木構造を示す木構造情報は、記憶部３５の木構造情報記憶部３５ａに記憶されている。木構造情報には、各ノードが、リーフノードであるか否かを示すリーフノードフラグが含まれている。

オンラインリメッシュが採用されている本実施形態では、物理計算処理部３２の変形計算処理部３２ａが、変形モデルに対する仮想ツールによる接触（インタラクション）などによって動的に生じる各四面体の歪の大きさに従って、変形操作が与えられた局所領域の四面体を再帰的に二分割又は二統合する。
例えば、図６に示すように、仮想ツール（鉗子など）によって引っ張り操作が与えられた変形モデルの局所領域のうち、仮想ツール側（図６に示す変形モデルの左上側）に近いほど、変形が大きい（歪が大きい）ため、二分割が多く行われて細分化され、より細かい四面体（例えば、分割Ｌｅｖｅｌ５の四面体）によって構成されることになる。逆に、引っ張り操作が解除されると、変形が小さいくなる（歪が小さくなる）ため、四面体の二統合が行われて、より粗い四面体（例えば、分割Ｌｅｖｅｌ３の四面体）によって構成されることになる。
このようにオンラインリメッシュでは、変形モデルに対して動的に生じる変形に応じて、メッシュ要素である四面体の二分割又は二統合が動的に行われるため、計算量を抑えつつ精度良い変形モデルが得られる。

［２．２ 共回転系変形モデル］
本実施形態では、変形モデルの変形計算の高速化を図るため、変形モデルとして共回転系変形モデルが採用されている（非特許文献２参照）。共回転系変形モデルは、線形有限要素法と座標変換を用いた手法であり、幾何学的非線形性を考慮可能である。本実施形態の変形計算部３２ａは、変形モデルを構成する四面体それぞれについて共回転系変形計算を行う。

図７は、共回転系変形計算の概要を示している。
図７のステップＳ１中の（ａ）は、バネＳの初期状態を示しており、このバネＳはその左端が固定端であり、右端が自由端である片持ち梁状態で存在しているものとする。この初期状態から、図７のステップ１中の（ｂ）に示すように、バネＳに回転を伴う変形が生じたものとする。
この場合、まず、変形状態にあるバネＳの回転成分（回転行列）Ｒを算出し、変形状態にあるバネＳから回転成分（回転行列）Ｒを除去し、図７のステップＳ２の（ｂ）に示すように、回転成分除去後の変形状態を求める。回転成分が除去されたバネＳの変形状態からバネＳの変位ｘと力（バネＳの回復力）Ｆとが、線形有限要素法によって計算される。

そして、図７のステップＳ３の（ｂ）に示すように、回転成分除去後の変形状態にあるバネＳを、抽出された回転成分に基づいて座標変換し、元の座標系に戻す。これにより、回転を伴う変形が生じたバネＳの弾性力ベクトルと変位ベクトルとを求めることができる。

本実施形態の変形モデルの場合、メッシュ要素として四面体が用いられているため、図７に示すバネＳに代えて、変形モデルを構成する複数の四面体それぞれが、共回転系変形計算の対象となる。四面体の共回転系変形計算によって、四面体の弾性力ベクトル及び四面体の４頂点それぞれの変位ベクトルなどを求めることができる。なお、線形有限要素法では、四面体の歪エネルギーなど様々な解析結果を得ることができる。

共回転系変形計算では、回転成分（回転行列）Ｒの抽出に多くの計算量を要する。このため、本実施形態のように多数の四面体からなる変形モデルにおいて、多数の四面体それぞれについて共回転系変形計算を行うと、変形モデルを構成する多数の四面体それぞれについて回転成分抽出の計算を行う必要があり、計算時間の増大を招く。
すなわち、多重解像度の変形モデルは、図５に示すような木構造において各リーフノードに対応する四面体によって構成される。そして、各リーフノードは、木構造情報記憶部３５ａに記憶された木構造情報に含まれるリーフノードフラグによって特定される。したがって、共回転系変形計算を行う変形計算部３２ａは、リーフノードフラグによって特定される多数のリーフノードに対応する四面体それぞれについて共回転系変形計算を行うことになり、多数のリーフノードに対応する多数の四面体の回転成分を抽出することで、計算時間が増大する。

そこで、本実施形態では、可能な限り、共回転系変形計算を行う対象となる四面体（リーフノード）よりも大きな四面体（図５に示すような木構造において、共回転系変形計算を行う四面体の上位階層にある四面体）を回転成分（回転行列）Ｒの抽出対象とする（図８参照）。
すなわち、共回転系変形計算に用いられる回転成分を抽出する回転抽出部３２ｂは、リーフノードフラグによって示されるリーフノードに対応する四面体の回転成分を抽出するのではなく、リーフノードフラグとは別の情報である「対象ノード情報」が示す対象ノードに対応する複数の四面体それぞれの回転成分を抽出するよう構成されている。
回転成分抽出の対象ノードは、例えば、図８において一点鎖線で囲まれた枠内のノードであり、回転成分抽出の対象ノードを示す対象ノード情報は、前述の対象ノード情報記憶部３５ｂに記憶されている。図８において、二点鎖線で囲まれた枠内のノードは、リーフノード（変形計算対象ノード）である。回転成分抽出対象のノードは、リーフノードと同じか、リーフノードよりも上位階層のノードとして設定されている。図８では、ルートノードは、１４個であるのに対し、回転成分抽出対象ノードは、８個であり、回転成分抽出計算の対象となるノード数が少なくなっており、計算時間を抑制することができる。なお、木構造においてリーフノードは、最下位のノードであるため、回転成分抽出対象ノードとして、リーフノードよりも下位のノードが設定されることはない。

そして、回転抽出部３２が、回転成分抽出対象の四面体（上位階層四面体）において抽出した回転成分（回転行列）Ｒは、回転成分抽出対象の四面体の子孫ノードである複数のリーフノードに対応する複数の四面体（変形計算対象四面体）それぞれの共回転系変形計算において用いられる回転成分Ｒとして用いられる。なお、回転成分抽出対象ノードとリーフノードとが一致している領域（図８の「一致領域」参照）においては、リーフノードに対応する四面体の回転成分が共回転系変形計算に用いられる。
このように、本実施形態においては、各リーフノードの共回転系変形計算は、可能な限り、各リーフノードの先祖ノードである対象ノードに対応する前記回転成分抽出対象メッシュ要素の回転成分を用いて行われる。

回転成分の抽出対象ノードをどの階層に設定するかによって、回転成分の計算量は変動するため、各四面体の回転に応じて、回転成分の抽出対象となる階層を適切に設定することが望まれる。
そこで、本実施形態の物理計算処理部３２は、対象ノードを動的に変更するため、対象ノード変更部３２ｃを備えている。対象ノード変更部３２ｃは、回転抽出対象の四面体にから抽出される回転成分の適切さを評価し、所定の近似精度を満たすように、回転抽出対象の階層を上位又は下位へ変更する。

図９は、対象ノード変更部３２ｃ及び回転抽出部３２ｂによる対象ノード変更処理手順を示している。まず、手術シミュレーション（変形計算）の初期状態において、対象ノード情報３５ｂに記憶されている対象ノード情報が示す対象ノードは、予め設定された初期ノードであり、初期ノードに対応する四面体が回転成分抽出対象四面体となる（ステップＳ１１）。初期ノードは、シミュレーション対象物（手術シミュレーションの場合は、軟組織）の局所特徴（内部の非一様性）に基づき、初期解像度を決める処理によって予め決定される。 そして、仮想ツールのインタラクションなどにより、変形モデルに変形が生じると、回転抽出部３２ｂは、対象ノード情報が示す対象ノード（初期状態では初期ノード）に対応する回転成分抽出対象四面体それぞれについて、当該変形に伴って生じた回転成分（回転行列）Ｒを抽出する計算を行う（ステップＳ１２）。

すると、対象ノード変更部３２ｃは、ステップＳ１２で抽出される回転成分の大きさの評価（回転成分抽出対象の適否の評価）を行う（ステップＳ１３）。各四面体から抽出される回転成分の大きさを示す評価値が、所定の閾値（第１閾値）よりも大きければ（ステップＳ１４）、その四面体は、回転成分抽出対象としては大きすぎることになるため、回転成分抽出対象ノードを、下位階層の四面体（２つの子ノード）に変更する（ステップＳ１５）。すなわち、回転成分抽出対象の四面体（対象ノード）を、その四面体の子である２分割四面体に変更する。
例えば、図８において、ルートノード５１が回転成分抽出対象ノードの一つであるとした場合において、このルートノード５１に対応する四面体から抽出される回転成分の評価値が、所定の第１閾値よりも大きければ、ルートノード５１は、対象ノードから除外され、新たにルートノード５１の２つの子ノード５２，５３が対象ノードとして設定され、対象ノード情報記憶部３５ｂに記憶される。

また、各四面体から抽出される回転成分の評価値が、所定の閾値（第２閾値；第１閾値＞第２閾値）よりも小さければ（ステップＳ１６）、その四面体は、回転成分抽出対象としては細かすぎ、上位階層の四面体で代用しても問題がないため、回転成分抽出対象ノードを、上位階層の四面体（親ノード）に変更する（ステップＳ１７）。すなわち、回転成分抽出対象の四面体（対象ノード）を、その四面体の親である四面体に変更する。
例えば、図８において、ノード５４が回転成分抽出対象ノードの一つであるとした場合において、このノード５４に対応する四面体から抽出される回転成分の評価値が、所定の第２閾値よりも小さければ、ノード５４及びその兄弟ノード（親が共通するノード）５５は、対象ノードから除外され、新たに、ノード５４の親ノード５２が対象ノードとして設定され、対象ノード記憶部３５ｂに記憶される。

ステップＳ１５及びステップＳ１６の後は、ステップＳ１２に戻り、対象ノード情報が示す対象ノードに対応する四面体における回転成分抽出が継続して行われる。なお、回転成分の評価値が、第１閾値と第２閾値との間にある場合には、対象ノードは適切であると評価され、対象ノードの変更は行われない。 なお、変形モデルに生じる変形が、回転対象成分抽出対象ノードを子ノードに変更する必要が生じるものの、親ノードに変更する必要がほとんどないような状況である場合（例えば、図６のように仮想ツールによってモデルをつまんだ場合）には、図９のフローチャートのステップＳ１６及びステップＳ１７を省略した第２の処理手順で、回転成分抽出対象ノードを変更することもできる。この第２の処理手順の場合、ステップＳ１４において、回転成分の評価値が所定の閾値（第１閾値）よりも大きい（ステップＳ１４において「ＹＥＳ」）と判断されると、回転成分抽出対象ノードが子ノードに変更されるが、ステップＳ１４において「ＮＯ」と判断されても、回転成分抽出対象ノードは変更されない。
また、変形モデルに生じる変形が、回転対象成分抽出対象ノードを親ノードに変更する必要が生じるものの、子ノードに変更する必要がほとんどないような状況である場合（例えば、図６の状態から仮想ツールによってモデルをつまんだ状態を解除して、モデルの形状が元に戻る場合）には、図９のフローチャートのステップＳ１４及びステップＳ１５を省略した第２の変形処理手順で、回転成分抽出対象ノードを変更することもできる。この第２の変形処理手順の場合、ステップＳ１３に続くステップＳ１６において、回転成分の評価値が所定の閾値（第２閾値）よりも小さい（ステップＳ１６において「ＹＥＳ」）と判断されると、回転成分抽出対象ノードが親ノードに変更される。
このように、図９の処理手順（第１の処理手順）に代えて、変形の仕方によって第２処理手順と第３処理手順を使い分けることもできる。第２処理手順と第３処理手順を使い分ける場合、ステップＳ１４の第１閾値とステップＳ１６の第２閾値は同じ値（第１閾値＝第２閾値）とすることができ、一つの閾値（一つの基準）を満たすように、ノードの選択が行われる。

本実施形態では、対象ノード変更部３２ｃを備えているため、動的に生じる変形に伴って各四面体に生じる回転成分を動的に抽出・評価し、回転成分抽出対象となる階層を変更するため、適切な階層で回転成分を抽出することができる。また、回転成分の大きさの評価は、回転成分抽出対象となる複数の四面体それぞれについて行われるため、変形モデルにおいて回転成分が大きい領域は下位階層の四面体で回転成分が計算され、回転成分が小さい領域は、上位階層の四面体で回転成分が計算される（図８参照）。したがって、変形計算と同様に、回転成分抽出においても、多重解像度となる。しかも、回転成分抽出の際は、変形計算よりは、対象となる四面体が少なくなるため、計算量を低減することができる。

［２．３ 回転成分の評価］
［２．３．１ 第１の評価方法］
第１の評価方法では、四面体から抽出される回転成分の評価は、回転成分抽出対象の四面体の回転成分と、当該回転成分抽出対象の四面体に隣接する他の回転成分抽出対象四面体の回転成分と、の差に基づいて行われる。
より具体的には、第１の評価方法では、四面体の回転成分をクォータニオンの回転量として求める。すなわち、第１の評価方法では、回転成分抽出対象の四面体におけるクォータニオンの回転量と、当該回転成分抽出対象の四面体に隣接する他の回転成分抽出対象四面体におけるクォータニオンの回転量と、の差分が求められる。
求めた差分を回転成分の評価値として扱い、当該評価値を、所定の閾値（第１閾値又は第２閾値）と比較することにより、対象ノード変更の要否を判定することができる。

［２．３．２ 第２の評価方法］
第２の評価方法においても、四面体から抽出される回転成分の評価は、回転成分抽出対象の四面体の回転成分と、当該回転成分抽出対象の四面体に隣接する他の回転成分抽出対象四面体の回転成分と、の差に基づいて行われる。
ただし、第２の評価方法では、四面体の回転成分の差を、回転行列の要素の差分の二乗和として求める。すなわち、第２の評価方法では、回転成分抽出対象の四面体の回転行列の要素それぞれと、当該回転成分抽出対象の四面体に隣接する他の回転成分抽出対象四面体の回転行列の要素それぞれと、の差分の二乗和が求められる。
求めた二乗和を回転成分の評価値として扱い、所定の閾値（第１閾値又は第２閾値）と比較することにより、対象ノード変更の要否を判定することができる。
なお、第２の評価方法で用いられる回転行列は、回転抽出部３２によって抽出された回転行列を用いることができる。
ここで、回転行列の要素の差分の二乗和ＳＳＤと回転角度θとの関係は、下記の式のようになる。
ＳＳＤ＝４（１−ｃｏｓθ）
したがって、例えば、隣接する四面体（要素）間での回転の差を３０度以下にするには、ＳＳＤ＝４（１−ｃｏｓ３０°）≒４×０．１３＝０．５２となり、回転行列の要素の差分の二乗和ＳＳＤを、０．５２以下にすればよい。

［２．３．３ 第３の評価方法］
第３の評価方法では、四面体から抽出される回転成分の評価は、回転成分抽出対象の四面体における（せん断）歪エネルギー（密度）に基づいて行われる。（せん断）歪エネルギーが大きい四面体では、比較的大きな回転が生じているものと推定でき、四面体から抽出される回転成分も大きくなる。
そこで、四面体の歪エネルギー（密度）の大きさを回転成分の評価値として扱い、所定の閾値（第１閾値又は第２閾値）と比較することにより、対象ノード変更の要否を判定することができる。
回転成分抽出対象の四面体の歪エネルギーは、回転成分抽出対象の四面体自体を対象に計算を行っても良いが、回転成分抽出対象の四面体の子孫にあたる複数のリーフ四面体（変形計算対象四面体）が存在する場合には、それらリーフ四面体について変形計算部３２ａによって有限要素法により計算される歪エネルギーを用いても良い。リーフ四面体（リーフノード）の歪エネルギーは、変形計算部３２ａによって計算されるため、回転成分抽出対象の四面体の歪エネルギーの大きさとして、回転成分抽出対象の四面体の子孫にあたる複数のリーフ四面体の歪エネルギーの値を利用することで、演算を簡素化できる。回転成分抽出対象の四面体の歪エネルギーの大きさとしては、回転成分抽出対象の四面体の子孫にあたる複数のリーフ四面体の歪エネルギー密度の平均を用いることができる。

［２．４ 隣接メッシュ情報］
ある回転成分抽出対象四面体における回転成分の評価のために、その回転成分抽出対象四面体だけでなく、その回転成分抽出対象四面体に隣接する他の回転成分抽出対象四面体を考慮する必要がある場合（例えば、第１の評価方法・第２の評価方法参照）、隣接する他の回転成分抽出対象四面体を容易に特定できるようにしておく必要がある。

ここで、変形計算の対象となるリーフ四面体については、オンラインリメッシュの処理の際に、共に分割すべき四面体（最長辺を共有する四面体）を効率よく探索するため、各リーフ四面体に隣接する他の隣接リーフ四面体を示す情報（第１隣接メッシュ情報）が必要となる。図１０（ａ）は、各リーフ四面体に隣接する他の隣接リーフ四面体を示す情報（第１隣接メッシュ情報）のデータ構造を示しており、四面体の第１面〜第４面それぞれについてその面に隣接する他のリーフ四面体を示すＩＤ（隣接メッシュＩＤ）が格納されている。この第１隣接メッシュ情報は、記憶部３５の第１隣接メッシュ情報記憶部３５ｃに記憶されている。

しかし、本実施形態では、変形計算の対象となる四面体と、回転成分抽出の対象となる四面体と、は異なるため、第１隣接メッシュ情報だけでは、回転成分抽出対象四面体に隣接する他の回転成分抽出対象四面体を特定することができない。
そこで、本実施形態では、回転成分抽出対象四面体に隣接する他の回転成分抽出対象四面体を示す情報（第２隣接メッシュ情報）が、第１隣接メッシュ情報とは別に、記憶部３５の第２隣接メッシュ情報記憶部３５ｃに記憶されている（図１０（ｂ）参照）。対象ノード変更部３２は、隣接四面体を考慮する必要がある場合、第２隣接メッシュ情報記憶部３５ｃに記憶されている第２隣接メッシュ情報を参照して、隣接四面体を特定する。
第２隣接メッシュ情報は、対象ノードが変更された場合には、それに伴って、更新される。

なお、本発明は、上記実施形態に限定されるものではなく、様々な変形が可能である。

１ 手術シミュレータ
２ ネットワーク
３Ａ，３Ｂ コンピュータ
４Ａ，４Ｂ 力覚入出力デバイス
４Ａ−１，４Ｂ−１ スタイラス
１０Ａ，１０Ｂ 訓練システム
３１通信処理部
３２ 物理計算処理部
３２ａ 変形計算部
３２ｂ 回転抽出部
３２ｃ 対象ノード変更部
３５ 記憶部
３５ａ 木構造情報記憶部
３５ｂ 対象ノード情報記憶部
３５ｃ 第１隣接メッシュ情報記憶部
３５ｄ 第２隣接メッシュ情報記憶部
３６ モニタ

Claims (6)

1. 動的に変形が生じる多重解像度の変形モデルの変形計算を行う非線形変形計算システムであって、
   前記変形モデルを構成する複数のメッシュ要素の木構造を示す木構造情報を記憶する第１記憶部と、
   前記木構造において、回転成分を抽出すべき複数の回転成分抽出対象メッシュ要素に対応する複数の対象ノードを示す対象ノード情報を記憶する第２記憶部と、
   前記木構造のリーフノードそれぞれに対応する複数のリーフメッシュ要素それぞれについての共回転系変形計算を、各リーフノードの先祖ノードである対象ノードに対応する前記回転成分抽出対象メッシュ要素の回転成分を用いて行う変形計算部と、
   前記回転成分抽出対象メッシュ要素の前記回転成分を抽出する回転抽出部と、
   前記対象ノード情報が示す前記対象ノードを動的に変更する対象ノード変更部と、
   を備え、
   前記対象ノード変更部は、回転成分抽出対象メッシュ要素から抽出される回転成分の大きさの評価を行い、その評価結果に基づいて、当該回転成分抽出対象メッシュ要素に対応する対象ノードの子ノード又は親ノードを新たな対象ノードに設定する
   非線形変形計算システム。
2. 前記評価は、前記回転成分抽出対象メッシュ要素における歪エネルギーに基づいて行われる
   請求項１記載の非線形変形計算システム。
3. 前記回転成分抽出対象メッシュ要素における歪エネルギーは、前記回転成分抽出対象メッシュ要素に対応する対象ノードの子孫ノードである複数のリーフノードに対応する複数のリーフメッシュ要素それぞれの歪エネルギーに基づいて算出される
   請求項２記載の非線形変形計算システム。
4. 前記評価は、当該回転成分抽出対象メッシュ要素の回転成分と当該回転成分抽出対象メッシュ要素に隣接する他の回転成分抽出対象メッシュ要素の回転成分との差に基づいて行われる
   請求項１記載の非線形変形計算システム。
5. 前記変形モデルが前記複数の回転成分抽出対象メッシュ要素によって構成されているとした場合における、各回転成分抽出対象メッシュそれぞれに隣接する他の回転成分抽出対象メッシュ要素を示す隣接メッシュ情報を記憶する第２記憶部を更に備える
   請求項１〜３のいずれか１項に記載の非線形変形計算システム。
6. 手術対象部位を示す変形モデルの変形計算を、請求項１〜５のいずれか１項に記載の非線形変形計算システムによって行うことを特徴とする手術シミュレータ。

Images