inopage

2020/04/11
本記事は学校の課題をプログラムで楽してこなそうとした際のメモを編集したものです。レポートの面影はありません。

目次

• 運動方程式を求める
  • モデル
  • 式変形用コード
  • 運動方程式
• シミュレーションしてみる
  • 順運動学、逆運動学
  • 順逆動力学計算用コード
  • 特異点
• 雑記

運動方程式を求める

※参考1-1のままやって一致すれば正しく出来ました、ということにしています。

モデル

上図のような２リンクマニピュレータを考えます。G_0=(x_0,y_0), G_1=(x_0,y_1), m_0, m_1, l_0, l_1をそれぞれ第一リンク、第二リンクの重心、質量、長さとします。また、リンク端からそれぞれの重心までの長さをそれぞれl_g0, l_g1としています。

$$l_{gi}=\frac{1}{2}l_i 　(i=0,1)$$

ラグランジュの運動方程式に当てはめて運動方程式を求めます。

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

$$ \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$$

あるいは、

$$ \frac{d}{dt}(\frac{\partial T}{\partial\dot q})-\frac{\partial T}{\partial q_i}+\frac{\partial U}{\partial q_i}+\frac{\partial D}{\partial \dot{q_i}}= Q_i　(i =1,2,..,f)$$

L:ラグランジアン

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

f:自由度

q:一般化座標

式変形用コード

ラグランジュ～でやると無限微分地獄で計算ミスのリスクが高いですね。少しの変更でものすごく時間がかかるのは嫌なので、Pythonで式変形する方法をやってみたいと思います。こういうときはmatlabっぽく書けるSimPyが便利 (雑記1.)。

• 式変形

simplify()

>from sympy import *
>x, y, z = symbols('x y z')
>init_printing(use_unicode=True)
>simplify(sin(x)**2 + cos(x)**2)
1
>simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1))
x - 1
>simplify(gamma(x)/gamma(x - 2))
(x - 2)⋅(x - 1)

式を単純化する関数で、計算量が不必要に大きくなることがあることに注意が必要です。神。

Simplification -SymPy1.5.1

• 解く

sp.solve (Eq)

引数に=0になる関数を入力します。

Solvers -SymPy 1.5.1

• 微分

diff(func, t)

第二引数に何で微分するのかを入れます。

Calculus -SymPy 1.5.1

参考=>1-2、1-3、1-4

• コード

環境:

• Python 3.7.3
• numpy 1.16.4
• sympy 1.4

https://github.com/inomatly/robot_arm/blob/master/solv_lagrange.ipynb

一応.pyも置いてありますが jupyter notebok で見ることを前提にしています。

式を単純化するだけのコードです。計算結果はTEXで出しています。

運動方程式

各リンクのx軸とのなす角τ_0, τ_1は

$$\tau_0=\theta_0 ...(1-1)$$

$$\tau_1=\theta_0+\theta_1 ...(1-2)$$

です。 また、第一リンクの重心位置(x_0, y_0)は

$$\begin{bmatrix}x_0 \cr y_0\end{bmatrix}=\begin{bmatrix}l_{g0}cos\theta_0\cr l_{g0}sin\theta_0\end{bmatrix} ...(1-3)$$

で表せます。また、第二リンクの重心位置(x_1,y_1)は

$$\begin{bmatrix}x_1 \cr y_1\end{bmatrix}=\begin{bmatrix}l_{g0}cos\theta_0+ l_{g1}cos(\theta_0+\theta_1)\cr l_{g0}sin\theta_0+ l_{g1}sin(\theta_0+\theta_1)\end{bmatrix} ...(1-4)$$

で表せる。各リンクの運動エネルギー、位置エネルギー、損失エネルギー、一般化力をT_0, T_1, U_0, U_1, D_0, D_1, Q_0, Q_1、一般化座標q_0, q_1として、ラグランジュの運動方程式

$$ \frac{d}{dt}(\frac{\partial T}{\partial\dot q})-\frac{\partial T}{\partial q_i}+\frac{\partial U}{\partial q_i}+\frac{\partial D}{\partial \dot{q_i}}= Q_i　(i =1,2) ...(1-5)$$

を用いて各リンクの運動方程式を求める。

ここで、一般化座標はq_i=θ_i (i=0,1)である。各リンクの慣性モーメントをI_0, I_1,関節の粘性摩擦係数をd_0, d_1としてT_0, T_1, U_0, U_1, D_0, D_1を求める。

$$T_0=\frac{1}{2} \left(I_{0} + l_{g0}^{2} m_{0}\right) \left(\frac{d}{d t} {q_{0}}{\left(t \right)}\right)^{2} ...(1-6)$$

$$U_0=g l_{g0} m_{0} \sin{\left({q_{0}}{\left(t \right)} \right)} ...(1-7)$$

$$D_0=\frac{d_{0} \left(\frac{d}{d t} {q_{0}}{\left(t \right)}\right)^{2}}{2} ...(1-8)$$

$$T_1=\frac{1}{2} I_{1} \left(\frac{d}{d t} {q_{0}}{\left(t \right)} + \frac{d}{d t} {q_{1}}{\left(t \right)}\right)^{2} \\+ \frac{1}{2} m_{1} \left(l_{0}^{2} \left(\frac{d}{d t} {q_{0}}{\left(t \right)}\right)^{2} + 2 l_{0} l_{g1} \cos{\left({q_{1}}{\left(t \right)} \right)} \left(\frac{d}{d t} {q_{0}}{\left(t \right)}\right)^{2} \\+ 2 l_{0} l_{g1} \cos{\left({q_{1}}{\left(t \right)} \right)} \frac{d}{d t} {q_{0}}{\left(t \right)} \frac{d}{d t} {q_{1}}{\left(t \right)} + l_{g1}^{2} \left(\frac{d}{d t} {q_{0}}{\left(t \right)}\right)^{2} \\+ 2 l_{g1}^{2} \frac{d}{d t} {q_{0}}{\left(t \right)} \frac{d}{d t} {q_{1}}{\left(t \right)} + l_{g1}^{2} \left(\frac{d}{d t} {q_{1}}{\left(t \right)}\right)^{2}\right) ...(1-9)$$

$$U_1=g m_{1} \left(l_{0} \sin{\left({q_{0}}{\left(t \right)} \right)} + l_{g1} \sin{\left({q_{0}}{\left(t \right)} + {q_{1}}{\left(t \right)} \right)}\right) ...(1-10)$$

$$D_1=\frac{d_{0} \left(\frac{d}{d t} {q_{1}}{\left(t \right)}\right)^{2}}{2} ...(1-11)$$

また、全体の運動エネルギー、位置エネルギー、損失エネルギーT, U, Dは

$$T=T_0+T_1 ...(1-12)$$

$$U=U_0+U_1 ...(1-13)$$

$$T=D_0+D_1 ...(1-14)$$

これより、

$$\frac{\partial T}{\partial \dot{q_0}}=I_{1} \left(\frac{d}{d t} {q_{0}}{\left(t \right)} + \frac{d}{d t} {q_{1}}{\left(t \right)}\right) + m_{1} \left(l_{0}^{2} \frac{d}{d t} {q_{0}}{\left(t \right)} \\+ 2 l_{0} l_{g1} \cos{\left({q_{1}}{\left(t \right)} \right)} \frac{d}{d t} {q_{0}}{\left(t \right)} + l_{0} l_{g1} \cos{\left({q_{1}}{\left(t \right)} \right)} \frac{d}{d t} {q_{1}}{\left(t \right)} + l_{g1}^{2} \frac{d}{d t} {q_{0}}{\left(t \right)} \\+ l_{g1}^{2} \frac{d}{d t} {q_{1}}{\left(t \right)}\right) +  \left(I_{0} + l_{g0}^{2} m_{0}\right) \frac{d}{d t} {q_{0}}{\left(t \right)}  ...(1-15)$$

$$\frac{d}{dt}(\frac{\partial T}{\partial \dot{q_0}})=I_{1} \left(\frac{d^{2}}{d t^{2}} {q_{0}}{\left(t \right)} + \frac{d^{2}}{d t^{2}} {q_{1}}{\left(t \right)}\right) + m_{1} \left(l_{0}^{2} \frac{d^{2}}{d t^{2}} {q_{0}}{\left(t \right)} \\- 2 l_{0} l_{g1} \sin{\left({q_{1}}{\left(t \right)} \right)} \frac{d}{d t} {q_{0}}{\left(t \right)} \frac{d}{d t} {q_{1}}{\left(t \right)} - l_{0} l_{g1} \sin{\left({q_{1}}{\left(t \right)} \right)} \left(\frac{d}{d t} {q_{1}}{\left(t \right)}\right)^{2} \\+ 2 l_{0} l_{g1} \cos{\left({q_{1}}{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} {q_{0}}{\left(t \right)} + l_{0} l_{g1} \cos{\left({q_{1}}{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} {q_{1}}{\left(t \right)} + l_{g1}^{2} \frac{d^{2}}{d t^{2}} {q_{0}}{\left(t \right)} \\+ l_{g1}^{2} \frac{d^{2}}{d t^{2}} {q_{1}}{\left(t \right)}\right) +  \left(I_{0} + l_{g0}^{2} m_{0}\right) \frac{d^{2}}{d t^{2}} {q_{0}}{\left(t \right)} ...(1-16)$$

$$\frac{\partial T}{\partial q_0}=0 ...(1-17)$$

$$\frac{\partial U}{\partial q_0}=g \left(l_{g0} m_{0} \cos{\left({q_{0}}{\left(t \right)} \right)} + m_{1} \left(l_{0} \cos{\left({q_{0}}{\left(t \right)} \right)} + l_{g1} \cos{\left({q_{0}}{\left(t \right)} + {q_{1}}{\left(t \right)} \right)}\right)\right) ..(1-18)$$

$$\frac{\partial D}{\partial \dot{q_0}}=d_{0} \frac{d}{d t} {q_{0}}{\left(t \right)} ...(1-19)$$

$$∴ Q_0=I_{1} \left(\frac{d^{2}}{d t^{2}} {q_{0}}{\left(t \right)} + \frac{d^{2}}{d t^{2}} {q_{1}}{\left(t \right)}\right) + d_{0} \frac{d}{d t} {q_{0}}{\left(t \right)} \\+ g \left(l_{g0} m_{0} \cos{\left({q_{0}}{\left(t \right)} \right)} + m_{1} \left(l_{0} \cos{\left({q_{0}}{\left(t \right)} \right)} + l_{g1} \cos{\left({q_{0}}{\left(t \right)} + {q_{1}}{\left(t \right)} \right)}\right)\right) \\+ m_{1} \left(l_{0}^{2} \frac{d^{2}}{d t^{2}} {q_{0}}{\left(t \right)} - 2 l_{0} l_{g1} \sin{\left({q_{1}}{\left(t \right)} \right)} \frac{d}{d t} {q_{0}}{\left(t \right)} \frac{d}{d t} {q_{1}}{\left(t \right)} - l_{0} l_{g1} \sin{\left({q_{1}}{\left(t \right)} \right)} \left(\frac{d}{d t} {q_{1}}{\left(t \right)}\right)^{2} \\+ 2 l_{0} l_{g1} \cos{\left({q_{1}}{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} {q_{0}}{\left(t \right)} + l_{0} l_{g1} \cos{\left({q_{1}}{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} {q_{1}}{\left(t \right)} + l_{g1}^{2} \frac{d^{2}}{d t^{2}} {q_{0}}{\left(t \right)} \\+ l_{g1}^{2} \frac{d^{2}}{d t^{2}} {q_{1}}{\left(t \right)}\right) +  \left(I_{0} + l_{g0}^{2} m_{0}\right) \frac{d^{2}}{d t^{2}} {q_{0}}{\left(t \right)} ...(1-20)$$

$$\frac{\partial T}{\partial \dot{q_1}}=I_{1} \left(\frac{d}{d t} {q_{0}}{\left(t \right)} + \frac{d}{d t} {q_{1}}{\left(t \right)}\right) \\+ l_{g1} m_{1} \left(l_{0} \cos{\left({q_{1}}{\left(t \right)} \right)} \frac{d}{d t} {q_{0}}{\left(t \right)} + l_{g1} \frac{d}{d t} {q_{0}}{\left(t \right)} + l_{g1} \frac{d}{d t} {q_{1}}{\left(t \right)}\right) ...(1-21)$$

$$\frac{d}{dt}(\frac{\partial T}{\partial \dot{q_1}})=I_{1} \left(\frac{d^{2}}{d t^{2}} {q_{0}}{\left(t \right)} + \frac{d^{2}}{d t^{2}} {q_{1}}{\left(t \right)}\right) + l_{g1} m_{1} \left(- l_{0} \sin{\left({q_{1}}{\left(t \right)} \right)} \frac{d}{d t} {q_{0}}{\left(t \right)} \frac{d}{d t} {q_{1}}{\left(t \right)} \\+ l_{0} \cos{\left({q_{1}}{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} {q_{0}}{\left(t \right)} + l_{g1} \frac{d^{2}}{d t^{2}} {q_{0}}{\left(t \right)} + l_{g1} \frac{d^{2}}{d t^{2}} {q_{1}}{\left(t \right)}\right) ...(1-22)$$

$$\frac{\partial T}{\partial q_1}=-  l_{0} l_{g1} m_{1} \left(\frac{d}{d t} {q_{0}}{\left(t \right)} + \frac{d}{d t} {q_{1}}{\left(t \right)}\right) \sin{\left({q_{1}}{\left(t \right)} \right)} \frac{d}{d t} {q_{0}}{\left(t \right)} ...(1-23)$$

$$\frac{\partial U}{\partial q_1}=g l_{g1} m_{1} \cos{\left({q_{0}}{\left(t \right)} + {q_{1}}{\left(t \right)} \right)} ...(1-24)$$

$$\frac{\partial D}{\partial \dot{q_1}}=d_{0} \frac{d}{d t} {q_{1}}{\left(t \right)} ...(1-25)$$

$$∴ Q_1= I_{1} \frac{d^{2}}{d t^{2}} {q_{0}}{\left(t \right)} + I_{1} \frac{d^{2}}{d t^{2}} {q_{1}}{\left(t \right)} + d_{0} \frac{d}{d t} {q_{1}}{\left(t \right)} +\\ g l_{g1} m_{1} \cos{\left({q_{0}}{\left(t \right)} + {q_{1}}{\left(t \right)} \right)} +  l_{0} l_{g1} m_{1} \sin{\left({q_{1}}{\left(t \right)} \right)} \left(\frac{d}{d t} {q_{0}}{\left(t \right)}\right)^{2} +\\ l_{0} l_{g1} m_{1} \cos{\left({q_{1}}{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} {q_{0}}{\left(t \right)} +  l_{g1}^{2} m_{1} \frac{d^{2}}{d t^{2}} {q_{0}}{\left(t \right)} + l_{g1}^{2} m_{1} \frac{d^{2}}{d t^{2}} {q_{1}}{\left(t \right)} ...(1-26)$$

以上を行列式の形にすると、

$$\begin{bmatrix}Q_0 \cr Q_1\end{bmatrix}=\\ \begin{bmatrix}I_0+m_0l_{g0}^2+I_1+m_1(l_0^2+l_{g1}^2+2l_0l_{g1}cos\theta_1)& I_1+m_1(l_{g1}^2+l_0l_{g1}cos\theta_1)\cr I_1+m_1(l_{g1}^2+l_0l_{g1}cos\theta_1)& I_1+m_1 l_{g1}^2\end{bmatrix}\begin{bmatrix}\ddot{\theta_0} \cr \ddot{\theta_1}\end{bmatrix}\\+\begin{bmatrix}d_0& 0\cr 0& d_1\end{bmatrix}\begin{bmatrix}\dot{\theta_0} \cr \dot{\theta_1}\end{bmatrix}\\+\begin{bmatrix}-m_1l_0l_{g1}(2\dot{\theta_0}\dot{\theta_1}+\dot{\theta_1}^2)sin\theta_1& m_0gl_{g0}cos(\theta_0)+m_1g(l_{0}cos\theta_0+l_{g1}cos(\theta_0+\theta_1))\cr m_1l_0l_{g1}\dot{\theta_0}^2sin\theta_1& m_1gl_{g1}cos(\theta_0+\theta_1)\end{bmatrix} \\...(1-27)$$

参考資料と一致しました。OK!! (雑記2.)

simplyfy式変形後は順番がめちゃくちゃなので綺麗に行列の形に直すのはかなり面倒です。気をつけてください。(雑記3.)

ここまでの参考

• 1-1."ロボット工学の基礎 ―マニピュレータの運動学・動力学―", 岡山大学 自然科学研究科 見浪 護

• 1-2.Sympyでオイラーラグランジュ方程式を解く -qiita

• 1-3.Jupyter Notebook上でSymPyの数式とLaTeXコマンドを組み合わせて表示 -qiita

• 1-4.Pythonで方程式を解く方法（SciPy、ニュートン法、二分法による計算）-Monte Carlo Note

シミュレーション

動く図を作っていきます。

順運動学と逆運動学

第二リンクの先端位置(x_0,y_0)は

$$\begin{bmatrix}x_0 \cr y_0\end{bmatrix}=\begin{bmatrix}l_{0}cos\theta_0\cr l_{0}sin\theta_0\end{bmatrix} ...(2-1)$$

第二リンクの先端位置(x_1,y_1)は

$$\begin{bmatrix}x \cr y\end{bmatrix}=\begin{bmatrix}x_1 \cr y_1\end{bmatrix}=\begin{bmatrix}l_{0}cos\theta_0+ l_{1}cos(\theta_0+\theta_1)\cr l_{0}sin\theta_0+ l_{1}sin(\theta_0+\theta_1)\end{bmatrix}...(2-2)$$

これを解いていきます。移項して,

$$x-l_{0}cos\theta_0= l_{1}cos(\theta_0+\theta_1)...(2-3)$$

$$y-l_{0}sin\theta_0= l_{1}sin(\theta_0+\theta_1) ...(2-4)$$

どちらも両辺2乗し、2式を足し合わせます。

$$x^2+y^2+l^2_0-2l_0(xcos\theta_0+ysin\theta_0)=l^2_1...(2-5)$$

cos,sinに関してまとめると,

$$\sqrt{x^2+y^2}cos(\theta_0+\alpha)=\frac{x^2+y^2+l^2_0-l^2_1}{2l_0}...(2-6)$$

ただし,

$$tan( \alpha)=\frac{y}{x}...(2-7)$$

$$∴ \theta_0=-\alpha\pm cos^{-1}(\frac{x^2+y^2+l^2_0-l^2_1}{2l^2_0\sqrt{x^2+y^2}})...(2-8)$$

次にθ_1を求めます。

(2-3)/(2-4)より,

$$\frac{y-l_{0}sin\theta_0}{x-l_{0}cos\theta_0}=tan(\theta_0+\theta_1)...(2-9)$$

$$∴\theta_1=-\theta_0+tan^{-1}(\frac{y-l_{0}sin\theta_0}{x-l_{0}cos\theta_0})...(2-10)$$

• 参考

2-1.2リンクマニピュレータの軌道追従制御 -qiita

2-2.【Python】2リンクマニピュレータの逆運動学シミュレーション -西住工房

2-3.大分大学工学部福祉環境工学科メカトロニクスコース松尾研究室ゼミ資料"MATLAB による 2 リンクロボットマニピュレータ制御のシミュレーション"

順逆運動学計算用コード

• そも

%matplotlib inline 
%matplotlib notebook

って何ですか、ということ。jupyter内で図を見えてくれるinlineは良く使います。使わなくても出るような気がするけど、出ない時にはとりあえず書き加える感じ。

一方notebookは編集可能らしい。アイコンをタッチしてから右クリックでドラッグしたり、いろんな操作ができ図を弄れる、上の電源ボタンみたいなもので画像を固定できる。

=>参考2-1、2-2

• wedget

jupyter widget

公式のUsing~を見てみます。

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

@interact(amp=(0.1, 4.0, 0.1), omega=(0.1, 4.0, 0.1), phase=(-np.pi, np.pi, 0.1), 
          fn = {'sin': np.sin, 'cos': np.cos, 'tan': np.tan})
def h(amp=1.0, omega=1.0, phase=0.0, fn=np.sin):
    domain=[-np.pi, np.pi]
    x = np.linspace(domain[0], domain[1], 100)
    y  = amp * fn(omega * x + phase)
    plt.plot(x, y)
    plt.plot(-phase/omega, 0, 'or')
    plt.xlim(domain)
    plt.ylim([-4, 4])

x=(最小,最大,ステップ)はinteractでも設定でき、初期設定としてinteractを使用したい場合はデコレータとして使用出来る。辞書型で与えると選択肢に出来る。

• コード 2

https://github.com/inomatly/robot_arm/blob/master/arm_sim.ipynb

一応.pyも置いてありますが jupyter notebok で見ることを前提にしています。

動いている様子:

楽しい。でもめっちゃ重い。(雑記3.)

このとき、(2-9)の±によって解（すなわち姿勢）は2通りあります。たとえば(1,1)は間接が(1,0)にある場合と(0,1)にある場合が考えられます。今回は計算を簡単にする為に1通りに統一しています。

特異点に関して

ヤコビ行列=0になるような姿勢は特異姿勢と呼ばれ、その点は特異点と呼ばれる。

(1)式を微分するとヤコビ行列Jは

$$\begin{bmatrix}\dot{x} \cr \dot{y}\end{bmatrix}=\begin{bmatrix}-l_{0}sin\theta_0 - l_{1}sin(\theta_0+\theta_1) & - l_{1}sin(\theta_0+\theta_1) \cr l_{0}cos\theta_0+ l_{1}cos(\theta_0+\theta_1) & l_{1}cos(\theta_0+\theta_1)\end{bmatrix}\begin{bmatrix}\dot{\theta_0} \cr \dot{\theta_1}\end{bmatrix}...(1)$$

$$∴ J=\begin{bmatrix}-l_{0}sin\theta_0 - l_{1}sin(\theta_0+\theta_1) & - l_{1}sin(\theta_0+\theta_1) \cr l_{0}cos\theta_0+ l_{1}cos(\theta_0+\theta_1) & l_{1}cos(\theta_0+\theta_1)\end{bmatrix}$$

特異点では

$$det J = l_0l_1sin\theta_1=0$$

$$∴ \theta_1=n\pi$$

n=0,1,2。この場合、たとえば(x,y) = (-2/√2,-2/√2)の伸びきった姿勢のような姿勢や、(0,0)のような縮こまった姿勢を示します。それ以上伸びたり縮んだり出来ない出来ませんから、一般に特異点では自由度が下がってることになります。

目標点が特異点を通る場合は問題ありませんが、特異点付近を通るときはアームが急激に動くため注意が必要です。(雑記5.)=>2-3

参考 2

2.1 Jupyter Notebookにおけるmatplotlib -Pythonオンライン学習サービスPyQ(パイキュー)

2-2 jupyter notebookにmatplotlibを使ってグラフを描画する -山本隆の開発日誌

2-3「ロボット工学　ー機械システムのベクトル解析ー, 広瀬茂男著 裳華房 2003年」の10章

雑記

1. コードに関して、「pythonじゃなくてmatlabとかmathematicaでいいんじゃないの」とも思いますが。新しく覚えるおが面倒なのと、学生のうちは安価ですが将来的には金がかかりますのでなるべくオープンソースを使っていきたいという意思表示です。オープンソース贔屓です。判官贔屓的な。
   
2. 「まだ微分で消耗してるの？」正月休みの1/3に手でゴリゴリ計算して、参考にしたのが正しいっぽいことを確認してからコード書いたので二度手間ですが、今後はラグランジュ方程式を使いまわせるのでＯＫということにしました。TEXコードまで出してくれるの神ですね。禁断の果実という感じがします。
   
3. ラグランジュ方程式は使えるようになったかな（わからないけど）と思いますが、お恥ずかしながらベクトル的な解法に苦手意識があります。これからやっていきたいと思います。
   
4. シミュレーションは、jupyter上の位置の問題なのか、かくかくしますね。難しい。動かしているときめっちゃチカチカするのは何故なのか。計算が重いんだろうか？だれかプログラミングおしえて。
   
5. ひよこの真似やワクワクっ！な動作のように脇を開閉する動きはこの類い...ではないですね。人間の腕は自由度のバケモノですから。ワクワクの時の手首が特異点かと言えばそんなわけない。

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

目次

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

目的

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

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

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

手順

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

エッジ処理

輪郭を取得します。

cv2.Canny

Canny -openCV

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 より

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を見るのが早いでしょう。