plant-raspberrypi3のブログ

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

M5Stick-CからMacにBluetoothで文字列を送信する

M5Stick-Cについて、一年ほど前に以下の記事を書いていました。

plant-raspberrypi3.hatenablog.com

今回、M5Stick-CのデータをMacに送信し、Macで（Pythonコードで）受け取る方法がわかったので、備忘録です。

Bluetooth経由でシリアル通信で送ります。

参考にさせていただいたのは以下のページです。

qiita.com

www.sejuku.net

1. M5Stick-CからBluetoothで文字列を送信

Arduino IDEでM5Stick-Cに以下のコードを書き込みます。

#include "BluetoothSerial.h"
#include <M5StickC.h>

BluetoothSerial bts;

void setup() {
  M5.begin();
  M5.Lcd.setRotation(3);
  M5.Lcd.fillScreen(BLACK);
  M5.Lcd.setCursor(0, 0, 2);
  M5.Lcd.println("Bluetooth TEST");
  bts.begin("M5StickC");
}

void loop() {
  M5.Lcd.setCursor(0, 20, 2);
  M5.Lcd.printf("Start Bluetooth");
  bts.println("Success");
  delay(1000);
}

bts.beginの引数はMac側でBluetoothを認識するときに表示される名前です。

bts.printlnの引数がMacに送信したい文字列です。

2. Mac側で文字列を受信（PythonのPySerial使用）

pyserialをインストール。

$ pip install pyserial

Macのシステム環境設定のBluetoothで「M5StickC」に接続する。

ターミナルで下記のコードを打ち込み、該当のBluetooth接続デバイスのパスを取得。

ls -l /dev/tty.*

今回は/dev/tty.M5StickC-ESP32SPPでした。

Pythonで以下のコードを動かすと、文字列が取得できるはずです。

import serial

m5data = serial.Serial('/dev/tty.M5StickC-ESP32SPP',9600, timeout=3)
line = m5data.readline()
print(line)
m5data.close()

b'Success\r\n'

これと同じ手順で、M5Stick-Cでセンサから取得した数値データなども送信できるはずです。

M5StickC ENV HatをArduino IDEで使ってみた

Raspberry Piと連携させて使いたい！ということでM5StickCを買ってみました。

www.switch-science.com

温度と湿度をモニターしてみよう、ということで、ENV Hatも。

www.switch-science.com

主に、Macでの設定方法の備忘録です。

この辺りを参考に。

make-muda.net

pages.switch-science.com

raspberrypi.mongonta.com

1. Arduino IDEのダウンロード

こちらのページから。

www.arduino.cc

Macの場合は、ダウンロード・解凍したら、Aruduino appをアプリケーションフォルダに移して完了。

2. ドライバのインストール

こちらのページから

jp.silabs.com

CP210X Driverをダウンロードし、インストールします。

3. ボードマネージャとライブラリの追加

「Arduino」->「Preferences」で環境設定画面を開く。

「追加のボードマネージャのURL」の右側のボタンをクリック。

https://dl.espressif.com/dl/package_esp32_index.json

を追加してOK。

「ツール」->「ボード」->「ボードマネージャ」

「esp32」で検索し、ESP32用のライブラリをインストール。

「ツール」 -> 「ボード」 -> 「M5Stick-C」

でM5Stick-Cを選択。

次に、

「スケッチ」->「ライブラリをインクルード」->「ライブラリを管理」

「M5StickC」で検索し、M5StickCのライブラリをインストール。

試しにHello Worldをしてみます。

M5Stick-CをUSBでMacに繋いで、

「ツール」 -> 「シリアルポート」

で新たに現れたシリアルポートを選択。

「ファイル」->「スケッチ例」->「M5StickC」 ->「Basics」 ->「HelloWorld」

でHelloWorldのコードを開き、「検証」と「マイコンボードに書き込む」を実行。

4. ENV Hatを動かしてみる

「ファイル」->「スケッチ例」->「M5StickC」->「Hat」->「ENV」

コードに

/*
    note: need add library Adafruit_BMP280 from library manage
    Github: https://github.com/adafruit/Adafruit_BMP280_Library
*/

とあるので、

「スケッチ」->「ライブラリをインクルード」->「ライブラリを管理」

BMP280で検索して、by Adafruitのライブラリをインストール。

「検証」と「マイコンボードに書き込む」を実行。

これで実行されるはずです。

Raspberry Pi 4&Raspbian BusterでOpenCV-Pythonをインストール

Raspberry Pi 4の4GBモデルがとうとう日本でも発売開始になりましたね！

jp.rs-online.com

www.switch-science.com

これまでのモデルとの違いなどはいろんなところで比較されているので、割愛します。

Macを使ったOSのインストール方法は以下のページを参考にしました。

qiita.com

このブログでは、Raspbian Busterを入れたRaspberry Pi 4でのOpenCV-Pythonのインストール方法をメモしていきます。

参考にしたのは以下のページ。

www.pyimagesearch.com

こちらのページではOpenCV-Python 4のインストール方法が紹介されていますが、今回は諸事情でversion 3をインストールしました。手順は全く同じです。

1. 依存パッケージのインストール

LXTerminalで下記のコードを順次実行。

sudo apt-get update && sudo apt-get upgrade
sudo apt-get install build-essential cmake pkg-config
sudo apt-get install libjpeg-dev libtiff5-dev libjasper-dev libpng-dev
sudo apt-get install libavcodec-dev libvformat-dev libswscale-dev libv4l-dev
sudo apt-get install libxvidcore-dev libx264-dev
sudo apt-get install libfontconfig1-dev libcairo2-dev
sudo apt-get install libgdk-pixbuf2.0-dev libpango1.0-dev
sudo apt-get install libgtk2.0-dev libgtk-3-dev
sudo apt-get install libatlas-base-dev gfortran
sudo apt-get install libhdf5-dev libhdf5-serial-dev libhdf5-103
sudo apt-get install libqtgui4 libqtwebkit4 libqt4-test python3-pyqt5
sudo apt-get install python3-dev

実際の作業では、上記を１つの.shファイルにして一括で実行しました。時々yes/no?と聞かれるので、その都度、yを押します。

2. pipとvirtualenvのインストール

次に、pipとvirtualenvをインストールしていきますが、楽するために、作業に入る前にパスの設定ファイルを準備しておきました。

nano opencvInstallpath

export WORKON_HOME=$HOME/.virtualenvs
export VIRTUALENVWRAPPER_PYTHON=/usr/bin/python3
source /usr/local/bin/virtualenvwrapper.sh

次に、以下を順次実行。

wget https://bootstrap.pypa.io/get-pip.py
sudo python get-pip.py
sudo python3 get-pip.py
sudo rm -rf ~/.cache/pip
sudo pip install virtualenv virtualenvwrapper
cat opencvInstallpath >> ~/.bashrc
pip install "picamera[array]" #if you use picamera

最後の行は、picameraを使うときだけ必要。

終わったら、

source ~/.bashrc

でパスを有効にします。

次に、virtualenv環境を構築します。

mkvirtualenv cv -p python3

cv環境が作られ、その中に移動します。

(cv) pi@raspberrypi:~ $

cv環境から抜けるには

deactivate

を実行します。

3. OpenCV-Pythonのインストール

cv環境中で、opencv-pythonをインストールします。

pip install opencv-contrib-python==3.4.4.19

これでうまくいきました。

ちなみに、いくつかのバージョンはエラーが出てしまってインストールできませんでした。（今回は試したのは3.3.0.10と3.4.7.28）

$ pip install opencv-contrib-python==3.4.7.28
Looking in indexes: https://pypi.org/simple, https://www.piwheels.org/simple
Collecting opencv-contrib-python==3.4.7.28
  Downloading https://www.piwheels.org/simple/opencv-contrib-python/opencv_contrib_python-3.4.7.28-cp37-cp37m-linux_armv7l.whl (15.3MB)
     |████████████▎                   | 5.9MB 42kB/s eta 0:03:43
ERROR: THESE PACKAGES DO NOT MATCH THE HASHES FROM THE REQUIREMENTS FILE. If you have updated the package versions, please update the hashes. Otherwise, examine the package contents carefully; someone may have tampered with them.

このあたり、詳しい情報をお持ちの方がいらっしゃいましたら、コメントいただけると嬉しいです。

最後に、opencv-pythonが正しくインストールされたか確認します。

$ python3
Python 3.7.3 (default, Apr  3 2019, 05:39:12) 
[GCC 8.2.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import cv2
>>>

エラーが出なければOKです。

以上です。

[2021/4/15追記] virtualenv環境をアクチベートするには

$ workon [env-name]

とします。

Raspberry Pi 7インチ タッチ・スクリーン ディスプレイ 設定覚書

公式のタッチ・スクリーン ディスプレイを購入したのですが、設定に苦戦したので覚書として。

raspberry-pi.ksyic.com

Rasbianのバージョン：Buster

主な問題点は、フレームを装着した状態だと画面が逆向きになる、太い黒縁が出てタッチした場所とカーソルが反応する場所が違う、など。

LXTerminalでconfig.txtの設定を変更します。

sudo nano /boot/config.txt

色々なページを参考に変えていって、最終的にうまくいったバージョンのメモ。

#hdmi_safe=1

disable_overscan=1

#次は全てコメントアウト
#overscan_left=32
#overscan_right=32
#overscan_top=16
#overscan_bottom=16

#framebuffer_width=1280
#framebuffer_height=960

#hdmi_force_hotplug=1

hdmi_group=2
hdmi_mode=87
hdmi_cvt=480 320 60 1 0 0 0

hdmi_drive=2

config_hdmi_boost=4

#中略

#最後に追記
lcd_rotate=2

以上です。

Raspberry Piストアができたらしい

TechCrunchの記事より。

jp.techcrunch.com

店構えがステキ。

記事によるとイギリスのケンブリッジだけらしいですが、日本にもできるといいなあ、、、

www.raspberrypi.org

The Raspberry Pi Store is a place where you can experience and buy Raspberry Pi products. Explore some of the things you can do with a Pi, discover our accessories and books, and get your hands on store-only exclusives.

Raspberry Pi Storeは、Raspberry Pi製品を体験して購入できる場所です。 Piを使ってできることのいくつかを探り、私たちのアクセサリーや本を見つけて、店内の限定品を手に入れよう。

Pythonで任意の四角形を正方形に変換する方法（射影変換・ホモグラフィ）

Python 画像処理

以前、台形補正についての記事を書きました。

plant-raspberrypi3.hatenablog.com

今回は、台形だけではなく、任意の四角形を正方形（または別の任意の四角形）に変換しよう、という試み。

射影変換の原理を勉強して、Pythonで実装してみました。

PythonライブラリのSympyとScikit-imageを使う方法です。

1. 射影変換の原理

クリックすると展開されます

射影変換（ホモグラフィ）の原理にはついて、たくさんの解説記事が見つかります。

shogo82148.github.io

yaju3d.hatenablog.jp

この記事では、私が自分なりに理解したことについて書いていきます。

・ステップ1

ベクトル (x,y)をベクトル (u,v)にうつす線形写像 Hを求める問題を考えます。

f:id:plant-raspberrypi3:20181118185617p:plain:w400

これは

{\displaystyle \mathbf{a} = \begin{pmatrix} x \\ y \\ 1 \end{pmatrix} \\ }
{\displaystyle \mathbf{b} = w \begin{pmatrix} u \\ v \\ 1 \end{pmatrix} = \begin{pmatrix} uw \\ vw \\ w \end{pmatrix} \\ }
{\displaystyle H = \begin{pmatrix} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & 1 \end{pmatrix} }

とおくと

\mathbf{b} = H \mathbf{a}

が成り立つ Hを求める問題と置き換えることができます。

・ステップ2

次に、四隅の座標を基準にして変形する場合を考えます。変換前の座標を

(x_{n}, y_{n})

変換後の座標を

(u_{n}, v_{n})

とします。ただし、 n = 1, 2, 3, 4 です。

具体的には、以下の図のような変換を考えます。

f:id:plant-raspberrypi3:20181118190217p:plain:w400

ここで、四隅の座標が与えられている場合、各 nについて以下の式が成り立ちます。

u_{n} w_{n} = h_{11} x_{n} + h_{12} y_{n} + h_{13}
v_{n} w_{n} = h_{21} x_{n} + h_{22} y_{n} + h_{23}

これは、 \mathbf{b} = H \mathbf{a}を計算することで得られます。

また、同様に

w_{n} = h_{31} x_{n} + h_{32} y_{n} + 1

も得られます。これを上の2つの式に代入して左辺を右辺に移動すると

h_{11} x_{n} + h_{12} y_{n} + h_{13} - u_{n} (h_{31} x_{n} + h_{32} y_{n} + 1) = 0
h_{21} x_{n} + h_{22} y_{n} + h_{23} - v_{n} (h_{31} x_{n} + h_{32} y_{n} + 1) = 0

となるので、これを n = 1, 2, 3, 4について解くことで h_{11}〜 h_{32}が得られ、 Hが求まります。

求めた Hを使い、

\mathbf{b} = H\mathbf{a}

の関係から画像の変換を行います。

2. Sympyで連立方程式を解きHを求める

クリックすると展開されます

前述の通り、 Hを求めるためには8元連立方程式を解く必要があります。

しかし、手計算で解くのはとても面倒なので、Sympyで解いていきます。

Sympyの使い方についてはこの辺りを参照。

3.2. Sympy : Python での代数計算 — Scipy lecture notes

環境は、MacOS Sierra, Python3.6です。 まずは必要なパッケージをインポートします。

from sympy import solve,symbols
import numpy as np

h_{11}〜 h_{32}を変数として設定します。

str_list = [f"h{i}" for i in [11,12,13,21,22,23,31,32]]
str_symbol = " ".join(str_list)

h11, h12, h13, h21, h22, h23, h31, h32 = symbols(str_symbol)

次に、 (x_{n}, y_{n})のセットをnumpyアレイxyに、 (u_{n}, v_{n})のセットをnumpyアレイuvに代入します。

Hを解くための連立方程式を以下のように記述します。

各 nごとに方程式を求め、リストexpに加えます。

exp = []
for i in range(4):
    exp_x = h11*xy[i,0]+h12*xy[i,1]+h13-(h31*xy[i,0]+h32*xy[i,1]+1)*uv[i,0]
    exp_y = h21*xy[i,0]+h22*xy[i,1]+h23-(h31*xy[i,0]+h32*xy[i,1]+1)*uv[i,1]
    exp.append(exp_x)
    exp.append(exp_y)

全ての方程式の入ったリストexpを引数としてSympyのsolve関数に渡し、連立方程式を解き、解をHに代入します。

result = solve(exp,dict=True)
r = result[0]

H = np.array([[r[h11],r[h12],r[h13]],
                           [r[h21],r[h22],r[h23]],
                           [r[h31],r[h32],1]],dtype=np.float64)

仮に以下のような図の変換の場合

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

方程式は

1: 10.0*h11 + 20.0*h12 + h13
2: 10.0*h21 + 20.0*h22 + h23
3: 20.0*h11 + 70.0*h12 + h13
4: 20.0*h21 + 70.0*h22 + h23 - 2000.0*h31 - 7000.0*h32 - 100.0
5: 90.0*h11 + 90.0*h12 + h13 - 9000.0*h31 - 9000.0*h32 - 100.0
6: 90.0*h21 + 90.0*h22 + h23 - 9000.0*h31 - 9000.0*h32 - 100.0
7: 90.0*h11 + 10.0*h12 + h13 - 9000.0*h31 - 1000.0*h32 - 100.0
8: 90.0*h21 + 10.0*h22 + h23

となり、Sympyによって Hを求めると

array([[ 2.06775593e+00, -4.13551186e-01, -1.24065356e+01],
       [ 2.13237330e-01,  1.70589864e+00, -3.62503462e+01],
       [ 8.18794424e-03, -4.13551186e-03,  1.00000000e+00]])

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

この Hが正しいかは、 (x_{n}, y_{n})の値を用いて H\mathbf{a}を計算して確かめることができます。

for i in range(4):
    xy1 = np.append(xy[i],1)
    uv1 = np.dot(H,xy1)
    uv1 = np.round(uv1/uv1[-1]).astype(int)
    uv1 = uv1[:-1]
    print(f"{xy[i]} -> {uv1}")

#出力
[20 10] -> [0 0]
[70 20] -> [100   0]
[90 90] -> [100 100]
[10 90] -> [  0 100]

3. 求めたHとScikit-imageを使って図形を変換する

クリックすると展開されます

Hを図形に適用して、実際に変換していきます。

変換がわかりやすいように多角形を描いてみました。

polygon = np.array([[20,10],[50,40],[70,20],[90,90],[50,60],[10,90]])

img = np.zeros((100,100),dtype=np.float64)
rr,cc = draw.polygon(polygon.T[0],polygon.T[1])
img[rr,cc] = 1

plt.imshow(img)
plt.show()

f:id:plant-raspberrypi3:20181118203324p:plain:w400

これを変換していきます。

前の項目で述べた、 Hを求める過程を関数にまとめます。

from sympy import solve,symbols
import numpy as np
from skimage import draw,transform
import matplotlib.pyplot as plt
from pprint import pprint

def solve_H(xy,uv):

    print("XY to UV")
    for i in range(4):
        print(f"{xy[i]} -> {uv[i]}")

    str_list = [f"h{i}" for i in [11,12,13,21,22,23,31,32]]
    str_symbol = " ".join(str_list)
    
    h11, h12, h13, h21, h22, h23, h31, h32 = symbols(str_symbol)

    exp = []

    for i in range(4):
        exp_x = h11*xy[i,0]+h12*xy[i,1]+h13-(h31*xy[i,0]+h32*xy[i,1]+1)*uv[i,0]
        exp_y = h21*xy[i,0]+h22*xy[i,1]+h23-(h31*xy[i,0]+h32*xy[i,1]+1)*uv[i,1]
        exp.append(exp_x)
        exp.append(exp_y)

    print("")
    print("Equations")
    for i,e in enumerate(exp):
        print(f"{i+1}: {e}")

    result = solve(exp,dict=True)
    r = result[0]

    H = np.array([[r[h11],r[h12],r[h13]],
                           [r[h21],r[h22],r[h23]],
                           [r[h31],r[h32],1]],dtype=np.float64)
    print("")
    print("Result_array")
    pprint(H)
    
    return H

この関数を使って解きます。四隅の変換前と変換後の座標を設定します。

xy = np.array([[20,10],[70,20],[90,90],[10,90]])
uv = np.array([[0,0],[100,0],[100,100],[0,100]])

早速、 Hを解いていきます。

が、そのままでやると何故かわからないですがうまくいかないことがわかったので、xとyを入れ替えたりゴソゴソやっています。numpyのxyと画像にした時のxyの方向が異なるせいかなぁ、と思っていますが、もし詳しい方がいらっしゃったらコメントいただけるとありがたいです。。。

xy_t = xy[:,::-1]
uv_t = uv[:,::-1]

H = solve_H(xy_t,uv_t)

#出力
XY to UV
[10 20] -> [0 0]
[20 70] -> [  0 100]
[90 90] -> [100 100]
[90 10] -> [100   0]

Equations
1: 10.0*h11 + 20.0*h12 + h13
2: 10.0*h21 + 20.0*h22 + h23
3: 20.0*h11 + 70.0*h12 + h13
4: 20.0*h21 + 70.0*h22 + h23 - 2000.0*h31 - 7000.0*h32 - 100.0
5: 90.0*h11 + 90.0*h12 + h13 - 9000.0*h31 - 9000.0*h32 - 100.0
6: 90.0*h21 + 90.0*h22 + h23 - 9000.0*h31 - 9000.0*h32 - 100.0
7: 90.0*h11 + 10.0*h12 + h13 - 9000.0*h31 - 1000.0*h32 - 100.0
8: 90.0*h21 + 10.0*h22 + h23

Result_array
array([[ 2.06775593e+00, -4.13551186e-01, -1.24065356e+01],
       [ 2.13237330e-01,  1.70589864e+00, -3.62503462e+01],
       [ 8.18794424e-03, -4.13551186e-03,  1.00000000e+00]])

最後に求めた Hを、skimageのtransformモジュールの機能を使って、画像に適用します。
使う機能はProjectiveTransformとwarpです。

ProjectiveTransformはinverseにすると想定通りの変換となります。

変換によって伸展した場所がモヤっとしてしまうので、最後に四捨五入によってバイナリ画像に修正してから表示します。

tform = transform.ProjectiveTransform(matrix=H).inverse
warped = transform.warp(img,tform,order=3,mode="edge")

plt.imshow(img)
plt.show()

plt.imshow(np.round(warped).astype(int))
plt.show()

変換前

f:id:plant-raspberrypi3:20181118200712p:plain:w400

変換後

f:id:plant-raspberrypi3:20181118200730p:plain:w400

うまくいきました！

原理的には正方形だけでなく、別の任意の四角形にも変換できるので、色々と面白いことができそうです。

以上です。

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

画像処理 数学 備忘録 scikit-image Python

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

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

長文です。

目次

  1. 前置き
  2. 二値画像で考えてみる
  3. さらに単純化してみる
  4. 「なぜ、M_{00}が領域の面積になるのか」の答え
  5. 「なぜ、(\frac{M_{10}}{M_{00}}, \frac{M_{01}}{M_{00}})が重心の座標になるのか」の答え
  6. numpyとscikit-imageで試す

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を見にいってみる。。。

Image moment - 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

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