(57)【特許請求の範囲】
【請求項１】
入力データを所定データサイズの複数の分割データに分割するデータの分割手段と、
前記所定データサイズに基づき、フィルタを複数の分割フィルタに分割するフィルタの分割手段と、
前記複数の分割フィルタごとにフィルタ係数の近似多項式を作成する作成手段と、
前記複数の分割データ、前記複数の分割フィルタおよび前記近似多項式を用いて、前記入力データに前記フィルタによるフィルタ処理を施す処理手段とを有することを特徴とするデータ処理装置。
【請求項２】
前記処理手段は、第一の分割データを選択する手段と、
前記第一の分割データを中心とする分割データの参照範囲を決定する手段と、
前記第一の分割データを前記フィルタ処理する際に参照すべき第二の分割データを前記参照範囲から選択する手段と、
前記第一の分割データに含まれる注目データの前記フィルタ処理において、前記第一の分割データと前記第二の分割データの位置関係に基づき前記近似多項式の係数を計算する手段と、
前記係数および前記近似多項式を用いて前記第二の分割データから前記注目データへ寄与するデータを計算する手段とを有することを特徴とする請求項1に記載されたデータ処理装置。
【請求項３】
前記フィルタの分割手段は、前記フィルタのフィルタサイズを拡大して、前記フィルタを奇数個の前記所定データサイズの分割フィルタに分割することを特徴とする請求項1または請求項2に記載されたデータ処理装置。
【請求項４】
画像データが表す画像を所定サイズの複数の分割領域に分割する画像の分割手段と、
前記所定サイズに基づき、二次元フィルタを複数の分割フィルタに分割するフィルタの分割手段と、
前記複数の分割フィルタごとにフィルタ係数の近似多項式を作成する作成手段と、
前記複数の分割領域、前記複数の分割フィルタおよび前記近似多項式を用いて、前記画像データに前記二次元フィルタによるフィルタ処理を施す処理手段とを有することを特徴とする画像処理装置。
【請求項５】
前記処理手段は、第一の分割領域を選択する手段と、
前記第一の分割領域を中心とする分割領域の参照範囲を決定する手段と、
前記第一の分割領域を前記フィルタ処理する際に参照すべき第二の分割領域を前記参照範囲から選択する手段と、
前記第一の分割領域に含まれる注目画素の前記フィルタ処理において、前記第一の分割領域と前記第二の分割領域の位置関係に基づき前記近似多項式の係数を計算する手段と、
前記係数および前記近似多項式を用いて前記第二の分割領域から前記注目画素へ寄与するデータを計算する手段とを有することを特徴とする請求項4に記載された画像処理装置。
【請求項６】
データの分割手段、フィルタの分割手段、作成手段、処理手段を有するデータ処理装置のデータ処理方法であって、
前記データの分割手段が、入力データを所定データサイズの複数の分割データに分割し、
前記フィルタの分割手段が、前記所定データサイズに基づき、フィルタを複数の分割フィルタに分割し、
前記作成手段が、前記複数の分割フィルタごとにフィルタ係数の近似多項式を作成し、
前記処理手段が、前記複数の分割データ、前記複数の分割フィルタおよび前記近似多項式を用いて、前記入力データに前記フィルタによるフィルタ処理を施すことを特徴とするデータ処理方法。
【請求項７】
データの分割手段、フィルタの分割手段、作成手段、処理手段を有する画像処理装置の画像処理方法であって、
前記データの分割手段が、画像データが表す画像を所定サイズの複数の分割領域に分割し、
前記フィルタの分割手段が、前記所定データサイズに基づき、二次元フィルタを複数の分割フィルタに分割し、
前記作成手段が、前記複数の分割フィルタごとにフィルタ係数の近似多項式を作成し、
前記処理手段が、前記複数の分割領域、前記複数の分割フィルタおよび前記近似多項式を用いて、前記画像データに前記二次元フィルタによるフィルタ処理を施すことを特徴とする画像処理方法。
【請求項８】
コンピュータを請求項1から請求項3の何れか一項に記載されたデータ処理装置の各手段として機能させることを特徴とするプログラム。
【請求項９】
コンピュータを請求項4または請求項5に記載された画像処理装置の各手段として機能させることを特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、ディジタルデータのフィルタ処理に関する。
【背景技術】
【０００２】
フィルタサイズが大きいフィルタ処理の高速化手法として、例えば、フィルタ処理を行う周波数域を制限する方法が知られている（例えば、特許文献1）。特許文献1は、二焦点レンズのピント暈けを補正する技術を開示する。つまり、暈けを表す関数の逆関数と撮影画像を畳込演算する際に、ピント暈けが低域成分に存在することを考慮して、帯域分割により撮影画像の低域成分を取り出して畳込演算の演算量を削減する。
【０００３】
また、フィルタによる畳込演算を、フィルタ処理する注目画素と、その周辺画素の相互作用と考えると、空間上に存在する多数の粒子の間において、何らかの力の相互作用が生じる場を解く多体問題(N-body problem)と見做すことができる。多体問題の高速化手法として、高速多重極法(fast multipole method)が知られている（例えば、非特許文献1）。高速多重極法は、重力またはクーロン力が相互作用する場において、遠く離れた物体間の作用力は小さいことに着目し、式展開とツリー構造による領域分割により、離れた粒子からの寄与分をまとめて効率よく計算する手法である。
【０００４】
しかし、特許文献1が開示する技術は、ある特定の帯域のみをフィルタ処理することはできるが、全周波数域をフィルタ処理する場合に適用することができない。
【０００５】
また、非特許文献1が開示する技術は、重力やクーロン力のような二点間の距離の二乗に反比例するような力に対して式変形を施して計算量を削減するが、相互作用力がより一般的な任意の関数で表される場合には適用することができない。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開2007-72558公報
【非特許文献】
【０００７】
【非特許文献１】L. Greengard and V. Rohklin「A Fast Algorithm for Particle Simulations」JOURNAL OF COMPUTATIONAL PHYSICS 73, 315-348(1987)
【発明の概要】
【発明が解決しようとする課題】
【０００８】
本発明は、ディジタルデータのフィルタ処理における計算量を低減して、フィルタ処理を高速化することを目的とする。
【課題を解決するための手段】
【０００９】
本発明は、前記の目的を達成する一手段として、以下の構成を備える。
【００１０】
本発明にかかるデータ処理は、入力データを所定データサイズの複数の分割データに分割し、前記所定データサイズに基づき、フィルタを複数の分割フィルタに分割し、前記複数の分割フィルタごとにフィルタ係数の近似多項式を作成し、複数の分割データ、前記複数の分割フィルタおよび前記近似多項式を用いて、前記入力データに前記フィルタによるフィルタ処理を施すことを特徴とする。
【発明の効果】
【００１１】
本発明によれば、ディジタルデータのフィルタ処理における計算量を低減して、フィルタ処理を高速化することができる。
【図面の簡単な説明】
【００１２】
【図１】実施例の画像処理装置の構成例を説明するブロック図。
【図２Ａ】画像処理装置が実行するフィルタ処理を説明するフローチャート。
【図２Ｂ】フィルタ処理の詳細を説明するフローチャート。
【図３】フィルタサイズの分割を説明する図。
【図４】畳込演算における画像データAとフィルタB'の関係を説明する図。
【図５】分割領域Rdi、Rdj周辺の拡大図を示す図。
【図６】フィルタ範囲の拡大を説明する図。
【図７】処理結果の確認に使用した画像Aの一例を示す図。
【図８】処理結果の確認に使用した画像Aの一例を示す図。
【図９】画像Aに対し、フィルタBによる処理を行って取得した画像Cと、フィルタB'による処理を行って取得した画像Dを示す図。
【図１０】フィルタBを用いる処理とフィルタB'を用いる処理の比較結果を説明する図。
【図１１】分割領域のサイズを20×20画素として、フィルタ処理した後の画像および空間周波数分布を示す図。
【図１２Ａ】一次元フィルタ処理を説明するフローチャート。
【図１２Ｂ】フィルタ処理の詳細を説明するフローチャート。
【発明を実施するための形態】
【００１３】
以下、本発明にかかる実施例のデータ処理および画像処理を図面を参照して詳細に説明する。
【実施例１】
【００１４】
［装置の構成］
図1のブロック図により実施例の画像処理装置の構成例を説明する。
【００１５】
CPU101は、RAM103をワークメモリとして、ROM102やハードディスクドライブ(HDD)108に格納されたプログラムを実行し、システムバス104を介して後述する各構成を制御する。これにより、後述する様々な処理が実行される。
【００１６】
汎用インタフェイス(I/F)105は、キーボードやマウス、ディジタルカメラ、スキャナ、メモリカードリーダ、ハードディスクドライブなど、様々なデバイス106を接続する、例えばUSBやIEEE1394などのシリアルバスインタフェイスである。CPU101は、汎用I/F105を介して、デバイス106から各種データを読み込み、各種データをデバイス106に書き込むことが可能である。
【００１７】
HDDインタフェイス(I/F)107は、HDD108や光ディスクドライブなどの二次記憶装置を接続する、例えばシリアルATA(SATA)などのインタフェイスである。CPU101は、HDD I/F107を介して、HDD108からデータを読み出し、HDD108にデータを書き込むことが可能である。さらに、CPU101は、HDD108に格納されたデータをRAM103に展開し、同様に、RAM103に展開されたデータをHDD108に保存することが可能である。そしてCPU101は、RAM103に展開したデータをプログラムと見做して実行することができる。
【００１８】
ビデオカード(VC)109は、モニタ110用のインタフェイスである。CPU101は、VC109を介して、後述するユーザインタフェイス(UI)、フィルタ処理すべき画像、フィルタ処理後の画像を含む各種画像をモニタ110に表示する。
【００１９】
ユーザは、デバイス106として汎用I/F105に接続されたキーボードやマウスにより、モニタ110に表示されたUIを操作して、例えば後述するフィルタ処理などの実行をCPU101に指示する。フィルタ処理の実行を指示されたCPU101は、ROM102またはHDD108に格納されたフィルタ処理のプログラムを実行する。そして、CPU101は、UIを介してユーザが指定する画像データをHDD108やデバイス106から読み込み、読み込んだ画像データにフィルタ処理を施す。さらに、フィルタ処理後の画像データを、UIを介してユーザが指定する、HDD108やデバイス106の領域に格納したり、フィルタ処理後の画像データが表す画像をモニタ110に表示したりする。
【００２０】
［フィルタ処理］
図2Aのフローチャートにより画像処理装置が実行するフィルタ処理を説明する。
【００２１】
CPU101は、ユーザが指定する画像データAを入力し(S101)、入力した画像データAが表す画像を横Nx画素×縦Ny画素（所定サイズ）の複数の領域に分割する(S102)。なお、画像を横方向（X方向）に分割する場合、画像の例えば右端において幅がNx画素に満たない領域が発生する可能性があるが、そのような領域も分割領域として扱う。同様に、画像を縦方向（Y方向）に分割する場合、画像の例えば下端において高さがNy画素に満たない領域が発生する可能性があるが、そのような領域も分割領域として扱う。また、以下の説明では、分割領域の数をMとする。
【００２２】
次に、CPU101は、フィルタ処理用のフィルタBを表すデータを例えばHDD108から入力する(S103)。フィルタBのフィルタサイズを横Nfx画素×縦Nfy画素とする。フィルタBは、注目画素からの距離の関数として表現されてもよいし、フィルタの各成分（フィルタ係数）を表すデータでもよい。ただし、後述するように、フィルタBのフィルタサイズを分割領域のサイズの自然数倍に拡大するため、フィルタ処理の裾に相当する部分の成分（フィルタ係数）がほぼ零であることが好ましい。
【００２３】
次に、CPU101は、分割領域のサイズNx×Nyに対応させて、フィルタBを分割する(S104)。その際、CPU101は、分割後のフィルタが縦横ともに奇数個の分割フィルタ（サイズNx×Ny）から構成されるように、フィルタサイズを拡大する。図3によりフィルタサイズの分割を説明する。つまり、分割フィルタの数を縦Lx×横Lyとすると、Lx、Lyが下式を満たすようにフィルタサイズを拡大する。なお、以下では、分割フィルタの数をLx×Ly=LLと表す。
Nx×Lx ≧ Nfx
Ny×Ly ≧ Nfy                    …(1)
ここで、Lx、Lyはそれぞれ奇数。
【００２４】
なお、拡大後のフィルタのサイズは下式で表される。
Nfx'×Nfy' = (Nx×Lx)×(Ny×Ly)      …(2)
【００２５】
次に、CPU101は、分割フィルタごとに、フィルタ係数をフィルタ中心からの距離の多項式関数で近似して、多項式の各項の係数を求める(S105)。なお、元のフィルタBを近似多項式で表したフィルタをフィルタB'とする。そして、CPU101は、詳細は後述するが、複数の分割領域、複数の分割フィルタおよび近似多項式を用いて、画像データAにフィルタB'によるフィルタ処理を施す(S106)。
【００２６】
●近似多項式
ここで、近似多項式の作成について説明する。
【００２７】
フィルタが二階微分が可能な関数で表される場合、各分割フィルタの中心周りでテイラ展開を行って近似多項式を作成すればよい。より一般的には、誤差の二乗和を最小にするように近似多項式の係数を決める最小二乗法を利用する。最小二乗法は、フィルタの関数が不連続または微分不可能な点を含む場合や、フィルタ係数が中心からの距離の関数ではなく予め画素に対応するセルごとに与えられている場合にも適用可能である。
【００２８】
ここで、テイラ展開を用いて近似多項式を求める例を説明する。例えば、フィルタBのフィルタ係数が、フィルタ中心からの距離(x, y)の関数f(x, y)として表されるとする。関数f(x, y)を二階微分が可能な関数とし、位置(a, b)周りで二次の項までテイラ展開すると、下式が得られる。
f(x, y) ≒ f(a, b) + fx(a, b)(x-a) + fy(a, b)(y-b)
+ fxx(a, b)(x-a)2/2 + fyy(a, b)(y-b)2/2
+ fxy(a, b)(x-a)(y-b)              …(3)
ここで、fx = ∂f/∂x、fy = ∂f/∂y、
fxx = ∂2f/∂x2、fyy = ∂2f/∂y2、
fxy = ∂2f/∂x/∂y
【００２９】
さらに、下記を定義して式(3)を変形すると式(4)が得られる。
axx = fxx(a, b)/2、ayy = fyy(a, b)/2、axy = fxy(a, b)、
ax = fx(a, b) - 2afxx(a, b) - bfxy(a, b)、ay = fy(a, b) - 2bfyy(a, b) - afxy(a,b)、
a0 = f(a, b) - afx(a, b)- bfy(a, b) + a2fxx(a, b)/2 + abfxy(a, b) + b2fyy(a, b)/2、
f(x, y) ≒ axxx2 + axyxy + ayyy2 + axx + ayy + a0      …(4)
【００３０】
なお、上記には二次多項式の例を示したが、これに限定されるものではなく、一次、あるいは、三次または四次以上の次数をもつ多項式でもよい。例えば一次の項までの近似多項式であれば、式(4)に代わり式(5)を用いることになる。式(5)を用いれば、項の数は半分になるが、近似精度（元の関数f(x, y)と(3)や(4)式で示す多項式との誤差）が低下する。
f(x, y) ≒ axx + ayy + a0      …(5)
【００３１】
他方、三次や四次以上の次数をもつ近似多項式を用いれば、近似精度は向上するが、例えば、項の数は三次で10、四次で15と増える分、計算量が増大する。一般にf(x, y)をn次多項式で近似すると式(6)のように表される。
f(x, y) ≒ ΣpΣqapqxpqj      …(6)
ここで、Σpの範囲はp=0からn、
Σqの範囲はq=0からn、
nは多項式の次数。
【００３２】
フィルタBをLL個の分割フィルタに分割した場合、式(6)の近似多項式もLL個生成され、係数{apq}の組もLL組得られる。なお、LL個の近似多項式のすべてが同じ次数をもつ必要はなく、例えば、フィルタ中心に近い分割フィルタは二次の近似多項式とし、フィルタ中心から遠い分割フィルタは一次の近似多項式としてもよい。
【００３３】
●フィルタ処理
図2Bのフローチャートによりフィルタ処理(S106)の詳細を説明する。
【００３４】
近似多項式の作成(S105)の後、CPU101は、領域分割した画像データAに対し、フィルタB'の近似多項式を使って各画素とフィルタの畳込演算を行う。CPU101は、カウンタiを1に初期化し(S111)、画像データAのi番目の分割領域Rdi（第一の分割領域）を選択する(S112)。
【００３５】
図4により畳込演算における画像データAとフィルタB'の関係を説明する。図4に示すように、CPU101は、分割領域Rdiとフィルタ中心が一致するように画像データAにフィルタB'を重ね、分割領域RdiとフィルタB'による畳込演算の対象領域Rfを決定する(S113)。そして、カウンタjを1に初期化し(S114)、対象領域Rfに含まれるj番目の分割領域Rdj（第二の分割領域）を選択する(S115)。
【００３６】
図4はカウンタiに対応する分割領域Rdiを中心とする対象領域Rfにおいて、カウンタjに対応する分割領域Rdjが選択された様子を示す。CPU101は、分割領域Rdiの位置に対する分割領域Rdjの位置に基づきフィルタB'の近似多項式の係数の組を選択する(S116)。
【００３７】
分割領域Rdi内の画素(x1, y1)を注目画素とすると、注目画素とフィルタB'の畳込演算は次のように表される。
h'(x1, y1) = h(x1, y1)＊f(x1, y1)
= Σx'Σy'h(x', y')f(x'-x1, y'-y1)      …(7)
ここで、h(x1, y1)は注目画素の画素値、
h'(x1, y1)はフィルタ処理後の画素値、
f(x1, y1)はフィルタ係数、
＊は畳込演算を表す。
【００３８】
式(7)において、分割領域Rdjに含まれる画素からの寄与分は、フィルタB'が二次多項式(4)で近似されているものとし、分割領域Rdjに含まれる画素の位置を(x2, y2)として、下式で表される。
h'(x1, y1)|Rdj = ΣΣh(x2, y2)f(x2-x1, y2-y1)
≒ ΣΣh(x2, y2){axx(x2-x1)2 + axy(x2-x1)(y2-y1)
+ ayy(y2-y1)2 + ax(x2-x1) + ay(y2-y1) + a0}
= ΣΣh(x2, y2){axxx12 + axyx1y1 + ayyy12
+ (-2axxx2 - axyy2 - ax)x1 + (-2ayyy2 - axyx2 - ay)y1
+ axxx22 + axyx2y2 + ayyy22 + axx2 + ayy2 + a0}    …(8)
ここで、Σの範囲は分割領域Rdjに含まれる画素(x2, y2)の範囲。
【００３９】
さらに、下記を定義して式(8)を変形すると式(9)が得られる。
Gxx = ΣΣh(x2, y2)x22、Gxy = ΣΣh(x2, y2)x2y2、Gyy = ΣΣh(x2, y2)y22、
Gx = ΣΣh(x2, y2)x2、Gy = ΣΣh(x2, y2)y2、G0 = ΣΣh(x2, y2)。
h'(x1, y1)|Rdj ≒ axxG0x12 + axyG0x1y1 + ayyG0y12
+ (-2axxGx - axyGy -axG0)x1 + (-2ayyGy - axyGx -ayG0)y1
+ (axxGxx + axyGxy + ayyGyy + axGx + ayGy + a0G0)
= φxxx12 + φxyx1y1 + φyyy12 + φxx1 + φyy1 + φ0      …(9)
ここで、φxx = axxG0、φxy = axyG0、φyy = ayyG0、
φx = -2axxGx - axyGy - axG0、φy = -2ayyGy - axyGx - ayG0、
φ0 = axxGxx + axyGxy + ayyGyy + axGx + ayGy + a0G0
【００４０】
式(9)はフィルタB'を、前述した二次多項式(4)で近似した式である。フィルタB'を、前述したより一般的なn次多項式(6)で近似した式は次のようになる。
h'(x1, y1)|Rdj = Σx2Σy2h(x2, y2)f(x2-x1, y2-y1)
≒ Σx2Σy2h(x2, y2){ΣpΣqapq(x2-x1)p(y2-y1)q}      …(10)
【００４１】
式(10)を、式(9)と同様にx1、y1について整理すると次式が得られる。
h'(x1, y1)|Rdj ≒ ΣpΣqφpqx1py1q      …(11)
【００４２】
式(11)の{φpq}(0≦p≦n、0≦q≦n)は、x1py1q項の係数であり、近似多項式(6)の係数{apq}と、分割領域Rdj内の画素の位置のみから得られる。図5により分割領域Rdi、Rdj周辺の拡大図を示す。図5示すように分割領域Rdi内の画素(x1, y1)とは異なる画素(x3, y3)を注目画素として、注目画素(x3, y3)とフィルタB'の畳込演算における分割領域Rdjからの寄与分を考えると、式(11)と同様に、次式が得られる。
h'(x3, y3)|Rdj ≒ ΣpΣqφpqx3py3q      …(12)
【００４３】
式(12)の{φpq}(0≦p≦n、0≦q≦n)の値は、前述したように、分割領域Rdiの画素(x1, y1)または画素(x3, y3)の値とは無関係であり、式(11)の{φpq}と同一の値を有する。つまり、分割領域Rdi内に含まれる各画素について、フィルタB'との畳込演算のうち、分割領域Rdjからの寄与分の計算は、以下の手順で行えばよい。
・分割領域Rdiと分割領域Rdjの位置関係からフィルタB'の近似多項式の係数の組を選択し(S116)、
・分割領域Rdjの係数{φpq}を計算し、得られた値を保存し(S117)、
・分割領域Rdiの各画素について、係数{φpq}および式(11)を用いて分割領域Rdjからの寄与分を計算する(S118)。
【００４４】
次に、CPU101は、カウンタjをインクリメントし(S119)、カウンタjのカウント値と分割フィルタの数LLを比較して(S120)、j≦LLならば処理をステップS115に戻す。また、j＞LLならばカウンタiをインクリメントし(S121)、カウンタjのカウント値と分割領域Rdの数Mを比較して(S122)、i≦Mならば処理をステップS112に戻す。
【００４５】
そして、j＞Mになると、CPU101は、式(11)によって演算した畳込演算の結果を画素ごとに加算して、画像データAとフィルタB'の畳込演算結果（フィルタ処理結果）の画像データA'として出力する(S123)。
【００４６】
［計算量］
フィルタを分割しない場合、注目画素一画素に対するフィルタ処理の計算回数はNfx×Nfy回であるから、一つの分割領域Rdが含む全画素の計算回数NRは次式によって表される。
NR = Nfx・Nfy・Nx・Ny      …(13)
【００４７】
一回の計算には、分割領域Rdが含む画素の値とフィルタ係数との乗算、乗算結果を新しい画素値として積算する処理が含まれる。また、本実施例のフィルタ処理における計算量は次のように見積られる。
【００４８】
●分割領域Rdjの係数{φpq}の計算(S117)
計算量は、分割領域の画素数とフィルタB'の分割フィルタの数の積に比例する。
α1・(Nx×Ny)(Nfx'/Nx)(Nfy'/Ny) = α1・Nfx'・Nfy'      …(14)
ここで、α1は一回当りの計算量を示す定数。
【００４９】
●分割領域Rdjからの寄与分の計算(S118)
つまり、係数{φpq}を用いた式(11)の計算であり、式(11)の計算を、一画素に付きフィルタB'に含まれる分割フィルタの数分繰り返し、さらに、分割領域の画素数分繰り返す。
α2・(Nfx'/Nx)(Nfy'/Ny)・Nx・Ny = α2・Nfx'・Nfy'      …(15)
ここで、α2は式(11)の一回当りの計算量を示す定数。
【００５０】
従って、式(14)と(15)を加算した結果が、本実施例のフィルタ処理における、分割領域とフィルタB'の畳込演算に要する計算量LFである。
LF = α・Nfx'・Nfy'      …(16)
ここで、α=α1・α2
【００５１】
また、フィルタを分割しない場合の計算数（式(13)）と式(16)から次式が得られる。
α・Nfx'・Nfy'/(Nfx・Nfy・Nx・Ny) = α'/(Nx・Ny)      …(17)
ここで、α' = α・(Nfx'/Nfx)・(Nfy'/Nfy)
【００５２】
フィルタB'を奇数個の分割フィルタによって構成するためにフィルタB'のサイズはフィルタBのサイズよりも大きくなる。フィルタサイズの拡大が計算量に影響する分としてα'を加味する。例えば、α'が9.0程度、Nx=Ny=8とすると、本実施例のフィルタ処理は、計算量を9/(8×8)≒1/7に削減することができる。
【００５３】
●計算量が削減される理由
式(17)に示すように計算量を削減できる理由は、畳込演算の式(7)または式(8)を、式(11)に示す近似式で置き換えたことにある。
【００５４】
ここで、フィルタ係数がフィルタ中心からの距離で決まるようなフィルタに対して、フィルタを分割せずに畳込演算を行う場合を考える。フィルタ処理する注目画素(x1, y1)が変わると、フィルタ範囲に含まれる周辺画素(x2, y2)のフィルタ中心からの距離が変わり、周辺画素(x2, y2)の寄与分も変化する。従って、周辺画素(x2, y2)からの寄与分Cは下式で表される。
C = h(x2, y2)×f(x2-x1, y2-y1)      …(18)
ここで、h(x2,y2)は周辺画素(x2,y2)の画素値、
f(x, y)はフィルタ関数、
(x1, y1)は注目画素の位置。
【００５５】
注目画素が(x3, y3)に移動すると、周辺画素(x2, y2)からの寄与分Cは下式で表される。
C = h(x2, y2)×f(x2-x3, y2-y3)      …(19)
【００５６】
つまり、同じ周辺画素(x2, y2)からの寄与分Cを計算するとしても、注目画素の移動に伴いフィルタ関数f(x, y)の値が変化するため、注目画素ごとに寄与分Cを計算する必要が生じる。
【００５７】
一方、式(11)においては、{φpq}が分割領域Rdj内の画素位置と画素値のみで決まり、フィルタ処理する分割領域Rdi内の注目画素の位置(x1, y1)と分離されている。従って、注目画素が移動した場合も、分割領域Rdj内の画素からの寄与分を表す{φpq}は同一（式(11)(12)参照）であり、{φpq}の再計算は不要で、{φpq}は定数と見做すことができる。また、式(11)に示す表現方法は、元のフィルタを二次多項式(4)で近似する場合の式の変形例（式(8)(9)）において示したように、多項式関数で置き換えることで、容易に実現することができる。
【００５８】
一般に、任意の連続関数は多項式で近似することが可能であるが、元の関数と近似関数の誤差を小さくすれば、多項式の次数を大きくする必要があり計算量が増大する。そこで、元のフィルタの範囲を領域分割し、一つの分割フィルタの範囲における近似多項式を求めるようにすれば、多項式の次数を大きくしなくても、分割フィルタの範囲内において充分に誤差が小さい近似多項式を得ることができる。さらに、予め近似多項式の形を決めておき、近似多項式の係数を最小二乗法を用いて求める手法を用いれば、元のフィルタは必ずしも連続関数である必要はない。
【００５９】
なお、式(13)から式(17)に示した計算量について、一回の計算量は、フィルタを分割しない場合は前述したように「一回の乗算」と「一回の加算」の和である。一方、本実施例における一回の計算量は、例えば二次の近似多項式を用いる式(9)においては、加算と乗算がそれぞれ五回ずつ含まれる。また、フィルタのサイズが分割フィルタ単位になるように、Nfx×NfyからNfx'×Nfy'に拡大する影響もある。そのため、式(17)の定数αは1よりも大きくなり、例えば二次の近似多項式を用いると9.0程度になる。
【００６０】
［フィルタ範囲の拡大］
本実施例におけるフィルタ処理は、奇数個の分割フィルタ単位に行うため、フィルタB'の範囲は、フィルタBの範囲に比べて拡大される。図6によりフィルタ範囲の拡大を説明する。
【００６１】
図6は、実線で示す72×80画素の画像を8×8画素に分割した画像Aに、破線で示すサイズが49×49画素のフィルタBを8×8画素のフィルタに分割した7×7の分割フィルタを有するフィルタB'を重ねた状態を示す。この場合、フィルタB'のサイズは56×56画素になる。
【００６２】
図6に示す分割領域の左上画素(32, 40)が注目画素501の場合、フィルタB（49×49画素）の範囲は、点線で示す左上(8, 16)、右下(56, 64)の範囲である。また、右下画素(39, 47)が注目画素502の場合、フィルタBの範囲は、点線で示す左上(15, 23)、右下(63, 71)の範囲である。一方、本実施例においては、同じ分割領域内の注目画素501、502に対してフィルタB'の範囲は変化せず、破線で示す左上(8, 16)、右下(63, 71)の範囲である。なお、本実施例において、注目画素501についてx=57より右、y=65より下の画素が参照範囲に追加され、注目画素502についてx=14より左、y=22より上の画素が参照画素に追加される。
【００６３】
このようなフィルタ範囲の違いは、畳込演算を分割フィルタ単位に実施するために生じる誤差であり、フィルタ範囲を拡大した部分のフィルタ係数がほぼ0であれば、誤差の影響は僅かである。従って、フィルタ中心から遠いフィルタ範囲の端部の係数がほぼ零になるフィルタに、本実施例を適用することが最も好ましい。
【００６４】
［処理結果］
●ガウシアンフィルタ
図7により処理結果の確認に使用した画像Aの一例を示す。図7に示す画像Aは、2048×2048画素の8ビットモノクロ画像であり、ランダムな位置に、ランダムな大きさの正方形や円が重畳し、さらにランダムなノイズが加わっている。また、フィルタBは、下式に示すガウシアンフィルタである。
f(x, y) = C1×exp{-(x2+y2)/2/σ2}      …(20)
ここで、C1は定数。
【００６５】
フィルタB（式(20)）のσは12画素とし、フィルタサイズは73×73画素である。フィルタBの係数f(x, y)の合計が1.0になるように定数C1=0.00249とする。また、分割領域の大きさは8×8画素である。フィルタB'の各分割フィルタを式(4)で示す二次多項式で近似し、フィルタBのサイズは11×11の分割フィルタである。つまり、フィルタBを73×73画素から88×88画素に拡大してフィルタB'にする。
【００６６】
画像Aに対し、フィルタBによる処理を行って画像Cを取得し、フィルタB'による処理を行って画像Dを取得した。そして、画像Cと画像Dの各画素ごとに画素値の差分を取り、全画像領域における差分値の最大、最小値を求めた。また、各画素の差分値の二乗平均平方根(RMS)を下式によって計算した。
RMS = √[ΣxΣy{hC'(x, y) - hD'(x, y)}2}/NA]      …(21)
ここで、hC(x, y)は画像Cの画素値、
hD(x, y)は画像Dの画素値、
NAは画素数。
【００６７】
その結果、全画像領域における差分値の最大は1、最小は-1であり、画素値の差は±1の範囲に収まり、RMS値は0.196であった。さらに、画像Dおよび差分値の分布に対して高速フーリエ変換(FFT)を施して空間周波数成分の分布を調べたところ、分割領域の大きさに関連するような空間周波数特性は見出されなかった。従って、画像Cと画像Dは殆ど同一の画像であり、フィルタ処理の違いによる差は殆ど生じないことを確認した。
【００６８】
また、フィルタB'による処理時間は、フィルタBによる処理時間の1/7.5であった。
【００６９】
●LoGフィルタ
図8により処理結果の確認に使用した画像Aの一例を示す。図8に示す画像Aは、1024×1024画素の8ビットモノクロ画像であり、ランダムな位置に、ランダムな大きさの正方形や円が重畳し、さらにランダムなノイズが加わっている。また、フィルタBは、下式に示すLoG (Laplacian of Gaussian)フィルタである。
f(x, y) = (x2 + y2 - 2σ2)／(2πσ6)×exp{-(x2+y2)/2/σ2}      …(22)
【００７０】
フィルタB（式(22)）のσは24画素とし、フィルタサイズは145×145画素である。また、分割領域の大きさは8×8画素である。フィルタB'の各分割フィルタを式(4)で示す二次多項式で近似し、フィルタBのサイズは19×19の分割フィルタである。つまり、フィルタBを145×145画素から152×152画素に拡大してフィルタB'にする。
【００７１】
図9により画像Aに対し、フィルタBによる処理を行って取得した画像C（図9(a)）と、フィルタB'による処理を行って取得した画像D（図9(b)）を示す。図8において、フィルタ処理によって得られた正の値の画素を黒色、負の値の画素を白色で表す。式(22)におけるσ値が大きいため、図9における正値の画素が構成する線の幅が太くなっているが、細線化すれば図8に示す図形のエッジに相当することは明らかである。また、図9(a)、図9(b)とも、同様に、エッジが抽出されている。
【００７２】
また、フィルタB'による処理時間は、フィルタBによる処理時間の1/6.2であった。
【００７３】
●分割領域のサイズの変更
以上では、分割領域のサイズが8×8画素の例を説明したが、以下では、分割領域のサイズを変化させた場合に、処理時間、および、フィルタを分割しない場合との誤差がどのように変化するかについて説明する。
【００７４】
画像Aは、図7に示す2048×2048画素の8ビットモノクロ画像である。また、フィルタBは式(20)に示すガウシアンフィルタ（73×73画素）である。そして、分割領域のサイズを6×6、8×8、10×10、12×12、14×14、16×16、18×18、20×20画素の八種類とする。分割領域のサイズを変更しながら、フィルタBを用いる処理とフィルタB'を用いる処理について、精度と処理時間の比較を行った。
【００７５】
図10によりフィルタBを用いる処理とフィルタB'を用いる処理の比較結果を説明する。図10において、横軸は「分割領域の辺の大きさ/フィルタのσ値」、左の縦軸は処理時間比（=フィルタB'を用いる処理/フィルタBを用いる処理）、右の縦軸は式(21)に示すRMS値をそれぞれ表す。
【００７６】
図10によれば、フィルタB'を用いる処理は、分割領域のサイズを大きくする程、処理時間が短くなり、分割領域のサイズを小さくする程、誤差が小さくなることがわかる。
【００７７】
図11により分割領域のサイズを20×20画素として、フィルタ処理した後の画像（図11(a)）および空間周波数分布（図11(b)）を示す。なお、図11(a)に示す画像は分割領域の20画素間隔の縞模様が見えるように部分的に拡大したものである。
【００７８】
図11(a)を参照すると、XY軸に対してそれぞれ±20画素、および、±20画素/整数に相当するピッチの模様が生じている。これは分割領域が元のフィルタのサイズに比べて大きいため、近似多項式の精度が低下し、誤差が生じているためと考えられる。
【００７９】
このように、フィルタを分割する本実施例のフィルタ処理は、分割領域のサイズを大きくする程、処理時間の短縮が期待できるが、その一方、誤差が大きくなり精度は低下する。従って、分割領域のサイズを適切に決めることにより、精度を低下させずに、フィルタ処理の処理時間を短縮することが可能になる。
【実施例２】
【００８０】
以下、本発明にかかる実施例2の画像処理を説明する。なお、実施例2において、実施例1と略同様の構成については、同一符号を付して、その詳細説明を省略する。
【００８１】
実施例1においては二次元フィルタ処理を説明したが、実施例2においては、図12Aのフローチャートにより一次元フィルタ処理を説明する。
【００８２】
CPU101は、一次元のデータA（例えばサウンドデータ）を入力し(S201)、入力データAをNxデータ（所定データサイズ）ごとに分割した分割データを生成する(S202)。なお、データAが画像データのような二次元データの場合、例えば画像の右端において、Nx画素に満たない分割データが発生する可能性があるが、そのような場合も分割データとして扱う。また、以下の説明では、分割データの数をMとする。
【００８３】
次に、CPU101は、一次元フィルタ処理用のフィルタBを表すデータを例えばHDD108から入力する(S203)。フィルタBのフィルタサイズをNfxデータとする。フィルタBは、注目データからの距離の関数として表現されてもよいし、フィルタの各成分（フィルタ係数）を表すデータでもよい。ただし、フィルタBのフィルタサイズを分割データの数の自然数倍に拡大するため、フィルタ処理の裾に相当する部分の成分（フィルタ係数）がほぼ零であることが好ましい。
【００８４】
次に、CPU101は、分割データの数Nxに対応させて、フィルタBを分割する(S204)。その際、CPU101は、分割後のフィルタが奇数個の分割フィルタ（サイズNx）から構成されるように、フィルタサイズを拡大する。つまり、分割フィルタの数をLxとすると、Lxが下式を満たすようにフィルタサイズを拡大する。
Nx×Lx ≧ Nfx            …(23)
ここで、Lxは奇数。
【００８５】
なお、拡大後のフィルタのサイズは下式で表される。
Nfx' = Nx×Lx            …(24)
【００８６】
次に、CPU101は、分割フィルタごとに、フィルタ係数をフィルタ中心からの距離の多項式関数で近似して、多項式の各項の係数を求める(S205)。なお、元のフィルタBを近似多項式で表したフィルタをフィルタB'とする。そして、CPU101は、詳細は後述するが、前記複数の分割データ、複数の分割フィルタおよび近似多項式を用いて、データAにフィルタBによるフィルタ処理を施す(S206)。
【００８７】
●近似多項式
フィルタBが二階微分が可能な関数として表記されている場合、例えば、各分割フィルタの中心に対してテイラ展開を行って近似多項式を求める。例えば、フィルタBのフィルタ係数が、フィルタ中心からの距離xの関数f(x)として表されるとする。関数f(x)を二階微分が可能な関数とし、位置aに対して二次の項までテイラ展開すると、下式が得られる。
f(x) ≒ f(a) + fx(a)(x-a) + fxx(a)(x-a)2/2      …(25)
ここで、fx = df/dx、
fxx = d2f/dx2
【００８８】
さらに、下記を定義して式(25)を変形すると式(26)が得られる。
axx = fxx(a)/2、
ax = fx(a) - 2afxx(a)、
a0 = f(a) - afx(a)- a2fxx(a)/2、
f(x, y) ≒ axxx2 + axx + a0      …(26)
【００８９】
なお、上記には二次多項式の例を示したが、これに限定されるものではなく、一次、あるいは、三次または四次以上の次数をもつ多項式でもよい。一般にf(x)をn次多項式で近似すると式(27)のように表される。
f(x) ≒ Σpapxp      …(27)
ここで、Σpの範囲はp=0からn、
nは多項式の次数。
【００９０】
フィルタBをLx個の分割フィルタに分割した場合、式(27)の近似多項式もLx個生成され、係数{ap}の組もLx組得られる。なお、Lx個の近似多項式のすべてが同じ次数をもつ必要はなく、例えば、フィルタ中心に近い分割フィルタは二次の近似多項式とし、フィルタ中心から遠い分割フィルタは一次の近似多項式としてもよい。
【００９１】
●フィルタ処理
図12Bのフローチャートによりフィルタ処理(S206)の詳細を説明する。
【００９２】
近似多項式の作成後、CPU101は、分割したデータAに対し、フィルタB'の近似多項式を使って各データとフィルタの畳込演算を行う。CPU101は、カウンタiを1に初期化し(S211)、データAのi番目の分割データDdi（第一の分割データ）を選択する(S212)。
【００９３】
CPU101は、分割データDdiとフィルタ中心が一致するようにデータAにフィルタB'を重ね、分割データDdiとフィルタB'による畳込演算の対象データ集合Sfを決定する(S213)。そして、カウンタjを1に初期化し(S214)、対象データ集合Sfに含まれるj番目の分割データDdj（第二の分割データ）を選択する(S215)。
【００９４】
CPU101は、分割データDdiの位置に対する分割データDdjの位置に基づきフィルタB'の近似多項式の係数の組を選択する(S216)。分割データDdi内のデータx1を注目データとすると、注目データとフィルタB'の畳込演算は次のように表される。
h'(x1) = h(x1)＊f(x1)
= Σx'h(x')f(x'-x1)      …(28)
ここで、h(x1)は注目データの値、
h'(x1)はフィルタ処理後の値、
f(x1)はフィルタ係数、
＊は畳込演算を表す。
【００９５】
式(28)において、分割データDdjに含まれるデータからの寄与分は、フィルタB'が二次多項式(26)で近似されているものとし、分割データDdjに含まれるデータの位置をx2として、下式で表される。
h'(x1)|Ddj = Σh(x2)f(x2-x1)
= Σh(x2){axx(x2-x1)2 + ax(x2-x1) + a0}
= Σh(x2){axxx12 + (-2axxx2 - ax)x1 + axxx22 + axx2 + a0}  …(29)
ここで、Σの範囲は分割データDdjに含まれるデータx2の範囲。
【００９６】
さらに、下記を定義して式(29)を変形すると式(30)が得られる。
Gxx = Σh(x2, y2)x22、Gx = Σh(x2, y2)x2、G0 = Σh(x2, y2)
h'(x1)|Ddj ≒ axxG0x12 + (-2axxGx - axG0)x1 + (axxGxx + axGx + a0G0)
≒ φxxx12 + φxx1 + φ0      …(30)
ここで、φxx = axxG0、φx = -2axxGx - axG0、φ0 = axxGxx + axGx + a0G0
【００９７】
式(30)はフィルタB'を、前述した二次多項式(26)で近似した式である。フィルタB'を、前述したより一般的なn次多項式(27)で近似した式は次のようになる。
h'(x1)|Ddj = Σx2h(x2)f(x2-x1)
≒ Σx2h(x2){Σpap(x2-x1)p}      …(31)
【００９８】
式(31)を、式(30)と同様にx1について整理すると次式が得られる。
h'(x1)|Ddj ≒ Σpφpx1p      …(32)
【００９９】
式(32)の{φp}(0≦p≦n)は、x1p項の係数であり、近似多項式(26)の係数{ap}と、分割データDdj内のデータの位置のみから得られる。つまり、分割データDdi内のデータx1とは異なるデータx3を注目データとして、注目データx3とフィルタB'の畳込演算における分割データDdjからの寄与分を考えると、式(32)と同様に、次式が得られる。
h'(x3)|Ddj ≒ Σpφpx3p      …(33)
【０１００】
式(33)の{φp}(0≦p≦n)の値は、前述したように、分割データDdiのデータx1またはデータx3の値とは無関係であり、式(31)の{φp}と同一の値を有する。つまり、分割データDdi内に含まれる各画素について、フィルタB'との畳込演算のうち、分割データDdjからの寄与分の計算は、以下の手順で行えばよい。
・分割データDdiと分割データDdjの位置関係からフィルタB'の近似多項式の係数の組を選択し(S216)、
・分割データDdjの係数{φp}を計算し、得られた値を保存し(S217)、
・分割データDdiの各データについて、係数{φp}および式(32)を用いて分割データDdjからの寄与分を計算する(S218)。
【０１０１】
次に、CPU101は、カウンタjをインクリメントし(S219)、カウンタjのカウント値と分割フィルタの数Lxを比較して(S220)、j≦Lxならば処理をステップS215に戻す。また、j＞Lxならばカウンタiをインクリメントし(S221)、カウンタjのカウント値と分割データRdの数Mを比較して(S222)、i≦Mならば処理をステップS212に戻す。
【０１０２】
そして、j＞Mになると、CPU101は、式(32)によって演算した畳込演算の結果をデータごとに加算して、データAとフィルタB'の畳込演算結果（フィルタ処理結果）のデータA'として出力する(S223)。
【０１０３】
［計算量］
フィルタを分割しない場合、注目データに対するフィルタ処理の計算回数はNfx回であるから、一つの分割データDdが含む全データの計算回数NRは次式によって表される。
NR = Nfx・Nx      …(34)
【０１０４】
一回の計算には、分割データDdが含むデータとフィルタ係数との乗算、乗算結果を新しいデータとして積算する処理が含まれる。また、本実施例のフィルタ処理における計算量LFは、実施例1で説明したように、次のように見積られる。
LF = α・Nfx'      …(35)
【０１０５】
また、フィルタを分割しない場合の計算数（式(34)）と式(35)から次式が得られる。
α・Nfx'/(Nfx・Nx) = α'/Nx      …(36)
ここで、α' = α・(Nfx'/Nfx)
【０１０６】
フィルタB'を奇数個の分割フィルタによって構成するためにフィルタB'のサイズはフィルタBのサイズよりも大きくなる。フィルタサイズの拡大が計算量に影響する分としてα'を加味する。例えば、α'が4程度、Nx=8とすると、本実施例のフィルタ処理は、計算量を4/8=1/2に削減することができる。
【０１０７】
［その他の実施例］
また、本発明は、以下の処理を実行することによっても実現される。即ち、上述した実施形態の機能を実現するソフトウェア（プログラム）を、ネットワーク又は各種記憶媒体を介してシステム或いは装置に供給し、そのシステムあるいは装置のコンピュータ（又はCPUやMPU等）がプログラムを読み出して実行する処理である。
【図２Ａ】
0005620829.tif 000002
【図２Ｂ】
0005620829.tif 000003
【図１０】
0005620829.tif 000004
【図１２Ａ】
0005620829.tif 000005
【図１２Ｂ】
0005620829.tif 000006
【図１】
0005620829.tif 000007
【図３】
0005620829.tif 000008
【図４】
0005620829.tif 000009
【図５】
0005620829.tif 000010
【図６】
0005620829.tif 000011
【図７】
0005620829.tif 000012
【図８】
0005620829.tif 000013
【図９】
0005620829.tif 000014
【図１１】
0005620829.tif 000015