Подгонка плоскости к двумерному массиву

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

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

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

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


python-3.x arrays regression plane
person Jacques    schedule 26.01.2016    source источник
comment
Этот вопрос слишком широк. Вам нужна помощь в чтении топоса в определенном формате? Вам нужна помощь с математикой? Вам нужна помощь с повторным рендерингом? Вы знаете, как сделать все это на языке X, но вам нужно перевести это на Python? Пожалуйста, разбейте свой вопрос на небольшие части, на которые можно ответить.   -  person bishop    schedule 26.01.2016

Ответы (1)


arrow_upward
8
arrow_downward

Сначала вам нужно правильно представить свои данные.

У вас есть два аргумента 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()
person Anton    schedule 27.01.2016