TOP > 国内特許検索 > 情報処理装置、情報処理方法及びプログラム > 明細書

明細書 :情報処理装置、情報処理方法及びプログラム

発行国 日本国特許庁(JP)
公報種別 再公表特許(A1)
発行日 令和元年11月7日(2019.11.7)
発明の名称または考案の名称 情報処理装置、情報処理方法及びプログラム
国際特許分類 G06T   7/00        (2017.01)
A61B   5/055       (2006.01)
FI G06T 7/00 612
A61B 5/055 380
国際予備審査の請求 未請求
全頁数 26
出願番号 特願2018-559002 (P2018-559002)
公序良俗違反の表示 1.WINDOWS
国際出願番号 PCT/JP2017/044527
国際公開番号 WO2018/123559
国際出願日 平成29年12月12日(2017.12.12)
国際公開日 平成30年7月5日(2018.7.5)
優先権出願番号 2016254752
優先日 平成28年12月28日(2016.12.28)
優先権主張国 日本国(JP)
指定国 AP(BW , GH , GM , KE , LR , LS , MW , MZ , NA , RW , SD , SL , ST , SZ , TZ , UG , ZM , ZW) , EA(AM , AZ , BY , KG , KZ , RU , TJ , TM) , EP(AL , AT , BE , BG , CH , CY , CZ , DE , DK , EE , ES , FI , FR , GB , GR , HR , HU , IE , IS , IT , LT , LU , LV , MC , MK , MT , NL , NO , PL , PT , RO , RS , SE , SI , SK , SM , TR) , OA(BF , BJ , CF , CG , CI , CM , GA , GN , GQ , GW , KM , ML , MR , NE , SN , TD , TG) , AE , AG , AL , AM , AO , AT , AU , AZ , BA , BB , BG , BH , BN , BR , BW , BY , BZ , CA , CH , CL , CN , CO , CR , CU , CZ , DE , DJ , DK , DM , DO , DZ , EC , EE , EG , ES , FI , GB , GD , GE , GH , GM , GT , HN , HR , HU , ID , IL , IN , IR , IS , JO , JP , KE , KG , KH , KN , KP , KR , KW , KZ , LA , LC , LK , LR , LS , LU , LY , MA , MD , ME , MG , MK , MN , MW , MX , MY , MZ , NA , NG , NI , NO , NZ , OM , PA , PE , PG , PH , PL , PT , QA , RO , RS , RU , RW , SA , SC , SD , SE , SG , SK , SL , SM , ST , SV , SY , TH , TJ , TM , TN , TR , TT
発明者または考案者 【氏名】中尾 恵
【氏名】松田 哲也
【氏名】森田 充樹
出願人 【識別番号】504132272
【氏名又は名称】国立大学法人京都大学
個別代理人の代理人 【識別番号】100104880、【弁理士】、【氏名又は名称】古部 次郎
【識別番号】100114546、【弁理士】、【氏名又は名称】頭師 教文
審査請求 未請求
テーマコード 4C096
5L096
Fターム 4C096AA18
4C096AB38
4C096DC21
4C096DC35
4C096DC40
5L096AA06
5L096AA09
5L096BA06
5L096BA13
5L096DA02
5L096HA09
要約 情報処理装置は、弾性体の一部分を撮像した画像から取得した観察変位量と弾性体の三次元モデルについて算出される推定変位量のうち観察領域に対応する領域部分の推定変位量との差分を最小化する最適解として、既知の作用点に作用する外力の大きさと当該三次元モデルにおける弾性率の分布とを求める演算部を有する。
特許請求の範囲 【請求項1】
弾性体の一部分を撮像した画像から取得した観察変位量と当該弾性体の三次元モデルについて算出される推定変位量のうち観察領域に対応する領域部分の推定変位量との差分を最小化する最適解として、既知の作用点に作用する外力の大きさと当該三次元モデルにおける弾性率の分布とを求める演算部
を有する情報処理装置。
【請求項2】
前記演算部は、次式に基づいて、前記外力の大きさと、前記弾性率の分布を求めることを特徴とする請求項1に記載の情報処理装置。
JP2018123559A1_000015t.gif ここで、
は、評価関数J(E,f)を最小にするEとfの集合であり、
Eは、弾性率であり、
は、弾性率の基準値であり、
fは、外力の大きさであり、
は、観察変位量であり、
(E)は、弾性率Eで表現される全体剛性マトリクスKの逆行列成分のうち観察できる部分であり、
(E)fは、推定変位量であり、
λは、係数であり、
|| ||は、1次元のノルムであり、
|| ||は、2次元のノルムである。
【請求項3】
弾性体の一部分を撮像した画像から観測変位量を取得する取得部と、
前記弾性体に作用する外力の作用点を受け付ける受付部と、
前記弾性体の三次元モデルにおける弾性率の分布と、前記外力の大きさとを変数として、当該三次元モデルの推定変位量を推定する推定部と、
前記観測変位量と当該観測変位量に対応する前記推定変位量との差分を最小化する最適解として、前記外力の大きさと前記三次元モデルにおける前記弾性率の分布を求める最適化部と
を有する情報処理装置。
【請求項4】
前記最適化部は、次式に基づいて、前記外力の大きさと、前記弾性率の分布を求めることを特徴とする請求項3に記載の情報処理装置。
JP2018123559A1_000016t.gif ここで、
は、評価関数J(E,f)を最小にするEとfの集合であり、
Eは、弾性率であり、
は、弾性率の基準値であり、
fは、外力の大きさであり、
は、観察変位量であり、
(E)は、弾性率Eで表現される全体剛性マトリクスKの逆行列成分のうち観察できる部分であり、
(E)fは、推定変位量であり、
λは、係数であり、
|| ||は、1次元のノルムであり、
|| ||は、2次元のノルムである。
【請求項5】
最適解として求められた前記弾性率の分布に基づいて前記三次元モデルを表現する表現部を更に有すること
を特徴とする請求項1又は3に記載の情報処理装置。
【請求項6】
前記画像は、単一の撮像カメラによって撮像された静止画像であること
を特徴とする請求項1又は3に記載の情報処理装置。
【請求項7】
前記最適解は、前記弾性率の分布の初期値に対する変化量の絶対値和と前記差分との総和を最小化する値として求められること
を特徴とする請求項1又は3に記載の情報処理装置。
【請求項8】
前記作用点の情報が複数与えられる場合、前記最適解は、個々の当該作用点について計算される前記差分の総和を最小化する値として求めること
を特徴とする請求項1又は3に記載の情報処理装置。
【請求項9】
前記作用点に対して複数の前記外力が与えられる場合、前記最適解は、個々の当該外力について計算される前記差分の総和を最小化する値として求めること
を特徴とする請求項1又は3に記載の情報処理装置。
【請求項10】
前記作用点は複数であり、当該作用点にはそれぞれ異なる前記外力が作用されること
を特徴とする請求項1又は3に記載の情報処理装置。
【請求項11】
前記最適解を求める空間分解能を低い状態から徐々に高めること
を特徴とする請求項1又は3に記載の情報処理装置。
【請求項12】
前記作用点に作用する前記外力の大きさが与えられる場合、前記最適解は、前記差分を最小化する前記弾性率の分布として求められること
を特徴とする請求項1又は3に記載の情報処理装置。
【請求項13】
弾性体の一部分を撮像した画像から観察変位量を取得する処理と、
前記弾性体に作用する外力の作用点を受け付ける処理と、
前記弾性体の三次元モデルにおける弾性率の分布と、前記外力の大きさとを変数として、当該三次元モデルの推定変位量を推定する処理と、
前記観測変位量と当該観測変位量に対応する前記推定変位量との差分を最小化する最適解として、前記外力の大きさと前記三次元モデルにおける前記弾性率の分布を求める処理と
を有する情報処理方法。
【請求項14】
コンピュータに、
弾性体の一部分を撮像した画像から観察変位量を取得する処理と、
前記弾性体に作用する外力の作用点を受け付ける処理と、
前記弾性体の三次元モデルにおける弾性率の分布と、前記外力の大きさとを変数として、当該三次元モデルの推定変位量を推定する処理と、
前記観測変位量と当該観測変位量に対応する前記推定変位量との差分を最小化する最適解として、前記外力の大きさと前記三次元モデルにおける前記弾性率の分布を求める処理と
を実行させるプログラム。
発明の詳細な説明 【技術分野】
【0001】
本発明は、情報処理装置、情報処理方法及びプログラムに関する。
【背景技術】
【0002】
癌病変に対する触診にみられるように、生体組織の硬さに関する情報は、医師が診断を下す上で有用である。近年では、臓器などの生体組織の弾性情報を定量的かつ客観的に評価し、その分布を画像表示するエラストグラフィ技術が実用化されている。エラストグラフィ技術には、例えば核磁気共鳴弾性計測法(Magnetic Resonance Elastography: MRE)や超音波エラストグラフィがある。
【0003】
MREは、磁気共鳴画像(Magnetic Resonanse Imaging: MRI)を用いる弾性率測定法であり、計測対象を外部から振動させた状態で撮影し、対象物体内を伝播する振動波を画像化する手法である。超音波エラストグラフィは、計測対象に外圧を加えたときに生じる組織のひずみ分布を利用する方法や対象組織内にせん断波を生じさせ、その伝播速度を利用する方法に代表される手法である。
【先行技術文献】
【0004】

【非特許文献1】O. Goksel, H. Eskandari, and S. E. Salcudean, “Mesh adaptation for improving elasticity reconstruction using the FEM inverse problem,”IEEE Trans. Med. Imaging、 vol. 32、no. 2, pp. 408-418, 2013.p21-23
【発明の概要】
【発明が解決しようとする課題】
【0005】
ところで、エラストグラフィ技術により弾性情報を取得できる領域範囲は、波が伝播する局所領域(対象組織の表面部分)や圧力を加えた領域の近傍範囲に限られる。
【0006】
そこで、本発明は、局所的な変位情報から弾性体全体における弾性率の分布を推測する手法の提供を目的とする。
【課題を解決するための手段】
【0007】
請求項1に記載の発明は、弾性体の一部分を撮像した画像から取得した観察変位量と当該弾性体の三次元モデルについて算出される推定変位量のうち観察領域に対応する領域部分の推定変位量との差分を最小化する最適解として、既知の作用点に作用する外力の大きさと当該三次元モデルにおける弾性率の分布とを求める演算部を有する情報処理装置である。
請求項2に記載の発明は、前記演算部は、次式に基づいて、前記外力の大きさと、前記弾性率の分布を求めることを特徴とする請求項1に記載の情報処理装置である。
JP2018123559A1_000003t.gif ここで、Eは、評価関数J(E,f)を最小にするEとfの集合であり、Eは、弾性率であり、Eは、弾性率の基準値であり、fは、外力の大きさであり、uは、観察変位量であり、L(E)は、弾性率Eで表現される全体剛性マトリクスKの逆行列成分のうち観察できる部分であり、L(E)fは、推定変位量であり、λは、係数であり、|| ||は、1次元のノルムであり、|| ||は、2次元のノルムである。
請求項3に記載の発明は、弾性体の一部分を撮像した画像から観測変位量を取得する取得部と、前記弾性体に作用する外力の作用点を受け付ける受付部と、前記弾性体の三次元モデルにおける弾性率の分布と、前記外力の大きさとを変数として、当該三次元モデルの推定変位量を推定する推定部と、前記観測変位量と当該観測変位量に対応する前記推定変位量との差分を最小化する最適解として、前記外力の大きさと前記三次元モデルにおける前記弾性率の分布を求める最適化部とを有する情報処理装置である。
請求項4に記載の発明は、前記最適化部は、次式に基づいて、前記外力の大きさと、前記弾性率の分布を求めることを特徴とする請求項3に記載の情報処理装置である。
JP2018123559A1_000004t.gif ここで、Eは、評価関数J(E,f)を最小にするEとfの集合であり、Eは、弾性率であり、Eは、弾性率の基準値であり、fは、外力の大きさであり、uは、観察変位量であり、L(E)は、弾性率Eで表現される全体剛性マトリクスKの逆行列成分のうち観察できる部分であり、L(E)fは、推定変位量であり、λは、係数であり、|| ||は、1次元のノルムであり、|| ||は、2次元のノルムである。
請求項5に記載の発明は、最適解として求められた前記弾性率の分布に基づいて前記三次元モデルを表現する表現部を更に有することを特徴とする請求項1又は3に記載の情報処理装置である。
請求項6に記載の発明は、前記画像は、単一の撮像カメラによって撮像された静止画像であることを特徴とする請求項1又は3に記載の情報処理装置である。
請求項7に記載の発明は、前記最適解は、前記弾性率の分布の初期値に対する変化量の絶対値和と前記差分との総和を最小化する値として求められることを特徴とする請求項1又は3に記載の情報処理装置である。
請求項8に記載の発明は、前記作用点の情報が複数与えられる場合、前記最適解は、個々の当該作用点について計算される前記差分の総和を最小化する値として求めることを特徴とする請求項1又は3に記載の情報処理装置である。
請求項9に記載の発明は、前記作用点に対して複数の前記外力が与えられる場合、前記最適解は、個々の当該外力について計算される前記差分の総和を最小化する値として求めることを特徴とする請求項1又は3に記載の情報処理装置である。
請求項10に記載の発明は、前記作用点は複数であり、当該作用点にはそれぞれ異なる前記外力が作用されることを特徴とする請求項1又は3に記載の情報処理装置である。
請求項11に記載の発明は、前記最適解を求める空間分解能を低い状態から徐々に高めることを特徴とする請求項1又は3に記載の情報処理装置である。
請求項12に記載の発明は、前記作用点に作用する前記外力の大きさが与えられる場合、前記最適解は、前記差分を最小化する前記弾性率の分布として求められることを特徴とする請求項1又は3に記載の情報処理装置である。
請求項13に記載の発明は、弾性体の一部分を撮像した画像から観察変位量を取得する処理と、前記弾性体に作用する外力の作用点を受け付ける処理と、前記弾性体の三次元モデルにおける弾性率の分布と、前記外力の大きさとを変数として、当該三次元モデルの推定変位量を推定する処理と、前記観測変位量と当該観測変位量に対応する前記推定変位量との差分を最小化する最適解として、前記外力の大きさと前記三次元モデルにおける前記弾性率の分布を求める処理とを有する情報処理方法である。
請求項14に記載の発明は、コンピュータに、弾性体の一部分を撮像した画像から観察変位量を取得する処理と、前記弾性体に作用する外力の作用点を受け付ける処理と、前記弾性体の三次元モデルにおける弾性率の分布と、前記外力の大きさとを変数として、当該三次元モデルの推定変位量を推定する処理と、前記観測変位量と当該観測変位量に対応する前記推定変位量との差分を最小化する最適解として、前記外力の大きさと前記三次元モデルにおける前記弾性率の分布を求める処理とを実行させるプログラムである。
【発明の効果】
【0008】
請求項1記載の発明によれば、変位量を取得できる領域が弾性体の一部分に限られる場合でも弾性体全域における弾性率の分布を求めることができる。
請求項2記載の発明によれば、変位量を取得できる領域が弾性体の一部分に限られる場合でも弾性体全域における弾性率の分布を求めることができる。
請求項3記載の発明によれば、変位量を取得できる領域が弾性体の一部分に限られる場合でも弾性体全域における弾性率の分布を求めることができる。
請求項4記載の発明によれば、変位量を取得できる領域が弾性体の一部分に限られる場合でも弾性体全域における弾性率の分布を求めることができる。
請求項5記載の発明によれば、弾性率の分布を三次元モデル上で視覚的に確認できる。
請求項6記載の発明によれば、弾性体の観察に使用する装置構成を簡素化できる。
請求項7記載の発明によれば、局在性が高い弾性率の分布を求めることができる。
請求項8記載の発明によれば、変位量の取得に用いた外力の作用点が1つだけの場合に比して推定精度を高めることができる。
請求項9記載の発明によれば、変位量の取得に用いた外力が1つだけの場合に比して推定精度を高めることができる。
請求項10記載の発明によれば、変位量の取得に用いた作用点と外力がそれぞれ1つだけの場合に比して推定精度を高めることができる。
請求項11記載の発明によれば、時間の経過とともに精度の高い推定結果を得ることができる。
請求項12記載の発明によれば、外力の大きさが与えられない場合と比して推定精度を高めることができる。
請求項13記載の発明によれば、変位量を取得できる領域が弾性体の一部分に限られる場合でも弾性体全域における弾性率の分布を求めることができる。
請求項14記載の発明によれば、変位量を取得できる領域が弾性体の一部分に限られる場合でも弾性体全域における弾性率の分布を求めることができる。
【図面の簡単な説明】
【0009】
【図1】実施の形態で説明する発明の概念イメージを説明する図である。
【図2】情報処理システムのハードウェア構成例を示す図である。
【図3】コンピュータプログラムの実行を通じて実現される機能構成例を説明する図である。
【図4】弾性体のメッシュモデル例を示す図である。
【図5】弾性体のメッシュモデルの変形例を説明する図である。
【図6】二次元のメッシュモデルを示す図である。
【図7】要素剛性マトリクスKeの重ね合わせ処理の概念を説明する図である。
【図8】シミュレーション実験の真値を与える弾性率Eの分布例を示す図である。
【図9】観察できる部分の割合を説明する図である。
【図10】推定ブロック数と弾性率Eが推定される単位区画(ブロック)の大きさとの関係を説明する図である。
【図11】推定ブロック数Nが10個の場合、20個の場合、40個の場合の実験結果を示す図である。
【図12】実施の形態2で使用する観察条件を説明する図である。
【図13】弾性率Eの分布の推定に用いたパターンが1組(実施の形態1)の場合と、パターンが2組(実施の形態2)の場合の実験結果を示す図である。
【図14】弾性率Eの分布の推定に用いたパターンが1組(実施の形態1)の場合と、パターンが2組(実施の形態2)の場合と、パターンが2組(実施の形態2)の場合に更に弾性率Eの変化のスパース性を加えた場合(実施の形態3)の実験結果を示す図である。
【図15】実施の形態3の場合における観測頂点数の割合と観測精度の関係を説明する図である。
【発明を実施するための形態】
【0010】
以下、添付図面を参照して、本発明の実施の形態について詳細に説明する。

【0011】
<概念イメージ>
図1は、後述する実施の形態で採用する推定手法の概念イメージを説明する図である。図1に示すように、各実施の形態は、弾性体1の全体を観察できない場合(すなわち、観察できる部分と観察できない部分がある場合)を前提条件とし、各実施の形態では、観察できる部分における局所的な変位情報を用いて弾性体全体における弾性率の分布を推定する手法について記述する。

【0012】
各実施の形態では、弾性体1を変形させる外力Fの作用点Paの位置情報が既知であることを利用して推定処理を実行する。
作用点Paの位置情報は、推定処理が始まる前に処理システムに与えられればよく、作用点Paの位置は任意である。従って、図1に示すように、作用点Paの位置が観察できる部分に存在する必要はない。
作用点Paの個数についての制限もなく、必ずしも1つである必要はない。従って、作用点Paの個数は、1個でも、複数個でも、弾性体1の表面全体でも良い。

【0013】
一方で、作用点Paに作用する外力Fの大きさfを与える情報や弾性体1が周囲に固定されている部位(固定点)Pfの位置情報は、既知であっても未知であっても良い。一般には、既知である情報が増えるほど、推定精度が高くなり、推定に要する計算時間も少なく済む。

【0014】
図1では、観察できる部分と観察できない部分の説明のため、弾性体1の一部分が遮蔽物2によって覆われている様子を表しているが、弾性体1の一部分が遮蔽物2で覆われていることを必須とするものではない。例えば弾性体1を一方向のみから観察する場合、観察視野に含まれる部分が観察できる部分であり、観察視野に含まれない部分が観察できない部分となる。従って、観察方向から見て弾性体1の陰になる部分(例えば裏側や背面側)は、観察できない部分の一例である。

【0015】
また、後述する実施の形態では、弾性体1のうち観察できる部分の変位量は、部位毎に特定できるものとする。すなわち、外力Fの作用による弾性体1の変位量は既知の情報として処理システムに与えられるものとする。
変位量は、観察手段の違いにより、弾性体1の表面領域についてのみ与えられる場合と弾性体1の内部領域も含めて与えられる場合がある。一般には、変位量が特定された領域が多いほど(観察できる部分の割合が増えるほど)、弾性率の推定精度が高くなる。

【0016】
<用語>
各実施の形態における弾性率には、押込・引張方向についての硬さを表すヤング率E、せん断方向についての硬さを表す剛性率G、ポアソン比ν、体積弾性率Kなどが含まれ、後述する実施の形態では、弾性率がヤング率であるものとして使用する。
各実施の形態において、弾性体1とは、外力Fの作用により歪が生じる(変形する)一方、外力Fの作用が停止されると元の寸法に戻る性質を有する物質をいうものとする。本実施の形態では、例えば生体組織における各種の臓器などの近似的にフックの法則を用いて表現される弾性体を想定する。
もっとも、各実施の形態では、弾性体1における弾性率の分布が一様でない場合を想定する。ここで、弾性率の分布が一様でない場合とは、例えば弾性体1の少なくとも一部の弾性率が周囲の弾性率と異なることをいう。

【0017】
<実施の形態1>
<システム構成>
図2は、情報処理システム100のハードウェア構成例を示す図である。
情報処理システム100は、弾性体1(図1参照)の弾性率の分布を推定する演算処理を実行する情報処理装置10と、弾性体1の撮像に用いる少なくとも1つの撮像カメラ20と、弾性体1のメッシュモデルの作成に使用する三次元画像データ(ボリュームデータ)を記憶する三次元画像データ記憶装置30と、推定された弾性率の分布の表示等に用いる表示装置40と、ユーザが操作するマウスやポインティングデバイス等で構成される入力装置50とを有している。

【0018】
情報処理装置10は、いわゆるコンピュータであり、プログラムを実行するCPU(Central Processing Unit)11と、BIOS(Basic Input / Output System)やファームウェアなどのプログラムやデータを記憶するROM(Read Only Memory)12と、プログラムに作業エリアを与えるRAM(Random Access Memory)13と、処理対象とする画像データを記憶する記憶部14とにより構成される。記憶部14は、ハードディスク装置や半導体メモリなどの記憶装置である。CPU11によるプログラムの実行を通じて実現される機能については後述する。

【0019】
撮像カメラ20は、弾性体1の観察に使用される撮像装置である。撮像カメラ20は、弾性体1の観察に適した撮像装置であればよい。例えば撮像カメラ20には、赤外線カメラ、内視鏡カメラ、超音波イメージング装置、光超音波イメージング装置などを使用する。
赤外線カメラや内視鏡カメラを用いる場合、情報処理装置10は、弾性体1の表面に存在する特徴点や赤外線マーカーの移動を通じて各部位の変位量を特定することができる。超音波イメージング装置を用いる場合、情報処理装置10は、弾性体1の内部も含めて各部位の変位量を特定することができる。なお、内視鏡カメラの撮像素子は、CCD(Charge Coupled Device)でもCMOS(Complementary MOS)でもよい。
本実施の形態では、単一の撮像カメラ20を使用し、弾性体1を一方向から撮影した静止画像を出力する。

【0020】
三次元画像データ記憶装置30は、弾性体1を断層撮影することにより取得した断層画像データを積み重ねるように再構成した画像(再構成画像)を三次元画像データとして記憶する大容量のハードディスク装置等で構成される。
三次元画像データ記憶装置30は、情報処理装置10に対して直接接続される場合だけでなく、ネットワーク上にサーバやネットワークストレージとして存在してもよい。ここで、断層画像データの取得には、例えばCT(Computed Tomography)、PET(Positron Emission Tomography)、MRI等を使用する。

【0021】
表示装置40は、例えば液晶ディスプレイで構成される。液晶ディスプレイは、液晶パネル、バックライトなどで構成される。なお、表示装置40は、有機EL(ElectroLuminescence)ディスプレイ等でもよい

【0022】
図3は、CPU11によるコンピュータプログラムの実行を通じて実現される機能構成例を説明する図である。
処理機能の観点から見た情報処理装置10は、三次元のメッシュモデルを生成するメッシュモデル生成部111と、撮像カメラ20から出力されるカメラ画像から観察部位についての変位量(観察変位量u)を取得する観察変位量取得部112と、作用点Paの位置情報を受け付ける作用点受付部113と、作用点情報をメッシュモデルに適用してメッシュモデルの各部位における変位量(推定変位量u’)を計算する変位量推定部114と、推定変位量u’と観察変位量uとの差分を最小化する弾性率E(結果として弾性率Eの関数である要素剛性マトリクスK)と外力Fの大きさfを最適解として求める最適化部115と、最適解とする要素剛性マトリクスKから特定される弾性率の分布をメッシュモデルに重ねて表現した画像を生成して表示装置40に出力する表現部116とを有している。これらの機能部は、請求の範囲における「演算部」の一例である。

【0023】
メッシュモデル生成部111は、三次元画像データ記憶装置30から観察対象とする弾性体1の三次元画像データを取り込むと、連続体として与えられた三次元画像データから単純な形状要素の集合体であるメッシュモデルを生成する機能部である。ここでのメッシュモデルは三次元モデルの一例であり、最小要素は立方体や四面体(三角錐)で表現される。

【0024】
図4は、弾性体1をメッシュ構造で表現したメッシュモデル200の例を示す図である。図4は、外力F(図5参照)が作用点Pa(図5参照)に加えられる前のメッシュモデル200を表している。図中の丸印は最小要素のノードを示す。なお、図5は、メッシュモデル200の右端下部を作用点Paとして外力Fを加えて弾性変形させた後のメッシュモデル200を表している。

【0025】
観察変位量取得部112は、外力Fが作用点Paに加えられる前のカメラ画像と外力Fが作用点Paに加えられた後のカメラ画像とを比較して画像内に現れる各部位の観察変位量uを取得する機能部である。
なお、観察変位量取得部112で取得される観察変位量uは、カメラ画像上で特定される各位置についての変位量であり、メッシュモデル200を構成するボクセルとの位置関係は不明である。

【0026】
そこで、観察変位量取得部112には、三次元画像データ又はメッシュモデル200(図4参照)とカメラ画像をマッチングし、観察変位量uをメッシュモデル200の個々の頂点に対応付ける処理機能が用意されている。ここでの対応付けの際には、三次元画像データ又はメッシュモデル200の寸法とカメラ画像の寸法とを揃える処理も実行される。

【0027】
一般には、三次元画像データ又はメッシュモデル200に含まれる画像情報が多いほどカメラ画像とのマッチング精度が高くなる。例えばメッシュモデル200をマッチングに使用する場合、メッシュモデル200の表面に三次元画像データの表面情報をレンダリングしたレンダリング画像をカメラ画像の輝度情報と比較すれば、メッシュモデル200の頂点に高いマッチング精度で観察変位量uを対応付けることができる。

【0028】
作用点受付部113は、作用点Paの位置情報を、入力装置50を通じて受け付ける機能部である。作用点受付部113は、例えばクリック操作が検出された際に、マウスカーソルが位置していたメッシュモデル200上の部位を外力Fの作用点Paとして受け付ける。
作用点受付部113は、受け付けた作用点Paの位置を、メッシュモデル200上の頂点の識別番号で特定し、変位量推定部114に出力する。

【0029】
変位量推定部114は、線形有限要素法で使用する力の釣り合い式(すなわち、フックの法則を表現する式)を用いて、メッシュモデル200を構成する全ての頂点について推定変位量u’を計算する機能部である。
具体的には、変位量推定部114は、メッシュモデル200の各頂点に加わる外力Fの大きさをf、各頂点における推定変位量をu’、全体剛性マトリクスをKとする力の釣り合い式(f=K・u’)を変形した式1に基づいて、推定変位量u’を計算する。
JP2018123559A1_000005t.gifなお、添え字のoは、観測できる部分に位置する頂点を意味し、Lは、全体剛性マトリクスKの逆行列(K-1)を意味する。

【0030】
本実施の形態の場合、弾性体1の弾性率Eの分布と弾性体1に加わる外力Fの大きさf(ベクトル)が最終的な推定対象となる。なお、弾性率Eは、全体剛性マトリクスKの内包パラメータであるので、全体剛性マトリクスKが特定されれば、弾性率Eも特定される。

【0031】
以下では、説明を簡単にするために、全体剛性マトリクスKと弾性率Eの分布との関係を、図6を用いて説明する。図6は、二次元のメッシュモデルである。言うまでもなく、情報処理システム100で扱うのは、三次元のメッシュモデル200である。

【0032】
図6に示すメッシュモデルの場合、三角形の要素Iは頂点番号1、2及び3で規定され、三角形の要素IIは頂点番号2、3及び4で規定される。このとき、頂点番号iに作用する力f(i=1~4)と頂点番号iにおける変位量u(i=1~4)は、全体剛性マトリクスKを用いた線形有限要素法での力の釣り合い式(式2)として与えられる。
JP2018123559A1_000006t.gif

【0033】
ここで、個々の頂点番号i(i=1~4)に作用する外力Fの大きさf(i=1~4)とその変位量u(i=1~4)との関係付けに用いる要素剛性マトリクスKは、以下の式3により導出される。
JP2018123559A1_000007t.gifここで、Vは三角形の要素の体積である。
また、Bマトリクスは、要素の形状関数から求まるマトリクスであり、Dマトリクスは、弾性率Eとポアソン比νから求まるマトリクスである。

【0034】
Dマトリクスは、次式(式4)により表すことができる。
JP2018123559A1_000008t.gif式4から分かるように、Dマトリクスは弾性率Eの関数である。従って、要素剛性マトリクスKも弾性率Eの関数として表現される。

【0035】
全体剛性マトリクスKは、三角形の要素Iの要素剛性マトリクスKと、三角形の要素IIの要素剛性マトリクスKとの和で与えられる。すなわち、全体剛性マトリクスKは、頂点番号1、2、3で構成される三角形の要素Iの要素剛性マトリクスKと、頂点番号2、3、4で構成される三角形の要素IIの要素剛性マトリクスKとの和で表される。

【0036】
図6に示すように、三角形の要素Iと三角形の要素IIとでは、互いに重なり合う頂点が存在する。そこで、重なり合う各頂点に対して要素剛性マトリクスKの重ね合わせを行う。
図7は、要素剛性マトリクスKの重ね合わせ処理の概念を説明する図である。図7は、全体剛性マトリクスKが、要素Iの要素剛性マトリクスKと要素IIの要素剛性マトリクスKとの和で表されることを示している。前述したように、要素剛性マトリクスKは弾性率Eの関数であったので、その和で表される全体剛性マトリクスKも弾性率Eの関数となる。
従って、本実施の形態におけるメッシュモデル200における全体剛性マトリクスKも弾性率Eの関数となる。

【0037】
ところで、本実施の形態のカメラ画像は、前述したように、作用点Paは既知であるが大きさfが未知の外力Fによって弾性変形された弾性体1の一部分を撮像した画像である。
従って、式2のうち未知のパラメータは、弾性率Eの分布(全体剛性マトリクスK)と外力Fの大きさf(ベクトル)となる。

【0038】
そこで、変位量推定部114では、これら未知のパラメータにそれぞれ任意の値を設定して推定変位量u’を計算する。なお、推定変位量u’の初回推定時には、変位量推定部114は、未知パラメータである弾性率Eの分布(全体剛性マトリクスK)と外力Fの大きさfのそれぞれに初期値を与えて推定変位量u’を計算する。

【0039】
式2をマトリクスで表現すると、次式(式5)となる。
JP2018123559A1_000009t.gifここで、u’は観察できる部分の推定変位量であり、u’は観察できない部分の推定変位量である。また、Lは弾性率Eで表現される全体剛性マトリクスKの逆行列成分のうち観察できる部分であり、Lは弾性率Eで表現される全体剛性マトリクスKの逆行列成分のうち観察できない部分である。fは観察できる部分の頂点に作用する外力Fの大きさ、fは観察できない部分の頂点に作用する外力Fの大きさである。

【0040】
観察できる部分の推定変位量u’に着目すると、式5は、次式(式6)として表現できる。
JP2018123559A1_000010t.gifここで、作用点Paが既知であるので、外力Fの大きさfのうち作用点Pa以外の頂点に対応する外力Fの大きさを0(ゼロ)で表すことができる。
変位量推定部114は、この式6を用いて弾性率Eで表される成分Lと外力Fの大きさfにより推定変位量u’(初期値を含む)を計算する。

【0041】
最適化部115は、メッシュモデル200の全ての頂点番号(観察できる部分だけでなく観察できない部分の頂点も含む)について推定された推定変位量u’と、撮像カメラ20によって撮像された部分画像から特定された観察変位量uとの差分を最小化する全体剛性マトリクスK(すなわち弾性率Eの分布)と外力fが得られるまで、未知パラメータの組み合わせを更新する。

【0042】
具体的には、最適化部115は、推定変位量u’と観察変位量uとの差分を最小化する(すなわち、観察結果に最も近い変形を与える)未知パラメータを求める最適化問題を解く処理を実行する。当該処理を式で表すと次式(式7)となる。
JP2018123559A1_000011t.gifここでのEは、評価関数J(E,f)を最小にするEとfの集合を意味する。
本実施の形態の場合、評価関数J(E,f)は、推定変位量u’と観察変位量uの差分のLノルムで与えられる。

【0043】
最適化部115は、評価関数Jの最小値が求まるまで、未知パラメータEとfを更新し、更新された値を変位量推定部114にフィードバックする。すなわち、未知パラメータEとfの更新値を用いた変位量推定部114による推定変位量u’の再計算と、最適化部115による推定変位量u’の評価とが、評価関数Jの最小値が求まるまで繰り返し実行される。

【0044】
本実施の形態では、最適化問題を解く手法として、共分散行列適応進化戦略(Covariance Matrix Adaptation Evolution Strategy:CMA-ES)を使用した。もっとも、最適化問題を解く手法には既知の他の手法を用いてもよい。例えば勾配法、モンテカルロ法、焼きなまし法などを用いてもよい。
なお、図3では、変位量推定部114と最適化部115を別の機能部として描画しているが、一つの機能部に含めてもよい。

【0045】
表現部116(図3参照)は、最適解として特定された弾性率Eの分布と外力Fの大きさのうち弾性率Eの分布を、メッシュモデル200に重ねて表現した画像を表示装置40(図2参照)に表示する。弾性率Eの分布がメッシュモデル200と関連付けて表現されることにより、ユーザは、観察できない部分も含めた弾性体1の弾性率Eの分布を視覚的に理解できる。
なお、表現部116は、弾性率Eの違いを色や輝度の違いとして表現する。また、表現部116は、メッシュモデル200の各部に対応する弾性率Eの数値を表示してもよい。なお、表現部116は、必要に応じ、最適解として求まった外力Fの大きさfを出力する。

【0046】
<シミュレーション実験による検証>
以下では、弾性体1について作成したメッシュモデル200(すなわち、図4に示すメッシュモデル)と、弾性体1を局所的に観察して得られる観察変位量uと、弾性体1に作用する外力Fの作用点Paの情報とを用いて推定した弾性率Eの分布の精度を検証した結果について説明する。

【0047】
以下では、シミュレーション実験では、弾性体1に対応するメッシュモデル200がN個の頂点を有するものとし、そのうちのM個(N>M)の頂点の観察変位量uを用いて弾性体1全体の弾性率Eの分布を推定する場合に、観測できる頂点の数が推定精度に与える影響について説明する。

【0048】
シミュレーション実験では、Visual C/C++を用いて、実施の形態に係る手法をプログラムとして記述し、作成されたプログラムを汎用計算機(OS:Windows 10 Pro、CPU:IntelCorei7-6700K、メモリ:16GB)で実行した。
なお、シミュレーションは、直方体形状のメッシュモデル200(図5参照)について実行した。このメッシュモデル200は、全体で99個の頂点を有し、最小要素を四面体とする。

【0049】
また、図5に示すように、メッシュモデル200の左端に位置する9個の頂点を固定端Pfとし、それ以外の点を自由点とした。
さらに、図5に示すように、外力Fの作用点Paとして、メッシュモデル200の右端下部に位置する3個の頂点を指定した。これら3個の頂点には、同じ大きさfの外力Fがいずれも下向き(-z軸方向)に作用するものとした。

【0050】
メッシュモデル200の外力Fによる変形は、弾性体1の内部構造を与える弾性率Eの分布状態に依存する。ここでは、本実施の形態で提案する手法によって推定された弾性率Eの分布の推定精度を評価するために使用する真値としての内部構造について説明する。図8は、真値として使用する弾性率Eの分布をメッシュモデル200に重ねて表した図である。図8に示すメッシュモデル200は、一部分の弾性率Eだけが大きく(硬く)、他の部分の弾性率Eは同じ値を採る。図中、白丸のノードは観察できる部分に位置することを表し、黒丸のノードは観察できない部分に位置することを表している。

【0051】
図8に示すメッシュモデル200は、40個の立方体ブロックのうち右端から左方向(-x方向)に7番目と8番目の左列に位置する4個の立方体ブロックの弾性率Eが2.0[MPa]と硬く、他の部分の弾性率Eが1.0[MPa]である場合を示している。
シミュレーション実験では、図8に示すメッシュモデル200に外力Fを作用させた場合に観測されたカメラ画像と作用点Paの位置に基づいて、メッシュモデル200の全ての頂点の弾性率Eの分布を推定した。

【0052】
なお、図8に示すメッシュモデル200は、図5の説明と同じく、弾性体1のメッシュモデル200の右端下部にある3個の頂点に下方向(-z軸方向)へ10Nの外力Fを作用させて変形した状態を表している。勿論、情報処理装置10にとっては、外力Fの大きさfも未知パラメータである。

【0053】
前述したように、シミュレーション実験の目的は、観測できる頂点の数が弾性率Eの推定精度に与える影響を評価することである。
そこで、シミュレーション実験では、観察できる頂点の数(観察変位量uが与えられる頂点の数)が全体の10%の場合、20%の場合、…90%の場合について評価した。図9は、メッシュモデル200のうち観察できる頂点の数の割合を説明する図である。図9は、メッシュモデル200を横方向(x-z平面)から観察した図であり、x軸方向について10個の小区分に分割されている。このうち1つの小区分を観察する場合が全体の10%を観察する場合に相当し、2つの小区分を観察する場合が全体の20%を観察する場合に相当する。以下同様に、9つの小区分を観察する場合が全体の90%を観察する場合に相当する。

【0054】
ところで、弾性率Eの推定精度は、推定空間の次元数の影響を受けると予想される。ここで、推定空間の次元数とは、弾性率Eの推定解像度で表現される。例えば推定空間の次元数が低いことは低解像度で推定することを意味し、次元数が高いことは高解像度で推定することを意味する。換言すると、推定空間の次元数は、1つのメッシュモデル200を何個の単位区画(ブロック)に分割するかを表している。

【0055】
以下のシミュレーション実験では弾性率Eの推定解像度を表す推定ブロック数Nを10個とする場合、N=20とする場合、N=40とする場合について検討した。
図10は、推定ブロック数Nと弾性率Eが推定される単位区画(ブロック)の大きさとの関係を説明する図である。

【0056】
例えばN=10の場合は、メッシュモデル200(図9参照)を10個の単位区画(ブロック)で表現することを意味し、N=20の場合は、メッシュモデル200を20個の単位区画(ブロック)で表現することを意味し、N=40の場合は、メッシュモデル200を40個の単位区画(ブロック)で表現することを意味する。
つまり、N=40の場合の単位区画(ブロック)は、図8において説明した立方体ブロックと同じ形状になる。当然ながら推定ブロック数Nの値が大きいほど(推定される単位区画の大きさが小さいほど)空間分解能が高くなり、弾性率Eの分布を高精細に推定できる。

【0057】
続いて、推定ブロック数Nと高い推定精度を得るために必要とされる観察点数の割合との関係を説明する。後述するように、推定ブロック数Nが大きいほど(推定される単位区画の大きさが小さいほど)、高い推定精度を得るにはより多くの割合で弾性体1を観察することが必要になる。
図11は、推定ブロック数Nが10個の場合、20個の場合、40個の場合の実験結果を示す図である。図11の縦軸は真の弾性率Eと推定された弾性率E’との平均二乗誤差(Root Mean Square Error: RMSE)であり、横軸は観察できる頂点数の全体に占める割合(図9参照)である。
勿論、弾性率E’の推定に際しては、前述したように、メッシュモデル200の全ての頂点に同じ値の弾性率Eを与えた状態で最適解を求める推定処理が開始される。
JP2018123559A1_000012t.gif

【0058】
図11によれば、メッシュモデル200を10個のブロックで表現する場合(低解像度で推定する場合)には、観測できる頂点の数がメッシュモデル200全体の15%程度でも高精度で弾性率Eを推定できることが分かる。また、この場合にはブロック数が少ないので、計算コストも低く済む。
また、メッシュモデル200を20個のブロックで表現する場合(中解像度で推定する場合)には、観測できる頂点の数がメッシュモデル200全体の50%程度になると高精度で弾性率Eを推定できることが分かる。
メッシュモデル200を40個のブロックで表現する場合(高解像度で推定する場合)には、観測できる頂点の数がメッシュモデル200全体の70%程度にならなければ最小二乗誤差が十分に小さくならないことが分かる。なお、ここでの低解像度、中解像度、高解像度は、単に3つの解像度間の相対的な関係を表したものにすぎない。

【0059】
<小括>
以上の通り、実施の形態1によれば、弾性体1の一部分(例えば全体の50%)を撮像した画像から取得した観察変位量uと弾性体1のメッシュモデル200について算出される推定変位量u’のうち観察領域に対応する領域部分の推定変位量との差分(u-u’)を最小化する最適解として、既知の作用点Paに作用する外力Fの大きさfとメッシュモデル200全域における弾性率Eの分布とを求めることにより、観察変位量uを取得できる領域が弾性体1の一部分に限られる場合でも弾性体1全域における弾性率Eの分布を求めることができる。

【0060】
また、推定精度は、推定空間の解像度と観察できる頂点の割合に依存し、もし高解像度で弾性体1の弾性率Eを推定したい場合には、撮像カメラ20によって撮像される弾性体1の面積を広くする必要があり、低解像で弾性体1の弾性率Eを推定したい場合には、撮像カメラ20によって撮像される弾性体1の面積は狭くてもよい。

【0061】
<実施の形態2>
本実施の形態では、前述した実施の形態1に比して、弾性率Eの分布の推定精度を高めるための手法について説明する。
図12は、実施の形態2で使用する観察条件を説明する図である。本実施の形態では、外力Fと作用点Paの組を2組用意する。図12(A)は、外力F1と作用点Pa1の第1の組(パターン1)を示し、図12(B)は、外力F2と作用点Pa2の第2の組(パターン2)を示す。
なお、図12に示すメッシュモデル200は、図8の場合と同じく、周囲より硬い弾性率Eを有する立方体ブロックを網掛けで表現している。

【0062】
図12(A)に示す第1の組(パターン1)は、メッシュモデル200の右端下部の3つの頂点を作用点Pa1として、下向き(-z軸方向)に第1の外力F1を作用させた場合の変形を表している。図12(B)に示す第2の組(パターン2)は、メッシュモデル200の右端のうち端面に向かって左辺に位置する3つの頂点を作用点Pa2として、左向きに(-y軸方向)に第2の外力F2を作用させた場合の変形を表している。

【0063】
なお、図12においては、作用点Paだけでなく外力Fが作用する方向も異なる例を示しているが、作用点Paは1つで外力Fの作用方向が複数でもよいし、作用点Paは2つで外力Fの作用方向は1つでもよい。
この場合も、実施の形態1の場合と同様、情報処理装置10は、各パターンについての観察変位量uを求め、その後、それぞれについて観察変位量uと推定変位量u’との差分(u-u’)を最小化する外力F1の大きさf1と外力F2の大きさf2と、メッシュモデル200全域における弾性率Eの分布を最適解として求める。

【0064】
すなわち、最適化部115は、次式(式9)を満たす未知パラメータを最適解として求める。
JP2018123559A1_000013t.gifここで、Uはu1oとu2oを列方向に並べた行列であり、Fはf1とf2を列方向に並べた行列である。また、3行目の式における2次元のノルムは行列のノルムであり、例えばフロベニアスノルムなどである。

【0065】
図13は、弾性率Eの分布の推定に用いたパターンが1組(実施の形態1)の場合と、パターンが2組(実施の形態2)の場合の実験結果を示す図である。図13の縦軸は真の弾性率Eと推定された弾性率E’との平均二乗誤差(Root Mean Square Error: RMSE)であり、横軸は観察できる頂点数の割合(図9参照)である。
図13からも分かるように、弾性率Eの分布の推定に使用する作用点Paと外力Fの大きさfとで規定される組数が「2」の場合には、実施の形態1の場合(組数が「1」の場合)に比して、観察できる部分の割合が同じであっても、推定精度を高めることができる。

【0066】
なお、本実施の形態では作点Paと外力Fの大きさfとで規定される組数を「2」としているが、組数は3以上でもよい。一般的に組数が増えるほど推定精度を高めることができる。例えば超音波などの振動を加えることによって評価対象物(例えば臓器)に生じるひずみ像を非侵襲的に評価するための手法の一例であるエラストグラフィ法では、評価対象物に作用させる外力Fが時間的に変化する。すなわち、エラストグラフィ法では、作用点Paと外力Fの大きさfとで規定される組数が複数である。従って、本実施の形態は、エラストグラフィ法との相性がよく、高い精度での推定が可能である。

【0067】
<実施の形態3>
本実施の形態でも、前述した実施の形態1に比して、弾性率Eの分布の推定精度を高めるための他の手法について説明する。
ここでは、弾性率Eの変化にスパース性を推定できる場合について説明する。スパース性があるとは、要素の大部分が0(ゼロ)を持つという性質のことである。具体的には、一部分だけが硬い又は柔らかいといった局在性があることが知られている弾性体1に適用できる手法である。例えば臓器は、健常組織の弾性率Eは均一であるのに対し、一部の腫瘍部位だけは健常組織に比して硬い(弾性率Eが高い)ことが想定される。従って、臓器の弾性率Eの分布の推定にはスパース性を利用することができる。

【0068】
最適化問題に弾性率Eの変化にスパース項を導入することにより、弾性率Eの推定値の変化が局所的であるという制約を与えることができ、その分、推定精度の向上が期待できる。
本実施の形態の場合は、最適化部115は、次式(式10)を満たす未知パラメータを最適解として求める。
JP2018123559A1_000014t.gif

【0069】
式10の第2項がスパース項である。ここで、Eは、ほとんどの頂点に適用される弾性率の基準値である。従って、スパース項は、各頂点について推定された弾性値Eと基準値との差分(すなわちΔE)について1次元のノルムLを求める式となる。
本実施の形態では、係数λを0.1とした。この式の場合、弾性率Eが基準値と異なる一部の頂点についてのみ差分ΔE(=E-E)が非0(ゼロ)となり、大部分の頂点においての差分ΔEは0(ゼロ)となる。このスパース項は、推定された弾性率Eの分布の局在性が高い場合にノルムが小さくなる。従って、スパース項のノルムが大きい弾性率Eの分布は評価が低くなる。なお、係数λの値は一例である。

【0070】
図14は、弾性率Eの分布の推定に用いたパターンが1組(実施の形態1)の場合と、パターンが2組(実施の形態2)の場合と、パターンが2組(実施の形態2)の場合に更に弾性率Eの変化のスパース性を加えた場合(実施の形態3)の実験結果を示す図である。図14の縦軸は真の弾性率Eと推定された弾性率E’との平均二乗誤差(Root Mean Square Error: RMSE)であり、横軸は観察できる頂点数の割合(図9参照)である。

【0071】
図14からも分かるように、実施の形態2に対して弾性率Eの変化のスパース性を組み合わせた場合には、実施の形態1の場合(組数が「1」の場合)及び実施の形態2の場合(組数が「2」の場合)のいずれに対しても、推定精度を高められることが分かる。特に、実施の形態2に対して弾性率Eの変化のスパース性を組み合わせた場合には、観察できる頂点数の割合が低い場合でも高い推定精度を得ることができる。

【0072】
図15は、本実施の形態の場合における観測頂点数の割合と観測精度の関係を説明する図である。図15の縦軸は弾性率であり、横軸はブロック番号である。図15は、推定ブロック数Nが40の場合について表している。図では、真の弾性率Eの変化が矩形パルスで表され、推定された弾性率Eの変化が折れ線で表されている。
図15(A)に示すように観測頂点数が全体の49%の場合には2つの波形がおおよそ一致し、図15(B)に示すように観測頂点数が全体の66%の場合には2つの波形の重なり具合が一段と高くなっている。従って、図15が前提とする観測条件の場合には、観察対象である弾性体1の約50%程度を観察できればどの部位が硬いかをおおよそ推定でき(RMSE=0.091907MPa)、約60%以上を観察できればどの部位が硬いかをより高い精度(RMSE=0.024834MPa)で推定できることが分かる。

【0073】
<他の実施の形態>
前述の実施の形態においては、外力Fの大きさfと弾性率Eの分布がいずれも未知である場合について説明したが、外力Fの大きさfが既知である場合には、未知パラメータの数が減るため、弾性率Eの分布の推定精度を実施の形態1に比して高めることができる。

【0074】
前述した各実施の形態は、様々な用途において活用できる。例えば医師による診断や手術を支援する医療支援システムに情報処理システム100を活用できる他、医療分野以外の産業分野にも、観察や分析の用途にも、生産又は製造の用途にも活用できる。

【0075】
例えば情報処理システム100を医療支援システムとして活用する場合には、前述したように、視野が限られるために臓器の全体像を一度に観察できなくても、観察できる部位に生じた臓器の変化の情報(観察変位量u)を用いて、被検体である患者の臓器に固有の弾性率Eの分布を高精度に推定することができる。

【0076】
前述の実施の形態においては、撮像カメラ20を通じて静止画像を撮像することを前提に説明したが、撮像カメラ20は弾性率Eの推定対象である弾性体1を動画像として撮像できる撮像手段であってもよい。この場合には、撮像された動画像からキャプチャされた静止画像を用いて前述した観察変位量uを取得してもよい。

【0077】
また、前述の実施の形態においては、撮像カメラ20が1台である場合について説明したが、撮像カメラ20は複数台あってもよい。また、撮影方向も1つに限らず、複数あってもよい。外力Fによる変形を複数の方向から撮影した画像を用いることによっても観察変位量uの情報を増やすことができ、推定精度を向上させることができる。

【0078】
また、前述の実施の形態(例えば実施の形態1)では、推定ブロック数Nが固定値として与えられる場合を前提に説明したが、推定ブロック数Nを増やしながら推定結果の空間分解能を徐々に高める手法を採用してもよい。すなわち、弾性率Eの分布の推定処理を複数ステップに分割してもよい。
例えば初回の推定ステップでは、推定ブロック数Nが10個の場合について推定処理を実行して短時間のうちに推定結果を得て推定結果を表示装置40に表示し、次の推定ステップでは、推定ブロック数Nを予め定めた値(例えば20個)に増加して推定処理の実行と推定結果の表示の更新を繰り返してもよい。このような手法によれば、ユーザは、弾性率Eの分布のおおよそのイメージを早期に把握しつつ、時間の経過につれてより高い空間分解能で取得できる。

【0079】
以上、本発明の実施の形態について説明したが、本発明の技術的範囲は上記実施の形態に記載の範囲には限定されない。上記実施の形態に、種々の変更又は改良を加えたものも、本発明の技術的範囲に含まれることは、請求の範囲の記載から明らかである。
【符号の説明】
【0080】
1…弾性体、10…情報処理装置、100…情報処理システム、111…メッシュモデル生成部、112…観察変位量取得部、113…作用点受付部、114…変位量推定部、115…最適化部、116…表現部、200…メッシュモデル、E…弾性率、F…外力、f…外力の大きさ、u…観察変位量、u’…推定変位量、Pa…作用点、Pf…固定点
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8
【図10】
9
【図11】
10
【図12】
11
【図13】
12
【図14】
13
【図15】
14