How to generate a procedural sphere.

I have received an email asking to explain how I generate a procedural sphere as per this post. Rather than respond to emails personally, I have decide in future to answer the questions on the blog, presenting the information to a wider audience for comment, correction, improvement and perhaps a little learning.

The base of my generated sphere is a cube. Depending on how smooth the generated sphere needs to be, I recursively split each side into 4 equally sized children. This is achieved my finding the center of the face, popping that point out so that it is the correct radius from the center, then making 4 faces using the original points and this new point. (There is some ASCII art showing this in the code below)

void CSphere::Initialize(float p_fRadius, int p_iMaxDepth)
{
        //         TLB----TRB   
        //        /|           /|
        //      /  |         /  |
        //    TLF----TRF    |
        //    |     BLB--|BRB
        //    |   /        |  /
        //    | /          |/
        //    BLF----BRF

        //Putting the Vertices of the initial Cube at p_fRadius is not correct as the distance of
        //p_fRadius, p_fRadius, p_fRadius from the Origin is greater than p_fRadius.

        CVector3 l_Vertices[8];
        l_Vertices[TOP_LEFT_FRONT] = MoveToRadiusDistance(CVector3(-p_fRadius,  p_fRadius, -p_fRadius), p_fRadius);
        l_Vertices[TOP_RIGHT_FRONT] = MoveToRadiusDistance(CVector3( p_fRadius,  p_fRadius, -p_fRadius), p_fRadius);
        l_Vertices[BOTTOM_RIGHT_FRONT] = MoveToRadiusDistance(CVector3( p_fRadius, -p_fRadius, -p_fRadius), p_fRadius);
        l_Vertices[BOTTOM_LEFT_FRONT] = MoveToRadiusDistance(CVector3(-p_fRadius, -p_fRadius, -p_fRadius), p_fRadius);
        l_Vertices[TOP_LEFT_BACK] = MoveToRadiusDistance(CVector3(-p_fRadius,  p_fRadius,  p_fRadius), p_fRadius);
        l_Vertices[TOP_RIGHT_BACK] = MoveToRadiusDistance(CVector3( p_fRadius,  p_fRadius,  p_fRadius), p_fRadius);
        l_Vertices[BOTTOM_RIGHT_BACK] = MoveToRadiusDistance(CVector3( p_fRadius, -p_fRadius,  p_fRadius), p_fRadius);
        l_Vertices[BOTTOM_LEFT_BACK] = MoveToRadiusDistance(CVector3(-p_fRadius, -p_fRadius,  p_fRadius), p_fRadius);

        //Initialize the faces of the cube (The face structure just stores the vertices for the face corners and has render functionality, not applicable to this explanation, and its depth in the face tree)
        m_pFaces[FRONT].Initialise(FRONT, l_Vertices[TOP_LEFT_FRONT], l_Vertices[TOP_RIGHT_FRONT], l_Vertices[BOTTOM_RIGHT_FRONT], l_Vertices[BOTTOM_LEFT_FRONT], DEPTH0);
        m_pFaces[RIGHT].Initialise(RIGHT, l_Vertices[TOP_RIGHT_FRONT], l_Vertices[TOP_RIGHT_BACK], l_Vertices[BOTTOM_RIGHT_BACK], l_Vertices[BOTTOM_RIGHT_FRONT], DEPTH0);
        m_pFaces[BACK].Initialise(BACK, l_Vertices[TOP_RIGHT_BACK], l_Vertices[TOP_LEFT_BACK], l_Vertices[BOTTOM_LEFT_BACK], l_Vertices[BOTTOM_RIGHT_BACK], DEPTH0);
        m_pFaces[LEFT].Initialise(LEFT, l_Vertices[TOP_LEFT_BACK], l_Vertices[TOP_LEFT_FRONT], l_Vertices[BOTTOM_LEFT_FRONT], l_Vertices[BOTTOM_LEFT_BACK], DEPTH0);
        m_pFaces[TOP].Initialise(TOP, l_Vertices[TOP_LEFT_BACK], l_Vertices[TOP_RIGHT_BACK], l_Vertices[TOP_RIGHT_FRONT], l_Vertices[TOP_LEFT_FRONT], DEPTH0);
        m_pFaces[BOTTOM].Initialise(BOTTOM, l_Vertices[BOTTOM_LEFT_FRONT], l_Vertices[BOTTOM_RIGHT_FRONT], l_Vertices[BOTTOM_RIGHT_BACK], l_Vertices[BOTTOM_LEFT_BACK], DEPTH0);

        //Subdivide each patch to the lowest resolution
        m_pPatches[FRONT].SubDivide(p_fRadius, p_iMaxDepth);
        m_pPatches[RIGHT].SubDivide(p_fRadius, p_iMaxDepth);
        m_pPatches[BACK].SubDivide(p_fRadius, p_iMaxDepth);
        m_pPatches[LEFT].SubDivide(p_fRadius, p_iMaxDepth);
        m_pPatches[TOP].SubDivide(p_fRadius, p_iMaxDepth);
        m_pPatches[BOTTOM].SubDivide(p_fRadius, p_iMaxDepth);
}

where

bool CSphereFace::SubDivide(float p_fRadius, int p_iMaxDepth)
{
        if (m_iDepth >= p_iMaxDepth)
            return false;

        //Create the Additional Vertices
        //                                             
        //    NW---------------D-------------NE  
        //    |                      |                  |   
        //    |                      |                  |    
        //    |                      |                  |     
        //    |                      |                  |    
        //    A--------------Center-----------C   
        //    |                      |                  |     
        //    |                      |                  |     
        //    |                      |                  |     
        //    |                      |                  |      
        //    SW----------------B-------------SE   
        //                                                 
        //

        g_vAdditionalVertices[A] = m_vBaseVertices[eNorthWest] + ((m_vBaseVertices[eSouthWest] - m_vBaseVertices[eNorthWest]) / 2.0f);
        g_vAdditionalVertices[B] = m_vBaseVertices[eSouthWest] + ((m_vBaseVertices[eSouthEast] - m_vBaseVertices[eSouthWest]) / 2.0f);
        g_vAdditionalVertices[C] = m_vBaseVertices[eNorthEast] + ((m_vBaseVertices[eSouthEast] - m_vBaseVertices[eNorthEast]) / 2.0f);
        g_vAdditionalVertices[D] = m_vBaseVertices[eNorthWest] + ((m_vBaseVertices[eNorthEast] - m_vBaseVertices[eNorthWest]) / 2.0f);
   
        //Create Child Nodes
        m_pChildren = new CSphereFace[4];
        m_pChildren[eNorthWest].Initialise(eNorthWest, m_vBaseVertices[eNorthWest], g_vAdditionalVertices[D], m_vBaseVertices[eCentre], g_vAdditionalVertices[A], m_iDepth + 1);
        m_pChildren[eNorthEast].Initialise(eNorthEast, g_vAdditionalVertices[D], m_vBaseVertices[eNorthEast], g_vAdditionalVertices[C], m_vBaseVertices[eCentre], m_iDepth + 1);
        m_pChildren[eSouthWest].Initialise(eSouthWest, g_vAdditionalVertices[A], m_vBaseVertices[eCentre], g_vAdditionalVertices[B], m_vBaseVertices[eSouthWest], m_iDepth + 1);
        m_pChildren[eSouthEast].Initialise(eSouthEast, m_vBaseVertices[eCentre], g_vAdditionalVertices[C], m_vBaseVertices[eSouthEast], g_vAdditionalVertices[B], m_iDepth + 1);   

        m_pChildren[eNorthWest].SubDivide(p_fRadius, p_iMaxDepth);
        m_pChildren[eNorthEast].SubDivide(p_fRadius, p_iMaxDepth);
        m_pChildren[eSouthWest].SubDivide(p_fRadius, p_iMaxDepth);
        m_pChildren[eSouthEast].SubDivide(p_fRadius, p_iMaxDepth);

        return true;
}

and

CVector3 MoveToRadiusDistance(CVector3 p_Vector, float p_fRadius)
{
       //Get the Normalized Vector, of this vertex from the origin (center of the sphere) and pop it out to the correct radius,
        return p_Vector.Normalize() * p_fRadius;
}

That is pretty much the basics of how I create the sphere from recursively subdividing a cube. If I overlooked anything please comment and I will correct the post to reflect any gap in information or any mistake.

Can I claim progress even if I am behind where I used be?

I am the first to admit annoyance at having to redo functionality which I have previously implemented, however with the experience and lessons learned when coding procedural planets the first time, my implementation this time is smaller (in code) and more efficient than before.

I shall try and highlight an area which is

1. Easier to implement
2. Less code
3. More efficient

From the centre of LOD (Level of Detail) every patch within a certain distance is rendered at the same LOD. The previous version of Decade calculated the distance from the centre of the patching being tested to the centre of LOD. The length between two positions in 3D space is calculated as

#define SQUARE(x)                ((x)*(x))

float Length(CVector3* p_pvOne, CVector3* p_pvTwo)
{
    return (float)(sqrt(SQUARE(pvTwo->X - p_pvOne->X) + SQUARE(pvTwo->Y - p_pvOne->Y) + SQUARE(pvTwo->Z - p_pvOne->Z)));
}

(The sqrt function has been traditionally considered slow and avoided if possible. I am unsure of its performance on modern hardware however the above could be optimised by not calculating the sqrt in the length function, and comparing it to the expected value squared.)

An better optimisation is to remove the requirement for calculating the distance from each patch to the centre of LOD completely. Instead of using the distance from the centre, I now recursively step from the centre to each neighbour, then onto each of their neighbours while incrementing the 'distance' with each step. When the 'distance' reaches a predefined value, the edge of the LOD area has been reached. In the following images 'distance limit' is set to 3.

The results above are not as desired. Some analysis showed that the east and west neighbours of the centre of LOD are within 3 recursive steps from the north neighbour (which is processed first). Because of this the east and west patches are marked as processed and will not be processed again when directly updated from the centre patch.

To overcome this, when processing a patch, if it is already flagged as processed, I compare its distance in steps from the centre of LOD. If the current distance is less than the stored distance I process again. Reading that sounds a little confusing so I shall try and explain in some steps. (only key steps are listed)

• Move from the centre patch to the north. Its distance is stored as 1
• Move from the current patch (north) to the east. Distance is stored as 2
• Move from the current patch (north->east) to the south (this is the centres east neighbour). Distance is stored as 3.
• Move from the centre patch to the east. This patch is flagged as processed at a distance of 3. The current distance is 1 therefore ignore previous processing and process again.
  

That looks better. It is a uniform area around the centre of LOD. However, it is still not correct. To render the next area of lower LOD I step 1 level up the patch tree and render outwards from here. Patches are not rendered if any of its children have been rendered as this would cause an unacceptable overlap and possible onscreen artifacts. This results in huge holes in the terrain.

A simple rule of "if one of my siblings is being rendered at a specific LOD, I must also render at that LOD even if I am not within range of the centre of LOD" fixes this problem.

The above example has resulted in code which is smaller, simpler to understand and faster to run than the previous version based on actual distance.

Planet rendered from orbit showing LOD

Same planet rendered from atmosphere, again showing LOD

Remove back facing patches before the render pipeline.

In a previous post, Procedural Planet - Subdividing Cube I mentioned how I remove complete patches which are facing away from the camera before the API culled their polygons in the render pipeline.

"Using frustum culling is not enough to remove any unrendered polygons from the planet. When close to a planet it can look like a flat terrain, just like the earth does for us as we stand on it, but from height it can be seen that the planet is in-fact spherical. With this knowledge it is possible to mathematically remove allot of the planet patches which are on the opposite side of the planet. With Back face culling the API would remove these anyway, however it would be very wasteful to pass these invisible faces down the render pipeline. By using a DotProduct with the LookAt vector of the camera and the Normal of the planet patch translated to model space, it is very simple to ignore these patches."

This code was part of what was lost from Decade with the recient SVN blooper, and therefore had to be rewritten. Despite having implemented the functionality about 18 months ago it took me some time to grasp the idea. I feel that a more technical post would be useful and will hopefully help anyone else when implementing similar functionality.

Anyone with experience in graphics programming will be familar with the concept of backface culling. As each polygon is passed down the rendering pipeline, it is tested to see if it is facing torwards or away from the camera. Since polygons facing away from the camera cannot be seen there is no need to render them. This concept can be applied when rendering a planet, however instead of just testing on a polygon by polygon basis, I test patch by patch. This allows me to ignore large collections of polygons with 1 test instead of the previously described polygon test.

How is this achieved?

Two pieces of information are required in order to test if a patch is facing towards or away from the camera.

1) The camera look-at vector. This is maintained by the camera object and updated as the camera moves and rotates around the world.
2) The normal vector of the patch. I calculate this when the patch is created by doing a CrossProduct of the vectors of 2 sides of the patch.

If the values of A,B,C and D in the above image are

l_vA: X=-11.547007 Y=11.547007 Z=-11.547007
l_vB: X=-6.6666670 Y=13.333334 Z=-13.333334
l_vC: X=-8.1649666 Y=8.1649666 Z=-16.329933
l_vD: X=-13.333334 Y=6.6666670 Z=-13.333334


The normal of the patch can be calculated using the following code.

m_vNormal = CalculateNormal(l_vB, l_vA, l_vC);

where

CVector3 CalculateNormal(CVector3 p_vOne, CVector3 p_vTwo, CVector3 p_vThree)
{
    CVector3 l_vA = p_vTwo - p_vOne;
    CVector3 l_vB = p_vThree - p_vOne;
    return Normalize(CrossProduct(l_vA, l_vB));
}

CVector3 CrossProduct(CVector3 p_vVectorOne, CVector3 p_vVectorTwo)
{
    CVector3 l_vResult;

    l_vResult.X = p_vVectorOne.Y * p_vVectorTwo.Z - p_vVectorOne.Z * p_vVectorTwo.Y;
    l_vResult.Y = p_vVectorOne.Z * p_vVectorTwo.X - p_vVectorOne.X * p_vVectorTwo.Z;
    l_vResult.Z = p_vVectorOne.X * p_vVectorTwo.Y - p_vVectorOne.Y * p_vVectorTwo.X;

    return l_vResult;
}

The calculated normal would be X=0.44721335 Y=-0.44721335 Z=0.77459687

With this information, as each patch is rendered, a simple DotProduct of these 2 pieces of information returns a floating point value. If this value is less than 0, the patch is facing away from the camera and therefore it and all its child patches can immediately be discarded.

float l_fDotProduct =  DotProduct(m_vNormal, p_pCamera->get_LookAt());
if (0.0f > l_fDotProduct)
    return;

where

float DotProduct(CVector3 p_vVectorOne, CVector3 p_vVectorTwo)
{
    return p_vVectorTwo.X * p_vVectorOne.X + p_vVectorTwo.Y * p_vVectorOne.Y + p_vVectorTwo.Z * p_vVectorOne.Z;
}

One more issue must be dealt with before we have a complete solution. The above works well for static objects where all vertices are relative to the origin. However what will happen if the object is rotating?
The answer of that question is shown in the following image.

It may not be obvious from an image instead of a realtime demo so I will try and explain. The above implemented patch culling is processed on
the raw sphere data. This is equilavent to removing back-facing patches (anything from the back of the sphere, relative to the camera, is removed), then rotating the sphere on the Y axis (in the image above
the sphere is rotated by 130 degrees) and then rendering with the API culling all back-facing polygons. This order is obviously incorrect.

The more correct sequence would be to rotate the sphere, remove all back-facing patches, then render remaining patches and allow the API to remove any back-facing polygons in the front-facing patches. Since rotation occurs in the render pipeline it isn't possible for us to rotate before we remove the back-facing patches.

The solution is to multiply the camera look-at vector by the modelview matrix. This is equilavent to transforming the camera by the same values that will be applied to the sphere, resulting in the correct back-facing patches being removed, regardless of what rotation/translation/scaling is applied to the sphere.


float l_fDotProduct =  DotProduct(m_vNormal, p_pCamera->get_LookAt() * p_pGraphics->get_Matrix(eModelView));
if (0.0f > l_fDotProduct)
    return;


(Note: Since p_pCamera->get_LookAt() * p_pGraphics->get_Matrix(eModelView) will yield the same value for every patch, it would be better to calculate this once per frame for each planet that is being rendered. This value can then be used within the test on each patch in the planet.)

Two steps forward, three steps back

In order to not admit what was probably a result of my own stupidity, I'm going to blame SVN. Regardless of what happened, it has now become apparent to me that the version of Decade Engine Source that I have isn't the latest version. It is missing my implementation of GPU planet generation. I've checked my online repository, backup disks etc .....

This means that I have to recode those sections. A chore, but on a positive side, I know the pit-falls and issues I encountered the last time, and can hopefully design around these and end up with better solution.

As per my previous post, I have also started IPhone and Android development. I shall be working on multiple projects at the same time, and rather than mix it all up on this blog, I have created a sister blog to this for Decade Mobile. Any updates which are specific to the mobile platforms shall be posted there.

Back online in Sydney

Hello again. Its been far to long! I am now settled and living in Sydney and have decided that the time to resume Decade is far overdue.

Development of Decade shall continue as before with procedural planery bodies, but over the past few months I have started to program IPhone/IPad and Android, therefore I think it would be fun to create a mobile Decade Engine and try to make some simple but fun phone games.

Let the adventure begin (yet again!)
Ciarán