特許 6799331 画像処理装置及び方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6799331

(24)【登録日】2020年11月25日

(45)【発行日】2020年12月16日

(54)【発明の名称】画像処理装置及び方法

(51)【国際特許分類】

   A61B 8/13 20060101AFI20201207BHJP

【ＦＩ】

   A61B8/13

【請求項の数】24

【全頁数】29

(21)【出願番号】特願2018-531909(P2018-531909)

(86)(22)【出願日】2017年8月1日

(86)【国際出願番号】JP2017027839

(87)【国際公開番号】WO2018025841

(87)【国際公開日】20180208

【審査請求日】2020年4月6日

(31)【優先権主張番号】特願2016-154658(P2016-154658)

(32)【優先日】2016年8月5日

(33)【優先権主張国】JP

【国等の委託研究の成果に係る記載事項】（出願人による申告）平成２６年度 独立行政法人科学技術振興機構、革新的研究開発推進プログラム（ＩｍＰＡＣＴ）「生体データ解析に基づく健康・医療リスク予測モデルの構築」に係る委託研究、産業技術力強化法第１９条の適用を受ける特許出願

(73)【特許権者】

【識別番号】504202472

【氏名又は名称】大学共同利用機関法人情報・システム研究機構

(74)【代理人】

【識別番号】100101454

【弁理士】

【氏名又は名称】山田 卓二

(72)【発明者】

【氏名】備瀬 竜馬

(72)【発明者】

【氏名】鄭 銀強

(72)【発明者】

【氏名】佐藤 いまり

【審査官】 姫島 あや乃

(56)【参考文献】

【文献】 特表２０１２−５３１２７９（ＪＰ，Ａ）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ａ６１Ｂ ８／００−８／１５

(57)【特許請求の範囲】

【請求項1】

所定の撮像対象物を含む所定領域が撮像された複数の画像データ内の撮像対象物を互いに位置合わせするように位置変換する制御手段を備えた画像処理装置であって、
前記制御手段は、前記複数の画像データをそれぞれ、前記撮像対象物を含む前景画像と、前記撮像対象物とは異なる背景画像であって画素値の平均値から離れるにつれて分布数が減少する所定の統計的分布を有する背景画像と、前記撮像対象物の画像上のデータ欠損を補完する補完誤差画像とに分解し、前記前景画像間の類似度を示す第１の評価関数を最小とするように、もしくは前記第１の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第１の位置合わせ処理を実行することを特徴とする画像処理装置。

【請求項2】

前記制御手段は、
（１）前記複数の前景画像の低ランク性を評価する評価関数であって、低ランクほど低い値を取る低ランク評価関数と、
（２）位置合わせ後の前記各背景画像の同一の画素の平均値が、前記各背景画像の平均値から離れるにつれてペナルティを与えるペナルティ関数と、
（３）前記補完誤差画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数と、
の和が最小とするように、もしくは当該和の最小値を含む所定範囲となるように位置合わせを行うことで、前記第１の位置合わせ処理を実行することを特徴とする請求項１記載の画像処理装置。

【請求項3】

前記制御手段は、所定の最適化アルゴリズム法を用いて、少なくとも、前記第１の評価関数が最小とするような位置合わせ、もしくは前記第１の評価関数の最小値を含む所定範囲となるような位置合わせ後の画像データ及び前記位置変換の関数の差分を求めることで、前記第１の位置合わせ処理を実行することを特徴とする請求項１又は２記載の画像処理装置。

【請求項4】

前記制御手段は、前記前景画像と、前記補完誤差画像と、前記背景画像とを出力して可視化することを特徴とする請求項１〜３のうちのいずれか１つに記載の画像処理装置。

【請求項5】

前記制御手段はさらに、前記第１の位置合わせ処理の前に、前記前景画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数が最小とすること、もしくは前記スパース性評価関数が前記スパース性評価関数の最小値を含む所定範囲となることを第２の評価関数とし、前記第２の評価関数が最小とするように、もしくは前記第２の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第２の位置合わせ処理を実行することを特徴とする請求項１〜４のうちのいずれか１つに記載の画像処理装置。

【請求項6】

前記制御手段は、所定の最適化アルゴリズム法を用いて、前記第２の評価関数が最小とするような位置合わせ、もしくは前記第２の評価関数の最小値を含む所定範囲となるような位置合わせ後の画像データ及び前記位置変換の関数の差分を求めることで、前記第２の位置合わせ処理を実行することを特徴とする請求項５記載の画像処理装置。

【請求項7】

前記撮像された複数の画像データは、光超音波画像の画像データであることを特徴とする請求項１〜６のうちのいずれか１つに記載の画像処理装置。

【請求項8】

前記撮像対象物は血管であり、前記背景画像はガウス分布に従うノイズであることを特徴とする請求項１〜７のうちのいずれか１つに記載の画像処理装置。

【請求項9】

前記第１の位置合わせ処理と前記第２の位置合わせ処理のうちの少なくとも一方により位置合わせした前記複数の画像データを平均化することで画像データを作成することを特徴とする請求項１〜８のうちのいずれか１つに記載の画像処理装置。

【請求項10】

前記第１の位置合わせ処理と前記第２の位置合わせ処理のうちの少なくとも一方により位置合わせした画像データのうち、前記各前景画像の画像データを平均化することで画像データを作成することを特徴とする請求項１〜８のうちのいずれか１つに記載の画像処理装置。

【請求項11】

所定の撮像対象物を含む所定領域が撮像された複数の画像データの撮像対象物を互いに位置合わせするように位置変換する制御手段を備えた画像処理装置のための画像処理方法であって、
前記制御手段が、前記複数の画像データをそれぞれ、前記撮像対象物を含む前景画像と、前記撮像対象物とは異なる背景画像であって画素値の平均値から離れるにつれて分布数が減少する所定の統計的分布を有する背景画像と、前記撮像対象物の画像上のデータ欠損を補完する補完誤差画像とに分解し、前記前景画像の類似度を示す第１の評価関数を最小とするように、もしくは前記第１の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第１の位置合わせ処理を実行することを特徴とする画像処理方法。

【請求項12】

前記第１の位置合わせ処理を実行することは、前記制御手段が、
（１）前記複数の前景画像の低ランク性を評価する評価関数であって、低ランクほど低い値を取る低ランク評価関数と、
（２）位置合わせ後の前記各背景画像の同一の画素の平均値が、前記各背景画像の平均値から離れるにつれてペナルティを与えるペナルティ関数と、
（３）前記補完誤差画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数と、
の和が最小とするように、もしくは前記和の最小値を含む所定範囲となるように位置合わせを行うことを含むことを特徴とする請求項１１記載の画像処理方法。

【請求項13】

前記制御手段が、前記第１の位置合わせ処理を実行するステップを実行する前に、前記前景画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数が最小とすること、もしくは当該スパース性評価関数が当該スパース性評価関数の最小値を含む所定範囲となることを第２の評価関数とし、前記第２の評価関数が最小とするように、もしくは前記第２の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第２の位置合わせ処理を実行するステップをさらに含むことを特徴とする請求項１１又は１２記載の画像処理方法。

【請求項14】

前記制御手段が、請求項１１又は１２記載の前記第１の位置合わせ処理を実行するステップを含むことを特徴とする画像処理装置の制御プログラム。

【請求項15】

前記制御手段がさらに、請求項１３記載の前記第２の位置合わせ処理を実行するステップを含むことを特徴とする請求項１４記載の画像処理装置の制御プログラム。

【請求項16】

前記制御手段が、請求項１１又は１２記載の前記第１の位置合わせ処理を実行するステップを実行する前記画像処理装置の制御プログラムを格納したことを特徴とする、コンピュータにより読取可能な記録媒体。

【請求項17】

前記制御手段がさらに、請求項１３記載の前記第２の位置合わせ処理を実行するステップを実行する前記画像処理装置の制御プログラムを格納したことを特徴とする、請求項１６記載の、コンピュータにより読取可能な記録媒体。

【請求項18】

所定の撮像対象物を含む所定領域が撮像された複数の画像データを互いに位置合わせするように位置変換する制御手段を備えた画像処理装置であって、
前記制御手段が、前記撮像対象物を含む前景画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数が最小とするように、もしくは前記スパース性評価関数が前記スパース性評価関数の最小値を含む所定範囲となることを第２の評価関数とし、前記第２の評価関数が最小とするように、もしくは前記第２の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第２の位置合わせ処理を実行することを特徴とする画像処理装置。

【請求項19】

前記制御手段は、所定の最適化アルゴリズム法を用いて、前記第２の評価関数が最小とするように、もしくは前記第２の評価関数の最小値を含む所定範囲となるような位置合わせ後の画像データ及び前記位置変換の関数の差分を求めることで、前記第２の位置合わせ処理を実行することを特徴とする請求項１８記載の画像処理装置。

【請求項20】

前記撮像された複数の画像データは、光超音波画像の画像データであることを特徴とする請求項１８又は１９記載の画像処理装置。

【請求項21】

前記撮像対象物は血管であることを特徴とする請求項１８〜２０のうちのいずれか１つに記載の画像処理装置。

【請求項22】

所定の撮像対象物を含む所定領域が撮像された複数の画像データを互いに位置合わせするように位置変換する制御手段を備えた画像処理装置のための画像処理方法であって、
前記制御手段が、前記撮像対象物を含む前景画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数が最小とすること、もしくは前記スパース性評価関数が前記スパース性評価関数の最小値を含む所定範囲となることを第２の評価関数とし、前記第２の評価関数が最小とすること、もしくは前記第２の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第２の位置合わせ処理を実行するステップを含むことを特徴とする画像処理方法。

【請求項23】

前記制御手段が、請求項２２記載の前記第２の位置合わせ処理を実行するステップを含むことを特徴とする画像処理装置の制御プログラム。

【請求項24】

前記制御手段が、請求項２２記載の前記第２の位置合わせ処理を実行するステップを実行する前記画像処理装置の制御プログラムを格納したことを特徴とする、コンピュータにより読取可能な記録媒体。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、複数の画像データの間で位置合わせを行う画像処理装置及び方法、画像処理装置の制御プログラム並びに当該制御プログラムを格納した記録媒体に関する。

【背景技術】

【0002】

光超音波（Photo-acoustic）画像化（イメージング）法は、癌、腫瘍、血管新生疾患など多くの病気の早期臨床診断のための有望な新技術である。光超音波画像化法は、光音響効果を利用している。すなわち、血管などの検査対象物は、短パルスの近赤外照射を吸収した後、超音波を放射する。検査対象物の３次元構造は熱音響波を検出することによって再構成することができる。この技術は、非侵襲的に任意の造影剤なしで、高空間分解能で生体内の血管を可視化することができる。

【0003】

従来の光超音波画像化法を用いた３次元イメージングにおいては、例えば胸や手などの検査対象物の全体画像を撮影するために、多数の局所領域をスキャンして複数のワンショットボリューム画像データを取得してそれらを合成している（例えば、特許文献１参照）。また、ワンショットボリューム画像データではノイズが強いため、スキャンの際に同じ場所を重複して撮影して、平均化することでノイズを軽減する手法が知られている。

【先行技術文献】

【特許文献】

【0004】

【特許文献1】特開２０１３−１５２４５２号公報

【非特許文献】

【0005】

【非特許文献1】Peng, Y. et al., "RASL: Robust alignment by sparse and low-rank decomposition for linearly correlated images", IEEE Transactions on Pattern Anal Match Intel, Vol. 34, No. 11, pp. 2233-2246, November 2012.

【非特許文献2】Lombaert, H. et al., "Spectral log-demons: Diffeomorphic image registration with very large deformations", International Journal of Computer Vision, Vol. 107, No. 3, pp. 254-271, May 2014.

【非特許文献3】Aylward, S. R. et al., "Registration and Analysis of Vascular Images", International Journal of Computer Vision, Vol. 55, No. 2-3, pp. 123-138, November 2003.

【非特許文献4】Frangi, A. F. et al., "Multiscale vessel enhancement filtering", Proceedings of MICCAI, pp. 130-137, 1998.

【非特許文献5】Lin, Z. et al., "The augmented Lagrange multiplier method for exact recovery of corrupted low-rank matrices", UIUC Technical Report UILU-ENG-09-2215, 2009.

【発明の概要】

【発明が解決しようとする課題】

【0006】

しかしながら、撮影中に体動による位置ずれが生じると、これらの方法では、画質が悪化するという問題があった。上述のように、既存の位置合わせ手法では、ワンショットボリューム画像データに強いノイズが多く含まれるため、局所領域間の類似度を用いた位置合わせ手法では、類似度の算出に対するノイズの影響が強く、位置合わせが難しい。

【0007】

また、ノイズの分離と位置合わせ対象物体行列の低ランク最適化を同時に解く方法（以下、ＲＡＳＬ法という。）が提案されているが（例えば、非特許文献１参照）、ノイズは疎なデータ欠損であるという仮定をもって問題を解いており、光超音波画像化法で取得した画像では仮定と異なり、うまく機能しないという問題点があった。

【0008】

本発明の目的は以上の問題点を解決し、従来技術に比較してノイズを除去することができ、しかも高精度で位置合わせを行うことができる画像処理装置及び方法、画像処理装置の制御プログラム並びにそれを格納した記録媒体を提供することにある。
本発明は、光超音波画像はノイズの影響が大きいものの、ノイズには統計的な特徴があり、このノイズの統計的な特徴を利用することによってノイズの影響を大きく低減できることを見いだしたことによる。

【課題を解決するための手段】

【0009】

第１の発明に係る画像処理装置は、所定の撮像対象物を含む所定領域が撮像された複数の画像データ内の撮像対象物を互いに位置合わせするように位置変換する制御手段を備えた画像処理装置であって、
前記制御手段は、前記複数の画像データをそれぞれ、前記撮像対象物を含む前景画像と、前記撮像対象物とは異なる背景画像であって画素値の平均値から離れるにつれて分布数が減少する所定の統計的分布を有する背景画像と、前記撮像対象物の画像上のデータ欠損を補完する補完誤差画像とに分解し、前記前景画像の類似度を示す第１の評価関数を最小とするように、もしくは前記第１の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第１の位置合わせ処理を実行することを特徴とする。

【0010】

前記画像処理装置において、前記制御手段は、
（１）前記複数の前景画像の低ランク性を評価する評価関数であって、低ランクほど低い値を取る低ランク評価関数と、
（２）位置合わせ後の前記各背景画像の同一の画素の平均値が、前記各背景画像の平均値から離れるにつれてペナルティを与えるペナルティ関数と、
（３）前記補完誤差画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数と、
の和が最小とするように、もしくは当該和の最小値を含む所定範囲となるように位置合わせを行うことで、前記第１の位置合わせ処理を実行することを特徴とする。

【0011】

また、前記画像処理装置において、前記制御手段は、所定の最適化アルゴリズム法を用いて、少なくとも、前記第１の評価関数が最小とするような位置合わせ、もしくは前記第１の評価関数の最小値を含む所定範囲となるような位置合わせ後の画像データ及び前記位置変換の関数の差分を求めることで、前記第１の位置合わせ処理を実行することを特徴とする。

【0012】

ここで、前記画像処理装置において、前記制御手段は、前記前景画像と、前記補完誤差画像と、前記背景画像とを出力して可視化することを特徴とする。

【0013】

さらに、前記画像処理装置において、前記制御手段はさらに、前記第１の位置合わせ処理の前に、前記前景画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数が最小とすること、もしくは前記スパース性評価関数の最小値を含む所定範囲となることを第２の評価関数とし、前記第２の評価関数が最小とすること、もしくは前記第２の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第２の位置合わせ処理を実行することを特徴とする。

【0014】

ここで、前記制御手段は、所定の最適化アルゴリズム法を用いて、前記第２の評価関数が最小とするような位置合わせ、もしくは前記第２の評価関数の最小値を含む所定範囲となるような位置合わせ後の画像データ及び前記位置変換の関数の差分を求めることで、前記第２の位置合わせ処理を実行することを特徴とする。

【0015】

前記画像処理装置において、前記撮像された複数の画像データは、光超音波画像の画像データであることを特徴とする。

【0016】

また、前記画像処理装置において、前記撮像対象物は血管であり、前記背景画像はガウス分布に従うノイズであることを特徴とする。

【0017】

さらに、前記画像処理装置において、前記第１の位置合わせ処理と前記第２の位置合わせ処理のうちの少なくとも一方により位置合わせした前記複数の画像データを平均化することで画像データを作成することを特徴とする。

【0018】

またさらに、前記画像処理装置において、前記第１の位置合わせ処理と前記第２の位置合わせ処理のうちの少なくとも一方により位置合わせした画像データのうち、前記各前景画像の画像データを平均化することで画像データを作成することを特徴とする。

【0019】

第２の発明に係る画像処理方法は、所定の撮像対象物を含む所定領域が撮像された複数の画像データの撮像対象物を互いに位置合わせするように位置変換する制御手段を備えた画像処理装置のための画像処理方法であって、
前記制御手段が、前記複数の画像データをそれぞれ、前記撮像対象物を含む前景画像と、前記撮像対象物とは異なる背景画像であって画素値の平均値から離れるにつれて分布数が減少する所定の統計的分布を有する背景画像と、前記撮像対象物の画像上のデータ欠損を補完する補完誤差画像とに分解し、前記前景画像の類似度を示す第１の評価関数を最小とするように、もしくは前記第１の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第１の位置合わせ処理を実行することを特徴とする。

【0020】

前記画像処理方法において、前記第１の位置合わせ処理を実行することは、前記制御手段が、
（１）前記複数の前景画像の低ランク性を評価する評価関数であって、低ランクほど低い値を取る低ランク評価関数と、
（２）位置合わせ後の前記各背景画像の同一の画素の平均値が、前記各背景画像の平均値から離れるにつれてペナルティを与えるペナルティ関数と、
（３）前記補完誤差画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数と、
の和が最小とするように、もしくは前記和の最小値を含む所定範囲となるように位置合わせを行うことを含むことを特徴とする。

【0021】

また、画像処理方法において、前記制御手段が、前記第１の位置合わせ処理を実行するステップを実行する前に、前記前景画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数が最小とすること、もしくは前記スパース性評価関数の最小値を含む所定範囲となることを第２の評価関数とし、前記第２の評価関数が最小とするように、もしくは前記第２の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第２の位置合わせ処理を実行するステップをさらに含むことを特徴とする。

【0022】

第３の発明に係る画像処理装置の制御プログラムは、前記制御手段が、前記第１の位置合わせ処理を実行するステップを含むことを特徴とする。

【0023】

前記画像処理装置の制御プログラムにおいて、前記制御手段がさらに、前記第２の位置合わせ処理を実行するステップを含むことを特徴とする。

【0024】

第４の発明に係るコンピュータにより読取可能な記録媒体は、前記制御手段が、前記第１の位置合わせ処理を実行するステップを実行する前記画像処理装置の制御プログラムを格納したことを特徴とする。

【0025】

前記コンピュータにより読取可能な記録媒体において、前記制御手段がさらに、前記第２の位置合わせ処理を実行するステップを実行する画像処理装置の制御プログラムを格納したことを特徴とする。

【0026】

第５の発明に係る画像処理装置は、所定の撮像対象物を含む所定領域が撮像された複数の画像データを互いに位置合わせするように位置変換する制御手段を備えた画像処理装置であって、
前記制御手段が、前記撮像対象物を含む前景画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数が最小とすること、もしくは前記スパース性評価関数の最小値を含む所定範囲となることを第２の評価関数とし、前記第２の評価関数が最小とするように、もしくは前記第２の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第２の位置合わせ処理を実行することを特徴とする。

【0027】

前記画像処理装置において、前記制御手段は、所定の最適化アルゴリズム法を用いて、前記第２の評価関数が最小とするような位置合わせ、もしくは前記第２の評価関数の最小値を含む所定範囲となるような位置合わせ後の画像データ及び前記位置変換の関数の差分を求めることで、前記第２の位置合わせ処理を実行することを特徴とする。

【0028】

また、前記画像処理装置において、前記撮像された複数の画像データは、光超音波画像の画像データであることを特徴とする。

【0029】

さらに、前記画像処理装置において、前記撮像対象物は血管であることを特徴とする。

【0030】

第６の発明に係る画像処理方法は、所定の撮像対象物を含む所定領域が撮像された複数の画像データを互いに位置合わせするように位置変換する制御手段を備えた画像処理装置のための画像処理方法であって、
前記制御手段が、前記撮像対象物を含む前景画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数が最小とすること、もしくは前記スパース性評価関数の最小値を含む所定範囲となることを第２の評価関数とし、前記第２の評価関数が最小とするように、もしくは前記第２の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第２の位置合わせ処理を実行するステップを含むことを特徴とする。

【0031】

第７の発明に係る画像処理装置の制御プログラムは、前記制御手段が、前記第２の位置合わせ処理を実行するステップを含むことを特徴とする。

【0032】

第８の発明に係るコンピュータにより読取可能な記録媒体は、前記制御手段が、前記第２の位置合わせ処理を実行するステップを実行する画像処理装置の制御プログラムを格納したことを特徴とする。

【発明の効果】

【0033】

従って、本発明に係る画像処理装置及び方法等によれば、従来技術に比較してノイズを除去することができ、しかも高精度で位置合わせを行うことができる。

【図面の簡単な説明】

【0034】

【図1】本発明の一実施形態に係る画像処理装置１０を含む画像処理システムの構成例を示すブロック図である。

【図2】従来技術に係る光超音波画像化法によって得られた２次元最大強度照射（ＭＩＰ）のワンショットボリューム画像データ（高輝度カーブは血管である）の複数の例を示す画像である。

【図3】（ａ）は図２のワンショットボリューム画像データの背景画像の強度分布を示すグラフであり、（ｂ）は図２の複数のワンショットボリューム画像データの同一画素位置での輝度平均の分布を示すグラフである。

【図4】図１の画像処理装置１０によって実行される位置合わせ処理を示すフローチャートである。

【図5】（ａ）は実画像の平均化画像であり、（ｂ）は位置ずれが発生した合成画像であり、（ｃ）はＢ−ＦＦＤ法により得られた平均化画像であり、（ｄ）はＬｏｇ−ｄｅｍｏｎ法により得られた平均化画像であり、（ｅ）はＲＡＳＬ法により得られた平均化画像であり、（ｆ）は本実施形態にかかる位置合わせ処理により得られた平均化画像である。

【図6】実施形態及び従来法による画像に対する評価結果であって、位置ずれレベルに対する位置ギャップのメトリックを示すグラフである。

【図7A】実施形態に係る画像処理装置１０に位置変換された画像データの画像である。

【図7B】図７Ａの位置変換された画像データから分離された前景画像データＡの画像である。

【図7C】図７Ａの位置変換された画像データから分離された後景画像データＢの画像である。

【図7D】図７Ａの位置変換された画像データから分離された補完画像データＣの画像である。

【図8A】ＭＩＰ法により得られた位置合わせなしの平均化画像データの画像である。

【図8B】本実施形態に係る位置合わせ処理で得られた平均化画像データの画像である。

【図8C】（ａ）及び（ｂ）は図８Ａの拡大図であり、（ｃ）及び（ｄ）はそれぞれ（ａ）及び（ｂ）の画像に対して本実施形態に係る位置合わせ処理を行って得られた画像である。

【発明を実施するための形態】

【0035】

以下、比較例及び本発明に係る実施形態について図面を参照して説明する。なお、以下の実施形態において、同様の構成要素については同一の符号を付している。

【0036】

１．比較例に係る背景技術及び実施形態の概要．
人体の全体の一部分をキャプチャする光超音波画像化法で撮像するために、検査対象物の局所領域をスキャンしてそれらをマージするマルチスキャン及びレジストレーションシステムが開発されてきた。光超音波画像化法によって再構成されたシングルショットボリューム画像データは通常、音と光の散乱及びセンサーレイアウトの制限によって引き起こされる多大なノイズを受けるという問題点があった。

【0037】

本発明の目的は、これらのノイズを有する多数のワンショットボリューム画像データから高品質の３次元画像データを生成することにある。これを達成することで、例えば血管がはっきりと見えるようになる。このようなノイズを低減するために、画像平均化技術がしばしば使用されている。背景画像でのランダムノイズの平均値が小さい定数となる場合には、直線相関を有する前景画像データが明確になる。このことは、画像データのサンプルが静的であるか、又は検査対象物の大きさが動きよりもはるかに大きいような画像平均化処理に限り有効である。

【0038】

しかしながら、局所的にスキャンされた複数のワンショットボリューム画像データをマージすることにより、再構成された光超音波画像化法による画像の場合には、患者は多くの場合、走査中に移動するので前記の仮定はほとんど適用できず、血管が薄い（０．５ｍｍ）であるという事実を考慮すれば、小さな動きが悪画質に影響を与えることは容易に想像できる。従って、そのような身体の動きに対応できる複数のワンショットボリューム画像データを正確に位置合わせする必要がある。

【0039】

低ランクベースの位置合わせ（アライメント）のフレームワークは、非常にデンスな（所定のしきい値以上の高密度な）ノイズや小さな前景画像を含んでいるような挑戦的な画像データの位置合わせを行うことができるという可能性を秘めている。低ランクベースの位置合わせのフレームワークでは、観測行列Ｄは、与えられた複数の画像からまず発生され、ここで、第ｉ番目の行はｉ番目のベクトル画像を示す。もし本発明者らが成功裏によく位置合わせできた画像を作成して位置合わせ画像間の差異を除去できれば、このことはデータ行列Ｄを低ランクにさせることができる。低ランクベースの位置合わせのフレームワークは、

【数1】

に示すように、観測行列Ｄを低ランク成分に変換する変換関数τと誤差成分αを検索することで実現する。前記式の演算子は、画像データに対して位置変換関数を用いて変換を行う演算子である。この最適化問題は基本的に不良設定問題であるため、通常の最適化の誤差成のいくつかの仮定を導入する必要がある。

【0040】

このフレームワークの主たる目的は、誤差成分αを設定することと、光超音波（ＰＡ）画像化法により得られたＰＡ画像データの特徴に基づいて、小さな前景と、非常にデンスなノイズと、破損した前景を含む複数のＰＡボリューム画像データを位置合わせするための最適化方法を設定することにある。この目的を達成するために、発明者らは極小に陥ることなく正確に、小さな血管領域を抽出する粗密（デンスおよびスパース）アプローチを提案する。粗い位置合わせのために、複数のＰＡボリューム画像データに対して、まずＦｒａｎｇｉフィルタ（例えば、非特許文献４参照）を適用することで背景ノイズを減少させ血管の強度を向上させる。次いで、マルチスケールのピラミッド法は、位置合わせされた画像データを低ランクの画像データに変換する位置変換関数を最適化するために使用される。前のステップで位置変換関数の初期推定値が与えられるとき、発明者らはさらに、ＰＡ画像データの特徴に基づいて、入力される画像データを、低ランクの前景画像（血管）と、デンスな（所定のしきい値以上の高密度な）背景画像（多大で厄介なノイズ）と、スパース（所定のしきい値未満の低密度な、まばらな）補完（誤差）画像成分（前景の相補的な部分）とに分解することによりの画像データを位置合わせする。本発明の実施形態における重要な新規性は、ＰＡ画像化の最適化において統計的な特徴（先験的な特徴情報）を有効的に用いていることであり、ここで、各空間ボクセルにおけるデンスなノイズの背景画像の平均値が一定になるように強制される。

【0041】

後述する実施例において、実際のＰＡ画像データを用いた実験結果は、本実施形態に係る提案手法が、非常に緻密なノイズ及び破損を受けることなく複数のボリューム画像データを位置合わせして血管を見出すために有効であることを実証し、クリアな血管構造を持つ高品質の３次元ボリューム画像データを得ることできることを示している。

【0042】

多数の公知方法は、医用画像化装置において複数の画像を位置合わせすることを提案してきた。しかし、ＰＡ画像において観測される血管は非常に密なノイズの背景にわたってまばらに分散している傾向があり、前景画像と背景画像の両方の画像から計算された画像の相似性に依存したこれらの位置合わせ方法を用いることでそのようなＰＡ画像を位置合わせすることを非常に困難にしている。画像の破損の影響を軽減するために、非特許文献１において開示されたロバスト位置合わせ法（ＲＡＳＬ法）が提案されている。このＲＡＳＬ法は医用画像化装置において成功裏にレジストレーションの問題に適用されている。ＲＡＳＬ法は、画像変換を最適化するが、位置合わせすべき複数の画像間の差がまばらな破損に起因していることを前提として、画像を低ランクの前景成分と、まばらな誤差成分とに分離している。ＲＡＳＬ法によって処理された複数の画像を比較すると、血管領域は、ＰＡ画像化の複数の血管画像において深刻なノイズを有して前景画像が小さい。この場合において、誤差成分はもはやスパースにならないので、ＲＡＳＬ法は正しく動作しない。すなわち、画像データを適当に変換することに代えて、前記誤差成分を調整することで低ランク成分が最適化される。

【0043】

２．実施形態に係る画像処理装置．
図１は本発明の一実施形態に係る画像処理装置１０を含む画像処理システムの構成例を示すブロック図である。実施形態に係る画像処理装置１０は、入力画像データ処理部２１と、粗い位置合わせ処理部２２と、細かい位置合わせ処理部２３と、出力画像データ処理部２４と、データ記憶部２５とを備えて構成される。画像処理装置１０には、画像処理装置１０の動作を指示する操作部１１と、撮像デバイス１２と、ディスプレイ１３とが接続される。ここで、画像処理装置１０は例えばディジタル計算機等のコンピュータにより構成される。

【0044】

図１において、撮像デバイス１２は例えば撮像センサとフロントエンド回路とを含んで構成され、例えば、ある局所領域を複数回撮像して、位置合わせすべき複数のボリューム画像データ（例えば２次元又は３次元）を生成して入力画像データ処理部２１に出力する。データ記憶部２５は画像処理装置１０で実行される位置合わせ処理を含む画像処理を実行するために必要なプログラム（図４の画像位置合わせ処理プログラム）及びデータ（例えば、初期位置変換行列τ_１，…，τ_ｍ）を予め格納するとともに、処理途中の画像データ及びパラメータセットを緩衝格納し、出力画像データ及び出力パラメータセットを格納する。なお、画像位置合わせ処理プログラムはデータ記憶部２５に予め格納してもよいし、もしくは、例えばＤＶＤ、ＣＤなどの光ディスクの記録媒体２７に格納された後、画像処理装置１０の動作開始時に光ディスクドライブ装置２６に当該記録媒体２７を挿入して光ディスクドライブ装置２６により当該画像位置合わせプログラムを読み出してデータ記憶部２５にロードして実行してもよい。

【0045】

入力画像データ処理部２１は入力される複数のボリューム画像データを所定の形式のベクトル化された３次元ボリューム画像データに変換して粗い位置合わせ処理部２２に出力する。粗い位置合わせ処理部２２は後述する図４の粗い位置合わせ処理（ステップＳ１〜Ｓ６）を実行することで所定のパラメータセット（Ａ^＊，Δτ^＊）及び位置変換された画像データ

【数2】

を出力する。細かい位置合わせ処理部２３は後述する図４の細かい位置合わせ処理（ステップＳ７〜Ｓ１２）を実行することで所定のパラメータセット（Ａ^＊，Ｂ^＊，Ｃ^＊，Δτ^＊）及び位置変換された画像データ

【数3】

を出力する。出力画像データ処理部２４は位置変換された画像データを表示用の所定の形式の画像データに変換してディスプレイ１３に出力して表示し、もしくは外部の画像データ記憶装置に出力する。

【0046】

以下、位置変換と三成分の分解の最適化推定処理を行う統一されたフレームワークを含む位置合わせ処理について以下詳述する。

【0047】

３．光超音波画像化法により得られた画像データの特徴．
図２は従来技術に係る光超音波画像化法によって得られた２次元最大強度照射（ＭＩＰ）のワンショットボリューム画像データ（高輝度カーブは血管である）の複数の例を示す画像である。ここで、各ボリューム画像データの強度範囲は−１から１までの範囲である。まず、これらの画像で以下の３成分を定義する。
（１）第１の成分は、高強度で血管を表す前景画像データである。
（２）第２の成分は、血管以外の部分にあるノイズを表す背景画像データである。
（３）第３の成分は、前景での破損が発生する複数のパーツを表す補完（誤差）画像データである。ここで、破損はＰＡ画像において、例えば血管の部分が生体内での不均一な光照射によって破壊されて発生する。例えば、図２において、破線の円の中に血管が表示されていない。対応する部分は左の画像に表示されている。この不足している部分は補完（誤差）画像データに相当する。

【0048】

図３（ａ）は図２のワンショットボリューム画像データの背景画像の強度分布を示すグラフであり、図３（ｂ）は図２の複数のワンショットボリューム画像データの同一画素位置での輝度平均の分布を示すグラフである。本発明者らは、実際のＰＡ画像データを調べることにより、図３に示すように、以下の統計的な特徴（先験的な特徴情報）を発見した。
（１）背景画像の強度分布が大きな標準偏差（０．２４９）を持つガウス分布に近いこと。
（２）各空間的なボクセルにおける位置合わせされた背景画像の平均値の分布は、所定の平均値（−０．００１）と、比較的小さな標準偏差（０．０８１）とを有するガウス分布に近いこと。

【0049】

ここで、先験的な特徴情報は、上記（２）の比較的小さな標準偏差を有するガウス分布に近いということを利用して、目的関数のペナルティ関数を設計している。標準偏差が小さいということは、平均値より離れる確率は低くなると考えるので、複数のワンショットボリューム画像データの同一画素位置での輝度平均値が、分布の平均値より離れると大きな値をとるようにしている。

【0050】

これとは対照的に、血管の平均強度は、血管が良好に位置合わせされたショットボリューム画像データの間で線形に相関するので、高い値を有する。特に、この統計的な特徴（先験的な特徴情報）はしばしば、本実施形態に係る画像処理装置１０等の画像化システムのみならず、従来暗黙のうちに使用する、ノイズを広く軽減するために従来技術に係る画像平均化技術が使用される多くのアプリケーションにおいて、高い確率で真実である。バックグラウンドにおいては、輝度値はガウス分布に従うので、無相関だが、血管の輝度強度は、バックグラウンドと比べ、高い値を取る。そのため、位置合わせに成功している場合は、血管上の同じ場所のボクセルにおいて、輝度値は、類似しており（相関があり）、高い値を取る。言い換えると、位置合わせがうまくいっていない場合は、ある画像ではあるボクセル値がバックグラウンドであるのに、ある画像ではそれに対応するボクセル値が血管領域で、輝度値が高いということがおき、この場合、そのボクセルの輝度平均値も高い値に引っ張られ、ペナルティ値が高くなる。また、ノイズの軽減に従来よく利用される従来技術に係る画像平均化手法は、位置合わせができていることが想定されている場合、前景領域の同じ場所は同等の輝度値を持っており（相関があり）、ノイズは無相関であることを暗黙に利用しており、同じ場所のボクセルが前景領域であると、複数の画像間で同様の値を持つので、その平均値は、その場所の本来の輝度値に近づき、ノイズ領域のボクセルは、無相関なので、その平均が小さい値になり、その結果、前景領域とノイズ領域のＳＮ比がよくなる。

【0051】

本実施形態では、本発明者らは、３つの成分の分解推定するために、これらの統計的な特徴（先験的な情報）を有効的に用いる。

【0052】

４．本実施形態に係る、粗い位置合わせから細かい位置合わせ（より高精度な位置合わせ）を行う低ランクベースの位置合わせ処理
緻密なノイズと小さな前景画像領域を含むＰＡ画像の場合には、現在の低ランクに基づく方法は、誤差成分を調整することにより、低ランクの成分を最適化する傾向がある。この問題を回避するために、発明者らは粗から密へのアプローチを用いる。

【0053】

粗い位置合わせのために、発明者らは、誤差成分を考慮せずに、ノイズ低減のデータを使用して変換を最適化するためのマルチスケールピラミッド法を用いる。これは、局所的な最小点に陥ることなく、適切な変換を見つけるのに役立つ。さらに高精度な位置合わせ画像を得るためには、本発明者らは、ＰＡ画像化の統計的特徴に基づいて、位置合わせすべき画像データを、低ランクの前景画像、高密度のノイズの背景画像、及びスパースな補完（誤差）画像データの３個の成分に適切に分解しながら、変換関数を最適化する。この手順の詳細は、次のセクションで説明する。

【0054】

４．１．ピラミッドを用いた低ランク最小化による粗い位置合わせ処理．
いまスキャン中に移動した可能性があるあるサンプルである、複数ｍ個のボリューム画像データ

【数4】

が与えていると仮定する。ここで、ｎ次元の空間Ｒ^ｎは例えば、幅、高さ、奥行きの３次元の空間である。このステップは概ね、前処理された前景画像データを使用して、これらのボリューム画像データを整列させて得られる。ここで、背景ノイズは、Ｆｒａｎｇｉフィルタ（例えば、非特許文献４参照）を通過させ、血管を増強することによって低減される。ここで、管状構造らしさを計算するために、ヘッセ行列の固有ベクトルを使用している。フィルタリングされた画像データは背景において偽陽性のデータを含み、血管が破損している。しかし、その画像データは粗い位置合わせのためには十分である。代わりに、誤差成分を調整する適切な変換関数を見つけるために、差分αを考慮せずに、前処理されたデータ行列Ｄをただ単に、低ランク成分の画像データ

【数5】

に変換する変換関数を最適化する。ここで、

【数6】

は観測データ行列であり、ｎは各ショットボリューム画像データ内のボクセル数である。また、τは例えばアフィン変換などの特定の位置線形変換（以下、位置変換という）関数を表す。例えば非特許文献１によれば、位置変換関数τの変化が小さいとき、

【数7】

は位置変換関数τについて線形化することで近似することができ、ランク最小化は次式（１）に示すように核ノルム

【数8】

を最小化するように緩和される。

【0055】

【数9】

（１）

【0056】

ここで、Δτは各反復において位置変換関数τの分散を表す。また「ｓ．ｔ．」は「ｓｕｃｈ ｔｈａｔ」の略であって「制約条件」を意味し、右側の条件式を満たすように、左側の式における前記核ノルムを最小化するときのパラメータＡ，τを求めることを意味する。以下、「ｓ．ｔ．」については以下同様である。また、左側の式の演算子はノルムを演算する演算子である。さらに、Ｊ_ｉは次式で定義されるｉ番目のデータの
アン行列である。

【0057】

【数10】

【0058】

そして、ε_ｉはｍ次元空間

【数11】

のための標準的な基底を表す。式（１）については、非特許文献１のＲＡＳＬ法と同様に、公知の拡張ラグランジュ未定乗数法（ＡＬＭ）アルゴリズムを用いて解くことができる（例えば、非特許文献６参照）。この最適化は、順次より細かいガウシアンピラミッドを有する画像データを成功裏に位置合わせするマルチスケールピラミッド法に適用される。

【0059】

４．２．三成分の分解を用いる細かい位置合わせ処理（例えば粗い位置合わせ処理の後に行う）．
前記前処理された画像データは、前のステップでの最適化はノイズを考慮していないので小さな位置合わせエラーにつながる、偽陽性と偽陰性のノイズを含むかもしれない。本実施形態では、ノイズを考慮して、元のボリュームを使用することにより粗い位置合わせの結果の画像データをさらに高精度な位置合わせ処理を行う。本実施形態に係る細かい位置合わせ処理では、密なノイズが含まれている画像を処理するために、位置合わせ処理で緻密なノイズ項を導入し、緻密なノイズの多い背景画像の特徴を用いて、低ランクの前景データＡ、高密度のノイズ背景データＢ、及びスパースな補完データＣに分解して、データの位置を調整することで位置合わせを行う。具体的には、同じ位置において位置合わせされた密なノイズの平均値の分布は、所定値の平均と小さな分散を有するガウス分布となる。本実施形態では、この統計的な特徴を次式のように、分解問題の定式化に導入する。

【0060】

【数12】

（２）

【0061】

ここで、演算子ｍｉｎの対象物は左の式の３つの項を含み、演算子ｍｉｎは左の式の３つの項の和が実質的に最小となるような前景画像データ行列Ａ及び位置変換行列τを求める。ここで、当該和がその最小値を含む所定範囲となるようにしてもよい。Ａ及びＢはそれぞれ、位置合わせされた、前景画像データ行列、背景画像データ行列を表し、Ｃは補完行列を表す。行列Ａ、Ｂ、及びＣの和は位置合わせされた画像データ

【数13】

に等しい。また、式（２）の「ベクトルａ」は、すべての要素がすべて１である１×ｍのベクトルである。さらに、式（２）の「ベクトルｂ」は、すべての要素がガウス分布の平均値であるすべて一定値（スカラー）ｂである１×ｎのベクトルである。ここで、実際のＰＡ画像データの、発明者らの観察に基づいてスカラーｂに０を設定する。式（２）の第２項

【数14】

は、同じ位置において位置合わせされた背景画像データの平均値がガウス分布の平均値から離れてとどまるときに、最適化関数にペナルティを課する項である。また、重み付けパラメータλ_１に次式を設定し、

【数15】

また、重み付けパラメータλ_２に次式を設定する。

【0062】

【数16】

【0063】

ここで、発明者らの実験において、κ_１＝１、κ_２＝３に設定した。

【0064】

４．３．全体の位置合わせ処理
図４は図１の画像処理装置１０によって実行される全体の位置合わせ処理を示すフローチャートである。全体の位置合わせ処理は、粗い位置合わせ処理と、細かい位置合わせ処理とから構成されるが、本発明はこれに限らず、細かい位置合わせ処理のみを実行してもよい。ここで、画像間の位置ずれが比較的大きな場合は、粗い位置合わせ処理が必要であるが、画像間の位置ずれが比較的小さい場合は、粗い位置合わせ処理を実行することなく、細かい位置合わせ処理を実行してもよい。

【0065】

図４において、ステップＳ１〜Ｓ６までの処理は、図１の粗い位置合わせ処理部２２によって実行される粗い位置合わせ処理であり、ステップＳ７〜Ｓ１２までの処理は、図１の細かい位置合わせ処理部２３によって実行される細かい位置合わせ処理である。

【0066】

図４のステップＳ１においてベクトル化された３次元ボリュームデータ

【数17】

及び初期線形変換行列τ_１，…，τ_ｍに基づいて、次式（３）を用いて
アン行列Ｊ_ｉを計算する。

【0067】

【数18】

（３）

【0068】

ここで、初期線形変換行列τ_ｉは次式で表される。

【0069】

【数19】

【0070】

次いで、ステップＳ２において、次式（４）を用いて画像データの変換及び正規化を行って線形変換されたデータ

【数20】

を計算する。

【0071】

【数21】

（４）

【0072】

次いで、ステップＳ３において、次式（５）を用いて、ＡＬＭアルゴリズムを用いた線形化された凸最適化を行ってパラメータセット（Ａ^＊，Δτ^＊，Ｙ^＊）を計算する。

【0073】

【数22】

（５）

【0074】

式（５）の右辺は評価関数である拡張ラグランジュ関数であり、その第１項は、前景のデータ行列Ａのノルムを最小化するとき（すなわち、位置合わせされたとき）の前景画像データを示す。第２項及び第３項はＡＬＭアルゴリズムを用いて線形化された凸最適化を行うための追加項である。ここで、Ｙは詳細後述するように、ＡＬＭアルゴリズムにおけるラグランジュ乗数行列である。また、右辺第２項の演算子＜Ｐ，Ｑ＞は行列Ｐ，Ｑ間の内積を演算する演算子であり、右辺第２項は前景画像データ行列Ａと位置変換後のラグランジュ乗数行列Ｙ^＊との内積を表す。さらに、右辺第３項は位置変換された画像データと前景画像データとの間の位置ずれの誤差の距離に対応する。式（５）の右辺では、演算子ｍｉｎの対象物は３つの項を含み、演算子ｍｉｎは３つの項の和を最小化するときの前景画像データ行列Ａ及び変換行列の差分Δτを求める。

【0075】

次いで、ステップＳ４において位置変換関数τを位置変換関数の差分Δτを用いて次式のように更新する。

【0076】

【数23】

【0077】

そして、ステップＳ５において、次式（６）の第１の終了条件（粗い位置合わせ誤差が所定値未満である）を満たすか否かが判断される。なお、ｔｏｌは終了条件判断用しきい値許容値（偏差）である。

【0078】

【数24】

（６）

【0079】

ステップＳ５において、ＹＥＳのときはステップＳ６に進む一方、ＮＯのときはステップＳ１に戻る。ステップＳ６ではパラメータセット（Ａ^＊，Δτ^＊）をデータ記憶部２５に出力して一時的に格納する。ここで、Ａ^＊は位置合わせ後の画像データであり、以上で粗い位置合わせ処理を終了する。

【0080】

次いで、以下の細かい位置合わせ処理を実行する。

【0081】

ステップＳ７において、ベクトル化された３次元ボリュームデータ

【数25】

及び初期線形変換行列τ_１，…，τ_ｍに基づいて、次式（７）を用いて
アン行列Ｊ_ｉを計算する。

【0082】

【数26】

（７）

【0083】

次いで、ステップＳ８において次式（８）を用いて画像データの位置変換及び正規化を行って線形変換されたデータ

【数27】

を計算する。

【0084】

【数28】

（８）

【0085】

次いで、ステップＳ９において、次式（９）を用いて、ＡＬＭアルゴリズムを用いた線形化された凸最適化を行ってパラメータセット（Ａ^＊，Ｂ^＊，Ｃ^＊，Δτ^＊，Ｙ^＊）を計算する。

【0086】

【数29】

（９）

【0087】

式（９）の右辺は評価関数である拡張ラグランジュ関数であり、その第１項は、前景のデータ行列Ａのノルムを最小化するとき（すなわち、位置合わせされたとき）の前景画像データを示す。第２項は背景画像データ行列Ｂが平均値０から離れたときにペナルティ値を付加する項である。第３項はスパース誤差画像データ行列Ｃの誤差に対応する。ここで、位置合わせ後の前景画像同士が類似していることを評価基準として低ランクを導入し、当該第３項の誤差が少なくなるようにＬ−１ノルムを導入して最適化問題を構成している。さらに、第４項及び第５項はＡＬＭアルゴリズムを用いて線形化された凸最適化を行うための追加項である。式（９）の右辺では、演算子ｍｉｎの対象物は５つの項を含み、演算子ｍｉｎは５つの項の和を最小化するときの前景画像データ行列Ａ、背景画像データ行列Ｂ、スパースな誤差画像データ行列Ｃ、及び変換行列の差分Δτを求める。

【0088】

次いで、ステップＳ１０において位置変換関数τを位置変換関数の差分Δτを用いて次式のように更新する。ここでの繰り返し処理により、評価関数を実質的に最小化し、もしくは評価関数の最小値を含む所定の範囲（例えば＋／−０％から〜＋／−５％ないし＋／−１０％）となるように行うことができる。

【0089】

【数30】

【0090】

次いで、ステップＳ１１において、式（１０）の第２の終了条件（細かい位置合わせ誤差が所定値未満である）を満たすか否かが判断される。なお、ｔｏｌは終了条件判断用しきい値許容値（偏差）である。

【0091】

【数31】

（１０）

【0092】

ステップＳ１１において、ＹＥＳのときはステップＳ１２に進む一方、ＮＯのときはステップＳ７に戻る。ステップＳ１２では、パラメータセット（Ａ^＊，Ｂ^＊，Ｃ^＊，Δτ^＊）を出力画像データ処理部２４に出力する。ここで、Ａ^＊は位置合わせ後の画像データであり、以上で細かい位置合わせ処理を終了し、当該画像位置合わせ処理を終了する。

【0093】

図４の２つの各反復では、位置変換関数τの現在の推定値で変換を線形化することにより、線形化及び最適化問題を解決した。この最適化問題は上述のＡＬＭアルゴリズムを用いて、線形化された凸最適化として解ける。

【0094】

本実施形態においては、次式の拡張増補ラグランジュ関数を最小化する最適化問題に帰着させて解を求めている。

【0095】

【数32】

（１１）

【0096】

ここで、

【数33】

はラグランジュ乗数行列であり、μは正のスカラー計数である。この問題は次式のように反復して最小化される。

【0097】

【数34】

（１２）

【0098】

ここで、ｓｖｄ（・）は特異値分解の演算子を表し、縮小演算子

【数35】

は核ノルムを最適化するための特異値閾値化アルゴリズムのためのものである。μ_ｋは単調に増加する正のシーケンスであり、本実施形態では、非特許文献１に従って、次式のように設定した。

【0099】

【数36】

【0100】

なお、提案手法の計算時間は、他の方法と同等であった。

【実施例】

【0101】

５．実験評価．
まず、本発明者らは実際の３次元ＰＡ画像データであるショットボリューム画像データに基づいて、本実施形態に係る提案された画像位置合わせ処理のロバスト性と精度を評価した。本発明者らは、スキャン中に被験者の身体を固定することにより、ピクセルレベルのずれが存在しないように、真実の画像データを５セット作成した。次いで、合成された真実のショットボリューム画像データに位置ギャップの異なるレベルを追加し、元のボリューム画像データをランダムに各軸（ｘ，ｙ，ｚ）に沿ったボクセルが均一に分布した擬似乱数［−ｄ ｄ］で置換した。ボクセルサイズは、［０．２５ ０．２５ ０．２５］（ｍｍ）で血管の厚さは０．５〜２ｍｍであった。ここで、発明者らは、１から１０のずれ量ｄを設定した。全体では、発明者らは１０×５個のデータセットを生成し、ここで、各データセットは、２１から３６個のショットボリューム画像データを含む。

【0102】

比較のために、発明者らは、以下の従来技術に係る方法により得られたものと発明者らの結果を比較した。
（１）類似機能を最小限に抑えながら、制御点を最適化するＢスプライン自由形状変形法（以下、Ｂ−ＦＦＤ法という）（例えば、非特許文献３参照）：。
（２）単純な最近傍検索で画像間の地点別の対応を発見するスペクトルＬｏｇ−ｄｅｍｏｎレジストレーション（以下、Ｌｏｇ−ｄｅｍｏｎ法という）（例えば、非特許文献４参照）；
（３）公平な比較のために血管拡張データに適用した、ピラミッド法を備えたＲＡＳＬ法（例えば、非特許文献１参照）。

【0103】

図６は実施形態及び従来法による画像に対する評価結果であって、位置ずれレベルに対する位置ギャップのメトリックを示すグラフである。図６において、縦軸の値が小さい場合は、各方法が正常にデータセットを位置合わせできたことを示す。これにより、本実施形態に係る提案手法はミスアライメントレベル７までうまく機能していることを示しており、このことは本当の身体の動きに対処するには十分であるというべきである。他の方法は、非常に低い位置ずれレベルで、良好に画像データの位置合わせを行うことができなかった。

【0104】

図５（ａ）は実画像の平均化画像であり、図５（ｂ）は位置ずれが発生した合成画像であり、図５（ｃ）はＢ−ＦＦＤ法により得られた平均化画像であり、図５（ｄ）はＬｏｇ−ｄｅｍｏｎ法により得られた平均化画像であり、図５（ｅ）はＲＡＳＬ法により得られた平均化画像であり、図５（ｆ）は本実施形態にかかる位置合わせ処理により得られた平均化画像である。すなわち、図５は、位置ずれレベル６で、それぞれの方法で得られた位置合わせ結果の平均画像の例を示している。

【0105】

図５（ａ）は真実の画像データを示し、ここで血管がはっきりと見える。図５（ｂ）から明らかなように、合成ずれが画質を大幅に悪化させることを示している。図５（ｆ）の提案方法の結果は、明確な血管を示しており、画像が真実の画像データ（図５（ａ））に類似している。このことは、提案手法が適切にショットボリューム画像データを位置合わせできることを示している。対照的に、血管は、さらに他の方法の結果で識別することが困難であった（図５（ｃ）、図５（ｄ）、図５（ｅ））。

【0106】

以下、本発明者らの提案方法は、実際のＰＡデータを３つの成分にどのように分解したかについて精査する。

【0107】

図７Ａは実施形態に係る画像処理装置１０により位置変換された画像データ

【数37】

の画像である。図７Ｂは図７Ａの位置変換された画像データ

【数38】

から分離された前景画像データＡの画像である。図７Ｃは図７Ａの位置変換された画像データ

【数39】

から分離された後景画像データＢの画像である。図７Ｄは図７Ａの位置変換された画像データ

【数40】

から分離された補完画像データＣの画像である。

【0108】

すなわち、図７Ａ〜図７Ｄはミスアラインメントレベル６での結果を示している。背景雑音は正しく後景画像データＢに分解されている。すべての血管は明らかに、低ランクの前景画像成分Ａに表示され、そして、血管のまばらな破損は、補完成分Ｃに表示されている。この結果は、提案方法が、元のノイズの多いＰＡ画像が、正常に３個の意味のある成分分解できたことを示している。

【0109】

次に、本発明者らは実際の体動に対する提案手法のロバスト性を精査した。実際のデータを取得するには、手は、螺旋パターンで局所領域を複数回走査する広視野範囲のＰＡ画像化システムを使用してスキャンした。そして、実際の身体の動きは、スキャン中に追加された。ショットボリューム画像データの総数は２０４８だった。総スキャン時間は約２分だった。連続したローカルショットボリューム画像データは８５％にわたって重畳された。すべての画像データは５１２×５１２×１００画素のデータであった。

【0110】

各ショットボリューム画像データは、局所エリアデータの一部であるので、発明者らはまず、１６個から５９個のショットボリューム画像データを含むローカルエリア（１６８×１６８×１００画素）のデータセットを生成した。次いで、発明者らは、各データセットにおいてショットボリューム画像データの位置合わせを行った。その後、全てのデータを生成するために位置合わせされた平均ボリューム画像データをともに縫い合わせた。比較のために、位置合わせなしの画像平均化技術が画像データデータに適用された。

【0111】

図８ＡはＭＩＰ法により得られた位置合わせなしの平均化画像データの画像である。また、図８Ｂは本実施形態に係る位置合わせ処理で得られた平均化画像データの画像である。さらに、図８Ｃ（ａ）及び（ｂ）は図８Ａの拡大図であり、図８Ｃ（ｃ）及び（ｄ）はそれぞれ図８Ｃ（ａ）及び（ｂ）の画像に対して本実施形態に係る位置合わせ処理を行って得られた画像である。

【0112】

図８Ａ、図８Ｃ（ａ）及び（ｂ）は位置合わせしない平均画像データを示す。これらの画像では、血管のコントラストが低く、いくつかの血管が体動によりぼやけている。これに対して、図８Ｂ、図８Ｃ（ｃ）及び（ｄ）は提案方法によるレジストレーション結果を示す。これらの画像では、血管は明らかに元の画像よりも鮮明であり、位置合わせしないときは平均画像データを見ることが困難であったいくつかの血管が見えるようになった。ここで、図８Ｂにおけるボックスに三回現れる血管について注意すべきであり、２本の静脈は動脈と並行して実行する場合にぼやけていない。これらの画像の特徴は、解剖学の医師によって確認され、その結果は、発明者らの提案方法が大幅にＰＡ画像化の画質を大幅に改善したことを示している。

【0113】

６．結論
発明者らは、ＰＡ画像の特徴を検討し、高品質の３次元ボリューム画像データを生成するために、ＰＡ画像の画像化のためのレジストレーション方法を提案した。提案方法を用いることで血管がはっきりと見えるようになった。位置合わせされた背景データの同じ位置の値の平均値が一定になる傾向を有するという統計的特徴を導入することで、発明者らの位置合わせ方法は、強力なノイズと大きなずれとを有する挑戦的なＰＡ画像データを処理することが可能である。実際のデータセットについての実験結果は、提案手法の有効性を実証しており、著しく画質が向上し、従来技術に比較して最良の位置合わせ精度を達成した。

【0114】

７．非特許文献１のＲＡＳＬ法と本実施形態の方法との相違点について（まとめ）．
以下の説明では、ＡＬＭアルゴリズムに係る項を省略して、ＲＡＳＬ法と本実施形態の方法の概要について以下に説明する。

【0115】

７−１．ＲＡＳＬ法
ＲＡＳＬ法はＲｏｂｕｓｔ−ＰＣＡを応用した手法であって、入力画像である位置合わせ前の入力画像データＤから、前景画像Ａと、ノイズ画像Ｅと、位置合わせのための位置変換τとを推定する。入力画像データＤを位置変換した位置合わせ後のデータを

【数41】

として、位置合わせ後の前景画像データ同士が類似していることの評価基準として、低ランクを導入し、差分誤差はなるべく少なくなるようにＬ−１ノルムを導入することで、最適化問題の評価関数を次式のように設定している。

【0116】

【数42】

【0117】

ここで、第１項は位置合わせ後の前景画像データであり、第２項はスパース（まばらな）誤差である。

【0118】

制約条件

【数43】

が非線形なので、位置変換関数である変数τのローカルな傾きとして、Δτに関して線形で表現することで問題を緩和し、解けるようにしている。ここで、変数τは繰り返し最適化の中で固定して、次式のごとくΔτを最適化の目的変数としていることを特徴としている。

【数44】

【0119】

【数45】

【0120】

７−２．３次元への適応におけるＲＡＳＬ法の問題点．
（問題点１）以下の非凸問題がある。
ＲＡＳＬ法では、位置ずれは小さいと仮定しており、３次元でかつ前景画像（血管）が小さいと、局所解に陥りやすいという問題点がった。具体的には、ＲＡＳＬ法をそのまま単純に３次元に適用すると正しく動作しないという問題点があった。
（問題点２）誤差画像領域はスパース誤差であると仮定しており、ＰＡ画像では、ノイズ領域はきわめて大きいという問題点があった。
７−３．本実施形態にかかる提案手法．
７−３−１．粗い位置合わせ処理．
ＲＡＳＬ法の第１の問題点を解決するために、位置合わせ処理により誤差成分で低ランクを実現することで、次式のごとく誤差成分なしで位置合わせの最適化のみで低ランクを検索している。

【0121】

【数46】

【0122】

前記の式では、位置変換してなる、位置合わせ後の前景画像同士が類似しているという探索基準で低ランクを導入し、前景画像Ａの核ノルムが最小となるように前景画像Ａ及び位置変換関数を探索する。
７−３−２．細かい位置合わせ処理．
ＲＡＳＬ法では、ノイズ領域は小さいと仮定しており、ＰＡ画像と異なるという問題点があった。本実施形態では、次式のごとく、前景画像Ａ、スパースな誤差画像Ｅに加え、背景画像Ｂを導入したことを特徴としている。背景画像Ｂはガウス分布などの統計的特徴（位置あわせ後の同じ座標の値の平均値は０に近づく）に従い、無相関である特徴を利用している。

【数47】

【0123】

前記の式では、第１項は前景画像Ａの核ノルムであり、第２項は背景画像Ｂがガウス分布に従って統計的に無相関であり、平均値が０から離れるとペナルティ値を付与する項であり、第３項はスパース誤差画像Ｅ（実施形態における補完画像Ｃと同義である）のＬ１ノルムである。

【0124】

前記細かい位置合わせ処理では、ノイズを強い分散を持つ高密度なガウス分布のノイズと、疎なデータ欠損との２つからなると仮定した問題を設定し、ガウス分布から得られるノイズの平均の分布は分散が弱いガウス分布に従うという統計上の特性を用いて、その問題を解く最適化手法を提案することで、ＰＡ画像の実態にあった仮定で位置合わせ問題を解くことができ、従来技術に比較して前景画像についての位置合わせ処理において良好な結果が得られる。

【0125】

８．変形例．
以上の実施形態において、図４に示すように、粗い位置合わせ処理を実行した後、細かい位置合わせ処理を実行しているが、本発明はこれに限らず、粗い位置合わせ処理を実行せずに、細かい位置合わせ処理のみを実行してもよい。このことは入力される複数のボリューム画像データ間での差分が比較的小さい場合に有効であり、処理時間の短縮を行うことができる。

【0126】

以上の実施形態において、光超音波画像（ＰＡ画像）に対する位置合わせ処理について説明しているが、本発明はこれに限らず、所定の対象物又は対象場所の複数の撮像画像に対して位置合わせ処理を行ってもよい。

【0127】

以上の実施形態において、前景画像として血管の一例を占めているが、本発明はこれに限らず、所定のしきい値以上の強度（画素値）を有する対象物の画像であってもよい。

【0128】

以上の実施形態において、位置変換関数として、アフィン変換を用いているが、本発明はこれに限らず、例えば平行移動変換、回転変換、拡大又は縮小変換、鏡映変換、射影変換などのグルーバルドメインでのパラメトリック線形変換関数であればよい。

【0129】

以上の実施形態において、ＰＡ画像の統計的特徴を有する分布としてガウス分布を用いているが、本発明はこれに限らず、以下の統計的分布を用いてもよい。
（１）ノイズの分布：位置合わせ対象画像が多い場合は、基本的にどんな分布でもよい。例えば一様分布、対数正規分布、カイ２乗分布、ワイブル分布、ガンマ分布などである。ただし、位置合わせ対象画像が少ない場合は、中心（平均）から離れるに従って、減少していく分布であり、例えば正規分布、ｔ分布、双曲線正割分布、ローレンツ分布、ラプラス分布などを用いることができる。
（２）平均化後の分布：中心（平均）から離れるに従って、減少していく分布ならばよく、例えば正規分布、ｔ分布、双曲線正割分布、ローレンツ分布、ラプラス分布などの分布を用いることができる。

【0130】

以上の実施形態においては、上述の評価関数を拡張ラグランジュ未定乗数法のＡＬＭアルゴリズムを用いて線形化された凸最適化を行って位置変換後の画像データ及び位置変換関数の差分などを求めているが、本発明はこれに限らず、拡張ラグランジュ未定乗数法に代えて所定の多変数探索法などの最適化アルゴリズム法を用いて位置変換後の画像データ及び位置変換関数の差分などを求めてもよい。例えば、以下の最適化アルゴリズム法を用いてもよい。
（１）ａｃｃｅｌｅｒａｔｅｄ ｐｒｏｘｉｍａｌ ｇｒａｄｉｅｎｔ（ＡＰＧ） ａｌｇｏｒｉｔｈｍ；
（２）Ｒｏｂｕｓｔ−ＰＣＡで用いられているＤｕａｌ Ｍｅｔｈｏｄ；
（３）Ｓｉｎｇｕｌａｒ Ｖａｌｕｅ Ｔｈｒｅｓｈｏｌｄｉｎｇ。

【0131】

以上の実施形態においては、図４の細かい位置合わせ処理において、パラメータセット（Ａ^＊，Ｂ^＊，Ｃ^＊，Δτ^＊）を出力しているが、本発明はこれに限らず、少なくともパラメータセット（Ａ^＊，Δτ^＊）を出力してもよい。

【0132】

以上の実施形態においては、低ランク化の方法として「核ノルム」を最小化しているが、本発明はこれに限らず、「核ノルム」に代えて、「行列のランクを評価し、低ランクほど小さい値を取る評価関数」を用いてもよい。

【0133】

以上の実施形態においては、補完誤差画像がスパースであることを評価する方法として、「Ｌ１ノルム」を利用しているが、本発明はこれに限らず、スパース性（行列やベクトルにおいて、非０の要素が少ないことをいう）の指標として、例えば「Ｌ０ノルム」などの「補完誤差画像のスパース性を評価し、スパースであるほど小さな値を取る評価関数」を用いてもよい。

【0134】

なお、撮像された複数の画像データを用いて画像を作成する画像処理装置において最終結果の画像データを得るために、前記位置合わせ処理により位置合わせした前記複数の画像データ（撮像された複数のオリジナル画像データ）を平均化することで、最終結果の画像データを作成してもよいし、もしくは、前記位置合わせ処理により位置合わせした画像のうち、前記複数の前景画像データを平均化することで、最終結果の画像データを作成してもよい。

【0135】

なお、画像処理装置１０は、前記分解された、前記前景画像と、前記補完誤差画像と、前記背景画像とを外部装置に出力して可視化してもよい。

【産業上の利用可能性】

【0136】

以上詳述したように、本発明に係る画像処理装置及び方法等によれば、従来技術に比較してノイズを除去することができ、しかも高精度で位置合わせを行うことができる。

【0137】

本発明に係る画像処理装置及び方法等は以下のために広く用いることができる。
（１）医療画像化装置：非侵襲で以下の血管系異常を、２次元又は３次元可視化することができ、早期診断、健康診断を行うことができる。例えば、皮膚毛細血管、がんの新生血管、関節の新生血管、動脈狭窄、閉塞、プラーク、脳酸素消費、血管走行異常などである。本実施形態では、特に、血管の自動抽出装置に有用である。
（２）製品検査装置：食品又は工業製品などの劣化、傷などの検査を行うことができる。

【符号の説明】

【0138】

１０…画像処理装置、
１１…操作部、
１２…撮像デバイス、
１３…ディスプレイ、
２１…入力画像データ処理部、
２２…粗い位置合わせ処理部、
２３…細かい位置合わせ処理部、
２４…出力画像データ処理部、
２５…データ記憶部、
２６…光ディスクドライブ装置、
２７…記録媒体。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7A】

【図7B】

【図7C】

【図7D】

【図8A】

【図8B】

【図8C】