最近、画像のモーメントについて（やっと）理解したので備忘録として。

「そもそもなんで画像のモーメントから重心座標が求められるの？」というところが主です。

長文です。

1. 前置き

OpenCV-Python Tutorialsでは、画像のモーメントは領域の面積や重心の求め方の項で紹介されています。

領域(輪郭)の特徴 — OpenCV-Python Tutorials 1 documentation

重心は $C_{x} =\frac{M_{10}}{M_{00}}$ と $C_{y} =\frac{M_{01}}{M_{00}}$ の関係から求められます。

モーメントについての詳細はWikipedia（英語版）に、、、ということなので、Wikipediaを見にいってみる。。。

In image processing, computer vision and related fields, an image moment is a certain particular weighted average (moment) of the image pixels' intensities, or a function of such moments, usually chosen to have some attractive property or interpretation.

Image moments are useful to describe objects after segmentation. Simple properties of the image which are found via image moments include area (or total intensity), its centroid, and information about its orientation.

画像処理、コンピュータビジョンおよび関連分野では、画像モーメントは、画像ピクセルの強度の特定の加重平均（モーメント）、または通常はある魅力的な特性または解釈を有するように選択されたそのようなモーメントの関数である。

画像モーメントは、分割後のオブジェクトを記述するのに便利です。画像モーメントによって見出される画像の単純な特性には、面積（または全強度）、重心、およびその向きに関する情報が含まれる。

なんのこっちゃΣ(￣ロ￣lll)

下の方の項目を見ても、積分記号がいっぱい出てきて、ますますわからない、、、orz

という経緯があり、原理については謎のまま天下り的に使っていました。

2. 二値画像で考えてみる

単純化のために二値画像の場合に絞って考えることにします。

画像のモーメント $M$ と添字( $n,m$ など)をWikipediaのような積分記号ではなく、総和（シグマ）で表すと次のようになります。

$M_{nm} = \sum^{画像全体} (x座標^{n} \times y座標^{m} \times 画素値)$

シグマの添字の付け方がイレギュラーな感じですが、画像中の全ピクセルについての総和と考えてください。

具体的な作業としては、次のような二値画像について

f:id:plant-raspberrypi3:20181113134231p:plain

1ピクセルずつ（x座標、y座標、画素値）を調べます。

f:id:plant-raspberrypi3:20181113134310g:plain

次に、全てのピクセルについて

x座標のn乗
y座標のm乗

を求めます。

さらに、上記の2つと画素値（0または1）を掛け合わせます。
例えば、(15,17,1)のピクセルについて、n=2、m=0の場合は

$15^{2} \times 17^{0} \times 1 = 225 \times 1 \times 1 = 225$

となり、このピクセルでの計算結果は225となります。

全てのピクセルでの計算結果を足し合わせたものが画像のモーメント $M_{nm}$ です。

3. さらに単純化してみる

二値画像の場合、画素値が0のピクセルの値は

$x座標^{n} \times y座標^{m} \times 0 = 0$

となって必ず0となります。一方、画素値が1のピクセルは

$x座標^{n} \times y座標^{m} \times 1 = x座標^{n} \times y座標^{m}$

となります。「画素値1の領域」を $領域A$ と名付けることにします。

f:id:plant-raspberrypi3:20181113154155p:plain

画素値が0のピクセルは計算結果も0となって消え去るので、最終的に

$M_{nm} = \sum^{領域A} (x座標^{n} \times y座標^{m})$

が画像のモーメントということになります。添字の意味は、領域A内の全ピクセルについての総和と考えてください。

ちょっとだけ単純になりました(^ ^)ノ

4. 「なぜ、 $M_{00}$ が領域の面積になるのか」の答え

$n = 0$ 、 $m = 0$ のときを考えます。その場合、各ピクセルの計算結果は

$x座標^{0} \times y座標^{0} = 1$

なので、画像のモーメントは

$M_{00} = \sum^{領域A} (1)$

となります。すなわち、 $M_{00}$ は領域A内のピクセルの総数であることがわかります。

領域A内のピクセル数は、二値画像の世界では領域Aの面積と同値なので

$M_{00} =領域Aの面積$

となることがわかりました。

5. 「なぜ、 $(\frac{M_{10}}{M_{00}}, \frac{M_{01}}{M_{00}})$ が重心の座標になるのか」の答え

いよいよ、本丸の重心座標について。先ほどの領域Aの重心座標を仮に $(C_{x}, C_{y})$ とおきます。

f:id:plant-raspberrypi3:20181113144826p:plain

ここで領域の重心とは、重心を境に上側と下側、右側と左側のバランスがとれている点です。

f:id:plant-raspberrypi3:20181113144702p:plain

つまり、重心を原点ととると、全てのx座標の合計、全てのy座標の合計のそれぞれが0となります。

具体的には、重心を原点とした時の各ピクセルの座標を(X, Y)とすると

$\sum^{領域A} X = 0\$

$\sum^{領域A} Y = 0$

と書けます。次に、Xの方だけ考えます。

f:id:plant-raspberrypi3:20181113150323p:plain

Xは元座標を基準とした $x$ と重心座標 $C_{x}$ で以下のように表せます。

$X = x - C_{x}$

この式を、前述の $\sum^{領域A} X = 0$ に代入して、

$\sum^{領域A} (x - C_{x}) = \sum^{領域A} x - C_{x} \sum^{領域A} (1) = 0$

ゆえに

$C_{x} = \frac{\sum^{領域A} x}{\sum^{領域A} (1)}$

となります。

ここで、

$\sum^{領域A} x = \sum^{領域A} (x^{1} \times y^{0}) = M_{10}$

また、

$\sum^{領域A} (1) = \sum^{領域A} (x^{0} \times y^{0}) = M_{00}$

なので、最終的に

$C_{x} = \frac{M_{10}}{M_{00}}$

となりました。

Yの方も同じ手順で求められるので、OpenCV-Python Tutorialsで説明されていた

重心は $C_{x} =\frac{M_{10}}{M_{00}}$ と $C_{y} =\frac{M_{01}}{M_{00}}$ の関係から求められます。

が成り立つことがわかりました。

理解してみると意外と難しくなかった。。。

6. numpyとscikit-imageで試す

・numpyで計算してみる

手を動かすほうが実感できるかなー、ということで。

まず、二値の円を描きます。

from skimage import draw, measure
import numpy as np
import matplotlib.pyplot as plt

img = np.zeros((10, 10), dtype=np.uint8)
rr, cc = draw.circle(4, 4, 5)
img[rr+1, cc] = 1

plt.imshow(img)
plt.show()

f:id:plant-raspberrypi3:20181113155523j:plain

黄色が1のピクセル、濃紺が0のピクセルを表しています。

次に、この円の画像のモーメントを求めてみます。

計算を簡単にするために、各座標のx、yの行列を作成しておきます。

y = np.zeros((10, 10), dtype=np.uint8)
y += np.arange(0,10,1, dtype=np.uint8)
x = y.T
print("x =")
print(x)
print("y =")
print(y)

#出力
x =
[[0 0 0 0 0 0 0 0 0 0]
 [1 1 1 1 1 1 1 1 1 1]
 [2 2 2 2 2 2 2 2 2 2]
 [3 3 3 3 3 3 3 3 3 3]
 [4 4 4 4 4 4 4 4 4 4]
 [5 5 5 5 5 5 5 5 5 5]
 [6 6 6 6 6 6 6 6 6 6]
 [7 7 7 7 7 7 7 7 7 7]
 [8 8 8 8 8 8 8 8 8 8]
 [9 9 9 9 9 9 9 9 9 9]]
y =
[[0 1 2 3 4 5 6 7 8 9]
 [0 1 2 3 4 5 6 7 8 9]
 [0 1 2 3 4 5 6 7 8 9]
 [0 1 2 3 4 5 6 7 8 9]
 [0 1 2 3 4 5 6 7 8 9]
 [0 1 2 3 4 5 6 7 8 9]
 [0 1 2 3 4 5 6 7 8 9]
 [0 1 2 3 4 5 6 7 8 9]
 [0 1 2 3 4 5 6 7 8 9]
 [0 1 2 3 4 5 6 7 8 9]]

それではいよいよ、 $M_{00}$ 、 $M_{10}$ 、 $M_{01}$ の計算。

for n,m in [[0,0],[1,0],[0,1]]:
    M_nm = np.sum(np.power(x,n) * np.power(y,m) * img)
    print(f"M{n}{m}: {M_nm}")

#出力
M00: 69
M10: 345
M01: 276

ここから、重心座標は

centroid = (345/69,276/69)
print(f"モーメントから求めた重心座標 : {centroid}")

#出力
モーメントから求めた重心座標 : (5.0, 4.0)

と計算できます。

重心を黒く塗りつぶして位置を確認してみます。

img[int(centroid[0]),int(centroid[1])] = 0

plt.imshow(img)
plt.show()

f:id:plant-raspberrypi3:20181113173155j:plain

大丈夫そうですね。

・scikit-imageでラクしてモーメントを求める

scikit-imageでは、画像のモーメントは
measureモジュールのmomentsという関数で求められます。

M = measure.moments(img)
print(M)

#出力
[[6.90000e+01 2.76000e+02 1.48000e+03 8.92800e+03]
 [3.45000e+02 1.38000e+03 7.40000e+03 4.46400e+04]
 [2.10100e+03 8.40400e+03 4.44400e+04 2.64352e+05]
 [1.42650e+04 5.70600e+04 2.96600e+05 1.73328e+06]]

画像のモーメントから面積と重心座標を計算してみます。

print(f"M00 (面積) : {M[0,0]}")
print(f"M10 : {M[1,0]}")
print(f"M01 : {M[0,1]}")

print(f"モーメントから求めた重心座標：{(M[1,0]/M[0,0],M[0,1]/M[0,0])}")

#出力
M00 (面積) : 69.0
M10 : 345.0
M01 : 276.0
モーメントから求めた重心座標：(5.0, 4.0)

numpyで求めた値と一致しています。よしよし。

・scikit-imageのさらに便利な関数たち

ちなみに、measureモジュールのregionprops関数を使うと
もっと簡単に面積・重心が求められます。

result_props = measure.regionprops(img)
print(f"面積: {result_props[0].area}")
print(f"重心座標: {result_props[0].centroid}")

#出力
面積: 69
重心座標: (5.0, 4.0)

measureモジュールにはmoments_centralという関数もあります。
これを使うと、重心を原点とした時の画像のモーメントが求められます。

M = measure.moments_central(img)

print(M[0,0])
print(M[1,0])
print(M[0,1])

#出力
69.0
0.0
0.0

これはラクチンですね(^ ^)ノ

plant-raspberrypi3のブログ

ラズベリーパイ3とPythonに挑戦して、植物を愛でたり画像を加工したりします。最近はscikit-imageの勉強してます。

画像のモーメントについての備忘録

目次

1. 前置き

2. 二値画像で考えてみる

3. さらに単純化してみる

4. 「なぜ、 $M_{00}$ が領域の面積になるのか」の答え

5. 「なぜ、 $(\frac{M_{10}}{M_{00}}, \frac{M_{01}}{M_{00}})$ が重心の座標になるのか」の答え

6. numpyとscikit-imageで試す

・numpyで計算してみる

・scikit-imageでラクしてモーメントを求める

・scikit-imageのさらに便利な関数たち

目次

1. 前置き

2. 二値画像で考えてみる

3. さらに単純化してみる

4. 「なぜ、が領域の面積になるのか」の答え

5. 「なぜ、が重心の座標になるのか」の答え

6. numpyとscikit-imageで試す

・numpyで計算してみる

・scikit-imageでラクしてモーメントを求める

・scikit-imageのさらに便利な関数たち

4. 「なぜ、 $M_{00}$ が領域の面積になるのか」の答え

5. 「なぜ、 $(\frac{M_{10}}{M_{00}}, \frac{M_{01}}{M_{00}})$ が重心の座標になるのか」の答え