Установка плоскости в 2D-массив

Question

У меня есть топологическое изображение, которое я пытаюсь выполнить вычитанием плоскости при использовании Python. Изображение представляет собой 256x256 двухмерный массив значений float32 между 0 и 1.

Что я хочу сделать, это использовать линейную регрессию, чтобы соответствовать плоскости этим данным и впоследствии вычесть эту плоскость из исходных значений.

Я не уверен, как это сделать.

Я новичок в языке Python и ценю любую помощь.

Answer 1

Сначала вам необходимо представить свои данные соответствующим образом.

У вас есть два аргумента X1 и X2, которые определяют координаты вашего топологического изображения и одно целевое значение Y, которое определяет высоту каждой точки. Для регрессионного анализа вам необходимо расширить список аргументов, добавив X0, который всегда равен единице.

Затем вам нужно развернуть параметры и цель в матрицах [m*m x 3] и [m*m x 1] соответственно. Вы хотите найти вектор theta, который будет описывать нужный вам самолет. Для этой цели можно использовать Нормальное уравнение:

Для демонстрации подхода я произвел некоторую топологическую поверхность. Вы можете увидеть поверхность, поверхность с подогнанной плоскости и поверхности после вычитания на картинке:

Вот код:

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()