python を使った楕円のプロットの仕方
別エントリで楕円の当てはめのスクリプトを書いているのですが，どうも楕円のプロットのさせ方がいまいちピンと来なかったのでまとめます．ここで最終的に解決しようとしている問題は
楕円の方程式として f(x, y) = 0 が与えられたときに，どうやって Python を使ってプロットするか？
です．
１．楕円の方程式
中心が (0, 0)，長軸の長さ 2a，短軸の長さ 2b の楕円の方程式は下記のように定義されます．
上記の式をみてもわかるように，y = f(x) の形ではなく，f(x, y) = const の形になっています．そのため，直線をプロットするときのように x を一定の範囲で変化させてその時の y の値を都度プロットするという形をとれません．これでは不便です．．．
２．楕円の中心が原点で，かつ長軸・短軸が x 軸と y 軸に平行な場合
このケースが一番単純な場合になります．この場合，数式は下記のような形になります．
媒介変数表示を用いたプロット
基本式のままではプロットするときに不便なので，媒介変数表示を使います．結局のところ方程式を満たすような変数表現を見つければいいだけなので，下記のように表現できます．
これで，t の値を 0 -> 2pi まで変化させたときの x と y の値を単純にプロットしていけば楕円のプロットができるようになります．下記，サンプルコードとその結果です．
import matplotlib.pyplot as plt
import numpy as np
import math
def generateVecFromEllipse(a, b):
t = np.linspace(0, 2 * math.pi, 101)
t = np.reshape(t, (t.shape[0], 1))
xVec = np.zeros((t.shape))
yVec = np.zeros((t.shape))
for i in range(t.shape[0]):
xVec[i] = a * math.cos(t[i])
yVec[i] = b * math.sin(t[i])
data = np.concatenate((xVec, yVec), axis=1)
return data
def plotData(dataOrg):
fig, ax = plt.subplots(ncols = 1, figsize=(10, 10))
plt.xlim(-10, 10)
plt.ylim(-10, 10)
ax.plot(dataOrg[:, 0], dataOrg[:, 1])
if __name__ == "__main__":
print("Ellipse Drawing Sample")
# 1. Generating Ellipse Point
dataEl = generateVecFromEllipse(3, 2)
# 2. Plot Data
plotData(dataEl)
プロット結果
３．長軸・短軸が座標軸に対して角度を持つ場合
で，問題はここからです．上記でプロットした楕円は中心が原点にピッタリそろっていて，長軸・短軸が座標２軸に対して平行ですが，こんな理想的な状況はほとんどありえません．．．よって，ここからは f(x, y) = const の楕円の方程式が与えられたときに，プロットする方法を説明します．まず，長軸・短軸に角度がある場合です．この場合，数式は下記のようになります．
この場合，楕円をプロットしたいときにはいくつかのステップを踏む必要があります．
１．楕円の長軸・短軸が何度傾いているか調べる．
２．楕円の長軸・短軸の長さを調べる．
３．一番基本となる楕円のプロット（２）で紹介した方法で，２で取得した長軸・短軸の楕円を描画．
４．３で作成した点群を１で取得した角度で回転させる．
上記の４ステップを一つづつ見ていきます．
１．楕円の長軸・短軸が何度傾いているか調べる．
ここでは，長軸・短軸の角度を求めます．一番簡単な方法は二次形式の係数行列 A をスペクトル分解する方法です．
直交変換行列 T が回転をあらわしているので，これが回転角になります．
２．楕円の長軸・短軸の長さを調べる．
１の結果より，下記のことが言えます．
３．一番基本となる楕円のプロット（２）で紹介した方法で，２で取得した長軸・短軸の楕円を描画．
基本となる楕円プロットと同じです．
４．３で作成した点群を１で取得した角度で回転させる．
１より求まったθを用いて，生成点群を回転させます．
プロット結果（最後のスクリプトを使って生成してます．）
４．長軸・短軸が座標軸に対して角度を持ち，かつ楕円の中心が原点でない場合．
最後のパターンです．二次元平面に書かれる楕円をすべて表現するにはこの方程式が必要です．導出は３とほぼ同じですが，点の変換に平行移動が入ってきます．つまり，，，
この場合も３の場合とほぼ同じですが，一次の項を平方完成する作業が増えます．
１．楕円の長軸・短軸が何度傾いているか調べる．
２．楕円の長軸・短軸の長さを調べる．
３．スペクトル分解で取得した回転行列を代入する．
４．一次の項を平方完成させ，原点からの平行移動ベクトル (a, b) を求める．
５．一番基本となる楕円のプロット（２）で紹介した方法で，２で取得した長軸・短軸の楕円を描画．
６．３で作成した点群を１で取得した角度で回転させる．
７．６で作成した点群を４で取得したベクトル分移動させる．
上記より，楕円の一般式が求まりました．下記スクリプトです．
# -*- coding: utf-8 -*-
from sympy import *
import sympy
import numpy as np
import matplotlib.pyplot as plt
import math
def generateVecFromEllipse(axis, center, T, rng=[0, 2*math.pi], num=101):
t = np.linspace(rng[0], rng[1], num)
t = np.reshape(t, (t.shape[0], 1))
xVec = np.zeros((t.shape))
yVec = np.zeros((t.shape))
for i in range(t.shape[0]):
xVec[i] = axis[0, 0] * math.cos(t[i])
yVec[i] = axis[1, 0] * math.sin(t[i])
dataTmp = np.concatenate((xVec, yVec), axis=1)
dataTmp = T * dataTmp.T
dataTmp = dataTmp.T
dataTmp[:, 0] = dataTmp[:, 0] + center[0, 0]
dataTmp[:, 1] = dataTmp[:, 1] + center[1, 0]
return dataTmp
def generateEllipseEqn(axis, center, T):
# 1. Creating Vector with Symbol
x_ = Symbol('x_')
y_ = Symbol('y_')
x_vec = np.matrix([[x_ - center[0, 0]], [y_ - center[1, 0]]])
# 2. Define Original Ellipse Matrix
# No tilt
ElMat = np.matrix([[1/axis[0, 0]**2, 0], [0, 1/axis[1, 0]**2]])
# 3. Rotate Coordinate
Tx_ = T.T * x_vec
# 4. Expand Ellipse Eqn.
fn = sympy.expand((Tx_.T * (ElMat * Tx_))[0, 0] - 1)
A = float(fn.coeff(x_, 2))
B = float(fn.coeff(x_*y_, 1))
C = float(fn.coeff(y_, 2))
D = float(fn.subs([(y_, 0)]).coeff(x_, 1))
E = float(fn.subs([(x_, 0)]).coeff(y_, 1))
F = float(fn.subs([(x_, 0), (y_, 0)]))
print("The ellipse eqn you are looking for is :")
print(fn)
return A, B, C, D, E, F
def getEllipseProperty(A, B, C, D, E, F):
if A < 0:
A = -A
B = -B
C = -C
D = -D
E = -E
F = -F
# Spectral Decomposition
M = np.matrix([[A, B/2], [B/2, C]])
lamdas, v = np.linalg.eigh(M)
# Diagonalize Coeffs Matrix
DiagA = v.T * M * v
# Apply translation for 1st order term.
tmp = np.matrix([D, E]) * v
# Calculate coefficient in rotated coords
AA = DiagA[0, 0]
BA = DiagA[0, 1] + DiagA[1, 0]
CA = DiagA[1, 1]
DA = tmp[0, 0]
EA = tmp[0, 1]
scale = F - DA**2 / (4*AA) - EA**2 / (4*CA)
# Normalize coeffs wrt to constant term.
AA = AA / abs(scale)
BA = BA / abs(scale)
CA = CA / abs(scale)
DA = DA / abs(scale)
EA = EA / abs(scale)
FA = abs(scale) / abs(scale)
# Ellipse Property Extraction
a = 1 / math.sqrt(abs(lamdas[0] / scale))
b = 1 / math.sqrt(abs(lamdas[1] / scale))
T = np.matrix([[v[0, 0], v[0, 1]], [v[1, 0], v[1, 1]]])
trans = v.T * np.matrix([[-DA/(2*AA)], [-EA/(2*CA)]])
valid = True
if AA * CA < 0:
valid = False
return valid, np.matrix([[a], [b]]), trans, T
def plotData(dataOrg, dataEst):
fig, ax = plt.subplots(ncols = 1, figsize=(10, 10))
plt.xlim(-10, 10)
plt.ylim(-10, 10)
ax.plot(dataOrg[:, 0], dataOrg[:, 1])
ax.plot(dataEst[:, 0], dataEst[:, 1])
if __name__ == "__main__":
print("Ellipse Drawing Sample")
# Define ellipse property.
axisXOrg = 3
axisYOrg = 1
centerX = 5
centerY = 3
rad = math.radians(120)
# Prepare arguments.
axisOrg = np.matrix([[axisXOrg], [axisYOrg]])
Rorg = np.matrix([[cos(rad), -sin(rad)], [sin(rad), cos(rad)]])
centOrg = np.matrix([[centerX], [centerY]])
# 0. Generating Data and Elilpse Equation Coefficient
dataOrg = generateVecFromEllipse(axisOrg, centOrg, Rorg)
A, B, C, D, E, F = generateEllipseEqn(axisOrg, centOrg, Rorg)
# 1. Extract Elilpse Property from Egn
valid, axis, centerEst, Rest = getEllipseProperty(A, B, C, D, E, F)
# 2. Generating Ellipse Point and Rotate
dataEst = generateVecFromEllipse(axis, centerEst, Rest)
# 3. Plot Data
plotData(dataOrg, dataEst)
実際に実行した結果が下記のようになります．青の楕円がオリジナルで，緑の楕円が楕円の方程式から上記スクリプトを使って生成した点群です．結局同じ楕円なので，きれいに重なります．
