(57)【要約】      （修正有）【課題】異なる姿勢で撮影された三次元画像間の領域の対応を提示する。【解決手段】第一の変形条件における第一の形状データ群と、第一の形状データ群に対応する第二の変形条件における第二の形状データ群の間の形状と変形の統計的な関係を表現する統計モデルを取得する手段と、対象症例の第一または第二の少なくとも一方の変形条件における第三の形状を取得する手段と、第三の形状と略一致するように統計モデルのパラメータを調整する手段と、パラメータに基づいて、対象症例の変形を推定する手段とを備える。【選択図】図１
2014018376.tif 000002
【特許請求の範囲】
【請求項１】
複数症例の各々について、第1の変形条件における形状を示す第1の形状データと第2の変形条件における形状を示す第2の形状データとを取得する取得手段と、
前記第1の形状データと前記第2の形状データとから、前記第1の変形条件と前記第2の変形条件に関する、症例の形状と変形の統計的な関係を示す統計モデルを生成する生成手段と、を備えることを特徴とする画像処理装置。
【請求項２】
前記第1の形状データおよび前記第2の形状データは、一の症例についてそれぞれの形状を表すノードの位置座標値を並べた第1のベクトルおよび第2のベクトルであり、前記第1のベクトルと前記第2のベクトルの対応する要素は対応するノードの位置座標値を表すことを特徴とする請求項1に記載の画像処理装置。
【請求項３】
前記生成手段は、同一の症例に関する前記第1のベクトルと前記第2のベクトルを並べたベクトル組について主成分分析を行ない、得られた主成分と前記ベクトル組の平均ベクトルにより、前記統計モデルを生成することを特徴とする請求項２に記載の画像処理装置。
【請求項４】
前記生成手段は、前記複数の症例を属性に従って複数のグループに分類し、該グループ毎に前記統計モデルを生成することを特徴とする請求項1乃至3のいずれか1項に記載の画像処理装置。
【請求項５】
前記生成手段は、異なる主成分数で主成分分析を行ない、主成分数毎に統計モデルを生成することを特徴とする請求項3に記載の画像処理装置。
【請求項６】
前記取得手段は、前記第1の変形条件で得られた画像と前記第2の変形条件で得られた画像から対応する対応点を抽出し、前記対応点の位置に基づいて、前記第1の形状データと前記第2の形状データを取得することを特徴とする請求項1乃至5のいずれか1項に記載の画像処理装置。
【請求項７】
前記取得手段は、前記第1の形状データと前記第2の形状データに対称性がある場合、一方の形状データを対称変換することにより他方の形状データを取得することを特徴とする請求項1乃至6のいずれか1項に記載の画像処理装置。
【請求項８】
前記取得手段は、前記第1の変形条件で得られた画像と前記第2の変形条件で得られた画像から、ユーザによる入力操作により対応する対応点を抽出し、前記対応点の位置に基づいて、前記第1の形状データと前記第2の形状データを取得することを特徴とする請求項1乃至5のいずれか1項に記載の画像処理装置。
【請求項９】
前記生成手段は、前記統計モデルとして、
前記複数症例から得られる複数の前記第1のベクトルを主成分分析し、得られた主成分と前記第1のベクトルの平均ベクトルにより第1の統計モデルを生成し、
前記複数の症例から得られる複数の前記第2のベクトルを主成分分析し、得られた主成分と前記第2のベクトルの平均ベクトルにより第2の統計モデルを生成し、
前記第1および第2の統計モデルの結合パラメータにより得られるベクトルを主成分分析して第3の統計モデルを生成する、ことを特徴とする請求項2に記載の画像処理装置。
【請求項１０】
前記生成手段は、
前記複数の症例の前記第1の変形条件で得られた複数の画像から得られる輝度値の分布を示す第1の輝度ベクトルを主成分分析し、得られた主成分と前記第1の輝度ベクトルの平均ベクトルにより第4の統計モデルを生成し、
更に、前記複数の症例の前記第2の変形条件で得られた複数の画像から得られる輝度値の分布を示す第2の輝度ベクトルを主成分分析し、得られた主成分と前記第2の輝度ベクトルの平均ベクトルにより第5の統計モデルを生成し、
前記第3の統計モデルを、前記第1、第2、第4および第5の統計モデルの結合パラメータにより得られるベクトルを主成分分析することにより取得することを特徴とする請求項９に記載の画像処理装置。
【請求項１１】
前記生成手段は、
前記複数の症例の臨床情報を数値化したベクトルを主成分分析し、得られた主成分と該ベクトルの平均ベクトルにより第4の統計モデルを生成し、
前記第3の統計モデルを、前記第1、第2および第4の統計モデルの結合パラメータにより得られるベクトルを主成分分析することにより取得することを特徴とする請求項９に記載の画像処理装置。
【請求項１２】
請求項1乃至11のいずれか1項に記載された画像処理装置により生成された統計モデルを取得する取得手段と、
対象について前記第1の変形条件における形状を示す第1の対象形状データと前記統計モデルが表す形状が略一致するように前記統計モデルのパラメータを算出する算出手段と、
前記算出手段により算出されたパラメータに基づいて、前記第1の対象形状データから前記第2の変形条件における形状を示す第2の対象形状データへの変形を推定する推定手段と、を備えることを特徴とする画像処理装置。
【請求項１３】
前記対象について前記第1の変形条件における形状を示す第1の対象形状データの座標系と、前記第1の対象形状データとは異なるモダリティを用いて得た前記第2の変形条件における形状を示す第2の対象形状データの座標系を、人体座標系に変換する変換手段を更に備え、
前記算出手段と前記推定手段は、前記変換手段により座標系が変換された画像を用いて行うことを特徴とする請求項12に記載の画像処理装置。
【請求項１４】
前記第2の対象形状データは超音波画像であり、
前記第2の対象形状データにおける前記対象の形状は、前記第2の対象形状データを取得する際の超音波プローブの位置を検出することにより取得されることを特徴とする請求項13に記載の画像処理装置。
【請求項１５】
複数症例の各々について、第1の変形条件における形状を示す第1の形状データと第2の変形条件における形状を示す第2の形状データとを取得する工程と、
前記第1の形状データと前記第2の形状データとから、前記第1の変形条件と前記第2の変形条件に関する、形状と変形の統計的な関係を示す統計モデルを生成する生成工程と、を備えることを特徴とする画像処理方法。 
【請求項１６】
請求項15に記載された画像処理方法により生成された統計モデルを取得する工程と、
対象について前記第1の変形条件における形状を示す第1の対象形状データと前記統計モデルが表す形状が略一致するように前記統計モデルのパラメータを算出する工程と、
前記算出する工程により算出されたパラメータに基づいて、前記第1の対象形状データから前記第2の変形条件における形状を示す第2の対象形状データへの変形を推定する工程と、を備えることを特徴とする画像処理方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、核磁気共鳴映像装置（MRI）、Ｘ線コンピュータ断層撮影装置（X線CT）、超音波画像診断装置（US）など、種々の医用画像収集装置（モダリティ）で撮影した医用画像の好適な画像処理装置及び方法に関する。
【背景技術】
【０００２】
医療の分野において、あるモダリティの画像上に注目部位があった場合に、その注目部位に対応する部位（対応部位）を他のモダリティの画像上で同定し、その二つの部位の対比によって診断を行う場合がある。モダリティ間で撮影体位が同じ場合には容易に同定・対比を行うことができる。しかしながら、モダリティ間で撮影体位が異なる場合には、撮影時の被検体の形状が異なるため、その同定・対比が難しくなる。そこで、双方の被検体の変形を推定すること（すなわち、変形を伴う画像間の位置合わせを行うこと）が試みられている。これにより、注目部位の位置情報に基づいて対応部位の位置を推定することや、一方の画像に変形を施して形状を他方と同一にした画像を生成することが可能となる。
【０００３】
例えば非特許文献1には、伏臥位で撮影されたMRI画像と、仰臥位で撮影される超音波画像とを位置を合わせて表示する技術が開示されている。具体的には、伏臥位で撮影されたMRI画像に基づいて有限要素モデル(FEM)を生成し、そのモデルを用いて伏臥位から仰臥位への変形シミュレーションを行う。そして、このシミュレーション結果に基づいて伏臥位で撮影されたMRI画像やMRI画像中に描出されている病変領域などを、仰臥位の超音波画像上に重畳して表示する。この表示を用いて、伏臥位と仰臥位との間の変形による位置の違いを補正できる。
【０００４】
また非特許文献2には、有限要素モデルによる変形シミュレーションを実行することにより、対象物の材料特性や境界条件が明らかでない場合に対処できる技術が開示されている。具体的には、予め様々な材料特性や境界条件を想定して変形シミュレーションを実行し、その結果得られる群から、対象物の変形を比較的に少ない数のパラメータで表現するモデルを生成する。そして、そのモデルを用いて対象物の変形を推定する。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】T. Carter, C. Tanner, N. Beechey-Newman, D. Barratt and D. Hawkes, “MR navigated breast surgery: Method and initial clinical experience,” MICCAI2008.
【非特許文献２】Y. Hu, D. Morgan, H. Ahmed, D. Pendse, M. Sahu, C. Allen, M. Emberton and D. Hawkes, "A Statistical Motion Model Based on Biomechanical Simulations," MICCAI2008.
【非特許文献３】Warren Cheung, Ghassan Hamarneh, “N-SIFT: N-Dimensional scale invariant feature transform for matching medical images,” 4th IEEE International Symposium on Biomedical Imaging: From Nano to Macro, 2007.（実施形態のみで参照される）
【発明の概要】
【発明が解決しようとする課題】
【０００６】
非特許文献1に記載の技術を用いた場合、伏臥位と仰臥位との間の高い精度の変形推定が期待できる。しかしながら、変形シミュレーションを繰り返して実行する必要があるため、変形推定の処理に多大な時間がかかるという課題があった。また、非特許文献2に記載の技術を用いた場合、変形推定の処理の時間が軽減されるという効果がある。しかしながら、対象症例に関して有限要素モデルを生成する必要があり、その作業が煩雑であるという課題があった。
【０００７】
本発明は、上記の課題に鑑みてなされたものであり、変形状態の異なる複数の医用画像間の位置合わせを、高速かつ簡便に行える仕組みを提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記の目的を達成するための、本発明の一実施形態による画像処理装置は以下の構成を備える。すなわち、
複数症例の各々について、第1の変形条件における形状を示す第1の形状データと第2の変形条件における形状を示す第2の形状データとを取得する取得手段と、
前記第1の形状データと前記第2の形状データとから、前記第1の変形条件と前記第2の変形条件に関する、症例の形状と変形の統計的な関係を示す統計モデルを生成する生成手段と、を有する。
【発明の効果】
【０００９】
本発明によれば、変形状態の異なる複数の医用画像間の位置合わせを、高速かつ簡便に行える仕組みを提供できる。
【図面の簡単な説明】
【００１０】
【図１】第1実施形態による処理システムの機能構成図。
【図２】第1実施形態による処理システムの装置構成図。
【図３】第1実施形態による学習フェーズの処理手順を示すフローチャート。
【図４】第1実施形態による変形推定フェーズの処理手順を示すフローチャート。
【図５】第1実施形態によるS310の処理手順を示すフローチャート。
【図６】第1実施形態によるS310の処理の説明図。
【図７】第1実施形態による処理装置のS320の処理手順を示すフローチャート。
【図８】第1実施形態による処理装置のS3200およびS3210の処理の説明図。
【図９】第1実施形態による処理装置のS3230の処理の説明図。
【図１０】第1実施形態による処理装置のS3240の処理の説明図。
【図１１】第1実施形態による処理装置のS420の処理手順を示すフローチャート。
【図１２】第2実施形態による処理装置のS4220の処理の説明図。
【図１３】第2実施形態による処理システムの機能構成図。
【図１４】第2実施形態による処理装置の変形推定フェーズの処理の説明図。
【図１５】第3実施形態による処理システムの機能構成図。
【図１６】第3実施形態による学習フェーズの処理手順を示すフローチャート。
【図１７】第3実施形態による処理装置のS720の処理手順を示すフローチャート。
【発明を実施するための形態】
【００１１】
以下、添付図面に従って本発明に係る処理装置及び方法の好ましい実施形態について詳説する。ただし、発明の範囲は図示例に限定されるものではない。
【００１２】
（第1実施形態）
［第1実施形態の概要］
本実施形態による処理装置は、異なる体位（撮影時の体位に応じた第1の変形条件および第2の変形条件）で撮影した医用画像（第1の形状データおよび第2の形状データ）が取得された場合に、その変形を推定し画像間の位置合わせを行うものである。なお、以下、仰臥位および伏臥位で撮像した乳房のMRI画像を処理対象とする場合について説明するが、被検体の体位やモダリティの種類はこれに限定されない。
【００１３】
本実施形態による処理装置は、多数の被検体における伏臥位から仰臥位への体位の変化に伴う乳房の変形に関する情報を、両体位で撮影した画像の対に基づいて取得する。ここで乳房の変形に関する情報とは、伏臥位における乳房の形状に関する情報と、仰臥位における乳房の形状の変位に関する情報とから成る。本実施形態では、これらの情報を事例データと称する。そして、取得した多数の事例データを主成分分析などの統計解析を施すことにより、乳房に関する変形の統計モデルを生成する。この処理を行うフェーズを学習フェーズと称する。本実施形態ではさらに、事例データには含まれない症例（未知症例）の乳房の変形を、学習フェーズで生成した統計モデルを用いて推定する処理を行う。この処理を行うフェーズを変形推定フェーズと称す。学習フェーズと変形推定フェーズにより、未知症例の伏臥位画像と仰臥位画像との間の対応関係を明らかにできる。また、その対応関係をユーザが視認できるように表示することができる。なお、本実施形態では、三次元医用画像の一例としてMRIによる画像を用いる場合を例として説明するが、本発明の実施形態はこれに限らず、X線CT画像やPET画像などであってもよい。なお、以下MRI画像を、単に画像と称する。
【００１４】
［機能構成］
図1は、本実施形態による処理システムの機能構成図である。本実施形態の処理装置100は、学習フェーズの処理を行う機能と、変形推定フェーズの処理を行う機能により構成される。学習フェーズの処理を行う機能として、事例データ取得部1000、変形取得部1010、メッシュ生成部1020、統計モデル生成部1030の要素がある。一方、変形推定フェーズの処理を行う機能として、対象画像取得部1100、対象形状取得部1110、変形パラメータ算出部1120、観察画像生成部1140の要素がある。また、処理装置100は、データサーバ120、モニタ160に接続される。
【００１５】
画像撮影装置110は、複数症例についての、人体である被検体の内部の3次元領域に関する情報を核磁気共鳴法により取得する。すなわち、画像撮影装置110は画像を取得する。画像撮影装置110は、取得した画像をデータサーバ120へ送信する。データサーバ120は、画像撮影装置110が撮影した画像を保持する。次に、処理装置100を構成する各要素について説明する。事例データ取得部1000は、画像撮影装置110が撮影した画像群（事例データ）を、データサーバ120を介して取得し、処理装置100へ取り込む。変形取得部1010は、事例データ取得部1000が取得した事例データから、各被検体の乳房の伏臥位および仰臥位の間の変形情報を取得する。メッシュ生成部1020は、事例データ取得部1000が取得した事例データと、変形取得部1010が取得した変形情報とに基づいて、各被検体の伏臥位および仰臥位の乳房のメッシュを生成する。
【００１６】
統計モデル生成部1030は、メッシュ生成部1020が生成した乳房のメッシュに基づき、多数の被検体について、乳房の形状と変形に関する統計モデルを導出する。対象画像取得部1100は、処理装置100が変形推定の対象とする被検体（対象症例）に関する画像をデータサーバ120から取得し、処理装置100へと取り込む。対象形状取得部1110は、対象画像取得部1100が取得した対象症例の画像に基づき、対象症例の伏臥位および仰臥位についての、体表形状や乳頭位置に関する情報を取得する。変形パラメータ算出部1120は、統計モデル生成部1030が導出した統計モデルと、対象形状取得部1110が取得した体表形状や乳頭位置に関する情報とに基づき、対象症例の伏臥位から仰臥位への変形を前述の統計モデルで近似するための変形パラメータを算出する。
【００１７】
観察画像生成部1140は、対象画像取得部1100が取得した伏臥位および仰臥位の乳房の画像と、統計モデル生成部1030が生成した統計モデル、および変形パラメータ算出部1120が算出した変形パラメータとに基づき、観察画像を生成する。そして、観察画像生成部1140は、その観察画像をユーザに提示するためにモニタ160へと出力する。モニタ160は観察画像生成部1140が生成した観察画像を表示する。
【００１８】
［装置構成］
図2は、本実施形態による処理システムの装置構成を示す図である。本実施形態の処理システムは、処理装置100、画像撮影装置110、データサーバ120、モニタ160、マウス170、キーボード180により構成される。処理装置100は、例えばパーソナルコンピュータ（PC）などで実現することができる。
【００１９】
処理装置100は、中央演算処理装置（CPU）211、主メモリ212、磁気ディスク213、表示メモリ214を有する。CPU211は、主として処理装置100の各構成要素の動作を制御する。主メモリ212は、CPU211が実行する制御プログラムを格納する。また、主メモリ212は、CPU211によるプログラム実行時の作業領域として使用される。磁気ディスク213は、オペレーティングシステム（OS）、周辺機器のデバイスドライブ、各処理等を行うためのプログラムを含む各種アプリケーションソフト等を格納する。表示メモリ214は、モニタ160のための表示用データを一時記憶する。モニタ160は、例えばCRTモニタや液晶モニタ等であり、表示メモリ214からのデータに基づいて画像を表示する。マウス170及びキーボード180は、ユーザによるポインティング入力及び文字やコマンド等の入力を夫々行う。そして、これら各構成要素は共通バス218により互いに通信可能に接続されている。
【００２０】
［処理フロー］
次に、本実施形態による処理装置100が行う処理に関して、図3および図4のフローチャートを用いて詳しく説明する。図3は本実施形態において処理装置100が実行する学習フェーズの処理のフローチャートである。また、図4は同様に変形推定フェーズの処理のフローチャートである。学習フェーズでは、伏臥位および仰臥位（第1の変形条件および第2の変形条件）の乳房の画像群（第1の形状データおよび第2の形状データ）に基づいて、乳房の形状と変形に関する統計モデルを生成する処理を行う。一方、変形推定フェーズでは、生成された統計モデルに基づいて、未知症例の変形を推定する処理を行う。なお、統計モデルの具体的な説明は、本実施形態の処理の説明の中で詳しく述べる。
【００２１】
本実施形態では、主メモリ212に格納されているプログラムをCPU211が実行することにより実現される。また、以下に説明する処理装置100が行う各処理の結果は、主メモリ212に格納され、記録される。まず、図3を用いて本実施形態の学習フェーズの処理について説明する。
【００２２】
(1)事例データの取得（S300）
本処理ステップの前に、画像撮影装置110が撮影した乳房の画像群（事例データ）がデータサーバ120に格納されているものとする。本処理ステップでは、事例データ取得部1000は、格納されている乳房の画像群（事例データ）を処理装置100に取り込む。事例データは、一乳房につき伏臥位および仰臥位の異なる2種類の体位で撮影した画像を含む。また、事例データの画像は、夫々3次元の座標系（画像座標系）によって表現される。本実施形態では、人体側方をX軸とし、人体の左側の方向をX軸の正方向とする。また、人体前後方向をY軸とし、人体の後ろ（背中）方向をY軸の正方向とする。そして、体軸方向をZ軸とし、人体の頭方向をZ軸の正方向とする。
【００２３】
伏臥位画像と仰臥位画像の間の座標は、乳房ごとに位置合わせが行われている。すなわち、乳房付近で、かつ変形による位置の変動の小さい胸骨体、肋骨、背骨、大胸筋面などの形状を基準として座標の調節が行われている。これらの形状は、同一乳房の伏臥位画像と仰臥位画像との間で略一致しているため、基準として適当である。なお、このような位置合わせを本処理ステップで実行するようにしても良い。いずれの場合でも、位置合わせは周知の方法（手動または自動）で実行できる。本実施形態では、このような事例データの取得処理をN個の乳房について実行する例を説明する。また、本実施形態における事例データは人体の右乳房に関するものとする。
【００２４】
(2)変形の取得（S310）
本処理ステップでは、変形取得部1010は、N個の事例データ（S300）に基づき、各乳房における伏臥位および仰臥位の間の変形情報を取得する。この情報は、伏臥位画像における任意の座標値と、これに対応する仰臥位画像の座標値との間の関係を意味する。例えば、伏臥位画像における病変部や血管の分岐などの部位の座標と、仰臥位画像におけるこれらの部位の座標との関係である。すなわち、数1に示すように伏臥位画像における任意の座標値xprone,iに基づき、それに対応する仰臥位画像の座標値xsupine,iを得る関数fiとして表現可能である。
【数１】
2014018376.tif 000003
ただし、添え字のiは事例データの乳房のインデックスであり、本実施形態では1≦i≦Nである。S310の処理について、図5のフローチャートと図6を用いて具体的に説明する。
【００２５】
(2-1) 対応点の取得（S3100）
本処理ステップでは、変形取得部1010は、伏臥位画像600および仰臥位画像610をモニタ160などに表示してユーザに提示する（図6）。ここで、夫々の画像には、乳腺の輪郭601や血管603が描出される。ユーザは、乳腺構造における特徴的な部位602や血管の分岐部604などを、伏臥位画像と仰臥位画像中での対応する領域（対応点）として視認できる。そして、ユーザは、マウス170やキーボード180を操作すること（以下「入力操作」という。）により、これらの位置座標を複数個取得する。
【００２６】
(2-2) 変形算出（S3110）
本処理ステップでは、変形取得部1010は、伏臥位画像と仰臥位画像における複数の対応点の位置座標（S3100）に基づいて、数1に示した関数fiを求める。この処理は、空間中に離散的に存在する複数の対応点に基づいて、連続的な変形を算出するものである。この処理は周知の補間方法を用いて行われる。例えば、放射基底関数(Radial Basis Function)を用いた方法や、FFD（Free Form Deformation）法と呼ばれるB-スプラインを用いた方法がある。この処理の結果、数1に示した関数fiが求められる。そして変形取得部1010は、求めた関数fiの情報を主メモリ212に記録する。
【００２７】
(2-3) 変形画像の生成（S3120）
本処理ステップでは、変形取得部1010は、関数fiを用いて、伏臥位画像を変形させた仮想的な仰臥位画像を生成する。この処理は伏臥位画像を構成する各画素の画像座標に対して、関数fiを適用することにより実行できる。
【００２８】
(2-4) 対応点取得の終了判定（S3130）
本処理ステップでは、変形取得部1010は、S310の処理を終了するか否かの判定を行う。終了すると判定した場合にはS310の処理を終了する。一方、終了しないと判定した場合には処理をS3100へと戻し、対応点の追加取得を行う。具体的には、まず変形取得部1010が、仮想的な仰臥位画像（S3120）と、被検体の仰臥位画像（S300）とをモニタ160に表示する。そして、ユーザによる入力操作により、処理終了か否かを判定する。ユーザは、仮想的な仰臥位画像が被検体の仰臥位画像に近くなるように、十分な対応点を取得したかを判断し、処理終了か否かを切り替えることができる。以上が本実施形態のステップ310の処理の説明である。なお、S310の処理はN個の事例データの全てに対して上記の処理を実行する。
【００２９】
(3) 伏臥位メッシュの生成（S320）
本処理ステップでは、メッシュ生成部1020は、事例データの伏臥位画像（S300）にメッシュ構造を生成する。図7はS320の処理を詳しく説明するフローチャートである。以下、この図を用いて説明する。
【００３０】
(3-1) 体表形状の取得（S3200）
本処理ステップでは、メッシュ生成部1020は、伏臥位画像に描出される被検体の体表の形状を取得する。この処理は、例えば、ユーザによる入力操作によって取得することで実行できる。この処理について図8を使って具体的に説明する。
【００３１】
本処理ステップでは、まずメッシュ生成部1020は、被検体の伏臥位の画像をモニタ160に表示してユーザに提示する。図8(a)は、モニタ160に表示した被検体の伏臥位画像800を示している。伏臥位画像800において、体表801は、乳房領域802と空気領域803との境界である。ユーザはモニタ160に表示された伏臥位画像800を目視することで、体表801を視認する。そして、メッシュ生成部1020は、ユーザによる入力操作によって、体表点群805を取得する。体表点群は、図8(b)に示すような体表上の点の集まりである。次に、メッシュ生成部1020は、この体表点群の位置をスプライン等の方法を用いた補間方法により連続的な曲面情報に変換する。そして、メッシュ生成部1020は、その情報を主メモリ212に記録する。
【００３２】
(3-2) 大胸筋面形状の取得（S3210）
本処理ステップでは、メッシュ生成部1020は、伏臥位画像に描出される被検体の大胸筋表面の形状を取得する。この処理は、S3200と同様な方法で実行できる。すなわち、メッシュ生成部1020は、ユーザの入力操作によって胸筋面804の点群806を取得し、補間処理によって連続的な曲面情報に変換する。そして、メッシュ生成部1020は、その情報を主メモリ212に記録する。
【００３３】
(3-3) 乳頭位置の取得（S3220）
本処理ステップでは、メッシュ生成部1020は、伏臥位画像に描出される被検体の乳頭の位置を取得する。この処理は、S3200やS3210と同様な方法で実行できる。すなわち、メッシュ生成部1020は、モニタ160に画像を表示し、ユーザの入力操作によって乳頭位置を取得する。そしてメッシュ生成部1020は、その情報を主メモリ212に記録する。
【００３４】
(3-4) バウンダリボックスを取得（S3230）
本処理ステップでは、メッシュ生成部1020は、S3200〜S3220で取得した情報に基づいて、伏臥位画像におけるバウンダリボックスを取得する。具体的には、事例データの乳房の上端、下端、内端、外端、前端、後端の位置を取得し、それらで囲まれた直方体領域を取得する。まず、乳房の上端、下端の位置の取得の処理を行う。この処理について図9を用いて説明する。
【００３５】
図9(a)は、伏臥位画像の乳頭位置（S3220）を含むサジタル面（Y-Z平面）900である。メッシュ生成部1020は、この画像をモニタ160に表示し、上端位置901、下端位置902をユーザによる入力操作によって取得する。取得した上端位置901の座標値をzUPPER,iとして主メモリ212に記録する。同様に下端位置902をzLOWER,iとして主メモリ212に記録する。ここで添え字のiは事例データのインデックスであり、zLOWER,i＜zUPPER,iである。
【００３６】
次に、乳房の内端、外端の位置の取得の処理を行う。図9(b)は、伏臥位画像の乳頭位置（S3220）を含むアキシャル面（X-Y平面）910である。メッシュ生成部1020は、この画像をモニタ160に表示し、内端位置911と外端位置912をユーザによる入力の操作によって取得する。取得した内端位置をxINNER,iとしてとして主メモリ212に記録する。同様に、外端位置をxEXTERIOR,iとして主メモリ212に記録する。人体の右乳房を事例データとする本実施形態では、xEXTERIOR,i＜xINNER,iである。
【００３７】
以上の処理で取得した事例データの上端、下端、内端、外端の位置に基づいて、X-Z平面上での2次元のバウンダリボックスの領域を取得する。この領域は、zUPPER,i≦z≦zUPPER,iかつ、xEXTERIOR,i≦x≦xINNER,iとなるx-z平面上の領域である。
【００３８】
次に、事例データの前端、後端の位置の取得の処理を行う。前端の位置は、体表点群（S3200）のうち、バウンダリボックス領域内におけるY軸上の最小値であり、周知の方法で取得する。取得した前端の位置をyANTERIOR,iとして、主メモリ212に記録する。また、後端の位置は、大胸筋面の点群（S3210）のうち、バウンダリボックスの領域内におけるY軸の最大値であり、周知の方法で取得する。取得した後端の位置をyPOSTERIOR,iとして、主メモリ212に記録する。以上の処理により取得した、事例データの上端、下端、内端、外端、前端、後端で囲まれた3次元の矩形領域が、本ステップで取得するバウンダリボックスとなる。
【００３９】
(3-5) 伏臥位のメッシュ生成（S3240）
本処理ステップでは、メッシュ生成部1020は、バウンダリボックス（S3230）の内部にメッシュを生成する。
【００４０】
図10はS3240の具体的な処理の一例について説明する図である。本処理ステップでは、バウンダリボックスを、乳頭位置951を中心として、X軸、Z軸の各方位に所定数で分割する。図10(a)に示す例では、乳頭位置を中心とし、X軸正方向およびX軸負方向に夫々4分割されている。また、同様にZ軸正方向およびZ軸負方向に夫々4分割されている。これにより、バウンダリボックス950は、X軸分割面952およびZ軸分割面953によって、X-Z平面上を16分割されることとなる。
【００４１】
図10(b)は、図10(a)の複数のZ軸分割面953の任意の一つの面(X-Y平面)を表す図である。この図に示すように、X-Y平面上において、X軸分割面952と体表形状954が交差する位置に体表ノード955を設定する。同様に、X軸分割面952と大胸筋面形状956が交差する位置に大胸筋面ノード957を設定する。そして、同じX軸分割面952上の体表ノード955と大胸筋面ノード957との間を内部ノード958により所定の分割数で分割する。図10(b)では、体表ノード955と大胸筋面ノード957との間を4分割した例を示している。全てのZ軸分割面において同様に処理を行う。そして、メッシュ生成部1020は、生成した体表ノード、大胸筋面ノード、内部ノードの位置を主メモリ212に記録する。
【００４２】
本実施形態では、図10のように分割すると、体表ノードが81個、大胸筋面ノードが81個、内部ノードが243個となる。全ノードの数をMと表記すると、本実施形態ではM=405である。これらのノード位置は、ラスタスキャンの順に番号（ノード番号）を付して記録する。本実施形態では、ノード位置を、xpnode,i,jと表記する。ここで、添え字のiは事例データのインデックス（1≦i≦N）、jはノード番号（1≦j≦M）である。また、これらノード間の隣接関係も記録する。以上の処理で記録したノードの位置とノード間の隣接関係の情報をメッシュ情報と称する。また、隣接関係にある8個のノード点を結んで作られる6面体をメッシュ要素と称する。すると、本実施形態では、64個のメッシュ要素が作られることになる。
【００４３】
本処理ステップでは、さらに、乳頭位置を表すノード番号をnnippleとして記録する。また、乳房の内端に位置する複数の体表ノードの番号をninner,kとして記録する。ここで添え字のkは複数のノード番号のインデックスであり、1≦k≦Kである。ただし、Kは乳房の内端に位置する体表ノードの数であり、本実施形態ではK=9とする。
【００４４】
以上が本実施形態のS320の処理の説明である。なお、本実施形態では、N個の事例データの全てに対して、本処理ステップの処理を実行する。この結果、事例データ毎に作られるメッシュの情報は、ノードの数やメッシュ要素の数が等しい。また、同じノード番号のノードの位置は、事例データの間で解剖学的にほぼ同じ位置を表すことになる。
【００４５】
(4) 仰臥位メッシュの算出（S330）
本処理ステップでは、メッシュ生成部1020は、仰臥位のメッシュのノード位置xsnode,i,jを算出する。この算出処理は、変形情報（関数fi）と、伏臥位のメッシュのノード位置情報xpnode,i,jとに基づく。具体的には、数2の計算をN個の事例データの全てについて全てのノードに対して実行される。
【数２】
2014018376.tif 000004
【００４６】
(5) 事例ベクトル群生成（S340）
本処理ステップでは、統計モデル生成部1030は、伏臥位および仰臥位のメッシュ情報に基づいて、事例ベクトルを生成する。本実施形態では、各事例の乳房の伏臥位および仰臥位のメッシュのノードの位置座標値を並べたベクトル組（第1のベクトルおよび第2のベクトル）を事例ベクトルと称する。本実施形態では、伏臥位および仰臥位のメッシュのノードの数は夫々M個であり、夫々のノードの位置座標は3次元ベクトルである。したがって、事例ベクトルは数3に示すような6×M次元のベクトルとなる。
【数３】
2014018376.tif 000005
統計モデル生成部1030は、N個の事例データの全てに対して上記の処理を実行し、N個の事例ベクトル群を生成する。
【００４７】
(6) 統計モデル導出（S350）
本処理ステップでは、統計モデル生成部1030は、事例ベクトル群xsample,iを統計的に解析して、統計モデルを導出する。事例ベクトル群を統計的に解析する方法としては、例えば主成分分析を用いる。この場合、事例ベクトル群xsample,iの平均ベクトルxaverageおよび主成分ベクトルed(1≦d≦D)を算出する。この２つのベクトルを統計モデルと称する。ここで、Dは主成分分析による算出する主成分の数であり、累積寄与率に閾値を設ける方法や、固定値を使うなどの方法により決定する。
【００４８】
ここで統計モデルについて説明を加える。本実施形態で導出した統計モデルは、乳房の形状と変形の統計的な記述の意味を持つ。例えば、事例データxsample,i(1≦i≦N)は、数4に示すように平均ベクトルxaverageと主成分ベクトルedの加重和（線型結合）により近似できる。
【数４】
2014018376.tif 000006
ここで、cd,i(1≦d≦D, 1≦i≦N)は主成分ベクトルの加重和に関する係数である。数4は数5に表記し直すことができる。
【数５】
2014018376.tif 000007
【００４９】
ここで、Eは主成分ベクトルを並べた行列、すなわちE={e1,e2,,,eD}Tであり、これを主成分行列と称する。またciは数5に示した係数cd,iを添え字のdについて並べたベクトル、すなわちci={c1,c2,,,,cD}Tであり、このベクトルを係数ベクトルと称する。本実施形態で導出する統計モデルは、この係数ベクトルを様々に変えることにより、事例データを含む様々な乳房の形状と変形を表現できる。
【００５０】
本処理ステップでは、次に、統計モデルで近似するための係数ベクトルciの算出を行う。主成分ベクトルcd,iは、互いに正規直交の関係にある。したがって、数6の計算により、事例ベクトルxsample,iを近似するための係数ベクトルを算出できる。
【数６】
2014018376.tif 000008
係数ベクトルは全ての事例データxsample,i (1≦i≦N)に対して算出される。そして、その結果として得た係数ベクトルci(1≦i≦N)を主メモリ212に記録される。
【００５１】
このように、S300〜S350の処理により、本実施形態の学習フェーズの処理が実行される。なお、上記の説明では、S300〜S330の各処理は全事例データを対象とする例を説明したが、本発明の実施形態はこれに限らない。例えば、S300〜S330までの各処理を一つの事例データを対象とした処理としてもよい。この場合、この一連の処理を事例データ数分繰り返す。繰り返しの回数は最大で事例データの数（N個）となるが、その途中段階で処理を打ち切り、S340の処理を実行しても良い。
【００５２】
例えば、後に説明するS400以降の処理を行うために十分な事例データが得られた場合に、処理を打ち切るようにできる。この時、打ち切りの判断は、ユーザによる入力操作や、導出された統計モデルの特性等に基づいて行うことができる。後者の場合は、導出された統計モデルの事例データに対する近似性能を検証すること等により行われる。例えば、Leave-one-out方式などを用いて検証し、その結果に基づいて自動的に判断を行うようにできる。次に図4を用いて、本実施形態における変形推定フェーズの処理について説明する。
【００５３】
(1) 対象症例データ読み込み（S400）
本処理ステップでは、対象画像取得部1100は、変形推定処理の対象とする症例の乳房の画像を、データサーバ120から取得する。このような症例を対象症例と称する。対象症例は、S300で画像を取得した症例とは異なり、変形の推定を行う対象となる症例である。
【００５４】
本処理ステップでは、対象症例の伏臥位および仰臥位の乳房の画像（第1の対象形状データおよび第2の対象形状データ）を取得する。これらの画像は、S300と同様に、胸骨などの変動の小さい部位等を基準として位置合わせが行われている。このような位置合わせを本処理ステップで実行し、その結果を取得してもよい。なお、本実施形態では、対象症例の右乳房を対象とする場合を例として説明する。
【００５５】
(2) 乳頭・体表形状データ抽出（S410）
本処理ステップでは、対象形状取得部1110は、対象症例の画像（S400）に基づき、乳頭位置および体表形状を取得する。この処理は、S3200とS3220と同様の方法で実現できる。この結果、対象症例の乳頭位置の3次元座標値と、体表形状を表す点群または陰関数などの情報を取得する。
【００５６】
本処理ステップでは、これ以外にも画像処理を行うことで、乳頭位置および体表形状を自動的に取得するようにしても良い。この場合、対象形状取得部1110は、まず、画像から周知の方法によってエッジ検出を行う。そして、3次元空間上の体表形状を曲面として抽出し、その体表形状の曲率を使って乳頭の位置を検出するようする。
【００５７】
本処理ステップにおいて、対象形状取得部1110はさらに、対象症例の伏臥位における正中線の位置の取得を行う。ここで、正中線とは、対象症例の左右の中央を通る体表上の曲線である。正中線は取得した対象症例の体表形状の一部であり、3次元空間における点群や関数などの形式で取得する。正中線の位置はユーザの入力操作に基づいて取得しても良いし、取得した体表形状や画像の輝度値のX軸（人体の左右）方向の対称性を利用して取得しても良い。この場合、対象形状取得部1110は、例えば、画像のX軸方向の対称面と交差する体表形状を正中線の位置として取得することができる。
【００５８】
(3) 最適パラメータの導出（S420）
本処理ステップでは、変形パラメータ算出部1120は、対象症例の形状に対して、統計モデルが表す形状が略一致するような統計モデルの係数と剛体変換を算出する。本実施形態では、この統計モデルの係数と剛体変換とを変形パラメータと称する。本処理ステップでは適切な変形パラメータを算出する。本処理ステップについて図11に示すフローチャートを使って説明する。
【００５９】
(3-1) パラメータの初期化（S4200）
本処理ステップでは、変形パラメータ算出部1120は、統計モデルの係数に関して初期値c’を設定する。初期値c’は、例えばゼロベクトルとすることができる。本実施形態では、S4210〜S4240までの処理を繰り返し実行することにより、初期値c’をより最適な値へと更新する。すなわち、本処理ステップで設定する初期値c’は統計モデルの係数の暫定値である。
【００６０】
(3-2) パラメータ候補の生成（S4210）
本処理ステップでは、変形パラメータ算出部1120は、係数の暫定値c’に基づいて、係数の変動の候補c”hを生成する。ただし、hは係数の変動の候補のインデックスであり、(1≦h≦H)である。ここで、Hは係数の変動の候補の総数である。c”hの具体な生成方法は、例えば数7に示すように、c’に基づき、各次元の要素に変動を与えるように生成できる。この場合、係数の変動の候補の総数は2D+1となる。Dは前述したように、主成分分析による算出する主成分の数である。
【数７】
2014018376.tif 000009
ここで、δは変動の幅であり、所定の定数が設定されるものとする。
【００６１】
次に、係数の変動の候補c”hの夫々に対し、数5に示した統計モデルを用いて、形状x”h(1≦h≦H)を数8により算出する。
【数８】
2014018376.tif 000010
本実施形態では、形状x”h(1≦h≦H)を暫定モデル形状群と称する。
【００６２】
(3-3) 剛体変換の算出（S4220）
本処理ステップでは、変形パラメータ算出部1120は、暫定モデル形状群（S4210）と対象症例の体表形状及び正中線（S410）に基づいて、剛体変換を表す行列を算出する。まず、変形パラメータ算出部1120は、伏臥位の乳頭位置と、暫定モデル形状群の夫々のnnipple番目のノード位置とに基づいて、これらを一致させる剛体変換（並進）を求める。その後、この剛体変換を暫定モデル形状群の夫々に適用する。この剛体変換を表す4行4列の行列を剛体変換行列m1,hと称する。そして、暫定モデル形状群x”hを剛体変換行列m1,hによって変換し、暫定モデル形状群x”m1,hを算出する。ただし、hは暫定モデル形状群のインデックスであり、1≦h≦Hである。
【００６３】
図12(a)は剛体変換行列m1,hによる剛体変換を適用した暫定モデル形状群x”m1,hと、対象症例の体表形状および正中線の位置の関係を説明する図である。この図において、暫定モデル形状の体表形状1200は、H個ある暫定モデル形状群x”m1,hの中の任意の一つに関する体表のノード位置を表している。図中では、ノード位置を分かりやすく示すために、隣接するノード間を線分で結んで図示している。対象症例の体表形状1210は、S410で取得した対象症例の体表形状である。暫定モデル形状の体表形状1200と対象症例の体表形状1210は、夫々3次元の曲面であるが、ここでは2次元の曲線として示している。
【００６４】
暫定モデル形状の正中線1220は、暫定モデル形状の体表形状1200の内端の位置である。また、対象症例の正中線1230は対象症例の体表形状1210の正中線の位置である。これら正中線は、3次元空間中における曲線であるが、ここでは点として表示している。
【００６５】
次に、変形パラメータ算出部1120は、図12(b)に示すように暫定モデル形状の正中線1220と対象症例の正中線1230を略一致させる剛体変換を求める。この剛体変換を表す行列を剛体変換行列m2,hと称する。そして、暫定モデル形状x”m1,hを剛体変換行列m2, hによって変換し、暫定モデル形状x”hを算出する。剛体変換行列m2,hは、乳頭位置を中心とした回転を表す剛体変換行列である。剛体変換行列m2,hは、ICP(Iterative Closest Point)法など周知の方法で算出できる。
【００６６】
本処理ステップでは、剛体変換m1,hとm2,hを用いて剛体変換後の暫定モデル形状x”hを算出する。また、暫定モデル形状x”hの正中線と対象症例の正中線との距離（残差）erigid,hを算出して、主メモリ212に記録する。
【００６７】
(3-4) 評価関数の計算（S4230）
本処理ステップでは、変形パラメータ算出部1120は、暫定モデル形状x”hと、対象症例の体表形状および乳頭位置（S410）とに基づいて、数9に示す評価関数を用いた評価値の算出を行う。ここで評価関数とは、暫定モデル形状x”hと対象症例の体表形状および乳頭位置との一致度を評価する関数である。
【００６８】
まず、暫定モデル形状x”hと対象症例の伏臥位の体表の一致度に関する評価値eprone_surface,hを算出する。この処理は、例えば、暫定モデル形状x”hの中の伏臥位の体表を表すノード位置の夫々から、対象症例の伏臥位の体表形状への最近傍探索を行い、その最近傍となる体表形状への距離の平均を算出することにより実現できる。
【００６９】
次に、暫定モデル形状x”hと対象症例の仰臥位の体表の一致度に関する評価値esupine_surface,hを算出する。この処理は、例えば、暫定モデル形状x”hの中の仰臥位の体表を表すノード位置の夫々から、対象症例の仰臥位の体表形状への最近傍探索を行い、その最近傍となる体表形状への距離の平均を算出することにより実現できる。
【００７０】
さらに、暫定モデル形状x”hと対象症例の仰臥位の乳頭位置の一致度に関する評価値esupine_nipple,hを算出する。この処理は、暫定モデル形状x”hの中の仰臥位の乳頭を表すノード位置と、対象症例の仰臥位の乳頭位置との距離を算出することにより実現できる。そして、数9のように、算出した各評価値と、正中線の位置に関する残差（S4220）により評価関数を生成し、評価値ehを算出する。
【数９】
2014018376.tif 000011
【００７１】
(3-5) 最適化終了判定（S4240）
本処理ステップでは、変形パラメータ算出部1120は、S420の処理を終了させるか否かの判定を行う。この判定は、H個の評価値ehに基づいて、例えば、以下の処理により実行される。まず、変形パラメータ算出部1120は、H個の評価値ehの中で最も値が小さいものを選択する。これをeoptとして主メモリ212に記録する。また、対応する統計モデルの係数を、係数の変動の候補c”hの中から選択し、coptとして主メモリ212に記録する。同様に、剛体変換についてもm1,opt,m2,optを選択して主メモリ212に記録する。eoptが所定の閾値よりも小さい場合には、S420の処理は終了する。一方、そうでない場合にはcoptを新たなc’として設定し、処理はS4210に戻る。
【００７２】
このように、変形パラメータ算出部1120は、最適な変形パラメータを算出した場合にS420の処理を終了する。しかし、最適な変形パラメータを算出していない場合には、より良い変形パラメータを求めるために、処理を反復する。
【００７３】
処理を反復するための判定方法は、上記の例に限らない。例えばS4210〜S4240のループ処理を複数回終えた段階において、前回までに選択された評価値より今回算出された評価値が所定の閾値よりも小さい場合に、S420の処理は終了しても良い。また、前回までに選択された評価値よりも今回算出された評価値への減少の割合が所定の閾値よりも小さい場合に、S420の処理は終了しても良い。また、ループ処理が所定の回数以上実行された場合には、S420の処理は終了しても良い。すなわち、非線形最適化に関する周知の打ち切り処理を用いることができる。
【００７４】
なお、上記に説明した変形パラメータの算出処理の方法は、本発明の実施形態の一例であり、これ以外の方法であっても良い。例えば、S420のパラメータの最適化の処理は、数9に示した評価関数に対する非線形問題と捉えることができる。したがって、これを他の周知な非線形最適化の手法に替えることができる。すなわち、ランダムサンプリング法や最急降下法、ニュートン法、シミュレーテッドアニーリング法、遺伝的プログラミング（GA）法などである。また、これらを組み合わせて用いるようにしても良い。例えば、ランダムサンプリング法で求めた変形パラメータを初期値とし、それに基づいて最急降下法を実行するようにしても良い。また、GA法の後に最急降下法やニュートン法を使うようにしても良い。また、変形パラメータを粗密探索するようにしても良い。
【００７５】
(4) 変形場の算出（S430）
本処理ステップでは、変形画像生成部1130は、統計モデル（S340）と変形パラメータ（S420）とに基づいて、対象症例の変形場を算出する。ここで変形場とは、対象症例の乳房における伏臥位から仰臥位への変形の場である。具体的には、まず、数10に示す計算を行う。
【数１０】
2014018376.tif 000012
【００７６】
次に、数10で算出したベクトルxoptを2×M個の3次元ベクトルに分解する。そして、その3次元座標ベクトルに対して剛体変換行列m1,opt,m2,opt（S420）による変換を行う。ここで1番目からM番目までのM個の3次元座標ベクトルは、対象症例の伏臥位の乳房に対応するメッシュのノード位置を表す。また、M+1番目から2×M番目までのM個の3次元座標ベクトルは、対象症例の仰臥位の乳房に対応するメッシュのノード位置を表す。そしてj(1≦j≦M)番目の伏臥位のノード位置とj+M番目の仰臥位のノード位置は、対象症例の同一部位に対応する推定を表すことになる。すなわち、これらの3次元ベクトルの差分が対象症例の伏臥位から仰臥位への変形の推定値を表すことになる。
【００７７】
本処理ステップでは、さらに、剛体変換を施した2×M個のノード位置の3次元座標ベクトルに基づいて、対象症例の伏臥位から仰臥位への変形場を算出する。この処理は、S320で行ったメッシュ生成の処理において記録したノードの接続に関する情報に基づいて実行する。すなわち、メッシュ要素毎に、その要素を構成する8個のノードの位置の情報を用いて、メッシュ要素の内部の変位を線形補間などにより求める。この方法は、メッシュモデルを用いた変形表現に関する周知の方法により実行できる。この結果、対象症例に関する変形を表す関数foptが求められる。ここで、関数foptはS310で説明した事例データに関する変形の表現と同様の形式の関数である。すなわち、対象症例の伏臥位の3次元の画像座標を引数として、対象症例の仰臥位の3次元の画像座標への変換を表す関数である。本処理ステップでは、算出した関数foptを主メモリ212に記録する。
【００７８】
(5) 観察画像の生成（S440）
本処理ステップでは、観察画像生成部1140は、対象症例の伏臥位および仰臥位の画像との間の対応関係を提示する観察画像を生成する。そして、その画像をモニタ160に表示する。
【００７９】
本処理ステップでは、観察画像生成部1140は、まず、伏臥位画像をモニタ160に表示する。そして、その画像中の注目領域をユーザの入力操作等から取得する。ここで、注目領域の伏臥位画像における画像座標値をxproiとする。そして、その座標値を関数fopt（S430）を用いて、対象症例の仰臥位画像における対応する領域の画像座標値xsroiに変換する。すなわち数11の計算を行う。
【数１１】
2014018376.tif 000013
【００８０】
そして、対象症例の仰臥位画像をモニタ160に表示するとともに、算出した画像座標値xsroiの位置にユーザが視認できるマーク等を重畳して表示する。これにより、伏臥位画像上でユーザが注目した領域に対応する領域が、仰臥位画像上に提示される。
【００８１】
なお、伏臥位画像における注目領域の画像座標値xproiは必ずしもユーザの入力操作によって取得する必要はない。例えば、ユーザが伏臥位画像を画像処理することにより病変の疑いのある位置を検出し、その位置をxproiとしても良い。また、上記には、伏臥位の画像の画像座標値から、これに対応する仰臥位の画像の画像座標値を算出する例を示したが、これと逆にしてもよい。すなわち、仰臥位画像の画像座標値から、これに対応する伏臥位画像の画像座標値を算出するようにしても良い。
【００８２】
このように、S400〜S440の処理により、本実施形態の変形推定フェーズの処理が実行される。以上が、第1実施形態における処理装置100の処理の説明である。第1実施形態によれば、対象症例の伏臥位画像と仰臥位画像との間の位置の対応を提示できる仕組みを提供できる。
【００８３】
（変形例1-1）
第1実施形態の説明では、学習フェーズで生成された統計モデルを用いて、変形推定フェーズの処理を行う場合を例として説明したが、本発明の実施形態はこれに限らない。例えば、処理装置100は学習フェーズの処理により統計モデルを導出し、これを主メモリ212等に記録するようにしても良い。このとき、処理装置100は必ずしも変形推定フェーズの処理を行う機能を有する必要は無い。そして、別の処理装置が、その記録された統計モデルを読み込んで変形推定フェーズの処理を実行する。この別の処理装置は、記録された統計モデルを読み込み、それに基づいて変形推定フェーズの処理を実行すれば良い。
【００８４】
（変形例1-2）
第1実施形態の説明では、学習フェーズの処理により単一の統計モデルを導出する場合を例として説明したが、本発明の実施形態はこれに限らない。例えば、乳房のサイズや形状の特徴などに基づいて症例を複数のグループに分類し、グループ毎の症例に基づいて複数の統計モデルを導出するようにしても良い。この場合、変形推定フェーズでは、対象症例の乳房のサイズや形状の特徴などの属性に基づいて、複数の統計モデルから適切な統計モデルを選択して利用するようにできる。この方法を用いれば、対象症例に乳房のサイズや形状の特徴が近い症例に基づいた統計モデルを利用できるため、変形推定の精度が向上する。または、これらの特徴が互いに近い事例データに基づいた統計モデルを利用できるため、より少ない主成分数Dで適切な変形推定を行うことができ、処理効率が向上する。
【００８５】
なお、症例のグループ分けの基準として説明した属性は、その一例に過ぎない。したがって、例えば、乳房に含まれる乳腺の割合や量、被検体の年齢、人種、住む地域、患っている病名や既往歴などといった属性に基づくグループ分けもできる。何れの場合も上記と同様の効果が期待できる。
【００８６】
また、グループ分けは、前述のように予め定めた方法に従っても良いし、他の方法に従っても良い。すなわち、複数の症例データを取得後に適切なグループ分けの方法を選択して、グループ分けを適用しても良い。また、事例データの総数に基づいてグループ数を可変にしたり、グループ分けの可否を決めたりしても良い。
【００８７】
また、事例データに含まれる症例の乳房のサイズなどの分布などに基づいてグループ分けの境界を定めても良い。このように、事例データの数や、事例データに含まれる乳房の特徴の分布などに応じて適応的にグループ分けができる。
【００８８】
従って、病院情報システムなどと接続して事例データを取得するような場合に、事例データの数などに応じた適切な統計モデルを導出できる。これにより、変形推定の精度や処理の効率化を図ることができる。
【００８９】
（変形例1-3）
第1実施形態の説明では、学習フェーズでは単一の主成分数Dを持つ統計モデルを導出する場合を例として説明したが、本発明の実施形態はこれに限らない。
【００９０】
例えば、学習フェーズで行う処理は、主成分数を様々に変えて、各主成分数毎に統計モデルを生成するようにしても良い。この場合、変形推定フェーズでは、複数の統計モデルのうち、まず低い主成分数の統計モデルを用いた処理を行う。そして、次により高い主成分数のモデルを用いて処理を行う。このとき、S4200で設置するパラメータの初期値は、低い主成分数の統計モデルで変形推定フェーズを実行した結果を用いるようにできる。
【００９１】
これによれば、高い主成分数の統計モデルを用いた変形推定フェーズの処理で局所解に陥りにくくできる。そのため、変形推定の高精度化と安定化を図れる効果がある。また、変形推定の処理を効率化できるという効果がある。
【００９２】
（変形例1-4）
第1実施形態の説明では、学習フェーズと変形推定フェーズで人体の右乳房を対象とする場合を例として説明したが、本発明の実施形態はこれに限らない。例えば、学習フェーズでは事例データの右乳房と左乳房とを夫々個別に処理をして夫々の統計モデルを導出するようにしても良い。この場合、変形推定フェーズで変形推定する対象の乳房が右乳房か左乳房かに基づいて統計モデルを切り替えて使うようにできる。
【００９３】
これにより、左右の乳房の夫々に対応した変形推定を実行できる。また、人体の左乳房と右乳房は、ほぼ対称な形状と変形の特性を有している。このため、学習フェーズで左右いずれかの乳房に基づいて生成された統計モデルを対称変換（鏡像）し、変形推定フェーズで用いることができる。すなわち、変形推定フェーズでは、対称変換した統計モデルを用いて、左右いずれの乳房の変形の推定を行うようにできる。これによれば、事例データとして左右いずれかのデータのみを用いて、左右両側の乳房を対象とした変形推定を実行できる。
【００９４】
また、事例データとして左右両側の乳房の事例データを用いる場合には、前述の対称性を利用することができる。すなわち、左右いずれかの事例データを対称変換することで、全事例データを用いて統計モデルを導出するようにできる。例えば、事例データに含まれる左乳房の事例データを対象変換して、右乳房のデータと見立てて使用できる。これによれば、少ない症例数から効率的に事例データを生成できる。また、事例データとして用いることができる症例数に限りがある場合に、より精度の高い変形推定を行うことができる。
【００９５】
（変形例1-5）
第1実施形態の説明では、統計モデルの係数の初期値としてゼロベクトルを設定する場合を例として説明したが、本発明の実施形態はこれに限らない。例えば、S350で算出した事例データの夫々に関する係数ci(1≦i≦N)のいずれか一つを用いるようにしても良い。
【００９６】
この場合、例えば、対象症例と形状が似た事例データを選択し、その事例データに対応する係数を初期値とすることができる。これによれば、変形パラメータの初期値の設定を、より正解に近い値とできることが期待できる。また、変形パラメータの推定処理を効率的に実行できるという効果がある。ここで、対象症例と形状が似た事例データを選択する方法として以下のような方法が考えられる。すなわち、対象症例の体表形状（S410）と、事例データの体表形状（S3200）とを比較して、類似性を評価して選択する方法である。
【００９７】
これ以外にも、事例データをユーザが選択し、その事例データに対応する係数を初期値とする方法もある。これは、対象症例の画像と事例データの複数の画像とをモニタ160に表示し、ユーザが、対象形状と形状が似た事例データを選択することにより実現できる。
【００９８】
（第2実施形態）
［第2実施形態の概要］
本実施形態による処理装置は、未知症例について、伏臥位で撮影したMRI画像と仰臥位で撮影した超音波画像との間の対応関係をユーザが視認できるように表示する。なお、第1実施形態と同様な学習フェーズで生成した統計モデルを利用する。
【００９９】
［機能構成］
図13は、本実施形態による処理システムの構成を示す図である。本実施形態において、第1実施形態と同様の機能を持つ構成要素には図1と同一の番号を付しており、説明は省略する。
【０１００】
本実施形態における処理装置200は、図1と比較してさらに位置姿勢計測装置130、超音波画像診断装置140と接続する。また処理装置200は、伏臥位画像取得部1170、変形画像生成部1130、計測値取得部1150、超音波画像取得部1160を持つ。超音波画像撮影装置140は不図示の超音波プローブを被検体に接触させて、被検体の内部を超音波撮影する。超音波画像撮影装置140は、被検体の断面領域を撮影した2次元のBモード超音波画像を撮影する。
【０１０１】
位置姿勢計測装置130は、3次元空間において、超音波プローブの位置と姿勢を計測する装置である。この装置は、例えば、磁気式や光学式の6自由度計測装置を超音波プローブに装着することで構成される。ここで、超音波プローブの位置と姿勢とは、プローブ座標系とセンサ座標系との間の相対的な位置と姿勢の関係を意味する。ここで、プローブ座標系とセンサ座標系について説明する。なお、以下において、超音波プローブの位置と姿勢の計測値をプローブ計測値と称する。
【０１０２】
プローブ座標系とは、超音波プローブを基準として定める直交座標系である。また、超音波プローブの撮影面（被検体と接する面）の中心を原点とする。そして、超音波ビームの放射方向をY軸、Y軸に直行する方向をX軸とする。ここでX軸は、超音波画像撮影装置140が撮影する超音波画像の撮影面に含まれるものとする。また、X軸とY軸の外積の方向をZ軸とする。
【０１０３】
センサ座標系は、位置姿勢計測装置130が基準として定める直交座標系である。例えば計測を行う3次元空間に固定して設置される不図示の計測ステーションを基準としている。計測値取得部1150は、位置姿勢計測装置130が計測したプローブ計測値を、並進と回転により表される4行4列の剛体変換行列として取得する。ただし、この計測値はセンサ座標系を基準とする。超音波画像取得部1160は、超音波画像撮影装置140が取得した被検体の内部を超音波撮影した超音波画像を取得する。
【０１０４】
［処理フロー］
次に、処理装置200が行う全体の動作に関して説明する。本実施形態では、CPU211が、主メモリ212に格納されているプログラムを実行することにより実現される。また以下に説明する処理装置100が行う各処理の結果は、主メモリ212に格納し、記録される。本実施形態における処理装置200の処理は、第1実施形態と同様に学習フェーズと変形推定フェーズの処理から成る。本実施形態による学習フェーズの処理は、第1実施形態のS300〜S350（図3）と同様であるので、説明は省略する。
【０１０５】
次に、本実施形態における変形推定フェーズの処理について、図14のフローチャートを用いて詳しく説明する。なお、計測値取得部1150は、処理のステップに関わらず、位置姿勢計測装置130が計測したプローブ計測値を逐次取得する。そしてその計測値を保持し、他の処理部からの要求に応じて最新の計測値を出力する。
【０１０６】
ここで、計測値取得部1150が取得するプローブ計測値は、直接的にはセンサ座標系における計測値である。しかし、本実施形態では、この計測値に所定の座標変換を施した値とする。すなわち、対象症例の体を基準とした座標系（人体座標系）への座標変換を施すものとする。
【０１０７】
基準座標系は、例えば対象症例の胸骨体の下端位置である剣状突起の先端位置を原点とし、対象症例の体側方向をX軸、前後方向をY軸、体軸方向をZ軸などと定めることができる。座標変換は、例えば、剣状突起の先端位置などを超音波画像診断装置140で描出した時のプローブ計測値などに基づいて実行できる。この処理は周知の技術で実現できるため詳細な説明は省略する。これにより、計測値取得部1150は、人体座標系を基準としたプローブ計測値として取得する。
【０１０８】
超音波画像取得部1160は、処理のステップに関わらず、超音波画像撮影装置140が撮影した超音波画像を逐次取得し、保持する。また、他の処理部からの要求に応じて最新の超音波画像を出力する。
【０１０９】
(1) 対象症例の伏臥位画像を取得（S500）
本処理ステップでは、対象画像取得部1100は、対象症例の伏臥位画像をデータサーバ120から取得する。ここで、取得する画像は、直接的には画像撮影装置が定める座標系を基準としている。しかし、本実施形態では、この計測値に上記と同様に人体座標系への座標変換を施した値とする。この座標変換は、例えば伏臥位画像に描出される被検体の剣状突起の先端などの位置などに基づいて実行できる。この処理は周知の技術で実現できるため詳細な説明は省略する。
【０１１０】
(2) 伏臥位乳頭位置・体表形状のデータを抽出（S510）
本処理ステップでは、形状取得部1110は、取得した対象症例の画像に基づき、伏臥位の乳頭位置と体表形状を抽出する。この処理は第1実施形態で説明したS410と同様の処理である。この結果、伏臥位の乳頭位置の3次元座標値と、体表形状を表す点群または陰関数などの情報を取得する。
【０１１１】
(3) 仰臥位乳頭位置・体表形状のデータを取得（S520）
本処理ステップでは、形状取得部1110は、計測値取得部1150による計測値に基づいて、対象症例の仰臥位における体表形状と乳頭位置とを取得する。この処理は例えば以下のように実行される。まず、ユーザは対象症例の乳頭位置に超音波プローブの撮影面中心が触れるように保持する。そして、形状取得部1110は、ユーザの実行命令を取得した際に、計測値取得部1150が取得した最新の計測値を取得する。そして、超音波プローブの撮影面の中心位置を対象症例の仰臥位の乳頭位置として記録する。
【０１１２】
同様にして、対象症例の仰臥位における体表形状も取得できる。すなわち、形状取得部1110は、超音波プローブを上記のように保持し、最新の計測値を取得することで、対象症例の仰臥位における体表の位置を取得する。そして、その処理を異なる複数の体表の位置に対して繰り返すことで、体表形状を表す複数の体表の位置を取得する。
【０１１３】
以上の説明では、超音波プローブを用いた場合を例として説明したが、本発明の実施形態はこれに限らない。例えば、位置姿勢計測装置130がスタイラスなど任意の物体の位置を計測し、その物体を対象症例の体表や乳頭の位置に接触させことにより、対象症例の体表形状や乳頭位置を取得するようにできる。
【０１１４】
(4) 最適な変形パラメータの算出（S530）
本処理ステップでは、処理装置200は、伏臥位の体表形状等（S510）と仰臥位の体表形状等（S520）に基づいて、最適な変形パラメータである統計モデルの係数と剛体変換を算出する。この処理は第1実施形態のS420と同様である。しかし、第1実施形態では、対象症例の体表形状等を、伏臥位と仰臥位ともに対象画像取得部1100（図1）において取得した画像に基づいて取得した値を用いるのに対し、本実施形態では以下の違いがある。すなわち、本実施形態では、伏臥位については第1実施形態と同じであるが、超音波プローブを接触させて計測した値を対象症例の仰臥位の体表形状等の情報として取得し、その値を用いる。
【０１１５】
(5) 変形場の算出（S540）
本処理ステップでは、処理装置200は、第1実施形態のS430と同様の処理を実行し、対象症例の変形場を算出する。すなわち、対象症例に関する変形を表す関数foptを算出する。
【０１１６】
(6) 変形画像の生成（S550）
本処理ステップでは、変形画像生成部1130は、伏臥位画像（S500）を、関数fopt（S550）により変形して、変形画像を生成する。ここで、変形画像とは、仰臥位における対象症例と空間的に整合した画像である。以下、変形画像を生成する具体的な処理について説明する。
【０１１７】
対象症例の伏臥位画像（S500）を、関数の形式で、Iprone(x,y,z)と表現する。この関数は、人体座標系の座標値を引数とし、画像の画素値（輝度値）を返す関数である。本処理ステップで生成する変形画像もこれと同様の形式で、Ideformed(x,y,z)と表記する。ここで、変形画像は、伏臥位画像と同様に人体座標系を基準とする。しかし、対象症例の伏臥位から仰臥位への姿勢の変形によって変形が生じる部位においては、座標の変換が必要となる。この座標の変換に関する情報は、S550で算出した関数foptに基づく。つまり、変形画像は、その画像を構成する画素について数12に示す計算を実行する。
【数１２】
2014018376.tif 000014
【０１１８】
ここで、f-1optは関数foptの逆関数であり、対象症例の仰臥位から伏臥位への変形の場である。f-1optはfoptに基づいて探索的に求めることができる。また、S550と同様に、伏臥位および仰臥位に関するノード位置を入れ替えて実行して求めるようにしても良い。以上に説明した処理により変形画像、Ideformed(x,y,z) が生成される。
【０１１９】
以下、本実施形態におけるS560〜S590の処理について説明する。なお、本実施形態におけるS560〜S590の処理は、S590の処理の分岐判断に基づき複数回繰り返して実行される場合がある。
【０１２０】
(7) 超音波画像の取得（S560）
本処理ステップでは、超音波画像取得部1160は、仰向け状態の被検体の乳房を超音波画像撮影装置140が撮影した超音波画像を取得する。超音波画像は、仰臥位の対象症例を超音波撮影した画像であり、本実施形態では二次元のBモード画像である。
【０１２１】
(8) 位置姿勢計測値取得（S570）
本処理ステップでは、計測値取得部1150は、超音波画像（S560）が撮影された時点における、プローブ計測値を取得する。本実施形態における計測値取得部1150は、人体座標系における計測値を取得する。なお、本実施形態において、S570で取得した超音波画像の画像座標系は、前述のプローブ座標系（のZ=0平面）と一致しているものとする。
【０１２２】
(9) 観察画像の生成（S580）
本処理ステップでは、観察画像生成部1140は、超音波画像（S560）、プローブ計測値（S570）、変形画像（S550）に基づき、観察画像を生成する。ここで観察画像とは、処理装置200がモニタ160を介してユーザに提示する画像である。すなわち、対象症例の伏臥位の画像（S500）と超音波画像（S560）とをユーザが容易に比較して観察できるようにした画像である。
【０１２３】
観察画像を生成する具体的な処理について説明する。まず、プローブ座標系を基準としている超音波画像について、プローブ座標系における撮影領域を取得する。例えば、超音波画像が矩形の場合には、その4隅の座標値を取得する。次に、プローブ計測値（S570）に基づいて、この4隅の座標値を、人体座標系を基準とした座標値に変換する。そして、その4隅の座標値を結ぶ人体座標系における矩形領域を算出する。そして、変形画像（S550）から、この矩形領域の画素を切り出すことにより対応断面画像を生成する。そして、この対応断面画像と、超音波画像（S560）とを比較して観察画像を生成し、モニタ160等に表示する。観察画像の具体例としては、対応断面画像と超音波画像とを並べた画像や、両画像の一方に他方を重畳した画像などである。
【０１２４】
これにより、超音波画像と、超音波画像を撮影した対象症例の略同一領域を撮影した画像とを比較して観察できる観察画像が生成される。以上の処理により生成した観察画像を、表示メモリ214への送信を介してモニタ160に表示する。
【０１２５】
(10) 終了判定（S590）
本処理ステップでは、ユーザは処理装置200の処理を終了させるか否かを判定する。終了と判定した場合には処理は終了し、そうでない場合には処理はS560に戻る。この判定は、ユーザによる入力操作に基づいて行うことができる。処理がS560に戻る場合には、ユーザの操作等に応じて、対象症例の異なる領域を撮影した超音波画像が撮影・取得される。そして、それに応じた観察画像が生成される。S560〜S590までの処理が高速に繰り返される場合には、ユーザの操作等に追随して観察画像が表示される。
【０１２６】
以上に説明した本発明の第2実施形態によれば、超音波撮影された仰臥位の対象症例と略一致するように、伏臥位の状態で撮影された対象症例の画像が変形される。そして、仰臥位の超音波画像と比較容易な状態で変形された伏臥位画像が提示される。したがって、対象症例の超音波画像と（MRI）画像とを、空間的な位置の対応付けが容易な形態で表示できる。
【０１２７】
（変形例2-1）
本実施形態では、ユーザが超音波プローブを対象症例に接触させ、その時の計測値を用いる場合を例として説明したが、本発明の実施形態はこれに限らない。例えば、処理装置200は、対象症例の体表形状等を非接触で計測できるレンジセンサと接続することができる。この場合、レンジセンサの計測値に基づいて対象症例の体表形状や乳頭位置を取得するようにできる。これによれば、超音波プローブの操作なく対象症例の体表形状と乳頭位置とを取得できるから、より簡便に処理装置200を使用できる効果がある。
【０１２８】
（第3実施形態）
[第3実施形態の概要]
本実施形態による処理装置は、第1実施形態に記載の処理装置が実行するS310の処理を自動的に実行する。これにより、ユーザが行う入力処理を低減し、より簡便に処理装置を使えるようにする。
【０１２９】
[機能構成]
図15は、本実施形態による処理システムの構成を示す図である。本実施形態において、第1実施形態と同様の機能を持つ構成要素には図1と同一の番号を付しており、説明は省略する。本実施形態における処理装置250は、伏臥位メッシュ生成部1170、仰臥位メッシュ生成部1180を持つ。
【０１３０】
伏臥位メッシュ生成部1170は、事例データ取得部1000が取得した事例データに基づいて、各症例の乳房に関する伏臥位のメッシュを生成する。仰臥位メッシュ生成部1180は、事例データ取得部1000が取得した事例データと、伏臥位メッシュ生成部1170が生成した伏臥位のメッシュに基づいて、各症例の乳房に関する仰臥位のメッシュを生成する。
【０１３１】
[処理フロー]
本実施形態による処理装置250は、第1実施形態の処理装置と同様に学習フェーズの処理と変形推定フェーズの処理により実行される。なお、本実施形態における変形推定フェーズの処理は、第1実施形態と同様であるので、ここでは説明を省略する。以下では、図16を用いて学習フェーズの処理についてのみ、説明する。
【０１３２】
(1) 事例データ取得（S700）
本処理ステップでは、事例データ取得部1000は第1実施形態のS300と同様の処理を実行する。
【０１３３】
(2) 伏臥位メッシュの生成（S710）
本処理ステップでは、伏臥位メッシュ生成部1170は、事例データの伏臥位画像（S700）にメッシュ構造を生成する。この処理は、第1実施形態においてメッシュ生成部1020が実行するS320の処理と同様の処理により実行される。
【０１３４】
(3) 仰臥位メッシュの生成（S720）
本処理ステップでは、仰臥位メッシュ生成部1180は、事例データ（S700）と、伏臥位メッシュ（S710）に基づいて、仰臥位メッシュを自動生成する。この仰臥位メッシュを生成する方法の一例について図17を用いて説明する。
【０１３５】
(3-1) 対応点の取得（S7200）
本処理ステップでは、仰臥位メッシュ生成部1180は、事例データ（S700）に含まれる伏臥位画像と仰臥位画像から複数の対応点の抽出を行い、抽出された複数の対応点の間の対応付けを行う。この方法は、例えば非特許文献3に記載のn-SIFT法など周知の技術により実行できる。
【０１３６】
(3-2) 仮の変形場の算出（S7210）
本処理ステップでは、仰臥位メッシュ生成部1180は、複数の対応点（S7200）の位置の情報に基づき、伏臥位画像から仰臥位画像への仮の変形場を生成する。この処理は、第1実施形態のS3110と同様の処理により実行する。
【０１３７】
(3-3) 仮の仰臥位メッシュの生成（S7220）
本処理ステップでは、仰臥位メッシュ生成部1180は、仮の変形場（S7220）に基づいて、伏臥位メッシュ（S710）を構成するノードの位置を変位させた仮の仰臥位メッシュのノードを生成する。この処理は、第1実施形態のS330の処理と同様の処理により実行する。
【０１３８】
(3-4) 仰臥位メッシュ群の最適化（S7230）
本処理ステップでは、仰臥位メッシュ生成部1180は、事例データ（S700）と、仮の仰臥位メッシュ（S7220）とに基づいて、仮の仰臥位メッシュをより適切なものに最適化する処理を実行する。
【０１３９】
この処理は、例えば以下のようにして実行される。まず、仰臥位メッシュ生成部1180は、仮の仰臥位メッシュを構成する各ノードの位置を微小に変化させた複数のメッシュを生成する。ここでは、これらのメッシュを仰臥位メッシュ候補群と称する。そして、この仰臥位メッシュ候補群の夫々について評価値を算出する。評価値の算出方法は後に詳述する。そして、複数の仰臥位メッシュ候補群の夫々について算出した評価値に基づき、その中で最も評価値が高いメッシュを選択する。そして、選択したメッシュを最適化された仰臥位メッシュとして記録する。
【０１４０】
また、本処理ステップでは、仰臥位メッシュ生成部1180は、上記の処理を繰り返して実行するようにしても良い。この時、2回目以降の処理では仮の仰臥位メッシュの代わりに、1回前の処理の結果として記録した仰臥位メッシュを用いるようにできる。すなわち、本処理ステップは、仮の仰臥位メッシュを、評価値が高くなるように逐次更新していくことができる。これにより、より精度良く仰臥位メッシュを生成することができる。
【０１４１】
ここで、評価値の算出方法の一例について説明する。まず、伏臥位メッシュ（S710）と仰臥位メッシュ候補群（S7230）の夫々のメッシュに基づいて、変形場の候補群を算出する。この処理は、S710で行ったメッシュ生成の処理において記録したノードの接続に関する情報に基づいて実行する。すなわち、メッシュ要素毎に、その要素を構成する8個のノードの位置の情報を用いて、メッシュ要素の内部の変位を線形補間などにより求める。この方法は、メッシュモデルを用いた変形表現に関する周知の方法により実行できる。この結果、仰臥位メッシュの候補群の夫々について、変形場が算出される。ここでは、これらの変形場を、変形場の候補群と称する。
【０１４２】
次に、変形場の候補群の夫々について以下の処理を実行する。すなわち、事例データ（S700）の伏臥位画像を、前述の変形場の候補群の夫々に基づいて変形させた、変形画像の候補群を生成する。この処理は、S560と同様の処理により実行する。
【０１４３】
そして、事例データの仰臥位画像（S700）と、生成した変形画像の候補群との間の画像類似度を算出する。画像類似度を算出する方法は、例えば両画像の画素値の差分の累積値に基づいて算出する方法や、両画像の画素値の相互情報量を用いて算出する方法など、周知の方法により実行できる。ここで、画像類似度は、スカラー値であり、その値が大きいほど両画像が類似していることを示す値である。以上の処理により算出した画像類似度を仰臥位メッシュの候補群の夫々に対する評価値とする。
【０１４４】
(4) 事例ベクトル群生成（S730）
本処理ステップでは、統計モデル生成部1030は、第1実施形態のS350と同様の処理を実行する。ここでは説明を省略する。
【０１４５】
以上に説明したS700〜S740の処理により、本実施形態における学習フェーズの処理が実行される。
【０１４６】
本実施形態による処理装置によれば、学習フェーズの処理において、伏臥位メッシュに対応する仰臥位メッシュの生成が、事例データの画像の画素値等に基づいて半自動的に実行される。したがって、ユーザによる入力操作の煩雑さを低減させることができ、ユーザにとってより簡便に対象症例の変形推定を行うことができる。
【０１４７】
（第4実施形態）
[第4実施形態の概要]
本実施形態による処理装置は、第1実施形態による処理装置が学習フェーズの処理として生成する統計モデルを多段階モデルとして生成する。そして、学習フェーズにおいてこの多段階の統計モデルを用いて対象症例の変形推定を行う。これにより、対象症例の変形推定の精度や計算効率をより高くできるようにする。
【０１４８】
[機能構成]
本実施形態による処理システムの構成は、図1に示した第1実施形態の構成と同一であるため、説明を省略する。
【０１４９】
[処理フロー]
本実施形態による処理装置が行う全体の動作は、第1実施形態の処理装置が行う全体の動作を説明する図3および図4を用いて説明できる。
【０１５０】
まず、本実施形態による処理装置が実行する学習フェーズの処理について図3を用いて説明する。S300〜S330の処理は、第1実施形態の同処理ステップと同一であるので、ここでは説明を省略する。
【０１５１】
(1) 事例ベクトル群生成（S340）
本処理ステップでは、統計モデル生成部1030は、伏臥位のメッシュ情報および仰臥位のメッシュの情報に基づいて事例ベクトルを生成する。この処理は、第1実施形態のS320の処理と同様である。ただし、本実施形態では、事例ベクトルとして、各事例の乳房の伏臥位および仰臥位のメッシュのノードの位置座標値を並べたベクトルを夫々生成する点で、第1実施形態とは異なる。本実施形態における伏臥位および仰臥位のメッシュのノードの数は夫々M個であり、夫々の位置座標は3次元ベクトルであるから、事例ベクトルとして、3×M次元の伏臥位の事例ベクトルxprone,iと、3×M次元の仰臥位の事例ベクトルxsupine,iとを夫々生成する。ただし、添え字のiは事例データの乳房のインデックスであり、N個の事例データを扱う本実施形態において、その範囲は1≦i≦Nである。
【０１５２】
(2) 統計モデル導出（S350）
本処理ステップでは、統計モデル生成部1030は、伏臥位の事例ベクトル群xprone,iおよび仰臥位の事例ベクトル群xsupine,i（S340）を統計的に解析することより、統計モデルを導出する。統計モデルを導出する具体的な一例について説明する。
【０１５３】
まず、伏臥位の事例ベクトル群xprone,iについて主成分分析を行い、伏臥位の平均ベクトルxpaverageと伏臥位の主成分ベクトルepd(1≦d≦Dp)を算出する。ここで、Dpは主成分分析により算出する主成分の数であり、例えば累積寄与率に閾値を設ける方法や、各固定値を使うなどの方法により決定する。さらに、数13の計算を実行することにより、事例データの伏臥位の1次パラメータベクトルcp,iを算出する。
【数１３】
2014018376.tif 000015
ここで、Epは伏臥位の主成分ベクトルを並べた行列、すなわちEp ={e p,1,e p,2,,,e p,Dp}Tであり、これを伏臥位の主成分行列と称する。
【０１５４】
次に、仰臥位の事例ベクトル群xsupine,iについても同様に主成分分析を行い、仰臥位の平均ベクトルxaaverageと仰臥位の主成分ベクトルesd(1≦d≦Ds)を算出する。ここで、Dsは主成分分析による算出する主成分の数であり、例えば累積寄与率に閾値を設ける方法や、固定値を使うなどの方法により決定する。さらに、数14の計算を実行することにより、各事例データの仰臥位の1次パラメータベクトルcs,iを算出する。
【数１４】
2014018376.tif 000016
ここで、Esは伏臥位の主成分ベクトルを並べた行列、すなわちEs ={e s,1,e s,2,,,e s,Ds }Tであり、これを仰臥位の主成分行列と称する。
【０１５５】
次に、伏臥位および仰臥位の1次パラメータを並べた結合1次パラメータベクトルca,iを数15に示すようにして生成する。
【数１５】
2014018376.tif 000017
そして、N個の事例データについて生成したN個の結合1次パラメータベクトルca,iを主成分分析して、結合1次パラメータベクトルの平均ベクトルcaaverageと主成分ベクトルead(1≦d≦Da)を算出する。ここで、このcaaverageを2次平均ベクトル、eadを2次主成分ベクトルと称する。また、Daは主成分分析による算出する主成分の数であり、例えば累積寄与率に閾値を設ける方法や、固定値を使うなどの方法により決定する。
【０１５６】
さらに、N個の事例データの夫々について、数16の計算を実行することにより、各事例データの2次パラメータベクトルca2,iを算出する。
【数１６】
2014018376.tif 000018
ここで、Eaは結合2次主成分ベクトルを並べた行列、すなわちEa ={ea,1,e a,2,,,ea,Da }Tであり、これを2次主成分行列と称する。
【０１５７】
本処理ステップで導出する統計モデルは、2次パラメータベクトルに基づいて伏臥位および仰臥位のメッシュを生成するモデルであり、その具体的な生成過程は数17から数19の順に計算を行うことによる。
【数１７】
2014018376.tif 000019
【数１８】
2014018376.tif 000020
【数１９】
2014018376.tif 000021
つまり、数17の右辺を計算することで、左辺値であるcp,csを算出し、それに基づいて数18および数19の右辺値を計算することで、左辺値であるxprone,xsupineを算出する。
【０１５８】
本処理ステップで導出する統計モデルは、具体的には、数17から数19に示したcaaverage、Ea、xpaverage、Ep、xsaverage、Esである。ここで、xpaverageとEpを第1の統計モデル、xsaverageとEsを第2の統計モデル、caaverageとEaを第3の統計モデルと称すると、数13と数14から第1の統計モデルおよび第2の統計モデルが求められ、これらの統計モデルにおける結合パラメータベクトルを用いて数16から第3の統計モデルが求められることがわかる。本処理ステップではこれらを主メモリ212に記録する。
【０１５９】
次に、本実施形態による処理装置が実行する変形推定フェーズの処理について図4を用いて説明する。なお、S400、S410、S440の処理は、第1実施形態の同処理ステップと同様であるので説明を省略する。
【０１６０】
(1) 最適パラメータの算出（S420）
本処理ステップでは、変形パラメータ算出部1120は、第1実施形態のS420と同様の処理を実行する。ただし、変形パラメータのうち、統計モデルの係数は数17に示したca2である。また、この係数から、伏臥位および仰臥位のメッシュのノード位置を算出する計算は、数17から数19の計算により実行する。すなわち、第1実施形態のS4210の処理の説明では数8により統計モデルの変動の候補から形状を算出するとしていたが、本実施形態では、数17から数19の計算により実行する。
【０１６１】
(2) 変形場の算出（S430）
本処理ステップでは、変形画像生成部1130は、第1実施形態のS430と同様の処理を実行する。ただし、第1実施形態のS430の説明では、S420で算出した統計モデルの係数から数10の計算により、伏臥位および仰臥位のメッシュのノードの位置を算出するとした。しかし、本実施形態では、それに代えて数17から数19の計算により実行する。
【０１６２】
以上に説明した方法により、第4実施形態における処理装置の処理が実施される。本実施形態では学習フェーズの処理で導出する統計モデルを、事例データの伏臥位および仰臥位の夫々に関する個別の統計モデルを導出し、それらのモデルの関係を記述するモデルとして統計モデルを導出している。これにより、伏臥位および仰臥位とを個別の統計モデルとして導出する際に、各モデルに異なる主成分数を設定できる。したがって、事例データの統計的な特性に応じて、より詳細に主成分数を定めることができる。そのため、変形推定の精度の向上や、計算効率の向上が行えるという効果がある。
【０１６３】
（変形例4-1）
本実施形態では、事例データの伏臥位および仰臥位のメッシュのノード位置に基づいて統計モデルを生成する場合を例として説明した。しかし、本発明の実施形態はこれに限らない。例えば、統計モデルは、事例データの伏臥位および仰臥位の画像の輝度値に基づくものであってよい。この場合の一例について説明する。
【０１６４】
例えば、本実施形態のS340で生成する事例ベクトル群は、前述の実施形態で説明した伏臥位の事例ベクトルxprone,iと仰臥位の事例ベクトルxsupine,iに加えて、伏臥位の事例画像ベクトル（第1の輝度ベクトル）xprone_img,iおよび仰臥位の事例画像ベクトル（第2の輝度ベクトル）xsupine_img,iも生成する。ここで伏臥位の事例画像ベクトルとは、例えば、S320で生成した伏臥位メッシュの複数のメッシュ要素の内部の領域における伏臥位画像の平均輝度値とすることができる。メッシュ要素の数が64個となる本実施形態では、xprone_img,iは64次元のベクトルとなる。仰臥位の事例画像ベクトルも同様に、S330で生成した仰臥位メッシュのメッシュ要素の内部の領域における仰臥位画像の平均輝度とすることができ、同様にxsupine_img,iも64次元のベクトルとなる。以上に説明した事例画像ベクトル群を生成して記録する。
【０１６５】
そして、S350では、伏臥位および仰臥位の事例ベクトルに対応する平均ベクトル、主成分行列、1次パラメータベクトルを算出する処理に加え以下の処理を行う。すなわち、前述の伏臥位および仰臥位の事例画像ベクトルについても同様に平均ベクトル、主成分行列、1次パラメータベクトルを算出する。ここでは、伏臥位の事例画像ベクトルに関する平均ベクトルをxpaverage_img、主成分行列をEp_img、1次パラメータベクトルをcp_img,iとする。また仰臥位の事例画像ベクトルに関する平均ベクトルをxsaverage_img、主成分行列をEs_imgとし、1次パラメータベクトルをcs_img,iとする。そして、数15で説明した結合1次パラメータベクトルを数20のように生成する。
【数２０】
2014018376.tif 000022
【０１６６】
そして、以下、第4実施形態のS350の処理と同様にして、1次パラメータベクトルを主成分分析して、2次平均ベクトルcaaverageと2次主成分ベクトルead(1≦d≦Da)を算出する。また、各事例データの2次パラメータベクトルca2,iは第4実施形態で説明した数16により計算することができる。
【０１６７】
また、第4実施形態で数17から数19を用いて説明した統計モデルの係数と、統計モデルが表現する情報との関係は、本変形形態では以下の数21から数25のようになる。また、第4実施形態で説明した統計モデルと対応させて、xpaverageとEpを第1の統計モデル、xsaverageとEsを第2の統計モデル、caaverageとEaを第3の統計モデル、xpaverage_imgとEp_imgを第4の統計モデル、xsaverage_imgとEs_imgを第5の統計モデルと称する。本処理ステップではこれらを主メモリ212に記録する。
【数２１】
2014018376.tif 000023
【数２２】
2014018376.tif 000024
【数２３】
2014018376.tif 000025
【数２４】
2014018376.tif 000026
【数２５】
2014018376.tif 000027
【０１６８】
すなわち、第4実施形態における統計モデルは、その係数に基づいて伏臥位および仰臥位のメッシュのノード位置を表現するが、本変形例の統計モデルはそれに加え、伏臥位および仰臥位の画像の輝度値（メッシュ要素内の輝度値の平均）も表現する。以上に説明した方法により、S350の処理を実行する。
【０１６９】
また、変形推定フェーズの処理の一部であるS420では、統計モデルが表現するメッシュのノード位置に関する情報に加えて、前述の画像の輝度値の情報も用いて統計モデルの係数と剛体変換を算出する。すなわち、S400で取得した対象症例の伏臥位および仰臥位の画像の輝度値と、統計モデルが表現する画像の輝度値との類似性を数9で説明した評価値の計算に加える。これにより、対象症例の伏臥位および仰臥位の体表形状や乳頭位置だけでなく、画像の輝度値にも基づいて変形パラメータの最適化が実行される。
【０１７０】
以上に説明した、本実施形態の変形形態によれば、事例データの乳房の形状だけでなく、その乳房の内部を撮影した画像の輝度値にも基づいて統計モデルを導出する。そのため、例えば乳房内部の乳腺の構造や脂肪の量など、その差異が画像の輝度値として表れ、かつ乳房の変形に関係（例えば相関性）がある特性を考慮した統計モデルが導出される。これにより、対象症例の画像から乳房の形状等の幾何的な特性だけでなく、乳房内部の画像の輝度値にも基づいて変形の推定を行うことができる。よって、より精度の高い変形推定を行える効果がある。
【０１７１】
また、以上の説明では伏臥位および仰臥位のメッシュ要素の内部の領域における画像の平均輝度値を要素値としたベクトルを用いて統計モデルを生成する場合を例として説明したが、本発明の実施形態はこれに限らない。たとえば、メッシュ要素の内部を、さらに細分化した小領域に分割し、その小領域毎に、その内部の領域における画像の平均輝度値を用いるようにできる。この場合、各症例の画像の輝度値を、より詳細に反映した統計モデルが構築できるため、さらに精度の高い変形推定を行える効果がある。
【０１７２】
（変形例4-2）
本実施形態では、事例データの伏臥位および仰臥位のメッシュのノード位置に基づいて統計モデルを生成する場合を例として説明した。しかし、本発明の実施形態はこれに限らない。例えば、統計モデルは、事例データの症例に関する臨床情報にも基づくものであってよい。ここで、変形形態の一例について説明する。
【０１７３】
例えば、本実施形態のS300では事例データである複数の症例の伏臥位および仰臥位の画像に加えて、症例の臨床情報を取得するようにできる。ここで臨床情報とは、例えば年齢、身長、体重、肥満度、BMI、人種、既往歴などである。
【０１７４】
そして、S340で生成する事例ベクトル群は、実施形態で説明した伏臥位の事例ベクトルxprone,iと仰臥位の事例ベクトルxsupine,iに加えて、事例データの症例の事例臨床情報ベクトルxinfo,iも生成する。ここで事例臨床情報ベクトルとは、S300で取得した臨床情報を数値化したベクトルである。臨床情報として例えば年齢、身長、体重、肥満度、BMIなどは、その数値の情報を並べたベクトルとすることができる。ここで、これらの各項目は予め事例データ内での分布等を調べて適切に正規化等の処理を施すのが望ましい。また、人種や既往歴など、元来数値情報でない情報は、数値情報へマッピング（ラベリング）したり、乳房の変形などと線形もしくは単調な関係に近づくように多次元ベクトル化などを施したりするようにできる。
【０１７５】
そして、S350では、本実施形態で説明した伏臥位および仰臥位の事例ベクトルに対応する平均ベクトル、主成分行列、1次パラメータベクトルを算出する処理に加え以下の処理を行う。すなわち、事例臨床情報ベクトルについても同様に平均ベクトル、主成分行列、1次パラメータベクトルを算出する。ここでは、事例臨床情報ベクトルに関する平均ベクトルをxaverage_info、主成分行列をEinfo、1次パラメータベクトルをcinfo,iとする。そして、数15で説明した統合1次パラメータベクトルを数26のように生成する。
【数２６】
2014018376.tif 000028
【０１７６】
そして、以下、第4実施形態のS350の処理と同様にして、1次パラメータベクトルを主成分分析して、2次平均ベクトルcaaverageと2次主成分ベクトルead(1≦d≦Da)を算出する。また、各事例データの2次パラメータベクトルca2,iは第4実施形態で説明した数16により計算することができる。
【０１７７】
また、第4実施形態で数17から数19を用いて説明した統計モデルの係数と、統計モデルが表現する情報との関係は、本変形形態では以下の数27から数30のようになる。
【数２７】
2014018376.tif 000029
【数２８】
2014018376.tif 000030
【数２９】
2014018376.tif 000031
【数３０】
2014018376.tif 000032
【０１７８】
すなわち、第4実施形態における統計モデルは、その係数に基づいて伏臥位および仰臥位のメッシュのノード位置を表現するが、本変形形態の統計モデルはそれに加え、臨床情報も表現する。以上に説明した方法により、S350の処理を実行する。
【０１７９】
また、変形推定フェーズの処理としては、まずS400では、対象症例の伏臥位および仰臥位画像に加えて、対象症例の臨床情報も取得する。そして、S420では、統計モデルが表現するメッシュのノード位置に関する情報に加えて、前述の臨床情報も用いて統計モデルの係数と剛体変換を算出する。すなわち、S400で取得した対象症例の臨床情報と、統計モデルが表現する臨床情報との類似性を数9で説明した評価値の計算に加える。これにより、対象症例の伏臥位および仰臥位の体表形状や乳頭位置だけでなく、臨床情報にも基づいて変形パラメータの最適化が実行される。
【０１８０】
以上に説明した、本実施形態の変形形態によれば、事例データの乳房の形状だけでなく、その事例データの症例に関する臨床情報にも基づいて統計モデルを導出する。そのため、例えば事例データの症例の年齢、身長、体重、肥満度、BMI、人種、既往歴など、乳房の変形に関係（例えば相関性）がある特性を考慮した統計モデルが導出される。これにより、対象症例の乳房の形状等の幾何的な特性だけでなく、臨床情報にも基づいて変形の推定を行うことができる。よって、より精度の高い変形推定を行える効果がある。
【０１８１】
（その他の実施形態）
また、本発明の目的は、前述した実施形態の機能を実現するソフトウェアのプログラムコードを記録した記録媒体（または記憶媒体）を、システムあるいは装置に供給し、そのシステムあるいは装置のコンピュータ（またはＣＰＵやＭＰＵ）が記録媒体に格納されたプログラムコードを読み出し実行することによっても、達成されることは言うまでもない。この場合、記録媒体から読み出されたプログラムコード自体が前述した実施形態の機能を実現することになり、そのプログラムコードを記録した記録媒体は本発明を構成することになる。
【０１８２】
また、コンピュータが読み出したプログラムコードを実行することにより、前述した実施形態の機能が実現されるだけでなく、そのプログラムコードの指示に基づき、コンピュータ上で稼働しているオペレーティングシステム（ＯＳ）などが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。
【０１８３】
さらに、記録媒体から読み出されたプログラムコードが、コンピュータに挿入された機能拡張カードやコンピュータに接続された機能拡張ユニットに備わるメモリに書き込まれた後、そのプログラムコードの指示に基づき、その機能拡張カードや機能拡張ユニットに備わるＣＰＵなどが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。
【０１８４】
本発明を上記記録媒体に適用する場合、その記録媒体には、先に説明したフローチャートに対応するプログラムコードが格納されることになる。
【０１８５】
なお、上述した本実施の形態における記述は、本発明に係る好適な画像処理装置の一例であり、本発明はこれに限定されるものではない。
【図１】
2014018376.tif 000033
【図２】
2014018376.tif 000034
【図３】
2014018376.tif 000035
【図４】
2014018376.tif 000036
【図５】
2014018376.tif 000037
【図６】
2014018376.tif 000038
【図７】
2014018376.tif 000039
【図８】
2014018376.tif 000040
【図９】
2014018376.tif 000041
【図１０】
2014018376.tif 000042
【図１１】
2014018376.tif 000043
【図１２】
2014018376.tif 000044
【図１３】
2014018376.tif 000045
【図１４】
2014018376.tif 000046
【図１５】
2014018376.tif 000047
【図１６】
2014018376.tif 000048
【図１７】
2014018376.tif 000049