(57)【要約】【課題】任意の信号位置に対し少ないメモリ使用量でフィルタ処理を施すことができる情報処理装置を提供すること。【解決手段】信号に線形フィルタを適用する情報処理装置１００であって、  初期値として与えられた１番目のベクトルを用いると共に、ｎ−１番目（ｎは２以上の整数）のベクトルと行列の積にｎ番目の信号値を加算してｎ番目のベクトルを求めるベクトル算出部１２と、前記ベクトル算出部が算出したｎ個のベクトルから取り出した２つのベクトルのうち、一方のベクトルと前記行列のべき乗との積から、他方のベクトルと前記行列のべき乗との積を減じた差分を線形フィルタの出力値として出力する出力値算出部１４と、を有することを特徴とする。【選択図】図４
2016103089.tif 000002
【特許請求の範囲】
【請求項１】
信号に線形フィルタを適用する情報処理装置であって、
初期値として与えられた１番目のベクトルを用いると共に、ｎ−１番目（ｎは２以上の整数）のベクトルと行列の積にｎ番目の信号値を加算してｎ番目のベクトルを求めるベクトル算出部と、
前記ベクトル算出部が算出したｎ個のベクトルから取り出した２つのベクトルのうち、一方のベクトルと前記行列のべき乗との積から、他方のベクトルと前記行列のべき乗との積を減じた差分を前記信号に対する線形フィルタの出力値として出力する出力値算出部と、を有することを特徴とする情報処理装置。
【請求項２】
前記ベクトル算出部が算出するｎ番目のベクトルの第i要素はｎ−１番目のベクトルの第i-1と等しく、前記ベクトル算出部は１番目のベクトルからｎ番目のベクトルの第一要素のみを記憶手段に記憶しておく請求項１に記載の情報処理装置。
【請求項３】
前記ベクトル算出部は、前記行列として固有値がexp-jω(ｊは虚数単位)の同伴行列を用いて、ｎ番目のベクトルの第一要素Fn,1を次式で算出する請求項２に記載の情報処理装置。
【数３２】
2016103089.tif 000039
【請求項４】
複数の周波数の線形結合で近似された線形フィルタの係数を取得するフィルタ係数取得部と、
前記出力値算出部が前記周波数ごとに算出した前記差分を加重加算して線形フィルタの出力値を算出する加重加算部と、を有する請求項１〜３いずれか１項に記載の情報処理装置。
【請求項５】
二分法で区切られた前記信号の両端を線形フィルタの適用位置に決定する適用位置決定部と、
テンプレートマッチングのためのテンプレートとして線形フィルタを用い、前記適用位置が決定されるごとに求められた複数の線形フィルタの出力値のうち出力値が最も大きい前記適用位置を出力する適用位置出力部と、を有する請求項４に記載の情報処理装置。
【請求項６】
線形フィルタを複数の周波数の線形結合で表されたガウシアンフィルタ又はガボールフィルタとした場合、
前記加重加算部は、ガウシアンフィルタ又はガボールフィルタを線形結合で表す前記周波数ごとに算出された前記差分を加重加算することで線形フィルタの出力値を算出する請求項４に記載の情報処理装置。
【請求項７】
前記信号が離散信号である場合、前記フィルタ係数取得部は複数の周波数の線形結合で近似された所定の関数を取得し、
前記信号から信号値を取得する位置が登録された位置登録情報を参照してフィルタ処理を施す前記信号の位置を決定する適用位置決定部を有し、
前記出力値算出部は前記適用位置決定部が決定した前記位置において前記差分を前記周波数ごとに算出し、
前記加重加算部は、前記周波数ごとに算出された前記差分を加重加算して前記離散信号が前記信号値を有さない前記位置における前記信号値を算出する請求項４に記載の情報処理装置。
【請求項８】
前記信号がＤ次元の信号（Ｄは２以上の整数）である場合、Ｄ−１次元までの信号に対応した線形フィルタの出力値に対して、Ｄ次元の信号に対応した線形フィルタを適用することを特徴とする請求項１に記載の情報処理装置。
【請求項９】
信号に線形フィルタを適用する情報処理装置によって行われる信号処理方法であって、
初期値として与えられた１番目のベクトルを用いると共に、ｎ−１番目（ｎは２以上の整数）のベクトルと行列の積にｎ番目の信号値を加算してｎ番目のベクトルを求めるベクトル算出ステップと、
前記ベクトル算出ステップにより算出されたｎ個のベクトルから取り出した２つのベクトルのうち、一方のベクトルと前記行列のべき乗との積から、他方のベクトルと前記行列のべき乗との積を減じた差分を線形フィルタの出力値として出力する出力値算出ステップと、を有することを特徴とする信号処理方法。
【請求項１０】
信号に線形フィルタを適用する情報処理装置に、
初期値として与えられた１番目のベクトルを用いると共に、ｎ−１番目（ｎは２以上の整数）のベクトルと行列の積にｎ番目の信号値を加算してｎ番目のベクトルを求めるベクトル算出ステップと、
前記ベクトル算出ステップにより算出されたｎ個のベクトルから取り出した２つのベクトルのうち、一方のベクトルと前記行列のべき乗との積から、他方のベクトルと前記行列のべき乗との積を減じた差分を線形フィルタの出力値として出力する出力値算出ステップと、を実行させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、情報処理装置、信号処理方法及びプログラムに関する。
【背景技術】
【０００２】
画像の平滑化、ノイズ除去、又は、濃度勾配の検出などにフィルタが用いられることは少なくない。フィルタを画像処理などの信号処理に適用する際に利用される演算方法の１つとして畳み込み演算が知られている。畳み込み演算は数学的には関数（信号）と関数（フィルタ）の積の積分で表されるが、画像データのようにデジタル信号が演算対象の場合、掛け算の結果を足し合わせることで求められる。
【０００３】
信号数Ｎの信号系列ｆ(f0, f1, …, fN-1)に対するカーネルサイズＫの線形フィルタg(g0, g1, …, gK-1)の信号位置xにおける出力値(g*f)xは次式で定義される。
【０００４】
【数１】
2016103089.tif 000003
右辺が線形フィルタgと信号系列fとの畳み込み演算を表している。計算機などが右辺の計算をそのまま計算すると、Ｋ回の乗算とＫ−１回の加算が必要なので計算量はカーネルサイズＫに比例する。計算量がＫに比例することをＯ（Ｋ）と称する。したがって、全信号位置{x=1, 2, …, N}において、線形フィルタの出力値を求める際の計算量はＯ(ＮＫ)となる（計算量がＮとＫの積に比例する）。このため、カーネルサイズが大きい線形フィルタでは計算量が膨大となることが知られている。
【０００５】
線形フィルタの畳み込み演算を高速化する方法の１つとして、高速フーリエ変換(Fast Fourier Transform, 以下ＦＦＴと称する。)を用いる方法が知られている。空間領域の線形フィルタには等価な周波数領域のフィルタが存在するため、ＦＦＴを用いることで空間領域での畳み込み演算を周波数領域の演算に変換することができる。ＦＦＴを用いた際の全信号位置での線形フィルタの計算量はＯ(ＮlogＮ)である。したがって、logＮ < Ｋ であれば、ＦＦＴを利用することで畳み込み演算をより高速に実行することが可能となる。
【０００６】
これに対し、畳み込み演算の計算量をＯ（Ｎ）程度に低減する方法が考案されている（例えば、非特許文献１、２参照）。非特許文献１には、矩形領域の画素値の累積値をテーブルとして保持することにより、カーネルサイズに非依存の計算量で畳み込み演算が可能な線形フィルタについて開示されている。非特許文献２には、矩形領域の画素値の累積値をインテグラルイメージと称して、インテグラルイメージを利用したリアルタイムの物体認識について開示されている。
【０００７】
図１を用いて非特許文献１，２における計算量の低減について説明する。図１（ａ）は積分計算の計算量を説明する図の一例である。f(x)を関数としf(x)とｘ軸の間のｘ＝ｂからｘ＝ａまでの面積Ａは「A＝f(b)＋f(b+1)＋f(b+2)＋…＋f(a)」で表すことができる（１，２…はサンプリング間隔の増分を表す）。また、ｘ＝０からｘ＝ａまでの面積をF(a)、ｘ＝０からｘ＝ｂまでの面積をF(b)とすれば面積Ａは「Ａ＝F(a)−F(b)」で表すことができる。したがって、F(a)とF(b)が予め計算されていれば、１回の減算で積分計算と同等の結果を得ることができる。
【０００８】
次に図１（ｂ）を用いてインテグラルイメージについて説明する。インテグラルイメージは上記一次元の計算を二次元で行うものである。点１のインテグラルイメージは領域Ａの画素値の合計であり、点２のインテグラルイメージは領域ＡとＢの画素値の合計であり、点３のインテグラルイメージは領域ＡとＣの画素値の合計であり、点４のインテグラルイメージは領域Ａ、Ｂ、Ｃ、Ｄの画素値の合計である。したがって、点１〜点４のインテグラルイメージが求められていれば、領域Ａ、Ｂ、Ｃ、Ｄのうち任意の領域の画素値の合計を求めることができる。例えば、領域Ｄの画素値の合計は、インテグラルイメージを使って「点４＋点１−（点２＋点３）」で求められる。
【０００９】
図１（ａ）や（ｂ）で求められる面積や画素値の合計の計算量は信号数Ｎに比例するため、カーネルサイズに依存しない。また、いったん面積や画素値の合計等を計算しておけば、決まった回数の加減算で積分計算と同等の結果を得ることができる。
【００１０】
ところで、図１で示した例は、信号であるｆ(x)や画素にフィルタ係数が１のフィルタによる畳み込み演算が施された例となっている。そこで、フィルタ係数を１以外の一般的な係数にも適用するための技術の１つとして三角関数を用いる計算方法が考案されている（例えば、非特許文献３、４参照。）。
【００１１】
非特許文献３には、フィルタ係数にＣｏｓ関数とＳｉｎ関数を適用するコサインインテグラルイメージ(Cosine Integral Image, 以下ＣＩＩと称す。)について開示されている。しかしながら、非特許文献３のＣＩＩは実部（Ｃｏｓ関数）と虚部（Ｓｉｎ関数）の２つの値を保持するためメモリ使用量が多い。また、非特許文献４には、ＣＩＩが必要とする大きなメモリ容量と頻繁なメモリアクセスを低減するため、画素値の合計を連続的に更新するＣＩＩについて開示されている。すなわち、画素値の合計が連続的に更新されるため画素値の合計を保持するためのメモリ容量を低減できる。
【発明の概要】
【発明が解決しようとする課題】
【００１２】
しかしながら、非特許文献４の計算方法ではメモリ使用量を抑制できるが、画素値の合計を連続的に合計する処理を含むため、非特許文献４の計算方法を適用できるのは全信号位置においてフィルタ処理を行う場合のみであるという問題がある。
【００１３】
本発明は、上記課題に鑑み、任意の信号位置に対し少ないメモリ使用量でフィルタ処理を施すことができる情報処理装置を提供することを目的とする。
【課題を解決するための手段】
【００１４】
本発明は、上記課題に鑑み、信号に線形フィルタを適用する情報処理装置であって、  初期値として与えられた１番目のベクトルを用いると共に、ｎ−１番目（ｎは２以上の整数）のベクトルと行列の積にｎ番目の信号値を加算してｎ番目のベクトルを求めるベクトル算出部と、前記ベクトル算出部が算出したｎ個のベクトルから取り出した２つのベクトルのうち、一方のベクトルと前記行列のべき乗との積から、他方のベクトルと前記行列のべき乗との積を減じた差分を前記信号に対する線形フィルタの出力値として出力する出力値算出部と、を有することを特徴とする。
【発明の効果】
【００１５】
任意の信号位置に対し少ないメモリ使用量でフィルタ処理を施すことができる情報処理装置を提供することができる。
【図面の簡単な説明】
【００１６】
【図１】積分計算の計算量を説明する図の一例である。
【図２】同伴行列べき乗型インテグラルイメージによる線形フィルタの計算方法を実行する情報処理装置のハードウェア構成図の一例である。
【図３】入力信号の一例及び標準のテンプレートの一例を示す図である。
【図４】情報処理装置の機能構成図の一例である。
【図５】情報処理装置がテンプレートマッチングを行う手順を示すフローチャート図の一例である。
【図６】二分法を模式的に説明する図の一例である。
【図７】ガウシアンフィルタの一例を示す図である。
【図８】情報処理装置が平滑化を行う手順を示すフローチャート図の一例である。
【図９】ガボールフィルタの一例を示す図である。
【図１０】情報処理装置がエッジ抽出を行う手順を示すフローチャート図の一例である。
【図１１】離散信号のリサンプリングの手順を模式的に示す図の一例である。
【図１２】情報処理装置がリサンプリングを行う手順を示すフローチャート図の一例である。
【発明を実施するための形態】
【００１７】
以下、本発明を実施するための形態について図面を参照しながら説明する。
【００１８】
＜インテグラルイメージの基本原理＞
本実施形態の線形フィルタについて説明するため、インテグラルイメージを用いたフィルタ処理の基本原理を説明する。背景技術で述べたようにインテグラルイメージを用いたフィルタ処理ではフィルタ係数に特徴があるが、以下ではこれらのインテグラルイメージを一般化したべき乗型インテグラルイメージが示される。後述するように、べき乗型インテグラルイメージには、フィルタ係数が実数のいわゆる通常のインテグラルイメージだけでなくＣＩＩ(Cosine Integral Image)などのインテグラルイメージが内包されている。そして、べき乗型インテグラルイメージから、本実施形態のインテグラルイメージ（後述する同伴行列べき乗型インテグラルイメージ）が新たに導出されることを説明する。なお、以下の説明で式の導出における動作の主体は出願人又は発明者である。
【００１９】
＜＜べき乗型インテグラルイメージ＞＞
信号数Ｎの信号系列ｆ(f0, f1, …, fN-1)に対するカーネルサイズＫの線形フィルタg(g0, g1, …, gK-1)の信号位置xにおける出力値(g*f)xを求める場合を例にして、べき乗型インテグラルイメージを導出する。なお、本実施形態の信号処理は代数学における環の定義を満たすことを前提とする。
【００２０】
【数２】
2016103089.tif 000004
式（１）において説明したように、右辺の畳み込み演算を式（２）のまま計算する際の計算量はＯ（Ｋ）である（Ｋ回の乗算とＫ−１回の加算）。
【００２１】
このフィルタ処理の高速化を図るため、Ｎ個の中間変数として(F0, F1, …, FN-1)を導入する（本実施形態では、この中間変数を「累積値」と称する。）。Fnは次の漸化式により算出される。
【００２２】
【数３】
2016103089.tif 000005
後述する情報処理装置は、この漸化式をn=0からN-1まで逐次計算することでｎ個のFnを求め、（F0, F1, …, FN-1)をメモリに記憶しておく。なお、ｎは２以上の整数であるとする。このように計算されたFnは、次のように変形することでfとαの積和で表現できることが分かる。
【００２３】
【数４】
2016103089.tif 000006
次に、Fn+KとαKFnを式（４）に代入して両者の差分を算出する。
【００２４】
【数５】
2016103089.tif 000007
さらに両辺にαmを乗ずると下式が得られる。
【００２５】
【数６】
2016103089.tif 000008
式（６）と式（２）を比較すると、式（２）の線形フィルタのフィルタ係数guを以下のように置き換えた場合に式（６）の左辺が線形フィルタの出力値(g*f)xと一致することが分かる。
【００２６】
【数７】
2016103089.tif 000009
したがって、後述する情報処理装置が式（６）の左辺に示されるαのべき乗及びFnを事前にテーブルで保持しておくことで、線形フィルタの出力値(g*f)xを得ることができる。式（６）の左辺は２回の乗算と１回の減算により求められるので、計算量のオーダーで言えば、Ｏ（１）である。すなわち、カーネルサイズに関わらず一定回数の計算量でフィルタ処理を実行できる。カーネルサイズＫに依存しないため、カーネルサイズが大きい線形フィルタでも固定時間で計算することができる。
【００２７】
なお、式（６）で必要な式（３）に基づくFnの計算の計算量（こちらも事前に算出される）は、１つのｎについて１回の乗算と１回の加算が必要なので、すべてのｎについてＯ（Ｎ）である。すなわち、式（６）により全信号位置で線形フィルタ処理を行う場合のフィルタ処理の計算量はカーネルサイズＫに依存せずＯ（Ｎ）である。畳み込み演算を行う場合のＯ（ＮＫ）及びＦＦＴを利用する場合のＯ（ＮlogＮ）よりも小さいことが分かる。
【００２８】
また、Fnはメモリに構築されたテーブル等に予め保持されることが可能なので、後述する情報処理装置などがいったんFnを計算した後は、任意のカーネルサイズＫの線形フィルタによる任意の信号位置における計算をＯ（１）の計算量で計算できる。
【００２９】
式（６）の左辺による線形フィルタの計算方法を「べき乗型インテグラルイメージ」と称することにする。
【００３０】
＜＜従来のインテグラルイメージの位置付け＞＞
従来のインテグラルイメージとは、背景技術などで述べた非特許文献１〜４のインテグラルイメージをいう。まず、一様係数の線形フィルタ（背景技術で述べたフィルタ係数が１の線形フィルタ）が用いられたインテグラルイメージに関し、べき乗型インテグラルイメージの式（６）においてα=1と置くと一様係数の線形フィルタに対する通常のインテグラルイメージが導かれる。
【００３１】
また、フィルタ係数を複素数に拡張しα=exp(-jω)＝cosω−jsinωと置けば、非特許文献３のＣＩＩと同等の処理が可能である（ｊは虚数単位）。
【００３２】
このように、べき乗型インテグラルイメージは従来のインテグラルイメージの一般化となっている。
【００３３】
＜＜同伴行列べき乗型インテグラルイメージ＞＞
これまではFnをスカラー値の信号系列としていたが、べき乗型インテグラルイメージはFnがＤ次元ベクトル（あるいは「環上の加群」の元）の場合、及び、αがＤ次元正方行列（正方行列とするのはべき乗するため）であっても成り立つ。
【００３４】
ここでαとして次の同伴行列Ａを用いる。
【００３５】
【数８】
2016103089.tif 000010
式（３）のFnに一次元の信号系列ｆを適用し、αに同伴行列Ａを適用すると式（９）のように書き直すことができる。式（９）の第一式は累積値Fnの初期値を与え、第二式はF0とｎ−１番目のFn-1からｎ番目のFnを求める漸化式である。すなわち、初期値として与えられた１番目のベクトルが与えられれば、ｎ−１番目（ｎは２以上の整数）のベクトルFn-1と同伴行列Ａの積に、ｎ番目の信号値fnを加算することでｎ番目のベクトルを求めることができる。
【００３６】
【数９】
2016103089.tif 000011
式（９）−１、式（９）−２、式（９）−３はそれぞれ説明のために示したF3までの計算結果である。式（９）の第二式（漸化式）と式（９）−１、式（９）−２、式（９）−３を見ると、Fnの第i要素はFn-1の第i-1要素と等しいことが分かる。すなわち、F１の第一要素はF2の第２要素に現れ、F2の第一要素はF3の第二要素に現れている。すなわちFn,i= Fn-1,i-1の関係がある。したがって、後述する情報処理装置はFnの第一要素のみをメモリに記録しておけば、すべてのnについてFnを再現することができる。
【００３７】
次に、同伴行列Ａを使って式（６）を書き直す。
【００３８】
【数１０】
2016103089.tif 000012
一方のベクトルFn+Kと行列のべき乗との積から、他方のベクトルFnと行列のべき乗との積を減じた差分により信号に対する線形フィルタの出力値が得られる。すなわち、式（１０）の左辺のベクトルの第一要素をメモリに記憶しておくことで、Ａのべき乗の１行１列の成分がフィルタ係数となる線形フィルタを実現できることが分かる。計算回数としては、２Ｄ回の乗算（第1項と第２項の行列計算における乗算）、２Ｄ−１回の加算又は減算（第１項の行列計算のＤ−１回の加算、第２項の行列計算のＤ−１回の加算、第１項と第２項の減算）であり、計算量はＯ（Ｄ）である。
【００３９】
なお、Ａのべき乗の１行１列の成分がフィルタ係数となる点について補足しておく。説明のためＡのべき乗を式（１０）−1のように置く（後述する式（１３）を簡易化したものである）。
【００４０】
【数１１】
2016103089.tif 000013
Ａ(m)ijは行列Ａをｍ乗した行列のｉ行ｊ列の成分を意味する。Ａが二次元ならば、式（１０）は式（１０）−２のように表せる。また、式（１０）−２の右辺は式（１０）−３のように展開できる。
【００４１】
【数１２】
2016103089.tif 000014
式（１０）−３に示すように、ベクトルの第一要素はＡのべき乗の行列の１行１列成分との畳み込みになっており、ベクトルの第二要素はＡのべき乗の行列の２行１列成分との畳み込みになっている。Ａのべき乗の行列の１行１列成分と２行１列成分は、後述する式（１３）で示すように位相しか違わない。また、位相は式（１０）のｍで調整できるため、ベクトルの第２要素は不要になる。したがって、Ａのべき乗の１行１列の成分がフィルタ係数となる。
【００４２】
続いて、Ｄ＝２の場合について具体的な例を用いて説明する。固有値がexp(-jω)の同伴行列Ａは式（１１）のように表現できる。
【００４３】
また、式（１１）の同伴行列Ａ及び前述の関係Fn,i= Fn-1,i-1により得られるFn-1を用いて式（９）の第二式を書き直すと式（９）の第二式は式（１１）−１のようになる。
【００４４】
【数１３】
2016103089.tif 000015
したがって、式（９）の更新式（第二式）からFnの第一要素は次の通りとなる。Fn,1はベクトルとして表された累積値Fnの第一要素である。
【００４５】
【数１４】
2016103089.tif 000016
次に線形フィルタの出力値の計算式を考える。Ａのべき乗は次の通りである。
【００４６】
【数１５】
2016103089.tif 000017
式（１３）を式（１０）に代入し両辺にsinωを乗じると式（１４）が得られる。式（１４）を導出する過程で、式の変形に前述の関係Fn,i= Fn-1,i-1を利用している（ベクトルの第一要素のみを取り出している）。
【００４７】
【数１６】
2016103089.tif 000018
式（１４）によれば、Fn+k,1、Fn+k-1,1、Fn,1、及び、Fn-1,1を式（１２）により求めれば、Ａのべき乗が式（１３）で表わされるフィルタ計算を実行できる。式（１４）を見ると明らかなように、右辺の信号値fxと三角関数との畳み込み演算が、左辺の累積値Fnに関する４回の乗算と３回の加算又は減算で実現できることが分かる。三角関数で計算される値は予めテーブルとして保持しておく。したがって、カーネルサイズＫに非依存で、式（１１）の同伴行列Ａによるフィルタ計算が可能であることを確認できる。
【００４８】
同伴行列Ａによって変わりうるが式（１４）又は式（１２）が本実施形態において線形フィルタの処理に用いられる式となる。式（１４）の左辺による線形フィルタの計算方法を「同伴行列べき乗型インテグラルイメージ」と称することにする。
【００４９】
＜＜ＣＩＩとの比較＞＞
式（１４）で示す線形フィルタの処理についてＣＩＩ等と比較する。上記のＣＩＩ又はべき乗型インテグラルイメージにおいてα=exp(-jω)とした場合のＣＩＩと同等のインテグラルイメージでも、畳み込み演算はカーネルサイズＫに非依存である。しかし、αとして複素数を用いた場合はFnも複素数であるので、実部と虚部の二つの値を保持しておく必要がある。
【００５０】
これに対して本実施形態のようにαに同伴行列を用いる場合、Fnの第一要素のみを保持しておけばよいので、ＣＩＩと比較してメモリ使用量を半分に低減することができる。
【００５１】
また、ＣＩＩを高速に計算する手法である非特許文献４はFnをメモリに保存する必要がない。しかし、非特許文献４では信号の合計を連続的に算出する処理を繰り返すため、非特許文献４の手法は全信号位置に対して線形フィルタを施す場合のみに効果的となる。すなわち、任意の信号位置に対し線形フィルタを適用することができない。
【００５２】
これに対して本実施形態の同伴行列べき乗型インテグラルイメージは、情報処理装置がフィルタを適用する位置に応じて累積値Fnを読み出すことで、任意の信号位置に対し線形フィルタを適用することができる。
【００５３】
＜＜複数周波数への対応＞＞
同伴行列べき乗型インテグラルイメージは、単一周波数の三角関数との畳み込み演算を高速化する。これに対し、任意係数の線形フィルタの処理を高速化するためには、予め線形フィルタの係数を異なる周波数に分解し、周波数ごとにインテグラルイメージを適用して出力値を得て、その出力値の加重和を取ることで目的の線形フィルタの出力値を得ることができる。周波数分解にはフーリエ変換、離散コサイン変換、離散サイン変換などを使うことができる、また、式（１４）においてmにより位相の調整ができる。
【００５４】
このように任意係数の線形フィルタの処理を高速化するためには、周波数ごとにインテグラルイメージを用意する必要があるため、本実施形態において高速化の効果が高いのは、線形フィルタが少数の周波数の三角関数の線形結合で近似できる場合である。後述するように、ガウシアンフィルタやガボールフィルタなど実用上重要な線形フィルタは少数の周波数の三角関数の線形結合で近似できる。また、周波数の数を制限しなければ任意の線形フィルタの係数を扱うことができる。
【００５５】
＜＜ 多次元への対応＞＞
これまでの説明では、一次元の信号系列f(f0, f1, …, fN-1)に対する線形フィルタについて説明したが、同伴行列べき乗型インテグラルイメージは静止画や動画などの多次元信号にも適用することができる。
【００５６】
一例としてi, jの二つのインデックスを持つ二次元信号{fij}を例にして説明する。まず、次の漸化式によりiに関する累積値Fijをすべてのj（０〜Ｊ−１）に対して算出する。
【００５７】
【数１７】
2016103089.tif 000019
次に、次の漸化式により累積値Fijを用いてjに関する累積値Gijをすべてのi（０〜Ｉ−１）に対して算出する。
【００５８】
【数１８】
2016103089.tif 000020
式（６）などこれまでの議論と同様に次式が成り立つ。
【００５９】
【数１９】
2016103089.tif 000021
式（１７）（１８）を利用すると、二次元信号{fij}とべき乗型係数との畳み込み演算が、次のようにGijの加減算で実行できることが分かる。
【００６０】
【数２０】
2016103089.tif 000022
なお、式（１９）の第二行から第三行への変形で、αとβの乗算の順序を入れ替えた。したがって、式（１９）が成り立つのはαとβについて乗法に対する可換性がある場合である。たとえば、αとβが実数体や複素数体の元であれば乗法の可換性があるので、本式を利用することができる。乗法の可換性がない場合は、各次元で線形フィルタを適用し、その線形フィルタの出力値に対して別の次元で線形フィルタを適用する、という逐次的なフィルタ処理を行うことで対応が可能である。
【００６１】
α、βが行列環の元である場合は一般には乗法の可換性がないので、逐次的なアプローチを取る必要がある。
【００６２】
式（１９）から明らかなように、α=β=1とおいた場合、式（１９）から一様係数の線形フィルタに対する通常のインテグラルイメージが導かれる。したがって、式（１９）はインテグラルイメージの一般化になっている。
【００６３】
また、式（１９）は二次元信号の例を説明しているが、一般に、Ｄ−１次元の累積値（二次元の例ではFij）の累積値（二次元の例ではGij）を取ることで、Ｄ次元での線形フィルタの計算が可能となる。したがって、多次元信号に対して高速かつ少ないメモリ使用量で線形フィルタの処理が可能となる。
【００６４】
以上説明したように、同伴行列べき乗型インテグラルイメージは、線形フィルタの処理を、漸化式による累積値の算出と累積値の差分演算で実現できるため、線形フィルタのカーネルサイズに依存しないで計算できる。また、同伴行列により累積値を求めるため累積値の第一成分のみをメモリに記憶しておけばよくメモリ使用量を低減することができる。
【実施例１】
【００６５】
以下の実施例では、上述した同伴行列べき乗型インテグラルイメージの適用例を説明する。本実施例では、一次元の離散信号の中から標準のテンプレートと類似した形状の信号列を検知するテンプレートマッチングへの適用例を説明する。
【００６６】
図２は、同伴行列べき乗型インテグラルイメージによる線形フィルタの計算方法を実行する情報処理装置１００のハードウェア構成図の一例を示す。情報処理装置１００は、汎用ＰＣ（Personal Computer）の他、プロジェクタなどの投影装置、電子黒板、スマートフォン、タブレット型端末、汎用ＰＣ（Personal Computer）、携帯電話、ウェアラブルＰＣ、ゲーム機器、カーナビゲーション端末、又は、ＰＤＡ（Personal Digital Assistant）等が挙げられる。
【００６７】
情報処理装置１００は、バス３１０に接続されたＣＰＵ３０１、ＲＯＭ３０２、ＲＡＭ３０３、ＨＤＤ３０５、ディスプレイ３０８、ネットワークＩ／Ｆ３０９、キーボード３１１、マウス３１２、メディアドライブ３０７、及び、光学ドライブ３１４を有する。ＣＰＵ３０１は、ＨＤ３０４に記憶されているプログラム３２０を実行して、情報処理装置１００の全体の動作を制御する。ＲＯＭ３０２はＩＰＬ等のＣＰＵ３０１の駆動に用いられるプログラムを記憶している。ＲＡＭ３０３はＣＰＵ３０１のワークエリアとして使用される。ＨＤ３０４は不揮発性メモリを搭載した記憶装置であり、同伴行列べき乗型インテグラルイメージによる線形フィルタの計算を行うプログラム３２０を記憶している。
【００６８】
ＨＤＤ(Hard Disk Drive)３０５はＣＰＵ３０１の制御にしたがってＨＤ３０４に対する各種データの読み出し又は書き込みを制御する。ディスプレイ３０８はカーソル、メニュー、ウィンドウ、文字、又は画像などの各種情報を表示する。ネットワークＩ／Ｆ３０９はＬＡＮやインターネットなどのネットワークとのインタフェースであり、他の機器との通信を行う。キーボード３１１及びマウス３１２は入出力装置であり、キーボード３１１は文字、数値、各種指示などの入力のための複数のキーを備えこれらからの入力を受け付ける。マウス３１２はマウスポインターの移動及び各種指示の選択や実行、処理対象の選択などを受け付ける。
【００６９】
メディアドライブ３０７はフラッシュメモリ等のメディア３０６に対するデータの読み出し又は書き込み（記憶）を制御する。光学ドライブ３１４は着脱可能な記録媒体の一例としてのＣＤ−ＲＯＭ(Compact Disc Read Only Memory)３１３等に対する各種データの読み出し又は書き込みを制御する。
【００７０】
なお、上記プログラム３２０は、インストール可能な形式又は実行可能な形式のファイルで、メディア３０６やＣＤ−ＲＯＭ３１３等のコンピュータで読み取り可能な記録媒体に記録して流通させるようにしてもよい。あるいは、プログラム３２０は、任意のサーバ型の情報処理装置からダウンロードされる形態で配布されてもよい。
【００７１】
＜信号波形の具体的な形状＞
図３を用いて本実施例で扱う信号波形等について説明する。図３（ａ）は入力信号の一例を、図３（ｂ）は標準のテンプレートの一例を示す。本実施例の情報処理装置１００は入力信号が標準のテンプレートと最も近い形状となる位置を探索する。
【００７２】
図３（ｂ）に示すように、情報処理装置１００は標準のテンプレートを異なるいくつかの角周波数で分解しておく。あるいは、既知の異なるいくつかの角周波数の波形を合成して標準のテンプレートを作成してもよい。図３（ｂ）では標準のテンプレートがsinω1xとsinω2xとの線形結合で表されている。
【００７３】
情報処理装置１００は、単一の周波数ごとに線形フィルタを入力信号に適用し、その出力値を加重加算することでテンプレートマッチングを実現する。加重加算にあたって、図２の例ではω1、ω2の正弦関数の重みをw(ω1)、w(ω1)として、各重みと単一周波数の線形結合で標準のテンプレートを表している。
【００７４】
＜情報処理装置の機能について＞
図４は、本実施例における情報処理装置１００の機能構成図の一例である。情報処理装置１００は信号入力部１１、累積値算出部１２、フィルタ位置決定部１３、累積差分算出部１４、加重加算部１５、終了判定部１６、データ出力部１７、複数周波数分解部１８、及び、記憶部１９を有している。
【００７５】
信号入力部１１はＲＡＭ３０３上に展開されたプログラム３２０に従ったＣＰＵ３０１からの命令によって動作することで実現される機能又は機能される手段である。この信号はＨＤ３０４、メディア３０６又はＣＤ−ＲＯＭ３１３に記憶されていてもよいし、ネットワークＩ／Ｆ３０９によりネットワークを介して取得されてもよい。信号入力部１１はテンプレートマッチングの対象となる信号の入力を受け付ける。
【００７６】
複数周波数分解部１８は、ＲＡＭ３０３上に展開されたプログラム３２０に従ったＣＰＵ３０１からの命令によって動作することで実現される機能又は機能される手段であり、線形フィルタの係数を複数の周波数に分解する。複数周波数分解部１８はフィルタ係数取得部の一例である。
【００７７】
累積値算出部１２は、ＲＡＭ３０３上に展開されたプログラム３２０に従ったＣＰＵ３０１からの命令によって動作することで実現される機能又は機能される手段であり、式（１２）により累積値Fnを算出して記憶部１９に記憶させる。累積値算出部１２はベクトル算出部の一例である。
【００７８】
フィルタ位置決定部１３は、ＲＡＭ３０３上に展開されたプログラム３２０に従ったＣＰＵ３０１からの命令によって動作することで実現される機能又は機能される手段であり、フィルタ処理を施す入力信号の位置n（フィルタの適用位置）を決定する。フィルタ位置決定部１３は適用位置決定部の一例である。
【００７９】
累積差分算出部１４はＲＡＭ３０３上に展開されたプログラム３２０に従ったＣＰＵ３０１からの命令によって動作することで実現される機能又は機能される手段であり、記憶部１９から累積値Fnを読み出して単一周波数のフィルタ出力値hnを算出する。累積差分算出部１４は出力値算出部の一例である。
【００８０】
加重加算部１５はフィルタ係数が複数の周波数の線形結合で表されている場合、累積差分算出部１４が算出した各周波数のフィルタ出力値hnを加重加算する。
【００８１】
終了判定部１６はＲＡＭ３０３上に展開されたプログラム３２０に従ったＣＰＵ３０１からの命令によって動作することで実現される機能又は機能される手段であり、フィルタによる計算（テンプレートマッチング）を終了するか否かを判定する。
【００８２】
データ出力部１７はＲＡＭ３０３上に展開されたプログラム３２０に従ったＣＰＵ３０１からの命令によって動作することで実現される機能又は機能される手段であり、計算結果をディスプレイ３０８等に出力したりＨＤ３０４に記憶したりする。データ出力部１７は適用位置出力部の一例である。なお、信号出力部１７１については実施例２、４にて説明し、エッジ強度出力部１７２については実施例３で説明する。
【００８３】
記憶部１９は、ＲＡＭ３０３、ＨＤ３０４及びメディア３０６などの記憶手段に構築され、累積値、フィルタ処理が施される入力信号の位置及びフィルタ出力値が記憶される。
【００８４】
図５は、本実施例において情報処理装置１００がテンプレートマッチングを行う手順を示すフローチャート図の一例である。
【００８５】
S101：まず、信号入力部１１が入力信号を取得する。また、複数周波数分解部１８はテンプレートとなるフィルタの係数を複数の周波数（例えば２〜３個）に分解する。あるいはテンプレートを線形結合で近似するための複数の周波数が与えられていてもよい。
S102：累積値算出部１２は入力信号に対して式（１２）により、全信号位置(n=0, 1, …, N-1)における累積値Fnを計算して記憶部１９に保存する。この時、式（１２）の角周波数ωは標準のテンプレートから分解された複数の周波数に対応して複数個用意される。図３の例では角周波数ωはω1とω2である。ωの集合をΩと称することにする。累積値はそれぞれの角周波数ωごとに算出される。ωに対応した累積値をF(ω)nと表記する。
【００８６】
S103：次に、フィルタ位置決定部１３がフィルタ処理を行う入力信号の位置ｎを決定する。これにより、全信号位置で式（１４）を算出する必要がなくなり計算量を低減できる。本実施例では、二分法でフィルタ処理を行う位置nを決定する。
【００８７】
図６は二分法を模式的に説明する図の一例である。ｈ(x)は信号の位置xにおけるフィルタの出力値を示す。図６（ａ）に示すように、フィルタ位置決定部１３は最初に一次元の信号系列の全体の両端x1とx2で式（１４）を用いてフィルタ出力値を求める。フィルタ位置決定部１３は両端のフィルタ出力値を比較する。次のサイクルでは（ステップＳ１０６でＮｏと判定される場合）、フィルタ出力値が大きい側の端点x2と（図６（ａ）では端点x2の方が大きいとした）、探索区間（x1とx2）の中点x3との区間が新たなフィルタ処理を行う位置nとなる。図６（ｂ）では中点x3のフィルタ出力値の方が端点x2のフィルタ出力値よりも大きいと判断されるものとする。このため、図６（ｃ）に示すように、次のサイクルでは中点x3のフィルタ出力値と、中点x3と端点x2との中点x4のフィルタ出力値とが比較される。このように、フィルタの適用位置が連続的でなくてもフィルタ処理が可能である。
【００８８】
フィルタ位置決定部１３はステップＳ１０６の判断で繰り返されるサイクルごとに、フィルタ処理を行う位置nを決定する。
【００８９】
S104：次に、累積差分算出部１４は、フィルタ位置決定部１３が決定した入力信号の位置ｎにおける累積値を記憶部１９から読み出し、式（１４）を用いて角周波数ω∈Ωに対するフィルタ出力値h(ω)nを算出する。
【００９０】
S105：次に、加重加算部１５は重みw(ω)の下で算出されたh(ω)nを加重加算し、ステップＳ１０２で求めた入力信号の位置nにおけるフィルタ出力値hnを算出する。
【００９１】
【数２１】
2016103089.tif 000023
S106：終了判定部１６はテンプレートマッチングの終了を判定する。終了基準としては、例えば、二分法の探索区間が信号の１サンプリング区間となり、これ以上探索区間を分割できなくなった場合に終了とする。あるいは、フィルタ出力値が特定閾値以上となった場合に終了としてもよい。終了と判定すればステップＳ１０７へ処理が進み、そうでなければ処理がステップＳ１０３に戻る。
【００９２】
S107：データ出力部１７は最終的に得られた、入力信号の位置を出力して処理を終了する。すなわち、二分法により最後にフィルタ出力値が大きい方であると判断された入力信号の位置が、線形フィルタのフィルタ出力値が最大となる位置なので、最もテンプレートと一致する入力信号の位置である。
【００９３】
以上の処理によって、入力信号中からテンプレートに類似する形状の部分の位置を高速に求めることができる。ステップＳ１０２の累積値の算出では計算量は信号数Ｎに比例するためＯ（Ｎ）で累積値を算出できる。ステップＳ１０４〜Ｓ１０５で累積値に対して決まった回数の計算（４回の乗算、３回の加算又は減算）で線形フィルタのフィルタ出力値hnが算出される。したがって、カーネルサイズ（テンプレートのサイズ）に非依存でテンプレートマッチングの処理が可能である。
【００９４】
本実施例の同伴行列べき乗型インテグラルイメージは信号位置に対して一つの累積値をメモリに保存すればよいので、２つの累積値（実部と虚部）の保存が必要なＣＩＩよりもメモリ使用量が少ない。また、非特許文献４の手法ではメモリ使用量は低減できるが、本実施例のステップＳ１０３に示したように任意のフィルタ位置を決定してフィルタ処理を行うことができない。
【００９５】
なお、本実施例の適用例には種々のバリエーションが含まれ得る。例えば、標準のテンプレートの周波数分解にはフーリエ変換、離散コサイン変換、離散サイン変換などを使うことができる。また、ＡＩＣ（赤池情報量基準）、ＴＩＣ（竹内情報量基準）、ＢＩＣ（ベイズ情報量基準）などの情報量基準を用いて、モデル化する周波数の数を最適化することもできる。
【００９６】
また、ステップＳ１０３では二分法によりフィルタ位置を決定したが、入力信号を数個に区切ってフィルタ出力値を求めてから、フィルタ出力値が大きい位置の近傍を細かく探索してもよい。また、ランダムサンプリングを用いてもよいし、模擬焼きなまし法（Simulated Annealing）を用いてもよい。この場合、ステップＳ１０６の終了判定は新しく探索する位置がなくなった時に終了とすればよい。
【実施例２】
【００９７】
本実施例ではガウシアンフィルタによる信号平滑化を同伴行列べき乗型インテグラルイメージにより実行する情報処理装置１００について説明する。本実施例では一次元の離散信号に対してガウシアンフィルタを適用し平滑化を行う。
【００９８】
ガウシアンフィルタのフィルタ係数guは、パラメータσを用いて次のように表される。
【００９９】
【数２２】
2016103089.tif 000024
図７はガウシアンフィルタの一例を示す図である。ガウシアンフィルタは低周波成分を多く含むため、信号の高周波成分を低減して平滑化することができる。平滑化効果を高めるためにはガウシアンフィルタのパラメータのσを大きくすることが好ましいが、これはカーネルサイズが大きくなることを意味する。従来の畳み込み演算でフィルタ処理を行うと処理時間が増大する。これに対し、本実施例ではカーネルサイズに非依存の計算量でガウシアンフィルタにより信号の平滑化を行うことができる。
【０１００】
ところで、ガウシアンフィルタのフィルタ係数guはＤＣＴを用いて、少数(2~3個)の異なる周波数のcosの線形結合で高精度に近似できることが知られている（非特許文献３参照）。本実施例においても少数の角周波数の集合Ωに対して、式（２２）のように余弦関数の線形結合でガウシアンフィルタのフィルタ係数guを近似する。
【０１０１】
【数２３】
2016103089.tif 000025
なお、式（１４）では正弦関数での畳み込み演算で定式化を行っているが、sin(x+π/2) = cos(x)の関係があるので、変数ｍを調整することで余弦関数の畳み込み演算が可能である。
【０１０２】
また、本実施例における情報処理装置１００の機能ブロック図は実施例１の図４と同じものを使用できる。しかしながら、フィルタ位置決定部１３は二分法ではなく記憶部１９に記憶されたフィルタ位置テーブルからフィルタ処理を行う入力信号の位置nを読み出す。フィルタ位置テーブルには予めフィルタ計算を行いたい位置nが登録されている。フィルタ位置テーブルは位置登録情報の一例である。
【０１０３】
図８は、本実施例において情報処理装置１００が平滑化を行う手順を示すフローチャート図の一例である。
S201：信号入力部１１が入力信号を取得する。また、複数周波数分解部１８はガウシアンフィルタの係数を複数の周波数（例えば２〜３個）に分解する。あるいはガウシアンフィルタの係数を線形結合で近似するための複数の周波数が与えられていてもよい。
S202：累積値算出部１２は全信号位置(n=0, 1, …, N-1)における角周波数ω∈Ωに対し式（１２）により累積値F(ω)nを算出し、記憶部１９に保存する。
【０１０４】
S203：フィルタ位置決定部１３はフィルタ処理を行う入力信号の位置nを決定する。フィルタ位置はフィルタ位置テーブルから読み出すことで決定する。
【０１０５】
S204：累積差分算出部１４が記憶部１９から累積値F(ω)nを読出し、式（１４）に基づき角周波数ω∈Ωに対するフィルタ出力値h(ω)nを算出する。
【０１０６】
S205：加重加算部１５は、算出されたh(ω)nを重みw(ω)の下で加重加算し、入力信号の位置nにおけるフィルタ出力値hnを算出する。
【０１０７】
S206：終了判定部１６は終了判定を行う。フィルタ位置テーブルに記録された全ての入力信号の位置でフィルタ処理がされていれば処理はステップＳ２０７に進む。ステップＳ２０６の判定がＮｏの場合、処理はステップＳ２０３に戻る。
【０１０８】
S207：信号出力部１７１は得られたフィルタ出力値hnを平滑化された信号として出力する。当然ながらフィルタ出力値hnは元の画像の画素の順番に並べられる。
【０１０９】
以上の処理によって、カーネルサイズに非依存の計算量でガウシアンフィルタによる入力信号の平滑化を行うことができる。ステップＳ２０２の計算量はＯ（Ｎ）であり、Ｓ２０４、Ｓ２０５では累積値に対して決まった計算量で線形フィルタのフィルタ出力値が算出されるため、計算量はカーネルサイズに非依存である。実施例１と同様にＣＩＩ及び非特許文献４に対する優位性も得られる。
【０１１０】
なお、本実施例のＳ２０３では予めフィルタ位置テーブルにフィルタ位置を記録していたが、平滑化を行いたい箇所をユーザがユーザインターフェースを通じてインタラクティブに指定してもよい。
【実施例３】
【０１１１】
本実施例ではガボールフィルタによるエッジ抽出を同伴行列べき乗型インテグラルイメージにより実行する情報処理装置１００について説明する。本実施例では一次元の離散信号に対してガボールフィルタを適用しエッジ抽出を行う。
【０１１２】
ガボールフィルタのフィルタ係数guは、パラメータσ, θ、及び
【０１１３】
【数２３Ｂ】
2016103089.tif 000026
を用いて次のように表される。
【０１１４】
【数２４】
2016103089.tif 000027
図９は、ガボールフィルタの一例を示す図である。形状から見て分かる通り、ガボールフィルタは信号値が大きく変化する位置でフィルタ出力値の絶対値が大きくなるので、エッジ抽出に利用することができる。
【０１１５】
次に、ガボールフィルタを同伴行列べき乗型インテグラルイメージに適用するために式（２３）を変形する。式（２３）は三角関数の加法定理を用いて次のように書き直すことができる。
【０１１６】
【数２５】
2016103089.tif 000028
実施例２のガウシアンフィルタと同様に余弦関数の線形結合でガボールフィルタのカーネルを近似する。
【０１１７】
【数２６】
2016103089.tif 000029
これにより式（２６）が得られる。
【０１１８】
【数２７】
2016103089.tif 000030
そして、Ω´＝｛ω±θ｜ω∈Ω｝とし以下のように置く
【０１１９】
【数２７Ｂ】
2016103089.tif 000031
ことで式（２６）は以下のように変形できる。
【０１２０】
【数２８】
2016103089.tif 000032
このようにガボールフィルタも少数の周波数の三角関数の線形結合として近似することができる。周波数にはωの他に
【０１２１】
【数２８Ｂ】
2016103089.tif 000033
が含まれるので、ガウシアンフィルタの二倍（
【０１２２】
【数２８Ｃ】
2016103089.tif 000034
の分だけ周波数が多い）の数の周波数により近似できる。
【０１２３】
ここで、ＦＦＴの計算量はＯ（ＮlogＮ）、同伴行列べき乗型インテグラルイメージの計算量は累積値Fnを算出しさえすれば、１つの周波数について決まった計算量である。よって、周波数の数をｓとすると下式が成り立つ場合にべき乗型インテグラルイメージを利用することに有効性がある。
ｓ＜Ｏ（ＮlogＮ）
なお、本実施例における情報処理装置１００の機能ブロック図は実施例１の図４と同じものを使用できる。フィルタ位置決定部１３がフィルタ処理を行う入力信号の位置を決定する方法は実施例２と同様である。
【０１２４】
図１０は、本実施例において情報処理装置１００がエッジ抽出を行う手順を示すフローチャート図の一例である。
S301：信号入力部１１が入力信号を取得する。また、複数周波数分解部１８はガボールフィルタの係数を複数の周波数（例えば４〜６個）に分解する。あるいはガボールフィルタの係数を線形結合で近似するための複数の周波数が与えられていてもよい。
S302：累積値算出部１２は全信号位置(n=0, 1, …, N-1)における角周波数ω∈Ω'に対し式（１２）により累積値F(ω)nを算出し、記憶部１９に保存する。
S303：フィルタ位置決定部１３がフィルタ処理を行う入力信号の位置を決定する。本実施例では、予めフィルタ計算を行いう入力信号の位置が記録されたフィルタ位置テーブルから読み出す。
S304：累積差分算出部１４は記憶部１９から累積値F(ω)nを読出し、式（１４）に基づき角周波数ω∈Ω'と
【０１２５】
【数２８Ｄ】
2016103089.tif 000035
に対するフィルタ出力値h(ω)nを算出する。
S305：加重加算部１５は算出されたh(ω)nを重みw(ω)の下で加重加算し、入力信号の位置ｎにおけるフィルタ出力値hnを算出する。
S306：終了判定部１６は終了判定を行う。フィルタ位置テーブルに記録された入力信号の全ての位置でフィルタ処理が行われていれば、処理はステップＳ３０７に進む、ステップＳ３０６の判定がＮｏであれば処理はステップＳ３０３に戻る。
S307：エッジ強度出力部１７２はステップＳ３０５で得られたフィルタ出力値hnを、フィルタ処理を行った入力信号の位置ｎに対応づけてエッジ強度信号として出力する。
【０１２６】
以上のように、本実施例の情報処理装置１００は、ステップＳ３０２においてＯ（Ｎ）で累積値を算出し、ステップＳ３０４、Ｓ３０５で累積値に対して決まった計算量の計算で線形フィルタのフィルタ出力値を算出できる。このため、カーネルサイズに非依存の計算量でガボールフィルタによる信号のエッジ抽出を行うこができる。実施例１、２と同様にＣＩＩ及び非特許文献４に対する優位性も得られる。
【０１２７】
なお、ステップＳ３０３ではフィルタ位置テーブルからフィルタ処理を行う入力信号の位置を記録しておくのでなく、エッジ抽出を行いたい箇所をユーザがユーザインターフェースを通じてインタラクティブに指定してもよい。
【実施例４】
【０１２８】
本実施例では同伴行列べき乗型インテグラルイメージを用いて１次元の離散信号から信号をリサンプリングする情報処理装置１００について説明する。入力信号のリサンプリングを動的に行う場合に有効である。
【０１２９】
図１１は、離散信号のリサンプリングの手順を模式的に示す図の一例である。図１１（ａ）は入力信号として入力された離散信号を示す。情報処理装置１００は、離散信号の間の値をいったんゼロ（0）で埋めて連続信号に変換する。図１１（ｂ）は離散信号がないｘの信号値がゼロにより埋められた離散信号を示す。
【０１３０】
次に、情報処理装置１００はsinc関数の線形フィルタ（以下、sincフィルタと呼ぶ）でエイリアス除去（アンチエイリアス処理）を行う。エイリアス除去はサンプリングされた位置で各値にsinc関数をそれぞれ適用して、適用結果を重ね合わせることで得られる。これにより、図１１（ｃ）に示すようになめらかに連続する入力信号が得られる。
【０１３１】
したがって、情報処理装置１００は、図１１（ｄ）に示すように信号の任意の位置で信号をリサンプリングすることができる。
【０１３２】
原信号がサンプリング周波数の半分以下の周波数成分のみを含む場合、この方法で完全に原信号を再現することができる（サンプリング定理）。
【０１３３】
実際のリサンプリングの計算では連続信号を再現する必要はないので、リサンプリングしたい位置で入力信号に対してsincフィルタを適用し、その出力値をリンサンプリング後の値とすれば良い。すなわち、一次元離散信号系(f0, f1, …, fN-1)から、位置x（実数）の信号値fxは次式で計算できる。
【０１３４】
【数２９】
2016103089.tif 000036
sinc(x)はsinc関数であり次のように定義される。
【０１３５】
【数３０】
2016103089.tif 000037
式（２８）を直接計算すると、Ｎ回の乗算とＮ−１回の加算が必要であり、計算量が大きい。そこで、同伴行列べき乗型インテグラルイメージを導入して式（２８）を計算することで高速化を行う。
【０１３６】
実施例１〜３と同様に、sinc関数を少数の周波数の三角関数の線形結合で近似する。
【０１３７】
【数３１】
2016103089.tif 000038
なお、本実施例における情報処理装置１００の機能ブロック図は実施例１の図４と同じものを使用できる。フィルタ位置決定部１３がフィルタ処理を行う入力信号の位置を決定する方法は実施例２と同様である。
【０１３８】
図１２は、本実施例において情報処理装置１００がリサンプリングを行う手順を示すフローチャート図の一例である。
S401：信号入力部１１が入力信号を取得する。また、複数周波数分解部１８はsinc関数を複数の周波数（例えば、２〜６個）の線形結合で近似する。あるいはsinc関数を線形結合で近似するための複数の周波数が与えられていてもよい。
S402：累積値算出部１２が全信号位置(n=0, 1, …, N-1)における角周波数ω∈Ωに対し式（１２）により累積値F(ω)nを算出し、記憶部１９に保存する。
S403：フィルタ位置決定部１３がフィルタ処理を行う入力信号の位置を決定する。フィルタ処理が行われる入力信号の位置はフィルタ位置テーブルから読み出す。
S404：累積差分算出部１４は記憶部１９から累積値F(ω)nを読み出し、式（１４）に基づき角周波数ω∈Ωに対するフィルタ出力値h(ω)nを算出する。
S405：加重加算部１５はステップＳ４０４で算出されたh(ω)nを重みw(ω)の下で加重加算し、位置nにおけるフィルタ出力値hnを算出する。
S406：終了判定部１６は終了判定を行う。フィルタ位置テーブルに記録された入力信号の全ての位置でフィルタ処理が行われていれば処理はステップＳ４０７に進む。ステップＳ４０６の判定がＮｏの場合、処理はステップＳ４０３に戻る。
S407：信号出力部１７１は得られたフィルタ出力値hnをリサンプリング信号として出力する。
【０１３９】
以上のように、本実施例の情報処理装置１００はステップＳ４０２のＯ（Ｎ）の計算量で累積値を算出し、ステップＳ４０４、Ｓ４０５で累積値に対して決まった計算量の計算で線形フィルタのフィルタ出力値を算出する。したがって、カーネルサイズに非依存の計算量でsincフィルタによる信号のリサンプリングを行うこができる。実施例１〜３と同様にＣＩＩ及び非特許文献４に対する優位性も得られる。
【０１４０】
なお、ステップＳ４０３ではフィルタ位置テーブルからフィルタ処理を行う入力信号の位置を記録しておくのでなく、エッジ抽出を行いたい箇所をユーザがユーザインターフェースを通じてインタラクティブに指定してもよい。
【０１４１】
また、情報処理装置１００は一度計算した累積値F(ω)nを使い続けることができるので、リサンプリングを繰り返し行う場合により高い計算効率を発揮する。たとえばディスプレイに拡大率を動的に変更しながら信号を表示する場合、動くＣＧ物体のテクスチャのレンダリングする場合などである。
【０１４２】
＜変形例＞
以上、本発明を実施するための最良の形態について実施例を用いて説明したが、本発明はこうした実施例に何等限定されるものではなく、本発明の要旨を逸脱しない範囲内において種々の変形及び置換を加えることができる。
【０１４３】
例えば、本実施形態ではべき乗型インテグラルイメージに同伴行列を適用したが、同伴行列と同等の性質を有する行列を適用することを制限するものではない。
【０１４４】
また、本実施形態では、同伴行列べき乗インテグラルイメージをテンプレートマッチング、ガウシアンフィルタ、ガボールフィルタ、及び、信号のリサンプリングに適用したが、同伴行列べき乗インテグラルイメージの適用範囲はこれらに限られない。
【符号の説明】
【０１４５】
１１    信号入力部
１２    累積値算出部
１３    フィルタ位置決定部
１４    累積差分算出部
１５    加重加算部
１６    終了判定部
１７    データ出力部
１８    複数周波数分解部
１９    記憶部
１００  情報処理装置
【先行技術文献】
【非特許文献】
【０１４６】
【非特許文献１】F.C. Crow, "Summed-area tables for texture mapping," in Proceedings of the 11th annual conference on Computer graphics and interactive techniques. ACM, 1984, pp. 207-212.
【非特許文献２】P. Viola and M. Jones, "Rapid object detection using a boosted cascade of simple features," in IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2001, vol. 1.
【非特許文献３】E. Elboher and M. Werman, "Cosine integral images for fast spatial and range filtering," in Image Processing  (ICIP), 2011 18th IEEE International Conference on. IEEE, pp. 89-92, 2011.
【非特許文献４】K. Sugimoto and S. Kamata, "Fast Image Filtering by DCT-Based Kernel Decomposition and Sequential Sum Update," in Image Processing  (ICIP), 2012 19th IEEE International Conference on. IEEE, pp.125-128, 2012.
【図１】
2016103089.tif 000040
【図２】
2016103089.tif 000041
【図３】
2016103089.tif 000042
【図４】
2016103089.tif 000043
【図５】
2016103089.tif 000044
【図６】
2016103089.tif 000045
【図７】
2016103089.tif 000046
【図８】
2016103089.tif 000047
【図９】
2016103089.tif 000048
【図１０】
2016103089.tif 000049
【図１１】
2016103089.tif 000050
【図１２】
2016103089.tif 000051