This is going to be a pretty short post because at this point in my learning experience I was still following a tutorial.
Reimer does a great job in these tutorials. XNA as a framework threaded the perfect gap for my level of knowledge. It abstracted
away enough of the low level boilerplate to keep progress fast whilst leaving enough intact to provide a great learning experience.
It did not take long to start showing results:
From there I decided to begin construction of the planet.
Taking inspiration from this fantastic blog i decided to follow his suggestion and implement a cubemap in order to project the terrain to the sphere.
All that was involved was to create 6 terrain patches and rotate them to form a cube.
I was able to port a terrain generator from my previous experiments with Roblox,
Then to project the cube to a sphere the cube verticies were normalised and then multiplied by the desired planet radius plus the height of the terrain at that point.
After some experimentation this is what I had to show for it:
In my first efforts, planets were round, oceans were blue and the land was green which had me pretty excited but failed to bring about any interest from anyone else. (I'm sure any programmer can relate to this!)
What it needed was more detail. I wanted to be able to fly in at orbital speeds, re-enter and land just like in Orbiter. The first complication was that representing a full sized planet all at once is impossible on current hardware and likely to be impossible for the foreseeable future, planets are just too big.
Taking inspiration from Steven Wittens fantastic blog i decided to have a go at "Chunked LOD"
The paper which describes it can be found here, However I will run through it and describe how I implemented it.
The concept:
The basic idea is that from a distance you can't see the fine details of the terrain anyway, so a simpler version is displayed instead. When the camera moves closer to the terrain and higher detail is needed, the terrain patch is split into 4 sub patches where each one is one quater the size of the original and those are displayed instead.
This fits nicely into a data structure called a "Quad tree."
Just like any other tree structure, the "Quad tree" is a hierarchy; A root node at the top which has child nodes who in turn have child nodes themselves.
What makes it a quad tree is that (as the name sugests) each node has 4 children.
This fits in perfectly with the previous cube approach to rendering the planet where the six sides of the cube are formed by six terrain parches representing six root nodes, for six independant quad trees.
At each draw call the algorithm starts from the top of the tree and ask the question;
Is the level of detail for this node sufficient?
If Yes: Render this chunk.
If No: Subdivide this node and recurse the test down to the children who may also need to be split.
All working properly it should look like this:
Another great advantage of this approach is
that it becomes trivial to implement highly effective culling of nodes
which lie outside of the camera's frustum - if a parent node's bounding
sphere lies outside of the camera's frustum, then all of it's children do
as well. We don't even need to visit them to know they will not need
to be rendered.
The implementation:
I'm just going to show stubs here because the implementation is dependent on the engine/platform so I'll get you to fill out the methods depending on what you need.
As you move further away from the surface, details get harder and harder to distinguish.
Eventually as you move further out into space, entire mountains can take up less than 1 pixel on your monitor and effectively become invisible.
Therefore if you had a low detail version of a terrain and a high detail version - at a certain distance the two meshes become impossible to tell apart.
So the plan is to determine when a terrain patch is far enough away for the difference between it, and a terrain patch of a different level of detail to be invisible, and to then swap it for the other chunk.
By taking into account how far away we are from the mesh, the estimated difference between the two meshes (more on this later) , the field of view of your camera and the size of your monitor it is possible to calculate how big the differences between 2 levels of detail appear on your monitor in screen space.
In the paper I referenced above, this is called the "Sceen space error metric"
In order to provide the position of the camera from the terrain chunk to calculate this metric it made sense for me to pass the camera object in itself since we also need it's bounding frustum for frustum culling.
The other thing necessary for the metric which isn't shown here is the width of the game window itself, how you want to go about providing that is up to you, XNA made this easy in my case by just querying the global GraphicsAdaptor.
The last step is to choose a tolerated error above which the higher quality mesh is rendered instead. This tolerated error needs to be
chosen so the detail swap happens before the difference between levels
becomes
obvious but not so early as to make the whole operation pointless.
class Camera
{
BoundingFrustum frustum;
Vector3 Position;
public BoundingFrustum GetFrustum()
{ return frustum; }
public Vector3 GetPosition()
{ return Position; }
}
class TerrainMesh
{
BoundingSphere sphere;
TerrainMesh()
{
GetHeightmap();
BuildMesh();
}
~TerrainMesh()
{ /*do stuff */}
private void BuildMesh()
{ /*do stuff */}
private void GetHeightmap()
{ /*do stuff */}
public void Render()
{ /*do stuff */}
public BoundingSphere GetBoundingSphere()
{ return sphere; }
}
As for the quad tree node class, I have filled out the basic recursive structure of the drawing here. Once again I have stubbed out the engine specific details, so fill out the methods as you need.
class QuadTreeNode
{
QuadTreeNode[] Children;
TerrainMesh Mesh;
// If Lod is sufficient draw this chunk otherwise test child nodes.
void Draw(Camera camera)
{
// if the chunk's bounding sphere is either intersecting or contained within the frustum:
if (camera.GetFrustum().Contains(Mesh.GetBoundingSphere()) != ContainmentType.Disjoint)
{
// if this chunk's detail is sufficient then draw it, also now
// is a good time to free up some memory by destroying child nodes
// since they won't be needed - unless you want to cache them.
if (LodIsSufficient(camera.GetPosition()))
{
Mesh.Render();
if (ChildrenAreValid()) DestroyChildNodes();
}
else
{
// if the children have already been created then pass the draw call down.
if (ChildrenAreValid())
foreach (QuadTreeNode chunk in Children) chunk.Draw(camera);
else
{ // otherwise we need to construct them.
// it's a good idea to do this on another thread to prevent halting the rendering.
// draw this chunk this frame - children should be ready next time, otherwise we will end up back here.
ConstructChildNodes();
Mesh.Render();
}
}
}
}
// Return true if chunk is detailed enough
bool LodIsSufficient(Vector3 CameraPosition)
{ /*do stuff */ return true; }
// return true if all children are constructed and ready to be rendered.
bool ChildrenAreValid()
{ /*do stuff */ return true; }
// Construct children if they do not already exist.
void ConstructChildNodes()
{ /*do stuff */}
// Release all the child's resources and pass the destroy call down to child's children.
void DestroyChildNodes()
{ /*do stuff */}
}
The complications to keep in mind:
Remember to define a maximum tree depth to prevent the draw call from creating new children if you are at the maximum detail you want to support.
When it comes to the estimated difference between the best quality mesh and the mesh currently being tested, it is sufficient to guess by starting with some big number ( I chose 2 number of quad tree levels ) and then just divide it by 2 each time you subdivide the mesh.
The reason this parameter exists, however, is to give you more control if you already have your height-maps in advance.
Imagine a terrain patch representing wide open plains, the high quality mesh is going to be very hard to distinguish from the low quality mesh very quickly because there is very little change in elevation happening.
By calculating a per-chunk difference estimation you can cause the LOD swap to happen much later when it is actually needed rather than at some global setting.
However if your height-maps are procedural like mine then this is not an option.
The final piece to my solution was to observe that the closer to the planet you get, the more of the terrain is hidden by the planet's curvature and can be culled from rendering.
This was harder in practice due to the fact that while a chunk may lie over the horizon on a perfectly flat sphere, should it contain hills or mountains which poke up above the horizon then it's possible that visible sections of terrain could be erroneously removed.
Solving this involved being more conservative in my culling where I would only prevent a chunk from being drawn if it was over the horizon by an amount more than the maximum possible height assigned by the terrain generator for that chunk.
This led to some chunks making it through when they should have been culled, but that was a better option than having chunks being culled when they should not have. This still feel hackish and it is an area that I will revisit.
Another problem which needed attention was filling cracks that appear in between two levels of detail.
I decided that this was a problem to be solved later, and was not in the scope of my tech demo. However a good looking solution at this point appears to be storing two levels of detail in each chunk and morphing between the current detail level and the next as the camera moves in.
This solves the ugly popping that happens when a chunk suddenly changed detail level as well as solving the cracking problem as neighboring chunks will also be morphing to fit this chunk.
Moving forward.
I added some vertex normals and this is the result.
The planet is still not textured - but regardless it is a massive leap in quality.
So at this point I turn my mind to the other half of my game - modifying my environment.
I
wanted free form building, similar to Minecraft, where the player places blocks to
create a structure instead of just plonking down prefabricated models on
the surface like in spore for example.
In order to support this gameplay, I required a way to to subtract the shape of a block from the terrain robustly without holes in order to support the digging of a foundation as well as tunneling into the ground.
Suddenly I am feeling very much out of my depth. This approach doesn't feel like it is going to work anymore.
Nothing for it but to dive in and find another way.
My first thoughts are - what if I just build the entire world out of cubes just like in Minecraft?
I understand distortion is inevitable but perhaps I could minimise it?
With no better plans in mind I got to work.
The plan:
The theory was to have six sides to a cube just as before. Each side is a cube itself, though tapering in towards a point like a pyramid.
These cubes slot together with each other to form one big cube.
Then I continued as before to normalize and multiply to project onto the sphere.
It took me a while to implement the LOD, I quickly ran into memory problems.
Rendered naively a 64^3 block volume contains 262144 blocks which requires 24 vertices each.
The reason it is 24 and not 8 is because each face of the cube may have a different texture coordinate and normal vector thus each of the 6 faces needs to store 4 vertices.
Suddenly the mesh is composed of 6.3 million vertices.
Each vertex needs a position, a normal and a texture coordinate which is a total of 8 floats.
At 4 bytes per float each just storing the mesh would take up 201mb of ram, and I wanted a planet much bigger than 64 meters across... I had at least 8 kilometers in mind which has a volume approximately 2 million times bigger! Impossible!
Clearly it was time to optimise.
Making the impossible possible:
The first step was to set up the generation of the cubes such that only faces exposed to a transparent block could generate vertices. The worst possible case now where the most faces would be exposed consisted of a 3d checkerboard pattern which held a volume of n^3 / 2 - 50% reduction.
This pattern is highly unnatural and would not be created by a terrain generator, the average case was now simply just the surface layer of the chunk which held 4096 blocks - sometimes slightly more because of overhangs or tunnels which meant an average case reduction of vertices by a factor of 64 times.
Another optimisation involved storing each vertex as a color instead of a vector, a color in xna is composed of 4 unsigned bytes which is exactly 1/3rd the size of a full vector.
A byte can only store 256 unique numbers which was sufficient to store which point in the volume the vertex was describing, a volume of 64 blocks required only 65 points on any axis - much less than 256.
Once in the shader i converted the vectors back to full sized vectors before offsetting them by a single global vector which the chunk being drawn had supplied.
A similar tactic was also used for the normal vector which combined resulted in reduced memory usage per vertex from 32 bytes to 16 bytes.
The final step was employing run length encoding of the volume which resulted in sometimes astonishing reductions in resource usage, a good introduction is outlined here.
Combined these optimisations allowed interactive frame rates rendering a planet 8192 blocks wide.
The results:
Then, with my improving shader programming skills, I decided to implement atmospheric scattering which is to this day one of the the most fiddly, frustrating things I have ever done. Regardless I managed to get it done and here is a video demonstrating my progress.
More complications:
However it became very clear that inevitable distortions were pretty bad where the cube's met on the surface.
"Perhaps I can distort it in reverse to smooth it out" I tell myself. I tried this and while it "kind of" worked (through a complicated process of projecting the cube to a cylinder, rotating it and then projecting it back to a sphere) in the sense that the cubes were now properly shaped cubes wherever iI stood, it resulted in landscapes bending and twisting awfully as you traveled while they morphed into shape - I would not recommend this approach to anyone without even considering what such an approach would mean to a physics engine.
I discovered the transvoxel algorithm here while surfing around looking for inspiration some weeks later which provided the spark of inspiration I needed.
It builds on an algorithm named "marching cubes" which is commonly found in medical scanners such as CT and MRI machines. These machines, as opposed to an xray for example, scan a volume as opposed to a single image.
The problem this algorithm solves is how to efficiently extract a mesh representing the isosurface of a three-dimensional scalar field (voxel volume) so that it can be 3d rendered on a computer.
Marching cubes:
When I first started the implementation, I must admit, it
stood outside of my abilities as a programmer at that point - however
that had been true for each step of the way so far already. I firmly
believe that one of the best ways to learn is to push yourself to achieve something that is just outside of your capabilities until you grow in understanding and
achieve it.
Anyway, the approach is in some ways similar to the previous blog post where I built a terrain out of cubes, except in this case, instead of just plonking down a cube when it is discovered that it lies on the border between empty and solid space, a sample is taken of the densities on the 8 corners of the cell and if one corner is solid and one corner is not - then the isosurface must travel through this cell.
Marching cubes is the implementation of a look up table which holds information on what triangles need to be created to triangulate the isosurface within that cell for all the possible combinations of the corners being in solid or empty space.
You can see from this image where the corners represented by black dots are in solid space, that for each case the triangles needed to form a watertight mesh separating the solid areas from the open corners is defined. Each of these cases can be rotated in order to satisfy any combination of empty / solid corners.
Once all the cells are processed in this way the triangles are combined to form a mesh.
This approach has a flaw however in that by taking each cell on it's own
there are a huge number of duplicated vertices which inflates memory usage in what is already an enormously memory intensive process. Eric Lengyel for the C4
game engine provided a fantastic and thorough explanation here in
his dissertation of how it could be done - I won't try to replicate his work here though I will post the result of my implementation and provide the actual source code here
You will notice in the previous image a kind of stair stepping going on. It is blurred because the vertices and the normals are being reused but the effect is still clear and by having a closer look at the geometry it becomes apparent that the mesh is in fact very angular.
Solving this was also described in Eric's paper, the key observation is that all vertices which need to be created for any given case are created in the center of one of the 12 edges of the cube. By interpolating the vertex position along the edge towards one corner or the other depending on how close the density at that corner was to being considered empty or solid as opposed to the other corner meant that the isosurface could be represented much more exactly.
There are still a few bumps here and there caused by my simplex noise generator which needs a bit of tweaking, however, the algorithm itself is working properly.
The transvoxel algorithm:
The inevitable problem, just as it was last time, is the planets are just too big. The level of detail adjustments are not optional - they are required.
The issue is, just as there were with the height map style terrain, two levels of detail next to each other,. This means there will be cracks forming in the mesh where they do not line up.
To solve this is not as simple as just sewing the chunks together or interpolating the levels of detail because the geometry is so much more dynamic.
The is what the transvoxel algorithm sets out to solve. It provides a whole host of additional configurations which extend the look up tables to support seamless sewing of adjacent terrain chunks of different detail levels.
This permits me to render a planet this way.
Since the terrain is no longer composed of cubes, having the voxels aligned to the planet surface is no longer required. The planet can be represented as one giant
volume rather than the 6 separate cube faces.
This eliminates the cause
of the awful distortion at the edges where the cube faces met and honestly I was glad to move away from the cube style terrain anyway.
Representing the terrain this way provides a pretty cool way to tunnel but finding a way to have the building blocks sit flush with the terrain sounded tricky.
Moving on.. again..
In the end what pushed me away from this approach was the sheer amount of data that needed to be generated in order to generate these terrains in the first place. If you had a pre computed world then I'm sure this approach would work just fine, however, for a procedural world like mine the generation times were just too long.
After some profiling I discovered some 98% of the time used to generate each chunk of mesh was used simply coming up with the density values in the first place, largely because of the requirement to fill n^3 data points instead of n^2 required for a height map.
All was not lost, however, I learned an enormous amount from these experiments and I felt that I had come out the other side a much more skilled programmer than the one who entered.
The solution I decided on was a height map terrain after all. The volumes were just not feasible on this scale. To support tunnelling, I would also have to be cutting holes in meshes.
Once again I was out of my depth but I am nothing, if not stubborn so I set off to learn how to accomplish this.
CSG is the technique of creating complex objects from simple shapes such as spheres and cubes through boolean set operations such as difference, intersection and union.
This concept can be extended, without too much trouble, to include freeform meshes and not just primitive shapes. To keep things simple I'll keep referring to this extended version as CSG.
It was at this point that I started off the beaten track. There was ample information on the concepts but very little in the way of concrete implementation which meant working a lot of this out myself.
So in an effort to smooth the way for those following me, grab a coffee and some snacks while I do my best to describe my approach as best as I can - it's going to be a long one.
The plan:
The approach is simple enough, these things often are - it is the implementation which is hard.
Here are 2 meshes, mesh A which is a red cube and mesh B which is a blue sphere.
I want to support the 3 main operations which are:
Difference: Subtract one meshes volume from the other, for example if mesh B was to subtract it's volume from mesh A, the algorithm must delete any portion of mesh A which lies inside of mesh B as well as deleting any portion of mesh B which does not lie inside of mesh A's volume.
Intersection: Return only the volume that both meshes share, any part of either mesh which does not lie inside the volume of the other mesh must be deleted.
Union: combine both meshes into one, deleteing any part of either mesh which lies inside of mesh A and mesh B.
To do this i will outline the general approach before going into detail.
1) Every polygon from either mesh which crosses inside the volume of the other mesh needs to be split along the intersection of the two meshes.
2) After splitting is complete all polygons should be classified as either inside or outside the volume of the other mesh and either kept or discarded depending on the operation.
To do this a data structure needs to be generated which represents the volume of the mesh in order to support classifying the polygons as inside or outside.
The data structure in use here is the "Binary space partition tree" (BSP tree) and revolves around classifying polygons as crossing a plane formed by another polygon or lieing to either side of that plane. A BSP tree will solve both step 1 and step 2.
I will go into detail and explain how to construct and use a BSP tree later, first I will describe how to classify a polygon to a plane and how to split a polygon against that plane which is a central concept used as part of constructing a BSP tree.
Classifying polygons:
At some point when one mesh intersects another, any polygons which lie inside or outside of the other mesh's volume need to be either kept or deleted depending on the operation.
This means a way to describe which side of a plane a polygon lies on is required.
The code to calculate this can be simplified greatly by observing that triangles are by definition planar.
By splitting a triangle any polygons created as a result of that split must also be planar.
Because I am working with meshes entirely composed of triangles then only planar polygons need ever be considered.
Because of this, determining which side of a plane a polygon lies on is reduced to simply calculating the distance to the plane for all points in the polygon and testing if they are all are at a positive distance from the plane or at a negative distance from the plane.
To calculate the distance of a point to a plane use the folowing formula.
DotProduct(point , PlaneNormal ) + D
where D is the dot product calculated between any of the points making up the polygon which the plane belongs to and the normal vector of the polygon being tested.
If the distance is positive or negative or 0 it can be classified as inside, outside or "on" the plane respectively.
A simple implementation:
//
public PointPlaneClassification ClassifyPointToPlane(Vertex Point)
{
double Distance = GetDistance(Point.Position);
if (IsEqual(Distance, 0.0, epsilon))
return PointPlaneClassification.Point_On_Plane;
if (IsGreaterThan(Distance, 0.0, epsilon))
return PointPlaneClassification.Point_Behind_Plane;
if (IsLessThan(Distance, 0.0, epsilon))
return PointPlaneClassification.Point_In_Front_Of_Plane;
return PointPlaneClassification.Point_On_Plane;
}
You might be wondering why I am using functions for the comparisons with 0 instead of the normal comparison operators - this is because CSG is an area where limited floating point precision is a really big problem. I recommend reading an introduction to managing finite floating point precision and how to work around it. (in this case by treating numbers within epsilon of 0 as 0).
It is then straight forward to extend this to be able to classify entire polygons.
//
public PolyIntersectType GetPolygonIntersectionType(Polygon polygon)
{
// short cut, if the two polygons are co-planar then only one point
// is required to know if the polygons coincident or disjoint on one side or the other.
double dot = Vector3Double.DotProduct(Normal, polygon.normal);
if (DoubleHelper.IsEqual(Math.Abs(dot), 1.0, DoubleHelper.epsilon))
{
for (int i = 0; i < polygon.Verticies.Length; i++)
{
PointPlaneClassification classification = ClassifyPointToPlane(polygon.Verticies[i]);
if (classification == PointPlaneClassification.Point_Behind_Plane)
return PolyIntersectType.DisjointPositive;
if (classification == PointPlaneClassification.Point_In_Front_Of_Plane)
return PolyIntersectType.DisjointNegative;
}
return PolyIntersectType.CoPlanar;
}
// not coplanar, test each point against the plane.
// if points exist on bothsides return Straddleing immediately.
int numFront = 0, numBack = 0;
for (int i = 0; i < polygon.Verticies.Length; i++)
{
PointPlaneClassification classification = ClassifyPointToPlane(polygon.Verticies[i]);
if (classification == PointPlaneClassification.Point_Behind_Plane)
numBack++;
if (classification == PointPlaneClassification.Point_In_Front_Of_Plane)
numFront++;
if (numFront != 0 && numBack != 0) return PolyIntersectType.Straddles;
}
// otherwise if all the points lie on the negative side, return disjoint negative.
if (numFront != 0) return PolyIntersectType.DisjointNegative;
// at this point the only possibility is that the polygon is disjoint on the positive side.
return PolyIntersectType.DisjointPositive;
}
Splitting polygons:
Now that it is possible to classify two polygons to a plane, it is time to describe how to split them by it.
In the image on the right the green triangle needs to be
classified as inside or outside the light blue plane which is formed by
the red triangle.
Since it crosses the plane this is impossible, the triangle is both inside and outside of the plane at the same time.
Therefore in order to resolve this the triangle is split into two peices where it crosses
the plane and one peice is classified as inside and the other as outside.
I had a go at it on my own and came up with an algorithm which walked around the edge of a polygon filling the vertex list of the first of 2 new polygons. The 2 new polygons represented the halves of the original polygon in the negative half space and the positive half space.
If an intersection was detected the intersection point was added to both polygon's vertex lists and the then the walking would continue except that it was the other polygon whose list was being filled.
Once the algorithm came back around and crossed the plane again once again the intersection point was added to both lists and the first polygon was once again the polygon whose list was being filled until the algorithm was back to the start.
This worked in the majority of cases, however, there was the odd edge case and handling them all was starting to become a pain so I decided to look around to if there was a better way.
Turns out I was pretty close, the overall concept was sound ,however, my approach was unnecessarily complicated.
One of my favourite references through this entire ordeal was "Real time collision detection" by Christer Ericson.
Basically, his approach was first of all to do away with filling one polygon, then the other and instead enumerate all the permutations that could happen when a vertex crossed a plane and insert verticies to either polygon depending on the combination of the current vertex and the next vertex being "On the plane" , "In front" or "Behind the plane". As a result the code came out much simpler and robust.
The final product was very similar to his example, for a thorough explanation I highly reccomend reading his book - it helped me immeasureably.
//
public Polygon[] SplitPolygon(Polygon polygon)
{
List<vertex> FrontVerticies = new List<vertex>(polygon.Verticies.Length + 1);
List<vertex> BackVerticies = new List<vertex>(polygon.Verticies.Length + 1);
Vertex A = polygon.Verticies[polygon.Verticies.Length-1];
PointPlaneClassification Aside = ClassifyPointToPlane(A);
for (int i = 0; i < polygon.Verticies.Length; i++)
{
Vertex B = polygon.Verticies[i];
PointPlaneClassification Bside = ClassifyPointToPlane(B);
switch (Bside)
{
case PointPlaneClassification.Point_In_Front_Of_Plane:
if (Aside == PointPlaneClassification.Point_Behind_Plane)
{
Vertex I = GetIntersection(B, A);
if (ClassifyPointToPlane(I) != PointPlaneClassification.Point_On_Plane) { }
FrontVerticies.Add(I);
BackVerticies.Add(I);
}
FrontVerticies.Add(B); // Output b to the front side
break;
case PointPlaneClassification.Point_Behind_Plane:
// Edge (a, b) straddles plane, output intersection point
if (Aside == PointPlaneClassification.Point_In_Front_Of_Plane)
{
Vertex I = GetIntersection(A, B);
if (ClassifyPointToPlane(I) != PointPlaneClassification.Point_On_Plane) { }
FrontVerticies.Add(I);
BackVerticies.Add(I);
}
else if (Aside == PointPlaneClassification.Point_On_Plane)
{ // Output a when edge (a, b) goes from ‘on’ to ‘behind’ plane
BackVerticies.Add(A);
}
// In all three cases, output b to the back side
BackVerticies.Add(B);
break;
case PointPlaneClassification.Point_On_Plane:
FrontVerticies.Add(B);
BackVerticies.Add(B);
break;
}
A = B;
Aside = Bside;
}
Polygon[] returnPolygons = new Polygon[2];
if (BackVerticies.Count > 2)
returnPolygons[0] = new Polygon(BackVerticies.ToArray(), polygon.Id , polygon.FlipNormals , polygon.normal);
if (FrontVerticies.Count > 2)
returnPolygons[1] = new Polygon(FrontVerticies.ToArray(), polygon.Id, polygon.FlipNormals, polygon.normal);
return returnPolygons;
}
You might be wondering why the polygon object has a "Id" and "FlipNormals" field - this all becomes important in the next section.
Binary space partition trees:
Bringing all these concepts together is the BSP tree.
It is all well and good to be able to decide if two polygons are intersecting and to have a way to split them, but it is not much use without the ability to determine if a given polygon lies inside a meshes volume.This is required so that a particular polygon can be kept or discarded as required.
A data structure which provides what is needed is the "Binary Space Partition" (BSP) tree.
I'm going to demonstrate how they work as a 2d diagram, to generalize into 3d space imagine that each line in the polygon is a triangle. The theory works exactly the same.
Each edge in this polygon is given a letter and has a normal vector indicated by a short line pointing out from the center of the edge, this points the way to the "outside" or "empty space".
I will refer to this polygon below simply as "the polygon"
The way a BSP tree is built is one edge at a time. Each time an edge is added it's plane is used to further divide the world into two equal halves, the inside half and the outside half.
I'm going to start with edge A and work through alphabetically.
Below is A's plane Indicated in red - demonstrating how it splits space into 2 halves, the upper half (outside) and the lower half (inside).
It is possible to be inside A's plane and outside of the polygon as a whole at the same time.
Now to add the next edge "B". Always start from the top of the tree.
Is this edge on the inside of A's plane? Or on the outside?
Since the outside is the side which the normal is pointing and since B lies on the opposite side, B is therefore on the "Inside" half of A's plane.
So B is inserted into the tree at A's inside node.
And the graphics representation:
The next step is to add C.
The process is the same, is C inside or outside of A?
C is inside.
Therefore the "inside" half of the tree is tested next.
Is C inside or outside of B?
C is outside.
Therefore C is inserted into the tree on B's outside.
Now it gets a little more tricky.
Edge D needs to be inserted, So begin at A as usual.
Is D inside or outside of A?
D is inside, So B is tested against next.
Is D inside or outside of B?
Trick question? D crosses over B, it is both inside and outside.
The solution is to split D into two halves and send each half down the appropriate side.
I'll name the inside half D2 to differentiate it from D which will go down the outside half.
D will then be tested against C, where it lies on the inside half and is added to the tree.
D2 is simply added to B's inside half.
And then last of all E.
E is inside A.
E is inside of B
E is inside of D2, so it is added to the tree on D2's inside half.
And there is the completed BSP tree!
The code:
I'll outline an example here of how this might work.
First some enums to describe the types of intersections.
class BspTree
{
Plane DividingPlane;
BspTree PositiveHalf;
BspTree NegativeHalf;
private BspTree(Polygon P)
{
DividingPlane = new Plane(P);
}
public static BspTree ConstructTree(List<Polygon> PolygonList)
{
BspTree root = new BspTree(PolygonList[0]);
for (int i = 1; i < PolygonList.Count; i++)
{
root.PartitionPolygon(PolygonList[i]);
}
return root;
}
public bool PolygonLiesInsideVolume(Polygon P)
{
PolyIntersectType Type = DividingPlane.GetPolygonIntersectionType(P);
switch (Type)
{
// to stop c# complaining about not all code paths returning a value.
default: throw new ArgumentOutOfRangeException();
case PolyIntersectType.DisjointNegative:
// if at leaf node polygon is in empty space.
if (NegativeHalf == null) return false;
// test negative half, return the result.
return NegativeHalf.PolygonLiesInsideVolume(P);
case PolyIntersectType.DisjointPositive:
// if at leaf node - polygon is in solid space.
if (PositiveHalf == null) return true;
// test positive half, return the result.
return PositiveHalf.PolygonLiesInsideVolume(P);
case PolyIntersectType.Straddles:
// partition polygon by plane.
Polygon[] Splits = DividingPlane.SplitPolygonByPlane(P);
bool negHalf = false; bool posHalf = false;
if (PositiveHalf != null && Splits[0] != null) posHalf = PositiveHalf.PolygonLiesInsideVolume(Splits[0]);
if (NegativeHalf != null && Splits[1] != null) negHalf = NegativeHalf.PolygonLiesInsideVolume(Splits[1]);
// if at a leaf node, if either half of the polygon is discovered
// to be inside the volume return true.
if (negHalf && PositiveHalf == null && Splits[0] != null) return true;
if (posHalf && NegativeHalf == null && Splits[1] != null) return true;
// return true if either halves are in solid space.
return negHalf && posHalf;
case PolyIntersectType.CoPlanar:
return false;
}
}
private void PartitionPolygon(Polygon P)
{
PolyIntersectType Type = DividingPlane.GetPolygonIntersectionType(P);
switch (Type)
{
case PolyIntersectType.DisjointNegative:
// partition to negative half space.
DividingPlane.GetPolygonIntersectionType(P);
PartitionToNegativeHalf(P);
break;
case PolyIntersectType.DisjointPositive:
// partition to positive half space.
PartitionToPositiveHalf(P);
break;
case PolyIntersectType.Straddles:
// split polygon and partition to both sides.
Polygon[] Splits = DividingPlane.SplitPolygonByPlane(P);
if (Splits[0] != null) PartitionToPositiveHalf(Splits[0]);
if (Splits[1] != null) PartitionToNegativeHalf(Splits[1]);
break;
case PolyIntersectType.CoPlanar:
// polygon shares the same plane as this node.
// no need to partition at all.
break;
}
}
private void PartitionToNegativeHalf(Polygon P)
{
// if the positive half has not been used yet - create it.
if (NegativeHalf == null)
{
NegativeHalf = new BspTree(P);
}
else // otherwise partition polygon to positive half.
{
NegativeHalf.PartitionPolygon(P);
}
}
private void PartitionToPositiveHalf(Polygon P)
{
// if the positive half has not been used yet - create it.
if (PositiveHalf == null)
{
PositiveHalf = new BspTree(P);
}
else // otherwise partition polygon to positive half.
{
PositiveHalf.PartitionPolygon(P);
}
}
}
Putting BSP trees to use: point in mesh test.
BSP trees are useful for a couple of reasons, the first is that they provide a really effective way of determining if a point is inside a mesh and secondly they as a side effect of the method of construction create all the splits we need to perform CSG on two meshes.
First I'll demonstrate how to use a BSP tree to decide if a point lies inside the mesh or not, it is a similar process to construction - just drop the point through the tree.
The red diamond is the point to be tested.
Starting from the top of the tree, the point is inside A, outside of B and outside of C. since it lands on an "outside" node it is therefore it is located outside of the mesh. In just 3 tests!
In this case the red diamond is inside A, outside of B, inside of C and inside of D - since the point lands on an "inside" node the point is inside the mesh.
Putting them to use 2: Constructive solid geometry.
The main purpose of using BSP trees, however, is that because splitting
polygons is a part of constructing them, merging two BSP trees created from two separate meshes will result in splits in all the
places needed to perform any CSG operation, by keeping a copy of the
original BSP trees from the mesh they can be used as demonstrated above
to classify the polygons which were split as either inside or outside
and then discard or keep them as required.
First off I am going to overlay a new triangle with sides X Y Z.
What I want to do is subtract this triangle from the first polygon so that it cuts a little triangle shaped hole - this is a "difference" operation.
First a new BSP tree needs to be built for the new triangle. If you have been following along the process should be pretty clear by now.
A convex mesh will always have all of it's edges as "inside" which effectively devolves the bsp tree to a linked list. There are ways to break this up a bit to improve performance but i will cover that at the end.
What needs to be done next is to take a copy of the polygon's BSP tree so that it can be merged with the triangle's BSP tree and still be able to use the original to classify polygons as inside or outside at the end.
Each edge of the triangle is inserted into the copy of the polygon's BSP tree one edge at a time just like before with one twist; if an edge from the triangle crosses an edge from the polygon it needs to be split just as normal but - it also must split the edge from the polygon that it crossed. ( This is only strictly required if the polygons, not their planes are intersecting - see optimisation 2 below.)
Edge X from the triangle is inserted first, X crosses edge A.
X is split the two halves are passed down both the inside and outside of A as normal.
However because A was crossed by X it also needs to be split - but it doesnt need to be partitioned. It is stored in the same node (A , A2)
Y is partitioned cleanly down the tree until edge E.
Y is split into two halves and passed down both sides, E is also split just as A was.
Z is able to be partitioned as normal since it doesnt intersect with anything, It is split by A but A is not split by Z.
Remember before when I commented on polygons having an "ID" member? Now that the two meshes polygons are all mixed up together it is important to be able to tell them apart in order to be able to separate them again later.
If you look at the last diagram you can see that all the splits that need to be made in order to cut a hole into the polygon have already been made - now all that remains is to decide which polygons need to be thrown out.
The way I did this was to start at the top of the tree and work my way down copying all the polygons from each mesh into separate lists. Once I had all the polygons I then, one at a time, tested them against the original BSP from the other mesh in order to determine if they lay inside or outside of the mesh's volume.
In order to create a triangle shaped hole, all the triangle edges that are located outside of the polygon are discarded and all of the polygon edges that are located inside of the triangle also need to be discarded.
I'll leave the implementation of that to the reader,
The result is this:
Improving performance:
My initial tests showed satisfactory performance when I intersected 2 small cubes together however once the size of the mesh was increased to what I expected to use in the game, the process just took too long and would cause frame rate skips.
The only way to improve performance in a program is to reduce the
amount of work the computer needs to do - by either working smarter or
reducing the size of the problem.
Optimisation 1: bounding boxes.
The first optimisation came about from realising that the only polygons which could be inside the other mesh or involved in an intersection must pass inside the the intersection of the two meshes bounding boxes. All other polygons could either be thrown out or kept depending on what the operation was right from the beginning.
The triangle's bounding box is marked in orange and the polygon's bounding box is maked in green.
The intersection (just the area shared by both boxes) is shaded in grey.
The edges highlighted in red are the only edges which need to be included in the bsp tree - all the others are automatically classed as "outside".
So in the above example where I wanted to cut a triangle shaped hole, this means edge A, B,C,D can immediately be put in the keep list, edge Z can be immediately discarded without any further processing.
This can result in massive performance gains when the meshes involved are composed of thousands of polygons. For my purposes, where the player would be making edits on the scale of 1 small block and because the size of the terrain was so big, the bounding box intersection only covered a handful of polygons from the terrain. This reduced the number of polygons involved massively and resulted in a big framerate improvement.
It is possible to eliminate even more polygons quickly by using the solid leaf BSP trees that were built before the merge to perform some extra culling.
The key is in realising that partitioning polygons is really quite a lightweight operation compared to splitting them.
Before creating the final merged BSP tree out of all the polygons which survived the bounding box intersection test, each list of polygons is dropped through the other mesh's BSP tree.
If a polygon is classified as either clearly inside or outside - then that polygon can be discarded or kept as required right there and then.
If it crosses a plane then return it back into the list to build the final merged BSP tree with.
Optimisation 2: To split or not to split.
The final tweak which I used came about from observing the possibility of two polygons to be aligned such that the plane of one polygon will cause a split in the second polygon even though the two polygons are not actually touching.
So when a polygon P from mesh A is being inserted into mesh B's bsp tree, P is split and partitioned just as normal - however the polygons in mesh B only need to be split as well if P is physically touching them. It is entirely correct to go ahead and split those polygons in B anyway - but a much higher mesh quality and increased performance can be gained by first running a strict polygon - polygon intersection test before deciding if to split or not.
No optimisation.
You can see how on the left, the result is absolutely valid - but there are massive amounts of unnecessary splits as showing in the wireframe on the right.
With optimisation:
The result is still just as valid but as shown on the wireframe on the right this mesh is much better. This runs much faster because all those extra polygons don't need to be processed.
The results:
Moving on:
It took me a long time to to get this far because it is an extensive topic and I had a lot of learning to do, however, I think the results speak for themselves. This is an excellent result and opens up many options for me going forward.
For example projectiles fired by the player creating damage in buildings or bombs / asteroids creating craters where they land.
Being able to tunnel into heightmap terrain was my main motivation however and having completed this experiment I decided to move back into planet rendering to provide an environment to use my new CSG engine within.
