特開 2023-61115 関節モデリング方法及びシミュレーション方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】公開特許公報(A)

(11)【公開番号】P2023061115

(43)【公開日】2023-05-01

(54)【発明の名称】関節モデリング方法及びシミュレーション方法

(51)【国際特許分類】

   A61B 34/10 20160101AFI20230424BHJP

【ＦＩ】

A61B34/10

【審査請求】未請求

【請求項の数】4

【出願形態】ＯＬ

(21)【出願番号】P 2021170917

(22)【出願日】2021-10-19

【国等の委託研究の成果に係る記載事項】（出願人による申告）令和３年度、国立研究開発法人日本医療研究開発機構、「医療分野研究成果展開事業 先端計測分析技術・機器開発プログラム」「生体関節コンピューテーショナルモデルの患者別迅速精密生成 － 関節外科手術の術前・術中支援 －」委託研究開発、産業技術力強化法第１７条の適用を受ける特許出願

(71)【出願人】

【識別番号】305027401

【氏名又は名称】東京都公立大学法人

(74)【代理人】

【識別番号】100165179

【弁理士】

【氏名又は名称】田▲崎▼ 聡

(74)【代理人】

【識別番号】100175824

【弁理士】

【氏名又は名称】小林 淳一

(74)【代理人】

【識別番号】100152272

【弁理士】

【氏名又は名称】川越 雄一郎

(74)【代理人】

【識別番号】100181722

【弁理士】

【氏名又は名称】春田 洋孝

(72)【発明者】

【氏名】伊井 仁志

(72)【発明者】

【氏名】藤江 裕道

(57)【要約】

【課題】シミュレーション時に計算量が少ない関節のモデルの関節モデリング方法を提供する。
【解決手段】第１の骨をモデリングし、第１の骨の形状を表す第１の骨モデルを生成する第１骨モデリングステップと、前記第１の骨と接続される第２の骨をモデリングし、第２の骨モデルを生成する第２骨モデリングステップと、前記第１の骨と前記第２の骨とを接続する靭帯を非線形な弾性体としてモデリングし、靭帯モデルを生成する靭帯モデリングステップと、を有する関節モデリング方法であって、前記第２の骨モデルは、前記第２の骨の形状の表面を基準とし、表面からの距離に依存する値を返すレベルセット関数である、関節モデリング方法。
【選択図】図１

【特許請求の範囲】

【請求項1】

第１の骨をモデリングし、第１の骨の形状を表す第１の骨モデルを生成する第１骨モデリングステップと、
前記第１の骨と接続される第２の骨をモデリングし、第２の骨モデルを生成する第２骨モデリングステップと、
前記第１の骨と前記第２の骨とを接続する靭帯を非線形な弾性体としてモデリングし、靭帯モデルを生成する靭帯モデリングステップと、
を有する関節モデリング方法であって、
前記第２の骨モデルは、前記第２の骨の形状の表面を基準とし、表面からの距離に依存する値を返すレベルセット関数である、
関節モデリング方法。

【請求項2】

前記靭帯モデリングステップは、
前記靭帯を含む関節に関する動揺性データに基づいて、前記靭帯の弾性率と自然長を推定するステップを含む、
請求項１に記載の関節モデリング方法。

【請求項3】

前記靭帯モデリングステップは、
マルコフ連鎖モンテカルロ法を使用して靭帯モデルを生成する、
請求項１又は２に記載の関節モデリング方法。

【請求項4】

請求項１から３のいずれか一項に記載の関節モデリング方法により生成された関節モデルを用いて、
前記第１の骨モデルが示す前記第１の骨の位置を、前記第２の骨モデルである前記レベルセット関数に入力して得られる前記第１の骨と前記第２の骨の距離に基づいて、前記第１の骨と前記第２の骨との間の反発力を計算し、
前記靭帯モデルによって前記第１の骨と前記第２の骨の位置関係から前記第１の骨と前記第２の骨とを接続する靭帯の復元力を計算し、
前記反発力及び前記復元力に陰解法を適用することで、関節の挙動をシミュレートする、
シミュレーション方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、関節モデリング方法及びシミュレーション方法に関する。

【背景技術】

【0002】

形態及び力学特性を有した関節のモデルを作成し、モデルを解析することで、当該関節に対して行う適切な手術の条件などを解析する技術が提案されている。関節のモデルを作成するために例えば有限要素法を用いて大腿骨、脛骨、靭帯、軟骨及び半月板を弾性体又は剛体として近似する手法がある（例えば、非特許文献１及び２）。

【先行技術文献】

【非特許文献】

【0003】

【非特許文献1】Beidokhti HN, et al., 2017, The influence of ligament modelling strategies on the predictive capability of finite element models of the human knee joint, Journal of Biomechanics, 65, 1-11.

【非特許文献2】Guess TM, et al., 2010, A subject specific multibody model of the knee with menisci, Medical Engineering & Physics, 32, 505-515.

【発明の概要】

【発明が解決しようとする課題】

【0004】

しかしながら、有限要素法を用いる場合、近似したそれぞれの部位同士の相互作用を扱う計算過程は複雑になり、計算負荷が高くなる。例えば、非特許文献２において、計算時間は１００時間と見積もられている。迅速性を求める臨床において、この計算負荷が高いことは数値シミュレーションを活用することを難しくする一因になっている。
本発明の目的は、上述した課題を解決する関節モデリング方法及びシミュレーション方法を提供することにある。

【課題を解決するための手段】

【0005】

本発明の一態様は、第１の骨をモデリングし、第１の骨の形状を表す第１の骨モデルを生成する第１骨モデリングステップと、前記第１の骨と接続される第２の骨をモデリングし、第２の骨モデルを生成する第２骨モデリングステップと、前記第１の骨と前記第２の骨とを接続する靭帯を非線形な弾性体としてモデリングし、靭帯モデルを生成する靭帯モデリングステップと、を有する関節モデリング方法であって、前記第２の骨モデルは、前記第２の骨の形状の表面を基準とし、表面からの距離に依存する値を返すレベルセット関数である、関節モデリング方法である。

【発明の効果】

【0006】

本発明によれば、シミュレーション時に計算量が少ない関節のモデルを作成することができる。

【図面の簡単な説明】

【0007】

【図1】第１の実施形態に係る関節モデリング装置の構成を示す図である。

【図2】関節モデルの一例である。

【図3】靭帯モデリング部が靭帯に係るパラメータを推定する方法を示すフローチャートである。

【発明を実施するための形態】

【0008】

以下、図面を参照しながら本発明の実施形態について詳しく説明する。
〈第１の実施形態〉
図１は、第１の実施形態に係る関節モデリング装置１の構成を示す図である。

【0009】

関節モデリング装置１は、形状データ取得部１０、第１骨モデリング部１１、第２骨モデリング部１２、動揺性データ取得部１３、靭帯モデリング部１４、モデル出力部１５を備える。関節モデリング装置１は、関節の形状データを取得し、関節モデルを出力する。関節の形状データは、関節により接続される２つの骨のデータを含み、例えば脛骨及び大腿骨である。関節モデルは、関節の挙動を再現するための関数や数値を含み、関節の挙動のシミュレーションに用いられる。図２は、関節モデルの一例である。図２に示す関節モデル２０は、第１の骨モデル２１、第２の骨モデル２２及び靭帯モデル２３を含む。第１の骨モデル２１が模擬する第１の骨は例えば、脛骨である。第２の骨モデル２２が模擬する第２の骨は第１の骨と接続する骨であって、第１の骨が脛骨である場合、第２の骨は例えば大腿骨である。靭帯モデル２３は、第１の骨と第２の骨とを接続する靭帯を模擬する。靭帯モデル２３が模擬する靭帯は、複数の靭帯を含んでもよい。例えば、第１の骨が脛骨、第２の骨が大腿骨であるとき、靭帯モデル２３は、外側側副靭帯モデル２３－１、前十字靭帯モデル２３－２、後十字靭帯モデル２３－３、後斜靭帯モデル２３－４、内側側副靭帯モデル２３－５を含む。関節モデル２０の形態は、対象となる関節を撮影したＣＴ画像に基づいて作成される。例えば、関節モデル２０の形態は、対象となる関節を撮影した画像の輝度に対してしきい値を用いて二値化を行い、さらにセグメンテーションを行い分割することで一続きの形状である脛骨及び大腿骨を得ることができる。脛骨と大腿骨とが一続きの形状になる場合、脛骨と大腿骨を分割された形状が得られるまで、二値化のしきい値の調整及びセグメンテーションを繰り返し実行してもよい。

【0010】

形状データ取得部１０は、関節を構成する第１の骨および第２の骨の形状を表すデータを取得する。第１骨モデリング部１１は第１の骨の形状を表すデータを用いて第１の骨モデル２１をモデリングする。第２骨モデリング部１２は第２の骨の形状を表すデータを用いて第２の骨モデル２２をモデリングする。動揺性データ取得部１３は、関節で測定された動揺性データを取得する。靭帯モデリング部１４は動揺性データに基づいて靭帯モデル２３に係るパラメータを推定しモデリングする。モデル出力部１５は、モデリングされた関節モデル２０を出力する。

【0011】

〈レベルセット関数〉
第１骨モデリング部１１は、第１の骨モデル２１の表面の形状データをポリゴンに分割する。ここで分割されるポリゴンは、有限要素法におけるメッシュと同等のものであり、３次元空間における長さを基準として分割されたものである。
第２骨モデリング部１２は、第２の骨をレベルセット関数としてモデリングし、第２の骨モデル２２を生成する。レベルセット関数は、例えば、符号付距離関数であって、第２の骨を基準とする座標空間上の点を入力とし、入力された点と第２の骨の表面との最短距離を返す。このとき、レベルセット関数は、第２の骨の内部では正の値を返し、第２の骨の外部では負の値を返す。つまり、レベルセット関数は、入力された点が第２の骨の表面にあるときは値０を返し、第２の骨の内部に入るほど大きい値を返し、第２の骨の外部で表面から離れるほど小さい値を返す。

【0012】

レベルセット関数は、第１の骨と第２の骨の接触判定に使用される。関節モデル２０において、第１の骨モデル２１のポリゴンの位置はレベルセット関数に代入され、返り値が正の値となる、つまり第２の骨の内部にあると判定されたポリゴンは反発力が働くようにシミュレートされる。ここで反発力はレベルセット関数により算出された値に比例し、かつ、第２の骨モデル２２の表面に垂直方向に働く。

【0013】

比較例として、第２の骨モデル２２を第１の骨モデル２１と同様にポリゴンに分割し、第２の骨モデル２２のポリゴンと第１の骨モデル２１とをそれぞれ接触判定させるモデルを考える。このようなモデルでは各々の骨のポリゴンの数の積と同じ回数接触判定を行う必要があるため、計算量が大きくなる。これに対して本実施形態の関節モデル２０においては、第２の骨モデル２２がレベルセット関数としてモデリングされていることから、各々の骨の接触判定には第１の骨のポリゴンの数と同じ回数の判定で接触判定を実現でき、計算量を削減することができる。

【0014】

〈靭帯モデリング部の推定方法〉
靭帯モデリング部１４は、靭帯を非線形な弾性体としてモデリングする。靭帯モデリング部１４は、膝関節に対して行った前後方動揺性の検査結果に基づいて靭帯に係るパラメータを推定する。例えば、靭帯モデリング部１４はマルコフ連鎖モンテカルロ法を使用することで、靭帯に係るパラメータである剛性及び自然長を推定する。

【0015】

図３は、靭帯モデリング部１４が靭帯に係るパラメータを推定する方法を示すフローチャートである。靭帯モデリング部１４は、初めに靭帯に係るパラメータkを設定する（ステップＳ３１）。パラメータkは例えば靭帯の剛性及び自然長を含むベクトルである。その後、靭帯モデリング部１４は、パラメータ候補k*を作成する（ステップＳ３２）。パラメータ候補k*は例えばパラメータkにランダムに与えられる背景雑音wを加えた値である。靭帯モデリング部１４は、パラメータk及びパラメータ候補k*に基づいて尤度関数を算出する（ステップＳ３３）。尤度関数は、パラメータk及びパラメータ候補k*が靭帯のパラメータとして設定されたときに膝関節に加えられる力とその変位を、動揺性データから見たときの尤もらしさを示す関数である。パラメータk及びパラメータ候補k*に基づき尤度関数を使用して算出される値をそれぞれp及びp*とする。kよりもk*の方がより尤もらしいパラメータであるとき、pよりもp*の方が値は大きくなる。

【0016】

靭帯モデリング部１４は遷移確率ζを算出する（ステップＳ３４）。ζは、例えば1とp*/pのうち、小さい値をとる。靭帯モデリング部１４は、パラメータkを更新する（ステップＳ３５）。靭帯モデリング部１４は、例えばζが０から１の一様乱数よりも大きいときにkをk*に更新し、それ以外の場合ではkを更新しない。靭帯モデリング部１４は、以上の動作を所定の回数繰り返した場合（ステップＳ３６：ＹＥＳ）、動作を終了し、所定の回数に達していない場合（ステップＳ３６：ＮＯ）、ステップＳ３２からの動作を再度行う。

【0017】

靭帯モデリング部１４が動作を終了する条件は回数を基準としなくてもよい。例えば、回数を基準として動作を終了する場合と異なりパラメータのばらつきを評価できない可能性があるものの、靭帯モデリング部１４は、尤度関数pとp*との値の差が一定値以下になったときに動作を終了してもよい。

【0018】

以上の手順を行うことで、靭帯モデリング部１４は、靭帯に係るパラメータを算出する。最適化に基づいて靭帯に係るパラメータを算出する際には、パラメータの解が局所最適解に陥る可能性があり、また、パラメータの感度に起因する解析結果のバラつきを評価することができなかった。本実施形態に係る靭帯モデリング部１４は、マルコフ連鎖モンテカルト法を使用することで、パラメータの解が局所最適解に陥ることがなく、パラメータが変化したことによる解析結果のバラつきを評価することができる。
また、従来、靭帯に係るパラメータを算出する際には、靭帯の剛性を膝関節の屈曲角度により変化させるという非物理的な仮定が採用されることがあった。しかしながら、本実施形態に係る靭帯モデリング部１４は、靭帯の剛性を屈曲角度に関わらず算出することで、より実際の膝関節に適合したモデルを作成することができる。

【0019】

〈関節モデルのシミュレーション〉
脛骨の表面形状の挙動は、例えば関節モデル２０の第１の骨モデル２１の表面の各々のポリゴンを有限要素法に基づく陽解法を適用することで、コンピュータであるシミュレータを用いてシミュレートすることができる。このとき、シミュレータは、各々のポリゴンにおいて、レベルセット関数による反発力及び靭帯モデル２３のパラメータや長さに依存する復元力を算出することで、関節の挙動をシミュレートする。具体的には、シミュレータは、第１の骨モデル２１が示す第１の骨の表面の各ポリゴンの位置を、第２の骨モデル２２であるレベルセット関数に入力することで、各ポリゴンと第２の骨の表面との距離を求める。シミュレータは、各ポリゴンと第２の骨の表面との距離から、第１の骨に生じる反発力の大きさと方向を計算する。次に、シミュレータは、靭帯モデル２３によって第１の骨と第２の骨の位置関係から靭帯の復元力を計算する。
次に、シミュレータは、第１の骨モデル２１の各々のポリゴンについて、陰解法によって、変位に係るパラメータを解く。陰解法を使用すると陽解法を使用するよりも短い時間で計算をすることができる。

【0020】

陰解法においては、次の２式（１）（２）に基づいて計算することができる。

【0021】

【数1】

【数2】

【0022】

ここで、ηは減衰係数、Ｖは速度、Ｆは力、Ｋは減衰行列、ωは角速度、Ｔはトルクである。また、上付き文字のｎやｎ＋１はステップｎ及びステップｎ＋１における各々の値を示す。式（１）（２）はまとめて式（３）のように書くことができる。

【0023】

【数3】

【0024】

ここで、Ａ^ｎ、Ｂ^ｎ、Ｃ^ｎ、Ｄ^ｎはそれぞれ、並進速度及び角速度に関する力及びトルクに対するヤコビアン行列である。式（３）によりＶ^ｎ＋１及びω^ｎ＋１を算出することで、変位Ｘ^ｎ＋１や姿勢クォータニオンＱ^ｎ＋１を算出することができる。

【0025】

〈他の実施形態〉
以上、図面を参照してこの発明の一実施形態について詳しく説明してきたが、具体的な構成は上述のものに限られることはなく、この発明の要旨を逸脱しない範囲内において様々な設計変更等をすることが可能である。
例えば、上述した実施形態では膝関節のモデリングについて特に説明したが、これに限られない。例えば、関節モデリング装置１は、肘関節など他の関節においても同様にモデリングを行うことができる。

【0026】

上述した実施形態における関節モデリング装置１の一部又は全部をコンピュータで実現するようにしてもよい。その場合、この機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することによって実現してもよい。なお、ここでいう「コンピュータシステム」とは、ＯＳや周辺機器のハードウェアを含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ－ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記録装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間の間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものを含んでもよい。また上記プログラムは、前述した機能の一部を実現するためのものであってもよく、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであってもよい。また、関節モデリング装置１の一部または全部は、ＦＰＧＡ（Field Programmable Gate Array）等のプログラマブルロジックデバイスを用いて実現されるものであってもよい。

【符号の説明】

【0027】

１ 関節モデリング装置、１０ 形状データ取得部、１１ 第１骨モデリング部、１２ 第２骨モデリング部、１３ 動揺性データ取得部、１４ 靭帯モデリング部、１５ モデル出力部、２０ 関節モデル、２１ 第１の骨モデル、２２ 第２の骨モデル、２３ 靭帯モデル

【図1】

【図2】

【図3】