(57)【要約】      （修正有）【課題】  Ｘ線画像における散乱防止格子の存在によって生ずる周期パターンを抑制する。【解決手段】  画像信号を周期的パターンに相当する空間周波数成分の表現に変換するために、画像信号に短時間ゲルツェル変換が適用される。この表現が処理される。次いで、変換の逆が計算され、そして逆変換の結果が画像信号から除去される。【選択図】  図２
2006289095.tif 000002
【特許請求の範囲】
【請求項１】
画像信号内に生ずる周期的パターンに相当する空間周波数成分が抑制される周期的パターン抑制方法であって、
−前記画像信号から、画像信号を空間周波数成分表現に変換するように画像信号に変換
を適用することにより、その中に含まれる周期的パターンに相当する空間周波数成分
を抽出し、
−空間周波数成分表現を処理し、
−空間周波数成分表現について逆変換を計算し、
−前記画像信号から前記逆変換の結果を除去することにより、画像信号内に生じている
空間周波数成分を抑制する
諸段階を含む方法。
【請求項２】
前記変換が短時間ゲルツェル変換である請求項１による方法。
【請求項３】
前記処理が、空間周波数成分表現への非線形高周波減衰フィルターの適用を含む請求項１又は２による方法。
【請求項４】
高周波減衰フィルターが中央値フィルターにより実行される請求項３により方法。
【請求項５】
高周波減衰フィルターが２次元フィルターである請求項３による方法。
【請求項６】
高周波減衰フィルターが、空間周波数成分と同じ方向の１次元フィルターである請求項３による方法。
【請求項７】
高周波減衰フィルターが、空間周波数成分に直角方向の１次元フィルターである請求項３による方法。
【請求項８】
高周波減衰フィルターが、いずれも１次元又は２次元である高周波減衰フィルターの連続である請求項３による方法。
【請求項９】
フィルターが複素表現の虚の部分に適用される請求項３による方法。
【請求項１０】
フィルターが複素表現の実の部分に適用される請求項３による方法。
【請求項１１】
フィルターが複素表現の大きさに適用される請求項３による方法。
【請求項１２】
フィルターが複素表現の位相に適用される請求項３による方法。
【請求項１３】
変換が、画像信号の濾波されたバージョンに適用される先行請求項のいずれかによる方法。
【請求項１４】
画像信号の濾波されたバージョンがＤＣ成分を含まない請求項１３による方法。
【請求項１５】
画像信号の濾波されたバージョンが、抑制すべき周波数のある帯域幅内の周波数のみを含む請求項１３による方法。
【請求項１６】
放射線撮影画像内の散乱防止格子線の画像に相当する空間周波数成分が除去される請求項１による方法。
【請求項１７】
コンピューター上で走るとき、先行請求項のいずれかの段階を実行するようにされたコンピュータープログラム製品。
【請求項１８】
先行請求項のいずれかの段階を実行するようにされたコンピューターで実行可能なプログラムコードを含むコンピューターによる読取り可能なキャリヤー媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、所与の周波数を有する周期的パターンの抑制方法に関する。より特別には、Ｘ線画像における散乱防止格子の存在より生ずる周期的パターンの抑制に関する。
【背景技術】
【０００２】
コンピューター放射線撮影、デジタル放射線撮影並びに古典的なフィルム使用のＸ線システムにおける散乱Ｘ線の量を減らすために通常使用される技法は、散乱防止格子の使用である。これら格子は、適切な視野確保のためある距離に離して置かれた鉛の箔のストリップである。
【０００３】
種々の形式の散乱防止格子がある。平行格子においては、鉛箔のストリップは平行であり、一方、ハニカム格子においてはストリップはハニカムパターンに置かれる。格子は静止し、又は動いている。これら格子の使用は放射線の散乱を効果的に減らすが、画像内に格子線のようなアーティファクトを作ることがある。
【０００４】
動いているシステムにおいては、格子の運動が画像内の格子線を取り除く。しかし、ある環境、例えば暴露時間が短い場合又はシステムの不具合の場合は、画像内にアーティファクトが残る。
【０００５】
画像がデジタルで形成された場合又は後でデジタル画像表現に変換された場合、ある尺度で画像を表示するときモアレアーティファクトが現れることがある。これら余り頻繁でないモアレアーティファクトは、妨害が最も大きくこれを無くす必要がある。画像内に格子線がある場合は、画像を表示する前にこれを除去すべきである。
【０００６】
診断情報を失うことなく、ズーミングされた画像の周期的アーティファクト、例えばモアレアーティファクトの負の影響を無くすことが本発明の目的である。
【０００７】
これまで特許文献１「放射線撮影の処理及び記憶のための方法、システム及びプログラム」、特許文献２「周期的パターンの抑制処理の方法及び装置」、及び非特許文献１「機能的構造を使用した離散コサイン変換の具体化のための直接公式化」が知られている。
【特許文献１】欧州特許第１  ２６５  １９４号  明細書
【特許文献２】米国特許第２００３／０９１２４３号  明細書
【特許文献３】欧州特許第０４  １０２  １８５．８号  明細書
【非特許文献１】チャウ・エル・ピー他（Chau L.P. et al）著「ダイレクト  フォーミュレーション  フォア  ザ  リアリゼーション  オブ  ディスクリート  コサイン  トランスフォーム  ユージング  リカーシブ  ストラクチュア（Direct formulation for the realization of descrete cosine transform using recursive structure)」、ＩＥＥＥ  トランズアクションズ  オン  サーキッツ  アンド  システムズ  ＩＩ：アナログ  アンド  デジタル  シグナル  プロセッシング（IEEE Transactions on circuits and systems II: analog and digital signal processing)，ＩＥＥＥ  インク.（IEEE Inc.）、ニューヨーク  米国（New York, US）
【発明の開示】
【課題を解決するための手段】
【０００８】
上述された態様は請求項１に説明された方法により実現される。
【０００９】
本発明の好ましい実施例についての特別な特徴は従属請求項において説明される。
【００１０】
本発明は、与えられた方向、画像の軸線と平行な方向における画像内の与えられた周波数を中心とする周期的変動を抑制するための方法を提供する。
【００１１】
この方法は、周期的変動に対して平行及び直角な方向の非線形高周波減衰フィルターを適用するための技法を提供する。
【００１２】
この非線形高周波減衰フィルターは中央値フィルターとすることができる。このフィルターは、オリジナル画像から最少量の関連情報が除去されることを確保する。このフィルターは、画像の画質を落とさない。
【００１３】
この方法により、変動周波数が正確には知られない周期的パターンの除去ができる。周期的変動の適切な表現を得るために短時間ゲルツェル（Geortzel)変換を行う。短時間ゲルツェル変換の結果は実部分又は虚部分として示すことができる。この表現又はこの表現からの誘導像について非線形高周波減衰濾波処理が行われる。短時間ゲルツェル変換から特徴的内容を除去した後、逆変換が計算され、そしてこの結果が入力データの補正用に使用される。
【００１４】
本発明の方法の具体例は、一般に、コンピューター上で走ったときに本発明の方法の諸段階を実行するようにされたコンピュータープログラム製品の形式で実施される。このコンピュータープログラム製品は、通常、ＣＤ−ＲＯＭ又はＤＶＤのようなコンピューター読出し可能なキャリヤー媒体に記憶される。或いは、コンピュータープログラム製品は電気信号の形式を取り、そして電子的通信を通して使用者に通信することができる。
【発明を実施するための最良の形態】
【００１５】
散乱防止格子の存在は、入力信号の周期的擾乱としてそれ自体を明らかにする。本発明による方法においては、この擾乱は純粋に周期的な性質のものであるとする。
【００１６】
周期的なアーティファクト、例えばズーミングされた画像のモアレアーティファクトの負の影響を、診断情報を除去することなしに無くすことが本発明の目的である。
【００１７】
周期的擾乱の主周波数を有するサイン波状信号の除去が、この目的に対して十分に適することが実験により示されている。
【００１８】
本発明は、加法的又は乗法的な性質の周期的アーティファクトを補正するために使用することができる。第２の事例においては、処理前に対数変換を使用し更に周期的アーティファクトの除去後に指数的変換を使用することにより、状況を加法的な性質に変換することができる（図１参照）。
【００１９】
本発明は、アーティファクトを含んだ画像に、中央値又は高周波を減衰するフィルターを周期的信号の方向で適用することにより画像を処理し、周期的アーティファクトの推定／抽出の後、周期的信号内にまだある診断情報を除去する諸段階を含む。通常、かかる高周波減衰フィルターは、信号の周期的な情報を除去する。
【００２０】
特許文献３において、この問題を克服するための技法が説明されている。この技法は、変動の周期が整数である場合にだけ作動する。
【００２１】
これは散乱防止格子の周期的擾乱に対する事例ではないため、周期的信号は与えられた周波数の複素表現又は振幅／位相表現に変換される。
【００２２】
かかる技法は一般にゲルツェル変換として知られる。
【００２３】
もし  s＝(s０,s１,…,sＮ）が入力信号であるならば、周波数ωについてのゲルツェル変換は次のように定義される。
【００２４】
【数１】
2006289095.tif 000003
簡単のため、これらの式に対して積分記号のみを使用することとする。これらの式は容易に我々の離散した問題に転用させることができる。しかし、散乱防止格子の存在のためＸ線画像内に生じた周期的パターンの抑制に対しては、周波数ωはある程度の精度で知られるだけである。更に、散乱防止格子の性質のため、この周波数は距離にわたって小ビット変化することが有り得る。本発明により、定位置に置かれたウインドウ付きのゲルツェル変換が計算される。この変換は、短時間ゲルツェル変換Γと呼ばれるであろう。
【００２５】
【数２】
2006289095.tif 000004
与えられた周波数ωに対する式１の結果は、各位置ｙに対する複素数である。複素数の逆変換が周期的擾乱の良好な推定である。
【００２６】
しかし、我々はこの複素表現（又はこの表現の適宜の部分又は誘導像、例えば虚の部分、振幅、位相…）について中央値又は高周波減衰フィルターを作用させる。このフィルターは、２次元フィルター、又は２個の分離したフィルターであって１個は格子の方向と平行であり１個は格子に直角なフィルターより構成することができる（図２及び図３参照）。
【００２７】
格子に平行な高周波減衰フィルターが既に適用された計画に続く場合は、周期的信号の確かな推定のためには格子に直角方向のフィルターだけで十分である（図４）。
【００２８】
入力信号ｓを
s(x)＝α(x)cos(ωx)+β(x)sin(ωx)
として定義できるとすれば、これは散乱防止格子のモデルのための合理的な仮定である。我々は、格子のゆっくりと変化している振幅α、βを抽出することを希望する。式（１）を次の２個の部分に分割する。
【００２９】
Γｃ(y,ω)＝∫w(x-y)s(x)cos(x)dx
Γｓ(y,ω)＝∫w(x-y)s(x)sin(x)dx
式（１）を  Γ(y、ω)＝Γｃ(y,ω)+iΓｓ(y,ω)  に変換する。我々のウインドウの広がりに亙って振幅がほぼ一致であるとすると、次の記号を使うことができる。
【００３０】
Γｃ(y,ω)≒α(y)Γｃｃ(y,ω)+β(y)Γｃｓ(y,ω)
Γｓ(y,ω)≒α(y)Γｃｓ(y,ω)+β(y)Γｓｓ(y,ω)
ここに、
Γｃｃ(y,ω)＝∫w(x-y)cos(x)cos(x)dx
Γｃｓ(y,ω)＝∫w(x-y)cos(x)sin(x)dx
Γｓｃｃ(y,ω)＝∫w(x-y)cos(x)cos(x)dx
である。この知識により、振幅を近似することができる。
【００３１】
【数３】
2006289095.tif 000005
振幅表現について何かの後処理が行われた場合は、新しい振幅表現（α’、β’）が得られ、再構築された信号ｓ’が次式で与えられる。
【００３２】
s’(x)＝α’(x)cos(ωx)+β’(x_sin(ω)
この変換の長所は倍増される。計算するためにただ１個のスペクトル成分を必要とするだけであり、更にこの変換はゆっくり変動する出力を作り、このことは、これをオリジナル入力ドメインに戻す変換前の処理にとって理想的である。
【００３３】
３．０１１画素の周期を有する合成格子により擾乱を受け小ビットのノイズを有する５０の周期の入力信号に対して（図５）、フーリエドメインにおける３個のスペクトル成分が抑制される方法と我々の方法とを比較する。周波数
【００３４】
【数４】
2006289095.tif 000006
に対する短時間ゲルツェル変換を計算するために使用されるウインドウが図６に示される。図７において、短時間ゲルツェル変換の出力が、複素記号の両部分として、及び振幅−位相記号として示される。両者の表現は極めて滑らかであり、これが残存している診断情報を除去するための高周波減衰フィルターの適用を容易にする。図７は、入力信号のフーリエスペクトルも表示する。図６は、両方法の結果が類似であることを示す。
【００３５】
信号の周期が３．０１１でなくて３．０２であり、同じ周波数の外側限度がフーリエドメインに入るとすると、別の結果が出る（図８参照）。短時間ゲルツェル変換はフーリエ法よりもうまく作動する。フーリエ技法により同じ結果を達成できるがより多くのスペクトル成分を抑制する必要がある。図９においては、７個のスペクトル成分が抑制される。これはフーリエ変換により多くの時間を消費する。
【００３６】
この方法は、周期が変化する信号を取るとき、より興味がある。図１０においては、格子の周期が２．９１１から３．１１１画素に変化する入力信号が作られる。７個のスペクトル成分を抑制することにより、短時間ゲルツェル変換と同じ結果が得られなかった。
【００３７】
ゲルツェル変換の代案は、短時間ゲルツェル変換以外の同じ濾波特性を有するノッチフィルターを設計することである。しかしゲルツェル変換は、その出力がゆっくりと変化するので、高周波の減衰、スパイク検出アルゴニズム…をゲルツェルドメインにおいて実行することが容易である長所を持つ（図７及び図１１を参照）。
【００３８】
非線形フィルター（この事例では図２の振幅表現に対する中央値フィルター）と組み合
った短時間ゲルツェル変換とノッチフィルターとの間の相違が図１２の入力信号に示される。入力信号は、与えれた周波数
s(x)＝1.5+UnitStep(x-128)+0.4cos(αx)+0.2sin(αx)
の信号で擾乱される単一の段階より構成される。
【００３９】
図１３において、振幅表現（α、β）が、この振幅表現の中央値フィルター（α’、β’）バージョンとともに示される。我々の非線形短時間ゲルツェル変換とノッチフィルターとの間の差が図１４において与えられる。我々のフィルターではほとんど完全にステップ関数を再構成できることが明らかに見られる。丸め付け誤差による小さい変動は、ノッチフィルターにより導かれた変動よりかなり小さい。
【図面の簡単な説明】
【００４０】
【図１】乗法問題を加法問題に変化するための流れ図を示す。
【図２】２次元の中央値フィルターによる本発明の実行を示す。
【図３】２個の分離可能な中央値フィルターによる本発明の実行を示す。
【図４】前処理段階として高域通過フィルターによる本発明の実行を示す。
【図５】コンピューター作成の入力信号についての本発明の結果を示す。本発明の結果は、フーリエドメインにおける３個のスペクトル成分が抑制される同様な技法の結果と比較される。
【図６】抑制を必要とする周期３．０１１を有する信号に対して短時間ゲルツェル変換を計算するために使用されるウインドウを示す。
【図７】図５の入力信号の短時間ゲルツェル変換の結果を示す。
【図８】コンピューター作成入力信号についての本発明の結果を示す。本発明の結果は、フーリエドメインにおける３個のスペクトル成分が抑制される同様な技法の結果と比較される。この場合、擾乱信号の周期はフーリエドメインにおけるスペクトル成分の中央周波数とは一致しない。
【図９】コンピューター作成の入力信号についての本発明の結果を示す。本発明の結果は、フーリエドメインにおける７個のスペクトル成分が抑制される同様な技法の結果と比較される。この場合、擾乱信号の周期はフーリエドメインにおけるスペクトル成分の中央周波数とは一致しない。
【図１０】コンピューター作成の入力信号についての本発明の結果を示す。本発明の結果は、フーリエドメインにおける７個のスペクトル成分が抑制される同様な技法の結果と比較される。この場合、擾乱信号の周期は２．９１１から３．１１１の間で変動する。
【図１１】図１０の入力信号の短時間ゲルツェル変換及び前記入力信号のフーリエ変換の結果及び前記入力信号のフーリエ変換の結果を示す。短時間ゲルツェル変換の結果はゆっくりと変化する。
【図１２】与えられた周波数のサイン及びコサイン関数の線形組合せにより分布された単一段階関数を示す。
【図１３】非線形後処理法、中央値フィルターと組み合わせられた本発明の方法の拡大表現を示す。
【図１４】図１２の信号に対するノッチフィルターの結果及び図１３に与えられた非線形後処理段階を使用した本発明の方法による補正の結果を示す。
【図１】
2006289095.tif 000007
【図２】
2006289095.tif 000008
【図３】
2006289095.tif 000009
【図４】
2006289095.tif 000010
【図５】
2006289095.tif 000011
【図６】
2006289095.tif 000012
【図７】
2006289095.tif 000013
【図８】
2006289095.tif 000014
【図９】
2006289095.tif 000015
【図１０】
2006289095.tif 000016
【図１１】
2006289095.tif 000017
【図１２】
2006289095.tif 000018
【図１３】
2006289095.tif 000019
【図１４】
2006289095.tif 000020