回答編集履歴

2020/01/22 12:11

投稿

スコア21956

test CHANGED Viewed

@@ -52,7 +52,9 @@
 ```python
+import cv2
-import numpy
+import numpy as np
 import matplotlib.pyplot as plt
@@ -60,385 +62,131 @@
-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.0
-        c = self.coef[2, 0]
-        d = self.coef[3, 0] / 2.0
-        f = self.coef[4, 0] / 2.0
-        g = self.coef[5, 0]
-        # finding center of ellipse [eqn.19 and 20] from (**)
-        x0 = (c * d - b * f) / (b ** 2.0 - a * c)
-        y0 = (a * f - b * d) / (b ** 2.0 - 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 = 0.5 * numpy.arctan((2.0 * 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=0.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]]
+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 で表される楕円を描画する。
     """
-    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.0
-    )
-    ellipse_y = (
-        center[1]
-        + width * numpy.cos(t) * numpy.sin(phi)
-        + height * numpy.sin(t) * numpy.cos(phi)
-        + y_noise / 2.0
-    )
-    return [ellipse_x, ellipse_y]
-import cv2
+    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])
+```
+## 結果
+![イメージ説明](6330b456050f23b7c2b6c56d9d8c4824.png)
+緑の点がデータ、赤の線が近似した楕円
+## 追記
+長軸 (major axis) と x 軸の角度、短軸 (minor axis) と x 軸の角度が考えられますが、楕円の式を行列表記して、固有値問題を問き、出てきた固有ベクトル2つが長軸、短軸になります。(小さい固有値の固有ベクトルが長軸、大きい固有値の固有ベクトルが短軸に対応する)
+[二次形式の意味，微分，標準形など | 高校数学の美しい物語](https://mathtrain.jp/quadraticform)
+あとは逆正弦関数で x 軸と長軸、または x 軸と短軸の角度を求めればよいです。
+```python
+import matplotlib.pyplot as plt
 import numpy as np
-import matplotlib.pyplot as plt
-from matplotlib.patches import Ellipse
-alpha = 5
-beta = 3
-N = 500
-DIM = 2
-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])
+fig, ax = plt.subplots(figsize=(6, 6))
-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 = 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(
+ax.grid()
-    "The ellipse is given by {0:.3}x^2 + {1:.3}xy+{2:.3}y^2 = 1".format(
-        x[0], x[1], x[2]
+a, b, c = 10, 12, 10
-    )
-)
@@ -450,9 +198,9 @@
     """
-    x = -np.linspace(-100, 100, 1000)
+    x = np.linspace(-1, 1, 100)
-    y = np.linspace(-100, 100, 1000)
+    y = np.linspace(-1, 1, 100)
     X, Y = np.meshgrid(x, y)
@@ -472,28 +220,60 @@
+# 固有値問題を解く
+Q = np.array([[a, b / 2], [b / 2, c]])
+eig_val, eig_vec = np.linalg.eig(Q)
-# Plot the noisy data
+print(eig_val, eig_vec)
+# 長軸、短軸を求める。
-fig, ax = plt.subplots()
+idx1, idx2 = eig_val.argsort()
+major_axis = eig_vec[:, idx1]  # 固有値が小さいほうの固有ベクトルが major axis
+minor_axis = eig_vec[:, idx2]  # 固有値が大きいほうの固有ベクトルが minor axis
-ax.grid()
+ax.annotate(
-ax.scatter(X, Y, label="Data Points", color="g")
+    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, x[0], x[1], x[2])
+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
 ```
-## 結果
-![イメージ説明](6330b456050f23b7c2b6c56d9d8c4824.png)
+![イメージ説明](b145bd531d08d205d36f4fde20ad4cfc.png)
-緑の点がデータ、赤の線が近似した楕円

2020/01/22 12:11

投稿

tiitoi

スコア21956

test CHANGED Viewed

@@ -485,3 +485,15 @@
 plot_ellipse_by_equation(ax, x[0], x[1], x[2])
 ```
+## 結果
+![イメージ説明](6330b456050f23b7c2b6c56d9d8c4824.png)
+緑の点がデータ、赤の線が近似した楕円

2020/01/22 06:27

投稿

tiitoi

スコア21956

test CHANGED Viewed

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