Применить матрицу вращения к координатам xy

Question

У меня есть xy координаты, которые представляют субъект в данном пространстве.На него ссылаются из другой точки и, следовательно, от центра .Как и в longitudinal axes не выравнивается вдоль x-axis.

Случайно сгенерированный ellipse ниже обеспечивает указание этого:

import numpy as np
from matplotlib.pyplot import scatter

xx = np.array([-0.51, 51.2])
yy = np.array([0.33, 51.6])
means = [xx.mean(), yy.mean()]  
stds = [xx.std() / 3, yy.std() / 3]
corr = 0.8         # correlation
covs = [[stds[0]**2          , stds[0]*stds[1]*corr], 
    [stds[0]*stds[1]*corr,           stds[1]**2]] 

m = np.random.multivariate_normal(means, covs, 1000).T
scatter(m[0], m[1])

К выпрямить координаты, которые я думал о применении вектора к rotation matrix.

Будет ли что-то подобное этой работе?

angle = 65.
theta = (angle/180.) * np.pi

rotMatrix = np.array([[np.cos(theta), -np.sin(theta)], 
                     [np.sin(theta),  np.cos(theta)]])

Это также может показаться глупым вопросом, но есть ли способ определить, перпендикулярны ли полученные vector из xy координат?Или вам просто нужно поиграться с rotation angle?

Answer 1

Вы можете использовать sklearn.decomposition.PCA (анализ главных компонентов) с n_components=2, чтобы извлечь наименьший угол, необходимый для поворота облака точек так, чтобы его главная ось была горизонтальной.

Runnableпример

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

np.random.seed(1)

xx = np.array([-0.51, 51.2])
yy = np.array([0.33, 51.6])
means = [xx.mean(), yy.mean()]  
stds = [xx.std() / 3, yy.std() / 3]
corr = 0.8         # correlation
covs = [[stds[0]**2,       stds[0]*stds[1]*corr], 
        [stds[0]*stds[1]*corr, stds[1]**2]]

m = np.random.multivariate_normal(means, covs, 1000)

pca = PCA(2)

# This was in my first answer attempt: fit_transform works fine, but it randomly 
# flips (mirrors) points across one of the principal axes.
# m2 = pca.fit_transform(m)

# Workaround: get the rotation angle from the PCA components and manually 
# build the rotation matrix.

# Fit the PCA object, but do not transform the data
pca.fit(m)

# pca.components_ : array, shape (n_components, n_features)
# cos theta
ct = pca.components_[0, 0]
# sin theta
st = pca.components_[0, 1]

# One possible value of theta that lies in [0, pi]
t = np.arccos(ct)

# If t is in quadrant 1, rotate CLOCKwise by t
if ct > 0 and st > 0:
    t *= -1
# If t is in Q2, rotate COUNTERclockwise by the complement of theta
elif ct < 0 and st > 0:
    t = np.pi - t
# If t is in Q3, rotate CLOCKwise by the complement of theta
elif ct < 0 and st < 0:
    t = -(np.pi - t)
# If t is in Q4, rotate COUNTERclockwise by theta, i.e., do nothing
elif ct > 0 and st < 0:
    pass

# Manually build the ccw rotation matrix
rotmat = np.array([[np.cos(t), -np.sin(t)], 
                   [np.sin(t),  np.cos(t)]])

# Apply rotation to each row of m
m2 = (rotmat @ m.T).T

# Center the rotated point cloud at (0, 0)
m2 -= m2.mean(axis=0)

fig, ax = plt.subplots()
plot_kws = {'alpha': '0.75',
            'edgecolor': 'white',
            'linewidths': 0.75}
ax.scatter(m[:, 0], m[:, 1], **plot_kws)
ax.scatter(m2[:, 0], m2[:, 1], **plot_kws)

Вывод

Предупреждение: pca.fit_transform() иногда переворачивает (отражает) облако точек

Основные компоненты могут случайным образом быть положительными или отрицательными.В некоторых случаях ваше облако точек может показаться перевернутым или даже зеркально отраженным по одной из его основных осей.(Чтобы проверить это, измените случайное начальное число и перезапустите код, пока не увидите переворачивание.) Здесь подробно обсуждается здесь (основано на R, но математика важна).Чтобы исправить это, вам нужно заменить строку fit_transform ручным перелистыванием знаков одного или обоих компонентов, а затем умножить матрицу перевернутых знаков на массив облаков точек.

Answer 2

Действительно, очень полезной концепцией здесь является линейное преобразование вектора v, выполняемое матрицей A.Если вы рассматриваете свои точки разброса как кончик векторов, начинающихся с (0,0), то очень легко повернуть их на любой угол theta.Матрица, которая выполняет такое вращение theta, будет иметь вид

A = [[cos(theta) -sin(theta]
     [sin(theta)  cos(theta)]]

Очевидно, что когда тета равна 90 градусам, это приводит к

A = [[ 0 1]
     [-1 0]]

И для применения вращения вам потребуется тольковыполнить умножение матриц w = A v

При этом текущая цель состоит в том, чтобы выполнить умножение матриц векторов, сохраненных в m, с x, y tips равными m[0],m[1].Вращенный вектор будет сохранен в m2.Ниже приведен соответствующий код для этого.Обратите внимание, что я транспонировал m для более простого вычисления умножения матриц (выполняется с @) и что угол поворота составляет 90 градусов против часовой стрелки.

import numpy as np
import matplotlib.pyplot as plt

xx = np.array([-0.51, 51.2])
yy = np.array([0.33, 51.6])
means = [xx.mean(), yy.mean()]  
stds = [xx.std() / 3, yy.std() / 3]
corr = 0.8         # correlation
covs = [[stds[0]**2          , stds[0]*stds[1]*corr], 
    [stds[0]*stds[1]*corr,           stds[1]**2]] 

m = np.random.multivariate_normal(means, covs, 1000).T
plt.scatter(m[0], m[1])

theta_deg = 90
theta_rad = np.deg2rad(theta_deg)
A = np.matrix([[np.cos(theta_rad), -np.sin(theta_rad)],
               [np.sin(theta_rad), np.cos(theta_rad)]])
m2 = np.zeros(m.T.shape)

for i,v in enumerate(m.T):
  w = A @ v.T
  m2[i] = w
m2 = m2.T

plt.scatter(m2[0], m2[1])

Это приводит к повернутому графику рассеяния: Вы можете быть уверены, что повернутая версия точно на 90 градусов против часовой стрелки с линейным преобразованием.

Редактировать

Чтобы найти вращениеугол, который необходимо применить для выравнивания графика рассеяния по оси x. Хороший подход - найти линейное приближение рассеянных данных с numpy.polyfit.Это дает линейную функцию, обеспечивая slope и точку пересечения оси y b.Затем получите угол поворота с помощью функции arctan наклона и вычислите матрицу преобразования, как и раньше.Вы можете сделать это, добавив следующую часть к коду

slope, b = np.polyfit(m[1], m[0], 1)
x = np.arange(min(m[0]), max(m[0]), 1)
y_line = slope*x + b
plt.plot(x, y_line, color='r')
theta_rad = -np.arctan(slope)

и результат к искомому графику

Редактировать 2

Поскольку @Peter Leimbigler указал, что numpy.polyfit не находит правильное глобальное направление рассеянных данных, я подумал, что вы можете получить средний наклон, усредняя значения x и yчасти данных.Это нужно для того, чтобы найти другой уклон, называемый slope2 (теперь обозначен зеленым цветом), чтобы применить вращение.Так просто,

slope, b = np.polyfit(m[1], m[0], 1)
x = np.arange(min(m[0]), max(m[0]), 1)
y_line = slope*x + b
slope2 = np.mean(m[1])/np.mean(m[0])
y_line2 = slope2*x + b
plt.plot(x, y_line, color='r')
plt.plot(x, y_line2, color='g')
theta_rad = -np.arctan(slope2)

И, применяя линейное преобразование с матрицей вращения, вы получаете

Answer 3

Если наклон двух линий, умноженных вместе, равен -1, то они перпендикулярны. В другом случае это верно, когда один наклон равен 0, а другой не определен (совершенно горизонтальная линия и совершенно вертикальная линия).