(57)【要約】【課題】従来よりも計算コストが小さく、かつ計算誤差の制御を自在に行うことが可能な画像フィルタ演算装置の提供。【解決手段】入力画像{x(p)}に対し、空間カーネルgs及びレンジ・カーネルgrのコンボリューション演算により画像平滑化を行う画像フィルタ演算装置であって、スケール・パラメータσｒ、許容誤差τ及びダイナミックレンジＲに対し、レンジ・カーネルgrをフーリエ級数近似した近似レンジ・カーネルg^rにおける打切項数Ｋ及び有意周期Ｔの組（Ｔ，Ｋ）を、近似誤差がτ以下となるように決定する近似パラメータ最適化手段と、近似レンジ・カーネルg^rの展開係数{a0,…,aKτ}を算出する展開係数演算手段と、展開係数{a0,…,aKτ}を用いコンボリューション演算を実行するフィルタ演算手段を備えた。【選択図】図１２
2016122430.tif 000002
【特許請求の範囲】
【請求項１】
入力画像{x(p)}（pは画素の位置，x(p)は位置pの画素値）に対して、標準偏差（以下「スケール・パラメータ」という。）σｓのガウス関数で定義される空間カーネルgs(q-p)（p，qは画素の位置）、及びスケール・パラメータσｒのガウス関数で定義されるレンジ・カーネルgr(x(q)-x(p))とのコンボリューション演算により画像の平滑化を行う画像フィルタ演算装置であって、
前記スケール・パラメータσｒ、所定の許容誤差τ及び前記入力画像{x(p)}のダイナミックレンジＲに対し、前記レンジ・カーネルgr(x(q)-x(p))を、有限項数のフーリエ級数により近似した近似レンジ・カーネルg^r(x(q)-x(p))における級数の打切項数Ｋ、及び前記近似レンジ・カーネルg^r(x(q)-x(p))における各項の展開係数ａｋ（ｋ＝０，…，Ｋ）を演算する際のフーリエ積分の積分領域［-T/2，T/2］を定める有意周期Ｔの組である近似パラメータ（Ｔ，Ｋ）を、前記近似レンジ・カーネルg^r(x(q)-x(p))の近似誤差が前記許容誤差τ以下となるように決定し出力する近似パラメータ最適化手段と、
前記近似パラメータ最適化手段が算出する前記近似パラメータ（Ｔ，Ｋ）に従い、前記近似レンジ・カーネルg^r(x(q)-x(p))の第０項から第Ｋ項までの展開係数{a0,…,aＫ}を算出する展開係数演算手段と、
前記展開係数演算手段が算出する展開係数{a0,…,aＫ}を用いて、前記入力画像{x(p)}に対し、前記空間カーネルgs(q-p)及び前記近似レンジ・カーネルg^r(x(q)-x(p))とのコンボリューション演算を実行するフィルタ演算手段と、を備えた画像フィルタ演算装置。
【請求項２】
前記近似パラメータ最適化手段は、前記打切項数Ｋを算出する打切項数演算手段を備え、
前記打切項数演算手段は、前記スケール・パラメータσｒ、所定の許容誤差τ及び前記入力画像{x(p)}のダイナミックレンジＲに対し、前記打切項数Ｋを
【数１】
2016122430.tif 000058
により算出することを特徴とする請求項１記載の画像フィルタ演算装置。
【請求項３】
前記近似パラメータ最適化手段は、前記有意周期Ｔを算出する有意周期最適化手段を備え、
前記展開係数演算手段は、前記レンジ・カーネルgr(x(q)-x(p))を積分領域［-T/2，T/2］でフーリエ積分した近似展開係数を、前記各展開係数{a0,…,aＫ}として算出するものであり、
前記有意周期最適化手段は、前記スケール・パラメータσｒ、前記打切項数Ｋ及び前記ダイナミックレンジＲに対する、変数Ｔの方程式
【数２】
2016122430.tif 000059
を満たすＴを前記有意周期Ｔとして算出するものであることを特徴とする請求項１又は２記載の画像フィルタ演算装置。
【請求項４】
前記近似パラメータ最適化手段は、
前記スケール・パラメータσｒ及び前記レンジ上限値Ｒ、並びに予め設定された有効幅パラメータα（ασｒ＜Ｒ）に基づき、前記有意周期Ｔを、
【数３】
2016122430.tif 000060
により算出する演算周期設定手段と、
前記許容誤差τ，前記有意周期Ｔ及び前記スケール・パラメータσｒに基づいて、前記打切項数Ｋを
【数４】
2016122430.tif 000061
により算出する打切項数決定手段と、を備えたことを特徴とする請求項１記載の画像フィルタ演算装置。
【請求項５】
前記展開係数演算手段は、前記各展開係数{a1,…,aＫ}を、
【数５】
2016122430.tif 000062
により算出するものであることを特徴とする請求項１乃至４の何れか一記載の画像フィルタ演算装置。
【請求項６】
コンピュータに読み込ませて実行することにより、前記コンピュータを請求項１乃至５の何れか一に記載の画像フィルタ演算装置として機能させることを特徴とするプログラム。
【請求項７】
標準偏差（以下「スケール・パラメータ」という。）σｒのガウス関数で定義されるガウシアン・カーネルgr(x)の近似値を算出するガウシアン・カーネル演算装置であって、
前記スケール・パラメータσｒ、所定の許容誤差τ及び前記入力画像{x(p)}のダイナミックレンジＲに対し、前記レンジ・カーネルgr(x(q)-x(p))を、有限項数のフーリエ級数により近似した近似レンジ・カーネルg^r(x(q)-x(p))における級数の打切項数Ｋ、及び前記近似レンジ・カーネルg^r(x(q)-x(p))における各項の展開係数ａｋ（ｋ＝０，…，Ｋ）を演算する際のフーリエ積分の積分領域［-T/2，T/2］を定める有意周期Ｔの組である近似パラメータ（Ｔ，Ｋ）を、前記近似レンジ・カーネルg^r(x(q)-x(p))の近似誤差が前記許容誤差τ以下となるように決定し出力する近似パラメータ最適化手段と、
前記近似パラメータ最適化手段が算出する前記近似パラメータ（Ｔ，Ｋ）に従い、前記近似ガウシアン・カーネルgr(x)の第０項から第Ｋ項までの展開係数{a0,…,aＫ}を算出する展開係数演算手段と、
前記展開係数演算手段が算出する展開係数{a0,…,aＫ}を用いて、前記近似ガウシアン・カーネルgr(x)を算出する近似関数演算手段と、を備えたガウシアン・カーネル演算装置。
【請求項８】
前記近似パラメータ最適化手段は、前記打切項数Ｋを算出する打切項数演算手段を備え、
前記打切項数演算手段は、前記スケール・パラメータσｒ、所定の許容誤差τ及び前記入力画像{x(p)}のダイナミックレンジＲに対し、前記打切項数Ｋを
【数６】
2016122430.tif 000063
により算出することを特徴とする請求項７記載の画像フィルタ演算装置。
【請求項９】
前記近似パラメータ最適化手段は、前記有意周期Ｔを算出する有意周期最適化手段を備え、
前記展開係数演算手段は、前記近似ガウシアン・カーネルgr(x)を積分領域［-T/2，T/2］でフーリエ積分した近似展開係数を、前記各展開係数{a0,…,aＫ}として算出するものであり、
前記有意周期最適化手段は、前記スケール・パラメータσｒ、前記打切項数Ｋ及び前記ダイナミックレンジＲに対する、変数Ｔの方程式
【数７】
2016122430.tif 000064
を満たすＴを前記有意周期Ｔとして算出するものであることを特徴とする請求項８記載のガウシアン・カーネル演算装置。
【請求項１０】
前記近似パラメータ最適化手段は、
前記スケール・パラメータσｒ及び前記レンジ上限値Ｒ、並びに予め設定された有効幅パラメータα（ασｒ＜Ｒ）に基づき、前記有意周期Ｔを、
【数８】
2016122430.tif 000065
により算出する演算周期設定手段と、
前記許容誤差τ，前記有意周期Ｔ及び前記スケール・パラメータσｒに基づいて、前記打切項数Ｋを
【数９】
2016122430.tif 000066
により算出する打切項数決定手段と、を備えたことを特徴とする請求項７記載のガウシアン・カーネル演算装置。
【請求項１１】
前記展開係数演算手段は、前記各展開係数{a1,…,aＫ}を、
【数１０】
2016122430.tif 000067
により算出するものであることを特徴とする請求項７乃至１０の何れか一記載のガウシアン・カーネル演算装置。
【請求項１２】
コンピュータに読み込ませて実行することにより、前記コンピュータを請求項７乃至１１の何れか一に記載のガウシアン・カーネル演算装置として機能させることを特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、入力されるサンプル値ｘ及び入力されるスケール・パラメータσに対し、ガウシアン・カーネルｇ（σ，ｘ）の演算を実行するガウシアン・カーネル演算装置と、それを用いた画像フィルタ演算装置に関する。
【背景技術】
【０００２】
画像処理技術の分野において、エッジを保存しつつノイズの低減を図るためのフィルタとして、バイラテラル・フィルタが知られている。一般に、ＢＬＦは次のように定義される。
【０００３】
（定義１）（バイラテラル・フィルタ）
Ｄ次元グレースケール画像の入力画像（target image）を考える。Ｚを整数全体の集合，Ω⊂ＺＤを画像領域（image domain），ｐ∈Ωを画素位置（pixel position）（ベクトル），Ｎｐ⊂Ωを画素位置ｐの近傍領域（フィルタ窓領域）とする。Ｒを実数全体の集合とし、ｘ（ｐ），ｚ（ｐ）∈Ｒを、それぞれ、入力画像，出力画像の位置ｐの画素（以下「画素ｐ」）の画素値（pixel intensity）とする。このとき、次式により定義されるフィルタを「バイラテラル・フィルタ（bilateral filter：ＢＬＦ）」という。
【０００４】
【数１】
2016122430.tif 000003
ここで、ｇ（σ，ｘ）はガウシアン・カーネル（Gaussian kernel）、σｓは空間領域の相関距離を決めるスケール・パラメータ（標準偏差）、σｒは画素値領域の相関距離を決めるスケール・パラメータ（標準偏差）、Ｄはベクトルｑの次元数、‖ｑ‖はベクトルｑのノルム（l2-ノルム）である。また、ｇｓを「空間カーネル（spatial kernel）」、ｇｒを「レンジ・カーネル（range kernel）」という。
（定義終わり）
【０００５】
ＢＬＦは、画像処理，コンピュータ・ビジョン，コンピュータ・グラフィックスなどにおいて、エッジ保存スムージングにおける基礎的なツールである。ＢＬＦは、フィルタ係数を２つの側向（lateral）（即ち、画素位置（pixel position）及び画素値（pixcel intensity））から定義するという明確な概念に基づき設計されている。また、ＢＬＦは、代表的な他のエッジ保存スムージング手法である異方性拡散法（anisotropic diffusion）（非特許文献６）や重み付き最小二乗法（weighted least square）（非特許文献７）とは異なり、反復演算プロセスを用いることなく、高いスムージング効果を得ることができる。そのため、ＢＬＦは、ノイズ除去，ハイダイナミックレンジ合成（high dynamic range imaging：ＨＤＲ），ステレオマッチングなどの様々な用途に於いて広く活用されている。
【０００６】
従来、ＢＬＦにおいては、スムージング効果の向上や計算コストの減少を図るために、様々な拡張が行われている。
【０００７】
まず、ＢＬＦの自然な拡張として、クロス／ジョイント拡張（クロスＢＬＦ／ジョイントＢＬＦ）について考案されている（非特許文献８，９）。ノイズがある入力画像のフィルタリングにおいて、ＢＬＦではノイズがある入力画像からフィルタ係数を定義するために、スムージング効果が一定程度相殺される。それに対して、クロス／ジョイントＢＬＦでは、入力画像とは異なる撮影条件で撮影されたガイダンス画像を用いてフィルタ係数を定義する。これにより、より高いスムージング効果が得られる。
【０００８】
式（１）のＢＬＦにおいては、ｇｒ（・）は、入力画像の画素値ｘ（・）によって変化するため、フィルタ係数は、それぞれの対象画素において動的に変化する。本来のＢＬＦは、フィルタ係数をｘ（・）から計算するが、それに対してクロス／ジョイントＢＬＦでは、フィルタ係数をｘ（・）に対応するガイダンス画像の画素値ｙ（・）から計算する。
【０００９】
（定義２）（クロス／ジョイント・バイラテラル・フィルタ）
ｙ（ｐ）∈Ｒをガイダンス画像の位置ｐにおける画素値とする。このとき、次式により定義されるフィルタを「クロス／ジョイント・バイラテラル・フィルタ（cross/joint extension of the BLF）」という。
【００１０】
【数２】
2016122430.tif 000004
（定義終わり）
【００１１】
一方、計算コストの減少を測るためのアルゴリズム的拡張として、定時間ＢＬＦ（Ｏ（１）ＢＬＦ）を含む高速化ＢＬＦがこれまでに幾つか考案されている（非特許文献１０，１１，１３−１６）。式（１ａ）に示したオリジナルのＢＬＦにおいては、１画素当たりの計算コストがフィルタ窓領域Ω内の画素数Ｍ（＝｜Ω｜）に依存する（即ち、計算量がＯ（Ｍ）となる）。そのため、大きな計算コストを要する。近年の高解像度化や高次元画像処理へ向かう傾向から考えて、ＢＬＦを低コスト化・高速化することは重要な技術課題である。高速化ＢＬＦは、この課題への対処から考案されたものであり、若干の精度的な犠牲の下で、Ｏ（log Ｍ）又はＯ（１）の計算量で実行できるように改良したものである。
【００１２】
従来の高速化ＢＬＦでは、ＢＬＦが空間軸とレンジ軸とにより測られる高次元空間（バイラテラル空間）におけるコンボリューション（convolution）として一般化できるという事実に基づいて、様々な改良を行っている。線型フィルタに関しては、過去に多くの高速化アルゴリズムが考案されているが、バイラテラル空間を直接取り扱うためには大きなメモリ・コストが必要となる。そこで、従来の高速化ＢＬＦアルゴリズムでは、この問題を様々な戦略により解決している。
【００１３】
非特許文献１０には、ボックス空間フィルタ（box spatial filter）を制限しヒストグラム法（histgram technique）を使用することによって、バイラテラル空間を単純化した高速化ＢＬＦが開示されている。
【００１４】
非特許文献１１では、ヒストグラム法（非特許文献１２）に類似のアプローチによる定時間Ｏ（１）ＢＬＦが開示されている。非特許文献１１においては、任意空間フィルタ（arbitrary spatial filter）とテイラー級数を用いた多項式レンジ・カーネル（polynomial range kernel）について記載されている。
【００１５】
非特許文献１３には、空間軸とレンジ軸との双方において、量子化と画素単位線型補間の近似を行い、高速フーリエ変換（ＦＦＴ）により線型フィルタリングを実行する高速化ＢＬＦアルゴリズムが記載されている。
【００１６】
非特許文献１４には、レンジ軸のみにおける量子化及び補間により、Ｏ（１）線型フィルタリングを使用することにより高速化した、高速化ＢＬＦアルゴリズム（Yang近似アルゴリズム）が記載されている。
【００１７】
非特許文献１５には、三次元ＦＦＴを用いた、ＢＬＦのより一般的なダウンサイジング・アプローチが記載されている。
【００１８】
非特許文献１６では、ＢＬＦをバイラテラル空間におけるガウス変換として捉え、高速ガウス変換（ＦＧＴ）により計算する手法が開示されている。
【００１９】
非特許文献１−３では、ガウシアン・レンジ・カーネルを三角函数と多項式函数で近似することにより高速化した、Ｏ（１）ＢＬＦが開示されている。この方法によれば、バイラテラル空間の量子化や補間などによって生じる並進変化問題（translation-variant problem）を解決することができる。このアイデアの基本は、函数の並進可能性（shiftability）（幾つかの基底函数の和で近似することにより、函数を任意量だけ平行移動しても函数形が不変であるという性質）にある（後述の式（６）参照）。
【００２０】
非特許文献１−３においては、ガウシアン・レンジ・カーネルｇ（σ，ｘ）を、次式（３ａ）又は（３ｂ）を用いて、余弦函数の冪乗又は多項式函数の冪乗により近似する（以下「Chaudhury近似」という）。
【００２１】
【数３】
2016122430.tif 000005
【００２２】
式（３ｂ）はネピア数の定義をそのまま適用すれば導かれる。式（３ａ）は、コサインをTaylor展開しＮ＞＞１の条件の下で式（３ｂ）に帰着させることによって導かれる。式（３ａ）におけるＮ＞Ｎ０の条件は、ｘの定義域［−T/2，T/2］内で近似式が正であり且つ単調であることを保証するための条件である。
【００２３】
また、非特許文献２では、式（３ａ）の収束をさらに速めた改良式として次の近似式（４）が記載されている。
【００２４】
【数４】
2016122430.tif 000006
【００２５】
さらに、非特許文献４，５では、他のガウシアン・カーネル近似式としては、ガウシアン・カーネルｇ（σ，ｘ）を、次式（５）により、余弦函数の多項式と指数函数との積和により近似する方法が記載されている（以下「Deriche近似」という）。
【００２６】
【数５】
2016122430.tif 000007
【００２７】
式（５）は、定義域ｘ∈［０，∞］におけるフィッティング式であり、ａ０，ａ１，ｂ０，ｂ１，ｃ０，ｃ１，ｗ０，ｗ１はフィッティング定数である。
【００２８】
ところで、ガウシアン・カーネル近似演算手法をバイラテラル・フィルタに適用する場合、定時間フィルタ（計算量がフィルタ窓のサイズＭに依存しないフィルタ）とするためには、ガウシアン・カーネルの近似式は並進可能（shiftable）でなければならない（非特許文献３）。ここで、函数φ（ｘ）（ｘ∈Ｒｄ）が並進可能であるとは、任意の並進移動τ（∈Ｒｄ）に対して函数φ（ｘ−τ）が有限個のある定められた基底函数φ１（ｚ），…，φＮ（ｚ）の線型和
【００２９】
【数６】
2016122430.tif 000008
により表されることをいう。ガウシアン・カーネルの近似式が並進可能であれば、バイラテラル・フィルタ（１ａ）の分子・分母の和分は、移動和（規格化していない移動平均）の式に帰着させることができ、周知の通り再帰を用いた計算によって定時間又はＯ（１）アルゴリズムによって計算が可能となるからである（非特許文献３）。
【００３０】
上記式（３ａ），（４），（５）の近似式は、いずれも並進可能であり、バイラテラル・フィルタに適用した場合には定時間フィルタとして構成することが可能である。
【先行技術文献】
【非特許文献】
【００３１】
【非特許文献１】K. N. Chaudhury, D. Sage, and M. Unser, "Fast O(1) bilateral filtering using trigonometric range kernels", IEEE Transactions on Image Processing, vol.20(12), pp.3376-3382, 2011.
【非特許文献２】K. N. Chaudhury, "Acceleration of the shiftable O(1) algorithm for bilateral filtering and non-local means", IEEE Transactions on Image Processing, vol. 22(4), pp. 1291- 1300, 2013.
【非特許文献３】K. N. Chaudhury, "Constant-time filtering using shiftable kernels", IEEE Signal Processing Letters, Vol.18, no.11, pp.651-654, 2011.
【非特許文献４】Rachid Deriche, "Recursively implementing the Gaussian and its derivatives", Technical report, INRIA, 1993.
【非特許文献５】Tan, Jason L. Dale, and Alan Johnston, ”Performance of three recursive algorithms for fast space-variant Gaussian filtering”. Real-Time Imaging 9, pp.251-228, 2003.
【非特許文献６】P. Perona, J. Malik, “Scale-space and edge detection using anisotropic diffusion”, IEEE Trans. on Pattern Recognition and Machine Intelligence, vol. 12, no. 7, pp. 629-639 , 1990.
【非特許文献７】R.L. Lagendijk, J. Biemond, and D.E. Boekee, "Regularized iterative image restoration with ringing reduction", IEEE Transactions on Acoustics, Vol.36, no.12, pp.1874-1888, 1998.
【非特許文献８】G. Petschnigg, et al., "Digital Photography with Flash and No-Flash Image Pairs", ACM Trans. on Graph. (Proc. of ACM SIGGRAPH 2004), Vol.23 no.3, pp.664-672, 2004.
【非特許文献９】E. Eisemann and F.Durand, "Flash Photography Enhancement via Intrinsic Relighting", ACM Trans. on Graph. (Proc. of ACM SIGGRAPH 2004), Vol.23, no.3, pp.673-678, 2004.
【非特許文献１０】B. Weiss, "Fast median and bilateral filtering", ACM Trans. on Graph. (Proc. of ACM SIGGRAPH 2006), Vol.25 no.3, pp.519-526, 2006.
【非特許文献１１】F. Porikli, "Constant time O(1) bilateral filtering", IEEE Conference on Computer Vision and Pattern Recognition 2008, pp.1-8, 2008.
【非特許文献１２】F. Porikli, "Integral histogram: a fast way to extract histograms in Cartesian spaces", IEEE Computer Soc. Conf. on Computer Vision and Pattern Recognition 2005, Vol.1, pp.829-836, 2005.
【非特許文献１３】F. Durand and J. Dorsey, "Fast Bilateral Filtering for the Display of High-Dynamic-Range Images", ACM Trans. on Graph. (Proc. of ACM SIGGRAPH 2002), Vol.21, no.3, pp.257-266, 2002.
【非特許文献１４】Q. Yang, K.-H. Tan, N. Ahuja, "Real-Time O(1) Bilateral Filtering", IEEE Conf. on Comput. Vis. and Pattern Recognit. 2009 (CVPR 2009), pp.557-564, 2009.
【非特許文献１５】S. Paris and F. Durand, "A Fast Approximation of the Bilateral Filter using a Signal Processing Approach", Int. J. of Comput. Vis. (IJCV), Vol.81 no.1, pp.24-52, 2009.
【非特許文献１６】S. Yoshizawa, A. Belyaev, and H. Yokota, "Fast Gauss Bilateral Filtering", Computer Graphics Forum, Vol.29, no.1, pp.60-74, 2010.
【非特許文献１７】W. Ray and R. Driver, "Further decomposition of the Karhunen-Loeve series representation of a stationary random process", IEEE Transactions on Information Theory, Vol.16 , no.6, pp.663-668, 1970.
【非特許文献１８】K. Sugimoto and S. Kamata, "Fast Gaussian filter with second-order shift property of DCT-5", Proc. IEEE Int. Coef. Image Process. (ICIP), pp.514-518, 2013.
【非特許文献１９】K. Sugimoto and S. Kamata, "Fast Gaussian filter using DCT-V", Proc. IEEE Int. Coef. Image Process. (ICIP), pp.514-518, 2013.
【非特許文献２０】R. Deriche, "Fast algorithms for low-level vision", IEEE 9th Int. Conf. on Pattern Recognition, Vol.1, pp.434-438, 1988.
【非特許文献２１】I. T. Young, L. J. Vliet, "Recursive implementation of the Gaussian filter", Signal Processing, Vol.44, no.2, pp.139-151, 1995.
【非特許文献２２】L. J. Vliet, I. T. Young and P. W. Verbeek, "Recursive Gaussian Derivative Filters",  IEEE Comp. Soc. Press, Vol.1, pp.509-514, 1998.
【非特許文献２３】G. Farnemack and C. F. Westin, "Improving Deriche-style Recursive Gaussian Filters", J. Math. Imaging Vol.26, pp.293-299, 2006.
【発明の概要】
【発明が解決しようとする課題】
【００３２】
非特許文献１５，１６に示されているように、Ｄ次元ＢＬＦは、新たな強度軸ｔを導入することによって、次のように（Ｄ＋１）次元コンボリューションの比に展開することができる。
【００３３】
【数７】
2016122430.tif 000009
ここで、δ（ｘ）はクロネッカ・デルタである。
【００３４】
座標ｐ∈ＲＤと強度ｔ∈Ｒとで張られた（Ｄ＋１）次元空間を「バイラテラル空間（bilateral space）」という。最初にｘ（ｐ）からξ（ｐ，ｔ），ξ０（ｐ，ｔ）が生成され、これがバイラテラル空間内の点とみなされる。ξ（ｐ，ｔ），ξ０（ｐ，ｔ）から、それらの模糊点群（blurred point clouds）β（ｐ，ｔ），β０（ｐ，ｔ）（次式参照）が式（７ａ）の分子，分母として計算される。
【００３５】
【数８】
2016122430.tif 000010
【００３６】
これらの操作は、（Ｄ＋１）次元のガウシアン・コンボリューションとみなすことができる。なぜなら、［ｐ，ｔ］Ｔ∈ＲＤ＋１は（Ｄ＋１）次元空間における対象画素の位置、［ｐ，ｘ（ｑ）］Ｔ∈ＲＤ＋１は（Ｄ＋１）次元空間における対象画素の近傍位置の一つとみなすことができるからである。ＢＬＦの出力は、これらの比から次のように計算される。
【００３７】
【数９】
2016122430.tif 000011
【００３８】
式（９ａ）のζ（ｐ，ｔ）はバイラテラル空間におけるＢＬＦ処理画像、式（９ｂ）はバイラテラル空間からもとの座標空間への投射（縮退）である。図１に、Ｄ次元ＢＬＦを（Ｄ＋１）次元のガウシアン・コンボリューションの比のＤ次元空間への投射として捉えた場合におけるＢＬＦ計算過程の概念図を示す。
【００３９】
以上のように、Ｄ次元ＢＬＦは、２つの（Ｄ＋１）次元ガウシアン・コンボリューションの比として計算される。しかしながら、この次元上昇によって計算コストが大幅に増加することになる。上述した従来のＢＬＦの高速化アルゴリズムにおいては、様々な手法によりβ（ｐ，ｔ），β０（ｐ，ｔ），ζ（ｐ，ｔ）を簡単化し、この計算困難性を減少させたものと捉えることができる。
【００４０】
非特許文献１３，１４，１５に記載の高速化ＢＬＦでは、バイラテラル空間における量子化と部分化線型補間（piece-wise linear interpolation）により、β（ｐ，ｔ），β０（ｐ，ｔ），ζ（ｐ，ｔ）のサイズを積極的に減少させている。量子化によって計算コストは大幅に減少するが、その過度な単純化によって近似精度が減少する。さらに、その精度は、空間軸方向又は強度軸方向のシフトに対してロバストではない（並進変化問題（translation-variant problem））。なぜなら、その精度は量子化閾値に強く依存するからである。
【００４１】
一方、非特許文献１−３では、β（ｐ，ｔ），β０（ｐ，ｔ）を並進可能カーネル（shiftable kernel）により近似した。即ち、ガウシアン・レンジ・カーネルを、三角函数カーネル（基底）や多項式カーネル（基底）のような並進可能カーネルにより置き換えた（式（３ａ），（４）参照）。この考え方により、並進変化問題が解決されるとともに、許容誤差のパラメータを用いて近似精度を制御するための合理的な方法が与えられる。しかしながら、実際の応用において許容できる程度の近似精度を得るためには、カーネルの演算に於いて、依然かなり大きな反復演算（iteration）回数が必要となる。
【００４２】
実際、Chaudhury近似（式（３ａ），（４））によりガウシアン・カーネル演算を行う場合、σの大きさにより収束速度が変化する（例えば、非特許文献２，Table.I, Table.III参照）。特に、σが小さくなるほど収束速度は遅くなり、計算打ち切り数Ｎが大きくなる。例えば、式（３ａ）の場合、ｘのレンジを［−T/2，T/2］とすると、近似式がレンジ内で正であり且つ単調であることを保証するためには、原則としてＮは少なくともＮ０＝４Ｔ２／π２σ２＝０．４０５（Ｔ／σ）２よりも大きくなければならない（非特許文献２）（尚、非特許文献２，３では、ガウシアン・カーネルは｜ｘ｜＞３σで殆ど０になることを利用してＮ０の条件をさらに緩和している）。一般にグレースケール画像ではＴ＝２５５が使用され、このときσ＝８ではＮ０＝４１３，σ＝５ではＮ０＝１０５３，σ＝３ではＮ０＝２９２９とされている。式（４）の近似では、式（３ａ）よりも収束は速いが、非特許文献２，Table.IIIによれば、σ＝８ではＮ０＝８８，σ＝５ではＮ０＝９２，σ＝３ではＮ０＝９６とされており、σが小さいと依然としてＮが大きくなり、１画素のＢＬＦ演算において１００オーダーの反復回数が必要となる。従って、依然として大きな計算コストが必要とされるという課題がある。エッジ保存性能を高めるには狭幅のレンジ・カーネルが必要とされるため、この問題の解決は、特にオンライン画像処理にＢＬＦを適用する場合などにおいて重要である。
【００４３】
また、Deriche近似（式（５））も並進可能であるため、式（３ａ），（４）の代わりに式（５）を使用しても並進変化問題を克服できる。しかしながら、Deriche近似は、ガウシアン・カーネルのフィッティング式であり、近似誤差が比較的大きく精度が悪い。また、σによって近似誤差が変化してしまうため、近似誤差の制御ができないという問題がある。
【００４４】
非特許文献１６に記載の高速化ＢＬＦでは、ξ（ｐ，ｔ），ξ０（ｐ，ｔ）を点群（point clouds）として取り扱い、Ｄ次元ガウシアン・コンボリューションの代わりに、点群に対して高速ガウス変換を行うことによりβ（ｐ，ｔ），β０（ｐ，ｔ）を求めている。しかし、このアルゴリズムは、σｓが小さくなると演算時間が長くなる傾向がある。
【００４５】
また、上記非特許文献１０，１１記載の高速化ＢＬＦは、ヒストグラムを使用するためにその管理コストが大きくなるという問題がある。
【００４６】
以上のように、上述した従来の高速化ＢＬＦは、特にオンライン・ビデオ処理への応用やさらなる高次元／高解像度画像処理への応用において、計算コスト及び精度の両方に於いて十分ではない。従って、Ｏ（１）ＢＬＦにおいて、計算コストと近似精度の間の権衡をさらに最適なものとすることが強く求められている。また、ＢＬＦの用途によって、要求される近似精度は様々であることから、計算誤差の制御を自由に行うことができることも重要である。
【００４７】
そこで、本発明の目的は、従来よりも計算コストが小さく、かつ計算誤差の制御を自在に行うことが可能なガウシアン・カーネル演算装置及びそれを用いた画像フィルタ演算装置を提供することにある。
【課題を解決するための手段】
【００４８】
式（１ａ），（２）から明らかな通り、ＢＬＦとクロス／ジョイントＢＬＦとは、定義式において些細な差違しかないので、以下、本明細書においては代表としてＢＬＦについてのみ説明する。以下に述べる本発明における解析やアプローチは、クロス／ジョイントＢＬＦに対しても全く同様に適用することが出来る。以下、〔１〕〜〔５〕で本発明における考え方と基本原理について説明し、〔６〕において本発明の構成及び作用について説明する。尚、以下の説明では、σは与えられているものとして、特に明示する必要のない限り、ガウシアン・カーネルｇ（σ，ｘ）はｇ（ｘ）と略記する。
【００４９】
〔１〕圧縮アプローチ
本発明においては、画像圧縮の観点から、上述のβ（ｐ，ｔ），β０（ｐ，ｔ），ζ（ｐ，ｔ）を効率的に圧縮することにより、計算コストの低減化を図ることを考える（以下これを「圧縮アプローチ」と呼ぶ）。図２は、バイラテラル空間へ展開した標準画像（“lenna”）に於いて、ｔ軸に沿って一定の間隔でスライスしたβ（ｐ，ｔ），β０（ｐ，ｔ），ζ（ｐ，ｔ）のレイヤ画像である。図２において、グレースケールの原画像のダイナミックレンジはｘ（ｐ）∈［０，１］に規格化している。各系列のレイヤ画像は、ｔに対して非常に滑らかに変化し、隣接するレイヤ画像間では高い相関を示していることがわかる。非可逆圧縮などでは、かかる冗長なデータは、僅かな精度犠牲のもとで、線型変換近似によって圧縮することができることが知られている。
【００５０】
そこで、図３に示すように、β（ｐ，ｔ）を非可逆圧縮することを考える。まず、β（ｐ，ｔ）をｔ軸に沿った線型変換Ｆｔにより変換し、次いで、変換領域に於いて幾つかの有意でないスペクトル成分を除去し、最後に、逆変換Ｆｔ−１により圧縮したβ（ｐ，ｔ）であるβ＾（ｐ，ｔ）を再構成する。β０（ｐ，ｔ），ζ（ｐ，ｔ）も同様にして圧縮できる。このようにして圧縮されたＢＬＦを、以下「圧縮バイラテラル・フィルタ（compressive bilateral filtering : CBLF）」と呼ぶ。今、この圧縮過程をＣｔと記す。Ｃｔはｔのみを変数として使用するので、圧縮過程Ｃｔは、次のようなガウシアン・カーネルの圧縮過程に帰着する。
【００５１】
【数１０】
2016122430.tif 000012
【００５２】
Ｃｔにより圧縮されたガウシアン・カーネルは、次のような線型変換Ｆｔの基底函数の和として与えられる。
【００５３】
【数１１】
2016122430.tif 000013
ここで、Ｎは有意なスペクトル成分の数、ｗｎ（・），φｎ（・）は、それぞれ、ｎ番目の変換係数，Ｆｔの基底函数である。式（１１）を式（１０）に代入すれば次のようになる。
【００５４】
【数１２】
2016122430.tif 000014
ここで、Ωｎ（ｐ）はｎ番目の構成画像である。
【００５５】
Ωｎ（ｐ）は、積画像φｎ（ｘ（ｑ））ｘ（ｑ）のガウシアン・コンボリューションと解される。従って、この圧縮アプローチは、ＢＬＦをＮ個のガウシアン・コンボリューションに分解する。Ｎはコンボリューション演算の反復回数に対応する。Ｎを小さくすれば、計算コストは減少するが、近似精度は劣化する。
【００５６】
〔２〕フーリエ基底によるガウシアン・カーネルの圧縮
画像圧縮において、ｌ２誤差を最小にするには、線型変換ＦｔとしてKarhunen-Loeve変換（ＫＬＴ）を用いればよいことが知られている。ＫＬＴでは、様々なｔに対してシフトされたｇｒ（ｘ（ｑ）−ｔ）の組から最適な基底函数φｎ（ｘ（ｑ））を求める。しかし、ＫＬＴでは、一般に導出される基底函数φｎ（ｘ（ｑ））及び変換係数ｗｎ（ｔ）が数学的に複雑な形となり、また、固有値問題を数値的に解くため計算コストが過剰に必要になるという問題がある。
【００５７】
ＫＬＴの適切な代替手段としては、フーリエ変換の族が考えられる。画像圧縮の分野においては、ターゲットが高い相関性を持つ定常１階マルコフ信号の場合、ＫＬＴは離散コサイン変換により近似できることが広く知られている（非特許文献１７）。そこで、本発明に於いては、ｔに対してシフトされたガウシアン・カーネルの組｛ｇｒ（ｘ（ｑ）−ｔ）｝が、高い相関性を持つ定常１階マルコフ信号であると仮定して（図２参照）、ガウシアン・レンジ・カーネルを、フーリエ級数（ＦＳ）を用いて圧縮する。
【００５８】
以下では、ＦＳにより圧縮されたガウシアン・レンジ・カーネルを考える。尚、ガウシアン・レンジ・カーネルであることが明らかである場合にはｇｒ（・）の添字ｒを省略し単にｇ（・）と記す。
【００５９】
定義域ｘ∈［−Ｒ，Ｒ］におけるガウシアン・カーネルｇ（ｘ）は、フーリエ基底を用いて、次のように圧縮することができる。
【００６０】
【数１３】
2016122430.tif 000015
ここで、Ｋ∈Ｎ（Ｎは自然数全体の集合）は、近似に使用する低周波成分の数（以下「打切項数（truncation term number）」という。）である。ａｋはフーリエ係数である。Ｔ／２はｇ（ｘ）が有意な値を採る範囲を表し、以下Ｔを「有意周期」と呼ぶ。離散ガウシアン・カーネルｇ（ｘ）が偶函数なので、cos項のみを考えればよい。従って、式（１１）は次のようになる。
【００６１】
【数１４】
2016122430.tif 000016
【００６２】
上式から、この圧縮アプローチでは、ガウシアン・コンボリューションにＮ＝２Ｋ＋２回の反復演算と１回の和演算が必要である。従って、ＢＬＦ演算を高速化するためには、Ｋをできる限り小さくすればよいことが分かる。
【００６３】
〔３〕周期的ガウシアン・カーネルの最適化
式（１３ａ）の近似により、定義域ｘ∈［−T/2，T/2］におけるガウシアン・カーネルｇ（ｘ）は、定義域ｘ∈（−∞，∞）の周期函数ｇ＾（ｘ；Ｋ，Ｔ）により近似される。以下、ｇ＾（ｘ；Ｋ，Ｔ）を「周期的ガウシアン・カーネル（periodic Gaussian kernel）」と呼ぶ。フーリエ変換による圧縮アプローチの圧縮レートを明らかにするため、まずａｋを次のような閉形式で近似する。
【００６４】
【数１５】
2016122430.tif 000017
ここで、｜ｘ｜＞Ｔ／２においてｇ（ｘ）≒０とした。
【００６５】
ａｋは指数函数的に減少することから、式（１３ａ）から、Ｋをある程度小さくしても十分な精度で近似することができ、σが与えられれば、ｋに対するａｋの減衰速度を決める制御因子はＴのみであることが分かる。Ｔが小さいほどａｋの減衰は速くなる。
【００６６】
図４に、３つの異なる周期条件における周期的ガウシアン・カーネルを示す。図４において、入力画像のダイナミックレンジを［０，Ｒ］としている。一般に、グレースケール画像ではＲ＝２５５，二値画像ではＲ＝１が広く用いられている。式（１１）に対してレンジ・カーネルがサポートしなければならない定義域は［−Ｒ，Ｒ］である。図４の各図の中央の両矢印は、各周期的ガウシアン・カーネルの中央の周期を表す。周期的ガウシアン・カーネルの有意周期ＴがＴ／２＝Ｒのとき、定義域全体に於いて周期的ガウシアン・カーネルは正確なガウシアン形状となる。ガウシアンは値がほぼ０の長いテールを持つため、有意周期ＴがＴ／２＝０．７Ｒのときも、周期的ガウシアン・カーネルはまだ定義域全域に於いてガウシアン形状を維持している。しかし、Ｔ／２＝０．５Ｒまで有意周期Ｔを縮めると、定義域［−Ｒ，Ｒ］の左右に隣接する周期のガウシアン・プロファイルが定義域［−Ｒ，Ｒ］内に侵入してくるため、周期的ガウシアン・カーネルはガウシアン形状から大きく崩れる。これは、式（１３ａ）において打切項数Ｋが与えられると、近似誤差の観点から最適な有意周期Ｔが存在することを示唆している。
【００６７】
打切項数Ｋと有意周期Ｔの最適な組合せ（Ｋ，Ｔ）を求めるため、まず、評価函数である近似誤差Ｅ（Ｋ，Ｔ）を次式の「カーネル誤差（kernel error）」により定義する。
【００６８】
【数１６】
2016122430.tif 000018
【００６９】
許容誤差τ（０＜τ＜＜１）が与えられると、組（Ｋ，Ｔ）は拘束条件Ｅ（Ｋ，Ｔ）≦τ２により拘束される。Ｅ（Ｋ，Ｔ）は、Ｋ軸に対して単調な函数であることから、拘束条件Ｅ（Ｋ，Ｔ）≦τ２を満たす最小のＫはＴの函数Ｋ（Ｔ）となる。すなわち、周期的ガウシアン・カーネルの最適化問題は、次のように定式化することができる。
【００７０】
【数１７】
2016122430.tif 000019
ここで、Ｎは自然数の集合、（Ｋτ，Ｔτ）は（Ｋ，Ｔ）の最適値を表す。
【００７１】
以下、式（１７ａ），（１７ｂ）を「周期的ガウシアン・カーネルの最適化」と呼ぶ。式（１７ａ），（１７ｂ）により、近似誤差を自在に制御しつつ（Ｋ，Ｔ）を合理的に決定する方法が与えられる。
【００７２】
またその一方、周期的ガウシアン・カーネルの最適化は、Ｏ（１）ＢＬＦの新たな手法でもある。従来の並進可能カーネルによるアプローチ（非特許文献１−３）ではＴ＝Ｒとした場合についてのみ注目しているのに対して、本発明の周期的ガウシアン・カーネルの最適化によるアプローチでは、図４のようにＴ≦Ｒの場合にまで一般化しているからである。これは、本発明における新しい着眼点である。従来の並進可能カーネルによるアプローチでは、実際にMAXFILTERと呼ばれる事前走査によってターゲット画像の有効なダイナミックレンジを事前に測定し、ダイナミックレンジＲ自体を減少しようとしている。しかしながら、この方法は、低コントラスト画像においてのみ有効であり、ターゲット画像のコンテンツに依存する。このコンテンツ依存性があるため、最悪ケース性能を改善することができないという問題がある。それに対して、本発明における周期的ガウシアン・カーネルの最適化によるアプローチでは、ガウシアン・カーネルの形状的特徴を十分に活用することにより、最悪ケース性能を確実に向上させることができる。
【００７３】
〔４〕近似誤差の評価
次に、式（１７ａ），（１７ｂ）の周期的ガウシアン・カーネルの最適化問題の具体的な解法について説明する。Ｔ＝Ｒの場合には、式（１３ａ）の直交性を利用して式（１７ａ）を解析的に解くことができるが、Ｔ＜Ｒの場合にはそれができない。また、式（１７ａ）を反復法により解くとすると、式（１６）の評価に大きな計算コストが必要とされるため、高速化の観点から問題がある。そこで、本発明においては、幾つかの適切な仮定のもとで、式（１６）の値を次の閉形式の近似式により推定する。
【００７４】
【数１８】
2016122430.tif 000020
ここで、erfc(・)は相補誤差函数である。
【００７５】
式（１８）のＥ＾（Ｋ，Ｔ）を「概算カーネル誤差（estimated kernel error）」と呼ぶ。式（１８）は式（１６）の右辺を区間［−T/2，T/2］の部分（式（１８）の右辺第１項）とそれ以外の部分（式（１８）の右辺第２項）とに分解して、上記仮定の下で近似したものである。尚、式（１８）の詳しい証明については、本発明の実施とは直接関係ないため省略する。
【００７６】
本発明に於いては、式（１６）の代わりに式（１８）を用いて近似誤差量の評価を行う（即ち、Ｅ（Ｋ，Ｔ）≒Ｅ＾（Ｋ，Ｔ）とする）。
【００７７】
式（１８）において、右辺の各項は非負である。従って、拘束条件Ｅ（Ｋ，Ｔ）≦τ２を満たすための必要条件は、erfc(πσ(2K+1)/T)∈(0, τ2]，erfc((T-R)/σ)∈[0,τ2]である。従って、Ｔの上限と下限が次のように定められる。
【００７８】
【数１９】
2016122430.tif 000021
【００７９】
式（１９ａ）より、σζτ＋R≦（πσ(2K+1)）／ζτであり、Ｋが整数であることを考慮すると、拘束条件Ｅ２（Ｋ，Ｔ）≦τ２を満たす最小のＫは次のように決定される。
【００８０】
【数２０】
2016122430.tif 000022
【００８１】
有意周期Ｔは式（１９ａ）の範囲内であることから、この範囲内に於いてエネルギー誤差Ｅ＾（Ｋ，Ｔ）を最小化するＴがＴτである。即ち、
【００８２】
【数２１】
2016122430.tif 000023
【００８３】
式（１８）より、式（２１）の目的関数Ｅ＾（Ｋ，Ｔ）の最適化問題は、次の方程式を解く問題に帰着する。
【００８４】
【数２２】
2016122430.tif 000024
【００８５】
上記式（２２）の解Ｔτは、二分法によって数値的に解くことができる。
【００８６】
図５に、Ｒ＝１，σ＝０．１の場合における幾つかのＫに対するカーネル誤差Ｅ（Ｋ，Ｔ）と概算カーネル誤差Ｅ＾（Ｋ，Ｔ）との計算値の比較を示す。図５において、横軸は周期Ｔ、縦軸はカーネル誤差Ｅ（Ｋ，Ｔ）又は概算カーネル誤差Ｅ＾（Ｋ，Ｔ）を示す。図５より、Ｋが十分に大きいとき、Ｅ（Ｋ，Ｔ）≒Ｅ＾（Ｋ，Ｔ）が充足されることが分かる。Ｒ＝１，σ＝０．１の条件では、許容誤差がτ２＝０．００５のときは、最小のＫは一定値Ｋ＝４となる。
【００８７】
〔５〕簡略化した周期的ガウシアン・カーネルの準最適化
最適な（Ｋ，Ｔ）の決定は、〔４〕において説明した手法を用いることによって実施することができるが、ＢＬＦの演算における多少の計算コストを犠牲にして準最適な（Ｋ，Ｔ）を求めるのであれば、以下に説明する、より簡略化した手法を採ることもできる。
【００８８】
図４より、周期的ガウシアン・カーネルｇ＾（ｘ；Ｋ，Ｔ）は、Ｋが十分に大きい場合には、閉区間［−T/2，T/2］のガウシアン・カーネルｇ（ｘ）をｘ軸方向に２ｎＴ（ｎは整数）だけ平行移動したピーク函数ｇｎ（ｘ）の和で表されることが分かる。従って、ｇ＾（ｘ；Ｋ，Ｔ）はHeavisideの階段函数Ｈ０（ｘ）を用いて次のように表すことができる。
【００８９】
【数２３】
2016122430.tif 000025
ここで、δ（ξ）はDiracのデルタ函数である。
【００９０】
Ｋが十分に大きい場合、Ｔ＝Ｒのときは、定義域［−Ｒ，Ｒ］においてｇ＾（ｘ；Ｋ，Ｔ）はｇ（ｘ）に一致する。Ｔ＜Ｒのときは、定義域［−Ｒ，Ｒ］の両端（ｘ＝±Ｒ）の近傍に於いて、ｇ０（ｘ）の隣接周期のピーク函数ｇ±１（ｘ）の影響を受けて誤差が生じる。従って、次のように、ｘ＝±Ｒにおいてｇ±１（±Ｒ）≒０を補償するようにＴを決定すればよい。
【００９１】
【数２４】
2016122430.tif 000026
【００９２】
ｘ≧ασにおいてｇ（ｘ）≒０が補償されるとすると、隣接周期のピーク函数ｇ±１（ｘ）の影響を受けないように補償するには、Ｔ−Ｒ＝ασであればよい。従って、有意周期Ｔの最適値Ｔτは、簡易的に次のように決定することができる。
【００９３】
【数２５】
2016122430.tif 000027
ここで、αは有意水準（関数値を有意とみなす割合）を表すパラメータである。例えば、有意水準を９９．７３％とすればα＝３となる。
【００９４】
式（２５）により有意周期Ｔの最適値Ｔτが決定されれば、上述の式（１７ｂ），（１９ａ）より打切項数Ｋの最適値Ｋτを決定することができる。例えば、最適値Ｋτを次式により計算できる。
【００９５】
【数２６】
2016122430.tif 000028
この方法は、ＴとＫとを同時に最適化しないため、一般に、〔４〕で説明した方法に比べて算出されるＫτが大きくなる傾向があり、算出される（Ｋτ，Ｔτ）は準最適値である。しかし、〔４〕に比べＴτを簡易に決定できるので、一つの有効な手法である。
【００９６】
〔６〕本発明の構成
本発明に係る画像フィルタ演算装置の第１の構成は、入力画像{x(p)}（pは画素の位置，x(p)は位置pの画素値）に対して、標準偏差（以下「スケール・パラメータ」という。）σｓのガウス関数で定義される空間カーネルgs(q-p)（p，qは画素の位置）、及びスケール・パラメータσｒのガウス関数で定義されるレンジ・カーネルgr(x(q)-x(p))とのコンボリューション演算により画像の平滑化を行う画像フィルタ演算装置であって、
前記スケール・パラメータσｒ、所定の許容誤差τ及び前記入力画像{x(p)}のダイナミックレンジＲに対し、前記レンジ・カーネルgr(x(q)-x(p))を、有限項数のフーリエ級数により近似した近似レンジ・カーネルg^r(x(q)-x(p))における級数の打切項数Ｋ、及び前記近似レンジ・カーネルg^r(x(q)-x(p))における各項の展開係数ａｋ（ｋ＝０，…，Ｋ）を演算する際のフーリエ積分の積分領域［-T/2，T/2］を定める有意周期Ｔの組である近似パラメータ（Ｔ，Ｋ）を、前記近似レンジ・カーネルg^r(x(q)-x(p))の近似誤差が前記許容誤差τ以下となるように決定し出力する近似パラメータ最適化手段と、
前記近似パラメータ最適化手段が算出する前記近似パラメータ（Ｔ，Ｋ）に従い、前記近似レンジ・カーネルg^r(x(q)-x(p))の第０項から第Ｋ項までの展開係数{a0,…,aＫ}を算出する展開係数演算手段と、
前記展開係数演算手段が算出する展開係数{a0,…,aＫ}を用いて、前記入力画像{x(p)}に対し、前記空間カーネルgs(q-p)及び前記近似レンジ・カーネルg^r(x(q)-x(p))とのコンボリューション演算を実行するフィルタ演算手段と、を備えたことを特徴とする。
【００９７】
この構成によれば、近似パラメータ最適化手段により、許容誤差τに応じて、打切項数Ｋのみならず、有意周期Ｔをも変更して最適化することで、従来よりも計算コストが小さく、かつ計算誤差の制御を自在に行うことが可能となる。
【００９８】
本発明に係る画像フィルタ演算装置の第２の構成は、前記第１の構成に於いて、前記近似パラメータ最適化手段は、前記打切項数Ｋを算出する打切項数演算手段を備え、
前記打切項数演算手段は、前記スケール・パラメータσｒ、所定の許容誤差τ及び前記入力画像{x(p)}のダイナミックレンジＲに対し、前記打切項数Ｋを
【００９９】
【数２７】
2016122430.tif 000029
により算出することを特徴とする。
【０１００】
本発明に係る画像フィルタ演算装置の第３の構成は、前記第１又は２の構成に於いて、前記近似パラメータ最適化手段は、前記有意周期Ｔを算出する有意周期最適化手段を備え、
前記展開係数演算手段は、前記レンジ・カーネルgr(x(q)-x(p))を積分領域［-T/2，T/2］でフーリエ積分した近似展開係数を、前記各展開係数{a0,…,aＫ}として算出するものであり、
前記有意周期最適化手段は、前記スケール・パラメータσｒ、前記打切項数Ｋ及び前記ダイナミックレンジＲに対する、変数Ｔの方程式
【０１０１】
【数２８】
2016122430.tif 000030
を満たすＴを前記有意周期Ｔとして算出するものであることを特徴とする。
【０１０２】
本発明に係る画像フィルタ演算装置の第４の構成は、前記第１の構成に於いて、前記近似パラメータ最適化手段は、
前記スケール・パラメータσｒ及び前記レンジ上限値Ｒ、並びに予め設定された有効幅パラメータα（ασｒ＜Ｒ）に基づき、前記有意周期Ｔを、
【０１０３】
【数２９】
2016122430.tif 000031
により算出する演算周期設定手段と、
前記許容誤差τ，前記有意周期Ｔ及び前記スケール・パラメータσｒに基づいて、前記打切項数Ｋを
【０１０４】
【数３０】
2016122430.tif 000032
により算出する打切項数決定手段と、を備えたことを特徴とする。
【０１０５】
本発明に係る画像フィルタ演算装置の第５の構成は、前記第１乃至４の何れか一の構成に於いて、前記展開係数演算手段は、前記各展開係数{a1,…,aＫ}を、
【０１０６】
【数３１】
2016122430.tif 000033
により算出するものであることを特徴とする。
【０１０７】
本発明に係る画像フィルタ演算プログラムは、コンピュータに読み込ませて実行することにより、前記コンピュータを前記第１乃至５の何れか一の構成の画像フィルタ演算装置として機能させることを特徴とする。
【０１０８】
本発明に係るガウシアン・カーネル演算装置の第１の構成は、標準偏差（以下「スケール・パラメータ」という。）σｒのガウス関数で定義されるガウシアン・カーネルgr(x)の近似値を算出するガウシアン・カーネル演算装置であって、
前記スケール・パラメータσｒ、所定の許容誤差τ及び前記入力画像{x(p)}のダイナミックレンジＲに対し、前記レンジ・カーネルgr(x(q)-x(p))を、有限項数のフーリエ級数により近似した近似レンジ・カーネルg^r(x(q)-x(p))における級数の打切項数Ｋ、及び前記近似レンジ・カーネルg^r(x(q)-x(p))における各項の展開係数ａｋ（ｋ＝０，…，Ｋ）を演算する際のフーリエ積分の積分領域［-T/2，T/2］を定める有意周期Ｔの組である近似パラメータ（Ｔ，Ｋ）を、前記近似レンジ・カーネルg^r(x(q)-x(p))の近似誤差が前記許容誤差τ以下となるように決定し出力する近似パラメータ最適化手段と、
前記近似パラメータ最適化手段が算出する前記近似パラメータ（Ｔ，Ｋ）に従い、前記近似ガウシアン・カーネルgr(x)の第０項から第Ｋ項までの展開係数{a0,…,aＫ}を算出する展開係数演算手段と、
前記展開係数演算手段が算出する展開係数{a0,…,aＫ}を用いて、前記近似ガウシアン・カーネルgr(x)を算出する近似関数演算手段と、を備えたことを特徴とする。
【０１０９】
本発明に係るガウシアン・カーネル演算装置の第２の構成は、前記第１の構成に於いて、前記近似パラメータ最適化手段は、
前記打切項数Ｋを算出する打切項数演算手段を備え、
前記打切項数演算手段は、前記スケール・パラメータσｒ、所定の許容誤差τ及び前記入力画像{x(p)}のダイナミックレンジＲに対し、前記打切項数Ｋを
【０１１０】
【数３２】
2016122430.tif 000034
により算出することを特徴とする。
【０１１１】
本発明に係るガウシアン・カーネル演算装置の第３の構成は、前記第１又は２の構成に於いて、前記近似パラメータ最適化手段は、前記有意周期Ｔを算出する有意周期最適化手段を備え、
前記展開係数演算手段は、前記近似ガウシアン・カーネルgr(x)を積分領域［-T/2，T/2］でフーリエ積分した近似展開係数を、前記各展開係数{a0,…,aＫ}として算出するものであり、
前記有意周期最適化手段は、前記スケール・パラメータσｒ、前記打切項数Ｋ及び前記ダイナミックレンジＲに対する、変数Ｔの方程式
【０１１２】
【数３３】
2016122430.tif 000035
を満たすＴを前記有意周期Ｔとして算出するものであることを特徴とする。
【０１１３】
本発明に係るガウシアン・カーネル演算装置の第４の構成は、前記第１の構成に於いて、前記近似パラメータ最適化手段は、
前記スケール・パラメータσｒ及び前記レンジ上限値Ｒ、並びに予め設定された有効幅パラメータα（ασｒ＜Ｒ）に基づき、前記有意周期Ｔを、
【０１１４】
【数３４】
2016122430.tif 000036
により算出する演算周期設定手段と、
前記許容誤差τ，前記有意周期Ｔ及び前記スケール・パラメータσｒに基づいて、前記打切項数Ｋを
【０１１５】
【数３５】
2016122430.tif 000037
により算出する打切項数決定手段と、を備えたことを特徴とする。
【０１１６】
本発明に係るガウシアン・カーネル演算装置の第５の構成は、前記第１乃至４の何れか一の構成に於いて、前記展開係数演算手段は、前記各展開係数{a1,…,aＫ}を、
【０１１７】
【数３６】
2016122430.tif 000038
により算出するものであることを特徴とする。
【０１１８】
本発明に係るガウシアン・カーネル演算プログラムは、コンピュータに読み込ませて実行することにより、前記コンピュータを第１乃至５の何れか一の構成のガウシアン・カーネル演算装置として機能させることを特徴とする。
【発明の効果】
【０１１９】
以上のように、本発明に係る画像フィルタ演算装置及びガウシアン・カーネル演算装置によれば、従来よりも計算コストが小さく、かつ計算誤差の制御を自在に行うことが可能となる。
【０１２０】
尚、本発明の画像フィルタ演算装置及びガウシアン・カーネル演算装置は、主としてＢＬＦへの適用を前提としているが、応用範囲はＢＬＦに限られず、上述したクロス／ジョイントＢＬＦなどにも適用可能である。
【図面の簡単な説明】
【０１２１】
【図１】Ｄ次元ＢＬＦを（Ｄ＋１）次元のガウシアン・コンボリューションの比のＤ次元空間への投射として捉えた場合におけるＢＬＦ計算過程の概念図である。
【図２】バイラテラル空間に於いて、ｔ軸に沿って一定の間隔でスライスしたβ（ｐ，ｔ），β０（ｐ，ｔ），ζ（ｐ，ｔ）のレイヤ画像である。
【図３】ｔ軸に沿った線型変換アプローチによる圧縮バイラテラル・フィルタの構成手順を示す図である。
【図４】３つの異なる周期条件における周期的ガウシアン・カーネルを示す図である。
【図５】Ｒ＝１，σ＝０．１の場合における幾つかのＫに対するカーネル誤差Ｅ（Ｋ，Ｔ）と概算カーネル誤差Ｅ＾（Ｋ，Ｔ）との計算値の比較図である。
【図６】実施例１に係るガウシアン・カーネル演算装置のハードウェア構成を表す図である。
【図７】実施例１に係るガウシアン・カーネル演算装置の機能構成を表すブロック図である。
【図８】実施例１のガウシアン・カーネル演算装置１０における、最適な近似パラメータ（Ｋτ，Ｔτ）及びに有意係数ａｋの演算処理の全体の流れを表すＰＡＤである。
【図９】実施例１のガウシアン・カーネル演算装置１０における、有意周期Ｔτの最適値演算処理を表すＰＡＤである。
【図１０】級数展開近似したガウシアン・カーネルの打切項数Ｋと近似誤差Ｅとの関係を示す図である。
【図１１】σｒとτとの関係を示す図である。
【図１２】実施例２に係るフィルタ演算装置の機能構成を表すブロック図である。
【図１３】実施例２のフィルタ演算装置２０における、バイラテラル演算処理の流れを表すＰＡＤである。
【図１４】ＧＦ演算器２９による２次元離散ガウシアン・フィルタ演算処理を表すＰＡＤである。
【図１５】非特許文献２のChaudhury近似アルゴリズム（ε＝0.8）及びCBLF（τ＝0.1）における近似精度[dB]と計算時間を比較した図である。
【図１６】非特許文献１４のYang近似アルゴリズム（Ｂ＝８）及びCBLFにおける近似精度[dB]と計算時間を比較した図である。
【図１７】原画像と、オリジナルのＢＬＦ、Chaudhury近似アルゴリズム、Yang近似アルゴリズム、及びＣＢＬＦの各アルゴリズムによって生成した平滑化画像と、１０倍に増幅した誤差画像である。
【図１８】実施例３に係るフィルタ演算装置の全体構成を示すブロック図である。
【図１９】図１８の演算パラメータメモリ５１及び近似パラメータ最適化部５２の構成を表すブロック図である。
【図２０】図１８のフィルタ演算回路５４の構成を表すブロック図である。
【図２１】図２０のＢ回路７８（Ｂ回路７８−１〜７８−Ｋ）の構成を表すブロック図である。
【発明を実施するための形態】
【０１２２】
以下、本発明を実施するための形態について、図面を参照しながら説明する。
【実施例１】
【０１２３】
図６は、本発明に係るガウシアン・カーネル演算装置のハードウェア構成を表す図である。
【０１２４】
図６において、本発明のガウシアン・カーネル演算装置は、コンピュータ１、入力装置２、出力装置３及び外部記憶装置４を備えている。入力装置２は、キーボードやマウス、ＣＤ−ＲＯＭ等の一般のコンピュータの入力装置である。出力装置３は、プリンタ、ディスプレイ等の一般のコンピュータの出力装置である。外部記憶装置４は、磁気ディスク等の一般の外部記憶装置である。本実施例に係るガウシアン・カーネル演算装置は、提供されるガウシアン・カーネル演算プログラムを、入力装置２からコンピュータ１に読み込ませて実行することにより実現される。
【０１２５】
図７は、本発明に係るガウシアン・カーネル演算装置１０の機能構成を表すブロック図である。図７の機能構成は、コンピュータ１において前記ガウシアン・カーネル演算プログラムを実行させることによって実現される。
【０１２６】
演算パラメータメモリ１１は、入力装置２から予め入力されるガウシアン・カーネルの演算パラメータσ，τ，Ｒなどを記憶する。ここで、σはガウシアン・カーネルのスケール・パラメータ、τは許容近似誤差、Ｒはサンプル値の採り得る最大値（レンジ上限値）である。
【０１２７】
近似パラメータ最適化手段１２は、上述した周期的ガウシアン・カーネルｇ＾（ｘ；Ｋ，Ｔ）の打切項数Ｋ及び有意周期Ｔを決定する。近似パラメータ最適化手段１２は、打切項数演算器１３，有意周期最適化器１４，及びＫＴレジスタ１４ａを備えている。打切項数演算器１３は、式（２０）により、許容誤差τに対して最小となる打切項数Ｋτを算出する。有意周期最適化器１４は、打切項数Ｋτに対して式（１８）の近似誤差を最小とする有意周期Ｔτの値を算出する。ＫＴレジスタ１４ａは、打切項数演算器１３及び有意周期最適化器１４が算出する最適化された近似パラメータ（Ｋτ，Ｔτ）を保持する。
【０１２８】
有意係数演算器１５は、近似パラメータ最適化手段１２により算出されＫＴレジスタ１４ａに保持された打切項数Ｋ及び有意周期Ｔの最適値（Ｋτ，Ｔτ）を用いて、式（１３ｂ）により周期的ガウシアン・カーネルｇ＾（ｘ；Ｋτ，Ｔτ）のフーリエ係数（以下「有意係数」という。）ａｋ（ｋ＝０，１，…，Ｋ）を算出する。有意係数メモリ１６は、有意係数演算器１５により算出された有意係数ａｋ（ｋ＝０，１，…，Ｋτ）を記憶する。
【０１２９】
近似関数演算器１７は、入力装置２から入力されるサンプル値ｘ及び有意係数メモリ１６に記憶された有意係数ａｋ（ｋ＝０，１，…，Ｋ）を用いて、式（１３ａ）により周期的ガウシアン・カーネルｇ＾（ｘ；Ｋ，Ｔ）の値を算出し、出力装置３へ出力する。
【０１３０】
以上のように構成された本実施例に係るガウシアン・カーネル演算装置１０について、以下その動作を説明する。本実施例に係るガウシアン・カーネル演算装置１０は、最初に入力装置２から演算パラメータ（Ｒ，σ，τ）が入力され、これらが演算パラメータメモリ１１に設定されると、まず、近似パラメータ最適化手段１２により打切項数Ｋ及び有意周期Ｔの最適値（Ｋτ，Ｔτ）が決定され、次に有意係数演算器１５により有意係数ａｋ（ｋ＝０，１，…，Ｋ）が算出され、これら有意係数ａｋが有意係数メモリ１６に保存される。その後、入力装置２からサンプル値ｘが入力される度に、近似関数演算器１７が、予め算出された打切項数Ｋ及び有意係数ａｋを用いて、周期的ガウシアン・カーネルｇ＾（ｘ；Ｋτ，Ｔτ）の値を式（１３ａ）により計算し、これを出力装置３へ出力する。
【０１３１】
〔１〕打切項数Ｋ及び有意係数ａｋの演算
図８，図９は、実施例１のガウシアン・カーネル演算装置１０における打切項数Ｋ及び有意周期Ｔ並びに有意係数ａｋの演算処理を表すフローチャートである。ここでは、予め入力装置２から演算パラメータ（Ｒ，σ，τ）が入力され、これらが演算パラメータメモリ１１に設定されているものとする。
【０１３２】
〔１．１〕全体の流れ
図８は、最適な近似パラメータ（Ｋτ，Ｔτ）及びに有意係数ａｋの演算処理の全体の流れを表すＰＡＤ（Problem Analysis Diagram）（JIS X 0128:1988／ISO 8631）である。ここでは、最初に式（１８）の近似誤差Ｅ＾が許容誤差τとなり且つＫが最小となるように（Ｋτ，Ｔτ）を設定し、その後、有意係数（フーリエ係数）ａｋの演算を行う。
【０１３３】
（ステップＳ１０１）
まず、打切項数演算器１３は、演算パラメータメモリ１１に設定された各パラメータ（Ｒ，σ，τ）を読み出し、式（２０）により、許容誤差τに対する最適な打切項数Ｋτを算出し、ＫＴレジスタ１４ａへ保存する。
【０１３４】
（ステップＳ１０２）
次に、有意周期最適化器１４は、演算パラメータメモリ１１に設定された各パラメータ（Ｒ，σ，τ）を読み出し、式（１９ａ）の閉区間において、Ｔの方程式（２２）を解くことにより、許容誤差τに対する最適な有意周期Ｔτを算出し、ＫＴレジスタ１４ａへ保存する。尚、方程式（２２）は解析的に解くことはできないため、反復法によって数値的に解く必要がある。しかし、Ｔの採り得る範囲は式（１９ａ）によって狭い範囲に制約されているため、二分法を用いて算出するのが好適である。この有意周期Ｔτの演算処理の具体例については後述する。
【０１３５】
（ステップＳ１０３）
次に、有意係数演算器１５は、０番目のフーリエ係数ａ０（＝１／Ｔτ）を算出し、有意係数メモリ１６へ保存する。
【０１３６】
（ステップＳ１０３−Ｓ１０４）
次に、有意係数演算器１５は、先にＫＴレジスタ１４ａに保持されている最適値（Ｋτ，Ｔτ）より、１番目〜Ｋτ番目のフーリエ係数ａｋ（ｋ＝１，…，Ｋτ）の近似値を式（１５）により算出し、有意係数メモリ１６へ保存する。
【０１３７】
以上のステップＳ１０１〜Ｓ１０４により、有意係数メモリ１６には、ガウシアン・カーネルの演算に必要なフーリエ係数ａｋ及び打切項数Ｋτが保存された状態となる。
【０１３８】
〔１．２〕有意周期Ｔτの演算処理
次に、図８のステップＳ１０２における有意周期Ｔτの演算処理の例について説明する。ここでは、Ｔの探索範囲が式（１９ａ）により狭い範囲に制約されていることから、二分法アルゴリズムを用いて最適値Ｔτの探索を行う例を示す。図９は、ステップＳ１０２における有意周期Ｔτの演算処理を表すＰＡＤである。
【０１３９】
（ステップＳ２０１）
まず、有意周期最適化器１４は、Ｔの探索範囲の初期値［Ｔ１，Ｔ２］，最大反復回数imax，方程式（２２）の根の精度εｘの設定を行う。探索範囲の初期値［Ｔ１，Ｔ２］は、式（１９ａ）により算出される。また、最大反復回数imaxは収束を考慮して適宜設定すればよいが、iter＝１２程度に設定しておけば十分である。
【０１４０】
（ステップＳ２０２）
次に、有意周期最適化器１４は、Ｔ＝Ｔ１及びＴ＝Ｔ２において、目的関数ｆ（Ｔ；Ｋτ，Ｒ，σ）の最初の評価を行い、ｚ１＝ｆ（Ｔ１；Ｋτ，Ｒ，σ），ｚmid＝ｆ（Ｔ２；Ｋτ，Ｒ，σ）に設定する。ここで、目的関数ｆ（Ｔ）は、方程式（２２）の左辺である。
【０１４１】
（ステップＳ２０３〜Ｓ２０５）
次に、有意周期最適化器１４は、探索方向の設定を行う。即ち、ｚ１＜０の場合には（Ｓ２０３）、探索方向ｄＴを正方向（Ｔ２−Ｔ１）とし、探索基準点ＴrtbをＴ１に設定する（Ｓ２０４）。ｚ１≧０の場合には（Ｓ２０３）、探索方向ｄＴを負方向（Ｔ１−Ｔ２）とし、探索基準点ＴrtbをＴ２に設定する（Ｓ２０５）。
【０１４２】
（ステップＳ２０６〜Ｓ２１１）
次に、有意周期最適化器１４は、二分探索によって方程式（２２）の解を算出する。即ち、以下Ｓ２０７〜Ｓ２１１の処理を反復実行する。
【０１４３】
まず、ｄＴを１／２倍し、現在の探索範囲の中間点ＴmidをＴrtb＋ｄＴにより算出し、中間点Ｔmidにおける目的関数値ｆ（Ｔmid；Ｋτ，Ｒ，σ）を算出しこれをｚmidとする（Ｓ２０７）。次いで、ｚmidが０以下の場合には探索基準点Ｔrtbを中間点Ｔmidにシフトさせる（Ｓ２０８，Ｓ２０９）。そして、｜ｄＴ｜が精度εｘ以下に達したか又は中間点Ｔmidの目的関数値ｚmidが０となった場合には、探索基準点Ｔrtbを解（最適値Ｔτ）として出力し（Ｓ２１０，Ｓ２１１）、そうでなければステップＳ２０７へ戻る。
【０１４４】
以上の処理を行うことによって、方程式（２２）の根である最適な有意周期Ｔτを容易に算出することができる。
【０１４５】
なお、本実施例に於いては、近似パラメータ最適化手段１２が式（１８），（１９ａ）に準じて（Ｋτ，Ｔτ）を決定する例について説明したが、より簡略化した方法として、近似パラメータ最適化手段１２が式（２２）によりＴτを算出し、このＴτを用いて、打切項数演算器１３が式（１９ａ）（具体的には図８のＫ（Ｔ）演算処理）により最適な打切項数Ｋτを算出するように構成することもできる。
【０１４６】
最後に、本実施例に係るガウシアン・カーネル演算装置１０の計算コストについて、従来のガウシアン・カーネル演算装置等と比較して説明する。図１０は、級数展開近似したガウシアン・カーネルの打切項数Ｋと近似誤差Ｅとの関係を示す図である。図１０において、（ａ），（ｂ），（ｃ）は、それぞれ、スケール・パラメータσが０．０５，０．１，０．２の場合を示している。尚、ダイナミックレンジＲは１に規格化している。横軸はスペクトル成分の打切項数（有意とみなすスペクトル成分の数）Ｎ、縦軸は近似誤差の平方根√Ｅ（式（１６））である。（ａ）〜（ｃ）において、「Chaudhury'13」は、背景技術で挙げたChaudhury近似の式（４）により計算した近似ガウシアン・カーネルの近似誤差を表す（尚、この場合、横軸の打切項数Ｋは式（４）の右辺の和の項数Ｎである）。「KLT」は、ガウシアン・カーネルをカルーネン・レーベ変換（ＫＬＴ）により近似した場合を表し、理論上最も最適とされる近似である。「Ours (T＝2R)」は、実施例１のガウシアン・カーネル演算装置において有意周期ＴをダイナミックレンジＲに固定した場合（即ち、通常のFourier級数近似）の近似誤差を表す。「Ours (T≦2R)」は、実施例１のガウシアン・カーネル演算装置において打切項数Ｋ及び有意周期Ｔの最適値（Ｋτ，Ｔτ）を図９のフローチャートに従った処理によって決定した場合の近似誤差を表す。また、（ａ）〜（ｃ）の各図においてＥ＝０．０５（５％）のラインが、この計算例における許容誤差τである。
【０１４７】
Ｅ＝０．５におけるＫの値を比較すると、明らかに、本実施例に係るガウシアン・カーネル演算装置は、従来のChaudhury近似に比べて少ない打切項数Ｋで所望の近似精度が得られることが分かる。スケール・パラメータσが小さくなるほどその傾向は顕著となる。また、Chaudhury近似や本実施例のガウシアン・カーネル演算装置の近似では、Ｋが大きくなると近似誤差Ｅが一定の値に漸近し、ある一定値（以下「漸近近似誤差値Ｅ∞」という。）よりも小さくならない。σが大きいほど漸近近似誤差値Ｅ∞は大きくなる傾向がある。しかし、図１０（ｃ）から明らかなように、本実施例のガウシアン・カーネル演算装置の近似では、Chaudhury近似よりも漸近近似誤差値Ｅ∞を小さくすることができ、大きなσにおいても精度のよい近似値を得ることができる。
【０１４８】
また、「Ours (T＝2R)」と「Ours (T≦2R)」を比較すると、本実施例のガウシアン・カーネル演算装置において、有意周期ＴをダイナミックレンジＲよりも小さい値に最適化することで、小さい打切項数Ｋで所望の近似精度を得ることができることが分かる。これは、スケール・パラメータσで規定されるガウシアン・カーネルの広がりに対して、有意周期ＴをＲから適度に縮小することにより、周期的ガウシアン・カーネルの形状がコサイン函数の形状に近づき、これにより高周波成分が減少することによる。すなわち、ダイナミックレンジＲがスケール・パラメータσに対して非常に大きい場合、周期的ガウシアン・カーネルの形状は、一定の周期で鋭いピークが並んだものとなり、これを近似するには多くの高周波成分を必要とするが、有意周期ＴをダイナミックレンジＲより小さくし各ピークの間隔を狭くして周期的ガウシアン・カーネルの形状を一定周期の滑らかな波の形状に近づけることで、近似に必要な高周波成分を減少させることができる。従って、本実施例のガウシアン・カーネル演算装置では、打切項数Ｋが小さいため、近似関数演算器１７において式（１３ａ）の演算を行う際に計算量（反復演算回数）を従来よりも減少させることができる。
【０１４９】
さらに、図１０は、レンジ・スケール・パラメータσｒが小さくなるほど、より大きなスペクトル成分の数Ｎ（すなわちＫ）が必要とされることが分かる。これは、σｒが小さくなるほど、圧縮性の側面からみて、冗長性の小さい模糊点群β（ｐ，ｔ），β０（ｐ，ｔ）が生成されることによる。さらに詳しくこの振る舞いを観察するため、図１１にσｒとτとの関係を示す。ガウシアン・レンジ・カーネルの近似効果はσｒ／Ｒとτのみに依存するため、図１１の等高線マップから空間フィルタの数Ｍ＝４Ｋ＋１が分かる。
【実施例２】
【０１５０】
〔１〕演算アルゴリズム
本実施例では、本発明に係るガウシアン・カーネル演算装置をバイラテラル・フィルタへ応用したフィルタ演算装置について説明する。式（１ａ）より、バイラテラル・フィルタは次のように表される。
【０１５１】
【数３７】
2016122430.tif 000039
ここで、ｘ（ｐ）は入力画像，ｚ（ｐ）は出力画像、β（ｐ），β０（ｐ）はそれぞれ分子画像，分母画像である。
【０１５２】
式（３３ｂ），（３３ｃ）に、本発明に係る周期的ガウシアン・カーネル（式（１３ａ））を代入して整理すると、分子画像β（ｐ），分母画像β０（ｐ）は、それぞれ次のように近似される。
【０１５３】
【数３８】
2016122430.tif 000040
【０１５４】
ｚ＾（ｐ）を「圧縮バイラテラル・フィルタ（compressive bilateral filter：CBLF）」、β＾（ｐ）を「近似分子画像」、β＾０（ｐ）を「近似分母画像」と呼ぶ。
【０１５５】
近似分子画像β＾（ｐ），近似分母画像β＾０（ｐ）において、ｑについての和（｛・｝で括られた部分）は、それぞれ、関数cos(πkx(p)/L)x(p)，sin(πkx(p)/L)x(p)，cos(πkx(p)/L)，sin(πkx(p)/L)をガウシアン・フィルタ処理した値となっている。そこで、ガウシアン・フィルタ関数GF（σ，ｐ，｛ｆ｝）（σはスケール・ファクタ，ｐは画素位置ベクトル、｛ｆ｝はスカラー関数ｆの関数値の集合）とおくと、近似分子画像β＾（ｐ），近似分母画像β＾０（ｐ）は次のように表すことができる。
【０１５６】
【数３９】
2016122430.tif 000041
【０１５７】
ガウシアン・フィルタ関数GF（σ，ｐ，｛ｆ｝）の計算については、既に多くのＯ（１）ガウシアン・フィルタ・アルゴリズムが知られているので（例えば、非特許文献１８−２３参照）、それを使用すればよい。以上の結果、本発明に係るガウシアン・カーネル演算手法を使用したバイラテラル・フィルタ演算アルゴリズムは次のようになる。
【０１５８】
【表１】
2016122430.tif 000042
【０１５９】
ここで、cos(・)，sin(・)は要素単位のベクトル演算として定義している。また、
【０１６０】
【数４０】
2016122430.tif 000043
は要素単位の和算，積算，除算を表す。函数 Calc_Significant_Period(・)は、適切な解法により式（２２）を解くことを表す。具体的には、実施例１の図９の処理を使用することができる。函数GF(・)は、空間ガウシアン・フィルタのＯ（１）アルゴリズムを表す。ここでは、非特許文献１８，１９に記載の効率的Ｏ（１）アルゴリズムを用いる。このアルゴリズムは、従来のＯ（１）ガウシアン・フィルタとして広く知られている再帰的ガウシアン・フィルタ（非特許文献２０−２３）に比べておよそ２倍程度速く動作し、広レンジのσｓにわたり安定して動作する。本実施例のCBLFの計算コストは、主に、これら（４Ｋ＋１）個の空間ガウシアン・フィルタの数（（２Ｋ＋１）個のβ（ｐ，ｔ）と２Ｋ個のβ０（ｐ，ｔ））に依存する。
【０１６１】
〔２〕装置構成と動作
本実施例に於いて、フィルタ演算装置のハードウェアの全体構成は図６と同様であるものとし、本実施例に係るフィルタ演算装置は、提供されるフィルタ演算プログラムを、入力装置２からコンピュータ１に読み込ませて実行することにより実現されるものとする。
【０１６２】
尚、以下の説明では、簡単のため、特に断らない限り、上述の近似分子画像β＾（ｐ），近似分母画像β＾０（ｐ）を、省略して、それぞれ単に分子画像β（ｐ），分母画像β０（ｐ）と呼ぶ。
【０１６３】
図１２は、実施例２に係るフィルタ演算装置の機能構成を表すブロック図である。図１２において、実施例１の図７と同様の部分については、同符号を付している。尚、本実施例に係るフィルタ演算装置は、マイコンなどを用いて組み込みボードとしてハードウェア的に構成してもよいが、提供されるフィルタ演算プログラムを、入力装置２からコンピュータ１に読み込ませて実行することにより実現するように構成してもよい。
【０１６４】
本実施例２のフィルタ演算装置２０は、実施例１で説明した、演算パラメータメモリ１１、近似パラメータ最適化手段１２、有意係数演算器１５、有意係数メモリ１６、近似函数演算器１７に加えて、入力画像メモリ２１、出力画像メモリ２２、項数カウンタ２３を備えている。近似パラメータ最適化手段１２、有意係数演算器１５、有意係数メモリ１６については、実施例１と同様なので説明は省略する。
【０１６５】
演算パラメータメモリ１１は、入力装置２から予め入力されるガウシアン・カーネルの演算パラメータσｒ，σｓ，τ，Ｒを記憶する。ここで、σｒ，σｓは空間次元及びレンジ次元のスケール・パラメータ（標準偏差）、τは許容近似誤差、Ｒは画素値の採り得る最大値（レンジ上限値）である。
【０１６６】
入力画像メモリ２１は、入力装置２から入力される入力画像ｘ＝（ｘｉｊ）を記憶する。尚、入力画像ｘはサイズＮ×Ｍ画素の画像データであり、ｘｉｊは位置（ｉ，ｊ）の画素の画素値を表す。各画素の画素値は０〜Ｒまでの値を採る。
【０１６７】
出力画像メモリ２２は、入力画像ｘに対してバイラテラル演算処理を行って得られる出力画像ｚ＝（ｚｉｊ）を記憶する。
【０１６８】
項数カウンタ２３は、近似演算を行うフーリエ級数の項番号ｋをカウントする。実施例１と同様、近似パラメータ最適化手段１２により、最適な打切項数Ｋτが演算され出力されるので、項数カウンタ２３は、項番号ｋを０からＫτまでカウントする。
【０１６９】
近似函数演算器１７は、近似パラメータ最適化手段１２により算出される（Ｋτ，Ｔτ）、及び有意係数演算器１５により算出される有意係数ａｋ（ｋ＝０，１，…，Ｋτ）を用いて、バイラテラル・フィルタ函数の演算を実行する。ここでは、前述したように、式（３１ａ），（３１ｂ）の分子画像β（ｘ），分母画像β０（ｘ）を個別に計算して、最後に式（３０ｃ）により出力画像ｚを算出する。
【０１７０】
近似函数演算器１７は、基底函数演算器２５、基底函数テーブル２６、分母画像演算器２７、分子画像演算器２８、ＧＦ演算器２９、分母画像メモリ３０、加算器３０ａ、分子画像メモリ３１、加算器３１ａ、画像除算器３２を備えている。
【０１７１】
基底函数演算器２５は、項数カウンタ２３が出力する項番号ｋにおいて入力画像ｘの各画素ｐ（ｉ，ｊ）について、式（３４ａ），（３４ｂ）に示したフーリエ展開式の基底函数
【０１７２】
【数４１】
2016122430.tif 000044
の演算を行う。
【０１７３】
基底函数テーブル２６は、基底函数演算器２５が算出する各基底函数の値を、ルックアップ・テーブル（ＬＵＴ）として記憶する。
【０１７４】
分母画像演算器２７は、式（３５ｂ）により、各画素ｐに対して分母画像β０の第ｋ項の画素値β(k)０（ｐ）＝β(k)０ｉｊの演算を行う。分子画像演算器２８は、式（３５ａ）により、各画素ｐに対して分子画像βの第ｋ項の画素値β(k)（ｐ）＝β(k)ｉｊの演算を行う。ＧＦ演算器２９は、式（３４ａ）（３４ｂ）に示したガウシアン・コンボリューション演算ＧＦ（・）を行う。
【０１７５】
分母画像メモリ３０は、分母画像演算器２７が出力する分母画像β０を記憶する。分子画像メモリ３１は、分子画像演算器２８が出力する分子画像βを記憶する。加算器３０ａ，３１ａは、それぞれ、分母画像β０，分子画像βの各画素の画素値β０ｉｊ，βｉｊを分母画像メモリ３０，画像メモリ３１から読み出し、分母画像演算器２７，分子画像演算器２８が出力する第ｋ項の値β(k)０ｉｊ，β(k)ｉｊを加えて分母画像メモリ３０，画像メモリ３１の画素値β０ｉｊ，βｉｊを更新する。
【０１７６】
画像除算器３２は、分母画像演算器２７，分子画像演算器２８による式（３４ｂ）（３４ｃ）の演算が完了した後、式（３３ａ）により出力画像ｚを算出し、出力画像メモリ２２へ保存する。出力画像メモリ２２へ保存された出力画像ｚは、出力装置３へ出力される。
【０１７７】
以上のように構成された本実施例に係るフィルタ演算装置２０について、以下その動作を説明する。図１３は、実施例２のフィルタ演算装置２０における、バイラテラル演算処理の流れを表すＰＡＤである。実施例１と同様、最初に式（１８）の近似誤差Ｅ＾が許容誤差τとなり且つＫが最小となるように（Ｋτ，Ｔτ）を設定し、その後、有意係数（フーリエ係数）ａｋの演算及びフィルタ函数の演算を行う。
【０１７８】
（ステップＳ４０１）
まず、打切項数演算器１３は、演算パラメータメモリ１１に設定された各パラメータ（Ｒ，σ，τ）を読み出し、式（２０）により、許容誤差τに対する最適な打切項数Ｋτを算出し、ＫＴレジスタ１４ａへ保存する。
【０１７９】
（ステップＳ４０２）
次に、有意周期最適化器１４は、演算パラメータメモリ１１に設定された各パラメータ（Ｒ，σ，τ）を読み出し、式（１９ａ）の閉区間において、Ｔの方程式（２２）を解くことにより、許容誤差τに対する最適な有意周期Ｔτを算出し、ＫＴレジスタ１４ａへ保存する。この（Ｋ，Ｔ）最適値探索処理の詳細については、実施例１で説明した通りである。
【０１８０】
（ステップＳ４０３）
次に、項数カウンタ２３は現在の項番号ｋを０に設定し、有意係数演算器１５は、０番目のフーリエ係数ａ０（＝１／２Ｔτ）を算出し、有意係数メモリ１６へ保存する。
【０１８１】
（ステップＳ４０４）
次に、分母画像β０，分子画像βのｋ＝０番目の項（ＤＣ項）の演算を行う。分母画像演算器２７は、分母画像β０のすべての画素値β０ｉｊをａ０に設定し、分母画像メモリ３０へ保存する。また、分子画像演算器２８は、分子画像βの各画素（ｉ，ｊ）に対する画素値βｉｊを、a0GF(σs, (i,j), {xkl})により算出し、分子画像メモリ３１へ保存する。ここで、ガウシアン・コンボリューション演算
【０１８２】
【数４２】
2016122430.tif 000045
は、ＧＦ演算器２９により実行される。
【０１８３】
（ステップＳ４０５）
次に、項数カウンタ２３が項番号ｋを１ずつインクリメントしながら、分母画像β０，分子画像βのｋ（＞０）番目の項（ＡＣ項）の演算を、以下のステップＳ４０６〜Ｓ４１１を順次繰り返すことにより実行する。
【０１８４】
（ステップＳ４０６）
まず、有意係数演算器１５が、ＫＴレジスタ１４ａのＴτ及び演算パラメータメモリ１１のスケール・パラメータσｒを用いて、ｋ番目のフーリエ係数ａｋを式（１５）により算出し、有意係数メモリ１６へ保存する。
【０１８５】
（ステップＳ４０７）
次に、基底函数演算器２５は、ＫＴレジスタ１４ａのＴτ及び入力画像メモリ２１の画素値（xij）を用いて、各画素（ｉ，ｊ）に対する項番号ｋの基底函数cos(kπxij/Ｔτ)，sin(kπxij/Ｔτ)の値を算出し、基底函数テーブル２６へ保存する。
【０１８６】
（ステップＳ４０８−Ｓ４０９）
次に、分母画像の第ｋ項の演算を行う。まず、分母画像演算器２７は、各画素（ｉ，ｊ）について、演算パラメータメモリ１１のスケール・パラメータσｓと基底函数テーブル２６のコサイン基底ｃｉｊ及びサイン基底ｓｉｊを用いて、コンボリューションΩcij＝GF(σｓ, (i,j), {ckl})，Ωsij＝GF(σｓ, (i,j), {skl})の演算を行う（Ｓ４０８）。ここで、ガウシアン・コンボリューション演算
【０１８７】
【数４３】
2016122430.tif 000046
は、ＧＦ演算器２９により実行される。次いで、分母画像演算器２７は、有意係数メモリ１６のフーリエ係数ａｋを用いて、分母画像の第ｋ項の値ａｋ（ｃｉｊΩcij＋ｓｉｊΩsij）を算出し、加算器３０ａはこれを分母画像メモリ３０の画素値β０ｉｊに加算し、分母画像メモリ３０の画素値β０ｉｊを更新する。これを、全画素について実行する。
【０１８８】
（ステップＳ４１０−Ｓ４１１）
次に、分子画像の第ｋ項の演算を行う。まず、分子画像演算器２８は、各画素（ｉ，ｊ）について、演算パラメータメモリ１１のスケール・パラメータσｓと基底函数テーブル２６のコサイン基底ｃｉｊ及びサイン基底ｓｉｊを用いて、コンボリューションΩcij＝GF(σｓ, (i,j), {cklxkl})，Ωsij＝GF(σｓ, (i,j), {sklxkl})の演算を行う（Ｓ４１０）。ここで、ガウシアン・コンボリューション演算
【０１８９】
【数４４】
2016122430.tif 000047
は、ＧＦ演算器２９により実行される。次いで、分子画像演算器２８は、有意係数メモリ１６のフーリエ係数ａｋを用いて、分子画像の第ｋ項の値ａｋ（ｃｉｊΩcij＋ｓｉｊΩsij）を算出し、加算器３１ａはこれを分子画像メモリ３１の画素値βｉｊに加算し、分子画像メモリ３１の画素値βｉｊを更新する。これを、全画素について実行する。
【０１９０】
（ステップＳ４０７）
以上のステップＳ４０５〜Ｓ４１１の処理が終了した後、画像除算器３２は、各画素（ｉ，ｊ）について、分子画像メモリ３１の画素値βｉｊを分母画像メモリ３０の画素値β０ｉｊで除算して出力画像ｚの画素値ｚｉｊを算出し、出力画像メモリ２２へ保存する。全画素について出力画像ｚの画素値が算出されると、出力画像ｚは出力装置３へ出力され、バイラテラル演算処理を終了する。
【０１９１】
以上の動作により、入力画像ｘに対するバイラテラル演算処理がされた出力画像ｚを得ることができる。次に、ＧＦ演算器２９によるガウシアン・コンボリューション演算GF（σ，ｐ，｛ｆ｝）について補足説明する。ガウシアン・コンボリューション演算については、ここではすでに公知のアルゴリズムを使用して行うが、本実施例では、一例として非特許文献１８，１９に記載のものを参考にして構成した例を示す。
【０１９２】
〔２．１〕ＧＦ演算器
ＧＦ演算器２９は、画素当たり定数時間アルゴリズムを用いてガウシアン・コンボリューション演算（ガウシアン・フィルタ演算）を行う。実際のガウシアン・コンボリューション演算では、離散値によって行われるため、以下ではpの連続関数f(p)の離散値をfpのように表記し、式（１ｂ）のガウシアン・カーネルg(σ,p)の離散値をgpと表記する（特に必要のない限りσは省略，pはスカラーｘ又はベクトル（ｘ，ｙ））。また、位置ｐにおけるガウシアン・コンボリューション演算GF(σ, p, {f(q)})はGF[fp]のように略記する。ＧＦ演算器２９は２次元のガウシアン・フィルタ演算を行うが、ガウシアン・カーネルｇ（σ，ｑ）は２次元座標（ｘ，ｙ）に対して変数分離できるので、ガウシアン・フィルタ演算はｘ方向とｙ方向に独立して実行することができる。すなわち、
【０１９３】
【数４５】
2016122430.tif 000048
ここで、GFx[・]，GFy[・]は、それぞれｘ座標，ｙ座標に対するガウシアン・フィルタ演算を表し、Ｗはガウシアン・フィルタ演算を行う窓領域のサイズを表す。通常、窓サイズＷはｇｘが有意な範囲（例えば、有意水準を９９．７３％とすればＷ＝３σ）に設定される。
【０１９４】
そこで、１次元のガウシアン・フィルタ演算について、１画素当たりの定数時間化を考えれば十分である。まず、ガウシアン・カーネルｇｘをDCT-5により展開すると次のようになる。
【０１９５】
【数４６】
2016122430.tif 000049
ここで、Ｇｋはｇｘのｋ番目のＤＣＴ係数、Ｋ１は打切項数である。Ｋ１はＤＣＴ係数Ｇｋが有意な範囲、例えば、有意水準を９９．７３％とすればＫ１＝３／（φσ）に設定すればよい。従って、函数ｆｘと上記ガウシアン・カーネルｇｘのコンボリューションは、次のように表される。
【０１９６】
【数４７】
2016122430.tif 000050
ここで、Ｆｋ（ｘ）は位置ｘにおける入力系列ｆｘのｋ番目の短時間ＤＣＴ係数、Ｃｋ（ｕ）はコサイン・カーネルである。短時間ＤＣＴ係数Ｆｋ（ｘ）については、次の初期値と漸化式によって表すことができ、この漸化式を用いて計算することにより、１画素当たりの式（３５）の計算量はＯ（１）となる。
【０１９７】
【数４８】
2016122430.tif 000051
【０１９８】
特に、ｋ＝０の場合は、漸化式は次のように簡単になる。
【０１９９】
【数４９】
2016122430.tif 000052
【０２００】
実際の計算に於いては、式（３６ｂ）において、Ｃｋ（１），Ｃｋ（Ｒ）における丸めの誤差が伝搬する。そこで、この誤差伝搬を抑圧するため、式（３５ａ）の右辺において、最大の重み係数Ｇ０を括り出す。
【０２０１】
【数５０】
2016122430.tif 000053
【０２０２】
このとき、Ｚｋ（ｘ）についての漸化式は次のようになる。
【０２０３】
【数５１】
2016122430.tif 000054
【０２０４】
以上より、１次元ガウシアン・フィルタ演算は次のようなアルゴリズムにより実行できる。
（１）最初に、Ｇ０，｛γｋＣｋ（ｕ）｜ｋ＝０，…，Ｋ１；ｕ＝−Ｗ，…，Ｗ｝，｛２Ｃｋ（１）｜ｋ＝０，…，Ｋ１｝を予め計算し、ルックアップ・テーブル（ＬＵＴ）として準備する。
（２）そしてまず、式（３６ａ），（３７ｂ）により、初期値としてＦ０（０），Ｚｋ（０），Ｚｋ（１）（ｋ＝１，…，Ｋ１）を計算し、式（３７ａ）により、まず（ｆ＊ｇ）０を計算する。
（３）次に、漸化式（３６ｂ）によりＦ０（ｘ）を更新する。
（４）次に、式（３７ａ）により（ｆ＊ｇ）ｘを計算し、漸化式（３６ｂ）によりＦ０（ｘ）を更新し、漸化式（３８）によりＺｋ（ｘ−１），Ｚｋ（ｘ）（ｋ＝０，…，Ｋ１）を更新するという演算を、ｘ＝１からＮ−１（Ｎは入力シーケンス（ｆ０，ｆ１，…，ｆＮ−１）の数）まで繰り返し実行し、すべてのｘについての（ｆ＊ｇ）ｘを求める。
【０２０５】
尚、ｘ＜０及びｘ≧Ｎにおけるｆｘは、例えば、次のような補外規則によって補充すればよい。
【０２０６】
【数５２】
2016122430.tif 000055
【０２０７】
以上から、ＧＦ演算器２９による２次元離散ガウシアン・フィルタ演算ＧＦ［ｆxy］は、図１４のＰＡＤに示した処理により実行することができる。
【０２０８】
（１）全体の処理（２次元離散ガウシアン・フィルタ処理）の流れ（図１４（ａ））
まず、ＧＦ演算器２９は、２次元座標（ｘ，ｙ）の函数である入力函数｛fxy｝（Ｎ×Ｍ行列として与えられる。）に対し、全座標位置（Ｎ×Ｍ点）において、ｙ座標に対する１次元ガウシアン・フィルタ処理を実行し、フィルタ処理後の函数を｛f(y) xy｝とする（Ｓ５０１）。函数｛f(y) xy｝はＮ×Ｍ行列である。次に、ＧＦ演算器２９は、函数｛f(y) xy｝に対し、全座標位置（Ｎ×Ｍ点）において、ｙ座標に対する１次元ガウシアン・フィルタ処理を実行し（Ｓ５０２）、フィルタ処理後の函数｛f(xy) xy｝（Ｎ×Ｍ行列）を、２次元離散ガウシアン・フィルタ処理後の函数として出力する（Ｓ５０３）。ここで、ステップＳ５０１，Ｓ５０２におけるｘ座標，ｙ座標に対する１次元ガウシアン・フィルタ処理は、図１４（ｂ）に示すように行われる。
【０２０９】
（２）１次元離散ガウシアン・フィルタ処理の流れ（図１４（ｂ））
図１４（ｂ）は、ステップＳ５０１，Ｓ５０２におけるＧＦｙ［・］，ＧＦｘ［・］の処理の詳細である（代表としてＧＦｘを示す）。入力函数｛fx｝は、座標ｘの函数（ｘ＝０，１，…，Ｎ−１）である。また、図１４（ｂ）において、｛  ｝で記された値については、ＧＦ演算器２９により予め計算され、内部メモリにＬＵＴとして格納されているものとする。
【０２１０】
（ステップＳ６０１）  まず、ＧＦ演算器２９は、漸化式（４４ｂ）の初期値Ｆ０（０）を式（４４ａ）により演算し、その結果を変数Ｆ０に格納する。ここで、Ｃ０（ｕ）＝１である。
【０２１１】
（ステップＳ６０２−Ｓ６０３）  次に、ＧＦ演算器２９は、漸化式（４７）の初期値Ｚｋ（０），Ｚｋ（１）（ｋ＝１，…，Ｋ１）を式（４６ｂ），（４４ａ）により演算し、その結果をそれぞれ変数Ｚｋ−，Ｚｋに格納する。
【０２１２】
（ステップＳ６０４）  次に、ＧＦ演算器２９は、現在、変数Ｆ０及び｛Ｚｋ−；ｋ＝１，…，Ｋ１｝に格納された値を用いて、式（４６ａ）により、まずｘ＝０における（ｆ＊ｇ）０を演算する。
【０２１３】
（ステップＳ６０５）  次に、ＧＦ演算器２９は、変数Ｆ０及び入力函数｛fx｝を用いて、漸化式（４５）によりＦ０（１）を演算し、その結果を変数Ｆ０に再格納する。ここで、ｘ＜０に対する関数値ｆｘは、式（４８）の補外規則によって補完される（後に必要となるｘ≧Ｎの関数値ｆｘの場合も同じ）。
【０２１４】
（ステップＳ６０６−Ｓ６１０）  次に、ｘ＝１，…，Ｎ−１について、以下のステップＳ６０６〜Ｓ６１０を、ｘをインクリメントしながら反復実行することにより、（ｆ＊ｇ）ｘを順次演算する。
【０２１５】
まず、ＧＦ演算器２９は、ｘの値を設定し（Ｓ６０６）、現在、変数Ｆ０及び｛Ｚｋ−；ｋ＝１，…，Ｋ１｝に格納された値を用いて、式（４６ａ）により、現在のｘにおける（ｆ＊ｇ）ｘを演算する（Ｓ６０７）。次いで、ＧＦ演算器２９は、現在の変数Ｆ０及び入力函数｛fx｝を用いて、漸化式（４５）によりＦ０（ｘ＋１）を演算し、その結果を変数Ｆ０に再格納するとともに、入力函数｛fx｝を用いて式（４４ｃ）によりΔ（ｘ）を演算する（Ｓ６０８）。最後に、ＧＦ演算器２９は、各々のｋ（ｋ＝１，…，Ｋ１）に対して、｛Ｚｋ−；ｋ＝１，…，Ｋ１｝，｛Ｚｋ；ｋ＝１，…，Ｋ１｝に格納された値を用いて、漸化式（４７）によりＺｋ（ｘ＋１）を演算し、現在Ｚｋに格納されているＺｋ（ｘ）をＺｋ−に再格納するとともに、新たに算出したＺｋ（ｘ＋１）をＺｋに再格納する（Ｓ６０９，Ｓ６１０）。
【０２１６】
（ステップＳ６１１）  以上の処理により算出された｛（ｆ＊ｇ）ｘ；ｘ＝０，…，Ｎ−１｝を、１次元離散ガウシアン・フィルタ処理後の函数として出力する。
【０２１７】
〔２．２〕空間カーネルにボックス・フィルタを用いた例
上述の実施例では、空間カーネルの演算器としてガウシアン・フィルタ（ＧＦ）演算器を使用した例を示したが、本発明に於いては、ＧＦ演算器に代えて、ボックス・フィルタ演算器（ＢＦ演算器）を使用することもできる。
【０２１８】
ボックス空間カーネル（Box spatial kernel）は次のように定義することができる。
【０２１９】
【数５３】
2016122430.tif 000056
【０２２０】
ここで、「‖・‖∞」はｌ∞−ノルムを表し、ベクトルｑの要素のうち最大のものの値をとる。
【０２２１】
上記のボックス空間カーネルを用いて、ＢＦ演算器は、空間領域の任意の関数ｆ（ｐ）とボックス空間カーネルとのコンボリューションとして、次のように定義される。
【０２２２】
【数５４】
2016122430.tif 000057
【０２２３】
ＢＦ演算器を用いる場合、式（３５ａ），（３５ｂ）におけるＧＦ（・）をＢＦ（・）に置換すればよい。
【０２２４】
〔２．３〕実験例
（１）フィルタリングにおける基本的性能
最後に、本実施例２に係るフィルタ演算装置２０のフィルタリングにおける基本的性能について、従来のバイラテラル・フィルタ演算アルゴリズムと比較して説明する。
【０２２５】
本実験では、２次元フィルタリングにおいて、近似精度と実行時間の間の性能のトレードオフを評価した。実験では、サイズが７６８×５１２画素または５１２×７６８画素の２４ＲＧＢ画像を含むコダックフォトＣＤを使用し、２４画像について平均近似精度と平均実行時間とを測定した。実行時間には、メモリ割り当て時間と演算時間が含まれる。対象アルゴリズムの近似精度は、原画像と、原画像から対象アルゴリズム及び本来のＢＬＦにより生成された平滑化画像との間の、ピーク信号対雑音比（ＰＳＮＲ）［dB］によって定量化した。また、±4σs又は±4σrをサポートするガウシアン空間カーネル又はガウシアン・レンジ・カーネルは、本来のＢＬＦにおいて正確な確度を実現するものと仮定する。一般に、約５０ｄＢであれば、十分な近似精度が実現できたものと評価される。
【０２２６】
最初に、ガウシアン空間／レンジ・カーネルを使用したChaudhury近似アルゴリズムに対するCBLFの優位性について確認する。図１５は、非特許文献２のChaudhury近似アルゴリズム（ε＝0.8）及びCBLF（τ＝0.1）における近似精度[dB]と計算時間を比較した図である。図１５（ａ）は、空間スケール・パラメータσｓが２及び６の場合における、レンジ・スケール・パラメータσｒに対する近似精度（ＰＳＮＲ）［dB］の関係を表す。図１５（ｂ）は、空間スケール・パラメータσｓが２及び６の場合における、レンジ・スケール・パラメータσｒに対する１Ｍ画素当たりの計算時間［s/Mpixels］の関係を表す。各図に於いて、「Chaudhury'13」がChaudhury近似アルゴリズムの結果を示し、「Ours」が本発明に係るフィルタ演算装置２０（CBLF）の結果を示す。また、Chaudhury近似アルゴリズムにおいては、許容誤差εは0.8とし、本発明に係るフィルタ演算装置２０（CBLF）においては許容誤差τは0.1とした。実験の結果、本発明に係るフィルタ演算装置２０は、Chaudhury近似アルゴリズムを用いた場合と比べ、同等の近似精度（約５０ｄＢ）において（図８（ａ））、５〜８倍程度計算速度が速いことが分かる（図８（ｂ））。尚、ここでは、本発明に係るフィルタ演算装置２０における空間フィルタの計算に於いては、〔２．１〕で説明した定数時間のＧＦ演算器を使用した。また、これらの結果は、CBLFの実行時間がσｓに依存しないことを実証している。すなわち、１画素当たりの計算量がO(1)であることを示している。本発明者は、様々なσｓに対して徹底的な実験を行い、この傾向について確認した。故に、本発明に係るフィルタ演算装置２０（ＣＢＬＦ）は、理論と実験において、従来のChaudhury近似アルゴリズムよりかなり優れていることが実証された。
【０２２７】
次に、ガウシアン・レンジ・カーネルとボックス空間カーネル（式（４９ａ））を用いた非特許文献１４のYang近似アルゴリズムに対する本発明に係るフィルタ演算装置２０（ＣＢＬＦ）の優位性について確認する。
非特許文献１４によれば、Yang近似アルゴリズムは、多くのケースに於いて、ボックス空間カーネルのサイズＢをＢ＝８とすることにより、十分な精度が実現される。
図１６は、非特許文献１４のYang近似アルゴリズム（Ｂ＝８）及びCBLFにおける近似精度[dB]と計算時間を比較した図である。図１６（ａ）は、空間スケール・パラメータσｓが２及び６の場合における、レンジ・スケール・パラメータσｒに対する近似精度（ＰＳＮＲ）［dB］の関係を表す。図１６（ｂ）は、空間スケール・パラメータσｓが２及び６の場合における、レンジ・スケール・パラメータσｒに対する１Ｍ画素当たりの計算時間［s/Mpixels］の関係を表す。各図に於いて、「Yang+'09（Ｂ＝８）」及び「Yang+'09（Ｂ＝４）」が、それぞれＢ＝８，４に対するYang近似アルゴリズムの結果を示し、「Ours」が本発明に係るフィルタ演算装置２０（CBLF）の結果を示す。
図１６（ａ）より、Yang近似アルゴリズム（Ｂ＝８）に匹敵する許容誤差とするためには、本発明に係るフィルタ演算装置２０（ＣＢＬＦ）においては、σｓ＝２に対してτ＝２、σｓ＝６に対してτ＝０．０７とすれば適当であることが分かる。図１６から明らかなように、ＣＢＬＦでは、近似精度におけるジグザグな変化と、実行時間の段階的な変化が見られる。これは、有意スペクトル項数Ｋの変化に起因して生じる現象である。与えられた許容誤差τに対する有意スペクトル項数Ｋは、図１１において、Ｋの等高線と与えられた許容誤差τの直線（図１１における点線）との交点から決まる。なお、空間フィルタの数Ｍは、Yang近似アルゴリズムではＭ＝２Ｂ、ＣＢＬＦではＭ＝４Ｋ＋１であることを注意しておく。さらに、図１６（ａ）より、Yang近似アルゴリズムにおいては、広範なσｒに対し、Ｂ＝８では十分な近似精度を得ることが可能であるが、Ｂ＝４では十分な近似精度を得ることができないことが分かる。Yang近似アルゴリズムのような量子化ベースのアルゴリズムでは、不十分な近似精度により、事実上のＭの下限が生じる。そして、これにより、σｒが大きくなるほど冗長性を排除することができなくなる。一方、ＣＢＬＦでは、およそ0.15≦σrのときＭを半分程度とすることができる。故に、ＣＢＬＦでは、特にσｒが大きい場合に、Yang近似アルゴリズムよりも性能トレードオフがより優れている（図１６（ｂ））。
【０２２８】
（２）ノイズの視覚的評価
最後に、オリジナルのＢＬＦ、Chaudhury近似アルゴリズム、Yang近似アルゴリズム、およびＣＢＬＦの近似誤差を目視により評価した結果について説明する。試験用の画像としては、標準画像Data-BAse3における256×256画素の３つのRGB画像を使用した。
【０２２９】
図１７は、原画像と、オリジナルのＢＬＦ、Chaudhury近似アルゴリズム、Yang近似アルゴリズム、及びＣＢＬＦの各アルゴリズムによって生成した平滑化画像と、１０倍に増幅した誤差画像である。すべてのＯ（１）ＢＬＦアルゴリズムは、オリジナルのＢＬＦに十分近似する平滑化画像を生成しているように見え、人間の視覚では、これらの差違は微小であるように見える。明らかに、Chaudhury近似アルゴリズムとＣＢＬＦは、殆ど同じ誤差画像を生成している。これは、双方ともコサイン項展開に基づいた同様のフレームワークによっているためである。一方、Yang近似アルゴリズムは、各セグメントにおいてある程度のオフセット誤差を示しているように見える。これらすべてのＯ（１）アルゴリズムにおいて、誤差は、主に鋭いエッジ部分に、高周波雑音として発生する。これらのエッジは視覚的に増幅されるため、視覚的に感じられる実誤差としては、十分小さく微細である。
【実施例３】
【０２３０】
次に、本発明に係るフィルタ演算装置をコプロセッサやＦＰＧＡなどの並列演算回路により構成した実施例について説明する。
【０２３１】
図１８は、実施例３に係るフィルタ演算装置の全体構成を示すブロック図である。図１８（ａ）は全体図、図１８（ｂ）は、入力ラインバッファ及び出力ラインバッファの構成図である。本実施例のフィルタ演算装置５０は、演算パラメータメモリ５１、近似パラメータ最適化部５２、入力ラインバッファアレイ５３、フィルタ演算回路５４、及び出力ラインバッファアレイ５５を備えている。
【０２３２】
演算パラメータメモリ５１は、フィルタ演算に必要な演算パラメータを記憶するメモリである。近似パラメータ最適化部５２は、演算パラメータメモリ５１に記憶された演算パラメータに基づいて、ガウシアン・カーネルの近似演算に必要な各パラメータを演算する回路ブロックである。入力ラインバッファアレイ５３は、原画像から切り出される部分画像が入力され一時的に記憶するメモリである。原画像のサイズをＷ×Ｈ画素（横Ｗ、縦Ｈ）とする。入力ラインバッファアレイ５３は、（２Ｍ＋１）行分の部分画像を記憶することができる（但し、２Ｍ＋１はフィルタの窓サイズであり、演算パラメータにより設定される）。入力ラインバッファアレイ５３は、図１８（ｂ）に示したように、Ｗ個のレジスタアレイを２Ｍ＋１行配列して構成されており、各列のレジスタ列はシフト・レジスタであり、各行毎にレジスタの値を読み出すことができるように構成されている。入力ラインバッファアレイ５３には、原画像の上端行から順に１行ずつ現画像データが入力される。フィルタ演算回路５４は、入力ラインバッファアレイ５３に保持された原画像に対してフィルタ演算を実行し、平滑化画像を出力する回路である。出力ラインバッファアレイ５５は、フィルタ演算回路５４が出力する平滑化画像を一時的に保持するメモリであり、入力ラインバッファアレイ５３と同様に構成されている。
【０２３３】
図１９は、図１８の演算パラメータメモリ５１及び近似パラメータ最適化部５２の構成を表すブロック図である。
【０２３４】
演算パラメータメモリ５１は、許容誤差τ、原画像のダイナミックレンジＲ、レンジ領域のスケール・パラメータσｒ、空間領域のスケール・パラメータσｓ、及びフィルタサイズ・パラメータＭの値を演算パラメータとしてそれぞれ記憶する。これらの各演算パラメータの値は、外部から入力される。
【０２３５】
近似パラメータ最適化部５２は、ζτ演算回路６１、ζτレジスタ６２、Ｋ演算回路６３、Ｋレジスタ６４、Ｄ演算回路６５、Ｄレジスタ６６、Ｔ演算回路６７、Ｔレジスタ６８、ａｋ演算回路６９、ａｋレジスタ７０、ｍ演算器７１、ｍレジスタ７２を備えている。
【０２３６】
ζτ演算回路６１は、演算パラメータメモリ５１に記憶された許容誤差τに基づき、式（１９ｂ）の逆相補誤差関数の演算を行うことによってζτを算出し、ζτレジスタ６２に保存する。ζτ演算回路６１は、二乗演算と逆相補誤差関数演算を行う回路であり、ルックアップ・テーブル（ＬＵＴ）により構成することができる。
【０２３７】
Ｋ演算回路６３は、演算パラメータメモリ５１に記憶されたダイナミックレンジＲ，スケール・パラメータσｒ、及びζτレジスタ６２に記憶されたζτに基づき、式（２０）の演算を行うことにより、打切項数Ｋを算出しＫレジスタ６４に保存する。
【０２３８】
Ｄ演算回路６５は、式（１９ａ）に示した、有意周期Ｔの範囲Ｄ＝［σｒζτ＋Ｒ，πσｒ（２Ｋ＋１）／ζτ］の下限値（σｒζτ＋Ｒ）及び上限値（πσｒ（２Ｋ＋１）／ζτ）を演算し、これらの値をＤレジスタ６６に保存する。
【０２３９】
Ｔ演算回路６７は、Ｄレジスタ６６に保存された有意周期Ｔの範囲Ｄにおいて、式（２２）を二分法により解くことによって有意周期Ｔを算出する回路である。Ｔ演算回路６７が算出する有意周期Ｔは、Ｔレジスタ６８に保存される。
【０２４０】
ａｋ演算回路６９は、Ｋレジスタ６４に保存された打切項数Ｋ、演算パラメータメモリ５１に記憶されたスケール・パラメータσｒ、及びＴレジスタ６８に保存された有意周期Ｔを読み出して、ｋ＝１〜Ｋまで、式（１５）を演算することにより各係数ａｋを演算子、ａｋレジスタ７０に保存する。ここで、Ｋは可変であるため、ａｋレジスタ７０はＫの変化幅に対応して十分な数のレジスタが必要とされるが、図１１より、Ｋの上限は高々４〜６程度である。
【０２４１】
ｍ演算器７１は、演算パラメータメモリ５１に記憶されたフィルタサイズ・パラメータＭから、フィルタの窓サイズｍ＝２Ｍ＋１を演算し、ｍレジスタ７２に保存する。
【０２４２】
図２０は、図１８のフィルタ演算回路５４の構成を表すブロック図である。フィルタ演算回路５４は、空間フィルタ演算器７５、ＤＣ成分ｂレジスタアレイ７６、ＤＣ成分ｂ０レジスタアレイ７７、Ｂ回路７８−ｋ（ｋ＝１〜Ｋ）、ｂレジスタアレイ７９−ｋ（ｋ＝１〜Ｋ）、ｂ０レジスタアレイ８０−ｋ（ｋ＝１〜Ｋ）、加算器８１，８２、及び除算器８３を備えている。
【０２４３】
空間フィルタ演算器７５は、実施例２で説明したＧＦ演算器であり、入力ラインバッファアレイ５３より、ｍ（＝２Ｍ＋１）行分の画素データ（ｘ［−Ｍ］〜ｘ［Ｍ］）を読み出し、中央行（０行目（図１８（ｂ）参照））の各画素に対して、順次式（３８）に示したＤＣ成分の空間フィルタ演算を行い、その結果得られるＷ個の近似分子画像β＾（式（３５ａ））のＤＣ成分（ｋ＝０）の値をＤＣ成分ｂレジスタアレイ７６に保存する。ＤＣ成分ｂレジスタアレイ７６は、１行分（Ｗ個）の近似分子画像β＾（式（３５ａ））のＣＤ成分の値を保持するレジスタである。
【０２４４】
ＤＣ成分ｂ０レジスタアレイ７７は、１行分（Ｗ個）の近似分母画像β０＾（式（３５ｂ））のＤＣ成分（ｋ＝０）の値を保持するレジスタであり、全て１が記憶されている。
【０２４５】
Ｂ回路７８−ｋ（ｋ＝１〜Ｋ）は、式（３５ａ）（３５ｂ）の近似分子画像β＾及び近似分母画像β０＾のｋ番目の周波数成分を演算する回路である。ｂレジスタアレイ７９−ｋ（ｋ＝１〜Ｋ）は、Ｂ回路７８−ｋが算出する１行分（Ｗ個）の近似分子画像β＾のｋ番目の周波数成分を保持するレジスタである。ｂ０レジスタアレイ８０−ｋ（ｋ＝１〜Ｋ）は、Ｂ回路７８−ｋが算出する１行分（Ｗ個）の近似分母画像β０＾のｋ番目の周波数成分を保持するレジスタである。
【０２４６】
加算器８１は、ｂレジスタアレイ７６，７９−１〜７９−Ｋに保持された各画素（Ｗ個）の近似分子画像β＾の周波数成分を、それぞれの画素毎に足しあせて、式（３５ａ）の右辺の和演算を行い、１行分の各画素に対する近似分子画像β＾の値を算出する。
【０２４７】
加算器８２は、ｂ０レジスタアレイ７７，８０−１〜８０−Ｋに保持された各画素（Ｗ個）の近似分母画像β０＾の周波数成分を、それぞれの画素毎に足しあせて、式（３５ｂ）の右辺の和演算を行い、１行分の各画素に対する近似分母画像β０＾の値を算出する。
【０２４８】
除算器８３は、加算器８１が出力する近似分子画像β＾の値を、加算器８２が出力する近似分母画像β０＾の値で、画素ごとに除算を行い、式（３３ａ）で表される最終的な平滑化画像を算出し、出力ラインバッファアレイ５５へ出力する。
【０２４９】
図２１は、図２０のＢ回路７８（Ｂ回路７８−１〜７８−Ｋ）の構成を表すブロック図である。
【０２５０】
Ｂ回路７８は、画素値セレクタ９１、項番号レジスタ９２、位相演算器９３、ｃｋ演算器９４、ｃｋレジスタアレイ９５、乗算器９６、ｃｆｋレジスタアレイ９７、ｓｋ演算器９８、ｓｋレジスタアレイ９９、乗算器１００、ｓｆｋレジスタアレイ１０１、空間フィルタ演算器１０２〜１０５、乗算器１０６，１０７，１１０，１１１、１０９，１１３、加算器１０８，１１２を備えている。
【０２５１】
画素値セレクタ９１は、入力ラインバッファアレイ５３に保持された画素の原画像データから、順次画素ｑを選択し、画素ｑの画素値ｘ（ｑ）を出力する。ここで、入力ラインバッファアレイ５３には、原画像の画素値が１行毎に追加入力されていくが、画素値セレクタ９１は、新たに追加された行（図１８（ｂ）におけるｘ［−Ｍ］の行）の各画素（Ｗ個）を順次選択するものとする。
【０２５２】
項番号レジスタ９２は、Ｂ回路７８の項番号ｋ（ｋ＝１，…，Ｋ）の値を保持するレジスタである。
【０２５３】
位相演算器９３は、画素値セレクタ９１により選択された画素ｑに対する位相値φk(q)＝２πkx(q)／Tを演算し出力する。
【０２５４】
ｃｋ演算器９４は、位相演算器９３が出力する位相値φk(q)のコサイン値ck(q)=cos(φk(q))を演算し出力する。ｃｋレジスタアレイ９５は、ｃｋ演算器９４が出力するコサイン値ck(q)を記憶するレジスタである。ｃｋレジスタアレイ９５は、図１８（ｂ）に示したラインバッファアレイ５３と同様の構造の（２Ｍ＋１）行Ｗ列のシフト・レジスタであり、新たに入力される行のコサイン値ck(q)は、−Ｍ番目の行に追加される。
【０２５５】
乗算器９６は、ｃｋ演算器９４が出力するコサイン値ck(q)と画素値セレクタ９１が出力する画素値x(q)との乗算を行い、乗算値ck(q)x(q)を出力する。ｃｆｋレジスタアレイ９７は、乗算器９６が出力する乗算値ck(q)x(q)を記憶するレジスタである。ｃｆｋレジスタアレイ９７も、図１８（ｂ）に示したラインバッファアレイ５３と同様の構造の（２Ｍ＋１）行Ｗ列のシフト・レジスタであり、新たに入力される行の乗算値ck(q)x(q)は、−Ｍ番目の行に追加される。
【０２５６】
ｓｋ演算器９８は、位相演算器９３が出力する位相値φk(q)のサイン値sk(q)=sin(φk(q))を演算し出力する。ｓｋレジスタアレイ９９は、ｓｋ演算器９８が出力するサイン値sk(q)を記憶するレジスタである。ｓｋレジスタアレイ９９も、図１８（ｂ）に示したラインバッファアレイ５３と同様の構造の（２Ｍ＋１）行Ｗ列のシフト・レジスタであり、新たに入力される行のサイン値sk(q)は、−Ｍ番目の行に追加される。
【０２５７】
乗算器１００は、ｓｋ演算器９８が出力するサイン値sk(q)と画素値セレクタ９１が出力する画素値x(q)との乗算を行い、乗算値sk(q)x(q)を出力する。ｓｆｋレジスタアレイ１０１は、乗算器１００が出力する乗算値sk(q)x(q)を記憶するレジスタである。ｓｆｋレジスタアレイ１０１も、図１８（ｂ）に示したラインバッファアレイ５３と同様の構造の（２Ｍ＋１）行Ｗ列のシフト・レジスタであり、新たに入力される行の乗算値sk(q)x(q)は、−Ｍ番目の行に追加される。
【０２５８】
空間フィルタ演算器１０２は、式（３５ｂ）の右辺第１項のガウシアン・フィルタ関数GF(σs，p，{ck(q)})（ｐは中央の行（図１８（ｂ）の第０行）のＷ個の各画素）の演算を行う空間フィルタである。σsは演算パラメータメモリ５１より入力され、{ck(q)}はｃｋレジスタアレイ９５より入力される。
【０２５９】
空間フィルタ演算器１０３は、式（３５ｂ）の右辺第２項のガウシアン・フィルタ関数GF(σs，p，{sk(q)})（ｐは中央の行（図１８（ｂ）の第０行）のＷ個の各画素）の演算を行う空間フィルタである。σsは演算パラメータメモリ５１より入力され、{sk(q)}はｓｋレジスタアレイ９９より入力される。
【０２６０】
乗算器１０６は、空間フィルタ演算器１０２が出力する値GF(σs，p，{ck(q)})と、ｃｋレジスタアレイ９５から読み出されるコサイン値ck(p)との乗算を行い、式（３５ｂ）の右辺第１項の乗算値ck(p)GF(σs，p，{ck(q)})を出力する。
【０２６１】
乗算器１０７は、空間フィルタ演算器１０３が出力する値GF(σs，p，{sk(q)})と、ｓｋレジスタアレイ９９から読み出されるサイン値sk(p)との乗算を行い、式（３５ｂ）の右辺第２項の乗算値sk(p)GF(σs，p，{sk(q)})を出力する。
【０２６２】
加算器１０８は、乗算器１０６の出力値と乗算器１０７の出力値とを加算し、式（３５ｂ）の右辺括弧内の和の演算を行う。乗算器１０９は、加算器１０８の出力値と、ａｋレジスタ７０から読み出される係数値akとの乗算を行い、各画素ｐに対する式（３５ｂ）の近似分母画像β＾０（ｐ）を算出し、ｂ０レジスタアレイ８０−ｋに保存する。
【０２６３】
空間フィルタ演算器１０４は、式（３５ａ）の右辺第１項のガウシアン・フィルタ関数GF(σs，p，{ck(q)x(q)})（ｐは中央の行（図１８（ｂ）の第０行）のＷ個の各画素）の演算を行う空間フィルタである。σsは演算パラメータメモリ５１より入力され、{ck(q)x(q)}はｃｆｋレジスタアレイ９７より入力される。
【０２６４】
空間フィルタ演算器１０５は、式（３５ａ）の右辺第２項のガウシアン・フィルタ関数GF(σs，p，{sk(q)x(q)})（ｐは中央の行（図１８（ｂ）の第０行）のＷ個の各画素）の演算を行う空間フィルタである。σsは演算パラメータメモリ５１より入力され、{sk(q)x(q)}はｓｆｋレジスタアレイ１０１より入力される。
【０２６５】
乗算器１１０は、空間フィルタ演算器１０４が出力する値GF(σs，p，{ck(q)x(q)})と、ｃｋレジスタアレイ９５から読み出されるコサイン値ck(p)との乗算を行い、式（３５ａ）の右辺第１項の乗算値ck(p)GF(σs，p，{ck(q)x(q)})を出力する。
【０２６６】
乗算器１１１は、空間フィルタ演算器１０５が出力する値GF(σs，p，{sk(q)x(q)})と、ｓｋレジスタアレイ９９から読み出されるサイン値sk(p)との乗算を行い、式（３５ａ）の右辺第２項の乗算値sk(p)GF(σs，p，{sk(q)x(q)})を出力する。
【０２６７】
加算器１１２は、乗算器１１０の出力値と乗算器１１１の出力値とを加算し、式（３５ａ）の右辺括弧内の和の演算を行う。乗算器１１３は、加算器１１２の出力値と、ａｋレジスタ７０から読み出される係数値akとの乗算を行い、各画素ｐに対する式（３５ａ）の近似分子画像β＾（ｐ）を算出し、ｂレジスタアレイ７９−ｋに保存する。
【０２６８】
以上のような構成により、実施例２と同様に、バイラテラル・フィルタの演算を行うことが可能となる。
【符号の説明】
【０２６９】
１  コンピュータ
２  入力装置
３  出力装置
４  外部記憶装置
１０  ガウシアン・カーネル演算装置
１１  演算パラメータメモリ
１２  近似パラメータ最適化手段
１３  打切項数演算器
１４  有意周期最適化器
１４ａ  ＫＴレジスタ
１５  有意係数演算器
１６  有意係数メモリ
１７  近似函数演算器
２０  フィルタ演算装置
２１  入力画像メモリ
２２  出力画像メモリ
２３  項数カウンタ
２５  基底函数演算器
２６  基底函数テーブル
２７  分母画像演算器
２８  分子画像演算器
２９  ＧＦ演算器
３０  分母画像メモリ
３０ａ  加算器
３１  分子画像メモリ
３１ａ  加算器
３２  画像除算器
５０  フィルタ演算装置
５１  演算パラメータメモリ
５２  近似パラメータ最適化部
５３  入力ラインバッファアレイ
５４  フィルタ演算回路
５５  出力ラインバッファアレイ
６１  ζτ演算回路
６２  ζτレジスタ
６３  Ｋ演算回路
６４  Ｋレジスタ
６５  Ｄ演算回路
６６  Ｄレジスタ
６７  Ｔ演算回路
６８  Ｔレジスタ
６９  ａｋ演算回路
７０  ａｋレジスタ
７１  ｍ演算器
７２  ｍレジスタ
７５  空間フィルタ演算器
７６  ＤＣ成分ｂレジスタアレイ
７７  ＤＣ成分ｂ０レジスタアレイ
７８，７８−ｋ（ｋ＝１〜Ｋ）  Ｂ回路
９−ｋ（ｋ＝１〜Ｋ）  ｂレジスタアレイ７
８０−ｋ（ｋ＝１〜Ｋ）  ｂ０レジスタアレイ
８１，８２  加算器
８３  除算器
９１  画素値セレクタ
９２  項番号レジスタ
９３  位相演算器
９４  ｃｋ演算器
９５  ｃｋレジスタアレイ
９６  乗算器
９７  ｃｆｋレジスタアレイ
９８  ｓｋ演算器
９９  ｓｋレジスタアレイ
１００  乗算器
１０１  ｓｆｋレジスタアレイ
１０２〜１０５  空間フィルタ演算器
１０６，１０７，１１０，１１１、１０９，１１３  乗算器
１０８，１１２  加算器
【図６】
2016122430.tif 000068
【図７】
2016122430.tif 000069
【図８】
2016122430.tif 000070
【図９】
2016122430.tif 000071
【図１２】
2016122430.tif 000072
【図１３】
2016122430.tif 000073
【図１４】
2016122430.tif 000074
【図１８】
2016122430.tif 000075
【図１９】
2016122430.tif 000076
【図２０】
2016122430.tif 000077
【図２１】
2016122430.tif 000078
【図１】
2016122430.tif 000079
【図２】
2016122430.tif 000080
【図３】
2016122430.tif 000081
【図４】
2016122430.tif 000082
【図５】
2016122430.tif 000083
【図１０】
2016122430.tif 000084
【図１１】
2016122430.tif 000085
【図１５】
2016122430.tif 000086
【図１６】
2016122430.tif 000087
【図１７】
2016122430.tif 000088