Surface Fairing

void TMesh::fairSurface(GLfloat lm, WeightType wt)
{
    int ii;
    deque<Vector> newVertices(numVertices);
    deque<bool> tagged(numVertices);

    for (ii = 0; ii < numVertices; ii++)
        tagged[ii] = false;

    // We loop through each corner
    for(unsigned int a = 0; a < numTriangles * 3; a++)
    {
        if(!tagged[cornerTable[a]])
        {
            // Create lists to keep track of our neighbors and their weights
            deque<Vector> neighborhood;
            deque<GLfloat> weights;
            
            // Visit the neighbors to calculate their weights
            int start = n(a);
            int current = start;
            do {
                // Add this vertex to our neighborhood list
                neighborhood.push_back(vt[cornerTable[current]]);
                
                // Calculate the weight for this edge
                GLfloat curWeight;
                if (wt == DISTANCE)
                    curWeight = ((vt[cornerTable[a]] - vt[cornerTable[current]]).length());
                else
                {
                    Vector v1 = vt[cornerTable[a]] - vt[cornerTable[n(a)]];
                    Vector v2 = vt[cornerTable[a]] - vt[cornerTable[p(a)]];
                    Vector v3 = vt[cornerTable[o(n(current))]] - vt[cornerTable[n(o(n(current)))]];
                    Vector v4 = vt[cornerTable[o(n(current))]] - vt[cornerTable[p(o(n(current)))]];
                    curWeight = (v1.cross(v2).length() + v3.cross(v4).length());
                }
                
                weights.push_back(curWeight);
                current = o(n(current));
            } while (current != start);

            // Normalize weights
            GLfloat totalWeight = 0.0;
            for (ii = 0; ii < weights.size(); ii++)
                totalWeight += weights[ii];
            for (ii = 0; ii < weights.size(); ii++)
                weights[ii] = weights[ii] / totalWeight;

            // Compute new vertex position based on the Gaussian smoothing step in Eq. 4
            Vector dx;
            for (ii = 0; ii < neighborhood.size(); ii++)
                dx += weights[ii] * (neighborhood[ii] - vt[cornerTable[a]]);
            newVertices[cornerTable[a]] = vt[cornerTable[a]] + dx * lm;
            
            // Tag this vertex so we don't recompute its location
            tagged[cornerTable[a]] = true;
        }
    }

    //replace old vertices
    for(unsigned int a = 0; a < newVertices.size(); a++)
        vt[a] = newVertices[a];
}
}

GLfloat lambda = 0.4;
GLfloat Kpb = 0.01;
GLfloat mu = 1.0 / (Kpb - (1.0 / lambda));
for (int ii = 0; ii < 100; ii++)
{
    tmesh->fairSurface(lambda, AREA);
    tmesh->fairSurface(mu, AREA);
}
tmesh->calculateNormals();