(57)【特許請求の範囲】
【請求項１】
対象物体に関する複数の変形条件下における変形形状群を生成する変形形状群生成手段と、
前記変形形状群の夫々における特徴点の位置に基づいて前記変形形状群を位置合わせする位置合わせ手段と、
前記位置合わせされた変形形状群に基づいて、前記複数の変形条件下における前記対象物体の複数の変形形状を表す変形形状モデルを生成する変形形状モデル生成手段と、
前記変形形状モデルを用いて、所望の変形条件下における前記対象物体の変形を推定する変形推定手段と、を有することを特徴とする情報処理装置。
【請求項２】
前記位置合わせ手段は、前記特徴点が二つである場合には、二つの特徴点の位置を結ぶ直線の中点を一致させることを特徴とする請求項１に記載の情報処理装置。
【請求項３】
前記位置合わせ手段は、前記特徴点が三つである場合には、三つの特徴点の重心を一致させることを特徴とする請求項１に記載の情報処理装置。
【請求項４】
前記変形形状モデル生成手段は、前記特徴点が複数である場合には、特徴点の位置の変形を変形形状モデルとして生成することを特徴とする請求項１に記載の情報処理装置。
【請求項５】
前記対象物体は人体の乳房であり、前記対象物体の特徴点は少なくとも人体の乳頭を含むことを特徴とする請求項１乃至請求項４のいずれか一項に記載の情報処理装置。
【請求項６】
前記対象物体はＭＲＩで撮像されたものであることを特徴とする請求項１乃至請求項５のいずれか一項に記載の情報処理装置。
【請求項７】
対象物体に関する複数の変形条件下における変形形状群を生成する変形形状群生成手段と、
対象物体に関する複数の変形形状群を、該変形形状群間で共通する解剖学的特徴の位置に基づいて位置合わせする位置合わせ手段と、
前記位置合わせされた変形形状群に基づいて、前記複数の変形条件下における前記対象物体の複数の変形形状を表す変形形状モデルを生成する変形形状モデル生成手段と、
前記変形形状モデルを用いて、所望の変形条件下における前記対象物体の変形を推定する変形推定手段とを有することを特徴とする情報処理装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、形状変形を行う情報処理装置に関し、特に対象となる形状にあうように物体の画像を変形する情報処理装置に関する。
【背景技術】
【０００２】
医療の分野において、医師は、患者を撮影した医用画像をモニタに表示し、表示された医用画像を読影して、病変部の状態や経時変化を観察する。この種の医用画像を生成する装置としては、単純Ｘ線撮影装置、Ｘ線コンピュータ断層撮影装置（Ｘ線ＣＴ）、核磁気共鳴映像装置（ＭＲＩ）、核医学診断装置（ＳＰＥＣＴ、ＰＥＴなど）、超音波画像診断装置（ＵＳ）等が挙げられる。
【０００３】
例えば、乳腺科領域においては、ＭＲＩで撮影した画像上で乳房内における病変部の位置を同定したうえで、超音波画像診断装置によって当該部位の状態を観察するという手順で画像診断を行う場合がある。ここで、乳腺科における一般的な撮影プロトコルでは、伏臥位（うつ伏せの体位）でＭＲＩの撮影を行い、仰臥位（あお向けの体位）で超音波撮影を行うことが多い。このとき医師は、撮影体位の差異に起因する乳房の変形を考慮して、伏臥位のＭＲＩ画像から得た病変部位置から仰臥位における病変部位置を推定したうえで、当該病変部の超音波撮影を行っていた。
【０００４】
しかし、撮影体位の差異に起因する乳房の変形は非常に大きいために、医師が推定する病変部の位置が実際とずれてしまい、本来観察したい病変部の超音波画像を描出できない、あるいは、その探索に長い時間を要してしまうという課題があった。超音波の撮影体位と同じ仰臥位でＭＲＩを撮影すればこの課題に対処できるが、仰臥位の撮影は被検者の呼吸の影響を受けてしまうため、読影に必要な鮮明なＭＲＩ画像を得ることができないという新たな課題が生じる。
【０００５】
もし仮に、伏臥位で撮影したＭＲＩ画像に画像処理による変形を施して、仮想的に仰臥位で撮影したＭＲＩ画像を生成できれば、変形後のＭＲＩ画像から病変部位置を同定することで、撮影体位の差を気にすることなく病変部の超音波撮影を行うことができる。例えば、伏臥位で撮影したＭＲＩ画像を読影して当該画像上における病変部位置を得た後に、伏臥位から仰臥位への変形情報に基づいて、仮想的な仰臥位のＭＲＩ画像上における病変部位置を算出することができる。あるいは、生成した仮想的な仰臥位のＭＲＩ画像を読影することで、当該画像上における病変部位置を直接的に求めることもできる。
【０００６】
これを実現する方法として、例えば、非特許文献１に示す方法を用いれば、伏臥位のＭＲＩ画像を仰臥位のＭＲＩ画像と同一の形状に変形させることができる。この方法では、まず、物理シミュレーションに基づいて伏臥位のＭＲＩ画像から仮想的な仰臥位のＭＲＩ画像が生成される。次に、仮想的な仰臥位のＭＲＩ画像と、実際に仰臥位で撮影したＭＲＩ画像との変形位置合わせが、画素値の類似度に基づいて実行される。そして、以上の処理で得られた対応関係に基づいて、伏臥位のＭＲＩ画像を仰臥位のＭＲＩ画像と同一の形状に変形させる処理が実行される。
【０００７】
また、非特許文献２には、統計動きモデル（以下、ＳＭＭ（Ｓｔａｔｉｓｔｉｃａｌ  Ｍｏｔｉｏｎ  Ｍｏｄｅｌ））を用いて変形前後の形状間の対応付けを高速に行う技術が開示されている。この技術ではまず、変形前の形状データに対する物理シミュレーションによって、対象物体の変形に係るパラメータ（以下変形パラメータと呼ぶ）を様々に設定した場合の変形形状群を取得して、その結果を主成分分析することでＳＭＭを生成する。そして、別途取得した変形後の表面形状データと、ＳＭＭの表面部分の形状とを比較して変形の推定を行うことにより、変形前後の形状間の対応付けを行う。
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】Ｔ．Ｊ．Ｃａｒｔｅｒ，Ｃ．Ｔａｎｎｅｒ，Ｗ．Ｒ．Ｃｒｕｍ  ａｎｄ  Ｄ．Ｊ．Ｈａｗｋｅｓ，“Ｂｉｏｍｅｃｈａｎｉｃａｌ  ｍｏｄｅｌ  ｉｎｉｔｉａｌｉｚｅｄ  ｎｏｎ−ｒｉｇｉｄ  ｒｅｇｉｓｔｒａｔｉｏｎ  ｆｏｒ  ｉｍａｇｅ−ｇｕｉｄｅｄ  ｂｒｅａｓｔ  ｓｕｒｇｅｒｙ，”９ｔｈ  Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｂｉｏｍｅｃｈａｎｉｃｓ  ｆｏｒ  Ｍｅｄｉｃｉｎｅ，９ｔｈ  ＭＩＣＣＡＩ  Ｃｏｎｆｅｒｅｎｃｅ  Ｗｏｒｋｓｈｏｐ．
【非特許文献２】Ｙ．Ｈｕ，Ｄ．Ｍｏｒｇａｎ，Ｈ．Ｕ．Ａｈｍｅｄ，Ｄ．Ｐｅｎｄｓｅ，Ｍ．Ｓａｈｕ，Ｃ．Ａｌｌｅｎ，Ｍ．Ｅｍｂｅｒｔｏｎ  ａｎｄ  Ｄ．Ｈａｗｋｅｓ，“Ａ  ｓｔａｔｉｓｔｉｃａｌ  ｍｏｔｉｏｎ  ｍｏｄｅｌ  ｂａｓｅｄ  ｏｎ  ｂｉｏｍｅｃｈａｎｉｃａｌ  ｓｉｍｕｌａｔｉｏｎｓ，”Ｐｒｏｃ．ＭＩＣＣＡＩ  ２００８，Ｐａｒｔ  Ｉ，ＬＮＣＳ  ５２４１，ｐｐ．７３７−７４４，２００８．
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかし、非特許文献１に記載の方法による処理を適正に作用させるためには、対象物体の変形パラメータの正確な値を事前に取得する必要があった。すなわち、変形パラメータが不明な場合には、非特許文献１に記載の方法を適用できないという課題があった。なお、変形パラメータが未知である場合に、変形パラメータのあらゆるパターンに基づく変形を試行するというアプローチも考えられるが、多数の変形の試行には膨大な時間を要するという課題があった。
【００１０】
また、非特許文献２に記載の方法は、対象物体の輪郭形状だけを用いて変形の推定を行っているため、人体の乳房表面のように滑らかな曲面形状を対象とすると推定に曖昧さが残り、精度の高い変形推定が実現できないという課題があった。
【００１１】
本発明は、上記の課題に鑑みてなされたものであり、対象物体の変形パラメータが未知の場合に、変形条件の違いによって生じる変形を高精度かつ高速に推定する仕組みを提供することを目的とする。
【課題を解決するための手段】
【００１２】
上記の目的を達成するための、本発明の一態様によるによる情報処理装置は以下の構成を備える。すなわち、
第１の変形条件下の対象物体における第１の形状の情報と特徴領域の第１の位置情報とから、前記特徴領域の位置を基準とする前記第１の形状の変形をモデル化する変形形状モデル生成手段と、
第２の変形条件下の前記対象物体における第２の形状の情報と前記特徴領域に対応する第２の位置情報に基づき、前記第１の位置と前記第２の位置を合わせ前記第１の形状を前記第２の形状に前記モデルを用いて変形推定する変形推定手段と、を備える。
【発明の効果】
【００１３】
対象物体の変形パラメータが未知の場合であっても、変形を高精度かつ高速に推定する仕組みを提供できる。
【図面の簡単な説明】
【００１４】
【図１】第１実施形態に係る情報処理装置の機能構成を示す図である。
【図２】第１実施形態に係る情報処理装置の装置構成を示す図である。
【図３】第１実施形態に係る情報処理装置の処理手順を示すフローチャートである。
【図４】第１実施形態に係る対象物体および対象物体の形状取得を説明する図である。
【図５】第１実施形態に係る重力加速度を説明する図である。
【図６】第１実施形態に係るメッシュモデルを説明する図である。
【図７】第１実施形態に係るステップＳ２０９の処理手順を示すフローチャートである。
【発明を実施するための形態】
【００１５】
以下、添付図面に従って本発明に係る情報処理装置及び方法の好ましい実施形態について詳説する。ただし、発明の範囲は図示例に限定されるものではない。
【００１６】
［第１実施形態］
図１は、本実施形態に係る情報処理装置１の機能構成を示している。本実施形態に係る情報処理装置は、画像撮影装置２としてのＭＲＩ装置、および形状測定装置３としてのレーザレンジセンサと接続される。情報処理装置１は、第１の変形条件下における対象物体を画像撮影装置２で撮影することによって得た第１の３次元画像データを取得する。また、第２の変形条件下における対象物体を形状測定装置３で計測することによって得た対象物体の表面形状（以下、第２の表面形状）を取得する。そして、前記第１の３次元画像に写る対象物体の形状（以下、第１の形状）が、前記第２の変形条件下における対象物体の形状（以下、第２の形状）と略一致するように、第１の３次元画像を変形させた変形画像を生成して表示するものである。
【００１７】
情報処理装置１は、さらに以下に説明する機能により構成されている。
第１画像取得部１００は、画像撮影装置２が対象物体を第１の変形条件下で撮影した第１の３次元画像（原画像）を取得し、第１形状取得部１０１、第１特徴領域位置取得部１０２、および変形画像生成部１１０へ送信する。
【００１８】
第１形状取得部１０１は、第１画像取得部１００が取得した第１の３次元画像を処理して第１の変形条件下における対象物体の形状（第１の形状）に関する情報を抽出し、第１の形状を記述する形状モデル（以下、第１の形状モデル）を生成する。そして、生成した第１の形状モデルを、第１特徴領域位置取得部１０２へ送信する。
【００１９】
第１特徴領域位置取得部１０２は、第１画像取得部１００が取得した第１の３次元画像を処理して、第１の変形条件下における対象物体に係る所定の特徴領域を１個抽出する。そして、その位置（以下、第１の位置）に関する第１の位置情報を、第１の形状モデルに組みこむ処理を実行する。そして、特徴領域に関する情報を付加した第１の形状モデルを、変形形状群生成部１０４、変位ベクトル算出部１０９、および変形画像生成部１１０へ送信する。
【００２０】
仮想変形パラメータ取得部１０３は、対象物体の変形パラメータがとり得る値を仮想的に設定した仮想変形パラメータを後に説明する方法でｎｐ組取得して、変形形状群生成部１０４へ送信する。ここで、対象物体の変形パラメータとは、例えば、対象物体の変形に係る性質を定める材料力学的な物理量（例えば弾性係数など）や、第１及び第２の変形条件下において対象物体に作用する外力に関する情報などである。
【００２１】
変形形状群生成部１０４は、仮想変形パラメータ取得部１０３から受信する複数の仮想変形パラメータの夫々に基づいて、第１の形状モデルに物理シミュレーションを施す。そして、第１の形状を変形させた複数の変形形状（以下、変形形状群）を生成し、関連付けられた特徴領域位置の変位を算出する。さらに、複数の変形形状の夫々を、特徴領域位置を基準とした座標系に変換する。すなわち変形形状間で特徴領域位置が一致するように座標変換する。そして、これらの情報を変形形状モデル生成部１０５へ送信する。
【００２２】
変形形状モデル生成部１０５は、前記変形形状群に基づいて、第２の変形条件下において対象物体がとり得る様々な変形形状を近似表現可能な変形形状モデルを生成し、変形成分推定部１０８へ送信する。
【００２３】
第２形状取得部１０６は、形状測定装置３から、対象物体の第２の形状の一部に関する情報として、第２の変形条件下における対象物体の表面形状（第２の表面形状）に関する情報を取得する。本実施形態では、物体表面上に密に配置した点群の位置を表す３次元座標の集合からなるレンジデータが、対象物体の表面形状として形状測定装置３から供給されるものとする。
【００２４】
第２特徴領域位置取得部１０７は、第２形状取得部１０６が取得した第２の表面形状（レンジデータ）に基づいて、対象物体の特徴領域の位置（第２の位置）を抽出し、第２の位置情報を変形推定部１０９へ送信する。
【００２５】
変形成分推定部１０８は、第２形状取得部１０６が取得した第２の表面形状および、第２特徴領域位置取得部１０７が抽出した第２の位置に基づいて、前記変形形状モデルで第２の形状を記述するための変形成分推定値を算出する。そして、その推定値を、変位ベクトル算出部１０９へ送信する。
【００２６】
変位ベクトル算出部１０９は、前記変形成分推定値に基づいて、前記第１の形状モデルを前記第２の形状へと変形させる変位ベクトルを算出し、変形画像生成部１１０へ送信する。
【００２７】
変形画像生成部１１０は、前記第１の形状モデルと、前記変位ベクトルに基づいて、前記第１の３次元画像を第２の形状へと変形させた第２の３次元画像（変形画像）を生成し、画像表示部１１１へ送信する。
【００２８】
画像表示部１１１は、前記第２の３次元画像を表示する。
【００２９】
図２は、第１実施形態に係る情報処理装置及びそれと接続される装置の構成を示す図である。情報処理装置１は、例えばパーソナルコンピュータ（ＰＣ）などで実現することができ、中央演算処理装置（ＣＰＵ）１０、主メモリ１１、磁気ディスク１２、表示メモリ１３、モニタ１４、マウス１５、キーボード１６を有する。
【００３０】
ＣＰＵ１０は、主として情報処理装置１の各構成要素の動作を制御する。主メモリ１１は、ＣＰＵ１０が実行する制御プログラムを格納したり、ＣＰＵ１０によるプログラム実行時の作業領域を提供したりする。磁気ディスク１２は、オペレーティングシステム（ＯＳ）、周辺機器のデバイスドライバ、後述する変形推定処理等を行うためのプログラムを含む各種アプリケーションソフト等を格納する。表示メモリ１３は、モニタ１４のための表示用データを一時記憶する。モニタ１４は、例えばＣＲＴモニタや液晶モニタ等であり、表示メモリ１３からのデータに基づいて画像を表示する。マウス１５及びキーボード１６はユーザによるポインティング入力及び文字等の入力をそれぞれ行う。上記各構成要素は共通バス１７により互いに通信可能に接続されている。
【００３１】
情報処理装置１は、イーサネット（登録商標）等によるＬＡＮ（ローカルエリアネットワーク）を介して形状測定装置３と接続されており、形状測定装置３から対象物体の表面形状を取得できる。また、情報処理装置１は、イーサネット（登録商標）等によるＬＡＮ（ローカルエリアネットワーク）を介して画像撮影装置２と接続されており、画像撮影装置２から画像データを取得できる。なお、本発明の形態はこれに限定されず、これらの機器との接続は、例えば、ＵＳＢやＩＥＥＥ１３９４等の他のインターフェイスを介して行ってもよい。また、これらのデータを管理するデータサーバから、ＬＡＮ等を介して必要なデータを読み込む構成であってもよい。また、情報処理装置１に記憶装置、例えばＦＤＤ、ＣＤ−ＲＷドライブ、ＭＯドライブ、ＺＩＰドライブ等を接続し、それらのドライブから必要なデータを読み込むようにしてもよい。
【００３２】
本実施形態に係る情報処理装置１は、第１の変形条件下における対象物体の形状から、第２の変形条件下における対象物体の形状への変形を推定し、それに基づいて第１の３次元画像を変形して表示する。すなわち、第１の３次元画像に変形処理を施して、第２の変形条件下における３次元画像を仮想的に生成して表示することに相当する。
【００３３】
本実施形態では、人体の乳房を対象物体とする場合を例として説明する。本実施形態における第１の変形条件は、重力方向に対して乳房がうつ伏せ（伏臥位）の状態であるものとする。また第２の変形条件は、重力方向に対して乳房があお向け（仰臥位）であるものとする。すなわち、第１の変形条件と第２の変形条件では乳房に作用する重力の方向が異なることを意味している。この条件が互いに異なるため、第１の形状と、第２の形状との間には変形が生じる。本実施形態に係る情報処理装置１は、第１の変形条件と第２の変形条件の夫々における対象物体に作用する重力方向の差異と、対象物体の弾性係数（ヤング率、ポアソン比）とが未知の変形パラメータであるとして、変形位置合わせを実行する。
【００３４】
次に、本実施形態に係る情報処理装置１が実行する処理を、図３のフローチャートを用いて詳しく説明する。
【００３５】
以下では、図４の（ａ）に示すような、変形しない胸壁面４０１に接続されている乳房４００を撮影する場合を例として、夫々の処理の説明を行う。なお乳房４００上には乳頭４０２があるものとする。図４においては、図示の都合上、乳房４００、胸壁面４０１、乳頭４０２が２次元平面上の物体であるかのように示しているが、本実施形態でこれらは３次元的な形状を持つものとし、同図はその断面を表しているものとする。
【００３６】
本実施形態では、ヤング率、ポアソン比はそれぞれスカラー量とし、それぞれｐｙ、ｐｐと表記する。また、重力加速度の差異は乳房４００に働く重力加速度の３次元のベクトル量ｐｇであり、その各成分をｐｇｘ、ｐｇｙ、ｐｇｚとする。すなわち、本実施形態における乳房４００の変形パラメータｐは、数１に示す５次元ベクトルによって記述される。
【００３７】
【数１】
0005546230.tif 000002
【００３８】
ここで、以下の説明を円滑に行うために、重力加速度の差異ｐｇの意味について、図５を用いて詳しく説明する。図５は、重力加速度により対象物体の任意の局所領域が受ける力のベクトルを示すベクトル図である。ここでは図示の都合上、２次元のベクトル図を用いて説明するが、本実施形態におけるｐｇは３次元のベクトルであるため、この図による説明を３次元に拡張して適用するものとする。
【００３９】
図５において、第１の重力条件５００は、第１の変形条件で対象物体に働く重力加速度のベクトルをあらわしている。図５では第１の重力条件５００をベクトル図の原点におき、後に説明する無重力条件５０１および第２の重力条件５０２を、第１の重力条件５００との相対ベクトルとして表記する。
【００４０】
無重力条件５０１は、第１の重力条件５００を基準としたときの、無重力状態の重力加速度のベクトルのとり得る値を示したものである。ここで、無重力条件５０１の重力加速度ベクトルの絶対値は、地球上での重力加速度の大きさとして一意に決められる。ここでは、その絶対値をｇとする。ただし、そのベクトルの方向は定まらない。したがって、図５における無重力状態の力のベクトルのとり得る値は、第１の重力条件５００を中心とした半径ｇの円周上のいずれかの点である。
【００４１】
第２の重力条件５０２は、第２の変形条件で対象物体に加わる重力加速度ベクトルのとり得る値を示したものである。第２の変形条件で対象物体に加わる重力加速度ベクトルは、無重力状態を基準とすると、その絶対値は前記同様にｇとなる。ただし、そのベクトルの方向は定まらない。したがって、無重力状態を基準としたベクトル図で表現したとすると、無重力状態を原点とした半径ｇの円上のいずれかのベクトルをとり得る。一方、図５のように第１の重力条件を基準とした場合には、第２の変形条件において対象物体に加わる重力加速度ベクトルは、無重力条件５０１の円周上の任意の点を中心として半径ｇの円を描いたときの、その円周上のベクトルをとり得る。したがって、図５において第２の重力条件は、第１の重力条件を中心とした半径２ｇの円内のいずれかのベクトルをとり得る。したがって、本実施形態における外力に関する３次元ベクトル量ｐｇは、絶対値が２ｇ以下の任意の３次元ベクトル量をとり得ることになる。
【００４２】
（Ｓ２００）
ステップＳ２００において、第１画像取得部１００は、画像撮影装置２が乳房４００を第１の変形条件下で撮影することによって得られるＭＲＩ画像を、第１の３次元画像（原画像）として取得する。ここで、図４の（ｂ）は、第１画像取得部１００が取得する第１の３次元画像の一例を示している。第１の３次元画像４０３の情報は撮影範囲の３次元空間内で定義される撮影画像の輝度の関数として数２のように記述されるものとする。
【００４３】
【数２】
0005546230.tif 000003
【００４４】
ここで、ｘ、ｙ、ｚは撮影範囲の３次元空間内における位置座標を意味し、撮影装置または撮影画像を基準とした直交する座標系（以下、ＭＲＩ座標系）において、原点からｘ［ｍｍ］、ｙ［ｍｍ］、ｚ［ｍｍ］だけ並進した位置を表すものとする。
【００４５】
（Ｓ２０１）
ステップＳ２０１において、第１形状取得部１０１は、ステップＳ２００で取得した第１の画像４０３を処理することで、第１の変形条件下における乳房４００の形状を表す第１の形状モデルを生成する。この処理について図４を用いて詳しく説明する。
【００４６】
まず、第１形状取得部１０１は、第１の３次元画像４０３（図４の（ｂ））に対して輪郭抽出の処理を行うことで、図４の（ｃ）に示す輪郭形状４０４を得る。輪郭抽出の処理とは、第１の３次元画像４０３における乳房４００の内部と外部の境界である輪郭を求める処理である。この処理の具体的な方法としては、第１の３次元画像４０３の輝度値の空間勾配を算出して、それに閾値処理を施すなどの方法で輪郭形状４０４を求めることができる。これ以外にも、第１の３次元画像４０３における乳房４００の内部と外部の輝度値のヒストグラムや輝度値のパターン（テクスチャ）の違いに基づいて画像の領域分割を行い、その領域の境界を輪郭形状４０４として求めることもできる。輪郭形状４０４を求める処理は上記のいずれかの方法、またはこれらを組み合わせた方法、またはこれ以外のいかなる方法で実行してもよい。
【００４７】
次に、第１形状取得部１０１は、輪郭形状４０４を適当な間隔で区切り、図４の（ｄ）に示す輪郭ノード群４０５を配置する。輪郭ノード群４０５は、個々の輪郭ノードが夫々、３次元位置座標の情報を持つ。ここではｍ１’個の輪郭ノードにより輪郭ノード群４０５が構成されるものとし、その夫々の位置座標を、３次元位置座標ベクトルｓ１ｉ（１≦ｉ≦ｍ１’）と表記する。
【００４８】
次に、第１形状取得部１０１は、この輪郭ノード群４０５の中で、乳房４００が胸壁面４０１と接する位置に対応するノード群を固定ノード群４０６とする。例えば、領域分割処理によって胸壁面４０１の領域を検出・認識して、当該領域の近傍のノード（当該領域からの距離が予め定めた閾値以内のノード）のうちの連続するノード群を固定ノード群４０６と判別する。そして、固定ノード群４０６とそれ以外のノード群（以下、表面ノード群）との識別が可能なように、必要な情報を記録する。
【００４９】
次に、第１形状取得部１０１は、上記の処理で得た輪郭ノード群４０５に関する情報を用いて、第１の形状モデルを生成する。本実施形態では、以下に説明するメッシュモデルを用いて、第１の形状モデルを表現する。以下、第１形状取得部１０１が行うメッシュモデルの生成処理を、図６を用いて説明する。図６において、輪郭ノード群４０５と固定ノード群４０６は、図４で説明したものと同一である。まず、第１形状取得部１０１は、輪郭ノード群４０５に囲まれた領域内に内部ノード群７００を生成する。このとき内部ノード群７００は、輪郭ノード群４０５で囲まれた領域内を、例えば等間隔に分割した分割位置に配置すればよい。ただし本発明の実施において、内部ノード群７００はいかなる方法で配置してもよい。さらに、第１形状取得部１０１は、輪郭ノード群４０５および内部ノード群４０６の接続に関する情報を生成する。ノード群の接続に関する情報の生成は、前記ノード群に対してドロネー分割を適用するなどの方法で実現できる。以上の処理により生成された輪郭ノード群および内部ノード群の接続に関する情報を表すデータのことを、以下ではメッシュモデルと呼ぶ。また、メッシュモデルを構成するノードの夫々に、上記で得た夫々のノードの位置情報を付与したもの（すなわち、第１の変形条件下における乳房４００の形状を表すように、該メッシュモデルを形成したもの）を、第１の形状モデルと呼ぶ。以下の説明では、生成した内部ノード群７００の夫々の位置座標を、３次元位置座標ベクトルｓ１ｉ（ｍ１’＋１≦ｉ≦ｍ１）によって表記する（ｍ１はメッシュモデルのノード総数）。そして、数３に示すように、輪郭ノード群４０５と内部ノード群７００を合わせた全てのノードの位置座標を縦に並べた３×ｍ１次元ベクトルｓ１によって、第１の形状モデルを構成するノード群の位置情報を表現する。
【００５０】
【数３】
0005546230.tif 000004
【００５１】
なお、以上の処理によって生成される第１の形状モデルは各部に送信され、以降の処理で利用されるものとする。
【００５２】
（Ｓ２０２）
ステップＳ２０２において、第１特徴領域位置取得部１０２は、ステップＳ２００で取得した第１の画像４０３を処理し、乳房４００の表面上に存在する所定の特徴領域を抽出する。この特徴領域としては、例えば乳頭４０２を用いることが好適である。また、第１の画像４０３を撮像する際に、ＭＲＩで撮影可能な複数のマーカ（図４には不図示）を被検体の表面に貼り付けておいて、これを特徴領域として用いてもよい。以下の説明では、この特徴領域の位置を表す座標値をｖ１＝（ｘ１、ｙ１、ｚ１）と表記し、以後、これを第１の位置と呼ぶ。
【００５３】
さらに、ステップＳ２０２において、第１特徴領域位置取得部１０２は、上記の処理で取得した第１の位置に関する情報を、第１の形状モデルに組みこむ処理を実行する。具体的には、第１の特徴領域に関して、第１の形状モデルを構成する表面ノード群の中から、その位置が特徴領域位置ｖ１に最も近接するノードを探索し、当該特徴領域を表すノード（以下、特徴領域ノード）として設定する。すなわち、特徴領域を表す特徴領域ノードとして当該ノードのインデックスｎを記録するとともに、その位置ｓ１ｎをｖ１によって置換する。つまり、数４に示す関係となる。
【００５４】
【数４】
0005546230.tif 000005
【００５５】
例えば、図４の（ｂ）に示す乳頭４０２を抽出した場合に、乳頭４０２の位置に最も近接するノードが乳頭ノード４０７（図４の（ｄ））として設定される。そして、乳頭ノード４０７の位置情報が、本ステップで抽出した乳頭４０２の位置によって置き換えられる。
【００５６】
（Ｓ２０３）
ステップＳ２０３において、仮想変形パラメータ取得部１０３は、変形パラメータがとり得る値を仮想的に組み合わせた仮想変形パラメータを複数取得する。本実施形態では、ｎｐ個の仮想変形パラメータｐｋ（１≦ｋ≦ｎｐ）を取得する場合を例として説明する。
【００５７】
仮想変形パラメータｐｋは、変形パラメータの各成分値がとり得る範囲を適当な間隔で分割して、それらの全ての組み合わせを得ることによって生成する。このとき、各成分が対象物体の変形に与える影響の大きさなどに応じて、分割の細かさを変更する。例えば、仮想変形パラメータｐｋのうちｐｙ、ｐｐに関しては、とりうる範囲を１０００＜ｐｙ＜４０００  ［ｋＰａ］，０＜ｐｐ＜０．５とし、ｐｇｘ，ｐｇｙ，ｐｇｚに関してはｐｇｘ２＋ｐｇｙ２＋ｐｇｚ２≦（２ｇ）２を満たす範囲とする。そして、例えばｐｙは変形に与える影響が大きいものとして前記範囲を３０［ｋＰａ］づつ１０分割し、ｐｐは変形に与える影響が比較的に小さいものとして前記範囲を０．１づつ５分割する。またｐｇｘ，ｐｇｙ，ｐｇｚ関しては、夫々−２ｇから＋２ｇまでの範囲を５分割した組み合わせのうち前記範囲に関する条件を満たす組み合わせを設定する。
（Ｓ２０４）
ステップＳ２０４において、変形形状群生成部１０４は、ステップＳ２０３で取得した変形パラメータの複数の仮説（仮想変形パラメータ）の夫々に基づいて、第１の形状モデルに変形を施した変形形状群を生成する処理を実行する。変形形状算出部１０４が行う前記処理は、例えば有限要素法を用いた物理シミュレーションによって実現できる。
【００５８】
まず、仮想変形パラメータｐｋ（１≦ｋ≦ｎｐ）の夫々を仮定して、第１の形状モデルに対して有限要素法に基づく物理シミュレーションを施すことで、メッシュモデルを構成する夫々のノードの変位ベクトルｄｋｉ（１≦ｋ≦ｎｐ，  １≦ｉ≦ｍ１）を算出する。そして、数５に示す計算を実行することで、第１の形状モデルのノードの位置ｓ１ｉ（１≦ｉ≦ｍ１）に、変位ベクトルｄｋｉ（１≦ｋ≦ｎｐ，  １≦ｉ≦ｍ１）に基づく変位を施す。これにより、夫々のノードの変位後の位置ｓｄｋｉ（１≦ｋ≦ｎｐ，１≦ｉ≦ｍ１）を算出する。
【００５９】
【数５】
0005546230.tif 000006
【００６０】
なお、ステップＳ２０２の処理で述べたように、メッシュモデルを構成するノードの中には特徴領域を表す特徴領域ノードが含まれている。したがって、上記の処理を行うことで、仮想変形パラメータｐｋ（１≦ｋ≦ｎｐ）の夫々を仮定した場合の、特徴領域の変位後の位置ｖｄｋ（＝ｓｄｋｎ）も推定される。
【００６１】
さらに、特徴領域の変位後の座標値が３次元空間の所定の基準点に移動するように、数５で得たｓｄｋｉを並進させる処理を施し、全てのｓｄｋｉを更新する。本実施形態では、具体的な一例として、前記所定の基準点として原点を用いる場合について説明する。この場合、数５で得た各ノードの位置ｓｄｋｉから、特徴領域の位置ｓｄｋｎを減ずる処理を実行する。この減算処理により、全てのノードは−ｓｄｋｎだけ平行移動して、特徴領域ノードは全ての変形形状で原点に位置するようになる。すなわち、全ての変形形状が、特徴領域位置を基準とした座標系に変換される。なお、本発明の実施における前記所定の基準点は必ずしも原点を用いる必要は無く、任意の点を基準点として用いることができる。
【００６２】
最後に、夫々の仮想変形パラメータｐｋ（１≦ｋ≦ｎｐ）に関して、全ノードの位置座標ｓｄｋｉ（１≦ｉ≦ｍ１）を縦に並べた３×ｍ１次元ベクトルｓｄｋを生成する。そして、ベクトルｓｄｋによって、当該仮想変形パラメータｐｋを仮定した場合に第１の形状が変形するであろう形状（すなわち変形形状）を表現する。
【００６３】
以上で示したステップＳ２０４の処理を行うことで、変形形状群生成部１０４は、特徴領域の位置が原点で一致するような変形形状群ｓｄｋ（１≦ｋ≦ｎｐ）を生成する。なお、本実施形態では有限要素法に基づく物理シミュレーションを用いて変形形状に関する情報を生成する一実施形態について説明したが、本発明の実施はこれに限らない。例えば差分法や有限差分法などに基づく物理シミュレーションを用いて対象物体の変形形状を算出することもできる。またＭＰＳ法などのメッシュフリー手法を用いれば、メッシュモデルを用いずに変形形状を算出することもできる。ステップＳ２０４の処理は、仮想変形パラメータの夫々に基づく変形形状の算出ができる方法であれば、前記に示した方法以外のいかなる方法を用いて行ってもよい。
【００６４】
（Ｓ２０５）
ステップＳ２０５において、変形形状モデル生成部１０５は、ステップＳ２０４で求めた複数の変形形状に関する情報ｓｄｋ（１≦ｋ≦ｎｐ）に基づいて、対象物体の変形を近似表現する変形形状モデルを生成する。
【００６５】
変形形状モデルの生成方法には様々な方法を用いることができるが、例えば非特許文献２に開示されているＳＭＭ（Ｓｔａｔｉｓｔｉｃａｌ  Ｍｏｔｉｏｎ  Ｍｏｄｅｌ）を用いることができる。この方法によれば、変形形状群ｓｄｋ（１≦ｋ≦ｎｐ）に対して主成分分析を施すことで、複数の固有変形成分を抽出し、その固有変形成分の線形和によって対象物体の変形を近似表現することができる。この方法の具体的な処理について説明する。
【００６６】
まず、ステップＳ２０４で求めた複数の変形形状に関する情報ｓｄｋ（１≦ｋ≦ｎｐ）から、平均形状ｓｄ＿ａｖｅを数６により算出する。
【００６７】
【数６】
0005546230.tif 000007
【００６８】
次に、ｓｄｋ（１≦ｋ≦ｎｐ）から平均形状ｓｄ＿ａｖｅを引いた正規化変形形状群ｓｄｋ’（１≦ｋ≦ｎｐ）を算出する。そして、ｓｄｋ’（１≦ｋ≦ｎｐ）の分散共分散行列を求め、その行列の固有値分解を行い、固有値λｉ（１≦ｉ≦ｎｅ）および固有ベクトルｅｉ（１≦ｉ≦ｎｅ）を得る。ここで、ｎｅは算出する固有ベクトルの数であり、固有値の累積寄与率がある閾値を超えるように選択する。なお、以下では必要に応じて、固有ベクトルｅｉを固有変形成分と呼ぶ。
【００６９】
なお、数７に示すように、以上の処理により得たｓｄ＿ａｖｅおよびｅｉを線形結合することで、夫々の変形形状ｓｄｋ（１≦ｋ≦ｎｐ）を近似表現することができる。
【００７０】
【数７】
0005546230.tif 000008
【００７１】
ここで、ｃｋｉ（１≦ｉ≦ｎｅ）は、ｋ番目の変形形状ｓｄｋを表すための線形結合の係数である。
【００７２】
以上の処理により取得した平均形状ｓｄ＿ａｖｅおよび固有変形成分ｅｉ（１≦ｉ≦ｎｅ）を、乳房４００の変形形状モデルと呼ぶ。変形形状モデルは、第２の変形条件下における乳房４００の形状ｓ２をｓｄ＿ａｖｅとｅｉの線形結合によって表現するためのものである。次式に示す係数ｃｉ（１≦ｉ≦ｎｅ）の値を調整することで、第２の変形条件下において乳房４００がとり得るであろう任意の形状ｓｄを記述できる。
【００７３】
【数８】
0005546230.tif 000009
【００７４】
なお、ステップＳ２０２の処理で述べたように、メッシュモデルを構成するノードの中には特徴領域を表す特徴領域ノードが含まれている。したがって、上記の線形結合により、第２の変形条件下において乳房４００上の特徴領域がとり得るであろう位置ｖｄ（＝ｓｄｎ）も表現される。しかし、ステップＳ２０４の処理の内容から明らかなように、変形形状群における特徴領域の位置は全て原点と一致し、各変形形状間において変動しない。したがって、特徴領域と対応する平均形状ｓｄ＿ａｖｅおよび固有変形成分ｅｉにおける要素の値は、共にゼロとなり、上記線形結合の係数の大小に関わらず、常に原点を示すことになる。なお、以下では必要に応じて、形状ｓｄを構成する各ノードの位置座標をｓｄｉ（１≦ｉ≦ｍ１）、特徴領域ノードの位置座標をｖｄと表記する。
【００７５】
（Ｓ２０７）
ステップＳ２０７において、第２形状取得部１０６は、形状計測装置３から、乳房４００の第２の形状の一部に関する情報として、第２の変形条件下における乳房４００の表面形状（第２の表面形状）を表すレンジデータを取得する処理を実行する。ここでレンジデータは、物体表面上に密に配置したｍ２個の点群の、形状計測装置３によって定義される座標系（以下、レンジセンサ座標系）における位置を表す３次元座標の集合ｓ２ｊ（１≦ｊ≦ｍ２）によって構成される。
【００７６】
（Ｓ２０８）
ステップＳ２０８において、第２特徴領域位置取得部１０７は、ステップＳ２０２において第１の位置を取得した乳房４００に係る所定の特徴領域について、第２の変形条件下における特徴領域の位置（第２の位置）を取得する処理を実行する。この処理は、例えば、ステップＳ２０７で得た第２の表面形状から、突起部などの特徴的な形状の部位を抽出することで実行される。以下では、第２の位置を表す座標値をｖ２＝（ｘ２、ｙ２、ｚ２）と表記する。
【００７７】
（Ｓ２０９）
ステップＳ２０９において、変形成分推定部１０８は、変形形状モデルによる形状表現（すなわち数８のｓｄ）が、乳房４００の第２の形状ｓ２を最も適切に表すような、線形結合の係数の組ｃｉ（１≦ｉ≦ｎｅ）を推定する。すなわち、線形結合の係数の組を推定することで、乳房４００の第２の形状ｓ２を推定する。なお、本ステップで求める線形結合の係数の組をｎｅ次元ベクトルｃｅｓｔで表記し、以下ではこれを変形成分推定値と呼ぶ。
【００７８】
この処理は以下の手順で実行される。まずステップＳ２０８で取得した第２の位置ｖ２が原点に位置するように、第２の表面形状を表す座標群ｓ２ｊ（１≦ｊ≦ｍ２）を並進させる。これにより、ステップＳ２０５で生成した変形形状モデルによって記述する特徴領域ノードと第２の位置ｖ２を一致させる。そして、上記特徴領域を中心とした回転の推定値Ｒと上記変形成分推定値ｃｅｓｔの最適化を行う。ここで回転の推定値Ｒとは回転による座標変換を表す３×３行列を意味する。以下、変形成分推定部１０８の具体的な処理を、図７に示すフローチャートに沿って詳しく説明する。
【００７９】
（Ｓ１１００）＜第２の表面形状の座標変換＞
ステップＳ１１００において、変形成分推定部１０８は、ステップＳ２０７で得た第２の表面形状を表す座標群ｓ２ｊ（１≦ｊ≦ｍ２）を、特徴領域位置を基準とした座標系に変換する。すなわち、ステップＳ２０８で取得した第２の位置ｖ２に基づいて、第２の表面形状を表す座標群ｓ２ｊ（１≦ｊ≦ｍ２）を並進した座標群ｓ２ｊ’（１≦ｊ≦ｍ２）を、数９に示す計算により算出する。
【００８０】
【数９】
0005546230.tif 000010
【００８１】
（Ｓ１１０１）＜推定値の初期化＞
ステップＳ１１００において、変形成分推定部１０８は、変形成分推定値ｃｅｓｔおよび回転の推定値Ｒの初期化を行う。この初期化は、例えばｃｅｓｔをゼロベクトルとし、Ｒを単位行列とすることができる。
【００８２】
（Ｓ１１０２）＜変形形状の生成＞
ステップＳ１１０１において、変形成分推定部１０８は、現在の変形成分推定値ｃｅｓｔに基づいて、数１０の計算を行うことで、推定変形形状ｓｄ＿ｅｓｔを生成する。
【００８３】
【数１０】
0005546230.tif 000011
【００８４】
ただし、ｃｉ＿ｅｓｔはｃｅｓｔのｉ番目の成分を表す。
【００８５】
次に、変形成分推定部１０８は、ステップＳ１１００で得た座標群ｓ２ｊ’（１≦ｊ≦ｍ２）とステップＳ１１０２で得た推定変形形状ｓｄ＿ｅｓｔとを最も整合させる回転Ｒを推定する。例えば、以下に示す一連の処理（ステップＳ１１０３からステップＳ１１０６）によって、ＩＣＰ（Ｉｔｅｒａｔｉｖｅ  Ｃｌｏｓｅｓｔ  Ｐｏｉｎｔ）アルゴリズム等によるＲの更新を行う。
【００８６】
（Ｓ１１０３）＜表面形状の座標変換＞
ステップＳ１１０３において、変形成分推定部１０８は、現在の回転の推定値Ｒに基づいて、ステップＳ１１００で得た座標群ｓ２ｊ’（１≦ｊ≦ｍ２）に回転を施した座標群ｓ２ｊ＿ｒｏｔａｔｅ（１≦ｊ≦ｍ２）を、数１１に示す計算により算出する。
【００８７】
【数１１】
0005546230.tif 000012
【００８８】
（Ｓ１１０４）＜最近傍探索＞
ステップＳ１１０４において、変形成分推定部１０８は、ステップＳ１１０２で得た推定変形形状ｓｄ＿ｅｓｔを構成する表面ノードの中から、ステップＳ１１０３で得た座標群ｓ２ｊ＿ｒｏｔａｔｅ（１≦ｊ≦ｍ２）の夫々に最も距離が近接するノードを探索して選択する。そして、選択したノード（以下、ｓｄｊ＿ｅｓｔ’（１≦ｊ≦ｍ２）と表記する）を、ｓ２ｊ’（１≦ｊ≦ｍ２）の夫々の対応点として記憶する。
【００８９】
（Ｓ１１０５）＜座標変換（回転Ｒ）の更新＞
ステップＳ１１０５において、変形成分推定部１０８は、座標群ｓ２ｊ’（１≦ｊ≦ｍ２）と推定変形形状ｓｄ＿ｅｓｔとを最も整合させる回転Ｒを推定する。すなわち、ステップＳ１１０４で得たｓ２ｊ’（１≦ｊ≦ｍ２）とｓｄｊ＿ｅｓｔ’（１≦ｊ≦ｍ２）との対応点間の距離の平均値を求める数１２を評価関数として、評価値ｄを最小化するＲを算出する。
【００９０】
【数１２】
0005546230.tif 000013
【００９１】
なお、複数の対応点を用いて座標系変換を求める方法は周知であるので、詳細な説明は省略する。
【００９２】
（Ｓ１１０６）＜収束判定＞
ステップＳ１１０６において、変形成分推定部１０８は、回転Ｒを更新する処理の収束判定を行う。例えば、まず始めに、ステップＳ１１０５で算出した評価値ｄと、前の反復処理の際にステップＳ１１０５で算出した評価値ｄとの差Δｄを算出する。そして、差Δｄが予め設定した閾値よりも小さい場合には、ステップＳ１１０７に処理を進め、そうでない場合には、ステップＳ１１０３に処理を戻して回転Ｒの更新を反復する。すなわち、ステップＳ１１０３からステップＳ１１０５までの処理を、ステップＳ１１０６での収束判定が肯定されるまで反復実行する。
【００９３】
上記の方法により回転の推定値Ｒは更新され、以後の処理は更新後のＲに基づいて実行される。
【００９４】
（Ｓ１１０７）＜第２の表面形状の剛体変換＞
ステップＳ１１０７において、変形成分推定部１０８は、ステップＳ１１０３と同様な処理を実行する。すなわち、ステップＳ１１０６で得た回転の推定値Ｒに基づいて、座標群ｓ２ｊ’（１≦ｊ≦ｍ２）に回転を施した座標群ｓ２ｊ＿ｒｏｔａｔｅ（１≦ｊ≦ｍ２）を、数１１に示す計算により算出する。
【００９５】
（Ｓ１１０８）＜最近傍探索＞
ステップＳ１１０８において、変形成分推定部１０８は、ステップＳ１１０４と同様な処理を実行する。すなわち、推定変形形状ｓｄ＿ｅｓｔを構成する表面ノードの中から、ステップＳ１１０７で得た座標群ｓ２ｊ＿ｒｏｔａｔｅ（１≦ｊ≦ｍ２）の夫々に最も距離が近接するノードを探索することで、ｓ２ｊ’（１≦ｊ≦ｍ２）とｓｄｊ＿ｅｓｔ’（１≦ｊ≦ｍ２）とを対応づける処理を実行する。
【００９６】
（Ｓ１１０９）＜推定値の更新＞
ステップＳ１１０９において、変形成分推定部１０８は、ステップＳ１１０８で対応付けられた点間の誤差評価値ｄ’が小さくなるように、変形成分推定値ｃｅｓｔを更新する処理を実行する。このとき、誤差評価値ｄ’としては、例えば対応付けられた各点間のユークリッド距離の平均値を用いることができ、その場合は数１３に示す計算により得ることができる。
【００９７】
【数１３】
0005546230.tif 000014
【００９８】
これ以外にも、例えば各形状の輪郭面における法線方向も考慮して算出するようにしてもよい。このｃｅｓｔの更新処理は、一般的に知られた非線形な最適化問題を解くことで実行することができる。その具体的な方法として例えば貪欲法を用いることができる。この場合、現在のｃｅｓｔのある要素を微小に増加あるいは減少させた新たな係数を生成し、その係数によって（数１０を用いて）変形形状を生成する。そして、その変形形状に関して数１３の評価値ｄ’を算出する。その値が元のｃｅｓｔに基づく評価値ｄ’よりも小さくなる場合には、ｃｅｓｔの当該要素の値を更新する。この処理を、ｃｅｓｔの各要素に対して独立に実行することで、対応点間の距離が小さくなるようにｃｅｓｔを更新できる。さらに、上記の処理を繰り返して実行するようにして、より最適に近いｃｅｓｔを求めるようにしてもよい。最適化の方法はこれ以外にも、一般的に知られた非線形最適化手法の如何なるものを用いてもよい。例えば、最急降下法やニュートン法などを用いることができる。
【００９９】
上記の方法によりｃｅｓｔは更新され、以後の処理は更新後のｃｅｓｔに基づいて実行される。
【０１００】
（Ｓ１１１０）＜終了判定＞
ステップＳ１１１０において、変形成分推定部１０８は、ステップＳ２０９の処理の終了判定を行う。例えば、まず始めに、ステップＳ１１０９で更新された変形成分推定値ｃｅｓｔに基づいて、数１３の評価値ｄ’を算出する。そして、評価値ｄ’に基づいて、以後に行う処理の切り変えを行う。例えば、評価値ｄ’が予め設定した閾値よりも小さい場合には、本処理（すなわち、ステップＳ２０９の処理）を終了するようにし、そうでない場合には、ステップＳ１１０２に処理を戻して、変形成分推定値ｃｅｓｔの更新処理を継続する。すなわち、ステップＳ１１０２からステップＳ１１０５までの処理を、ステップＳ１１０５での終了判定が否定される限り繰り返す。
【０１０１】
以上に説明したステップＳ２０９の処理により、変形成分推定部１０８は、特徴領域を中心とした回転の推定値Ｒおよび変形成分推定値ｃｅｓｔを算出する処理を実行する。
【０１０２】
（Ｓ２１０）
ステップＳ２１０において、変位ベクトル算出部１０９は、ステップＳ２０９で算出したｃｅｓｔに基づいて数１０の計算を行い、乳房４００の第２の形状の推定値ｓ２＿ｅｓｔを得る。そして、第１の形状モデルの各ノードを第２の形状に変形させるための変位ベクトルｄｉ（１≦ｉ≦ｍ１）を、数１４を用いて算出する。
【０１０３】
【数１４】
0005546230.tif 000015
【０１０４】
ここで、ｓ１ｉは、第１の形状モデルにおけるｉ番目のノードの３次元位置座標ベクトルを示している。また、ｓ２ｉ＿ｅｓｔは、第２の形状の推定値ｓ２＿ｅｓｔが示すｉ番目のノードの３次元位置座標ベクトルであり、ｓ２＿ｅｓｔの３（ｉ−１）＋１、３（ｉ−１）＋２、３（ｉ−１）＋３番目の要素がこれに相当する。
【０１０５】
（Ｓ２１１）
ステップＳ２１１において、変形画像生成部１１０は、ステップＳ２００で取得した第１の３次元画像に変形を施し、変形後の乳房４００の形状が第２の形状と同様であるような第２の３次元画像（変形画像）を生成する。この変形は、ステップＳ２０１で生成した第１の形状モデルと、ステップＳ２１０で算出した各ノードの変位ベクトルｄｉ（１≦ｉ≦ｍ１）に基づいて、公知の画像変形手法によって実行される。
【０１０６】
（Ｓ２１２）
ステップＳ２１２において、画像表示部１１１は、ステップＳ２１１で生成した第２の３次元画像をモニタ１４などに表示する。
【０１０７】
以上に説明したように、本実施形態における情報処理装置１によれば、特徴領域位置を基準とした座標系に基づいて変形の推定を行うことで、特徴領域位置が一致することを保証しながら、簡便な計算で変形パラメータを推定できる。その結果、対象物体の変形パラメータが未知の場合に、第２の変形条件下における対象物体の形状に略一致するように変形させた変形画像を高速かつ高精度に生成して表示することができる。
【０１０８】
（変形例１−１）＜ＳＭＭの代わりに正規化変形形状を基底とするモデル＞
本実施形態では、ステップＳ２０５で変形形状モデル生成部１０５が行う処理として、ＳＭＭを用いる方法を例として説明したが、本発明の実施はこれに限らない。例えば、正規化変形形状群ｓｄｋ’をそのままｅｉとして用いて、以後の処理を実行するようにできる。この方法によれば、より簡易な処理で変形形状モデルを生成できる効果がある。またステップＳ２０５の処理としては、この方法に限らず、ステップＳ２０４で求めた複数の変形形状に関する情報ｓｄｋ（１≦ｋ≦ｎｐ）に基づいて、対象物体の変形を近似表現する変形形状モデルを生成する処理であれば、どのような処理であってもよい。
【０１０９】
（変形例１−２）＜両方の入力が３Ｄ画像の場合＞
本実施形態では形状測定装置３としてレンジセンサを用いていたが、形状測定装置３として画像撮影装置２を用いてもよい。
【０１１０】
この場合、第２の条件下における対象物体を画像撮影装置２によって撮影して３次元画像データ（第２の３次元画像データ）を得る。そして、第２形状取得部１０６は、ステップＳ２０７の処理において、画像撮影装置２から第２の３次元画像データを取得し、この画像から第２の変形条件下における対象物体の表面形状（すなわち、座標群ｓ２ｊ）を抽出する。また、第２特徴領域位置取得部１０７は、ステップＳ２０８の処理において、第１特徴領域位置取得部１０２と同様な処理により、第２の３次元画像データから第２の位置（すなわち、座標値ｖ２）を取得する。
【０１１１】
この場合、対象物体の特徴領域は対象物体の内部にあってもよい。例えば、物体内部において解剖学的に特徴を有する部位や腫瘤等を特徴領域としてもよい。このとき、第１特徴領域位置取得部１０２は、ステップＳ２０２の処理において、物体内部の特徴領域を抽出する処理を実行する。例えば、第１の３次元画像４０３の輝度値が周囲よりも高いような特徴を持つ領域を抽出するようにし、その位置を第１の位置とする。なお、物体内部の特徴領域に対応する特徴領域ノードは、内部ノード群４０６の中から選択する。ステップＳ２０８における第２特徴領域位置取得部１０７の処理も同様である。
【０１１２】
これらの方法によれば、第１の変形条件下で取得した３次元画像と第２の変形条件下で取得した３次元画像との位置合わせを行うことができる。なお、形状測定装置３として、画像撮影装置２とは異なるＭＲＩ装置や、Ｘ線ＣＴ装置等の他の３次元画像撮影装置を用いてもよいことはいうまでもない。
【０１１３】
（変形例１−３）＜ＭＲＩ／ＵＳ位置合わせ＞
また、第１の変形条件下で取得したＭＲＩ画像と第２の変形条件下で取得した超音波画像との位置合わせを行うことが目的の場合には、形状測定装置３として、位置センサを取り付けた超音波プローブを持つ超音波撮影装置を用いてもよい。
【０１１４】
特徴領域が物体表面上に存在する場合、第２形状取得部１０６は、ステップＳ２０７の処理において、被検体に接触させた超音波プローブを操作者が操作して得たプローブ先端位置の集合を、第２の変形条件下における物体の表面形状（座標群ｓ２ｊ）として取得する。また、操作者が超音波プローブを特徴領域に接触さた際のプローブ先端位置を、第２の位置（座標値ｖ２）として取得する。
【０１１５】
一方、特徴領域が物体内部に存在する場合には、第２形状取得部１０６は、ステップＳ２０７の処理において、被検体に接触させた超音波プローブを操作者が操作して撮影した超音波画像群と、夫々の超音波画像を撮影した際のプローブ位置を取得する。また、このときのプローブ先端位置の集合を、第２の変形条件下における物体の表面形状（座標群ｓ２ｊ）として取得する。そして、第２特徴領域位置取得部１０７は、ステップＳ２０８の処理において、ステップＳ２０７で得た超音波画像群の中から、物体内部の特徴領域を抽出する。そして、超音波画像上における該対応領域の座標と、該超音波画像を撮像した際のプローブ位置とを利用して、形状測定装置３が基準とする座標系における特徴領域の３次元位置を算出し、その位置を第２の位置（座標値ｖ２）とする。
【０１１６】
（変形例１−４）＜特徴領域の取得方法のバリエーション＞
本実施形態では、形状測定装置３で計測した被検体のレンジデータを用いて第２の位置を取得していたが、第２の位置を取得可能な方法であれば、本発明の実施はこれに限らない。例えば、磁気式の位置センサなどを備えることで先端部位置を計測可能なペン型の指示デバイス（スタイラス）を乳頭やマーカ等の特徴領域に接触させ、この接触点の位置を第２の位置として直接的に計測してもよい。
【０１１７】
また、本実施形態では、ステップＳ２０２で第１特徴領域位置取得部１０２が行う処理として、第１の３次元画像を処理することで第１の位置を求める場合を例として説明したが、本発明の実施はこれに限らない。例えば、マウス１５やキーボード１６などによってユーザの操作に基づいて、第１の位置を取得するようにしてもよい。これによれば、第１の３次元画像からの特徴領域を抽出する処理を省略できるため、より確実に（画像処理の安定性の制約を受けずに）本発明を実施することができる効果がある。この場合、第１特徴領域位置取得部１０２は、第１の３次元画像をモニタ１４などに表示してユーザに提示し、ユーザがその画像を視認しながら、その画像上に特徴領域に関する情報を設定できるようにすることが望ましい。
【０１１８】
また、第２特徴領域位置取得部１０７も同様に、マウス１５やキーボード１６などによってユーザの操作に基づいて、第２の位置を取得するようにしてもよく、前記同様の効果が期待できる。
【０１１９】
（変形例１−５）＜仮想変形パラメータの設定方法＞
ステップＳ２０３の処理における仮想変形パラメータに関する前記の範囲、分割数などは具体的な実施の一例に過ぎず、本発明の実施はこれに限らない。また、ステップＳ２０３の処理は、後段の処理に必要な仮想変形パラメータｐｋを取得できる方法であれば、いかなる方法で行ってもよい。例えば、マウス１５やキーボード１６等からなるユーザインタフェースを介してユーザがｐｋの値を入力し、仮想変形パラメータ取得部１０３がそれを取得する構成であってもよい。これ以外にも、各パラメータの範囲や分割の細かさなどをユーザが入力し、仮想変形パラメータ取得部１０３がその指示に基づいてｐｋの値を自動的に生成する構成であってもよい。また、磁気ディスク１２などにｐｋに関する情報をあらかじめ記録しておいて、仮想変形パラメータ取得部１０３がそれを取得できるようにしてもよい。
【０１２０】
（変形例１−６）＜ＭＲＩ以外＞
本実施形態では画像撮影装置２としてＭＲＩ装置を用いる場合を例として説明したが、本発明の実施はこれに限らない。例えば、Ｘ線ＣＴ装置、超音波画像診断装置、核医学装置などを用いることができる。
【０１２１】
［第２実施形態］
第１実施形態では、対象物体から１個の特徴領域を抽出して、それに基づいて処理を実行する場合について説明した。それに対し、第２の実施形態では複数の特徴領域を抽出することで、より高効率・高精度に処理を実行する場合について説明する。なお、本実施形態は、第１の実施形態の第１特徴領域位置取得部１０２，変形形状群生成部１０４，第２特徴領域位置取得部１０６および変形成分推定部１０８の処理を一部変更したものである。それ以外の機能については第１の実施形態と同様の内容であるため、説明を省略する。
【０１２２】
本実施形態に係る情報処理装置の機能構成は、第１の実施形態に係る情報処理装置１の機能構成を示す図１と同様である。
【０１２３】
本実施形態に係る情報処理装置の全体の処理フローは図３で説明した第１の実施形態の処理フローと同様である。ただし、本実施形態では、ステップＳ２０２，ステップＳ２０４，ステップＳ２０８およびステップＳ２０９で実行する処理の内容の一部が第１の実施形態とは異なるため、以下ではこれについて詳しく説明する。以下、具体的に２組の特徴領域を用いる場合、例えば、ＭＲＩで撮影された両乳房が対象物体であって、左右の乳頭を特徴領域として用いる場合を例として説明する。
【０１２４】
（Ｓ２０２）
ステップＳ２０２において、第１特徴領域位置取得部１０２は、対象物体の表面上に存在する複数の特徴領域（すなわち左右の乳頭）を、第１の実施形態と同様にして抽出し、その位置を取得する処理を実行する。以下では、取得した特徴領域の位置を表す座標値をｖ１ｈ（１≦ｈ≦２）と表記し、以後、これを第１の位置と呼ぶ。
【０１２５】
さらに、上記の処理で取得した第１の位置に関する情報を、第１の形状モデルに組みこむ処理を実行する。具体的には、第１の特徴領域の夫々に関して、第１の形状モデルを構成する表面ノード群の中から、その位置が特徴領域位置ｖ１ｈ（１≦ｈ≦２）に最も近接するノードを探索し、当該特徴領域を表すノード（以下、特徴領域ノード）として設定する。すなわち、特徴領域を表す特徴領域ノードとして当該ノードのインデックスｎｈ（１≦ｈ≦２）を記録するとともに、その位置ｓ１ｎｈをｖ１ｈによって置換する。
【０１２６】
（Ｓ２０４）
ステップＳ２０４において、変形形状群生成部１０４は、第１の実施形態と同様の物理シミュレーションによる変形を施した夫々の変形形状に対して、変位後の２個の特徴領域位置（ｓｄｋｎ１，ｓｄｋｎ２と表記する）を基準とした座標系へと座標変換を施す。ここで、変位後の２個の特徴領域位置を基準とした座標系とは、例えば、当該特徴領域の中点ｓ（すなわち、（ｓｄｋｎ１＋ｓｄｋｎ２）／２）を原点として、さらに、特徴領域間を結ぶ直線（ｓｄｋｎ１からｓｄｋｎ２に向かう直線ｌ）を所定の空間軸（例えばｘ軸）と定義するような座標系である。この座標系への座標変換は、例えば、前記中点ｓが原点となるように全てのノードを−ｓだけ平行移動した後に、直線ｌが所定の空間軸と一致するように、直線ｌと所定の空間軸とがなす平面と直交して原点を通る軸まわりの回転を施せばよい。
【０１２７】
（Ｓ２０８）
ステップＳ２０８において、第２特徴領域位置取得部１０７は、特徴領域の夫々に関して第２の位置を取得する処理を実行する。ここで、第２の位置を表す座標値をｖ２ｈ（１≦ｈ≦２）と表記する。
【０１２８】
（Ｓ２０９）
ステップＳ２０９における処理は、第１の実施形態と同様に図７を用いて説明できる。この図においてステップＳ１１００，Ｓ１１０５およびＳ１１０９の処理は第１の実施形態と一部処理の内容が異なるため、その差異を説明する。それ以外の処理の内容は第１の実施形態と同様の処理の内容であるため説明を省略する。
【０１２９】
（Ｓ１１００）
ステップＳ１１００において、変形成分推定部１０８は、ステップＳ２０７で得た第２の表面形状を表す座標群ｓ２ｊを、ステップＳ２０８で取得した２個の第２の位置ｖ２ｈを基準とした座標系に変換する。この座標変換は、上記特徴領域を結ぶ直線の中点が原点に移動する並進ｔを施した後に、特徴領域間を結ぶ直線が前記所定の空間軸と一致するような回転Ｒ’を施す剛体変換である。また、第２の位置ｖ２ｈの夫々にも同様の座標変換を施し、座標変換後の第２の位置ｖ２ｈ’を算出する。
【０１３０】
（Ｓ１１０５）
ステップＳ１１０３において、変形成分推定部１０８は、第１の実施形態と同様の処理を実行して、座標群ｓ２ｊ’と推定変形形状ｓｄ＿ｅｓｔとを最も整合させる回転Ｒを推定する。ただし、回転Ｒの推定は、前記所定の空間軸（すなわち、２個の特徴領域間を結ぶ直線）を回転軸とした回転成分に限定して推定を行う。
【０１３１】
（Ｓ１１０９）
ステップＳ１１０９において、変形成分推定部１０８は、第１の実施形態と同様の処理を実行する。なお、誤差評価値ｄ’の算出方法として、数１３に示した計算結果に、特徴領域の対応点間における距離（すなわち、Σ｜ｖｄｈ＿ｅｓｔ−Ｒ・ｖ２ｈ’｜）を所定の重みで誤差評価値ｄ’に加えてもよい。ここで、ｖｄｈ＿ｅｓｔは、ステップＳ１１０２で生成した推定変形形状ｓｄ＿ｅｓｔ中における、夫々の特徴領域ノードの位置座標を表している。
【０１３２】
以上に説明した第２の実施形態における情報処理装置１によれば、２組の特徴領域が略一致することを拘束条件とすることができるため、第１の実施形態に係る効果に加えて、さらに効率的に処理を実行することができる効果がある。具体的には第１の実施形態のステップＳ２０９における回転Ｒの最適化処理の探索空間が３自由度であるのに対して、本変形例によれば１自由度の探索空間に限定できる効果がある。
【０１３３】
（変形例２−１）＜特徴領域が３個＞
本実施形態では、２組の特徴領域を取得する場合を例として説明したが、本発明の実施はこれに限らない。例えば、３組の特徴領域を基準とした座標系で変形形状モデルを生成して、ステップＳ２０９ではそれらの特徴領域を基準とした座標系で変形位置合わせを実行するようにできる。
【０１３４】
この場合、ステップＳ２０４の処理において、変形形状群生成部１０４は、物理シミュレーションによる変形を施した夫々の変形形状に対して、変位後の３個の特徴領域位置を基準とした座標系へと座標変換を施す。この座標変換は、例えば当該特徴領域の重心位置が原点に、全特徴領域が所定の平面（例えばｘ−ｙ平面）上に、さらに前記特徴領域の中の所定の１つが所定の軸（例えばｘ軸）上に移動するような剛体変換である。
【０１３５】
また、ステップＳ１１００の処理において、変形成分推定部１０８は、ステップＳ２０７で得た第２の表面形状を表す座標群ｓ２ｊを、ステップＳ２０８で取得した３個の第２の位置を基準とした座標系に変換する。この座標変換は、上記３個の特徴領域の重心位置が原点に移動する並進ｔを施した後に、特徴領域を互いに結んで作られる平面が前記所定の空間平面と一致し、前記所定の特徴領域が前記所定の軸上に位置するような回転Ｒ’を施す剛体変換である。
【０１３６】
また、ステップＳ１１０３からステップＳ１１０６までの処理による回転Ｒの推定は、不要であるので省略できる。
【０１３７】
以上に説明した方法によれば、３組の特徴領域が略一致することを拘束条件とすることができるため、第２の実施形態に係る効果に加えて、さらに効率的に処理を実行することができる効果がある。具体的には第１の実施形態のステップＳ２０９における回転Ｒの最適化処理を省略できる効果がある。
【０１３８】
（変形例２−２）＜特徴領域が４個以上＞
本実施形態では、２組の特徴領域を用いる場合を例として説明し、またその変形例２−１では３組の特徴領域を用いる場合を例として説明したが、本発明の実施はこれに限らない。例えば、４組以上の特徴領域を基準とした座標系で変形形状モデルを生成して、ステップＳ２０９ではそれらの特徴領域を基準とした座標系で変形位置合わせを実行するようにできる。
【０１３９】
この場合、ステップＳ２０４の処理において、変形形状群生成部１０４は、物理シミュレーションによる変形を施した夫々の変形形状に対して、変位後の４個以上の特徴領域位置を基準とした座標系へと座標変換を施す。この座標変換は、例えば、第１の変形条件下における前記４個以上の特徴領域位置（ｖ１ｈ）を基準位置とし、その基準位置の夫々に対して前記変形を施した変形形状における対応する夫々の前記特徴領域の位置を略一致の状態にする剛体変換とすることができる。このとき、第１の変形条件下における特徴領域位置と変形を施した変形形状における対応する夫々の特徴領域との間の距離の総和が最小となることをもって前記略一致の状態とすることもできる。また、前記基準位置の重心と各特徴領域とを結ぶベクトルの角度の差の総和が最小となることをもって前記略一致の状態とすることもできる。
【０１４０】
また、ステップＳ１１００の処理において、変形成分推定部１０８は、ステップＳ２０７で得た第２の表面形状を表す座標群ｓ２ｊを、ステップＳ２０８で取得した４個以上の第２の特徴領域位置を基準とした座標系に変換する。この座標変換は例えば、第１の変形条件下における前記４個以上の特徴領域の位置を基準位置とし、その基準位置に対して前記第２の位置を略一致させる剛体変換とすることができる。
【０１４１】
また、ステップＳ１１０３の処理における回転Ｒの推定は、不要であるので省略できる。
【０１４２】
以上に説明した方法によれば、４組以上の特徴領域が略一致することを拘束条件とすることができるため、第２の実施形態に係る効果に加えて、さらに効率的に処理を実行することができる効果がある。具体的には第１の実施形態のステップＳ２０７における回転Ｒの最適化処理を省略することができる効果がある。また、上記方法は４組以上の特徴領域を取得した場合に限らず、例えば特徴領域が２組および３組の場合においても同様にして実行することができる。
【０１４３】
（変形例２−３）＜特徴領域が複数の場合には「第２の表面形状」を使わなくてもよい＞
本実施形態では特徴領域の組を取得して、こらの特徴領域位置を基準とした座標系を利用して、第２の表面形状と前記変形形状モデルとの位置合わせを行う場合を例として説明したが、本発明の実施はこれに限らない。例えば、第２の表面形状を用いずに、第２の位置だけを利用して前記変形形状モデルとの位置合わせを行うようにできる。
【０１４４】
この場合、ステップＳ１１０３からステップＳ１１０８の処理を省略することができる。またステップＳ１１０９では、第２の表面形状を用いずに、特徴領域の対応点間における距離の総和だけに基づいて、変形成分推定値を更新するようにできる。またステップＳ１１１０では、前記距離の総和だけに基づいて繰り返し処理の終了判定を行うようにできる。
【０１４５】
以上に説明した方法によれば、複数の特徴領域の組を取得する本発明の実施において、変形成分推定値の算出の処理を、より簡易に実行でき、処理を効率化・高速化できる効果がある。
【０１４６】
（変形例２−４）＜特徴領域が１〜４個以上＞
本実施形態およびその変形例２−１から変形例２−３および第１の実施形態では、予め定めされた特徴領域の組の数に応じた好適な処理方法の例について説明したが、本発明の実施はこれらの何れかの形態であることに限定されない。例えばステップＳ２０２およびステップＳ２０８で取得する特徴領域の組の数をユーザの指示によって切り替えられるようにし、その数に応じて前記の処理を切り替えて実行する仕組みを持つようにしても良い。これによれば、取得した特徴領域の組の数に基づいて好適な処理方法を選択して実行することが可能となり、ユーザの利便性を向上させる効果がある。
【０１４７】
［その他の実施形態］
また、本発明の目的は、前述した実施形態の機能を実現するソフトウェアのプログラムコードを記録した記録媒体（または記憶媒体）を、システムあるいは装置に供給し、そのシステムあるいは装置のコンピュータ（またはＣＰＵやＭＰＵ）が記録媒体に格納されたプログラムコードを読み出し実行することによっても、達成されることは言うまでもない。この場合、記録媒体から読み出されたプログラムコード自体が前述した実施形態の機能を実現することになり、そのプログラムコードを記録した記録媒体は本発明を構成することになる。
【０１４８】
また、コンピュータが読み出したプログラムコードを実行することにより、前述した実施形態の機能が実現されるだけでなく、そのプログラムコードの指示に基づき、コンピュータ上で稼働しているオペレーティングシステム（ＯＳ）などが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。
【０１４９】
さらに、記録媒体から読み出されたプログラムコードが、コンピュータに挿入された機能拡張カードやコンピュータに接続された機能拡張ユニットに備わるメモリに書き込まれた後、そのプログラムコードの指示に基づき、その機能拡張カードや機能拡張ユニットに備わるＣＰＵなどが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。
【０１５０】
本発明を上記記録媒体に適用する場合、その記録媒体には、先に説明したフローチャートに対応するプログラムコードが格納されることになる。
【０１５１】
なお、上述した本実施の形態における記述は、本発明に係る好適な情報処理装置の一例であり、本発明はこれに限定されるものではない。
【符号の説明】
【０１５２】
１  情報処理装置
１００  第１画像取得部
１０１  第１形状取得部
１０２  第１特徴領域位置取得部
１０３  仮想変形パラメータ取得部
１０４  変形形状群生成部
１０５  変形形状モデル生成部
１０６  第２形状取得部
１０７  第２特徴領域位置取得部
１０８  変形成分推定部
１０９  変位ベクトル算出部
１１０  変形画像生成部
１１１  画像表示部
２  画像撮影装置
３  形状測定装置
【図１】
0005546230.tif 000016
【図２】
0005546230.tif 000017
【図３】
0005546230.tif 000018
【図４】
0005546230.tif 000019
【図５】
0005546230.tif 000020
【図６】
0005546230.tif 000021
【図７】
0005546230.tif 000022