将平面拟合到二维数组答案

【问题标题】：Fitting a plane to a 2D array将平面拟合到二维数组
【发布时间】：2016-01-26 01:04:22
【问题描述】：

我有一个拓扑图像，我正在尝试使用 Python 执行平面减法。该图像是一个 256x256 二维数组，由介于 0 和 1 之间的 float32 值组成。

我想做的是使用线性回归来拟合这个数据的平面，然后从原始值中减去这个平面。

我不确定如何实现这一目标。

我是 Python 语言的新手，感谢任何帮助。

【问题讨论】：

这个问题太笼统了。您需要帮助阅读特定格式的拓扑吗？你需要数学方面的帮助吗？您需要重新渲染方面的帮助吗？你知道如何在语言 X 中完成所有这些工作，但需要将其翻译成 Python 吗？请将您的问题分解成可回答的小块。

标签： arrays python-3.x regression plane

【解决方案1】：

首先，您需要以正确的方式表示您的数据。

您有两个参数X1 和X2，它们定义了拓扑图像的坐标，还有一个目标值Y，它定义了每个点的高度。对于回归分析，您需要通过添加 X0 来扩展参数列表，它始终等于 1。

然后你需要将参数和目标分别展开成矩阵[m*m x 3]和[m*m x 1]。您想找到矢量theta，它将描述所需的平面。为此，您可以使用Normal Equation：

为了演示该方法，我生成了一些拓扑表面。可以在图片上看到曲面、拟合平面的曲面和减法后的曲面：

代码如下：

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np

m = 256 #size of the matrix

X1, X2 = np.mgrid[:m, :m]

fig = plt.figure()
ax = fig.add_subplot(3,1,1, projection='3d')
jet = plt.get_cmap('jet')

#generation of the surface
F = 3        
i = np.minimum(X1, m-X1-1)
j = np.minimum(X2, m-X2-1)
H = np.exp(-.5*(np.power(i, 2)  +  np.power(j, 2)   )/(F*F))
Y = np.real(  np.fft.ifft2   (H  *  np.fft.fft2(  np.random.randn(m, m))))
a = 0.0005; b = 0.0002; #parameters of the tilted plane
Y = Y + (a*X1 + b*X2); #adding the plane
Y = (Y - np.min(Y)) / (np.max(Y) - np.min(Y)) #data scaling

#plot the initial topological surface
ax.plot_surface(X1,X2,Y, rstride = 1, cstride = 1, cmap = jet, linewidth = 0)


#Regression
X = np.hstack(   ( np.reshape(X1, (m*m, 1)) , np.reshape(X2, (m*m, 1)) ) )
X = np.hstack(   ( np.ones((m*m, 1)) , X ))
YY = np.reshape(Y, (m*m, 1))

theta = np.dot(np.dot( np.linalg.pinv(np.dot(X.transpose(), X)), X.transpose()), YY)

plane = np.reshape(np.dot(X, theta), (m, m));

ax = fig.add_subplot(3,1,2, projection='3d')
ax.plot_surface(X1,X2,plane)
ax.plot_surface(X1,X2,Y, rstride = 1, cstride = 1, cmap = jet, linewidth = 0)


#Subtraction
Y_sub = Y - plane
ax = fig.add_subplot(3,1,3, projection='3d')
ax.plot_surface(X1,X2,Y_sub, rstride = 1, cstride = 1, cmap = jet, linewidth = 0)

plt.show()

【讨论】：