My Bug, My Bad #1: Fractal Spheres

posted by Craig Gidney on November 27, 2012

In this series of posts, the dual of the “Unfathomable Bugs” series, I describe interesting bugs that are my own fault.

Today’s bug comes courtesy of me. Thanks, self, this series wouldn’t exist without the generous contributions of entities like you.

The Setup

Recently, I had my first foray into tessellating a sphere (computing an approximation of its surface with a set of triangles). Considering I work at a company known for expertise in a development platform for creating 3d games, you might find this surprising. However, unlike my co-workers, my professional background is actually in application development and, although I did create games as a hobby, they were mainly maps/mods for StarCraft (e.g. Storm the Fort) and WarCraft 3 (e.g. Power Towers) instead of stand-alone programs with custom graphics engines. As a result, even though I understand the math behind 3d graphics very well, I haven’t had to solve a lot of the common problems like sphere tessellation.

I encountered the sphere tessellation problem as part of optimizing and augmenting bubble-rendering code to allow deformations of the bubbles. After implementing some obvious optimizations, like caching the vertices making up the tessellation instead of constantly recomputing them, I noticed that the tessellation itself could be improved. It was constructed out of points evenly spaced by longitude and latitude. As a result, there were significantly more points near the “poles” than there were near the “equator”. You can see this for yourself if you have a globe of the earth: the areas separated by the latitude and longitude lines are smaller near the poles. I figured that, by constructing a more uniform tessellation, I could cut the number of points by about a quarter while maintaining an equivalent level of detail.

After thinking about how to tessellate uniformly for all of a minute, I realized the problem had probably already been solved a thousand times. A quick googling gave an answer: refine a simple starting shape by placing new points at the center of faces. “Of course!” I thought, “What an obvious strategy (in hindsight)!”. I quickly implemented it: one method to create a tetrahedron, one method to refine a shape to be more sphere-like, and one method to iteratively refine the tetrahedron to a desired level of detail. Package it into a class and voila! Bubbles that look… erm… hmmm:

Incorrect bubble spheres

The shader makes it hard to tell what exactly is wrong here, but it looks sort of like some of the triangles are missing or overlapping. I stared at this, and the code I wrote, for quite some time without being able to find the bug. Maybe you’re a bit more clever than I was, and you don’t even need to see the code to know what’s wrong. Do you know what the problem is? Here’s the code, including the refining method that works by placing a new vertex in the center of each triangular face and connecting it to the corners of that face:

//WARNING: THE FOLLOWING CODE CONTAINS A BUG

class UnitSphereApproximation {
public:
    shared_ptr<std::vector<Vector3>> m_vertices;
    shared_ptr<std::vector<uint16>> m_triangles;
    UnitSphereApproximation(shared_ptr<std::vector<Vector3>> vertices, shared_ptr<std::vector<uint16>> triangles)
        : m_vertices(vertices),
          m_triangles(triangles) {
    }
    static UnitSphereApproximation Tetrahedron() {
        // vertices
        auto tetrahedronVertices = std::make_shared<std::vector<Vector3>>();
        auto e = sqrtf(0.5f);
        tetrahedronVertices->push_back(Vector3(+1, 0,-e).normalize());
        tetrahedronVertices->push_back(Vector3(-1, 0,-e).normalize());
        tetrahedronVertices->push_back(Vector3( 0,+1,+e).normalize());
        tetrahedronVertices->push_back(Vector3( 0,-1,+e).normalize());

        // triangles
        auto tetrahedronTriangles = std::make_shared<std::vector<uint16>>();
        tetrahedronTriangles->push_back(0);
        tetrahedronTriangles->push_back(2);
        tetrahedronTriangles->push_back(1);

        tetrahedronTriangles->push_back(0);
        tetrahedronTriangles->push_back(3);
        tetrahedronTriangles->push_back(2);

        tetrahedronTriangles->push_back(0);
        tetrahedronTriangles->push_back(1);
        tetrahedronTriangles->push_back(3);

        tetrahedronTriangles->push_back(3);
        tetrahedronTriangles->push_back(1);
        tetrahedronTriangles->push_back(2);

        return UnitSphereApproximation(tetrahedronVertices, tetrahedronTriangles);
    }
    static UnitSphereApproximation Refine(UnitSphereApproximation sphere) {
        auto refinedVertices = std::make_shared<std::vector<Vector3>>();
        auto refinedTriangles = std::make_shared<std::vector<uint16>>();
        for (auto e : *sphere.m_vertices) { // re-use existing vertices, maintaining indexes
            refinedVertices->push_back(e);
        }
        // for each triangle
        for (auto tri = sphere.m_triangles->begin(); tri != sphere.m_triangles->end(); tri += 3) {
            // compute center vertex of face
            auto ic = refinedVertices->size();
            auto vc = Vector3::ZERO;
            for (auto j = 0; j < 3; j++) {
                vc += (*sphere.m_vertices)[tri[j]];
            }
            refinedVertices->push_back(vc.normalize());
        
            // replace triangle with three triangles surrounding the center vertex
            for (auto j = 0; j < 3; j++) {
                refinedTriangles->push_back(tri[(j + 0) % 3]);
                refinedTriangles->push_back(tri[(j + 1) % 3]);
                refinedTriangles->push_back(ic);
            }
        }
        return UnitSphereApproximation(refinedVertices, refinedTriangles);
    }
    static UnitSphereApproximation WithMinimumPointCount(int minimumPointCount) {
        auto sphere = Tetrahedron();
        while (sphere.m_vertices->size() < minimumPointCount) {
            sphere = Refine(sphere);
        }
        return sphere;
    }
};

(My apologies to veteran C++ programmers, who can probably recognize that I'm used to writing in C# based on the presence of some unnecessary copies and the gratuitous use of smart pointers.)

Don't spend too long scouring the code for misplaced vertices, reversed orderings, or memory corruption. It's not that type of bug.

The Bug

Eventually, I mentioned my predicament to a coworker. Intrigued, they reviewed the code and confirmed that it looked right. However, having lots of experience in 3d graphics, they knew about the graphics diagnostics functionality in visual studio 2012 that lets you see what's being drawn in stages.

One good look at my "spheres", without a shader obscuring the problem, and I knew what was wrong (view the image directly for a larger version):

Fractal tessellation created by refining faces of tetrahedron

The problem is the refining strategy. Not the code implementing the strategy, but the idea itself. In particular, it never makes any edges shorter. It adds shorter and shorter edges, but the original edges are never replaced. This creates a very cool shape, with a volume that does approximate the volume of a sphere, but with a fractal surface instead of a smooth spherical surface. The computer was doing exactly what I said, but what I said did not imply what I really wanted. While looking for the bug I might as well have been trying to figure out why 2+2 didn't return 5 by confirming there were no arithmetic overflows in the code.

Fixing the bug is straightforward: place new vertices at the center of each edge instead of at the center of each face. This ensures edges get shorter which, in combination with ensuring vertices are on the unit sphere, ensures the approximation gets smoother. I could also have continued placing new vertices at the center of faces, with a more intelligent method to choose neighbors, but I expected that to be significantly more tricky to get right. In any case the edge refining strategy works:

Sphere tessellation created by refining edges of tetrahedron

With the benefit of a surface that isn't being fractalled, the bubbles even end up looking like bubbles (note: using a different shader than the shot showing the incorrect bubbles):

Correct bubble spheres

The fixed code for refining a shape is a bit more complicated but, on the other hand, is actually correct. Here it is:

static UnitSphereApproximation Refine(UnitSphereApproximation sphere) {
    auto refinedVertices = std::make_shared<std::vector<Vector3>>();
    // re-use existing vertices, maintaining their indexes
    for (auto e : *sphere.m_vertices) { 
        refinedVertices->push_back(e);
    }
    // place new vertices at centers of spherical edges between existing vertices
    auto edgeVertexMap = std::map<int, int>();
    // for each triangle
    for (auto tri = sphere.m_triangles->begin(); tri != sphere.m_triangles->end(); tri += 3) {
        for (auto j = 0; j < 3; j++) {
            auto i1 = tri[j];
            auto i2 = tri[(j + 1) % 3];
            if (i1 >= i2) continue; // avoid adding the same edge vertex twice (once from X to Y and once from Y and X)
            auto undirectedEdgeId = i1*sphere.m_vertices->size() + i2;
            
            edgeVertexMap[undirectedEdgeId] = refinedVertices->size();
            auto edgeCenter = ((*sphere.m_vertices)[i1] + (*sphere.m_vertices)[i2]).normalize();
            refinedVertices->push_back(edgeCenter);
        }
    }

    auto refinedTriangles = std::make_shared<std::vector<uint16>>();
    // for each triangle: create new triangles, using old and new vertices
    for (auto tri = sphere.m_triangles->begin(); tri != sphere.m_triangles->end(); tri += 3) {
        // find vertices at the center of each of the triangle's edges
        int edgeCenterVertices[3];
        for (auto j = 0; j < 3; j++) {
            auto i1 = tri[j];
            auto i2 = tri[(j + 1) % 3];
            auto undirectedEdgeId = min(i1, i2)*sphere.m_vertices->size() + max(i1, i2);
            edgeCenterVertices[j] = edgeVertexMap[undirectedEdgeId];
        }
    
        // create a triangle covering the center, touching the three edges
        refinedTriangles->push_back(edgeCenterVertices[0]);
        refinedTriangles->push_back(edgeCenterVertices[1]);
        refinedTriangles->push_back(edgeCenterVertices[2]);
        // create a triangle for each corner of the existing triangle
        for (auto j = 0; j < 3; j++) {
            refinedTriangles->push_back(tri[j]);
            refinedTriangles->push_back(edgeCenterVertices[(j + 0) % 3]);
            refinedTriangles->push_back(edgeCenterVertices[(j + 2) % 3]);
        }
    }
    return UnitSphereApproximation(refinedVertices, refinedTriangles);
}

(Again, my apologies to veteran C++-ers.)

Summary

This bug was particularly subtle to me, which is a bit embarrassing considering how obvious it is in hindsight. In my defense the bug is not only hiding in the spec, instead of the code, but it only appears after the second iteration of refinement. I simulated the first iteration with pen and paper, but the result is just a cube (a perfectly fine low-quality approximation of a sphere). Also, it's hard to feel bad about making a stupid mistake when you get a neat looking fractal out of it.

Extra: Icosahedron

One thing I noticed, while making the animation of a tetrahedron being refined, is that the vertices of the original tetrahedron remain quite obvious. I decided to try starting from an icosahedron:

Sphere tessellation created by refining edges of icosahedron

This looks better initially but after shading my eyes can't see a difference, and we ended up just using the tetrahedron. In any case, here's the code I wrote to generate an icosahedron:

UnitSphereApproximation UnitSphereApproximation::Icosahedron() {
    auto vertices = std::make_shared<std::vector<Vector3>>();
    auto e = 0.5f + sqrtf(1.25f);
    // create vertices (4 per axis plane)
    for (auto d = 0; d < 3; d++) {
        for (auto s1 = -1; s1 <= +1; s1 += 2) {
            for (auto s2 = -1; s2 <= +1; s2 += 2) {
                auto v = Vector3::ZERO;
                v[(d + 1) % 3] = e * s1;
                v[(d + 2) % 3] = 1 * s2;
                vertices->push_back(v.normalize());
            }
        }
    }

    auto triangles = std::make_shared<std::vector<uint16>>();
    // create the triangles that have a point on each of the three axis planes
    auto id = [](int d, int s1, int s2){ return d * 4 + (s1 + 1) + ((s2 + 1) >> 1); };
    for (auto s1 = -1; s1 <= +1; s1 += 2) {
        for (auto s2 = -1; s2 <= +1; s2 += 2) {
            for (auto s3 = -1; s3 <= +1; s3 += 2) {
                auto rev = s1*s2*s3 == -1;
                auto i1 = id(0,s1,s2);
                auto i2 = id(1,s2,s3);
                auto i3 = id(2,s3,s1);
                triangles->push_back(i1);
                triangles->push_back(rev ? i3 : i2);
                triangles->push_back(rev ? i2 : i3);
            }
        }
    }
    // create the triangles that have two points on one axis plane and one point on another
    for (auto d = 0; d < 3; d++) {
        for (auto s1 = -1; s1 <= +1; s1 += 2) {
            for (auto s2 = -1; s2 <= +1; s2 += 2) {
                auto rev = s1*s2 == +1;
                auto i2 = id(d,s1,-1);
                auto i1 = id(d,s1,+1);
                auto i3 = id((d + 2) % 3,s2,s1);
                triangles->push_back(i1);
                triangles->push_back(rev ? i3 : i2);
                triangles->push_back(rev ? i2 : i3);
            }
        }
    }
        
    return UnitSphereApproximation(vertices, triangles);
}

Maybe someone will find that useful or at least interesting.

---

Discuss on Reddit

---


Twisted Oak Studios offers consulting and development on high-tech interactive projects. Check out our portfolio, or Give us a shout if you have anything you think some really rad engineers should help you with.

Archive

More interesting posts (18 of 42 articles)

Or check out our Portfolio.