JP2017012341A - 画像処理方法および装置並びにプログラム - Google Patents

Abstract

【課題】複数の画像間で共通する血管分岐をより高い精度で特定できるようにする。【解決手段】第１の血管分岐を含む第１の画像における前記第１の血管分岐の分岐点の周辺に位置する複数の画像領域の各々について、該画像領域の画像の特徴を表す指標を算出する第１の指標算出手段と、第２の血管分岐を含む第２の画像における前記複数の画像領域と実質的に対応関係にある複数の対応領域の各々について、前記指標を算出する第２の指標算出手段と、前記複数の画像領域の各々について算出された前記指標と、前記複数の対応領域の各々について算出された前記指標との類似度を算出する類似度算出手段とを備えた画像処理装置を提供する。【選択図】図１

Description

本発明は、複数の画像間で共通する血管分岐を特定するのに有用な技術に関する。

近年、特に医療分野においては、画像診断の精度向上を目的として、被検体の同一の部位を互いに異なる複数のモダリティ（modality）で撮像し、得られた複数のボリューム（volume）画像（３次元画像）またはその断層像を同時に表示して参照することが行われている。この際、一般的には、同一の部位が互いに重なるように、これら複数のボリューム画像の位置合せが行われる。

ボリューム画像の位置合せの方法は種々存在するが、そのうちの一つとして血管分岐の特徴を用いる方法が提案されている（特許文献１、要約等参照）。この方法では、複数のボリューム画像間で共通する同一の血管分岐を特定し、共通する血管分岐が互いに重なるようにボリューム画像同士を位置合せしている。共通する血管分岐の特定は、例えば、血管分岐の周辺において抽出された血管画像の相互の相関係数に基づいて行われる。

特開２０１５−３９５７８号公報

上記の方法による位置合せの成功率は、位置合せをするボリューム画像間での共通の血管分岐におけるセグメンテーション（segmentation）の精度に依存する。セグメンテーションの誤差は、血管分岐の平面を求める演算における回転、移動の誤差につながる。比較する血管のセグメンテーション同士において相互の相関係数が偶然にも類似してしまった場合には、真に共通する血管分岐の特定に失敗する。

このような事情により、複数の画像間で共通する血管分岐を、より高い精度で特定することが可能な技術が望まれている。

第１の観点発明は、
第１の血管分岐を含む第１の画像における前記第１の血管分岐の分岐点の周辺に位置する複数の画像領域の各々について、該画像領域の画像の特徴を表す指標を算出する処理と、
第２の血管分岐を含む第２の画像における前記複数の画像領域と実質的に対応関係にある複数の対応領域の各々について、該対応領域の画像の特徴を表す指標を算出する処理と、
前記複数の画像領域の各々について算出された指標と、前記複数の対応領域の各々について算出された指標との類似度を算出する処理とをコンピュータ（computer）に実行させる画像処理方法を提供する。

第２の観点の発明は、
第１の血管分岐を含む第１の画像における前記第１の血管分岐の分岐点の周辺に位置する複数の画像領域の各々について、該画像領域の画像の特徴を表す指標を算出する第１の指標算出手段と、
第２の血管分岐を含む第２の画像における前記複数の画像領域と実質的に対応関係にある複数の対応領域の各々について、該対応領域の画像の特徴を表す指標を算出する第２の指標算出手段と、
前記複数の画像領域の各々について算出された指標と、前記複数の対応領域の各々について算出された指標との類似度を算出する類似度算出手段とを備えた画像処理装置を提供する。

第３の観点の発明は、
前記第２の指標算出手段が、前記複数の画像領域の各々について、前記第２の画像における該画像領域と対応する領域及び該領域の近傍の領域に複数の暫定領域を設定し、該複数の暫定領域の各々について該前提領域の画像の特徴を表す指標を算出し、該画像領域の指標に最も近い指標が算出された前記暫定領域を該画像領域の対応領域として特定する上記第２の観点の画像処理装置を提供する。

第４の観点の発明は、
前記複数の画像領域の位置と前記複数の対応領域の位置との近似度を算出する近似度算出手段と、
前記類似度及び前記近似度の大きさに基づいて前記第１の血管分岐と前記第２の血管分岐とが同一の血管分岐を表している蓋然性の程度を示す指標値を算出する蓋然性指標値算出手段とをさらに備えた上記第３の観点の画像処理装置を提供する。

第５の観点の発明は、
前記第１の指標算出手段が、互いに異なる複数の第２の画像の各々について前記類似度を算出し、
前記第２の指標算出手段が、前記複数の第２の画像の各々について前記近似度を算出し、
前記蓋然性指標値算出手段が、前記第１の画像と前記複数の第２の画像の各々とからなる複数の組合せについて前記蓋然性指標値を算出し、
前記複数の組合せの中で最も高い前記蓋然性指標値が算出された組合せを特定する特定手段をさらに備えた上記第４の観点の画像処理装置を提供する。

第６の観点の発明は、
第１のモダリティのボリューム画像から前記第１の画像を生成する第１の生成手段と、
第２のモダリティのボリューム画像から前記複数の第２の画像を生成する第２の生成手段とをさらに備えた上記第５の観点の画像処理装置を提供する。

第７の観点の発明は、
前記第１の生成手段が、前記第１のモダリティのボリューム画像における血管分岐を含むスラブ（slab）に対して画素値の最大値投影または最小値投影を行って前記第１の画像を生成し、
前記第２の生成手段が、前記第２のモダリティのボリューム画像における血管分岐を含むスラブに対して画素値の最大値投影または最小値投影を行って前記第２の画像を生成する上記第６の観点の画像処理装置を提供する。

第８の観点の発明は、
前記第１の生成手段が、前記第１のモダリティの種類に応じて前記最大値投影及び最小値投影のいずれかを行い、
前記第２の生成手段が、前記第２のモダリティの種類に応じて前記最大値投影及び最小値投影のいずれかを行う上記第６の観点または第７の観点の画像処理装置を提供する。

第９の観点の発明は、
前記特定手段により特定された組合せに対応する血管分岐同士が重なるように、前記第１のモダリティのボリューム画像と前記第２のモダリティのボリューム画像との位置合せを行う位置合せ手段をさらに備えた上記第６の観点から第８の観点のいずれか一つの観点の画像処理装置を提供する。

第１０の観点の発明は、
前記第１のモダリティのボリューム画像から生成された血管分岐を含む複数の画像の中から１つ画像を前記第１の画像として指定する指定手段をさらに備えた上記第６の観点から第９の観点のいずれか一つの観点の画像処理装置を提供する。

第１１の観点の発明は、
前記指定手段が、操作者の操作に応じて前記第１の画像を指定する上記第１０の観点の画像処理装置を提供する。

第１２の観点の発明は、
前記第１のモダリティが、超音波撮像装置であり、
前記第２のモダリティは、磁気共鳴撮像装置である上記第６の観点から第１１の観点のいずれか一つの観点の画像処理装置を提供する。

第１３の観点の発明は、
前記第１のモダリティが、超音波撮像装置であり、
前記第２のモダリティが、放射線断層撮像装置である上記第６の観点から第１１の観点のいずれか一つの観点の画像処理装置を提供する。

第１４の観点の発明は、
前記ボリューム画像が、肝臓を含む画像である上記第６の観点から第１３の観点のいずれか一つの観点の画像処理装置を提供する。

第１５の観点の発明は、
前記第１の画像における画素値に基づいて前記複数の画像領域を設定する設定手段をさらに備えた上記第２の観点から第１４の観点のいずれか一つの観点の画像処理装置を提供する。

第１６の観点の発明は、
前記設定手段が、前記第１の画像を複数の部分画像に分割し、前記複数の部分画像の各々について該部分画像の特徴度を算出し、該特徴度の大きさに基づいて前記複数の画像領域を設定する上記第１５の観点の画像処理装置を提供する。

第１７の観点の発明は、
前記設定手段が、前記特徴度として、前記部分画像内の複数の画素における角らしさ（corner-ness）の程度が反映された値を算出する上記第１６の観点の画像処理装置を提供する。

第１８の観点の発明は、
前記設定手段が、前記複数の画像領域を前記第１の画像における前記第１の血管分岐の分岐点の周りに分散して設定する上記第１５の観点から第１７の観点のいずれか一つの観点の画像処理装置を提供する。

第１９の観点の発明は、
前記第１及び第２の指標算出手段が、画像の特徴を表す指標として複数の画素における自己相似記述子（self-similarity descriptor）を算出する上記第２の観点から第１８の観点のいずれか一つの観点の画像処理装置を提供する。

第２０の観点の発明は、
コンピュータを、上記第２の観点から第１９の観点のいずれか一つの観点の画像処理装置として機能させるためのプログラム（program）を提供する。

上記観点の発明によれば、上記構成により、画像における血管分岐の分岐点の周辺に位置する複数の領域の各々について画像の特徴を表す指標を算出し、複数の画像間でのこれら指標の類似度を算出するので、血管分岐に係るより多くの情報を利用して、複数の画像にそれぞれ含まれる血管分岐が、共通する同一の血管分岐を表している蓋然性の程度をより高い精度で求めることができる。その結果、複数の画像間で共通する血管分岐を、より高い精度で特定することが可能になる。

本実施形態に係る画像処理装置の構成を概略的に示す機能ブロック図である。 本実施形態に係る画像処理装置における処理の流れを示すフロー図である。 分岐平面画像を生成する方法を説明するための図である。 クエリパッチを設定する方法を説明するための図である。 自己相似記述子を計算する方法を説明するための図である。 類似パッチを探索する方法を説明するための図である。 パッチの記述子の類似度を算出する例を示す図である。 パッチのジオメトリの近似度を算出する例を示す図である。 血管分岐のベストマッチングを特定する例と、ボリューム画像の位置合せの例とを示す図である。

以下、発明の実施形態について説明する。なお、これにより発明は限定されない。

まず、本実施形態に係る画像処理装置の概要について説明する。

本実施形態に係る画像処理装置は、同一の被検者における同一の部位を含む２つのボリューム画像を位置合せし、その後、任意のスライス（slice）を表す断層画像を再構成して表示する装置である。位置合せは、次のような手順で行う。上記２つのボリューム画像の各々について血管を抽出して血管分岐を検出する。上記２つのボリューム画像間で共通する同一の血管分岐を特定する。これらの血管分岐が互いに重なるようにボリューム画像の座標変換を行う。

また本実施形態では、この位置合せにおいて、共通する同一の血管分岐を特定する際に、比較対象となる血管分岐の分岐点の周辺における組織構造の類似性を参照する。この方法は、人間が血管分岐の画像を目視して血管分岐の同一性を判断するときの思考に似ており、精度の高い特定が可能になる。具体的には、次のような処理を行う。

まず、比較対象となる一方のボリューム画像から検出された血管分岐の分岐点の周辺に複数の画像領域を設定し、画像領域ごとに画像の特徴を表す指標を算出する。また、比較対象となる他方の血管分岐の周辺に上記複数の画像領域に対応した位置の近傍で、算出された指標に類似した指標を持つ対応領域を探索する。そして、一方の血管分岐側の各画像領域と他方の血管分岐側の各対応領域との間で、指標の類似度とジオメトリ（geometry）の近似度とを算出し、これら類似度及び近似度から比較対象の血管分岐が共通する同一の血管分岐である蓋然性の程度を示す蓋然性指標値を求める。このような処理を、比較対象となる他方の血管分岐を変えながら行い、その都度、蓋然性指標値を算出する。例えば、一方のボリューム画像から検出された血管分岐がｍ個、他方のボリューム画像から検出された血管分岐がｎ個の場合、ｍ×ｎ通りの組合せに対して蓋然性指標値を算出する。そして、最も高い蓋然性指標値が算出された血管分岐の組合せに対して、これら血管分岐は共通する同一の血管分岐であるとの判定を行う。

画像処理装置は、この判定結果を利用することで、ボリューム画像間に共通する同一の血管分岐を高い精度で特定することができ、ボリューム画像間の位置合せを高い精度で行うことができる。

次に、本実施形態に係る画像処理装置の詳細について説明する。

図１は、本実施形態に係る画像処理装置１の構成を概略的に示す機能ブロック（block）図である。画像処理装置１は、例えば、コンピュータ（computer）ＣＰに所定のプログラム（program）を実行させることにより実現させることができる。コンピュータは、一般的に、中央演算処理装置、記憶装置、入力装置、画像出力装置、入出力インタフェース（interface）等を備えている。

図１に示すように、画像処理装置１は、画像取得部２、血管分岐検出部３、分岐平面画像生成部４、平滑化処理部５、分岐平面画像表示部６、血管分岐指定部７、クエリパッチ（query patch）設定部８、クエリパッチ指標算出部９、類似パッチ探索部１０、蓋然性指標値算出部１１、ベストマッチ（best match）分岐特定部１２、位置合せ部１３、及び断層像再構成表示部１４を有している。

なお、血管分岐検出部３、分岐平面画像生成部４、及び平滑化処理部５は、発明における第１及び第２の生成手段の一例である。血管分岐指定部７は、発明における指定手段の一例である。クエリパッチ設定部８は、発明における設定手段の一例である。クエリパッチ指標算出部９は、発明における第１の指標算出手段の一例である。類似パッチ探索部１０は、発明における第２の指標算出手段の一例である。蓋然性指標値算出部１１は、発明における類似度算出手段、近似度算出手段、及び蓋然性指標値算出手段の一例である。ベストマッチ分岐特定部１２は、発明における特定手段の一例である。位置合せ部１３は、発明における位置合せ手段の一例である。

画像取得部２は、位置合せの対象となる複数のボリューム画像を取得する。

血管分岐検出部３は、取得された複数のボリューム画像の各々について、血管分岐を検出する。血管分岐は、その分岐点の位置と分岐点から枝分かれする血管部分に対応した分岐ベクトル（vector）とで規定する。

分岐平面画像生成部４は、ボリューム画像において検出された個々の血管分岐について、分岐平面画像を生成する。

平滑化処理部５は、同一の血管分岐を表す複数の分岐平面画像を加算したものに対して平滑化処理を施し、分岐平面画像の画質を向上させる。

分岐平面画像表示部６は、生成された分岐平面画像を、画像出力装置の画面に表示させる。

血管分岐指定部７は、位置合せの対象である複数のボリューム画像間で共通する血管分岐を特定する際の当該血管分岐として、一方のボリューム画像で検出された血管分岐の中から一つを指定する。

クエリパッチ設定部８は、一方のボリューム画像側で選択された分岐平面画像に対して、複数のクエリパッチを設定する。

クエリパッチ指標算出部９は、複数のクエリパッチの各々について、クエリパッチの画像の特徴を表す指標を算出する。

類似パッチ探索部１０は、他方のボリューム画像側における複数の分岐平面画像の各々において、クエリパッチごとにクエリパッチの指標と最も類似する指標を持つ類似パッチを探索する。

蓋然性指標値算出部１１は、他方のボリューム画像側における複数の分岐平面画像の各々について、同一分岐蓋然性指標値を算出する。同一分岐蓋然性指標値は、一方のボリューム画像側で指定された分岐平面画像の血管分岐と他方のボリューム画像側のある分岐平面画像の血管分岐とが共通する同一の血管分岐である蓋然性の程度を示す指標値である。

ベストマッチ分岐特定部１２は、他方のボリューム画像側の血管分岐の中で同一分岐蓋然度が最も高い血管分岐を、一方のボリューム画像側で指定された血管分岐と共通する同一の血管分岐であると特定する。

位置合せ部１３は、一方のボリューム画像と他方のボリューム画像とを、ベストマッチの血管分岐同士が重なるように位置合せを行う。

断層像再構成表示部１４は、位置合せされた複数のボリューム画像において任意の断層像を再構成し、画面に表示する。

これより、本実施形態に係る画像処理装置１における処理の流れについて説明する。

図２は、本実施形態に係る画像処理装置１における処理の流れを示すフロー（flow）図である。

ステップ（step）Ｓ１では、ボリューム画像を取得する。具体的には、画像取得部２が、同一の被検者における同一の部位を含むＵＳボリューム画像Ｖ_US及びＭＲボリューム画像Ｖ_MRを取得する。ＵＳボリューム画像は、超音波撮像装置により得られたボリューム画像である。また、ＭＲボリューム画像は、磁気共鳴撮像装置により得られたボリューム画像である。ボリューム画像の取得は、例えば、内部または外部の記憶装置や記憶媒体等からの読み込みにより行われる。また、ボリューム画像に含まれる部位は、血管を含む変形性の高い臓器等が考えられる。本例では、この部位として肝臓を想定する。

ステップＳ２では、血管分岐を検出する。具体的には、血管分岐検出部３が、ＵＳボリューム画像及びＭＲボリューム画像の各々について、次のような処理を順次行い、血管分岐を検出する。まず、ボリューム画像において既知の方法により血管を表す画像を抽出する。次に、抽出された血管の画像に対してスムージング（smoothing）処理、細線化処理等を行って血管の略中心線で形成されるツリー（tree）構造を求める。次いで、このツリー構造の枝を上流から下流へと追跡することにより、血管の分岐点を検出する。血管分岐は、その分岐点の位置と分岐点から下流側に枝分かれする血管部分に対応した分岐ベクトルとで規定する。

これにより、描画の態様、すなわち、物質とその画素値との対応関係が互いに異なる複数のボリューム画像において検出された血管分岐を、同一種類の画像として扱うことができる。すなわち、モダリティの違いを気にせずに血管分岐同士を比較することができる。例えば、超音波撮像装置により得られた超音波画像における血管分岐と、Ｘ線ＣＴ装置により得られたＣＴ画像における血管分岐とを区別して扱う必要がなくなる。

ステップＳ３では、分岐平面画像を生成する。

図３は、分岐平面画像を生成する方法を説明するための図である。具体的には、分岐平面画像生成部４が、ＵＳボリューム画像及びＭＲボリューム画像の各々について、次のような処理を順次行い、分岐平面画像を生成する。図３に示すように、まず、ボリューム画像において検出された個々の血管分岐について、分岐点ＢＰ及び分岐ベクトルＢＶが含まれる部分平面ＢＨを規定する。次に、その部分平面ＢＨに所定の厚みΔｔを持たせたスラブＳＬＢを設定する。次いで、そのスラブＳＬＢ内の画素値（輝度値）を厚み方向ｔに最大値投影（ＭＩＰ）もしくは最小値投影（ＭｉｎＩＰ）して、分岐平面画像ＢＪを生成する。なお、造影剤が注入された被検体の肝臓を表すＸ線ＣＴ画像など、血管がその周辺の組織の平均的な画素値（輝度値）よりも高い画素値で表されるような場合には、投影処理の種類として最大値投影処理を用いる。一方、造影剤が注入された被検体の肝臓を表すＭＲ画像や超音波ｂモード画像など、血管がその周辺の組織の平均的な画素値（輝度値）よりも低い画素値で表されるような場合には、投影処理の種類として最小値投影処理を用いる。厚みΔｔは、例えば、血管の直径と同程度、具体的には実空間上で２〜７ｍｍに対応する程度とする。

これにより、血管抽出の精度の問題で、検出された分岐点及び分岐ベクトルに誤差が含まれる場合にも、ある程度の抽出のずれを補うことが可能になる。

ステップＳ４では、平滑化処理を行う。具体的には、平滑化処理部５が、同一の血管分岐を表す複数の分岐平面画像を加算したものに対して平滑化処理を施し、分岐平面画像の画質を向上させる。モダリティによっては、同一の撮像領域に対して複数の画像が取得される場合がある。そこで、これをうまく利用し、上述のように画像の加算及び平滑化処理を行うことにより、分岐平面画像におけるコントラスト（contrast）に対するノイズ（noise）を抑え、画質を向上させることができる。平滑化処理には、例えば、ガウシアンフィルタ（Gaussian filter）や移動平均フィルタ等を用いる。なお、この画像加算及び平滑化処理は、必須の処理ではなく、必要に応じて実施すればよい。

ステップＳ５では、分岐平面画像を表示する。具体的には、分岐平面画像表示部６が、ＵＳボリューム画像側の分岐平面画像とＭＲボリューム画像側の分岐平面画像とを、画像出力装置の画面に表示させる。

これにより、操作者は、ＵＳボリューム画像とＭＲボリューム画像のそれぞれに、どのような血管分岐が検出されたかを目視で確認することができる。

ステップＳ６では、血管分岐を指定する。具体的には、血管分岐指定部７が、位置合せの対象であるＵＳボリューム画像とＭＲボリューム画像との間で共通する血管分岐を特定する際の当該血管分岐として、ＵＳボリューム画像で検出された血管分岐の中から一つを指定する。指定は、例えば、操作者が入力装置を操作して、表示された分岐平面画像の中から所望の一つを選択することにより行われる。

ステップＳ７では、クエリパッチを設定する。

図４は、クエリパッチを設定する方法を説明するための図である。具体的には、クエリパッチ設定部８が、ＵＳボリューム画像側で選択された分岐平面画像に対して、複数のクエリパッチを設定する。複数のクエリパッチは、指定された血管分岐を特徴づける部分的な画像領域（パッチ）の組合せであり、ＭＲボリューム画像側で当該血管分岐と共通するもう一方の血管分岐を特定する際に用いる。複数のクエリパッチは、次のような手順により設定される。

まず、図４の右側の分岐平面画像ＢＪ_US1に示すように、選択された分岐平面画像を複数の矩形領域に分割する。次に、分割された矩形領域ごとに、その矩形領域の画像が特徴的（複雑またはユニーク（unique））であればあるほど値が大きくなるような特徴度を算出する。このような特徴度としては、例えば、矩形領域の画像内の複数の画素における角らしさ（corner-ness）の程度が反映された値を考えることができる。より詳しくは、例えば、矩形領域内の画像を構成する各画素に対して、その画素が描画上の角を表している蓋然性の程度を表す特徴量を算出し、それを平均化したものを考えることができる。この特徴量は、コーナーネス（Corner-ness）と呼ばれる。コーナーネスは、その画素が描画上の角を表している蓋然性、すなわち角らしさの程度を示す特徴量であり、角らしさの程度が大きいほど大きな値をとる。コーナーネスは、次のように定義される。

コーナーネスの値は、明示的に算出された固有値なしの主曲率（principal curvatures）のヘシアン行列Ｈ（Hessian）のトレース（trace；Tr）と行列式（determinant；Det）との比により算出される。

ここで、Ｄxx，Ｄyy，及びＤxyは、ヘシアン行列Ｈのエレメント（要素）であり、導関数演算子（derivative operator）としてのソーベル・フィルタ（Sobel filter）を用いて算出される。ソーベル・フィルタは、エッジ（edge）検出用として一般的なフィルタである。

図４の左側の画像Ｃ_mapは、分岐平面画像ＢＪ_US1の各画素に対するコーナーネスの値をその画素の輝度値としたときのマップ（map）である。

次いで、分岐平面画像において、コーナーネスの平均値が比較的高い複数の矩形領域を、クエリパッチに設定する。このとき、クエリパッチは、その位置に偏りが生じないよう、分岐点ＢＰの周りに分散して配置される。本例では、分岐平面画像を第１象限Ａ１〜第４象限Ａ４に分割し、各象限において、コーナーネスの平均値が閾値以上であって上位から２番目までの矩形領域を、クエリパッチに設定する。すなわち、分岐点の周辺において最大で８個のクエリパッチが設定される。

図４の例では、第１象限Ａ１で１つ、第２象限Ａ２で２つ、第３象限Ａ３で１つ、第４象限Ａ４で２つの合計６つのクエリパッチＱＰが設定されている。

なお、コーナーネスの詳細については、非特許文献：Carlos R del-Blancoa, et.al: Motion estimation through efficient matching of a reduced number of reliable singular points SPIE-IS&T/ Vol. 6811 68110N-1を参照されたい。

ステップＳ８では、クエリパッチの画像の特徴を表す指標を算出する。具体的には、クエリパッチ指標算出部９が、複数のクエリパッチの各々について、クエリパッチの画像を構成する画素ごとにその画素の特徴を表す記述子（descriptor）を算出する。本例では、記述子として、自己相似記述子（self-similarity descriptor）を用いる。自己相似記述子は、例えば、次のように定義される。

図５は、自己相似記述子を説明するための図である。図５の左上図は、自己相似記述子の計算方法の概念図を示している。

今、座標（ｕ，ｖ）の画素Ｏについて自己相似記述子ｄ（ｕ，ｖ）を計算することを考える。自己相似記述子ｄ（ｕ，ｖ）は、画素Ｏを中心とした相対的に小さい幅ａ×高さｂの第１の領域ＲＣの画像と、画素Ｏを中心とし第１の領域ＲＣを囲む相対的に大きい、幅ｗ×高さｈの第２の領域ＲＳの画像との比較によって求められる。第１の領域ＲＣのサイズ（size）ａ×ｂは、典型的には、５×５画素である。また、第２の領域ＲＳのサイズｗ×ｈは、典型的には、４０×４０画素である。上記の比較には、第１の領域ＲＣの画像と第２の領域ＲＳの画像との間の画素値（輝度値）のＳＳＤ（Sum of Squared Difference)が用いられる。ＳＳＤを求める際には、第１の領域ＲＣの画像をテンプレートとして第２の領域ＲＳの画像上をラスタスキャン（raster scan）させる。そして、テンプレート（template）の位置ごとに、テンプレート内の各位置において同じ位置の画素間で画素値（輝度値）の差の２乗を求めてその和を取る。これにより、ＳＳＤは、テンプレートの位置の数分の要素を持つ配列として得られる。最終的に、画素（ｕ，ｖ）におけるＳＳＤ（ｕ，ｖ）は標準化され、“相関表面（correlation surface）”Ｓ（ｕ，ｖ）に変換される。これを自己相似記述子ｄ（ｕ，ｖ）とする。

ここで、σ²はローカル（local）またはグローバルノイズ（global noise）の見積値、Ｎは自己相似記述子が算出される領域である。

図５の下列は、分岐画像平面ＢＪ_US1と、その画素ごとの自己相似記述子ｄの絶対値を輝度値にして表したマップｄ_mapとを示している。この図において、自己相似記述子ｄは、対数極座標プロット（log polar plot）の形式で可視化されている。

なお、自己相似記述子の詳細については、非特許文献：Shechtman，E.Irani，M.et .al: Matching local self-similarities across images and videos. In: CVPR 2007，pp .1-8 IEEE (2007)を参照されたい。

ステップＳ９では、類似パッチを探索する。具体的には、類似パッチ探索部１０が、ＭＲボリューム画像側における複数の分岐平面画像の各々において、クエリパッチごとにクエリパッチの記述子と最も類似する記述子を持つ類似パッチを探索する。類似パッチの探索範囲は、探索ミスを抑え、効率を高めるため、血管分岐との位置関係がクエリパッチと同じになる対応位置の近傍（周辺）に限定する。

図６は、類似パッチを散策する様子を模式的に示す図である。例えば、図６に示すように、分岐平面画像ＢＪ_US1に複数のクエリパッチＱＰ₁，ＱＰ₂，ＱＰ₃があったとして、そのうちの１つのクエリパッチＱＰ₁に着目する。このクエリパッチＱＰ₁は、Ｏを中心とした幅ｗｐ×高さｈｐの矩形領域である。ここで、ＭＲボリューム画像側の分岐平面画像ＢＪ_MR1において、血管分岐を基準としてクエリパッチＱＰ₁と位置的に対応する対応領域ＱＰ₁’を規定する。さらに、この対応領域ＱＰ₁’の中心Ｏ‘と同じ中心を有し、幅及び高さがそれぞれ２×ｗｐ、２×ｈｐである拡張領域ＥＰ₁を規定する。次に、サイズがこのクエリパッチＱＰ₁と同じである幅ｗｐ×高さｈｐの暫定パッチＴＰを、その中心が拡張領域ＥＰ₁内で走査するように順次設定し、その都度、暫定パッチＴＰの各画素の記述子ｄの配列を算出する。クエリパッチＱＰ₁の各画素の記述子ｄの配列に最も類似する記述子の配列が算出された暫定パッチＴＰを、このクエリパッチＱＰ₁の類似パッチに決定する。このような処理を、すべてのクエリパッチについて行う。

ステップＳ１０では、同一分岐蓋然性指標値を算出する。具体的には、蓋然性指標値算出部１１が、ＭＲボリューム画像側における複数の分岐平面画像の各々について、同一分岐蓋然性指標値を算出する。

同一分岐蓋然性指標値は、ＵＳボリューム画像側で指定された分岐平面画像の血管分岐とＭＲボリューム画像側のある分岐平面画像の血管分岐とが共通する同一の血管分岐である蓋然性の程度を示す指標値である。同一分岐蓋然性指標値は、複数のクエリパッチの群ｙとこれらに対応する複数の類似パッチの群ｘとを比較した結果から算出する。具体的には、複数のクエリパッチｙと複数の類似パッチの群ｘとの間で記述子類似度とジオメトリ近似度とを算出し、これらの類似度及び近似度に基づいて同一分岐蓋然性指標値Ｐ（ｘ，ｙ）を算出する。なお、ジオメトリ近似度は、幾何学的な類似度ということもできる。

ここで、ｄⁱ _yは、複数のクエリパッチの群ｙにおけるｉ番目のクエリパッチの記述子を示し、ｌⁱ _yは、このクエリパッチの絶対座標における位置を示すものとする。同様に、ｄⁱ _xは、複数の類似パッチの群ｘにおけるｉ番目の類似パッチの記述子を示し、ｌⁱ _xは、この類似パッチの絶対座標における位置を示すものとする。また、ｃ_yは複数のクエリパッチの群ｙの中心点（“origin”point）、ｃ_xは複数の類似パッチの群ｘの中心点を示すものとする。

すると、複数のクエリパッチの群ｙと複数の類似パッチの群ｘとの間での同一分岐蓋然性指標値Ｐ（ｘ，ｙ）は、次式によって表すことができる。

ここで、Ｐ（ｄⁱ _y｜ｄⁱ _x）は、ガウシアン分布を用いた記述子間の類似度である（記述子類似度という）。また、Ｐ（ｌⁱ _y｜ｌⁱ _x，ｃ_y，ｃ_x）は、パッチの群の相対的なジオメトリ間の近似度である（以下、ジオメトリ近似度という）。

本方法において、すべてのパッチに対するジオメトリ近似度Ｐ（ｌⁱ _y｜ｌⁱ _x，ｃ_y，ｃ_x）の直接的な計算は、共通する分岐平面画像における回転及び平行移動によるミスアライメントをそれぞれ±５度以下、±５ｍｍ以下と仮定して行われる。

ジオメトリ近似度Ｐ（ｌⁱ _y｜ｌⁱ _x，ｃ_y，ｃ_x）の値は、パッチの群における各パッチの中心の周辺におけるサイズ固定の矩形領域には１が設定され、それ以外の領域には０が設定される。

最終的な同一分岐蓋然性指標値Ｐ（ｘ，ｙ）は、ジオメトリ近似度Ｐ（ｌⁱ _y｜ｌⁱ _x，ｃ_y，ｃ_x）が１となる画素についての記述子類似度Ｐ（ｄⁱ _y｜ｄⁱ _x）の総和として求められる。

図７は、記述子類似度のマッピングの例を示す図である。本例は、ＵＳボリューム画像側の分岐平面画像ＢＪ_US1にテンプレートパッチＰ₁を設定し、ＭＲボリューム画像側の分岐平面画像ＢＪ_MR1，ＢＪ_MR2における各領域について、このテンプレートパッチＰ₁の記述子類似度Ｐ（ｄⁱ _y｜ｄⁱ _x）を求め、マッピングしたものである。なお、分岐平面画像ＢＪ_MR1は、分岐平面画像ＢＪ_US1と共通する同一の血管分岐を含んでいる。

この例においては、共通する分岐平面画像における回転（rotation）または平行移動（shift）によるミスマッチ（mismatch）をある程度補うため、ガウシアン分布パラメータの分散または平方偏差は調整されている。

図７の右側の画像は、分岐平面画像ＢＪ_MR1，ＢＪ_MR2に記述子類似度のマップＬＤ_map1，ＬＤ_map2を重ね合せたものである。マップＬＤ_map1において枠で囲んだ部分が、特に高い類似度を持つことが分かる。

図８は、同一分岐蓋然性指標値の算出例を示す図である。図８の上図は、ＵＳボリューム画像側の分岐平面画像ＢＪ_US1において血管分岐の分岐点ＢＰの周りに設定されたクエリパッチＱＰを示している。図８の中列は、ＭＲボリューム画像側の分岐平面画像ＢＪ_MR1，ＢＪ_MR2，ＢＪ_MR3を示している。図８の下列は、各パッチの中心周辺の画素位置での同一分岐蓋然性指標値を輝度値で表したマップＬＧ_map1，ＬＧ_map2，ＬＧ_map3を示している。

矩形領域のサイズは、自動抽出による誤差によって見積もられる回転または平行移動の誤差量に基づいて決定される。

ステップＳ１１では、ベストマッチの血管分岐を特定する。具体的には、ベストマッチ血管分岐特定部１２が、ＭＲボリューム画像側の血管分岐の中で同一分岐蓋然度が最も高い血管分岐を、ＵＳボリューム画像側で指定された血管分岐と共通する同一の血管分岐であると特定する。

図９は、共通するまたはベストマッチングの血管分岐の自動探索の様子を示している。ＵＳボリューム画像側で検出されたすべての血管分岐の分岐平面画像ＢＪ_US,ALLの中から画像ＢＪ_US1が指定されている。ＭＲボリューム画像側で検出されたすべての血管分岐の分岐平面画像ＢＪ_MR,ALLの中から最も高い同一分岐蓋然性指標値が算出された画像ＢＪ_MR1が特定されている。

ステップＳ１２では、ボリューム画像の位置合せを行う。具体的には、位置合せ部１３が、ＵＳボリューム画像とＭＲボリューム画像とを、ベストマッチの血管分岐同士が重なるように位置合せを行う。位置合せは、一方または両方のボリューム画像の座標を変換することにより行われる。

ステップＳ１３では、任意の断層像を表示する。具体的には、断層像再構成表示部１４が、位置合せされたＵＳボリューム画像とＭＲボリューム画像とにおいて、任意の断層像を再構成し、画面に表示する。

図９の右下に、位置合せ後の画像Ｇ_US,MRの例を示す。

以上、本実施形態によれば、互いに異なる２つの画像それぞれに含まれる血管分岐が同一である蓋然性の程度を求める際に、画像に含まれる血管分岐の分岐点の周辺に位置する複数の部分的な画像領域の特徴を用いるので、血管分岐に係るより多くの情報を利用することができ、上記蓋然性の程度をより高い精度で求めることができる。その結果、２つの画像間で共通する同一の血管分岐が重なるように画像の位置合せを行う方法において、より高い精度で位置合せを行うことができる。

また、本実施形態の位置合せ方法は、特許文献１（特開２０１５−３９５７８号公報）による位置合せ方法、すなわち、抽出された血管の相互相関類似測定に基づいた位置合せ方法のパフォーマンスを改善させることができる。

また、本実施形態の位置合せ方法では、血管分岐の検出ではある程度のレベル（level）で偽の血管分岐を誤検出するにもかかわらず、位置合せの成功率は改善されることが確認された。つまり、本提案の方法は、血管の抽出（セグメンテーション）が完璧でない場合においても、位置合せをよりロバスト（robust）に行うことができる。

また、本実施形態の位置合せ方法は、血管を画素の輝度値で差別化できるすべてのモダリティの画像に対して適用可能である。すなわち、位置合せを行う２つの画像の組合せとしては、超音波撮像装置によるＵＳ画像と磁気共鳴撮像装置によるＭＲ画像との組合せだけでなく、ＵＳ画像と放射線断層撮像装置によるＣＴ画像や、ＭＲ画像とＣＴ画像の組合せなど、あらゆるモダリティの画像に適用できる。

また、本実施形態の方法は、構成がシンプル（simple）であり、計算コスト（cost）を抑えることができる。

また、本実施形態は、画像処理装置であるが、コンピュータをこのような画像処理装置として機能させるためのプログラムや、当該プログラムが記憶されたコンピュータ読み取り可能な記憶媒体もまた発明の実施形態の一例である。なお、当該記憶媒体としては、一過性のものだけでなく、非一過性のものも考えることができる。

なお、発明の実施形態は、その趣旨を逸脱しない範囲において、種々の変更・追加等が可能である。

１ 画像処理装置
２ 画像取得部
３ 血管分岐検出部
４ 分岐平面画像生成部
５ 平滑化部
６ 分岐平面画像表示部
７ 血管分岐指定部
８ クエリパッチ設定部
９ クエリパッチ指標算出部
１０ 類似パッチ探索部
１１ 蓋然性指標値算出部
１２ ベストマッチ分岐特定部
１３ 位置合せ部
１４ 断層像再構成・表示部

Claims (20)

1. 第１の血管分岐を含む第１の画像における前記第１の血管分岐の分岐点の周辺に位置する複数の画像領域の各々について、該画像領域の画像の特徴を表す指標を算出する処理と、
   第２の血管分岐を含む第２の画像における前記複数の画像領域と実質的に対応関係にある複数の対応領域の各々について、該対応領域の画像の特徴を表す指標を算出する処理と、
   前記複数の画像領域の各々について算出された指標と、前記複数の対応領域の各々について算出された指標との類似度を算出する処理とをコンピュータに実行させる画像処理方法。
2. 第１の血管分岐を含む第１の画像における前記第１の血管分岐の分岐点の周辺に位置する複数の画像領域の各々について、該画像領域の画像の特徴を表す指標を算出する第１の指標算出手段と、
   第２の血管分岐を含む第２の画像における前記複数の画像領域と実質的に対応関係にある複数の対応領域の各々について、該対応領域の画像の特徴を表す指標を算出する第２の指標算出手段と、
   前記複数の画像領域の各々について算出された指標と、前記複数の対応領域の各々について算出された指標との類似度を算出する類似度算出手段とを備えた画像処理装置。
3. 前記第２の指標算出手段は、前記複数の画像領域の各々について、前記第２の画像における該画像領域と対応する領域及び該領域の近傍の領域に複数の暫定領域を設定し、該複数の暫定領域の各々について該暫定領域の画像の特徴を表す指標を算出し、該画像領域の指標に最も近い指標が算出された前記暫定領域を該画像領域の対応領域として特定する請求項２に記載の画像処理装置。
4. 前記複数の画像領域の位置と前記複数の対応領域の位置との近似度を算出する近似度算出手段と、
   前記類似度及び前記近似度の大きさに基づいて前記第１の血管分岐と前記第２の血管分岐とが同一の血管分岐を表している蓋然性の程度を示す指標値を算出する蓋然性指標値算出手段とをさらに備えた請求項３に記載の画像処理装置。
5. 前記第１の指標算出手段は、互いに異なる複数の第２の画像の各々について前記類似度を算出し、
   前記第２の指標算出手段は、前記複数の第２の画像の各々について前記近似度を算出し、
   前記蓋然性指標値算出手段は、前記第１の画像と前記複数の第２の画像の各々とからなる複数の組合せについて前記蓋然性指標値を算出し、
   前記複数の組合せの中で最も高い前記蓋然性指標値が算出された組合せを特定する特定手段をさらに備えた請求項４に記載の画像処理装置。
6. 第１のモダリティのボリューム画像から前記第１の画像を生成する第１の生成手段と、
   第２のモダリティのボリューム画像から前記複数の第２の画像を生成する第２の生成手段とをさらに備えた請求項５に記載の画像処理装置。
7. 前記第１の生成手段は、前記第１のモダリティのボリューム画像における血管分岐を含むスラブに対して画素値の最大値投影または最小値投影を行って前記第１の画像を生成し、
   前記第２の生成手段は、前記第２のモダリティのボリューム画像における血管分岐を含むスラブに対して画素値の最大値投影または最小値投影を行って前記第２の画像を生成する請求項６に記載の画像処理装置。
8. 前記第１の生成手段は、前記第１のモダリティの種類に応じて前記最大値投影及び最小値投影のいずれかを行い、
   前記第２の生成手段は、前記第２のモダリティの種類に応じて前記最大値投影及び最小値投影のいずれかを行う請求項６または請求項７に記載の画像処理装置。
9. 前記特定手段により特定された組合せに対応する血管分岐同士が重なるように、前記第１のモダリティのボリューム画像と前記第２のモダリティのボリューム画像との位置合せを行う位置合せ手段をさらに備えた請求項６から請求項８のいずれか一項に記載の画像処理装置。
10. 前記第１のモダリティのボリューム画像から生成された血管分岐を含む複数の画像の中から１つ画像を前記第１の画像として指定する指定手段をさらに備えた請求項６から請求項９のいずれか一項に記載の画像処理装置。
11. 前記指定手段は、操作者の操作に応じて前記第１の画像を指定する請求項１０に記載の画像処理装置。
12. 前記第１のモダリティは、超音波撮像装置であり、
    前記第２のモダリティは、磁気共鳴撮像装置である請求項６から請求項１１のいずれか一項に記載の画像処理装置。
13. 前記第１のモダリティは、超音波撮像装置であり、
    前記第２のモダリティは、放射線断層撮像装置である、請求項６から請求項１１のいずれか一項に記載の画像処理装置。
14. 前記ボリューム画像は、肝臓を含む画像である請求項６から請求項１３のいずれか一項に記載の画像処理装置。
15. 前記第１の画像における画素値に基づいて前記複数の画像領域を設定する設定手段をさらに備えた請求項２から請求項１４のいずれか一項に記載の画像処理装置。
16. 前記設定手段は、前記第１の画像を複数の部分画像に分割し、前記複数の部分画像の各々について該部分画像の特徴度を算出し、該特徴度の大きさに基づいて前記複数の画像領域を設定する請求項１５に記載の画像処理装置。
17. 前記設定手段は、前記特徴度として、前記部分画像内の複数の画素における角らしさ（corner-ness）の程度が反映された値を算出する請求項１６に記載の画像処理装置。
18. 前記設定手段は、前記複数の画像領域を前記第１の画像における前記第１の血管分岐の分岐点の周りに分散して設定する請求項１５から請求項１７のいずれか一項に記載の画像処理装置。
19. 前記第１及び第２の指標算出手段は、画像の特徴をを表す指標として複数の画素における自己相似記述子（self-similarity descriptor）を算出する請求項２から請求項１８のいずれか一項に記載の画像処理装置。
20. コンピュータを、請求項２から請求項１９のいずれか一項に記載の画像処理装置として機能させるためのプログラム。