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 – out of memory 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

def main():
    # 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.

verts.extend(list(None for i in range(old_edges)))

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.

old_verts = len(verts)

# 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).

# add new vertex only if non-existent
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):

# new external 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.

# new internal edges
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(( 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))

Subdivision code-snippet

The complete code of the subdivision can be seen here:

def subdivide_ply(verts, faces, level):
    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

What do you do when you have 2 idle projectors by your computer? The answer is obviously a high definition projection area to be filled with lo.v.e. (lots of valuable experiments).

Two short throw projectors in one seamless desktop

I’ve been following the work of the Vision3D since 2009. This lab in Montreal is specialized in computer vision (recherche fondamentale et appliquée sur les aspects tridimensionnels de la vision par ordinateur). Lead by Sébastien Roy they have been producing (and sharing!) on calibration of projection surface (e.g. domes o/), multiple projector systems, and content toolsets.

lt-align manual calibration process

The Vision3D lab main tool in that area is Light Twist. This tool was presented in the LGM2009 with a live showcase of the system in a cylinder. In the last week I tried to have light twist going with a multi projector system (aiming to use this for a dome later on) but so far I’m stuck in the playback of content (and I suspect the calibration stage is wrong). Anyways, light twist will be a topic of another post, once I get it up and running.

Plugin enabled – video in the middle of the screens, desktop working normally

Since 2009 the light twist project shifted its focus from labs to end users. In 2011 they finally presented a new project called lt-align and lt-compiz-plugin. The lt-align is a software to quickly calibrate the screens alignment, very easy to use.


The Compiz plugin requires some fooling around with ubuntu settings, but once things are in place it works like a charm. I’m yet to make it work with Unity, so I can have real fullscreen across the desktops.

Recording of the alignment process and video playback

Elephants Dream – Stitched Edition 😉

Note: there is an extra package you need to compile the lt-compiz-plugin:`sudo apt-get install compiz-plugins-main-dev. And I didn’t have to restart compiz with ccp to make it work. Also I changed the shortcuts to start the plugin because Alt+F* were taken by other OS commands.

Time to make it real and project in a large wall

In this picture you can see Djalma Lucio, sys admin that oversees all the computer installations at Visgraf on IMPA. A great professional and a very funny guy to work with. Think about someone that actually enjoys opening a xorg.conf file. And you can also see in the right Aldo Zang. Check it out his ARLuxRender project – a plugin system for LuxRender “which allows to render scenes with mixtures of real and virtual objects directly, without post-processing”.

I hope to post more in the coming months in domes, projections, a special video project … 😉 I went on a 3-month leave of my work at UBC to join the research lab at Visgraf/IMPA, under the coordination of prof. Luiz Velho. This is the second week only, but it’s been already a great experience. And above all, it’s nice to be back home (Rio de Janeiro, Brazil).

Happy Twisting,
Dalai