LOD层次细节算法-大规模实时地形的绘制

背景

在介绍层次细节算法之前，先来看两幅图片。

图一

图二

这两幅图片是用层次细节算法也即LOD算法绘制的地形网格。为了更清晰的看清地形网格的结构，我没有给其贴上纹理。这两幅图片看上去给人第一感觉就是分辨率不同，图一分辨率较低，图二分辨率很高。图一图二是由同一个程序生成的，图一是在调节系数为1的情况下生成的，图二是在调节系数为25的情况下生成的。为了增加对比度，我故意把两幅图片的分辨率调节的差别很大。这个地形网格如果达到全分辨率的话将会是513像素*513像素。然而读者看到的图片并没有达到全分辨率。为什么呢？大家想一下，在现实世界中，人眼的视角是有一个范围的并不能看到360度范围的场景；随着视线的往远处移动，看到的东西越来越模糊；大家在想一下另外一个问题，把图片贴到一个物体上，如果这个物体的表面是平的，那么任务肯定非常容易完成，如果物体表面凸凹不平，那么你就不能把图片直接贴在上面，你必须把图片剪成小的片段然后再艰难的贴在物体表面。同理层次细节算法就是模拟上述现实场景的技术，只不过我们将换一套专有名词来描述。上述地形网格在满足下面三个条件时候被绘制，一，不在照相机视景体内的网格部分将不会被绘制；二，距离相机视点远的地方网格以低分辨率来绘制，近的地方以高分辨率来绘制；三，粗糙的部分以高分辨率绘制，平坦的部分以较低的分辨率来绘制。这三个条件合称节点评价系统。

地形高程图

在层次细节地形绘制的过程中，每一个顶点坐标(x,y,z)被分为两部分处理,(x,z)与y两个部分。至于为什么这样做，当你看完本篇文章之后你就会明白。在OpenGL里面Y轴是垂直向上的，因此顶点的y坐标值代表顶点的高度。高程图就是存储y坐标值的，其中一种方法就是使用raw格式的文件。raw格式文件是8位的，也就是说把raw格式的图片看成width*hight大小的矩形的话，其中每个元素表示一个8位的数据，范围是0到256，而这个范围在乘以相应的比例正好可以表示高度值。在层次细节算法中要求地形的大小必须是正方形，而且必须满足边长的像素数为(2的n次方)+1;后文大家会明白为什么会有这两个限制。本文我们使用.raw格式的高程图，其制作方法有很多，本文介绍一种最简单的方法，用photoshop生成，如果要想生成自己想要的那种地形需要使用专业的方法来生成高程图，本文使用的是随机生成的高程图。方法很简单，打开PS快捷键crtl+N新建一个513*513像素大小的项目，确定后从工具栏里选择滤镜->渲染->分层云彩,然后存储为选择.raw格式，至此高程图就完成了。下面给出加载高程图的程序代码：

 

01.voidGLLod::loadRawFile(LPSTRstrName, intnSize) 

02.{ 

03.   FILE*pFile=NULL; 

04.    pFile=fopen(strName,"rb"); 

05.   if(pFile==NULL) 

06.    { 

07.        MessageBox(NULL,TEXT("不能打开高度图文件"),TEXT("错误"),MB_OK); 

08.       return; 

09.    } 

10.   fread(pHeightMap,1,nSize,pFile); 

11.   intresult=ferror(pFile); 

12.   if(result) 

13.    { 

14.        MessageBox(NULL,TEXT("读取数据失败"),TEXT("错误"),MB_OK); 

15.    } 

16.   fclose(pFile); 

17.}

从高程图中得到坐标(x,z)处的高度值的程序代码：

?

01.intGLLod::getRawHeight(intx, intz) 

02.{ 

03.   intxx=x%(map_size+1); 

04.   intzz=z%(map_size+1); 

05.   if(!pHeightMap)return0; 

06.   returnpHeightMap[xx+(zz*(map_size+1))]; 

07.}

层次细节算法(LOD算法) 
现在要步入正题了，我将会向大家介绍层次细节算法的原理。LOD算法采用的是四叉树的结构来处理的。

我们看上图，网格一是最初始的没有分割的正方形，其边长是如果以像素表示边长的话那么边上有个像素(假设每两个最近的像素间的距离为1)。那么在这里我们可以看出每一个网格的边长是5个像素，若每个像素间的距离为一的话，那么每个网格的边长是4，下文中我们默认每个像素间的距离是1。将网格一均分为4个小网格，那么这四个小网格是原网格的四个儿子节点，继续划分，图二中的每个儿子网格继续划分为四个儿子。我们可以有选择的划分某些网格，每个网格一旦划分的话，就必需划分为四个相等的小网格。这种划分方法很明显是符合四叉树的，只不过这个四叉树每个节点要么有四个儿子，要么没有儿子。上面描述的网格是在x-z平面上，至此还没有考虑y坐标呢。等划分到一定标准的时候就会在此基础上考虑y坐标，这个时候就可以渲染了，不过距离这一步还有很远的距离。这个划分结束的标准就是前文中所说的节点评价系统(如果满足节点评价系统，但是网格已经达到最大分辨率的话也不能再划分了)。下面我们就来详细讲述节点评价系统。

节点评价系统

相机裁剪  上面我们划分了网格一，直至网格三，但是那个划分是盲目的，接下来我们就必需按照一定的标准有目的的划分，来看几副图片。

图一图二

图三

图中的红色线代表照相机(也可以理解为人的眼睛)的视野范围，即使只有一小部分在相机视野内，我们就需要对其进行划分，显然图一的正方形在视野内，对其划分得到图二中的网格，这个时候可以发现图二中的最右边的两个儿子网格仍然在视野内，因此对其继续划分，最左边两个儿子网格不在视野内，因此不必划分也即被裁剪掉了，于是得到了图三中的网格。现在呈现给读者的是平面的，实际上考虑y轴坐标的话，网格节点是一个三维的，为了方便我们以平面的方式来阐述，实际上是三维的。那么如何进行三维的裁剪呢？在笔者的前一篇文章《3D坐标系、矩阵运算、视景体与裁剪》中对其算法进行了描述，这里只给出代码。

?

01.intGLFrustum::isAabbInFrustum( AABB &aabb) 

02.{ 

03.       //aabb是一个AABB包围盒      

04.    calculateFrustumPlanes();//计算平截头体的六个面的方程  

05.   boolinsect=false;//相机裁剪的标志  

06.   for(inti=0;i<6;i++) 

07.    { 

08.               //接下来3个if语句是轴分离的方法调整aabb包围盒  

09.       if(g_frustumPlanes[i][0]<0.0f) 

10.        { 

11.           inttemp=aabb.min[0]; 

12.            aabb.min[0]=aabb.max[0]; 

13.            aabb.max[0]=temp; 

14.        } 

15.       if(g_frustumPlanes[i][1]<0.0f) 

16.        { 

17.           inttemp=aabb.min[1]; 

18.            aabb.min[1]=aabb.max[1]; 

19.            aabb.max[1]=temp; 

20.        } 

21.       if(g_frustumPlanes[i][2]<0) 

22.        { 

23.           inttemp=aabb.min[2]; 

24.            aabb.min[2]=aabb.max[2]; 

25.            aabb.max[2]=temp; 

26.        } 

27.         

28.       if((g_frustumPlanes[i][0]*aabb.min[0]+g_frustumPlanes[i][1]*aabb.min[1]+g_frustumPlanes[i][2]*aabb.min[2]+g_frustumPlanes[i][3])>0.0f) 

29.        { 

30.             

31.           return0;//不可见  

32.        } 

33. 

34.       if((g_frustumPlanes[i][0]*aabb.max[0]+g_frustumPlanes[i][1]*aabb.max[1]+g_frustumPlanes[i][2]*aabb.max[2]+g_frustumPlanes[i][3])>=0.0f) 

35.        { 

36.            insect=true;//裁剪  

37.        } 

38.             

39.    } 

40.   if(insect)return1;//裁剪  

41.   return2;//完全可见  

42.}

这其中的calculateFrustumPlanes();是用来计算平截头体的六个面的方程读者可以点击这里了解详情。

视点距离上面的相机裁剪部分说在视野范围内的继续划分，那么视野内的那部分是不是一直继续划分呢？当然不是，过度的划分并不会带来视觉上的改观而会增加GPU的处理负担，试想如果有1000万个顶点的话，那处理起来是非常消耗GPU的。根据常识，在现实世界中我们看远处的物体会感觉到模糊，看近处的会感觉比较清晰，同样LOD算法里也是采用这个理念，距离视点远的网格节点就不必要继续划分，而近处的网格则要继续划分。

如上图，一个小兔子的眼睛看着右上角的儿子网格，我们定义距离L是眼睛到网格中心的距离，d是目标网格的边长。满足条件的时候网格继续划分，否则不需要继续划分。C1是一个可以根据实际渲染情况调节的因子。图一 图二就是不同调节因子下所形成的两幅地形网格。

粗糙程度

在物体粗糙的部分需要继续划分以更好的显示，而平坦的部分则不需要做过多的划分，比如一个平面直接贴纹理就行了，做过多的划分是没有意义的。

那么怎么定义网格的粗糙程度呢？如上图所示，如果在xz平面考虑问题的话，那么网格的每个节点的边都是在一个平面上的，如果考虑y坐标的话原来是直线的边会变成曲线，我们以每个网格四条边的起伏程度和中心点的起伏程度的最大值来定义网格的粗糙度。举个例子在上图中dh4的值是那条边的两个顶点的y坐标值相加再除以2以后减去边的中点的y坐标值得到dh4,这个dh4就是所在边的起伏度，同样方法计算dh1 dh2 dh3。中心点需要计算两次，因为中心点所在的边有两条，分别是对应的两条对角线。这6个值计算出来后选择其中最大的作为粗糙度。设粗糙度为DHmax,那么满足条件时继续划分，否则不划分。其中C2是可以调节的因子。这个条件可以和上一个条件合并得到。

总结来说if(相机裁剪通过&&视点-粗糙值合适&&没有达到最大分辨率) then 划分节点，否则不划分。

消除裂缝如果仅仅做到上面所说的是不是就可以渲染地形了呢？答案是否定的，因为这样会产生裂缝。什么是裂缝呢？我们来看两幅图片就可以知道了。

地形一

地形二

为了方便大家观察，每一副图片部分网格以线框的模式渲染另一部分以填充的方式渲染。图一是正常的，图二却有许多列缝。这是什么原因呢？为了弄明白这一点我们先放一放，看另外一个问题，网格是怎么渲染的呢？

如图是一个即将送往3D API渲染的网格节点。共有九个顶点，中心点是0点，另外还有8个点。在OpenGL里这个网格将会以三角形扇的方式进行渲染，比如032是一个三角形，021又是一个三角形。但是这样渲染是有问题的，看接下来的图片：

在上面这个图中左右两个网格的划分层次相差为1，也就是右侧的网格比左侧的多划分一次。这时候问题出现了，当以顶点1 2 3渲染三角形和以顶点2 4 5渲染三角形以及以顶点5 4 3渲染三角形时候就出现了裂缝。图上面是在xz平面上看不出问题，但是当给顶点赋予y坐标值的时候问题就出现了，因为点4可能和点2 点3不在一个高度，所以点2 点4点3组成的可能是一条折线。假如点2 点3的高度相同，点4比点2 点3 高，那么点2 点3 点4便组成了一个三角形，这个三角形就是地形二上面的裂缝。大家可以看下比较大的裂缝会发现正好是个三角形，这就是很多诸如点2 点3 点4组成的三角形裂缝。那么怎么解决这个问题呢？可以在点1 点4之间增加一条线段，这个处理起来比较麻烦，本文采取的是将点4点5组成的线段取消，也即删除点4，这个时候裂缝便会消失，如地形一那样。但是我们忽略了一个问题，这个还是比较棘手的问题，我们再看一张图片：

在这幅图里面，右侧的网格比左侧的网格多划分了2次，这个时候顶点3 顶点6 顶点5 顶点4组成的边比先前那个例子更加复杂了。即使删除点5，点3 点6 点4仍然可能组成一个三角形裂缝。点6不能删除，因为点6是正方形网格的顶点，删除它就等于删除网格了，所以只能删除边的中点。这下怎么办呢？如果有一种方法保证左右两侧的两个正方形网格划分层次相差小于等于1，那么就可以按照前文所说的那样通过删除点来消除三角形裂缝。能做到这一点吗？答案是肯定的。假设左侧的网格划分值为f2，右侧的网格的父亲节点划分值为f1，当f1网格需要划分的时候必须保证f2也划分，如何做到这一点呢？通过观察如果满足表达式;f2<f1的情况下f2会在f1划分时也被划分，因为f1如果被划分的话说明f1<1，而f2<f1从而f2<1也即f1满足节点评价系统，所以也会被划分。也即

因为d2=2d1所以化简以后得到

现在问题集中在L2和L1的比值。看图：

从视点做垂直于xz平面的直线交予xz面与O点。此时有代入后得到

也就是说DHmax2>DHmax1的话就可以满足上式，从而消除三角形裂缝。如何使得DHmax2>DHmax1呢？可以这样做，对于左边的正方形网格求紧贴其四条边的边长为其一半的8个小正方形的DHmax值，再与其自身的DHmax比较，并将最大值当初该网格的DHmax值。这样就可以保证右侧的小正方形父亲被划分时候左侧的大的正方形网格也被划分，这样就可以保证他们的划分层次相差小于等于1。如此便不会出现三角形裂缝。调整DHmax的函数代码是：

?

01.voidGLLod::modifyDHMatrix() 

02.{ 

03.   intedgeLength=2; 

04.   while(edgeLength<=map_size) 

05.    { 

06.       inthalfEdgeLength=edgeLength>>1; 

07.       inthalfChildEdgeLength=edgeLength>>2; 

08.       for(intz=halfEdgeLength;z<map_size;z+=edgeLength) 

09.             

10.        { 

11.           for(intx=halfEdgeLength;x<map_size;x+=edgeLength) 

12.             

13.           if(edgeLength==2) 

14.            { 

15.               intDH6[6]; 

16.                DH6[0]=abs(((getRawHeight(x-halfEdgeLength,z+halfEdgeLength)+getRawHeight(x+halfEdgeLength,z+halfEdgeLength))>>1)-getRawHeight(x,z+halfEdgeLength)); 

17.                DH6[1]=abs(((getRawHeight(x+halfEdgeLength,z+halfEdgeLength)+getRawHeight(x+halfEdgeLength,z-halfEdgeLength))>>1)-getRawHeight(x+halfEdgeLength,z)); 

18.                DH6[2]=abs(((getRawHeight(x-halfEdgeLength,z-halfEdgeLength)+getRawHeight(x+halfEdgeLength,z-halfEdgeLength))>>1)-getRawHeight(x,z-halfEdgeLength)); 

19.                DH6[3]=abs(((getRawHeight(x-halfEdgeLength,z+halfEdgeLength)+getRawHeight(x-halfEdgeLength,z-halfEdgeLength))>>1)-getRawHeight(x-halfEdgeLength,z)); 

20.                DH6[4]=abs(((getRawHeight(x-halfEdgeLength,z-halfEdgeLength)+getRawHeight(x+halfEdgeLength,z+halfEdgeLength))>>1)-getRawHeight(x,z)); 

21.                DH6[5]=abs(((getRawHeight(x+halfEdgeLength,z-halfEdgeLength)+getRawHeight(x-halfEdgeLength,z+halfEdgeLength))>>1)-getRawHeight(x,z)); 

22.               intDHMax=DH6[0]; 

23.               for(inti=1;i<6;i++) 

24.                { 

25.                   if(DHMax<DH6[i]) 

26.                        DHMax=DH6[i]; 

27.                     

28.                } 

29. 

30.                setDHMatrix(x,z,DHMax); 

31.            } 

32.           else 

33.            { 

34.               intDH14[14]; 

35.               intnumDH=0; 

36. 

37.               intneighborX; 

38.               intneighborZ; 

39.                neighborX=x-edgeLength; 

40.                neighborZ=z; 

41.               if(neighborX>0) 

42.                { 

43.                    DH14[numDH]=getDHMatrix(neighborX+halfChildEdgeLength,neighborZ-halfChildEdgeLength); 

44.                    numDH++; 

45.                    DH14[numDH]=getDHMatrix(neighborX+halfChildEdgeLength,neighborZ+halfChildEdgeLength); 

46.                    numDH++; 

47.                } 

48.                neighborX=x; 

49.                neighborZ=z-edgeLength; 

50.               if(neighborZ>0) 

51.                { 

52.                    DH14[numDH]=getDHMatrix(neighborX-halfChildEdgeLength,neighborZ+halfChildEdgeLength); 

53.                    numDH++; 

54.                    DH14[numDH]=getDHMatrix(neighborX+halfChildEdgeLength,neighborZ+halfChildEdgeLength); 

55.                    numDH++; 

56.                } 

57.                neighborX=x+edgeLength; 

58.                neighborZ=z; 

59.               if(neighborX<map_size) 

60.                { 

61.                    DH14[numDH]=getDHMatrix(neighborX-halfChildEdgeLength,neighborZ-halfChildEdgeLength); 

62.                    numDH++; 

63.                    DH14[numDH]=getDHMatrix(neighborX-halfChildEdgeLength,neighborZ+halfChildEdgeLength); 

64.                    numDH++; 

65.                } 

66.                neighborX=x; 

67.                neighborZ=z+edgeLength; 

68.               if(neighborZ<map_size) 

69.                { 

70.                    DH14[numDH]=getDHMatrix(neighborX-halfChildEdgeLength,neighborZ-halfChildEdgeLength); 

71.                    numDH++; 

72.                    DH14[numDH]=getDHMatrix(neighborX+halfChildEdgeLength,neighborZ-halfChildEdgeLength); 

73.                    numDH++; 

74.                } 

75.                DH14[numDH]=abs(((getRawHeight(x-halfEdgeLength,z+halfEdgeLength)+getRawHeight(x+halfEdgeLength,z+halfEdgeLength))>>1)-getRawHeight(x,z+halfEdgeLength)); 

76.                numDH++; 

77.                DH14[numDH]=abs(((getRawHeight(x+halfEdgeLength,z+halfEdgeLength)+getRawHeight(x+halfEdgeLength,z-halfEdgeLength))>>1)-getRawHeight(x+halfEdgeLength,z)); 

78.                numDH++; 

79.                DH14[numDH]=abs(((getRawHeight(x-halfEdgeLength,z-halfEdgeLength)+getRawHeight(x+halfEdgeLength,z-halfEdgeLength))>>1)-getRawHeight(x,z-halfEdgeLength)); 

80.                numDH++; 

81.                DH14[numDH]=abs(((getRawHeight(x-halfEdgeLength,z+halfEdgeLength)+getRawHeight(x-halfEdgeLength,z-halfEdgeLength))>>1)-getRawHeight(x-halfEdgeLength,z)); 

82.                numDH++; 

83.                DH14[numDH]=abs(((getRawHeight(x-halfEdgeLength,z-halfEdgeLength)+getRawHeight(x+halfEdgeLength,z+halfEdgeLength))>>1)-getRawHeight(x,z)); 

84.                numDH++; 

85.                DH14[numDH]=abs(((getRawHeight(x+halfEdgeLength,z-halfEdgeLength)+getRawHeight(x-halfEdgeLength,z+halfEdgeLength))>>1)-getRawHeight(x,z)); 

86.                numDH++; 

87.               intDHMax=DH14[0]; 

88.               for(inti=1;i<14;i++) 

89.                { 

90.                   if(DHMax<DH14[i]) 

91.                        DHMax=DH14[i]; 

92.                } 

93.                setDHMatrix(x,z,DHMax); 

94.         

95.                 

96.            } 

97.        } 

98.        edgeLength=edgeLength<<1; 

99.         

100.    } 

101.}

生成地形网格的代码如下：

?

01.voidGLLod::updateQuadTreeNode(intcenterX, intcenterZ,intedgeLength,intchild) 

02.{ 

03.     

04.   if(edgeLength>2) 

05.    { 

06.         

07.       if(isObjectCulled(centerX,centerZ,edgeLength)) 

08.        { 

09.             

10.            quadNode[centerX+centerZ*(map_size+1)].blend=2; 

11. 

12.        } 

13.       else 

14.        { 

15.             

16.                 

17.             

18.           floatfViewDistance,f; 

19.           inthalfChildEdgeLength; 

20.           intchildEdgeLength; 

21.           intblend; 

22.           intcenterQuad[3]; 

23.            centerQuad[0]=centerX; 

24.            centerQuad[2]=centerZ; 

25.            centerQuad[1]=getRawHeight(centerX,centerZ); 

26.            fViewDistance=frustum.distanceOfTwoPoints(centerQuad); 

27.             

28.            f=fViewDistance/(edgeLength*mfMinResolution*(max(mfDetailLevel*getDHMatrix(centerX,centerZ),1.0f))); 

29.         

30.           if(f<1.0f) 

31.                blend=1; 

32.           else 

33.                blend=0; 

34.             

35.           inttemp=centerX+centerZ*(map_size+1); 

36.            quadNode[temp].blend=blend; 

37.            quadNode[temp].centerX=centerX; 

38.            quadNode[temp].centerY=centerQuad[1]; 

39.            quadNode[temp].centerZ=centerZ; 

40.            quadNode[temp].child=child; 

41.            quadNode[temp].edgeLength=edgeLength; 

42.           if(blend==1) 

43.            { 

44.               inthalfChildEdgeLength=edgeLength>>2; 

45.               intchildEdgeLength=edgeLength>>1; 

46.             

47.                updateQuadTreeNode(centerX-halfChildEdgeLength,centerZ-halfChildEdgeLength,childEdgeLength,1); 

48.                updateQuadTreeNode(centerX+halfChildEdgeLength,centerZ-halfChildEdgeLength,childEdgeLength,2); 

49.                updateQuadTreeNode(centerX-halfChildEdgeLength,centerZ+halfChildEdgeLength,childEdgeLength,3); 

50.                updateQuadTreeNode(centerX+halfChildEdgeLength,centerZ+halfChildEdgeLength,childEdgeLength,4); 

51.            } 

52.             

53.        } 

54.    } 

55.}

原文来自：http://my.oschina.net/u/997451/blog/122649