一、 实习目的
掌握移动曲面法数字高程模型内插原理及其内插子程序的设计方法，了解其它逐点高程内插方法的基本原理。
二、 实习内容
根据提供的10个数据点的坐标（Xn,Yn,Zn）和待求点的平面坐标（Xp,Yp），要求利用移动二次曲面拟合法，由格网点P（Xp,Yp）周围的10个已知点内插出待求格网点P的高程，编制相应的程序并进行调试，最后解算出格网点P的高程并提交源程序代码。
三、资料准备
已知数据点坐标
编程计算点（110,110）上的高程。
四. 操作步骤
1、读入已知点的坐标，建立以待定点为原点的局部坐标系；
2、建立误差方程式：
3、组成法方程，解算六个系数：
4、整理计算结果，以实习报告的形式递交成果。
程序：
// 摄影测量实习2.cpp : 定义控制台应用程序的入口点。
//
#include "stdafx.h"#include "iostream"using namespace std;void InverseMatrix(double a[],int n){int *is,*js,i,j,k,l,u,v;double d,p;is = (int*)malloc(n * sizeof(int));js = (int*)malloc(n * sizeof(int));for(k = 0; k <= n - 1;k++){d = 0.0;for(i = k; i <= n - 1; i++)for(j = k; j <= n - 1; j++){ l = i * n + j; p = fabs(a[l]);if(p > d){ d = p; is[k] = i; js[k] = j;}}if(d + 1.0 == 1.0){ free(is); free(js); printf("err**not inv\n");return ;}if(is[k] != k)for(j = 0; j <= n - 1;j++){ u = k * n + j; v = is[k] * n + j;p = a[u]; a[u] = a[v]; a[v] = p;}if(js[k] != k)for(i = 0; i <= n - 1; i++){ u = i * n + k; v = i * n + js[k];p = a[u]; a[u] = a[v]; a[v] = p;}l = k * n + k;a[l] = 1.0/a[l];for(j = 0; j <= n - 1;j++)if(j != k){ u = k * n + j; a[u] = a[u] * a[l];}for(i = 0; i <= n - 1;i++)if(i != k)for(j = 0; j <= n - 1;j++)if(j != k){ u = i * n + j;a[u] = a[u] - a[i * n + k] * a[k * n + j];}for(i = 0; i <= n - 1;i++)if(i != k){u = i * n + k; a[u] = -a[u] * a[l];}}for(k = n - 1;k >= 0;k--){ if (js[k] != k)for (j = 0;j <= n - 1;j++){ u = k * n + j; v = js[k] * n + j;p = a[u]; a[u] = a[v]; a[v] = p;}if(is[k] != k)for(i = 0; i <= n - 1; i++){ u = i * n + k; v = i * n + is[k];p = a[u]; a[u] = a[v]; a[v] = p;}}free(is); free(js);}typedef struct {double X;double Y;double Z;}DATA;int _tmain(int argc, _TCHAR* argv[]){FILE *fp1; DATA data[10]; fp1=fopen("Table.txt","r"); if(fp1==0) { printf("Error!Can't open it!\n"); return 0;} for(int i=0;!feof(fp1);i++){ fscanf(fp1,"%lf %lf %lf",&data[i].X,&data[i].Y,&data[i].Z); data[i].X=data[i].X-110;data[i].Y=data[i].Y-110;} fclose(fp1);double L1[10];double result[6]={0};for(int i=0;i<10;i++){L1[i]=data[i].Z;}double A[10][6];for(int i=0;i<10;i++){A[i][0]=data[i].X*data[i].X;A[i][1]=data[i].X*data[i].Y;A[i][2]=data[i].Y*data[i].Y;A[i][3]=data[i].X;A[i][4]=data[i].Y;A[i][5]=1;}double AT[6][10];for(int i=0;i<10;i++){for(int j=0;j<6;j++){AT[j][i]=A[i][j];}}double t;double ATA[6][6];for(int i=0;i<6;i++){for (int m=0;m<6;m++){t=0;for (int j=0;j<10;j++){t=AT[i][j]*A[j][m]+t; }ATA[i][m]=t;}} InverseMatrix(*ATA,6);double A1[6][10];double tt;for(int i=0;i<6;i++){for (int m=0;m<10;m++){tt=0;for (int j=0;j<6;j++){tt=ATA[i][j]*AT[j][m]+tt; }A1[i][m]=tt;}}for(int i=0;i<6;i++){ double o=0;for( int j=0;j<10;j++){o=A1[i][j]*L1[j]+o;}result[i]=o;}printf("%f %f %f %f %f %f ",result[0],result[1],result[2],result[3],result[4],result[5]);system("pause");return 0;}
