1) 针对格网中每个基本单元，创建一个参考变量，即YuanSu[i][j](0<=i<=Row; 0<=j<=Column)同时分配内存。此参考变量赋值0或1，其中1代表此位置处的邻域高程数值将发生变化，那么此位置将绘制等高线。
col = (char **)calloc(m_iRows, sizeof(char *)); //指向指针的指针初始化
row = (char **)calloc(m_iRows, sizeof(char *)); //指向指针的指针初始化
for(y=0; y<m_iRows; y++) //m_iRows:格网的行数
{
col[y] = (char *)calloc(m_iColumns, sizeof(char)); //每个单元分配内存，m_iColumns : 格网的列数
row[y] = (char *)calloc(m_iColumns, sizeof(char)); //每个单元分配内存
}
2）在以高程数值从最小值开始依次递增的最外部一层循开始，接着对整个格网数据进行循环，在给参考变量YuanSu[i][j]赋数值0或1之后，然后以YuanSu[i][j] == 1为判断条件，寻找并记录等高线。
for(zValue=zMin, ID=0; zValue<=zMax; zValue+=zStep) //从最小Z数值计算起，知道最大Z数值
{
//-------------------------------------------------
// Step1: Find Border Cells //找到边界单元，边界：高度趋势发生改变的地方。
for(y=0; y<m_iRows-1; y++)
for(x=0; x<m_iColumns-1; x++)
if( m_p3dDEMPoints[y*m_iColumns+x][2] >= zValue )
{
row[y][x] = m_p3dDEMPoints[y*m_iColumns+(x+1)][2] < zValue ? 1 : 0; //‘’就代表边界
col[y][x] = m_p3dDEMPoints[(y+1)*m_iColumns+x][2] < zValue ? 1 : 0;
}
else
{
row[y][x] = m_p3dDEMPoints[y*m_iColumns+(x+1)][2] >= zValue ? 1 : 0; //‘’就代表边界
col[y][x] = m_p3dDEMPoints[(y+1)*m_iColumns+x][2] >= zValue ? 1 : 0;
}
//-------------------------------------------------
// Step2: Interpolation + Delineation
for(y=0; y<m_iRows-1; y++)
for(x=0; x<m_iColumns-1; x++)
{
if(row[y][x]) //如果数值为
{
for(i=0; i<2; i++) //两次调用函数，ID值增加
{
if(i==0) m_iCountourfind = 0;
else m_iCountourfind = 1;
ContourFind(x, y, zValue, true, ID++);
}
row[y][x] = 0; //代表已经处理过
}
if(col[y][x]) //如果数值为
{
for(i=0; i<2; i++) //两次调用函数，ID值继续增加
{
if(i==0) m_iCountourfind = 0;
else m_iCountourfind = 1;
ContourFind(x, y, zValue, false, ID++);
}
col[y][x] = 0; //代表已经处理过
}
}
}
3）寻找等高线的过程，是在本位置的高程数值，本条等高线的高程ZValue和邻近位置的高程数值做插值处理后，得到X、Y位置并做记录。
do
{
//-------------------------------------------------
// Interpolation...
d = m_p3dDEMPoints[y*m_iColumns+x][2]; //取本位置的高程数值
d = (d - z) / (d - m_p3dDEMPoints[zy*m_iColumns+zx][2]); //以本位置的高程数值，ZValue和邻近位置的高程数值做插值处理
xPos = m_dXMin + m_iXCellSize * (x + d * (zx - x)); //通过高度差的比例，线性插值求得X、Y位置
yPos = m_dYMin + m_iYCellSize * (y + d * (zy - y)); //
if(m_iCountourfind == 0)
{
if (m_iPointNum_one >= m_iBufferSize_one)
{
m_iBufferSize_one += 1024;
m_point3dtmp_one = (POINT3d *)realloc(m_point3dtmp_one,m_iBufferSize_one * sizeof(POINT3d)); //为顶点增加新的内存分配
}
if(m_point3dtmp_one!=NULL)
{
m_point3dtmp_one[m_iPointNum_one][0] = xPos; //记录下X位置
m_point3dtmp_one[m_iPointNum_one][1] = yPos; //记录下Y位置
m_point3dtmp_one[m_iPointNum_one][2] = z+0.2; //此处高度为ZValue,增加.2个单位的目的是为了显示时抬高一些，便于观察。
m_iPointNum_one++;
}
}
else
{
if (m_iPointNum_two >= m_iBufferSize_two)
{
m_iBufferSize_two += 1024;
m_point3dtmp_two = (POINT3d *)realloc(m_point3dtmp_two,m_iBufferSize_two * sizeof(POINT3d)); //为顶点增加新的内存分配
}
if(m_point3dtmp_two!=NULL)
{
m_point3dtmp_two[m_iPointNum_two][0] = xPos; //记录下X位置
m_point3dtmp_two[m_iPointNum_two][1] = yPos; //记录下Y位置
m_point3dtmp_two[m_iPointNum_two][2] = z+0.2; //此处高度为ZValue,增加.2个单位的目的是为了显示时抬高一些，便于观察。
m_iPointNum_two++;
}
}
//-------------------------------------------------
// Naechstes (x/y) (col/row) finden...
if( !ContourFindNext(Dir, x ,y ,doRow) ) //x、y将发生增或减的变化，doRow也改变数值
doContinue = ContourFindNext(Dir, x, y, doRow);
Dir = (Dir + 5) % 8;
//-------------------------------------------------
// Loeschen und initialisieren...
if(doRow)
{
row[y][x] = 0; //表示已经处理过
zx = x + 1; //重新计算数值
zy = y;
}
else
{
col[y][x] = 0; //表示已经处理过
zx = x; //重新计算数值
zy = y + 1;
}
}
while(doContinue);
4）寻找等高线过程结束的条件：
bool doContinue;
if(doRow)
{
switch(Dir) //邻域上个方向的判断
{
case 0:// Norden
if(row[y+1][x ])
{ y++; Dir=0; doContinue=true; break; }
case 1:// Nord-Ost
if(col[y ][x+1])
{ x++; doRow=false;Dir=1; doContinue=true; break; }
case 2:// Osten ist nicht...
case 3:// Sued-Ost
if(y-1>=0) if(col[y-1][x+1])
{ x++; y--; doRow=false;Dir=3; doContinue=true; break; }
case 4:// Sueden
if(y-1>=0) if(row[y-1][x ])
{ y--; Dir=4; doContinue=true; break; }
case 5:// Sued-West
if(y-1>=0) if(col[y-1][x ])
{ y--; doRow=false;Dir=5; doContinue=true; break; }
case 6:// Westen ist nicht...
case 7:// Nord-West
if(col[y ][x ])
{ doRow=false;Dir=7; doContinue=true; break; }
default:
Dir=0; doContinue=false;
}
}
else
{
switch(Dir) //邻域上个方向的判断
{
case 0:// Norden ist nicht...
case 1:// Nord-Ost
if(row[y+1][x ])
{ y++; doRow=true; Dir=1; doContinue=true; break; }
case 2:// Osten
if(col[y ][x+1])
{ x++; Dir=2; doContinue=true; break; }
case 3:// Sued-Ost
if(row[y ][x ])
{ doRow=true; Dir=3; doContinue=true; break; }
case 4:// Sueden ist nicht...
case 5:// Sued-West
if(x-1>=0) if(row[y ][x-1])
{ x--; doRow=true; Dir=5; doContinue=true; break; }
case 6:// Westen
if(x-1>=0) if(col[y ][x-1])
{ x--; Dir=6; doContinue=true; break; }
case 7:// Nord-West
if(x-1>=0) if(row[y+1][x-1])
{ x--; y++; doRow=true; Dir=7; doContinue=true; break; }
default:
Dir=0; doContinue=false;
}
}
return(doContinue);
原始DEM显示图:
生成的等高线图: