ＩＶＲ画像の画質改善法
【課題】ＩＶＲ画像について、ノイズを除去すると同時に、着目対象であるガイドワイヤや血管系等の線状の陰影を強調することができる方法を提供する。
【解決手段】画質改善すべき元のＩＶＲ画像（Ｓ１）に対しＩＣＡ・Shrinkageフィルタ処理がなされ（Ｓ２）、ノイズ信号が除去された第１のＩＶＲ画像が生成される（Ｓ３）。第１のＩＶＲ画像に対し多重スケールフィルタ処理がなされ（Ｓ４）、その画像中の線状陰影の信号を強調した第２のＩＶＲ画像が生成される（Ｓ５）。第２のＩＶＲ画像のコントラストが反転せしめられた後、第１および第２のＩＶＲ画像が合成され（Ｓ６）、画質の改善されたＩＶＲ画像が生成される（Ｓ７）。
【解決手段】画質改善すべき元のＩＶＲ画像（Ｓ１）に対しＩＣＡ・Shrinkageフィルタ処理がなされ（Ｓ２）、ノイズ信号が除去された第１のＩＶＲ画像が生成される（Ｓ３）。第１のＩＶＲ画像に対し多重スケールフィルタ処理がなされ（Ｓ４）、その画像中の線状陰影の信号を強調した第２のＩＶＲ画像が生成される（Ｓ５）。第２のＩＶＲ画像のコントラストが反転せしめられた後、第１および第２のＩＶＲ画像が合成され（Ｓ６）、画質の改善されたＩＶＲ画像が生成される（Ｓ７）。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、ＩＶＲ（Interventional Radiology）画像の画質を改善する方法に関するものである。
【背景技術】
【０００２】
ＩＶＲは、画像診断的介入治療法または放射線診断技術を応用した治療法とも呼ばれ、近年、血管系病変の治療等において果たす役割は増大しており、今後も、非常に有効な治療手段としてより一層発展することが期待されている。そして、この治療法を支える重要な要素の１つとして、放射線の透視の際の画質向上に関する技術がある。
【０００３】
ＩＶＲにおいては、長時間の治療による患者の放射線被爆を最小限に抑えるため、通常の放射線撮影時の数十分の一程度の強度の放射線が照射され、透視画像が取得される。このため、放射線量の不足によって、粒子性のノイズ（量子ノイズ）が、着目対象であるガイドワイヤ、カテーテルおよび血管系の視認性を悪くしている。
そして、ＩＶＲでは、放射線による透視画像中においてガイドワイヤ、カテーテルおよび血管系等の線状の陰影が、良好な信号対雑音比（Ｓ／N比）で観察できることが極めて重要となっている。
【０００４】
そこで、ガイドワイヤ、カテーテルおよび血管系等の線状の陰影を強調すると同時に、背景のノイズを除去する手法が必要とされる。この場合、従来のノイズ除去法においては、フーリエ変換やwavelet変換等が用いられており、フーリエ変換やwavelet変換では、三角関数やwavelet関数が基底関数として用いられ、ノイズを除去すべき対象画像は、これらの基底関数に基づいて成分分解される。そして、各成分についてフィルタリングがなされ、対象画像からノイズが除去される。ところで、三角関数およびwavelet関数は、いずれも、その関数形が予め数学的に決められていて、画像に基づいて導出されたものではなく、画像の処理に適した基底関数とはいえない。このため、フーリエ変換やwavelet変換等によっては、元の画像の信号成分を失うことなくノイズのみを除去することが困難であり、十分に満足のいくノイズ除去処理ができなかった（例えば、非特許文献１参照）。
また、血管等の画像信号を強調するフィルタとして、多重スケール（Multi-scale）フィルタがあるが、背景にある骨等の画像が見えにくくなったり、ノイズが強調されてしまったりするという問題があった（例えば、非特許文献２参照）。
すなわち、従来技術においては、ＩＶＲ画像中のノイズ除去と線状陰影の強調とを同時に達成し得る手法はこれまでに存在しなかった。
【０００５】
【非特許文献１】R.D. Nowak and R.Baraniuk, "Wavelet Domain Filtering for Photon Imaging Systems", IEEE Transactions on Image Processing, Volume: 8 Issue: 5 , May 1999 Page(s): 666-678
【非特許文献２】Frangi A. et al.,"Multiscale vessel enhancement filtering", Lecture Notes in Computer Science, Springer, LNCS 1496, pp.130-137 (1998).
【発明の開示】
【発明が解決しようとする課題】
【０００６】
したがって、本発明の課題は、ＩＶＲ画像について、ノイズを除去すると同時に、着目対象であるガイドワイヤや血管系等の線状の陰影を強調することができる方法を提供することにある。
【課題を解決するための手段】
【０００７】
上記課題を解決するため、本発明は、（ａ）画質改善すべき元のＩＶＲ画像に対しＩＣＡ・Shrinkageフィルタ処理を行い、ノイズ信号を除去した第１のＩＶＲ画像を生成するステップと、（ｂ）前記第１のＩＶＲ画像に対し多重スケールフィルタ処理を行い、その画像中の線状陰影の信号を強調した第２のＩＶＲ画像を生成するステップと、（ｃ）前記第２のＩＶＲ画像のコントラストを反転させた後、前記第１および第２のＩＶＲ画像を合成して画質の改善されたＩＶＲ画像を生成するステップと、を有することを特徴とするＩＶＲ画像の画質改善法を構成したものである。
【０００８】
上記構成において、必要に応じて、前記ステップ（ｃ）において前記第１および第２の画像を合成するとき、それぞれの画像に重み係数を付与することもできる。
【０００９】
上記課題を解決するため、また、本発明は、（ａ’）画質改善すべき同一の元のＩＶＲ画像を２つ用意し、（ｂ’）一方の前記元のＩＶＲ画像に対しＩＣＡ・Shrinkageフィルタ処理を行い、ノイズ信号を除去した第１のＩＶＲ画像を生成するとともに、他方の前記元のＩＶＲ画像に対し多重スケールフィルタ処理を行い、その画像中の線状陰影の信号を強調した第２のＩＶＲ画像を生成するステップと、（ｃ’）前記第２のＩＶＲ画像のコントラストを反転させた後、前記第１および第２のＩＶＲ画像を合成して、画質の改善されたＩＶＲ画像を生成するステップと、を有することを特徴とするＩＶＲ画像の画質改善法を構成したものである。
【００１０】
上記構成において、必要に応じて、前記ステップ（ｃ’）において前記第１および第２の画像を合成するとき、それぞれの画像に重み係数を付与することもできる。
【００１１】
また、本発明の好ましい実施例によれば、前記ＩＣＡ・Shrinkageフィルタ処理は、（i）予め与えられた複数のサンプル画像から、独立成分分析（ＩＣＡ）によって、前記サンプル画像の特徴を統計的に表す基底関数の組を抽出するステップと、（ii）前記ＩＶＲ画像をＩＣＡ変換して、前記ＩＶＲ画像を前記基底関数の線形和の形に成分分解するステップと、（iii）前記ＩＶＲ画像を前記ＩＣＡ変換したときの各変換係数を、Shrinkage関数によってフィルタ処理して前記ＩＶＲ画像に含まれるノイズ成分を除去するステップと、（iv）前記フィルタ処理した後の各変換係数を用いてＩＣＡ逆変換を行い、ノイズの除去されたＩＶＲ画像を得るステップと、からなっている。
【００１２】
この場合、前記ステップ（iii）において、前記Shrinkage関数として、ポアッソンノイズに対するクロスバリデーション（cross-validation）法に基づくShrinkage関数を用いることが好ましく、また、前記ステップ（i）において、前記サンプル画像を、前記ＩＶＲ画像に類似する対象物について前記ＩＶＲ画像の取得時と近似した条件下で取得したデジタル画像に基づいて作成することが好ましい。
【００１３】
本発明の別の好ましい実施例によれば、前記多重スケールフィルタ処理は、
（i）強調すべき線状陰影の幅に対応するスケールｓをもつガウス関数Ｇ（ｘ，ｙ）を生成するステップと、
（ii）Ｇ（ｘ，ｙ）のx方向の微分∂Ｇ（ｘ，ｙ）／∂ｘ（＝Ｇｘ（ｘ，ｙ））およびｙ方向の微分∂Ｇ（ｘ，ｙ）／∂ｙ（＝Ｇｙ（ｘ，ｙ））を求め、さらに、Ｇｘ（ｘ，ｙ）のｘ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｘ２（＝Ｇｘｘ（ｘ，ｙ））およびｙ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｘ∂ｙ（＝Ｇｘｙ（ｘ，ｙ））と、Ｇｙ（ｘ，ｙ）のｘ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｙ∂ｘ（＝Ｇｙｘ（ｘ，ｙ））およびｙ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｙ２（＝Ｇｙｙ（ｘ，ｙ））とを求めるステップと、
（iii）前記ＩＶＲ画像の各画素値Ｌ（ｘ，ｙ）の微分∂２Ｌ（ｘ，ｙ）／∂ｘ２（＝Ｌｘｘ（ｘ，ｙ））、∂２Ｌ（ｘ，ｙ）／∂ｘ∂ｙ（＝Ｌｘｙ（ｘ，ｙ））、∂２Ｌ（ｘ，ｙ）／∂ｙ∂ｘ（＝Ｌｙｘ（ｘ，ｙ））および∂２Ｌ（ｘ，ｙ）／∂ｙ２（＝Ｌｙｙ（ｘ，ｙ））を、
Ｌｘｘ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｘｘ（ｘ，ｙ）、
Ｌｘｙ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｘｙ（ｘ，ｙ）、
Ｌｙｘ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｙｘ（ｘ，ｙ）、
Ｌｙｙ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｙｙ（ｘ，ｙ）、
ここで、γは任意の定数である、
によって求めるステップと、
（iv）前記ＩＶＲ画像の各画素の位置におけるヘッセ行列Ｈを
【数１】
により求めるステップと、
（v）前記ヘッセ行列Ｈの２つの固有値λ１、λ２（｜λ１｜＜｜λ２｜）および対応する固有ベクトルν１、ν２を求め、
【数２】
を求めるステップと、
（vi）前記ＩＶＲ画像の各画素値を、
【数３】
ここで、ｃおよびβは任意の定数である、
によって求めるステップと、
（vii）前記ガウス関数Ｇ（ｘ，ｙ）のスケールｓの値を変更して、前記ステップ（ii）〜（vi）を繰り返すステップと、
（viii）前記ステップ（ii）〜（vii）を所定の回数繰り返した後、得られたＩｓ（ｘ，ｙ）を用いて、
【数４】
に従って、前記ＩＶＲ画像の各画素値を求めるステップと、からなっている。
【発明の効果】
【００１４】
本発明によれば、ＩＣＡ・Shrinkageフィルタ処理と多重スケールフィルタ処理を組み合わせたことにより、ＩＶＲ画像中のエッジ部分などの細部がぼけることなく、画像全体に分散する量子ノイズの低減を図ることができるとともに、画像中の血管、ガイドワイヤやカテーテル等の直線状の陰影部分を強調することができる。
【発明を実施するための最良の形態】
【００１５】
以下、添付図面を参照して、本発明の好ましい実施例について説明する。図１は、本発明の１実施例によるＩＶＲ画像の画質改善法のフロー図である。図１を参照して、本発明によれば、まず、画質改善すべき元のＩＶＲ画像（図１のＳ１）に対しＩＣＡ・Shrinkageフィルタ処理がなされ（図１のＳ２）、ノイズ信号が除去された第１のＩＶＲ画像が生成される（図１のＳ３）。
【００１６】
図２は、ＩＣＡ・Shrinkageフィルタ処理のフロー図である。以下、図２を参照して、ＩＣＡ・Shrinkageフィルタ処理について説明する。まず、予め与えられたサンプル画像から、独立成分分析（ＩＣＡ）法によって、サンプル画像の特徴を統計的に表す基底関数の組が抽出される（図２のＳ２０、Ｓ２１）。サンプル画像としては、画質を改善すべきＩＶＲ画像により適合した基底関数の組を得るべく、当該ＩＶＲ画像に類似する別のＩＶＲ画像が用いられる。また、ＩＣＡ法によるサンプル画像からの基底関数の抽出は統計的になされるので、できるだけ多数のサンプル画像を用いることが、ＩＶＲ画像により適合した基底関数の組を得るうえで好ましい。
【００１７】
ＩＣＡ法を用いてサンプル画像から基底関数の組を抽出する方法を次に説明する。
ＩＣＡ法によれば、サンプル画像［ｘ１↑，ｘ２↑，・・・，ｘＮ↑］（文字の後ろに「↑」の記号を付した場合はベクトル表示であることを表す。以下同様。）は、基底関数［ａ１↑，ａ２↑，・・・，ａＭ↑］の線形和で表され、その結合係数ｓは統計的に独立である。これを数式で表すと次のようになる。
【００１８】
【数５】
ここで、ｘｉ↑＝Σｓｉｊａｊ↑（ｊについての総和）はｉ番目のサンプル画像であり、行列Ｓの列ベクトルｓｉ↑はその画像の独立な結合係数であり、特徴ベクトルと呼ばれる。また、ａｊ↑はｊ番目の基底関数である。
【００１９】
（１）式を行列式として表すと、
Ｘ↑＝Ａ↑・Ｓ↑ （２）
となる。
サンプル画像セット行列Ｘ↑のみから、基底関数と特徴ベクトルを求めることは、一種のブラインドソース分離問題とみなせるので、ＩＣＡ法を用いて解くことができる。
ＩＣＡでは、あるｎ×ｍの分離行列W↑を用いて、
Ｙ↑＝Ｗ↑・Ｘ↑ （３）
においてＹ↑を計算し、Ｙ↑の各成分が互いに統計的に独立となるようなＷ↑を求める。
Ｗ↑を求めるに当たり、Ｙ↑の各成分が独立か否かを判定する評価基準として、相互情報量（Kullback-Leibler divergence）または高次統計量（kurtosis）等が既に提案されている。この実施例では、Bell & Senjowskyが提案した相互情報量の最大化による手法を用いる。
【００２０】
今、分離信号［ｙ１，ｙ２，・・・，ｙｎ］が互いに独立であるとき、各成分ｙｉの同時分布密度関数ｐ（ｙ↑）＝Ｐ（ｙ１，ｙ２，・・・，ｙｎ）は、各成分ｙｉの周辺密度関数ｐ（ｙｉ）の積で表される。すなわち、
【００２１】
【数６】
となる。評価関数として、ｐ（ｙ↑）とΠｐ（ｙｉ）（ｉについての積）との距離を表すKullback-Leibler Divergence （Ｄ）が用いられる。Kullback-Leibler Divergence （Ｄ）は、エントロピーＨと分離行列Ｗ↑を用いて以下のように表される。
【００２２】
【数７】
ここで、
【００２３】
【数８】
Ｈ（ｙ↑）およびＨ（ｙｉ）は、それぞれ、同時分布のエントロピーと周辺分布のエントロピーである。（３）式より、Ｐ（Ｙ↑）＝Ｐ（Ｖ↑）／｜Ｗ↑｜であるので、同時分布エントロピーＨ（ｙ↑）は、次のように書き換えられる。
【００２４】
【数９】
よって、（５）式のKullback-Leibler
divergence （Ｄ）は次のように書き換えられる。
【００２５】
【数１０】
ＩＣＡ法では、ｐ（ｙ↑）とΠｐ（ｙｉ）（ｉについての積）が一致するように分離行列Ｗ↑を求めるが、これは、（９）式のＤを最小にするＷ↑を求める問題となる。
一般に、Ｗ↑を解析的に求めることはできず、（９）式のＤを最小にするＷ↑を求めるには、下記の自然勾配法が用いられる。
Ｗ↑の変化分は、（７）式で与えられる。通常の最急降下法に比べると、収束がよく、計算も簡単である。
【００２６】
【数１１】
こうして、Ｗ↑の更新則は次式で与えられる。
【００２７】
【数１２】
ここで、μは学習率であり、Ｉ↑は単位行列であり、アンサンブル平均は時間平均に置き換えている。なお、ｐ（ｙｉ）が未知であるから、ｇ（ｙ↑）は近似で求めることになる。この実施例では、ｇ（ｙ↑）はシグモイド関数である。（１２）式により、Ｙが独立となるようなＷ↑が求まり、また、（３）式よりＹ↑が求まる。
【００２８】
こうして、この実施例では、対象となるＩＶＲ画像とは別の４枚のＩＶＲ画像を用い、ＩＣＡ法によって基底関数の組を抽出した。
まず、これら４枚のＩＶＲ画像（各々１２８画素×１２８画素のサイズを有する）から任意に８画素×８画素の大きさのパッチ１２,０００個を切り出し、基底関数の組を学習によって抽出するためのサンプル画像とした。この場合、サンプル画像を形成するパッチのサイズは特に限定されないが、抽出される基底関数の組が、画質を改善すべきＩＶＲ画像の空間的特徴を最もよく反映し得るようなパッチのサイズを設定することが好ましい。
【００２９】
そして、各サンプル画像を６４×１の１次元ベクトルに展開し、行列Ｘ↑の列ベクトルとする。それによって、６４×１２,０００のサンプル画像セット行列Ｘ↑が得られる。行列Ｘ↑の行の数はサンプル画像のパッチの個数を表し、列の数はサンプル画像の数を表す。
（１２）式を用いてＷ↑を学習すると、６４×６４の行列Ｗ↑が求められ、このときのＷ↑の各列が基底関数となる。この実施例では、基底関数の数は６４となる。
【００３０】
再び図２を参照して、次に、ＩＶＲ画像が、上で求めた基底関数の組を用いてＩＣＡ変換され、基底関数の線形和の形に成分分解される（図２のＳ２２、Ｓ２３）。これを数式で表すと次のようになる。
【００３１】
【数１３】
【００３２】
ここで、Ｐ↑はＩＶＲ画像であり、［ａ０↑，ａ１↑，・・・，ａＮ−１↑］は基底関数であり、［ｓ０，ｓ１，・・・，ｓＮ−１］はＩＣＡ変換の変換係数である。
【００３３】
そして、ＩＶＲ画像をＩＣＡ変換したときの各変換係数が、Shrinkage関数によってフィルタ処理され、ＩＶＲ画像に含まれるノイズ成分が除去される（図２のＳ２４）。この実施例では、Shrinkage関数として、
【００３４】
【数１４】
ここで、δ２は、ポアッソンノイズの累乗であり、
【数１５】
【００３５】
この関数ｇ（ｓ↑）は、ポアッソンノイズに対するクロスバリデーション（cross-validation ）法に基づくShrinkage関数であり、ポアッソンノイズの除去に最適である。なお、Shrinkage関数の関数形はこれに限定されるものではなく、例えば、次式で表されるような公知の適当な関数形を用いることができる。
【数１６】
こうして、ノイズ分散δ２の増加に伴って、係数ｓの縮小される割合が大きくなる。
【００３６】
さらに、フィルタ処理後の各変換係数を用いてＩＣＡ逆変換が行われ、ノイズ成分の除去されたＩＶＲ画像（第１のＩＶＲ画像）が得られる（図２のＳ２５、Ｓ２６）。これを式で表すと次のようになる。
【００３７】
【数１７】
【００３８】
ここで、Ｐ（ハット）は第１のＩＶＲ画像であり、［ａ０↑，ａ１↑，・・・，ａＮ−１↑］は基底関数であり、［ｓ０（ハット），ｓ１（ハット），・・・，ｓＮ−１（ハット）］はフィルタ処理後の変換係数である。
【００３９】
図５は、この実施例で用いた元のＩＶＲ画像を示したものであり、図６は、図５のＩＶＲ画像に対しＩＣＡ・Shrinkageフィルタ処理を行って得られた第１のＩＶＲ画像を示したものである。図５および図６のＩＶＲ画像の比較から、ＩＣＡ・Shrinkageフィルタ処理によって、画像中の縞模様のような構造的なノイズがかなり除去されており、処理前のＩＶＲ画像を劣化させることなく、ノイズ成分を除去できることがわかる。
【００４０】
なお、上記実施例では、対象ＩＶＲ画像をＩＣＡ変換した後、変換係数をShrinkage関数によってフィルタ処理することでノイズ成分を除去するようにしたが、このフィルタ処理を行うとともに、ＩＶＲ画像における分解された画像部分のうち、目視によりノイズを含んでいると認識された画像部分を除去するようにしてもよい。
【００４１】
ＩＣＡ・Shrinkageフィルタ処理によれば、複数のサンプル画像から独立成分分析（ＩＣＡ）によって学習、抽出した基底関数に基づき、元のＩＶＲ画像をＩＣＡ変換し、変換係数をShrinkage関数によってフィルタ処理することでノイズ成分を除去した後、ＩＣＡ逆変換によって第１のＩＶＲ画像を得るようにしたので、元の画像の信号成分を失うことなくノイズ成分を除去することができる。
【００４２】
再び図１を参照して、本発明によれば、次に、第１のＩＶＲ画像に対し多重スケール（Multi-scale）フィルタ処理がなされ（図１のＳ４）、その画像中の線状陰影の信号を強調した第２のＩＶＲ画像が生成される（図１のＳ５）。図３には、多重スケールフィルタ処理のフロー図が示してある。
【００４３】
図３を参照して、多重スケールフィルタ処理において、まず、強調すべき線状陰影の幅に対応するスケールｓをもつガウス関数Ｇ（ｘ，ｙ）が生成される（図３のＳ３０）。すなわち、ガウス関数Ｇのスケールｓの大きさに応じて、ＩＶＲ画像中の強調される線状の陰影部分の幅が決定される。
【００４４】
次に、Ｇ（ｘ，ｙ）のx方向の微分∂Ｇ（ｘ，ｙ）／∂ｘ（＝Ｇｘ（ｘ，ｙ））およびｙ方向の微分∂Ｇ（ｘ，ｙ）／∂ｙ（＝Ｇｙ（ｘ，ｙ））が求められ、さらに、Ｇｘ（ｘ，ｙ）のｘ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｘ２（＝Ｇｘｘ（ｘ，ｙ））およびｙ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｘ∂ｙ（＝Ｇｘｙ（ｘ，ｙ））と、Ｇｙ（ｘ，ｙ）のｘ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｙ∂ｘ（＝Ｇｙｘ（ｘ，ｙ））およびｙ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｙ２（＝Ｇｙｙ（ｘ，ｙ））が求められる（図３のＳ３１）。
【００４５】
そして、第１のＩＶＲ画像の各画素値Ｌ（ｘ，ｙ）の微分∂２Ｌ（ｘ，ｙ）／∂ｘ２（＝Ｌｘｘ（ｘ，ｙ））、∂２Ｌ（ｘ，ｙ）／∂ｘ∂ｙ（＝Ｌｘｙ（ｘ，ｙ））、∂２Ｌ（ｘ，ｙ）／∂ｙ∂ｘ（＝Ｌｙｘ（ｘ，ｙ））および∂２Ｌ（ｘ，ｙ）／∂ｙ２（＝Ｌｙｙ（ｘ，ｙ））が、
Ｌｘｘ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｘｘ（ｘ，ｙ）、
Ｌｘｙ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｘｙ（ｘ，ｙ）、
Ｌｙｘ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｙｘ（ｘ，ｙ）、
Ｌｙｙ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｙｙ（ｘ，ｙ）、
ここで、γは任意の定数である、
によって求められる（図３のＳ３２）。
【００４６】
さらに、第１のＩＶＲ画像の各画素の位置におけるヘッセ行列Ｈが
【数１８】
によって求められ（図３のＳ３３）、ヘッセ行列Ｈの２つの固有値λ１、λ２（｜λ１｜＜｜λ２｜）および対応する固有ベクトルν１、ν２が求められる（図３のＳ３４）。この場合、λ１の固有ベクトルν１は、関係する線状陰影の軸方向を表し、λ２の固有ベクトルν２は、当該線状陰影の軸に対し垂直な方向を表している。
【００４７】
そして、
【数１９】
が計算され（図３のＳ３５）、ＩＶＲ画像の各画素値が、
【数２０】
ここで、ｃおよびβはいずれも任意の定数である、
によって求められる（図３のＳ３５）。この式において、βの値を大きくすればするほど画像は明るくなり、ｃの値を大きくすればするほど画像は暗くなる。
【００４８】
次に、ガウス関数Ｇ（ｘ，ｙ）のスケールｓが変更され、図３のステップＳ３０〜Ｓ３５が繰り返される。そして、すべてのｓについて計算が完了したとき、得られたＩｓ（ｘ，ｙ）を用いて、
【数２１】
に従って、ＩＶＲ画像の各画素値が求められる（図３のＳ３７）。この場合、計算に使用されるスケールｓの値の組は、対象となるＩＶＲ画像中に現われた血管、ガイドワイヤおよびカテーテル等の線状の陰影部分の各幅に対応して、予め決定される。
【００４９】
再び図１を参照して、本発明の方法によれば、さらに、得られた第２のＩＶＲ画像のコントラストが反転せしめられ、その画像と第１のＩＶＲ画像が合成され（図１のＳ６）、画質の改善されたＩＶＲ画像が生成される（図１のＳ７）。なお、第１および第２のＩＶＲ画像の合成に際し、必要に応じて、各画像に重み係数が付与される。
【００５０】
図７は、図６のＩＶＲ画像（第１のＩＶＲ画像）に対し多重スケールフィルタ処理を行って得られた第２のＩＶＲ画像を示したものである。また、図８は、図７のＩＶＲ画像（第２のＩＶＲ画像）のコントラストを反転させた画像と、図６のＩＶＲ画像（第１のＩＶＲ画像）を合成して得られた、画質の改善されたＩＶＲ画像を示したものである。この場合、第１のＩＶＲ画像の割合を１とし、コントラスト反転した第２のＩＶＲ画像の割合を０．１として両画像を合成し、得られた画像を１．１で割っている。
図５および図８のＩＶＲ画像の比較から、本発明による画質改善法によれば、ＩＶＲ画像におけるノイズ除去と線状陰影の強調を同時に達成することができることがわかる。
【００５１】
上記実施例では、ＩＶＲ画像が静止画像である場合について説明してきたが、本発明の方法は、ＩＶＲ画像が動画像である場合にも適用可能である。この場合には、動画像を１フレーム毎の画像に分解し、分離した画像毎に上記と同様の処理を行い、処理後の各画像を再合成すればよい。
【００５２】
図４は、本発明の別の実施例によるＩＶＲ画像の画質改善法のフロー図である。図４を参照して、この実施例では、画質改善すべき同一の元のＩＶＲ画像が２つ用意される（図４のＳ１０）。そして、一方の元のＩＶＲ画像に対しＩＣＡ・Shrinkageフィルタ処理がなされ、ノイズ信号が除去された第１のＩＶＲ画像が生成される一方（図４のＳ１１、Ｓ１２）、他方の元のＩＶＲ画像に対し多重スケールフィルタ処理がなされ、その画像中の線状陰影の信号が強調された第２のＩＶＲ画像が生成される（図４のＳ１３、Ｓ１４）。
【００５３】
そして、第２のＩＶＲ画像のコントラストが反転せしめられ（図４のＳ１５）、その画像と第１のＩＶＲ画像とが合成され（図４のＳ１６）、画質の改善されたＩＶＲ画像が生成される（図４のＳ１７）。なお、第１および第２のＩＶＲ画像の合成の際に、必要に応じて、各画像に重み係数が付与され得る。
この実施例の場合にも、図１の実施例の場合と同様の効果が得られる。
【図面の簡単な説明】
【００５４】
【図１】本発明の１実施例によるＩＶＲ画像の画質改善法のフロー図である。
【図２】図１中のＩＣＡ・Shrinkageフィルタ処理のフロー図である。
【図３】図１中の多重スケールフィルタ処理のフロー図である。
【図４】本発明の別の実施例によるＩＶＲ画像の画像改善法のフロー図である。
【図５】画質を改善すべきＩＶＲ画像を示したものである。
【図６】図５のＩＶＲ画像に対しＩＣＡ・Shrinkageフィルタ処理を行って得られた第１のＩＶＲ画像を示したものである。
【図７】図６のＩＶＲ画像に対し多重スケールフィルタ処理を行って得られた第２のＩＶＲ画像を示したものである。
【図８】図７のＩＶＲ画像（第２のＩＶＲ画像）のコントラストを反転させ、その画像と図６のＩＶＲ画像（第１のＩＶＲ画像）を合成して得られた、画質の改善されたＩＶＲ画像を示したものである。
【技術分野】
【０００１】
本発明は、ＩＶＲ（Interventional Radiology）画像の画質を改善する方法に関するものである。
【背景技術】
【０００２】
ＩＶＲは、画像診断的介入治療法または放射線診断技術を応用した治療法とも呼ばれ、近年、血管系病変の治療等において果たす役割は増大しており、今後も、非常に有効な治療手段としてより一層発展することが期待されている。そして、この治療法を支える重要な要素の１つとして、放射線の透視の際の画質向上に関する技術がある。
【０００３】
ＩＶＲにおいては、長時間の治療による患者の放射線被爆を最小限に抑えるため、通常の放射線撮影時の数十分の一程度の強度の放射線が照射され、透視画像が取得される。このため、放射線量の不足によって、粒子性のノイズ（量子ノイズ）が、着目対象であるガイドワイヤ、カテーテルおよび血管系の視認性を悪くしている。
そして、ＩＶＲでは、放射線による透視画像中においてガイドワイヤ、カテーテルおよび血管系等の線状の陰影が、良好な信号対雑音比（Ｓ／N比）で観察できることが極めて重要となっている。
【０００４】
そこで、ガイドワイヤ、カテーテルおよび血管系等の線状の陰影を強調すると同時に、背景のノイズを除去する手法が必要とされる。この場合、従来のノイズ除去法においては、フーリエ変換やwavelet変換等が用いられており、フーリエ変換やwavelet変換では、三角関数やwavelet関数が基底関数として用いられ、ノイズを除去すべき対象画像は、これらの基底関数に基づいて成分分解される。そして、各成分についてフィルタリングがなされ、対象画像からノイズが除去される。ところで、三角関数およびwavelet関数は、いずれも、その関数形が予め数学的に決められていて、画像に基づいて導出されたものではなく、画像の処理に適した基底関数とはいえない。このため、フーリエ変換やwavelet変換等によっては、元の画像の信号成分を失うことなくノイズのみを除去することが困難であり、十分に満足のいくノイズ除去処理ができなかった（例えば、非特許文献１参照）。
また、血管等の画像信号を強調するフィルタとして、多重スケール（Multi-scale）フィルタがあるが、背景にある骨等の画像が見えにくくなったり、ノイズが強調されてしまったりするという問題があった（例えば、非特許文献２参照）。
すなわち、従来技術においては、ＩＶＲ画像中のノイズ除去と線状陰影の強調とを同時に達成し得る手法はこれまでに存在しなかった。
【０００５】
【非特許文献１】R.D. Nowak and R.Baraniuk, "Wavelet Domain Filtering for Photon Imaging Systems", IEEE Transactions on Image Processing, Volume: 8 Issue: 5 , May 1999 Page(s): 666-678
【非特許文献２】Frangi A. et al.,"Multiscale vessel enhancement filtering", Lecture Notes in Computer Science, Springer, LNCS 1496, pp.130-137 (1998).
【発明の開示】
【発明が解決しようとする課題】
【０００６】
したがって、本発明の課題は、ＩＶＲ画像について、ノイズを除去すると同時に、着目対象であるガイドワイヤや血管系等の線状の陰影を強調することができる方法を提供することにある。
【課題を解決するための手段】
【０００７】
上記課題を解決するため、本発明は、（ａ）画質改善すべき元のＩＶＲ画像に対しＩＣＡ・Shrinkageフィルタ処理を行い、ノイズ信号を除去した第１のＩＶＲ画像を生成するステップと、（ｂ）前記第１のＩＶＲ画像に対し多重スケールフィルタ処理を行い、その画像中の線状陰影の信号を強調した第２のＩＶＲ画像を生成するステップと、（ｃ）前記第２のＩＶＲ画像のコントラストを反転させた後、前記第１および第２のＩＶＲ画像を合成して画質の改善されたＩＶＲ画像を生成するステップと、を有することを特徴とするＩＶＲ画像の画質改善法を構成したものである。
【０００８】
上記構成において、必要に応じて、前記ステップ（ｃ）において前記第１および第２の画像を合成するとき、それぞれの画像に重み係数を付与することもできる。
【０００９】
上記課題を解決するため、また、本発明は、（ａ’）画質改善すべき同一の元のＩＶＲ画像を２つ用意し、（ｂ’）一方の前記元のＩＶＲ画像に対しＩＣＡ・Shrinkageフィルタ処理を行い、ノイズ信号を除去した第１のＩＶＲ画像を生成するとともに、他方の前記元のＩＶＲ画像に対し多重スケールフィルタ処理を行い、その画像中の線状陰影の信号を強調した第２のＩＶＲ画像を生成するステップと、（ｃ’）前記第２のＩＶＲ画像のコントラストを反転させた後、前記第１および第２のＩＶＲ画像を合成して、画質の改善されたＩＶＲ画像を生成するステップと、を有することを特徴とするＩＶＲ画像の画質改善法を構成したものである。
【００１０】
上記構成において、必要に応じて、前記ステップ（ｃ’）において前記第１および第２の画像を合成するとき、それぞれの画像に重み係数を付与することもできる。
【００１１】
また、本発明の好ましい実施例によれば、前記ＩＣＡ・Shrinkageフィルタ処理は、（i）予め与えられた複数のサンプル画像から、独立成分分析（ＩＣＡ）によって、前記サンプル画像の特徴を統計的に表す基底関数の組を抽出するステップと、（ii）前記ＩＶＲ画像をＩＣＡ変換して、前記ＩＶＲ画像を前記基底関数の線形和の形に成分分解するステップと、（iii）前記ＩＶＲ画像を前記ＩＣＡ変換したときの各変換係数を、Shrinkage関数によってフィルタ処理して前記ＩＶＲ画像に含まれるノイズ成分を除去するステップと、（iv）前記フィルタ処理した後の各変換係数を用いてＩＣＡ逆変換を行い、ノイズの除去されたＩＶＲ画像を得るステップと、からなっている。
【００１２】
この場合、前記ステップ（iii）において、前記Shrinkage関数として、ポアッソンノイズに対するクロスバリデーション（cross-validation）法に基づくShrinkage関数を用いることが好ましく、また、前記ステップ（i）において、前記サンプル画像を、前記ＩＶＲ画像に類似する対象物について前記ＩＶＲ画像の取得時と近似した条件下で取得したデジタル画像に基づいて作成することが好ましい。
【００１３】
本発明の別の好ましい実施例によれば、前記多重スケールフィルタ処理は、
（i）強調すべき線状陰影の幅に対応するスケールｓをもつガウス関数Ｇ（ｘ，ｙ）を生成するステップと、
（ii）Ｇ（ｘ，ｙ）のx方向の微分∂Ｇ（ｘ，ｙ）／∂ｘ（＝Ｇｘ（ｘ，ｙ））およびｙ方向の微分∂Ｇ（ｘ，ｙ）／∂ｙ（＝Ｇｙ（ｘ，ｙ））を求め、さらに、Ｇｘ（ｘ，ｙ）のｘ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｘ２（＝Ｇｘｘ（ｘ，ｙ））およびｙ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｘ∂ｙ（＝Ｇｘｙ（ｘ，ｙ））と、Ｇｙ（ｘ，ｙ）のｘ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｙ∂ｘ（＝Ｇｙｘ（ｘ，ｙ））およびｙ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｙ２（＝Ｇｙｙ（ｘ，ｙ））とを求めるステップと、
（iii）前記ＩＶＲ画像の各画素値Ｌ（ｘ，ｙ）の微分∂２Ｌ（ｘ，ｙ）／∂ｘ２（＝Ｌｘｘ（ｘ，ｙ））、∂２Ｌ（ｘ，ｙ）／∂ｘ∂ｙ（＝Ｌｘｙ（ｘ，ｙ））、∂２Ｌ（ｘ，ｙ）／∂ｙ∂ｘ（＝Ｌｙｘ（ｘ，ｙ））および∂２Ｌ（ｘ，ｙ）／∂ｙ２（＝Ｌｙｙ（ｘ，ｙ））を、
Ｌｘｘ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｘｘ（ｘ，ｙ）、
Ｌｘｙ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｘｙ（ｘ，ｙ）、
Ｌｙｘ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｙｘ（ｘ，ｙ）、
Ｌｙｙ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｙｙ（ｘ，ｙ）、
ここで、γは任意の定数である、
によって求めるステップと、
（iv）前記ＩＶＲ画像の各画素の位置におけるヘッセ行列Ｈを
【数１】
により求めるステップと、
（v）前記ヘッセ行列Ｈの２つの固有値λ１、λ２（｜λ１｜＜｜λ２｜）および対応する固有ベクトルν１、ν２を求め、
【数２】
を求めるステップと、
（vi）前記ＩＶＲ画像の各画素値を、
【数３】
ここで、ｃおよびβは任意の定数である、
によって求めるステップと、
（vii）前記ガウス関数Ｇ（ｘ，ｙ）のスケールｓの値を変更して、前記ステップ（ii）〜（vi）を繰り返すステップと、
（viii）前記ステップ（ii）〜（vii）を所定の回数繰り返した後、得られたＩｓ（ｘ，ｙ）を用いて、
【数４】
に従って、前記ＩＶＲ画像の各画素値を求めるステップと、からなっている。
【発明の効果】
【００１４】
本発明によれば、ＩＣＡ・Shrinkageフィルタ処理と多重スケールフィルタ処理を組み合わせたことにより、ＩＶＲ画像中のエッジ部分などの細部がぼけることなく、画像全体に分散する量子ノイズの低減を図ることができるとともに、画像中の血管、ガイドワイヤやカテーテル等の直線状の陰影部分を強調することができる。
【発明を実施するための最良の形態】
【００１５】
以下、添付図面を参照して、本発明の好ましい実施例について説明する。図１は、本発明の１実施例によるＩＶＲ画像の画質改善法のフロー図である。図１を参照して、本発明によれば、まず、画質改善すべき元のＩＶＲ画像（図１のＳ１）に対しＩＣＡ・Shrinkageフィルタ処理がなされ（図１のＳ２）、ノイズ信号が除去された第１のＩＶＲ画像が生成される（図１のＳ３）。
【００１６】
図２は、ＩＣＡ・Shrinkageフィルタ処理のフロー図である。以下、図２を参照して、ＩＣＡ・Shrinkageフィルタ処理について説明する。まず、予め与えられたサンプル画像から、独立成分分析（ＩＣＡ）法によって、サンプル画像の特徴を統計的に表す基底関数の組が抽出される（図２のＳ２０、Ｓ２１）。サンプル画像としては、画質を改善すべきＩＶＲ画像により適合した基底関数の組を得るべく、当該ＩＶＲ画像に類似する別のＩＶＲ画像が用いられる。また、ＩＣＡ法によるサンプル画像からの基底関数の抽出は統計的になされるので、できるだけ多数のサンプル画像を用いることが、ＩＶＲ画像により適合した基底関数の組を得るうえで好ましい。
【００１７】
ＩＣＡ法を用いてサンプル画像から基底関数の組を抽出する方法を次に説明する。
ＩＣＡ法によれば、サンプル画像［ｘ１↑，ｘ２↑，・・・，ｘＮ↑］（文字の後ろに「↑」の記号を付した場合はベクトル表示であることを表す。以下同様。）は、基底関数［ａ１↑，ａ２↑，・・・，ａＭ↑］の線形和で表され、その結合係数ｓは統計的に独立である。これを数式で表すと次のようになる。
【００１８】
【数５】
ここで、ｘｉ↑＝Σｓｉｊａｊ↑（ｊについての総和）はｉ番目のサンプル画像であり、行列Ｓの列ベクトルｓｉ↑はその画像の独立な結合係数であり、特徴ベクトルと呼ばれる。また、ａｊ↑はｊ番目の基底関数である。
【００１９】
（１）式を行列式として表すと、
Ｘ↑＝Ａ↑・Ｓ↑ （２）
となる。
サンプル画像セット行列Ｘ↑のみから、基底関数と特徴ベクトルを求めることは、一種のブラインドソース分離問題とみなせるので、ＩＣＡ法を用いて解くことができる。
ＩＣＡでは、あるｎ×ｍの分離行列W↑を用いて、
Ｙ↑＝Ｗ↑・Ｘ↑ （３）
においてＹ↑を計算し、Ｙ↑の各成分が互いに統計的に独立となるようなＷ↑を求める。
Ｗ↑を求めるに当たり、Ｙ↑の各成分が独立か否かを判定する評価基準として、相互情報量（Kullback-Leibler divergence）または高次統計量（kurtosis）等が既に提案されている。この実施例では、Bell & Senjowskyが提案した相互情報量の最大化による手法を用いる。
【００２０】
今、分離信号［ｙ１，ｙ２，・・・，ｙｎ］が互いに独立であるとき、各成分ｙｉの同時分布密度関数ｐ（ｙ↑）＝Ｐ（ｙ１，ｙ２，・・・，ｙｎ）は、各成分ｙｉの周辺密度関数ｐ（ｙｉ）の積で表される。すなわち、
【００２１】
【数６】
となる。評価関数として、ｐ（ｙ↑）とΠｐ（ｙｉ）（ｉについての積）との距離を表すKullback-Leibler Divergence （Ｄ）が用いられる。Kullback-Leibler Divergence （Ｄ）は、エントロピーＨと分離行列Ｗ↑を用いて以下のように表される。
【００２２】
【数７】
ここで、
【００２３】
【数８】
Ｈ（ｙ↑）およびＨ（ｙｉ）は、それぞれ、同時分布のエントロピーと周辺分布のエントロピーである。（３）式より、Ｐ（Ｙ↑）＝Ｐ（Ｖ↑）／｜Ｗ↑｜であるので、同時分布エントロピーＨ（ｙ↑）は、次のように書き換えられる。
【００２４】
【数９】
よって、（５）式のKullback-Leibler
divergence （Ｄ）は次のように書き換えられる。
【００２５】
【数１０】
ＩＣＡ法では、ｐ（ｙ↑）とΠｐ（ｙｉ）（ｉについての積）が一致するように分離行列Ｗ↑を求めるが、これは、（９）式のＤを最小にするＷ↑を求める問題となる。
一般に、Ｗ↑を解析的に求めることはできず、（９）式のＤを最小にするＷ↑を求めるには、下記の自然勾配法が用いられる。
Ｗ↑の変化分は、（７）式で与えられる。通常の最急降下法に比べると、収束がよく、計算も簡単である。
【００２６】
【数１１】
こうして、Ｗ↑の更新則は次式で与えられる。
【００２７】
【数１２】
ここで、μは学習率であり、Ｉ↑は単位行列であり、アンサンブル平均は時間平均に置き換えている。なお、ｐ（ｙｉ）が未知であるから、ｇ（ｙ↑）は近似で求めることになる。この実施例では、ｇ（ｙ↑）はシグモイド関数である。（１２）式により、Ｙが独立となるようなＷ↑が求まり、また、（３）式よりＹ↑が求まる。
【００２８】
こうして、この実施例では、対象となるＩＶＲ画像とは別の４枚のＩＶＲ画像を用い、ＩＣＡ法によって基底関数の組を抽出した。
まず、これら４枚のＩＶＲ画像（各々１２８画素×１２８画素のサイズを有する）から任意に８画素×８画素の大きさのパッチ１２,０００個を切り出し、基底関数の組を学習によって抽出するためのサンプル画像とした。この場合、サンプル画像を形成するパッチのサイズは特に限定されないが、抽出される基底関数の組が、画質を改善すべきＩＶＲ画像の空間的特徴を最もよく反映し得るようなパッチのサイズを設定することが好ましい。
【００２９】
そして、各サンプル画像を６４×１の１次元ベクトルに展開し、行列Ｘ↑の列ベクトルとする。それによって、６４×１２,０００のサンプル画像セット行列Ｘ↑が得られる。行列Ｘ↑の行の数はサンプル画像のパッチの個数を表し、列の数はサンプル画像の数を表す。
（１２）式を用いてＷ↑を学習すると、６４×６４の行列Ｗ↑が求められ、このときのＷ↑の各列が基底関数となる。この実施例では、基底関数の数は６４となる。
【００３０】
再び図２を参照して、次に、ＩＶＲ画像が、上で求めた基底関数の組を用いてＩＣＡ変換され、基底関数の線形和の形に成分分解される（図２のＳ２２、Ｓ２３）。これを数式で表すと次のようになる。
【００３１】
【数１３】
【００３２】
ここで、Ｐ↑はＩＶＲ画像であり、［ａ０↑，ａ１↑，・・・，ａＮ−１↑］は基底関数であり、［ｓ０，ｓ１，・・・，ｓＮ−１］はＩＣＡ変換の変換係数である。
【００３３】
そして、ＩＶＲ画像をＩＣＡ変換したときの各変換係数が、Shrinkage関数によってフィルタ処理され、ＩＶＲ画像に含まれるノイズ成分が除去される（図２のＳ２４）。この実施例では、Shrinkage関数として、
【００３４】
【数１４】
ここで、δ２は、ポアッソンノイズの累乗であり、
【数１５】
【００３５】
この関数ｇ（ｓ↑）は、ポアッソンノイズに対するクロスバリデーション（cross-validation ）法に基づくShrinkage関数であり、ポアッソンノイズの除去に最適である。なお、Shrinkage関数の関数形はこれに限定されるものではなく、例えば、次式で表されるような公知の適当な関数形を用いることができる。
【数１６】
こうして、ノイズ分散δ２の増加に伴って、係数ｓの縮小される割合が大きくなる。
【００３６】
さらに、フィルタ処理後の各変換係数を用いてＩＣＡ逆変換が行われ、ノイズ成分の除去されたＩＶＲ画像（第１のＩＶＲ画像）が得られる（図２のＳ２５、Ｓ２６）。これを式で表すと次のようになる。
【００３７】
【数１７】
【００３８】
ここで、Ｐ（ハット）は第１のＩＶＲ画像であり、［ａ０↑，ａ１↑，・・・，ａＮ−１↑］は基底関数であり、［ｓ０（ハット），ｓ１（ハット），・・・，ｓＮ−１（ハット）］はフィルタ処理後の変換係数である。
【００３９】
図５は、この実施例で用いた元のＩＶＲ画像を示したものであり、図６は、図５のＩＶＲ画像に対しＩＣＡ・Shrinkageフィルタ処理を行って得られた第１のＩＶＲ画像を示したものである。図５および図６のＩＶＲ画像の比較から、ＩＣＡ・Shrinkageフィルタ処理によって、画像中の縞模様のような構造的なノイズがかなり除去されており、処理前のＩＶＲ画像を劣化させることなく、ノイズ成分を除去できることがわかる。
【００４０】
なお、上記実施例では、対象ＩＶＲ画像をＩＣＡ変換した後、変換係数をShrinkage関数によってフィルタ処理することでノイズ成分を除去するようにしたが、このフィルタ処理を行うとともに、ＩＶＲ画像における分解された画像部分のうち、目視によりノイズを含んでいると認識された画像部分を除去するようにしてもよい。
【００４１】
ＩＣＡ・Shrinkageフィルタ処理によれば、複数のサンプル画像から独立成分分析（ＩＣＡ）によって学習、抽出した基底関数に基づき、元のＩＶＲ画像をＩＣＡ変換し、変換係数をShrinkage関数によってフィルタ処理することでノイズ成分を除去した後、ＩＣＡ逆変換によって第１のＩＶＲ画像を得るようにしたので、元の画像の信号成分を失うことなくノイズ成分を除去することができる。
【００４２】
再び図１を参照して、本発明によれば、次に、第１のＩＶＲ画像に対し多重スケール（Multi-scale）フィルタ処理がなされ（図１のＳ４）、その画像中の線状陰影の信号を強調した第２のＩＶＲ画像が生成される（図１のＳ５）。図３には、多重スケールフィルタ処理のフロー図が示してある。
【００４３】
図３を参照して、多重スケールフィルタ処理において、まず、強調すべき線状陰影の幅に対応するスケールｓをもつガウス関数Ｇ（ｘ，ｙ）が生成される（図３のＳ３０）。すなわち、ガウス関数Ｇのスケールｓの大きさに応じて、ＩＶＲ画像中の強調される線状の陰影部分の幅が決定される。
【００４４】
次に、Ｇ（ｘ，ｙ）のx方向の微分∂Ｇ（ｘ，ｙ）／∂ｘ（＝Ｇｘ（ｘ，ｙ））およびｙ方向の微分∂Ｇ（ｘ，ｙ）／∂ｙ（＝Ｇｙ（ｘ，ｙ））が求められ、さらに、Ｇｘ（ｘ，ｙ）のｘ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｘ２（＝Ｇｘｘ（ｘ，ｙ））およびｙ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｘ∂ｙ（＝Ｇｘｙ（ｘ，ｙ））と、Ｇｙ（ｘ，ｙ）のｘ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｙ∂ｘ（＝Ｇｙｘ（ｘ，ｙ））およびｙ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｙ２（＝Ｇｙｙ（ｘ，ｙ））が求められる（図３のＳ３１）。
【００４５】
そして、第１のＩＶＲ画像の各画素値Ｌ（ｘ，ｙ）の微分∂２Ｌ（ｘ，ｙ）／∂ｘ２（＝Ｌｘｘ（ｘ，ｙ））、∂２Ｌ（ｘ，ｙ）／∂ｘ∂ｙ（＝Ｌｘｙ（ｘ，ｙ））、∂２Ｌ（ｘ，ｙ）／∂ｙ∂ｘ（＝Ｌｙｘ（ｘ，ｙ））および∂２Ｌ（ｘ，ｙ）／∂ｙ２（＝Ｌｙｙ（ｘ，ｙ））が、
Ｌｘｘ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｘｘ（ｘ，ｙ）、
Ｌｘｙ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｘｙ（ｘ，ｙ）、
Ｌｙｘ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｙｘ（ｘ，ｙ）、
Ｌｙｙ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｙｙ（ｘ，ｙ）、
ここで、γは任意の定数である、
によって求められる（図３のＳ３２）。
【００４６】
さらに、第１のＩＶＲ画像の各画素の位置におけるヘッセ行列Ｈが
【数１８】
によって求められ（図３のＳ３３）、ヘッセ行列Ｈの２つの固有値λ１、λ２（｜λ１｜＜｜λ２｜）および対応する固有ベクトルν１、ν２が求められる（図３のＳ３４）。この場合、λ１の固有ベクトルν１は、関係する線状陰影の軸方向を表し、λ２の固有ベクトルν２は、当該線状陰影の軸に対し垂直な方向を表している。
【００４７】
そして、
【数１９】
が計算され（図３のＳ３５）、ＩＶＲ画像の各画素値が、
【数２０】
ここで、ｃおよびβはいずれも任意の定数である、
によって求められる（図３のＳ３５）。この式において、βの値を大きくすればするほど画像は明るくなり、ｃの値を大きくすればするほど画像は暗くなる。
【００４８】
次に、ガウス関数Ｇ（ｘ，ｙ）のスケールｓが変更され、図３のステップＳ３０〜Ｓ３５が繰り返される。そして、すべてのｓについて計算が完了したとき、得られたＩｓ（ｘ，ｙ）を用いて、
【数２１】
に従って、ＩＶＲ画像の各画素値が求められる（図３のＳ３７）。この場合、計算に使用されるスケールｓの値の組は、対象となるＩＶＲ画像中に現われた血管、ガイドワイヤおよびカテーテル等の線状の陰影部分の各幅に対応して、予め決定される。
【００４９】
再び図１を参照して、本発明の方法によれば、さらに、得られた第２のＩＶＲ画像のコントラストが反転せしめられ、その画像と第１のＩＶＲ画像が合成され（図１のＳ６）、画質の改善されたＩＶＲ画像が生成される（図１のＳ７）。なお、第１および第２のＩＶＲ画像の合成に際し、必要に応じて、各画像に重み係数が付与される。
【００５０】
図７は、図６のＩＶＲ画像（第１のＩＶＲ画像）に対し多重スケールフィルタ処理を行って得られた第２のＩＶＲ画像を示したものである。また、図８は、図７のＩＶＲ画像（第２のＩＶＲ画像）のコントラストを反転させた画像と、図６のＩＶＲ画像（第１のＩＶＲ画像）を合成して得られた、画質の改善されたＩＶＲ画像を示したものである。この場合、第１のＩＶＲ画像の割合を１とし、コントラスト反転した第２のＩＶＲ画像の割合を０．１として両画像を合成し、得られた画像を１．１で割っている。
図５および図８のＩＶＲ画像の比較から、本発明による画質改善法によれば、ＩＶＲ画像におけるノイズ除去と線状陰影の強調を同時に達成することができることがわかる。
【００５１】
上記実施例では、ＩＶＲ画像が静止画像である場合について説明してきたが、本発明の方法は、ＩＶＲ画像が動画像である場合にも適用可能である。この場合には、動画像を１フレーム毎の画像に分解し、分離した画像毎に上記と同様の処理を行い、処理後の各画像を再合成すればよい。
【００５２】
図４は、本発明の別の実施例によるＩＶＲ画像の画質改善法のフロー図である。図４を参照して、この実施例では、画質改善すべき同一の元のＩＶＲ画像が２つ用意される（図４のＳ１０）。そして、一方の元のＩＶＲ画像に対しＩＣＡ・Shrinkageフィルタ処理がなされ、ノイズ信号が除去された第１のＩＶＲ画像が生成される一方（図４のＳ１１、Ｓ１２）、他方の元のＩＶＲ画像に対し多重スケールフィルタ処理がなされ、その画像中の線状陰影の信号が強調された第２のＩＶＲ画像が生成される（図４のＳ１３、Ｓ１４）。
【００５３】
そして、第２のＩＶＲ画像のコントラストが反転せしめられ（図４のＳ１５）、その画像と第１のＩＶＲ画像とが合成され（図４のＳ１６）、画質の改善されたＩＶＲ画像が生成される（図４のＳ１７）。なお、第１および第２のＩＶＲ画像の合成の際に、必要に応じて、各画像に重み係数が付与され得る。
この実施例の場合にも、図１の実施例の場合と同様の効果が得られる。
【図面の簡単な説明】
【００５４】
【図１】本発明の１実施例によるＩＶＲ画像の画質改善法のフロー図である。
【図２】図１中のＩＣＡ・Shrinkageフィルタ処理のフロー図である。
【図３】図１中の多重スケールフィルタ処理のフロー図である。
【図４】本発明の別の実施例によるＩＶＲ画像の画像改善法のフロー図である。
【図５】画質を改善すべきＩＶＲ画像を示したものである。
【図６】図５のＩＶＲ画像に対しＩＣＡ・Shrinkageフィルタ処理を行って得られた第１のＩＶＲ画像を示したものである。
【図７】図６のＩＶＲ画像に対し多重スケールフィルタ処理を行って得られた第２のＩＶＲ画像を示したものである。
【図８】図７のＩＶＲ画像（第２のＩＶＲ画像）のコントラストを反転させ、その画像と図６のＩＶＲ画像（第１のＩＶＲ画像）を合成して得られた、画質の改善されたＩＶＲ画像を示したものである。
【特許請求の範囲】
【請求項１】
（ａ）画質改善すべき元のＩＶＲ画像に対しＩＣＡ・Shrinkageフィルタ処理を行い、ノイズ信号を除去した第１のＩＶＲ画像を生成するステップと、
（ｂ）前記第１のＩＶＲ画像に対し多重スケールフィルタ処理を行い、その画像中の線状陰影の信号を強調した第２のＩＶＲ画像を生成するステップと、
（ｃ）前記第２のＩＶＲ画像のコントラストを反転させた後、前記第１および第２のＩＶＲ画像を合成して、画質の改善されたＩＶＲ画像を生成するステップと、を有することを特徴とするＩＶＲ画像の画質改善法。
【請求項２】
前記ステップ（ｃ）において前記第１および第２の画像を合成するとき、それぞれの画像に重み係数を付与することを特徴とする請求項１に記載のＩＶＲ画像の画質改善法。
【請求項３】
（ａ’）画質改善すべき同一の元のＩＶＲ画像を２つ用意し、
（ｂ’）一方の前記元のＩＶＲ画像に対しＩＣＡ・Shrinkageフィルタ処理を行い、ノイズ信号を除去した第１のＩＶＲ画像を生成するとともに、他方の前記元のＩＶＲ画像に対し多重スケールフィルタ処理を行い、その画像中の線状陰影の信号を強調した第２のＩＶＲ画像を生成するステップと、
（ｃ’）前記第２のＩＶＲ画像のコントラストを反転させた後、前記第１および第２のＩＶＲ画像を合成して画質の改善されたＩＶＲ画像を生成するステップと、を有することを特徴とするＩＶＲ画像の画質改善法。
【請求項４】
前記ステップ（ｃ’）において前記第１および第２の画像を合成するとき、それぞれの画像に重み係数を付与することを特徴とする請求項３に記載のＩＶＲ画像の画質改善法。
【請求項５】
前記ＩＣＡ・Shrinkageフィルタ処理は、
（i）予め与えられた複数のサンプル画像から、独立成分分析（ＩＣＡ）によって、前記サンプル画像の特徴を統計的に表す基底関数の組を抽出するステップと、
（ii）前記ＩＶＲ画像をＩＣＡ変換して、前記ＩＶＲ画像を前記基底関数の線形和の形に成分分解するステップと、
（iii）前記ＩＶＲ画像を前記ＩＣＡ変換したときの各変換係数を、Shrinkage関数によってフィルタ処理して前記ＩＶＲ画像に含まれるノイズ成分を除去するステップと、
（iv）前記フィルタ処理した後の各変換係数を用いてＩＣＡ逆変換を行い、ノイズの除去されたＩＶＲ画像を得るステップと、からなることを特徴とする請求項１〜請求項４のいずれかに記載のＩＶＲ画像の画質改善法。
【請求項６】
前記ステップ（iii）において、前記Shrinkage関数として、ポアッソンノイズに対するクロスバリデーション（cross-validation）法に基づくShrinkage関数を用いることを特徴とする請求項５に記載のＩＶＲ画像の画質改善法。
【請求項７】
前記ステップ（i）において、前記サンプル画像を、前記ＩＶＲ画像に類似する対象物について前記ＩＶＲ画像の取得時と近似した条件下で取得したデジタル画像に基づいて作成することを特徴とする請求項５または請求項６に記載の方法。
【請求項８】
前記多重スケールフィルタ処理は、
（i）強調すべき線状陰影の幅に対応するスケールｓをもつガウス関数Ｇ（ｘ，ｙ）を生成するステップと、
（ii）Ｇ（ｘ，ｙ）のx方向の微分∂Ｇ（ｘ，ｙ）／∂ｘ（＝Ｇｘ（ｘ，ｙ））およびｙ方向の微分∂Ｇ（ｘ，ｙ）／∂ｙ（＝Ｇｙ（ｘ，ｙ））を求め、さらに、Ｇｘ（ｘ，ｙ）のｘ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｘ２（＝Ｇｘｘ（ｘ，ｙ））およびｙ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｘ∂ｙ（＝Ｇｘｙ（ｘ，ｙ））と、Ｇｙ（ｘ，ｙ）のｘ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｙ∂ｘ（＝Ｇｙｘ（ｘ，ｙ））およびｙ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｙ２（＝Ｇｙｙ（ｘ，ｙ））とを求めるステップと、
（iii）前記ＩＶＲ画像の各画素値Ｌ（ｘ，ｙ）の微分∂２Ｌ（ｘ，ｙ）／∂ｘ２（＝Ｌｘｘ（ｘ，ｙ））、∂２Ｌ（ｘ，ｙ）／∂ｘ∂ｙ（＝Ｌｘｙ（ｘ，ｙ））、∂２Ｌ（ｘ，ｙ）／∂ｙ∂ｘ（＝Ｌｙｘ（ｘ，ｙ））および∂２Ｌ（ｘ，ｙ）／∂ｙ２（＝Ｌｙｙ（ｘ，ｙ））を、
Ｌｘｘ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｘｘ（ｘ，ｙ）、
Ｌｘｙ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｘｙ（ｘ，ｙ）、
Ｌｙｘ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｙｘ（ｘ，ｙ）、
Ｌｙｙ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｙｙ（ｘ，ｙ）、
ここで、γは任意の定数である、
によって求めるステップと、
（iv）前記ＩＶＲ画像の各画素の位置におけるヘッセ行列Ｈを
【数１】
により求めるステップと、
（v）前記ヘッセ行列Ｈの２つの固有値λ１、λ２（｜λ１｜＜｜λ２｜）および対応する固有ベクトルν１、ν２を求め、
【数２】
を求めるステップと、
（vi）前記ＩＶＲ画像の各画素値を、
【数３】
ここで、ｃおよびβは任意の定数である、
によって求めるステップと、
（vii）前記ガウス関数Ｇ（ｘ，ｙ）のスケールｓの値を変更して、前記ステップ（ii）〜（vi）を繰り返すステップと、
（viii）前記ステップ（ii）〜（vii）を所定の回数繰り返した後、得られたＩｓ（ｘ，ｙ）を用いて、
【数４】
に従って、前記ＩＶＲ画像の各画素値を求めるステップと、からなっていることを特徴とする請求項１〜請求項７のいずれかに記載のＩＶＲ画像の画質改善法。
【請求項１】
（ａ）画質改善すべき元のＩＶＲ画像に対しＩＣＡ・Shrinkageフィルタ処理を行い、ノイズ信号を除去した第１のＩＶＲ画像を生成するステップと、
（ｂ）前記第１のＩＶＲ画像に対し多重スケールフィルタ処理を行い、その画像中の線状陰影の信号を強調した第２のＩＶＲ画像を生成するステップと、
（ｃ）前記第２のＩＶＲ画像のコントラストを反転させた後、前記第１および第２のＩＶＲ画像を合成して、画質の改善されたＩＶＲ画像を生成するステップと、を有することを特徴とするＩＶＲ画像の画質改善法。
【請求項２】
前記ステップ（ｃ）において前記第１および第２の画像を合成するとき、それぞれの画像に重み係数を付与することを特徴とする請求項１に記載のＩＶＲ画像の画質改善法。
【請求項３】
（ａ’）画質改善すべき同一の元のＩＶＲ画像を２つ用意し、
（ｂ’）一方の前記元のＩＶＲ画像に対しＩＣＡ・Shrinkageフィルタ処理を行い、ノイズ信号を除去した第１のＩＶＲ画像を生成するとともに、他方の前記元のＩＶＲ画像に対し多重スケールフィルタ処理を行い、その画像中の線状陰影の信号を強調した第２のＩＶＲ画像を生成するステップと、
（ｃ’）前記第２のＩＶＲ画像のコントラストを反転させた後、前記第１および第２のＩＶＲ画像を合成して画質の改善されたＩＶＲ画像を生成するステップと、を有することを特徴とするＩＶＲ画像の画質改善法。
【請求項４】
前記ステップ（ｃ’）において前記第１および第２の画像を合成するとき、それぞれの画像に重み係数を付与することを特徴とする請求項３に記載のＩＶＲ画像の画質改善法。
【請求項５】
前記ＩＣＡ・Shrinkageフィルタ処理は、
（i）予め与えられた複数のサンプル画像から、独立成分分析（ＩＣＡ）によって、前記サンプル画像の特徴を統計的に表す基底関数の組を抽出するステップと、
（ii）前記ＩＶＲ画像をＩＣＡ変換して、前記ＩＶＲ画像を前記基底関数の線形和の形に成分分解するステップと、
（iii）前記ＩＶＲ画像を前記ＩＣＡ変換したときの各変換係数を、Shrinkage関数によってフィルタ処理して前記ＩＶＲ画像に含まれるノイズ成分を除去するステップと、
（iv）前記フィルタ処理した後の各変換係数を用いてＩＣＡ逆変換を行い、ノイズの除去されたＩＶＲ画像を得るステップと、からなることを特徴とする請求項１〜請求項４のいずれかに記載のＩＶＲ画像の画質改善法。
【請求項６】
前記ステップ（iii）において、前記Shrinkage関数として、ポアッソンノイズに対するクロスバリデーション（cross-validation）法に基づくShrinkage関数を用いることを特徴とする請求項５に記載のＩＶＲ画像の画質改善法。
【請求項７】
前記ステップ（i）において、前記サンプル画像を、前記ＩＶＲ画像に類似する対象物について前記ＩＶＲ画像の取得時と近似した条件下で取得したデジタル画像に基づいて作成することを特徴とする請求項５または請求項６に記載の方法。
【請求項８】
前記多重スケールフィルタ処理は、
（i）強調すべき線状陰影の幅に対応するスケールｓをもつガウス関数Ｇ（ｘ，ｙ）を生成するステップと、
（ii）Ｇ（ｘ，ｙ）のx方向の微分∂Ｇ（ｘ，ｙ）／∂ｘ（＝Ｇｘ（ｘ，ｙ））およびｙ方向の微分∂Ｇ（ｘ，ｙ）／∂ｙ（＝Ｇｙ（ｘ，ｙ））を求め、さらに、Ｇｘ（ｘ，ｙ）のｘ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｘ２（＝Ｇｘｘ（ｘ，ｙ））およびｙ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｘ∂ｙ（＝Ｇｘｙ（ｘ，ｙ））と、Ｇｙ（ｘ，ｙ）のｘ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｙ∂ｘ（＝Ｇｙｘ（ｘ，ｙ））およびｙ方向の微分∂２Ｇ（ｘ，ｙ）／∂ｙ２（＝Ｇｙｙ（ｘ，ｙ））とを求めるステップと、
（iii）前記ＩＶＲ画像の各画素値Ｌ（ｘ，ｙ）の微分∂２Ｌ（ｘ，ｙ）／∂ｘ２（＝Ｌｘｘ（ｘ，ｙ））、∂２Ｌ（ｘ，ｙ）／∂ｘ∂ｙ（＝Ｌｘｙ（ｘ，ｙ））、∂２Ｌ（ｘ，ｙ）／∂ｙ∂ｘ（＝Ｌｙｘ（ｘ，ｙ））および∂２Ｌ（ｘ，ｙ）／∂ｙ２（＝Ｌｙｙ（ｘ，ｙ））を、
Ｌｘｘ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｘｘ（ｘ，ｙ）、
Ｌｘｙ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｘｙ（ｘ，ｙ）、
Ｌｙｘ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｙｘ（ｘ，ｙ）、
Ｌｙｙ（ｘ，ｙ）＝ｓγＬ（ｘ，ｙ）×Ｇｙｙ（ｘ，ｙ）、
ここで、γは任意の定数である、
によって求めるステップと、
（iv）前記ＩＶＲ画像の各画素の位置におけるヘッセ行列Ｈを
【数１】
により求めるステップと、
（v）前記ヘッセ行列Ｈの２つの固有値λ１、λ２（｜λ１｜＜｜λ２｜）および対応する固有ベクトルν１、ν２を求め、
【数２】
を求めるステップと、
（vi）前記ＩＶＲ画像の各画素値を、
【数３】
ここで、ｃおよびβは任意の定数である、
によって求めるステップと、
（vii）前記ガウス関数Ｇ（ｘ，ｙ）のスケールｓの値を変更して、前記ステップ（ii）〜（vi）を繰り返すステップと、
（viii）前記ステップ（ii）〜（vii）を所定の回数繰り返した後、得られたＩｓ（ｘ，ｙ）を用いて、
【数４】
に従って、前記ＩＶＲ画像の各画素値を求めるステップと、からなっていることを特徴とする請求項１〜請求項７のいずれかに記載のＩＶＲ画像の画質改善法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【公開番号】特開２００９−２８３１９（Ｐ２００９−２８３１９Ａ）
【公開日】平成２１年２月１２日（２００９．２．１２）
【国際特許分類】
【出願番号】特願２００７−１９５６９３（Ｐ２００７−１９５６９３）
【出願日】平成１９年７月２７日（２００７．７．２７）
【出願人】（５９３００６６３０）学校法人立命館 (359)
【出願人】（０００００１９９３）株式会社島津製作所 (3,708)
【Ｆターム（参考）】
【公開日】平成２１年２月１２日（２００９．２．１２）
【国際特許分類】
【出願日】平成１９年７月２７日（２００７．７．２７）
【出願人】（５９３００６６３０）学校法人立命館 (359)
【出願人】（０００００１９９３）株式会社島津製作所 (3,708)
【Ｆターム（参考）】
[ Back to top ]