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

平面方程拟合计算
平面方程的一般表达式为：
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 <math.h>
#include <stdlib.h>
#include <Windows.h>
#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<n;i++)
matrix2[i]=(double *)malloc(n*sizeof(double));
for(j=0;j<n;j++)
*(matrix2[j]+i)=(AlCo(matrix1,n,i,j)/d);
}
double Determinant(double* matrix[],int n)
{
double result=0,temp;
int i;
if(n==1)
result=(*matrix[0]);
else
{
for(i=0;i<n;i++)
{
temp=AlCo(matrix,n,n-1,i);
result+=(*(matrix[n-1]+i))*temp;
}
}
return result;
}
double AlCo(double* matrix[],int jie,int row,int column) {
double result;
if((row+column)%2 == 0)
result = Cofactor(matrix,jie,row,column);
else result=(-1)*Cofactor(matrix,jie,row,column);
return result;
}
double Cofactor(double* matrix[],int jie,int row,int column) {
double result;
int i,j;
double* smallmatr[MAX-1];
for(i=0;i<jie-1;i++)
smallmatr[i]= new double[jie - 1];
for(i=0;i<row;i++)
for(j=0;j<column;j++)
*(smallmatr[i]+j)=*(matrix[i]+j);
for(i=row;i<jie-1;i++)
for(j=0;j<column;j++)
*(smallmatr[i]+j)=*(matrix[i+1]+j);
for(j=column;j<jie-1;j++)
*(smallmatr[i]+j)=*(matrix[i]+j+1);
for(i=row;i<jie-1;i++)
for(j=column;j<jie-1;j++)
*(smallmatr[i]+j)=*(matrix[i+1]+j+1);
result = Determinant(smallmatr,jie-1);
for(i=0;i<jie-1;i++)
delete[] smallmatr[i];
return result;
}。