3D-плоскость с наименьшими квадратами

Ответ 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-пространстве, вы можете узнать о "геометрических" наименьших квадратах.

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

Ответ 2

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

сначала, заданные точки r_i\n\R, я = 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.

Ответ 3

Как и в любом методе наименьших квадратов, вы выполняете следующее:

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

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

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

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

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

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

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

Кодирование

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

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

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

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

Ответ 4

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

Ответ 5

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

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-разрядная программа)

Ответ 6

Уравнение для плоскости: 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()

Ответ 7

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

Если это ваши очки: (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

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

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

Если я неправильно понял ваш вопрос, дайте нам знать.

ИЗМЕНИТЬ

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

Ответ 8

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

Ответ 9

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

Код С++ с использованием 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;
}