Not only of domes can an artist leave of. The fisheye mode shown in the previous post is sufficient for planetarium productions, but may not satisfy artists looking for ‘real’ fisheye lens. The most common lens found in today’s market are ‘equisolid’ lens (ref: HDRI for CGI).

After a productive holiday we now have something new to play with 😉

model by Jar-Artist

model by Jar-Artist

Models/scene gently provided by Jar-Artist


The proof is in the pudin

It’s “easy” to put up a nice artistic effect together and simply assume everything is working as it should. However, I wanted to make a system that could match the effect produced by real lens. For that I had to build a fail proof experiment. People not into coding may not know, but building a reliable testing and debugging setup is one of the keys for efficient coding. In my opinion it’s always worthy to put time into this. Note: this is also why Blender users can help a lot with bug fixing by simply building proper test files for the (also carefully/methodologically) reported bugs).

1 – Photoshooting

Take some pictures with a tripod rotating the camera on its center (the focal centre actually). We have been doing this for the past two weeks so it was smooth. Those pictures were taken by Aldo Zang in the Visgraf Lab at IMPA.

2 – Stiching

I don’t get tired of recommending Hugin for stitching and panorama making – This open source project sometimes works better even than autopano pro (a pretty good commercial alternative).


3 – Rendering

I put the panorama as a background plate, calibrated the aligment, added a simple floor + spheres. This was done with the (yet to be released) IBL Toolkit. Apart from that my Blender camera needs to match the settings of the real camera+lens I’m aiming at.

In this case all the pictures for the stitching were taken with a Nikon DX2S and a fisheye 10.5mm lens. I created new presets for the sensor dimensions and the render output.

4 – Results

I was quite pleased when I compared the rendered output with the original image. The aspects we should be looking at are only field of view and line distortion across the lens:


Also note the bottom left corner of the photo. This subtle shadowing is due to the vignetting of the lens. This is not present in the cycles render because I’m not implementing a real camera model (as shown here and here).


The complete patch can be found here. The core is the following function (simplified here). I elaborated this from the ‘classic’ fisheye equisolid formula: radius = 2 * focallens * sin ( angle / 2):

__device float3 fisheye_equisolid_to_direction(
float u, float v, float lens, float width, float height)
    u = (u - 0.5) * width;
    v = (v - 0.5) * height;

    float r = sqrt(u*u + v*v);
    float phi = acos(u/r);
    float theta = 2 * asin(r/(2 * lens));

    if (v < 0) phi = -phi;

    return make_float3(

I hope you like it. If you can share any work you did with this feature I would love to see it.

There is a new addon landing. If you work with IBL in Blender for modeling come by soon. In the mean time enjoy the teaser (or poke me to provide some feedback and perhaps even join the alpha testing period).

IBL autosetup for Cycles and Luxrender

IBL Re-Aligment

Home made panoramas also work 😉

Point projection for calibration tweaking

Ideas, questions, comments, feel free to drop a line 😉

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.

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 –
1.2) Download a sample ply file: icosahedron.ply.

2) Unzip and run the command:
./ 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.


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

    return verts, faces

I hope the explanation was clear.
It’s always fun to reinvent the wheel 😉