WO2018123559A1

WO2018123559A1 - 情報処理装置、情報処理方法及びプログラム

Info

Publication number: WO2018123559A1
Application number: PCT/JP2017/044527
Authority: WO
Inventors: 恵中尾; 松田　哲也; 充樹森田
Original assignee: Kyoto University NUC
Current assignee: Kyoto University NUC
Priority date: 2016-12-28
Filing date: 2017-12-12
Publication date: 2018-07-05
Anticipated expiration: 2019-06-28
Also published as: JPWO2018123559A1; EP3564901A4; US20200121225A1; EP3564901A1

Abstract

情報処理装置は、弾性体の一部分を撮像した画像から取得した観察変位量と弾性体の三次元モデルについて算出される推定変位量のうち観察領域に対応する領域部分の推定変位量との差分を最小化する最適解として、既知の作用点に作用する外力の大きさと当該三次元モデルにおける弾性率の分布とを求める演算部を有する。

Description

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

　本発明は、情報処理装置、情報処理方法及びプログラムに関する。

　癌病変に対する触診にみられるように、生体組織の硬さに関する情報は、医師が診断を下す上で有用である。近年では、臓器などの生体組織の弾性情報を定量的かつ客観的に評価し、その分布を画像表示するエラストグラフィ技術が実用化されている。エラストグラフィ技術には、例えば核磁気共鳴弾性計測法（Magnetic Resonance Elastography: MRE）や超音波エラストグラフィがある。

　ＭＲＥは、磁気共鳴画像（Magnetic Resonanse Imaging: MRI）を用いる弾性率測定法であり、計測対象を外部から振動させた状態で撮影し、対象物体内を伝播する振動波を画像化する手法である。超音波エラストグラフィは、計測対象に外圧を加えたときに生じる組織のひずみ分布を利用する方法や対象組織内にせん断波を生じさせ、その伝播速度を利用する方法に代表される手法である。

O. Goksel, H. Eskandari, and S. E. Salcudean, "Mesh adaptation ｆor improving elasticity reconstruction using the FEM inverse problem,"IEEE Trans. Med. Imaging、 vol. 32、no. 2, pp. 408-418, 2013.p21-23

　ところで、エラストグラフィ技術により弾性情報を取得できる領域範囲は、波が伝播する局所領域（対象組織の表面部分）や圧力を加えた領域の近傍範囲に限られる。

　そこで、本発明は、局所的な変位情報から弾性体全体における弾性率の分布を推測する手法の提供を目的とする。

　請求項１に記載の発明は、弾性体の一部分を撮像した画像から取得した観察変位量と当該弾性体の三次元モデルについて算出される推定変位量のうち観察領域に対応する領域部分の推定変位量との差分を最小化する最適解として、既知の作用点に作用する外力の大きさと当該三次元モデルにおける弾性率の分布とを求める演算部を有する情報処理装置である。
　請求項２に記載の発明は、前記演算部は、次式に基づいて、前記外力の大きさと、前記弾性率の分布を求めることを特徴とする請求項１に記載の情報処理装置である。

　ここで、Ｅ^＊は、評価関数Ｊ（Ｅ，ｆ）を最小にするＥとｆの集合であり、Ｅは、弾性率であり、Ｅ_０は、弾性率の基準値であり、ｆは、外力の大きさであり、ｕ_０は、観察変位量であり、Ｌ_０（Ｅ）は、弾性率Ｅで表現される全体剛性マトリクスＫの逆行列成分のうち観察できる部分であり、Ｌ_０（Ｅ）ｆは、推定変位量であり、λは、係数であり、||　||_１は、１次元のノルムであり、||　||_２は、２次元のノルムである。
　請求項３に記載の発明は、弾性体の一部分を撮像した画像から観測変位量を取得する取得部と、前記弾性体に作用する外力の作用点を受け付ける受付部と、前記弾性体の三次元モデルにおける弾性率の分布と、前記外力の大きさとを変数として、当該三次元モデルの推定変位量を推定する推定部と、前記観測変位量と当該観測変位量に対応する前記推定変位量との差分を最小化する最適解として、前記外力の大きさと前記三次元モデルにおける前記弾性率の分布を求める最適化部とを有する情報処理装置である。
　請求項４に記載の発明は、前記最適化部は、次式に基づいて、前記外力の大きさと、前記弾性率の分布を求めることを特徴とする請求項３に記載の情報処理装置である。

　ここで、Ｅ^＊は、評価関数Ｊ（Ｅ，ｆ）を最小にするＥとｆの集合であり、Ｅは、弾性率であり、Ｅ_０は、弾性率の基準値であり、ｆは、外力の大きさであり、ｕ_０は、観察変位量であり、Ｌ_０（Ｅ）は、弾性率Ｅで表現される全体剛性マトリクスＫの逆行列成分のうち観察できる部分であり、Ｌ_０（Ｅ）ｆは、推定変位量であり、λは、係数であり、||　||_１は、１次元のノルムであり、||　||_２は、２次元のノルムである。
　請求項５に記載の発明は、最適解として求められた前記弾性率の分布に基づいて前記三次元モデルを表現する表現部を更に有することを特徴とする請求項１又は３に記載の情報処理装置である。
　請求項６に記載の発明は、前記画像は、単一の撮像カメラによって撮像された静止画像であることを特徴とする請求項１又は３に記載の情報処理装置である。
　請求項７に記載の発明は、前記最適解は、前記弾性率の分布の初期値に対する変化量の絶対値和と前記差分との総和を最小化する値として求められることを特徴とする請求項１又は３に記載の情報処理装置である。
　請求項８に記載の発明は、前記作用点の情報が複数与えられる場合、前記最適解は、個々の当該作用点について計算される前記差分の総和を最小化する値として求めることを特徴とする請求項１又は３に記載の情報処理装置である。
　請求項９に記載の発明は、前記作用点に対して複数の前記外力が与えられる場合、前記最適解は、個々の当該外力について計算される前記差分の総和を最小化する値として求めることを特徴とする請求項１又は３に記載の情報処理装置である。
　請求項１０に記載の発明は、前記作用点は複数であり、当該作用点にはそれぞれ異なる前記外力が作用されることを特徴とする請求項１又は３に記載の情報処理装置である。
　請求項１１に記載の発明は、前記最適解を求める空間分解能を低い状態から徐々に高めることを特徴とする請求項１又は３に記載の情報処理装置である。
　請求項１２に記載の発明は、前記作用点に作用する前記外力の大きさが与えられる場合、前記最適解は、前記差分を最小化する前記弾性率の分布として求められることを特徴とする請求項１又は３に記載の情報処理装置である。
　請求項１３に記載の発明は、弾性体の一部分を撮像した画像から観察変位量を取得する処理と、前記弾性体に作用する外力の作用点を受け付ける処理と、前記弾性体の三次元モデルにおける弾性率の分布と、前記外力の大きさとを変数として、当該三次元モデルの推定変位量を推定する処理と、前記観測変位量と当該観測変位量に対応する前記推定変位量との差分を最小化する最適解として、前記外力の大きさと前記三次元モデルにおける前記弾性率の分布を求める処理とを有する情報処理方法である。
　請求項１４に記載の発明は、コンピュータに、弾性体の一部分を撮像した画像から観察変位量を取得する処理と、前記弾性体に作用する外力の作用点を受け付ける処理と、前記弾性体の三次元モデルにおける弾性率の分布と、前記外力の大きさとを変数として、当該三次元モデルの推定変位量を推定する処理と、前記観測変位量と当該観測変位量に対応する前記推定変位量との差分を最小化する最適解として、前記外力の大きさと前記三次元モデルにおける前記弾性率の分布を求める処理とを実行させるプログラムである。

　請求項１記載の発明によれば、変位量を取得できる領域が弾性体の一部分に限られる場合でも弾性体全域における弾性率の分布を求めることができる。
　請求項２記載の発明によれば、変位量を取得できる領域が弾性体の一部分に限られる場合でも弾性体全域における弾性率の分布を求めることができる。
　請求項３記載の発明によれば、変位量を取得できる領域が弾性体の一部分に限られる場合でも弾性体全域における弾性率の分布を求めることができる。
　請求項４記載の発明によれば、変位量を取得できる領域が弾性体の一部分に限られる場合でも弾性体全域における弾性率の分布を求めることができる。
　請求項５記載の発明によれば、弾性率の分布を三次元モデル上で視覚的に確認できる。
　請求項６記載の発明によれば、弾性体の観察に使用する装置構成を簡素化できる。
　請求項７記載の発明によれば、局在性が高い弾性率の分布を求めることができる。
　請求項８記載の発明によれば、変位量の取得に用いた外力の作用点が１つだけの場合に比して推定精度を高めることができる。
　請求項９記載の発明によれば、変位量の取得に用いた外力が１つだけの場合に比して推定精度を高めることができる。
　請求項１０記載の発明によれば、変位量の取得に用いた作用点と外力がそれぞれ１つだけの場合に比して推定精度を高めることができる。
　請求項１１記載の発明によれば、時間の経過とともに精度の高い推定結果を得ることができる。
　請求項１２記載の発明によれば、外力の大きさが与えられない場合と比して推定精度を高めることができる。
　請求項１３記載の発明によれば、変位量を取得できる領域が弾性体の一部分に限られる場合でも弾性体全域における弾性率の分布を求めることができる。
　請求項１４記載の発明によれば、変位量を取得できる領域が弾性体の一部分に限られる場合でも弾性体全域における弾性率の分布を求めることができる。

実施の形態で説明する発明の概念イメージを説明する図である。情報処理システムのハードウェア構成例を示す図である。コンピュータプログラムの実行を通じて実現される機能構成例を説明する図である。弾性体のメッシュモデル例を示す図である。弾性体のメッシュモデルの変形例を説明する図である。二次元のメッシュモデルを示す図である。要素剛性マトリクスＫｅの重ね合わせ処理の概念を説明する図である。シミュレーション実験の真値を与える弾性率Ｅの分布例を示す図である。観察できる部分の割合を説明する図である。推定ブロック数と弾性率Ｅが推定される単位区画（ブロック）の大きさとの関係を説明する図である。推定ブロック数Ｎが１０個の場合、２０個の場合、４０個の場合の実験結果を示す図である。実施の形態２で使用する観察条件を説明する図である。弾性率Ｅの分布の推定に用いたパターンが１組（実施の形態１）の場合と、パターンが２組（実施の形態２）の場合の実験結果を示す図である。弾性率Ｅの分布の推定に用いたパターンが１組（実施の形態１）の場合と、パターンが２組（実施の形態２）の場合と、パターンが２組（実施の形態２）の場合に更に弾性率Ｅの変化のスパース性を加えた場合（実施の形態３）の実験結果を示す図である。実施の形態３の場合における観測頂点数の割合と観測精度の関係を説明する図である。

　以下、添付図面を参照して、本発明の実施の形態について詳細に説明する。

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

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

　一方で、作用点Ｐａに作用する外力Ｆの大きさｆを与える情報や弾性体１が周囲に固定されている部位（固定点）Ｐｆの位置情報は、既知であっても未知であっても良い。一般には、既知である情報が増えるほど、推定精度が高くなり、推定に要する計算時間も少なく済む。

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

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

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

＜実施の形態１＞
＜システム構成＞
　図２は、情報処理システム１００のハードウェア構成例を示す図である。
　情報処理システム１００は、弾性体１（図１参照）の弾性率の分布を推定する演算処理を実行する情報処理装置１０と、弾性体１の撮像に用いる少なくとも１つの撮像カメラ２０と、弾性体１のメッシュモデルの作成に使用する三次元画像データ（ボリュームデータ）を記憶する三次元画像データ記憶装置３０と、推定された弾性率の分布の表示等に用いる表示装置４０と、ユーザが操作するマウスやポインティングデバイス等で構成される入力装置５０とを有している。

　情報処理装置１０は、いわゆるコンピュータであり、プログラムを実行するＣＰＵ（Central Processing Unit）１１と、ＢＩＯＳ（Basic Input / Output System）やファームウェアなどのプログラムやデータを記憶するＲＯＭ（Read Only Memory）１２と、プログラムに作業エリアを与えるＲＡＭ（Random Access Memory）１３と、処理対象とする画像データを記憶する記憶部１４とにより構成される。記憶部１４は、ハードディスク装置や半導体メモリなどの記憶装置である。ＣＰＵ１１によるプログラムの実行を通じて実現される機能については後述する。

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

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

　表示装置４０は、例えば液晶ディスプレイで構成される。液晶ディスプレイは、液晶パネル、バックライトなどで構成される。なお、表示装置４０は、有機ＥＬ（ElectroLuminescence）ディスプレイ等でもよい

　図３は、ＣＰＵ１１によるコンピュータプログラムの実行を通じて実現される機能構成例を説明する図である。
　処理機能の観点から見た情報処理装置１０は、三次元のメッシュモデルを生成するメッシュモデル生成部１１１と、撮像カメラ２０から出力されるカメラ画像から観察部位についての変位量（観察変位量ｕ_ｏ）を取得する観察変位量取得部１１２と、作用点Ｐａの位置情報を受け付ける作用点受付部１１３と、作用点情報をメッシュモデルに適用してメッシュモデルの各部位における変位量（推定変位量ｕ_ｏ’）を計算する変位量推定部１１４と、推定変位量ｕ_ｏ’と観察変位量ｕ_ｏとの差分を最小化する弾性率Ｅ（結果として弾性率Ｅの関数である要素剛性マトリクスＫ）と外力Ｆの大きさｆを最適解として求める最適化部１１５と、最適解とする要素剛性マトリクスＫから特定される弾性率の分布をメッシュモデルに重ねて表現した画像を生成して表示装置４０に出力する表現部１１６とを有している。これらの機能部は、請求の範囲における「演算部」の一例である。

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

　図４は、弾性体１をメッシュ構造で表現したメッシュモデル２００の例を示す図である。図４は、外力Ｆ（図５参照）が作用点Ｐａ（図５参照）に加えられる前のメッシュモデル２００を表している。図中の丸印は最小要素のノードを示す。なお、図５は、メッシュモデル２００の右端下部を作用点Ｐａとして外力Ｆを加えて弾性変形させた後のメッシュモデル２００を表している。

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

　そこで、観察変位量取得部１１２には、三次元画像データ又はメッシュモデル２００（図４参照）とカメラ画像をマッチングし、観察変位量ｕ_ｏをメッシュモデル２００の個々の頂点に対応付ける処理機能が用意されている。ここでの対応付けの際には、三次元画像データ又はメッシュモデル２００の寸法とカメラ画像の寸法とを揃える処理も実行される。

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

　作用点受付部１１３は、作用点Ｐａの位置情報を、入力装置５０を通じて受け付ける機能部である。作用点受付部１１３は、例えばクリック操作が検出された際に、マウスカーソルが位置していたメッシュモデル２００上の部位を外力Ｆの作用点Ｐａとして受け付ける。
　作用点受付部１１３は、受け付けた作用点Ｐａの位置を、メッシュモデル２００上の頂点の識別番号で特定し、変位量推定部１１４に出力する。

　変位量推定部１１４は、線形有限要素法で使用する力の釣り合い式（すなわち、フックの法則を表現する式）を用いて、メッシュモデル２００を構成する全ての頂点について推定変位量ｕ_ｏ’を計算する機能部である。
　具体的には、変位量推定部１１４は、メッシュモデル２００の各頂点に加わる外力Ｆの大きさをｆ、各頂点における推定変位量をｕ_ｏ’、全体剛性マトリクスをＫとする力の釣り合い式（ｆ＝Ｋ・ｕ_ｏ’）を変形した式１に基づいて、推定変位量ｕ_ｏ’を計算する。

なお、添え字のｏは、観測できる部分に位置する頂点を意味し、Ｌは、全体剛性マトリクスＫの逆行列（Ｋ^－１）を意味する。

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

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

　図６に示すメッシュモデルの場合、三角形の要素Iは頂点番号１、２及び３で規定され、三角形の要素ＩＩは頂点番号２、３及び４で規定される。このとき、頂点番号ｉに作用する力ｆ_ｉ（ｉ＝１～４）と頂点番号ｉにおける変位量ｕ_ｉ（ｉ＝１～４）は、全体剛性マトリクスＫを用いた線形有限要素法での力の釣り合い式（式２）として与えられる。

　ここで、個々の頂点番号ｉ（ｉ＝１～４）に作用する外力Ｆの大きさｆ_ｉ（ｉ＝１～４）とその変位量ｕ_ｉ（ｉ＝１～４）との関係付けに用いる要素剛性マトリクスＫ_ｅは、以下の式３により導出される。

ここで、Ｖは三角形の要素の体積である。
　また、Ｂマトリクスは、要素の形状関数から求まるマトリクスであり、Ｄマトリクスは、弾性率Ｅとポアソン比νから求まるマトリクスである。

　Ｄマトリクスは、次式（式４）により表すことができる。

式４から分かるように、Ｄマトリクスは弾性率Ｅの関数である。従って、要素剛性マトリクスＫ_ｅも弾性率Ｅの関数として表現される。

　全体剛性マトリクスＫは、三角形の要素Ｉの要素剛性マトリクスＫ_ｅと、三角形の要素ＩＩの要素剛性マトリクスＫ_ｅとの和で与えられる。すなわち、全体剛性マトリクスＫは、頂点番号１、２、３で構成される三角形の要素Ｉの要素剛性マトリクスＫ_ｅと、頂点番号２、３、４で構成される三角形の要素ＩＩの要素剛性マトリクスＫ_ｅとの和で表される。

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

　ところで、本実施の形態のカメラ画像は、前述したように、作用点Ｐａは既知であるが大きさｆが未知の外力Ｆによって弾性変形された弾性体１の一部分を撮像した画像である。
　従って、式２のうち未知のパラメータは、弾性率Ｅの分布（全体剛性マトリクスＫ）と外力Ｆの大きさf（ベクトル）となる。

　そこで、変位量推定部１１４では、これら未知のパラメータにそれぞれ任意の値を設定して推定変位量ｕ_ｏ’を計算する。なお、推定変位量ｕ_ｏ’の初回推定時には、変位量推定部１１４は、未知パラメータである弾性率Ｅの分布（全体剛性マトリクスＫ）と外力Ｆの大きさfのそれぞれに初期値を与えて推定変位量ｕ_ｏ’を計算する。

　式２をマトリクスで表現すると、次式（式５）となる。

ここで、ｕ_ｏ’は観察できる部分の推定変位量であり、ｕ_ｉ’は観察できない部分の推定変位量である。また、Ｌ_ｏは弾性率Ｅで表現される全体剛性マトリクスＫの逆行列成分のうち観察できる部分であり、Ｌ_ｉは弾性率Ｅで表現される全体剛性マトリクスＫの逆行列成分のうち観察できない部分である。ｆ_ｏは観察できる部分の頂点に作用する外力Ｆの大きさ、ｆ_ｉは観察できない部分の頂点に作用する外力Ｆの大きさである。

　観察できる部分の推定変位量ｕ_ｏ’に着目すると、式５は、次式（式６）として表現できる。

ここで、作用点Ｐａが既知であるので、外力Ｆの大きさｆのうち作用点Ｐａ以外の頂点に対応する外力Ｆの大きさを０（ゼロ）で表すことができる。
　変位量推定部１１４は、この式６を用いて弾性率Ｅで表される成分Ｌ_ｏと外力Ｆの大きさｆにより推定変位量ｕ_ｏ’（初期値を含む）を計算する。

　最適化部１１５は、メッシュモデル２００の全ての頂点番号（観察できる部分だけでなく観察できない部分の頂点も含む）について推定された推定変位量ｕ_ｏ’と、撮像カメラ２０によって撮像された部分画像から特定された観察変位量ｕ_ｏとの差分を最小化する全体剛性マトリクスＫ（すなわち弾性率Ｅの分布）と外力ｆが得られるまで、未知パラメータの組み合わせを更新する。

　具体的には、最適化部１１５は、推定変位量ｕ_ｏ’と観察変位量ｕ_ｏとの差分を最小化する（すなわち、観察結果に最も近い変形を与える）未知パラメータを求める最適化問題を解く処理を実行する。当該処理を式で表すと次式（式７）となる。

ここでのＥ^＊は、評価関数Ｊ（Ｅ，ｆ）を最小にするＥとｆの集合を意味する。
　本実施の形態の場合、評価関数Ｊ（Ｅ，ｆ）は、推定変位量ｕ_ｏ’と観察変位量ｕ_ｏの差分のＬ^２ノルムで与えられる。

　最適化部１１５は、評価関数Ｊの最小値が求まるまで、未知パラメータＥとｆを更新し、更新された値を変位量推定部１１４にフィードバックする。すなわち、未知パラメータＥとｆの更新値を用いた変位量推定部１１４による推定変位量ｕ_ｏ’の再計算と、最適化部１１５による推定変位量ｕ_ｏ’の評価とが、評価関数Ｊの最小値が求まるまで繰り返し実行される。

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

　表現部１１６（図３参照）は、最適解として特定された弾性率Ｅの分布と外力Ｆの大きさのうち弾性率Ｅの分布を、メッシュモデル２００に重ねて表現した画像を表示装置４０（図２参照）に表示する。弾性率Ｅの分布がメッシュモデル２００と関連付けて表現されることにより、ユーザは、観察できない部分も含めた弾性体１の弾性率Ｅの分布を視覚的に理解できる。
　なお、表現部１１６は、弾性率Ｅの違いを色や輝度の違いとして表現する。また、表現部１１６は、メッシュモデル２００の各部に対応する弾性率Ｅの数値を表示してもよい。なお、表現部１１６は、必要に応じ、最適解として求まった外力Ｆの大きさｆを出力する。

＜シミュレーション実験による検証＞
　以下では、弾性体１について作成したメッシュモデル２００（すなわち、図４に示すメッシュモデル）と、弾性体１を局所的に観察して得られる観察変位量ｕ_ｏと、弾性体１に作用する外力Ｆの作用点Ｐａの情報とを用いて推定した弾性率Ｅの分布の精度を検証した結果について説明する。

　以下では、シミュレーション実験では、弾性体１に対応するメッシュモデル２００がＮ個の頂点を有するものとし、そのうちのＭ個（Ｎ＞Ｍ)の頂点の観察変位量ｕ_ｏを用いて弾性体１全体の弾性率Ｅの分布を推定する場合に、観測できる頂点の数が推定精度に与える影響について説明する。

　シミュレーション実験では、Visual C/C++を用いて、実施の形態に係る手法をプログラムとして記述し、作成されたプログラムを汎用計算機（ＯＳ：Windows 10 Pro、ＣＰＵ：IntelCorei7-6700K、メモリ：16GB）で実行した。
　なお、シミュレーションは、直方体形状のメッシュモデル２００（図５参照）について実行した。このメッシュモデル２００は、全体で９９個の頂点を有し、最小要素を四面体とする。

　また、図５に示すように、メッシュモデル２００の左端に位置する９個の頂点を固定端Ｐｆとし、それ以外の点を自由点とした。
　さらに、図５に示すように、外力Ｆの作用点Ｐａとして、メッシュモデル２００の右端下部に位置する３個の頂点を指定した。これら３個の頂点には、同じ大きさｆの外力Ｆがいずれも下向き（－ｚ軸方向）に作用するものとした。

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

　図８に示すメッシュモデル２００は、４０個の立方体ブロックのうち右端から左方向（－ｘ方向）に７番目と８番目の左列に位置する４個の立方体ブロックの弾性率Ｅが２．０［ＭＰａ］と硬く、他の部分の弾性率Ｅが１．０［ＭＰａ］である場合を示している。
　シミュレーション実験では、図８に示すメッシュモデル２００に外力Ｆを作用させた場合に観測されたカメラ画像と作用点Ｐａの位置に基づいて、メッシュモデル２００の全ての頂点の弾性率Ｅの分布を推定した。

　なお、図８に示すメッシュモデル２００は、図５の説明と同じく、弾性体１のメッシュモデル２００の右端下部にある３個の頂点に下方向（－ｚ軸方向）へ１０Ｎの外力Ｆを作用させて変形した状態を表している。勿論、情報処理装置１０にとっては、外力Ｆの大きさｆも未知パラメータである。

　前述したように、シミュレーション実験の目的は、観測できる頂点の数が弾性率Ｅの推定精度に与える影響を評価することである。
　そこで、シミュレーション実験では、観察できる頂点の数（観察変位量ｕ_０が与えられる頂点の数）が全体の１０％の場合、２０％の場合、…９０％の場合について評価した。図９は、メッシュモデル２００のうち観察できる頂点の数の割合を説明する図である。図９は、メッシュモデル２００を横方向（ｘ－ｚ平面）から観察した図であり、ｘ軸方向について１０個の小区分に分割されている。このうち１つの小区分を観察する場合が全体の１０％を観察する場合に相当し、２つの小区分を観察する場合が全体の２０％を観察する場合に相当する。以下同様に、９つの小区分を観察する場合が全体の９０％を観察する場合に相当する。

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

　以下のシミュレーション実験では弾性率Ｅの推定解像度を表す推定ブロック数Ｎを１０個とする場合、Ｎ＝２０とする場合、Ｎ＝４０とする場合について検討した。
　図１０は、推定ブロック数Ｎと弾性率Ｅが推定される単位区画（ブロック）の大きさとの関係を説明する図である。

　例えばＮ＝１０の場合は、メッシュモデル２００（図９参照）を１０個の単位区画（ブロック）で表現することを意味し、Ｎ＝２０の場合は、メッシュモデル２００を２０個の単位区画（ブロック）で表現することを意味し、Ｎ＝４０の場合は、メッシュモデル２００を４０個の単位区画（ブロック）で表現することを意味する。
　つまり、Ｎ＝４０の場合の単位区画（ブロック）は、図８において説明した立方体ブロックと同じ形状になる。当然ながら推定ブロック数Ｎの値が大きいほど（推定される単位区画の大きさが小さいほど）空間分解能が高くなり、弾性率Ｅの分布を高精細に推定できる。

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

　図１１によれば、メッシュモデル２００を１０個のブロックで表現する場合（低解像度で推定する場合）には、観測できる頂点の数がメッシュモデル２００全体の１５％程度でも高精度で弾性率Ｅを推定できることが分かる。また、この場合にはブロック数が少ないので、計算コストも低く済む。
　また、メッシュモデル２００を２０個のブロックで表現する場合（中解像度で推定する場合）には、観測できる頂点の数がメッシュモデル２００全体の５０％程度になると高精度で弾性率Ｅを推定できることが分かる。
　メッシュモデル２００を４０個のブロックで表現する場合（高解像度で推定する場合）には、観測できる頂点の数がメッシュモデル２００全体の７０％程度にならなければ最小二乗誤差が十分に小さくならないことが分かる。なお、ここでの低解像度、中解像度、高解像度は、単に３つの解像度間の相対的な関係を表したものにすぎない。

＜小括＞
　以上の通り、実施の形態１によれば、弾性体１の一部分（例えば全体の５０％）を撮像した画像から取得した観察変位量ｕ_ｏと弾性体１のメッシュモデル２００について算出される推定変位量ｕ_ｏ’のうち観察領域に対応する領域部分の推定変位量との差分（ｕ_ｏ－ｕ_ｏ’）を最小化する最適解として、既知の作用点Ｐａに作用する外力Ｆの大きさｆとメッシュモデル２００全域における弾性率Ｅの分布とを求めることにより、観察変位量ｕ_ｏを取得できる領域が弾性体１の一部分に限られる場合でも弾性体１全域における弾性率Ｅの分布を求めることができる。

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

＜実施の形態２＞
　本実施の形態では、前述した実施の形態１に比して、弾性率Ｅの分布の推定精度を高めるための手法について説明する。
　図１２は、実施の形態２で使用する観察条件を説明する図である。本実施の形態では、外力Ｆと作用点Ｐａの組を２組用意する。図１２（Ａ）は、外力Ｆ１と作用点Ｐａ１の第１の組（パターン１）を示し、図１２（Ｂ）は、外力Ｆ２と作用点Ｐａ２の第２の組（パターン２）を示す。
　なお、図１２に示すメッシュモデル２００は、図８の場合と同じく、周囲より硬い弾性率Ｅを有する立方体ブロックを網掛けで表現している。

　図１２（Ａ）に示す第１の組（パターン１）は、メッシュモデル２００の右端下部の３つの頂点を作用点Ｐａ１として、下向き（－ｚ軸方向）に第１の外力Ｆ１を作用させた場合の変形を表している。図１２（Ｂ）に示す第２の組（パターン２）は、メッシュモデル２００の右端のうち端面に向かって左辺に位置する３つの頂点を作用点Ｐａ２として、左向きに（－ｙ軸方向）に第２の外力Ｆ２を作用させた場合の変形を表している。

　なお、図１２においては、作用点Ｐａだけでなく外力Ｆが作用する方向も異なる例を示しているが、作用点Ｐａは１つで外力Ｆの作用方向が複数でもよいし、作用点Ｐａは２つで外力Ｆの作用方向は１つでもよい。
　この場合も、実施の形態１の場合と同様、情報処理装置１０は、各パターンについての観察変位量ｕ_ｏを求め、その後、それぞれについて観察変位量ｕ_ｏと推定変位量ｕ_ｏ’との差分（ｕ_ｏ－ｕ_ｏ’）を最小化する外力Ｆ１の大きさｆ１と外力Ｆ２の大きさｆ２と、メッシュモデル２００全域における弾性率Ｅの分布を最適解として求める。

　すなわち、最適化部１１５は、次式（式９）を満たす未知パラメータを最適解として求める。

ここで、Ｕ_０はｕ_1ｏとｕ_2ｏを列方向に並べた行列であり、Ｆはｆ₁とｆ₂を列方向に並べた行列である。また、３行目の式における２次元のノルムは行列のノルムであり、例えばフロベニアスノルムなどである。

　図１３は、弾性率Ｅの分布の推定に用いたパターンが１組（実施の形態１）の場合と、パターンが２組（実施の形態２）の場合の実験結果を示す図である。図１３の縦軸は真の弾性率Ｅと推定された弾性率Ｅ’との平均二乗誤差（Root Mean Square Error: RMSE）であり、横軸は観察できる頂点数の割合（図９参照）である。
　図１３からも分かるように、弾性率Ｅの分布の推定に使用する作用点Ｐａと外力Ｆの大きさｆとで規定される組数が「２」の場合には、実施の形態１の場合（組数が「１」の場合）に比して、観察できる部分の割合が同じであっても、推定精度を高めることができる。

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

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

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

　式１０の第２項がスパース項である。ここで、Ｅ_０は、ほとんどの頂点に適用される弾性率の基準値である。従って、スパース項は、各頂点について推定された弾性値Ｅと基準値との差分（すなわちΔＥ）について１次元のノルムＬ^１を求める式となる。
　本実施の形態では、係数λを０．１とした。この式の場合、弾性率Ｅが基準値と異なる一部の頂点についてのみ差分ΔＥ（＝Ｅ－Ｅ_０）が非０（ゼロ）となり、大部分の頂点においての差分ΔＥは０（ゼロ）となる。このスパース項は、推定された弾性率Ｅの分布の局在性が高い場合にノルムが小さくなる。従って、スパース項のノルムが大きい弾性率Ｅの分布は評価が低くなる。なお、係数λの値は一例である。

　図１４は、弾性率Ｅの分布の推定に用いたパターンが１組（実施の形態１）の場合と、パターンが２組（実施の形態２）の場合と、パターンが２組（実施の形態２）の場合に更に弾性率Ｅの変化のスパース性を加えた場合（実施の形態３）の実験結果を示す図である。図１４の縦軸は真の弾性率Ｅと推定された弾性率Ｅ’との平均二乗誤差（Root Mean Square Error: RMSE）であり、横軸は観察できる頂点数の割合（図９参照）である。

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

　図１５は、本実施の形態の場合における観測頂点数の割合と観測精度の関係を説明する図である。図１５の縦軸は弾性率であり、横軸はブロック番号である。図１５は、推定ブロック数Ｎが４０の場合について表している。図では、真の弾性率Ｅの変化が矩形パルスで表され、推定された弾性率Ｅの変化が折れ線で表されている。
　図１５（Ａ）に示すように観測頂点数が全体の４９％の場合には２つの波形がおおよそ一致し、図１５（Ｂ）に示すように観測頂点数が全体の６６％の場合には２つの波形の重なり具合が一段と高くなっている。従って、図１５が前提とする観測条件の場合には、観察対象である弾性体１の約５０％程度を観察できればどの部位が硬いかをおおよそ推定でき（ＲＭＳＥ＝０．０９１９０７ＭＰａ）、約６０％以上を観察できればどの部位が硬いかをより高い精度（ＲＭＳＥ＝０．０２４８３４ＭＰａ）で推定できることが分かる。

＜他の実施の形態＞
　前述の実施の形態においては、外力Ｆの大きさｆと弾性率Ｅの分布がいずれも未知である場合について説明したが、外力Ｆの大きさｆが既知である場合には、未知パラメータの数が減るため、弾性率Ｅの分布の推定精度を実施の形態１に比して高めることができる。

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

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

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

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

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

　以上、本発明の実施の形態について説明したが、本発明の技術的範囲は上記実施の形態に記載の範囲には限定されない。上記実施の形態に、種々の変更又は改良を加えたものも、本発明の技術的範囲に含まれることは、請求の範囲の記載から明らかである。

　１…弾性体、１０…情報処理装置、１００…情報処理システム、１１１…メッシュモデル生成部、１１２…観察変位量取得部、１１３…作用点受付部、１１４…変位量推定部、１１５…最適化部、１１６…表現部、２００…メッシュモデル、Ｅ…弾性率、Ｆ…外力、ｆ…外力の大きさ、ｕ_ｏ…観察変位量、ｕ_ｏ’…推定変位量、Ｐａ…作用点、Ｐｆ…固定点

Claims

　弾性体の一部分を撮像した画像から取得した観察変位量と当該弾性体の三次元モデルについて算出される推定変位量のうち観察領域に対応する領域部分の推定変位量との差分を最小化する最適解として、既知の作用点に作用する外力の大きさと当該三次元モデルにおける弾性率の分布とを求める演算部
　を有する情報処理装置。
　前記演算部は、次式に基づいて、前記外力の大きさと、前記弾性率の分布を求めることを特徴とする請求項１に記載の情報処理装置。

　ここで、
　Ｅ^＊は、評価関数Ｊ（Ｅ，ｆ）を最小にするＥとｆの集合であり、
　Ｅは、弾性率であり、
　Ｅ_０は、弾性率の基準値であり、
　ｆは、外力の大きさであり、
　ｕ_０は、観察変位量であり、
　Ｌ_０（Ｅ）は、弾性率Ｅで表現される全体剛性マトリクスＫの逆行列成分のうち観察できる部分であり、
　Ｌ_０（Ｅ）ｆは、推定変位量であり、
　λは、係数であり、
　||　||_１は、１次元のノルムであり、
　||　||_２は、２次元のノルムである。
　弾性体の一部分を撮像した画像から観測変位量を取得する取得部と、
　前記弾性体に作用する外力の作用点を受け付ける受付部と、
　前記弾性体の三次元モデルにおける弾性率の分布と、前記外力の大きさとを変数として、当該三次元モデルの推定変位量を推定する推定部と、
　前記観測変位量と当該観測変位量に対応する前記推定変位量との差分を最小化する最適解として、前記外力の大きさと前記三次元モデルにおける前記弾性率の分布を求める最適化部と
　を有する情報処理装置。
　前記最適化部は、次式に基づいて、前記外力の大きさと、前記弾性率の分布を求めることを特徴とする請求項３に記載の情報処理装置。

　ここで、
　Ｅ^＊は、評価関数Ｊ（Ｅ，ｆ）を最小にするＥとｆの集合であり、
　Ｅは、弾性率であり、
　Ｅ_０は、弾性率の基準値であり、
　ｆは、外力の大きさであり、
　ｕ_０は、観察変位量であり、
　Ｌ_０（Ｅ）は、弾性率Ｅで表現される全体剛性マトリクスＫの逆行列成分のうち観察できる部分であり、
　Ｌ_０（Ｅ）ｆは、推定変位量であり、
　λは、係数であり、
　||　||_１は、１次元のノルムであり、
　||　||_２は、２次元のノルムである。
　最適解として求められた前記弾性率の分布に基づいて前記三次元モデルを表現する表現部を更に有すること
　を特徴とする請求項１又は３に記載の情報処理装置。
　前記画像は、単一の撮像カメラによって撮像された静止画像であること
　を特徴とする請求項１又は３に記載の情報処理装置。
　前記最適解は、前記弾性率の分布の初期値に対する変化量の絶対値和と前記差分との総和を最小化する値として求められること
　を特徴とする請求項１又は３に記載の情報処理装置。
　前記作用点の情報が複数与えられる場合、前記最適解は、個々の当該作用点について計算される前記差分の総和を最小化する値として求めること
　を特徴とする請求項１又は３に記載の情報処理装置。
　前記作用点に対して複数の前記外力が与えられる場合、前記最適解は、個々の当該外力について計算される前記差分の総和を最小化する値として求めること
　を特徴とする請求項１又は３に記載の情報処理装置。
　前記作用点は複数であり、当該作用点にはそれぞれ異なる前記外力が作用されること
　を特徴とする請求項１又は３に記載の情報処理装置。
　前記最適解を求める空間分解能を低い状態から徐々に高めること
　を特徴とする請求項１又は３に記載の情報処理装置。
　前記作用点に作用する前記外力の大きさが与えられる場合、前記最適解は、前記差分を最小化する前記弾性率の分布として求められること
　を特徴とする請求項１又は３に記載の情報処理装置。
　弾性体の一部分を撮像した画像から観察変位量を取得する処理と、
　前記弾性体に作用する外力の作用点を受け付ける処理と、
　前記弾性体の三次元モデルにおける弾性率の分布と、前記外力の大きさとを変数として、当該三次元モデルの推定変位量を推定する処理と、
　前記観測変位量と当該観測変位量に対応する前記推定変位量との差分を最小化する最適解として、前記外力の大きさと前記三次元モデルにおける前記弾性率の分布を求める処理と
　を有する情報処理方法。
　コンピュータに、
　弾性体の一部分を撮像した画像から観察変位量を取得する処理と、
　前記弾性体に作用する外力の作用点を受け付ける処理と、
　前記弾性体の三次元モデルにおける弾性率の分布と、前記外力の大きさとを変数として、当該三次元モデルの推定変位量を推定する処理と、
　前記観測変位量と当該観測変位量に対応する前記推定変位量との差分を最小化する最適解として、前記外力の大きさと前記三次元モデルにおける前記弾性率の分布を求める処理と
　を実行させるプログラム。