This post is my presentation for the “Geometry Processing” course I’m attending in IMPA, the National Institute of Pure and Applied Math in Rio de Janeiro, Brazil. http://lhf.impa.br/cursos/pg/
The course is ministred by the professors:
Luiz Henrique de Figueiredo
Luiz Velho
What’s up doc?
The problem proposed was to create an algorithm to subdivide a triangular mesh inscripted in a sphere (e.g. an icosahedron) multiple times in order to approximate the mesh to a complete sphere. We need to solve two different mesh structures:
- STL-like files – where triangles are listed by the coordinates of their vertices.
- OFF-like files – where triangles are listed by the indices of their vertices, and the vertices are listed by their coordinates.
I’m using only .ply files for both approaches. If you choice -m stl and the faces are indexed, they are internally converted to a storage format that treat them as OFF files.
The delivery of the project is the following program/script. It requires python to be installed in the computer.
How to run it
1.1) Download the final script here – subdivider.zip
1.2) Download a sample ply file: icosahedron.ply.
2) Unzip and run the command:
./subdivider.py icosahedron.ply -l 5 -m ply -o icosahedron_5.ply
* in MS Windows you need to call it with the python console or set it to automatically run using the Python installation available.
* if you want to force the program to use the stl method use -m stl in the arguments.
3) The output file is ready. You can check the result in Blender or MeshLab
STL output file | PLY output file
For the rest of this document I’m explaining only the ‘off’ (ply actually) approach. The ‘stl’ one is more a warm up, given that it’s no challenge compared to what comes next. Both algorithms can be found in the complete script above.
Result and performance
Result rendered in Blender, levels 0, 1, 2, 3 and 4
Performance Comparison
| Subdivision | Python | PyPy |
| level 1 | 0.0001 | 0.0002 |
| level 2 | 0.0003 | 0.0005 |
| level 3 | 0.0011 | 0.0019 |
| level 4 | 0.0046 | 0.0074 |
| level 5 | 0.0186 | 0.0881 |
| level 6 | 0.0746 | 0.2395 |
| level 7 | 0.3336 | 0.2790 |
| level 8 | 1.5101 | 0.8054 |
| level 9 | 6.6167 | 3.0026 |
| level 10 | |
11.3281 |
Without any changes in the code, the execution with PyPy showed an improvement of 100%. It still doesn’t compare to native/compiled code, but those are satisfactory results given the simplicity of the implementation and the flexibility of using a script language.
Rational
The core of the code is the subdivision algorithm. I decided to look for a solution where a face could be seen as as much independent as possible. That means no pointer to neighbour triangles and no need to navigate in the triangles in a fancy particular order. The key I pursued is to create unique indices in a ‘universal’ way. A system that could be used by neighbour triangles unaware of each other and that could result in the same values.
This may not be the most efficient algorithm, but it’s very simple to implement (the code in Python ended very compact).
Code overview
# 1) parse args
input, output, method,level = parse_args()
# 2) import ply
raw_verts, raw_faces = import_ply(input)
# 3) convert to internal structure
verts, faces = convert(raw_verts, raw_faces, method)
# 4) subvidide
verts, faces = subdivide(verts, faces, level, method)
# 5) export ply
export_ply(output, verts, faces)
Apart from the subdivision, they are other small steps needed: (1) parse the arguments (level, method, file, …); (2) import the ply file; (3) convert the raw data into an internal structure (more on that soon); (4) subdivide; (5) dump the new data into another ply file.
The conversion (3) rearrange the faces to store not only the vertices indices but also the indices of the edges. It’s a simple code that creates a key for a pair of vertex indices that will always return the same regardless of their order:
key = "%020d%020d" % (min(v1,v2), max(v1,v2)).
If the key was already created we get the id from a dictionary. Otherwise creates a new id and assign it to this key.
Subdivision code explanation
* take a look at the code in the end of this page, or in the download link in the top.
First I allocate new vertices in the verts array. We are going to create one vertex per edge, thus we know from beginning the size of the new array. With the array pre-created it’s easy to see if a specific vertex has been created (by the other triangle that shares this edge) or if it needs to be initialized in the list.
Next I create a unique vertex index that can be guessed/retrieved from both triangles that share the edge where the vertex is generated. Taking advantage of the fact that the number of new vertices is the same as the number of edges in the previous interaction, I use the edge index to directly assign a value for the vertex. old_verts works as an offset to the new vertices indices.
# new vertex indices
v4 = old_verts + ea
v5 = old_verts + eb
v6 = old_verts + ec
Now that we know the exactly id of the vertex, we can see if it was created already. There is no risk of array overflow since the new verts array was already expanded to accommodate all the vertices (filled with None objects).
if verts[v4] == None:
verts[v4] = (verts[v2]+verts[v3]).normalize()
For the new edges, 6 are external and come from the existent edges. It’s crucial to also find a unique index for them. My first attempt was to use an orientation system, always naming the edges clock-wise. This is not reliable though, I found it impossible to implement correctly. The final solution was to name an edge based on its opposite vertex. For example, in the initial triangle the vertex 1 faces the edge 1, the vertex 2 faces the edge 2 and the vertex 3 faces the edge 3.
Also part of the challenge is to split the edge in two new edges with universal indices. I used arbitrary comparison (vertex a > vertex b) to determine which edge will inherit the original index of the edge, and which one will be allocated passed the initial edge offset (old_edges):
e1 = ea + (old_edges if v2 > v3 else 0)
e4 = ea + (old_edges if v2 < v3 else 0)
The internal edges (edge 7, 8 and 9) are brand new and all they need is to get an index from the edge pile and increment the counter. This is the equivalent in C for e7 = edge_inc ++. edge_inc is initialized as an offset to avoid overriding of edge indices.
e7 = edge_inc; edge_inc += 1
Finally it’s time to create the faces following the original schema.
new vertices and edges
new_faces.append(( v6, v2, v4, e1, e8, e6))
new_faces.append(( v5, v4, v3, e4, e2, e9))
new_faces.append(( v6, v4, v5, e9, e7, e8))
Subdivision code-snippet
The complete code of the subdivision can be seen here:
for i in range(level):
old_verts = len(verts)
old_edges = int(len(faces) * 1.5)
edge_inc = old_edges * 2
verts.extend(list(None for i in range(old_edges)))
new_faces = []
for v1, v2, v3, ea, eb, ec in faces:
# new vertex indices
v4 = old_verts + ea
v5 = old_verts + eb
v6 = old_verts + ec
# add new vertex only if non-existent
if verts[v4] == None:
verts[v4] = (verts[v2]+verts[v3]).normalize()
if verts[v5] == None:
verts[v5] = (verts[v3]+verts[v1]).normalize()
if verts[v6] == None:
verts[v6] = (verts[v1]+verts[v2]).normalize()
# note, the following could be 'optmized' to run the
# if only once . the gain is of 3%
# new external edges
e1 = ea + (old_edges if v2 > v3 else 0)
e4 = ea + (old_edges if v2 < v3 else 0)
e2 = eb + (old_edges if v3 > v1 else 0)
e5 = eb + (old_edges if v3 < v1 else 0)
e3 = ec + (old_edges if v1 > v2 else 0)
e6 = ec + (old_edges if v1 < v2 else 0)
# new internal edges
e7 = edge_inc; edge_inc += 1
e8 = edge_inc; edge_inc += 1
e9 = edge_inc; edge_inc += 1
new_faces.append(( v1, v6, v5, e7, e5, e3))
new_faces.append(( v6, v2, v4, e1, e8, e6))
new_faces.append(( v5, v4, v3, e4, e2, e9))
new_faces.append(( v6, v4, v5, e9, e7, e8))
faces=new_faces
return verts, faces
I hope the explanation was clear.
It’s always fun to reinvent the wheel 😉
Dalai