• entries
232
1463
• views
960211

# Horizon culling

2039 views

Status:

First of all, a little status check. I know I haven't posted a lot of dev journals in the past 2 months, but there's a good reason for that. And it's not because I've been inactive, quite the contrary. I easily spent 80+ hours of work ( split 50/50 between normal job and Infinity work ) in the past 3 weeks.

I feel extremely tired.

Other than the main reason, work, I must also say ( even if I've already complained about it before ) that the pressure is increasing. By pressure, I mean community / management pressure, of course. Private messages, people disturbing me on IRC to ask simple questions ( killing my train-of-thought - if you're no programmer, don't try to understand ), e-mails, requests to get informations on tons of topics..

I think I will soon arrive to a point where, to be able to continue to get actual work done, I'll have to stop coming on IRC or answering emails/PMs at all.

Anyway, back on topic:

I had to re-prioritize a bit the development due to "business opportunities". Nothing that I can announce yet ( if it even gets done ) though. In any case, other than the priorities shift, it shouldn't affect Infinity or its development in any way.

In the past 2 months, I've mostly finished the development of the terrain engine. It's coming along very well. I'm also integrating a lot of older prototypes into the game client. Not particularly hard, but it can be quite long, and a lot of shaders had to be rewritten in GLSL, and tested..

I've added a lot of new features, but I will post more detailed journals about each of them in the coming days.

For now, I just want to make a quick journal about horizon culling.

Horizon culling

If you're no programmer/mathematician, you can stop reading here: you will not understand.

I've been using HC for a while, but the algorithm I've been using wasn't mathematically correct; it was also computationally expensive.

The goal of HC is to hide terrain patches that are hidden by the curvature of the planet's surface.

In maths term, the idea is to find whether a sphere ( approximating the object / terrain patch ) is totally inside the "shadow cone" of another sphere ( the whole planet ) from a camera viewpoint.

Surprisingly, it looks like it's very hard to find such an algorithm on the web. So with my coworker Inigo Quilez, who is better at maths than me, we came up with the following equations. First, a diagram showing the setup in 2D:

- the red circle represents the first sphere ( planet ), at an origin O with a radius R.
- the blue frustum represents the shadow cone. The camera viewpoint is called C and is located at the apex. This cone is tangent to the red sphere. In 3D, it forms an ellipsis, but in 2D, I've represented it as a point Q.
- the yellow sphere represents the object, at an origin O' with a radius R'.

Some other definitions:
- D is the distance between the camera C and the red sphere's origin O;
- D' is the distance between C and the yellow sphere's origin O'.
- T is the distance between C and Q ( Q being a tangent point on the red sphere ).

In order to determine whether the object ( yellow sphere ) is occluded by the red sphere ( planet ), there are two conditions:
- it must be in the shadow cone ( blue frustum ) AND
- it must be behind the capping plane drawn in brown in the diagram, defined by the plane at point P with normal N.

Those are two different problems with their own equations. Let's have a look at the first one.

Cone test

- alpha is the angle for the planet's cone.
- beta is the angle for the object's cone.
- gamma is the angle between the camera and the two sphere's origins.

Calculating gamma is easy:
- normalize OC = |OC|
- normalize O'C = |O'C|
- perform a dot product between the two to get the cosine:
cos(gamma) = dot(|OC|,|O'C|)

The object is visible iif gamma < alpha - beta

Unfortunately, you should notice that we need to subtract two angles. We cannot do that directly with dot-products ( cosine of angles ), but we can use a little transformation of the original formula, as cos(a + b)=cos(a).cos(b) - sin(a). sin(b):

cos(gamma) < cos(alpha).cos(beta) + sin(alpha).sin(beta)

From now on we'll refer to cos(gamma) as K.

Next step: let's find the values of cos(alpha), cos(beta), sin(alpha) and sin(beta):

In triangle COQ we get:
cos(alpha) = T / D
Using pythagoras theorem on COQ we also get:
D.D = T.T + R.R
T = sqrt(D.D - R.R)

so
cos(alpha) = sqrt(D.D - R.D) / D

In the same way,
cos(beta) = sqrt(D'.D' - R'.R') / D'

For the sine it's even more easy: from triangle COQ again we get:
sin(alpha) = R / D
and
sin(beta) = R' / D'

Are we done ? Not yet. Notice that there are some square roots and slow divisions. Those can be optimized. The equation now is:

cos(gamma) < cos(alpha).cos(beta) + sin(alpha).sin(beta)
K < (sqrt(D.D - R.D) / D).(sqrt(D'.D' - R'.R') / D') + (R.R')/(D.D')

D.D' can be put in common and moved to the left of the equation:
K.D.D' < sqrt(D.D - R.R).sqrt(D'.D' - R'.R') + R.R'
K.D.D' - R.R' < sqrt(D.D - R.R).sqrt(D'.D' - R'.R')

We square both sides:
K.K.D.D.D'.D' - 2.K.D.D'.R.R' + R.R.R'.R' < (D.D - R.R).(D'.D' - R'.R')
NOTE: we are only allowing to do this when the left term is positive. Fortunately, the right term is always positive, so if the left term is negative, the result of the equation is TRUE and we can stop here.

We develop the right side:
K.K.D.D.D'.D' - 2.K.D.D'.R.R' + R.R.R'.R' < D.D.D'.D' - R.R.D'.D' - D.D.R'.R' + R.R.R'.R'

After simplifying and re-arranging:
-2.K.R.R'.D.D' + R.R.D'.D' + R'.R'.D.D < (1 - K.K).D.D.D'.D'

Dividing by D.D.D'.D' twice gives:
-2.K.(R / D).(R' / D') + (R / D).(R / D) + (R' / D').(R' / D') < 1 - K.K

The fantastic thing is that if you set K1 = R / D and K2 = R' / D', you get an equation without any square root nor division:
-2.K.K1.K2 + K1.K1 + K2.K2 < 1 - K.K

Finding P and N for the clipping plane:

N is easy: it equals -|CO|

To create the clipping plane, we also need a point on this plane, which we called P:
P = O + N.y

Let's find the value of y:

In triangle CQP we have:
T.T = x.x + h.h ( Pythagoras theorem )
In triangle POQ:
R.R = h.h + y.y

So h.h = R.R - y.y
We replace in the first equation:
T.T = x.x + R.R - y.y
But we also know that T.T = D.D - R.R and x = D - y (by definition)
So
D.D - R.R = (D - y).(D - y) + R.R - y.y
D.D - 2.R.R = D.D + y.y - 2.D.y - y.y
-2.R.R = -2.D.y
y = R.R / D

Final algorithm:

So many formulas.. the code must be complex, right ? Not really. Let's have a look:

//////		Performs horizon culling with an object's bounding sphere, given a view point.///		This function checks whether the object's sphere is inside the view cone formed///		by the view point and this sphere. The view cone is capped so that the visibility///		is false only in the 'shadow' of this sphere.///		@param view Position of view point in world space///		@param obj Bounding sphere of another object.///		@return true if the object's bounding sphere is visible from the viewpoint, false if the///		sphere is in the shadow cone AND not in front of the capping plane.///TBool SSphere3DF::horizonCulling(const SVec3D& view, const SSphere3DF& obj) const{	SVec3D O1C = m_center - view;	SVec3D O2C = obj.m_center - view;	const TFloat D1 = O1C.getLength();	const TFloat D2 = O2C.getLength();	const TFloat R1 = m_radius;	const TFloat R2 = obj.m_radius;	const TFloat iD1 = 1.0f / D1;	const TFloat iD2 = 1.0f / D2;		O1C *= iD1;	O2C *= iD2;	const TFloat K = O1C.dot(O2C);	const TFloat K1 = R1 * iD1;	const TFloat K2 = R2 * iD2;	TBool status = true;	if ( K > K1 * K2 )	{		status = (-2.0f * K * K1 * K2 + K1 * K1 + K2 * K2 < 1.0f - K * K);	}	TFloat y = R1 * R1 * iD1;	SVec3D P = m_center - y * O1C;	SVec3D N = -O1C;	SPlane plane(N, P);	status = status || (plane.getDistance(obj.m_center) > obj.m_radius);	return status;}

The algorithm has of course been extensively tested and validated, so you can safely use it in your own projects ( if you find a use for it ) :)

Awesome! The end result you achieved is very succinct. It looks rather inexpensive as well, considering the amount of geometry it is capable of culling.

Your formula will undoubtedly have a wide range of uses, and is especially invaluable for other projects simulating space environments.

Quote:
 The algorithm has of course been extensively tested and validated, so you can safely use it in your own projects ( if you find a use for it ) :)

If you only knew how many hours I spent working up that functionality... (and you have shown me why mine didn't work). Thanks majorly for that!

I didn't work myself through all those formulas but when it comes to culling terrain patches, why don't you use some sort of "backface culling".
This is just a thought but lets assume you build your planets through reccursive subdivision of a tetrahedron (which would lead to an icosphere with easy LODing because of the hierarchy) and a terrain patch would be a bunch of triangles you get from the subdivision of a single triangle way up in your hierarchy. Now as mountains are relativly small compared to the size of a planet, this trianglesoup closely resembles the single parent triangle and if this parent triangle is backfacing (add a few degrees to be sure) then the whole terrainpatch is likely to be occluded. The backfacing test would boil down to a single dotproduct of view-direction and parent-normal. If its positive (or greater then some error-threshold) then its backfacing.
In addition you could speed things up by making use of the hierarchy and checking the top levels first. If they are "extremely" backfacing (eg. the backside of the planet) discard all following levels and if they are "slightly" backfacing (eg. parrallel to the view direction) go down a few levels and check those again.
However for other geometry like ships, moons, ... the above equations are still necessery, so as I said, just a thought...

Quote:
 Original post by Ohforf sake I didn't work myself through all those formulas but when it comes to culling terrain patches, why don't you use some sort of "backface culling".

That's a good question. This is the "old algorithm" I mentionned at the beginning of the post. I was doing it that way, but it's really not as easy as you think :)

Quote:
 Original post by Ohforf sake Now as mountains are relativly small compared to the size of a planet,

A dangerous assumption. Remember that in my case, I have a whole planet down to a meter above the ground. This means that a single terrain patch can look extremely tiny compared to a mountain..

Quote:
 Original post by Ohforf sake The backfacing test would boil down to a single dotproduct of view-direction and parent-normal. If its positive (or greater then some error-threshold) then its backfacing.

Of course, this is the first idea that would come to anybody's mind. And that was my first implementation.

The funny thing is, the code I posted above is a generalization of that. Think about it: it checks whether a sphere is inside the "shadow cone" of the planet sphere. And how is a cone defined ? By an apex ( the camera ) and an angle.

So checking if a point is inside the cone comes down to comparing a dot-product, which is exactly what your "backfacing" idea is.

The problem start when you try to take the terrain patch (object)'s size into account. A point-inside-cone isn't good enough, especially for the higher levels of details, as a single terrain patch can be hundreds or thousands of kilometers big.

You can add a small epsilon or play with threshold angles depending on the object size, but you'll always find a case where it doesn't work, and a terrain patch gets culled while it should be visible on screen. The reason for that is that dot-products don't add up linearly, so when you're adding two angles:
angle1 < angle2 + angle3
is not equivalent to doing:
cos(angle1) < cos(angle2) + cos(angle3)

So the basic idea of "correcting" the threshold by a function of the object size doesn't work, because you're working on the cosine and not the angles themselves.

Of course, you can do it on angles, but it requires some expensive acos calls, which my HC code above avoids.

Quote:
 Original post by Ohforf sake In addition you could speed things up by making use of the hierarchy and checking the top levels first.

Yes, that's an obvious optimization. I already implemented it, and fortunately the bounding spheres of the patches of the quadtree are conservative, which means that any child's bounding sphere is inside its parent's, so it works perfectly fine.

Didn't help much for performance though, as the bottleneck wasn't in horizon culling in the first place :)

Whoa, most crucial entry.

I don't know what horizon culling is, but suddenly I want to implement it into my engine lol.

Keep up the good work man, moooreee screeennssss /end zombie voice

Always interesting to see what you are working on, and thanks for sharing the fruits of your hard work.

About being tired, take it easy and don't work too hard. Programming can really burn you out if you take a too heavy load on. Our minds are not designed to be inside a computer all day long, although that would have been cool... ;)

First, your engine screen shots are incredible. Wow.

Some time ago at work, we realized that our horizon culling code was lacking. Your post got me thinking about the problem again. I have a different approach and thought that it might be useful to you.

Consider the following image:

E - earth center
P - terrain patch center
H - horizon
V - viewer position

The circle E of center E is the earth, the occluder. The circle P of center P is the bounding circle of the terrain patch.

Imagine this: if P moved slightly closer to the viewer along the surface of the earth, circle P would be visible over the horizon. Circle P would be below the horizon if P moved away from the viewer. If you knew the distance V to P where circle P is at the edge of being visible, you could could compare that to the actual distance of V to P. If the actual distance is less, circle P is visible.

Here is the main monkey business from a glut program I wrote for testing:

static bool Visible(const sVec2d& viewerPosition, double horizonDistance,
{
const sVec2d pToEVec = {patchPosition.x - earthPosition.x, patchPosition.y - earthPosition.y};
double  temp0 = pToEVec.MagnitudeSqrd() - (temp1 * temp1);
if (temp0 > 0.0)
{
temp0 = sqrt(temp0) + horizonDistance;
const sVec2d vToPVec = {patchPosition.x - gViewerPosition.x, patchPosition.y - gViewerPosition.y};
}
return false;
}



Some notes:

1. The horizon distance is input to this function since it only has to be calculated once and can be used over and over again for each terrain patch. While I framed the problem in terms of a terrain patch, the code works fine for a bounding circle off of circle P. For example, you might use this code to determine if a planet occludes a spacecraft.

2. This should translate over to 3D with minor modifications.

3. I really wish I could get rid of that square root.

4. I do not do a plane test as you do. Circle E should be the minimum size of the planet, i.e., no terrain geometry is inside the circle E. If circle P is ever completely inside the circle E, circle P is not visible. In my code, such a case would cause the function to return false.

5. I can make the glut test program available if there is any interest.

I wrote the algorithm and code last night, so there is a good chance that there are code bugs and deficiencies with my thought process. Anyway, I hope this helps. As Eek! the cat says, "It never hurts to help."

Keep the screen shots coming.

Deron
Point Break

Does anybody have a copy of the diagram in Ysanya's post? I'm trying to understand the code but without the diagram it's too hard.

## Create an account

Register a new account