(57)【特許請求の範囲】
【請求項１】
Ｘ線検出器により取得された蛍光透視画像のシーケンス内のノイズを減らす方法であって、画像が１個以上の特徴を含み、画像がカーブレット係数ｃｔ，ｓ，θｓｂ，ｉ，ｊにより表され、ここで、ｓは前記画像のスケール、θｓｂは前記画像内の特徴の角度に紐付けられた前記特徴の副帯域、および（ｉ，ｊ）は前記特徴の相対位置である、方法において、少なくとも以下のステップ、すなわち、
・副帯域毎にカーブレット係数を用いて、且つ副帯域θｓｂに関連付けられたフィルタの最大応答を判定することにより、カーブレット領域内の画像内の前記特徴を検出（３６）し、前記検出された特徴をメモリに保存するステップと、
・前記検出された特徴を局所的に照合するステップであって、
・カーブレット係数の集合を対象として、対象とするカーブレット係数が時点ｔ以前に検出された特徴に属するか否かを検証し、
・前記カーブレット係数が検出された特徴に属する場合、前記係数が属する前記特徴が、より早い時点で検出された特徴と一致するか否かを検証（４４）し、一致する場合、前記カーブレット係数を第１の集合に分類し、
・前記係数が検出されない（４２）場合、または前記カーブレット係数が検出された特徴に属するがより早い時点で検出された特徴と一致しない（４５）場合、少なくとも第２の集合を形成する、ステップと、
・前記第１の集合に属する前記カーブレット係数の時間領域フィルタリング用に第１の係数を、および第２の集合に属する前記カーブレット係数の時間領域フィルタリング用に少なくとも第２の係数を適用するステップと、
・前記時間領域フィルタリングされたカーブレット係数を用いてノイズ除去済み画像を決定するステップと
を組み合わせて含むことを特徴とする方法。
【請求項２】
以下のステップ、すなわち、
・副帯域毎に前記カーブレット係数を連結する前記カーブレット領域内の処理により前記画像内の前記１個または複数の特徴を検出するステップであって、
・分解において最高周波数カーブレットｓｍａｘにおけるスケールを採用し、
・各カーブレット係数ｃｔ，ｓ，θｓｂ，ｉ，ｊについて、高周波情報を含む低スケールでその親係数を計算し、
・次いで前記カーブレット係数を合算（３６）することにより、連結された係数ｃｃｏｎｃａｔｅｎａｔｅｄのマトリクスを取得し、
・前記特徴の前記検出を前記連結された係数のマトリクスに適用し、
・各々の連結された係数について、前記副帯域に関連付けられた有向フィルタの最大応答を計算（３７）し、
【数１】
0006375373.tif 000027
ここで、ｐ、ｑは向きを決定する位置であり、
・最大応答の結果を、所与の画像について維持および記録される前記特徴を与えるヒステリシスにより閾値化（３８）して、係数が高閾値より大きい場合、低閾値よりも高いその拡張隣接係数の全てを維持し、
・最高周波数のスケールにおける前記特徴の位置を選択することにより基準特徴を画定する、ステップ
を含む、前記特徴を検出するステップを含むことを特徴とする、請求項１に記載の方法。
【請求項３】
以下のステップ、すなわち、
・前記時点ｔより早い時点ｔ−１で決定されたカーブレット係数の集合を対象とするステップであって、
・カーブレット係数が前記時点ｔまたは前記時点ｔ−１で検出された特徴に属するか否かを検証（４１）し、
・前記検証の結果が否定的（４２）であれば、前記カーブレット係数は被写体に属していないものと分類し、
・前記検証の結果が肯定的（４３）であれば、
・前記係数が属する前記特徴が前記時点ｔ−１または時点ｔ−２で検出された特徴に一致するか否かを検証（４４）し、
・前記検証の結果が否定的（４５）であれば、前記カーブレット係数は被写体に属していないものと分類し、
・前記検証の結果が肯定的（４６）であれば、前記カーブレット係数は被写体に属しているものと分類する、ステップと、
・少なくとも２種類のカーブレット係数に対応する少なくとも２組の係数を取得するステップと
を含む、前記特徴を局所的に照合するステップ（４０）を含むことを特徴とする、請求項１に記載の方法。
【請求項４】
以下のステップ、すなわち、
・被写体に属していない前記カーブレット係数について、前記時点ｔで検出された前記係数と時点ｔ−１で検出された前記係数とを平均化するステップと、
・被写体に属する前記カーブレット係数について、より弱いフィルタリングを適用するステップと
を含む、時間領域フィルタリングのステップ（４７）を含むことを特徴とする、請求項１に記載の方法。
【請求項５】
前記ノイズ除去済み画像が、前記時間領域フィルタリング済み係数を用いて決定（４９）される逆変換ステップ（４８）を含むことを特徴とする、請求項１に記載の方法。
【請求項６】
前記時間領域フィルタリングステップが、２個の利得係数αおよびβ（ここで、０＜α≦β）を定義し、
・前記カーブレット係数が検出されなかった場合、または検出されたが一致しなかった場合、前記時点ｔで前記利得αを前記係数に適用し、および時点ｔ−１で利得（１−α）を同一係数に適用し、
【数２】
0006375373.tif 000028
ここで、
【数３】
0006375373.tif 000029
は、前記時点ｔでのノイズ除去済み係数であり、
【数４】
0006375373.tif 000030
は、前記時点ｔでのノイズ係数であり、
【数５】
0006375373.tif 000031
は、前記時点ｔ−１でのノイズ除去済み係数であり、
・前記係数が前記時点ｔまたは前記時点ｔ−１で検出された場合、利得βを前記係数に適用すると共に、その隣接係数から前記検出された係数までの距離の関数として所与の減少率でβからαまで減少する利得を前記係数の前記隣接係数に適用することを特徴とする、請求項１に記載の方法。
【請求項７】
前記利得の前記減少率が線形または正規分布に従い適用されることを特徴とする、請求項６に記載の方法。
【請求項８】
前記時点ｔで検出された前記カーブレット係数が、時点ｔ−１で検出され且つそれらと一致する前記カーブレット係数の集合と平均化されることを特徴とする、請求項１に記載の方法。
【請求項９】
前の時点で特徴に属していた前記係数、または現時点で一致した特徴に属する前記係数、あるいは前記２種類の係数のいずれかの平均化において、前記時点ｔでノイズ除去された前記カーブレット係数をノイズカーブレット係数の代わりに用いることを特徴とする、請求項８に記載の方法。
【請求項１０】
同一被写体に属する前記特徴に照合すべく、以下のステップ、すなわち、
・コスト関数を定義し、ならびに前記特徴の長さ、強度、変位、および角度を類似度指標として使用するステップと、
・全エネルギーを次式で定義するステップであって、
【数６】
0006375373.tif 000032
ｗｌ、ｗｉ、ｗｄおよびｗａは前記指標の信頼度に対応し、ここで、ｗｌ、ｗｉ、ｗｄおよびｗａ＝１である、ステップと、  ・前記全エネルギーを、一方の特徴がｔで、他方がｔ−１で検出された特徴の全てのペアについて計算するステップと、  ・エネルギーの重みを最小化する照合リストを取得すべく適応された方法を適用するステップとが実行され、  ・前記カーブレット係数が特徴に対応していない場合、前記方法が時間平均化ステップを実行することを特徴とする、請求項１に記載の方法。
【請求項１１】
比較される最小の特徴の長さＬｍｉｎおよび比較される最大の特徴の長さＬｍａｘが定義され、
【数７】
0006375373.tif 000033
ここで、Ｉｍｉｎは前記２個の特徴のうち低い方の強度であり、Ｉｍａｘは他方の特徴の強度であり、
【数８】
0006375373.tif 000034
指数ＩおよびＥが前記時点ｔおよびｔ−１における前記特徴の極値を表し、μが受容可能な動きを固定させるパラメータであり、ｍが寸法であり、ｎがｍでないケースではμを除算する寸法であり、
【数９】
0006375373.tif 000035
ここで、θは前記時点ｔおよびｔ−１で検出された前記特徴の向きを表し、νは向きの変動が大き過ぎる場合にＥａｎｇｌｅを１に近くするためのパラメータであり、［π］はθの単位がラジアンであることを示すことを特徴とする、請求項１０に記載の方法。
【請求項１２】
｜ξ×β｜より大きい利得を有する前記カーブレット係数に係数ρを乗算することにより、画像のコントラストを高めるステップを含み、１＜ρであり、ここで、ξはα＜ξβ＜βであるように選択されることを特徴とする、請求項６に記載の方法。
【請求項１３】
乗数係数ρを、スケール内近傍でρから１まで減らしながら、前記検出されたカーブレット係数に適用することを特徴とする、請求項１〜１２のいずれか一項に記載の方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明の主題は、カーブレット領域における特徴検出、および検出された特徴と特徴に関連付けられた係数の種類に応じた時間領域フィルタリングとの照合により、蛍光透視画像シーケンスのノイズを減らす方法に関する。
【０００２】
本発明は特に、Ｘ線撮影の分野、より具体的には、医療画像分野で用いられる。
【背景技術】
【０００３】
医療画像の分野において、Ｘ線蛍光透視法による撮影は、「最小限に侵襲的」と呼ばれる外科的処置を行う間、患者の画像ストリームを医師に提供できるようにする。これらの「最小限に侵襲的」な手術は、患者の手術後の痛みを和らげると共に感染のリスクを減らすため、利用が普及しつつある。これらの手術は、例えばＸ線により得られた画像シーケンスを用いて、患者の動脈系内で医療器具を追跡する能力を要する。このような処置は長時間かかり、患者が被曝する放射線量は患者の健康に有害な影響を及ぼす恐れがある。従って患者だけでなく医療従事者の健康のためにも、受けるＸ線量を可能な限り減少させる必要がある。しかし、Ｘ線量を減らすことはノイズの増大および画像コントラストの低下につながる。従って、１画像当たりのＸ線量を減らせることが望まれる場合、高画質の画像を表示するために利用可能なノイズ除去またはノイズ減少アルゴリズムを有していることが重要である。
【０００４】
画像内のノイズを減らすべく従来技術で各種の方法が開発されてきた。これらの方法は、例えば、ウェーブレット変換を利用する。ウェーブレット変換の主な短所は、画像の輪郭を正確に考慮することが困難な点である。ウェーブレットが取り得る向きの数は利用する種類に依存する。いずれの場合も、基本機能が異方性を示さず、このため曲線の表現の最適特性が得られず、真のウェーブレットは多くても６個の異なる向きしか取り得ない。この個数は、自然の画像に含まれる様々な向きを最適に表現するには不十分である。
【０００５】
Ｐａｔｔｅｒｎ  Ｒｅｃｏｇｎｉｔｉｏｎ，２００８，ＩＣＰＲ，２００８，１９ｔｈ  Ｉｎｔｅｒｎａｔｉｏｎａｌ  Ｃｏｎｆｅｒｅｎｃｅ  ｏｎ，ｐｐ．１−４，ＩＥＥＥ，２００８，Ｂｉｓｍｕｔｈ  ｅｔ  ａｌ．に公開された“Ａ  ｄｅｖｉｃｅ  ｅｎｈａｎｃｉｎｇ  ａｎｄ  ｄｅｎｏｉｓｉｎｇ  ａｌｇｏｒｉｔｈｍ  ｆｏｒ  ｘ−ｒａｙ  ｃａｒｄｉａｃ  ｆｌｕｏｒｏｓｃｏｐｙ”という名称の文献では、ノイズ除去済み背景を得るためにウェーブレット変換を用いるが、これは輪郭に干渉する傾向がある。
【０００６】
Ｂｉｏｍｅｄｉｃａｌ  Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｉｎｇ  ａｎｄ  ｃｏｎｔｒｏｌ，ｖｏｌ．７，ｎ°２，ｐｐ  １７３−１７９，２０１２，Ｔｏｍｉｃ  ｅｔ  ａｌ．に公開された“Ａｄａｐｔｉｖｅ  ｓｐａｔｉｏ−ｔｅｍｐｏｒａｌ  ｄｅｎｏｉｓｉｎｇ  ｏｆ  ｆｌｕｏｒｏｓｃｏｐｉｃ  ｘ−ｒａｙ  ｓｅｑｕｅｎｃｅｓ”という名称の文献では、２次元＋時間信号を１次元信号へと索引を再構築するため、時空的アプローチを進めている。従って、輪郭を維持しながらウェーブレット変換を適用した後で、空間および時間領域の両方でノイズ除去済み信号が得られる。しかし、この方法は、索引の再構築が原因で、心臓の鼓動等の顕著な動きを考慮することは空間位置の変更を要するため不可能である。
【０００７】
他の方法は、シーケンス内の画像の冗長性を利用して、ノイズによる空間彩色を回避すべく、時間挙動の側面だけに注目してきた。Ａｕｖｒａｙ  ｅｔ  ａｌ．による文献“Ｊｏｉｎｔ  ｍｏｔｉｏｎ  ｅｓｔｉｍａｔｉｏｎ  ａｎｄ  ｌａｙｅｒ  ｓｅｇｍｅｎｔａｔｉｏｎ  ｉｎ  ｔｒａｎｓｐａｒｅｎｔ  ｉｍａｇｅ  ｓｅｑｕｅｎｃｅｓ  −  Ａｐｐｌｉｃａｔｉｏｎ  ｔｏ  ｔｒａｎｓｐａｒｅｎｔ−ｍｏｔｉｏｎ−ｃｏｍｐｅｎｓａｔｅｄ  ｎｏｉｓｅ  ｒｅｄｕｃｔｉｏｎ  ｉｎ  ｘ−ｒａｙ  ｉｍａｇｅ  ｓｅｑｕｅｎｃｅｓ”，２００８に記述されている一方法は、動きを補償しながら時間領域フィルタリングを行う。１個のピクセルが単一の動きベクトルに関連付けられているとする近似を行うことが可能である。これにより方法は簡素化されるが誤差が生じる。蛍光透視画像は増倍的であるため、身体部位別の層が重なり合っている。しかし、これらは異なる動きに従属するため、異なる経路を辿る複数の被写体が同一ピクセルに出現する可能性がある。Ａｕｖｒａｙらの発想は従って、透明層を分離して、これらの層の各々での動きベクトルを評価するものである。この方法には、取得が極めて困難である画像の各部の透明層の個数を評価する必要があるという短所がある。
【０００８】
図１に、カーブレット法による空間的ノイズ除去を実行する従来技術で公知の方法を概略的に示す。画像に２次元空間変換のステップが施され、次いで係数の相関特性が検証される。これらの特性に準拠していない場合、閾値化関数を適用する前に対象係数の振幅を修正する。方法は、例えば、閾値化関数で４個の領域を囲む３個の閾値を用いる。これらの領域を用いて係数同士のコヒーレンスを検証し、次いでノイズに起因し得る顕著な変化を検出する。現在開発中のアルゴリズムにより、ノイズ除去を行う間に輪郭を維持することができる。
【０００９】
上述の空間処理により、ノイズ除去済み画像に存在するアーチファクトの個数を減らせるが、完全に除去することができない。特に、Ｘ線量が少ないほど、アーチファクトの個数が増大する。空間冗長性が低いこと、且つ時間的制約により、処理の性能特性を向上させることはできない。にもかかわらず、処理を、時間経過に伴い画像のピクセル値を平均化する時間領域フィルタリングと合わせて用いることが可能である。フィルタは空間フィルタリングの前に適用され、両者を組み合せて用いることによりアーチファクトのないノイズ除去済み画像を得ることができる。しかし、動いている被写体を維持するために、時間領域フィルタのノイズ除去能力が低下する。従って、これを用いることでＸ線量の削減に関する可能性が制約される。
【００１０】
Ｓ．Ｎｉｇａｍ  ａｎｄ  Ａ．Ｋｈａｒｅ，Ｃｏｍｐｕｔｅｒ  Ｖｉｓｉｏｎ，ＩＥＴ，ｖｏｌ．６，ｎｏ．３，ｐｐ  ２３１−２５１，２０１２による“Ｃｕｒｖｅｌｅｔ  ｔｒａｎｓｆｏｒｍ−ｂａｓｅｄ  ｔｅｃｈｎｉｑｕｅ  ｆｏｒ  ｔｒａｃｋｉｎｇ  ｏｆ  ｍｏｖｉｎｇ  ｏｂｊｅｃｔｓ”という名称の刊行物は、カーブレット内の被写体を追跡する方法を記述している。この方法は、追跡している被写体の位置を更新するために（スケール間およびスケール内）全エネルギーの保存を利用する。無ノイズビデオについて結果が公表されており、空間的方法またはウェーブレット追跡方法よりも良好な追跡結果を示している。
【先行技術文献】
【非特許文献】
【００１１】
【非特許文献１】Ｐａｔｔｅｒｎ  Ｒｅｃｏｇｎｉｔｉｏｎ，２００８，ＩＣＰＲ，２００８，１９ｔｈ  Ｉｎｔｅｒｎａｔｉｏｎａｌ  Ｃｏｎｆｅｒｅｎｃｅ  ｏｎ，ｐｐ．１−４，ＩＥＥＥ，２００８，Ｂｉｓｍｕｔｈ  ｅｔ  ａｌ．に公開された“Ａ  ｄｅｖｉｃｅ  ｅｎｈａｎｃｉｎｇ  ａｎｄ  ｄｅｎｏｉｓｉｎｇ  ａｌｇｏｒｉｔｈｍ  ｆｏｒ  ｘ−ｒａｙ  ｃａｒｄｉａｃ  ｆｌｕｏｒｏｓｃｏｐｙ”
【非特許文献２】Ｂｉｏｍｅｄｉｃａｌ  Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｉｎｇ  ａｎｄ  ｃｏｎｔｒｏｌ，ｖｏｌ．７，ｎ°２，ｐｐ  １７３−１７９，２０１２，Ｔｏｍｉｃ  ｅｔ  ａｌ．に公開された“Ａｄａｐｔｉｖｅ  ｓｐａｔｉｏ−ｔｅｍｐｏｒａｌ  ｄｅｎｏｉｓｉｎｇ  ｏｆ  ｆｌｕｏｒｏｓｃｏｐｉｃ  ｘ−ｒａｙ  ｓｅｑｕｅｎｃｅｓ”
【非特許文献３】Ａｕｖｒａｙ  ｅｔ  ａｌ．による文献“Ｊｏｉｎｔ  ｍｏｔｉｏｎ  ｅｓｔｉｍａｔｉｏｎ  ａｎｄ  ｌａｙｅｒ  ｓｅｇｍｅｎｔａｔｉｏｎ  ｉｎ  ｔｒａｎｓｐａｒｅｎｔ  ｉｍａｇｅ  ｓｅｑｕｅｎｃｅｓ  −  Ａｐｐｌｉｃａｔｉｏｎ  ｔｏ  ｔｒａｎｓｐａｒｅｎｔ−ｍｏｔｉｏｎ−ｃｏｍｐｅｎｓａｔｅｄ  ｎｏｉｓｅ  ｒｅｄｕｃｔｉｏｎ  ｉｎ  ｘ−ｒａｙ  ｉｍａｇｅ  ｓｅｑｕｅｎｃｅｓ”，２００８
【非特許文献４】Ｓ．Ｎｉｇａｍ  ａｎｄ  Ａ．Ｋｈａｒｅ，Ｃｏｍｐｕｔｅｒ  Ｖｉｓｉｏｎ，ＩＥＴ，ｖｏｌ．６，ｎｏ．３，ｐｐ  ２３１−２５１，２０１２による“Ｃｕｒｖｅｌｅｔ  ｔｒａｎｓｆｏｒｍ−ｂａｓｅｄ  ｔｅｃｈｎｉｑｕｅ  ｆｏｒ  ｔｒａｃｋｉｎｇ  ｏｆ  ｍｏｖｉｎｇ  ｏｂｊｅｃｔｓ”
【発明の概要】
【課題を解決するための手段】
【００１２】
本発明の主題は、ノイズ減少に関して時間領域フィルタリングの性能を向上させると同時に、動いている被写体を維持することが可能な新たなアプローチに基づいている。本発明による方法は、シーケンス内の画像の離散的カーブレット変換を用いる。
【００１３】
本発明は、Ｘ線検出器により取得された蛍光透視画像のシーケンス内のノイズを減らす方法であって、画像が１個以上の特徴を含み、画像がカーブレット係数ｃｔ，ｓ，θｓｂ，ｉ，ｊにより表され、ここで、ｓは画像のスケール、θｓｂは画像内の特徴の角度に紐付けられた特徴の副帯域、および（ｉ，ｊ）は特徴の相対位置である、方法において、少なくとも以下のステップ、すなわち、
・副帯域毎にカーブレット係数を用いて、且つ副帯域θｓｂに関連付けられたフィルタの最大応答を判定することにより、カーブレット領域内の画像内の特徴を検出し、検出された特徴をメモリに保存するステップと、
・検出された特徴を局所的に照合するステップであって、
・カーブレット係数の集合を対象として、対象とするカーブレット係数が時点ｔ以前に検出された特徴に属するか否かを検証し、
・カーブレット係数が検出された特徴に属する場合、係数が属する特徴が、より早い時点で検出された特徴と一致するか否かを検証し、一致する場合、カーブレット係数を第１の集合に分類し、
・係数が検出されない場合、またはカーブレット係数が検出された特徴に属するがより早い時点で検出された特徴と一致しない場合、少なくとも第２の集合を形成する、ステップと、
・第１の集合に属するカーブレット係数の時間領域フィルタリング用に第１の係数を、および第２の集合に属するカーブレット係数の時間領域フィルタリング用に少なくとも第２の係数を適用するステップと、
・時間領域フィルタリングされたカーブレット係数を用いてノイズ除去済み画像を決定するステップと
を組み合わせて含むことを特徴とする。
【００１４】
別の変型実施形態によれば、本方法は、例えば以下のステップ、すなわち、
・副帯域毎にカーブレット係数を連結するカーブレット領域内の処理により画像内の１個または複数の特徴を検出するステップであって、
・分解において最高周波数カーブレットｓｍａｘにおけるスケールを採用し、
・各カーブレット係数ｃｔ，ｓ，θｓｂ，ｉ，ｊについて、高周波情報を含む低スケールでその親係数を計算し、
・次いでこれらのカーブレット係数を合算することにより、連結された係数ｃｃｏｎｃａｔｅｎａｔｅｄのマトリクスを取得し、
・特徴の検出結果を連結された係数のマトリクスに適用し、
・各々の連結された係数について、副帯域に関連付けられた有向フィルタの最大応答を計算し、
【数１】
0006375373.tif 000002
ここで、ｐ、ｑは向きを決定する位置であり、
・最大応答の結果を、所与の画像について維持および記録される特徴を与えるヒステリシスにより閾値化して、係数が高閾値より大きい場合、低閾値よりも高いその拡張隣接係数の全てを維持し、
・最高周波数のスケールにおける特徴の位置を選択することにより基準特徴を画定する、ステップ
を含む、特徴を検出する少なくとも１個のステップを含む。
【００１５】
一変型実施形態によれば、本方法は、例えば以下のステップ、すなわち、
・時点ｔより早い時点ｔ−１で決定されたカーブレット係数の集合を対象とするステップであって、
・カーブレット係数が時点ｔまたは時点ｔ−１で検出された特徴に属するか否かを検証し、
・検証の結果が否定的であれば、カーブレット係数は被写体に属していないものと分類し、
・検証の結果は肯定的であれば、
・係数が属する特徴が時点ｔ−１または時点ｔ−２で検出された特徴に一致するか否かを検証し、
・検証の結果が否定的であれば、カーブレット係数は被写体に属していないものと分類し、
・検証の結果が肯定的であれば、カーブレット係数は被写体に属しているものと分類する、ステップと、
・少なくとも２種類のカーブレット係数に対応する少なくとも２組の係数を取得するステップと
を含む、特徴を局所的に照合するステップを含む。
【００１６】
本方法はまた、例えば以下のステップ、すなわち、
・被写体に属していないカーブレット係数について、時点ｔで検出された係数と時点ｔ−１で検出された係数とを平均化するステップと、
・被写体に属するカーブレット係数について、より弱いフィルタリングを適用するステップと
を含む、時間領域フィルタリングのステップを含んでいてよい。
【００１７】
本方法はまた、ノイズ除去済み画像が、時間領域フィルタリング済み係数を用いて決定される逆変換ステップを含んでいてよい。
【００１８】
一実施形態によれば、時間領域フィルタリングステップは、２個の利得係数αおよびβ（ここで、０＜α＜β≦１）を定義する。
・カーブレット係数が検出されなかった場合、または検出されたが一致しなかった場合、時点ｔで利得αを係数に適用し、および時点ｔ−１で利得（１−α）を同一係数に適用し、
【数２】
0006375373.tif 000003
・係数が時点ｔまたは時点ｔ−１で検出された場合、利得βを係数に適用すると共に、その隣接係数から検出された係数までの距離の関数として所与の減少率でβからαまで減少する利得を係数の隣接係数に適用する。
【００１９】
一実施形態によれば、利得の減少率は線形または正規分布に従い適用される。
【００２０】
時点ｔで検出されたカーブレット係数が、例えば時点ｔ−１で検出され且つそれらと一致するカーブレット係数の集合と平均化される。
【００２１】
一変型実施形態によれば、前の時点で特徴に属していた係数、または現時点で一致した特徴に属する係数、あるいはこれら２種類の係数のいずれかの平均化において、時点ｔでノイズ除去されたカーブレット係数をノイズカーブレット係数の代わりに用いる。
【００２２】
同一被写体に属する特徴を照合すべく、本方法は例えば以下のステップ、すなわち、
・コスト関数を定義し、ならびに特徴の長さ、強度、変位、および角度を類似度指標として使用するステップと、
・全エネルギーを次式で定義するステップであって、
【数３】
0006375373.tif 000004
ｗｌ、ｗｉ、ｗｄおよびｗａは指標の信頼度に対応し、ここで、ｗｌ、ｗｉ、ｗｄおよびｗａ＝１である、ステップと、
・全エネルギーを、一方の特徴が時点ｔで、他方が時点ｔ−１で検出された特徴の全てのペアについて計算するステップと、
・エネルギーの重みを最小化する照合リストを取得すべく適応された方法を適用するステップと
が実行され、
・カーブレット係数が特徴に対応していない場合、方法は時間平均化ステップを実行する。
【００２３】
比較される最小の特徴の長さＬｍｉｎおよび比較される最大の特徴の長さＬｍａｘが例えば以下のように定義され、
【数４】
0006375373.tif 000005
ここで、Ｉｍｉｎは２個の特徴のうち低い方の強度であり、Ｉｍａｘは他方の特徴の強度であり、
【数５】
0006375373.tif 000006
指数ＩおよびＥは時点ｔおよびｔ−１における特徴の極値を表し、
【数６】
0006375373.tif 000007
ここで、θは時点ｔおよびｔ−１で検出された特徴の向きを表す。
【００２４】
本方法は、｜ξ×β｜より大きい利得を有するカーブレット係数に係数ρを乗算することにより、画像のコントラストを高めるステップを含んでいてよく、１＜ρであり、ここで、ξはα＜ξβ＜βであるように選択される。
【００２５】
一変型実施形態によれば、乗数係数ρを、スケール内近傍でρから１まで減らしながら、検出されたカーブレット係数に適用する。
【００２６】
本発明による方法の他の特徴および利点は、添付図面に非限定的に示す例示的な一実施形態に続く記述を精査すればより明確に明らかになろう。
【図面の簡単な説明】
【００２７】
【図１】従来技術による方法の概略的なレイアウトを示す。
【図２】本発明の方法の概略的なレイアウトを示す。
【図３】本方法を実装可能なステップを示す。
【図４】画像の取得および取得された画像の処理を可能にするシステムの説明図を示す。
【発明を実施するための形態】
【００２８】
以下の記述は、例示的であって一切限定的ではない。図２に、図１に示すような従来技術で公知の原理に従い空間ノイズ除去により得られたノイズ除去済み係数を用いて適用した時間領域フィルタリングの原理を概略的に示す。
【００２９】
本方法は、シーケンス内の画像の離散的カーブレット変換を用いる。ノイズ除去済み係数が検出され（２０）、次いで局所的に特徴を照合するステップへ送信され（２１）、以下にその詳細を述べる。一致したデータは次いで、画像内で被写体の動きを維持可能にしながら時間領域フィルタリングのステップ２２および最終画像を生成すべく２次元逆離散的カーブレット変換（ＩＤＣＴ）ステップ２３を施される。
【００３０】
上記カーブレットにおいて、画像内に存在する２次元不連続箇所等の全ての曲線状の要素は区間の形状で表され、これらの要素により画像の再構築が可能になる。
【００３１】
以下の説明において、Ｘ線検出器により取得された一連の蛍光透視画像ｘｔを対象とし、ｔはシーケンス内での取得時点を表し、Ｘ線検出器は取得されたデータを適当なデータ処理装置に送信する。各画像は、各々が信号レベルを表す値ｘｔ（ｐ，ｑ）を有するピクセルのマトリクスからなり、ｐおよびｑは各々マトリクスの行および列を表す。値は、例えば１６ビットに符号化され、０〜６５５３５の範囲にある。コントラスト対ノイズ比は以下のように定義される。
ＣＮＲ＝（被写体のグレイレベル−背景のグレイレベル）／σ
σは標準偏差として定義され、ノイズは好適に、関心対象の被写体が存在しない画像の領域にわたり決定され、コントラストは当業者に公知の方法を用いて関心対象の被写体に対応する点に基づいて決定される。
【００３２】
基準カーブレット係数とそれに紐付けされた係数との間で確立可能な親子関係と呼ばれる紐付けを以下に定義する。このカテゴリにおいて、基準係数に関して定義される親係数、隣接係数および従弟係数を区別することができる。親係数は、基準カーブレット係数と同一相対位置にあるが、周波数が低い次スケールに位置する係数として定義される。隣接係数は、カーブレット係数の位置に接続された位置にある係数に対応する。カーブレットの方向にある隣接係数は、整列隣接係数と呼ばれ、非整列隣接係数と呼ばれる他の係数と区別される。従弟係数は、基準カーブレット係数と同一相対位置にあるが、向きが異なる。向きが隣接している場合、用語「隣接従弟」を用いる。向きが係数の向きと反対である場合、用語「逆向き従弟」を用いる。
【００３３】
図３に、本発明による方法が実行する一連のステップの一例を示す。
【００３４】
本発明による方法の第１のステップにおいて、Ｘ線検出器により生成されて画像処理装置へ送信されたノイズ画像ｘｔは、ｔ（３０）で、空間ノイズ除去により得られるノイズ除去済み係数を取得すべく空間フィルタリングされる（３１）。
【００３５】
【数７】
0006375373.tif 000008
を時点ｔで画像のカーブレット変換の離散的係数とする。以下において、表記を煩雑にしないよう指数Ｄは含まれない。カーブレット係数は、画像の解像度を表す係数のスケールｓ（低周波数対高周波）、画像内の特徴の角度に紐付けられた係数の副帯域θｓｂ、および係数が属する画像のカーブレット表現の副帯域内の係数の座標であるトルク（ｉ，ｊ）より表される係数の相対位置により識別される。そのような係数またはカーブレット係数をｃｔ，ｓ，θｓｂ，ｉ，ｊと表記する。
【００３６】
図３に示すステップは、親３２および相関する隣接係数３３を発見して、これらの親または相関する隣接係数が同一領域に存在し、且つ潜在的にステップ３４で係数Ｃｔを修正すべきか否かを判定可能にするものであり、本出願人により出願されて当業者には公知の仏国特許出願公開第２９７８２７３号明細書に記述されている。
【００３７】
閾値化ステップ３５は、取得された画像について、空間的にノイズ除去済みカーブレット係数
【数８】
0006375373.tif 000009
を生成する。閾値化関数３５は、例えば２個のステップからなる。第１段階として、カーブレット係数の振幅とカーブレット係数に紐付けられた係数とのコヒーレンスが上で定義されたように検証され、必要であればコヒーレンスは再構築される。第２段階において、閾値化関数が適用される。
【００３８】
閾値化関数は、例えば、（初期ノイズ画像について計算された）標準偏差に依存する３個の閾値Ｔ１、Ｔ２およびＴ３（Ｔ１＜Ｔ２＜Ｔ３）の関数として定義され、以下のように書くことができる。
【数９】
0006375373.tif 000010
γおよびλは、関数がそれ自体の定義域にわたり連続的であるように計算されたパラメータである。
【００３９】
隣接および親が対象カーブレット係数と相関するため、後述するように情報のコヒーレンスを検証することが可能である。
【００４０】
空間閾値化ステップ３５は、例えば、対象係数と同一の関数領域で発見された係数に対して実行される。
【００４１】
親係数ｃｐａｒｅｎｔは、最低周波数のスケールで、同じ向きを表す副帯域内のカーブレット係数ｃｔ，ｓ，θｓｂ，ｉ，ｊと同一相対位置にあるカーブレット係数である。ここでｓ−１は、カーブレットのスケールｓに関して低周波数のスケールを表す。向きの個数が少なければ、
【数１０】
0006375373.tif 000011
さもなければ同一副帯域
【数１１】
0006375373.tif 000012
が保たれる。
【００４２】
同様に、ｍｓ，θｓｂおよびｎｓ，θｓｂでスケールｓにおける副帯域の寸法を表せば、親係数の新しい位置は以下のように導かれる。
【数１２】
0006375373.tif 000013
【００４３】
次いで、カーブレット係数｜ｃｔ，ｓθｓｂ’ｉ，ｊ｜が隣接カーブレット係数｜ｃｎｅｉｇｈｂｏｒｓ｜および親カーブレット係数
【数１３】
0006375373.tif 000014
と同一の閾値の間に位置していることが、異なる領域に属する全てのケースを処理しながら検証される。以下の部分において、表記を簡素化すべく、明らかに読み取れる場合は添え字を割愛し、親の表記または隣接の指示だけに添え字を付与する。
【００４４】
空間閾値化ステップ３５の後で、所与の画像に対応するノイズ除去済みカーブレット係数が得られる。
【００４５】
本発明による方法は、画像をカーブレットに空間フィルタリング分解した後で特徴を検出するステップを実行する。従来のフィルタリングによる分解では特に、検出実行時の偽陽性を減らすことができる。本発明の範囲から逸脱することなく、ノイズ分解を利用することも可能であり、後述するステップは、非フィルタリングカーブレット係数にも適用されてよい。
【００４６】
画像内の特徴の検出は、例えば、カーブレットへの分解のスケールの連結３６を用いて実行される。これは、画像のカーブレット係数を副帯域毎に連結することから始まる。換言すれば、最高周波数カーブレットでのスケールが分解において採用される。これをｓｍａｘと呼ぶ。各カーブレット係数ｃｔ，ｓ、θｓｂ，ｉ，ｊについて、その親は、関心対象の（換言すれば、高周波情報を含む）低周波数スケールで計算される。本例では、親は２個のスケールｕＭａｘ＝２だけについて計算される。これらの親係数は次いで合算され、従って次式が得られる。
【数１４】
0006375373.tif 000015
ここで、ｕＭａｘは、対象とするスケールの個数、ｓｍａｘ−ｕは、表現のスケールの添え字である。
【００４７】
連結された係数は従って時点ｔおよび副帯域θｓｂ（添え字がｓｍａｘの）だけに依存する。連結された係数のこれらのマトリクスに行検出が適用される。第１のステップは、各連結された係数について、副帯域θｓｂに関連付けられた有向フィルタの最大応答３７を、有向線形ベクトルとの畳み込みにより計算することである。
【数１５】
0006375373.tif 000016
フィルタは、有向フィルタであって、当業者に公知の方法に応じて、例えば方法の最初で計算されるか、または事前に知られている。副帯域に関連付けられた有向ベクトルの個数は、分解のパラメータ、およびベクトルのサイズに依存する。典型的に、低周波数スケールで６個のスケールおよび１２個の角度を有するカーブレットに分解する場合、小さい被写体を検出するために、例えば長さが５であるベクトルを用いる。これらの条件下で、全ての向きをカバーするのに１副帯域当たり１または２個のベクトルで十分である。
【００４８】
畳み込みの結果は次いで、例えば閾値化ヒステリシス３８により、行を識別するために分割される。このように得られた特徴は、画像同士での比較用に処理装置のメモリに保存される。ヒステリシスによる閾値化の利点は、偽陽性を制限して区間の細切れ検出を回避することである。
【００４９】
本方法は、連結された係数の全てのマトリクスに対して全体的な閾値、すなわち高閾値Ｔｈｉｇｈ＝αＴｏｔｓｕおよび低閾値Ｔｌｏｗ＝βＴｏｔｓｕを（α＞β）選択する。この目的のため、連結副帯域の最大値のヒストグラムが計算される。いくつかの副帯域がノイズしか含んでいない場合、ヒストグラムは低い値で「ピーク」を示す。これらの値では更に、他の最大値が示されている。全体的な閾値Ｔｏｔｓｕは、当業者には公知のＯｔｓｕの方法をヒストグラムに適用して得られた値により与えられ、この値によりヒストグラムを２分割することができる。これらの閾値を選択する別の方法は、ノイズの標準偏差の関数として定義されるものである。従って、Ｔｈｉｇｈ＝ασｓｂおよびＴｌｏｗ＝βσｓｂが得られよう。
【００５０】
次のステップは、特徴を照合するステップ４０である。特徴の検出により、各スケールでの位置が公知である特徴の集合を取得することができる。これらの位置がスケール間で紐付けられている前提で、基準対象と呼ばれる特徴の位置、換言すれば、最高周波数スケールでの特徴の位置だけがここで関心対象である。照合は、２個のステップを含んでいる。すなわちコスト関数の定義および関数の最適化である。当業者には公知であるいくつかの最適化方法が存在する。
【００５１】
後述する例において、選択される最適化方法は、Ｋｕｈｎ−Ｍｕｎｋｒｅｓ法である。コスト関数として、カーブレット領域で重要と判断される量を考慮する関数を用いる。
【００５２】
画像が実時間（毎秒１５〜３０画像）で取得され、且つこれらが身体部位別の画像である前提で、検出される輪郭は２個の画像間でほとんど修正されないものと思われる。例えば、類似度指標として、以下の指標を選択する。すなわち特徴の長さ、強度、動き、角度である。ある時間にわたり特徴を照合することを想定しているため、時点ｔ−１で検出された特徴（３９）は従って、時点ｔ−１で検出されて例えばシステムに保存された係数の値を用いて、時点ｔで検出された特徴と比較される（４０）。「時点ｔ−１」という表現は、本明細書でより広義に、対象時点ｔよりも早い時点（ｔ−１）、（ｔ−２）、（ｔ−３）等を表す。
【００５３】
特徴の照合を行うために実行されるステップは以下の通りである。カーブレット係数が、時点ｔで検出された特徴、または時点ｔ−１で検出された特徴に属する否かを検証する（４１）。カーブレット係数が検出されたいずれの特徴にも属していない（４２）場合、本方法は次いで、係数が被写体ではなく画像の背景に属すると宣言する。カーブレット係数が検出された特徴に属する（４３）場合、本方法は次いで、係数が属する特徴が、前の時点ｔ−１または時点ｔ−２で検出された特徴と一致するか否かを検証（４４）する。応答が否定的（４５）である場合、本方法は特徴を画像の背景に属すると見なす。その逆の場合、検出されて一致した特徴を構成する係数は、被写体の一部を形成すると見なされる（４６）。これらのステップの終了時点で、本方法は、カーブレット係数および一致した特徴の複数の集合を決定して、カーブレット係数の種類に応じた時間領域フィルタリングのステップ４７を適用する。第１の集合は、例えば、検出されて一致したカーブレット係数を含み、第２の集合は、検知されなかったか、または検出されたが一致しないカーブレット係数を含んでいてよい。第２の集合を２個の部分集合に分割して、第１の部分集合は検知されなかった係数を、第２の部分集合は検出されたが一致しない係数を含むようにしてもよい。この場合、時間領域フィルタリングステップの係数は異なるであろう。これは、例えば検知されたが一致しない係数に対し、検出されなかった係数よりも強いフィルタリングを適用するとの発想である。
【００５４】
特徴に対応するエネルギーを以下で定義する。Ｌｍｉｎは、比較される最小の特徴の長さ、Ｌｍａｘは比較される最大の特徴の長さを表す。
【数１６】
0006375373.tif 000017
同様に、Ｉｍｉｎは２個の特徴のうち弱い方の強度、Ｉｍａｘは他方の強度を表す。
【数１７】
0006375373.tif 000018
全エネルギーの最小化、従って全て同程度の大きさである最小化対象エネルギーの決定が求められている。特徴が伸縮してもエネルギーは一定であるように特徴の長さが順序付けられる。上の例において、これらの可能性は共に均等である。強度についても同じことが言える。従って、変化がなければこれらのエネルギーは０に近く、特徴の特性に顕著な変化があれば１に近い。
【００５５】
特徴の動きに関して、区間の両端間の、例えばマンハッタンの意味での距離を測定する。換言すれば、時点ｔと時点ｔ−１で検出された特徴の第１終端の間、および時点ｔと時点ｔ−１で検出された特徴の他方の端の間の距離を測定する。これにより、中心回りの回転の場合に距離がゼロになることが回避される。更に、特徴の動きが大きいと判断される場合に、距離を値１近くまで戻すために、各々結果を副帯域の寸法により正規化する。このように、１本の軸に沿った距離にμ／ｍを乗算し、ここで、ｍは寸法（他方のケースではｎで除算する）、μは受容可能な動きを固定させるパラメータである。例えば、パラメータが４に等しい場合、寸法の１／４の動きが許される。この値を超えると、エネルギーは１より大きくなり、全エネルギーに対して特徴が一致しないという顕著な悪影響を及ぼす。これを保証すべく、動きのエネルギーが１より大きい場合、一致を不可能にすべく全エネルギーに対し最大値が設定される。
【数１８】
0006375373.tif 000019
指数ＩおよびＥは、時点ｔおよびｔ−１での特徴の極値を表す。最後に、特徴の角度に紐付けされたエネルギーを次式で定義する。
【数１９】
0006375373.tif 000020
但し「θ」が時点ｔおよびｔ−１での特徴の向きを表す。これは従って、悪影響を受ける、πを法とする変動である。向きの変動が大き過ぎる場合、このエネルギーにパラメータνを乗算することで１に近くなる。ここではラジアンを用いているため、１０分の１のオーダーの変動は大きいと判断される。従って、このパラメータの値が１０であれば変動が０．１ラジアンに等しい場合にエネルギーを１に近くすることができる。パラメータは、維持しようとする動きの振幅の関数として調整することができる。
【００５６】
従って、全エネルギーは次式のように定義される。
【数２０】
0006375373.tif 000021
係数ｗｌ、ｗｉ、ｗｄ、ｗαは、指標に対する信頼度に対応する重みである。これらの４個の重みの合計は１に等しくなければならない。例えば、長さに紐付けされたエネルギーに関連付けられた重みは、角度に関連付けられた重みよりも低いが、その理由は、（例えば検出が僅かに動くため）特徴の極値位置が僅かに変動すると、特定のケースにおいて長さ（例えば、特徴の長さが短い場合）の相当大きい変動が結果的に生じる一方、角度はこれらの小さい動きにさほど影響を受けないからである。しかし、長さは情報を提供するため、長さを維持することは有利である。
【００５７】
この全エネルギーＥｔｏｔａｌは、特徴のペアの一方が時点ｔで検出され他方が時点ｔ−１で検出された特徴の全てのペアについて計算される。不要な計算を抑制すべく、本方法は隣接副帯域内の特徴のエネルギー計算に限定される。これは、特徴が広く分離された副帯域に存在する場合、それらの角度が大幅に異なって、極めて高いエネルギーが生じるからである。
【００５８】
最後に、上述のように、最大エネルギーＥｍａｘが決定される。特徴が最大エネルギーに達するかまたは超えた場合、照合は不可能である。
【００５９】
全てのエネルギーが計算されると、Ｋｕｈｎ−Ｍｕｎｋｒｅｓ法が適用されて、エネルギーの重みを最小化する照合リストを得ることができる。一致した特徴が実際の信号であると見なされる。不一致特徴に関する限り、これらは値が異常である係数の集合を含んでいる。これらは一般に「外れ値」と呼ばれる。これは用語の誤用であるが、これらの不一致特徴もまた「外れ値」と表記するものとする。
【００６０】
以下のステップは、係数の種類に応じて適合された時間領域フィルタリング４７である。時間領域フィルタリング４７は、係数の種類に応じて異なるフィルタリングを適用すべく、カーブレット係数の種類、すなわち検出されて一致した係数、またはノイズと見なされる係数に関して以前に得られた情報を用いる。この方法により、動きがない領域内で強いフィルタリングを適用し、関心対象の特徴を表す領域内で弱いフィルタリングを適用することができる。
【００６１】
【数２１】
0006375373.tif 000022
を、係数のクラスに依存する時点ｔでの時間領域フィルタリングの利得とする。時間領域フィルタリングは次式のように書くことができる。
【数２２】
0006375373.tif 000023
【００６２】
一実装方法では、例えば、２個の利得係数αおよびβ（０＜α＜β≦１）が関連付けられた２個のクラスの係数を使用する。
【００６３】
カーブレット係数が検出されなかったか、または検出されたが一致しなかった場合、第２のクラスに属する。係数は異常な特徴の一部をなすと見なされ、利得αが係数に適用され、同一係数に利得（１−α）が時点ｔ−１で適用される。
【数２３】
0006375373.tif 000024
【００６４】
係数が時点ｔまたは時点ｔ−１で検出され、且つ係数が属する特徴が一致した場合、利得βが係数に適用され、その隣接係数には隣接係数から検出された係数までの距離の関数としてβからαまで減少する利得が適用され、減少率は例えば線形または正規分布に従うものであってよい。
【数２４】
0006375373.tif 000025
【００６５】
減少率はスケール内近傍で適用される。利得β’（α＜β’≦β）から出発して同一方法が同一スケールの全ての副帯域に対して適用される。再構築中の反跳アーチファクトを回避すべくランプを用いる。検出がなかった低周波数のスケールであっても利得の適用を選択してよい。この場合、それ自体の子係数のうち１個が検出されて一致した係数のコントラストが高まる。衝突（ランプおよび間引きに起因して同一係数に各種の利得が適用され得るケース）の場合、最も高い利得が維持される。
【００６６】
β＝１の場合、本方法により、関心対象の被写体のコントラストを維持しながらノイズを減らすことができ、且つ遅延アーチファクトの生成が回避される。
【００６７】
第２の実施形態において、時点ｔで検出されたカーブレット係数は、時点ｔ−１で検出され且つそれ自体と一致するカーブレット係数の集合と平均化される。
【００６８】
第３の方法において、前の時点である特徴に属していた係数、または現時点で一致した特徴に属する係数、あるいはこれら２種類の係数の平均化において、ノイズ係数｜ｃｔ｜の代わりに時点ｔでのノイズ除去済み係数
【数２５】
0006375373.tif 000026
を用いる。
【００６９】
一変型実施形態によれば、検出された要素を用いてコントラストを高めることが可能である。このステップは時間領域フィルタリングステップの前後で実行されてよい。
【００７０】
一方法は、例えば、利得が｜ξ×β｜よりも高い全てのカーブレット係数に係数ρ（ρ＞１）を乗算するものである。これにより、検出された係数を隣接副帯域に存在する係数と共に増大させることができ、この場合も再構築中に輪郭の付近での反跳の問題が避けられる。係数ξはα＜ξβ＜βであるように選択されなければならない。「外れ値」は照合ステップにより除外されており、従ってコントラストの増大は輪郭だけに適用される。
【００７１】
時間領域フィルタリングの係数を考慮しない別の実施形態は、時間領域フィルタリング方法と同様の仕方でコントラスト強調係数を適用するものである。乗数係数ρを、スケール内近傍でρから１まで減らしながら、検出された要素に適用する。同一スケールの副帯域に利得ρ’（１＜ρ’≦ρ）を用いて同一方法が適用される。本方法はスケール全体にわたり適用されてもよい。スケールに依存する利得ρを用いることも可能であり、これにより局所的にノイズを増やすことなくコントラストを更に高めることができる。
【００７２】
本方法の最後のステップ４８は、ノイズ除去済み画像４９を再構成すべく２次元ＩＤＣＴ変換を適用するものである。
【００７３】
図４に、本発明による方法を実装可能にするシステムの一例を示す。システムは、画像を検出するＸ線検出器５０を含み、検出器は、本発明による方法のステップを実行すべく設計されているプロセッサ５２を含む処理装置５１、および変換されたカーブレットを保存するメモリ５３に接続されている。
【００７４】
本発明による方法は特に、画像の透明度の問題を特に解決可能にするカーブレットの照合に基づいている。実際、Ｘ線物理学により、各種の身体部位別の層が画像内の同一平面に表され、１個のピクセルで様々な層を表すことができる。これに関連して、１個の動きベクトルに各ピクセルを関連付けることができないため、且つこれらの層が異なる動きを示すため再位置決めを正確に行うことができず、１個のピクセルに存在する被写体は異なる動きを示す。カーブレットはそれ自体の向きに応じて要素を分離し、微細要素から背景を分離するため、微細な係数に対して、第１の近似として、各係数に単一の動きベクトルを関連付けることができると仮定することが可能である。
【図１】
0006375373.tif 000036
【図２】
0006375373.tif 000037
【図３】
0006375373.tif 000038
【図４】
0006375373.tif 000039