最近需要编一个利用程序，给一组三维点，拟合空间平面。找了很久，没有找到如何利用最小二乘拟合空间平面的纯数学理论知识，因为对这块不是很熟悉，有哪位大侠手里有这方面的数学资料，请发给我好么，我想知道很具体的理论算法，这样才能编程序。Thank you !!
    我的邮箱：yudandan.1985.hi@

[ 本帖最后由 丹85丹 于 2009-10-5 16:02 编辑 ]

百度上有，我是不会。给你个地址。http://blog.
有很详尽的理论推导过程，及现成的完整的程序。

谢谢你，这个我已经看到了。但是这个拟合的是一条空间的直线，我要拟合的是一个空间的平面，即a*x+b*y+c*z+m=0.呵呵。在百度上找了很久了，还是没有找到。

题目的意思是求一个平面P，使得 ∑(Distance(pi,P)^2)最小化, pi是给定点集里的点。 Right？
对于a*x+b*y+c*z+m=0,m!=0的情况,写成a*x+b*y+c*z+1=0
(x,y,z)记做p,
(a,b,c)记做n,
则上述平面可记做n·p+1=0
点pi到平面的有向距离di=n·pi+1

那么∑(Distance(pi,P)^2) = ∑(n·pi + 1)^2  := M(n)
其实就是确定适当的n使得M(n)最小化

联列
∂M/∂a = 0
∂M/∂b = 0
∂M/∂c = 0

推到后面似乎就是解一个3元一次线性方程组，用一下克拉默法则就直接得解了

补充，对于m=0的情况，你不妨也列出来去算，算出来一个m=1时候的最优(a,b,c)，一个m=0的最优(a,b,c)然后检查一下哪个产生的∑Distance最小就取哪个。

Thank you ！！！！
按你说的解下。太感谢了

觉得这个关于最小二乘拟合空间平面的也很好，只是我怎么上传呢，貌似不可以呢

楼主你好，麻烦问下你找到了最小二乘拟合平面的程序了吗？可不可以给我发一份，我急着用，谢谢了

#include "stdafx.h"
#include <math.h>
#include <stdio.h>

#define m_col 3
#define m_dim 3

void InvMatrix(const double M[][m_col], const int n, double invM[][m_col]);
void LinearFit(const double Coef[][m_dim], const int n, double *Para);

void main()
{
    double A[4][3] = {0,1,7,1,0,6,2,1,11,-1,-1,-1};
    double Para[3];
    LinearFit(A, 4, Para);
    printf("%lf %lf %lf",Para[0],Para[1],Para[2]);
}

void LinearFit(const double Coef[][m_dim], const int n, double *Para)
{
    double (*A)[m_col] = new double[m_dim][m_col];
    for(int i=0;i<m_dim-1;i++)
    {
        for(int j=0;j<m_dim-1;j++)
        {
            A[i][j] = 0;
            for(int k=0;k<n;k++)
                A[i][j]+= Coef[k][i]*Coef[k][j];
        }
    }
    for(int i=0;i<m_dim-1;i++)
    {
        A[i][m_dim-1] = A[m_dim-1][i] = 0;
        for(int j=0;j<n;j++)
            A[i][m_dim-1]+= Coef[j][i];
        A[m_dim-1][i] = A[i][m_dim-1];
    }
    A[m_dim-1][m_dim-1] = n;

    double *b = new double[m_dim];
    for(int i=0;i<m_dim-1;i++)
    {
        b[i] = 0;
        for(int j=0;j<n;j++)
            b[i]+= Coef[j][m_dim-1]*Coef[j][i];
    }
    b[m_dim-1] = 0;
    for(int i=0;i<n;i++)
        b[m_dim-1]+= Coef[i][m_dim-1];

    double (*invA)[m_col] = new double[m_dim][m_col];
    InvMatrix(A, m_dim, invA);

    for(int i=0;i<m_dim;i++)
    {
        Para[i] = 0;
        for(int j=0;j<m_dim;j++)
            Para[i]+= invA[i][j]*b[j];
    }

    delete []A;
    delete []invA;
    delete b;
}

void InvMatrix(const double M[][m_col], const int n, double invM[][m_col])
{
    if(n==1)
        invM[0][0] = 1/M[0][0];
    else
    {
        double cu, cd, normb;
        double (*b)[m_col] = new double[n][m_col];
        double (*invv)[m_col] = new double[n][m_col];
        for(int j=0;j<n;j++)
        {
            double *schu = new double[j];
            double *schd = new double[j];
            if(j>0)
            {
                for(int k=0;k<j;k++)
                {
                    schu[k] = schd[k] = 0;
                    for(int i=0;i<n;i++)
                    {
                        schu[k]+= b[i][k]*M[i][j];
                        schd[k]+= b[i][k]*b[i][k];
                    }
                }
            }

            normb = 0;
            for(int i=0;i<n;i++)
            {
                b[i][j] = M[i][j];
                if(j>0)
                {
                    for(int k=0;k<j;k++)
                        b[i][j]-= b[i][k]*schu[k]/schd[k];
                }
                normb+= b[i][j]*b[i][j];
            }
            normb = sqrt(normb);
            for(int i=0;i<n;i++)
                invv[j][i] = b[i][j]/normb;

            delete schu;
            delete schd;
        }

        double (*c)[m_col] = new double[n][m_col];
        for(int i=0;i<n;i++)
        {
            for(int j=0;j<=i;j++)
            {
                cu = cd = 0;
                if(j<i)
                {
                    for(int k=0;k<n;k++)
                    {
                        cu+= M[k][i]*b[k][j];
                        cd+= b[k][j]*b[k][j];
                    }
                    c[j][i] = cu/sqrt(cd);
                }
                else
                {
                    for(int k=0;k<n;k++)
                        cd+= b[k][j]*b[k][j];
                    c[j][i] = sqrt(cd);
                }
            }
        }

        double (*invc)[m_col] = new double[n][m_col];
        for(int j=0;j<n;j++)
        {
            for(int i=n-1;i>=0;i--)
            {
                if(i>j)
                    invc[i][j] = 0;
                else if(i==j)
                    invc[i][j] = 1/c[i][j];
                else
                {
                    invc[i][j] = 0;
                    for (int k=i+1;k<=j;k++)
                        invc[i][j]-= c[i][k]*invc[k][j];
                    invc[i][j]/= c[i][i];
                }
            }
        }

        for(int i=0;i<n;i++)
        {
            for(int j=0;j<n;j++)
            {
                invM[i][j] = 0;
                for(int k=0;k<n;k++)
                    invM[i][j]+= invc[i][k]*invv[k][j];
            }
        }
        delete []b;
        delete []c;
        delete []invv;
        delete []invc;
    }
}

最近做项目恰好遇到这个问题，编了个小程序。

linerfit函数为矩阵法的多元线性最小二乘求解函数， 当m_dim设置为3时即可求平面的最小二乘，当然m_dim设置为2可求平面直线最小二乘。注意的是m_col应定义为大于等于m_dim。若目标函数为z=ax+by+c;那么coef传参格式为:
x1 y1 z1
x2 y2 z2
....
求解结果para = [a b c].

当然你也可以用InvMatrix子函数做其他用途，其为QR法快速求逆矩阵。

[ 本帖最后由 eastlife0108 于 2014-4-7 15:21 编辑 ]