3D Самолет наименьших квадратов

Question

Каков алгоритм для вычисления плоскости наименьших квадратов в (x, y, z) пространстве с учетом набора трехмерных точек данных? Другими словами, если бы у меня было несколько точек, таких как (1, 2, 3), (4, 5, 6), (7, 8, 9) и т. Д., Как можно было бы рассчитать плоскость наилучшего соответствия f (x, y) = ax + by + c? Каков алгоритм для получения a, b и c из набора трехмерных точек?

Answer 1

Если у вас есть n точек данных (x [i], y [i], z [i]), вычислите симметричную матрицу 3x3 A, записи которой:

sum_i x[i]*x[i],    sum_i x[i]*y[i],    sum_i x[i]
sum_i x[i]*y[i],    sum_i y[i]*y[i],    sum_i y[i]
sum_i x[i],         sum_i y[i],         n

Также вычислите вектор 3 элементов b:

{sum_i x[i]*z[i],   sum_i y[i]*z[i],    sum_i z[i]}

Тогда решите Ax = b для данных A и b. Три компонента вектора решения являются коэффициентами для плоскости подгонки наименьших квадратов {a, b, c}.

Обратите внимание, что это подгонка "обычных наименьших квадратов", которая подходит только тогда, когда ожидается, что z будет линейной функцией от x и y. Если вы ищете в общем случае «плоскость наилучшего соответствия» в 3-мерном пространстве, вы можете узнать о «геометрических» наименьших квадратах.

Обратите внимание, что это не удастся, если ваши точки расположены на одной линии, как ваши примеры точек.

Answer 2

, если кто-то не скажет мне, как набирать уравнения здесь, позвольте мне записать окончательные вычисления, которые вы должны сделать:

сначала, с учетом точек r_i \ n \ R, i = 1..N, вычислить центр масс всех точек:

r_G = \frac{\sum_{i=1}^N r_i}{N}

затем вычислите вектор нормали n, который вместе с базовым вектором r_G определяет плоскость, вычисляя матрицу 3x3 A как

A = \sum_{i=1}^N (r_i - r_G)(r_i - r_G)^T

с этой матрицей вектор нормали n теперь задается собственным вектором A, соответствующим минимальному собственному значению A.

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

Это решение основано на теореме Релея-Ритца для эрмитовой матрицы A.

Answer 3

См. «Подгонка данных наименьших квадратов» Дэвида Эберли, чтобы узнать, как я придумал этот способ, чтобы минимизировать геометрическое соответствие (ортогональное расстояние от точек до плоскости).

bool Geom_utils::Fit_plane_direct(const arma::mat& pts_in, Plane& plane_out)
{
    bool success(false);
    int K(pts_in.n_cols);
    if(pts_in.n_rows == 3 && K > 2)  // check for bad sizing and indeterminate case
    {
        plane_out._p_3 = (1.0/static_cast<double>(K))*arma::sum(pts_in,1);
        arma::mat A(pts_in);
        A.each_col() -= plane_out._p_3; //[x1-p, x2-p, ..., xk-p]
        arma::mat33 M(A*A.t());
        arma::vec3 D;
        arma::mat33 V;
        if(arma::eig_sym(D,V,M))
        {
            // diagonalization succeeded
            plane_out._n_3 = V.col(0); // in ascending order by default
            if(plane_out._n_3(2) < 0)
            {
                plane_out._n_3 = -plane_out._n_3; // upward pointing
            }
            success = true;
        }
    }
    return success;
}

По времени 37 микросекунды подгонка плоскости к 1000 точкам (Windows 7, i7, 32-битная программа)

Answer 4

Как и в случае любого метода наименьших квадратов, вы действуете так:

Перед началом кодирования

1. Запишите уравнение для плоскости при некоторой параметризации, скажем 0 = ax + by + z + d в параметрах (a, b, d).

2. Найдите выражение D(\vec{v};a, b, d) для расстояния от произвольной точки \vec{v}.

3. Запишите сумму S = \sigma_i=0,n D^2(\vec{x}_i) и упростите ее, пока она не будет выражена в виде простых сумм компонентов v, таких как \sigma v_x, \sigma v_y^2, \sigma v_x*v_z ...

4. Запишите выражения минимизации для каждого параметра dS/dx_0 = 0, dS/dy_0 = 0 ..., которые дают вам набор из трех уравнений в трех параметрах и суммы из предыдущего шага.

5. Решите этот набор уравнений для параметров.

(или для простых случаев просто посмотрите на форму). Использование пакета символической алгебры (например, Mathematica) может значительно облегчить вам жизнь.

Кодировка

• Введите код для формирования необходимых сумм и найдите параметры из последнего набора, указанного выше.

Альтернативы

Обратите внимание: если бы у вас было всего три точки, вам лучше было бы найти самолет, который проходит через них.

Также, если аналитическое решение неосуществимо (не в случае с плоскостью, но возможно вообще), вы можете выполнить шаги 1 и 2 и использовать минимизатор Монте-Карло на сумму в шаге 3 .

Answer 5

Уравнение для плоскости: ax + by + c = z. Поэтому настройте такие матрицы со всеми вашими данными:

    x_0   y_0   1  
A = x_1   y_1   1  
          ... 
    x_n   y_n   1  

И

    a  
x = b  
    c

И

    z_0   
B = z_1   
    ...   
    z_n

Другими словами: Ax = B. Теперь решите для x, какие у вас коэффициенты. Но поскольку (я полагаю) у вас более 3 баллов, система переопределена, поэтому вам нужно использовать левое псевдообращение. Итак, ответ:

a 
b = (A^T A)^-1 A^T B
c

А вот простой код Python с примером:

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

N_POINTS = 10
TARGET_X_SLOPE = 2
TARGET_y_SLOPE = 3
TARGET_OFFSET  = 5
EXTENTS = 5
NOISE = 5

# create random data
xs = [np.random.uniform(2*EXTENTS)-EXTENTS for i in range(N_POINTS)]
ys = [np.random.uniform(2*EXTENTS)-EXTENTS for i in range(N_POINTS)]
zs = []
for i in range(N_POINTS):
    zs.append(xs[i]*TARGET_X_SLOPE + \
              ys[i]*TARGET_y_SLOPE + \
              TARGET_OFFSET + np.random.normal(scale=NOISE))

# plot raw data
plt.figure()
ax = plt.subplot(111, projection='3d')
ax.scatter(xs, ys, zs, color='b')

# do fit
tmp_A = []
tmp_b = []
for i in range(len(xs)):
    tmp_A.append([xs[i], ys[i], 1])
    tmp_b.append(zs[i])
b = np.matrix(tmp_b).T
A = np.matrix(tmp_A)
fit = (A.T * A).I * A.T * b
errors = b - A * fit
residual = np.linalg.norm(errors)

print "solution:"
print "%f x + %f y + %f = z" % (fit[0], fit[1], fit[2])
print "errors:"
print errors
print "residual:"
print residual

# plot plane
xlim = ax.get_xlim()
ylim = ax.get_ylim()
X,Y = np.meshgrid(np.arange(xlim[0], xlim[1]),
                  np.arange(ylim[0], ylim[1]))
Z = np.zeros(X.shape)
for r in range(X.shape[0]):
    for c in range(X.shape[1]):
        Z[r,c] = fit[0] * X[r,c] + fit[1] * Y[r,c] + fit[2]
ax.plot_wireframe(X,Y,Z, color='k')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()

Answer 6

CGAL::linear_least_squares_fitting_3

Функция linear_least_squares_fitting_3 вычисляет наиболее подходящее 3D линия или плоскость (в смысле наименьших квадратов) набора трехмерных объектов, таких как как точки, сегменты, треугольники, сферы, шары, кубоиды или тетраэдры.

http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Principal_component_analysis_ref/Function_linear_least_squares_fitting_3.html

Answer 7

Это сводится к проблеме Total Least Squares , которая может быть решена с помощью SVD разложения.

C ++ код с использованием OpenCV:

float fitPlaneToSetOfPoints(const std::vector<cv::Point3f> &pts, cv::Point3f &p0, cv::Vec3f &nml) {
    const int SCALAR_TYPE = CV_32F;
    typedef float ScalarType;

    // Calculate centroid
    p0 = cv::Point3f(0,0,0);
    for (int i = 0; i < pts.size(); ++i)
        p0 = p0 + conv<cv::Vec3f>(pts[i]);
    p0 *= 1.0/pts.size();

    // Compose data matrix subtracting the centroid from each point
    cv::Mat Q(pts.size(), 3, SCALAR_TYPE);
    for (int i = 0; i < pts.size(); ++i) {
        Q.at<ScalarType>(i,0) = pts[i].x - p0.x;
        Q.at<ScalarType>(i,1) = pts[i].y - p0.y;
        Q.at<ScalarType>(i,2) = pts[i].z - p0.z;
    }

    // Compute SVD decomposition and the Total Least Squares solution, which is the eigenvector corresponding to the least eigenvalue
    cv::SVD svd(Q, cv::SVD::MODIFY_A|cv::SVD::FULL_UV);
    nml = svd.vt.row(2);

    // Calculate the actual RMS error
    float err = 0;
    for (int i = 0; i < pts.size(); ++i)
        err += powf(nml.dot(pts[i] - p0), 2);
    err = sqrtf(err / pts.size());

    return err;
}

Answer 8

Похоже, все, что вы хотите сделать, это линейная регрессия с двумя регрессорами. Страница википедии на эту тему должна рассказать вам все, что вам нужно знать, а затем немного.

Answer 9

Все, что вам нужно сделать, это решить систему уравнений.

Если это ваши очки: (1, 2, 3), (4, 5, 6), (7, 8, 9)

Это дает вам уравнения:

3=a*1 + b*2 + c
6=a*4 + b*5 + c
9=a*7 + b*8 + c

Итак, ваш вопрос на самом деле должен звучать так: как мне решить систему уравнений?

Поэтому я рекомендую прочитать этот ТАК вопрос.

Если я неправильно понял ваш вопрос, сообщите нам.

EDIT

Не обращайте внимания на мой ответ, поскольку вы, вероятно, имели в виду что-то другое.