سؤال

ما هي الخوارزمية المستخدمة لحساب مستوى المربعات الصغرى في الفضاء (x، y، z)، في ضوء مجموعة من نقاط البيانات ثلاثية الأبعاد؟بمعنى آخر، إذا كان لدي مجموعة من النقاط مثل (1، 2، 3)، (4، 5، 6)، (7، 8، 9)، وما إلى ذلك، فكيف يمكن حساب أفضل مستوى مناسب لـ f؟ (س، ص) = الفأس + بواسطة + ج؟ما هي الخوارزمية للحصول على a وb وc من مجموعة من النقاط ثلاثية الأبعاد؟

هل كانت مفيدة؟

المحلول

إذا كان لديك 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

قم أيضًا بحساب ناقل العناصر الثلاثة 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، فقد ترغب في التعرف على المربعات الصغرى "الهندسية".

لاحظ أيضًا أن هذا سيفشل إذا كانت نقاطك في سطر، كما هو الحال مع نقاط المثال الخاصة بك.

نصائح أخرى

ما لم يخبرني شخص ما كيفية كتابة المعادلات هنا، اسمحوا لي أن أكتب الحسابات النهائية ما عليك القيام به:

أولا، نظرا نقاط r_i \ ن \ R، ط = 1..N، وحساب مركز الكتلة من جميع نقاط:

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

وبعد ذلك، حساب ناقلات ن العادي، الذي جنبا إلى جنب مع r_G ناقلات قاعدة تحدد الطائرة عن طريق حساب 3X3 المصفوفة A ك

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

ومع هذه المصفوفة، والآن نظرا لناقلات ن العادي من بالمتجه الذاتي من A المقابلة لالقيمة الذاتية الحد الأدنى من A.

لمعرفة المزيد عن أزواج بالمتجه الذاتي / القيمة الذاتية، واستخدام أي مكتبة الجبر الخطي الذي تختاره.

ويستند هذا الحل على Rayleight ريتز نظرية لالهرمتري مصفوفة A.

وانظر "المربعات تركيب البيانات عن طريق ديفيد Eberly لكيفية خطرت لي هذه واحدة للحد من صالح الهندسي (المسافة المتعامدة من نقطة إلى الطائرة).

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 نقطة (ويندوز 7، I7، برنامج 32BIT و)

كما هو الحال مع أي نهج المربعات الصغرى، يمكنك المتابعة على النحو التالي:

قبل البدء في الترميز

  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.

CGAL::linear_least_squares_fitting_3
<اقتباس فقرة>   

وظيفة linear_least_squares_fitting_3 يحسب أفضل 3D المناسب   خط أو الطائرة (في أقل شعور الساحات) من مجموعة من الكائنات 3D مثل   كنقاط، قطاعات، مثلثات، دوائر، والكرات، أشباه المكعبات أو رباعي الأسطح.

http://www.cgal.org/Manual/ احدث / doc_html / cgal_manual / Principal_component_analysis_ref / Function_linear_least_squares_fitting_3.html

ومعادلة طائرة هي: الفأس + من + ج = ض. لذلك انشاء المصفوفات مثل هذا مع جميع البيانات الخاصة بك:

    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

وبعبارة أخرى: فأس = B. الآن حل ل x والتي هي معاملات الخاصة بك. ولكن منذ (افترض) كان لديك أكثر من 3 نقاط، وهذا النظام هو أكثر من التي تحددها لذلك تحتاج إلى استخدام معكوس الزائفة الأيسر. هكذا كان الجواب:

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

وهنا بعض التعليمات البرمجية بايثون بسيطة مع مثال:

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

وهذا يقلل إلى إجمالي المربعات المشكلة، يمكن حلها باستخدام <أ href ل = "https://en.wikipedia.org/wiki/Singular-value_decomposition" يختلط = "نوفولو noreferrer"> SVD التحلل.

وC ++ التعليمات البرمجية باستخدام مكتبة برمجية مفتوحة للرؤية الحاسوبية:

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;
}

وكل ما عليك القيام به هو حل نظام المعادلات.

وإذا كانت هذه هي النقاط: (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 السؤال.

إذا كنت أسيء فهمها سؤالك إعلامنا.

تعديل : ل

وتجاهل جوابي وربما كنت تعني شيء آخر.

وهذا يبدو وكأنه كل ما تريد القيام به هو الانحدار الخطي مع 2 regressors. و يكيبيديا الصفحة على هذا الموضوع ينبغي أن أقول لك كل ما تحتاج إلى معرفته وبعد ذلك بعض.

مرخصة بموجب: CC-BY-SA مع الإسناد
لا تنتمي إلى StackOverflow
scroll top