现阶段生成数字高程模型（DEM）的方法较多，如以摄影测量得到的像对为数据源跟踪生成等高线及DEM，由机载激光测距仪记录规则点集后生产数据，也可采用传统的地形图扫描后跟踪等高线，记录一连串离散点集，接着运用各类算法进行处理，最后生成不规则三角网（TIN）与规则格网（GRID）DEM的方法。本文主要介绍的就是以等高线（参考图一）和离散点集为数据源，产生TIN与GRID DEM的技术路线。具体步骤如下：
1）跟踪等高线生成离散点集，记录在文本文件中。参考图二和图三。
2）读取文本文件中的数据，进行预处理。主要工作是找到XY轴方向上最小最大数值，压缩数据范围，避免数据范围跨度太大或太小，即出现数据分布稠密或稀疏的情况。
while (!_demfile.eof()) {
_demfile >> point3dXYZ[i][0]>> point3dXYZ[i][1]>> point3dXYZ[i][2];
point3dXYZ[i][2] = point3dXYZ[i][2] / 2; //因为XY轴在随后调整，因此相应调Z轴数值
if(xMin>point3dXYZ[i][0]) xMin = point3dXYZ[i][0];//得到整个范围的最大与最小数值
if(xMax<point3dXYZ[i][0]) xMax = point3dXYZ[i][0];
if(yMin>point3dXYZ[i][1]) yMin = point3dXYZ[i][1];
if(yMax<point3dXYZ[i][1]) yMax = point3dXYZ[i][1];
if(zMin>point3dXYZ[i][2]) zMin = point3dXYZ[i][2];
if(zMax<point3dXYZ[i][2]) zMax = point3dXYZ[i][2];
i++;
}
_demfile.close(); //文件流读取完毕，关闭
_demfile.clear(); //文件流清除
3）完成TIN中点，线，三角形和网的定义。
class AFX_EXT_CLASS TIN_Point
{ //TIN中的点
friend class TIN;
public:
int Get_ID (void) { return( m_ID );} //ID数值
const POINT3d & Get_Point (void){ return( m_Point ); } //普通的点
double Get_X(void) {return( m_Point[0] ); } //点上的X数值
double Get_Y(void) {return( m_Point[1] ); } //点上的Y数值
double Get_Z(void) {return( m_Point[2] ); }
int Get_Neighbor_Count (void){ return( m_nNeighbors ); } //邻接点的个数
TIN_Point *Get_Neighbor (int iNeighbor) { return( iNeighbor >= 0 && iNeighbor < m_nNeighbors ? m_Neighbors[iNeighbor] : NULL ); } //得到某一邻接点
int Get_Triangle_Count (void){return( m_nTriangles );} //得到邻接三角形个数
TIN_Triangle * Get_Triangle(int iTriangle) { return( iTriangle >= 0 && iTriangle < m_nTriangles ? m_Triangles[iTriangle] : NULL ); }
private:
TIN_Point(void); //构造函数
TIN_Point(int ID, POINT3d pptxyz);
virtual ~TIN_Point(void);
int m_ID, m_nNeighbors, m_nTriangles;//本身ID数值，邻接点个数，邻接三角形个数
POINT3d m_Point;//自身存储着XY数值
TIN_Point **m_Neighbors;//邻接顶点
TIN_Triangle **m_Triangles;//邻接三角形
bool _Add_Neighbor (TIN_Point *pNeighbor); //增加顶点
bool _Add_Triangle (TIN_Triangle *pTriangle); //增加三角形
bool _Del_Relations (void);
};
class AFX_EXT_CLASS TIN_Edge //TIN中的边
{ friend class TIN; //申明友元类
public: //得到边的端点
TIN_Point * Get_Point (int iPoint){ return( m_Points[iPoint % 2] ); }
private:
TIN_Edge(TIN_Point *a, TIN_Point *b); //构造函数
virtual ~TIN_Edge(void);
TIN_Point *m_Points[2]; //边的两个端点
};
4）以找到预处理后离散点集的最大外包三角形作为队列中第一个三角形开始，依次加入点集中的每个顶点，接着判断此新加顶点与已经存在的每个三角形的外接圆是否存在包含关系，如果包含，则改变此关联三角形的边信息后，然后以新加入的顶点和关联三角形的顶点为参数创建新的三角形，最后收尾工作是，以完成的三角化后每个三角形中的顶点为参数，建立TIN中的顶点、边和三角形之间的关系。
bool TIN::_Triangulate(TIN_Point **Points, int nPoints, TTIN_Triangle *Triangles, int &nTriangles)
{
int i, j, k, inside, trimax,
nedge = 0,
emax = 200,
status = 0,
*complete = NULL;
double xmid, ymid, dmax,
xp, yp, x1, y1, x2, y2, x3, y3, xc, yc, r;
TTIN_Edge *edges = NULL; //TTIN_Edge:边的结构
//-----------------------------------------------------
// Allocate memory for the completeness list, flag for each triangle //
trimax = 4 * nPoints;
if( (complete = (int *)malloc(trimax * sizeof(int))) == NULL ) //给即将生成的三角形的标志分配内存
{
status = 1;
if( edges ) free(edges); //如果不为空，首先释放内存
if( complete ) free(complete);
}
//-----------------------------------------------------
// Allocate memory for the edge list
if( (edges = (TTIN_Edge *)malloc(emax * sizeof(TTIN_Edge))) == NULL ) //给边的标志分配内存
{
status = 2;
if( edges ) free(edges);
if( complete ) free(complete);
}
//准备找到最大外包三角形作为三角形队列中第一个三角形,并准备不断的添加新顶点
_Extent_Update(); //更新整个离散顶点的矩形范围
dmax= m_Extent.Get_XRange() > m_Extent.Get_YRange() ? m_Extent.Get_XRange() : m_Extent.Get_YRange(); //取范围中XY两轴的较大数值
xmid= m_Extent.Get_XCenter(); //中心点
ymid= m_Extent.Get_YCenter();
Points[nPoints + 0]->m_Point[0] = xmid - 20 * dmax;//第一个顶点的X数值
Points[nPoints + 0]->m_Point[1] = ymid - dmax;//第一个顶点的Y数值，此数值取得相当小
Points[nPoints + 1]->m_Point[0] = xmid; //第二个顶点的X数值
Points[nPoints + 1]->m_Point[1] = ymid + 20 * dmax; //第二个顶点的Y数值，居中，但偏高
Points[nPoints + 2]->m_Point[0] = xmid + 20 * dmax; //第三个顶点的X数值
Points[nPoints + 2]->m_Point[1] = ymid - dmax;//第三个顶点的Y数值，第三个顶点与第一个顶点高度相同，在X轴上对称
Triangles[0].p1 = nPoints + 0; //三角形的顶点记录的是顶点数组中的索引
Triangles[0].p2 = nPoints + 1;
Triangles[0].p3 = nPoints + 2;
complete [0]= false;
nTriangles = 1; //初始化三角形个数为1
//-----------------------------------------------------
//每次增加一个顶点到已存在的TIN格网中
for(i=0; i<nPoints; i++) //对于所有顶点进行循环
{
xp = Points[i]->m_Point[0]; //使用浮点数进行记录
yp = Points[i]->m_Point[1]; //使用浮点数进行记录
nedge = 0; //对于每个顶点，初始化边的个数为0
//如果顶点位于三角形的外接圆内，则原三角形的三条边被加到一个缓存中，并移走此三角形
for(j=0; j<nTriangles; j++) //对于所有刚刚建立的三角形进行循环
{
if( complete[j] )
{
continue; //此三角形的完整性不用破坏，继续下次循环
}
x1 = Points[Triangles[j].p1]->m_Point[0];//三角形的三个顶点
y1 = Points[Triangles[j].p1]->m_Point[1];
x2 = Points[Triangles[j].p2]->m_Point[0];
y2 = Points[Triangles[j].p2]->m_Point[1];
x3 = Points[Triangles[j].p3]->m_Point[0];
y3 = Points[Triangles[j].p3]->m_Point[1];
inside = _CircumCircle(xp, yp, x1, y1, x2, y2, x3, y3, &xc, &yc, &r); //判断顶点与外接圆的关系
if( xc + r < xp ) //顶点落在外接圆之外
{
complete[j] = true; //标记此三角形为完整性
}
if( inside ) //顶点落在外接圆里面
{
//检查是否超出边链表的大小
if( nedge + 3 >= emax )
{
emax+= 100; //增加边的个数
if( (edges = (TTIN_Edge *)realloc(edges, emax * sizeof(TTIN_Edge))) == NULL ) //增加新的内存分配
{
status = 3;
if( edges ) free(edges);
if( complete ) free(complete);
}
}
edges[nedge + 0].p1 = Triangles[j].p1; //存储边的信息
edges[nedge + 0].p2 = Triangles[j].p2; //顶点的索引：1,2; 2,3; 3,1
edges[nedge + 1].p1 = Triangles[j].p2;
edges[nedge + 1].p2 = Triangles[j].p3;
edges[nedge + 2].p1 = Triangles[j].p3;
edges[nedge + 2].p2 = Triangles[j].p1;
nedge += 3; //对于最外层循环的顶点来说，相关边的个数加上3
Triangles[j]= Triangles[nTriangles - 1]; /将链表中最后一个三角形赋值过来
complete [j]= complete [nTriangles - 1];
nTriangles--; //三角形个数减去
j--;
} //对应：if( inside )
} //结束对所有三角形的循环
//标记被重复使用的边
//若三角形以逆时针方向存储，则内部被多次使用的边存储的顶点顺序是相反的
for(j=0; j<nedge-1; j++)
{
for(k=j+1; k<nedge; k++)
{
if( (edges[j].p1 == edges[k].p2) && (edges[j].p2 == edges[k].p1) ) //同样的一条边存储的方向刚刚相反
{
edges[j].p1 = -1; //做标记
edges[j].p2 = -1;
edges[k].p1 = -1;
edges[k].p2 = -1;
}
if( (edges[j].p1 == edges[k].p1) && (edges[j].p2 == edges[k].p2) ) //同样的一条边
{
edges[j].p1 = -1;
edges[j].p2 = -1;
edges[k].p1 = -1;
edges[k].p2 = -1;
}
}
}
//组成新三角形
//忽略所有已经打过标志的边
//所有的边以顺时针方向安排
for(j=0; j<nedge; j++)
{
if( edges[j].p1 < 0 || edges[j].p2 < 0 ) //比如等于-1的情况
{
continue;
}
if( nTriangles >= trimax )
{
status = 4;
if( edges ) free(edges);
if( complete ) free(complete);
}
Triangles[nTriangles].p1 = edges[j].p1; //记录关联的边的一个顶点索引
Triangles[nTriangles].p2 = edges[j].p2;//记记录关联的边的另一个顶点索引
Triangles[nTriangles].p3 = i;//对于最外层循环的顶点来说，记录此顶点的索引
complete [nTriangles] = false;
nTriangles++;
}
} //对所有的顶点循环完毕 //对应着for(i=0;i<;i++)
//-----------------------------------------------------
// Remove triangles with supertriangle vertices //移走初始定义的大外接三角形
// These are triangles which have a vertex number greater than nPoints //这些三角形的顶点索引超出了顶点个数
for(i=0; i<nTriangles; i++)
{
if(Triangles[i].p1 >= nPoints
|| Triangles[i].p2 >= nPoints
|| Triangles[i].p3 >= nPoints )
{
Triangles[i] = Triangles[nTriangles - 1]; //不断的把最后一个三角形赋值过来
nTriangles--;
i--;
}
}
return( status == 0 );
}
5）对TIN中的三角形进行循环，以每个三角形的顶点为参数，插值生成GRID DEM。参见图四.
图一:等高线图
图二:离散点
图三:导出的数据存储在文本中
图四:生成的GRID DEM