Geosphere generation

In one of the projects I’ve been working on recently, I needed to programmatically generate a geosphere primitive.Photo courtesy of Global Nerdy  I had some old almost-complete code for this purpose, but as I dug it up I discovered – like many of us do – that 3 years ago my coding style was horrible, and the code itself is cryptic.  

So instead of modifying my own code, I went online to try to find the solution I would have based the code upon (and then not placed in the comments *smack*).  But little luck, all my Google searches for geosphere algorithm, geodesic sphere generation and the like didn’t find me what I wanted.  So I eventually decrypted my own code, and I’ll post it here for fellow lost travellers.  I believe most of this comes from Paul Bourke, but I’m not sure of the origin of the icosahedron code.

The basic idea of a geosphere is quite simple: Take a 3D primitive, subdivide each of its faces over a number of repetitions to increase the resolution, then push every vertex out to the radius of the sphere you’re trying to approximate.  Almost any 3D primitive is good for the job, but not all will ensure the resulting faces are equal sizes.  I believe the best 2 primitives to use are either the Octahedron (8-sided) or Icosahedron (20-sided).  While they won’t result in exactly equal faces, the difference is so small it’s never really noticeable.  So, on to the code:

public Geosphere(float radius, Vector3 position, int depth)
{
GenerateIcosohedron();
SubdivideToDepth(depth);
PushVerticesOutToRadius(radius);
}

This is all pretty straightforward.  Let’s look closer at each method:

private void GenerateIcosohedron()
{
//Calculate variables for algorithm
float piOnFive = (float)Math.PI / 5.0f;
float piOnTen = (float)Math.PI / 10.0f;

float unitRadius = 1.0f;

float insideExtent = (float)Math.Cos(piOnFive);
float sideLength = 2 * (float)Math.Sin(piOnFive);
float Cx = (float)Math.Cos(piOnTen);
float Cz = (float)Math.Sin(piOnTen);
float H1 = (float)Math.Sqrt((sideLength) * (sideLength) - unitRadius);
float H2 = (float)Math.Sqrt((insideExtent + unitRadius) * (insideExtent + unitRadius) - insideExtent * insideExtent);
float Y2 = 0.5f * (H2 - H1);
float Y1 = Y2 + H1;
float r = unitRadius;
float s = sideLength;
float h = insideExtent;

//create the icosahedron
Vertices = new List();
Vertices.Add(new Vector3(0, Y1, 0)); //a
Vertices.Add(new Vector3(0, Y2, r)); //b
Vertices.Add(new Vector3(Cx, Y2, Cz)); //c
Vertices.Add(new Vector3(0.5f * s, Y2, -h)); //d
Vertices.Add(new Vector3(-0.5f * s, Y2, -h)); //e
Vertices.Add(new Vector3(-Cx, Y2, Cz)); //f
Vertices.Add(new Vector3(0, -Y2, -r)); //g
Vertices.Add(new Vector3(-Cx, -Y2, -Cz)); //h
Vertices.Add(new Vector3(-0.5f * s, -Y2, h)); //i
Vertices.Add(new Vector3(0.5f * s, -Y2, h)); //j
Vertices.Add(new Vector3(Cx, -Y2, -Cz)); //k
Vertices.Add(new Vector3(0, -Y1, 0)); //l

//create the indices list
_IndexTriangles.Add(new IndexTriangle(0, 1, 2));
_IndexTriangles.Add(new IndexTriangle(0, 2, 3));
_IndexTriangles.Add(new IndexTriangle(0, 3, 4));
_IndexTriangles.Add(new IndexTriangle(0, 4, 5));
_IndexTriangles.Add(new IndexTriangle(0, 5, 1));
_IndexTriangles.Add(new IndexTriangle(1, 8, 9));
_IndexTriangles.Add(new IndexTriangle(9, 2, 1));
_IndexTriangles.Add(new IndexTriangle(2, 9, 10));
_IndexTriangles.Add(new IndexTriangle(10, 3, 2));
_IndexTriangles.Add(new IndexTriangle(3, 10, 6));
_IndexTriangles.Add(new IndexTriangle(6, 4, 3));
_IndexTriangles.Add(new IndexTriangle(4, 6, 7));
_IndexTriangles.Add(new IndexTriangle(7, 5, 4));
_IndexTriangles.Add(new IndexTriangle(5, 7, 8));
_IndexTriangles.Add(new IndexTriangle(8, 1, 5));
_IndexTriangles.Add(new IndexTriangle(11, 6, 10));
_IndexTriangles.Add(new IndexTriangle(11, 10, 9));
_IndexTriangles.Add(new IndexTriangle(11, 9, 8));
_IndexTriangles.Add(new IndexTriangle(11, 8, 7));
_IndexTriangles.Add(new IndexTriangle(11, 7, 6));
}

I wish I knew where I got this code from.  I believe the vertex generation stuff is from Paul Bourke, while I calculated the indices myself.  Note that I’m using an IndexTriangle class to store the indices.  This makes it easier to subdivide the triangles in the next step.  I’m also assuming the icosohedron has a unit radius (1.0f), because I’ll be pushing vertices out to the radius lateron.

private struct IndexTriangle
{
public IndexTriangle(int a, int b, int c)
{
this.a = a;
this.b = b;
this.c = c;
}
public int a;
public int b;
public int c;
};

The next step is to subdivide each triangle and repeat to a desired depth.  The subdivision I’m using splits each triangle into 4 smaller ones like so:

subdivide

private void SubdivideToDepth(int depth)
{
for (int i = 0; i < depth; ++i)
{
List newIndexTriangles = new List();
foreach (IndexTriangle indexTriangle in _IndexTriangles)
{
Vector3 newVectorOne = Vector3.Lerp(Vertices[indexTriangle.a], Vertices[indexTriangle.b], 0.5f);
Vector3 newVectorTwo = Vector3.Lerp(Vertices[indexTriangle.b], Vertices[indexTriangle.c], 0.5f);
Vector3 newVectorThree = Vector3.Lerp(Vertices[indexTriangle.c], Vertices[indexTriangle.a], 0.5f);
Vertices.Add(newVectorOne);
Vertices.Add(newVectorTwo);
Vertices.Add(newVectorThree);

IndexTriangle newTriOne = indexTriangle;
newTriOne.b = Vertices.IndexOf(newVectorOne);
newTriOne.c = Vertices.IndexOf(newVectorThree);
newIndexTriangles.Add(newTriOne);
newIndexTriangles.Add(new IndexTriangle(Vertices.IndexOf(newVectorOne),
indexTriangle.b, Vertices.IndexOf(newVectorTwo)));
newIndexTriangles.Add(new IndexTriangle(Vertices.IndexOf(newVectorTwo),
indexTriangle.c, Vertices.IndexOf(newVectorThree)));
newIndexTriangles.Add(new IndexTriangle(Vertices.IndexOf(newVectorOne),
Vertices.IndexOf(newVectorTwo), Vertices.IndexOf(newVectorThree)));


}
_IndexTriangles = newIndexTriangles;
}
}

And finally, push all the created vertices out to the radius, so we end up with the approximation of a sphere

private void PushVerticesOutToRadius(float radius)
{
List newVertices = new List();
foreach (Vector3 vertex in Vertices)
{
float rootrad = (float)Math.Sqrt(vertex.X * vertex.X +
vertex.Y * vertex.Y +
vertex.Z * vertex.Z);
newVertices.Add(new Vector3(vertex.X * (radius / rootrad),
vertex.Y * (radius / rootrad),
vertex.Z * (radius / rootrad)));
}
Vertices = newVertices;
}

The code is looking pretty ugly posted up at the moment, hopefully I’ll have time to tidy it up tonight.

One thought on “Geosphere generation

Leave a Reply

Your email address will not be published. Required fields are marked *