(57)【要約】【課題】窓関数を用いる画像処理に際し、少ない時間および労力で最適な窓関数を決定するのに有利な画像処理方法を提供する。【解決手段】この画像処理方法は、画像を窓関数を用いて処理する窓関数処理工程（Ｓ３００）と、窓関数処理工程で処理された空間領域の画像を、フーリエ変換により空間周波数領域の成分に変換するフーリエ変換処理工程（Ｓ４００）と、空間周波数領域における周波数軸の軸上成分とローパスフィルタのカットオフ周波数とに基づいて、窓関数を決定する窓関数決定工程（Ｓ４１０）と、窓関数決定工程で決定された窓関数を用いて処理された空間周波数領域の成分を、ローパスフィルタを用いて処理するローパスフィルタ処理工程（Ｓ５００）と、ローパスフィルタ処理工程で処理された空間周波数領域の成分を、逆フーリエ変換により空間領域の画像に変換する逆フーリエ変換処理工程（Ｓ６００）とを含む。【選択図】図１
2015108872.tif 000002
【特許請求の範囲】
【請求項１】
画像を処理する方法であって、
前記画像を窓関数を用いて処理する窓関数処理工程と、
前記窓関数処理工程で処理された空間領域の画像を、フーリエ変換により空間周波数領域の成分に変換するフーリエ変換処理工程と、
前記空間周波数領域における周波数軸の軸上成分と、ローパスフィルタのカットオフ周波数とに基づいて、前記窓関数を決定する窓関数決定工程と、
前記窓関数決定工程で決定された前記窓関数を用いて処理された前記空間周波数領域の成分を、前記ローパスフィルタを用いて処理するローパスフィルタ処理工程と、
前記ローパスフィルタ処理工程で処理された前記空間周波数領域の成分を、逆フーリエ変換により空間領域の画像に変換する逆フーリエ変換処理工程と、
を含むことを特徴とする画像処理方法。
【請求項２】
前記窓関数処理工程の前に、前記画像の周辺部のデータを伸長するデータ伸長工程と、
前記逆フーリエ変換処理工程の後に、逆フーリエ変換された前記画像の周辺部からデータを切り出すデータ切り出し工程と、
を含むことを特徴とする請求項１に記載の画像処理方法。
【請求項３】
前記窓関数決定工程では、前記空間周波数領域において、前記周波数軸の軸上成分が閾値を横切る位置の空間周波数が、前記カットオフ周波数よりも小さくなる窓関数を選択することを特徴とする請求項１または２に記載の画像処理方法。
【請求項４】
前記窓関数決定工程では、前記空間周波数領域において、前記カットオフ周波数での前記周波数軸の軸上成分が、閾値よりも小さくなる窓関数を選択することを特徴とする請求項１または２に記載の画像処理方法。
【請求項５】
前記閾値は、前記画像の少なくとも１つについて、前記窓関数処理工程で処理された画像に対する前記空間周波数領域の成分における、前記カットオフ周波数での周波数軸の軸上成分と、前記逆フーリエ変換処理工程で変換された前記空間領域の画像とに基づいて決定されることを特徴とする請求項３または４に記載の画像処理方法。
【請求項６】
前記逆フーリエ変換処理工程で変換された前記空間領域の画像には、該空間領域の画像の周辺部を含むことを特徴とする請求項５に記載の画像処理方法。
【請求項７】
前記閾値は、前記画像の少なくとも１つについて、前記フーリエ変換処理工程で変換された前記空間周波数領域における前記周波数軸の軸上成分の振幅スペクトルの原点を除く最大振幅に基づいて決定されることを特徴とする請求項３または４に記載の画像処理方法。
【請求項８】
前記画像を処理する前に前記閾値が予め存在する場合に、
前記窓関数処理工程の前に、該閾値を取得する閾値取得工程を含むことを特徴とする請求項３ないし７のいずれか１項に記載の画像処理方法。
【請求項９】
前記窓関数処理工程の前に、前記閾値が存在するかどうかを判断する閾値設定判断工程と、
前記閾値設定判断工程で前記閾値が存在しないと判断した場合に前記閾値を設定する閾値設定工程と、
を含むことを特徴とする請求項３ないし７のいずれか１項に記載の画像処理方法。
【請求項１０】
前記カットオフ周波数は、前記空間周波数領域において、前記周波数軸の原点と、該原点から所定の半径におけるリング状の周波数成分の総和から極大値を示す半径との間の極小値に基づいて決定されることを特徴とする請求項１ないし９のいずれか１項に記載の画像処理方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、画像処理方法に関する。
【背景技術】
【０００２】
例えば撮像して得られた画像（画像データ）を処理するに際し、特にフーリエ変換を用いた空間周波数領域でローパスフィルタなどの各種処理を施す画像処理においては、一般に、有限の大きさを有する画像を周期画像として取り扱う。そのため、元の画像の四端の領域では、不連続部が生じ得る。この不連続性に起因するギブスの現象により、不連続部には高周波成分のノイズ、いわゆるリンギングノイズが発生し、画質を低下させる可能性がある。そこで、ギブスの現象を抑えるために、元の画像に対して窓関数処理を実行し、画像の四隅の領域に連続性を持たせる提案がなされている。
【０００３】
特許文献１は、デジタルカメラの合焦点を特定するために、撮像された２次元画像の特定エリアに対して、レンズ情報などに基づいて複数の窓関数の中から適切な窓関数を選択して窓関数処理を行う制御装置を開示している。この制御装置によれば、画像の端の不連続点に起因して不要な高周波成分が混入することによる合焦点特定の精度の低下を抑えることができる。また、特許文献２は、信号波形の処理区間について、例えばハニング窓関数に対して振幅ピーク点で引き伸ばした窓関数を生成し、結果的にウエイト関数を生成する方法を開示している。この方法によれば、波形の中心部はそのままで、波形の両端のみの影響を抑制するので、信号波形の処理区間においては、信号処理ロスが改善されて、信号解析処理の精度を向上させることができる。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開平１１−２９５５８９号公報
【特許文献２】特開２００６−２１１１２７号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、特許文献１では複数の窓関数から適切な窓関数を選択する、または特許文献２では窓関数の波形を引き伸ばすなどの記載はあるものの、実際に採用する窓関数が最適であるかどうかを判断する明確な基準は記載されていない。これに対して、現状、最終的に窓関数を決定するためには、基本的に以下のような処理が実行される。例えば、ローパスフィルタ処理を行うような画像処理の場合、その処理を実行する処理部は、まず、対象となる計測画像を取得し、窓関数処理を行って画像の周辺部の不連続性を抑える。次に、処理部は、フーリエ変換を行った後に、空間周波数領域でローパスフィルタ処理を行い、さらに逆フーリエ変換を行うことで、高周波成分を除去した画像を求める。ここで、窓関数が適切でなかった場合には、ローパスフィルタにより失われた高周波成分に起因して、画像にリンギングノイズが残留する。そこで、処理部は、その後、窓関数が適切であるかどうかを、画像にリンギングノイズが残存していないかで判断し、適切でなければ、窓関数を変更した上で、上記判断を繰り返す。すなわち、処理部は、ある窓関数が適切であるかを判断する前に、その都度、フーリエ変換および逆フーリエ変換の２回のフーリエ変換を行う必要があり、さらに、これを複数の窓関数で評価しなければならない。その結果、最適な窓関数の決定までに、多大な時間と労力がかかる。
【０００６】
本発明は、このような状況に鑑みてなされたものであり、例えば、窓関数を用いる画像処理に際し、少ない時間および労力で最適な窓関数を決定するのに有利な画像処理方法を提供することを目的とする。
【課題を解決するための手段】
【０００７】
上記課題を解決するために、本発明は、画像を処理する方法であって、画像を窓関数を用いて処理する窓関数処理工程と、窓関数処理工程で処理された空間領域の画像を、フーリエ変換により空間周波数領域の成分に変換するフーリエ変換処理工程と、空間周波数領域における周波数軸の軸上成分と、ローパスフィルタのカットオフ周波数とに基づいて、窓関数を決定する窓関数決定工程と、窓関数決定工程で決定された窓関数を用いて処理された空間周波数領域の成分を、ローパスフィルタを用いて処理するローパスフィルタ処理工程と、ローパスフィルタ処理工程で処理された空間周波数領域の成分を、逆フーリエ変換により空間領域の画像に変換する逆フーリエ変換処理工程とを含むことを特徴とする。
【発明の効果】
【０００８】
本発明によれば、例えば、窓関数を用いる画像処理に際し、少ない時間および労力で最適な窓関数を決定するのに有利な画像処理方法を提供することができる。
【図面の簡単な説明】
【０００９】
【図１】第１実施形態における画像処理工程を示すフローチャートである。
【図２】第１実施形態におけるデータ伸長工程を説明するための図である。
【図３】第１実施形態における窓関数を説明するためのグラフである。
【図４】第１実施形態におけるカットオフ周波数を説明するためのグラフである。
【図５】第１実施形態における振幅スペクトルを示すグラフである。
【図６】第１実施形態におけるローパスフィルタ処理後の画像を示す図である。
【図７】第１実施形態における閾値の設定方法を説明するためのグラフである。
【図８】第２実施形態における画像処理工程を示すフローチャートである。
【図９】第４実施形態における閾値の設定方法を説明するためのグラフである。
【発明を実施するための形態】
【００１０】
以下、本発明を実施するための形態について図面などを参照して説明する。
【００１１】
（第１実施形態）
まず、本発明の第１実施形態に係る画像処理方法について説明する。本実施形態に係る画像処理方法は、例えば、被検体を撮像した画像を処理し、その結果から被検体の状態（形状、構造や収差など）を検査する画像検査装置などの画像処理装置や、独立して画像を処理する情報処理装置などで実行され得る。以下、これらの装置に含まれ得る処理部（制御部）が画像処理方法を実行するものとして説明する。
【００１２】
図１は、本実施形態における画像処理工程の流れを示すフローチャートである。まず、処理部は、画像処理の対象となる計測画像の原画像Ｉ０（ｘ，ｙ）を取得する（ステップＳ１００）。次に、処理部は、データ伸長工程として、ステップＳ１００で取得した原画像Ｉ０（ｘ，ｙ）の端辺（四辺）に画像データを付加し、データ付加後の画像Ｉ（ｘ，ｙ）を求める（ステップＳ２００）。
【００１３】
図２は、データ伸長工程を説明するための図であり、図２（ａ）は、原画像Ｉ０（ｘ，ｙ）の一例を示す図である。図２（ｂ）は、データ付加後の画像Ｉ（ｘ，ｙ）の一例を示す図であり、それに含まれる原画像Ｉ０（ｘ，ｙ）を破線で示している。まず、原画像Ｉ０（ｘ，ｙ）は、図２（ａ）に示すように、ｘ軸方向の画素数がＮ１、ｙ軸方向の画素数がＮ２の２次元画像である。また、原画像Ｉ０（ｘ，ｙ）には、例えば、必要な周波数成分としての図中ｃ１を中心とする低周波成分と、不要な周波数成分としての図中ｃ２、ｃ３を中心とするリング状の高周波成分とが存在する。ステップＳ２００での画像データの付加は、周波数分解能を上げるための処理であり、この処理により、必要な低周波成分と不要な高周波成分とを分離させておくほど、以下のローパスフィルタ処理においてカットオフ周波数を設定するのが容易となる。画像データの付加方法としては、原画像Ｉ０（ｘ，ｙ）の端辺の画像データを、ｘ軸方向にＭ１画素、ｙ軸方向にＭ２画素だけ繰り返して引き延ばす方法がある。この場合、図２（ｂ）に示す画像Ｉ（ｘ，ｙ）のｘ軸方向の画素数Ｎｘと、ｙ軸方向の画素数Ｎｙとは、式（１）で表される。
【００１４】
【数１】
2015108872.tif 000003
【００１５】
また、図２（ｂ）に示す画像Ｉ（ｘ，ｙ）の四隅の各領域Ｇ１、Ｇ２、Ｇ３、Ｇ４は、それぞれＭ１画素×Ｎ１画素の領域を有し、例えば領域Ｇ１は、その値が図２（ａ）における点ｇ１と同一の値を持つ平坦部となっている。
【００１６】
次に、処理部は、閾値取得工程として、以下のステップＳ４１０での窓関数の最適化処理で用いる閾値Ｔｈを取得する（ステップＳ２１０）。なお、本実施形態では、画像処理工程が実行される前に、既に閾値Ｔｈが存在するものと想定する。この場合、閾値Ｔｈは、例えば、処理部に回線を介して接続されている不図示の記憶装置に予め保存されている。また、予め閾値Ｔｈを準備するに際し、閾値Ｔｈを設定する具体的な方法については、ステップＳ４１０での処理に関連するため、後述する。
【００１７】
次に、処理部は、窓関数処理工程として、ステップＳ２００で得られたデータ付加後の画像Ｉ（ｘ，ｙ）に対して窓関数処理を実行し、画像Ｊ（ｘ，ｙ）を求める（ステップＳ３００）。ステップＳ３００における窓関数としては、例えば、式（２）に示す２次元のガウシアンフィルタｈ（ｘ，ｙ）を採用可能である。図３は、本実施形態における各工程での窓関数を説明するためのグラフである。このうち、図３（ａ）は、ガウシアンフィルタｈ（ｘ，ｙ）の一例を示すグラフである。
【００１８】
【数２】
2015108872.tif 000004
【００１９】
ただし、ｎは窓の長さ、σは窓の形状を特定する実数パラメータである。以下、ステップＳ３００での窓関数処理に関し、式（２）に示すガウシアンフィルタｈ（ｘ，ｙ）から窓関数処理済みの画像Ｊ（ｘ，ｙ）を求める具体的な方法について説明する。ここでは、ガウシアンフィルタｈ（ｘ，ｙ）に対して、式（３）に示すように、ｘ軸方向、ｙ軸方向に投影し、ｘ軸に投影した信号ｈｙ(ｘ)と、ｙ軸に投影した信号ｈｘ(ｙ)とそれぞれを求める。図３（ｂ）は、ｘ軸に投影したｈｙ(ｘ)の一例を示すグラフである。
【００２０】
【数３】
2015108872.tif 000005
【００２１】
ただし、信号ｈｙ（ｘ）、ｈｘ（ｙ）の長さは、それぞれｎである。また、ステップＳ２００でのデータ付加後の画像Ｉ（ｘ，ｙ）のｘ軸方向、ｙ軸方向の各画素数は、式（１）に示すように、Ｎｘ、Ｎｙである。そこで、処理部は、ｎ＜Ｎｘ、ｎ＜Ｎｙの関係を考慮し、長さｎをＮｘ、Ｎｙにともに拡張した後、自己相関関数β（ｘ）、β（ｙ）をそれぞれ求める。信号ｈｙ（ｘ）、ｈｘ（ｙ）を拡張する方法は、例えば、式（４）に示すように、端の値をＮｘ、Ｎｙまで拡張すればよい。または、式（４）に換えて、端の値は十分ゼロに収束している仮定して、ゼロ値を拡張してもよい。
【００２２】
【数４】
2015108872.tif 000006
【００２３】
式（３）、（４）により拡張された信号ｈｙ（ｘ）、ｈｘ（ｙ）に対する自己相関関数β（ｘ）、β（ｙ）は、式（５）で表される。
【００２４】
【数５】
2015108872.tif 000007
【００２５】
ただし、ＦＴはフーリエ変換、ＩＦＴは逆フーリエ変換、ｃｏｎｊは複素共役である。本実施形態では、相関計算は、直接空間軸で行われるのではなく、式（５）に示すようにフーリエ変換を利用して空間周波数領域としてから行われる。特に、信号ｈｘ（ｙ）が長く、すべての位置ずれにおける相関値を求めるときには、直接計算よりも、このようにフーリエ変換を利用した相関値計算の方が、処理速度の面で有利である。図３（ｃ）は、拡張したｈｙ(ｘ)に対する自己相関関数β（ｘ）の一例を示すグラフである。ここで、自己相関関数β（ｘ）は、ピークで正規化し、そのピークが原点となるようにシフト処理されている。次に、処理部は、自己相関関数β（ｘ）、β（ｙ）を用いて、式（６）より、窓関数Ｗ（ｘ，ｙ）を求める。
【００２６】
【数６】
2015108872.tif 000008
【００２７】
図３（ｄ）は、窓関数Ｗ（ｘ，ｙ）の一例を示すグラフである。本実施形態における窓関数Ｗ（ｘ，ｙ）は、図３（ｄ）に示すように、画像の中心部では重みが１で、周辺部にいくほど重みが滑らかにゼロになる特徴がある。そして、最終的にステップＳ３００で得られる窓関数処理済みの画像Ｊ（ｘ，ｙ）は、式（７）に示すような窓関数Ｗ（ｘ，ｙ）に関する重みつき和により求められる。
【００２８】
【数７】
2015108872.tif 000009
【００２９】
ただし、Ｉｃ（ｘ，ｙ）は、ステップＳ２００で得られたデータ付加後の画像Ｉ（ｘ，ｙ）とガウシアンフィルタｈ（ｘ，ｙ）とのコンボリューションにより求められる画像（ボケ画像）を表す。なお、画像Ｉｃ（ｘ，ｙ）を求めるときには巡回型畳み込みを行うため、画像Ｉｃ（ｘ，ｙ）の四隅の部分には、連続性を持たせている。このように、画像Ｊ（ｘ，ｙ）は、画像の中心部では、画像Ｉ（ｘ，ｙ）に一致し、周辺部では、ボケ画像Ｉｃ（ｘ，ｙ）に一致すると同時に、上記のとおり四隅の部分では連続性を持つ。
【００３０】
次に、処理部は、フーリエ変換処理工程として、ステップＳ３００で得られた窓関数処理済みの画像Ｊ（ｘ，ｙ）に対してフーリエ変換を実行し、空間周波数領域の成分に変換する（ステップＳ４００）。
【００３１】
次に、処理部は、窓関数決定工程として、空間周波数軸の軸上成分に着目し、窓関数が最適であるかどうかを判断する（ステップＳ４１０）。なお、窓関数が最適であるかどうかを判断する具体的な方法については、このときに用いる閾値Ｔｈと合わせ、後述する。ここで、処理部は、窓関数が最適でないと判断したら（ＮＯ）、次に式（２）中の窓関数のパラメータσの値を変更し（ステップＳ４２０）、その後、ステップＳ３００に戻る。すなわち、窓関数が最適でない場合には、窓関数のパラメータσの値を変更してガウシアンフィルタｈ（ｘ，ｙ）を再生成し、窓関数処理済みの画像Ｊ（ｘ，ｙ）を繰り返し求めることになる。一方、ステップＳ４１０にて、処理部は、窓関数が最適であると判断したら（ＹＥＳ）、ステップＳ５００に移行する。
【００３２】
次に、処理部は、ローパスフィルタ処理工程として、空間周波数領域において、カットオフ周波数ｆｃに基づいてローパスフィルタ処理を実行し、不要な高周波成分を除去する（ステップＳ５００）。以下、カットオフ周波数ｆｃの決定方法について具体的に説明する。先のステップＳ４００で、窓関数処理済みの画像Ｊ（ｘ，ｙ）に対してフーリエ変換したときの振幅スペクトルをＰ（ｕ，ｖ）とする。図４は、カットオフ周波数ｆｃを説明するためのグラフである。このうち、図４（ａ）は、振幅スペクトルＰ（ｕ，ｖ）を模式的に示すグラフである。なお、振幅スペクトルに換えて、振幅の二乗であるパワースペクトルを用いてもよい。図４（ａ）において、画像Ｊ（ｘ，ｙ）のｘ軸に対応する空間周波数軸がｕ軸であり、ｙ軸の対応する空間周波数軸がｖ軸に対応する。また、画像Ｊ（ｘ，ｙ）のフーリエ変換後の振幅スペクトルＰ（ｕ、ｖ）を極座標形式で表し、Ｐ（ｒ、θ）とする。まず、空間周波数領域において、原点からの半径ｒで所定の幅±Δｒのリング状の領域に着目し、処理部は、式（８）より、ｒ±Δｒ内の周波数成分の総和Ｐ（ｒ）を求める。
【００３３】
【数８】
2015108872.tif 000010
【００３４】
図４（ｂ）は、半径ｒに対する周波数成分の総和Ｐ（ｒ）の一例を示すグラフである。不要な高周波成分は、図４（ａ）に示すようなリング状の特性を持つ。そこで、この特性を利用し、処理部は、図４（ｂ）に示すような極大値を持つ周波数での半径ｒ０を求め、このときの周波数を不要な高周波成分と特定する。次に、処理部は、原点から半径ｒ０までの間で極小値をもつ周波数をローパスフィルタのカットオフ周波数ｆｃと決定する。なお、原点から半径ｒ０までの間の極小値付近は、Ｓ／Ｎ比が小さい領域でもあるため、極小値が複数存在し、一意にカットオフ周波数ｆｃを決定することができない場合もあり得る。この場合、制御部は、例えば、原点と半径ｒ０との距離の１／２に相当する周波数をローパスフィルタのカットオフ周波数ｆｃと決定してもよい。
【００３５】
次に、処理部は、逆フーリエ変換処理工程として、逆フーリエ変換処理を実行し、空間領域に戻された画像Ｉ１（ｘ，ｙ）を求める（ステップＳ６００）。なお、画像Ｉ１（ｘ，ｙ）では、ｘ軸方向の画素数がＮｘ、ｙ軸方向の画素数がＮｙである。次に、処理部は、データ切り出し工程として、ステップＳ２００で付加したデータを切り出す（ステップＳ７００）。具体的には、処理部は、画像Ｉ１（ｘ，ｙ）に対して、ｘ軸方向の両側からＭ１、ｙ軸方向の両側からＭ２だけデータを切り出すことで、最終的に、原画像Ｉ０（ｘ，ｙ）と同じ画素数の画像Ｉ'（ｘ，ｙ）を求める。そして、処理部は、ステップＳ７００の後、画像処理工程を終了する。
【００３６】
次に、ステップＳ４１０の窓関数決定工程において、窓関数のパラメータσが最適なものであるかを判断する（最適なパラメータσを決定する）方法について具体的に説明する。図５は、パラメータσが、ある値のときに、先のステップＳ４００で、窓関数処理済みの画像Ｊ（ｘ，ｙ）に対してフーリエ変換したときの振幅スペクトルＰ（ｕ，ｖ）に関するグラフである。このうち、図５（ａ）は、振幅スペクトルＰ（ｕ，ｖ）の一例を３次元状に示すグラフである。これに対して、図５（ｂ）は、パラメータσが小さい値のときの、一方、図５（ｃ）は、パラメータσが大きい値のときの各振幅スペクトルＰ（ｕ，ｖ）についての空間周波数軸のｕ、ｖ軸上の値を示すグラフである。図５（ｂ）および図５（ｃ）を参照すると、パラメータσが小さい値ではｕ、ｖ軸上の値は大きいが、パラメータσが大きい値ではｕ、ｖ軸上の値が小さくなる。
【００３７】
図６は、ステップＳ５００でのローパスフィルタ処理の後の画像Ｉ１（ｘ，ｙ）に関する図である。このうち、図６（ａ)は、パラメータσが、ある値のときに、ステップＳ５００でローパスフィルタ処理後、ステップＳ６００で逆フーリエ変換して空間領域に戻した画像Ｉ１（ｘ，ｙ）を示す図である。なお、画像Ｉ１（ｘ，ｙ）は、ステップＳ２００での画像データを付加したままであり、ｘ、ｙ軸方向の各画素数は、Ｎｘ，Ｎｙである。リンギングノイズは、図６（ａ）に示すような２次元画像からでも、わずかに確認されるが、本実施形態では、処理部は、図６（ａ）中の一点鎖線Ａ−Ａ’に対応する１次元信号を抜き出し、リンギングノイズを判断する。図６（ｂ）は、パラメータσが小さい値のときの、一方、図６（ｃ）は、パラメータσが大きい値のときの画像Ｉ１（ｘ，ｙ）に対するＡ−Ａ’に対応する１次元信号を抽出したグラフである。
【００３８】
ここで、図５（ｂ）を参照すると、空間周波数軸のｕ、ｖ軸上成分が所定の閾値Ｔｈを横切る周波数ｆ１が、カットオフ周波数ｆｃよりも大きい。この場合、図６（ｂ）に示すように、リンギングノイズが残存する。これに対して、図５（ｃ）を参照すると、ｕ、ｖ軸上の成分が閾値Ｔｈを横切る位置の周波数ｆ２が、カットオフ周波数ｆｃよりも小さく、この場合には、図６（ｃ）に示すように、リンギングノイズが抑制されていることがわかる。すなわち、処理部は、最適なパラメータσを決定する際には、空間周波数領域における周波数軸の軸上成分が閾値Ｔｈを横切る位置の周波数が、カットオフ周波数ｆｃよりも小さくなるときのパラメータσを選択すればよい。一方、図５（ｂ）を参照すると、カットオフ周波数ｆｃにおける空間周波数軸の軸上成分Ｔｂは、閾値Ｔｈよりも大きいが、図５（ｃ）を参照すると、カットオフ周波数ｆｃにおける空間周波数軸の軸上成分Ｔｃは、閾値Ｔｈよりも小さい。したがって、上記の決定方法に換えて、処理部は、最適なパラメータσを決定する際には、空間周波数領域において、カットオフ周波数ｆｃでの周波数軸の軸上成分が、閾値Ｔｈよりも小さくなるときのパラメータσを選択するものとしてもよい。
【００３９】
次に、閾値Ｔｈの設定方法について具体的に説明する。この設定方法で設定された閾値Ｔｈは、図１に示す画像処理工程のステップＳ２１０にて処理部に取得される。図７は、窓関数のパラメータσがσ１、σ２、σ３のときに、一例として、それぞれのフーリエ変換の振幅スペクトルＰ（ｕ）の空間周波数軸のｕ軸上成分を示すグラフである。なお、図７では、各パラメータの値に対応し、図６（ｂ）および図６（ｃ）と同様、逆フーリエ変換して空間領域に戻した画像Ｉ１（ｘ，ｙ）におけるＡ−Ａ’に対応する１次元信号を抽出したグラフも併記している。図７に示すように、パラメータσ＝σ１のときは、空間画像のリンギングノイズが十分抑制されている。また、パラメータσ＝σ３のときは、リンギングノイズが残存している。そして、パラメータσ＝σ２のときは、リンギングノイズが許容範囲内にあるとした場合に、振幅スペクトルＰ（ｕ）がσ２より下方に位置していることが必要である。そこで、処理部は、パラメータσ＝σ２の振幅スペクトルＰ（ｕ）の空間周波数軸のｕ軸上成分で、ローパスフィルタのカットオフ周波数ｆｃにおける値Ｔｈを閾値に設定する。ここで、リンギングノイズが許容範囲内であるか否かを判断する基準としては、例えば、図７に示す画像の平坦部における振動成分の振幅ΔＰが、所定の値を超えているかどうかで判断すればよい。
【００４０】
ローパスフィルタ処理では、ローパスフィルタのカットオフ周波数ｆｃにおける空間周波数軸の軸上成分が小さいほど、失われた高周波成分に起因するリンギングノイズが小さくなる傾向がある。そして、窓関数のパラメータσの値を大きくすれば、カットオフ周波数ｆｃにおける空間周波数軸の軸上成分が小さくなる。そこで、本実施形態では、上記のように、窓関数のパラメータσの値を変更させながらリンギングノイズを抑制させる際に空間周波数軸の軸上成分に着目することで、窓関数を最適化させる。したがって、本実施形態では、リンギングノイズを抑制する最適な窓関数のパラメータσの値を、従来の空間領域での画像で判断しなくとも、空間周波数領域で（空間周波数軸の軸上だけで）求めることができる。すなわち、従来、ある窓関数が最適であるかを判断するに際し、フーリエ変換および逆フーリエ変換の２回のフーリエ変換を行い、かつその窓関数に対して評価することを繰り返していた。これに対して、本実施形態に係る画像処理方法では、ある窓関数が最適であるかを判断するに際しては、１回のフーリエ変換（図１のステップＳ４００）後に評価すること（図１のステップＳ４１０）を繰り返すだけでよく、逆フーリエ変換を要しない。これにより、最適な窓関数の決定までに要する時間と労力とを抑えることができる。
【００４１】
以上のように、本実施形態によれば、窓関数を用いる画像処理に際し、少ない時間および労力で最適な窓関数を決定するのに有利な画像処理方法を提供することができる。
【００４２】
（第２実施形態）
次に、本発明の第２実施形態に係る画像処理方法について説明する。上記の第１実施形態では、図１に示す画像処理工程内のステップＳ４２０において、窓関数のパラメータとして、ガウシアンフィルタのパラメータσを用いた。これに対して、本実施形態に係る画像処理方法の特徴は、ステップＳ４２０における窓関数のパラメータとして、カイザー窓のパラメータαを用いる点にある。ここで、カイザー窓γ（ｘ）とは、式（９）で定義される窓関数である。
【００４３】
【数９】
2015108872.tif 000011
【００４４】
ただし、Ｉ０は第１種の０次の変形ベッセル関数であり、αは窓の形状を特定する実数パラメータであり、整数ｎは窓の長さである。カイザー窓γ（ｘ）は、パラメータαの値を大きくしていくと、広がりが狭くなり、一方、パラメータαの値を小さくしていくと、広がりが大きくなる性質がある。このカイザー窓γ（ｘ）を、第１の実施形態における式（４）に示すｘ軸に投影した信号ｈｙ（ｘ）、ｙ軸に投影した信号ｈｘ（ｙ）に適用する。すなわち、式（３）（図３（ｂ））を式（９）に示すカイザー窓に変更し、式（５）および式（６）から窓関数処理済みの画像Ｊ（ｘ，ｙ）を求める。
【００４５】
本実施形態によれば、第１実施形態と同様の効果を奏し、特に、ガウシアンフィルタよりもカイザー窓の方が窓の形状のダイナミックレンジが広いので、窓関数のパラメータをより広範に振ることができる。
【００４６】
（第３実施形態）
次に、本発明の第３実施形態に係る画像処理方法について説明する。上記の第１実施形態では、図１に示す画像処理工程内のステップＳ４１０における窓関数の最適化処理で用いられる閾値Ｔｈは、画像処理工程の開始前に予め存在しており、処理部は、閾値取得工程であるステップＳ２１０において、その閾値Ｔｈを取得する。これに対して、本実施形態に係る画像処理方法の特徴は、予め閾値Ｔｈが存在していないものと想定し、画像処理工程内に、新たに閾値Ｔｈを設定する工程を追加している点にある。
【００４７】
図８は、本実施形態における画像処理工程の流れを示すフローチャートである。図８に示す画像処理工程は、第１実施形態における図１に示す画像処理工程と比較すると、ステップＳ２１０での閾値取得工程に換えて、予め閾値Ｔｈを設定するかどうかを判断する閾値設定判断工程など（ステップＳ８００以降）を含む。なお、その他のステップＳ１００からステップＳ７００までの各工程は、図１と同一であるので、説明を省略する。この場合、処理部は、ステップＳ２００でのデータ伸長工程の後で、かつ、ステップＳ３００での窓関数処理工程の前に、閾値Ｔｈを設定するかどうかを判断する（ステップＳ８００：閾値設定判断工程）。ここで、処理部は、例えば、計測対象の画像を今回初めて処理するときに、閾値Ｔｈを設定する（ＹＥＳ）と判断し、その最初の計測画像（少なくとも１つの画像）を用いて設定する。
【００４８】
処理部は、ステップＳ８００で閾値Ｔｈを設定すると判断した場合、まず、ステップＳ３００と同様の窓関数処理を実行する（ステップＳ８１０）。次に、処理部は、ステップＳ４００と同様に、窓関数処理後の画像Ｊ（ｘ，ｙ）にフーリエ変換を実行し、空間周波数領域での振幅スペクトルＰ（ｕ，ｖ）に変換する（ステップＳ８２０）。次に、処理部は、ステップＳ５００〜Ｓ７００の流れと同様に、振幅スペクトルＰ（ｕ、ｖ）に対してローパスフィルタ処理（ステップＳ８３０）、逆フーリエ変換（ステップＳ８４０）、およびデータの切り出し（ステップＳ８５０）を順次実行する。次に、処理部は、最終的に得られる空間領域での画像Ｉ’（ｘ，ｙ）を評価し、リンギングノイズが残存していないかどうかを判断する（ステップＳ８６０）。具体的には、処理部は、画像Ｉ’（ｘ，ｙ）の平坦部やエッジ近傍などに着目して、直接的にリンギングノイズが残存していないか確認することで評価を行う。ここで、処理部は、リンギングノイズが残存していない（ＮＯ）、すなわち、このときの窓関数は最適であると判断したら、以下のステップＳ８８０に移行する。そして、処理部は、ステップＳ８２０で得られた空間周波数領域における振幅スペクトルＰ（ｕ，ｖ）に対して、ローパスフィルタのカットオフ周波数ｆｃにおける空間周波数軸のｕ、ｖ軸上の値を閾値Ｔｈとして設定する（ステップＳ８８０）。一方、ステップＳ８６０にて、処理部は、リンギングノイズが残存している（ＹＥＳ）と判断したら、ステップＳ８７０で窓関数のパラメータ（例えば式（２）におけるパラメータσ）を変更した後、ステップＳ８１０に戻り、窓関数処理以下の工程を繰り返す。
【００４９】
そして、ステップＳ８８０にて閾値Ｔｈが設定された後は、処理部は、ステップＳ１００に戻り、次に２回目以降の計測画像を取得して、以後、ステップＳ８００では閾値Ｔｈを設定しない（ＮＯ）と判断し、ステップＳ３００に移行する。
【００５０】
このように、画像処理工程の前に予め閾値Ｔｈが存在せずとも、最初の計測画像を用いて１回だけ閾値Ｔｈを設定すれば、それ以降の計測画像に対しては、その閾値Ｔｈに基づいて窓関数の最適化を行うことができる。結果的に、本実施形態によれば、第１実施形態と同様の効果を奏する。
【００５１】
なお、図８に示す画像処理工程では、ステップＳ８５０としてデータ切り出し工程を含むものとしているが、このデータ切り出し工程を含まないものとしてもよい。この場合、処理部は、ステップＳ８６０でのリンギングノイズが残存しているかどうかの判断を、データ切り出し後の空間領域での画像Ｉ’（ｘ，ｙ）に対してではなく、データ切り出し前の空間領域での画像Ｉ１（ｘ，ｙ）に対して実行する。データ切り出し前の画像Ｉ１（ｘ，ｙ）は、図２（ｂ）に示す画像Ｉ（ｘ，ｙ）と同じ大きさ（Ｎｘ画素×Ｎｙ画素）を有する。また、図２（ｂ）の四隅（破線で示される画像Ｉ（ｘ，ｙ）と同じ大きさの領域の周辺部）に示される領域Ｇ１、Ｇ２、Ｇ３、Ｇ４は、上記のとおり、それぞれＭ１画素×Ｎ１画素の大きさを有する平坦部である。すなわち、画像Ｉ１（ｘ，ｙ）においても、各領域Ｇ１〜Ｇ４に相当する四隅の領域に、もともと平坦部が存在している。そこで、これらの四隅の領域の平坦部に着目すれば、処理部は、ステップＳ８６０での判断をより効果的に実行できるので、結果的に、閾値Ｔｈをより迅速に設定することができる。
【００５２】
（第４実施形態）
次に、本発明の第４実施形態に係る画像処理方法について説明する。上記の第１実施形態では、予め設定しておく閾値Ｔｈの設定方法として、画像Ｉ１（ｘ，ｙ）に基づいて、リンギングノイズが許容範囲内にあるときの空間周波数軸の軸上成分と、ローパスフィルタのカットオフ周波数ｆｃとの交点から求める方法を例示した。これに対して、本実施形態に係る画像処理方法の特徴は、空間周波数領域での軸上成分の振幅スペクトルＰ（ｕ，ｖ）の値に基づいて、予め閾値Ｔｈを設定する点にある。
【００５３】
図９は、窓関数のパラメータσを変化させたときの振幅スペクトルＰ（ｕ，ｖ）における空間周波数軸のｕ、ｖ軸上の値を示すグラフである。このうち、図９（ａ）は、パラメータσが小さい値のときで、リンギングノイズが残存している状態である。一方、図９（ｂ）は、図９（ａ）よりパラメータσが大きい値のときで、リンギングノイズが許容範囲内にある状態である。また、図９（ｂ）において、カットオフ周波数ｆｃとの交点（周波数領域Ｄの基準点）の振幅値である閾値を図中Ｔｈで示し、図９（ａ）にも図示した。
【００５４】
まず、図９（ａ）において、Ａで示される値は、振幅スペクトルの原点における値で、画像のＤＣ成分を表し、画像の平均値や中間値で引いたりする処理により、容易に変動し得る。そのため、振幅スペクトルの原点のように上記処理によって左右される値に基づいて閾値Ｔｈを設定するのは望ましくない。そこで、原点を除く振幅スペクトルの最大振幅（図９（ａ）中のＢで示される位置）で正規化し、さらに、縦軸の表示を、式（１０）に示すようなデシベル表示とする。
【００５５】
【数１０】
2015108872.tif 000012
【００５６】
このとき、閾値Ｔｈは、−４６ｄＢにあり、原点を除く振幅スペクトルの最大振幅を１とすると、０．５％に相当する。ここで、閾値が図９（ｂ）中のＴｈｕに示すようにＴｈよりも大きいと、リンギングノイズが残存する程パラメータσが小さいときでも、軸上成分が閾値を横切る位置の空間周波数が、カットオフ周波数ｆｃよりも小さくなる条件を満たしてしまう場合がある。そこで、図９（ｂ）中に示すＴｈを閾値の上限とする。
【００５７】
一方、カットオフ周波数ｆｃ近傍の周波数領域Ｄでは、振幅スペクトルは、パラメータσの変化に伴い大きく変化する。しかしながら、カットオフ周波数ｆｃよりも大きな周波数領域Ｅでは、振幅スペクトルは、パラメータσを変化させても大きく変化せず、振幅自体も小さいため、Ｓ／Ｎ比が小さい。そのため、閾値を図９（ｂ）中のＴｈｄに示すようにＳ／Ｎ比が小さい周波数領域Ｅと同等の振幅レベルまで小さくすると、図中Ｆに示すように軸上成分が閾値を横切る位置の空間周波数がカットオフ周波数ｆｃよりも大きくなる可能性がある。つまり、図９（ｂ）のように、リンギングノイズが許容範囲内にあるにも関わらず、条件を満たさない場合があり望ましくない。そこで、閾値の下限は、カットオフ周波数ｆｃより大きな周波数領域Ｅの最大振幅Ｆ（−５０ｄＢ）よりも大きくする。なお、このときの値は、原点を除く振幅スペクトルの最大振幅を１とすると、０．３％に相当する。
【００５８】
このように、閾値Ｔｈを空間周波数領域での振幅スペクトルの値に基づいて設定する場合には、その閾値を振幅スペクトルの原点を除く最大振幅の０．３％から０．５％までの間にあるものとすればよい。結果的に、本実施形態によっても、第１実施形態と同様の効果を奏する。
【００５９】
以上、本発明の好ましい実施形態について説明したが、本発明は、これらの実施形態に限定されず、その要旨の範囲内で種々の変形および変更が可能である。
【符号の説明】
【００６０】
Ｓ３００  窓関数処理工程
Ｓ４００  フーリエ変換処理工程
Ｓ４１０  窓関数決定工程
Ｓ５００  ローパスフィルタ処理工程
Ｓ６００  逆フーリエ変換処理工程
【図１】
2015108872.tif 000013
【図４】
2015108872.tif 000014
【図７】
2015108872.tif 000015
【図８】
2015108872.tif 000016
【図９】
2015108872.tif 000017
【図２】
2015108872.tif 000018
【図３】
2015108872.tif 000019
【図５】
2015108872.tif 000020
【図６】
2015108872.tif 000021