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.

The Agile Cowboy Programmer

I was recently informed that within my development team I’m often thought of as the ‘rebel’ or cowboy coder. Being, as I am, very attached to programming in rigid methodologies with well-defined practices in place – this revelation surprised me. When I pushed my informant for further information, it seems that I’m viewed this way because I’m the guy writing Agile code and refactoring other classes in order to get them into tests!

This led me to wonder: is a coder truly a rebel if what they are doing is trying to improve the testing systems and reducing software entropy? At what point should respecting the team take precedence over making and keeping the code-base maintainable? Is there such as thing as a ‘good’ cowboy coder?

It’s a hard line for me to draw: on the one hand I’ve a driving need to ensure that the code I create is the best it can be. On the other I feel it’s important for me to try and explain to my fellow teamsters exactly *why* I’m changing their classes, but this would require me to pull back on my own tests and wait for them to begin to understand and write their own tests.

It’s a difficult decision, especially when every agile developer I’ve talked to has had little or no success convincing others to get ‘test infected’ without higher-level support.