JP6626344B2

JP6626344B2 - 画像処理装置、画像処理装置の制御方法およびプログラム

Info

Publication number: JP6626344B2
Application number: JP2015257328A
Authority: JP
Inventors: 中野　雄太; 雄太中野
Original assignee: Canon Inc
Current assignee: Canon Inc
Priority date: 2015-09-29
Filing date: 2015-12-28
Publication date: 2019-12-25
Anticipated expiration: 2035-12-28
Also published as: US20180276799A1; US10672111B2; US20170109871A1; JP2017064370A; US10007973B2

Description

本発明は、画像処理装置、画像処理装置の制御方法およびプログラムに関し、特に、Ｘ線コンピュータ断層撮影装置（Ｘ線ＣＴ）などの医用画像収集装置（モダリティ）で撮像した医用画像を処理する技術に関する。

医用画像（例えば、被検体内部の情報を表す３次元断層画像）を用いた画像診断において、病変部や血管、特定の臓器などの解剖学的部位を強調するために造影剤を用いて撮像を行うことがある。造影剤の投与によって、隣り合う臓器や臓器の中に存在する病変部・血管などのコントラストが高くなり、その形状を把握しやすくなる。しかし、造影剤を投与してからの経過時間が異なると、同じ解剖学的部位であっても染まり方が異なる。そこで、投与からの時間経過を同定するための情報が、画像やカルテなどに記録されている。

例えば、造影ＣＴにおいては、時相と呼ばれる造影剤投与からどれくらい時間が経過してから撮像したかを表すパラメータが定義されており、画像のヘッダ部分に含まれている。時相は、４つのカテゴリに分類して記録されることが多い。まず1つ目は、造影剤を投与しないで撮像された「非造影」と呼ばれるカテゴリである。次に、造影剤を投与した画像については、投与からの経過時間が早い順に、「早期相」、「門脈相」、「晩期相」の３つのカテゴリに分類される。

ところで、臓器などの解剖学的部位の濃度値の分布は、同じ部位であっても時相によって大きく異なる。そのため、例えば各部位の領域を抽出する画像処理において、様々な時相の画像を処理対象とすると、一つの画像処理パラメータでは正確な抽出が行えない場合がある。医用画像処理の代表的な領域抽出手法の一つとして最大事後確率推定法（以下、ＭＡＰ（ＭａｘｉｍｕｍａＰｏｓｔｅｒｉｏｒｉ）法）というものがある。この手法では、予め用意した各部位の存在確率を表す確率アトラスと、各部位の濃度値の事前分布を用いて、画像中の各ボクセルがどの部位に所属するかを示す尤度を算出する。そして、その尤度に基づき、画像から各部位の領域を抽出する。

非特許文献１では、ＭＡＰ法を用いて肝臓や腎臓領域を抽出している。ＭＡＰ法では、各部位の濃度値の事前分布をなるべく精度よく与えることが重要である。

国際公開第２０１１／１２５３６９号米国特許第６９３７２１２号明細書

ＨｙｕｎｊｉｎＰａｒｋ，ＰｅｙｔｏｎＨ. Ｂｌａｎｄ，ＣｈａｒｌｅｓＲ. Ｍｅｙｅｒ"ＣｏｎｓｔｒｕｃｔｉｏｎｏｆａｎＡｂｄｏｍｉｎａｌＰｒｏｂａｂｉｌｉｓｔｉｃＡｔｌａｓａｎｄｉｔｓＡｐｐｌｉｃａｔｉｏｎｉｎＳｅｇｍｅｎｔａｔｉｏｎ"，ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓｏｎＭｅｄｉｃａｌＩｍａｇｉｎｇ，Ｖｏｌ．２２，Ｎｏ．４，Ａｐｒｉｌ２００３．

しかしながら、非特許文献１に記載の技術では、各部位の濃度値の分布は時相によって大きく異なるため、画像データの時相情報が不明な場合には濃度値の事前分布を精度よく与えることができず、領域抽出を安定的に行うことができない。また、画像データにおいて時相情報の欠損または誤りがあると、時相情報を用いた適切な画像処理パラメータの選択（例えば、ＭＡＰ法における各部位の濃度値の事前分布の選択）ができない。すなわち、画像データに時相情報の欠損または誤りがあった場合、後段の画像処理の精度が低下するという課題がある。

本発明は、上記の課題に鑑みてなされたものであり、画像データに時相情報の欠損または誤りがあった場合でも、後段の画像処理の精度を向上させることを目的とする。

上記の目的を達成するために、本発明の一態様による画像処理装置は以下の構成を備える。即ち、
被検体を造影して撮像した画像を取得する取得手段と、
前記被検体において造影される前記被検体の第一の解剖学的部位を表す第一の領域を前記画像から抽出する第一の領域抽出手段と、
前記第一の領域の濃度値に関する特徴量と前記第一の解剖学的部位の複数の時相における濃度値に関する統計データとの比較結果に基づいて、前記画像の時相を推定する推定手段と、
前記推定手段による推定結果に基づいて、前記被検体の第二の解剖学的部位を表す第二の領域を前記画像から抽出する第二の領域抽出手段と
を備えることを特徴とする。

本発明によれば、画像データに時相情報の欠損または誤りがあった場合でも、後段の画像処理の精度を向上させることが可能となる。

第１の実施形態に係る画像処理装置の構成例を示す図。第１の実施形態に係る画像処理装置が実施する全体の処理手順を示すフローチャート。第１の実施形態の位置合わせ処理で用いるランドマークを画像上に示した図。第２の実施形態に係る画像処理装置の構成例を示す図。第２の実施形態に係る画像処理装置が実施する全体の処理手順を示すフローチャート。第３の実施形態に係る画像処理装置の構成例を示す図。第３の実施形態に係る画像処理装置が実施する全体の処理手順を示すフローチャート。第３の実施形態における時相情報取得部の構成例を示す図。第３の実施形態における時相情報取得部が実施する処理手順を示すフローチャート。グラフカット法におけるグラフを説明する模式図。第４の実施形態に係る画像処理装置の構成例を示す図。第４の実施形態に係る画像処理装置が実施する全体の処理手順を示すフローチャート。第４の実施形態における第二の領域抽出部が実施する処理手順を示すフローチャート。

以下、添付図面を参照しながら本発明に係る画像処理装置の実施形態を説明する。なお、以下の実施形態において示す構成は一例に過ぎず、本発明は図示された構成に限定されるものではない。

（第１の実施形態）
＜概要＞
本実施形態に係る画像処理装置は、３次元医用画像から、臓器や病変部など被検体の所定の解剖学的部位の領域の抽出を行う。このとき、入力された画像から時相（例えば、「非造影」、「早期相」、「門脈相」、「晩期相」）を推定し、その推定結果に基づいて後段の処理のパラメータを設定する。これにより、入力画像のヘッダ部分に時相情報が含まれていなくても、もしくは時相情報が誤っていても、適切なパラメータが選ばれ、良好な領域抽出結果を得ることができる。

より具体的には、本実施形態に係る画像処理装置は、時相情報を用いずに行う第一の解剖学的部位の領域抽出（第一の領域抽出）と、推定された時相情報を用いて行う第二の解剖学的部位の領域抽出（第二の領域抽出）とを行う。第一の領域抽出の時点では時相情報が使えないため、抽出対象となる部位に関する濃度値の事前分布を高い精度で得ることができない。そこで、本実施形態の画像処理装置は、統計的に患者間で位置変動の少ない解剖学的部位だけを第一の解剖学的部位として選択する。そして、第一の領域抽出として、位置情報に重みをおいた（濃度値の事前分布が曖昧でも実施可能な）領域抽出を実行する。さらに、この第一の領域抽出の結果を用いて、時相の推定に用いる画像特徴量として、第一の領域内の濃度値の分布情報を算出する。そして、この濃度値の分布情報に基づいて、画像の時相を推定する。

一方、第二の領域抽出は、推定された時相情報を用いて行われる。すなわち、本実施形態の画像処理装置は、推定された時相に基づいて、第二の領域抽出に用いる画像処理のパラメータを選択する。より具体的には、時相情報に基づいて、第二の解剖学的部位の濃度値の事前分布を高い精度で推定し、当該事前分布を用いた領域抽出を実行する。その結果、時相情報を用いずに領域抽出するよりも良好な結果が得られる。

なお、以下の説明において、本実施形態に係る画像処理装置は、被検体の胸腹部の３次元画像を処理対象とする場合を例として説明する。また、解剖学的部位として、肝臓、心臓、左腎臓、右腎臓、および脾臓の各領域を画像から抽出する場合を例として説明する。ここで、肝臓と心臓は、患者間での位置変動が少ない部位であるという性質を有している。すなわち、時相が不明であっても抽出が容易な部位である。また、これらは、濃染のパターンが互いに異なる（心臓は肝臓よりも早くコントラストのピークを迎える）部位であるという性質を有している。すなわち、二つの領域内の濃度値の分布情報を組み合わせることで、時相の推定を安定して行うことができる。以上の性質を考慮して、本実施形態では、肝臓および心臓を第一の解剖学的部位として選択し、他の臓器を第二の解剖学的部位としている。以下、図１乃至図３を用いて、本実施形態の画像処理システムの構成及び処理を説明する。

＜画像処理システムの構成＞
図１は、本実施形態に係る画像処理システムの構成例を示す。同図に示すように、本実施形態における画像処理装置１００は、データサーバ１１０、データベース１２０および表示部１３０と接続されている。

データサーバ１１０は、画像データとして、ある条件（モダリティ、撮像モード、日時、体位等）で被検体を予め撮像して得られた３次元断層画像（ボリュームデータ）を保持している。３次元断層画像を撮像するモダリティは、ＭＲＩ（magnetic resonance imaging）装置、Ｘ線ＣＴ（Computed Tomography）装置、３次元超音波撮影装置、光音響トモグラフィ装置、ＰＥＴ／ＳＰＥＣＴ装置などであってもよい。ここでＰＥＴとはpositron emission tomographyの略であり、ＳＰＥＣＴとはSingle photon emission computed tomographyの略である。画像データは、画像データ取得部１０１０を介して画像処理装置１００に入力される。

データベース１２０は、抽出処理に用いる統計データを保持している。ここでの統計データとは、各部位の確率アトラスや、各部位における濃度値の事前分布（平均値や分散値）のことである。統計データは、統計データ取得部１０２０を介して画像処理装置１００に入力される。

表示部１３０は液晶モニタ等であり、画像処理装置１００が生成する表示画像等の各種の情報を表示する。また、表示部１３０には、ユーザからの指示を取得するためのＧＵＩ（graphical user interface）も配置されている。

＜画像処理装置の構成＞
画像処理装置１００は、以下に説明する構成要素により構成されている。画像データ取得部１０１０は、画像処理装置１００へ入力される画像データ（入力画像）をデータサーバ１１０から取得する。そして、取得した入力画像を、第一の領域抽出部１０４０、特徴量算出部１０５０、第二の領域抽出部１０７０、および表示制御部１０８０へと出力する。

統計データ取得部１０２０は、画像処理装置１００へ入力される統計データをデータベース１２０から取得する。そして、取得した統計データを、第一の領域抽出部１０４０、時相推定部１０６０、第二の領域抽出部１０７０へと出力する。部位選択部１０３０は、第一の解剖学的部位を選択する。そして、選択した部位の情報を、特徴量算出部１０５０および時相推定部１０６０へと出力する。第一の領域抽出部１０４０は、第一の解剖学的部位の領域（第一の領域）を入力画像から抽出する。この処理には、統計データ取得部１０２０が取得した統計データを用いる。そして、抽出した第一の領域の情報を、特徴量算出部１０５０、表示制御部１０８０へと出力する。

特徴量算出部１０５０は、第一の領域抽出部１０４０が抽出した第一の領域の情報を用いて、第一の解剖学的部位の濃度値に関する特徴量を入力画像から算出する。そして、算出した特徴量を時相推定部１０６０へと出力する。時相推定部１０６０は、特徴量算出部１０５０が算出した濃度値に関する特徴量を用いて、時系列での濃度値変化に関する状態、例えば入力画像の造影時の時相（例えば「非造影」、「早期相」、「門脈相」、「晩期相」）を推定する。この処理には、統計データ取得部１０２０が取得した統計データを用いる。そして、推定した時相の情報を第二の領域抽出部１０７０へと出力する。

第二の領域抽出部１０７０は、時相推定部１０６０が推定した時相の情報を用いて、第二の解剖学的部位の領域（第二の領域）を入力画像から抽出する。この処理には、統計データ取得部１０２０が取得した統計データを用いる。そして、抽出した第二の領域の情報を表示制御部１０８０へと出力する。表示制御部１０８０は、入力画像、および領域抽出結果を、ユーザによる不図示の操作部からの入力に応じて表示部１３０に表示させる制御を行う。

ここで、データベース１２０が保持する確率アトラスについて説明する。確率アトラスは、標準化された被検体の画像内において、夫々の画素の位置が何れの部位であるかの確率を示すものである。確率アトラスは、以下の手順で作成される。まず、事前に撮像した多数の症例の夫々の画像データに対して、解剖学的部位を表すラベル画像を手動で作成する。本実施形態では、解剖学的部位のラベルとして、心臓、肝臓、右腎臓、左腎臓、脾臓の５つのラベルを用いる。次に、症例間で画像空間の標準化を行い、ラベル画像を重ね合わせる。この処理を行うと、各画素において、重ね合わされたラベルの数から各部位の存在確率が算出できる。算出された存在確率を画像値として持たせた画像が確率アトラスとなる。なお、空間の標準化は、後述するステップＳ２０１０およびＳ２０２０と同様な処理によって行うことができる。これによって、アトラス作成の際の空間標準化だけでなく、入力画像と確率アトラスの空間標準化も同じアルゴリズムで行うことできる。

また、データベース１２０は、各部位における濃度値の事前分布を、時相ごとに格納している。また、時相が不明な場合に用いる濃度値の事前分布を、各部位ごとに格納している。本実施形態では、濃度値の事前分布としてガウス分布を仮定し、平均値と分散値を保持する場合について説明する。以下の説明では、時相tにおける解剖学的部位lの濃度値の事前分布を、平均値I_{l_t_DB}と分散値σ_{l_t_DB} ²で表す。また、時相が不明な場合に用いる解剖学的部位lの濃度値の事前分布を、平均値I_{l_DB}と分散値σ_{l_DB} ²で表す。ここで、濃度値の事前分布は、確率アトラスの作成に用いた複数症例のラベル画像から算出する。すなわち、夫々の部位について、同じ時相の全て症例における当該部位の領域内の濃度値の平均値と分散値を求め、これを当該時相における当該部位の濃度値の事前分布とする。また、夫々の部位について、夫々の時相の事前分布の平均（平均値の平均と分散値の平均）を求め、これを、時相が不明な場合の当該部位の濃度値の事前分布とする。なお、時相が不明な場合の事前分布は、何れかの時相の事前分布で代用してもよい。この場合、何れの時相の事前分布を用いるかは、データベース１２０に含まれる複数症例の画像データを用いた予備実験で定めることができる。本実施形態では、時相が不明な場合には門脈相の事前分布を適用する。

なお、時相が不明な場合の事前分布は、他の方法で求めてもよい。例えば、夫々の部位について、時相を限定しない全ての症例における当該部位の領域内の濃度値の平均値と分散値を求め、これを、時相が不明な場合の当該部位の濃度値の事前分布と定義してもよい。この場合、全ての時相での濃度値の統計を求めることになるので、時相ごとの統計よりも分散が大きくなる。その結果として、時相が不明な場合には、濃度値の事前分布よりも確率アトラスによる位置情報に重みをおいた推定が実施可能となる。

＜画像処理装置の処理＞
図２は、本実施形態に係る画像処理装置１００が実施する全体の処理手順を示すフローチャートである。

（Ｓ２０００：データの取得）
ステップＳ２０００において、画像データ取得部１０１０は、データサーバ１１０から処理対象である入力画像を取得する。また、統計データ取得部１０２０は、データベース１２０から統計データを取得する。

（Ｓ２０１０：ランドマークの配置）
ステップＳ２０１０において、第一の領域抽出部１０４０は、ステップＳ２０００で取得した入力画像上にランドマークを配置する。このランドマークは後段の確率アトラスと入力画像を位置合わせする際に用いられる。本実施形態では、入力画像中から肺領域と骨盤領域をそれぞれ閾値処理により抽出し、その領域の重心位置に基づいてランドマークを配置する。

具体的には、まず、抽出した肺領域の重心を含むａｘｉａｌ断面と、骨盤領域の重心を含むａｘｉａｌ断面を見つける。次に、そのａｘｉａｌ断面に対して閾値処理により体表を抽出する。最後に、肺と骨盤それぞれの重心から、ａｘｉａｌ断面上で上下左右に探索点を動かし、体表にぶつかったところをランドマークとする。

図３にランドマークを配置した入力画像の模式図を示す。Ｃ３１はｃｏｒｏｎａｌ断面、Ａ３１は肺領域の重心を通るａｘｉａｌ断面、Ａ３２は骨盤領域の重心を通るａｘｉａｌ断面である。Ｃ３１上の波線の矢印がＡ３１とＡ３２の位置を表している。Ａ３１とＡ３２に配置されている白い点が本実施形態に係るランドマークに該当する。なお、ランドマークは必ずしも自動配置される必要はなく、ユーザがユーザインタフェースを用いて入力してもよい。

（Ｓ２０２０：確率アトラスとの位置合わせ）
ステップＳ２０２０において、第一の領域抽出部１０４０は、ステップＳ２０１０で配置されたランドマークを用いて、データベース１２０に格納されている確率アトラスと入力画像との位置合わせを行う。本実施形態では、ステップＳ２０１０で配置されたランドマークと、確率アトラス上のランドマークとを対応させる。そして代表的な変形モデルの一つであるＲａｄｉａｌＢａｓｉｓＦｕｎｃｔｉｏｎ（ＲＢＦ）で変形場を作成し、確率アトラスを入力画像に合うように変形させる。これにより、患者による各解剖学的部位の位置のバイアスを取り除くことができる。なお、画像変形に用いる手法は必ずしもＲＢＦである必要はなく、代表的な変形モデルの一つであるＦｒｅｅＦｏｒｍＤｅｆｏｒｍａｔｉｏｎ（ＦＦＤ）などを用いてもよい。

（Ｓ２０３０：第一の解剖学的部位の選択）
ステップＳ２０３０において、部位選択部１０３０は、画像処理装置１００が抽出対象としている解剖学的部位の中から第一の解剖学的部位を選択する。本実施形態では、先に述べた理由により、心臓および肝臓を第一の解剖学的部位として選択する。なお、第一の解剖学的部位の選択は、入力画像に応じて適応的に行うようにしてもよい。例えば、ステップＳ２０２０で得た確率アトラスとの位置合わせ結果に基づき、候補部位である心臓と肝臓の夫々が入力画像に含まれているか否かを判定し、入力画像に含まれている部位のみを第一の解剖学的部位として選択してもよい。

（Ｓ２０４０：第一の領域抽出）
ステップＳ２０４０において、第一の領域抽出部１０４０は、ステップＳ２０３０で選択された第一の解剖学的部位の領域を入力画像から抽出する。本実施形態では、領域抽出手法として、非特許文献１に記載されているＭＡＰ法を採用する。ＭＡＰ法は、各画素において観測される特徴量v（ここでは濃度値）に基づいて、事後確率が最大になる部位ラベルを夫々の画素に割り当てる手法である。数式で表すと以下のようになる。

ここで、Lは、抽出対象である解剖学的部位に割り当てられたラベルの集合、lはその中の何れか一つのラベルである。また、vは濃度値、p_p(l|v)は位置pにおけるlの事後確率である。第一の領域抽出部１０４０は、夫々の画素について、夫々のlの事後確率p_p (l|v)を求め、式（１）によって夫々の画素のラベルを決定する。なお、第一の領域抽出部１０４０が本ステップで行う処理の目的は第一の解剖学的部位の領域抽出であるが、実際の処理としては、第二の解剖学的部位も含む全ての部位の領域抽出が実行される。

次に、上述した式（１）の演算に必要な、夫々のlの事後確率p_p (l|v)の導出処理について説明する。周知のように、事後確率p_p (l|v)は、ベイズの定理を用いて以下のように書き換えられる。

すなわち、事後確率p_p(l|v)は、位置pにおいて部位ラベルがlであるときの特徴量vの尤度p_p(v|l)と、位置pにおける部位ラベルlの事前確率p_p(l)から算出できる。そこで、第一の領域抽出部１０４０は、夫々の画素について、夫々のlのp_p(l)とp_p(v|l)を求め、式（２）によってp_p(l|v)を算出する。ここで、夫々の画素におけるlの事前確率p_p(l)は、ステップＳ２０００で統計データの一部として取得した確率アトラスから取得できる。

最後に、上述した式（２）の演算に必要な、夫々のlにおける尤度p_p(v|l)の導出処理について説明する。ステップＳ２０２０で説明したとおり、入力画像と確率アトラスとは位置合わせされている。これを用いて、第一の領域抽出部１０４０は、部位ラベルがlであるときの特徴量vの尤度p_p(v|l)、すなわち、各部位の濃度値の分布を推定する。すなわち、次式を用いてp_p(v|l)を算出する。

ここで、pは画像中の画素の位置を表し、Ipは位置pにおける入力画像の濃度値を表す。また、Lはラベルの集合を表す。また、パラメータI_lとσ_l ²は、部位lの領域内における濃度値の分布情報（平均値と分散値）である。この分布情報は、統計データから得た事前分布を初期値として、期待値最大化（ＥＭ：ＥｘｐｅｃｔａｔｉｏｎＭａｘｉｍｉｚａｔｉｏｎ）アルゴリズム（以後、ＥＭアルゴリズムとする）によって反復的に推定される。このとき、分布情報の初期値が入力画像における実際の分布情報に近いほど推定精度が改善される。本ステップの時点では時相情報は使えないため、時相が不明な場合の部位lの事前分布（すなわち、平均値I_{l_DB}と分散値σ_{l_DB} ²）を初期値として用いる。しかし、位置情報である確率アトラスを考慮するＭＡＰ法を採用しているため、患者間の位置変動が少ない心臓と肝臓であれば、分布情報が曖昧な条件下であっても良好な抽出を行うことができる。

以上の処理により、第一の解剖学的部位の領域が抽出される。なお、領域抽出に用いる手法は必ずしもＭＡＰ法である必要はなく、グラフカット法やＬｅｖｅｌｓｅｔ法などを用いてもよい。その際はグラフカット法における各ノードの重みや、Ｌｅｖｅｌｓｅｔ法におけるｚｅｒｏｌｅｖｅｌｓｅｔ（ｆｒｏｎｔ）の設定に、確率アトラスなどの位置情報を付与することで、本実施形態と同等の結果が期待できる。

（Ｓ２０５０：特徴量算出）
ステップＳ２０５０において、特徴量算出部１０５０は、ステップＳ２０４０で抽出された第一の解剖学的部位の領域から、当該領域の濃度値に関する特徴量を算出する。本実施形態では、第一の解剖学的部位の夫々について、領域内の平均濃度値を算出する。なお、算出する特徴量は必ずしも平均濃度値である必要はなく、領域内における濃度値の分散や最大値・最小値などを、当該領域の濃度値に関する特徴量として算出してもよい。また、これらの組み合わせを特徴量としてもよい。

（Ｓ２０６０：時相推定）
ステップＳ２０６０において、時相推定部１０６０は、ステップＳ２０５０で算出した濃度値に関する特徴量を用いて、入力画像の時相を推定する。ここでは、ステップＳ２０００で統計データとして取得した、第一の解剖学的部位の各時相における平均濃度値を用いる。すなわち、入力画像の特徴量と統計データとを比較することで、時相の推定を行う。

具体的には、次式を用いて各時相の評価値D_tを算出する。

ここで、tは時相を表す。L'は第一の解剖学的部位の集合であり、lはその何れか一つ（本実施形態では肝臓か心臓）を表す。また、I_{l_t_DB}は、統計データとして取得した、時相tにおける部位ｌの濃度値の事前分布における平均値を表す。また、I_{l_in}は、入力画像における部位lの平均濃度値であり、ステップＳ２０５０で得た特徴量を表す。すなわち、評価値D_tは、第一の解剖学的部位の平均濃度値の、時相をtと仮定した際の事前分布との差異であり、これをL1ノルムとして求めている。時相推定部１０６０は、夫々の時相に対して評価値D_tの値を算出する。そして、最小のD_tを与える時相tを時相の推定結果t*とする。

なお、評価値D_tの定義は式（４）に限定されるものではなく、例えば、事前分布と平均濃度値のL2ノルム（二乗和の平方根）を用いてもよい。また、事前分布における分散値を考慮した事前分布と平均濃度値のマハラノビス距離をD_tとして用いてもよい。なお、ステップＳ２０５０で平均濃度値以外の特徴量を求めている場合には、統計データとしても対応する値を予め時相ごとに用意し、その差異を評価すればよい。また、複数の特徴量を併用する場合には、それらの評価値の和を用いることができる。

なお、時相の推定方法は上のような算出方法である必要はない。様々な時相を含む多数の症例の第一の解剖学的部位の濃度値特性を事前に学習しておいて、ステップＳ２０５０で得た特徴量（I_{l_in}）を入力として時相を推定する任意の識別器を用いるようにしてもよい。例えば、ＮｅｕｒａｌＮｅｔｗｏｒｋやＳｕｐｐｏｒｔＶｅｃｔｏｒＭａｃｈｉｎｅのような識別器に特徴量を入力して時相を識別してもよい。

なお、本ステップで推定した時相情報は、データサーバ１１０が保持する入力画像の付帯情報として保存するようにしてもよい。特に、入力画像のヘッダ部分に時相情報を保持する項目がある場合には、当該項目に時相情報を記録するようにしてもよい。

（Ｓ２０７０：第二の領域抽出）
ステップＳ２０７０において、第二の領域抽出部１０７０は、ステップＳ２０６０で推定した時相情報を用いて、第二の解剖学的部位の領域を入力画像から抽出する。まず、ステップＳ２０６０で推定した時相情報t*に基づいて、適切な統計データとして、各部位の濃度値の事前分布（すなわち、平均値I_{l_t*_DB}と分散値σ_{l_t*_DB} ²）を選択する。そして、選択した事前分布に基づき、ステップＳ２０４０と同様の処理により第二の領域を抽出する。なお、本ステップの実際の処理では、ステップＳ２０４０と同様に、第一の解剖学的部位と同一の部位も含む全ての部位の領域抽出を行う。なお、第一の解剖学的部位の領域については、本ステップの結果を用いてもよいし、ステップＳ２０４０で得た結果に差し替えて用いてもよい。最後に、第二の領域抽出部１０７０は、全領域の抽出結果を、３次元のラベル画像の形態でデータサーバ１１０へ出力して保存する。

（Ｓ２０８０：領域抽出結果の表示）
ステップＳ２０８０において、表示制御部１０８０は、ステップＳ２０７０で領域抽出された各部位の領域を入力画像上に表示する。本実施形態では、不図示のＵＩ（user interface）を介してユーザが指定した入力画像の断面画像に、ステップＳ２０７０で得たラベル画像の当該断面を重畳して表示する。

以上説明したように、本実施形態によれば、時相情報が欠損、または誤って記載されている画像に対しても、画像から時相を推定することで、高精度な領域抽出ができる。

（変形例１−１）
上記の実施形態では、時相が「非造影」、「早期相」、「門脈相」、「晩期相」の何れの状態であるかを推定していた。しかし、事前分布の選択に十分な分類であれば、画像の時系列での濃度値変化に関する他の状態の分類でもよい。例えば、濃度値の分布が類似する複数の時相を統合して用いてもよい。例えば、「早期相」と「晩期相」を一つにまとめて、「非造影」、「門脈相」、「早期相または晩期相」の何れであるかを推定するようにしてもよい。また、推定する状態は医学の分野で分類される「時相」とは必ずしも一致していなくてもよく、事前分布のバリエーションに基づいて適宜定義してもよい。

また、推定するのは「時相」のような離散的な状態ではなく、造影剤投与からの経過時間のような連続的な状態であってもよい。この場合、データベース１２０には、統計データとして、夫々の解剖学的部位における濃度値の事前分布を、造影剤投与からの経過時間の関数（事前分布関数）の形式で保持しておく。また、第一の解剖学的部位における濃度値の特徴量を引数として、造影剤投与からの経過時間を推定する関数（経過時間関数）を保持しておく。

ここで、事前分布関数は、事前に撮像された複数症例の画像データの夫々の解剖学的部位に関して、夫々の症例における経過時間と濃度値の分布情報の組を得て、時間を説明変数、濃度値の平均値と分散値の夫々を目的変数とした関数を当てはめることで生成する。また、経過時間関数は、事前に撮像された複数症例の画像データの第一の解剖学的部位に関して、夫々の症例における経過時間と夫々の（肝臓と心臓の）濃度値の特徴量の組を得て、夫々の特徴量を説明変数、時間を目的変数とした関数を当てはめることで生成する。

そして、ステップＳ２０６０において、時相推定部１０６０は、時相推定の代わりに以下の処理を行う。すなわち、上記の経過時間関数を用いて、第一の領域抽出結果から算出される濃度値の特徴量から経過時間を推定する。また、ステップＳ２０７０において、第二の領域抽出部１０７０は、ステップＳ２０６０得た経過時間に基づき、上記の（夫々の第二の解剖学的部位の）事前分布関数を用いて、第二の解剖学的部位の夫々の濃度値の事前分布を推定する。そして、該事前分布を初期値として用いて、第二の領域の抽出処理を行うようにすることができる。

（変形例１−２）
上記の実施形態におけるステップＳ２０７０の処理では、時相ごとに定める濃度値の事前分布として、平均値と分散値を得ていた。しかし、平均値と分散値の両方を必ずしも時相に応じて求めなくてもよい。例えば、分散値には共通の固定値を用いることにして、平均値だけを時相に応じて得るようにしてもよい。

（第２の実施形態）
第１の実施形態では、第一の領域抽出の結果から推定した時相の情報を利用して、第二の領域抽出を行っていた。一方、本実施形態における画像処理装置は、時相の推定を介さずに同様の効果を得る。以下、本実施形態に係る画像処理装置について、第１の実施形態との相違部分を中心に説明する。

＜画像処理装置の構成＞
図４は、本実施形態に係る画像処理システムの構成例を示す。第１の実施形態と同じ部分は同じ番号で示している。すなわち、本実施形態に係る画像処理装置４００は、第１の実施形態における時相推定部１０６０の代わりに事前分布推定部４０６０を有する。また、第二の領域抽出部４０７０の処理が、第１の実施形態における第二の領域抽出部１０７０の処理と異なっている。その他の構成については第１の実施形態と同様であるため説明を省略する。

事前分布推定部４０６０は、特徴量算出部１０５０が算出した第一の解剖学的部位の濃度値に関する特徴量を用いて、第二の解剖学的部位の濃度値の事前分布を推定する。この処理には、統計データ取得部１０２０が取得した統計データを用いる。そして、推定した事前分布の情報を第二の領域抽出部４０７０へと出力する。

本実施形態における統計データには、第二の解剖学的部位の夫々についての、当該部位の濃度値の事前分布（平均値と分散値）と、第一の解剖学的部位の濃度値に関する特徴量とを対応付けた関数が含まれている。この関数は、第一の解剖学的部位の濃度値に関する特徴量を引数として与えると、第二の解剖学的部位の濃度値の事前分布を返す対応表のような役割を果たす。この関数は、確率アトラスと同様に、事前に撮像された多数の症例の夫々の画像データと手動で作成された解剖学的部位を表すラベル画像とを用いて作成される。

本実施形態では、まず、画像データとラベル画像を用いて各解剖学的部位の濃度値の分布情報（平均値と分散値）を算出する。次に、第一の解剖学的部位である肝臓と心臓の濃度値の特徴量をそれぞれ軸に割り当てた特徴量空間を定義する。最後に、第二の解剖学的部位の濃度値の分布情報を目的変数、第一の解剖学的部位の濃度値の特徴量を説明変数とした関数を当てはめることで作成する。

第二の領域抽出部４０７０は、事前分布推定部４０６０が推定した第二の解剖学的部位の濃度値の事前分布を用いて、第二の解剖学的部位の領域を入力画像から抽出する。そして、抽出した第二の領域の情報を表示制御部１０８０へと出力する。

＜画像処理装置の処理＞
図５は、本実施形態に係る画像処理装置４００が実施する全体の処理手順を示すフローチャートである。なお、ステップＳ２０００〜Ｓ２０５０およびＳ２０８０で各部が行う処理は、第１の実施形態と同様である。

（Ｓ５０６０：事前分布の推定）
ステップＳ５０６０において、事前分布推定部４０６０は、ステップＳ２０５０で算出した第一の解剖学的部位の濃度値に関する特徴量を用いて、第二の解剖学的部位の濃度値の事前分布を推定する。具体的には、第二の解剖学的部位の事前分布を表す濃度値関数の夫々に第一の解剖学的部位の濃度値に関する特徴量を入力して、事前分布の値を得る。

（Ｓ５０７０：第二の領域抽出）
ステップＳ５０７０において、第二の領域抽出部４０７０は、ステップＳ５０６０で推定した第二の解剖学的部位の濃度値の事前分布を用いて、第二の解剖学的部位の領域を入力画像から抽出する。本ステップの処理は、ステップＳ２０４０と同様である。ただし、ステップＳ５０６０で推定した事前分布を用いる点だけが、ステップＳ２０４０の処理とは異なっている。

以上説明したように、本実施形態によれば、特徴量空間を用いることで、時相よりも細かい状態で分類することが可能となり、より高精度なパラメータ選択を行うことができる。

このように、第１の実施形態、第２の実施形態によれば、画像データに時相情報の欠損または誤りがあった場合でも、画像の時相、または、解剖学的部位の濃度値の事前分布を推定できる。そして、その推定結果を用いることで、後段の画像処理の精度を向上させることができる。

（第３の実施形態）
背景技術において説明した通り、造影剤を投与してからの経過時間が異なると、同じ解剖学的部位であっても染まり方が異なる。そのため、例えば解剖学的部位の領域を抽出する画像処理において、さまざまな時相の画像を処理対象とすると、一つの画像処理パラメータでは正確な抽出が行えない場合がある。特許文献１では、ＤＩＣＯＭ（ＤｉｇｉｔａｌＩｍａｇｉｎｇａｎｄＣｏｍｍｕｎｉｃａｔｉｏｎｉｎＭｅｄｉｃｉｎｅ）規格に従って保存された医用画像のヘッダ部から撮影時刻を取得し、取得した撮影時刻に基づいて領域拡張法や動的輪郭法で使用するパラメータを設定する技術が開示されている。

しかしながら、実臨床で使用する医用画像では、領域拡張法や動的輪郭法といった簡単な領域抽出手法では領域を高精度に抽出することができない。そこで、第３の実施形態、第４の実施形態では、さまざまな時相の医用画像中に存在する注目領域を抽出する領域抽出処理の精度を向上させる技術を提供する。

＜概要＞
本実施形態に係る画像処理装置は、当該装置へ入力された医用画像（入力画像）から臓器や病変部など被検体の所定の解剖学的部位（抽出対象の解剖学的部位）の領域を抽出する。このとき、入力画像から時相の情報（例えば、「非造影」、「早期相」、「門脈相」、「晩期相」）を推定し、その後、入力画像の各画素について、推定した時相の情報に基づき帰属度を算出することを特徴とする。ここで帰属度とは、抽出対象の解剖学的部位のうちのいずれに属しているかを表す値のことである。後段の領域抽出処理では、時相情報に基づいて算出された帰属度を利用して、入力画像から抽出対象の解剖学的部位に対応する領域を抽出する。このような方法により、入力画像がどのような時相であっても、抽出対象の解剖学的部位を良好に抽出することが可能となる。

解剖学的部位に対応する入力画像中の領域は、造影剤投与からの経過時間に応じて特徴的な濃染パターンを有する。本実施形態に係る画像処理装置では、この濃染パターンに基づいて入力画像の時相の情報を推定する。そのために、はじめに解剖学的部位の抽出を行う。しかしながら入力画像から時相情報を取得する段階では、帰属度も（当然のことながら時相情報も）使用できないため、解剖学的部位に対応する領域を高い精度で抽出することができない。そこで、本実施形態の画像処理装置は、統計的に患者間で位置変動の少ない解剖学的部位を、時相推定のための解剖学的部位として選択する。そして、当該部位を、位置情報に重みをおいた（濃度値の事前分布が曖昧でも実施可能な）領域抽出手法により抽出する。このような抽出手法で得られた時相推定のための解剖学的部位の領域から、当該領域内の濃度値の分布情報を時相推定のための画像特徴量として取得する。そして、この濃度値の分布情報に基づいて、画像の時相を推定する。

次に推定した時相の情報に基づいて、各画素の帰属度を算出する。本実施形態における帰属度とは、入力画像のそれぞれの画素（もしくは複数の画素からなる領域）が、抽出対象の解剖学的部位のうちのいずれの部位の一部であるかを表す値であり、各画素（もしくは領域）の画素値を基準となる値と比較することで算出される。このような画素値に基づく帰属度の算出方法では、抽出対象の解剖学的部位のそれぞれについて、比較対象の基準値を入力画像の時相に応じて、適切に選択する必要がある。そこで本実施形態では、入力画像の時相の情報に基づいて、帰属度算出のための基準値を選択する。そして選択された基準値に基づいて、各画素の帰属度を算出する。

抽出対象の解剖学的部位に対応する領域を入力画像から抽出する処理では、算出された帰属度に基づいて領域抽出を実行する。領域抽出に用いる手法は、帰属度を利用する手法であればどのような方法でもよい。時相に基づいて算出された帰属度を利用することで、抽出対象の解剖学的部位に対応する領域を、入力画像の時相の状態によらず良好に抽出することができる。

以下の説明において、本実施形態に係る画像処理装置は、被検体の胸腹部の３次元画像を処理対象とする場合を例として説明する。そして解剖学的部位として、肝臓、心臓、左腎臓、右腎臓、および脾臓の各領域を画像から抽出する場合を例として説明する。また現在注目している解剖学的部位（すなわち肝臓、心臓、左腎臓、右腎臓、脾臓）以外の臓器や軟部組織、空気等の領域を、便宜的に一つの解剖学的部位と見なし、以降、「その他の解剖学的部位」と記述する。本実施形態では、「その他の解剖学的部位」も抽出対象の解剖学的部位の一つに加える。

上述の解剖学的部位のうち、肝臓と心臓は患者間での位置変動が少ない部位であるという性質を有している。すなわち、時相が不明であっても抽出が容易な部位である。また、これら２つの部位は、濃染のパターンが互いに異なる（心臓は肝臓よりも早くコントラストのピークを迎える）部位であるという性質を有している。すなわち、二つの領域内の濃度値の分布情報を組み合わせることで、時相の推定を安定して行うことができる。以上の性質を考慮して、本実施形態では、肝臓および心臓を時相推定のための解剖学的部位として選択する。

本明細書では、以降、抽出対象の解剖学的部位を集合Lで、時相情報を取得する際に参照する解剖学的部位をL_tで表す。また肝臓、心臓、左腎臓、右腎臓、脾臓、「その他の解剖学的部位」をそれぞれ、liver, heart, lkidney, rkidney, spleen, othersというラベルで表す。すなわち、L = {liver, heart, lkidney, rkidney, spleen, others}であり、L_t= {liver, heart }である。なお本明細書では実人体中の解剖学的部位（肝臓、心臓、左腎臓、右腎臓、脾臓）と画像処理装置内で考慮する解剖学的部位(ラベル値でいうとliver, heart, lkidney, rkidney, spleen)を１対１で対応させている。しかしながら、１つの解剖学的部位が同じ時相であっても複数の異なる画像特徴を持つ場合は、該解剖学的部位を複数の仮想的な該解剖学的部位に対応させることができる。例えば、肝臓中に濃度値の明るい部分と暗い部分が存在する場合や、肝臓内に腫瘍などの異常領域が存在する場合は、肝臓を肝臓１と肝臓２(ラベル値はliver1, liver2)という２つの仮想的な解剖学的部位として表す。このような場合はliver1とliver2の両方をLに加えたうえで、以下に説明する処理を実行する。それにより、肝臓領域の明るい部分と暗い部分の双方を高い精度で抽出することが可能となる。他の解剖学的部位についても同様である。

本明細書では、帰属度を利用した領域抽出手法として、特許文献２に記載のグラフカット法による領域抽出手法を説明する。ただし、上述の通り、領域抽出手法はグラフカット法に限定されるものではなく、帰属度を利用する領域抽出手法であれば他の領域抽出手法、例えばＬｅｖｅｌｓｅｔ法やＳｕｐｅｒｐｉｘｅｌ法であってもよい。

以下、図６乃至図１０を用いて、本実施形態の画像処理システムの構成及び処理を説明する。

＜画像処理システムの構成＞
本実施形態に係る画像処理システムの構成例を図６に示す。同図に示すように、本実施形態における画像処理装置６００は、データサーバ６１０および表示部６２０と接続されている。

データサーバ６１０は、画像データとして、ある条件（モダリティ、撮像モード、日時、体位等）で被検体を予め撮像して得られた３次元断層画像（ボリュームデータ）を保持している。３次元断層画像を撮像するモダリティは、ＭＲＩ（magnetic resonance imaging）装置、Ｘ線ＣＴ（Computed Tomography）装置、３次元超音波撮影装置、光音響トモグラフィ装置、ＰＥＴ／ＳＰＥＣＴ装置などであってもよい。ここでＰＥＴとはpositron emission tomographyの略であり、ＳＰＥＣＴとはSingle photon emission computed tomographyの略である。これらの３次元断層画像にはスライス枚数が１枚の画像、すなわち２次元画像も含まれる。画像データは、画像データ取得部６０１０を介して画像処理装置６００に入力される。

またデータサーバ６１０は、抽出処理に用いる統計データを保持している。統計データは、統計データ取得部６０２０を介して画像処理装置６００に入力される。統計データの詳細については後述する。

表示部６２０は液晶モニタ等であり、画像処理装置６００が生成する表示画像等の各種の情報を表示する。また表示部６２０には、ユーザからの指示を取得するためのＧＵＩ（graphical user interface）も配置されている。

＜画像処理装置６００の構成＞
画像処理装置６００は、以下に説明する構成要素により構成されている。画像データ取得部６０１０は、画像処理装置６００へ入力される画像データ（入力画像）をデータサーバ６１０から取得する。そして取得した入力画像を、時相情報取得部６０３０、第一の帰属度算出部６０４０、第一の領域抽出部６０５０、および表示制御部６０６０へと出力する。

統計データ取得部６０２０は、画像処理装置６００へ入力される統計データをデータサーバ６１０から取得する。そして、取得した統計データを、時相情報取得部６０３０、第一の帰属度算出部６０４０、および第一の領域抽出部６０５０へと出力する。

時相情報取得部６０３０は、入力画像について時相の情報を推定する。本処理部は、画像データ取得部６０１０から入力画像を受け取る画像取得処理を行う。また統計データ取得部６０２０から統計データを受け取る。そして取得した統計データに基づいて、入力画像の時相の情報を推定する。この時、入力画像から時相推定のための解剖学的部位の領域を抽出し、その領域の画像情報を参照して、時相の情報を推定する。本処理部で得られた時相の情報は、第一の帰属度算出部６０４０に出力される。

第一の帰属度算出部６０４０は、統計データ取得部６０２０から取得した統計データと、時相情報取得部６０３０で推定した時相の情報を利用して、入力画像中の各画素について帰属度を算出する。ここで帰属度は、抽出対象の部位ごとに計算される。すなわち、入力画像中の１つの画素について、抽出対象の解剖学的部位の数だけ帰属度が計算される。算出した帰属度は、第一の領域抽出部６０５０に出力される。

第一の領域抽出部６０５０は、統計データ取得部６０２０が取得した統計データと、第一の帰属度算出部６０４０で算出した帰属度を参照して、入力画像から抽出対象の解剖学的部位に属する領域を抽出する。すなわち、入力画像中のそれぞれの画素について、その画素が属する部位を決定する。第一の領域抽出部６０５０は、抽出した領域の情報を表示制御部６０６０へと出力する。表示制御部６０６０は、入力画像、および領域抽出結果を、ユーザによる不図示の操作部からの入力に応じて表示部６２０に表示させる制御を行う。

図８は時相情報取得部６０３０の機能構成を表す図である。時相情報取得部６０３０は、時相推定のための領域抽出部８０１０、部位選択部８０２０、特徴量算出部８０３０、時相推定部８０４０を備えている。

時相推定のための領域抽出部８０１０は、時相推定のための解剖学的部位の領域（時相推定のための領域）を入力画像から抽出する。この処理には、統計データ取得部６０２０が取得した統計データを用いる。部位選択部８０２０は、時相推定のための解剖学的部位を選択する。そして選択した部位の情報を、特徴量算出部８０３０および時相推定部８０４０へと出力する。特徴量算出部８０３０は、時相推定のための領域抽出部８０１０が抽出した時相推定のための領域の情報を用いて、時相推定のための解剖学的部位の濃度値に関する特徴量（時相推定のための特徴量）を入力画像から算出する（特徴取得処理）。そして、算出した特徴量を時相推定部８０４０へと出力する。

時相推定部８０４０は、特徴量算出部８０３０が算出した濃度値に関する特徴量を用いて、時系列での濃度値変化に関する状態、例えば入力画像の造影時の時相（例えば「非造影」、「早期相」、「門脈相」、「晩期相」）を推定する。この処理には、統計データ取得部６０２０が取得した統計データを用いる。そして、推定した時相の情報を第一の帰属度算出部６０４０へと出力する。

ここで、データサーバ６１０が保持する統計データについて説明する。統計データには、存在確率アトラス、形状モデル、抽出対象の解剖学的部位の濃度値の事前分布が含まれる。

存在確率アトラスは、標準化された被検体の画像内において、夫々の画素の位置が何れの解剖学的部位であるかの確率（存在確率）を示すものである。存在確率アトラスは、以下の手順で事前に作成される。まず、事前に撮像した多数の症例の夫々の画像データに対して、解剖学的部位の領域を表すマスク画像を手動で作成する。ここでマスク画像とは、画像内の解剖学的部位の領域に含まれる画素とそうでない画素で異なる画素値を取る２値画像のことである。本実施形態では、心臓、肝臓、左腎臓、右腎臓、脾臓、「その他の解剖学的部位」の６つのマスク画像を用いる。ただし「その他の解剖学的部位」については、他の５つのマスク画像より作成可能である。

次に、症例間で画像空間の標準化を行い、マスク画像を重ね合わせる。この処理を行うと、各画素において、重ね合わされたマスクの数から各部位の存在確率が算出できる。算出された存在確率を画像値として持たせた画像が存在確率アトラスとなる。なお、空間の標準化は、後述するステップＳ９０１０およびＳ９０２０と同様の処理によって行うことができる。これによって、アトラス作成の際の空間標準化だけでなく、入力画像と存在確率アトラスの空間標準化も同じアルゴリズムで行うことできる。

形状モデルは、抽出対象の解剖学的部位の形状を表すモデルである。本実施形態で用いる形状モデルは、レベルセット分布モデルを用いて表現される。レベルセット分布モデルを用いた形状表現は、画像処理の分野において公知の手法であるため、本詳細については説明を省略する。存在確率アトラスと同様、形状モデルも事前に作成される。作成手順は以下の通りである。なお、形状モデルは抽出対象のそれぞれの解剖学的部位について１つ作成する。

以下の説明は、抽出対象の１つの解剖学的部位について、対応する形状モデルを作成する手順である。形状モデルの作成には、複数の被検体のマスク画像を用いる。これには存在確率アトラスの作成時に使用したマスク画像と同じもの（空間標準化されたマスク画像）を使用する。はじめに空間標準化された複数個のマスク画像のそれぞれに対して符号付き距離変換を適用して、距離画像を生成する。ここで符号付き距離変換は、マスク画像中の解剖学的部位の領域に対応する画素に領域の境界からの距離値が正の値で、逆に領域外の画素に境界からの距離値が負の値で格納されている画像である。次に生成された個々の距離画像をそれぞれ１つの多次元ベクトルに変換する。

結果として、それぞれのマスク画像に１対１で対応する多次元ベクトルが生成される。こうして得られた多次元ベクトルの集合に対して主成分分析を適用して、１つの平均値ベクトルと複数の固有ベクトルを得る。この平均値ベクトルとそれぞれの固有ベクトルを、レベルセット分布モデルで表現された形状モデルの平均形状（画像）とそれぞれを固有形状（画像）とする。このような手順を抽出対象の解剖学的部位のそれぞれについて実施することで、すべての解剖学的部位の形状モデルを取得できる。

抽出対象の解剖学的部位の濃度値の事前分布は、時相情報取得部６０３０が推定する時相に１対１で対応している。また、時相が不明な場合に用いる濃度値の事前分布を、各部位ごとに格納している。本実施形態では、濃度値の事前分布としてガウス分布を仮定し、平均値と分散値を保持する場合について説明する。

以下の説明では、時相tにおける解剖学的部位lの濃度値の事前分布を、平均値I_{l_t_DB}と分散値σ_{l_t_DB} ²で表す。また、時相が不明な場合に用いる解剖学的部位lの濃度値の事前分布を、平均値I_{l_DB}と分散値σ_{l_DB} ²で表す。ここで、濃度値の事前分布は、存在確率アトラスの作成に用いた複数症例のラベル画像から事前に算出する。すなわち、夫々の部位について、同じ時相の全ての症例における当該部位の領域内の濃度値の平均値と分散値を求め、これを当該時相における当該部位の濃度値の事前分布とする。

また、夫々の部位について、夫々の時相の事前分布の平均（平均値の平均と分散値の平均）を求め、これを、時相が不明な場合の当該部位の濃度値の事前分布とする。なお、時相が不明な場合の事前分布は、何れかの時相の事前分布で代用してもよい。この場合、何れの時相の事前分布を用いるかは、データサーバ６１０に含まれる複数症例の画像データを用いた予備実験で定めることができる。本実施形態では、時相が不明な場合には門脈相の事前分布を適用する。

なお、時相が不明な場合の事前分布は、他の方法で求めてもよい。例えば、夫々の部位について、時相を限定しない全ての症例における当該部位の領域内の濃度値の平均値と分散値を求め、これを、時相が不明な場合の当該部位の濃度値の事前分布と定義してもよい。この場合、全ての時相での濃度値の統計を求めることになるので、時相ごとの統計よりも分散が大きくなる。その結果として、時相が不明な場合には、濃度値の事前分布よりも存在確率アトラスによる位置情報に重みをおいた推定が実施可能となる。

＜画像処理装置６００の処理＞
図７は、本実施形態に係る画像処理装置６００が実施する全体の処理手順を示すフローチャートである。

（Ｓ７０００：データの取得）
ステップＳ７０００において、画像データ取得部６０１０は、データサーバ６１０から処理対象である入力画像を取得する。また、統計データ取得部６０２０は、データサーバ６１０から統計データを取得する。

（Ｓ７０１０：時相情報の取得）
ステップＳ７０１０において、時相情報取得部６０３０は、入力画像の時相の情報を推定する。ここで図９のフローチャートを参照して、時相情報取得部６０３０が行う処理の詳細を説明する。

（Ｓ９０１０：ランドマークの配置）
ステップＳ９０１０において、時相推定のための領域抽出部８０１０は、ステップＳ７０００で取得した入力画像上にランドマークを配置する。このランドマークは後段の存在確率アトラスと入力画像を位置合わせする際に用いられる。本実施形態では、入力画像中から肺領域と骨盤領域をそれぞれ閾値処理により抽出し、その領域の重心位置に基づいてランドマークを配置する。

ランドマークを配置した入力画像の模式図については、第１の実施形態で図３を参照して説明した通りである。

（Ｓ９０２０：存在確率アトラスとの位置合わせ）
ステップＳ９０２０において、時相推定のための領域抽出部８０１０は、ステップＳ９０１０で配置されたランドマークを用いて、データサーバ６１０に格納されている存在確率アトラスと、入力画像との位置合わせを行う。本実施形態では、ステップＳ９０１０で配置されたランドマークと、存在確率アトラス上のランドマークとを対応させる。そして代表的な変形モデルの一つであるＲａｄｉａｌＢａｓｉｓＦｕｎｃｔｉｏｎ（ＲＢＦ）で変形場を作成し、存在確率アトラスを入力画像に合うように変形させる。これにより、患者による各解剖学的部位の位置のバイアスを取り除くことができる。なお、画像変形に用いる手法は必ずしもＲＢＦである必要はなく、代表的な変形モデルの一つであるＦｒｅｅＦｏｒｍＤｅｆｏｒｍａｔｉｏｎ（ＦＦＤ）などを用いてもよい。

（Ｓ９０３０：時相推定のための解剖学的部位の選択）
ステップＳ９０３０において、部位選択部８０２０は、画像処理装置６００が抽出対象としている解剖学的部位の中から時相推定のための解剖学的部位を選択する。本実施形態では、先に述べた理由により、心臓および肝臓を時相推定のための解剖学的部位として選択する。なお、時相推定のための解剖学的部位の選択は、入力画像に応じて適応的に行うようにしてもよい。例えば、ステップＳ９０２０で得た存在確率アトラスとの位置合わせ結果に基づき、候補部位である心臓と肝臓の夫々が入力画像に含まれているか否かを判定し、入力画像に含まれている部位のみを時相推定のための解剖学的部位として選択してもよい。

（Ｓ９０４０：時相推定のための領域の抽出）
ステップＳ９０４０において、時相推定のための領域抽出部８０１０は、ステップＳ９０３０で選択された時相推定のための解剖学的部位の領域を入力画像から抽出する。本実施形態では、領域抽出手法として、非特許文献１に記載されている最大事後確率法（ＭａｘｉｍｕｍａＰｏｓｔｅｒｉｏｒＰｒｏｂａｂｉｌｉｔｙ，ＭＡＰ法）を採用する。ＭＡＰ法は、各画素において観測される特徴量v（ここでは濃度値）に基づいて、事後確率が最大になる部位ラベルを夫々の画素に割り当てる手法である。数式で表すと以下のようになる。

ここで、Lは抽出対象である解剖学的部位に割り当てられたラベルの集合、lはその中の何れか一つのラベルである。また、vは濃度値、p_p(l|v)は位置pにおけるlの事後確率である。時相推定のための領域抽出部８０１０は、夫々の画素について、夫々のlの事後確率p_p (l|v)を求め、式（５）によって夫々の画素のラベルを決定する。なお、時相推定のための領域抽出部８０１０が本ステップで行う処理の目的は時相推定のための解剖学的部位の領域抽出であるが、実際の処理としては、時相推定のための解剖学的部位も含む全ての解剖学的部位の領域抽出が実行される。

次に、上述した式（５）の演算に必要な、夫々のlの事後確率p_p (l|v)の導出処理について説明する。周知のように、事後確率p_p (l|v)は、ベイズの定理を用いて以下のように書き換えられる。

すなわち、事後確率p_p(l|v)は、位置pにおいて部位ラベルがlであるときの特徴量vの尤度p_p(v|l)と、位置pにおける部位ラベルlの事前確率p_p(l)から算出できる。そこで、第一の領域抽出部６０５０は、夫々の画素について、夫々のlのp_p(l)とp_p(v|l)を求め、式（６）によってp_p(l|v)を算出する。ここで、夫々の画素におけるlの事前確率p_p(l)は、ステップＳ７０００で統計データの一部として取得した存在確率アトラスから取得できる。

最後に、上述した式（６）の演算に必要な、夫々のlにおける尤度p_p(v|l)の導出処理について説明する。ステップＳ９０２０で説明したとおり、入力画像と存在確率アトラスとは位置合わせされている。これを用いて、時相推定のための領域抽出部８０１０は、部位ラベルがlであるときの特徴量vの尤度p_p(v|l)、すなわち、各部位の濃度値の分布を推定する。すなわち、次式を用いてp_p(v|l)を算出する。

ここで、pは画像中の画素の位置を表し、I_pは位置pにおける入力画像の濃度値を表す。また、Lはラベルの集合を表す。また、パラメータI_lとσ_l ²は、部位lの領域内における濃度値の分布情報（平均値と分散値）である。この分布情報は、統計データから得た事前分布を初期値として、期待値最大化（ＥＭ：ＥｘｐｅｃｔａｔｉｏｎＭａｘｉｍｉｚａｔｉｏｎ）アルゴリズム（以後、ＥＭアルゴリズムとする）によって反復的に推定される。このとき、分布情報の初期値が入力画像における実際の分布情報に近いほど推定精度が改善される。本ステップの時点では時相情報は使えないため、時相が不明な場合の部位lの事前分布（すなわち、平均値I_{l_DB}と分散値σ_{l_DB} ²）を初期値として用いる。しかし、位置情報である存在確率アトラスを考慮するＭＡＰ法を採用しているため、患者間の位置変動が少ない心臓と肝臓であれば、分布情報が曖昧な条件下であっても良好な抽出を行うことができる。

以上の処理により、時相推定のための解剖学的部位の領域が抽出される。なお、領域抽出に用いる手法は必ずしもＭＡＰ法である必要はなく、グラフカット法やＬｅｖｅｌｓｅｔ法などを用いてもよい。その際はグラフカット法における各ノードの重みや、Ｌｅｖｅｌｓｅｔにおけるｚｅｒｏｌｅｖｅｌｓｅｔ（ｆｒｏｎｔ）の設定に、存在確率アトラスなどの位置情報を付与することで、本実施形態と同等の結果が期待できる。

（Ｓ９０５０：時相推定のための特徴量の算出）
ステップＳ９０５０において、特徴量算出部８０３０は、ステップＳ９０４０で抽出された時相推定のための解剖学的部位の領域から、当該領域の濃度値に関する特徴量（時相推定のための特徴量）を算出する。本実施形態では、時相推定のための解剖学的部位の夫々について、領域内の平均濃度値を算出する。なお、算出する特徴量は必ずしも平均濃度値である必要はなく、領域内における濃度値の分散や最大値・最小値などを、当該領域の濃度値に関する特徴量として算出してもよい。また、これらの組み合わせを特徴量としてもよい。

（Ｓ９０６０：時相推定）
ステップＳ９０６０において、時相推定部８０４０は、ステップＳ９０５０で算出した時相推定のための特徴量を用いて、入力画像の時相を推定する。ここでは、ステップＳ７０００で統計データとして取得した、時相推定のための解剖学的部位の各時相における平均濃度値を用いる。すなわち、入力画像の特徴量と統計データとを比較することで、時相の推定を行う。具体的には、次式を用いて各時相の評価値D_tを算出する。

ここで、tは時相を表す。L_tは時相推定のための解剖学的部位の集合であり、l_tはその何れか一つ（本実施形態では肝臓か心臓）を表す。また、I_{l１_t_DB}は、統計データとして取得した、時相tにおける部位l_tの濃度値の事前分布における平均値を表す。また、I_{lt_in}は、入力画像における部位l_tの平均濃度値であり、ステップＳ９０５０で得た特徴量を表す。すなわち、評価値D_tは、時相推定のための解剖学的部位の平均濃度値の、時相をtと仮定した際の事前分布との差異であり、これをL1ノルムとして求めている。時相推定部８０４０は、夫々の時相に対して評価値D_tの値を算出する。そして、最小のD_tを与える時相tを時相の推定結果t*とする。

なお、評価値D_tの定義は式（８）に限定されるものではなく、例えば、事前分布と平均濃度値のL2ノルム（二乗和の平方根）を用いてもよい。また、事前分布における分散値を考慮した事前分布と平均濃度値のマハラノビス距離をD_tとして用いてもよい。なお、ステップＳ９０５０で平均濃度値以外の特徴量を求めている場合には、統計データとしても対応する値を予め時相ごとに用意し、その差異を評価すればよい。また、複数の特徴量を併用する場合には、それらの評価値の和を用いることができる。

なお、時相の推定方法は上のような算出方法である必要はない。様々な時相を含む多数の症例の時相推定のための解剖学的部位の濃度値特性を事前に学習しておいて、ステップＳ９０５０で得た特徴量（I_{l１_in}）を入力として時相を推定する任意の識別器を用いるようにしてもよい。例えば、ＮｅｕｒａｌＮｅｔｗｏｒｋやＳｕｐｐｏｒｔＶｅｃｔｏｒＭａｃｈｉｎｅのような識別器に特徴量を入力して時相を識別してもよい。

なお、本ステップで推定した時相情報は、データサーバ６１０が保持する入力画像の付帯情報として保存するようにしてもよい。特に、入力画像のヘッダ部分に時相情報を保持する項目がある場合には、当該項目に時相情報を記録するようにしてもよい。

以上でステップＳ７０１０の説明を終える。

（Ｓ７０２０：各臓器への帰属度（第一の帰属度）の算出）
ステップＳ７０２０において、第一の帰属度算出部６０４０は、抽出対象の各々の解剖学的部位について第一の帰属度を算出する。ここで解剖学的部位の集合Lに含まれる一つの部位をlで表す。またLからlを除いた集合をL'_lとし（すなわちL'_l = L − {l}）、L'_lに含まれる個々の解剖学的部位をl'_lとする。すると本ステップでは、入力画像のそれぞれの画素pについて、lの事後確率p_p(l|v)と、L'_lに含まれる解剖学的部位l'_lの事後確率p_p(l'_l|v)の最大値p_p(l'_{l_max}|v)を決定する。そして、それらの事後確率の組(p_p(l|v), p_p(l'_{l_max}|v))を画素pにおける解剖学的部位lの第一の帰属度として算出する。本ステップでは、入力画像のすべての画素について、この第一の帰属度（事後確率の組）を計算する。

まずステップＳ７０１０で推定した時相情報t*に基づいて、解剖学的部位lの濃度値の事前分布（すなわち平均値I_{l_t*_DB}と分散値σ_{l_t*_DB} ²）を統計データ取得部６０２０より取得する。次にステップＳ９０４０で実行したＭＡＰ法を再度実行することで、Lに含まれるすべて解剖学的部位について事後確率p_p (l|v)を計算する。ただし本ステップでのＭＡＰ法では、式（７）の濃度値の分布情報（平均値l_lと分散値σ_l ²）の初期値として時相情報t*に基づいて選択された平均値I_{l_t*_DB}と分散値σ_{l_t*_DB} ²を使用する。本処理により、抽出対象のそれぞれの解剖学的部位について、事後確率p_p(liver|v), p_p(heart|v), ..., p_p(others|v)を得る。

本実施形態では、注目部位の事後確率と、注目部位以外の部位の事後確率の最大値からなる１組の数値を第一の帰属度としたが、第一の帰属度の表現はこれに限定されるものではない。例えば、注目部位の事後確率のみを第一の帰属度としてもよい。また、第一の注目部位の事後確率、第二の注目部位の事後確率、それら以外の部位の事後確率の最大値というように、複数の事後確率を１組の数値として第一の帰属度としてもよい。さらに本実施形態では事後確率の最大値を第一の帰属度の一方の数値としたが、事後確率の平均値を用いてもよい。

（Ｓ７０３０：第一の帰属度を用いて得られる領域に対する形状モデルの当てはめ）
ステップＳ７０３０において、第一の領域抽出部６０５０は、第一の帰属度を用いて領域を抽出した後、抽出した領域に対して形状モデルを当てはめる。本ステップの処理は、抽出対象のそれぞれの解剖学的部位について、それぞれ独立に実行される。以下では、解剖学的部位lについて形状エネルギーを計算する例を説明する。

はじめに第一の領域抽出部６０５０は、データサーバ６１０に格納されている形状モデルを、統計データ取得部６０２０を介して取得する。以下、解剖学的部位lの形状モデルのうち、平均形状画像をφ_l0と、N_shape個存在している固有形状画像のうちk番目の固有形状画像をφ_lk（k = 1, ..., N_shape）と記述する。

次に前ステップで算出した第一の帰属度に基づいて、抽出対象の解剖学的部位lの領域を抽出する。これは入力画像の各画素で、第一の帰属度((p_p(l|v), p_p(l'_max|v))に対して閾値処理p_p(l|v)−p_p(l'_max|v)＞T_lを適用することで得られる。本処理により、抽出対象の解剖学的部位lの領域を表すマスク画像M_ltargetを得る。

次にマスク画像M_ltargetに対して、形状モデルを当てはめる。レベルセット分布モデルで表現された形状モデルでは、形状を表すマスク画像は次式で表現される画像φ_lを値０で閾値処理することで得られる。

ここで式（１０）中の(c_l1, c_l2, ..., c_lNshape)は形状モデルのパラメータである。このパラメータにさまざまな値を設定することで、形状モデルは解剖学的部位lの取りうる様々な形状を表現することができる。マスク画像M_ltargetに対する形状モデルの当てはめでは、次の処理を行う。事前に、パラメータ(c_l1, c_l2, ..., c_lNshape)にさまざまな値を設定することで、複数のパラメータからなるパラメータの集合を定義する。次にパラメータの集合に含まれるそれぞれのパラメータについて式（１０）でφ_lを計算し、さらにφ_lを値０で閾値処理をする。この閾値処理で取得したマスク画像M_{lshape(cl1, cl2, ..., clNshape)}とM_ltargetを比較して、その差を求める。このような処理を事前に定義したすべてのパラメータについて行うことで、差がもっとも小さくなるようなパラメータ(c'_l1, c'_l2, ..., c'_lNshape)を取得する。さらに取得したパラメータを、再度、式（１０）に代入することで、対応する形状モデルの当てはめ結果M'_{lshape(c'l1, c'l2, ..., c'lNshape)}を取得する。

なお、形状モデルのパラメータの集合はさまざまな方法で決定することができる。もっとも簡単な方法は、ある一定範囲の値を網羅的に割り当てることである。具体的にはc_{li_min} ＜ c_{li_max} (ただしi = 1, 2, ..., N_shape)なる値をあらかじめ決め、その２つの数から設定される値の範囲[c_{li_min}, c_{li_max}]を一定数Zで分割して、c_liに割り当てる。すなわちc_li = k * (c_{li_max} - c_{li_min}) / Z + c_{li_min} (ただしk = 0, 1, ..., Z)のように割り当てる。ここでc_{li_min}とc_{li_max}は、形状モデルを構築する際に得られる固有値に基づいて設定することができる。例えば、固有形状画像φ_lkに対応する固有値をλ_lkとする時、c_{li_min} = −Y *λ_lk、c_{li_max} = Y *λ_lkとする。ここでYは、事前に決めた定数である。また固有値λ_lkをそのまま使う代わりに、その逆数1/λ_lkを使用してもよいし、平方根を使用してもよい。

（Ｓ７０４０：第一の帰属度を用いたグラフの生成）
ステップＳ７０４０において、第一の領域抽出部６０５０は、特許文献２で開示されているグラフカット法で使用するグラフを生成する。

ここで図１０を参照して、グラフカット法におけるグラフについて、本実施形態を理解するために必要な事項に限定して説明する。図１０では図面の見やすさを考慮して、一辺が３画素の２次元画像１６００（図１０（ａ））と、それに対応するグラフ１６１０（図１０（ｂ））を用いて、グラフを説明する。

グラフカット法で生成されるグラフ１６１０は、画像１６００の各画素に１対１、もしくは画像中で隣接する複数個の画素に対して１つ対応する画像ノードから構成される。グラフはさらに、領域の前景ラベルと背景ラベルに対応するターミナル・ノード１６１２、１６１３から構成される。隣接する複数の画像ノード間は、エッジで互いに連結される。例えば、画像中のある１つの画素１６０１に対応する画像ノード１６１４と、その画素の周辺４近傍にある画素１６０２に対応する画像ノード１６１５は、エッジ１６１６で互いに連結されている。グラフの連結近傍数は、領域抽出の精度と計算リソースのバランスを考慮して決定されるが、一般的には２次元画像の場合は４または８、３次元画像の場合は６または２６が利用される。画像ノード間を連結するエッジに加えて、グラフ内のすべての画像ノードはそれぞれ、２つのターミナル・ノード１６１２、１６１３とエッジで連結される。例えば図１０（ｂ）の画像ノード１６１５は、ターミナル・ノード１６１２、１６１３とエッジ１６１８、１６１９で連結されている。以下、画像ノード間に張られるエッジをｎ−ｌｉｎｋ、画像ノードとターミナル・ノードの間に張られる２本のエッジをｔ−ｌｉｎｋと記述する。

すべてのｎ−ｌｉｎｋとｔ−ｌｉｎｋは、エネルギーと呼ばれる正の数を持つ。一般にはｎ−ｌｉｎｋには２つの画像ノードに対応する画素（または領域）間の類似性に基づいて、エネルギーを割り当てる。そしてｔ−ｌｉｎｋには、対応する画素の前景らしさと背景らしさの尤度に基づいて、エネルギーを割り当てる。

図７におけるステップＳ７０４０の説明に戻る。本ステップでは、ステップＳ７０２０で算出した第一の帰属度(p_p(l|v), p_p(l'_max|v))とステップＳ７０３０で得られた形状モデルの当てはめ結果に基づいて、グラフカット法で利用するグラフを生成する。グラフは抽出対象のそれぞれの解剖学的部位に対して１つ作成されることに注意されたい。以下では解剖学的部位lに対応するグラフG_lを作成する例を説明する。

まず、入力画像の各画素に対応する画像ノードを作成する。そして画像ノードとターミナル・ノード間をｔ−ｌｉｎｋで、画像ノード間をｎ−ｌｉｎｋで連結する。

次にステップＳ７０２０で算出した第一の帰属度(p_p(l|v), p_p(l'_max|v))に基づいて、グラフG_l中のすべてのｔ−ｌｉｎｋにエネルギー値を設定する。今、画素pに対応する画像ノードをN_pとする。そして画像ノードN_pと２つのターミナル・ノードN_s、N_tを連結するｔ−ｌｉｎｋに割り当てるエネルギーを、それぞれE_sp、E_ptとする。するとE_spとE_ptは、次式のように設定される。

式（１１）にしたがって、グラフG_l中のすべてのｔ−ｌｉｎｋにエネルギーを割り当てる。

最後にグラフ中のすべてのｎ−ｌｉｎｋにエネルギーを設定する。ｎ−ｌｉｎｋに割り当てるエネルギーの計算は、公知の方法にしたがって行う。入力画像中の画素Ａに対応する画像ノードをＮ_Ａと、また画素Ａの近傍画素Ｂに対応する画像ノードをＮ_Ｂとする。そして画素Ａの画素値をＩ_Ａと、画素Ｂの画素値をＩ_Ｂとする。さらに画素Ａと画素Ｂの画像上での距離を返す関数をｄｉｓｔ（Ａ，Ｂ）とする。また、形状モデルの当てはめ結果を参照する関数をＳ（Ａ，Ｂ）とする（詳細は後述）。すると画像ノードＮ_ＡとＮ_Ｂを連結するｎ−ｌｉｎｋに割り当てるエネルギーＥ_ＮＡＮＢは、次式で計算される。

式（１２）においてw_intensityとw_modelは、それぞれ第１項（画像濃度値に関する項）と第２項（形状モデルの当てはめ結果に関する項）の寄与を決定する重み係数である。またσは画素値Ｉ_ＡとＩ_Ｂの間の近接性を評価するための定数である。この定数はステップＳ９０４０の説明で述べた濃度値の分布情報σ_ｌ ^２等とは異なる種類の数値であることに注意されたい。

式（１２）の第２項は、形状モデルの当てはめ結果を参照して値を返すエネルギー項である。本項は、本ステップで得られる領域の境界が、前ステップで取得したマスク画像M'_{lshape(c'l1, c'l2, ..., c'lNshape)}中の解剖学的部位の領域の境界と似た形状になるようにするために導入される。このような効果を持つ項はさまざまな形式で表現可能であるが、本実施形態では以下の項を用いる。

式（１３）において、φ_ｌはM'_{lshape(c'l1, c'l2, ..., c'lNshape)}に対して符号付き距離変換を適用して得られた距離画像であり、ωは正の定数である。ここで符号付き距離変換は、M'_{lshape(c'l1, c'l2, ..., c'lNshape)}中の画素のうち、解剖学的部位の領域に属する画素に正の距離値を、そうでない画素に負の距離値を与える。以上で説明した式（１２）と式（１３）にしたがって、グラフG_l中のすべてのｎ−ｌｉｎｋにエネルギーを割り当てる。

以上で述べた方法で、抽出対象のそれぞれの解剖学的部位について、対応するグラフを生成する。

（Ｓ７０５０：グラフカット法による領域抽出）
ステップＳ７０５０において、第一の領域抽出部６０５０は抽出対象のすべての解剖学的部位の領域を入力画像から抽出する。本ステップでは、特許文献２で開示されているグラフカット法を用いて領域抽出を行う。具体的には、ステップＳ７０４０で作成したグラフG_lをＭａｘ−ｆｌｏｗ／ｍｉｎ−ｃｕｔアルゴリズムを用いて２つの部分グラフに分割し、その分割結果に基づき解剖学的部位lの領域を特定する。この処理を抽出対象のすべての解剖学的部位について行うことで、すべての解剖学的部位の領域を取得する。

ステップＳ７０５０の処理を終えたのち、第一の領域抽出部６０５０は、全領域の抽出結果を３次元のラベル画像の形態でデータサーバ６１０へ出力して保存する。

（Ｓ７０６０：領域抽出結果の表示）
ステップＳ７０６０において、表示制御部６０６０は、ステップＳ７０５０で領域抽出された各部位の領域を入力画像上に表示する。本実施形態では、不図示のＵＩ（user interface）を介してユーザが指定した入力画像の断面画像に、ステップＳ７０４０で得たラベル画像の当該断面を重畳して表示する。

ここまでで説明したように、本実施形態によれば、さまざまな時相の医用画像中に存在する注目領域を、領域拡張法や動的輪郭法といった簡便な領域抽出手法よりも高い精度で抽出することが可能となる。

（変形例３−１）
上記の実施形態では、時相が「非造影」、「早期相」、「門脈相」、「晩期相」の何れの状態であるかを推定し、その結果に基づいて第一の帰属度を算出する際に必要な事前分布を選択していた。しかし、事前分布の選択に十分な分類であれば、画像の時系列での濃度値変化に関する他の状態の分類でもよい。例えば、濃度値の分布が類似する複数の時相を統合して用いてもよい。例えば、「早期相」と「晩期相」を一つにまとめて、「非造影」、「門脈相」、「早期相または晩期相」の何れであるかを推定するようにしてもよい。また、推定する状態は医学の分野で分類される「時相」とは必ずしも一致していなくてもよく、事前分布のバリエーションに基づいて適宜定義してもよい。

また、推定するのは「時相」のような離散的な状態ではなく、造影剤投与からの経過時間のような連続的な状態であってもよい。この場合、データサーバ６１０には、統計データとして、造影剤投与からの経過時間を引数として受け取り、その時間における抽出対象の解剖学的部位の濃度値分布を返す関数（事前分布関数）を保持しておく。また、時相推定のための特徴量を引数として受け取り、造影剤投与からの経過時間を推定する関数（経過時間関数）を保持しておく。

ここで事前分布関数と経過時間関数は、事前に撮像された複数症例の画像から取得する。具体的には、はじめに、それぞれの画像から、造影剤投与からの経過時間、抽出対象のそれぞれの解剖学的部位についての濃度値の分布情報（平均値と分散値）、時相推定のための特徴量を取得する。そして、造影剤投与からの経過時間と各々の濃度値の分布情報の組に対して、経過時間を説明変数、濃度値の分布情報を目的変数とした関数を当てはめることで、事前分布関数を取得する。また、造影剤投与からの経過時間と時相推定のための特徴量の組に対して、時相推定のための特徴量を説明変数、経過時間を目的変数とした関数を当てはめることで、経過時間関数を取得する。

本変形例における処理について説明する。ステップＳ９０６０において時相情報取得部６０３０は、時相推定の代わりに以下の処理を行う。すなわち、時相推定のための特徴量を算出したのち、上記の経過時間関数を用いて時相推定のための特徴量から経過時間を推定する。そしてステップＳ７０２０において第一の帰属度算出部６０４０は、推定した経過時間と事前分布関数を利用して、抽出対象のそれぞれの解剖学的部位に対応する領域について、濃度値の事前分布を推定する。

ここまでで説明したように、本変形例における画像処理装置では、時相の情報（すなわち「非造影」、「早期相」、「門脈相」、「晩期相」など）の代わりとなる情報を推定する。このような方法により、さまざまな時相の医用画像中に存在する注目領域を、領域拡張法や動的輪郭法といった簡便な領域抽出手法よりも高い精度で抽出することが可能となる。

（変形例３−２）
第３の実施形態では、時相が「非造影」、「早期相」、「門脈相」、「晩期相」の何れの状態であるかを推定し、その結果に基づいて第一の帰属度を算出する際に必要な事前分布を選択していた。しかし、第一の帰属度の算出に必要な事前分布を選択することができれば、時相以外の画像の時系列での濃度値変化に関する情報を用いてもよい。

本変形例では、画像の時系列での濃度値変化を表す特徴として、ステップＳ９０５０で算出した時相推定のための特徴量を用いる。そして統計データ取得部６０２０が取得した統計データを用いて、この特徴量から直接（すなわち時相（「非造影」、「早期相」、「門脈相」、「晩期相」）の情報を介することなく）、解剖学的部位の濃度値の事前分布を推定する。

本変形例における統計データには、抽出対象のそれぞれの解剖学的部位についての濃度値の事前分布（平均値と分散値）と、時相推定のための特徴量とを対応付けた関数が含まれている。この関数は、時相推定のための特徴量を引数として与えると、抽出対象のそれぞれの解剖学的部位についての濃度値の事前分布を返す対応表のような役割を果たす。この関数は、存在確率アトラスと同様に、事前に撮像された多数の症例の夫々の画像データと手動で作成された解剖学的部位を表すラベル画像とを用いて作成される。具体的には、まず、画像データとラベル画像を用いて各解剖学的部位の濃度値の分布情報（平均値と分散値）を算出する。次に、時相推定のための解剖学的部位である肝臓と心臓の濃度値の特徴量をそれぞれ軸に割り当てた特徴量空間を定義する。最後に、抽出対象のそれぞれの解剖学的部位についての濃度値の分布情報を目的変数、時相推定のための特徴量を説明変数とした関数を当てはめることで作成する。

本変形例における処理について説明する。ステップＳ７０１０で時相情報取得部６０３０は、ステップＳ９０１０からＳ９０５０までの処理を行う。ステップＳ７０２０において第一の帰属度算出部６０４０は、ステップＳ９０５０で算出した時相推定のための特徴量を用いて、抽出対象のそれぞれの解剖学的部位について濃度値の事前分布を推定する。具体的には、それぞれの解剖学的部位の事前分布を返す濃度値関数に、時相推定のための特徴量を入力して、事前分布の値を得る。そして推定した濃度値の事前分布を用いて、第一の帰属度を算出する。

ここまでで説明したように、本変形例における画像処理装置では、特徴量空間を用いることで、時相の情報（すなわち「非造影」、「早期相」、「門脈相」、「晩期相」など）よりも細かい状態を取得することが可能となる。その結果、さまざまな時相の医用画像中に存在する注目領域を、領域拡張法や動的輪郭法といった簡便な領域抽出手法よりも高い精度で抽出することが可能となる。

（変形例３−３）
上記の実施形態におけるステップＳ７０２０の処理では、時相ごとに定める濃度値の事前分布として、平均値と分散値を得ていた。しかし、平均値と分散値の両方を必ずしも時相に応じて求めなくてもよい。例えば、分散値には共通の固定値を用いることにして、平均値だけを時相に応じて得るようにしてもよい。

（変形例３−４）
上記の実施形態におけるステップＳ７０１０において時相情報取得部６０３０は、入力画像から時相推定のための特徴量を算出して、その結果を利用して時相の情報を取得していた。しかしながら入力画像に時相の情報が付帯している場合には、その情報を利用してもよい。例えば、入力画像がＤＩＣＯＭ形式で保存されており、そのヘッダ部（ＤＩＣＯＭヘッダ）に撮影時刻や時相の情報（「非造影」、「早期相」、「門脈相」、「晩期相」など）が格納されている場合は、それらの情報を取得するようにしてもよい。また入力画像がＤＩＣＯＭ形式で保存されていない場合であっても、入力画像に対応する撮影時刻や時相の情報が何らかの形式でデータサーバ６１０に格納されている場合は、それらの情報を取得するようにしてもよい。

（第４の実施形態）
第３の実施形態では、入力画像の時相の情報を推定し、推定した時相の情報に基づいて、第一の帰属度を算出していた。そして算出した第一の帰属度に基づいて、入力画像から抽出対象の解剖学的部位に対応する領域を抽出していた。本実施形態における画像処理装置も、第３の実施形態と同様の方法で抽出対象の解剖学的部位に対応する領域を抽出する。しかしながら、本実施形態では、第３の実施形態で領域を抽出した後、抽出した領域を参照して、第３の実施形態で使用したものとは異なる帰属度（第二の帰属度）を計算する。そして第二の帰属度に基づいて、再度、抽出対象の解剖学的部位に対応する領域を抽出する。これにより、第３の実施形態における画像処理装置よりも、高い精度で解剖学的部位の領域を抽出する。以下、図１１から図１３までを参照して、本実施形態に係る画像処理装置について、第３の実施形態との相違部分を中心に説明する。

図１１は、本実施形態に係る画像処理システムの構成例を示している。第３の実施形態と同じ部分は同じ番号で示している。図１１から分かる通り、本実施形態に係る画像処理装置７００は、第３の実施形態における画像処理装置６００の構成要素に加えて、第二の帰属度算出部７０６０、第二の領域抽出部７０７０を有する。ここで、画像データ取得部６０１０、統計データ取得部６０２０、時相情報取得部６０３０、第一の帰属度算出部６０４０、第一の領域抽出部６０５０、表示制御部６０６０、データサーバ６１０、表示部６２０については、その処理も第３の実施形態に係る画像処理装置６００と同様である。そのため、説明を省略する。

第一の領域抽出部６０５０は、第３の実施形態に係る画像処理装置６００における第一の領域抽出部と同様の処理を実行する。しかしながら、抽出対象の解剖学的部位に対応する領域を入力画像から抽出した後、その抽出結果を、表示制御部６０６０ではなく、第二の帰属度算出部７０６０に出力する。

第二の帰属度算出部７０６０は、第一の領域抽出部６０５０が抽出した領域を参照して、第一の帰属度算出部６０４０で算出する帰属度とは異なる帰属度（第二の帰属度）を算出する。算出した帰属度は、第二の領域抽出部７０７０に出力される。

第二の領域抽出部７０７０は、第二の帰属度算出部７０６０で算出した帰属度に基づいて、入力画像から、再度、抽出対象の解剖学的領域に対応する領域を抽出する。本処理部では、第一の領域抽出部６０５０とは異なり、反復的に領域抽出処理を実行する。抽出した領域の情報は表示制御部６０６０に出力される。

＜画像処理装置の処理＞
図１２は、本実施形態に係る画像処理装置７００が実施する全体の処理手順を示すフローチャートである。本図において、ステップＳ７０００〜Ｓ７０６０で各処理部が行う処理は第３の実施形態と同様である。そのため、説明を省略する。

（Ｓ１２０６０：各臓器への帰属度（第二の帰属度）の算出）
ステップＳ１２０６０において第二の帰属度算出部７０６０は、抽出対象のそれぞれの解剖学的部位について、第二の帰属度を算出する。第二の帰属度は、ＭＡＰ法を適用して得られる尤度に基づいて算出される。第一の帰属度と同様、第二の帰属度も入力画像の各画素について計算される。すなわち、入力画像の各画素pについて、解剖学的部位lの尤度p_p(v|l)と、「その他の解剖学的部位」の尤度p_p(v|others)の組(p_p(v|l), p_p(v|others))を解剖学的部位lの第二の帰属度として算出する。本ステップでは、入力画像のすべての画素について、この第二の帰属度（尤度の組）を計算する。

抽出対象のそれぞれの解剖学的部位lの尤度p_p(v|l)は、ステップＳ９０４０で実行したＭＡＰ法を再度実行することで計算できる。ＭＡＰ法を実行する際には、式（７）の濃度値の分布情報（平均値l_lと分散値σ_l ²）が初期値として必要となる。本ステップでは、この初期値をステップＳ７０５０で取得した解剖学的部位の領域から計算する。すなわち、入力画像の中で部位lとして抽出された領域内の画素の画素値を参照して、その平均値と分散値を計算する。このようにして得られる値を初期値としてＭＡＰ法を実行することで、Lに含まれるそれぞれの部位について尤度p_p(v|l)を取得する。最後に解剖学的部位lについて尤度の組(p_p(v|l), p_p(v|others))を作成し、これを解剖学的部位lの第二の帰属度とする。

本実施形態では、注目部位の尤度と、「その他の解剖学的部位」の尤度からなる１組の数値を第二の帰属度としたが、第二の帰属度の表現はこれに限定されるものではない。例えば、注目部位の尤度のみを第二の帰属度としてもよい。また、第一の注目部位の尤度、第二の注目部位の尤度、「その他の解剖学的部位」の尤度というように、複数の尤度を１組として第二の帰属度としてもよい。さらに第３の実施形態における第一の帰属度と同様、複数の尤度の最大値や平均値を用いてもよい。

なお本ステップにおける濃度値の分布情報の初期値（式（７）における平均値と分散値の初期値）は、ステップＳ７０２０で取得した濃度値の分布情報、すなわちＥＭアルゴリズムで推定した平均値と分散値を使用してもよい。

（Ｓ１２０７０：グラフカット法の繰り返しによる領域抽出）
ステップＳ１２０７０において、第二の領域抽出部７０７０は抽出対象の解剖学的部位に対応する領域を抽出する。ここで図１３に示したフローチャートを参照して、本ステップにおける処理の詳細を説明する。

（Ｓ１３０１０：抽出した領域の補正）
図１３から分かる通り、ステップＳ１２０７０において第二の領域抽出部７０７０は、ステップＳ１３０１０、Ｓ１３０２０、Ｓ１３０３０、Ｓ１３０４０、Ｓ１３０５０の処理を反復的に繰り返すことで、抽出対象の解剖学的部位に対応する領域を抽出する。以下では、n回目の抽出処理について説明する。また、前回の抽出処理で抽出した解剖学的部位lの領域をマスク画像M^[n] _lに格納した場合について説明する。ここでM^[n] _lは２値画像であり、その画素値M^[n] _l(p)は画素pが部位lの領域に属するかどうかで異なる値をとる。例えば、画素pが部位lの領域に属する場合はM^[n] _l(p) = 1（もしくは0以外の何らかの定数）で、そうでなければM^[n] _l(p) = 0である。なおM^[0] _l、すなわち反復処理の開始時のマスク画像には、図１２のステップＳ７０５０で抽出した解剖学的部位lの領域をマスク画像としたものを使用する。

ステップＳ１３０１０、Ｓ１３０２０、Ｓ１３０３０、Ｓ１３０４０、Ｓ１３０５０の反復処理は、しばしば計算時間を費やす処理になる。そこで計算時間を削減するために、反復処理を実施する前に、反復処理への入力、具体的には入力画像、前ステップで算出した第二の帰属度等の画像サイズを小さくして実施することが可能である。このような画像サイズの変更（画像の縮小）は、公知の画像処理であるサンプリング法、例えば最近傍補間法、線形補間法、三次補間法などで実行できる。

ステップＳ１３０１０において、第二の領域抽出部７０７０は、解剖学的部位の領域についての前回の抽出結果を補正する。本ステップでは、抽出する解剖学的部位に応じてさまざまな処理を実施する。本実施形態では、抽出対象の解剖学的部位の一つである肝臓を例に挙げて説明する。

肝臓には腫瘍や壊死領域など、解剖学的な異常領域が存在することがある。これらの異常領域はしばしば、入力画像中の肝臓領域内に、明るい（または暗い）球形領域として存在していることが多い。このような異常領域が肝臓の辺縁部に存在している場合、グラフカット法ではこれらの異常領域を抽出できないことがある。そこで本ステップにおいて、これらの異常領域を抽出し、グラフカット法で抽出した領域に加える。

本処理は公知の画像処理方法の組み合わせで実施する。最初に肝臓の前回の抽出結果を格納したマスク画像M^[n-1] _liverに対して、図形拡大と収縮処理（モルフォロジ演算であるdilationとerosion処理）を適用する。これにより前回の抽出処理で抽出した肝臓領域に、前回の抽出処理で抽出できなかった肝臓辺縁部の明るい（または暗い）球形領域を加えた、新たなマスク画像を取得できる。このマスク画像をM^[n] _tmpと記述する。

次にM^[n] _tmpとM^[n-1] _liverとの間で画像間差分を行うことで、領域の辺縁部に存在する、明るい（または暗い）球形領域のみを抽出する。具体的には、明るい球形領域を抽出する場合には、M^[n] _tmp(p) = 1かつM^[n-1] _liver = 0かつI(p) ＞ T_tmpを満たす画素pの画素値M^[n] _tmp(p)を新たにM^[n] _tmp(p) = 1にする。逆に、この条件を満たさない画素の画素値をM^[n] _tmp(p) = 0にする。もし暗い球形領域を抽出する場合は、I(p) ＞ T_tmp のかわりにI(p) ＜ T_tmpを条件として同様の閾値処理を行う。こうして得たM^[n] _tmpに対してラベリング処理を適用することで、M^[n] _tmp中の個々の領域を取得する。

最後にM^[n] _tmp中の個々の領域について、領域の妥当性について判定処理を行う。現在、抽出しようとしている領域は球形領域である。そこでM^[n] _tmp中の個々の領域について球形度を計算し、その値が閾値T_sphericityよりも小さい領域をM^[n] _tmp中で削除する。球形度の計算方法も、公知の画像処理手法を利用すればよい。また、この処理で領域の画素数を取捨選択の条件の一つに加えてもよい。すなわち、個々の領域について、その領域を構成する連結画素の画素数を数え、その値がある閾値T_volumeよりも小さい領域を削除する。これにより球形度では球形であると判断される領域（特に小領域）をM^[n] _tmp中で削除することができる。

以上の処理で得たマスク画像M^[n] _tmpを、前回の抽出結果M^[n-1] _liverに統合する。

なお本ステップは選択的に実施することが可能である。すなわち、入力画像中に異常領域が存在しないことが予め分かっている場合は、本ステップの実行を省略することができる。また、本ステップを後述のステップＳ１３０３０の後に実行しても、同様の効果を得ることができる。

本ステップは抽出対象の解剖学的部位に合わせて、適切な処理を実行してよい。さらに本ステップの処理は時相に応じて変更してもよい。例えば、心臓、右腎臓、左腎臓、脾臓を抽出する場合に、特に明るい領域を画像処理手法で選択的に抽出してもよい。

（Ｓ１３０２０：第二の帰属度を用いたグラフの生成）
ステップＳ１３０２０において、第二の領域抽出部７０７０は、ステップＳ１２０６０で算出した第二の帰属度に基づいて、グラフカット法で使用するグラフを生成する。グラフは抽出対象のそれぞれの解剖学的部位に対して１つ作成されることに注意されたい。以下では、解剖学的部位lに対応するグラフG^[n] _lを作成する例を説明する。

グラフは、ステップＳ７０４０と同様の処理で作成される。本ステップとステップＳ７０４０でのグラフ作成の違いは、ｔ−ｌｉｎｋとｎ−ｌｉｎｋに割り当てられるリンクのエネルギーの計算方法である。

ｔ−ｌｉｎｋに割り当てられるエネルギーの計算方法について説明する。今、画素pに対応する画像ノードをN_pとする。そして画像ノードN_pと２つのターミナル・ノードN_s、N_tを連結するｔ−ｌｉｎｋに割り当てるエネルギーを、それぞれE_sp、E_ptとする。するとE_spとE_ptは、前ステップで算出した第二の帰属度(p_p(v|l), p_p(v|others))と、マスク画像M^[n-1] _lから生成される距離画像D^[n] _lの画素値（距離値）D^[n] _l(p)に基づいて計算される。計算の詳細を説明する。

はじめに距離値D^[n] _l(p)を求める。距離値D^[n] _l(p)は公知の画像処理手法である符号付き距離変換をマスク画像M^[n-1] _lに適用することで得られる。ここで使用する符号付き距離変換は、M^[n] _l中の画素のうち部位lの領域に含まれる画素に対して正の距離値を、また含まれない画素に対して負の距離値を与える。

次に第二の帰属度(p_p(v|l), p_p(v|others))と距離値D^[n] _l(p)から、E_spとE_ptを算出する際に参照する値B^[n] _l(p)を計算する。この値は次式により計算される。

ここでD_constは距離値を表す数である。この値を正の値とすることにより、今ステップで抽出される解剖学的部位の領域は、前ステップで得られた領域よりも大きな領域になることが期待できる。また逆に負の値とすることにより、今ステップで得られる領域は前ステップで得られた領域よりも小さな領域になることが期待できる。本実施形態ではD_constを正の数とする。

式（１４）のw_likelihoodとw_distanceはそれぞれ、第二の帰属度に関係する項と距離値に関係する項にかける重み定数である。これらは実験的に決定する。

最後に式（１４）で計算されるB^[n] _l(p)に基づいて、E_spとE_ptを算出する。E_spとE_ptの値の算出式は、B^[n] _l(p)の符号により異なる。B^[n] _l(p)が正の数である場合は、次式に基づいてE_spとE_ptの値を算出する。

逆に、B^[n] _l(p)が負の数である場合は、次式に基づいてE_spとE_ptの値を算出する。

式（１５）または式（１６）にしたがって、グラフG^[n] _l中のすべてのｔ−ｌｉｎｋにエネルギーを割り当てる。

次にｎ−ｌｉｎｋに割り当てるエネルギーの計算方法について説明する。本ステップでもステップＳ７０４０と同様、式（１２）と式（１３）を用いてｎ−ｌｉｎｋに割り当てるエネルギーを計算する。ステップＳ７０４０での計算との違いは、式（１３）中のθ_lの計算方法である。本ステップでは、θ_lをｔ−ｌｉｎｋのエネルギーを計算する際に取得した距離画像D^[n] _lとする。すなわちθ_l= D^[n] _lとする。このような計算方法により、グラフG^[n] _l中のすべてのｎ−ｌｉｎｋにエネルギーを割り当てる。

以上で述べた方法で、抽出対象のそれぞれの解剖学的部位について、対応するグラフG^[n] _lを生成する。

（Ｓ１３０３０：グラフカット法による領域抽出）
ステップＳ１３０３０において、第二の領域抽出部７０７０は、前ステップで算出したそれぞれのグラフG^[n] _lに対してグラフカット法を適用し、入力画像から抽出対象の解剖学的部位に対応する領域を抽出する。本ステップの処理は、ステップＳ７０５０と同じである。

（Ｓ１３０４０：抽出した領域の評価）
ステップＳ１３０４０において、第二の領域抽出部７０７０は、前ステップで抽出した領域を所定の基準で評価する。この評価結果は、次ステップで反復処理の終了判定に使用される。

本ステップでは、３つの基準で抽出した領域を評価する。以下では、n回目の領域抽出処理で取得したすべての部位の領域を格納したラベル画像を使用して、本ステップで評価する基準について説明する。このラベル画像では、画素pが部位lの領域に属する場合は、対応する画素値R_p ^[n]がlとなる。すなわちR_p ^[n] = lである。

本ステップで評価する１つめの基準は、前回の領域抽出で取得した部位lの領域と、今回の領域抽出で取得した領域の差である。具体的には、部位lに関してR_p ^[n-1] = lかつR_p ^[n-1] != l、またはR_p ^[n-1] != lかつR_p ^[n-1] = lとなる画素の数である。なお、"a != b"は２つの数値ａ，ｂが異なることを表す。

２つ目の基準は、前回の領域抽出では部位lの領域として抽出されなかった領域のうち、今回の領域抽出で部位lの領域として抽出された部分である。具体的には、部位lに関してR_p ^[n-1] != lかつR_p ^[n-1] = lとなる画素の数である。

３つ目の基準は、２つ目の基準とは逆に、前回の領域抽出で部位lの領域として抽出された部分のうち、今回の領域抽出では部位lの領域として抽出されなかった部分である。具体的には、部位lに関してR_p ^[n-1] = lかつR_p ^[n-1] != lとなる画素の数である。

解剖学的部位の集合Lに含まれるそれぞれの部位の領域を、上述の３つの基準で評価する。ｎ回目の領域抽出処理で取得した部位lの評価結果をV_{l_1} ^[n] _、V_{l_2} ^[n] _、V_{l_3} ^[n]で表す。

（Ｓ１３０５０：評価値の判定とフローの制御）
ステップＳ１３０５０において、第二の領域抽出部７０７０は前ステップでの評価結果に基づいて、反復的な領域抽出処理の終了を判定する。本ステップでは評価結果を、条件１：V_{l_1} ^[n-1] = 0、条件２：V_{l_2} ^[n-1] = V_{l_2} ^[n] かつV_{l_3} ^[n-1] = V_{l_3} ^[n]、という２つの条件で判定する。そして２つの条件のうち、少なくとも一方の条件が満たされていることを判定する。もし、少なくとも一方の条件が満たされている場合は、ステップＳ１３０６０の処理を実行する。もし条件が満たされていない場合は、ステップＳ１３０１０の処理に戻る。

（Ｓ１３０６０：グラフカット法の反復で得られた領域に基づくグラフの生成）
本ステップと次ステップで、抽出対象の解剖学的部位に対応する領域を入力画像から再度抽出する。ステップＳ１３０６０とＳ１３０７０は、ステップＳ１３０１０、Ｓ１３０２０、Ｓ１３０３０、Ｓ１３０４０、Ｓ１３０５０の反復処理で抽出した領域の精度をさらに高めることを目的として実施する。そのため、この２つのステップは入力画像と同じ画像サイズ、または少なくとも反復処理で利用していた画像サイズよりも大きな画像サイズで実施することが望ましい。もし反復処理において画像サイズを縮小させていた場合は、反復処理の出力、例えば解剖学的部位の抽出結果であるマスク画像の画像サイズを、入力画像と同等の画像サイズに変換する処理を行ったうえで、実施するとよい。このような画像サイズの変換は、公知の画像処理であるサンプリング法、例えば最近傍補間法、線形補間法、三次補間法などで実行できる。

ステップＳ１３０６０において第二の領域抽出部７０７０は、ステップＳ１３０１０からＳ１３０５０までの反復処理で抽出した解剖学的部位の領域に基づいて、グラフカット法で使用するグラフを生成する。以下では、ステップＳ１３０１０からＳ１３０５０までの反復計算がm回目で終了した場合について説明する。これまでと同様、グラフは抽出対象のそれぞれの解剖学的部位に対して１つ作成されることに注意されたい。以下では、解剖学的部位lに対応するグラフG^[m] _lを作成する例を説明する。また、解剖学的部位lに対応する領域がマスク画像M^[m] _lに格納されているとする。

グラフG^[m] _lは、ステップＳ１３０２０と同様の処理で生成される。本ステップとステップＳ１３０１０でのグラフ作成との違いは、ｔ−ｌｉｎｋに割り当てられるエネルギーの計算方法である。ステップＳ１３０２０では式（１３）にしたがって、E_spとE_ptを算出する際に参照する値B^[n] _l(p)を計算していた。本ステップでは、次式に示すように、マスク画像M^[m] _lを符号付き距離変換して得られる距離画像D^[m] _lの距離値D^[m] _l(p)に基づいて計算される。ここで使用する符号付き距離変換は、M^[m] _l中の画素のうち部位lの領域に含まれる画素に対して正の距離値を、また含まれない画素に対して負の距離値を与える。

式（１７）のw'_distanceは距離値に関係する項にかける重み定数である。この値は実験的に決定する。

最後に式（１７）で計算されるB^[m] _l(p)に基づいて、E_spとE_ptを算出する。算出式は、ステップＳ１３０２０での式（１５）と式（１６）と同様の式を用いる。すなわち、B^[m] _l(p)が正の数である場合は式（１５）を、逆にB^[m] _l(p)が負の数である場合は、式（１６）を用いる。ただしいずれの式でも、B^[n] _l(p)ではなくB^[m] _l(p)とする。このような算出を行うことで、グラフG^[m] _l中のすべてのｔ−ｌｉｎｋにエネルギーを割り当てる。

以上で述べた方法で、抽出対象のそれぞれの解剖学的部位について、対応するグラフG^[m] _lを生成する。

（Ｓ１３０７０：グラフカット法）
ステップＳ１３０７０において第二の領域抽出部７０７０は、前ステップで算出したそれぞれのグラフG^[m] _lに対してグラフカット法を適用し、入力画像から抽出対象の解剖学的部位に対応する領域を抽出する。本ステップの処理は、ステップＳ７０５０と同じである。

以上で、ステップＳ１２０７０において第二の画像処理部７０７０が実行する処理の詳細についての説明を終える。

ステップＳ１２０７０の処理を終えた後、第二の領域抽出部７０７０は、抽出した領域の情報を表示制御部６０６０に出力する。そして表示制御部６０６０は、ステップＳ７０６０の処理を実行する。本処理は第３の実施形態に係る画像処理装置６００と同じである。

ここまでで説明したように、本実施形態によれば、さまざまな時相の医用画像中に存在する解剖学的部位を、第３の実施形態よりも高い精度で抽出することが可能となる。

以上のように、第３の実施形態、第４の実施形態によれば、さまざまな時相の医用画像中に存在する注目領域を高い精度で抽出することが可能となる。また、さまざまな時相の医用画像中に存在する解剖学的部位を、領域拡張法や動的輪郭法といった簡便な領域抽出手法よりも高い精度で抽出することができる。

（その他の実施形態）
本発明は、上述の実施形態の１以上の機能を実現するプログラムを、ネットワーク又は記憶媒体を介してシステム又は装置に供給し、そのシステム又は装置のコンピュータにおける１つ以上のプロセッサがプログラムを読出し実行する処理でも実現可能である。また、１以上の機能を実現する回路（例えば、ＡＳＩＣ）によっても実現可能である。

１００：画像処理装置、１１０：データサーバ、１２０：データベース、１３０：表示部、４００：画像処理装置、１０１０：画像データ取得部、１０２０：統計データ取得部、１０３０：部位選択部、１０４０：領域抽出部、１０５０：特徴量算出部、１０６０：時相推定部、１０７０：領域抽出部、１０８０：表示制御部、４０６０：事前分布推定部、４０７０：領域抽出部、６００：画像処理装置、６０１０：画像データ取得部、６０２０：統計データ取得部、６０３０：時相情報取得部、６０４０：第一の帰属度算出部、６０５０：第一の領域抽出部、６０６０：表示制御部、７０６０：第二の帰属度算出部、７０７０：第二の領域抽出部、８０１０：時相推定のための領域抽出部、８０２０：部位選択部、８０３０：特徴量算出部、８０４０：時相推定部

Claims

被検体を造影して撮像した画像を取得する取得手段と、
前記被検体において造影される前記被検体の第一の解剖学的部位を表す第一の領域を前記画像から抽出する第一の領域抽出手段と、
前記第一の領域の濃度値に関する特徴量と前記第一の解剖学的部位の複数の時相における濃度値に関する統計データとの比較結果に基づいて、前記画像の時相を推定する推定手段と、
前記推定手段による推定結果に基づいて、前記被検体の第二の解剖学的部位を表す第二の領域を前記画像から抽出する第二の領域抽出手段と
を備えることを特徴とする画像処理装置。
前記第一の解剖学的部位は、患者間で位置変動が少ない部位であることを特徴とする請求項１に記載の画像処理装置。
前記第一の解剖学的部位は、濃染のパターンが互いに異なる複数の部位であることを特徴とする請求項１又は２に記載の画像処理装置。
前記第一の解剖学的部位は、肝臓および心臓であることを特徴とする請求項３に記載の画像処理装置。
前記第二の解剖学的部位は、前記第一の解剖学的部位と同一の部位を含み、
前記第二の領域抽出手段は、前記第一の領域抽出手段により抽出された第一の解剖学的部位の領域を、より高精度に抽出する処理を含むことを特徴とする請求項１乃至４の何れか１項に記載の画像処理装置。
被検体を造影して撮像した画像を取得する取得手段と、
前記被検体の第一の解剖学的部位を表す第一の領域を前記画像から抽出する第一の領域抽出手段と、
前記第一の領域内における濃度値に関する特徴量に基づいて、前記被検体の前記第一の解剖学的部位とは解剖学的に異なる第二の解剖学的部位を表す第二の領域内における濃度値の事前分布を推定する推定手段と、
前記事前分布を用いて前記第二の領域を前記画像から抽出する第二の領域抽出手段と
を備えることを特徴とする画像処理装置。
被検体を造影して撮像した画像を取得する画像取得手段と、
前記画像から該画像の時系列の濃度値変化を表す特徴を取得する特徴取得手段と、
前記特徴を用いて前記画像から抽出対象の解剖学的部位への帰属度を算出する算出手段と、
前記帰属度を用いて前記画像中に存在する前記解剖学的部位を抽出する抽出手段と、
グラフカット法におけるグラフを生成する生成手段と、を備え、
前記生成手段は、前記グラフの夫々の画像ノードとターミナル・ノードの間のリンクのエネルギーを、夫々の画像ノードに対応する画素に関する前記帰属度に基づいて割り当て、
前記抽出手段は、前記グラフに基づいて前記解剖学的部位を抽出することを特徴とする画像処理装置。
被検体を造影して撮像した画像を取得する画像取得手段と、
前記画像から該画像の時系列の濃度値変化を表す特徴を取得する特徴取得手段と、
前記特徴を用いて前記画像から抽出対象の解剖学的部位への帰属度を算出する算出手段と、
前記帰属度を用いて前記画像中に存在する前記解剖学的部位を抽出する抽出手段と、を備え、
前記算出手段は、前記帰属度を、前記画像内での前記解剖学的部位の存在確率に基づいて算出することを特徴とする画像処理装置。
前記特徴は、時相であることを特徴とする請求項７又は８に記載の画像処理装置。
前記算出手段は、前記帰属度を、前記画像の画像濃度値に基づいて算出することを特徴とする請求項７乃至９の何れか１項に記載の画像処理装置。
被検体を造影して撮像した画像を取得する取得手段と、
前記被検体において造影される前記被検体の所定の解剖学的部位を表す領域を前記画像から抽出する領域抽出手段と、
前記領域の濃度値に関する特徴量と前記所定の解剖学的部位の複数の時相における濃度値に関する統計データとの比較結果に基づいて、前記画像の時相を推定する推定手段と、
を備えることを特徴とする画像処理装置。
画像処理装置の制御方法であって、
取得手段が、被検体を造影して撮像した画像を取得する取得工程と、
第一の領域抽出手段が、前記被検体において造影される前記被検体の第一の解剖学的部位を表す第一の領域を前記画像から抽出する第一の領域抽出工程と、
推定手段が、前記第一の領域の濃度値に関する特徴量と前記第一の解剖学的部位の複数の時相における濃度値に関する統計データとの比較結果に基づいて、前記画像の時相を推定する推定工程と、
第二の領域抽出手段が、前記推定工程による推定結果に基づいて、前記被検体の第二の解剖学的部位を表す第二の領域を前記画像から抽出する第二の領域抽出工程と
を有することを特徴とする画像処理装置の制御方法。
画像処理装置の制御方法であって、
取得手段が、被検体を造影して撮像した画像を取得する取得工程と、
第一の領域抽出手段が、前記被検体の第一の解剖学的部位を表す第一の領域を前記画像から抽出する第一の領域抽出工程と、
推定手段が、前記第一の領域内における濃度値に関する特徴量に基づいて、前記被検体の前記第一の解剖学的部位とは解剖学的に異なる第二の解剖学的部位を表す第二の領域内における濃度値の事前分布を推定する推定工程と、
第二の領域抽出手段が、前記事前分布を用いて前記第二の領域を前記画像から抽出する第二の領域抽出工程と
を有することを特徴とする画像処理装置の制御方法。
画像処理装置の制御方法であって、
画像取得手段が、被検体を造影して撮像した画像を取得する画像取得工程と、
特徴取得手段が、前記画像から該画像の時系列の濃度値変化を表す特徴を取得する特徴取得工程と、
算出手段が、前記特徴を用いて前記画像から抽出対象の解剖学的部位への帰属度を算出する算出工程と、
生成手段が、グラフカット法におけるグラフを生成する生成工程と、
抽出手段が、前記帰属度を用いて前記画像中に存在する前記解剖学的部位を抽出する抽出工程と、を有し、
前記生成工程では、前記グラフの夫々の画像ノードとターミナル・ノードの間のリンクのエネルギーを、夫々の画像ノードに対応する画素に関する前記帰属度に基づいて割り当て、
前記抽出工程では、前記グラフに基づいて前記解剖学的部位を抽出することを特徴とする画像処理装置の制御方法。
画像処理装置の制御方法であって、
画像取得手段が、被検体を造影して撮像した画像を取得する画像取得工程と、
特徴取得手段が、前記画像から該画像の時系列の濃度値変化を表す特徴を取得する特徴取得工程と、
算出手段が、前記特徴を用いて前記画像から抽出対象の解剖学的部位への帰属度を算出する算出工程と、
抽出手段が、前記帰属度を用いて前記画像中に存在する前記解剖学的部位を抽出する抽出工程と、を有し、
前記算出工程では、前記帰属度を、前記画像内での前記解剖学的部位の存在確率に基づいて算出することを特徴とする画像処理装置の制御方法。
画像処理装置の制御方法であって、
取得手段が、被検体を造影して撮像した画像を取得する取得工程と、
領域抽出手段が、前記被検体において造影される前記被検体の所定の解剖学的部位を表す領域を前記画像から抽出する領域抽出工程と、
推定手段が、前記領域の濃度値に関する特徴量と前記所定の解剖学的部位の複数の時相における濃度値に関する統計データとの比較結果に基づいて、前記画像の時相を推定する推定工程と、
を有することを特徴とする画像処理装置の制御方法。
請求項１３乃至１６の何れか１項に記載の画像処理装置の制御方法の各工程をコンピュータに実行させるためのプログラム。