测绘-编写数字高程模型（DEM）内插程序

一、 实习目的

掌握移动曲面法数字高程模型内插原理及其内插子程序的设计方法，了解其它逐点高程内插方法的基本原理。

二、 实习内容

根据提供的10个数据点的坐标（X_n,Y_n,Z_n）和待求点的平面坐标（X_p,Y_p），要求利用移动二次曲面拟合法，由格网点P（X_p,Y_p）周围的10个已知点内插出待求格网点P的高程，编制相应的程序并进行调试，最后解算出格网点P的高程并提交源程序代码。

三、资料准备

已知数据点坐标

点号	X	Y	Z
1	102	110	15
2	109	113	18
3	105	115	19
4	103	103	17
5	108	105	21
6	105	108	15
7	115	104	20
8	118	108	15
9	116	113	17
10	113	118	22

 

编程计算点（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;
}