通过椭圆的线性回归显示了意外的行为

Question

我在一张空图像上画了一个二维椭圆。现在，我想通过椭圆拟合一条直线来获得长轴。知道有很多选择(PCA，图像时刻等)，我认为线性回归应该可以完成这项工作。然而，只有当椭圆的旋转平行于x轴时，它才“有效”。为什么会这样呢？难道任何均匀分布的对称点云都不应该给出中线吗？

这是我使用的代码：

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from skimage.draw import ellipse
from ipywidgets import interact
from sklearn.linear_model import LinearRegression

@interact
def rotateAndFit(rot:(-90,90)=-90):
    im = np.zeros((300,300), dtype=np.float64)
    im[ellipse(im.shape[0]//2, # center x
               im.shape[1]//2-10, # center y
               120,            # radius major axis
               40,             # radius minor axis
               im.shape,       # image shape
               rot/180*np.pi)] = 1  # rotation angle in degree,


    # Get corresponding x and y values
    y, x = np.where(im)

    # Do Linear Regression
    lr = LinearRegression()
    lr.fit(x[None].T,y)

    plt.imshow(im)
    plt.plot([0, 300], [lr.intercept_, lr.coef_[0]*300+lr.intercept_])
    plt.axis([0,300,300,0])
    plt.title('rotation $r = {}°$'.format(rot))

代码提供了以下输出：

​

​

我真的很困惑，有什么想法吗？我使用脊线和套索回归来调整权重，但它们降低了权重，但看起来权重，即斜率必须更陡峭，我认为线性回归低估了斜率。有趣的是，线性回归通常是“点”对称的，但不是跨线对称的……我理解，接近0°的行为，斜率不可能是无穷大的。但它至少应该在低旋转度下工作。

Answer 1

回归线不能与主轴重合，因为回归在y方向上最小化，而不是垂直于回归线。下面的示例使用正交距离回归，而不是y中的线性回归，它给出了所需的结果：

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from skimage.draw import ellipse
from ipywidgets import interact
from sklearn.linear_model import LinearRegression
from scipy.odr import ODR, Model, Data

def lin(beta, x):
    a,b = beta
    return a*x+b

@interact(rot=(-90,90))
def rotateAndFit(rot=-90):
    im = np.zeros((300,300), dtype=np.float64)
    im[ellipse(im.shape[0]//2, # center x
               im.shape[1]//2-10, # center y
               120,            # radius major axis
               40,             # radius minor axis
               im.shape,       # image shape
               rot/180*np.pi)] = 1  # rotation angle in pi (40°),


    y, x = np.where(im)

    d = Data(x,y)
    m = Model(lin)
    o = ODR(d, m, [0,0])
    out = o.run()
    a,b = out.beta

    plt.imshow(im)
    plt.plot([0, 300], [b, a*300+b])
    plt.axis([0,300,300,0])
    plt.title('rotation $r = {}°$'.format(rot))

Answer 2

为了补充上面的答案，这里有一个示意图，展示了你使用OLS和椭圆的点所做的事情。因为你在每个x处有两个点，所以你基本上是在每个x处找到y的平均值。对于未旋转的椭圆(图中的顶部椭圆)，这是你所期望的。但是当你旋转椭圆(底部椭圆)时，在给定的x处，椭圆上两点中间的点远离椭圆的主轴。所以，你必须最小化与旋转轴正交的方向。

​

​