level of detail algorithmrendering of massive realtime terrain

BackGround:

Before I introduce the algorithm of level of detail，let's look at two pictures below:

The two pictures are grid terrian which is rendered using the level of detail algorithm.In order to display the structure in detail i render them in the mode line without attching texture.The first impression we get from the two pictures is the distinction of the resolution.a lower resolution of picture one and a higher resolution of picture two.The two grid terrian are generated from the identical program.The adjustment coefficent of the first terrian is 1 and 25 of the second one.I deliberately exggerating the diffenerce of the resolution between the two grid terrian in order to increase the contrast.If the grid terrian reach the full resolution ,it will be 513pixel*513pixel.Howerer as you see the two picture above has not reached the full resolution.why?let's think,one person can't get a 360 degree scene in the real world.As the view point moves father ,the scene becomes blur.let's consider anther question about mapping a texture on a plane,it's esay to map a texture on a plat suface,but will be harder to map a texture on a no-plat surface,in that case,we may cut the texture into som samller pieces,and carefully place one piece on some area of the surface and then one by one the same processing.And the same way LOD algorithm works,we just uses different terminoloy to discrible in LOD.The LOD terrian will be rendered if three conditions meets below.One:the section of the terrion will not be processed if it is not in the view area of the camera,Second:the futher part of the terrian won't be rendered in high resolution and the nearer part will be.Third:The rough part of the terrian will be rendered in high resolution and the plat part low resolution.we call the three conditons node evalution system.

Height Map

when we process LOD,we will seperate every vertice into coordinate(x,z) and (y) which y is got just from the Height Map.In OpenGL the y axis is straight up,so the y coordinate represents the height of the vertice.The Height Map just stores the Y which wanted.One of method to create the height map is making a .raw file by photoshop.Every pixel of the raw file within 8bits which range from 0 to 255,we can multiply a proportion when we really draw a LOD terrian.In LOD the terrian must be a square who has a edge length and if we use the number of pixels to reprensent the edge length.a function to load a raw file show as below:

void GLLod::loadRawFile(LPSTR strName, int nSize)
{
	FILE *pFile=NULL;
	pFile=fopen(strName,"rb");
	if(pFile==NULL)
	{
		MessageBox(NULL,TEXT("不能打开高度图文件"),TEXT("错误"),MB_OK);
		return;
	}
	fread(pHeightMap,1,nSize,pFile);
	int result=ferror(pFile);
	if(result)
	{
		MessageBox(NULL,TEXT("读取数据失败"),TEXT("错误"),MB_OK);
	}
	fclose(pFile);
}

we can use the below function to get the Y values from the raw height map:

int GLLod::getRawHeight(int x, int z)
{
	int xx=x%(map_size+1);
	int zz=z%(map_size+1);
	if(!pHeightMap) return 0;
	return pHeightMap[xx+(zz*(map_size+1))];
}

Level Of Detail algorithm

Down to the bussiness of LOD algorithm now,let's begin.LOD uses quad-tree structure ;

we evenly divide grid one into four little grids which in the same size,and we call them the child grid of the origin big grid,contiue we will reach the grid three ,ofcourse we can choose the one we want to divile and levae alone the other grid.Once a grid divided ,it must evenly divided into four same child grid.obviouly this process is just alike Quad-Tree structure.By now we deal the grid in X-Z plane and let out the y values which we should take account of when nesscery.

Node evalution system

we have divide the grid one along the way to grid third,however it is a blind way to divide grid without considering some other conditions which i will tell you int the following section;

图二图三

the two red line represent the camera's viewing range.The grid should be divided as long as even a very little part of the grid within the viewing range .Obvisouly grid one is within the camera's viewing range,so it is divided into grid two,the same way we find that the top-right and bottom-right child grid of the grid two should be divided,the we geg grid thrid.Now what you see is a x-z plane,in the practice we must consider y axis.how can we deal with a 3D scissors?please see my other blog:http://blog.csdn.net/cs_huster/article/details/8794731 I just give the source code :

int GLFrustum::isAabbInFrustum( AABB &aabb)
{
        //aabb是一个AABB包围盒	
	calculateFrustumPlanes();//计算平截头体的六个面的方程
	bool insect=false;//相机裁剪的标志
	for(int i=0;i<6;i++)
	{
                //接下来3个if语句是轴分离的方法调整aabb包围盒
		if(g_frustumPlanes[i][0]<0.0f)
		{
			int temp=aabb.min[0];
			aabb.min[0]=aabb.max[0];
			aabb.max[0]=temp;
		}
		if(g_frustumPlanes[i][1]<0.0f)
		{
			int temp=aabb.min[1];
			aabb.min[1]=aabb.max[1];
			aabb.max[1]=temp;
		}
		if(g_frustumPlanes[i][2]<0)
		{
			int temp=aabb.min[2];
			aabb.min[2]=aabb.max[2];
			aabb.max[2]=temp;
		}
		
		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)
		{
			
			return 0;//不可见
		}

		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)
		{
			insect=true;//裁剪
		}
			
	}
	if(insect) return 1;//裁剪
	return 2;//完全可见
}

calculateFrustumPlanes() is used to caculate the equation of the six planes of the frustum.

Eye distance The part of the grid which is far from our eyes should be rendered in low resolution and the near part the high resolution.let's see a picture;

the rabbit looks at the top right child grid of the big grid.we define l as the distance between the top-right child grid and the ribbit's eyes.when the condition is met we will continue divide the grid into four identical child grid,otherwise we should't do the divide.The value of C1 can be adjust to the appropriate according to the actual situation.picture one and picture is just generated from different value of C1.

The roughness we should render the rough part in high resolution and the plat part low resolution.just as the plat part we can only draw a rectangle mapped with a texture.

How can we define the roughness?Just as the picture above shows if we assume the grid be within the x-z plane,then it is a 2D grid,actually it is a 3D grid,now we must consider the Y axis.We define every edge's fluctuation as (one point's y value+another point's y value )/2-the y value of middle point between them. then we got six value as dh1-dh6 and choose the max as DHmax.Then we use DHmax reprensent the roughness of the grid.If the condition meets we will continue the grid,otherwise won't divide.If we consider this condition and the condition in View range we will get a new conditon which replace the two old condition.and the new condition is

As a summery if(grid in the camera's view&&viewpoint-roughness allows&&not reach the full resolution) we will divide the grid,otherwise won't divide.

Remove the cracks:

can we render the terrian grid if we meet the condition above?the answer is no because of the cracks.

terrian one

terrian two.

Terrian seems normal but terrian get some defect,what causes it?In order to grab this let's look at the picture below to recognise how a grid be rendered.

The grid above will soon deliver to 3D api to render,it has 9 vertices range from 0 to 8.we will render it as a trangle fan.But it will cause crack defect,if we look at another picture we will understand.

as shown above,the devide level between left and right grid is different,right grid has a lever higher than the left one.The problem just arises when we render the two grid.Let's see,if we render a triangle with point 2、4、5 and anthoer triangle with point 3 4 5 and one more triangle with point 1 2 3.Then the problem just be caused by point 1 2 3.So far,we couldn't find the reason,because I just show the picture in a 2D plane,let's take account for a 3D space,in this case,point 2 point 3 and point 4 are not int one line,let's assume that point 4 are higher than ponint 2 and point 3,then the three point(2 4 3)just compose a triangle which we didn't render?DO you remember we just render one triangle with point 1 2 3,another one with point 2 4 5,and the third with point 3 4 5,just without a  triangle  which composed with point 2 4 3,so the crack triange will occur at this time just as shown in "terrian two" above.How can we solve this problem?We can add a segment between point 1 and point 4,also we can delete segment which composed with point 5 and point 4.In this blog I try to solve the crack problem with the latter idea for the purpose of simplified.But we ignore a problem,Let's see another picture.

The picture above shows：the right grid divided two more time than the left one,in other words,the right grid is two lever larger than the left one.In this case a more complex broken line formed with point 3 6 5 and 4,we just cann't solve this problem with deleting  the segment composed with point 5 and 2,becase after ding that there are also three points:4 6 3,these three points still can construct a crack triangle.Can we delete one more segment which composed with point 2 and point 6,ofcourse not,the right grid would disappear shall we do that,since point 6 is a prerequisite part of the right grid.If point 6 disaperas then the grid no longer exists.If we can ensure the right grid divided just one more time than the left one or the same time,we can avoid this bad case.Can we do this? The answer is affirmatory.Let's assum that the divided equation for the left is f2 and the right grid's father grid f1.If f1 is divided we ensure f2 also divided,how can we do this?If we ensure f2<f1,then we can do it.That is

because d2=2d1;then we can get:

Now the problem concentrated in l2 and l1.

As show above,we construct a line perpendicular to the plane d1-d-d2.Then we can get

after substitution we get:

That is to say if we can guarantee DHmax2>DHmax1,then we can achieve our intent.for each grid when we intent to get it's Dhmax,we try to get another 8 smaller grids' Dhmax besides the grid it'self,the we choose the max Dhmax and assinment the max Dhmax to the origin grid.The 8 grids are along the origin grid's four edges,each edge 2 smaller grids.And the code is:

void GLLod::modifyDHMatrix()
{
	int edgeLength=2;
	while(edgeLength<=map_size)
	{
		int halfEdgeLength=edgeLength>>1;
		int halfChildEdgeLength=edgeLength>>2;
		for(int z=halfEdgeLength;z<map_size;z+=edgeLength)
			
		{
			for(int x=halfEdgeLength;x<map_size;x+=edgeLength)
			
			if(edgeLength==2)
			{
				int DH6[6];
				DH6[0]=abs(((getRawHeight(x-halfEdgeLength,z+halfEdgeLength)+getRawHeight(x+halfEdgeLength,z+halfEdgeLength))>>1)-getRawHeight(x,z+halfEdgeLength));
				DH6[1]=abs(((getRawHeight(x+halfEdgeLength,z+halfEdgeLength)+getRawHeight(x+halfEdgeLength,z-halfEdgeLength))>>1)-getRawHeight(x+halfEdgeLength,z));
				DH6[2]=abs(((getRawHeight(x-halfEdgeLength,z-halfEdgeLength)+getRawHeight(x+halfEdgeLength,z-halfEdgeLength))>>1)-getRawHeight(x,z-halfEdgeLength));
				DH6[3]=abs(((getRawHeight(x-halfEdgeLength,z+halfEdgeLength)+getRawHeight(x-halfEdgeLength,z-halfEdgeLength))>>1)-getRawHeight(x-halfEdgeLength,z));
				DH6[4]=abs(((getRawHeight(x-halfEdgeLength,z-halfEdgeLength)+getRawHeight(x+halfEdgeLength,z+halfEdgeLength))>>1)-getRawHeight(x,z));
				DH6[5]=abs(((getRawHeight(x+halfEdgeLength,z-halfEdgeLength)+getRawHeight(x-halfEdgeLength,z+halfEdgeLength))>>1)-getRawHeight(x,z));
				int DHMax=DH6[0];
				for(int i=1;i<6;i++)
				{
					if(DHMax<DH6[i])
						DHMax=DH6[i];
					
				}

				setDHMatrix(x,z,DHMax);
			}
			else
			{
				int DH14[14];
				int numDH=0;

				int neighborX;
				int neighborZ;
				neighborX=x-edgeLength;
				neighborZ=z;
				if(neighborX>0)
				{
					DH14[numDH]=getDHMatrix(neighborX+halfChildEdgeLength,neighborZ-halfChildEdgeLength);
					numDH++;
					DH14[numDH]=getDHMatrix(neighborX+halfChildEdgeLength,neighborZ+halfChildEdgeLength);
					numDH++;
				}
				neighborX=x;
				neighborZ=z-edgeLength;
				if(neighborZ>0)
				{
					DH14[numDH]=getDHMatrix(neighborX-halfChildEdgeLength,neighborZ+halfChildEdgeLength);
					numDH++;
					DH14[numDH]=getDHMatrix(neighborX+halfChildEdgeLength,neighborZ+halfChildEdgeLength);
					numDH++;
				}
				neighborX=x+edgeLength;
				neighborZ=z;
				if(neighborX<map_size)
				{
					DH14[numDH]=getDHMatrix(neighborX-halfChildEdgeLength,neighborZ-halfChildEdgeLength);
					numDH++;
					DH14[numDH]=getDHMatrix(neighborX-halfChildEdgeLength,neighborZ+halfChildEdgeLength);
					numDH++;
				}
				neighborX=x;
				neighborZ=z+edgeLength;
				if(neighborZ<map_size)
				{
					DH14[numDH]=getDHMatrix(neighborX-halfChildEdgeLength,neighborZ-halfChildEdgeLength);
					numDH++;
					DH14[numDH]=getDHMatrix(neighborX+halfChildEdgeLength,neighborZ-halfChildEdgeLength);
					numDH++;
				}
				DH14[numDH]=abs(((getRawHeight(x-halfEdgeLength,z+halfEdgeLength)+getRawHeight(x+halfEdgeLength,z+halfEdgeLength))>>1)-getRawHeight(x,z+halfEdgeLength));
				numDH++;
				DH14[numDH]=abs(((getRawHeight(x+halfEdgeLength,z+halfEdgeLength)+getRawHeight(x+halfEdgeLength,z-halfEdgeLength))>>1)-getRawHeight(x+halfEdgeLength,z));
				numDH++;
				DH14[numDH]=abs(((getRawHeight(x-halfEdgeLength,z-halfEdgeLength)+getRawHeight(x+halfEdgeLength,z-halfEdgeLength))>>1)-getRawHeight(x,z-halfEdgeLength));
				numDH++;
				DH14[numDH]=abs(((getRawHeight(x-halfEdgeLength,z+halfEdgeLength)+getRawHeight(x-halfEdgeLength,z-halfEdgeLength))>>1)-getRawHeight(x-halfEdgeLength,z));
				numDH++;
				DH14[numDH]=abs(((getRawHeight(x-halfEdgeLength,z-halfEdgeLength)+getRawHeight(x+halfEdgeLength,z+halfEdgeLength))>>1)-getRawHeight(x,z));
				numDH++;
				DH14[numDH]=abs(((getRawHeight(x+halfEdgeLength,z-halfEdgeLength)+getRawHeight(x-halfEdgeLength,z+halfEdgeLength))>>1)-getRawHeight(x,z));
				numDH++;
				int DHMax=DH14[0];
				for(int i=1;i<14;i++)
				{
					if(DHMax<DH14[i])
						DHMax=DH14[i];
				}
				setDHMatrix(x,z,DHMax);
		
				
			}
		}
		edgeLength=edgeLength<<1;
		
	}
}

The algorithm for generating Lod Terrian is:

void GLLod::updateQuadTreeNode(int centerX, int centerZ,int edgeLength,int child)
{
	
	if(edgeLength>2)
	{
		
		if(isObjectCulled(centerX,centerZ,edgeLength))
		{
			
			quadNode[centerX+centerZ*(map_size+1)].blend=2;

		}
		else
		{
			
				
			
			float fViewDistance,f;
			int halfChildEdgeLength;
			int childEdgeLength;
			int blend;
			int centerQuad[3];
			centerQuad[0]=centerX;
			centerQuad[2]=centerZ;
			centerQuad[1]=getRawHeight(centerX,centerZ);
			fViewDistance=frustum.distanceOfTwoPoints(centerQuad);
			
			f=fViewDistance/(edgeLength*mfMinResolution*(max(mfDetailLevel*getDHMatrix(centerX,centerZ),1.0f)));
		
			if(f<1.0f)
				blend=1;
			else
				blend=0;
			
			int temp=centerX+centerZ*(map_size+1);
			quadNode[temp].blend=blend;
			quadNode[temp].centerX=centerX;
			quadNode[temp].centerY=centerQuad[1];
			quadNode[temp].centerZ=centerZ;
			quadNode[temp].child=child;
			quadNode[temp].edgeLength=edgeLength;
			if(blend==1)
			{
				int halfChildEdgeLength=edgeLength>>2;
				int childEdgeLength=edgeLength>>1;
			
				updateQuadTreeNode(centerX-halfChildEdgeLength,centerZ-halfChildEdgeLength,childEdgeLength,1);
				updateQuadTreeNode(centerX+halfChildEdgeLength,centerZ-halfChildEdgeLength,childEdgeLength,2);
				updateQuadTreeNode(centerX-halfChildEdgeLength,centerZ+halfChildEdgeLength,childEdgeLength,3);
				updateQuadTreeNode(centerX+halfChildEdgeLength,centerZ+halfChildEdgeLength,childEdgeLength,4);
			}
			
		}
	}
}

If you have find some error in this blog,you can contact me by:microland@126.com