基于OpenGL的三维地形的构建

qq_46029700

已于 2023-07-23 19:20:04 修改

阅读量6k

点赞数 14

分类专栏：笔记文章标签： arcgis c++

于 2023-07-23 18:41:01 首次发布

本文链接：https://blog.csdn.net/qq_46029700/article/details/131882344

版权

笔记专栏收录该内容

3 篇文章 0 订阅

订阅专栏

（1）在arcmap里面加载DEM数据后，使用【创建渔网】工具，设置模板范围为DEM范围，像元宽高均设置为20m。勾选创建标注点。如图1所示。

（2）使用【值提取至点】工具，输入点要素为刚刚创建的渔网点，提取栅格值为DEM。如图2

（3）对提取后的点数据【添加XY坐标】，如图3

（4）去掉边界的点数据，首先选中不是边界的点数据，然后重新导出，因为边界点数据在【值提取至点】时可能提取不到。如图4

（5）遥感影像纹理准备，这里添加在线地图影像，使用【绘图】工具绘制出高程点的矩形框范围，然后截图。如图5所示

（6）另外计算出的x,y投影坐标较大，这里进行中心化。中心化公式如下：

（7）最后导出该属性表为excel，并拷贝最后的结构体形式字段到工程头文件中。数据预处理完成

1、目的要求

（1）地形高程的随机生成

（2）地形平滑技术

（3）根据高程确定地形的颜色

（4）地形的动态生成

（5）地形的绘制

2、实验内容

三维地形动态生成及模拟程序设计

3、实验原理及核心代码

本次实验思路框架如下：

3.1 地形数据的预处理

3.1.2原理

我们知道在opengl里面是可以进行三维图形的绘制的，三维图形是由一个个顶点和边组成的，顶点之间相连便可组成边，进而形成面。所以顶点数据尤为重要。这里最终需要数据形式为三维点数组来存储点数据。数据形式如下：

三维地形原始数据结构

PointZ terri[ROW*COL]

={PoinZ00,PointZ01,PoinZ02,…PoinZij,…PointZrow-1,col-1}；

其中ROW和COL为块状栅格的行列数目，在c++里面这里是采用一维数组来存储二维数组数据，第i行j列寻址结果为terri[i*COL+j]。另外PointZ为三维点结构体，其包含x,y,h三个坐标:

三维点数据结构

struct PointZ{

GLfloat x;//x坐标

GLfloat y;// y坐标

GLfloat h;//高程h

};

3.1.2数据处理操作

（1）在arcmap里面加载DEM数据后，使用【创建渔网】工具，设置模板范围为DEM范围，像元宽高均设置为20m。勾选创建标注点。如图1所示。

图1 创建渔网

（2）使用【值提取至点】工具，输入点要素为刚刚创建的渔网点，提取栅格值为DEM。如图2

图2 值提取至点

（3）对提取后的点数据【添加XY坐标】，如图3

图3 添加XY坐标

（4）去掉边界的点数据，首先选中不是边界的点数据，然后重新导出，因为边界点数据在【值提取至点】时可能提取不到。如图4

图4 去掉外围点数据

（5）遥感影像纹理准备，这里添加在线地图影像，使用【绘图】工具绘制出高程点的矩形框范围，然后截图。如图5所示

图5 截图在线地图

（6）另外计算出的x,y投影坐标较大，这里进行中心化。中心化公式如下：

式中xi为数据组{x1,x2,x3,…,xi,…xn}中的数据，通过其减去数据组数据的平均值，保留其相对波动是构造地形一种比较好的方式，本身在地理学中，高程通常讲究的是在某一基准面下的高程——即相对高差。在arcmap中通过【添加字段】和【字段计算器】进行计算。计算过程如图6，7所示。

图6 统计平均值

图7 字段计算器输入计算公式

（7）注意这里得到的x,y,h其实与实际opengl里面绘制时有些许区别，如下图7所示，虚线的四边形为栅格，可以看到预处理后的高程点数据的x,y,h与opengl里面的对应关系为：x->x y->z h->y，此外还需要对y值进行取反操作。对字段y取反后进行最后输入之前所讲结构体的形式的字段，如图8所示

图7 arcmap与opengl绘制对应起来

图8 计算最后的x,y,h

（7）最后导出该属性表为excel，并拷贝最后的结构体形式字段到工程头文件中。数据预处理完成

3.2 地形的生成

3.2.1原理

虽然我们存储二维高程点数据采用的是一维数组数据结构，但是我们进行分析时可以从其逻辑结构进行分析。如图9所示，每个格网的折点就是高程点数据所在位置，取其中任意一个i行j列的点进行地形的构造(这里i从0开始,j从1开始，可能在实际生活中我们称为i+1行,j+1列，问题不大，注意一下就行)，每个格网被分成了2个三角形，三角形绘制方向在图9中已经给出。

图9 地形生成原理示意图

写成算法伪代码的形式如下：

三维地形生成算法

for(int i=0;i<ROW-1;i++){

for(int j=0;j<COL-1;j++){

glBegin(GL_TRIANGLES);

//第一个三角形

glVertex3fv(PointZij);//第I 行j 列点

glVertex3fv(PointZi+1,j);//第I＋1行j 列点

glVertex3fv(PointZi+1,j+1);//第I＋1行j＋1列点

//第二个三角形

glVertex3fv(PointZi+1,j+1);//第I＋1行j＋1列点

glVertex3fv(PointZi,j+1);//第I 行j＋1列点

glVertex3fv(PointZi,j);//第I 行j 列点

glEnd();

}

3.2.2核心代码

根据上述思路，这里绘制线和三角面两种情况。核心代码如下：

三维地形生成核心代码

case WIRE:

{

glColor4f(0,1,0,1);

for(int i=0;i<ROW-1;i++){

for(int j=0;j<COL-1;j++){

glBegin(GL_LINE_LOOP);

glVertex3f(terri[i*COL+j].x,terri[i*COL+j].h,terri[i*COL+j].z);//i,j

glVertex3f(terri[(i+1)*COL+j].x,terri[(i+1)*COL+j].h,terri[(i+1)*COL+j].z);//i+1,j glVertex3f(terri[(i+1)*COL+j+1].x,terri[(i+1)*COL+j+1].h,terri[(i+1)*COL+j+1].z);//i+1,j+1

glEnd();

glBegin(GL_LINE_LOOP);

glVertex3f(terri[(i+1)*COL+j+1].x,terri[(i+1)*COL+j+1].h,terri[(i+1)*COL+j+1].z);//i+1,j+1

glVertex3f(terri[i*COL+j+1].x,terri[i*COL+j+1].h,terri[i*COL+j+1].z);//i,j+1

glVertex3f(terri[i*COL+j].x,terri[i*COL+j].h,terri[i*COL+j].z);//i,j

glEnd();

}

break;

case SOLID:

{

glColor3f(0.5,0.5,0.5);

for(int i=0;i<ROW-1;i++){

for(int j=0;j<COL-1;j++){

glBegin(GL_TRIANGLES);

//第一个三角形

glVertex3f(terri[i*COL+j].x,terri[i*COL+j].h,terri[i*COL+j].z);//第I 行j 列点 glVertex3f(terri[(i+1)*COL+j].x,terri[(i+1)*COL+j].h,terri[(i+1)*COL+j].z);//第I＋1行j 列点 glVertex3f(terri[(i+1)*COL+j+1].x,terri[(i+1)*COL+j+1].h,terri[(i+1)*COL+j+1].z);//第I＋1行j＋1列点

//第二个三角形

glVertex3f(terri[(i+1)*COL+j+1].x,terri[(i+1)*COL+j+1].h,terri[(i+1)*COL+j+1].z);//第I＋1行j＋1列点 glVertex3f(terri[i*COL+j+1].x,terri[i*COL+j+1].h,terri[i*COL+j+1].z);//第I 行j＋1列点

glVertex3f(terri[i*COL+j].x,terri[i*COL+j].h,terri[i*COL+j].z);//第I 行j 列点

glEnd();

}

break;

3.2.3效果展示

图10 三角网线框地形

图11 三角网平面地形

3.3 真实感地形与光照模型

3.2节在展示地形时，三角网平面展示效果不是很好，这是因为缺乏光照的效果。本节将进行光照的添加。

3.3.1原理

opengl是采用的phone光照模型，然后设置光照一般需要设置光源、物体材质、物体各平面的法向量。

（1）光源：光源的设置主要通过glLightfv( GL_LIGHT0, prama, prama2 )进行设置，然后启动就行，必要时也可以通过glLightModelfv(GL_LIGHT_MODE-L_AMBIENT,lmodel_ambilent)设置全局环境光(不过我们本次不进行设置)。为达到比较好的光照效果，本次采用光源衰减效果。光源衰减方程如下：

我们需要设置Kc,Kl和Kq3个参数的值，使用的函数为glLightf(GL_LIGHT1, parama, value);具体设置多大，我们可以根据Ogre3D的Wiki给出的参考[1]。如下表1.本实验中Kc=1.0，Kl=0.0007，Kq=0.0002。

表1 光源衰减参考

另外本次实验中，光源会进行移动，照射角度会发生改变，一般在进行平移旋转后调用glLightfv(GL_LIGHT0,GL_POSITION,value)就行，本次实验通过旋转角度变量来改变光源的位置。

（2）材质：在光照条件下，一般会掩盖绘制物体的色彩，比如glColor4f会失效，不过可以通过设置glEnable(GL_COLOR_MATERIAL)，来使物体的绘制色彩代替物体的材质，一般紧跟上glLightModeli( GL_FRONT, GL_AMBIENT-_AND_D-IFFUSE )。此外还可以通过传统方式glMaterialfv(GL_FRONT, prama, prama2);来进行设置。在创建材质时为使后续调用方便，也可以创建材质列表。

（3）法向量：下面介绍地形表面法向量的计算，因为是三角平面，可以通过向量叉积来得到平面的法向量，但是要注意方向。如下图12所示三角平面ABC，三个顶点p1,p2,p3。

图12 三角平面

向量AC、AB，AB×AC = 朝向屏幕外的法向量；AC×AB = 朝向屏幕内的法向量；如果p1～p3属于3维空间，则具体计算步骤为：

用a表示向量AB，用b表示向量AC，则a=(x1,y1,z1)；b=(x2,y2,z2);

3.3.2核心代码

光源设置核心代码

//白天阳光

GLfloat white[] = { 1, 1, 1, 1.0 };

GLfloat gray[]={0.1,0.1,0.1,1.0};

glLightfv( GL_LIGHT0, GL_POSITION, lightPos );

glLightfv( GL_LIGHT0, GL_SPECULAR, white );

glLightfv( GL_LIGHT0, GL_DIFFUSE, white );

glLightfv(GL_LIGHT0, GL_AMBIENT, gray);

glLightf(GL_LIGHT0, GL_CONSTANT_ATTENUATION, 1.0); // c 系数

glLightf(GL_LIGHT0, GL_LINEAR_ATTENUATION, 0.0007); // l 系数

glLightf(GL_LIGHT0, GL_QUADRATIC_ATTENUATION, 0.0002); // q 系数

//夜晚月光

GLfloat white2[] = { 0.4, 0.4, 0.4, 1.0 };

GLfloat gray2[]={0.2,0.2,0.2,1.0};

glLightfv( GL_LIGHT1, GL_POSITION, lightPos );

glLightfv( GL_LIGHT1, GL_SPECULAR, white2 );

glLightfv( GL_LIGHT1, GL_DIFFUSE, white2 );

glLightfv(GL_LIGHT1, GL_AMBIENT, gray2);

//衰减系数参考https://learnopengl-cn.github.io/02%20Lighting/05%20Light%20casters/#_4

glLightf(GL_LIGHT1, GL_CONSTANT_ATTENUATION, 1.0); // c 系数

glLightf(GL_LIGHT1, GL_LINEAR_ATTENUATION, 0.0007); // l 系数

glLightf(GL_LIGHT1, GL_QUADRATIC_ATTENUATION, 0.0002); // q 系数

glEnable( GL_DEPTH_TEST );

光源移动

glRotatef(angle, 0.0f, 0.0f, 1.0f);

glTranslatef(lightPos[0],lightPos[1],lightPos[2]);

glLightfv(GL_LIGHT0,GL_POSITION,lightPos);

材质设置核心代码

//常规方法

GLfloat terri_mat_ambient[] = {0.1f, 0.1f, 0.1f, 1.0f};

GLfloat terri_mat_diffuse[] = {0.5f, 0.5f, 0.5f, 1.0f};

GLfloat terri_mat_specular[] = {0.0f, 0.0f, 0.0f, 1.0f};

GLfloat terri_mat_emission[] = {0.f, 0.0f, 0.0f, 1.0f};

GLfloat terri_mat_shininess = 0.0f;

TerriMat=glGenLists(1);

glNewList(TerriMat, GL_COMPILE);

glMaterialfv(GL_FRONT, GL_AMBIENT, terri_mat_ambient);

glMaterialfv(GL_FRONT, GL_DIFFUSE, terri_mat_diffuse);

glMaterialfv(GL_FRONT, GL_SPECULAR, terri_mat_specular);

glMaterialfv(GL_FRONT, GL_EMISSION, terri_mat_emission);

glMaterialf (GL_FRONT, GL_SHININESS, terri_mat_shininess);

glEndList();

//用绘制色彩代替

glCallList(ColorMat);//回归默认材质

glEnable(GL_COLOR_MATERIAL);

glLightModeli( GL_FRONT, GL_AMBIENT_AND_DIFFUSE );

glColor4f(0,1,0,1);

计算平面法向量核心代码

//2.2计算三角平面的法向量，返回一个GLfloat数组

GLfloat* calNormal(PointZ p1,PointZ p2,PointZ p3){

double x1,y1,z1,x2,y2,z2;

x1=p2.x-p1.x;

y1=p2.h-p1.h;

z1=p2.z-p1.z;

x2=p3.x-p1.x;

y2=p3.h-p1.h;

z2=p3.z-p1.z;

GLfloat a[3]={x1,y1,z1};

GLfloat b[3]={x2,y2,z2};

GLfloat result[3];

result[0]=y1*z2-y2*z1;

result[1]=-(x1*z2-x2*z1);

result[2]=x1*y2-x2*y1;

return result;

}

//在绘制三角形面前使用glNormal3fv调用上述函数就行

3.3.3效果展示

图13三角网线框(夜晚)

图14三角网平面(白天)

3.4 地形的分级色彩

地形有高有低，如果用不同色彩代表不同的高程值，会让人更加直观的体验到高程的变化。

3.4.1原理

地形分级色彩的实现不是很难，其实思路就是创建色带——>然后根据三角面平均高程获取色带索引号。

（1）色带的创建：本次色带没有写死，自己写的一个很简单的色带变化，给定分类数目Level，然后色彩(RGBA)由低高程绿(0,1,0,1)过渡到红(1,0,0,1)，R值线性增加,G值线性减少。每次增加或减少量为1/Level。最后生成一个色带数组，后面高程根据自己的范围选择色带索引。

（2）如何根据高程选择色带索引：首先根据该三角面的三个顶点计算他的平均高程，然后这里举个例子来说明如何得到色带索引：假设最高高程为99，最低高程为-62。假设分5类时，颜色高程分度值为(99-62)/5=32.2；

如果某一面计算出平均高程-60，与最低高程插值为2，2与颜色分度值32.2相除取整为0，取色带编号0。

如果某一面计算出平均高程-29.7，与最低高程差值为32.3，32.3/32.2=1，取色带编号1。

如果某一面计算出平均高程99,与最低高差161，161/32.2=5.0，这里可能会出现浮点数的问题，这里只能去4，所以统一在得到色带索引编号时,如果大于色带最大索引，就取色带最大索引(这里为Level-1),如果小于色带最小索引，就取色带最小索引，这里为0；

（3）小结：之前尝试过写死色带，和用if-else去取色带编号，感觉代码扩展性不强，采用上述方法可扩展性较好，并且后面可直接把重心放到色带变化的函数上去。

3.4.2核心代码

色带初始化核心代码

//5.2色带初始化

void initColorRamp(){

//从(0,1,0)过渡到(1,0,0)吧，就使用最简单的方式，根据LEVEL数目，G每次减少1/LEVEL,R每次增加1/LEVEL

//GLfloat colorDivision=1/LEVEL;

for(int i=0;i<LEVEL;i++){

colorRamp[i][0]=i/(GLfloat)LEVEL;

colorRamp[i][1]=1-i/(GLfloat)LEVEL;

colorRamp[i][2]=0;

colorRamp[i][3]=1;

}

根据高程获取色带索引编号核心代码

//2.4根据平均高程设置渲染色彩

GLfloat* setLevelColor(GLfloat h){

/*if(h<(MIN_HEIGHT+colorDivision)) return colorRamp[0];

else if(h<(MIN_HEIGHT+colorDivision*1)) return colorRamp[0];*///这种效率太低，根据与最低高程的差值与颜色分度值的除运算进行

//例如分5类时，颜色分度值为(99--62)/5=32.2；高程-60，与最低高程插值为2，2与颜色分度值32.2，2/32.2取色带编号0

//高程-29.7，与最低高程差值为32.3,32.3/32.2=1,取色带编号1

//高程99,与最低高差161，161/32.2=5.0，这里可能会出现浮点数的问题，这里只能去4，所以统一在得到index时,如果大于LEVEL-1，就取LEVEL-1,如果小于0，就取0；

int index=((h-MIN_HEIGHT)/heightColorDivision);

if(index>LEVEL-1) index=LEVEL-1;

if(index<0) index=0;

return colorRamp[index];

}

3.4.3效果展示

图15 三角网线框分级色彩(夜晚)

图16 三角网平面分级色彩(白天)

3.5 纹理映射之遥感影像融合

利用opengl的纹理映射技术，可把研究区域的遥感影象与地形融合在一起，取得更为逼真的三维视觉效果[2]。

3.5.1原理

（1）纹理映射的一般步骤：

纹理定义(glTexImage2D说明所映射为例的内容，包括纹理数据的指针、纹理大小、纹理类别等)

纹理控制(glTexParameteri，说明纹理以何种方式映射到模型表面上，opengl提供有有纹理滤波、重复和伸缩等等)

纹理映射方式说明(glTexEnvf，用来调整三维模型色彩或者将纹理与三维模型原来的色彩进行融合)

启用关闭纹理映射(glEnable(GL_TEXTURE_2D)，glDiable(GL_TEXTU-RE_2D))

三维地形高程点纹理坐标与几何坐标的对应(纹理坐标通过glTexCoord指定)[3]

本次重点对纹理坐标的计算和纹理映射方式进行说明。

（2）纹理坐标的计算：在 OpenGL 中的纹理坐标空间（ u，v ）为u：0.0～1.0；v ：0.0～1.0，则数字地形与影象融合的纹理映射过程，其关键就是计算出高程点对应的纹理空间（ u，v）中的坐标。根据实验情况，纹理映射的过程可以由下图17表示

图17 纹理映射的关键

该图需要注意地形高程点的最低值与遥感影像的最低值位置不一样，主要是纵方向，横方向是一致的。了解这个情况后，我们可以用比例来计算。计算步骤如下：

a. 地形高程点横向共分为COL-1等份，所以遥感影像纹理也要分为COL-1等分，每一等份分度值为1.0/(COL-1)。同理，纵向遥感影像纹理分度值为1.0/(ROW-1)等份。

b. 因为对于i行j列的高程点，其向下移动了i行,i个分度值，因此纹理的纵坐标为1.0-i/(ROW-1)；横坐标为i/(COL-1)

（2）纹理光照的融合：纹理和光照的融合函数默认是glTexEnvf(GL_TEX-TURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE)，该效果会将光照的效果给隐藏掉，我们为了看到遥感影像叠加上去后能反映光源的位置影像，这里可以将GL_REPLACE替换为GL_BLEND。

3.5.2核心代码

纹理载入代码较模式化，这里不再引入。

遥感影像纹理融合核心代码

case RS:

{

glEnable(GL_TEXTURE_2D);

glCallList(RSMat);

glBindTexture(GL_TEXTURE_2D, texObj[RSIMG]);//贴遥感影像时，需要将遥感影像根据DEM行(ROW-1)列(COL-1)数进行分割

//glColor3f(0.5,0.5,0.5);

for(int i=0;i<ROW-1;i++){

for(int j=0;j<COL-1;j++){

glBegin(GL_TRIANGLES);

//第一个三角形

//glNormal3fv(calNormal(terri[i*COL+j],terri[(i+1)*COL+j],terri[i*COL+j+1])); //法向量1

glTexCoord2f(j*xFactor,1-i*yFactor);

glVertex3f(terri[i*COL+j].x,terri[i*COL+j].h,terri[i*COL+j].z);//第I 行j 列点

glTexCoord2f(j*xFactor,1-(i+1)*yFactor);

glVertex3f(terri[(i+1)*COL+j].x,terri[(i+1)*COL+j].h,terri[(i+1)*COL+j].z);//第I＋1行j 列点

glTexCoord2f((j+1)*xFactor,1-i*yFactor); glVertex3f(terri[i*COL+j+1].x,terri[i*COL+j+1].h,terri[i*COL+j+1].z);//第I 行j＋1列点

//第二个三角形 //glNormal3fv(calNormal(terri[i*COL+j],terri[(i+1)*COL+j],terri[i*COL+j+1])); //法向量2

glTexCoord2f(j*xFactor,1-(i+1)*yFactor); glVertex3f(terri[(i+1)*COL+j].x,terri[(i+1)*COL+j].h,terri[(i+1)*COL+j].z);//第I＋1行j 列点

glTexCoord2f((j+1)*xFactor,1-(i+1)*yFactor); glVertex3f(terri[(i+1)*COL+j+1].x,terri[(i+1)*COL+j+1].h,terri[(i+1)*COL+j+1].z);//第I＋1行j＋1列点

glTexCoord2f((j+1)*xFactor,1-i*yFactor);

glVertex3f(terri[i*COL+j+1].x,terri[i*COL+j+1].h,terri[i*COL+j+1].z);//第I 行j＋1列点

glEnd();

}

glDisable(GL_TEXTURE_2D);

}

break;

3.5.3效果展示

图18 地形融合遥感影像

3.6 动画虚拟之洪水模拟

3.6.1原理

洪水模拟主要是利用opengl的双缓存技术进行动画的效果，利用<time.h>时间库函数来控制时间进而控制洪水立方体的绘制高度。具体实现示意如下图19，图20。

图19 洪水开始模拟前

图20 洪水模拟某一时刻洪水立方体的绘制

另外洪水在地形上需要展示一种半透明的效果。我们需要根据混合的知识来调整绘制顺序——就是：首先绘制所有不透明的物体。如果两个物体都是不透明的，则谁先谁后都没有关系。然后，将深度缓冲区设置为只读。接下来，绘制所有半透明的物体。

3.6.2核心代码

洪水模拟核心代码

//绘制洪水立方体

void drawWater(){

if(height!=MIN_HEIGHT){

GLdouble h=height-MIN_HEIGHT;//水深

glPushAttrib(GL_CURRENT_BIT);

glEnable(GL_COLOR_MATERIAL);

glLightModeli( GL_FRONT, GL_AMBIENT_AND_DIFFUSE );

glDepthMask(GL_FALSE);//禁用深度缓存

glColor4f(0,0,1,0.5);

glTranslatef(0,MIN_HEIGHT+h/2,0);

glPushMatrix();

glScalef((MAX_DEMX-MIN_DEMX-1)/h,1,(MAX_DEMY-MIN_DEMY-1)/h);

glutSolidCube(h);

glPopMatrix();

glDepthMask(GL_TRUE);//启用深度缓存

glDisable(GL_COLOR_MATERIAL);

glPopAttrib();

}

//洪水模拟菜单

case START:

glEnable(GL_BLEND);

glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

height=MIN_HEIGHT;

while(height<MAX_HEIGHT){

height+=SPEED;

display();

Sleep(SPEED);

}

break;

case CLEAR:

glDisable(GL_BLEND);

height=MIN_HEIGHT;

glutPostRedisplay();

break;

3.6.3效果展示

图21 洪水模拟中间和结束状态

3.7 地形的动态浏览

3.7.1原理

动态浏览包括视角的切换和视口的变换，其中视角的切换包括改变观察方向和观察距离。通过按键或者鼠标来改变他们的参数然后重绘场景以达到地形动态浏览的效果。鼠标相比键盘实现稍微要动点脑筋，下面依次来阐述。

（1）视角的切换：即变换视角矩阵glLookAt函数。本次观察场景采用球体坐标系进行观察，如下图22所示。有3个重要的变量，观察半径g_r，与y轴形成的夹角g_phai，相机投影在xoz平面与x轴正向的夹角g_theta。观察半径与观察距离有关，另外两个角度与观察方向有关。可以通过设置键盘事件改变其值来进行地形的动态绘制。这里面具体glLookAt参数为： gluLookAt((g_r) * sin(g_phai) * cos(g_theta), (g_r) * cos(g_phai), (g_r) * sin(g_phai) * sin(g_theta), 0, 0, 0, 0, sin(g_phai) >= 0 ? 1 : -1, 0);

图22 视角变换的参考示意图

（2）视口的变换:视口的变换是投影坐标系到窗体坐标系的变换。这里主要是通过改变glViewport(viewx, viewy, w, h)的前两个参数来进行改变，viewx和viey默认值是0，0.通过改变其值，从而调整投影到视窗的位置。具体原理可以见图23.上图最开始是0,0下图是改成viewx和viewy后的结果。注意：一般情况下场景是占据整个窗体的，这里绘制起时是一般情形下。

图23 视口变换原理

（3）opengl的鼠标事件：上述变换变换参数通过键盘事件是很简单的，但通过鼠标事件稍微有点复杂，此外opengl默认是没有鼠标滚轮事件的。这里最后需要实现鼠标中键按下能改变观察半径，鼠标左键移动能改成观察方向的效果。

其中鼠标中键实现原理如图24所示，鼠标左键实现原理如图25所示。

图24 鼠标中键改变观察半径原理

图25 鼠标左键改变观察角度原理

3.7.2核心代码

按键内代码比较简单，这里只给下视口变换时的核心代码和鼠标事件的代码。

视口变换键盘核心代码

//7.1 通过特殊按键控制场景的进退平移

void keyboard2( int key, int x, int y )

{

switch( key ) {

case GLUT_KEY_LEFT:

{

viewx-=viewStep;

int w = glutGet (GLUT_WINDOW_WIDTH); int h = glutGet (GLUT_WINDOW_HEIGHT);

glViewport(viewx, viewy, w, h);//规定视口的大小（视口中图形全部填充满）

break;

}

case GLUT_KEY_RIGHT:

{

viewx+=viewStep;

int w = glutGet (GLUT_WINDOW_WIDTH); int h = glutGet (GLUT_WINDOW_HEIGHT);

glViewport(viewx, viewy, w, h);

break;

}

case GLUT_KEY_DOWN:

{

viewy-=viewStep;

int w = glutGet (GLUT_WINDOW_WIDTH); int h = glutGet (GLUT_WINDOW_HEIGHT);

glViewport(viewx, viewy, w, h);//规定视口的大小（视口中图形全部填充满）

break;

}

case GLUT_KEY_UP:

{

viewy+=viewStep;

int w = glutGet (GLUT_WINDOW_WIDTH); int h = glutGet (GLUT_WINDOW_HEIGHT);

glViewport(viewx, viewy, w, h);//规定视口的大小（视口中图形全部填充满）

break;

}

glutPostRedisplay();

}

鼠标事件核心代码

//6.1鼠标按下事件

void myMouse(int button,int state,int x,int y){

if(button == GLUT_LEFT_BUTTON && state == GLUT_DOWN)//鼠标左键按下

{

startPoint.x=x;startPoint.y=y;

ButtonDown=LeftButtonDown;

}

if(button==GLUT_MIDDLE_BUTTON &&state==GLUT_DOWN){

startPoint.x=x;startPoint.y=y;

ButtonDown=MiddleButtonDown;

}

//6.2鼠标移动事件

void myMotion(int x,int y){

endPoint.x=x;endPoint.y=y;

int dx=endPoint.x-startPoint.x,dy=endPoint.y-startPoint.y;

switch (ButtonDown)

{

case LeftButtonDown:

{

double k1,k2;//各自所占比例

if(dx==0&&dy==0) return;

k1=(double)dx/(double)(abs(dx)+abs(dy));

k2=(double)dy/(double)(abs(dx)+abs(dy));

g_theta+= k1*move_angle*PI / 180.0;

g_phai-= k2*move_angle*PI / 180.0;

glutPostRedisplay();//刷新

}

break;

case MiddleButtonDown:

{

g_r+=dy;

glutPostRedisplay();//刷新

}

break;

default:

break;

}

startPoint.x=endPoint.x;startPoint.y=endPoint.y;

}

3.7.3效果展示

图26 键盘左箭头移动后的场景

图27 键盘上箭头和左箭头移动后的场景

图28 视口变换后再视角变换后的场景

3.8 地形场景美化之天空盒与地形基座

3.8.1原理

（1）天空盒：天空盒贴图即立方体贴图，按照纹理映射的过程贴图即可，只不过需要注意一下绘制的立法体要完整的包裹整个场景，立方体各个面与纹理坐标的对应关系。这里天空盒子共分为两个场景下——白天和夜晚。为提高效率，可以采用调用显示列表形式，注意天空盒子不受场景光照的效果影像，可以通过设置其材质为自发光白光就行。

（2）地形基座：地形基座即通过绘制地形面周围的5个多边形面即可，除了最底下的面外，其余周围的多边形的构成均是地下一条边，外加该边垂直上升下的高程点。但绘制不能采用GL_POLYGON的形式，不然会出现有些多边形如下图29的效果。

图29 使用GL_POLYGON绘制的错误结果

经过观察与测试，OpenGL内部绘制面GL_POLYGON应该是封装了GL_TRANGLE_FAN的形式，而采用GL_QUAD_STRIP应该是封装了GL_TRANGLE_STRIP的形式。上述两者对顶点的顺序有很高的要求，所以上述几种绘制我均没有得到理想的效果。最后我采取的是GL_QUADS的形式，如图30所示。

图30 使用GL_QUADS绘制折线多边形

3.8.2核心代码

天空盒子(白天)绘制核心代码

//1.3用列表创建天空盒子

void CreateSkyBoxList()

{

SkyBoxObj = glGenLists(1); //2

glNewList(SkyBoxObj, GL_COMPILE);

glPushAttrib(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_CURRENT_BIT);//这里将颜色等属性进行压栈,这句话必须要，是那个画球体的问题，不然进行遥感影像的时候没效果

glEnable(GL_TEXTURE_2D);

glCallList(SunMat);

//第一个面，后面

glBindTexture(GL_TEXTURE_2D, texObj[BACKIMG]);

glBegin(GL_POLYGON);

glNormal3f(0.0,0.0,1.0);

glTexCoord2f(1.0f, 0.0f);glVertex3f(-SkyBoxWidth,-SkyBoxWidth,-SkyBoxWidth);//顶点坐标

glTexCoord2f(0.0f, 0.0); glVertex3f(SkyBoxWidth,-SkyBoxWidth,-SkyBoxWidth);

glTexCoord2f(0,1); glVertex3f(SkyBoxWidth,SkyBoxWidth,-SkyBoxWidth);

glTexCoord2f(1,1); glVertex3f(-SkyBoxWidth,SkyBoxWidth,-SkyBoxWidth);

glEnd();

//第二个面，下面

glBindTexture(GL_TEXTURE_2D, texObj[BOTTOMIMG]);

glBegin(GL_POLYGON);

glNormal3f(0.0,1.0,0.0);

glTexCoord2f(0.0f, 0.0f);glVertex3f(-SkyBoxWidth,-SkyBoxWidth,SkyBoxWidth);//顶点坐标

glTexCoord2f(1.0f, 0.0); glVertex3f(SkyBoxWidth,-SkyBoxWidth,SkyBoxWidth);

glTexCoord2f(1,1); glVertex3f(SkyBoxWidth,-SkyBoxWidth,-SkyBoxWidth);

glTexCoord2f(0,1); glVertex3f(-SkyBoxWidth,-SkyBoxWidth,-SkyBoxWidth);

glEnd();

//第三个面，前面

glBindTexture(GL_TEXTURE_2D, texObj[FRONTIMG]);

glBegin(GL_POLYGON);

glNormal3f(0.0,0.0,-1.0);

glTexCoord2f(0.0f, 0.0f);glVertex3f(-SkyBoxWidth,-SkyBoxWidth,SkyBoxWidth);//顶点坐标

glTexCoord2f(1.0f, 0.0); glVertex3f(SkyBoxWidth,-SkyBoxWidth,SkyBoxWidth);

glTexCoord2f(1,1); glVertex3f(SkyBoxWidth,SkyBoxWidth,SkyBoxWidth);

glTexCoord2f(0,1); glVertex3f(-SkyBoxWidth,SkyBoxWidth,SkyBoxWidth);

glEnd();

//第四个面，左面

glBindTexture(GL_TEXTURE_2D, texObj[LEFTIMG]);

glBegin(GL_POLYGON);

glNormal3f(1.0,0.0,0.0);

glTexCoord2f(1.0f, 0.0f);glVertex3f(-SkyBoxWidth,-SkyBoxWidth,SkyBoxWidth);//顶点坐标

glTexCoord2f(0.0f, 0.0); glVertex3f(-SkyBoxWidth,-SkyBoxWidth,-SkyBoxWidth);

glTexCoord2f(0,1); glVertex3f(-SkyBoxWidth,SkyBoxWidth,-SkyBoxWidth);

glTexCoord2f(1,1); glVertex3f(-SkyBoxWidth,SkyBoxWidth,SkyBoxWidth);

glEnd();

//第五个面，右面

glBindTexture(GL_TEXTURE_2D, texObj[RIGHTIMG]);

glBegin(GL_POLYGON);

glNormal3f(-1,0.0,0.0);

glTexCoord2f(0.0f, 0.0f);glVertex3f(SkyBoxWidth,-SkyBoxWidth,SkyBoxWidth);//顶点坐标

glTexCoord2f(1.0f,0.0); glVertex3f(SkyBoxWidth,-SkyBoxWidth,-SkyBoxWidth);

glTexCoord2f(1,1); glVertex3f(SkyBoxWidth,SkyBoxWidth,-SkyBoxWidth);

glTexCoord2f(0,1); glVertex3f(SkyBoxWidth,SkyBoxWidth,SkyBoxWidth);

glEnd();

//第六个面，上面

glBindTexture(GL_TEXTURE_2D, texObj[TOPIMG]);

glBegin(GL_POLYGON);

glNormal3f(0,-1.0,0.0);

glTexCoord2f(0.0f, 0.0); glVertex3f(SkyBoxWidth,SkyBoxWidth,SkyBoxWidth);

glTexCoord2f(1,0); glVertex3f(SkyBoxWidth,SkyBoxWidth,-SkyBoxWidth);

glTexCoord2f(1,1); glVertex3f(-SkyBoxWidth,SkyBoxWidth,-SkyBoxWidth);

glTexCoord2f(0.0f, 1.0f);glVertex3f(-SkyBoxWidth,SkyBoxWidth,SkyBoxWidth);//顶点坐标

glEnd();

glDisable(GL_TEXTURE_2D);

glPopAttrib();

glEndList();

}

地形基座绘制核心代码

//4.6绘制地形底座

void drawBase(){

switch (Base)

{

case BASE:

{

glPushAttrib(GL_CURRENT_BIT);

glEnable(GL_COLOR_MATERIAL);

glLightModeli( GL_FRONT, GL_AMBIENT_AND_DIFFUSE );

glColor3f(0.5,0.5,0.5);

//共5个面，最低高程-1的目的是解决深度冲突的问题

//绘制底面

glBegin(GL_QUADS);

glVertex3f(MIN_DEMX,MIN_HEIGHT-1,MIN_DEMY);

glVertex3f(MIN_DEMX,MIN_HEIGHT-1,MAX_DEMY);

glVertex3f(MAX_DEMX,MIN_HEIGHT-1,MAX_DEMY);

glVertex3f(MAX_DEMX,MIN_HEIGHT-1,MIN_DEMY);

glEnd();

//绘制左面,第i~0行，第0列的点

glBegin(GL_QUADS);//注意逆时针方向

for(int i=0;i<ROW-1;i++){

glNormal3f(-1,0,0);

glVertex3f(terri[i*COL].x,terri[i*COL].h,terri[i*COL].z);

glVertex3f(terri[i*COL].x,MIN_HEIGHT-1,terri[i*COL].z);

glVertex3f(terri[(i+1)*COL].x,MIN_HEIGHT-1,terri[(i+1)*COL].z);

glVertex3f(terri[(i+1)*COL].x,terri[(i+1)*COL].h,terri[(i+1)*COL].z);

}

glEnd();

//绘制前面,第ROW-1行，第COL-1~0列的点

glBegin(GL_QUADS);//注意逆时针方向

for(int i=0;i<COL-1;i++){

glNormal3f(0,0,1);

glVertex3f(terri[(ROW-1)*COL+i].x,terri[(ROW-1)*COL+i].h,terri[(ROW-1)*COL+i].z);

glVertex3f(terri[(ROW-1)*COL+i].x,MIN_HEIGHT-1,terri[(ROW-1)*COL+i].z);

glVertex3f(terri[(ROW-1)*COL+i+1].x,MIN_HEIGHT-1,terri[(ROW-1)*COL+i+1].z);

glVertex3f(terri[(ROW-1)*COL+i+1].x,terri[(ROW-1)*COL+i+1].h,terri[(ROW-1)*COL+i+1].z);

}

glEnd();

//绘制右面,第0~i行，第COL-1列的点

glBegin(GL_QUADS);//注意逆时针方向

for(int i=0;i<ROW-1;i++){

glNormal3f(1,0,0);

glVertex3f(terri[i*COL+COL-1].x,terri[i*COL+COL-1].h,terri[i*COL+COL-1].z);

glVertex3f(terri[i*COL+COL-1].x,MIN_HEIGHT-1,terri[i*COL+COL-1].z);

glVertex3f(terri[(i+1)*COL+COL-1].x,MIN_HEIGHT-1,terri[(i+1)*COL+COL-1].z);

glVertex3f(terri[(i+1)*COL+COL-1].x,terri[(i+1)*COL+COL-1].h,terri[(i+1)*COL+COL-1].z);

}

glEnd();

//绘制后面,第0行，第0~i列的点

glBegin(GL_QUADS);//注意逆时针方向

for(int i=0;i<COL-1;i++){

glNormal3f(0,0,-1);

glVertex3f(terri[i].x,terri[i].h,terri[i].z);

glVertex3f(terri[i].x,MIN_HEIGHT-1,terri[i].z);

glVertex3f(terri[i+1].x,MIN_HEIGHT-1,terri[i+1].z);

glVertex3f(terri[i+1].x,terri[i+1].h,terri[i+1].z);

}

glEnd();

glDisable(GL_COLOR_MATERIAL);

glPopAttrib();

}

default:

break;

}

3.8.3 效果展示

图31 白天天空盒子

图32 带底座的夜晚天空盒子

3.9地形场景图片导出

3.9.1原理

图形截图的过程其实就是将屏幕上的像素读取到内存了，然后将内存中的数据保存到文件。为使保存的文件名具有唯一性，这里改进老师的代码，加入时间参数，以截图的时间命名文件名，主要也是调用<time.h>库函数。

3.9.2核心代码

老师的截图函数这里就不再展示了，这里只展示返回当前时间字符串以作为截图后的文件名。

截图文件命名核心代码

//2.5获取时间字符串，方便截图

string getTime() //2020-09-11 22:00:49 这个只能到秒

{

struct tm t; //tm结构指针

time_t timep;

time (&timep);

char tmp[64];

localtime_s(&t,&timep);

//const struct tm t2=t;

strftime(tmp, sizeof(tmp), "%Y%m%d%H%M%S", &t);//注意存储时不要出现特殊字符

return tmp;

}

3.9.3 效果展示

图33 截图文件

4、实验总结

本次实验在创建地形这个骨架上我感觉难度不是很大，另外其他的有关图形显示骨架代码的撰写上面也还勉强，自己在编写代码时主要是受光照、纹理混合作用效果的影像，感觉有点难度，很多时候得排查修改蛮久才能得到自己想要的效果，属性堆栈还有有点作用的——那是我卡的最久的一次。不过整体做下来，还是很有意思的。

色带的变化函数自己写的比较简陋，后面可以再研究一下，不过感觉很多时候自己都不用去研究色带的变化函数，就像arcgis二次开发一样，直接调用色带API就行了。另外分级色彩本次只实现了最简单的等距分级色彩，其实还有几何间断点分级、自然间断点分级等等。

本次三维地形的建模，自己4月底做过初步后由于其他事情时间原因后面没有动了，再次做已经是6月中旬了，第二次开工时自己也尝试查阅了很多文献，发现2000年左右的前人的思路很多跟我差不多，还是挺高兴，也借鉴了他们的一些思路。他们提到的多级分辨率模型[4]在渲染场景不同深度的地形时具有加速渲染的效果，其实跟我之前想的大数据如何去进行图形渲染的思想差不多，类似栅格金字塔和矢量切片的思想。

最后本实验已经录制成视频上传至B站，视频地址：三维GIS作业——基于OpenGL的地形模拟_哔哩哔哩_bilibili