JPWO2018025841A1 - 画像処理装置及び方法 - Google Patents

画像処理装置及び方法 Download PDF

Info

Publication number
JPWO2018025841A1
JPWO2018025841A1 JP2018531909A JP2018531909A JPWO2018025841A1 JP WO2018025841 A1 JPWO2018025841 A1 JP WO2018025841A1 JP 2018531909 A JP2018531909 A JP 2018531909A JP 2018531909 A JP2018531909 A JP 2018531909A JP WO2018025841 A1 JPWO2018025841 A1 JP WO2018025841A1
Authority
JP
Japan
Prior art keywords
image
evaluation function
image data
processing apparatus
alignment
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
JP2018531909A
Other languages
English (en)
Other versions
JP6799331B2 (ja
Inventor
竜馬 備瀬
竜馬 備瀬
銀強 鄭
銀強 鄭
いまり 佐藤
いまり 佐藤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Inter University Research Institute Corp Research Organization of Information and Systems
Original Assignee
Inter University Research Institute Corp Research Organization of Information and Systems
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Inter University Research Institute Corp Research Organization of Information and Systems filed Critical Inter University Research Institute Corp Research Organization of Information and Systems
Publication of JPWO2018025841A1 publication Critical patent/JPWO2018025841A1/ja
Application granted granted Critical
Publication of JP6799331B2 publication Critical patent/JP6799331B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/13Tomography
    • GPHYSICS
    • G02OPTICS
    • G02BOPTICAL ELEMENTS, SYSTEMS OR APPARATUS
    • G02B21/00Microscopes
    • G02B21/36Microscopes arranged for photographic purposes or projection purposes or digital imaging or video purposes including associated control and data processing arrangements
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00General purpose image data processing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00Geometric image transformations in the plane of the image

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • Medical Informatics (AREA)
  • Animal Behavior & Ethology (AREA)
  • Pathology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Biophysics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • General Health & Medical Sciences (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Multimedia (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Optics & Photonics (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)
  • Image Processing (AREA)

Abstract

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

Description

本発明は、複数の画像データの間で位置合わせを行う画像処理装置及び方法、画像処理装置の制御プログラム並びに当該制御プログラムを格納した記録媒体に関する。
光超音波(Photo-acoustic)画像化(イメージング)法は、癌、腫瘍、血管新生疾患など多くの病気の早期臨床診断のための有望な新技術である。光超音波画像化法は、光音響効果を利用している。すなわち、血管などの検査対象物は、短パルスの近赤外照射を吸収した後、超音波を放射する。検査対象物の3次元構造は熱音響波を検出することによって再構成することができる。この技術は、非侵襲的に任意の造影剤なしで、高空間分解能で生体内の血管を可視化することができる。
従来の光超音波画像化法を用いた3次元イメージングにおいては、例えば胸や手などの検査対象物の全体画像を撮影するために、多数の局所領域をスキャンして複数のワンショットボリューム画像データを取得してそれらを合成している(例えば、特許文献1参照)。また、ワンショットボリューム画像データではノイズが強いため、スキャンの際に同じ場所を重複して撮影して、平均化することでノイズを軽減する手法が知られている。
特開2013−152452号公報
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. 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. 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. Frangi, A. F. et al., "Multiscale vessel enhancement filtering", Proceedings of MICCAI, pp. 130-137, 1998. 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.
しかしながら、撮影中に体動による位置ずれが生じると、これらの方法では、画質が悪化するという問題があった。上述のように、既存の位置合わせ手法では、ワンショットボリューム画像データに強いノイズが多く含まれるため、局所領域間の類似度を用いた位置合わせ手法では、類似度の算出に対するノイズの影響が強く、位置合わせが難しい。
また、ノイズの分離と位置合わせ対象物体行列の低ランク最適化を同時に解く方法(以下、RASL法という。)が提案されているが(例えば、非特許文献1参照)、ノイズは疎なデータ欠損であるという仮定をもって問題を解いており、光超音波画像化法で取得した画像では仮定と異なり、うまく機能しないという問題点があった。
本発明の目的は以上の問題点を解決し、従来技術に比較してノイズを除去することができ、しかも高精度で位置合わせを行うことができる画像処理装置及び方法、画像処理装置の制御プログラム並びにそれを格納した記録媒体を提供することにある。
本発明は、光超音波画像はノイズの影響が大きいものの、ノイズには統計的な特徴があり、このノイズの統計的な特徴を利用することによってノイズの影響を大きく低減できることを見いだしたことによる。
第1の発明に係る画像処理装置は、所定の撮像対象物を含む所定領域が撮像された複数の画像データ内の撮像対象物を互いに位置合わせするように位置変換する制御手段を備えた画像処理装置であって、
前記制御手段は、前記複数の画像データをそれぞれ、前記撮像対象物を含む前景画像と、前記撮像対象物とは異なる背景画像であって画素値の平均値から離れるにつれて分布数が減少する所定の統計的分布を有する背景画像と、前記撮像対象物の画像上のデータ欠損を補完する補完誤差画像とに分解し、前記前景画像の類似度を示す第1の評価関数を実質的に最小とするように、もしくは前記第1の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第1の位置合わせ処理を実行することを特徴とする。
前記画像処理装置において、前記制御手段は、
(1)前記複数の前景画像の低ランク性を評価する評価関数であって、低ランクほど低い値を取る低ランク評価関数と、
(2)位置合わせ後の前記各背景画像の同一の画素の平均値が、前記各背景画像の平均値から離れるにつれてペナルティを与えるペナルティ関数、
(3)前記補完誤差画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数
の和が実質的に最小とするように、もしくは当該和の最小値を含む所定範囲となるように位置合わせを行うことで、前記第1の位置合わせ処理を実行することを特徴とする。
また、前記画像処理装置において、前記制御手段は、所定の最適化アルゴリズム法を用いて、少なくとも、前記第1の評価関数が実質的に最小とするような位置合わせ、もしくは前記第1の評価関数の最小値を含む所定範囲となるような位置合わせ後の画像データ及び前記位置変換の関数の差分を求めることで、前記第1の位置合わせ処理を実行することを特徴とする。
ここで、前記画像処理装置において、前記制御手段は、前記前景画像と、前記補完誤差画像と、前記背景画像とを出力して可視化することを特徴とする。
さらに、前記画像処理装置において、前記制御手段はさらに、前記第1の位置合わせ処理の前に、前記前景画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数のみが実質的に最小とすること、もしくは前記スパース性評価関数の最小値を含む所定範囲となることを第2の評価関数とし、前記第2の評価関数が実質的に最小とすること、もしくは前記第2の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第2の位置合わせ処理を実行することを特徴とする。
ここで、前記制御手段は、所定の最適化アルゴリズム法を用いて、前記第2の評価関数が実質的に最小とするような位置合わせ、もしくは前記第2の評価関数の最小値を含む所定範囲となるような位置合わせ後の画像データ及び前記位置変換の関数の差分を求めることで、前記第2の位置合わせ処理を実行することを特徴とする。
前記画像処理装置において、前記撮像された複数の画像データは、光超音波画像の画像データであることを特徴とする。
また、前記画像処理装置において、前記撮像対象物は血管であり、前記背景画像はガウス分布に従うノイズであることを特徴とする。
さらに、前記画像処理装置において、前記第1の位置合わせ処理と前記第2の位置合わせ処理のうちの少なくとも一方により位置合わせした前記複数の画像データを平均化することで画像データを作成することを特徴とする。
またさらに、前記画像処理装置において、前記第1の位置合わせ処理と前記第2の位置合わせ処理のうちの少なくとも一方により位置合わせした画像データのうち、前記各前景画像の画像データを平均化することで画像データを作成することを特徴とする。
第2の発明に係る画像処理方法は、所定の撮像対象物を含む所定領域が撮像された複数の画像データの撮像対象物を互いに位置合わせするように位置変換する制御手段を備えた画像処理装置のための画像処理方法であって、
前記制御手段が、前記複数の画像データをそれぞれ、前記撮像対象物を含む前景画像と、前記撮像対象物とは異なる背景画像であって画素値の平均値から離れるにつれて分布数が減少する所定の統計的分布を有する背景画像と、前記撮像対象物の画像上のデータ欠損を補完する補完誤差画像とに分解し、前記前景画像の類似度を示す第1の評価関数を実質的に最小とするように、もしくは前記第1の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第1の位置合わせ処理を実行することを特徴とする。
前記画像処理方法において、前記第1の位置合わせ処理を実行することは、前記制御手段が、
(1)前記複数の前景画像の低ランク性を評価する評価関数であって、低ランクほど低い値を取る低ランク評価関数と、
(2)位置合わせ後の前記各背景画像の同一の画素の平均値が、前記各背景画像の平均値から離れるにつれてペナルティを与えるペナルティ関数、
(3)前記補完誤差画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数
の和が実質的に最小とするように、もしくは前記和の最小値を含む所定範囲となるように位置合わせを行うことを含むことを特徴とする。
また、画像処理方法において、前記制御手段が、前記第1の位置合わせ処理を実行するステップを実行する前に、前記前景画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数のみが実質的に最小とすること、もしくは前記スパース性評価関数の最小値を含む所定範囲となることを第2の評価関数とし、前記第2の評価関数が実質的に最小とするように、もしくは前記第2の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第2の位置合わせ処理を実行するステップをさらに含むことを特徴とする。
第3の発明に係る画像処理装置の制御プログラムは、前記制御手段が、前記第1の位置合わせ処理を実行するステップを含むことを特徴とする。
前記画像処理装置の制御プログラムにおいて、前記制御手段がさらに、前記第2の位置合わせ処理を実行するステップを含むことを特徴とする。
第4の発明に係るコンピュータにより読取可能な記録媒体は、前記制御手段が、前記第1の位置合わせ処理を実行するステップを実行する前記画像処理装置の制御プログラムを格納したことを特徴とする。
前記コンピュータにより読取可能な記録媒体において、前記制御手段がさらに、前記第2の位置合わせ処理を実行するステップを含むことを特徴とする。
第5の発明に係る画像処理装置は、所定の撮像対象物を含む所定領域が撮像された複数の画像データを互いに位置合わせするように位置変換する制御手段を備えた画像処理装置であって、
前記制御手段が、前記撮像対象物を含む前景画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数のみが実質的に最小とすること、もしくは前記スパース性評価関数の最小値を含む所定範囲となることを第2の評価関数とし、前記第2の評価関数が実質的に最小とするように、もしくは前記第2の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第2の位置合わせ処理を実行することを特徴とする。
前記画像処理装置において、前記制御手段は、所定の最適化アルゴリズム法を用いて、前記第2の評価関数が実質的に最小とするような位置合わせ、もしくは前記第2の評価関数の最小値を含む所定範囲となるような位置合わせ後の画像データ及び前記位置変換の関数の差分を求めることで、前記第2の位置合わせ処理を実行することを特徴とする。
また、前記画像処理装置において、前記撮像された複数の画像データは、光超音波画像の画像データであることを特徴とする。
さらに、前記画像処理装置において、前記撮像対象物は血管であることを特徴とする。
第6の発明に係る画像処理方法は、所定の撮像対象物を含む所定領域が撮像された複数の画像データを互いに位置合わせするように位置変換する制御手段を備えた画像処理装置のための画像処理方法であって、
前記制御手段が、前記撮像対象物を含む前景画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数のみが実質的に最小とすること、もしくは前記スパース性評価関数の最小値を含む所定範囲となることを第2の評価関数とし、前記第2の評価関数が実質的に最小とするように、もしくは前記第2の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第2の位置合わせ処理を実行するステップを含むことを特徴とする。
第7の発明に係る画像処理装置の制御プログラムは、前記制御手段が、前記第2の位置合わせ処理を実行するステップを含むことを特徴とする。
第8の発明に係るコンピュータにより読取可能な記録媒体は、前記制御手段が、前記第2の位置合わせ処理を実行するステップを含むことを特徴とする。
従って、本発明に係る画像処理装置及び方法等によれば、従来技術に比較してノイズを除去することができ、しかも高精度で位置合わせを行うことができる。
本発明の一実施形態に係る画像処理装置10を含む画像処理システムの構成例を示すブロック図である。 従来技術に係る光超音波画像化法によって得られた2次元最大強度照射(MIP)のワンショットボリューム画像データ(高輝度カーブは血管である)の複数の例を示す画像である。 (a)は図2のワンショットボリューム画像データの背景画像の強度分布を示すグラフであり、(b)は図2の複数のワンショットボリューム画像データの同一画素位置での輝度平均の分布を示すグラフである。 図1の画像処理装置10によって実行される位置合わせ処理を示すフローチャートである。 (a)は実画像の平均化画像であり、(b)は位置ずれが発生した合成画像であり、(c)はB−FFD法により得られた平均化画像であり、(d)はLog−demon法により得られた平均化画像であり、(e)はRASL法により得られた平均化画像であり、(f)は本実施形態にかかる位置合わせ処理により得られた平均化画像である。 実施形態及び従来法による画像に対する評価結果であって、位置ずれレベルに対する位置ギャップのメトリックを示すグラフである。 実施形態に係る画像処理装置10に位置変換された画像データの画像である。 図7Aの位置変換された画像データから分離された前景画像データAの画像である。 図7Aの位置変換された画像データから分離された後景画像データBの画像である。 図7Aの位置変換された画像データから分離された補完画像データCの画像である。 MIP法により得られた位置合わせなしの平均化画像データの画像である。 本実施形態に係る位置合わせ処理で得られた平均化画像データの画像である。 (a)及び(b)は図8Aの拡大図であり、(c)及び(d)はそれぞれ(a)及び(b)の画像に対して本実施形態に係る位置合わせ処理を行って得られた画像である。
以下、比較例及び本発明に係る実施形態について図面を参照して説明する。なお、以下の実施形態において、同様の構成要素については同一の符号を付している。
1.比較例に係る背景技術及び実施形態の概要.
人体の全体の一部分をキャプチャする光超音波画像化法で撮像するために、検査対象物の局所領域をスキャンしてそれらをマージするマルチスキャン及びレジストレーションシステムが開発されてきた。光超音波画像化法によって再構成されたシングルショットボリューム画像データは通常、音と光の散乱及びセンサーレイアウトの制限によって引き起こされる多大なノイズを受けるという問題点があった。
本発明の目的は、これらのノイズを有する多数のワンショットボリューム画像データから高品質の3次元画像データを生成することにある。これを達成することで、例えば血管がはっきりと見えるようになる。このようなノイズを低減するために、画像平均化技術がしばしば使用されている。背景画像でのランダムノイズの平均値が小さい定数となる場合には、直線相関を有する前景画像データが明確になる。このことは、画像データのサンプルが静的であるか、又は検査対象物の大きさが動きよりもはるかに大きいような画像平均化処理に限り有効である。
しかしながら、局所的にスキャンされた複数のワンショットボリューム画像データをマージすることにより、再構成された光超音波画像化法による画像の場合には、患者は多くの場合、走査中に移動するので前記の仮定はほとんど適用できず、血管が薄い(0.5mm)であるという事実を考慮すれば、小さな動きが悪画質に影響を与えることは容易に想像できる。従って、そのような身体の動きに対応できる複数のワンショットボリューム画像データを正確に位置合わせする必要がある。
低ランクベースの位置合わせ(アライメント)のフレームワークは、非常にデンスな(所定のしきい値以上の高密度な)ノイズや小さな前景画像を含んでいるような挑戦的な画像データの位置合わせを行うことができるという可能性を秘めている。低ランクベースの位置合わせのフレームワークでは、観測行列Dは、与えられた複数の画像からまず発生され、ここで、第i番目の行はi番目のベクトル画像を示す。もし本発明者らが成功裏によく位置合わせできた画像を作成して位置合わせ画像間の差異を除去できれば、このことはデータ行列Dを低ランクにさせることができる。低ランクベースの位置合わせのフレームワークは、
Figure 2018025841
に示すように、観測行列Dを低ランク成分に変換する変換関数τと誤差成分αを検索することで実現する。前記式の演算子は、画像データに対して位置変換関数を用いて変換を行う演算子である。この最適化問題は基本的に不良設定問題であるため、通常の最適化の誤差成のいくつかの仮定を導入する必要がある。
このフレームワークの主たる目的は、誤差成分αを設定することと、光超音波(PA)画像化法により得られたPA画像データの特徴に基づいて、小さな前景と、非常にデンスなノイズと、破損した前景を含む複数のPAボリューム画像データを位置合わせするための最適化方法を設定することにある。この目的を達成するために、発明者らは極小に陥ることなく正確に、小さな血管領域を抽出する粗密(デンスおよびスパース)アプローチを提案する。粗い位置合わせのために、複数のPAボリューム画像データに対して、まずFrangiフィルタ(例えば、非特許文献4参照)を適用することで背景ノイズを減少させ血管の強度を向上させる。次いで、マルチスケールのピラミッド法は、位置合わせされた画像データを低ランクの画像データに変換する位置変換関数を最適化するために使用される。前のステップで位置変換関数の初期推定値が与えられるとき、発明者らはさらに、PA画像データの特徴に基づいて、入力される画像データを、低ランクの前景画像(血管)と、デンスな(所定のしきい値以上の高密度な)背景画像(多大で厄介なノイズ)と、スパース(所定のしきい値未満の低密度な、まばらな)補完(誤差)画像成分(前景の相補的な部分)とに分解することによりの画像データを位置合わせする。本発明の実施形態における重要な新規性は、PA画像化の最適化において統計的な特徴(先験的な特徴情報)を有効的に用いていることであり、ここで、各空間ボクセルにおけるデンスなノイズの背景画像の平均値が一定になるように強制される。
後述する実施例において、実際のPA画像データを用いた実験結果は、本実施形態に係る提案手法が、非常に緻密なノイズ及び破損を受けることなく複数のボリューム画像データを位置合わせして血管を見出すために有効であることを実証し、クリアな血管構造を持つ高品質の3次元ボリューム画像データを得ることできることを示している。
多数の公知方法は、医用画像化装置において複数の画像を位置合わせすることを提案してきた。しかし、PA画像において観測される血管は非常に密なノイズの背景にわたってまばらに分散している傾向があり、前景画像と背景画像の両方の画像から計算された画像の相似性に依存したこれらの位置合わせ方法を用いることでそのようなPA画像を位置合わせすることを非常に困難にしている。画像の破損の影響を軽減するために、非特許文献1において開示されたロバスト位置合わせ法(RASL法)が提案されている。このRASL法は医用画像化装置において成功裏にレジストレーションの問題に適用されている。RASL法は、画像変換を最適化するが、位置合わせすべき複数の画像間の差がまばらな破損に起因していることを前提として、画像を低ランクの前景成分と、まばらな誤差成分とに分離している。RASL法によって処理された複数の画像を比較すると、血管領域は、PA画像化の複数の血管画像において深刻なノイズを有して前景画像が小さい。この場合において、誤差成分はもはやスパースにならないので、RASL法は正しく動作しない。すなわち、画像データを適当に変換することに代えて、前記誤差成分を調整することで低ランク成分が最適化される。
2.実施形態に係る画像処理装置.
図1は本発明の一実施形態に係る画像処理装置10を含む画像処理システムの構成例を示すブロック図である。実施形態に係る画像処理装置10は、入力画像データ処理部21と、粗い位置合わせ処理部22と、細かい位置合わせ処理部23と、出力画像データ処理部24と、データ記憶部25とを備えて構成される。画像処理装置10には、画像処理装置10の動作を指示する操作部11と、撮像デバイス12と、ディスプレイ13とが接続される。ここで、画像処理装置10は例えばディジタル計算機等のコンピュータにより構成される。
図1において、撮像デバイス12は例えば撮像センサとフロントエンド回路とを含んで構成され、例えば、ある局所領域を複数回撮像して、位置合わせすべき複数のボリューム画像データ(例えば2次元又は3次元)を生成して入力画像データ処理部21に出力する。データ記憶部25は画像処理装置10で実行される位置合わせ処理を含む画像処理を実行するために必要なプログラム(図4の画像位置合わせ処理プログラム)及びデータ(例えば、初期位置変換行列τ,…,τ)を予め格納するとともに、処理途中の画像データ及びパラメータセットを緩衝格納し、出力画像データ及び出力パラメータセットを格納する。なお、画像位置合わせ処理プログラムはデータ記憶部25に予め格納してもよいし、もしくは、例えばDVD、CDなどの光ディスクの記録媒体27に格納された後、画像処理装置10の動作開始時に光ディスクドライブ装置26に当該記録媒体27を挿入して光ディスクドライブ装置26により当該画像位置合わせプログラムを読み出してデータ記憶部25にロードして実行してもよい。
入力画像データ処理部21は入力される複数のボリューム画像データを所定の形式のベクトル化された3次元ボリューム画像データに変換して粗い位置合わせ処理部22に出力する。粗い位置合わせ処理部22は後述する図4の粗い位置合わせ処理(ステップS1〜S6)を実行することで所定のパラメータセット(A,Δτ)及び位置変換された画像データ
Figure 2018025841
を出力する。細かい位置合わせ処理部23は後述する図4の細かい位置合わせ処理(ステップS7〜S12)を実行することで所定のパラメータセット(A,B,C,Δτ)及び位置変換された画像データ
Figure 2018025841
を出力する。出力画像データ処理部24は位置変換された画像データを表示用の所定の形式の画像データに変換してディスプレイ13に出力して表示し、もしくは外部の画像データ記憶装置に出力する。
以下、位置変換と三成分の分解の最適化推定処理を行う統一されたフレームワークを含む位置合わせ処理について以下詳述する。
3.光超音波画像化法により得られた画像データの特徴.
図2は従来技術に係る光超音波画像化法によって得られた2次元最大強度照射(MIP)のワンショットボリューム画像データ(高輝度カーブは血管である)の複数の例を示す画像である。ここで、各ボリューム画像データの強度範囲は−1から1までの範囲である。まず、これらの画像で以下の3成分を定義する。
(1)第1の成分は、高強度で血管を表す前景画像データである。
(2)第2の成分は、血管以外の部分にあるノイズを表す背景画像データである。
(3)第3の成分は、前景での破損が発生する複数のパーツを表す補完(誤差)画像データである。ここで、破損はPA画像において、例えば血管の部分が生体内での不均一な光照射によって破壊されて発生する。例えば、図2において、破線の円の中に血管が表示されていない。対応する部分は左の画像に表示されている。この不足している部分は補完(誤差)画像データに相当する。
図3(a)は図2のワンショットボリューム画像データの背景画像の強度分布を示すグラフであり、図3(b)は図2の複数のワンショットボリューム画像データの同一画素位置での輝度平均の分布を示すグラフである。本発明者らは、実際のPA画像データを調べることにより、図3に示すように、以下の統計的な特徴(先験的な特徴情報)を発見した。
(1)背景画像の強度分布が大きな標準偏差(0.249)を持つガウス分布に近いこと。
(2)各空間的なボクセルにおける位置合わせされた背景画像の平均値の分布は、所定の平均値(−0.001)と、比較的小さな標準偏差(0.081)とを有するガウス分布に近いこと。
ここで、先験的な特徴情報は、上記(2)の比較的小さな標準偏差を有するガウス分布に近いということを利用して、目的関数のペナルティ関数を設計している。標準偏差が小さいということは、平均値より離れる確率は低くなると考えるので、複数のワンショットボリューム画像データの同一画素位置での輝度平均値が、分布の平均値より離れると大きな値をとるようにしている。
これとは対照的に、血管の平均強度は、血管が良好に位置合わせされたショットボリューム画像データの間で線形に相関するので、高い値を有する。特に、この統計的な特徴(先験的な特徴情報)はしばしば、本実施形態に係る画像処理装置10等の画像化システムのみならず、従来暗黙のうちに使用する、ノイズを広く軽減するために従来技術に係る画像平均化技術が使用される多くのアプリケーションにおいて、高い確率で真実である。バックグラウンドにおいては、輝度値はガウス分布に従うので、無相関だが、血管の輝度強度は、バックグラウンドと比べ、高い値を取る。そのため、位置合わせに成功している場合は、血管上の同じ場所のボクセルにおいて、輝度値は、類似しており(相関があり)、高い値を取る。言い換えると、位置合わせがうまくいっていない場合は、ある画像ではあるボクセル値がバックグラウンドであるのに、ある画像ではそれに対応するボクセル値が血管領域で、輝度値が高いということがおき、この場合、そのボクセルの輝度平均値も高い値に引っ張られ、ペナルティ値が高くなる。また、ノイズの軽減に従来よく利用される従来技術に係る画像平均化手法は、位置合わせができていることが想定されている場合、前景領域の同じ場所は同等の輝度値を持っており(相関があり)、ノイズは無相関であることを暗黙に利用しており、同じ場所のボクセルが前景領域であると、複数の画像間で同様の値を持つので、その平均値は、その場所の本来の輝度値に近づき、ノイズ領域のボクセルは、無相関なので、その平均が小さい値になり、その結果、前景領域とノイズ領域のSN比がよくなる。
本実施形態では、本発明者らは、3つの成分の分解推定するために、これらの統計的な特徴(先験的な情報)を有効的に用いる。
4.本実施形態に係る、粗い位置合わせから細かい位置合わせ(より高精度な位置合わせ)を行う低ランクベースの位置合わせ処理
緻密なノイズと小さな前景画像領域を含むPA画像の場合には、現在の低ランクに基づく方法は、誤差成分を調整することにより、低ランクの成分を最適化する傾向がある。この問題を回避するために、発明者らは粗から密へのアプローチを用いる。
粗い位置合わせのために、発明者らは、誤差成分を考慮せずに、ノイズ低減のデータを使用して変換を最適化するためのマルチスケールピラミッド法を用いる。これは、局所的な最小点に陥ることなく、適切な変換を見つけるのに役立つ。さらに高精度な位置合わせ画像を得るためには、本発明者らは、PA画像化の統計的特徴に基づいて、位置合わせすべき画像データを、低ランクの前景画像、高密度のノイズの背景画像、及びスパースな補完(誤差)画像データの3個の成分に適切に分解しながら、変換関数を最適化する。この手順の詳細は、次のセクションで説明する。
4.1.ピラミッドを用いた低ランク最小化による粗い位置合わせ処理.
いまスキャン中に移動した可能性があるあるサンプルである、複数m個のボリューム画像データ
Figure 2018025841
が与えていると仮定する。ここで、n次元の空間Rは例えば、幅、高さ、奥行きの3次元の空間である。このステップは概ね、前処理された前景画像データを使用して、これらのボリューム画像データを整列させて得られる。ここで、背景ノイズは、Frangiフィルタ(例えば、非特許文献4参照)を通過させ、血管を増強することによって低減される。ここで、管状構造らしさを計算するために、ヘッセ行列の固有ベクトルを使用している。フィルタリングされた画像データは背景において偽陽性のデータを含み、血管が破損している。しかし、その画像データは粗い位置合わせのためには十分である。代わりに、誤差成分を調整する適切な変換関数を見つけるために、差分αを考慮せずに、前処理されたデータ行列Dをただ単に、低ランク成分の画像データ
Figure 2018025841
に変換する変換関数を最適化する。ここで、
Figure 2018025841
は観測データ行列であり、nは各ショットボリューム画像データ内のボクセル数である。また、τは例えばアフィン変換などの特定の位置線形変換(以下、位置変換という)関数を表す。例えば非特許文献1によれば、位置変換関数τの変化が小さいとき、
Figure 2018025841
は位置変換関数τについて線形化することで近似することができ、ランク最小化は次式(1)に示すように核ノルム
Figure 2018025841
を最小化するように緩和される。
Figure 2018025841
(1)
ここで、Δτは各反復において位置変換関数τの分散を表す。また「s.t.」は「such that」の略であって「制約条件」を意味し、右側の条件式を満たすように、左側の式における前記核ノルムを最小化するときのパラメータA,τを求めることを意味する。以下、「s.t.」については以下同様である。また、左側の式の演算子はノルムを演算する演算子である。さらに、Jは次式で定義されるi番目のデータの
アン行列である。
Figure 2018025841
そして、εはm次元空間
Figure 2018025841
のための標準的な基底を表す。式(1)については、非特許文献1のRASL法と同様に、公知の拡張ラグランジュ未定乗数法(ALM)アルゴリズムを用いて解くことができる(例えば、非特許文献6参照)。この最適化は、順次より細かいガウシアンピラミッドを有する画像データを成功裏に位置合わせするマルチスケールピラミッド法に適用される。
4.2.三成分の分解を用いる細かい位置合わせ処理(例えば粗い位置合わせ処理の後に行う).
前記前処理された画像データは、前のステップでの最適化はノイズを考慮していないので小さな位置合わせエラーにつながる、偽陽性と偽陰性のノイズを含むかもしれない。本実施形態では、ノイズを考慮して、元のボリュームを使用することにより粗い位置合わせの結果の画像データをさらに高精度な位置合わせ処理を行う。本実施形態に係る細かい位置合わせ処理では、密なノイズが含まれている画像を処理するために、位置合わせ処理で緻密なノイズ項を導入し、緻密なノイズの多い背景画像の特徴を用いて、低ランクの前景データA、高密度のノイズ背景データB、及びスパースな補完データCに分解して、データの位置を調整することで位置合わせを行う。具体的には、同じ位置において位置合わせされた密なノイズの平均値の分布は、所定値の平均と小さな分散を有するガウス分布となる。本実施形態では、この統計的な特徴を次式のように、分解問題の定式化に導入する。
Figure 2018025841
(2)
ここで、演算子minの対象物は左の式の3つの項を含み、演算子minは左の式の3つの項の和が実質的に最小となるような前景画像データ行列A及び位置変換行列τを求める。ここで、当該和がその最小値を含む所定範囲となるようにしてもよい。A及びBはそれぞれ、位置合わせされた、前景画像データ行列、背景画像データ行列を表し、Cは補完行列を表す。行列A、B、及びCの和は位置合わせされた画像データ
Figure 2018025841
に等しい。また、式(2)の「ベクトルa」は、すべての要素がすべて1である1×mのベクトルである。さらに、式(2)の「ベクトルb」は、すべての要素がガウス分布の平均値であるすべて一定値(スカラー)bである1×nのベクトルである。ここで、実際のPA画像データの、発明者らの観察に基づいてスカラーbに0を設定する。式(2)の第2項
Figure 2018025841
は、同じ位置において位置合わせされた背景画像データの平均値がガウス分布の平均値から離れてとどまるときに、最適化関数にペナルティを課する項である。また、重み付けパラメータλに次式を設定し、
Figure 2018025841
また、重み付けパラメータλに次式を設定する。
Figure 2018025841
ここで、発明者らの実験において、κ=1、κ=3に設定した。
4.3.全体の位置合わせ処理
図4は図1の画像処理装置10によって実行される全体の位置合わせ処理を示すフローチャートである。全体の位置合わせ処理は、粗い位置合わせ処理と、細かい位置合わせ処理とから構成されるが、本発明はこれに限らず、細かい位置合わせ処理のみを実行してもよい。ここで、画像間の位置ずれが比較的大きな場合は、粗い位置合わせ処理が必要であるが、画像間の位置ずれが比較的小さい場合は、粗い位置合わせ処理を実行することなく、細かい位置合わせ処理を実行してもよい。
図4において、ステップS1〜S6までの処理は、図1の粗い位置合わせ処理部22によって実行される粗い位置合わせ処理であり、ステップS7〜S12までの処理は、図1の細かい位置合わせ処理部23によって実行される細かい位置合わせ処理である。
図4のステップS1においてベクトル化された3次元ボリュームデータ
Figure 2018025841
及び初期線形変換行列τ,…,τに基づいて、次式(3)を用いて
アン行列Jを計算する。
Figure 2018025841
(3)
ここで、初期線形変換行列τは次式で表される。
Figure 2018025841
次いで、ステップS2において、次式(4)を用いて画像データの変換及び正規化を行って線形変換されたデータ
Figure 2018025841
を計算する。
Figure 2018025841
(4)
次いで、ステップS3において、次式(5)を用いて、ALMアルゴリズムを用いた線形化された凸最適化を行ってパラメータセット(A,Δτ,Y)を計算する。
Figure 2018025841
(5)
式(5)の右辺は評価関数である拡張ラグランジュ関数であり、その第1項は、前景のデータ行列Aのノルムを最小化するとき(すなわち、位置合わせされたとき)の前景画像データを示す。第2項及び第3項はALMアルゴリズムを用いて線形化された凸最適化を行うための追加項である。ここで、Yは詳細後述するように、ALMアルゴリズムにおけるラグランジュ乗数行列である。また、右辺第2項の演算子<P,Q>は行列P,Q間の内積を演算する演算子であり、右辺第2項は前景画像データ行列Aと位置変換後のラグランジュ乗数行列Yとの内積を表す。さらに、右辺第3項は位置変換された画像データと前景画像データとの間の位置ずれの誤差の距離に対応する。式(5)の右辺では、演算子minの対象物は3つの項を含み、演算子minは3つの項の和を最小化するときの前景画像データ行列A及び変換行列の差分Δτを求める。
次いで、ステップS4において位置変換関数τを位置変換関数の差分Δτを用いて次式のように更新する。
Figure 2018025841
そして、ステップS5において、次式(6)の第1の終了条件(粗い位置合わせ誤差が所定値未満である)を満たすか否かが判断される。なお、tolは終了条件判断用しきい値許容値(偏差)である。
Figure 2018025841
(6)
ステップS5において、YESのときはステップS6に進む一方、NOのときはステップS1に戻る。ステップS6ではパラメータセット(A,Δτ)をデータ記憶部25に出力して一時的に格納する。ここで、Aは位置合わせ後の画像データであり、以上で粗い位置合わせ処理を終了する。
次いで、以下の細かい位置合わせ処理を実行する。
ステップS7において、ベクトル化された3次元ボリュームデータ
Figure 2018025841
及び初期線形変換行列τ,…,τに基づいて、次式(7)を用いて
アン行列Jを計算する。
Figure 2018025841
(7)
次いで、ステップS8において次式(8)を用いて画像データの位置変換及び正規化を行って線形変換されたデータ
Figure 2018025841
を計算する。
Figure 2018025841
(8)
次いで、ステップS9において、次式(9)を用いて、ALMアルゴリズムを用いた線形化された凸最適化を行ってパラメータセット(A,B,C,Δτ,Y)を計算する。
Figure 2018025841
(9)
式(9)の右辺は評価関数である拡張ラグランジュ関数であり、その第1項は、前景のデータ行列Aのノルムを最小化するとき(すなわち、位置合わせされたとき)の前景画像データを示す。第2項は背景画像データ行列Bが平均値0から離れたときにペナルティ値を付加する項である。第3項はスパース誤差画像データ行列Cの誤差に対応する。ここで、位置合わせ後の前景画像同士が類似していることを評価基準として低ランクを導入し、当該第3項の誤差が少なくなるようにL−1ノルムを導入して最適化問題を構成している。さらに、第4項及び第5項はALMアルゴリズムを用いて線形化された凸最適化を行うための追加項である。式(9)の右辺では、演算子minの対象物は5つの項を含み、演算子minは5つの項の和を最小化するときの前景画像データ行列A、背景画像データ行列B、スパースな誤差画像データ行列C、及び変換行列の差分Δτを求める。
次いで、ステップS10において位置変換関数τを位置変換関数の差分Δτを用いて次式のように更新する。ここでの繰り返し処理により、評価関数を実質的に最小化し、もしくは評価関数の最小値を含む所定の範囲(例えば+/−0%から〜+/−5%ないし+/−10%)となるように行うことができる。
Figure 2018025841
次いで、ステップS11において、式(10)の第2の終了条件(細かい位置合わせ誤差が所定値未満である)を満たすか否かが判断される。なお、tolは終了条件判断用しきい値許容値(偏差)である。
Figure 2018025841
(10)
ステップS11において、YESのときはステップS12に進む一方、NOのときはステップS7に戻る。ステップS12では、パラメータセット(A,B,C,Δτ)を出力画像データ処理部24に出力する。ここで、Aは位置合わせ後の画像データであり、以上で細かい位置合わせ処理を終了し、当該画像位置合わせ処理を終了する。
図4の2つの各反復では、位置変換関数τの現在の推定値で変換を線形化することにより、線形化及び最適化問題を解決した。この最適化問題は上述のALMアルゴリズムを用いて、線形化された凸最適化として解ける。
本実施形態においては、次式の拡張増補ラグランジュ関数を最小化する最適化問題に帰着させて解を求めている。
Figure 2018025841
(11)
ここで、
Figure 2018025841
はラグランジュ乗数行列であり、μは正のスカラー計数である。この問題は次式のように反復して最小化される。
Figure 2018025841
(12)
ここで、svd(・)は特異値分解の演算子を表し、縮小演算子
Figure 2018025841
は核ノルムを最適化するための特異値閾値化アルゴリズムのためのものである。μは単調に増加する正のシーケンスであり、本実施形態では、非特許文献1に従って、次式のように設定した。
Figure 2018025841
なお、提案手法の計算時間は、他の方法と同等であった。
5.実験評価.
まず、本発明者らは実際の3次元PA画像データであるショットボリューム画像データに基づいて、本実施形態に係る提案された画像位置合わせ処理のロバスト性と精度を評価した。本発明者らは、スキャン中に被験者の身体を固定することにより、ピクセルレベルのずれが存在しないように、真実の画像データを5セット作成した。次いで、合成された真実のショットボリューム画像データに位置ギャップの異なるレベルを追加し、元のボリューム画像データをランダムに各軸(x,y,z)に沿ったボクセルが均一に分布した擬似乱数[−d d]で置換した。ボクセルサイズは、[0.25 0.25 0.25](mm)で血管の厚さは0.5〜2mmであった。ここで、発明者らは、1から10のずれ量dを設定した。全体では、発明者らは10×5個のデータセットを生成し、ここで、各データセットは、21から36個のショットボリューム画像データを含む。
比較のために、発明者らは、以下の従来技術に係る方法により得られたものと発明者らの結果を比較した。
(1)類似機能を最小限に抑えながら、制御点を最適化するBスプライン自由形状変形法(以下、B−FFD法という)(例えば、非特許文献3参照):。
(2)単純な最近傍検索で画像間の地点別の対応を発見するスペクトルLog−demonレジストレーション(以下、Log−demon法という)(例えば、非特許文献4参照);
(3)公平な比較のために血管拡張データに適用した、ピラミッド法を備えたRASL法(例えば、非特許文献1参照)。
図6は実施形態及び従来法による画像に対する評価結果であって、位置ずれレベルに対する位置ギャップのメトリックを示すグラフである。図6において、縦軸の値が小さい場合は、各方法が正常にデータセットを位置合わせできたことを示す。これにより、本実施形態に係る提案手法はミスアライメントレベル7までうまく機能していることを示しており、このことは本当の身体の動きに対処するには十分であるというべきである。他の方法は、非常に低い位置ずれレベルで、良好に画像データの位置合わせを行うことができなかった。
図5(a)は実画像の平均化画像であり、図5(b)は位置ずれが発生した合成画像であり、図5(c)はB−FFD法により得られた平均化画像であり、図5(d)はLog−demon法により得られた平均化画像であり、図5(e)はRASL法により得られた平均化画像であり、図5(f)は本実施形態にかかる位置合わせ処理により得られた平均化画像である。すなわち、図5は、位置ずれレベル6で、それぞれの方法で得られた位置合わせ結果の平均画像の例を示している。
図5(a)は真実の画像データを示し、ここで血管がはっきりと見える。図5(b)から明らかなように、合成ずれが画質を大幅に悪化させることを示している。図5(f)の提案方法の結果は、明確な血管を示しており、画像が真実の画像データ(図5(a))に類似している。このことは、提案手法が適切にショットボリューム画像データを位置合わせできることを示している。対照的に、血管は、さらに他の方法の結果で識別することが困難であった(図5(c)、図5(d)、図5(e))。
以下、本発明者らの提案方法は、実際のPAデータを3つの成分にどのように分解したかについて精査する。
図7Aは実施形態に係る画像処理装置10により位置変換された画像データ
Figure 2018025841
の画像である。図7Bは図7Aの位置変換された画像データ
Figure 2018025841
から分離された前景画像データAの画像である。図7Cは図7Aの位置変換された画像データ
Figure 2018025841
から分離された後景画像データBの画像である。図7Dは図7Aの位置変換された画像データ
Figure 2018025841
から分離された補完画像データCの画像である。
すなわち、図7A〜図7Dはミスアラインメントレベル6での結果を示している。背景雑音は正しく後景画像データBに分解されている。すべての血管は明らかに、低ランクの前景画像成分Aに表示され、そして、血管のまばらな破損は、補完成分Cに表示されている。この結果は、提案方法が、元のノイズの多いPA画像が、正常に3個の意味のある成分分解できたことを示している。
次に、本発明者らは実際の体動に対する提案手法のロバスト性を精査した。実際のデータを取得するには、手は、螺旋パターンで局所領域を複数回走査する広視野範囲のPA画像化システムを使用してスキャンした。そして、実際の身体の動きは、スキャン中に追加された。ショットボリューム画像データの総数は2048だった。総スキャン時間は約2分だった。連続したローカルショットボリューム画像データは85%にわたって重畳された。すべての画像データは512×512×100画素のデータであった。
各ショットボリューム画像データは、局所エリアデータの一部であるので、発明者らはまず、16個から59個のショットボリューム画像データを含むローカルエリア(168×168×100画素)のデータセットを生成した。次いで、発明者らは、各データセットにおいてショットボリューム画像データの位置合わせを行った。その後、全てのデータを生成するために位置合わせされた平均ボリューム画像データをともに縫い合わせた。比較のために、位置合わせなしの画像平均化技術が画像データデータに適用された。
図8AはMIP法により得られた位置合わせなしの平均化画像データの画像である。また、図8Bは本実施形態に係る位置合わせ処理で得られた平均化画像データの画像である。さらに、図8C(a)及び(b)は図8Aの拡大図であり、図8C(c)及び(d)はそれぞれ図8C(a)及び(b)の画像に対して本実施形態に係る位置合わせ処理を行って得られた画像である。
図8A、図8C(a)及び(b)は位置合わせしない平均画像データを示す。これらの画像では、血管のコントラストが低く、いくつかの血管が体動によりぼやけている。これに対して、図8B、図8C(c)及び(d)は提案方法によるレジストレーション結果を示す。これらの画像では、血管は明らかに元の画像よりも鮮明であり、位置合わせしないときは平均画像データを見ることが困難であったいくつかの血管が見えるようになった。ここで、図8Bにおけるボックスに三回現れる血管について注意すべきであり、2本の静脈は動脈と並行して実行する場合にぼやけていない。これらの画像の特徴は、解剖学の医師によって確認され、その結果は、発明者らの提案方法が大幅にPA画像化の画質を大幅に改善したことを示している。
6.結論
発明者らは、PA画像の特徴を検討し、高品質の3次元ボリューム画像データを生成するために、PA画像の画像化のためのレジストレーション方法を提案した。提案方法を用いることで血管がはっきりと見えるようになった。位置合わせされた背景データの同じ位置の値の平均値が一定になる傾向を有するという統計的特徴を導入することで、発明者らの位置合わせ方法は、強力なノイズと大きなずれとを有する挑戦的なPA画像データを処理することが可能である。実際のデータセットについての実験結果は、提案手法の有効性を実証しており、著しく画質が向上し、従来技術に比較して最良の位置合わせ精度を達成した。
7.非特許文献1のRASL法と本実施形態の方法との相違点について(まとめ).
以下の説明では、ALMアルゴリズムに係る項を省略して、RASL法と本実施形態の方法の概要について以下に説明する。
7−1.RASL法
RASL法はRobust−PCAを応用した手法であって、入力画像である位置合わせ前の入力画像データDから、前景画像Aと、ノイズ画像Eと、位置合わせのための位置変換τとを推定する。入力画像データDを位置変換した位置合わせ後のデータを
Figure 2018025841
として、位置合わせ後の前景画像データ同士が類似していることの評価基準として、低ランクを導入し、差分誤差はなるべく少なくなるようにL−1ノルムを導入することで、最適化問題の評価関数を次式のように設定している。
Figure 2018025841
ここで、第1項は位置合わせ後の前景画像データであり、第2項はスパース(まばらな)誤差である。
制約条件
Figure 2018025841
が非線形なので、位置変換関数である変数τのローカルな傾きとして、Δτに関して線形で表現することで問題を緩和し、解けるようにしている。ここで、変数τは繰り返し最適化の中で固定して、次式のごとくΔτを最適化の目的変数としていることを特徴としている。
Figure 2018025841
Figure 2018025841
7−2.3次元への適応におけるRASL法の問題点.
(問題点1)以下の非凸問題がある。
RASL法では、位置ずれは小さいと仮定しており、3次元でかつ前景画像(血管)が小さいと、局所解に陥りやすいという問題点がった。具体的には、RASL法をそのまま単純に3次元に適用すると正しく動作しないという問題点があった。
(問題点2)誤差画像領域はスパース誤差であると仮定しており、PA画像では、ノイズ領域はきわめて大きいという問題点があった。
7−3.本実施形態にかかる提案手法.
7−3−1.粗い位置合わせ処理.
RASL法の第1の問題点を解決するために、位置合わせ処理により誤差成分で低ランクを実現することで、次式のごとく誤差成分なしで位置合わせの最適化のみで低ランクを検索している。
Figure 2018025841
前記の式では、位置変換してなる、位置合わせ後の前景画像同士が類似しているという探索基準で低ランクを導入し、前景画像Aの核ノルムが最小となるように前景画像A及び位置変換関数を探索する。
7−3−2.細かい位置合わせ処理.
RASL法では、ノイズ領域は小さいと仮定しており、PA画像と異なるという問題点があった。本実施形態では、次式のごとく、前景画像A、スパースな誤差画像Eに加え、背景画像Bを導入したことを特徴としている。背景画像Bはガウス分布などの統計的特徴(位置あわせ後の同じ座標の値の平均値は0に近づく)に従い、無相関である特徴を利用している。
Figure 2018025841
前記の式では、第1項は前景画像Aの核ノルムであり、第2項は背景画像Bがガウス分布に従って統計的に無相関であり、平均値が0から離れるとペナルティ値を付与する項であり、第3項はスパース誤差画像E(実施形態における補完画像Cと同義である)のL1ノルムである。
前記細かい位置合わせ処理では、ノイズを強い分散を持つ高密度なガウス分布のノイズと、疎なデータ欠損との2つからなると仮定した問題を設定し、ガウス分布から得られるノイズの平均の分布は分散が弱いガウス分布に従うという統計上の特性を用いて、その問題を解く最適化手法を提案することで、PA画像の実態にあった仮定で位置合わせ問題を解くことができ、従来技術に比較して前景画像についての位置合わせ処理において良好な結果が得られる。
8.変形例.
以上の実施形態において、図4に示すように、粗い位置合わせ処理を実行した後、細かい位置合わせ処理を実行しているが、本発明はこれに限らず、粗い位置合わせ処理を実行せずに、細かい位置合わせ処理のみを実行してもよい。このことは入力される複数のボリューム画像データ間での差分が比較的小さい場合に有効であり、処理時間の短縮を行うことができる。
以上の実施形態において、光超音波画像(PA画像)に対する位置合わせ処理について説明しているが、本発明はこれに限らず、所定の対象物又は対象場所の複数の撮像画像に対して位置合わせ処理を行ってもよい。
以上の実施形態において、前景画像として血管の一例を占めているが、本発明はこれに限らず、所定のしきい値以上の強度(画素値)を有する対象物の画像であってもよい。
以上の実施形態において、位置変換関数として、アフィン変換を用いているが、本発明はこれに限らず、例えば平行移動変換、回転変換、拡大又は縮小変換、鏡映変換、射影変換などのグルーバルドメインでのパラメトリック線形変換関数であればよい。
以上の実施形態において、PA画像の統計的特徴を有する分布としてガウス分布を用いているが、本発明はこれに限らず、以下の統計的分布を用いてもよい。
(1)ノイズの分布:位置合わせ対象画像が多い場合は、基本的にどんな分布でもよい。例えば一様分布、対数正規分布、カイ2乗分布、ワイブル分布、ガンマ分布などである。ただし、位置合わせ対象画像が少ない場合は、中心(平均)から離れるに従って、減少していく分布であり、例えば正規分布、t分布、双曲線正割分布、ローレンツ分布、ラプラス分布などを用いることができる。
(2)平均化後の分布:中心(平均)から離れるに従って、減少していく分布ならばよく、例えば正規分布、t分布、双曲線正割分布、ローレンツ分布、ラプラス分布などの分布を用いることができる。
以上の実施形態においては、上述の評価関数を拡張ラグランジュ未定乗数法のALMアルゴリズムを用いて線形化された凸最適化を行って位置変換後の画像データ及び位置変換関数の差分などを求めているが、本発明はこれに限らず、拡張ラグランジュ未定乗数法に代えて所定の多変数探索法などの最適化アルゴリズム法を用いて位置変換後の画像データ及び位置変換関数の差分などを求めてもよい。例えば、以下の最適化アルゴリズム法を用いてもよい。
(1)accelerated proximal gradient(APG) algorithm;
(2)Robust−PCAで用いられているDual Method;
(3)Singular Value Thresholding。
以上の実施形態においては、図4の細かい位置合わせ処理において、パラメータセット(A,B,C,Δτ)を出力しているが、本発明はこれに限らず、少なくともパラメータセット(A,Δτ)を出力してもよい。
以上の実施形態においては、低ランク化の方法として「核ノルム」を最小化しているが、本発明はこれに限らず、「核ノルム」に代えて、「行列のランクを評価し、低ランクほど小さい値を取る評価関数」を用いてもよい。
以上の実施形態においては、補完誤差画像がスパースであることを評価する方法として、「L1ノルム」を利用しているが、本発明はこれに限らず、スパース性(行列やベクトルにおいて、非0の要素が少ないことをいう)の指標として、例えば「L0ノルム」などの「補完誤差画像のスパース性を評価し、スパースであるほど小さな値を取る評価関数」を用いてもよい。
なお、撮像された複数の画像データを用いて画像を作成する画像処理装置において最終結果の画像データを得るために、前記位置合わせ処理により位置合わせした前記複数の画像データ(撮像された複数のオリジナル画像データ)を平均化することで、最終結果の画像データを作成してもよいし、もしくは、前記位置合わせ処理により位置合わせした画像のうち、前記複数の前景画像データを平均化することで、最終結果の画像データを作成してもよい。
なお、画像処理装置10は、前記分解された、前記前景画像と、前記補完誤差画像と、前記背景画像とを外部装置に出力して可視化してもよい。
以上詳述したように、本発明に係る画像処理装置及び方法等によれば、従来技術に比較してノイズを除去することができ、しかも高精度で位置合わせを行うことができる。
本発明に係る画像処理装置及び方法等は以下のために広く用いることができる。
(1)医療画像化装置:非侵襲で以下の血管系異常を、2次元又は3次元可視化することができ、早期診断、健康診断を行うことができる。例えば、皮膚毛細血管、がんの新生血管、関節の新生血管、動脈狭窄、閉塞、プラーク、脳酸素消費、血管走行異常などである。本実施形態では、特に、血管の自動抽出装置に有用である。
(2)製品検査装置:食品又は工業製品などの劣化、傷などの検査を行うことができる。
10…画像処理装置、
11…操作部、
12…撮像デバイス、
13…ディスプレイ、
21…入力画像データ処理部、
22…粗い位置合わせ処理部、
23…細かい位置合わせ処理部、
24…出力画像データ処理部、
25…データ記憶部、
26…光ディスクドライブ装置、
27…記録媒体。

Claims (24)

  1. 所定の撮像対象物を含む所定領域が撮像された複数の画像データ内の撮像対象物を互いに位置合わせするように位置変換する制御手段を備えた画像処理装置であって、
    前記制御手段は、前記複数の画像データをそれぞれ、前記撮像対象物を含む前景画像と、前記撮像対象物とは異なる背景画像であって画素値の平均値から離れるにつれて分布数が減少する所定の統計的分布を有する背景画像と、前記撮像対象物の画像上のデータ欠損を補完する補完誤差画像とに分解し、前記前景画像間の類似度を示す第1の評価関数を実質的に最小とするように、もしくは前記第1の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第1の位置合わせ処理を実行することを特徴とする画像処理装置。
  2. 前記制御手段は、
    (1)前記複数の前景画像の低ランク性を評価する評価関数であって、低ランクほど低い値を取る低ランク評価関数と、
    (2)位置合わせ後の前記各背景画像の同一の画素の平均値が、前記各背景画像の平均値から離れるにつれてペナルティを与えるペナルティ関数、
    (3)前記補完誤差画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数
    補完の和が実質的に最小とするように、もしくは当該和の最小値を含む所定範囲となるように行うことで、前記第1の位置合わせ処理を実行することを特徴とする請求項1記載の画像処理装置。
  3. 前記制御手段は、所定の最適化アルゴリズム法を用いて、少なくとも、前記第1の評価関数が実質的に最小とするような位置合わせ、もしくは前記第1の評価関数の最小値を含む所定範囲となるような位置合わせ後の画像データ及び前記位置変換の関数の差分を求めることで、前記第1の位置合わせ処理を実行することを特徴とする請求項1又は2記載の画像処理装置。
  4. 前記制御手段は、前記前景画像と、前記補完誤差画像と、前記背景画像とを出力して可視化することを特徴とする請求項1〜3のうちのいずれか1つに記載の画像処理装置。
  5. 前記制御手段はさらに、前記第1の位置合わせ処理の前に、前記前景画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数のみが実質的に最小とすること、もしくは前記スパース性評価関数のみが前記スパース性評価関数の最小値を含む所定範囲となることを第2の評価関数とし、前記第2の評価関数が実質的に最小とするように、もしくは前記第2の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第2の位置合わせ処理を実行することを特徴とする請求項1〜4のうちのいずれか1つに記載の画像処理装置。
  6. 前記制御手段は、所定の最適化アルゴリズム法を用いて、前記第2の評価関数が実質的に最小とするような位置合わせ、もしくは前記第2の評価関数の最小値を含む所定範囲となるような位置合わせ後の画像データ及び前記位置変換の関数の差分を求めることで、前記第2の位置合わせ処理を実行することを特徴とする請求項5記載の画像処理装置。
  7. 前記撮像された複数の画像データは、光超音波画像の画像データであることを特徴とする請求項1〜6のうちのいずれか1つに記載の画像処理装置。
  8. 前記撮像対象物は血管であり、前記背景画像はガウス分布に従うノイズであることを特徴とする請求項1〜7のうちのいずれか1つに記載の画像処理装置。
  9. 前記第1の位置合わせ処理と前記第2の位置合わせ処理のうちの少なくとも一方により位置合わせした前記複数の画像データを平均化することで画像データを作成することを特徴とする請求項1〜8のうちのいずれか1つに記載の画像処理装置。
  10. 前記第1の位置合わせ処理と前記第2の位置合わせ処理のうちの少なくとも一方により位置合わせした画像データのうち、前記各前景画像の画像データを平均化することで画像データを作成することを特徴とする請求項1〜8のうちのいずれか1つに記載の画像処理装置。
  11. 所定の撮像対象物を含む所定領域が撮像された複数の画像データの撮像対象物を互いに位置合わせするように位置変換する制御手段を備えた画像処理装置のための画像処理方法であって、
    前記制御手段が、前記複数の画像データをそれぞれ、前記撮像対象物を含む前景画像と、前記撮像対象物とは異なる背景画像であって画素値の平均値から離れるにつれて分布数が減少する所定の統計的分布を有する背景画像と、前記撮像対象物の画像上のデータ欠損を補完する補完誤差画像とに分解し、前記前景画像の類似度を示す第1の評価関数を実質的に最小とするように、もしくは前記第1の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第1の位置合わせ処理を実行することを特徴とする画像処理方法。
  12. 前記第1の位置合わせ処理を実行することは、前記制御手段が、
    (1)前記複数の前景画像の低ランク性を評価する評価関数であって、低ランクほど低い値を取る低ランク評価関数と、
    (2)位置合わせ後の前記各背景画像の同一の画素の平均値が、前記各背景画像の平均値から離れるにつれてペナルティを与えるペナルティ関数、
    (3)前記補完誤差画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数
    の和が実質的に最小とするように、もしくは前記和の最小値を含む所定範囲となるように位置合わせを行うことを含むことを特徴とする請求項11記載の画像処理方法。
  13. 前記制御手段が、前記第1の位置合わせ処理を実行するステップを実行する前に、前記前景画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数のみが実質的に最小とすること、もしくは当該スパース性評価関数のみが当該スパース性評価関数の最小値を含む所定範囲となることを第2の評価関数とし、前記第2の評価関数が実質的に最小とするように、もしくは前記第2の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第2の位置合わせ処理を実行するステップをさらに含むことを特徴とする請求項11又は12記載の画像処理方法。
  14. 前記制御手段が、請求項11又は12記載の前記第1の位置合わせ処理を実行するステップを含むことを特徴とする画像処理装置の制御プログラム。
  15. 前記制御手段がさらに、請求項13記載の前記第2の位置合わせ処理を実行するステップを含むことを特徴とする請求項14記載の画像処理装置の制御プログラム。
  16. 前記制御手段が、請求項11又は12記載の前記第1の位置合わせ処理を実行するステップを実行する前記画像処理装置の制御プログラムを格納したことを特徴とする、コンピュータにより読取可能な記録媒体。
  17. 前記制御手段がさらに、請求項13記載の前記第2の位置合わせ処理を実行するステップを含むことを特徴とする、請求項14記載のコンピュータにより読取可能な記録媒体。
  18. 所定の撮像対象物を含む所定領域が撮像された複数の画像データを互いに位置合わせするように位置変換する制御手段を備えた画像処理装置であって、
    前記制御手段が、前記撮像対象物を含む前景画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数のみが実質的に最小とするように、もしくは前記スパース性評価関数のみが前記スパース性評価関数の最小値を含む所定範囲となることを第2の評価関数とし、前記第2の評価関数が実質的に最小とするように、もしくは前記第2の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第2の位置合わせ処理を実行することを特徴とする画像処理装置。
  19. 前記制御手段は、所定の最適化アルゴリズム法を用いて、前記第2の評価関数が実質的に最小とするように、もしくは前記第2の評価関数の最小値を含む所定範囲となるような位置合わせ後の画像データ及び前記位置変換の関数の差分を求めることで、前記第2の位置合わせ処理を実行することを特徴とする請求項18記載の画像処理装置。
  20. 前記撮像された複数の画像データは、光超音波画像の画像データであることを特徴とする請求項18又は19記載の画像処理装置。
  21. 前記撮像対象物は血管であることを特徴とする請求項18〜20のうちのいずれか1つに記載の画像処理装置。
  22. 所定の撮像対象物を含む所定領域が撮像された複数の画像データを互いに位置合わせするように位置変換する制御手段を備えた画像処理装置のための画像処理方法であって、
    前記制御手段が、前記撮像対象物を含む前景画像のスパース性を評価する評価関数であって、スパースなほど低い値を取るスパース性評価関数のみが実質的に最小とすること、もしくは前記スパース性評価関数のみが前記スパース性評価関数の最小値を含む所定範囲となることを第2の評価関数とし、前記第2の評価関数が実質的に最小とすること、もしくは前記第2の評価関数の最小値を含む所定範囲となるように位置合わせされた前景画像の画像データを求める第2の位置合わせ処理を実行するステップを含むことを特徴とする画像処理方法。
  23. 前記制御手段が、請求項22記載の前記第2の位置合わせ処理を実行するステップを含むことを特徴とする画像処理装置の制御プログラム。
  24. 前記制御手段が、請求項22記載の前記第2の位置合わせ処理を実行するステップを含むことを特徴とする、コンピュータにより読取可能な記録媒体。
JP2018531909A 2016-08-05 2017-08-01 画像処理装置及び方法 Active JP6799331B2 (ja)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
JP2016154658 2016-08-05
JP2016154658 2016-08-05
PCT/JP2017/027839 WO2018025841A1 (ja) 2016-08-05 2017-08-01 画像処理装置及び方法

Publications (2)

Publication Number Publication Date
JPWO2018025841A1 true JPWO2018025841A1 (ja) 2019-06-06
JP6799331B2 JP6799331B2 (ja) 2020-12-16

Family

ID=61073599

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2018531909A Active JP6799331B2 (ja) 2016-08-05 2017-08-01 画像処理装置及び方法

Country Status (2)

Country Link
JP (1) JP6799331B2 (ja)
WO (1) WO2018025841A1 (ja)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111316323B (zh) * 2019-01-31 2023-06-06 深圳市瑞立视多媒体科技有限公司 一种三维轨迹数据的异常值处理方法及装置

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2012531279A (ja) * 2009-06-29 2012-12-10 ヘルムホルツ・ツェントルム・ミュンヒェン・ドイチェス・フォルシュンクスツェントルム・フューア・ゲズントハイト・ウント・ウムベルト(ゲーエムベーハー) 吸収マップを定量的に抽出するための熱音響撮像

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2012531279A (ja) * 2009-06-29 2012-12-10 ヘルムホルツ・ツェントルム・ミュンヒェン・ドイチェス・フォルシュンクスツェントルム・フューア・ゲズントハイト・ウント・ウムベルト(ゲーエムベーハー) 吸収マップを定量的に抽出するための熱音響撮像

Also Published As

Publication number Publication date
WO2018025841A1 (ja) 2018-02-08
JP6799331B2 (ja) 2020-12-16

Similar Documents

Publication Publication Date Title
US8842936B2 (en) Method, apparatus, and program for aligning images
JP2016041247A (ja) 医療用画像のコンピュータ支援による解析方法
JP2003153082A (ja) 画像の位置合わせ装置および画像処理装置
Ko et al. Rigid and non-rigid motion artifact reduction in X-ray CT using attention module
JP2002157593A (ja) 異常陰影検出方法および装置
US20220138936A1 (en) Systems and methods for calcium-free computed tomography angiography
Išgum et al. Automated aortic calcium scoring on low‐dose chest computed tomography
Xu et al. Cardiac CT motion artifact grading via semi-automatic labeling and vessel tracking using synthetic image-augmented training data
Jeon et al. Mm-net: Multiframe and multimask-based unsupervised deep denoising for low-dose computed tomography
JP6799331B2 (ja) 画像処理装置及び方法
Sigit et al. Development of Healthcare Kiosk for Checking Heart Health
Carminati et al. Reconstruction of the descending thoracic aorta by multiview compounding of 3-d transesophageal echocardiographic aortic data sets for improved examination and quantification of atheroma burden
Yan et al. Calcium removal from cardiac CT images using deep convolutional neural network
Huang et al. MLNAN: Multi-level noise-aware network for low-dose CT imaging implemented with constrained cycle Wasserstein generative adversarial networks
Habijan et al. Centerline tracking of the single coronary artery from x-ray angiograms
Razeto et al. Accurate registration of coronary arteries for volumetric CT digital subtraction angiography
Shamonin et al. Automatic lung lobe segmentation of COPD patients using iterative B-spline fitting
Yamada et al. ROI extraction of chest CT images using adaptive opening filter
Karmakar et al. Detailed investigation of lumen-based tomographic co-registration
Ali et al. Elastic registration of high-resolution 3D PLI data of the human brain
Hai Wavelet-based image fusion for enhancement of ROI in CT image
Luengo-Oroz et al. Extraction of the coronary artery tree in cardiac computer tomographic images using morphological operators
Cao et al. A two-stage framework for optical coherence tomography angiography image quality improvement
Fazlali et al. Robust catheter identification and tracking in X-ray angiographic sequences
Razeto et al. Accurate, fully-automated registration of coronary arteries for volumetric CT digital subtraction angiography

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20200406

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20200825

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20200930

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20201020

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20201113

R150 Certificate of patent or registration of utility model

Ref document number: 6799331

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250