多点最小二乘法平面方程拟合计算

平面方程拟合计算

平面方程的一般表达式为：

0=+++D Cz By Ax , (0≠C )

C

D y C B x C A z ---

= 记：C D a C B a C A a -=-=-=210,, 则：210a y a x a z ++=

平面方程拟合:

对于一系列的n 个点)3(≥n ：

1,,1,0),,,(-=n i z y x i i i

要用点1,,1,0),,,(-=n i z y x i i i 拟合计算上述平面方程，则使： ()∑-=-++=1

02

210n i z a y a x a S 最小。

要使得S 最小，应满足： 2,1,0,0==∂∂k a S k

即：⎪⎩⎪⎨⎧=-++=-++=-++∑

∑∑0)(20)(20)(2210210210i i i i i i i i i i i z a y a x a y z a y a x a x z a y a x a

有，⎪⎩⎪⎨⎧=++=++=++∑∑∑∑∑∑∑∑∑∑∑i

i i i i i i i i i i i i i i z n a y a x a z y y a y a y x a z x x a y x a x a 21022102120 或，⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛∑∑∑∑∑∑∑∑∑∑∑i i i i i i i i i

i i i i i i z z y z x a a a n y x y

y y x x

y x x 21022 解上述线形方程组，得：210,,a a a

即：210a y a x a z ++=

下面程序实际求得的是以下的参数：

01=+++Z D

C Y

D B X D A 即:AX+BY+CZ+1=0

其程序代码如下：

#include "stdafx.h"

#include

#include

#include

#define MAX 10

void Inverse(double *matrix1[],double *matrix2[],int n,double d); double Determinant(double* matrix[],int n);

double AlCo(double* matrix[],int jie,int row,int column); double Cofactor(double* matrix[],int jie,int row,int column);

int _tmain(int argc, _TCHAR* argv[])

{

double array[12][3],Y[3];

double A,B,C;

A =

B =

C = 0.0;

ZeroMemory(array,sizeof(array));

ZeroMemory(Y ,sizeof(Y));

for (int i = 0;i < 12;i++)

{

for (int j = 0;j < 3;j++)

{

array[i][j] = (double)rand();

}

}

for (int i = 0; i < 12;i++)

{

array[i][0] = 1.0;

}//设计了12个最简单的数据点，x = 1平面上的点， double *Matrix[3],*IMatrix[3];

for (int i = 0;i < 3;i++)

{

Matrix[i] = new double[3];

IMatrix[i] = new double[3];

}

for (int i = 0;i < 3;i++)

{

for (int j = 0;j < 3;j++)

{

*(Matrix[i] + j) = 0.0;

}

}

for (int j = 0;j < 3;j++)

{

for (int i = 0;i < 12;i++)

{

*(Matrix[0] + j) += array[i][0]*array[i][j];

*(Matrix[1] + j) += array[i][1]*array[i][j];

*(Matrix[2] + j) += array[i][2]*array[i][j];

Y[j] -= array[i][j];

}

}

double d = Determinant(Matrix,3);

if (abs(d) < 0.0001)

{

printf("\n矩阵奇异");

getchar();

return -1;

}

Inverse(Matrix,IMatrix,3,d);

for (int i = 0;i < 3;i++)

{

A += *(IMatrix[0] + i)*Y[i];

B += *(IMatrix[1] + i)*Y[i];

C += *(IMatrix[2] + i)*Y[i];

}

printf("\n A = %5.3f,B = %5.3f,C= %5.3f",A,B,C);

for (int i = 0;i < 3;i++)

{

delete[] Matrix[i];

delete[] IMatrix[i];

}

getchar();

return 0;

}

void Inverse(double *matrix1[],double *matrix2[],int n,double d) {

int i,j;

for(i=0;i

matrix2[i]=(double *)malloc(n*sizeof(double));