• Create Account

### #Actualtaby

Posted 24 April 2013 - 09:36 AM

Nice work! This could surely come in handy at some point.

Thank you. I hope it comes it handy.

The paper 'Geometric Signal Processing on Polygonal Meshes' by G. Taubin mentions a few edge weight schemes:

- constant unit weight: w_ij = 1

- inverse edge length weight: w_ij = 1 / |e_ij|

- curvature normal-esque (cotan) weight: w_ij = 1 / tan(theta_ij0) + 1 / tan(theta_ij1)

The code in the zip file from the first post uses the constant unit weight scheme (see indexed_mesh::laplace_smooth()).

I didn't put the other schemes into the code zip file because I'm not sure if I'm implementing them correctly, but I'll include them here in case anyone wishes to play with them to see if/where I'm going wrong.

Here's the code to replace the laplace_smooth function with the inverse edge length weight scheme:

void indexed_mesh::inverse_edge_length_smooth(const float scale)
{
vector<vertex_3> displacements(vertices.size(), vertex_3(0, 0, 0));

// Get per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
{
// Skip rogue vertices (which were probably made rogue during a previous
// attempt to fix mesh cracks).
if(0 == vertex_to_vertex_indices[i].size())
continue;

vector<float> weights(vertex_to_vertex_indices[i].size(), 0.0f);

// Calculate weights based on inverse edge lengths.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];

float edge_length = vertices[i].distance(vertices[neighbour_j]);

if(0 == edge_length)
edge_length = numeric_limits<float>::epsilon();

weights[j] = 1.0f / edge_length;
}

// Normalize the weights so that they sum up to 1.
float s = 0;

for(size_t j = 0; j < weights.size(); j++)
s += weights[j];

if(0 == s)
s = numeric_limits<float>::epsilon();

for(size_t j = 0; j < weights.size(); j++)
weights[j] /= s;

// Sum the displacements.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];
displacements[i] += (vertices[neighbour_j] - vertices[i])*weights[j];
}
}

// Apply per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
vertices[i] += displacements[i]*scale;
}

Here's the code to replace the laplace_smooth function with the curvature normal-esque weight scheme (it works most times, but no guarantees):

void indexed_mesh::curvature_normal_smooth(const float scale)
{
vector<vertex_3> displacements(vertices.size(), vertex_3(0, 0, 0));

// Get per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
{
if(0 == vertex_to_vertex_indices[i].size())
continue;

vector<float> weights(vertex_to_vertex_indices[i].size(), 0.0f);

size_t angle_error = 0;

// For each vertex pair (ie. each edge),
// calculate weight based on the two opposing angles (ie. curvature normal scheme).
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t angle_count = 0;

size_t neighbour_j = vertex_to_vertex_indices[i][j];

// Find out which two triangles are shared by the edge.
for(size_t k = 0; k < vertex_to_triangle_indices[i].size(); k++)
{
for(size_t l = 0; l < vertex_to_triangle_indices[neighbour_j].size(); l++)
{
size_t tri0_index = vertex_to_triangle_indices[i][k];
size_t tri1_index = vertex_to_triangle_indices[neighbour_j][l];

// This will occur twice per edge.
if(tri0_index == tri1_index)
{
// Find the third vertex in this triangle (the vertex that doesn't belong to the edge).
for(size_t m = 0; m < 3; m++)
{
// This will occur once per triangle.
if(triangles[tri0_index].vertex_indices[m] != i && triangles[tri0_index].vertex_indices[m] != neighbour_j)
{
size_t opp_vert_index = triangles[tri0_index].vertex_indices[m];

// Get the angle opposite of the edge.
vertex_3 a = vertices[i] - vertices[opp_vert_index];
vertex_3 b = vertices[neighbour_j] - vertices[opp_vert_index];
a.normalize();
b.normalize();

float dot = a.dot(b);

if(-1 > dot)
dot = -1;
else if(1 < dot)
dot = 1;

float angle = acosf(dot);

// Curvature normal weighting.
float slope = tanf(angle);

if(0 == slope)
slope = numeric_limits<float>::epsilon();

// Note: Some weights will be negative, due to obtuse triangles.
// You may wish to do weights[j] += fabsf(1.0f / slope); here.
weights[j] += 1.0f / slope;

angle_count++;

break;
}
}

// Since we found a triangle match, we can skip to the first vertex's next triangle.
break;
}
}
} // End of: Find out which two triangles are shared by the vertex pair.

if(angle_count != 2)
angle_error++;

} // End of: For each vertex pair (ie. each edge).

if(angle_error != 0)
{
cout << "Warning: Vertex " << i << " belongs to " << angle_error << " edges that do not belong to two triangles (" << vertex_to_vertex_indices[i].size() - angle_error << " edges were OK)." << endl;
cout << "Your mesh probably has cracks or holes in it." << endl;
}

// Normalize the weights so that they sum up to 1.
float s = 0;

// Note: Some weights will be negative, due to obtuse triangles.
// You may wish to do s += fabsf(weights[j]); here.
for(size_t j = 0; j < weights.size(); j++)
s += weights[j];

if(0 == s)
s = numeric_limits<float>::epsilon();

for(size_t j = 0; j < weights.size(); j++)
weights[j] /= s;

// Sum the displacements.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];

displacements[i] += (vertices[neighbour_j] - vertices[i])*weights[j];
}
}

// To do: Find out why there are cases where displacement is much, much, much larger than all edge lengths put together.

// Apply per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
vertices[i] += displacements[i]*scale;
}

### #9taby

Posted 24 April 2013 - 09:35 AM

Nice work! This could surely come in handy at some point.

Thank you. I hope it comes it handy.

The paper 'Geometric Signal Processing on Polygonal Meshes' by G. Taubin mentions a few edge weight schemes:

- constant unit weight: w_ij = 1

- inverse edge length weight: w_ij = 1 / |e_ij|

- curvature normal-esque (cotan) weight: w_ij = 1 / tan(theta_ij0) + 1 / tan(theta_ij1)

The code in the zip file from the first post uses the constant unit weight scheme (see indexed_mesh::laplace_smooth()).

I didn't put the other schemes into the code zip file because I'm not sure if I'm implementing them correctly, but I'll include them here in case anyone wishes to play with them to see if/where I'm going wrong.

Here's the code to replace the laplace_smooth function with the inverse edge length weight scheme:

void indexed_mesh::inverse_edge_length_smooth(const float scale)
{
vector<vertex_3> displacements(vertices.size(), vertex_3(0, 0, 0));

// Get per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
{
// Skip rogue vertices (which were probably made rogue during a previous
// attempt to fix mesh cracks).
if(0 == vertex_to_vertex_indices[i].size())
continue;

vector<float> weights(vertex_to_vertex_indices[i].size(), 0.0f);

// Calculate weights based on inverse edge lengths.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];

float edge_length = vertices[i].distance(vertices[neighbour_j]);

if(0 == edge_length)
edge_length = numeric_limits<float>::epsilon();

weights[j] = 1.0f / edge_length;
}

// Normalize the weights so that they sum up to 1.
float s = 0;

for(size_t j = 0; j < weights.size(); j++)
s += weights[j];

if(0 == s)
s = numeric_limits<float>::epsilon();

for(size_t j = 0; j < weights.size(); j++)
weights[j] /= s;

// Sum the displacements.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];
displacements[i] += (vertices[neighbour_j] - vertices[i])*weights[j];
}
}

// Apply per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
vertices[i] += displacements[i]*scale;
}

Here's the code to replace the laplace_smooth function with the curvature normal-esque weight scheme (it works most times, but no guarantees:

void indexed_mesh::curvature_normal_smooth(const float scale)
{
vector<vertex_3> displacements(vertices.size(), vertex_3(0, 0, 0));

// Get per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
{
if(0 == vertex_to_vertex_indices[i].size())
continue;

vector<float> weights(vertex_to_vertex_indices[i].size(), 0.0f);

size_t angle_error = 0;

// For each vertex pair (ie. each edge),
// calculate weight based on the two opposing angles (ie. curvature normal scheme).
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t angle_count = 0;

size_t neighbour_j = vertex_to_vertex_indices[i][j];

// Find out which two triangles are shared by the edge.
for(size_t k = 0; k < vertex_to_triangle_indices[i].size(); k++)
{
for(size_t l = 0; l < vertex_to_triangle_indices[neighbour_j].size(); l++)
{
size_t tri0_index = vertex_to_triangle_indices[i][k];
size_t tri1_index = vertex_to_triangle_indices[neighbour_j][l];

// This will occur twice per edge.
if(tri0_index == tri1_index)
{
// Find the third vertex in this triangle (the vertex that doesn't belong to the edge).
for(size_t m = 0; m < 3; m++)
{
// This will occur once per triangle.
if(triangles[tri0_index].vertex_indices[m] != i && triangles[tri0_index].vertex_indices[m] != neighbour_j)
{
size_t opp_vert_index = triangles[tri0_index].vertex_indices[m];

// Get the angle opposite of the edge.
vertex_3 a = vertices[i] - vertices[opp_vert_index];
vertex_3 b = vertices[neighbour_j] - vertices[opp_vert_index];
a.normalize();
b.normalize();

float dot = a.dot(b);

if(-1 > dot)
dot = -1;
else if(1 < dot)
dot = 1;

float angle = acosf(dot);

// Curvature normal weighting.
float slope = tanf(angle);

if(0 == slope)
slope = numeric_limits<float>::epsilon();

// Note: Some weights will be negative, due to obtuse triangles.
// You may wish to do weights[j] += fabsf(1.0f / slope); here.
weights[j] += 1.0f / slope;

angle_count++;

break;
}
}

// Since we found a triangle match, we can skip to the first vertex's next triangle.
break;
}
}
} // End of: Find out which two triangles are shared by the vertex pair.

if(angle_count != 2)
angle_error++;

} // End of: For each vertex pair (ie. each edge).

if(angle_error != 0)
{
cout << "Warning: Vertex " << i << " belongs to " << angle_error << " edges that do not belong to two triangles (" << vertex_to_vertex_indices[i].size() - angle_error << " edges were OK)." << endl;
cout << "Your mesh probably has cracks or holes in it." << endl;
}

// Normalize the weights so that they sum up to 1.
float s = 0;

// Note: Some weights will be negative, due to obtuse triangles.
// You may wish to do s += fabsf(weights[j]); here.
for(size_t j = 0; j < weights.size(); j++)
s += weights[j];

if(0 == s)
s = numeric_limits<float>::epsilon();

for(size_t j = 0; j < weights.size(); j++)
weights[j] /= s;

// Sum the displacements.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];

displacements[i] += (vertices[neighbour_j] - vertices[i])*weights[j];
}
}

// To do: Find out why there are cases where displacement is much, much, much larger than all edge lengths put together.

// Apply per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
vertices[i] += displacements[i]*scale;
}

### #8taby

Posted 24 April 2013 - 09:35 AM

Nice work! This could surely come in handy at some point.

Thank you. I hope it comes it handy.

The paper 'Geometric Signal Processing on Polygonal Meshes' by G. Taubin mentions a few edge weight schemes:

- constant unit weight: w_ij = 1

- inverse edge length weight: w_ij = 1 / |e_ij|

- curvature normal-esque (cotan) weight: w_ij = 1 / tan(theta_ij0) + 1 / tan(theta_ij1)

The code in the first post uses the constant unit weight scheme (see indexed_mesh::laplace_smooth()).

I didn't put the other schemes into the code zip file because I'm not sure if I'm implementing them correctly, but I'll include them here in case anyone wishes to play with them to see if/where I'm going wrong.

Here's the code to replace the laplace_smooth function with the inverse edge length weight scheme:

void indexed_mesh::inverse_edge_length_smooth(const float scale)
{
vector<vertex_3> displacements(vertices.size(), vertex_3(0, 0, 0));

// Get per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
{
// Skip rogue vertices (which were probably made rogue during a previous
// attempt to fix mesh cracks).
if(0 == vertex_to_vertex_indices[i].size())
continue;

vector<float> weights(vertex_to_vertex_indices[i].size(), 0.0f);

// Calculate weights based on inverse edge lengths.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];

float edge_length = vertices[i].distance(vertices[neighbour_j]);

if(0 == edge_length)
edge_length = numeric_limits<float>::epsilon();

weights[j] = 1.0f / edge_length;
}

// Normalize the weights so that they sum up to 1.
float s = 0;

for(size_t j = 0; j < weights.size(); j++)
s += weights[j];

if(0 == s)
s = numeric_limits<float>::epsilon();

for(size_t j = 0; j < weights.size(); j++)
weights[j] /= s;

// Sum the displacements.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];
displacements[i] += (vertices[neighbour_j] - vertices[i])*weights[j];
}
}

// Apply per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
vertices[i] += displacements[i]*scale;
}

Here's the code to replace the laplace_smooth function with the curvature normal-esque weight scheme (it works most times, but no guarantees:

void indexed_mesh::curvature_normal_smooth(const float scale)
{
vector<vertex_3> displacements(vertices.size(), vertex_3(0, 0, 0));

// Get per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
{
if(0 == vertex_to_vertex_indices[i].size())
continue;

vector<float> weights(vertex_to_vertex_indices[i].size(), 0.0f);

size_t angle_error = 0;

// For each vertex pair (ie. each edge),
// calculate weight based on the two opposing angles (ie. curvature normal scheme).
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t angle_count = 0;

size_t neighbour_j = vertex_to_vertex_indices[i][j];

// Find out which two triangles are shared by the edge.
for(size_t k = 0; k < vertex_to_triangle_indices[i].size(); k++)
{
for(size_t l = 0; l < vertex_to_triangle_indices[neighbour_j].size(); l++)
{
size_t tri0_index = vertex_to_triangle_indices[i][k];
size_t tri1_index = vertex_to_triangle_indices[neighbour_j][l];

// This will occur twice per edge.
if(tri0_index == tri1_index)
{
// Find the third vertex in this triangle (the vertex that doesn't belong to the edge).
for(size_t m = 0; m < 3; m++)
{
// This will occur once per triangle.
if(triangles[tri0_index].vertex_indices[m] != i && triangles[tri0_index].vertex_indices[m] != neighbour_j)
{
size_t opp_vert_index = triangles[tri0_index].vertex_indices[m];

// Get the angle opposite of the edge.
vertex_3 a = vertices[i] - vertices[opp_vert_index];
vertex_3 b = vertices[neighbour_j] - vertices[opp_vert_index];
a.normalize();
b.normalize();

float dot = a.dot(b);

if(-1 > dot)
dot = -1;
else if(1 < dot)
dot = 1;

float angle = acosf(dot);

// Curvature normal weighting.
float slope = tanf(angle);

if(0 == slope)
slope = numeric_limits<float>::epsilon();

// Note: Some weights will be negative, due to obtuse triangles.
// You may wish to do weights[j] += fabsf(1.0f / slope); here.
weights[j] += 1.0f / slope;

angle_count++;

break;
}
}

// Since we found a triangle match, we can skip to the first vertex's next triangle.
break;
}
}
} // End of: Find out which two triangles are shared by the vertex pair.

if(angle_count != 2)
angle_error++;

} // End of: For each vertex pair (ie. each edge).

if(angle_error != 0)
{
cout << "Warning: Vertex " << i << " belongs to " << angle_error << " edges that do not belong to two triangles (" << vertex_to_vertex_indices[i].size() - angle_error << " edges were OK)." << endl;
cout << "Your mesh probably has cracks or holes in it." << endl;
}

// Normalize the weights so that they sum up to 1.
float s = 0;

// Note: Some weights will be negative, due to obtuse triangles.
// You may wish to do s += fabsf(weights[j]); here.
for(size_t j = 0; j < weights.size(); j++)
s += weights[j];

if(0 == s)
s = numeric_limits<float>::epsilon();

for(size_t j = 0; j < weights.size(); j++)
weights[j] /= s;

// Sum the displacements.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];

displacements[i] += (vertices[neighbour_j] - vertices[i])*weights[j];
}
}

// To do: Find out why there are cases where displacement is much, much, much larger than all edge lengths put together.

// Apply per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
vertices[i] += displacements[i]*scale;
}

### #7taby

Posted 24 April 2013 - 09:34 AM

Nice work! This could surely come in handy at some point.

Thank you. I hope it comes it handy.

The paper 'Geometric Signal Processing on Polygonal Meshes' by G. Taubin mentions a few edge weight schemes:

- constant unit weight: w_ij = 1

- inverse edge length weight: w_ij = 1 / |e_ij|

- curvature normal-esque (cotan) weight: w_ij = 1 / tan(theta_ij0) + 1 / tan(theta_ij1)

The code in the first post uses the constant unit weight scheme (see indexed_mesh::laplace_smooth()). I didn't put the other schemes into the code zip file because I'm not sure if I'm implementing them correctly, but I'll include them here in case anyone wishes to play with them to see if/where I'm going wrong.

Here's the code to replace the laplace_smooth function with the inverse edge length weight scheme:

void indexed_mesh::inverse_edge_length_smooth(const float scale)
{
vector<vertex_3> displacements(vertices.size(), vertex_3(0, 0, 0));

// Get per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
{
// Skip rogue vertices (which were probably made rogue during a previous
// attempt to fix mesh cracks).
if(0 == vertex_to_vertex_indices[i].size())
continue;

vector<float> weights(vertex_to_vertex_indices[i].size(), 0.0f);

// Calculate weights based on inverse edge lengths.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];

float edge_length = vertices[i].distance(vertices[neighbour_j]);

if(0 == edge_length)
edge_length = numeric_limits<float>::epsilon();

weights[j] = 1.0f / edge_length;
}

// Normalize the weights so that they sum up to 1.
float s = 0;

for(size_t j = 0; j < weights.size(); j++)
s += weights[j];

if(0 == s)
s = numeric_limits<float>::epsilon();

for(size_t j = 0; j < weights.size(); j++)
weights[j] /= s;

// Sum the displacements.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];
displacements[i] += (vertices[neighbour_j] - vertices[i])*weights[j];
}
}

// Apply per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
vertices[i] += displacements[i]*scale;
}

Here's the code to replace the laplace_smooth function with the curvature normal-esque weight scheme (it works most times, but no guarantees:

void indexed_mesh::curvature_normal_smooth(const float scale)
{
vector<vertex_3> displacements(vertices.size(), vertex_3(0, 0, 0));

// Get per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
{
if(0 == vertex_to_vertex_indices[i].size())
continue;

vector<float> weights(vertex_to_vertex_indices[i].size(), 0.0f);

size_t angle_error = 0;

// For each vertex pair (ie. each edge),
// calculate weight based on the two opposing angles (ie. curvature normal scheme).
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t angle_count = 0;

size_t neighbour_j = vertex_to_vertex_indices[i][j];

// Find out which two triangles are shared by the edge.
for(size_t k = 0; k < vertex_to_triangle_indices[i].size(); k++)
{
for(size_t l = 0; l < vertex_to_triangle_indices[neighbour_j].size(); l++)
{
size_t tri0_index = vertex_to_triangle_indices[i][k];
size_t tri1_index = vertex_to_triangle_indices[neighbour_j][l];

// This will occur twice per edge.
if(tri0_index == tri1_index)
{
// Find the third vertex in this triangle (the vertex that doesn't belong to the edge).
for(size_t m = 0; m < 3; m++)
{
// This will occur once per triangle.
if(triangles[tri0_index].vertex_indices[m] != i && triangles[tri0_index].vertex_indices[m] != neighbour_j)
{
size_t opp_vert_index = triangles[tri0_index].vertex_indices[m];

// Get the angle opposite of the edge.
vertex_3 a = vertices[i] - vertices[opp_vert_index];
vertex_3 b = vertices[neighbour_j] - vertices[opp_vert_index];
a.normalize();
b.normalize();

float dot = a.dot(b);

if(-1 > dot)
dot = -1;
else if(1 < dot)
dot = 1;

float angle = acosf(dot);

// Curvature normal weighting.
float slope = tanf(angle);

if(0 == slope)
slope = numeric_limits<float>::epsilon();

// Note: Some weights will be negative, due to obtuse triangles.
// You may wish to do weights[j] += fabsf(1.0f / slope); here.
weights[j] += 1.0f / slope;

angle_count++;

break;
}
}

// Since we found a triangle match, we can skip to the first vertex's next triangle.
break;
}
}
} // End of: Find out which two triangles are shared by the vertex pair.

if(angle_count != 2)
angle_error++;

} // End of: For each vertex pair (ie. each edge).

if(angle_error != 0)
{
cout << "Warning: Vertex " << i << " belongs to " << angle_error << " edges that do not belong to two triangles (" << vertex_to_vertex_indices[i].size() - angle_error << " edges were OK)." << endl;
cout << "Your mesh probably has cracks or holes in it." << endl;
}

// Normalize the weights so that they sum up to 1.
float s = 0;

// Note: Some weights will be negative, due to obtuse triangles.
// You may wish to do s += fabsf(weights[j]); here.
for(size_t j = 0; j < weights.size(); j++)
s += weights[j];

if(0 == s)
s = numeric_limits<float>::epsilon();

for(size_t j = 0; j < weights.size(); j++)
weights[j] /= s;

// Sum the displacements.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];

displacements[i] += (vertices[neighbour_j] - vertices[i])*weights[j];
}
}

// To do: Find out why there are cases where displacement is much, much, much larger than all edge lengths put together.

// Apply per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
vertices[i] += displacements[i]*scale;
}

### #6taby

Posted 24 April 2013 - 09:34 AM

Nice work! This could surely come in handy at some point.

Thank you. I hope it comes it handy.

The paper 'Geometric Signal Processing on Polygonal Meshes' by G. Taubin mentions a few edge weight schemes:

- constant unit weight: w_ij = 1

- inverse edge length weight: w_ij = 1/|e_ij|

- curvature normal-esque (cotan) weight: w_ij = 1/tan(theta_ij0) + 1/tan(theta_ij1)

The code in the first post uses the constant unit weight scheme (see indexed_mesh::laplace_smooth()). I didn't put the other schemes into the code zip file because I'm not sure if I'm implementing them correctly, but I'll include them here in case anyone wishes to play with them to see if/where I'm going wrong.

Here's the code to replace the laplace_smooth function with the inverse edge length weight scheme:

void indexed_mesh::inverse_edge_length_smooth(const float scale)
{
vector<vertex_3> displacements(vertices.size(), vertex_3(0, 0, 0));

// Get per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
{
// Skip rogue vertices (which were probably made rogue during a previous
// attempt to fix mesh cracks).
if(0 == vertex_to_vertex_indices[i].size())
continue;

vector<float> weights(vertex_to_vertex_indices[i].size(), 0.0f);

// Calculate weights based on inverse edge lengths.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];

float edge_length = vertices[i].distance(vertices[neighbour_j]);

if(0 == edge_length)
edge_length = numeric_limits<float>::epsilon();

weights[j] = 1.0f / edge_length;
}

// Normalize the weights so that they sum up to 1.
float s = 0;

for(size_t j = 0; j < weights.size(); j++)
s += weights[j];

if(0 == s)
s = numeric_limits<float>::epsilon();

for(size_t j = 0; j < weights.size(); j++)
weights[j] /= s;

// Sum the displacements.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];
displacements[i] += (vertices[neighbour_j] - vertices[i])*weights[j];
}
}

// Apply per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
vertices[i] += displacements[i]*scale;
}

Here's the code to replace the laplace_smooth function with the curvature normal-esque weight scheme (it works most times, but no guarantees:

void indexed_mesh::curvature_normal_smooth(const float scale)
{
vector<vertex_3> displacements(vertices.size(), vertex_3(0, 0, 0));

// Get per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
{
if(0 == vertex_to_vertex_indices[i].size())
continue;

vector<float> weights(vertex_to_vertex_indices[i].size(), 0.0f);

size_t angle_error = 0;

// For each vertex pair (ie. each edge),
// calculate weight based on the two opposing angles (ie. curvature normal scheme).
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t angle_count = 0;

size_t neighbour_j = vertex_to_vertex_indices[i][j];

// Find out which two triangles are shared by the edge.
for(size_t k = 0; k < vertex_to_triangle_indices[i].size(); k++)
{
for(size_t l = 0; l < vertex_to_triangle_indices[neighbour_j].size(); l++)
{
size_t tri0_index = vertex_to_triangle_indices[i][k];
size_t tri1_index = vertex_to_triangle_indices[neighbour_j][l];

// This will occur twice per edge.
if(tri0_index == tri1_index)
{
// Find the third vertex in this triangle (the vertex that doesn't belong to the edge).
for(size_t m = 0; m < 3; m++)
{
// This will occur once per triangle.
if(triangles[tri0_index].vertex_indices[m] != i && triangles[tri0_index].vertex_indices[m] != neighbour_j)
{
size_t opp_vert_index = triangles[tri0_index].vertex_indices[m];

// Get the angle opposite of the edge.
vertex_3 a = vertices[i] - vertices[opp_vert_index];
vertex_3 b = vertices[neighbour_j] - vertices[opp_vert_index];
a.normalize();
b.normalize();

float dot = a.dot(b);

if(-1 > dot)
dot = -1;
else if(1 < dot)
dot = 1;

float angle = acosf(dot);

// Curvature normal weighting.
float slope = tanf(angle);

if(0 == slope)
slope = numeric_limits<float>::epsilon();

// Note: Some weights will be negative, due to obtuse triangles.
// You may wish to do weights[j] += fabsf(1.0f / slope); here.
weights[j] += 1.0f / slope;

angle_count++;

break;
}
}

// Since we found a triangle match, we can skip to the first vertex's next triangle.
break;
}
}
} // End of: Find out which two triangles are shared by the vertex pair.

if(angle_count != 2)
angle_error++;

} // End of: For each vertex pair (ie. each edge).

if(angle_error != 0)
{
cout << "Warning: Vertex " << i << " belongs to " << angle_error << " edges that do not belong to two triangles (" << vertex_to_vertex_indices[i].size() - angle_error << " edges were OK)." << endl;
cout << "Your mesh probably has cracks or holes in it." << endl;
}

// Normalize the weights so that they sum up to 1.
float s = 0;

// Note: Some weights will be negative, due to obtuse triangles.
// You may wish to do s += fabsf(weights[j]); here.
for(size_t j = 0; j < weights.size(); j++)
s += weights[j];

if(0 == s)
s = numeric_limits<float>::epsilon();

for(size_t j = 0; j < weights.size(); j++)
weights[j] /= s;

// Sum the displacements.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];

displacements[i] += (vertices[neighbour_j] - vertices[i])*weights[j];
}
}

// To do: Find out why there are cases where displacement is much, much, much larger than all edge lengths put together.

// Apply per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
vertices[i] += displacements[i]*scale;
}

### #5taby

Posted 24 April 2013 - 09:29 AM

Nice work! This could surely come in handy at some point.

Thank you. I hope it comes it handy.

The paper 'Geometric Signal Processing on Polygonal Meshes' by G. Taubin mentions a few edge weight schemes:

- constant unit weight: w_ij = 1

- inverse edge length weight: w_ij = 1/|e_ij|

- curvature normal-esque (cotan) weight: w_ij = cotan(theta_ij0) + cotan(theta_ij1)

The code in the first post uses the constant unit weight scheme (see indexed_mesh::laplace_smooth()). I didn't put the other schemes into the code zip file because I'm not sure if I'm implementing them correctly, but I'll include them here in case anyone wishes to play with them to see if/where I'm going wrong.

Here's the code to replace the laplace_smooth function with the inverse edge length weight scheme:

void indexed_mesh::inverse_edge_length_smooth(const float scale)
{
vector<vertex_3> displacements(vertices.size(), vertex_3(0, 0, 0));

// Get per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
{
// Skip rogue vertices (which were probably made rogue during a previous
// attempt to fix mesh cracks).
if(0 == vertex_to_vertex_indices[i].size())
continue;

vector<float> weights(vertex_to_vertex_indices[i].size(), 0.0f);

// Calculate weights based on inverse edge lengths.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];

float edge_length = vertices[i].distance(vertices[neighbour_j]);

if(0 == edge_length)
edge_length = numeric_limits<float>::epsilon();

weights[j] = 1.0f / edge_length;
}

// Normalize the weights so that they sum up to 1.
float s = 0;

for(size_t j = 0; j < weights.size(); j++)
s += weights[j];

if(0 == s)
s = numeric_limits<float>::epsilon();

for(size_t j = 0; j < weights.size(); j++)
weights[j] /= s;

// Sum the displacements.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];
displacements[i] += (vertices[neighbour_j] - vertices[i])*weights[j];
}
}

// Apply per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
vertices[i] += displacements[i]*scale;
}

Here's the code to replace the laplace_smooth function with the curvature normal-esque (cotan) weight scheme (it works most times, but no guarantees:

void indexed_mesh::curvature_normal_smooth(const float scale)
{
vector<vertex_3> displacements(vertices.size(), vertex_3(0, 0, 0));

// Get per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
{
if(0 == vertex_to_vertex_indices[i].size())
continue;

vector<float> weights(vertex_to_vertex_indices[i].size(), 0.0f);

size_t angle_error = 0;

// For each vertex pair (ie. each edge),
// calculate weight based on the two opposing angles (ie. curvature normal scheme).
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t angle_count = 0;

size_t neighbour_j = vertex_to_vertex_indices[i][j];

// Find out which two triangles are shared by the edge.
for(size_t k = 0; k < vertex_to_triangle_indices[i].size(); k++)
{
for(size_t l = 0; l < vertex_to_triangle_indices[neighbour_j].size(); l++)
{
size_t tri0_index = vertex_to_triangle_indices[i][k];
size_t tri1_index = vertex_to_triangle_indices[neighbour_j][l];

// This will occur twice per edge.
if(tri0_index == tri1_index)
{
// Find the third vertex in this triangle (the vertex that doesn't belong to the edge).
for(size_t m = 0; m < 3; m++)
{
// This will occur once per triangle.
if(triangles[tri0_index].vertex_indices[m] != i && triangles[tri0_index].vertex_indices[m] != neighbour_j)
{
size_t opp_vert_index = triangles[tri0_index].vertex_indices[m];

// Get the angle opposite of the edge.
vertex_3 a = vertices[i] - vertices[opp_vert_index];
vertex_3 b = vertices[neighbour_j] - vertices[opp_vert_index];
a.normalize();
b.normalize();

float dot = a.dot(b);

if(-1 > dot)
dot = -1;
else if(1 < dot)
dot = 1;

float angle = acosf(dot);

// Curvature normal weighting.
float slope = tanf(angle);

if(0 == slope)
slope = numeric_limits<float>::epsilon();

// Note: Some weights will be negative, due to obtuse triangles.
// You may wish to do weights[j] += fabsf(1.0f / slope); here.
weights[j] += 1.0f / slope;

angle_count++;

break;
}
}

// Since we found a triangle match, we can skip to the first vertex's next triangle.
break;
}
}
} // End of: Find out which two triangles are shared by the vertex pair.

if(angle_count != 2)
angle_error++;

} // End of: For each vertex pair (ie. each edge).

if(angle_error != 0)
{
cout << "Warning: Vertex " << i << " belongs to " << angle_error << " edges that do not belong to two triangles (" << vertex_to_vertex_indices[i].size() - angle_error << " edges were OK)." << endl;
cout << "Your mesh probably has cracks or holes in it." << endl;
}

// Normalize the weights so that they sum up to 1.
float s = 0;

// Note: Some weights will be negative, due to obtuse triangles.
// You may wish to do s += fabsf(weights[j]); here.
for(size_t j = 0; j < weights.size(); j++)
s += weights[j];

if(0 == s)
s = numeric_limits<float>::epsilon();

for(size_t j = 0; j < weights.size(); j++)
weights[j] /= s;

// Sum the displacements.
for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
{
size_t neighbour_j = vertex_to_vertex_indices[i][j];

displacements[i] += (vertices[neighbour_j] - vertices[i])*weights[j];
}
}

// To do: Find out why there are cases where displacement is much, much, much larger than all edge lengths put together.

// Apply per-vertex displacement.
for(size_t i = 0; i < vertices.size(); i++)
vertices[i] += displacements[i]*scale;
}

PARTNERS