Innovotionブログ - Kinectのための平滑化フィルタ(2) のためのノートです。

平滑化フィルタ その1

テストデータの準備

平滑化フィルタの性質を見るためにまず1次元のテストデータを用意し、後にそのテストデータが各フィルタによってどのように平滑化されるかを見ていきます。Kinectで取得した各関節データのx,y,z座標を想定していますが、当面はある程度の正確度と精度があることを仮定します。つまり異常な外れ値がある場合や値が取得できないときのことは考えません。

In [1]:
import numpy as np
In [2]:
import matplotlib.pyplot as plt
In [3]:
%matplotlib inline
plt.rcParams['figure.figsize'] = (12.0, 9.0)
In [4]:
import random

ここまではおまじないです。

In [5]:
t=np.arange(0.0,10.0,1.0/15)

x軸とする時間軸データです。0秒から10秒まで、毎秒15フレームだと思ってください。

In [6]:
tData=np.zeros(t.size,np.float32)
In [7]:
for i in range(t.size):
    if t[i]<2:
        tData[i]=0
    elif t[i]<3:
        tData[i]=t[i]-2
    elif t[i]<5:
        tData[i]=1
    elif t[i]<8:
        tData[i]=np.cos(0.5*np.pi*(t[i]-5));
    else:
        tData[i]=0
In [8]:
plt.plot(t,tData)
plt.ylim(-1.0,1.2)
Out[8]:
(-1.0, 1.2)

真の値としてはこのグラフのようなデータを使います。 静止→等速運動→静止→滑らかに非等速運動→静止 という動きを想定して定義しました。

In [9]:
oData=np.zeros(t.size,np.float32)
In [10]:
for i in range(oData.size):
    oData[i] = tData[i] + random.gauss(0,0.05)
In [11]:
plt.plot(t,tData)
plt.plot(t,oData,'.')
Out[11]:
[<matplotlib.lines.Line2D at 0x10ec55c10>]

誤差を含む観測データとして、真のデータに平均値0,分散0.05の正規分布を加えました。ばらつきが欲しいだけなのでここでは正規分布の定義などには立ち入りませんが、もっと分散の値を大きくすればばらつきも大きくできます。 さて、これでテストデータの準備ができました。この観測データ(緑の点)を各平滑化フィルタにかけ、その結果が真値データ(青い線)にどれくらい近くなっているかを見ていきます。

記号の準備

フィルタの説明のため、少しだけ記号を準備します。Microsoft社の"Skeletal Joint Smoothing White Paper"と同様です。以下、$i$フレーム目の観測データ、すなわちフィルターへの入力を$X_i$、フィルターの出力を$\widehat{X}_i$と書くことにします。

単純移動平均(Simple Moving Average)フィルタ

In [12]:
def SimpleAverage(data,N):
    filtered = np.zeros(data.size)
    for i in range(N):
        filtered[i]=data[i]
    for i in range(data.size):
        if i>N:
            sum=0
            for j in range(N):
                sum = sum + data[i-j]
            filtered[i] = sum/N
    return filtered

単純移動平均フィルタは非常に単純で、自然数のパラメータ$N$に対して$i$フレーム目の出力を$i-N+1$フレームから$i$フレームの$N$フレーム分の入力データを平均したものです。

$$\widehat{X}_i=\frac{1}{N}\sum _{k=0}^{N-1} X_{i-k} $$

In [13]:
pData1=SimpleAverage(oData,5)
plt.plot(t,tData,"b")
plt.plot(t,oData,'.g')
plt.plot(t,pData1,"r")
Out[13]:
[<matplotlib.lines.Line2D at 0x10f053f10>]

N=5としてフィルターにかけた結果が赤い線です。$0<t<2 $ や $3<t<5$ の静止状態では悪くありませんが、動いているときには青い線からはっきり遅れているのが読み取れると思います。Nの値をさらに大きくすればより滑らかになりますが、より遅延がひどくなりますので試してみてください。遅延が問題とならない場合や、プレイヤーが動かないときには使えそうなフィルターだということが言えそうです。

二重移動平均(Double Moving Average)フィルター

In [14]:
def DoubleMovingAverage(data,N):
    filtered1 = np.zeros(data.size)
    filtered2 = np.zeros(data.size)
    for i in range(N):
        filtered1[i]=data[i]
    for i in range(N):
        filtered2[i]=filtered1[i]
    for i in range(data.size):
        if i>N:
            sum=0
            for j in range(N):
                sum = sum + data[i-j]
            filtered1[i] = sum/N
    for i in range(data.size):
        if i>N:
            sum=0
            for j in range(N):
                sum = sum + filtered1[i-j]
            filtered2[i] = sum/N
    #return filtered1*2-filtered2
    return filtered2

二重移動平均フィルタは単純移動平均フィルタを2重にかけたものです。 $$Y_{i}=\frac{1}{N}\sum _{k=0}^{N-1} X_{i-k} $$ $$\widehat{X}_i = \frac{1}{N}\sum _{k=0}^{N-1} Y_{i-k}$$ 必ずしも$N$は同一に取る必要はありませんが、ここでは同じ$N$だけの平均をくり返しています。

In [15]:
pData2=DoubleMovingAverage(oData,4)
pData1=SimpleAverage(oData,4)
plt.plot(t,tData,'b')
plt.plot(t,oData,'.g')
plt.plot(t,pData2,'r')
plt.plot(t,pData1,'y')
Out[15]:
[<matplotlib.lines.Line2D at 0x10f080f90>]

赤い線が$N=4$として二重移動平均フィルタをかけたもの、黄色い線が同じく$N=4$として単純移動平均フィルタをかけたものです。赤い線の方が黄色い線より滑らかになっていますが、動いているときの遅延はさらに大きくなっています。ここで、線の間隔に注目すると、赤い線、黄色い線、青い線の間隔が概ね等しくなっています。実際、真の値が線型であるとき、すなわち一次関数で表されるときはこれは厳密な意味で成り立ちます。すなわち、真の値 $$T_i=ai+b$$ に各フィルタをかけたものをそれぞれ$SA_i, DA_i$とおくと、 $$\begin{eqnarray}SA_i&=&\frac{1}{N}\sum_{k=0}^{N-1} T_{i-k} = \frac{1}{N}\sum_{k=0}^{N-1} \left\{a(i-k)+b\right\}=ai+b-\frac{N-1}{2}a \\ DA_i&=&\frac{1}{N}\sum_{k=0}^{N-1} SA_{i-k} = \frac{1}{N}\sum_{k=0}^{N-1}\left\{ a(i-k)+b-\frac{N-1}{2}a \right\}=ai+b-\frac{N-1}{2}a-\frac{N-1}{2}a \end{eqnarray}$$ なので $$ T_i-SA_i=SA_i-DA_i$$ が成り立っています。そこで、真の値が線型であると予想されるとき、つまり動作が等速であると予想されるときは、二重移動平均フィルタと単純移動平均フィルタを組み合わせて遅延をキャンセルするという手があります。

In [16]:
pData3=2*SimpleAverage(oData,4)-DoubleMovingAverage(oData,4)
plt.plot(t,tData,'b')
plt.plot(t,oData,'.g')
plt.plot(t,pData3,'r')
Out[16]:
[<matplotlib.lines.Line2D at 0x10edc5dd0>]

滑らかさは単純移動平均フィルタと同程度ですが、等速運動時の遅延はきれいにキャンセルできているのが分かります。その一方で$t=3$など速度が変化した点の直後にこれまで見られなかったような真の値からの外れが見えます。等速運動時の遅延がないことは大きなメリットですが、速度変化には弱いフィルタであると言えます。これをマウスポインタの座標に直接使うとすれば、手を動かしているときは快適だが思ったところでぴたりと止まらない、という動きになるかもしれません。

指数平滑化(Exponential Smoothing)フィルタ

これはパラメータ$0\leq \alpha\leq 1$に対して次のように定義されます。 $$ \begin{eqnarray} \widehat{X}_i&=&\alpha X_i+(1-\alpha)\widehat{X}_{i-1},\\ \widehat{X}_0&=&\alpha X_{0}. \end{eqnarray}$$これは新たな入力値$X_{i}$と直前の出力値$\widehat{X}_{i-1}$とを$1-\alpha : \alpha$の割合で混ぜ合わせることで滑らかさを出していると考えられます。また、書き換えると $$\widehat{X}_i=\alpha \sum_{k=0}^{i} (1-\alpha)^{i-k}X_{k}$$ となり、大ざっぱに言えば指数的に重みをつけて平均を取っていると見ることもできるので指数平滑化フィルタという名前がついています。ただし、重みの総和が $$\alpha \sum_{k=0}^{i} (1-\alpha)^{i-k}=\alpha \frac{1-(1-\alpha)^{i+1}}{1-(1-\alpha)}=1-(1-\alpha)^{i+1}$$ と厳密に言えば1より小さくなっていますので、$i$が小さいときには注意が必要なケースもあるかもしれません。

In [17]:
def ExponentialSmoothing(data,alpha):
    filtered = np.zeros(data.size)
    filtered[0]=alpha*data[0]
    for i in range(data.size):
        if i>0:
            filtered[i]=alpha*data[i]+(1-alpha)*filtered[i-1]
    return filtered
In [18]:
pData4=ExponentialSmoothing(oData,0.9)
pData5=ExponentialSmoothing(oData,0.5)
pData6=ExponentialSmoothing(oData,0.1)
plt.plot(t,tData,'b')
plt.plot(t,oData,'.g')
plt.plot(t,pData4,'r')
plt.plot(t,pData5,'k')
plt.plot(t,pData6,'y')
Out[18]:
[<matplotlib.lines.Line2D at 0x10f0b1610>]

パラメータ$\alpha=0.9,0.5,0.1$についてそれぞれ赤、黒、黄の線でプロットしました。$\alpha$の値が小さくなるにつれて滑らかさが増す一方で遅延がひどくなっているのが分かります。滑らかさと遅延の大きさがきれいにトレードオフの関係になっていますが、目的に応じた微妙なパラメータの調整がやりやすいのが指数平滑化フィルタの利点です。

二重指数平滑化(Double Exponential Smoothing)フィルタ

二重指数平滑化フィルタは、パラメータ$0\leq \alpha ,\gamma\leq 1$に対して $$ \begin{eqnarray} \widehat{X}_i&=&\alpha X_{i}+(1-\alpha )\left( \widehat{X}_{i-1} +b_{i-1} \right) ,\\ b_i&=&\gamma \left( \widehat{X}_{i}-\widehat{X}_{i-1} \right) +(1-\gamma)b_{i-1},\\ \widehat{X}_{0}&=&X_{0},\\ b_0&=&X_{0}, \end{eqnarray}$$ で定義されます。ここで2つ目の漸化式に注目すると、$b_{i-1}$は$\left\{ \widehat{X}_{k}-\widehat{X}_{k-1}\right\}$というデータを指数平滑化したものになっており、$i-1$フレームにおける1フレームあたりの速度と見ることができます。したがって、その速度が続くとすると次フレームでの出力値は$\widehat{X}_{i-1}+b_{i-1}$になりそうです。その予測値が一つ目の漸化式の第二項に入っています。したがって、このフィルタは1つ前のフレームの速度から予測された値と新たな入力データを混合したものだと見ることができます。

In [19]:
def DoubleExponentialSmoothing(data,alpha,gamma):
    filtered=np.zeros(data.size)
    b=np.zeros(data.size)
    b[0]=data[1]-data[0]
    filtered[0]=data[0]
    for i in range(data.size):
        filtered[i]=alpha*data[i]+(1-alpha)*(filtered[i-1]+b[i-1])
        b[i]=gamma*(filtered[i]-filtered[i-1])+(1-gamma)*b[i-1]
    return filtered
In [20]:
pData5=ExponentialSmoothing(oData,0.4)
pData6=DoubleExponentialSmoothing(oData,0.4,0.5)
plt.plot(t,tData,'b')
plt.plot(t,oData,'.g')
plt.plot(t,pData6,'r')
plt.plot(t,pData5,'y')
Out[20]:
[<matplotlib.lines.Line2D at 0x10f0d9e50>]

パラメータ$\alpha=0.4$の指数平滑化フィルタの出力を黄色で、$\alpha=0.4,\gamma=0.5$の二重指数平滑化フィルタの出力を赤でプロットしています。滑らかさは変わらないようですが、等速運動している区間においては遅延がかなり解消されているのが分かります。