random vertex colors using Numpy

The bit of research presented in a previous article was not completely academical: speed improvements using Numpy are very welcome when for example generating random vertex colors for very large meshes.

With the results of the measurements in the back of my mind and spurred by a discussion on BlenderArtists I created a new version of the random vertex colors addon. (available on GitHub), which reduced the time to generate random vertex colors for approximately 1 million faces from 3.2 seconds to 0.8 seconds on my machine (or even 0.4s for a 32 bit variant), i.e. up to 8 times faster which is not bad.

Some explanation

Assiging random vertex colors to faces means that you have to assign colors to all the loops of a polygon the same random color (for the difference between vertices and loops, check the BMesh design document. The same applies to the regular Mesh objects we use in the code below). This can be done in a straight forward manner once you have located the vertex color layer:
mesh = context.scene.objects.active.data
vertex_colors = mesh.vertex_colors.active.data
polygons = mesh.polygons

for poly in polygons:
 color = [random(), random(), random()]
 for loop_index in range(poly.loop_start, poly.loop_start + poly.loop_total):
  vertex_colors[loop_index].color = color
Straight forward as this may be, a loop inside a loop is time consuming and so is generating lots of random numbers, especially because Python does not do this in parallel even if you have a processor with multiple cores.

Fortunately for us Blender comes bundled with Numpy, which is a library that can manipulate huge arrays of numbers in a very efficient manner. This allows for a much more efficient approach (although as shown previously a significant speed increase is only noticeable for large meshes).

startloop = np.empty(npolygons, dtype=np.int)
numloops = np.empty(npolygons, dtype=np.int)
polygon_indices = np.empty(npolygons, dtype=np.int)

polygons.foreach_get('index', polygon_indices)
polygons.foreach_get('loop_start', startloop)
polygons.foreach_get('loop_total', numloops)

colors = np.random.random_sample((npolygons,3))
We can even reduce storage (and get an additional speedup) if we change the types to 32 bit variants. There will be no loss of accuracy as these are the sizes used by Blender internally. (Would you do a lot of additional calculations this might be different of course). The change would only alter the declarations of the arrays:
startloop = np.empty(npolygons, dtype=np.int32)
numloops = np.empty(npolygons, dtype=np.int32)
polygon_indices = np.empty(npolygons, dtype=np.int32)

polygons.foreach_get('index', polygon_indices)
polygons.foreach_get('loop_start', startloop)
polygons.foreach_get('loop_total', numloops)

colors = np.random.random_sample((npolygons,3)).astype(np.float32)
As shown above we start out by creating Numpy array that will hold the startloop indices and the number of of loops in each polygon as well as an array that will hold the polygon indices itself. This last one isn't strictly needed for assigning random values because we don't care which random color we assign to which polygon but for other scenarios it might make more sense so we keep it here. We get all these indices from the Mesh object using the fast foreach_get method. We then use Numpy buil-in random_sample function the generate the random colors (3 random floats between 0 and 1) for all polygons.
loopcolors = np.empty((nloops,3)) # or loopcolors = np.empty((nloops,3), dtype=np.float32)
loopcolors[startloop] = colors[polygon_indices]
numloops -= 1
nz = np.flatnonzero(numloops)
while len(nz):
 startloop[nz] += 1
 loopcolors[startloop[nz]] = colors[polygon_indices[nz]]
 numloops[nz] -= 1
 nz = np.flatnonzero(numloops)
The real work is done in the code above: we first create an empty array to hold the colors for all individual loops. Then we assign all the loops with index startloop with the colors that corresponds to the color of the polygon it belongs to. Note that startloop and polygon_indices are arrays with the same length, i.e. the number of polygons.

Now we also have an array numloops which holds for each polygon the number of loops. We decrement this number of loops by one for all polygons in one go (line 3) and create an array of indices of those elements in numloops that are still greater than zero (line 4). If we are still left with one or more indices (line 5) we increment the index of the startloop for all those nonzero indices (line 6).

Line 7 then assigns again in one go a polygon color to a loop at a certain index for all polygons where there are still a non zero number of loops. And finally we again reduce our numloop counters. Note that this all works because the loop indices of all loops of a polygon are consecutive.

loopcolors = loopcolors.flatten()
vertex_colors.foreach_set("color", loopcolors)
The final lines flatten the array of loop colors to a 1-D array and write back the colors to the Mesh object with the fast foreach_set method.

Discussion

Now even though we end up with about four times as much code, it is much faster while still quite readable, as long as you get your head around using index arrays. Basically our single while loops executes the inner loop of assigning colors to loops for all polygons in parallel. The penalty is for this speed increase is memory consumption: We use 5 large numpy arrays here.

Copying vertices to Numpy arrays in Blender

If you want to do a lot of mathematical work on very large Blender meshes using plain Python might be too slow. Fortunately Blender comes bundled with Numpy which allows for fast calculations on vast arrays of numerical data.
To use Numpy you would need to copy for example vertex coordinates first before performing any calculations and then copy the results back again. This consumes extra memory and time. You wouldn't do this for a single operation like determining the average position of all vertices because Python can perform this so efficiently that it wouldn't compensate the extra setup time. However, when calculating many iterations of some erosion algorithm this might be well worth the cost.
In order to help in deciding whether this setup cost is worth the effort I decided to produce some benchmark values. My aim was to find the fastest method from several alternatives and to produce actual timings. Now your timings will be different of course but the trends will be comparable and I have provided the source code for my benchmark programs on GitHub.

results, copying to a Numpy array

First we have a look at how much time it takes to copy the vertex coordinates from a mesh to a Numpy array. Given a the array of vertices me.vertices and a Numpy array verts we have several options:
# blue line
verts = np.array([v.co for v in me.vertices])

# red line
verts = np.fromiter((x for v in me.vertices for x in v.co),
                     dtype=np.float64)
verts.shape = (len(me.vertices), 3)

# orange line
verts = np.fromiter((x for v in me.vertices for x in v.co),
                     dtype=np.float64,
                     count=len(me.vertices)*3)
verts.shape = (len(me.vertices), 3)

# green line
verts = np.empty(count*3, dtype=np.float64)
me.vertices.foreach_get('co', verts)
verts.shape = shape
The blue line represents the creating of a Numpy array from a list. Even with list comprehension this has the disadvantage of creating a intermediate copy, which is also a Python list. Apparently this a costly operation as this naive method performes the worst of all.
The red line shows a huge increase in performance as we use a generator expression to access the coordinates. This will not create an intermediate list but the result is a one dimensional array. Changing the shape of an array in Numpy involves no copying so is very cheap.
We can still improve a little bit on our timing by preallocating the size for the Numpy array. The result is shown with the orange line.
The green line shows the fastest option. Blender provides a fast method to access attributes in a property collection and this method wins hands down. (thanks Linus Yng for pointing this out)

results, copying from a Numpy array

Copying data back from a Numpy array can also be done in more than one way, the most interesting ones are shown below:
# blue line
for i,v in enumerate(me.vertices):
    v.co = verts[i]

# green line
for n,v in zip(verts, me.vertices):
    v.co = n

# light blue line
def f2(n,v):
    v.co = n

for v in itertools.starmap(f2, zip(verts, me.vertices)): pass

# orange line
verts.shape = count*3
me.vertices.foreach_set('co', verts)
verts.shape = shape

The conventional methods perform not significantly different, but again Blenders built-in foreach_set method outperforms all other other options by a mile, just like its foreach_get counterpart.


Discussion

The actual timings on my machine (a intel core i7 running at 4GHz) indicate that we can copy about 1.87 Million vertex coordinates per second from a Blender mesh to a Numpy array and about 1.12 Million vertex coordinates per second from a Numpy array to a Blender mesh using 'conventional methods'. The 50% difference might be due to the fact that the coordinate data in a Blender mesh is scattered over the internal memory (as compared to Numpy's contiguous block of ram) and writing to scattered memory incurs a bigger hit form cache misses and the like than reading from it. Note that a vertex location consists of three 64 bit floating point numbers.

Using foreach_get and foreach_set this performance is greatly improved to 13.8 Million vertex coordinates per second and 10.8 Million vertex coordinates per second respectively (on my machine).

Overall, with this performance this means that on my machine even a simple scaling operation on a million vertices might already be faster with Numpy, like shown in the first example below, which on my machine is about 7x faster (not that scaling is such a perfect example, chances are Blender's scale operator is faster but it gives you an idea how much there is to gain if even for simple operations that are not natively available in Blender.

    # numpy code with fast Blender foreach
    fac = np.array([1.1, 1.1, 1.1])
    count = len(me.vertices)
    shape = (count, 3)
    
    fac = np.array([1.1, 1.1, 1.1])

    verts = np.empty(count*3, dtype=np.float64)
    me.vertices.foreach_get('co', verts)
    verts.shape = shape
    np.multiply(verts,fac,out=verts)
    verts.shape = count*3
    me.vertices.foreach_set('co', verts)
    verts.shape = shape

    # 'classic code'
    fac = Vector((1.1, 1.1, 1.1))

    for v in me.vertices:
        v.co[0] *= fac[0]
        v.co[1] *= fac[1]
        v.co[2] *= fac[2]

Benchmark availability

The .blend file with code and instructions can be downloaded from my GitHub repository.





Blender add-on: create a camera view filling backdrop

whether for product presentations or outdoor scenes, often you need a backdrop that nicely fills the camera view, preferably with a nice curve upwards to hide the horizon and also not larger than necessary to reduce the number of vertices and particles.

Creating this by hand is not massively time consuming but tedious enough to benefit from a small add-on: BackDrop.

Features

By default it creates a flat horizontal mesh at z=0 that fills the camera view exactly. By working with the margin option you can create a backdrop that is slightly bigger than the exact view and unchecking the Zero level option will allow you to position the backdrop on another z-coordinate.

The Lift option will elevate the far end of the backdrop resulting in a curved mesh. The curvature can be adjusted but that doesn't work well yet.

The backdrop mesh is parented to the camera so it will move with it if you rotate or translate the camera.

Availability

The plugin is available on my blenderaddons project on GitHub. Just download backdrop.py and install it from File -> user preferences -> Add-ons -> install from file

Bugs and limitations

The plugin cannot in all circumstances create a suitable backdrop, for example if the camera is pointing upward. If this is the case, a suitable error is displayed in the properties in the toolbar.

Also, the curvature control sucks, and in fact create a curves surface using bezier interpolation is not optiomal since the points are not evenly spaced. For this we need to convert it to an arc length parameterization but I don't have time to do that now.

DumpMesh: clean, complete, sexy

The previous bugfix release of DumpMesh was already able to dump quite a number attributes and attribute layers but in the latest version I went the whole way and dumped every attribute layer that happens to be defined. And I really mean every layer: if you have added a skin modifier the BMSkinVert layer will be dumped for example and even if you assign values to the generic float, string, etc. attribute layers programmatically the will be dumped as well.

Note that vertex groups will not be dumped: vertex groups are not part of a mesh but of an object. (yes, if you didn't know that already, you can have more than one object pointing to the same mesh data, each with its own vertex groups)

This is done by code that is a lot more Pythonic (whatever that means, for me pragmatism is more important than style, as long as it is readable and documented). It uses a lot of introspection to keep the code short and sweet as well. Of course it helps that Blender attribute layers (aka CustomData) is very flexible and well designed (well,... mostly, see Bugs below)

Functionally the add-on hasn't changed much: the interface is cleaner (no more separate selection of what will be dumped, you either choose just geometry or everything), and now all objects that are selected are dumped. This might give a slightly faster workflow.

Availability

The updated code for DumpMesh and CreateMesh can be found on my GitHub repository. DumpMesh can be downloaded directly from this link, as can CreateMesh if your are interested to look at the code.

Bugs

Blenders Python API is generally well designed and adequately documented but there are still some bugs and strange design issues:
BMVertSkin
not documented at all but it is used by the skin modifier. And the bm.verts.layers.skin layer the skin modifier creates has a name that is an empty string
BMTexPoly
has a documentation entry but is marked as todo. Fair enough, because this might be for Ptex stuff that is sort of shelved for now. In the add-on it is mostly ignored.
bm.loops.index_update()
If you add vertices, edges or faces to a bmesh you must ensure that the sequences they are members of can be indexed. You can do this with bm.verts.ensure_lookup_table() or equivalent. If you not only need to be able to do bm.verts[i] but need the actual index as well, like bm.verts[i].index you need to call bm.verts.index_update() otherwise all indices are -1. So far so good, however for loops this appears to be either unfinished or I don't understand it: there is no ensure_lookup_table() for loops, apparently this is taken care of by the faces, but any index attribute of a BMLoop is -1. now bm.faces[i].loops does have an index_update() but it will generate only indices for the loops of that particular face and will start at 0. However when transfering a bmesh to a regular mesh the indices of all loops in the whole mesh are properly initialized, so there is apparently hidden functionality to take care of this. Of course most of the time you can do without but in the addon I dump loop layers as dictionaries of dictionaries that can be indexed as d[faceindex][loopindex] where loopindex is initialy derived from the complete mesh and therefore numbered as one unbroken sequence of all the loops.
factories instead of instantiation from a class
BMesh, BMVert, BMSinVert, etc. are not fully fledged classes. You cannot subclass them and you cannot instantiate them with BMesh(). also most have no __repr__() function and also cannot be pickled. This means that especially for layer attributes like BMSkinVert with more than one attribute you have to know what to assign to which attribute, which is cumbersome. introspection is no help here because you cannot know which attribute is to be saved and which one is just for temporary use.
These aren't deal breakers, just things to keep in mind when you are designing an add-on.

Blender add-on: DumpMesh [Update]

In a previous article I presented a small add-on that could dump all sorts of mesh characteristics like vertex coordinates, faces, edges, seams, etc. as Python code. This code could then be included in other add-ons for example to be used as basic building blocks for parameterized objects.

I have now added the edge selection status and the active uv-map to the list of items that can be dumped. I also fixed a couple of bugs and made dumping most characteristics optional. DumpMesh not only dumps mesh information into a text block but optionally also generates complete add-on code to regenerate the mesh. It serves both as a way to test the generated Python data and might serve as an example of how to create a bmesh object from scratch, including how to add data layers for edges (crease) and loops (uvs).

From a Python point of view DumpMesh.py might be interesting because I had to find a way to include a large chunk of Python code as a string (the source code for the CreateMesh add-on). I solved this by including this source code as a base64 encoded string that is decoded on the fly. This way any quotes or other nasty stuff won't cause the same challenge as trying to create a string literal (which would contain quotes but also escaped quotes etc...).

CreateMesh was even more challenging in a way: DumpMesh generates source code that defines all sorts of list. However, many of these list are optional and all might have a suffix chosen by the user. This means that the CreateMesh add-on must have a way to check if the lists are defined. Fortunately this information is readily accessible in Python because any globally defined variabele is an entry in the dict returned by globals(). Anyway, some annotations can be found at the end of this article.

Availability

The updated code for DumpMesh and CreateMesh can be found on my GitHub repository. DumpMesh can be downloaded directly from this link, as can CreateMesh if your are interested to look at the code.

CreateMesh Python tricks

The lists and dicts with data created by DumpMesh look like this:
suffix = ""

verts = [(-1,-1,-1),(-1.0365,-1,1.0122), ... ]

faces = [(1, 3, 2, 0),(3, 7, 6, 2), ... ]

edges = [(0, 1),(1, 3),(3, 2),(2, 0), ... ]

seams = {0: False,1: False,2: False,3: ... }

crease = {0: 0.0,1: 0.0,2: 0.0,3: 0.0,4: ... }

selected = {0: True,1: True,2: True,3: ... }

uv = {0: {1: (0,0),3: (1,0),2: (1,1),0: (0,1),}, 1: ... }
With a suffix the code would look something like this:
suffix = "_cube"

verts_cube = [(-1,-1,-1),(-1.0365,-1,1.0122), ... ]

faces_cube = [(1, 3, 2, 0),(3, 7, 6, 2), ... ]

edges_cube = [(0, 1),(1, 3),(3, 2),(2, 0), ... ]

...
The geometry() function that creates a bmesh from this data has be clever because for example seams might not be defined at all or might not be called seams but seams_cube for example. This is solved by checking the dict of global symbols returned by globals():
def geometry():

 # we check if certain lists and dicts are defined
 have_seams  = 'seams' + suffix in globals()
 have_crease  = 'crease' + suffix in globals()
 have_selected  = 'selected' + suffix in globals()
 have_uv  = 'uv' + suffix in globals()

 # we deliberately shadow the global entries so we don't have to deal
 # with the suffix if it's there
 verts = globals()['verts' + suffix]
 edges = globals()['edges' + suffix]
 faces = globals()['faces' + suffix]
 if have_seams: seams = globals()['seams' + suffix]
 if have_crease: crease = globals()['crease' + suffix]
 if have_selected: selected = globals()['selected' + suffix]
 if have_uv: uv = globals()['uv' + suffix]

BMesh pitfalls

A peculiarity of the way BMesh works is that adding elements like verts or edges to their respective lists does not guarantee that they are indexable! That means that after adding a number of vertices you must call bm.verts.ensure_lookup_table() in order to access a vertex as for example bm.verts[3]. Now they rela tricky part is that being able to index the list of verts is not sufficient to get the index of a vertex! That is, if you want to retrieve bm.verts[3].index you need to call bm.verts.index_update() first. Even though this is documented here and here, it baffled me at first. (note: what goes for verts, goes for edges and faces as well).

So for verts the complete code looks like this:

 bm = bmesh.new()

 for v in verts:
  bm.verts.new(v)
        # ensure bm.verts can be indexed
 bm.verts.ensure_lookup_table()
        # ensures all bm.verts have an index (= different thing!)
 bm.verts.index_update()

Blender add-on: DumpMesh

DumpMesh is an add-on that will write a lot of information about a currently selected mesh object in the form of Python code. The resulting code is stored as a text object in the Text Editor.

Of course it is far more effective to store an object by saving the .blend file and reuse it afterward but the purpose of this add-on is to get access to mesh data in a form that can easily be used in other add-ons, for example as the basis of some parameterized object, where creating lists of vertex coordinates by hand is very tedious and error prone. It therefore produces Python definitions for

  • a list of vertex coordinates
  • a list of edges
  • a list of faces
  • a dictionary with edge seams
  • a dictionary with edge creases

A small, shortened, example of what is produced is shown below:

verts = [
 (-1.0, -1.0, -1.0),
 (-1.0, -1.0, 1.0),
...
 (1.0, 1.0, 1.0),
]

faces = [
 (1, 3, 2, 0),
 (3, 7, 6, 2),
...
 (5, 7, 3, 1),
]

edges = [
 (0, 1),
 (1, 3),
...
 (0, 4),
]

seams = {
 0: True,
 1: True,
...
 11: False,
}

crease = {
 0: 0.0,
 1: 1.0,
...    
 11: 0.0,
}

By default it also produces code for an operator that recreates a mesh object based on the values in the lists and dictionaries mentioned above.

The latter may sound a bit strange, an add-on that produces code for an add-on, but this way it is very simple to verify that the dumped data indeed produces the desired object when used. Your work flow may look like this:

  • Select the object to dump in the 3D view,
  • Select Object -> DumpMesh, the code will show up in the text editor,
  • Run the script in the text editor by clicking Run script, it will create a new menu entry
  • Select Add -> Mesh -> CreateMesh, and a duplicate mesh should be added to your scene.

If satisfied with the result you can save the generated code for reuse in your own add-on.

Availability

The code for DumpMesh is available on GitHub. Just download the file and use File -> User preferences -> Add-ons -> Install from file in the usual manner. In the code I have embedded the code for the CreateMesh operator as base64 encoded strings. This is in my opinion the simplest way to store a chunk of Python in antother Python program without bothering about quotes of all kinds. The original non-encoded version is on GitHub as well.

Noise experiment with Open Shading Language

I am still looking for good procedural noise to be used as the basis for bark materials. Such noise would need some anisotropy because the cracks in bark are almast always oriented: either along the direction of the trunk or perpendicular to it, depending on the stretching characteristics during growth. Now you could scale for example Perlin noise or use Gabor noise (the latter is available in OSL as well, check for example this article), but we could also create it from scratch.

Theory

In the shader we randomly distribute a number of line segments through space and calculate the distance to these line segments. If we are close enough to a segment we add a value to our sum and in the end the value for a pixel is the sum of all these contributions. Now the essence is in the distribution of those line segments because we do not only control their length but also how far they may deviate from a preferred direction.
The code for the shader itself is quite simple (I have omitted all auxiliary functions):
shader dtnoise(
 point  Pos = P,            // texture position and scale
 float  Scale = 1,
 
 int    Np = 15,            // number of ticks (impulse) per cell
 
 int    Seed = 42,          // random seed
 
 float  Radius = 0.3,       // distance to impulse that will add to result
 float  Size = 0.7,         // length of an impulse
 vector Dir = vector(1,0,0),// direction of anisotropy
 float  Theta = M_PI/2,     // cone (half)angle. The default = 90% which means no anisotropy
 
 float  Step = 0.3,         // cutoff value for Fac output
 float  Gain = 0.3,         // brightness of the result
 
 output float Sum = 0,      // the sum of all the tick values
 output float Fac = 0       // the 
){
 point p = Pos * Scale;
 point f = floor(p);

 vector an= normalize(Dir);
 vector r = vector(an[2],an[1],an[0]); //no need to norm again
 vector uv = cross(an,r);
 vector vv = cross(an,u);

 
 int xx,yy,zz, np;
 
 int nn[3];
 
 for( xx=-1; xx<=1; xx++){
  for( yy=-1; yy<=1; yy++){
   for( zz=-1; zz<=1; zz++){
    point ff = f + vector(xx,yy,zz);
    int s=Seed;
    
    nn[0] = int(ff[0]);
    nn[1] = int(ff[1]);
    nn[2] = int(ff[2]);
    
    for( np=0; np < Np; np++){
     vector pd1 = noise("cell",ff,s);
     vector pd2 = randomcone(uv, vv, an, Theta, nn, s+1);
    
     point p1 = ff + pd1;
     point p2 = p1 + Size * pd2;
    
     float d = distance(p1,p2,p);     
     Sum += (1 - smoothstep(0, Radius, d));
     
     s+=3;
    }
   }
  }
 }
 
 Sum *= Gain;
 Fac = smoothstep(0, Step, Sum);
}
The essence is in line 44 and 45: here we pick a random point and another point along a random direction, where this random direction is restricted to a cone whose narrowness we can control with the Theta parameter. Those two points make up the start and end points of a line segment and if the distance to this line segment is small enough we a a value. (make sure that the sum of Size and Length < 1 to prevent artifacts.

Samples

With few points per cell the anisotropy is very clear (left: Theta = 0.05, right Theta = 1.51 (approx. 90 degrees):
   
 Adding points increases the complexity of the noise (we reduced the gain, the number of points is 40):
 
If we compare our noise to gabor and perlin noise (each thresholded and suitably scaled in a preferred direction) we see that the patterns are similar but each with its own character.

How useful is all is, it up to you of course :-) It is nice to have a different procedural noise in your toolbox but on the other hand, it ties you to OSL (i.e. you can't use it on the GPU) and it is relatively slow.

Code availability

As usual the full code is available on GitHub.

Mmmm, maybe this advertising widget below is a bit over sized. Still, if you would like to learn some more about Open Shading Language for Blender, you might check it out.