質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

89.09%

原点を中心とした楕円の近似をしたいです.

解決済

回答 1

投稿

  • 評価
  • クリップ 0
  • VIEW 592
退会済みユーザー

退会済みユーザー

今あるプログラムは“x”と“y”の項があるため,原点を通らない近似となってしまいます.

近似式を
ax^2+bxy+cy^2=1
となるようにしたいです.

以下にプログラムを載せておきます.
わかるかたおねがいします.

入力データは【data.txt】の列を読み取るようになっています.

お願いします.

import numpy
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

"""Demonstration of least-squares fitting of ellipses
    __author__ = "Ben Hammel, Nick Sullivan-Molina"
    __credits__ = ["Ben Hammel", "Nick Sullivan-Molina"]
    __maintainer__ = "Ben Hammel"
    __email__ = "bdhammel@gmail.com"
    __status__ = "Development"
    Requirements
    ------------
    Python 2.X or 3.X
    numpy
    matplotlib
    References
    ----------
    (*) Halir, R., Flusser, J.: 'Numerically Stable Direct Least Squares
        Fitting of Ellipses'
    (**) http://mathworld.wolfram.com/Ellipse.html
    (***) White, A. McHale, B. 'Faraday rotation data analysis with least-squares
        elliptical fitting'
"""

class LSqEllipse:

    def fit(self, data):
        """Lest Squares fitting algorithm
        Theory taken from (*)
        Solving equation Sa=lCa. with a = |a b c d f g> and a1 = |a b c>
            a2 = |d f g>
        Args
        ----
        data (list:list:float): list of two lists containing the x and y data of the
            ellipse. of the form [[x1, x2, ..., xi],[y1, y2, ..., yi]]
        Returns
        ------
        coef (list): list of the coefficients describing an ellipse
           [a,b,c,d,f,g] corresponding to ax**2+2bxy+cy**2+2dx+2fy+g
        """
        x, y = numpy.asarray(data, dtype=float)

        #Quadratic part of design matrix [eqn. 15] from (*)
        D1 = numpy.mat(numpy.vstack([x**2, x*y, y**2])).T
        #Linear part of design matrix [eqn. 16] from (*)
        D2 = numpy.mat(numpy.vstack([x, y, numpy.ones(len(x))])).T

        #forming scatter matrix [eqn. 17] from (*)
        S1 = D1.T*D1
        S2 = D1.T*D2
        S3 = D2.T*D2

        #Constraint matrix [eqn. 18]
        C1 = numpy.mat('0. 0. 2.; 0. -1. 0.; 2. 0. 0.')

        #Reduced scatter matrix [eqn. 29]
        M=C1.I*(S1-S2*S3.I*S2.T)

        #M*|a b c >=l|a b c >. Find eigenvalues and eigenvectors from this equation [eqn. 28]
        eval, evec = numpy.linalg.eig(M)

        # eigenvector must meet constraint 4ac - b^2 to be valid.
        cond = 4*numpy.multiply(evec[0, :], evec[2, :]) - numpy.power(evec[1, :], 2)
        a1 = evec[:, numpy.nonzero(cond.A > 0)[1]]

        #|d f g> = -S3^(-1)*S2^(T)*|a b c> [eqn. 24]
        a2 = -S3.I*S2.T*a1

        # eigenvectors |a b c d f g>
        self.coef = numpy.vstack([a1, a2])
        self._save_parameters()

    def _save_parameters(self):
        """finds the important parameters of the fitted ellipse

        Theory taken form http://mathworld.wolfram
        Args
        -----
        coef (list): list of the coefficients describing an ellipse
           [a,b,c,d,f,g] corresponding to ax**2+2bxy+cy**2+2dx+2fy+g
        Returns
        _______
        center (List): of the form [x0, y0]
        width (float): major axis
        height (float): minor axis
        phi (float): rotation of major axis form the x-axis in radians
        """

        #eigenvectors are the coefficients of an ellipse in general form
        #a*x^2 + 2*b*x*y + c*y^2 + 2*d*x + 2*f*y + g = 0 [eqn. 15) from (**) or (***)
        a = self.coef[0,0]
        b = self.coef[1,0]/2.
        c = self.coef[2,0]
        d = self.coef[3,0]/2.
        f = self.coef[4,0]/2.
        g = self.coef[5,0]

        #finding center of ellipse [eqn.19 and 20] from (**)
        x0 = (c*d-b*f)/(b**2.-a*c)
        y0 = (a*f-b*d)/(b**2.-a*c)

        #Find the semi-axes lengths [eqn. 21 and 22] from (**)
        numerator = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
        denominator1 = (b*b-a*c)*( (c-a)*numpy.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
        denominator2 = (b*b-a*c)*( (a-c)*numpy.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
        width = numpy.sqrt(numerator/denominator1)
        height = numpy.sqrt(numerator/denominator2)

        # angle of counterclockwise rotation of major-axis of ellipse to x-axis [eqn. 23] from (**)
        # or [eqn. 26] from (***).
        phi = .5*numpy.arctan((2.*b)/(a-c))

        self._center = [x0, y0]
        self._width = width
        self._height = height
        self._phi = phi

    @property
    def center(self):
        return self._center

    @property
    def width(self):
        return self._width

    @property
    def height(self):
        return self._height

    @property
    def phi(self):
        """angle of counterclockwise rotation of major-axis of ellipse to x-axis
        [eqn. 23] from (**)
        """
        return self._phi

    def parameters(self):
        return self.center, self.width, self.height, self.phi


def make_test_ellipse(center=[1,1], width=1, height=.6, phi=3.14/5):
    """Generate Elliptical data with noise

    Args
    ----
    center (list:float): (<x_location>, <y_location>)
    width (float): semimajor axis. Horizontal dimension of the ellipse (**)
    height (float): semiminor axis. Vertical dimension of the ellipse (**)
    phi (float:radians): tilt of the ellipse, the angle the semimajor axis
        makes with the x-axis
    Returns
    -------
    data (list:list:float): list of two lists containing the x and y data of the
        ellipse. of the form [[x1, x2, ..., xi],[y1, y2, ..., yi]]
    """
    t = numpy.linspace(0, 2*numpy.pi, 1000)
    x_noise, y_noise = numpy.random.rand(2, len(t))

    ellipse_x = center[0] + width*numpy.cos(t)*numpy.cos(phi)-height*numpy.sin(t)*numpy.sin(phi) + x_noise/2.
    ellipse_y = center[1] + width*numpy.cos(t)*numpy.sin(phi)+height*numpy.sin(t)*numpy.cos(phi) + y_noise/2.

    return [ellipse_x, ellipse_y]

import numpy as np
import matplotlib.pyplot as plt
import ellipses as el
from matplotlib.patches import Ellipse

alpha = 5
beta = 3
N = 500
DIM = 2

fig = plt.figure(1)
ax = fig.add_subplot(111)




np.random.seed(2)

# Generate random points on the unit circle by sampling uniform angles
theta = np.random.uniform(0, 2*np.pi, (N,1))
eps_noise = 0.2 * np.random.normal(size=[N,1])
circle = np.hstack([np.cos(theta), np.sin(theta)])

# Stretch and rotate circle to an ellipse with random linear tranformation
B = np.random.randint(-3, 3, (DIM, DIM))
noisy_ellipse = np.loadtxt("data.txt")

# Extract x coords and y coords of the ellipse as column vectors
X = noisy_ellipse[:,0:1]
Y = noisy_ellipse[:,1:]

# Formulate and solve the least squares problem ||Ax - b ||^2
A = np.hstack([X**2, X * Y, Y**2, X, Y])
b = np.ones_like(X)
x = np.linalg.lstsq(A, b)[0].squeeze()

# Print the equation of the ellipse in standard form
print('The ellipse is given by {0:.3}x^2 + {1:.3}xy+{2:.3}y^2+{3:.3}x+{4:.3}y = 1'.format(x[0], x[1],x[2],x[3],x[4]))

# Plot the noisy data
plt.scatter(X, Y, label='Data Points', color='g')
  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

回答 1

checkベストアンサー

0

今あるプログラムは“x”と“y”の項があるため,原点を通らない近似となってしまいます.

単純に x と y の項を消して、近似すればいいのではないでしょうか。

質問のコードからの変更点

a x^2 + b xy + c y^2 = 1 のモデルとして近似する。

A = np.hstack([X ** 2, X * Y, Y ** 2])
b = np.ones_like(X)
x = np.linalg.lstsq(A, b)[0].squeeze()

# Print the equation of the ellipse in standard form
print(
    "The ellipse is given by {0:.3}x^2 + {1:.3}xy+{2:.3}y^2 = 1".format(
        x[0], x[1], x[2]
    )
)

コード全体

楕円の点を作成するのに OpenCV を使っていますが、本題には関係ないので、その部分は CSV から読み込むように置き換えてください。

import cv2
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

noisy_ellipse = cv2.ellipse2Poly((0, 0), (50, 70), 45, 0, 360, 5)

# Extract x coords and y coords of the ellipse as column vectors
X = noisy_ellipse[:, 0:1]
Y = noisy_ellipse[:, 1:]

# Formulate and solve the least squares problem ||Ax - b ||^2
A = np.hstack([X ** 2, X * Y, Y ** 2])
b = np.ones_like(X)
x = np.linalg.lstsq(A, b)[0].squeeze()

# Print the equation of the ellipse in standard form
print(
    "The ellipse is given by {0:.3}x^2 + {1:.3}xy+{2:.3}y^2 = 1".format(
        x[0], x[1], x[2]
    )
)


def plot_ellipse_by_equation(ax, a, b, c):
    """a x^2 + b xy + c y^2 = 1 で表される楕円を描画する。
    """
    x = -np.linspace(-100, 100, 1000)
    y = np.linspace(-100, 100, 1000)
    X, Y = np.meshgrid(x, y)
    eqn = a * X ** 2 + b * X * Y + c * Y ** 2
    Z = 1

    # 楕円の式をそのまま描画できる関数は matplotlib にないので、
    # f = a * X ** 2 + b * X * Y + c * Y ** 2 の f = 1 の等高線を描画する形で
    # 楕円を描画する。
    ax.contour(X, Y, eqn, levels=[Z], colors=["r"])


# Plot the noisy data
fig, ax = plt.subplots()
ax.grid()
ax.scatter(X, Y, label="Data Points", color="g")
plot_ellipse_by_equation(ax, x[0], x[1], x[2])

結果

イメージ説明

緑の点がデータ、赤の線が近似した楕円

追記

長軸 (major axis) と x 軸の角度、短軸 (minor axis) と x 軸の角度が考えられますが、楕円の式を行列表記して、固有値問題を問き、出てきた固有ベクトル2つが長軸、短軸になります。(小さい固有値の固有ベクトルが長軸、大きい固有値の固有ベクトルが短軸に対応する)

二次形式の意味,微分,標準形など | 高校数学の美しい物語

あとは逆正弦関数で x 軸と長軸、または x 軸と短軸の角度を求めればよいです。

import matplotlib.pyplot as plt
import numpy as np

fig, ax = plt.subplots(figsize=(6, 6))
ax.grid()

a, b, c = 10, 12, 10


def plot_ellipse_by_equation(ax, a, b, c):
    """a x^2 + b xy + c y^2 = 1 で表される楕円を描画する。
    """
    x = np.linspace(-1, 1, 100)
    y = np.linspace(-1, 1, 100)
    X, Y = np.meshgrid(x, y)
    eqn = a * X ** 2 + b * X * Y + c * Y ** 2
    Z = 1

    # 楕円の式をそのまま描画できる関数は matplotlib にないので、
    # f = a * X ** 2 + b * X * Y + c * Y ** 2 の f = 1 の等高線を描画する形で
    # 楕円を描画する。
    ax.contour(X, Y, eqn, levels=[Z], colors=["r"])

# 固有値問題を解く
Q = np.array([[a, b / 2], [b / 2, c]])
eig_val, eig_vec = np.linalg.eig(Q)
print(eig_val, eig_vec)

# 長軸、短軸を求める。
idx1, idx2 = eig_val.argsort()
major_axis = eig_vec[:, idx1]  # 固有値が小さいほうの固有ベクトルが major axis
minor_axis = eig_vec[:, idx2]  # 固有値が大きいほうの固有ベクトルが minor axis

ax.annotate(
    s="", xy=major_axis, xytext=(0, 0), arrowprops=dict(facecolor="r"),
)
ax.annotate(
    s="", xy=minor_axis, xytext=(0, 0), arrowprops=dict(facecolor="b"),
)

plot_ellipse_by_equation(ax, a, b, c)

# 角度を求める
theta1 = np.rad2deg(np.arctan(major_axis[1] / major_axis[0]))
theta2 = np.rad2deg(np.arctan(minor_axis[1] / minor_axis[0]))
print(f"angle between major axis and x axis= {theta1:.2f}")
print(f"angle between minor axis and x axis= {theta2:.2f}")
# angle between major axis and x axis= -45.00
# angle between minor axis and x axis= 45.00

イメージ説明

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2020/01/22 16:55

    できました!
    ありがとうございます.

    あと,楕円の回転角の大きさを求める方法ありますか??
    わかるなら教えていただきたいです.

    キャンセル

  • 2020/01/22 21:12

    追記しました。

    キャンセル

  • 2020/01/22 22:15

    神様ですか?
    本当にありがとうございます.

    まだ,3次元楕円の問題があるのですが,少しわかった部分あるので自分で頑張りたいと思います.

    また,わからなかったら質問を投稿するので,お願いします.

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 89.09%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる