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 ) :)

static bool Visible(const sVec2d& viewerPosition, double horizonDistance,
    const sVec2d& patchPosition, double patchRadius,
    const sVec2d& earthPosition, double earthRadius)
{
    const sVec2d pToEVec = {patchPosition.x - earthPosition.x, patchPosition.y - earthPosition.y};
    const double temp1 = earthRadius - patchRadius;
    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 ((temp0 * temp0) + (patchRadius * patchRadius)) > vToPVec.MagnitudeSqrd();
    }
    return false;
}