(57)【特許請求の範囲】
【請求項１】
画像に含まれる線状もしくはエッジ状の陰影の方向とその強さを検出する手段と、
前記検出手段により検出された方向と強さに応じた指向性を有する空間フィルタを前記画像に施すフィルタ手段と、を具備し、
前記画像の陰影の画素レベルが背景の画素より大きい場合、前記検出手段は画像の関心領域の共分散行列の２つの固有値のうち、大きい方の固有値に対応する固有ベクトルの向きに基づいて方向を算出し、
前記フィルタ手段は、該方向の周波数特性がその他の方向の周波数特性と異なるディジタル空間フィルタを前記画像に施すことを特徴とする画像処理装置。
【請求項２】
画像に含まれる線状もしくはエッジ状の陰影の方向とその強さを検出する手段と、
前記検出手段により検出された方向と強さに応じた指向性を有する空間フィルタを前記画像に施すフィルタ手段と、を具備し、
前記画像の陰影の画素レベルが背景の画素より小さい場合、前記検出手段は画像の関心領域の共分散行列の２つの固有値のうち、小さい方の固有値に対応する固有ベクトルの向きに基づいて方向を算出し、
前記フィルタ手段は、該方向の周波数特性がその他の方向の周波数特性と異なるディジタル空間フィルタを前記画像に施すことを特徴とする画像処理装置。
【請求項３】
前記検出手段は、関心領域を水平方向と垂直方向にシフトしながら画像全体に対する関心領域のサイズの局所毎の陰影の方向とその強さを逐次求め、
前記フィルタ手段は、方向と強さによって決定される２次元空間フィルタ係数を関心領域のサイズ毎に変更しながら画像全体に空間フィルタを施すことを特徴とする請求項１、または請求項２に記載の画像処理装置。
【請求項４】
前記検出手段は、関心領域内の画素配列に対して、各行毎に所定の行方向係数と画素値の間の積和演算を行い、求められた各行の演算結果に対して所定の列方向係数との間で積和演算を行った結果の数値を、複数の行方向係数と列方向係数の組み合わせに対して少なくとも５個算出し、その数値を基にして所定の演算を行って関心領域の共分散行列を求めることを特徴とする請求項１、または請求項２に記載の画像処理装置。
【請求項５】
画像に含まれる線状もしくはエッジ状の陰影の方向とその強さを検出する手段と、
前記検出手段により検出された方向と強さに応じた指向性を有する空間フィルタを前記画像に施すフィルタ手段と、を具備し、
前記検出手段は、関心領域内の画素分布の平均値に基づいて関心領域内のランダムノイズの分散を推定し、その分散と関心領域のマトリクスサイズで定まる確率分布を強さを表す指標について予測しておき、関心領域内の画素値を基に計算された強さが予測された確率分布の下で取り得る確率を計算し、この確率を強さの指標とすることを特徴とする画像処理装置。
【請求項６】
前記検出手段は、共分散行列の２つの固有値の比を方向の強さの指標とすることを特徴とする請求項１、または請求項２に記載の画像処理装置。
【請求項７】
画像に含まれる線状もしくはエッジ状の陰影の方向とその強さを検出する手段と、
前記検出手段により検出された方向と強さに応じた指向性を有する空間フィルタを前記画像に施すフィルタ手段と、を具備し、
前記検出手段は、設定した関心領域の共分散行列の要素と、関心領域の周辺の領域の共分散行列の要素とを重み付け平均して新しく共分散行列を求め、この新しく求められた共分散行列の固有ベクトルの方向と固有値とを、それぞれ関心領域のパターン走行方向とパターン走行方向の強さの指標として、パターン方向ベクトルを決め、この方向ベクトルによって決定される空間フィルタを施すことを特徴とする画像処理装置。
【請求項８】
前記検出手段は、設定した関心領域と、その周辺の複数の領域について新しく共分散行列を求め、複数の領域についての固有ベクトルの方向を求め、最も多くの領域で求められた方向を関心領域のパターン走行方向として求めることを特徴とする請求項７に記載の画像処理装置。
【請求項９】
前記検出手段は、入力画像に対して周波数低減フィルタを施し、その出力画像と入力画像との差画像に対して方向検出処理を行うことを特徴とする請求項１、または請求項２に記載の画像処理装置。
【請求項１０】
前記検出手段は、入力画像の画素値の符号を反転して方向検出処理を行う請求項１、または請求項２に記載の画像処理装置。
【請求項１１】
前記検出手段は、入力画像の最大ビット数で決定される最大値から入力画像を減じた画像に対して方向検出処理を行うことを特徴とする請求項１、または請求項２に記載の画像処理装置。
【請求項１２】
前記検出手段は、入力画像の信号スペクトルのある帯域を強調しつつ、ノイズスペクトルを抑えるプリプロセス処理を行うことを特徴とする請求項１、または請求項２のいずれか一項に記載の画像処理装置。
【請求項１３】
予測しておく強さを表す指標の分布は、画像の関心領域の中心位置と、関心領域内の画素分布の重心位置に基づいて決定された位置のまわりの共分散行列から求められる固有値に対して行われることを特徴とする請求項５に記載の画像処理装置。
【請求項１４】
予測しておく強さを表す指標の分布は、画像の関心領域の中心位置と、関心領域内の画素分布の重心位置に基づいて決定された位置のまわりの共分散行列から求められる固有値の差または比に対して行われることを特徴とする請求項５に記載の画像処理装置。
【請求項１５】
直接変換方式のＸ線検出器を用いて被検体のＸ線画像を得る手段と、
前記Ｘ線画像が入力される請求項１乃至請求項１４のいずれか一項に記載の画像処理装置と、
を具備することを特徴とするＸ線診断装置。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は画像処理装置に関し、特に、医用画像、例えばＸ線画像上に描写された血管影やガイドワイヤなどの線状陰影やエッジ状の陰影を含む画像において、陰影と背景を正確に分離し、陰影のコントラストを損なわずに、背景のノイズ成分を選択的に抑圧する画像処理装置に関するものである。
【０００２】
【従来の技術】
Ｘ線を利用した医療技術として、透視下におけるカテーテル治療が盛んに行われている。しかし、Ｘ線透視は患者の被曝量を低減するために、Ｘ線撮影に比べてＸ線の線量を低減しているので、画像上に重畳するノイズが大きく、ガイドワイヤやカテーテル、もしくは造影された血管影などが背景ノイズに邪魔されて見にくいという問題がある。また、線量を増加して透視を行うと、背景ノイズは相対的に小さくなるが、患者の被曝や医療技術者の被曝量が大きくなり、現実的ではない。そのため、画像処理によってノイズを低減させたり、コントラストを向上させる技術が現在に至っても研究され続けている。
【０００３】
最も一般に用いられるノイズ低減のための画像処理技術は、時間的に連続して入力される画像の加算平均をとる手法である。この技術はすでに公知であり、その応用技術としてのリカーシブフィルタと呼ばれる回帰フィルタも広く一般に使用されている。しかし、時間的な加算は動きの無い領域に対してはノイズ低減作用があるが、心臓血管に挿入されたガイドワイヤのような動きのある物体に対してはそのコントラストを薄めたり、残像を生じたりする欠点がある。そのため、被写体の動きを検出して画像加算効果が変わるように加算平均の重みを画像毎に変更することが行われている（例えば，特許Ｐ２５０８０７８、特公平６−６９４４７号公報、特開平３−１９８８３６号公報、特開平６−４７０３５号公報、特開平７−７９９５６号公報、特開平８−２５５２３８号公報参照）。
【０００４】
しかしながら、被写体の動きの判断は基本的に現在の画像とそれまでの画像の差分に基づいた手法で行われるため、被写体の動きと同時に背景ノイズによる画素濃度の時間変更も動きと誤認される。
【０００５】
したがって、医療の進歩とともに細径化されているカテーテルやコントラストの小さいガイドワイヤに対して、それらの動きを検出するために検出精度を増加させる場合には、ノイズ低減の効果に限界が生じている。ノイズ低減効果を大きくするような処理パラメータを設定すると、残像が生じてかえって見にくい画像になる。
【０００６】
一方、ノイズ低減やコントラストを向上させる目的で、線形空間フィルタも一般に使用されている。これは空間的に分布する画素の重み付き平均を得て、スムージングさせたり、空間的な微分作用を持たせることによってコントラストを強調する処理である。しかし、スムージングによってノイズだけでなくガイドワイヤやカテーテルのような線状陰影がぼやけてコントラストが損なわれる欠点や、コントラスト強調の場合にはノイズも強調されたりする欠点がある。このため、空間的に血管影をパターン認識する装置（特開平４−１２２３５５号公報）、あるいは被写体のエッジ領域を検出してエッジ領域にはコントラスト強調を施し、それ以外の領域にはスムージングを施す装置（特許ＰＳ６０２４５０８４）などが開発された。前者は１画素単位に多くの平均化回路や比較回路を必要とするため、透視処理などのリアルタイム演算を行うには非常に多くの回路を必要とするうえ、ノイズを低減させる効果はない。また、後者は単に１ビットの情報によってエッジ領域と背景領域を区別するため、誤認識があった場合には、極端にスムージングを施された部分に隣接して誤認識されたノイズがコントラスト強調されて出力されることなどがあり、誤認識が敏感に出力画像のアーチファクトとなって画質を損ないやすい。
【０００７】
また、特許ＰＨ０９２４８９２０には、画像の線状陰影を２ビットの情報によって線状陰影と背景を識別する画像処理方法の記載がある。これにより、識別精度は向上するものの誤認識によって多少のアーチファクトが生じるという問題がある。また、線状陰影そのものに対してノイズ低減を行うことが有効にできないとう問題がある。
【０００８】
一方、その他の方式の画像処理方法として、方向性を求めてその方向成分を強調したり、あるいは方向成分を有する領域とそれ以外の領域に分けて別々のフィルタを施す手法が報告されている。これらはDirectional  Adaptive  filterと称され、いくつかの技術が公知である。
【０００９】
ＵＳＰ第５，７９９，１００号明細書、またはＵＳＰ第５，７１７，７９１号明細書には、Ｘ線画像に対してウェーブレット変換を用いて方向を検出する方式についての記述がある。また、ＭＲＩ画像に対して４方向のフィルタをかけて原画との差を分析することにより、エッジのない場所ではノイズを低減させる作用の強いフィルタを施す手法がＵＳＰ第５，６０２，９３４号明細書に記載されている。しかしながら、これらの手法はリアルタイムに処理を行うにはアルゴリズムや演算が複雑である。特にウェーブレット変換を用いる手法はその変換を行うだけでもリアルタイム処理を行うことが困難になる。また、ＵＳＰ第５，６０２，９３４号明細書に記載された手法は４方向の角度に指向性を持つフィルタを実際に施すが、この角度以外の角度で線状陰影やエッジが存在した場合の精度に難点があると同時に、フィルタ演算回数や画像間演算の量が多く、さらなる高速性、精度向上が望まれる。
【００１０】
このように、従来の画像処理技術では、線状陰影もしくはエッジ状の陰影を含む画像において、線状陰影やエッジ状陰影のコントラストを損なわずに、背景のノイズ成分を選択的に、しかもリアルタイムで抑圧することが困難であった。なお、この問題は、Ｘ線画像に限らず、ＣＴ画像やＭＲＩ画像等、ランダムノイズが重畳するような画像全般に対しても同様に生じ得る。
【００１１】
【発明が解決しようとする課題】
本発明の目的は、線状もしくはエッジ状の陰影と背景ノイズとの分離の精度を向上させ、陰影のコントラストを損なわずに背景ノイズを抑圧でき、しかもリアルタイム処理が可能な画像処理装置を提供することである。
【００１２】
【課題を解決するための手段】
前記課題を解決し目的を達成するために、本発明は以下に示す手段を用いている。
【００１３】
本発明の一態様による画像処理装置は、画像に含まれる線状もしくはエッジ状の陰影の方向とその強さを検出する手段と、前記検出手段により検出された方向と強さに応じた指向性を有する空間フィルタを前記画像に施すフィルタ手段と、を具備し、
前記画像の陰影の画素レベルが背景の画素より大きい場合、前記検出手段は画像の関心領域の共分散行列の２つの固有値のうち、大きい方の固有値に対応する固有ベクトルの向きに基づいて方向を算出し、
前記フィルタ手段は、該方向の周波数特性がその他の方向の周波数特性と異なるディジタル空間フィルタを前記画像に施すことを特徴とする。
【００１４】
本発明の他の態様による画像処理装置は、画像に含まれる線状もしくはエッジ状の陰影の方向とその強さを検出する手段と、前記検出手段により検出された方向と強さに応じた指向性を有する空間フィルタを前記画像に施すフィルタ手段と、を具備し、
前記画像の陰影の画素レベルが背景の画素より小さい場合、前記検出手段は画像の関心領域の共分散行列の２つの固有値のうち、小さい方の固有値に対応する固有ベクトルの向きに基づいて方向を算出し、
前記フィルタ手段は、該方向の周波数特性がその他の方向の周波数特性と異なるディジタル空間フィルタを前記画像に施すことを特徴とする。
【００１５】
本発明の他の態様による画像処理装置は、画像に含まれる線状もしくはエッジ状の陰影の方向とその強さを検出する手段と、前記検出手段により検出された方向と強さに応じた指向性を有する空間フィルタを前記画像に施すフィルタ手段と、を具備し、
前記検出手段は、関心領域内の画素分布の平均値に基づいて関心領域内のランダムノイズの分散を推定し、その分散と関心領域のマトリクスサイズで定まる確率分布を強さを表す指標について予測しておき、関心領域内の画素値を基に計算された強さが予測された確率分布の下で取り得る確率を計算し、この確率を強さの指標とすることを特徴とする。
【００１６】
本発明の他の態様による画像処理装置は、画像に含まれる線状もしくはエッジ状の陰影の方向とその強さを検出する手段と、前記検出手段により検出された方向と強さに応じた指向性を有する空間フィルタを前記画像に施すフィルタ手段と、を具備し、
前記検出手段は、設定した関心領域の共分散行列の要素と、関心領域の周辺の領域の共分散行列の要素とを重み付け平均して新しく共分散行列を求め、この新しく求められた共分散行列の固有ベクトルの方向と固有値とを、それぞれ関心領域のパターン走行方向とパターン走行方向の強さの指標として、パターン方向ベクトルを決め、この方向ベクトルによって決定される空間フィルタを施すことを特徴とする。
【００３６】
【発明の実施の形態】
以下、図面を参照して本発明による画像処理装置の実施形態を説明する。
【００３７】
第１実施形態
ここでは、処理対象画像としてはＸ線画像を例にとり説明する。一般に、ディジタルＸ線画像に含まれるノイズは１画素に入射するＸ線フォトンの個数Ｎの平方根に比例することが知られている。例えば、１画素に１００個のフォトン数が入射すれば、出力信号の期待値Ｌとノイズの標準偏差σは次のように表される。
【００３８】
Ｌ＝ａ×１００×Ｅ
ここで、ＥはＸ線フォトンの平均エネルギー、ａはＸ線エネルギーを出力信号に変換する際のゲインである。
【００３９】
σ＝ａ×１０×Ｅ
したがって均一な被写体にＸ線を照射して得られた画像は平均的にはＬの画像レベルの上にσの標準偏差を持つようなノイズが分布したものとなる。インターベンショナルラジオロジー（ＩＶＲ）で用いるような細いガイドワイヤやカテーテルはそのノイズの中で線状の陰影として描出されている。本発明の目的はこのようなノイズに埋もれた線状の陰影に対し、シャープネスを低減することなくノイズを低減することである。
【００４０】
なお、Ｘ線以外の医用画像においても、血管影、腫瘍陰影や臓器のエッジなどの線状陰影を鮮明に観測する要求があり、本発明はこれらの画像にも適用可能であるし、さらに、医用画像に限らず、画像全般にも適用可能である。
【００４１】
図１は本発明の第１実施形態に係る画像処理装置の概略構成を示す図である。本実施形態の基本的な手法は線状陰影の走行方向に沿ってノイズ低減フィルタを施すものである。入力画像がプリプロセス回路１０、空間フィルタ回路１２に入力される。プリプロセス回路１０は線状陰影、または画像のエッジを多少強調するか、もしくはノイズを全体的に低減して、走行方向を検知しやすくするためのものである。あるいは、Ｘ線画像では観測対象となるのはカテーテルやガイドワイヤなどの医用機具、もしくは造影血管であるが、通常、この陰影は背景に対して低い画素値である。これはＸ線阻止能の大きな物体、例えば比重が水に比べて大きい金属で構成されるためであるが、後段の処理においては、この陰影を逆に背景に対して大きな画素値とするほうが都合が良い。そのため、プリプロセス回路１０ではこの観察対象となる陰影を背景に対して大きな画素値となるような操作を加えることも行われる。
【００４２】
プリプロセス処理の後、画像信号は走行方向検出回路１４に供給され、走行方向とその強さが検出される。この処理は画像上の関心領域毎に行われ、図２に示すような線状陰影（ａ）またはエッジ（ｂ）を有するパターンにおいて、その線状陰影もしくはエッジがどの方向を向いて走行しているかを検出する。この検出の結果、方向とその強さが出力される。走行方向は関心領域内のパターンの走行方向であるから、複雑に走行している陰影であっても関心領域の大きさを適切に定めると関心領域内ではほぼ直線の走行とみなせ、容易に検出できる。また、走行方向の強さとは、はっきりした陰影では強く、ぼやけたような陰影では弱くなるような指標である。例えば、図３（ａ），（ｂ）は走行方向の強さが図２（ａ），（ｂ）に比べて小さい例である。
【００４３】
この２つの出力はベクトルフィルタ回路１６により修正処理される。この修正処理の目的は、ノイズに埋もれた信号がノイズによって誤った方向や強さで求められても、周辺の関心領域の出力結果とある程度相関を持たせることにより、誤検出の影響を小さくするものである。例えば、ノイズパターンの中で偶然に線状陰影のようなパターンが形成され、方向を強く検出したとしても周辺の方向がバラバラな場合は方向の強さを弱める作用がある。また、逆にノイズが偶然に線状陰影の一部の領域でコントラストを失わさせ、あたかも線が切れたようになった場合は、その領域で間違った方向を検出するが、周辺の出力が一定の方向で検出されればそのような領域でも方向を修正する。
【００４４】
ベクトルフィルタ回路１６の出力が空間フィルタ係数発生回路１８に供給され、空間フィルタ回路１２のフィルタ係数が決定される。これにより、ベクトルフィルタ回路１６の出力に応じて、関心領域内部で陰影の走行方向と強さに応じた空間フィルタが施される。
【００４５】
図４に空間フィルタ回路１２の係数の一例を示す。（ａ）は水平方向に強い方向性を持つような陰影を検出したとき、（ｂ）は水平方向に弱い方向性を持つような陰影を検出したとき、（ｃ）は水平方向から約３０度斜めに傾いているような少し強めの方向性を検出したときに、それぞれ施すべき空間フィルタである。例示した空間フィルタは係数の和が１０になるようにしているが、実際の画像に対してはフィルタ後の画像のレベルが変わらないように、各係数を１／１０倍にして和が１になるように調整される。これらの空間フィルタは方向に応じて周波数低減効果が異なる。例えば、図４（ａ）は水平方向には周波数低減効果があるが、垂直方向には全く何のフィルタ操作をしないため、水平方向に走行する線状陰影の周辺をぼかすことがない。
【００４６】
このように画像の各関心領域毎に陰影の方向とその強さに応じて、線状またはエッジの走行方向にノイズ低減を行うことが本発明の基本的な原理である。
【００４７】
以下、図１の各部の具体的な作用を説明する。
【００４８】
プリプロセス回路１０はエッジを検出しやすくしたり、ノイズを低減させて、後段の処理においてノイズによる誤認識を減らすことが目的である。例えば、一般に使用されている通常の３×３マトリクスの係数をすべて１／９に設定した線形フィルタは周波数低減作用があることが広く知られているが、多少そのような操作をしてノイズを低減すれば、非常にノイズの多い画像に対しては後段の処理において誤認識が減る。また、辺縁がぼけた画像が入力された場合は、エッジ強調作用をもつようなバンドパスフィルタ処理等を施せば、線状陰影をよりはっきり描出せしめるため、後段の走行方向検出処理の精度が向上する。また、Ｘ線画像ではカテーテルやガイドワイヤは金属やＸ線阻止能の大きい物体で構成されることが多いため、背景に対して小さい画素値で陰影が形成される。プリプロセス回路１０では、この陰影が背景よりも大きな値を取るような操作を行う。これは後段の走行方向検出回路１４は、このような陰影の空間的広がりを計算し、広がりの大きい方向を走行方法とする処理を行うため、画素値を大きくしておく必要があるためである。なお、この処理をしない場合は走行方向検出回路１４において空間的な広がりの小さい方向を走行方向とすればよい。しかし、ここでは、この観察対象となる陰影が背景よりも大きな画素値をとる状況に対して説明する。観察対象となる陰影の濃度が背景よりも小さな場合のプリプロセスにおいて具体的に最も簡単な方法は画素値の符号を反転するだけの処理である。なお、符号付き演算を行わないような処理装置では、入力画像全面の画素値の最大値から入力画像を減じる処理でもよい。例えば、入力画像が１２ビットで表現されている場合には、入力画像の最大値は２１２−１＝４０９５である。したがって、プリプロセスでは４０９５−（入力画像）を演算すると入力画像に含まれている線状陰影の画素値は背景よりも大きな値となる。
【００４９】
本発明の中で最も重要な部分は走行方向やその強さを検出する走行方向検出回路１４である。この検出手段の一つとして、主成分分析法に準じた方式をまず以下に説明する。
【００５０】
図５（ａ）のような画像の関心領域をまず仮定する。これは、例えばＸ線透視下においてガイドワイヤが人体に挿入されたとき、そのガイドワイヤが描出されている部分画像である。この関心領域内の画素濃度分布は中心を（０，０）として固定された座標点（ｉ，ｊ）に画素値Ｄ（ｉ，ｊ）の個数だけデータが散布されたものと考えると、線状パターンの走行方向とはデータの空間的分布がもっとも広がった方向と考えることができる。図５（ａ）の例では画素濃度０の上に濃度１０の線状パターンが走行している例であるが、この画像は図５（ｂ）の黒丸の座標位置、例えばアドレス（ｉ，ｊ）＝（２，−１）の位置にＤ（２，−１）＝１０個のデータが存在しているものと解釈する。画像をこのようにアドレス（ｉ，ｊ）にＤ（ｉ，ｊ）個のデータが分布したものとしてとらえると、主成分分析手法が応用可能であり、共分散行列を次のように表すと、Ｓ１１、Ｓ１２、Ｓ２２は次のようになる。
【００５１】
【数１】
0004112762.tif 000002
【００５２】
ここで、Ｍは関心領域の範囲を表す。図５に例示した関心領域ではＭ＝３である。また、Ｘｇ，Ｙｇ，Ｚはそれぞれ水平方向、垂直方向の重心とデータ個数の総和を示しており、それぞれ下式で表される。Ｓ’１１、Ｓ’１２、Ｓ’２２は原点の周りの２次モーメントである。
【００５３】
【数２】
0004112762.tif 000003
【００５４】
この演算において、Ｚは関心領域内の画素値の総和であるが、Ｚ／（２Ｍ＋１）２は関心領域の画素値の平均値を表す。したがって、仮に入力画像の平均値が別の手法で求められている場合は、Ｚを求める必要がない。例えば、後述する図１３に示されているような変形例は、プリプロセス回路１０で画像平均値を求めるため、そこで求められた結果からＺを推定することが可能である。
【００５５】
そこで、共分散行列Ａの固有値は次のように求められる。
【００５６】
【数３】
0004112762.tif 000004
【００５７】
固有ベクトルの向きをそれぞれの固有値に対してθ１、θ２とすると
ｔａｎ（θ１）＝（λ１−Ｓ１１）／Ｓ１２            （ｅ−１２）
ｔａｎ（θ２）＝（λ２−Ｓ１１）／Ｓ１２            （ｅ−１３）
となる。つまり、固有値λ１、λ２のうち値の大きいほうに対応する固有ベクトルの向きが分散を最大にする主軸の向き（空間的広がりの大きい方向）であるから、λ１のほうを採用してθ１を求めることができる。なお、θ１は主軸の向きであるから０から１８０度の範囲で設定することができる。
【００５８】
また、方向性の強さの定義の方法でもっとも一般的なものは
Ｃ＝λ１／（λ１＋λ２）                          （ｅ−１４）
として寄与率を求める方法であり、Ｃ＝０〜１の間で求めることができる。つまり画像の関心領域の画素値Ｄ（ｉ，ｊ）、（ｉ＝−Ｍ，−（Ｍ−１），．．．，Ｍ；ｊ＝−Ｍ，−（Ｍ−１），．．．，Ｍ））に対してＳ１１，Ｓ２２，Ｓ１２を演算し、θ１とＣを求めれば走行パターンの方向と強さの指標をそれぞれ求めることができる。
【００５９】
以上、走行方向と強さを求めるための原理を説明したが、以下に実際の演算を行う場合の実施例を示す。例えば、水平方向の重心Ｘｇを求めるには、図６に示すように、先ず、画像の垂直方向の総和を各列毎に計算し、２Ｍ＋１個の数値を求めた後、水平方向に積和演算すればよい。この演算は一般の線形一次元フィルタと同じ演算方式でもよく、フィルタ係数は原点からの距離に応じて定められる。つまり、水平方向の重心Ｘｇは関心領域内の画素に対して以下の操作で求められる。
【００６０】
（１）フィルタ係数をすべて１として垂直方向に１次元フィルタを施すのと同じ演算で出力を取り出す。
【００６１】
（２）この操作を全ての列に対して行い、水平方向の数値の並び（２Ｍ＋１個）を求める
（３）２Ｍ＋１個の数値に対してフィルタ係数＝（−Ｍ，−（Ｍ−１），−（Ｍ−２），．．．，（Ｍ−１），Ｍ）として水平方向に１次元フィルタをかけ、この出力をＸｇとする。
【００６２】
同様にして垂直、水平方向のフィルタ係数を入れ換えると、垂直方向の重心Ｙｇを求めることができる。
【００６３】
データ個数の総和Ｚはフィルタ係数をすべて１．０とすれば、同じハードウェアで求めることが出来る。
【００６４】
さらにＳ’１１，Ｓ’１２，Ｓ’２２などもフィルタ係数を変えて積和演算を行えばよく、Ｓ’１１は次のように求めることが出来る。
【００６５】
（１）フィルタ係数をすべて１として垂直方向に1次元フフィルタを施すのと同じ演算で出力を取り出す。
【００６６】
（２）この操作を全ての列に対して行い、水平方向の数値の並び（２Ｍ＋１個）を求める。
【００６７】
（３）２Ｍ＋１個の数値に対して、フィルタ係数＝（Ｍ２，（Ｍ−１）２，（Ｍ−２）２，．．．，（Ｍ−１）２，Ｍ２）として水平方向に1次元フィルタをかけ、この出力をＳ’１１とする。
【００６８】
また、Ｓ’２２は水平と垂直のフィルタ係数を入れ替えれば求めることができる。Ｓ’１２に対しては、次のように求める。
【００６９】
（１）フィルタ係数を（−Ｍ，−（Ｍ−1），−（Ｍ−2），．．．，（Ｍ−1），Ｍ）として垂直方向に1次元フィルタを施すのと同じ演算で出力を取り出す。
【００７０】
（２）この操作を全ての列に対して行い、水平方向の数値の並び（２Ｍ＋１個）を求める
（３）２Ｍ＋１個の数値に対してフィルタ係数＝（−Ｍ，−（Ｍ−1），−（Ｍ−２），．．．，（Ｍ−1），Ｍ）として水平方向に1次元フフィルタをかけ、この出力をＳ’１２とする。
【００７１】
上記の演算をまとめると、図7に示すような構成で実施される。垂直方向の積和演算によって３つのベクトルＶ０（ｉ），Ｖ１（ｉ），Ｖ２（ｉ）を求め、これらを水平方向の積和演算することにより、（ｅ−４）式〜（ｅ−９）式の各値を求めることが出来る。
【００７２】
これらの数値が求まった後は、単なる数値の演算であるからＲＯＭや高速ＣＰＵなどを用いて高速演算することも可能である。
【００７３】
関心領域のサイズはいかなるサイズでも計算が可能であるが、７×７以上が望ましい。線状のパターンが含まれない領域や、ノイズだけの領域では、走行方向の強さが非常に小さくなることが期待されるが、これはノイズに対して水平、垂直方向の分散が等しいこと、つまりＳ１１＝Ｓ２２と近似できることが必要条件である。このためには充分な統計量が必要であり、一般的には５０個以上の統計量が望ましく、それゆえに７×７以上のサイズが望ましい。
【００７４】
さて、画像全体に対してはこの関心領域の中心位置を少しずつ変えて局所パターンの走行方向と強さを求めることが行われる。例えば５００×５００マトリクスの画素から成る画像に対し、関心領域の中心位置を２画素ずつ水平方向にシフトしながら走行方向と強さを求めていくと水平方向を２５０個に分割した領域毎の走行方向とその強さの情報が得られる。同様に、垂直方向にも２画素ずつシフトしながら画像全体にわたって走査をすれば２５０×２５０の情報が得られる。この関心領域のシフト量が小さい場合は、あるひとつの関心領域とシフトした関心領域の間で重なった画素情報が多いため、求まった走行方向と強さは互いに高い相関がある。つまり、冗長度が大きい。例えば、１画素ずつシフトした場合は、５００×５００の局所情報を得るが、画像全体の演算回数が大きい割には互いに冗長性の大きい結果であるため、多少もったいない演算をしていることになる。逆に、あまりシフト量が大きいと、演算回数は少ないが、局所情報の独立性が大きくなり、画像上のノイズに影響されやすくなる。隣接した関心領域とある程度の相関を保つことは次のベクトルフィルタを施す場合にも必要であり、シフト量は関心領域のサイズの半分以下が望ましい。例えば、好適な実施例としては、７×７のサイズを持つ関心領域を２画素シフトしながら画像全体を走査して局所的な走行方向と強さを求める手法が推奨される。
【００７５】
走行方向検出回路１４により走行方向θ１と強さＣを求めた後、ベクトルフィルタ回路１６により近傍の情報に基づいて、より確かな方向と強さを求める処理を行う。図８に示すのはそのベクトルフィルタ１６の一例である。画像の関心領域の走行方向が場所毎に図のように求まったとき、３×３の画素の中央の画素の走行方向がノイズによって間違った方向として求まったとしても、ベクトルフィルタ１６を施すことにより、近傍の画素の走行方向の情報から修正される様子を示している。
【００７６】
このベクトルフィルタ１６は以下のように構成される。
【００７７】
１つ目の好適な構成例は、近傍の小領域の方向ベクトルの角度を分類して多数決決定する方法である。図９には２２．５毎に分類した概念図を示している。例えば、この例の場合は、４５度を中心として±１１．２５度の範囲に入る方向が最も多いため、ベクトルフィルタの出力は４５度という結果になる。一方、走行方向の強さは全方向ベクトルの強さの総和に対して多数決で選択されたベクトルの強さだけの総和で決めることができる。例えば、対象となる小領域とその近傍の８つの領域の走行方向θと強さＣが次のような場合を考える。
【００７８】
小領域位置              方向θ（度）        強さＣ
＃０  対象部              １３１              ０．３
＃１  対象部の右            ８０              ０．５
＃２  対象部の右上          ４４              ０．１
＃３  対象部の上            ４２              ０．９
＃４  対象部の左上          ７０              ０．７
＃５  対象部の左            ４５              ０．８
＃６  対象部の左下          ５０              ０．８
＃７  対象部の下            ４０              ０．５
＃８  対象部の右下          ４９              ０．６
多数決で４５±１１．２５度の範囲に入る方向を持つ小領域は＃２、＃３、＃５、＃６、＃７、＃８の６つであるから、次のように決定する。
【００７９】
Ｃ＝（０．１＋０．９＋０．８＋０．８＋０．５＋０．６）／
（０．３＋０．５＋０．１＋０．９＋０．７＋０．８＋０．８＋０．５＋０．６）
＝３．７／５．２
＝０．７
ベクトルフィルタの２つ目の好適な構成例は、近傍の小領域のＳ１１、Ｓ２２、Ｓ１２を全て平均化した行列Ａｇを計算し、Ａｇ行列に基づいて方向ベクトルとその強さを求める方法である。
【００８０】
行列Ａｇは次のように与えられる。
【００８１】
【数４】
0004112762.tif 000005
【００８２】
近傍小領域のＳ１１等の和を算出するとき、単純な平均ではなく、重み付き平均を利用してもよい。
【００８３】
以上の２つの好適な構成例により説明したベクトルフィルタの目的は対象とした小領域の周辺まで含めてある程度広い範囲内でパターン走行方向とその強さを推定するためであり、ノイズなどで特異的な方向にパターン走行方向が誤認されても、周辺の領域から修正を加えることができる。
【００８４】
ベクトルフィルタ回路１６により近傍の関心領域での走行方向検出結果に基づいてベクトルの修正を行った結果、この走行方向に特有の周波数低減フィルタを施す。このフィルタの例は、既に図４で示されているが、このフィルタ係数は局所的に求められた走行方向と強さの情報に基づいて空間フィルタ係数発生回路１８によって画素毎に変更される。図１０には入力画像、ベクトルフィルタ回路１６の出力例、およびその情報に基づいて施すべき５×５の空間フィルタを３つの局所に対して例示している。入力画像の中の線状陰影やエッジ部は局所的に強い方向成分を持っているが、そういう部分にはその方向に空間周波数を低減する空間フィルタをかけ、背景部には全方向性の空間フィルタが施されることが説明されている。
【００８５】
このように本発明では、画素単位に係数が設定可能な空間フィルタを施すわけであるが、この空間フィルタが大型であればあるほど実時間演算がしにくくなる。そのため、１つの大型なフィルタを施す代わりに、小さいサイズの空間フィルタを連続して施すことで演算回路を簡易化させることもできる。例えば５×５の空間フィルタは２５個の積和演算を行う必要があるが、３×３の空間フィルタを２回連続して施すようにすれば１８個の積和演算で済み、回路規模を小さくすることができる。図１１にはその一例を示す。フィルタ出力１とフィルタ出力２は同一の結果となる。
【００８６】
以上説明したように、本実施形態によれば、画像に含まれる線状もしくはエッジ部のパターンの走行方向、およびその強さに応じた空間フィルタをその走行方向に沿って施すことで、線状陰影やエッジ部の辺縁をぼかしたり、コントラストを低下することなく、ノイズ低減を達成することができる。
【００８７】
以下、本発明による画像処理装置の他の実施形態を説明する。他の実施形態の説明において第１の実施形態と同一部分は同一参照数字を付してその詳細な説明は省略する。
【００８８】
第２実施形態
第２の実施形態のブロック図は第１の実施形態のブロック図と同一であるので、図示省略する。第２実施形態は、走行方向の検出処理の変形に関する実施形態である。
【００８９】
第１実施形態では、Ｓ１１，Ｓ１２，Ｓ２２は（ｅ−１）式から（ｅ−３）式に示した式により求めた。しかし、一般に重心は画像の画素値の平均値によって変動する。したがって、画像の平均レベルが大きいときには重心は関心領域の原点に近いところで求まり、小さいときには重心が原点から離れる。また、ある２つの画像の差分をとったような画像、例えばＸ線診断装置で得られるＤＳＡ（digital subtraction angiography）画像のように画素平均値がゼロに近い場合は、データ個数の総和Ｚがゼロになることがあり、（ｅ−１）式等の演算がエラーを起こすことがある。したがって、本発明の耐用性を高める手段として画素値の平均値に依存しない手法を第２実施形態として示す。
【００９０】
第２実施形態では、（ｅ−１）式、（ｅ−２）式、（ｅ−３）式に代えて下記の式によりＳ１１，Ｓ１２，Ｓ２２を求める。
【００９１】
【数５】
0004112762.tif 000006
【００９２】
この演算の意義は関心領域内の画素値の平均値からの差の成分に対する重心をＸｇ’，Ｙｇ’として定義したものである。
【００９３】
このような定義で重心を求めることにより、平均値Ｄavがいかなる値をとっても重心位置は平均値からの差分の成分だけを対象にして求めるわけであるから、平均値には依存しない値として求められる。
【００９４】
この演算を実際のハードウェアで行うには、プリプロセス回路１０で次の演算を行えばよい。
【００９５】
プロプロセス回路出力
＝入力画像−入力画像のローパスフィルタ像（＝Ｄav）    （ｅ−２６）
この演算によれば、走行方向検出回路１４に入力される画像は局所的に平均値がゼロとなった画像となるため、Ｄav＝０として（ｅ−２４）式を演算することができ、ハードウェア構成が簡素化される。
【００９６】
第３実施形態
第３実施形態も、走行方向の検出処理の変形に関する実施形態である。本実施形態では、Ｓ１１，Ｓ１２，Ｓ２２を単純に原点の周りの単純な２次モーメントで表すことであり、次のように求められる。
【００９７】
Ｓ１１＝Ｓ’１１                                        （ｅ−２７）
Ｓ１２＝Ｓ’１２                                        （ｅ−２８）
Ｓ２２＝Ｓ’２２                                        （ｅ−２９）
この式を用いて行列Ａの固有値を求めると、線状陰影から関心領域の中心が離れていくにしたがって、方向性の強さが次第に弱まって求められ、さらに離れていくと、線状陰影の走行方向とは別の走行方向を検出することがあり得るが、線状陰影の付近の関心領域の方向性は正しく求められるので、実用に供せられる。ただし、このときには同時にＸｇ，Ｙｇの演算も行い、Ｘｇ，Ｙｇの計算結果からあまりに大きな重心のシフトがあることが分かる場合は、関心領域の中心位置は線状陰影から離れているわけであるから、強制的に方向性の強さをゼロとするなどの操作を行う必要がある。
【００９８】
第４実施形態
第４実施形態も、走行方向の検出処理の変形に関する実施形態である。走行方向θ１は数学的に画像の平均値に依存せず求まるが、強さを示す寄与率Ｃは平均値が大きいときには次第に小さくなるという数学的性質を持つ。また、画像の平均値をゼロにした場合、λ１＋λ２＝０となり方向性の強さを正しく求めることができない場合もある。
【００９９】
これを解決するために画像の平均値に依存しない量Ｐを次のように定義し、
Ｐ＝λ１−λ２                                        （ｅ−３０）
いろいろな条件下で予め測定しておいたλ１−λ２と比較して、走行方向の強さを表す指標Ｃ’を得ることができる。
【０１００】
Ｐは画像上に線状もしくはエッジ状のパターンがなくノイズだけが分布している場合は、関心領域のサイズ、および画像上のノイズの標準偏差σによって定まる確率分布に従う。この定式化は複雑であるが、ほぼ正規分布に従う。一方、画像上のノイズの標準偏差σは画像の撮影条件などから知ることができる。すでに述べたように、Ｘ線画像の場合は、画像のノイズの標準偏差σはＸ線フォトン数Ｎの平方根に比例するため、画像の平均値をＬとするとＸ線フォトン数を推定することができ、σも求められる。例えば、１つのＸ線フォトンに対してＬ０という画素値が生じるような場合、画像の平均値Ｌを観測すれば、
Ｎ＝Ｌ／Ｌ０
という関係があるため、σは
σ＝Ｌ０（Ｌ／Ｌ０）1/2
となる。つまり、画像の平均値を得ることにより、ノイズだけが存在する場合のＰの確率分布を計算することができる。この分布をｆ（Ｐ）（ｆ（Ｐ）は積分すると１．０となるような分布）とすれば、実際に関心領域内の走行方向検出処理で求まったＰｘ＝λ１−λ２がこの分布のどこに位置するかを求めて、０〜１までの間の確からしさを算出することができる。図１２はあらかじめ求められているｆ（Ｐ）、および実際に得られたＰｘの例であるが図中の斜線部の積分値ηを用いて
Ｃ’＝１−η
とすれば、走行方向の強さを表す指標Ｃ’を得ることができる。
【０１０１】
第４実施形態の信号の流れを図１３に示す。プリプロセス回路１０では画像の平均値を求める処理を行い、その平均値に対応するノイズの標準偏差σを基に、ｆ（Ｐ）を発生する。一方、走行方向検出回路１４はＰｘ＝λ１−λ２を求める。ｆ（Ｐ）分布のＰｘの位置からηを算出して、Ｃ’を求める。
【０１０２】
第５実施形態
次に、第５実施形態として、本発明の画像処理装置の具体的な応用例としてのＸ線診断装置を説明する。従来からＸ線診断装置は検体を非破壊で検査できる有用な装置であり、高精細で動画が得られる事が望まれている。さらに、装置の大型、重厚化を避けるために、イメージインテンシファイヤの代わりに平面型検出器が必要とされている。これらの要請を達成する為に、アクティブマトリックス型の平面検出器が提案されている。この様な平面検出器には間接変換方式と直接変換方式とがあるが、Ｘ線照射により発生する電子・正孔対を直接検出する直接変換方式は、間接変換方式と比べて高精細な画像を得る事ができる。この直接変換方式は、例えば、米国特許第５，３１９，２０６号明細書（E.I.DuPont）や、ＳＰＩＥ論文誌，vol. 243，P237，1995年（D.L.Lee et al.）や、特開平１１−４８２１号に記載されている。図１４に特開平１１−４８２１号公報に記載の検出器を示す。
【０１０３】
図１４（ａ）、（ｂ）に沿って説明する。まず、光導電体３０１に入射したＸ線Ｘは、電子・正孔対３０９を発生させ、この発生した電荷は信号容量３０５と、画素電極３０４と対向する電極３００とで形成される容量とで分配される。Ｘ線入射の電荷は、トランジスタのゲートＧがゲート信号Ｇｌからの信号によって開く事によりドレインＤ１〜Ｄ４を通して図示しない積分器に集められる。直接変換方式では、光導電体３０１に入射した情報が全て電子・正孔対３０９として収集されるという点で、画素電極形状を精細化することで、高精細な画像を実現する事が可能である。ここで、３０８は電源である。
【０１０４】
ところで、Ｘ線を利用した医療技術として透視下におけるカテーテル治療が盛んに行われているが、Ｘ線透視は撮影に比べて線量を低減して患者被爆を低減していることから、画像上に重畳するノイズが大きく、カテーテル、その他のガイドワイヤや血管等の線状陰影が背景ノイズに邪魔されて、見えにくくなるという問題がある。これは、上述したような直接変換方式の平面検出器を備えたＸ線診断システムにて行われることを考慮した場合、その原理上高精細な画像を形成し得ることから、低線量時にはＸ線の量子ノイズがより顕著に現れることに起因する。
【０１０５】
しかしながら、直接変換方式の平面検出器を備えたＸ線診断システムに本発明の画像処理装置を組み合わせることにより、本来、高分解能、かつ高精細の面で優れた直接変換方式の平面検出器の特徴を生かしつつ、線状陰影やエッジ状陰影のコントラストを損なわずに、背景のノイズを選択的に抑圧でき、直接変換方式の平面検出器にとって、特によい結果をもたらすことができる。すなわち、被検体を透過したＸ線を図１４に示すような直接変換方式の平面検出器により検出して得た被検体のＸ線像を図１の入力画像として本発明の画像処理装置に入力し、線状陰影やエッジ状陰影の方向とその強さに応じた指向性を有する空間フィルタ処理を施すことにより、高精細の画像が得られるＸ線診断装置が実現される。そのため、低線量のカテーテル治療時においても、線状陰影のコントラストを損なわずに、背景ノイズを選択的に抑圧できる。第５実施形態のＸ線診断装置は、上述した第１〜第４実施形態のいずれも適用可能である。
【０１０６】
本発明は上述した実施形態に限定されず、種々変形して実施可能である。例えば、上述の説明では、Ｘ線画像を例にとったが、これに限らず、ＣＴ画像やＭＲＩ画像等、ランダムノイズが重畳するような画像全般に対しても同様に適用可能である。
【０１０７】
【発明の効果】
以上説明したように本発明によれば、線状陰影もしくはエッジ状の陰影を含む画像において、線状陰影やエッジ状陰影のコントラストを損なわずに、背景のノイズ成分を選択的に、しかもリアルタイムで抑圧することができる画像処理装置を提供することができる。
【図面の簡単な説明】
【図１】本発明による画像処理装置の第１の実施形態の構成を示すブロック図。
【図２】処理対象の画像の一例を示す図。
【図３】図２よりも方向の強さが弱い画像の例を示す図。
【図４】第１実施形態の空間フィルタの例を示す図である。
【図５】第１の実施形態の走行方向検出回路が主成分分析法を適用して走行方向、強さをするための原理を示す図。
【図６】第１の実施形態の走行方向検出回路が走行方向、強さをするための積和演算を実行する部分の回路図。
【図７】第１の実施形態の走行方向検出回路の具体的な回路図。
【図８】第１実施形態のベクトルフィルタ回路の動作の概略を示す図。
【図９】第１実施形態のベクトルフィルタ回路の動作の詳細を示す図。
【図１０】第１実施形態の空間フィルタ回路の動作を示す図。
【図１１】第１実施形態において、１つの大型の空間フィルタを施す代わりに、複数の小型の空間フィルタを施す動作を示す図。
【図１２】本発明による画像処理装置の第２実施形態の走行方向検出の原理を示す図。
【図１３】本発明による画像処理装置の第４実施形態の構成を示すブロック図。
【図１４】本発明による画像処理装置の第５実施形態に用いられる直接変換方式の平面型検出器の構成を示す図。
【符号の説明】
１０…プリプロセス回路
１２…空間フィルタ回路
１４…走行方向検出回路
１６…ベクトルフィルタ回路
１８…空間フィルタ係数発生回路
【図１】
0004112762.tif 000007
【図２】
0004112762.tif 000008
【図３】
0004112762.tif 000009
【図４】
0004112762.tif 000010
【図５】
0004112762.tif 000011
【図６】
0004112762.tif 000012
【図７】
0004112762.tif 000013
【図８】
0004112762.tif 000014
【図９】
0004112762.tif 000015
【図１０】
0004112762.tif 000016
【図１１】
0004112762.tif 000017
【図１２】
0004112762.tif 000018
【図１３】
0004112762.tif 000019
【図１４】
0004112762.tif 000020