numpyで積分画像を作る - カオスノカケラ

numpyで積分画像を作る際は img.cumsum(axis=0).cumsum(axis=1) が良いよという話。

積分画像 (integral image) とは

数式的表記

画像 ${I(x,y)}$ があるとき、その積分画像 ${ii(x,y)}$ は次のように表される。

${ \displaystyle ii(x,y) = \sum_{w=1}^{x} \sum_{h=1}^{y} I(w,h) }$

となる。ここで画像 ${I(x,y)}$ は左上の画素を原点 ${(1,1)}$ として右方向に ${x}$ 軸を、下方向に ${y}$ 軸をとったものであることに注意。

numpy的表記

わかりやすいように画像を白黒画像(1チャンネル)として、 ${(H,W)-}$ 型ndarrayとする。すなわち img[0,0] は画像左上の画素、 img[H-1,W-1] は右下の画素、img[0,W-1] は右上の画素を表す。

積分画像は次の図のように定義され、積分画像における各画素の値は、元の画像で左上の画素からその画素までの区間の総和となる。つまり ii[y,x] == img[:y+1,:x+1].sum() の関係が成り立つ。

f:id:yuuhopro:20200616205555j:plain — 積分画像の説明図

応用

これがいったい何の役に立つかというと、主に画像中の矩形領域の画素の総和を計算するのに使える。 (矩形はここでは全ての辺が ${x}$ 軸または ${y}$ 軸に平行な長方形(正方形)のことをいう。 AABB: Axis Aligned Bounding Box )

具体的に考えてみる。

f:id:yuuhopro:20200616203207j:plain — 画像中の矩形

この画像中の矩形箇所(図中斜線領域)の画素値の総和を求める際に、通常であれば二重のfor文で全ての画素値を足し合わせることになる。

つまり

def area_sum(img,x1,x2,y1,y2):
    result = 0
    for h in range(y1,y2):
        for w in range(x1,x2):
            result += img[h,w]
    return result

または、numpyの機能を使うと

def area_sum(img,x1,x2,y1,y2):
    return img[y1:y2,x1:x2].sum()

となる。いずれにしてもアルゴリズム上では同じ計算量である。 ${x_2-x_1=N}$ 、 ${y_2-y_1=M}$ とすると計算量は ${O(NM)}$ である。

一方で積分画像を予め用意しておいた場合は、

def area_sum(ii, x1,x2,y1,y2):
    return ii[y2,x2] + ii[y1,x1] - ii[y2,x1] - ii[y1,x2]

という ${O(1)}$ のアルゴリズムで求まる。これは、図形的に表すと以下のようになる。

f:id:yuuhopro:20200616212435j:plain — 矩形領域の計算方法

計算で使用している全ての項が、左上から数えた画素値の和つまり積分画像に記憶されている値になっていることがミソである。これにより簡単に画像中の任意の矩形範囲の画素値の和を高速に求めることができるようになった。これは求めたい範囲が複数個ある場合に非常に有効になる。一箇所しか求めない場合は積分画像を作成する際の計算量でメリットが打ち消されてしまうが、調べたい矩形が何千箇所もの数になる場合は圧倒的な差である。

積分画像の作成

積分画像作成の愚直な実装では以下のようになる。

ii = np.empty((H,W))
for h in range(H):
    for w in range(W):
        ii[h,w] = img[:h+1,:w+1].sum()

当然このような素人のアルゴリズムは採用すべきではない。通常であれば累積和により求めることが妥当である。つまり以下のようにすると .sum() の部分の二乗オーダーの計算が省略でき、高速である。

ii = np.empty((H,W))
ii[0,0] = img[0,0]
for w in range(1,W):
    ii[0,w] = ii[0,w-1] + img[0,w]
for h in range(1,H):
    ii[h,0] = ii[h-1,0] + img[h,0]
    for w in range(1,W):
        ii[h,w] = ii[h,w-1] + img[h,w]

しかしながら、pythonを使う以上はfor文に依る低速化は無視できない。そこで numpy の cumsum() 関数により、各列での累積和を計算し、その値を用いて各行の累積和を計算する。

まず内部的な動作を示す。以下のように各列ごとの累積和を求め、その後に各行ごとの累積和を求めるようにすると自動的に積分画像ができる。

ii = np.empty((H,W))

# cumsum(axis=0)
for w in range(W):
    ii[0,w] = img[0,w]
    for h in range(1,H):
        ii[h,w] = ii[h-1,w] + img[h,w]

# cumsum(axis=1)
for h in range(H):
    ii[h,0] = ii[h,0] + img[h,0]
    for w in range(1,W):
        ii[h,w] = ii[h,w-1] + img[h,w]

これは単純に

ii = img.cumsum(axis=0).cumsum(axis=1)

によって計算可能であり、完全にfor文を排して高速に計算できる。

おまけ cv2.integral, cv2.integral2 について

画像処理、画像認識に使えるライブラリの一つである cv2(OpenCV) は積分画像計算用の関数を備えている。

ii = cv2.integral(img)

で積分画像が求まるが、numpyの cumsum() を用いたものと違い、最初に 0 が挿入される。

つまり [ 1,1,1,1,1 ] という配列に対してそれぞれの関数を適用すると、

cumsum()適用: 出力 [1,2,3,4,5]
integral()適用: 出力[0,1,2,3,4,5]

となる。したがって integral() では入力サイズと出力サイズは異なる。

また、 cv2.integral2() を利用すれば、

ii, ii2 = cv2.integral2(img)

のように戻り値は二つであり、それぞれは

ii : cv2.integral( img )
ii2 : cv2.integral( img**2 )

となっている。つまり通常の積分画像に加えて、同時に画素値を2乗した場合の積分画像も計算してくれる。

これは、おそらく矩形範囲の画素値の分散などの計算に便利である。