コンテンツへスキップ

inomatyについて

電気電子を専攻している大学生。趣味は電子工作、プログラミング、合気道、ランニング、読書。 脳科学と心理学が好き。肉ばかり食べている。 twitter => @oritatamigasa1

画像から、全体の中心点と模様展の半径、角度を取得してファイルで出力する方法です。opencvを使用します。

画像から、全体の中心点と模様(濃い青丸)の半径、角度を取得したので、メモを修正して公開します。最終的に極座標グラフにします。

目次

  • 目的
  • 手順
  • 環境
  • コード
  • 解説

目的

最近、画像から半径と角度をたくさん取得する必要に駆られました。その作業はかなり手間なので自動化してやろうというのが目的です。

私の使用する画像をそのまま使うわけにはいかないので、似たような画像をいらすとやで拾ってきました。輪郭があんまり丸くなかったので加工しています。

螺旋の貝殻のイラスト

この画像から全体の中心点と模様(濃い青丸)の半径、角度を取得します。

手順

  1. グレースケール&位置座標を取得
  2. 円の取得の中心点の座標を取得
  3. 中心点を基準にして (rcosθ,rsinθ)で表現
  4. 各点(0~i)を(r, θ)で表記
  5. csvで出力

環境

  • Python 3.8.0
  • conda 4.8.0
  • opencv 4.1.2
  • matplotlib 3.1.2
  • numpy 1.17.3

コード

import cv2
import matplotlib.pyplot as plt
import numpy as np
import csv
import math
import sys

"""
あらかじめ同じフォルダに「white.png」と「読み込み対象.png」をおいておく必要がある。
芯の中心を原点にして各点(小円)の半径、角度を取得する。
"""
#input
input='rasen'
img = cv2.imread('input'+input+'.png',cv2.IMREAD_COLOR)#ndarray.グレイスケールで読み込みとエラー。
white = cv2.imread('white.png',0)#ndarray

#事前処理
img = cv2.medianBlur(img, 5) #メディアンフィルタを用いて画像を平滑化
cv2.imwrite("output/img.png",img)#debag1。outputフォルダに作成される。
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
cv2.imwrite("output/gray.png",gray)#debag2
edge = cv2.Canny(gray, 20, 40)
cv2.imwrite("output/edge.png",edge)#debag3
#中心円の取得
circles = cv2.HoughCircles(edge, cv2.HOUGH_GRADIENT, 1, 1200, param1=60, param2=20, minRadius = 800, maxRadius =1000)
#print(circles)
circles = np.uint16(np.around(circles))
x_c=np.int64(circles[0][0][0])
y_c=-np.int64(circles[0][0][1])
r_c=np.int64(circles[0][0][2])
center =[x_c,y_c]
norm_c= np.linalg.norm(center, ord=2)
for i in circles[0,:]:
    # 円を描く
    cv2.circle(img,(i[0],i[1]),i[2],(0,0,0),4)
#小円の取得
circles = cv2.HoughCircles(edge, cv2.HOUGH_GRADIENT, 1, 200,param1=60,param2=20, minRadius = 15 ,maxRadius = 50) #param1は色の違い、param2は形の雑さ
circles = np.uint16(np.around(circles))#整数になり、16bitに変換された配列になる。
j=0
for i in circles[0,:]:
    # 小円を描く
    j+=1
    cv2.circle(img,(i[0],i[1]),i[2],(0,0,0),4)
print(str(j)+"points")

#データの整形
#radiusは不要
radius=circles[:,:,2]
radius=np.array(radius[0,:])
x=circles[:,:,0]
x=x[0,:]
y=circles[:,:,1]
y=y[0,:]

#csvファイルの作成
xs=[]
ys=[]
points=[]
with open('output/'+input+'_x-y.csv', 'w') as f:
    writer = csv.writer(f, lineterminator='\n')
    writer.writerow(center)
    for i in range(len(x)):
        xs.append(np.int64(x[i]))
        ys.append(-np.int64(y[i]))
        points.append([xs[i],ys[i]])# points=[x,y]
        writer.writerow([xs[i],ys[i]])

# L2ノルム(距離)の取得とCSVファイルの生成
ds = []
Theta=[]
cos=[]
with open('output/'+input+'_d-Theta.csv', 'w') as f:
    writer = csv.writer(f, lineterminator='\n')
    writer.writerow(['Radius','Deg[deg]'])#ラベル
    for i in range(len(xs)):
        X = np.array([xs[i]-x_c,ys[i]-y_c])#中心からのベクトル
        norm = np.linalg.norm(X, ord=2)#中心からの距離
        ds.append(norm)
        cos.append(X[0]/norm)#cos = X[0]/SQRT(X[0]^2+X[1]^2)
        if ys[i]>y_c:#芯の中心より上側。
            Theta.append(360-math.degrees(math.acos(cos[i])))
        else:        #芯の中心より下側。
            Theta.append(math.degrees(math.acos(cos[i])))
        writer.writerow([ds[i],Theta[i]])

cv2.imwrite('output/'+input+'circle.png',img)# debag4

なお、角度は時計回りに作っています。

メモ:

画像は、色情報の2次元マトリックスになっており、画像処理の基本は左上から右上まで走査したあと一個下の行を走査するイメージです。取得した各点の順番は左上が先だと思います。

解説

1. グレースケール&xy座標で画像を取得

読み込み

cv2.imread(img,flag)
  • - cv2.IMREAD_COLOR : カラー画像として読み込む.画像の透明度は無視される.デフォルト値
  • - cv2.IMREAD_GRAYSCALE : グレースケール画像として読み込む
  • - cv2.IMREAD_UNCHANGED : アルファチャンネルも含めた画像として読み込む 上記のフラグを使う代わりに,単に1, 0, -1 の整数値を与えて指定することもできます.

画像を扱う -OpenCV より

フィルタリング

ここではあまり働いてませんが、やると精度が上がります。下の関数はC++(のはず)の記法ですが、引数の参考まで。

medianBlur(const Mat& src, Mat& dst, int ksize)

画像フィルタリング -openCV

グレイスケール

cv2.cvtColor(img, code)
  • code – 色空間の変換コード.説明を参照してください

色空間の変換 -OpenCV

debag2: グレイスケールまでの画像

エッジ処理

輪郭を取得します。

cv2.Canny

Canny -openCV

debag3 :エッジ検出までの画像

2. ハフ変換で円を検出する

ハフ変換は画像中の直線や円などを検出する操作です。ここでもpython用のサイトより見やすかったのでC++(のはず)へのリンクを張っておきます。

HoughCircles(Mat& image, vector<Vec3f>& circles, int method, double dp, double minDist, double param1=100, double param2=100, int minRadius=0, int maxRadius=0)¶

検出された円の (中心の x 座標, 中心の y 座標, 半径) のタプルを返します。

  • - image – 8ビット,シングルチャンネル,グレースケールの入力画像.
  • - circles – 検出された円を出力するベクトル.各ベクトルは,3要素の浮動小数点型ベクトル (x, y, radius) としてエンコードされます.
  • - method – 現在のところ, CV_HOUGH_GRADIENT メソッドのみが実装されています.基本的には 2段階ハフ変換 で,これについては Yuen90 で述べられています.
  • - dp – 画像分解能に対する投票分解能の比率の逆数.例えば, dp=1 の場合は,投票空間は入力画像と同じ分解能をもちます.また dp=2 の場合は,投票空間の幅と高さは半分になります.
  • - minDist – 検出される円の中心同士の最小距離.このパラメータが小さすぎると,正しい円の周辺に別の円が複数誤って検出されることになります.逆に大きすぎると,検出できない円がでてくる可能性があります.
  • - param1 – 手法依存の 1 番目のパラメータ. CV_HOUGH_GRADIENT の場合は, Canny() エッジ検出器に渡される2つの閾値の内,大きい方の閾値を表します(小さい閾値は,この値の半分になります).
  • - param2 – 手法依存の 2 番目のパラメータ. CV_HOUGH_GRADIENT の場合は,円の中心を検出する際の投票数の閾値を表します.これが小さくなるほど,より多くの誤検出が起こる可能性があります.より多くの投票を獲得した円が,最初に出力されます.
  • - minRadius – 円の半径の最小値.
  • - maxRadius – 円の半径の最大値.

cv::HoughCircles¶ -OpenCV より

出力は中心位置と半径です。単位は画素だと思われます。

メモ:

上記コードでは以下のように使用しています。

circles = cv2.HoughCircles(gray, cv2.HOUGH_GRADIENT, 1, 180, param1=60, param2=40, minRadius = 800, maxRadius =1100)

上のコードでは入力画像は2504*2607です.minRadius = 800なので、画像の半分以上を覆うような大きい円だけを検出しています。 param2は経験上、円の雑さに関するパラメータで、小さいほど歪な円を検出します。

円を描画

void circle(Mat& img, Point center, int radius, const Scalar& color, int thickness=1, int lineType=8, int shift=0)¶
  • - img – 円を描画する画像.
  • - center – 円の中心座標.
  • - radius – 円の半径.
  • - color – 円の色.
  • - thickness – 円の枠線の太さ.負の値の場合,円が塗りつぶされます.
  • - lineType – 円の枠線の種類, Line() の説明を参照してください.
  • - shift – 中心点の座標と半径の値において,小数点以下の桁を表すビット数

cv::circle -OpenCV より

debag4: ハフ変換までの画像

3. 円の中心点の座標を取得

各点は(x,y,r)を取得します。芯の位置ベクトル($x_c,y_c$)から、中心点からの各点の位置ベクトル(x_i,y_i)=(x-x_c,y-y_c)を求めます。x軸単位ベクトルを(x_e,y_e)=(1,0)とすると、内積から、

$$cos\theta = \frac{(x-x_c)x_e+(y-y_c)y_e}{\sqrt{x_c^2+y_c^2}\sqrt{(x-x_c)^2+(y-y_c)^2}}$$

$$cos\theta = \frac{x-x_c}{\sqrt{(x-x_c)^2+(y-y_c)^2}}$$

$$\theta = deg(arccos(\frac{x-x_c}{\sqrt{(x-x_c)^2+(y-y_c)^2}}))[deg]$$

【NumPy入門】ベクトルの大きさ(ノルム)を計算するnp.linalg.norm

np.linalg.norm(X, ord=2)

L2ノルム(ユークリッドノルム)は一般的な"長さ"です。

ノルムの意味とL1,L2,L∞ノルム -高校数学の美しい物語

グラフにして確認

結果が正しいかを確認します。アバウトでOKなので、画像とグラフを重ねていきます。

細かいことはかつて書きました=>[python]データのインポートと極座標系のグラフを表示- イノマタの趣味部屋

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

"""
CSVファイルから極座標系のグラフをプロットする
入力:(CSVファイル名, 出力画像ファイル名,半径,角度)
出力:出力画像ファイル名.png
"""
def polar(in_file,out_file,param_r,param_theta):    
    df = pd.read_csv("output/"+in_file)#ファイル読み込み
    r = df[param_r]
    deg = df[param_theta]
    radian = np.deg2rad(deg)#
    Len=len(r)

    fig = plt.figure(figsize=(10, 10))#新しい図面
    ax = fig.add_subplot(111, projection='polar')#軸の表示
    ax.set_theta_direction(-1)#時計回り
    #ax.plot(radian,r, '-o')
    ax.plot(radian,r, 'or')
    
    i=2
    #while i < Len:
    #    ax.plot([radian[i],radian[i-1]],[r[i], r[i-1]], 'xr-')
    #    i +=2
    plt.savefig(out_file)
    plt.show()

polar('0%_d-Theta.csv','0%.png','Radius','Deg[deg]')
取得した数値からグラフを作成

なお、画像は下記サイトで背景を透明にしました。

画像に透明色を設定する -peko step

はじめの画像と重ね合わせると、うまくグラフにできていることが確認できます。 取得したCSVは正確みたいですね。

グラフと元々の画像を重ねた画像

もともとの画像を修正したりエッジを強調したりすれはもっといい結果が得られるはずです。

雑記

パラメータの変更も自動化できれば良かったですが、できませんでした。

目的は達成したものの、手で測定するのと同じくらい時間がかかってしまったのも悲しいです。追加があれば効果が発揮できるんですが。

参考

画像処理の基本は、古い本ですがこの本が良いと思います。

コンピュータ画像処理 田村 秀行 (著) オーム社

他に、最新のアルゴリズムとか実用的な部分はOpenCVを見るのが早いでしょう。

「pythonによる制御工学入門」を読ませて頂きました。面白かったのでご紹介します。
pythonによる制御工学入門

pythonによる制御工学入門」を読ませて頂きました。面白かったのでご紹介します。

サポートページ =>『Pythonによる制御工学入門』サポートページ

期間:12日間程度(1日2~3時間)

特徴

内容

  • pythonの基礎
  • 制御モデル
    • 伝達関数モデル
    • 状態空間モデル
  • 制御対称のふるまい
    • 時間応答
    • 周波数応答
  • 制御系設計
    • PID制御
    • 2自由度制御
    • ゲインチューニング(限界感度法)
    • ゲインチューニング(モデルマッチング)
    • 状態フィードバック制御(極配置法)
    • 状態フィードバック制御(最適レギュレータ)
  • 安定性
    • 閉ループ
    • 開ループ
  • アドバンスドな制御系設計
    • オブザーバ
    • ロバスト制御

感想

私の事前知識:

古典制御の伝達関数やPID制御は授業で出てきたので何となくイメージ出来ていました(わかるとは言ってない)。現代制御は全くわかりませんでした。

pythonに関してはゆるふわで常用しています。

全体の感想

表題からわかるようにpythonですぐ動かせます。パラメータを弄って理解を深めたり、自分で設計するうえで真似すればに使えるため便利です。

実用や、制御を勉強するための入門として全体を俯瞰する目的で書かれていると思われます。そのため、後半は特にあっさりした説明になっています。

pythonの記法に関しては、(私は読み飛ばしましたがパッと見た感じ)必要な部分は書いてありそうです。pythonを使ったことのない方には多少大変かもしれません。

モジュールはNumpy,Matplotli,Scipy,Sympy,Python-control(, slycotもインストールする必要あり)。 Python-control (リンクは2019/12/25時点で最新バージョンのマニュアルへ)はMATLABっぽく書けるパッケージなので、MATLABを使ったことがある方に良いかもしれません。

モデル

最適レギュレータを下で紹介していますので、状態空間モデルに関して少し書きます。

伝達関数モデルではフィードバック系全体について、出力y(t)、入力u(t)に対してP(t)を

$$\frac{y}{u}=P(t)$$

の形で表します。それに対して状態空間モデルでは、内部の状態を全て観測して、

$$\dot{ \textbf{x}}(t)= \textbf{Ax}(t)+ \textbf{Bu}(t)$$ $$\textbf{y}(t)= \textbf{Cx}(t)+ \textbf{Du}(t)$$

の形で表します。例えば、このようなアームを考えてみます。

アームのモデル

トルクを$$\tau$$

粘性摩擦による力を$$\mu \dot{\theta}$$

重力加速度をgとすると、運動方程式は

$$\tau(t)=J\ddot{\theta}+\mu \dot{\theta}+Mglsin\theta$$

と書けます。入力u(t)、出力y(t)をそれぞれ

$$u(t)=\tau(t)$$

$$y(t)=\theta(t)$$

とし、角度 $\theta$ を十分小さいとして $sin(\theta)$ を

$$sin(y(t))=y(t)$$

と近似すると、

$$u(t)=J\ddot{y}(t)+\mu \dot{y}(t)+Mgly(t)$$

と書けます。これをもとに伝達関数が作れます。

さらに、状態空間モデルでは、

$$\textbf{x}(t)= \begin{bmatrix} \theta & \dot{\theta} \end{bmatrix}^T$$

と置いて

$$\dot{\textbf{x}}(t)= \begin{bmatrix} \dot{\theta} \cr \ddot{\theta} \end{bmatrix} =\begin{bmatrix} 0 & 1\cr -\frac{Mgl}{J} & -\frac{\mu}{J} \end{bmatrix}\textbf{x}(t) + \begin{bmatrix} 0 \cr \frac{1}{J} \end{bmatrix}u(t) $$

$$y(t)=\theta(t)=\begin{bmatrix} 1 & 0 \end{bmatrix} \textbf{x}(t) $$

となります。状態空間モデルの形にまとめると、

$$\dot{ \textbf{x}}(t)= \textbf{Ax}(t)+ \textbf{Bu}(t)$$ $$\textbf{y}(t)= \textbf{Cx}(t)+ \textbf{Du}(t)$$

$$\textbf{A}=\begin{bmatrix} 0 & 1\cr -\frac{Mgl}{J} & -\frac{\mu}{J} \end{bmatrix},\textbf{B}=\begin{bmatrix} 0 \cr \frac{1}{J} \end{bmatrix}, \textbf{C}=\begin{bmatrix} 1 & 0 \end{bmatrix},\textbf{D}=0$$

のように表せます。

ただし、状態空間モデルでは各パラメータが観測出来る必要があります。

序盤はかなりわかりやすく書かれています。古典制御と現代制御の基本的な部分はわかった気がします。

制御系設計

pythonで自分で動かせるという点で最もよかったのが設計の部分でした。教科書のグラフをサラッと見るよりわかった感がありました。

授業で使った教科書ではPID制御の理論をやったあとで、「なんで制御できるのかはわかったけど、kp,ki,kdの具体的な調べ方がわからん」と思ってたので、ゲインチューニングがコードと一緒に書いてあるのはかなり有り難かったです。

MATLABのようにマクローリン展開などが一瞬で出来ます。

(ここでちょっと思ったのが、MATLABのような有料ツールに近いことが無料で出来るのだから計算速度等が落ちているのではということ。実際のところはどうなのかは調べてません。)

個人的ハイライト:最適レギュレータ

p.176~状態空間モデルでのフィードバックゲインFを下のように与えます。

$$u=\textbf{Fx}$$

この F を"いい感じ"(雑)に選ぶ為に使われる手法の1つが最適レギュレータらしいです。

最適レギュレータでは、最適な値を求めるために、評価関数を最小にすることを考えます。

$$\textbf{Q}=\textbf{Q}^T>0,\textbf{R}=\textbf{R}^T>0$$

に対して評価関数

$$J=\int^\infty_0\textbf{x}(t)^T\textbf{Qx}(t)+u^T(t)\textbf{R}u(t)dt$$

を最小化するコントローラは

$$u=\textbf{F}_{opt}\textbf{x}$$

$$\textbf{F}_{opt}=\textbf{R}^{-1}\textbf{B}^T\textbf{P}$$

の形で得られます。ただし、

$$\textbf{P}=\textbf{P}^T>0$$

はリカッチ方程式

$$\textbf{A}^T\textbf{P}+\textbf{A}\textbf{P}+\textbf{P}\textbf{B}\textbf{R}^{-1}\textbf{B}\textbf{P}+\textbf{Q}=0$$

を満たす唯一の正定対象行列です。 なんでこうなるのかは詳しくは書いてませんが今後調べます。 また、Jの最小値は

$$\textbf{x}(0)^T\textbf{Px}(0)$$

です。

python上でMATLABの関数に似たものを使用できるpython_controlというものがあります。これで、一見難しそうな計算が一瞬の内に出来ます!!これはやばい。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline 
from control.matlab import *
# グラフの設定
def plot_set(fig_ax,*args):
    # x軸のラベルを1つ目を引数に
    fig_ax.set_xlabel(args[0])
    # y軸のラベルを2つ目を引数に
    fig_ax.set_ylabel(args[1])
    fig_ax.grid(ls=':')
    if len(args)==3:
        fig_ax.legend(loc=args[2])

A = [[0,1],[-4,-5]]
B = [[0],[1]]
C= np.eye(2)
D = np.zeros([2,1])
P= ss(A,B,C,D)#状態空間モデル
#最適レギュレータによるFの計算
Q = [[100,0],[0,1]]
R=1
X,E,F=care(P.A, P.B, Q, R)#フィードバックゲインF,リカッチ方程式の解X,閉ループ極Eを返す。
F=-F
#状態空間モデルでの応答
Acl = P.A+P.B*F
Pfb = ss(Acl, P.B, P.C, P.D)#状態空間モデル
Td = np.arange(0, 5, 0.01)#時間
Ud = 1*(Td>0)#時間Tdが正の範囲では1,負では0となるステップ信号。
X0 = [-0.3,0.4]#初期値
xst,t=step(Pfb,Td)#ステップ応答
xin,_=initial(Pfb,Td,X0)
x,_,_=lsim(Pfb,Ud,Td,X0)

fig, ax = plt.subplots(1,2,figsize=(6,2.3))
for i in [0,1]:
    ax[i].plot(t,x[:,i])
    ax[i].plot(t,xst[:,i], ls='-')
    ax[i].plot(t,xin[:,i], ls='-.')

plot_set(ax[0],'t','$x_1$')
plot_set(ax[1],'t','$x_2$')
plt.savefig('最適レギュレータによるフィードバック制御.png')
最適レギュレータのステップ応答。
( 青が応答、オレンジがゼロ状態応答、点線がゼロ入力応答。)

最適レギュレータの関数:

現代文明に感動しました。伝達関数モデルでのゲインチューニングの方法もすぐできて衝撃的でした。

次に読みたい

最後に参考文献としておすすめ本が載っているのもポイント高い。

解析力学キャンパス・ゼミ (馬場 敬之 (著), マセマ出版社)を読んだ感想です

解析力学キャンパス・ゼミ 馬場 敬之 (著), マセマ出版社

期間: 15日程度

特徴

単純な例を多く挙げつつきちんと導出や計算過程を書いてくれている

内容

  • 基礎事項
    • 回転行列
    • ラグランジュの運動方程式、ハミルトンの正準方程式を使用した練習問題
  • ラグランジュの運動方程式
    • オイラー角
    • オイラーの方程式
    • 変分原理
    • 仮想仕事とダランベールの原理
  • ハミルトンの正準方程式
    • 位相空間とトラジェクトリー
    • 正準変換
    • リウビルの定理
    • ポアソン括弧
    • 無限小変換

感想

事前知識:

力学の大学院授業を受講して何となく一般化座標のイメージが出来ていた。ラグランジュの運動方程式はわっぱりわかっていなかった。

全体

さすがマセマ。数学的に難しい部分をわかり安く説明していて、自分で導出できる。全体的にわかったような気がする(気がするというのは重要)。

練習問題はほぼなく、導出と例題がほとんど。

他の難しい顔した本を読む前に一回読むといいと思う。

ラグランジュの運動方程式

一般化座標で表記できるため、足が固定されていない移動ロボットの運動を考える上で重要。

ラグランジュの運動方程式:

$$ \frac{d}{dt}(\frac{\partial L}{\partial\dot q_i})-\frac{\partial L}{\partial q_i}= Q_i (i =1,2,..,f)$$

$$L=T-U$$

L:ラグランジアン

Q:非保存力による一般化力

f:自由度

q:一般化座標

一般化座標を用いて記述することで、ニュートンの運動方程式で考えなければならなかった外力や束縛力を無視して記述できる。これは素晴らしい!!

デカルト座標、極座標、球座標で考えた後、座標を一般化する。力学の記憶が定かではなくてもわかるように、自由度やコリオリの力など出てくる専門用語には必ず説明があった。

オイラー角がちょいイメージしにくかったが、他のどの本より丁寧に書いていたと思う。

かなり丁寧かつ簡単に書いている。最適降下線問題やサイクロイド曲線の問題程度までの難易度で、わかりやすい。

ハミルトンの正準方程式

ハミルトンの正準方程式:

$$ \frac{d q_i}{d t}=\frac{\partial H}{\partial p_i} $$ $$ \frac{d p_i}{d t}= - \frac{\partial H}{\partial q_i} $$

$$H=\Sigma ^f_{i=1}p_i \dot q_i - L$$

$$p_i=\frac{\partial L}{\partial \dot q_i} (i =1,2,..,f)$$

H:ハミルトニアン

q_i:一般化座標

p_i:一般化運動量

トラジェクトリー(位相空間内の代表点が描く軌跡みたいなもの)がイメージしにくく、大変だった。

個人的ハイライト

p164~ :単振動のトラジェクトリーは楕円形、という話から E=$h\nu$ ,光量子のエネルギーに触れたところが面白かった。計算量は増えるが、この考えから統計力学や量子力学に発展すると言うのが面白い。

関連

他のいろんな力学の本と合わせて読みたい。

次に読みたい

統計力学キャンパス・ゼミ馬場 敬之 (著) マセマ出版社

解析力学 細谷曉夫著 -やまなみ書房 <=無料で公開してくださっている神様ですありがとうございます!!!!