6491 and 7491 PROJECTS
Alexander Powell

Surface Subdivision

CS6491: Computer Graphics
Alexander Powell - [ alex at powelltown dot com ]
October 2, 2003

This project is an implementation of the Loop, Butterfly, and Jarek's subdivision methods. I implemented the project on top of the TMesh class in my group's Project 1 source. The particular method which performs the subdivision is as follows:


void TMesh::subdivide(int mode)
{
    int numCrn = numTriangles * 3;
    int newNumVertices = (numVertices + (numCrn / 2) + 1);
    
    cout << "Old number of vertices: " << numVertices << endl;
    
	GLfloat time = glutGet(GLUT_ELAPSED_TIME);
    
    int m[numCrn];
    bool tagged[numVertices];
    Vector displacement[numVertices];
    for (int ii = 0; ii < numCrn; ii++)
        m[ii] = -1;
    for (ii = 0; ii < numVertices; ii++)
        tagged[ii] = false;
    
    std::deque<Vector> newVertices(numVertices);
    for (int c = 0; c < numCrn; c++)
    {
        //printf("Computing new vertices... %f\r", (double)jj / (numTriangles * 3));

        if ((mode == LOOP) || (mode == JAREK))
        {
            if (tagged[cornerTable[c]] == false)
            {
                tagged[cornerTable[c]] = true;
                std::deque<int> neighbors;
                
                int cur = c;
                do {
                    neighbors.push_back(n(cur));
                    cur = n(o(n(cur)));
                } while (cur != c);
                
                float numNeighbors = neighbors.size();
                Vector sum = Vector(0, 0, 0);
                for (int ii = 0; ii < numNeighbors; ii++)
                    sum += v(neighbors[ii]);
                
                if (mode == JAREK)
                {
                    sum /= numNeighbors;
                    displacement[cornerTable[c]] = (sum - v(c)) * 0.125f;
                    newVertices[cornerTable[c]] = v(c) + displacement[cornerTable[c]];
                }
                else // if (mode == LOOP)
                {
                    float pp = cos((2.0f * M_PI) / numNeighbors);
                    float alpha = (1.0f / numNeighbors) * (0.625f - (0.375f - 0.25f * pp * pp));
                    
                    newVertices[cornerTable[c]] = ((1.0f - numNeighbors * alpha) * v(c) + alpha * sum);
                }
            }
        }

        if (m[c] == -1)
        {
            Vector newV;
            
            if ((mode == MIDPOINT) || (mode == BUTTERFLY))
            {
                Vector cv = v(c);
                Vector cnv = v(n(c));
                Vector cpv = v(p(c));
                Vector cov = v(o(c));
                Vector clv = v(l(c));
                Vector crv = v(r(c));
                Vector colv = v(l(o(c)));
                Vector corv = v(r(o(c)));
                
                
                if (mode == MIDPOINT)
                    newV = (cnv + cpv) * 0.5f;
                else if (mode == BUTTERFLY)
                    newV = (cnv + cpv) * 0.5f + (cv + cov) * 0.125f - (clv + crv + colv + corv) * 0.0625f;
            }
            else if (mode == LOOP)
            {
                Vector cv = v(c);
                Vector cnv = v(n(c));
                Vector cpv = v(p(c));
                Vector cov = v(o(c));
                
                newV = (cnv + cpv) * 0.375f + (cv + cov) * 0.125f;
            }
            else if (mode == JAREK)
            {
                newV = (v(n(c)) + v(p(c))) * 0.5f - (displacement[cornerTable[n(c)]] + displacement[cornerTable[p(c)]]) * 0.5f;
            }
            
            int newVertexIdx = addVertex(newV);
            
            m[c] = newVertexIdx;
            m[o(c)] = newVertexIdx;
        }
    }
    
    if (mode == LOOP)
    {
        for (int ii = 0; ii < newVertices.size(); ii++)
        {
            vt[ii] = newVertices[ii];
        }
    }
    
    printf("New number of vertices: %d\n", numVertices);
    
	cout<<glutGet(GLUT_ELAPSED_TIME) - time<<" ms to calculate new vertices"<<endl;
	time = glutGet(GLUT_ELAPSED_TIME);
    
    // Create space for our new triangles
    int *newCorners = new int[numTriangles * 12];
       
    int cnt = 0;

    printf("Old number of triangles: %d\n", numTriangles);

    for (int jj = 0; jj < numCrn; jj = jj + 3)
    {
        int c = jj;
        
        //printf("Computing new triangles... %f\r", (double)jj / numCrn);
        
        // Triangle 1
        newCorners[cnt] = cornerTable[c];
        cnt++;
        newCorners[cnt] = m[p(c)];
        cnt++;
        newCorners[cnt] = m[n(c)];
        cnt++;
        
        // Triangle 2
        newCorners[cnt] = m[p(c)];
        cnt++;
        newCorners[cnt] = cornerTable[n(c)];
        cnt++;
        newCorners[cnt] = m[c];
        cnt++;
        
        // Triangle 3
        newCorners[cnt] = m[n(c)];
        cnt++;
        newCorners[cnt] = m[c];
        cnt++;
        newCorners[cnt] = cornerTable[p(c)];
        cnt++;
        
        // Triangle 4
        newCorners[cnt] = m[c];
        cnt++;
        newCorners[cnt] = m[n(c)];
        cnt++;
        newCorners[cnt] = m[p(c)];
        cnt++;
    }
    
	cout<<glutGet(GLUT_ELAPSED_TIME) - time<<" ms to calculate new triangles"<<endl;
	time = glutGet(GLUT_ELAPSED_TIME);
    
    delete []cornerTable;
    cornerTable = newCorners;
    numTriangles *= 4;
    
    printf("New number of triangles: %d %d\n", cnt / 3, numTriangles);
    
	cout<<glutGet(GLUT_ELAPSED_TIME) - time<<" ms to clean up after subdivision"<<endl;
    
	time = glutGet(GLUT_ELAPSED_TIME);
    createOppositeTable();
	cout<<glutGet(GLUT_ELAPSED_TIME) - time<<" ms to recalculate opposites"<<endl;
    
	time = glutGet(GLUT_ELAPSED_TIME);
    calculateNormals();
	cout<<glutGet(GLUT_ELAPSED_TIME) - time<<" ms to recalculate normals"<<endl;

	time = glutGet(GLUT_ELAPSED_TIME);
	colorShells();
	cout<<glutGet(GLUT_ELAPSED_TIME) - time<<" ms to recolor the shells"<<endl;
    }
}


The method can be called with modes MIDPOINT (which only subdivides into strict midpoints, i.e. it doesn't change the surface of the mesh at all-- just the triangle count), BUTTERFLY, LOOP, or JAREK. The particular mode you choose doesn't change the triangle splitting code at all-- it just changes the vertices of the final mesh.

The full source can be retrieved here: P2.tar.gz

Pictures

Click on a picture for a bigger (but likely of worse quality) version.





squares.csg - From left to right, top to bottom: Original isosurface
1 iteration of Butterfly :: 2 iterations of Butterfly :: 1 iteration of Loop :: 2 iterations of Loop
1 iteration of Jarek's :: 2 iterations of Jarek's :: 2 iterations of Jarek's without the wireframe, demonstrating the oscillation resulting in a "rougher" surface.

chair.csg - From left to right: Original isosurface :: 2 iterations of Butterfly :: 2 iterations of Loop (note the smoother joints on the Loop subdivision compared with the "creases" that are still evident in Butterfly).

spheres.csg - From left to right: Original isosurface :: 1 iteration of Butterfly :: 1 iteration of Loop :: A demonstration of the combination of Taubin's surface fairing (10 iterations) and Butterfly subdivision (1 iteration). The combination of surface fairing and subdivision provides a pleasing shape which retains the overall size of the original and is computationally faster than just using subdivision to get to such a smooth surface. My opinion is that this combination provides the best results using the smallest amount of CPU time.

Running times

The subdivision methods ran rather slow on my laptop. I found that the biggest problem was with memory reallocation, and since this particular portion of the project was not heavily optimized, there is probably *lots* of room for improvement in these times:

  • Subdivision of ~300 triangles: 4.3 seconds, average
  • Subdivision of ~1000 triangles: 22.7 seconds, average
  • Subdivision of ~3000 triangles: 74.2 seconds, average

Note that I didn't separate the times into the categories of which subdivision technique was used-- all of the subdivision times were bounded mostly by the memory allocation/reallocation, and therefore all of them ran equally as fast/slow, within only 50-100 ms.

about me
navigate
links