画像処理装置、画像処理方法並びに画像処理プログラム
【課題】グリッド像などの周期的パターンを含む放射線画像から周期的パターンを抑制するための通過帯域特性を適切に決定する。
【解決手段】周期的パターンを含む画像を表す画像信号から少なくとも１つの方向について周波数解析を行うことにより、周期的パターンに対応する周波数スペクトルのピーク位置ｐｆとピーク幅ｐｗを検出し、検出されたピーク位置ｐｆおよびピーク幅ｐｗに基づいて、画像信号から周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性Ｆを決定する。
【解決手段】周期的パターンを含む画像を表す画像信号から少なくとも１つの方向について周波数解析を行うことにより、周期的パターンに対応する周波数スペクトルのピーク位置ｐｆとピーク幅ｐｗを検出し、検出されたピーク位置ｐｆおよびピーク幅ｐｗに基づいて、画像信号から周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性Ｆを決定する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、画像信号における周期的パターンに応じた空間周波数成分を抑制する画像処理装置、方法およびプログラム、特に放射線画像において、撮影に使用されたグリッドに起因する周期的パターンを抑制する画像処理装置、方法およびプログラムに関するものである。
【背景技術】
【０００２】
従来より、放射線（Ｘ線、α線、β線、γ線、電子線、紫外線等）を照射するとこの放射線エネルギーの一部が蓄積され、その後、可視光やレーザ光などの励起光を照射すると、蓄積された放射線エネルギーに応じて輝尽発光光が発光される蓄積性蛍光体（輝尽性蛍光体）を利用して、例えば支持体上にこの蓄積性蛍光体を積層した蓄積性蛍光体シートに人体などの被写体を透過した放射線を照射することにより、放射線画像情報を一旦蓄積記録し、この蓄積性蛍光体シートにレーザ光などの励起光を照射して、輝尽発光光を生じさせ、この輝尽発光光を光電変換して画像信号を得る放射線画像読取装置がCR（Computed Radiography）として、広く実用に供されている。
【０００３】
ここで、上述した蓄積性蛍光体シートなどに被写体の放射線画像を撮影記録する際に、被写体により散乱された放射線がシートに照射されないように例えば４本／mm程度の細かなピッチで放射線の透過しない鉛等と透過しやすいアルミニウムや木材等とが交互に配置されたグリッドを被写体とシートとの間に配置して撮影を行うことがある。グリッドを用いて撮影を行うと被写体により散乱された放射線がシートに照射されにくくなるため、被写体の放射線画像のコントラストを向上させることができるが、このグリッド像が画像に含まれるため、画質を低下させるという問題がある。また、グリッド像が含まれた画像を拡大縮小すると、拡大縮小率に応じて折り返しによるエリアジングが生じる。さらに、エリアジングがグリッド像等の空間周波数と重なると、細かな縞模様（モアレ）が生じ、再生画像が見難いものとなってしまう。
【０００４】
特許文献１には、このような周期的パターンを画像から除去する方法として、グリッド像の縞模様の並び方向に一次元ハイパスフィルタ処理を施し、更にグリッド像の縞模様と平行な方向に一次元ローパスフィルタ処理を施してグリッド像に応じた空間周波数成分を原画像から抽出し、抽出された空間周波数成分を原画像から差し引くグリッド像除去方法が開示されている。
【０００５】
また、特許文献２のように、画像データを２次元または１次元フーリエ変換することにより周波数解析を行ってモアレ縞のピーク位置（周波数）、ピークの高さ、半値全幅（ＨＷＦＭ）、ピークの総エネルギー、グリッドの向きなどを取得し、画像のピクセルサイズ、ピーク位置、ピークのエネルギーに基づいて決定したカーネルサイズにより平滑化フィルタリングを行うことによりグリッドによるアーチファクトを抑制した画像を得る方法が提案されている。
【０００６】
また、特許文献３のように、画像データに２次元フーリエ変換を施し、周期的パターンに対応する空間周波数成分を検出し、検出された空間周波数成分のみを除去するフィルタによりフィルタリング処理を行う画像処理方法が開示されている。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開２００３−１５０９５４号公報
【特許文献２】米国特許６２６９１７６号公報
【特許文献３】特開２００９−１９５５１２号公報
【発明の概要】
【発明が解決しようとする課題】
【０００８】
ここで、画像診断に適した高品質な画像を得るためには、被写体を表す成分をできるだけ除去しないで、モアレ等の周期的パターンを除去したいという要求がある。しかし、特許文献１に記載された方法は、グリッド像の縞模様のように一次元ハイパスフィルタ処理を施すことにより周期的パターンに対応する空間周波数成分を検出するものであるため、周期的パターンが高周波数帯域ではなく低周波数帯域に存在した場合には意図したとおりに周期的パターンを検出することができず、特許文献１に記載された方法を適用することができなかった。
【０００９】
また、特許文献２に記載された方法は、平滑化フィルタリング処理が施されることにより局所的なノイズ成分が抑制されるものの、検出された周期的パターンに対応する空間周波数成分のみを画像から除去するものではないため、周期的パターンを含む空間周波数成分を集中的に除去することはできなかった。
【００１０】
また、特許文献３に記載された方法によれば、２次元フーリエ変換により周期的パターンに対応する空間周波数成分を検出するため計算負荷が高く、より少ない計算負荷で周期的パターンに対応する空間周波数成分の通過帯域特性を決定する方法が求められていた。
【００１１】
本発明は上記の事情に鑑みてなされたものであり、周期的パターンに対応する周波数成分が中周波数帯域または低周波数帯域に存在する場合にも、少ない計算負荷で適切に通過帯域特性を決定することができる画像処理装置および画像処理方法ならびに画像処理プログラムを提供することを目的とする。
【課題を解決するための手段】
【００１２】
本発明による画像処理装置は、周期的パターンを含む画像を表す画像信号から少なくとも１つの方向について周波数解析を行うことにより、前記周期的パターンに対応する周波数スペクトルのピーク位置とピーク幅を検出する周期的パターン検出手段と、前記検出されたピーク位置およびピーク幅に基づいて、前記画像信号から前記周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定する通過帯域特性決定手段とを備えたことを特徴とするものである。
【００１３】
本発明による画像処理方法は、周期的パターンを含む画像を表す画像信号から少なくとも１つの方向について周波数解析を行うことにより、前記周期的パターンに対応する周波数スペクトルのピーク位置とピーク幅を検出し、前記検出されたピーク位置およびピーク幅に基づいて、前記画像信号から前記周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定することを特徴とするものである。
【００１４】
本発明による画像処理プログラムは、コンピュータを、周期的パターンを含む画像を表す画像信号から少なくとも１つの方向について周波数解析を行うことにより、前記周期的パターンに対応する周波数スペクトルのピーク位置とピーク幅を検出する周期的パターン検出手段と、前記検出されたピーク位置およびピーク幅に基づいて、前記画像信号から前記周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定する通過帯域特性決定手段として機能させることを特徴とするものである。
【００１５】
上記「周期的パターン」とは、原画像が含む周期的なパターンを持ったノイズを意味する。例えば、上記グリッドを利用して放射線画像を蓄積性蛍光体シートに撮影したときに原画像に含まれるグリッド像やモアレ等を意味する。
【００１６】
また、「通過帯域特性」は、所定の帯域の通過または除去を行うフィルタの特性を意味し、帯域通過特性（バンドパス特性）と帯域除去特性（バンドストップ特性）の両方を含む。
【００１７】
また、ピーク位置は、周波数スペクトルのピーク周波数の代表的な値を表すものであれば任意に定義してよい。同様に、ピーク幅も、周波数スペクトルのピークの幅の代表的な値を表すものであれば任意に定義してよい。
【００１８】
上記通過帯域特性決定手段は、周期的パターンに対応する周波数成分のピーク位置を含み、ピーク位置を通過する帯域に含み、ピークの位置とピーク幅に基づいて決定可能なものであればいかなる通過帯域特性を決定してもよい。
【００１９】
なお、周期的パターンに対応する周波数成分を除去するために、帯域通過の幅はピーク幅以上であることが好ましい。例えば、前記通過帯域特性決定手段は、前記ピーク位置を含み、前記ピーク幅が大きくなるほど通過帯域幅が大きくなるように前記通過帯域特性を特定するものであることが好ましい。
【００２０】
また、通過帯域特性において、低周波数帯域は被写体成分を多く含むため、低周波数帯域においてはフィルタの帯域通過が広いほど周期的パターンに対応する周波数成分と一緒に被写体成分が除去される可能性が高い。このことにより、低周波帯域においては、周期的パターンのピークを含み、かつ幅が狭い通過帯域特性が好ましい。例えば、通過帯域特性決定手段は、前記ピーク位置が低周波数帯域に位置するほど通過帯域の幅が狭くなるように前記通過帯域特性を決定するものであってもよい。
【００２１】
また、通過帯域特性決定手段は、前記ピーク位置と前記ピーク幅により規定される正規分布関数により前記通過帯域特性を特定するものであってもよい。例えば、正規分布関数を、前記ピーク位置を平均とし、前記ピーク幅の半値を標準偏差として規定することができ、かかる正規分布関数のピークを周波数レスポンスの１に対応付けるように正規分布曲線を拡大縮小したものを通過帯域特性として決定してもよい。
【００２２】
また、前記通過帯域特性決定手段は、前記通過帯域特性を、前記ピーク幅以上の幅を有し、前記ピーク位置を中心とする矩形関数により決定してもよい。
【００２３】
また、前記通過帯域特性決定手段は、前記ピークの高さにも基づいて前記通過帯域特性を決定するものであってもよい。スペクトルの大きさが小さい場合などには、周期的パターンに対応する周波数成分を含む周波数帯域を全て抑制しなくとも（周期的パターンに対応する周波数成分を一部画像信号に残留させても）、診断画像として十分使用可能な品質まで周期的パターンを抑制することができるためである。例えば、ピークの高さが大きくなる程帯域通過フィルタによる帯域通過の周波数レスポンスが大きくなり、ピークの高さが小さくなる程帯域通過フィルタによる帯域通過の周波数レスポンスが小さくなるように通過帯域特性を決定することが考えられる。
【００２４】
また、本画像処理装置において、前記画像信号から前記決定された通過帯域特性に基づく１次元フィルタリング処理を行うことにより前記周期的パターンを抑制する抑制手段をさらに備えることが好ましい。
【００２５】
また、本画像処理装置における前記通過帯域特性決定手段は、前記通過帯域特性における通過帯域の幅に上限または下限を設けたものであってもよい。非常に狭い帯域のみを通過（または除去）する帯域通過（または帯域除去）フィルタによりフィルタリング処理を行うことは、画像のアーチファクト等の発生原因となるため好ましくなく、一方、あまりに広帯域を通過する帯域通過（帯域除去）フィルタによりフィルタリング処理を行うことは、被写体成分等の画像診断上重要な情報までも周期的パターンに対応する成分と一緒に除去してしまう可能性があるため好ましくないためである。
【００２６】
前記通過帯域特性決定手段は、前記画像の撮影条件に応じて前記通過帯域特性を異ならせて決定するものであってもよい。なお、画像の撮影部位や撮影条件は、原画像の撮影装置において指定される上記撮影条件や撮影部位の情報などの撮影メニューを利用して取得してもよく、ユーザのマニュアル操作により入力された被写体の臓器や病変などの撮影部位情報を利用してもよい。また、本発明による画像処理装置は、前記画像から該画像の撮影部位を抽出する部位抽出手段をさらに備え、抽出された撮影部位情報を利用してもよい。例えば、本出願人により提案されている、特開２００２−１０９５４８号公報、特開２００３−６６６１号公報に記載された基準となる平均的な心胸郭の輪郭と略相似形のテンプレートを用いてテンプレートマッチング処理を行なって胸郭を自動検出する方法を用いてもよい。
【発明の効果】
【００２７】
本発明の画像処理装置、画像処理方法および画像処理プログラムによれば、周期的パターンを含む画像を表す画像信号から少なくとも１つの方向について周波数解析を行うことにより、周期的パターンに対応する周波数スペクトルのピーク位置とピーク幅を検出し、検出されたピーク位置およびピーク幅に基づいて、画像信号から周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定するため、ピーク位置およびピーク幅に基づいて、画像の被写体成分を過剰に抑制しない範囲で周波数成分を適切に抑制できる１次元フィルタの通過帯域特性を好適に決定できる。この結果、比較的小さい計算負荷で高品質な画像を得ることができる。
【図面の簡単な説明】
【００２８】
【図１】放射線画像撮影装置の概略を示した図
【図２】グリッド撮影により得られた放射線画像を示す図
【図３】放射線画像読取装置の一例を示す斜視図
【図４】走査方向と読取画像との関係を示した図
【図５】本実施形態における画像処理装置の概略構成図
【図６】本実施形態における画像処理装置の処理の流れを示す図
【図７】本実施形態における領域抽出処理を説明する概略図
【図８】本実施形態における規格化処理を説明する図
【図９Ａ】本実施形態における通過帯域特性決定処理を説明する図
【図９Ｂ】本実施形態における通過帯域特性決定処理の変形例を説明する図
【図１０Ａ】本実施形態における抑制処理を説明する図
【図１０Ｂ】本実施形態における抑制処理の変形例を説明する図
【図１１】本実施形態における通過帯域特性の変形例を示す図（その１）
【図１２】本実施形態における通過帯域特性の変形例を示す図（その２）
【図１３】本実施形態における通過帯域特性の変形例を示す図（その３）
【発明を実施するための形態】
【００２９】
以下、図面を参照して本発明の実施の形態を説明する。尚、以下の実施の形態では、本発明における周期的パターン抑制処理装置を放射線画像読取装置に用いた場合について説明するが、デジタルカメラ等を用いた通常撮影において網戸やブラインド等を介して撮影した際に撮影画像に含まれる周期的パターンを抑制するための画像処理装置等に用いられてもよい。尚、放射線画像読取装置は、蓄積性蛍光体シートに記録された人体の放射線画像をレーザビーム走査によりデジタル画像信号として読み取るものである。
【００３０】
図１は放射線画像撮影装置の概略図である。放射線源１から放射された放射線２は、被写体３を透過して静止グリッド（以下、単に「グリッド」と表記する）４に到達する。グリッド４は放射線２を吸収する鉛４ａと、放射線２を透過するアルミニウム４ｂとが例えば４本／ｍｍ程度のピッチで交互に配置されているものである。また、放射線２がアルミニウム４ｂを透過して蓄積性蛍光体シート１１に入射するように鉛４ａは位置に応じて多少傾きを変化させて設置されている。
【００３１】
従って、被写体３を透過した放射線２は、鉛４ａに吸収されて蓄積性蛍光体シート１１への到達を遮られる一方で、アルミニウム４ｂを透過して蓄積性蛍光体シート１１に到達する。蓄積性蛍光体シート１１には被写体像と共に４本／ｍｍの縞模様状のグリッド像が記録される。一方、被写体３内で散乱された散乱放射線２ａは、位置に応じて傾きを持って設置された鉛４ａによって吸収され、又はグリッド４の表面で反射するため、蓄積性蛍光体シート１１には到達しない。従って、蓄積性蛍光体シート１１は散乱放射線２ａの照射の少ない鮮明な放射線画像を記録することができる。図２は、図１に示した放射線画像撮影装置を用いて撮影した際に蓄積性蛍光体シート１１に蓄積記録された被写体像５とグリッド像に起因するモアレパターン６の放射線画像の一例である。
【００３２】
図３は、放射線画像読取装置の斜視図と機能ブロック図を組み合わせた図である。読取部１０の所定位置にセットされた蓄積性蛍光体シート１１は、駆動手段（不図示）によって駆動されるエンドレスベルト等のシート搬送手段１５によって、例えば走査ピッチ１０本／ｍｍで矢印Ｙ方向に搬送（副走査）される。一方、レーザ光源１６から発せられた光ビーム１７は、モータ２４により駆動され矢印方向に高速回転する回転多面鏡１８によってｆθレンズ等の集束レンズ１９に向けて反射される。光ビーム１７は集束レンズ１９を通過した後、ミラー２０によって蓄積性蛍光体シート１１に向けて反射される。蓄積性蛍光体シート１１は光ビーム１７によって副走査方向（矢印Ｙ方向）と略直角な矢印Ｘ方向に主走査される。
【００３３】
蓄積性蛍光体シート１１において光ビーム１７が照射された箇所からは蓄積記録されている放射線画像情報に応じた光量の輝尽発光光２１が発散される。輝尽発光光２１は光ガイド２２の入射端面２２ａから入射し、光ガイド２２内の内部を全反射を繰り返して出射端面２２ｂから出射する。出射された輝尽発光光２２はフォトマルチプライヤ２３によって受光されて光電変換され、アナログ画像信号Ｓａに変換される。
【００３４】
アナログ画像信号Ｓａはログアンプ２６によって対数的に増幅された後、Ａ／Ｄ変換器２８において空間周波数ｆｓ＝１０ｃｙｃｌｅ／ｍｍに対応するサンプリング間隔でサンプリングされてデジタル化され、デジタル画像信号Ｓｄ（以下、単に「画像信号Ｓｄ」と表記する）を出力する。この画像信号Ｓｄは、図４に示すように、蓄積性蛍光体シート１１に対して主走査方向（横方向）に光ビーム１７を走査させながら蓄積性蛍光体シート１１を副走査方向（縦方向）に移動させることにより蓄積性蛍光体シート１１を２次元走査して得られた放射線画像情報を示す。尚、このようにして得られた画像信号Ｓｄには、被写体３に対応する放射線画像５の情報の他に、グリッド４に対応するグリッド像（不図示）に起因するモアレパターン６の情報も含まれている。
【００３５】
画像信号Ｓｄは一旦記憶部２９に記憶された後、画像信号処理部３０に入力される。画像信号処理部３０は本発明における画像処理方法を実行するための画像処理装置４０を備えている。
【００３６】
図５は本実施の形態の画像処理装置４０の概略構成を表す図である。図５を用いて、本実施形態による画像処理装置４０について以下説明する。
【００３７】
本実施形態の画像処理装置４０は、周期的パターンを含む画像を表す画像信号Ｓｄから少なくとも１つの方向について周波数解析を行うことにより、周期的パターンに対応する周波数スペクトルのピーク位置とピーク幅を検出する周期的パターン検出手段４１と、検出されたピーク位置およびピーク幅に基づいて、画像信号から周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定する通過帯域特性決定手段４２と、画像信号から決定された通過帯域特性に基づく１次元フィルタリング処理を行うことにより周期的パターンを抑制する抑制手段４３とを備えるものである。
【００３８】
本実施形態における画像処理プログラムと画像処理プログラムが参照するデータは、インストール時に記憶部２９に記憶され、起動時に記憶部２９に含まれるメモリにロードされる。画像処理プログラムは、画像処理装置４０を構成する画像信号処理部３０の中央処理装置に実行させる処理として、周期的パターン検出処理と、通過帯域特性決定処理と、抑制処理を規定している。そして、プログラムの規定にしたがって、中央処理装置が上記各処理を実行することにより、画像信号処理部３０の中央処理装置は、周期的パターン検出手段４１、通過帯域特性決定手段４２、抑制手段４３として機能する。
【００３９】
図６は、本画像処理装置４０の画像処理の流れを表すフローチャートである。図７は周期的パターン検出手段４１によるライン状のサンプル領域抽出処理を説明するための図である。図６および図７に従って、本実施形態の画像処理装置４０の画像処理について詳細に説明する。
【００４０】
まず、図６に示すように、画像データＳｄが周期的パターン検出手段４１に入力される（Ｓ０１）。すると、図７に示すように、周期的パターン検出手段４１は、画像データＳｄにより表される画像におけるｘ方向およびｙ方向のそれぞれについて、ｎｘｍ本のライン状の領域を抽出する。なお、ｘ方向とｙ方向とで領域の抽出の仕方は同一であるため、ここでは、ｘ方向の抽出についてのみ説明する。
【００４１】
周期的パターン検出手段４１は、まず画像データＳｄにより表される画像Ｇ１上にｘ方向に延在するｎ（本実施形態ではｎ＝９）の矩形領域Ａ１〜Ａｎを抽出する（Ｓ０２）。そして、ｎ個の領域画像データＳＲｘｉのそれぞれに、ライン状のｍ個（本実施形態ではｍ＝９）の領域画像データＳＲｘｉｊ（ｉ＝１〜ｎ,ｊ＝１〜ｍ）を得る。ここで、領域Ａ１のｙ方向の中心は画像上縁から画像高さの１／３０の位置に配置され、領域Ａ９のＹ方向の中心は画像上縁から画像高さの１／３０の位置に配置され、領域Ａ２乃至Ａ７は領域Ａ１およびＡ９間に等間隔に位置するように配置される。また、ライン状のｍ個（本実施形態ではｍ＝９）の各領域画像データＳＲｘｉｊ（ｉ＝１〜ｎ,ｊ＝１〜ｍ）は、ｙ方向に３画素間隔で配置され、ｘ方向への画素数は１０２４である。
【００４２】
周期的パターン検出手段４１は、ＳＲｘｉｊ（ｉ＝１〜ｎ,ｊ＝１〜ｍ）のそれぞれに高速フーリエ変換を施して周波数解析を行う（Ｓ０３）。具体的には各領域画像データＳＲｘｉｊの画素配列ＳＲｘｉｊ（ｋ）（ｋ＝０〜１０２３）に対して高速フーリエ変換が施される。すなわち、スペクトルをＰｉｊ（ｋ）（ｋ＝０〜１０２３）、高速フーリエ変換により得られる演算結果の実部をＲｅ（ｋ）、虚部をＩｍ（ｋ）とすると、下記の式（１）によりスペクトルＰｉｊ（ｋ）が算出される。
【００４３】
【数１】
なお、スペクトルは実際にはナイキスト周波数以下の周波数についてのみ算出されるため、データ数としては１０２４の１／２すなわち５１２となる。したがって、以降の説明においては、スペクトルの参照符号としてＰｉｊ（ｋ′）（ｋ′＝０〜５１１）を用いるものとする。そして、周期的パターン検出手段４１は、ｎｘｍ個の領域画像データＳＲｘｉｊについてそれぞれ算出された、ｍｘｎ個のスペクトルＰｉｊ（ｋ′）（ｉ＝１〜ｎ,ｊ＝１〜ｍ）を加重平均処理したスペクトルＰａ（ｋ’）を算出する（Ｓ０４）。
【００４４】
図８左図はスペクトルＰａ（ｋ’）の例を示す図である。図８左図に示すように求められたスペクトルは低周波数から高周波数となるにつれて小さくなるとともに、所々に細かなピークを有するものとなっている。
【００４５】
次いで、周期的パターン検出手段４１は、加重平均処理したスペクトルＰａ（ｋ’）を規格化する（Ｓ０５）。周期的パターン検出手段４１は、スペクトルＰａ（ｋ’）を平滑化してスペクトルＰａ（ｋ’）の細かなピークが平滑化により除去されたスペクトルＰｂ（ｋ’）を算出する。そして算出されたスペクトルＰｂ（ｋ’）を図８左図に示す画像信号Ｓｄから算出されたスペクトルＰａ（ｋ’）から減算することにより規格化し、図８右図の規格化されたスペクトルＰ（ｋ’）を取得する。規格化処理により、ピークをより明確にして精密に評価するためである。
【００４６】
次いで、周期的パターン検出手段４１は、規格化されたスペクトルＰからスペクトルの最大値における周波数である最大値周波数（ピーク位置）とピーク幅を取得する（Ｓ０６）。図９Ａは、通過帯域特性を決定する処理を説明する図であり、図８右の規格化されたスペクトルＰを拡大したものである。まず、周期的パターン検出手段４１は、図８右および図９Ａに示すように、スペクトルＰについて、誤検出を避けるため所定の帯域（ｆＴｈ１≦ｆ≦ｆＴｈ２）を評価帯域とし、最大スペクトルを示した周波数（ピーク位置）ｐｆとピークの最大値（ピーク高さｐｈ）を取得する。そして、評価帯域においてピーク位置を含み所定のしきい値Ｔｈ以上の周波数帯域の幅をピーク幅ｐｗとして取得し、取得されたピーク位置ｐｆ、ピーク幅ｐｗ、ピーク高さｐｈをメモリに保存する。なお、本実施形態では、図８右図に示すように、ｆＴｈ１＝０．３（ｃｙｃｌｅ／ｍｍ）、ｆＴｈ２＝３．３（ｃｙｃｌｅ／ｍｍ）である。かかる評価帯域は、被写体成分を多く含む低周波数を除き、グリッド像の折り返し成分であるモアレパターン６に対応する周波数成分が存在すると予想される、ナイキスト周波数以下の周波数帯域を含むように設定されている。
【００４７】
次いで、通過帯域特性決定手段４２は、取得したピーク位置ｐｆおよびピーク幅ｐｗに基づいて通過帯域特性を決定する（Ｓ０７）。本実施形態においては、平均μをｐｆ、標準偏差σをｐｗ／２とすることにより規定される、下記の式（２）の正規分布関数Ｎに基づいて通過帯域特性（フィルタ特性）を特定する。すなわち、図９Ａに示されるように、スペクトルＰにおける周期的パターン成分に対応するピーク位置ｐｆで周波数レスポンスが１となるように正規分布関数Ｎを拡大縮小した関数が通過帯域特性Ｆとして決定される。本実施形態では、通過帯域特性決定手段４２は同様にｙ方向についても通過帯域特性を決定する。
【数２】
【００４８】
次いで、抑制手段４３は、決定されたｘ、ｙ両方向の通過帯域特性に基づいて、それぞれ周知の手法を用いて１次元の帯域通過フィルタ４４を作成する（Ｓ０８）。例えば、フィルタ特性に基づくディジタルフィルタの設計については、高橋・池原、「ディジタルフィルタ」、培風館（特に第6.3節、pp84-86）、岩田・実践、「ディジタル・フィルタ設計入門」、ＣＱ出版（特に第8章）などに記載された方法が適用できる。
【００４９】
図１０Ａは、本実施形態の抑制手段４３による周期的パターンの抑制処理を説明する図である。本実施形態による抑制手段４３は、図９Ａに示すように、決定した通過帯域特性Ｆに基づいて、決定された通過帯域を除去する通過帯域特性Ｆ’を有する帯域除去フィルタ（バンドストップフィルタ）ＢＳＦを１次元の帯域通過フィルタ４４として生成する。
【００５０】
そして、抑制手段４３は、画像信号から、決定された通過帯域特性に基づく１次元フィルタリング処理を行うことにより周期的パターンを抑制する（Ｓ０９）。すなわち、本実施形態では、ｘ、ｙ方向ごとに作成された１次元フィルタ（ＢＳＰ）によるｘ、ｙ方向に対するフィルタリング処理を施す。
【００５１】
ここでは、図１０Ａに示すように、抑制手段４３は、画像信号Ｓｄから生成した帯域除去フィルタＢＳＦによるフィルタリング処理を行い、周期的パターンを抑制した画像である、第２処理済信号Ｓｐを抽出する。そして、この第２処理済信号Ｓｐに対応する画像が生成されて（Ｓ１０）、生成された画像が出力される（Ｓ１１）。その後、出力された画像は、モニタ等の表示手段（不図示）に表示され、又はプリンタ等の出力手段（不図示）に出力される。
【００５２】
本実施形態によれば、検出されたピーク位置ｐｆおよびピーク幅ｐｗに基づいて、画像信号から周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定するため、ピーク位置ｐｆおよびピーク幅ｐｗに基づいて、画像の被写体成分を過剰に抑制しない範囲で周波数成分を適切に抑制する１次元フィルタの通過帯域特性を好適に決定できる。このため、本実施形態による画像処理を低周波数帯域に存在する周期的パターンにも好適に適用できる。また、本実施形態では、周期的パターンに対応する空間周波数成分の検出処理および周期的パターンに対応する空間周波数成分の抑制処理を１次元の処理を行って、周波数成分を適切に抑制する１次元フィルタの通過帯域特性を決定する。このため、比較的小さい計算負荷にも関わらず、決定された通過帯域特性に基づいて好適に周期的パターンの抑制処理を行うことができる。また、集中的に周期的パターンに対応する空間周波数のみを抑制できるため高品質な画像を得ることができる。
【００５３】
また、本実施形態による通過帯域特性決定手段４２が、ピーク位置ｐｆを平均とし、ピーク幅ｐｗの半値を標準偏差とする正規分布関数に基づいて、通過帯域特性Ｆを特定するものであるため、決定された通過帯域特性Ｆはピーク位置ｐｆを中心としピーク幅ｐｗを適切に反映して滑らかに分布する特性を示す。すなわち、本実施形態の通過帯域特性Ｆは、ピーク位置を含み、ピーク幅が大きくなるほど通過帯域幅が大きくなるものであるため、周期的パターンに対応する空間周波数成分のピーク位置ｐｆを包含するような通過帯域特性Ｆを適切に決定することができる。さらに、本実施形態によれば、ピーク位置を中心にピーク幅ｐｗを含むような通過帯域特性Ｆが決定される可能性が高く好適に周期的パターンに対応する空間周波数成分を含む通過帯域特性Ｆを決定できる。
【００５４】
また、本実施形態による画像処理装置４０は、画像信号から決定された通過帯域特性に基づく１次元フィルタリング処理を行うことにより周期的パターンを抑制する抑制手段４３をさらに備えているため、決定された通過帯域特性Ｆに基づいて、好適に周期的パターンに対応する空間周波数成分を抑制できる。
【００５５】
以下、本実施形態の変形例について記す。
【００５６】
上記実施形態において周期的パターン検出手段４１は、ピーク位置、ピーク幅、ピーク高さについて、ピークの位置、幅、高さのそれぞれの代表的な値を定義するものであれば、任意の定義により特定してよい。例えば、ピーク幅を所定値以上の周波数レスポンスを有する帯域の数倍（２乃至３倍）に設定してもよい。
【００５７】
また、通過帯域特性決定手段４２は、周期的パターンに対応する周波数成分のピーク位置を含み、ピークの位置とピーク幅に基づいて決定可能な通過帯域特性を決定するものであればいかなる通過帯域特性を決定してもよい。図９Ｂは、通過帯域特性を矩形関数で規定する例を示す図である。例えば、図９Ｂに示すように、ピーク幅ｐｗと同じ周波数帯域を通過させ、ピーク位置ｐｆを中心とする矩形関数により通過帯域特性Ｆ１を決定してもよい。
【００５８】
また、本発明により決定された通過帯域特性は、ピーク位置ｐｆとピーク幅ｐｗに基づいて決定された通過帯域特性に種々の補正を行ったものであってもよい。以下に、本発明による通過帯域特性の変形例を説明する。
【００５９】
一般に、低周波数帯域は被写体成分を多く含むため、低周波数帯域においてはフィルタの帯域通過が広いほど周期的パターンに対応する周波数成分と一緒に被写体成分が除去される可能性が高い。このことにより、低周波帯域においては、周期的パターンのピークを含み、かつ幅が狭い通過帯域特性が好ましい。図１１は、通過帯域特性の変形例を示す図である。本実施形態の変形例として、図１１に示すように、上記実施形態における通過帯域特性決定手段４２は、上記通過帯域特性Ｆが周期的パターンのピーク位置ｐｆを含み、ピーク位置ｐｆが低周波数帯域に位置するほど通過帯域の幅が狭くなるように決定してもよい。例えば、ピーク位置ｐｆに通過帯域特性Ｆの通過帯域幅を対応付けて記憶部２９に記憶しておき、周期的パターン検出手段４１がピーク位置ｐｆを検出すると、検出されたピーク位置ｐｆ対応付けられた通過帯域幅になるように上記実施形態の通過帯域特性Ｆを幅方向に拡大縮小する補正を行い、その結果得られた関数を通過帯域特性Ｆ２として決定することが考えられる。この場合には、従来技術に比べて低周波成分に多く含まれる被写体成分を維持しつつモアレ等の周期的パターンに対応する空間周波数成分のみを好適に除去できる。
【００６０】
また、非常に尖鋭な帯域のみを通過（または除去）する帯域通過（帯域除去）フィルタによりフィルタリング処理を行うことは、画像のアーチファクト等の発生原因となるため好ましくない。そこで、通過帯域特性決定手段４２は、通過帯域の幅に下限を設けてもよい。図１２は、上記実施形態の通過帯域特性の変形例を示す図である。例えば、周期的パターンがあまりに鋭いピークである場合など、上記通過帯域特性Ｆの幅が所定値より小さい場合に、通過帯域の幅が広くなるように上記実施形態の通過帯域特性Ｆを補正した関数を通過帯域特性として決定してもよい。例えば、図１２に示すように、ピーク位置ｐｆおよびピーク幅ｐｗに基づいて決定された通過帯域特性Ｆ３ａの所定の周波数レスポンスＲＴｈにおける通過帯域幅ｆｗが所定値より小さい場合、通過帯域特性Ｆの幅を所定の割合で広くなるように補正したものを通過帯域特性Ｆ３として決定してもよい。または、ピーク幅ｐｗを所定の割合で大きくなるように補正し、補正したピーク幅に基づいて通過帯域特性を決定するものであってもよい。
【００６１】
上記場合には、通過帯域が狭い通過帯域特性に基づいてフィルタリング処理をおこなうことによるアーチファクトの発生を抑制でき、より品質の高い画像を生成することができる。
【００６２】
また、あまりに広帯域を通過する帯域通過（帯域除去）フィルタによりフィルタリング処理を行うことは、被写体成分等の画像診断上重要な情報までも周期的パターンに対応する成分と一緒に除去してしまう可能性があるため好ましくない。そこで、通過帯域特性決定手段４２は、通過帯域の幅に上限を設けてもよい。例えば、通過帯域特性Ｆの幅ｆｗが所定値より大きい場合に、通過帯域特性の幅が狭くなるように通過帯域特性Ｆを補正したものを通過帯域特性として決定してもよい。例えば、図１２の通過帯域特性Ｆ４に示すように、ピーク位置ｐｆおよびピーク幅ｐｗに基づいて決定された通過帯域特性Ｆ４ａの所定の周波数レスポンスＲｔｈにおける通過帯域幅ｆｗが所定値より大きい場合、通過帯域特性Ｆ４ａの幅を所定の割合で狭くなるように補正したものを通過帯域特性Ｆ４として決定してもよい。または、ピーク幅ｐｗが所定値よりも大きい場合に、ピーク幅ｐｗを所定の割合で小さくした補正ピーク幅ｐｗにより通過帯域特性を決定してもよい。また、ピーク幅ｐｗを用いて算出された通過帯域特性Ｆ４ａに対して、ピーク位置ｐｆから所定値ｆａだけ低い周波数より小さい周波数帯域における周波数レスポンスを０にするように通過帯域特性Ｆ４ａを補正したものを通過帯域特性Ｆ４ｂとして決定してもよい。すなわち、図１２に示すように、周波数（ｐｆ−ｆａ）より小さい周波数帯域における周波数レスポンスを０、周波数（ｐｆ−ｆａ）以上の周波数レスポンスをＦ４ａで規定するような通過帯域特性Ｆ４ｂにしてもよい。また、同様に低周波数側だけでなく高周波数側に対しても通過帯域を狭くする補正を行ってもよい。例えば、上記通過帯域特性Ｆ４ａまたはＦ４ｂをピーク位置から所定値ｆｂだけ高い周波数（ｐｆ＋ｆｂ）より大きい周波数帯域における周波数レスポンスを０とし、周波数（ｐｆ＋ｆｂ）以下の周波数レスポンスを通過帯域特性Ｆ４ａにより規定する通過帯域特性Ｆ４ｃに補正してもよい。
【００６３】
上記場合には、通過帯域が広い通過帯域特性に基づいてフィルタリング処理をおこなうことによる被写体成分等の画像診断上重要な情報までも周期的パターンに対応する成分と一緒に除去してしまう可能性を低減でき、より品質の高い画像を生成することができる。
【００６４】
また、通過帯域特性決定手段４２は、ピークの高さ（スペクトルの大きさ）にも基づいて通過帯域特性を決定するものであってもよい。図１３は、本実施形態の変形例を示す図であり、ピークの高さｐｈ（スペクトルの大きさ）が小さいほど通過帯域特性の周波数レスポンスが小さくなるようにピーク位置ｐｆとピーク幅ｐｗに基づいて決定された通過帯域特性Ｆ５ａを補正して通過帯域特性Ｆ５を決定した例を説明するための図である。図１３の通過帯域特性Ｆ５に示すように、スペクトルの高さｐｈが小さい場合などには、周期的パターンに対応する周波数成分を含む周波数帯域を全て抑制しなくとも（周期的パターンに対応する周波数成分を一部画像信号に残留させても）、診断画像として十分使用可能な品質まで周期的パターンを抑制することができるためである。この場合、ピークの高さｐｈが大きくなる程帯域通過フィルタによる帯域通過の周波数レスポンスが大きくなり、ピークの高さｐｈが小さくなる程帯域通過フィルタによる帯域通過の周波数レスポンスが小さくなるものであれば、任意の方法により通過帯域特性を決定してよい。
【００６５】
また、通過帯域特性決定手段４２は、画像の撮影条件に応じて通過帯域特性を異ならせて決定するものであってもよい。例えば、記憶部２９に、撮影部位や撮影条件に適応する通過帯域特性を予め決定して撮影部位や撮影条件ごとに対応付けて記憶しておき、通過帯域特性決定手段４２が、処理対象である画像の撮影部位や撮影条件に対応付けられた通過帯域特性を特定し、特定された通過帯域特性を用いて周期的パターンの抑制処理を行うことが考えられる。
【００６６】
上記場合には、撮影部位や撮影条件に応じた通過帯域特性に応じて適切に抑制処理を施すことができるため、好適に周期的パターンに対応する周波成分を抑制した第２処理済信号を抽出することができる。結果として、より高品質に周期的パターンに対応する周波成分を抑制した画像を生成することができる。
【００６７】
なお、画像の撮影部位や撮影条件は、原画像の撮影装置において指定される上記撮影条件や撮影部位の情報などの撮影メニューを利用して取得してもよく、ユーザのマニュアル操作により入力された被写体の臓器や病変などの撮影部位情報を利用してもよい。また、上記画像処理装置４０は、不図示の画像の撮影部位を抽出する部位抽出手段をさらに備え、抽出された撮影部位情報を利用してもよい。撮影部位情報の抽出方法として、部位又は病変などが抽出できるものであればいかなる方法を適用してもよい。例えば、本出願人により提案されている、特開２００２−１０９５４８号公報、特開２００３−６６６１号公報に記載された基準となる平均的な心胸郭の輪郭と略相似形のテンプレートを用いてテンプレートマッチング処理を行なって胸郭を自動検出する方法を用いてもよい。
【００６８】
なお、かかる場合には、通過帯域特性決定手段４２は、ユーザによる撮影部位の入力操作を必要とせず、自動抽出された結果に基づいて撮影部位を抽出でき、抽出された撮影部位に基づいて適切に通過帯域特性を決定することができる。このため、ユーザは撮影部位に応じて周期的パターンに対応する周波成分を抑制した画像を容易かつ高品質に生成することができる。また、通過帯域特性決定手段４２が撮影メニューから撮影部位または撮影条件を自動抽出して、自動抽出された結果に基づいて通過帯域特性を決定した場合にも同じ効果が得られる。
【００６９】
また、同一の装置で続けて同じ部位の周期的パターンの抑制処理を行う場合など、ほぼ同じような撮影条件の場合には、最初の被写体の画像信号Ｓｄに対して本実施形態の画像処理により帯域通過フィルタ（帯域除去フィルタ）を決定し、かかるフィルタを後続の画像信号にも適用してもよい。ほぼ同じような撮影条件の場合には、同じような周期的パターンを有する画像が取得されると考えられるため、代表的な画像について決定した帯域通過フィルタ（帯域除去フィルタ）を用いることで十分周期的パターンの抑制効果が得られるとともに画像ごとに通過帯域特性を決定する時間と労力を省略できるため効率がよい。
【００７０】
なお、抑制手段４３は、最終的に画像信号Ｓｄに対して決定した通過帯域特性に基づいて周期的パターンの抑制処理を行った画像信号（本実施形態におけるＳｐ）が得られるものであれば、水平方向および垂直方向に帯域通過フィルタおよび帯域除去フィルタのいかなる組合せおよびフィルタリング処理の順序により抑制処理を行ってもよい。例えば、決定した通過帯域特性Ｆに基づいて、決定した通過帯域を通過させる通過帯域特性Ｆを有する帯域通過フィルタ（バンドパスフィルタ）ＢＰＦを生成してもよい。図１０Ｂは、上記実施形態の抑制処理の変形例を示す図である。図１０Ｂに示すように、画像信号Ｓｄを帯域通過フィルタＢＰＦによるフィルタリング処理した画像信号Ｓｅ（画像信号Ｓｄから抽出した周期的パターン画像）を生成し、画像信号Ｓｄから減算することにより上記実施形態と同様に周期的パターンを抑制できる。この場合、抑制手段４３は帯域通過フィルタＢＰＦと減算器４５を備え、画像信号Ｓｄに対して帯域通過フィルタＢＰＦによるフィルタリング処理を施して、周期的パターン成分を抽出した周期的パターン画像信号Ｓｅを算出し、その後、減算器４５は周期的パターン画像信号Ｓｅを減算して、周期的パターンを抑制した画像である、処理済信号Ｓｐを抽出する。
【００７１】
上記実施形態において周期的パターン検出手段４１は、複数のサンプル領域に変えて１つの代表的な領域についてのみ複数の線上の画像信号を取得して、周波数解析を行って加重平均をとってスペクトルＰａを取得してもよく、１つの線上の画像信号を取得して、これについてのみ周波数解析を行ったスペクトルを規格化されたスペクトルＰａの代わりに用いてもよい。
【００７２】
上記実施形態において周期的パターン検出手段４１は、規格化処理を省略してもよい。ただし、規格化処理を行った場合には、ピーク幅やピーク値やピーク高さなどの情報を精度よく検出でき好ましい。
【００７３】
また、本実施形態において、絶縁基板上に夫々が画素に対応する複数個の光電変換素子を２次元状に形成した放射線画像読取装置であって、この放射線画像読取装置上に形成された画像情報を担持する放射線が照射されると照射された放射線を可視光に変換する蛍光体層（シンチレータ）を積層して成る放射線固体検出器（以下、「光変換方式」の放射線固体検出器という）を用いた放射線画像読取装置に本発明を適用した例について説明したが、本発明は、上記実施形態に限定されず種々のタイプの放射線固体検出器を用いた放射線画像読取装置に適用可能である。例えば、絶縁基板上に夫々が画素に対応する複数個の電荷収集電極を２次元状に形成した放射線画像読取装置であって、この２次元画像読取装置上に形成された、画像情報を担持する放射線が照射されると前記画像情報を担持する電荷を発生する放射線導電体とを積層して成る放射線固体検出器（以下、「直接変換方式」の放射線固体検出器という）を用いた放射線画像撮影装置に適用してもよい。
【００７４】
また、光変換方式の放射線固体検出器としては、例えば特開昭59-211263 号、特開平2-164067号、ＰＣＴ国際公開番号ＷＯ92/06501号、Signal,noise,and read out considerations in the development of amorphous silicon photodiode arrays for radiotherapy and diagnostic x-ray imaging，L.E.Antonuk et.al ，University of Michigan，R.A.Street Xerox，PARC，SPIE Vol.1443 Medical Imaging V;Image Physics(1991) ，p.108-119 等に記載された放射線固体検出器が適用できる。
【００７５】
一方、直接変換方式の放射線固体検出器としては、例えば、(i) 放射線の透過方向の厚さが通常のものより10倍程度厚く設定された放射線固体検出器（MATERIAL PARAMETERS IN THICK HYDROGENATED AMORPHOUS SILICON RADIATION DETECTORS,Lawrence Berkeley Laboratory.University of California,Berkeley.CA 94720 XeroxParc.Palo Alto.CA 94304)、あるいは(ii)放射線の透過方向に、金属板を介して２つ以上積層された放射線固体検出器(Metal/Amorphous Silicon Multilayer Radiation Detectors,IEE TRANSACTIONS ON NUCLEAR SCIENCE.VOL.36.NO.2.APRIL 1989) 、あるいは(iii) ＣｄＴｅ等を使用した放射線固体検出器（特開平1-216290号）等に記載された放射線固体検出器が適用できる。
【００７６】
上記の各実施形態はあくまでも例示であり、上記のすべての説明が本発明の技術的範囲を限定的に解釈するために利用されるべきものではない。この他、上記の実施形態におけるシステム構成、ハードウェア構成、処理フロー、モジュール構成や具体的処理内容等に対して、本発明の趣旨から逸脱しない範囲で様々な改変を行ったものも、本発明の技術的範囲に含まれる。
【符号の説明】
【００７７】
１ 放射線源
２ 放射線
３ 被写体
４ グリッド
５ 被写体像
６ グリッド線
１０ 読取部
２９ 記憶部
３０ 画像処理部
４０ 画像処理装置
４１ 周期的パターン検出手段
４２ 通過帯域特性決定手段
４３ 抑制手段
４４ 帯域通過フィルタ
４５ 減算器
ＢＰＦ 帯域通過フィルタ
ＢＳＦ 帯域除去フィルタ
Ｆ，Ｆ１，Ｆ２，Ｆ３，Ｆ４，Ｆ５，Ｆ６ 通過帯域特性
Ｐ，Ｐａ，Ｐｂ 周波数スペクトル
ｐｆ ピーク位置（周波数）
ｐｗ ピーク幅
ｐｈ ピークの高さ
【技術分野】
【０００１】
本発明は、画像信号における周期的パターンに応じた空間周波数成分を抑制する画像処理装置、方法およびプログラム、特に放射線画像において、撮影に使用されたグリッドに起因する周期的パターンを抑制する画像処理装置、方法およびプログラムに関するものである。
【背景技術】
【０００２】
従来より、放射線（Ｘ線、α線、β線、γ線、電子線、紫外線等）を照射するとこの放射線エネルギーの一部が蓄積され、その後、可視光やレーザ光などの励起光を照射すると、蓄積された放射線エネルギーに応じて輝尽発光光が発光される蓄積性蛍光体（輝尽性蛍光体）を利用して、例えば支持体上にこの蓄積性蛍光体を積層した蓄積性蛍光体シートに人体などの被写体を透過した放射線を照射することにより、放射線画像情報を一旦蓄積記録し、この蓄積性蛍光体シートにレーザ光などの励起光を照射して、輝尽発光光を生じさせ、この輝尽発光光を光電変換して画像信号を得る放射線画像読取装置がCR（Computed Radiography）として、広く実用に供されている。
【０００３】
ここで、上述した蓄積性蛍光体シートなどに被写体の放射線画像を撮影記録する際に、被写体により散乱された放射線がシートに照射されないように例えば４本／mm程度の細かなピッチで放射線の透過しない鉛等と透過しやすいアルミニウムや木材等とが交互に配置されたグリッドを被写体とシートとの間に配置して撮影を行うことがある。グリッドを用いて撮影を行うと被写体により散乱された放射線がシートに照射されにくくなるため、被写体の放射線画像のコントラストを向上させることができるが、このグリッド像が画像に含まれるため、画質を低下させるという問題がある。また、グリッド像が含まれた画像を拡大縮小すると、拡大縮小率に応じて折り返しによるエリアジングが生じる。さらに、エリアジングがグリッド像等の空間周波数と重なると、細かな縞模様（モアレ）が生じ、再生画像が見難いものとなってしまう。
【０００４】
特許文献１には、このような周期的パターンを画像から除去する方法として、グリッド像の縞模様の並び方向に一次元ハイパスフィルタ処理を施し、更にグリッド像の縞模様と平行な方向に一次元ローパスフィルタ処理を施してグリッド像に応じた空間周波数成分を原画像から抽出し、抽出された空間周波数成分を原画像から差し引くグリッド像除去方法が開示されている。
【０００５】
また、特許文献２のように、画像データを２次元または１次元フーリエ変換することにより周波数解析を行ってモアレ縞のピーク位置（周波数）、ピークの高さ、半値全幅（ＨＷＦＭ）、ピークの総エネルギー、グリッドの向きなどを取得し、画像のピクセルサイズ、ピーク位置、ピークのエネルギーに基づいて決定したカーネルサイズにより平滑化フィルタリングを行うことによりグリッドによるアーチファクトを抑制した画像を得る方法が提案されている。
【０００６】
また、特許文献３のように、画像データに２次元フーリエ変換を施し、周期的パターンに対応する空間周波数成分を検出し、検出された空間周波数成分のみを除去するフィルタによりフィルタリング処理を行う画像処理方法が開示されている。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開２００３−１５０９５４号公報
【特許文献２】米国特許６２６９１７６号公報
【特許文献３】特開２００９−１９５５１２号公報
【発明の概要】
【発明が解決しようとする課題】
【０００８】
ここで、画像診断に適した高品質な画像を得るためには、被写体を表す成分をできるだけ除去しないで、モアレ等の周期的パターンを除去したいという要求がある。しかし、特許文献１に記載された方法は、グリッド像の縞模様のように一次元ハイパスフィルタ処理を施すことにより周期的パターンに対応する空間周波数成分を検出するものであるため、周期的パターンが高周波数帯域ではなく低周波数帯域に存在した場合には意図したとおりに周期的パターンを検出することができず、特許文献１に記載された方法を適用することができなかった。
【０００９】
また、特許文献２に記載された方法は、平滑化フィルタリング処理が施されることにより局所的なノイズ成分が抑制されるものの、検出された周期的パターンに対応する空間周波数成分のみを画像から除去するものではないため、周期的パターンを含む空間周波数成分を集中的に除去することはできなかった。
【００１０】
また、特許文献３に記載された方法によれば、２次元フーリエ変換により周期的パターンに対応する空間周波数成分を検出するため計算負荷が高く、より少ない計算負荷で周期的パターンに対応する空間周波数成分の通過帯域特性を決定する方法が求められていた。
【００１１】
本発明は上記の事情に鑑みてなされたものであり、周期的パターンに対応する周波数成分が中周波数帯域または低周波数帯域に存在する場合にも、少ない計算負荷で適切に通過帯域特性を決定することができる画像処理装置および画像処理方法ならびに画像処理プログラムを提供することを目的とする。
【課題を解決するための手段】
【００１２】
本発明による画像処理装置は、周期的パターンを含む画像を表す画像信号から少なくとも１つの方向について周波数解析を行うことにより、前記周期的パターンに対応する周波数スペクトルのピーク位置とピーク幅を検出する周期的パターン検出手段と、前記検出されたピーク位置およびピーク幅に基づいて、前記画像信号から前記周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定する通過帯域特性決定手段とを備えたことを特徴とするものである。
【００１３】
本発明による画像処理方法は、周期的パターンを含む画像を表す画像信号から少なくとも１つの方向について周波数解析を行うことにより、前記周期的パターンに対応する周波数スペクトルのピーク位置とピーク幅を検出し、前記検出されたピーク位置およびピーク幅に基づいて、前記画像信号から前記周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定することを特徴とするものである。
【００１４】
本発明による画像処理プログラムは、コンピュータを、周期的パターンを含む画像を表す画像信号から少なくとも１つの方向について周波数解析を行うことにより、前記周期的パターンに対応する周波数スペクトルのピーク位置とピーク幅を検出する周期的パターン検出手段と、前記検出されたピーク位置およびピーク幅に基づいて、前記画像信号から前記周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定する通過帯域特性決定手段として機能させることを特徴とするものである。
【００１５】
上記「周期的パターン」とは、原画像が含む周期的なパターンを持ったノイズを意味する。例えば、上記グリッドを利用して放射線画像を蓄積性蛍光体シートに撮影したときに原画像に含まれるグリッド像やモアレ等を意味する。
【００１６】
また、「通過帯域特性」は、所定の帯域の通過または除去を行うフィルタの特性を意味し、帯域通過特性（バンドパス特性）と帯域除去特性（バンドストップ特性）の両方を含む。
【００１７】
また、ピーク位置は、周波数スペクトルのピーク周波数の代表的な値を表すものであれば任意に定義してよい。同様に、ピーク幅も、周波数スペクトルのピークの幅の代表的な値を表すものであれば任意に定義してよい。
【００１８】
上記通過帯域特性決定手段は、周期的パターンに対応する周波数成分のピーク位置を含み、ピーク位置を通過する帯域に含み、ピークの位置とピーク幅に基づいて決定可能なものであればいかなる通過帯域特性を決定してもよい。
【００１９】
なお、周期的パターンに対応する周波数成分を除去するために、帯域通過の幅はピーク幅以上であることが好ましい。例えば、前記通過帯域特性決定手段は、前記ピーク位置を含み、前記ピーク幅が大きくなるほど通過帯域幅が大きくなるように前記通過帯域特性を特定するものであることが好ましい。
【００２０】
また、通過帯域特性において、低周波数帯域は被写体成分を多く含むため、低周波数帯域においてはフィルタの帯域通過が広いほど周期的パターンに対応する周波数成分と一緒に被写体成分が除去される可能性が高い。このことにより、低周波帯域においては、周期的パターンのピークを含み、かつ幅が狭い通過帯域特性が好ましい。例えば、通過帯域特性決定手段は、前記ピーク位置が低周波数帯域に位置するほど通過帯域の幅が狭くなるように前記通過帯域特性を決定するものであってもよい。
【００２１】
また、通過帯域特性決定手段は、前記ピーク位置と前記ピーク幅により規定される正規分布関数により前記通過帯域特性を特定するものであってもよい。例えば、正規分布関数を、前記ピーク位置を平均とし、前記ピーク幅の半値を標準偏差として規定することができ、かかる正規分布関数のピークを周波数レスポンスの１に対応付けるように正規分布曲線を拡大縮小したものを通過帯域特性として決定してもよい。
【００２２】
また、前記通過帯域特性決定手段は、前記通過帯域特性を、前記ピーク幅以上の幅を有し、前記ピーク位置を中心とする矩形関数により決定してもよい。
【００２３】
また、前記通過帯域特性決定手段は、前記ピークの高さにも基づいて前記通過帯域特性を決定するものであってもよい。スペクトルの大きさが小さい場合などには、周期的パターンに対応する周波数成分を含む周波数帯域を全て抑制しなくとも（周期的パターンに対応する周波数成分を一部画像信号に残留させても）、診断画像として十分使用可能な品質まで周期的パターンを抑制することができるためである。例えば、ピークの高さが大きくなる程帯域通過フィルタによる帯域通過の周波数レスポンスが大きくなり、ピークの高さが小さくなる程帯域通過フィルタによる帯域通過の周波数レスポンスが小さくなるように通過帯域特性を決定することが考えられる。
【００２４】
また、本画像処理装置において、前記画像信号から前記決定された通過帯域特性に基づく１次元フィルタリング処理を行うことにより前記周期的パターンを抑制する抑制手段をさらに備えることが好ましい。
【００２５】
また、本画像処理装置における前記通過帯域特性決定手段は、前記通過帯域特性における通過帯域の幅に上限または下限を設けたものであってもよい。非常に狭い帯域のみを通過（または除去）する帯域通過（または帯域除去）フィルタによりフィルタリング処理を行うことは、画像のアーチファクト等の発生原因となるため好ましくなく、一方、あまりに広帯域を通過する帯域通過（帯域除去）フィルタによりフィルタリング処理を行うことは、被写体成分等の画像診断上重要な情報までも周期的パターンに対応する成分と一緒に除去してしまう可能性があるため好ましくないためである。
【００２６】
前記通過帯域特性決定手段は、前記画像の撮影条件に応じて前記通過帯域特性を異ならせて決定するものであってもよい。なお、画像の撮影部位や撮影条件は、原画像の撮影装置において指定される上記撮影条件や撮影部位の情報などの撮影メニューを利用して取得してもよく、ユーザのマニュアル操作により入力された被写体の臓器や病変などの撮影部位情報を利用してもよい。また、本発明による画像処理装置は、前記画像から該画像の撮影部位を抽出する部位抽出手段をさらに備え、抽出された撮影部位情報を利用してもよい。例えば、本出願人により提案されている、特開２００２−１０９５４８号公報、特開２００３−６６６１号公報に記載された基準となる平均的な心胸郭の輪郭と略相似形のテンプレートを用いてテンプレートマッチング処理を行なって胸郭を自動検出する方法を用いてもよい。
【発明の効果】
【００２７】
本発明の画像処理装置、画像処理方法および画像処理プログラムによれば、周期的パターンを含む画像を表す画像信号から少なくとも１つの方向について周波数解析を行うことにより、周期的パターンに対応する周波数スペクトルのピーク位置とピーク幅を検出し、検出されたピーク位置およびピーク幅に基づいて、画像信号から周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定するため、ピーク位置およびピーク幅に基づいて、画像の被写体成分を過剰に抑制しない範囲で周波数成分を適切に抑制できる１次元フィルタの通過帯域特性を好適に決定できる。この結果、比較的小さい計算負荷で高品質な画像を得ることができる。
【図面の簡単な説明】
【００２８】
【図１】放射線画像撮影装置の概略を示した図
【図２】グリッド撮影により得られた放射線画像を示す図
【図３】放射線画像読取装置の一例を示す斜視図
【図４】走査方向と読取画像との関係を示した図
【図５】本実施形態における画像処理装置の概略構成図
【図６】本実施形態における画像処理装置の処理の流れを示す図
【図７】本実施形態における領域抽出処理を説明する概略図
【図８】本実施形態における規格化処理を説明する図
【図９Ａ】本実施形態における通過帯域特性決定処理を説明する図
【図９Ｂ】本実施形態における通過帯域特性決定処理の変形例を説明する図
【図１０Ａ】本実施形態における抑制処理を説明する図
【図１０Ｂ】本実施形態における抑制処理の変形例を説明する図
【図１１】本実施形態における通過帯域特性の変形例を示す図（その１）
【図１２】本実施形態における通過帯域特性の変形例を示す図（その２）
【図１３】本実施形態における通過帯域特性の変形例を示す図（その３）
【発明を実施するための形態】
【００２９】
以下、図面を参照して本発明の実施の形態を説明する。尚、以下の実施の形態では、本発明における周期的パターン抑制処理装置を放射線画像読取装置に用いた場合について説明するが、デジタルカメラ等を用いた通常撮影において網戸やブラインド等を介して撮影した際に撮影画像に含まれる周期的パターンを抑制するための画像処理装置等に用いられてもよい。尚、放射線画像読取装置は、蓄積性蛍光体シートに記録された人体の放射線画像をレーザビーム走査によりデジタル画像信号として読み取るものである。
【００３０】
図１は放射線画像撮影装置の概略図である。放射線源１から放射された放射線２は、被写体３を透過して静止グリッド（以下、単に「グリッド」と表記する）４に到達する。グリッド４は放射線２を吸収する鉛４ａと、放射線２を透過するアルミニウム４ｂとが例えば４本／ｍｍ程度のピッチで交互に配置されているものである。また、放射線２がアルミニウム４ｂを透過して蓄積性蛍光体シート１１に入射するように鉛４ａは位置に応じて多少傾きを変化させて設置されている。
【００３１】
従って、被写体３を透過した放射線２は、鉛４ａに吸収されて蓄積性蛍光体シート１１への到達を遮られる一方で、アルミニウム４ｂを透過して蓄積性蛍光体シート１１に到達する。蓄積性蛍光体シート１１には被写体像と共に４本／ｍｍの縞模様状のグリッド像が記録される。一方、被写体３内で散乱された散乱放射線２ａは、位置に応じて傾きを持って設置された鉛４ａによって吸収され、又はグリッド４の表面で反射するため、蓄積性蛍光体シート１１には到達しない。従って、蓄積性蛍光体シート１１は散乱放射線２ａの照射の少ない鮮明な放射線画像を記録することができる。図２は、図１に示した放射線画像撮影装置を用いて撮影した際に蓄積性蛍光体シート１１に蓄積記録された被写体像５とグリッド像に起因するモアレパターン６の放射線画像の一例である。
【００３２】
図３は、放射線画像読取装置の斜視図と機能ブロック図を組み合わせた図である。読取部１０の所定位置にセットされた蓄積性蛍光体シート１１は、駆動手段（不図示）によって駆動されるエンドレスベルト等のシート搬送手段１５によって、例えば走査ピッチ１０本／ｍｍで矢印Ｙ方向に搬送（副走査）される。一方、レーザ光源１６から発せられた光ビーム１７は、モータ２４により駆動され矢印方向に高速回転する回転多面鏡１８によってｆθレンズ等の集束レンズ１９に向けて反射される。光ビーム１７は集束レンズ１９を通過した後、ミラー２０によって蓄積性蛍光体シート１１に向けて反射される。蓄積性蛍光体シート１１は光ビーム１７によって副走査方向（矢印Ｙ方向）と略直角な矢印Ｘ方向に主走査される。
【００３３】
蓄積性蛍光体シート１１において光ビーム１７が照射された箇所からは蓄積記録されている放射線画像情報に応じた光量の輝尽発光光２１が発散される。輝尽発光光２１は光ガイド２２の入射端面２２ａから入射し、光ガイド２２内の内部を全反射を繰り返して出射端面２２ｂから出射する。出射された輝尽発光光２２はフォトマルチプライヤ２３によって受光されて光電変換され、アナログ画像信号Ｓａに変換される。
【００３４】
アナログ画像信号Ｓａはログアンプ２６によって対数的に増幅された後、Ａ／Ｄ変換器２８において空間周波数ｆｓ＝１０ｃｙｃｌｅ／ｍｍに対応するサンプリング間隔でサンプリングされてデジタル化され、デジタル画像信号Ｓｄ（以下、単に「画像信号Ｓｄ」と表記する）を出力する。この画像信号Ｓｄは、図４に示すように、蓄積性蛍光体シート１１に対して主走査方向（横方向）に光ビーム１７を走査させながら蓄積性蛍光体シート１１を副走査方向（縦方向）に移動させることにより蓄積性蛍光体シート１１を２次元走査して得られた放射線画像情報を示す。尚、このようにして得られた画像信号Ｓｄには、被写体３に対応する放射線画像５の情報の他に、グリッド４に対応するグリッド像（不図示）に起因するモアレパターン６の情報も含まれている。
【００３５】
画像信号Ｓｄは一旦記憶部２９に記憶された後、画像信号処理部３０に入力される。画像信号処理部３０は本発明における画像処理方法を実行するための画像処理装置４０を備えている。
【００３６】
図５は本実施の形態の画像処理装置４０の概略構成を表す図である。図５を用いて、本実施形態による画像処理装置４０について以下説明する。
【００３７】
本実施形態の画像処理装置４０は、周期的パターンを含む画像を表す画像信号Ｓｄから少なくとも１つの方向について周波数解析を行うことにより、周期的パターンに対応する周波数スペクトルのピーク位置とピーク幅を検出する周期的パターン検出手段４１と、検出されたピーク位置およびピーク幅に基づいて、画像信号から周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定する通過帯域特性決定手段４２と、画像信号から決定された通過帯域特性に基づく１次元フィルタリング処理を行うことにより周期的パターンを抑制する抑制手段４３とを備えるものである。
【００３８】
本実施形態における画像処理プログラムと画像処理プログラムが参照するデータは、インストール時に記憶部２９に記憶され、起動時に記憶部２９に含まれるメモリにロードされる。画像処理プログラムは、画像処理装置４０を構成する画像信号処理部３０の中央処理装置に実行させる処理として、周期的パターン検出処理と、通過帯域特性決定処理と、抑制処理を規定している。そして、プログラムの規定にしたがって、中央処理装置が上記各処理を実行することにより、画像信号処理部３０の中央処理装置は、周期的パターン検出手段４１、通過帯域特性決定手段４２、抑制手段４３として機能する。
【００３９】
図６は、本画像処理装置４０の画像処理の流れを表すフローチャートである。図７は周期的パターン検出手段４１によるライン状のサンプル領域抽出処理を説明するための図である。図６および図７に従って、本実施形態の画像処理装置４０の画像処理について詳細に説明する。
【００４０】
まず、図６に示すように、画像データＳｄが周期的パターン検出手段４１に入力される（Ｓ０１）。すると、図７に示すように、周期的パターン検出手段４１は、画像データＳｄにより表される画像におけるｘ方向およびｙ方向のそれぞれについて、ｎｘｍ本のライン状の領域を抽出する。なお、ｘ方向とｙ方向とで領域の抽出の仕方は同一であるため、ここでは、ｘ方向の抽出についてのみ説明する。
【００４１】
周期的パターン検出手段４１は、まず画像データＳｄにより表される画像Ｇ１上にｘ方向に延在するｎ（本実施形態ではｎ＝９）の矩形領域Ａ１〜Ａｎを抽出する（Ｓ０２）。そして、ｎ個の領域画像データＳＲｘｉのそれぞれに、ライン状のｍ個（本実施形態ではｍ＝９）の領域画像データＳＲｘｉｊ（ｉ＝１〜ｎ,ｊ＝１〜ｍ）を得る。ここで、領域Ａ１のｙ方向の中心は画像上縁から画像高さの１／３０の位置に配置され、領域Ａ９のＹ方向の中心は画像上縁から画像高さの１／３０の位置に配置され、領域Ａ２乃至Ａ７は領域Ａ１およびＡ９間に等間隔に位置するように配置される。また、ライン状のｍ個（本実施形態ではｍ＝９）の各領域画像データＳＲｘｉｊ（ｉ＝１〜ｎ,ｊ＝１〜ｍ）は、ｙ方向に３画素間隔で配置され、ｘ方向への画素数は１０２４である。
【００４２】
周期的パターン検出手段４１は、ＳＲｘｉｊ（ｉ＝１〜ｎ,ｊ＝１〜ｍ）のそれぞれに高速フーリエ変換を施して周波数解析を行う（Ｓ０３）。具体的には各領域画像データＳＲｘｉｊの画素配列ＳＲｘｉｊ（ｋ）（ｋ＝０〜１０２３）に対して高速フーリエ変換が施される。すなわち、スペクトルをＰｉｊ（ｋ）（ｋ＝０〜１０２３）、高速フーリエ変換により得られる演算結果の実部をＲｅ（ｋ）、虚部をＩｍ（ｋ）とすると、下記の式（１）によりスペクトルＰｉｊ（ｋ）が算出される。
【００４３】
【数１】
なお、スペクトルは実際にはナイキスト周波数以下の周波数についてのみ算出されるため、データ数としては１０２４の１／２すなわち５１２となる。したがって、以降の説明においては、スペクトルの参照符号としてＰｉｊ（ｋ′）（ｋ′＝０〜５１１）を用いるものとする。そして、周期的パターン検出手段４１は、ｎｘｍ個の領域画像データＳＲｘｉｊについてそれぞれ算出された、ｍｘｎ個のスペクトルＰｉｊ（ｋ′）（ｉ＝１〜ｎ,ｊ＝１〜ｍ）を加重平均処理したスペクトルＰａ（ｋ’）を算出する（Ｓ０４）。
【００４４】
図８左図はスペクトルＰａ（ｋ’）の例を示す図である。図８左図に示すように求められたスペクトルは低周波数から高周波数となるにつれて小さくなるとともに、所々に細かなピークを有するものとなっている。
【００４５】
次いで、周期的パターン検出手段４１は、加重平均処理したスペクトルＰａ（ｋ’）を規格化する（Ｓ０５）。周期的パターン検出手段４１は、スペクトルＰａ（ｋ’）を平滑化してスペクトルＰａ（ｋ’）の細かなピークが平滑化により除去されたスペクトルＰｂ（ｋ’）を算出する。そして算出されたスペクトルＰｂ（ｋ’）を図８左図に示す画像信号Ｓｄから算出されたスペクトルＰａ（ｋ’）から減算することにより規格化し、図８右図の規格化されたスペクトルＰ（ｋ’）を取得する。規格化処理により、ピークをより明確にして精密に評価するためである。
【００４６】
次いで、周期的パターン検出手段４１は、規格化されたスペクトルＰからスペクトルの最大値における周波数である最大値周波数（ピーク位置）とピーク幅を取得する（Ｓ０６）。図９Ａは、通過帯域特性を決定する処理を説明する図であり、図８右の規格化されたスペクトルＰを拡大したものである。まず、周期的パターン検出手段４１は、図８右および図９Ａに示すように、スペクトルＰについて、誤検出を避けるため所定の帯域（ｆＴｈ１≦ｆ≦ｆＴｈ２）を評価帯域とし、最大スペクトルを示した周波数（ピーク位置）ｐｆとピークの最大値（ピーク高さｐｈ）を取得する。そして、評価帯域においてピーク位置を含み所定のしきい値Ｔｈ以上の周波数帯域の幅をピーク幅ｐｗとして取得し、取得されたピーク位置ｐｆ、ピーク幅ｐｗ、ピーク高さｐｈをメモリに保存する。なお、本実施形態では、図８右図に示すように、ｆＴｈ１＝０．３（ｃｙｃｌｅ／ｍｍ）、ｆＴｈ２＝３．３（ｃｙｃｌｅ／ｍｍ）である。かかる評価帯域は、被写体成分を多く含む低周波数を除き、グリッド像の折り返し成分であるモアレパターン６に対応する周波数成分が存在すると予想される、ナイキスト周波数以下の周波数帯域を含むように設定されている。
【００４７】
次いで、通過帯域特性決定手段４２は、取得したピーク位置ｐｆおよびピーク幅ｐｗに基づいて通過帯域特性を決定する（Ｓ０７）。本実施形態においては、平均μをｐｆ、標準偏差σをｐｗ／２とすることにより規定される、下記の式（２）の正規分布関数Ｎに基づいて通過帯域特性（フィルタ特性）を特定する。すなわち、図９Ａに示されるように、スペクトルＰにおける周期的パターン成分に対応するピーク位置ｐｆで周波数レスポンスが１となるように正規分布関数Ｎを拡大縮小した関数が通過帯域特性Ｆとして決定される。本実施形態では、通過帯域特性決定手段４２は同様にｙ方向についても通過帯域特性を決定する。
【数２】
【００４８】
次いで、抑制手段４３は、決定されたｘ、ｙ両方向の通過帯域特性に基づいて、それぞれ周知の手法を用いて１次元の帯域通過フィルタ４４を作成する（Ｓ０８）。例えば、フィルタ特性に基づくディジタルフィルタの設計については、高橋・池原、「ディジタルフィルタ」、培風館（特に第6.3節、pp84-86）、岩田・実践、「ディジタル・フィルタ設計入門」、ＣＱ出版（特に第8章）などに記載された方法が適用できる。
【００４９】
図１０Ａは、本実施形態の抑制手段４３による周期的パターンの抑制処理を説明する図である。本実施形態による抑制手段４３は、図９Ａに示すように、決定した通過帯域特性Ｆに基づいて、決定された通過帯域を除去する通過帯域特性Ｆ’を有する帯域除去フィルタ（バンドストップフィルタ）ＢＳＦを１次元の帯域通過フィルタ４４として生成する。
【００５０】
そして、抑制手段４３は、画像信号から、決定された通過帯域特性に基づく１次元フィルタリング処理を行うことにより周期的パターンを抑制する（Ｓ０９）。すなわち、本実施形態では、ｘ、ｙ方向ごとに作成された１次元フィルタ（ＢＳＰ）によるｘ、ｙ方向に対するフィルタリング処理を施す。
【００５１】
ここでは、図１０Ａに示すように、抑制手段４３は、画像信号Ｓｄから生成した帯域除去フィルタＢＳＦによるフィルタリング処理を行い、周期的パターンを抑制した画像である、第２処理済信号Ｓｐを抽出する。そして、この第２処理済信号Ｓｐに対応する画像が生成されて（Ｓ１０）、生成された画像が出力される（Ｓ１１）。その後、出力された画像は、モニタ等の表示手段（不図示）に表示され、又はプリンタ等の出力手段（不図示）に出力される。
【００５２】
本実施形態によれば、検出されたピーク位置ｐｆおよびピーク幅ｐｗに基づいて、画像信号から周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定するため、ピーク位置ｐｆおよびピーク幅ｐｗに基づいて、画像の被写体成分を過剰に抑制しない範囲で周波数成分を適切に抑制する１次元フィルタの通過帯域特性を好適に決定できる。このため、本実施形態による画像処理を低周波数帯域に存在する周期的パターンにも好適に適用できる。また、本実施形態では、周期的パターンに対応する空間周波数成分の検出処理および周期的パターンに対応する空間周波数成分の抑制処理を１次元の処理を行って、周波数成分を適切に抑制する１次元フィルタの通過帯域特性を決定する。このため、比較的小さい計算負荷にも関わらず、決定された通過帯域特性に基づいて好適に周期的パターンの抑制処理を行うことができる。また、集中的に周期的パターンに対応する空間周波数のみを抑制できるため高品質な画像を得ることができる。
【００５３】
また、本実施形態による通過帯域特性決定手段４２が、ピーク位置ｐｆを平均とし、ピーク幅ｐｗの半値を標準偏差とする正規分布関数に基づいて、通過帯域特性Ｆを特定するものであるため、決定された通過帯域特性Ｆはピーク位置ｐｆを中心としピーク幅ｐｗを適切に反映して滑らかに分布する特性を示す。すなわち、本実施形態の通過帯域特性Ｆは、ピーク位置を含み、ピーク幅が大きくなるほど通過帯域幅が大きくなるものであるため、周期的パターンに対応する空間周波数成分のピーク位置ｐｆを包含するような通過帯域特性Ｆを適切に決定することができる。さらに、本実施形態によれば、ピーク位置を中心にピーク幅ｐｗを含むような通過帯域特性Ｆが決定される可能性が高く好適に周期的パターンに対応する空間周波数成分を含む通過帯域特性Ｆを決定できる。
【００５４】
また、本実施形態による画像処理装置４０は、画像信号から決定された通過帯域特性に基づく１次元フィルタリング処理を行うことにより周期的パターンを抑制する抑制手段４３をさらに備えているため、決定された通過帯域特性Ｆに基づいて、好適に周期的パターンに対応する空間周波数成分を抑制できる。
【００５５】
以下、本実施形態の変形例について記す。
【００５６】
上記実施形態において周期的パターン検出手段４１は、ピーク位置、ピーク幅、ピーク高さについて、ピークの位置、幅、高さのそれぞれの代表的な値を定義するものであれば、任意の定義により特定してよい。例えば、ピーク幅を所定値以上の周波数レスポンスを有する帯域の数倍（２乃至３倍）に設定してもよい。
【００５７】
また、通過帯域特性決定手段４２は、周期的パターンに対応する周波数成分のピーク位置を含み、ピークの位置とピーク幅に基づいて決定可能な通過帯域特性を決定するものであればいかなる通過帯域特性を決定してもよい。図９Ｂは、通過帯域特性を矩形関数で規定する例を示す図である。例えば、図９Ｂに示すように、ピーク幅ｐｗと同じ周波数帯域を通過させ、ピーク位置ｐｆを中心とする矩形関数により通過帯域特性Ｆ１を決定してもよい。
【００５８】
また、本発明により決定された通過帯域特性は、ピーク位置ｐｆとピーク幅ｐｗに基づいて決定された通過帯域特性に種々の補正を行ったものであってもよい。以下に、本発明による通過帯域特性の変形例を説明する。
【００５９】
一般に、低周波数帯域は被写体成分を多く含むため、低周波数帯域においてはフィルタの帯域通過が広いほど周期的パターンに対応する周波数成分と一緒に被写体成分が除去される可能性が高い。このことにより、低周波帯域においては、周期的パターンのピークを含み、かつ幅が狭い通過帯域特性が好ましい。図１１は、通過帯域特性の変形例を示す図である。本実施形態の変形例として、図１１に示すように、上記実施形態における通過帯域特性決定手段４２は、上記通過帯域特性Ｆが周期的パターンのピーク位置ｐｆを含み、ピーク位置ｐｆが低周波数帯域に位置するほど通過帯域の幅が狭くなるように決定してもよい。例えば、ピーク位置ｐｆに通過帯域特性Ｆの通過帯域幅を対応付けて記憶部２９に記憶しておき、周期的パターン検出手段４１がピーク位置ｐｆを検出すると、検出されたピーク位置ｐｆ対応付けられた通過帯域幅になるように上記実施形態の通過帯域特性Ｆを幅方向に拡大縮小する補正を行い、その結果得られた関数を通過帯域特性Ｆ２として決定することが考えられる。この場合には、従来技術に比べて低周波成分に多く含まれる被写体成分を維持しつつモアレ等の周期的パターンに対応する空間周波数成分のみを好適に除去できる。
【００６０】
また、非常に尖鋭な帯域のみを通過（または除去）する帯域通過（帯域除去）フィルタによりフィルタリング処理を行うことは、画像のアーチファクト等の発生原因となるため好ましくない。そこで、通過帯域特性決定手段４２は、通過帯域の幅に下限を設けてもよい。図１２は、上記実施形態の通過帯域特性の変形例を示す図である。例えば、周期的パターンがあまりに鋭いピークである場合など、上記通過帯域特性Ｆの幅が所定値より小さい場合に、通過帯域の幅が広くなるように上記実施形態の通過帯域特性Ｆを補正した関数を通過帯域特性として決定してもよい。例えば、図１２に示すように、ピーク位置ｐｆおよびピーク幅ｐｗに基づいて決定された通過帯域特性Ｆ３ａの所定の周波数レスポンスＲＴｈにおける通過帯域幅ｆｗが所定値より小さい場合、通過帯域特性Ｆの幅を所定の割合で広くなるように補正したものを通過帯域特性Ｆ３として決定してもよい。または、ピーク幅ｐｗを所定の割合で大きくなるように補正し、補正したピーク幅に基づいて通過帯域特性を決定するものであってもよい。
【００６１】
上記場合には、通過帯域が狭い通過帯域特性に基づいてフィルタリング処理をおこなうことによるアーチファクトの発生を抑制でき、より品質の高い画像を生成することができる。
【００６２】
また、あまりに広帯域を通過する帯域通過（帯域除去）フィルタによりフィルタリング処理を行うことは、被写体成分等の画像診断上重要な情報までも周期的パターンに対応する成分と一緒に除去してしまう可能性があるため好ましくない。そこで、通過帯域特性決定手段４２は、通過帯域の幅に上限を設けてもよい。例えば、通過帯域特性Ｆの幅ｆｗが所定値より大きい場合に、通過帯域特性の幅が狭くなるように通過帯域特性Ｆを補正したものを通過帯域特性として決定してもよい。例えば、図１２の通過帯域特性Ｆ４に示すように、ピーク位置ｐｆおよびピーク幅ｐｗに基づいて決定された通過帯域特性Ｆ４ａの所定の周波数レスポンスＲｔｈにおける通過帯域幅ｆｗが所定値より大きい場合、通過帯域特性Ｆ４ａの幅を所定の割合で狭くなるように補正したものを通過帯域特性Ｆ４として決定してもよい。または、ピーク幅ｐｗが所定値よりも大きい場合に、ピーク幅ｐｗを所定の割合で小さくした補正ピーク幅ｐｗにより通過帯域特性を決定してもよい。また、ピーク幅ｐｗを用いて算出された通過帯域特性Ｆ４ａに対して、ピーク位置ｐｆから所定値ｆａだけ低い周波数より小さい周波数帯域における周波数レスポンスを０にするように通過帯域特性Ｆ４ａを補正したものを通過帯域特性Ｆ４ｂとして決定してもよい。すなわち、図１２に示すように、周波数（ｐｆ−ｆａ）より小さい周波数帯域における周波数レスポンスを０、周波数（ｐｆ−ｆａ）以上の周波数レスポンスをＦ４ａで規定するような通過帯域特性Ｆ４ｂにしてもよい。また、同様に低周波数側だけでなく高周波数側に対しても通過帯域を狭くする補正を行ってもよい。例えば、上記通過帯域特性Ｆ４ａまたはＦ４ｂをピーク位置から所定値ｆｂだけ高い周波数（ｐｆ＋ｆｂ）より大きい周波数帯域における周波数レスポンスを０とし、周波数（ｐｆ＋ｆｂ）以下の周波数レスポンスを通過帯域特性Ｆ４ａにより規定する通過帯域特性Ｆ４ｃに補正してもよい。
【００６３】
上記場合には、通過帯域が広い通過帯域特性に基づいてフィルタリング処理をおこなうことによる被写体成分等の画像診断上重要な情報までも周期的パターンに対応する成分と一緒に除去してしまう可能性を低減でき、より品質の高い画像を生成することができる。
【００６４】
また、通過帯域特性決定手段４２は、ピークの高さ（スペクトルの大きさ）にも基づいて通過帯域特性を決定するものであってもよい。図１３は、本実施形態の変形例を示す図であり、ピークの高さｐｈ（スペクトルの大きさ）が小さいほど通過帯域特性の周波数レスポンスが小さくなるようにピーク位置ｐｆとピーク幅ｐｗに基づいて決定された通過帯域特性Ｆ５ａを補正して通過帯域特性Ｆ５を決定した例を説明するための図である。図１３の通過帯域特性Ｆ５に示すように、スペクトルの高さｐｈが小さい場合などには、周期的パターンに対応する周波数成分を含む周波数帯域を全て抑制しなくとも（周期的パターンに対応する周波数成分を一部画像信号に残留させても）、診断画像として十分使用可能な品質まで周期的パターンを抑制することができるためである。この場合、ピークの高さｐｈが大きくなる程帯域通過フィルタによる帯域通過の周波数レスポンスが大きくなり、ピークの高さｐｈが小さくなる程帯域通過フィルタによる帯域通過の周波数レスポンスが小さくなるものであれば、任意の方法により通過帯域特性を決定してよい。
【００６５】
また、通過帯域特性決定手段４２は、画像の撮影条件に応じて通過帯域特性を異ならせて決定するものであってもよい。例えば、記憶部２９に、撮影部位や撮影条件に適応する通過帯域特性を予め決定して撮影部位や撮影条件ごとに対応付けて記憶しておき、通過帯域特性決定手段４２が、処理対象である画像の撮影部位や撮影条件に対応付けられた通過帯域特性を特定し、特定された通過帯域特性を用いて周期的パターンの抑制処理を行うことが考えられる。
【００６６】
上記場合には、撮影部位や撮影条件に応じた通過帯域特性に応じて適切に抑制処理を施すことができるため、好適に周期的パターンに対応する周波成分を抑制した第２処理済信号を抽出することができる。結果として、より高品質に周期的パターンに対応する周波成分を抑制した画像を生成することができる。
【００６７】
なお、画像の撮影部位や撮影条件は、原画像の撮影装置において指定される上記撮影条件や撮影部位の情報などの撮影メニューを利用して取得してもよく、ユーザのマニュアル操作により入力された被写体の臓器や病変などの撮影部位情報を利用してもよい。また、上記画像処理装置４０は、不図示の画像の撮影部位を抽出する部位抽出手段をさらに備え、抽出された撮影部位情報を利用してもよい。撮影部位情報の抽出方法として、部位又は病変などが抽出できるものであればいかなる方法を適用してもよい。例えば、本出願人により提案されている、特開２００２−１０９５４８号公報、特開２００３−６６６１号公報に記載された基準となる平均的な心胸郭の輪郭と略相似形のテンプレートを用いてテンプレートマッチング処理を行なって胸郭を自動検出する方法を用いてもよい。
【００６８】
なお、かかる場合には、通過帯域特性決定手段４２は、ユーザによる撮影部位の入力操作を必要とせず、自動抽出された結果に基づいて撮影部位を抽出でき、抽出された撮影部位に基づいて適切に通過帯域特性を決定することができる。このため、ユーザは撮影部位に応じて周期的パターンに対応する周波成分を抑制した画像を容易かつ高品質に生成することができる。また、通過帯域特性決定手段４２が撮影メニューから撮影部位または撮影条件を自動抽出して、自動抽出された結果に基づいて通過帯域特性を決定した場合にも同じ効果が得られる。
【００６９】
また、同一の装置で続けて同じ部位の周期的パターンの抑制処理を行う場合など、ほぼ同じような撮影条件の場合には、最初の被写体の画像信号Ｓｄに対して本実施形態の画像処理により帯域通過フィルタ（帯域除去フィルタ）を決定し、かかるフィルタを後続の画像信号にも適用してもよい。ほぼ同じような撮影条件の場合には、同じような周期的パターンを有する画像が取得されると考えられるため、代表的な画像について決定した帯域通過フィルタ（帯域除去フィルタ）を用いることで十分周期的パターンの抑制効果が得られるとともに画像ごとに通過帯域特性を決定する時間と労力を省略できるため効率がよい。
【００７０】
なお、抑制手段４３は、最終的に画像信号Ｓｄに対して決定した通過帯域特性に基づいて周期的パターンの抑制処理を行った画像信号（本実施形態におけるＳｐ）が得られるものであれば、水平方向および垂直方向に帯域通過フィルタおよび帯域除去フィルタのいかなる組合せおよびフィルタリング処理の順序により抑制処理を行ってもよい。例えば、決定した通過帯域特性Ｆに基づいて、決定した通過帯域を通過させる通過帯域特性Ｆを有する帯域通過フィルタ（バンドパスフィルタ）ＢＰＦを生成してもよい。図１０Ｂは、上記実施形態の抑制処理の変形例を示す図である。図１０Ｂに示すように、画像信号Ｓｄを帯域通過フィルタＢＰＦによるフィルタリング処理した画像信号Ｓｅ（画像信号Ｓｄから抽出した周期的パターン画像）を生成し、画像信号Ｓｄから減算することにより上記実施形態と同様に周期的パターンを抑制できる。この場合、抑制手段４３は帯域通過フィルタＢＰＦと減算器４５を備え、画像信号Ｓｄに対して帯域通過フィルタＢＰＦによるフィルタリング処理を施して、周期的パターン成分を抽出した周期的パターン画像信号Ｓｅを算出し、その後、減算器４５は周期的パターン画像信号Ｓｅを減算して、周期的パターンを抑制した画像である、処理済信号Ｓｐを抽出する。
【００７１】
上記実施形態において周期的パターン検出手段４１は、複数のサンプル領域に変えて１つの代表的な領域についてのみ複数の線上の画像信号を取得して、周波数解析を行って加重平均をとってスペクトルＰａを取得してもよく、１つの線上の画像信号を取得して、これについてのみ周波数解析を行ったスペクトルを規格化されたスペクトルＰａの代わりに用いてもよい。
【００７２】
上記実施形態において周期的パターン検出手段４１は、規格化処理を省略してもよい。ただし、規格化処理を行った場合には、ピーク幅やピーク値やピーク高さなどの情報を精度よく検出でき好ましい。
【００７３】
また、本実施形態において、絶縁基板上に夫々が画素に対応する複数個の光電変換素子を２次元状に形成した放射線画像読取装置であって、この放射線画像読取装置上に形成された画像情報を担持する放射線が照射されると照射された放射線を可視光に変換する蛍光体層（シンチレータ）を積層して成る放射線固体検出器（以下、「光変換方式」の放射線固体検出器という）を用いた放射線画像読取装置に本発明を適用した例について説明したが、本発明は、上記実施形態に限定されず種々のタイプの放射線固体検出器を用いた放射線画像読取装置に適用可能である。例えば、絶縁基板上に夫々が画素に対応する複数個の電荷収集電極を２次元状に形成した放射線画像読取装置であって、この２次元画像読取装置上に形成された、画像情報を担持する放射線が照射されると前記画像情報を担持する電荷を発生する放射線導電体とを積層して成る放射線固体検出器（以下、「直接変換方式」の放射線固体検出器という）を用いた放射線画像撮影装置に適用してもよい。
【００７４】
また、光変換方式の放射線固体検出器としては、例えば特開昭59-211263 号、特開平2-164067号、ＰＣＴ国際公開番号ＷＯ92/06501号、Signal,noise,and read out considerations in the development of amorphous silicon photodiode arrays for radiotherapy and diagnostic x-ray imaging，L.E.Antonuk et.al ，University of Michigan，R.A.Street Xerox，PARC，SPIE Vol.1443 Medical Imaging V;Image Physics(1991) ，p.108-119 等に記載された放射線固体検出器が適用できる。
【００７５】
一方、直接変換方式の放射線固体検出器としては、例えば、(i) 放射線の透過方向の厚さが通常のものより10倍程度厚く設定された放射線固体検出器（MATERIAL PARAMETERS IN THICK HYDROGENATED AMORPHOUS SILICON RADIATION DETECTORS,Lawrence Berkeley Laboratory.University of California,Berkeley.CA 94720 XeroxParc.Palo Alto.CA 94304)、あるいは(ii)放射線の透過方向に、金属板を介して２つ以上積層された放射線固体検出器(Metal/Amorphous Silicon Multilayer Radiation Detectors,IEE TRANSACTIONS ON NUCLEAR SCIENCE.VOL.36.NO.2.APRIL 1989) 、あるいは(iii) ＣｄＴｅ等を使用した放射線固体検出器（特開平1-216290号）等に記載された放射線固体検出器が適用できる。
【００７６】
上記の各実施形態はあくまでも例示であり、上記のすべての説明が本発明の技術的範囲を限定的に解釈するために利用されるべきものではない。この他、上記の実施形態におけるシステム構成、ハードウェア構成、処理フロー、モジュール構成や具体的処理内容等に対して、本発明の趣旨から逸脱しない範囲で様々な改変を行ったものも、本発明の技術的範囲に含まれる。
【符号の説明】
【００７７】
１ 放射線源
２ 放射線
３ 被写体
４ グリッド
５ 被写体像
６ グリッド線
１０ 読取部
２９ 記憶部
３０ 画像処理部
４０ 画像処理装置
４１ 周期的パターン検出手段
４２ 通過帯域特性決定手段
４３ 抑制手段
４４ 帯域通過フィルタ
４５ 減算器
ＢＰＦ 帯域通過フィルタ
ＢＳＦ 帯域除去フィルタ
Ｆ，Ｆ１，Ｆ２，Ｆ３，Ｆ４，Ｆ５，Ｆ６ 通過帯域特性
Ｐ，Ｐａ，Ｐｂ 周波数スペクトル
ｐｆ ピーク位置（周波数）
ｐｗ ピーク幅
ｐｈ ピークの高さ
【特許請求の範囲】
【請求項１】
周期的パターンを含む画像を表す画像信号から少なくとも１つの方向について周波数解析を行うことにより、前記周期的パターンに対応する周波数スペクトルのピーク位置とピーク幅を検出する周期的パターン検出手段と、
前記検出されたピーク位置およびピーク幅に基づいて、前記画像信号から前記周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定する通過帯域特性決定手段とを備えたことを特徴とする画像処理装置。
【請求項２】
前記通過帯域特性決定手段は、前記ピーク位置を含み、前記ピーク幅が大きくなるほど通過帯域幅が大きくなるように前記通過帯域特性を特定するものであることを特徴とする請求項１記載の画像処理装置。
【請求項３】
前記通過帯域特性決定手段は、前記ピーク位置と前記ピーク幅により規定される正規分布関数により前記通過帯域特性を特定するものであることを特徴とする請求項２記載の画像処理装置。
【請求項４】
前記通過帯域特性決定手段は、前記前記ピーク位置が低周波数帯域に位置するほど通過帯域の幅が狭くなるように前記通過帯域特性を決定するものであることを特徴とする請求項１から３のいずれか１項記載の画像処理装置。
【請求項５】
前記通過帯域特性決定手段は、前記ピークの高さにも基づいて前記通過帯域特性を決定するものであることを特徴とする請求項１から４いずれか１項記載の画像処理装置。
【請求項６】
前記画像信号から前記決定された通過帯域特性に基づく１次元フィルタリング処理を行うことにより前記周期的パターンを抑制する抑制手段をさらに備えたことを特徴とする請求項１から５いずれか１項記載の画像処理装置。
【請求項７】
前記通過帯域特性決定手段は、前記通過帯域特性における通過帯域の幅に上限または下限を設けたものであることを特徴とする請求項１から６いずれか１項記載の画像処理装置。
【請求項８】
前記通過帯域特性決定手段は、前記画像の撮影条件に応じて前記通過帯域特性を異ならせて決定するものであることを特徴とする請求項１から７いずれか１項記載の画像処理装置。
【請求項９】
前記通過帯域特性決定手段は、前記通過帯域特性を、前記ピーク幅以上の幅を有し、前記ピーク位置を中心とする矩形関数により決定するものであることを特徴とする請求項１から８いずれか１項記載の画像処理装置。
【請求項１０】
周期的パターンを含む画像を表す画像信号から少なくとも１つの方向について周波数解析を行うことにより、前記周期的パターンに対応する周波数スペクトルのピーク位置とピーク幅を検出し、
前記検出されたピーク位置およびピーク幅に基づいて、前記画像信号から前記周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定することを特徴とする画像処理方法。
【請求項１１】
コンピュータを、
周期的パターンを含む画像を表す画像信号から少なくとも１つの方向について周波数解析を行うことにより、前記周期的パターンに対応する周波数スペクトルのピーク位置とピーク幅を検出する周期的パターン検出手段と、
前記検出されたピーク位置およびピーク幅に基づいて、前記画像信号から前記周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定する通過帯域特性決定手段として機能させることを特徴とする画像処理プログラム。
【請求項１】
周期的パターンを含む画像を表す画像信号から少なくとも１つの方向について周波数解析を行うことにより、前記周期的パターンに対応する周波数スペクトルのピーク位置とピーク幅を検出する周期的パターン検出手段と、
前記検出されたピーク位置およびピーク幅に基づいて、前記画像信号から前記周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定する通過帯域特性決定手段とを備えたことを特徴とする画像処理装置。
【請求項２】
前記通過帯域特性決定手段は、前記ピーク位置を含み、前記ピーク幅が大きくなるほど通過帯域幅が大きくなるように前記通過帯域特性を特定するものであることを特徴とする請求項１記載の画像処理装置。
【請求項３】
前記通過帯域特性決定手段は、前記ピーク位置と前記ピーク幅により規定される正規分布関数により前記通過帯域特性を特定するものであることを特徴とする請求項２記載の画像処理装置。
【請求項４】
前記通過帯域特性決定手段は、前記前記ピーク位置が低周波数帯域に位置するほど通過帯域の幅が狭くなるように前記通過帯域特性を決定するものであることを特徴とする請求項１から３のいずれか１項記載の画像処理装置。
【請求項５】
前記通過帯域特性決定手段は、前記ピークの高さにも基づいて前記通過帯域特性を決定するものであることを特徴とする請求項１から４いずれか１項記載の画像処理装置。
【請求項６】
前記画像信号から前記決定された通過帯域特性に基づく１次元フィルタリング処理を行うことにより前記周期的パターンを抑制する抑制手段をさらに備えたことを特徴とする請求項１から５いずれか１項記載の画像処理装置。
【請求項７】
前記通過帯域特性決定手段は、前記通過帯域特性における通過帯域の幅に上限または下限を設けたものであることを特徴とする請求項１から６いずれか１項記載の画像処理装置。
【請求項８】
前記通過帯域特性決定手段は、前記画像の撮影条件に応じて前記通過帯域特性を異ならせて決定するものであることを特徴とする請求項１から７いずれか１項記載の画像処理装置。
【請求項９】
前記通過帯域特性決定手段は、前記通過帯域特性を、前記ピーク幅以上の幅を有し、前記ピーク位置を中心とする矩形関数により決定するものであることを特徴とする請求項１から８いずれか１項記載の画像処理装置。
【請求項１０】
周期的パターンを含む画像を表す画像信号から少なくとも１つの方向について周波数解析を行うことにより、前記周期的パターンに対応する周波数スペクトルのピーク位置とピーク幅を検出し、
前記検出されたピーク位置およびピーク幅に基づいて、前記画像信号から前記周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定することを特徴とする画像処理方法。
【請求項１１】
コンピュータを、
周期的パターンを含む画像を表す画像信号から少なくとも１つの方向について周波数解析を行うことにより、前記周期的パターンに対応する周波数スペクトルのピーク位置とピーク幅を検出する周期的パターン検出手段と、
前記検出されたピーク位置およびピーク幅に基づいて、前記画像信号から前記周期的パターンに対応する空間周波数成分を抑制する１次元フィルタの通過帯域特性を決定する通過帯域特性決定手段として機能させることを特徴とする画像処理プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９Ａ】
【図９Ｂ】
【図１０Ａ】
【図１０Ｂ】
【図１１】
【図１２】
【図１３】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９Ａ】
【図９Ｂ】
【図１０Ａ】
【図１０Ｂ】
【図１１】
【図１２】
【図１３】
【公開番号】特開２０１２−２０３５０３（Ｐ２０１２−２０３５０３Ａ）
【公開日】平成２４年１０月２２日（２０１２．１０．２２）
【国際特許分類】
【出願番号】特願２０１１−６５５１０（Ｐ２０１１−６５５１０）
【出願日】平成２３年３月２４日（２０１１．３．２４）
【出願人】（３０６０３７３１１）富士フイルム株式会社 (25,513)
【Ｆターム（参考）】
【公開日】平成２４年１０月２２日（２０１２．１０．２２）
【国際特許分類】
【出願日】平成２３年３月２４日（２０１１．３．２４）
【出願人】（３０６０３７３１１）富士フイルム株式会社 (25,513)
【Ｆターム（参考）】
[ Back to top ]