Stippling applying the Hard Disk Sampling method with a CMYK color basis. Original image from pennlive.com. Click for full resolution.

At the end of writing my post on Stippling and Blue Noise, I left it wondering how would this apply if we wanted to use colors. Just to recapitulate, the idea of stippling is to generate a reproduction of an image by halftoning, making use of fixed-size points distributed on a plane. The points need to be distributed in a way which is both aesthetically pleasing (for which we place the points carefully in a way which shows Blue Noise properties) and be spaced according to the local density of the image. If we use black dots over a white background, this simply means we need to add more dots on the darker areas of the image.

This is all better summarized in this awesome video:

It also illustrates very well why we would want to use a computer program to carry out this task

More than one sample class

So the principle looks simple enough for black and white. Now, if we want to add color, the first immediate question is: Can we not just do the same for each individual R, G, B or C, M, Y, K channel? As it turns out, that would more or less work, but it will look worse than the black and white image. The reason for it is that combining individual Blue Noise distributions does not lead to a result with the same propertie. This is easy to reason about by looking at the first row on Figure 1 below: if we place samples carefully by enforcing a minimum distance between each one of them, but do so on a per-class (color) basis, nothing is preventing points belonging to different classes from ending up close to each other. Similarly, as shown in the second row, merging all samples into a single uniform distribution, and then sepparating them again into individual classes will cause noticeable gaps in the individual classes’ distribution.

Figure 1. Top row: Uniform individual classes does not lead to uniform combined result.
Bottom row: Uniform result does not lead to uniform individual classes.
Image by [Wei 10], click for a better view.

Therefore the Hard-Disk Sampling method boils down to an easy extension of the original dart-throwing approach, where we check this carefully-built conflict matrix to decide whether or not to accept each randomly-generated color sample, producing results that sit inbetween the two extreme cases depicted in Figure 1. This sounds easy-enough to give it a shot with a quick prototype!

Hard-Disk Sampling implementation

Figure 2. Stippling result using a CMYK color basis with subtractive blending. The source image is “Portrait of Jack Nicholson” by Martin Schoeller. Click for full resolution version.

The first detail to bear in mind is that the minimum distance values for each class vary in space: in our stippling case, each individual pixel of the original image has potentially a different color distribution which translates into different importance for each component of the color basis we choose. The distance matrix for a 2D image will thus be 2N dimensional, where N is the number of point classes we consider.

How to choose the point classes depends on the look we’re after: we can use for instance Red, Green and Blue and combine them with additive blending over a black background to emulate a screen. Or we could use Cyan, Magenta, Yellow and Black over a white background to emulate print. Any combination will do, but obviously the more classes and the closer their respective colors are to the original image, the better match. I personally find the results of restrictive color basis like the two aforementioned to be more interesting, as they let you better appreciate how the result emerges from the combination of all the individual samples rather than from a perfect match of each individual sample.

The UI in the prototype provided will let you pick a color basis, which has a blending mode associated: (+) for additive, and (-) for subtractive. Then the density for each individual color class is extracted from each pixel on the source image by projecting the original RGB values to the elements of the chosen basis. I’ve tried different approaches to generate the per-channel densities, and what seems to be working fairly well is to take a projection of each source color into the color basis by doing the dot-product of their resulting vectors in RGB space.

= source pixel

= desired N-class color basis

= densites for , where:

is the projection of the input sample into the i-th color basis vector.
for additive basis, and
for subtractive basis

Note this is a non-linear response to intensity, which I found to work a bit better than the linear one.

Going back to the sampling issue shown in Figure 1, I was wondering: if for our particular case of stippling where we are interested in the combined result and not so much in the individual color channels, could we not get away with the simpler method of a uniform combined result, as shown in the second row, despite of the less uniform distribution of the individual color pigments?

To test that out it was easy to simplify the algorithm to calculate the distance matrix based solely on the brightness of the input pixels, disregarding class information. Then, when once a random sample had been accepted based exclusively on the spatial density, its class is chosen randomly from the available color basis according to a probability distribution matching the underlying pixel color.

Figure 3. Comparison between single class uniform distribution (left) taking only distances between samples into account, and hard disk sampling (right) producing a uniform distribution both per class and on the combined result. Note how the sharp edges around the eyes are defined similarly, but the results differ in the smooth areas. Click for a larger view.

With this test I was expecting to still get the good spatial distribution of a single-class dart throwing, and an on-average correct color: each individual sample would be off, but on average the sum of the samples should more or less match the tone of the underlying region. Figure 3 shows the comparison between the two methods, which differ much more than I was expecting. Truth is that in high-frequency regions the adaptive sampling predominates over the more even distribution of the former method – the difference is more evident in the smoother gradients on the image such as the clothing, where you can see that my naive simplification tends to produce noisier results. Having a uniform distribution for the combined result helps achieving good gradients of density, and having good per-class distribution produces more uniform, smoother tones, with less color clumps.

Retries

An interesting aspect of the Hard Disk Sampling method as proposed in Wei’s paper is that at the core of it, we still rely on random sampling to place points in the image. This is very fast at the beginning but once the image starts to fill up, it becomes increasingly unlikely for a sample to be accepted. More so, since for a given position on the image, we’ve prioritized the classes (the densities of each class), it may be extremely hard to place samples of faint tones in the image by just random chance. Recall from the original stippling post that this is the reason why Dart-throwing is so easy to implement, but also so slow for most practical uses.

Figure 4: A blue sample couldn’t find a suitable place in the image by random chance, so it is ‘forced’ by removing much more likely green samples which will be reintroduced eventually.

There are ways to be smarter about the way we place samples (other than rely on random chance), but the alternative method proposed in the paper is to discard existing samples after we’ve given up trying to place a new one. This is simple and ingenuous: we can think of this as having samples that settle down int he image, and from time to time we need to “shake” portions of the image to dust them off and allow room for improbable samples to settle. More rigurously though what we do is rely on the fact that the “probable” samples we remove to make room have a better chance to settle on the final image again, and the result will converge faster than if we expect the unlikely samples to find a suitable spot by randomness alone. This also leads to a better color rendition with less noise as you can see in Figure 3.

Downloads

I wrote a small prototype application which does color stippling of an input image using the Hard Disk Sampling method described in this post. It will generate both images and XML data files which should be easy to parse and render with other applications.

The code is written in C# and is available as usual on GitHub at: https://github.com/joesfer/ColorStippling
And you can also download the precompiled binaries directly from here.

References

[Wei 10] Li-Yi Wei. 2010. Multi-Class Blue Noise Sampling. Siggraph.

[Cook 86]  Cook, R. L. 1986. Stochastic sampling in computer graphics. ACM Transactions on Graphics 5, 1, 51–72.

Scenes from Wrath of the Titans involve extremely dense volumetric objects. We adapt the Transmittance Function Mapping algorithm for high quality interactive previsualization and tuning of those media. (c) 2012 Warner Bros. Entertainment

It is not very often that I get a chance to talk about what I do in my day job. I am privileged to work with very smart people in great projects, but it is a fact that these projects are normally protected with teeth and claws, and that includes their technology. In fact I would normally keep the topics on this sketchbook deliberately away from possible overlaps to avoid breaking any confidentiality arrangement and also because there’s simply lots of other interesting topics I want to still look into and share!

I was fortunate to collaborate on a project with some colleagues at MPC, which ended up being published at DigiPro 2012, and therefore is public domain now. The text is on a bit more technical tone than I would normally like to have in these posts, but it has some pretty pictures which I think are worth skimming through.

The topic is about rendering of volumetric effects in real-time, including lighting and shadowing. We wanted to give our artists a tool to work with their fluid simulations without always requiring slow offline renders, thus reducing their iteration times.

The technical paper can be found here: Bringing Transmittance Function Maps to the Screen: Wrath of the Titans and Prometheus

Just a quick post to write about the ongoing migration to GitHub and source code revamping I’m currently carrying out.

Since the beginning of this sketchbook it’s been my intention to provide with the source code of the projects with the hope that anyone finds it valuable as a complement for the articles. So far this had been done in a less-than-ideal way by attaching zip/rar files with the contents of a Visual Studio solution. Some other times this hasn’t been even possible since the code depended on my internal libraries and therefore it was hard to isolate. My intention is therefore to slowly migrate the existing code from SVN to git, reorganize the library dependencies as git submodules so this no longer becomes a limitation to distribute code, and upload the suitable projects to GitHub.

Another change will be moving to cross-platform development. So far I’ve been developing primary under Windows/VS, but this has become a limitation the times I’ve been asked about, for example, a Mac versions of a given plugin. Instead of having to compile the code for each version of Windows/Mac/Linux available out there, I’m removing the dependency on the IDE so that the interested people are able to build the code themselves (this will hopefully help finding bugs as well). For this I’m using CMake, which is a cross-platform build system which from a description file will be able to generate Visual Studio solutions under windows, GCC makefiles under Linux, and so on so forth. One exception to this will be the C# projects, for which I’m sticking to Windows for now.

So the GitHub repository can be found at https://github.com/joesfer – I’ll be fleshing it out in the next weeks, please let me know of any issue you may find.

Please bear with me as these are quite big changes and I’m also learning as I go hopefully it’ll be much better in the long run.

It’s been a while since the last post, and truth is I’ve spent most of my spare time lately devoted to my photography. Back however to more technical matters, one of the topics I find fascinating is NPR techniques, and particularly Stippling as an interesting example of them. I’ve been reading about it for some time now, and learnt a couple of interesting things along the way, so I thought I would share them in a post.

When stippling an image, we’re basically interested in generating a distribution of points with a density distribution matching the tones of the original image, the size of the points is generally constant, and we add more or less of them to create darker or lighter tones. In these areas it is important to generate a distribution free of sampling artifacts, which reveal as patterns in the image –the eye is really good picking up these patterns-.

So essentially all we need is a big number of points distributed randomly, and a way of quantizing their required density when we place them over the image. To generate the points, the first approach would be using purely random points. It turns out, however that pure randomness does not lead to aesthetically pleasing results, as points tend to clutter in areas and leave empty spaces. On the other extreme, using a perfectly homogeneous scheme –say a grid- would lead to artifacts in the form of aliasing. We thus need something still random but still somewhat a bit more even.

Random pattern on the left, Blue noise on the right. Image from |Kopf 06|

In the in-between land there’s a flavor of noise called Blue Noise, which is the name given to a particular distribution of frequencies with “minimal low frequency components and no concentrated spikes in energy”. This means that we have no structural aliasing given by the low frequencies, and even spacing with minimal cluttering. It is generally accepted that blue noise properties make for some of the best sampling patterns for a wide range of applications (see [Lopez 10]) .

Noise patterns: the top image shows structural repetition. The middle one displays radial patterns and clustering. The third one displays blue noise characteristics. Images from |Kopf 06|

One way of generating point distributions with blue noise properties is by means of a Poisson-Disk Distribution [Cook 86]. One of the most straightforward ways of generating such distribution is by using a so called Dart-Throwing Algorithm [Cook 86], which generates candidate points in a random fashion over the plane, but will only accept each candidate if it lies at a minimum distance from every previous accepted sample. McCool and Fiume [McCool 82] describe a more practical variant
of dart throwing, where the dart radius is gradually decreased as more samples are placed. The Dart-Throwing algorithm is therefore inherently sequential and slow due to its quadratic complexity, but the quality of the results is nevertheless high.

An important property of any distributions generated with the aforementioned variation proposed by McCool and Fiume is its progressiveness: Similarly to some low-discrepancy sequences, if we put all the generated points into a list of N elements, and draw only a portion of that list [1, N < M], the resulting point distribution also covers the entire plane (as opposed to sweeping it), but simply does it with a lower density. This is important because it means that if we can stipple a spatial region with a given density by plotting all the points of a tile T, we will also be able to generate any lighter density by using the same T by plotting a smaller set of its points.

The amount of points we typically need to stipple an image, which is in the range of thousands in the best of cases, makes the straightforward implementation of the dart-throwing algorithm impracticable. A big part of the problem we’re facing with stippling is how to efficiently generate lots of points in a blue noise distribution. There is a significant amount of research on the matter, and the rest of the post will deal with a couple of approaches I liked and implemented myself, along with notes on related work.

Efficient generation of Poisson-Disk Point Distributions

To generate PD point distributions with blue noise properties there are several families of approaches, among which we have:

• Incremental methods, which place samples sequentially. The Dart-Throwing algorithm described above is a basic example of them, but one could easily come up with more sophisticated implementations which made use of auxiliary structures such as hashing grids or trees to reduce the cost of the nearest-neighbor lookups. A proposal I found particularly interesting is described in the paper [Lopez 10], further details on the implementation of this method are provided below.
• High-throughput methods: which usually rely a heavy precomputation stage but a very fast execution time. A prime example of this kind of algorithms is [Kopf 06], also implemented and described in this post. Also impressive results are obtained by [Wei 08] with a GPU-based Poisson-Disk distribution generation algorithm, orders of magnitude faster than the sequential approach.
• Approaches based on Constrained Voronoi Diagrams (CVD), where each sample corresponds to a cell. These approaches, despite being extremely expensive in time, generate probably the most pleasing distributions. [Hongwei 09] report gains of 10x in execution time in their paper.

Constrained Voronoi Tessellation to generate samples. Image from the paper |Hongwei 09|

The differences between types of approaches are both quality (depicted as even, perceptually pleasing distributions) and execution time. The best results seem to be obtained with CVDs or sequential algorithms, but their execution times are orders of magnitude slower than high-throughput approaches such as the one based in Wang Tiles, which run on real-time. These later generate a more noisy result however, which as we’ll see, may not be so much aesthetically pleasing for the particular case of stippling – Nevertheless blue noise sampling may be used for a variety of applications such as object or texture placement as the authors show in the [Cohen 03], where the density of the distribution will likely be much less dense and speed is crucial.

Adaptive Incremental Sampling

The first one of the approaches I tried is the one described in the paper [Lopez 10], which is a more efficient extension of a sequential algorithm with good quality results and decent execution time.

The key idea behind the approach is simple: having a sample point on the plane, and a minimum acceptance distance limiting how close any other sample can get to it, we build a disk and generate all subsequent samples on its boundary*, instead of randomly on the plane. For each accepted new sample, we calculate the intersection between its disk and the currently active one, and subtract a portion of the active boundary area –so no new samples will overlap-. When there’s no remaining boundary on the active disk, we de-queue it and proceed with the next one.

This defines a pattern in how the samples are generated, as new samples are only generated on the currently active disk, and their corresponding disks queued for later processing. The samples then grow across the image from an initial seed (which it’s quite pretty to watch!). To speed-up disk overlapping checks, a secondary image buffer is carried, where disks are rasterized. In order to accept a sample we have to first check whether the corresponding rasterized disk would overlap a non-empty pixel in the disk buffer.

The acceptance distance for any sample is given by the density distribution, i.e. the amount of black, in the region covered by the sample on the source image

* More precisely, an angle is sampled on the remaining boundary of the active disk. The new disk will be placed on the line that parts from the center along the direction defined by this angle, at a distance given by the image densities of both the active and the new disk.

Stippled image generated with the AIS algorithm. Click for full resolution image.

Implementation notes

As part of the test application provided at the end of this post, I worked on an implementation of the AIS algorithm using C#. Despite the running times I’m getting are not nearly as fast as the ones reported in the paper, they’re not too bad for the quality. Although this is conveniently omitted in the paper, care has to be taken when implementing the disk adjustment heuristic to deal with sub-pixel disks, especially as they tend to make the process diverge, and  may even get the algorithm stuck due to the inability of the disk image buffer to handle disks smaller than 1 pixel. I found using a correct disk area calculation (as opposed to a squared window approximation) to produce less artifacts albeit being slightly slower. Dark images such as the zebra displayed in the screenshot below are remarkably difficult to handle robustly due to the high contrast and sample densities required

This being said, the results look pretty nice in most cases. There’s a free parameter which controls the overall density of the stippled image – not every image is suited for a fixed value of this parameter, and tinkering with it helps achieving optimal results.

Recursive Wang-Tiles

The second approach I tried belongs to the high-throughput family described above. This algorithm relies on a heavy precomputation stage to generate an auxiliary structure called Wang Tiles used in combination of PD distributions to achieve fast running times.

Wang Tiles [Cohen 03] are a construction which enables arranging a finite set of elements (the tiles) in a non-periodic (pseudorandom) sequence to cover the plane. The tiles are squares which have a color assigned to each one of its four edges. The algorithm states that a tile can only be placed if the color of its edges matches those of the surrounding tiles that have been previously placed. For any set of tiles which elements fulfill the required restrictions, an element is chosen randomly.

Set of Wang Tiles arranged in a valid configuration. Image from |Cohen 03|

We start with a tile and sequentially insert random tiles from left to right, simply ensuring that the left edge color on each newly placed tile matches the right edge color of the last one. Once we’ve reached the right bound of the area to cover, we continue building a new row, now matching the left and top edge with the tiles from the previous row. This is righteously called scanline algorithm.

We start with a fixed set of tiles, and place them satisfying 1 or 2 restrictions depending on how many neighbor edges we need to match at any given time. For C possible edge colors, we’ll need at least N*C², N>=1 tiles available to satisfy every restriction on the top and left edges. N determines the number of available choices for any placement, from where we chose using a uniformly distributed random number.

What it’s interesting from the sequence generated by the aforementioned algorithm is its aperiodicity, its lack of structure or repetitive patterns.

Aperiodic sequence generated from a finite set of tiles. Image from |Cohen 03|

Now what do we put in the tiles? Wang Tiles have been used to synthesize textures from samples; we’ll use them to cover the plane with a noise pattern, generated from a finite, relatively small, set of PD distributions. This process is described in the paper [Kopf 06]. Having this, all we need to do at runtime is execute the scanline algorithm to cover the desired image area with tiles from the precomputed set.

Texture synthesis from samples. Image from |Cohen 03|

Noise synthesis from samples. Image from |Cohen 03|

Before that, there are a couple more details to worry about: Edge seams and recursiveness.

Edge Seams

By simply choosing random squared samples of textures or noise patterns and then putting those together next to each, we won’t guarantee that the seams match. We need to build the tiles in such a way that, whenever the edge colors on the opposite edges match, so do the image on the tile.

To do this, the authors propose choosing square samples for the colors themselves, and then combining those sample colors to build each tile (4 source images for each tile) as depicted in the figure below:

Tile merging: a source distribution is generated for each tile, and then merged with 4 additional distribution corresponding to each edge color. This will ensure the contents of adjacent tiles match on the boundaries. Image from |Kopf 06|

This way, since each matching edge comes from the same image, they will align properly.

For the case of blue noise samples however, things get a bit more involved: we not only need the noise patterns to match on the edges, but we also need the resulting merged tiles to preserve their PD properties: we don’t want merged samples to form clusters. The authors of [Kopf 06] propose making use of a Voronoi diagram containing all the samples to be merged, and then build a seam from the path that minimizes a cost function which penalizes nearby samples. Put in plain words, this means merging the source PD distributions by cutting along an irregular seam passing by those areas where points are more far away on average.

Minimum cost seam generated from traversing the Voronoi boundaries. Edge cost is shown in red.

After this stage we’ll have a set of tiles that we can blindly arrange according to the scanline algorithm to generate a sequence which will also preserve PD distribution properties. Note that generating the Wang Tiles is the bulk of the work in this algorithm, but it’s a one off process. After that, we’ll be able to cover as much area as we want with blue noise without having to run a sequential dart throwing algorithm at runtime.

Recursiveness

The second issue we have to worry about is: we know our source tiles are progressive, and therefore plotting a subset of them produces a less dense, but equally filling noise distribution. We use this property to generate light shades on our stippled image – but if we happen to require more density than the maximum our tiles’ distribution can provide, we won’t have enough points!

We thus need to make the tiles recursive: allowing any tile to be subdivided into a set of NxN smaller tiles (built following the same edge-matching rules) which adds up to N² more points than the source tile.

Tile subdivision. Image from |Kopf 06|

Recursiveness along with Progressiveness are the properties that will allow us achieving arbitrarily dark or light shades in the stippled image.

Putting it all together

The stippling algorithm based on Recursive Wang Tiles is based in the following steps:

Preprocessing:

1. Let C be the number of edge colors, generate a complete set of N*C² source PD point distributions.
2. Generate C more PD distributions, one for each edge
3. Build the complete set of unique N*C² Wang Tiles by merging each source distribution with the 4 edge distributions. Note that the original source distribution won’t match on the seams, but these merged ones will.

At runtime: Cover the source image with a non-periodic sequence of the pre-generated tiles; for each tile, decide whether to insert each point based on the required image density. If we’ve inserted all the points on a tile distribution and we still require more density, subdivide the tile.

Stippled image generated using the Recursive Wang Tiles method. Click for full size.

Note that the image generated with this method is somewhat noisier than the result achieved with AIS, however if you run the provided software you’ll notice how much faster this method is (excluding the preprocessing, that is).

Color?

To finish the post, I wanted to provide a hint of what I consider a failed test either as a note to self, or in case anybody reading this is willing to come up with an idea to make it work in the future

AIS extended to color. Pointillism?

Despite personally finding the monochrome stippled images to be beautiful, one trivial extension I thought of was extending the AIS algorithm to plot colored disks by averaging the underlying image pixels. I was aiming for a somewhat pointillist effect, which kind of reminds of. However the image doesn’t completely make it for me because the approach is basing the stipples (the strokes in this color version) exclusively on the darkness of the image, completely ignoring the homogeneity on the tone – this will prevent dark broad strokes from appearing on the result no matter what. Obviously that’s what the stippling algorithm is designed to do: generate dark shades from dense clustering of points, but I do think something like this could be achieved:

Downloads

An application implementing the Adaptive Incremental Sampling and Recursive Wang Tiles algorithms described in this post is provided. The source code is written in C# under Visual Studio 2010. If you just want to run it, there’s binaries already built, along with a pre-generated  set of Wang Tiles to use included in the compressed file.

Precompiled binaries: stippling.zip

Source code repository on GitHub: https://github.com/joesfer/Stippling

Stippling application: AIS and Wang Tiles approaches available

Stippling application: Wang Tiles generator, with (optionally) detailed step breakdown.

References

[Cook 86]  Cook, R. L. 1986. Stochastic sampling in computer graphics. ACM Transactions on Graphics 5, 1, 51–72.

[McCool 82] McCool, M., and Fiume, E. 1992. Hierarchical poisson disk sampling distributions. In Proc. Graphics interface ’92, 94–105.

[Kopf 06] Johannes Kopf and Daniel Cohen-Or and Oliver Deussen and Dani Lischinski. Recursive Wang Tiles for Real-Time Blue Noise. ACM Transactions on Graphics (Proceedings of SIGGRAPH 2006).

[Hongwei 09] Hongwei Li, Diego Nehab, Li-Yi Wei, Pedro Sander, and Chi-Wing Fu. Fast Capacity Constrained Voronoi Tessellation.

[Cohen 03] Michael F Cohen, Jonathan Shade, Stefan Hiller, Oliver Deussen. Wang Tiles for Image and Texture Generation. ACM Transactions on Graphics 22 (3) p. 287.

[Lopez 10] Ignacio Ascencio-Lopez and Oscar Meruvia-Pastor and Hugo Hidalgo-Silva. Adaptive Incremental Stippling using the Poisson-Disk Distribution. Journal of graphics, gpu, and game tools.

[Wei 08] Li-Yi Wei. Parallel Poisson Disk Sampling, in SIGGRAPH 2008, Association for Computing Machinery, Inc., March 2008

Fig 1. Sampling the surface

It’s not unusual that a given algorithm requires values or positions to be sampled from a geometric mesh – some of the projects described in past entries of this blog are examples of it: the Weathering Simulation and Space Colonization Algorithm required points to be placed on the surface of a mesh, wereas the Voronoi Shattering seeds the fragments with positions obtained along the volume of a source object.

In this post I’ll be providing more detail on the implementation of those two kinds of samplers, and at the end of the post you can find the link to some source code which should hopefully help complementing the description.

Sampling the surface of a mesh

Placing random samples along the triangles of a polygonal mesh is the easiest case. The process consists in simply chosing a random triangle, and then generate a point within its bounds until we’ve reached a desired number of samples.

Sampling a triangle can be done using barycentric coordinates as depicted in Figure 1. We start with two variables u and v with random values between 0 and 1 such that u + v ≤ 1 (for this we can keep generating values until the condition is met, or even quicker, just pick the complementary value u = 1 – u, v = 1 – v when u + v > 1). The resulting point is given by the weighed sum of the triangle vertices such as:

s = u*A + v*B + w*C = u*A + v*B + (1 - (u+v)) * C

In order to guarantee a more uniform coverage of the samples, it’s better to weight the probability for a triangle to be sampled according to its area, so that larger triangles will be sampled more often and hold more samples along their surface. We can use the cumulative distribution to pick elements from the triangle array with the following simple steps:

1. Calculate the area of every triangle and store it into an array.

E.g. [ 0.1, 0.4, 0.5, 0.2 ]

2. Find the smaller value on the array and divide every element by it, rounding to the next integer. The array is now holding how proportionally larger each triangle is compared to the smallest one.

ceil( [ 0.1, 0.4, 0.5, 0.2 ] / 0.1) = [ 1, 4, 5, 2 ]

3. Add up the elements of the array and initialize a second array of that length, so that the index i of a triangle is repeated as many times as it’s corresponding values on the first array:

A1: [ 1, 4, 5, 2 ]

A2: [ 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3 ]

4. All that’s left is picking a random element within A2, t = A2[ random(0,1) * length(A2) ] as our triangle index.

Sampling the volume of a mesh

Things get a bit more involved when we want to place samples not just on the surface but also in the volume enclosed by the object. I’ve personally been using a couple of alternatives with different pros and cons that I found useful, an implementation as maya nodes can be found in the provided source code

Using ray intersections

A very straightforward way of sampling a volume is the one illustrated in Figure 2: we begin by obtaining the bounds of the mesh and then iterativelly shoot rays joining two random points on opposite faces.

Fig 2. Sampler by Raymarching

Assuming the mesh is closed, each ray will intersect the geometry in a number of segments (red dashed lines) each defined by an entry and an exit point. We’ll generate samples for each segment, putting care in trying to preserve an overall sampling density #samples on segment = length segment * density, which in the implementation is defined as a “linear density factor” obtained from the bounding box volume and the number of samples requested. Rinse and repeat until we have enough samples.

The pros of this approach is that it fits the shape perfectly since we can precisely obtain the entry and exit points for each ray. Maya provides built-in facilities for retrieving all the intersections on the MFnMesh class by using an internal grid acceleration structure (which is fair enough for moderate meshes). The downsides of this technique are that it may be a bit slow, depending on our requirements, and that the distribution is far from ideal for few samples – although it tends to converge decently when the number of requested samples increases. I found it useful to limit the number of samples a given segment can contribute with to enforce the algorithm shooting more rays and therefore increasing the density of segments instead of the number of samples along each segment.

Sampling voxels

A different approach to sample volumes is to first voxelize the source geometry and then obtain the samples as random points within the voxels. If our voxels are the same size, we just need to pick a random voxel and a random point within that voxel and repeat. If the voxels have different sizes though, we can use the same cumulative distribution concept described above and sample according their volume for a more uniform coverage.

Fig 3. Sampling a voxelized version of the geometry

There are several ways of generating voxels from a mesh (Fig. 3), but it can be done particularly efficiently by using the GPU: the idea is to rasterize the object in a fairly low-res image, and use the pixels of those image as a grid to tell which cells are occupied. This can be done in a multi-pass algorithm by playing with the camera clip planes to “peel” the object into slices rendered to texture. Alternatively, we can do it in a single pass as described in [Eisemann 08] by encoding the depth information in the bits of the color components of each pixel.

The idea is simple but powerful: let’s assume we’re using a render target using 32 bits 8+8+8+8 RGBA channels. For each pixel on the texture, we can set treat the color bits as cells in depth, and then set each one of the 32 different bits to 1 if there is a filled voxel occupying the corresponding cell. We do this by first disabling the depth mask, so no fragments are culled, and rendering the object once from an orthographic projection fitted around the bounding box. For each fragment rendered, we find it’s corresponding bit corresponding to it’s normalized depth in the bounding box, and output a bitmask which is combined in the resulting texture using OpenGL’s bit operator blend mode.

Fig 4. Packing voxel information into color channels

To generate the voxels from the image, we check the value for every pixel, and instead of treating it as color, we treat it as a bit array and output a voxel for each active bit. Even better, by just changing the bitmasks used on each fragment, you can trivially change the behavior of the algorithm from volume to surface sampling

The main benefit of this technique by far is speed: you can voxelize complex objects in real time, and generate the samples from those voxels is cheap as well. On the other hand the number of voxels is limited by the bit depth used in the target texture. With 32 bits per component you can go up to 128 voxels in depth (this is how the provided implementation works), and to go beyond those numbers we’d have to make use of a multipass algorithm as described above. Another drawback is lack of precission: since the voxels won’t match the surface of the object perfectly, the samples we obtain from them may potentially be slithgly outside the original mesh. The problem is of course minimized as we increase the resolution of the voxelization.

Fig 5. Voxel Sampler

Quality of random numbers

Just to make a quick note about this – we’ve been using random uniform variables in these methods to pick elements, or sample points within a triangle or a voxel. A straightforward way of doing this is to use the default rand method in C. For the sake of clarity that’s the method provided in the implementation, but one must be aware that the quality of these pseudo-random number sequences is generally very poor. Even with a better random number generator you may find that chosing points randomly along a surface tends to form clusters and gaps. If that’s the case, better sampling methods can be used instead to generate sequences of values and improve the distribution. I find the illustrations on this paper particularly clear to understand the differences visually.

Downloads

Precompiled plugins for Maya 2011 32 and 64 bit under Windows: sampler.zip
Source code available on GitHub: https://github.com/joesfer/Sampler

The plugin also includes auxiliary preview nodes to help visualizing the samples and the voxels. Since I thought the voxelizer may be useful for purposes other than sampling (or just because producing Lego versions of meshes is fun!) it is encapsulated it in it’s own node for easy reuse.

References

Stack Overflow. Picking a random number based on probabilities.

[Eisemann 08] Elmar Eisemann, Xavier D´ecoret. “Single-Pass GPU Solid Voxelization for Real-Time Applications”. GI ’08: Proceedings of Graphics Interface 2008 322 (2008) 73-80.

Fragment carving

On the previous post we left it at the point where we had a Voronoi Diagram generated from a set of points sampled from the source geometry’s volume. We’ll now take it from there and, in this post, I’ll be talking about the mesh slicing.

Fig. 1. Voronoi cell intersecting source geometry to form fragment

Given the voronoi diagram and the input geometry, the way to obtain the fragments for our shatter algorithm is by computing the intersection of the mesh with each one of the voronoi cells (Fig 1). Since the VD define a convex hull, it means that we can calculate the intersection with the simple following method:

 For each voronoi cell C do For each face F in C do Obtain the plane P containing F Split the geometry in two using P, and discard all the polygons on the front side of P Identify the set of vertices lying exactly on P, and generate  a new set of faces filling the hole. end end 

Fig 2. Slicing the geometry with the faces of a VD cell

The process is depicted in Fig. 2. Note that by filling the holes left from each cut before we move on to the next one -marked as red faces in the picture-, we solve the problem of keeping the shape closed defining a volume.

Hole filling

Fig. 3. First attempt - Naive tessellation

Filling the holes left by each cut while we’re carving the shape is a matter of triangulating a flat set of vertices lying exactly over the cut plane. The vertices will be a combination of new vertices generated by splitting existing faces, and a few more previously existing ones, which happened to be coplanar already. Linking the vertices we’ll have edges describing outlines and holes (e.g. when we slice a torus in two as in Figure 3), and since everything happens on a plane we can just project the coordinates from 3D to 2D, triangulate, and then reproject the result back to 3D.

I initially made the mistake of going for a overly simplistic approach that led to more problems than actual benefits – I’m describing it here anyway possibly as an example of something that didn’t go well!

Version 1 of my hole filling algorithm was trying to identify sequences of linked segments describing loops. Then the relative position of those loops was tested to determine which of them were outlines and which were holes. Next, for every set of outline and its contained holes, every pair of vertices was linked with a new edge unless that edge would intersect with an existing one. The result is depicted in Figure 3.

Unfortunatelly there were several problems with this approach:

Fig. 4. Finding loops

• First of all, identifying the closed loops from a pool of edges turned out to be harder than expected: ideally it would be a matter of finding a sequence of edges v1-v2, v2-v3, v3-v4… eventually finding edge vn-v1 closing the loop. However, tiny imperfections on the mesh slicing such as parallel faces perpendicular to the plane or just tiny or degenerated faces being discarded would cause missing edges -open loops-, inner loops or dangling segments, and an overall actual input being far from the ideal expected, see Figure 4.
• Missing edges were particularly harmful because they would cause loop sequences not to be found, this then resulting in entire holes not being filled. When the shape (now left open) was sliced again, even more edges would be missing, just making things worse and worse.
• Even when things didn’t go wrong, the quality of the triangulation was not very good, being very prone to stretched triangles that would just cause further trouble when being sliced.

The overall method, despite being simple in principle, turned out to be not so simple in practice when having to deal with real-world details, and above all, it happened to be just too fragile for it to work for the general case.

Fig. 5. Higher quality tessellation

Since the quality of the tessellation proved to be crucial for the success of the algorithm, I decided to throw away the first attempt and try a different approach which luckily ended up working much better: in Version 2, given the set of coplanar vertices, we would then perform a Constrained Delaunay Triangulation over it, generating an initial set of triangles filling both the inside of the outline and the holes. A Constrained Delaunay triangulation (CDT) is a generalization of the Delaunay triangulation that forces certain required segments into the triangulation. In our case, we need to ensure the triangulation conforms the shape by including every segment coplanar to the cut plane, with no need for a prior classification to tell whether it belongs to an outline or a hole. Fewer steps involved, and a lot more robust.

Calculating a CDT involves generating a regular DT first, and then progressively insert the required edges by removing all the intersecting triangles and re-tessellating the affected area (Figure 6). I based my implementation on [Domiter 04], which is quite nicely explained. There’s a detail that is not covered in that text, which is how to get rid of the triangles that fill the hole loops (inner circle in Fig. 5). The author proposes using a seed point provided by the user to start a triangle-killing greedy traversal on neighbor triangles. What I did instead in order to provide a (rather ineficient) automatic alternative was to start on each triangle centroid and trace a semi infinite line to count the number of intersections with each one of the source outline/hole edges (this excludes the new edges generated from the triangles) and this way determining whether the triangle is inside or outside the shape -being discarded in the later case-. You can find some explanatory diagrams in [Rupert].

This method vastly improved the quality of the tessellation, and the number of degenerate triangles that came out of further slicing it. Had it not been enough, though, I would probably had required to go one step further and implement a Conforming Delaunay Triangulation which would further reduce the distortion of the resulting triangles by splitting long edges.

Fig. 6. Triangles affected by an edge a-b inserted in a DT

Numerical Robustness

A final note about implementation details: the main issue with the method described above is that it performs many geometrc operations on a relatively small region of the mesh – A given set of faces can be sliced and reconstructed many times before the final shape of a fragment is achieved, and sooner than later this leads to problems in the form of really small, and/or almost degenerate triangles. Think of how wrong things can go when a split plane just slightly touches the geometry, or it cuts it in an odd angle leading very tiny holes which have to be filled just to be cut again on further steps.

Just like any other numeric algorithm, a good deal of care has to be put on the implementation of these operations. A big number of nasty situations can be alleviated by avoiding exact tests and using tolerances or epsilons instead. To give a simple example, when cutting out the mesh with a plane it turned out better to be slightly permissive when testing which side of the plane a vertex was on, just to avoid slicing a triangle over and over again infinitesimally away from an existing vertex. Therefore, instead of doing something like:

If (distance( Vertex, Plane ) < 0) then discard it

You would do something like:

If (distance( Vertex, Plane ) < -epsilon) then discard it

Where epsilon is a small number, so that we avoid incurring in a false positive case where a coplanar vertex decides to return a distance of -0.000001 from the plane. So far so good, but things start smelling fishy again when we find ourselves with a growing number of epsilons, magic numbers, with arbitrary small values. How should we choose them? should they all have the same value?

While it’s true that different magnitudes of “small” may suit each algorithm, and some testing will help finding sensible ranges, instead of using absolute epsilon values, it is better to make the tolerance dependent on the magnitude of the numbers you’re testing.

Fig. 7. Representable numbers using floating point

The reason for this is that real numbers are approximated in a computer internally by using a fixed number of bits providing a finite precision. That precision, however, is not homogeneous along the line of numbers – it converges at it’s best around 0 but then it starts degrading with larger magnitudes (Fig. 7). Taking single-precision floats, the gap between two consecutive representable numbers close to 0 is as little as 0.000001, but only 0.0625 when our numbers are in the range of millions. That’s sixty two thousand times less precise! hence it makes sense to adjust our epsilon to grow it larger with the magnitude of the inputs. If you’re interested in further details, this is all explained in better detail in Christer Ericson’s presentation for GDC ’07, and his blog and related book Real-Time Collision Detection is an invaluable reference for this and many other related topics.

The takeaway of all this is that, when comparing numbers, instead of doing:

Equal( A, B ) { return abs( A - B ) == 0; }

We would then be adding a tolerance value as discussed above:

Equal( A, B, epsilon ) { return abs( A - B ) < epsilon; }

But even better, we’d make the epsilon value depend on the magnitude of our input values to make up for the non-homogeneous precision:

Equal( A, B, epsilon ) { magnitude = max( abs( A ), abs( B ) ); // you may want to clamp magnitude to 1, to avoid further // decreasing epsilon return abs( A - B ) < magnitude * epsilon; }

References

[Domiter 04] Vid Domiter. Constrained Delaunay Triangulation using Plane Subdivision.

[Rupert] Ruppert’s Delaunay Refinement Algorithm.

[Ericson 07] Christer Ericson. Numerical Robustness for Geometric Calculations. GDC 2007.

The project I’ve been working on for some time now is an implementation of a shattering algorithm used to break objects into pieces which could then be plugged to a physics simulation. The technique I went for is so-called Voronoi Shattering because it applies Voronoi Diagrams to generate the fragments, being fairly easy for the artist to define the distribution of those cracks and the overall look.

Shattering a mesh into pieces is an interesting problem that involved a lot more work than I initially thought! – as it is often the case with this kind of geometric algorithms, getting the basic prototype running for easy cases is not too bad, but turning that into a more robust algorithm became 80% of the work. I think it will still be interesting to describe the things that worked along those that didn’t, so I’m splitting this into a couple of posts to keep them manageable.

The problem

Generally speaking, we want to split an object into many solid fragments given by closed surfaces. Provided that we know how to partition the source object, the most immediate problem that arises is that we only have information about the shell of the object, and therefore we must reconstruct the interior volume to figure out how the cracks propagate through that object.

About the breakage process, it’s worth noting that the simulation is itself a very broad subject, and it largely depends on the physical properties of the material we’re simulating: for instance, wereas hard minerals and crystals will generally shatter in clean fractures, other materials such as wood will show some flexibility before breaking, and then splinter off  producing a very different kind of fragments. The problem we’re tackling here lies within the first kind of clean-fracturing materials.

Fig 1. Materials influence the appearance of the fractures

The approach: Voronoi Shattering

Purely procedural simulations have little value if they can’t be controlled. I therefore wanted to work on a solution which, while being automatic, would also give the artist some degree of control about the way the fragments are generated without having to describe each one of them individually.

Before getting into details, let’s have a look at the overall picture of the process to better understand each stage (see Fig. 2): we start with a source polygonal mesh which is sampled to produce a set of seed points directing the slicing. Each seed is turned into a set of planes describing it’s corresponding Voronoi Cell (more on this next) which is used to iteratively slice the mesh. After each split, a hole filling stage adds new faces to the mesh to cover the resulting gap. The process is repeated for all the planes of a seed, for each seed, resulting in a set of convex forms corresponding to the fragments.

Fig 2. Process Overview

I’ll skip the details about the Volume Sampling for now to avoid sidetracking too much, and leave that (interesting) subject for a further post. For now all that matters is that we start with a set of points related somehow to the mesh we’re breaking – they could be obtained by sampling as it is the case here, but also by using other simpler means such as particle systems. Given those seed points telling us where to break our object, the complex task of splitting the mesh is decomposed into an iterative process of bulding each fragment by splitting the source geometry by a different set of planes each time. Let’s move on to see how we obtain the splitting planes from a seed point by using Voronoi Diagrams.

Voronoi Diagrams

Fig 3. Vornoi Diagram in 2D

A Voronoi Diagram (VD) is a way of decompositing the space which takes a set of points and produces an equal number of disjoint convex areas (called cells), one for each point. A cell encloses the space which is closer to its associated point Si than to any other point of S.

Voronoi Diagrams have many applications, and for our purposes we’re particularly interested in the fact that visually they look a lot like cracks – by tweaking the distribution of the initial set of seed points, we can somewhat control the look of how we’re going to break our object. Furthermore, since we’ll be using the cells to slice the mesh, it is very convenient that there will be no intersections among them, and that each cell is guaranteed to be convex; this simplifies the mesh cutting by allowing us to use infinite split planes.

Delaunay Tetrahedralization

There are algorithms that can generate a Voronoi Diagram directly from a set of points, but a more common way of doing so is by using an intermediate Delaunay Triangulation (DT) of the source points. The graph representing the information of a VD is said to be dual of the corresponding graph of a DT, and vice versa. This means that having one, we have all the information we need to reconstruct the other (for example, the vertices of a DT are the centers of each cell in a VD). Therefore, in order to obtain the VD we can first generate a Delaunay Triangulation from the seed points and obtain the dual graph from it. See Figure 4 for a 2D representation.

Fig 4. DT (black) / VD (red) equivalence for the set of yellow points

Things get a bit trickier because we’re working on a 3D domain (our points describe a volume), and therefore instead of a 2D triangulation, we require extending the DT to 3 dimensions; this is called a tetrahedralization, since the primitives are no longer triangles but tetrahedra (pyramids with a triangular base). As you probably guessed, an extra dimension to the problem causes an nice evil explosion of new cases to be taken care of!

There is a lot of literature about DT for the 2-dimensional case,  being a very well studied problem. However the tetrahedralization doesn’t seem to be so popular, and some of the algorithms available for 2D DT don’t extend to higher dimensions. A good solution would be making use of 3rd party libraries such as CGAL, which it’s really complete and free (in fact it’d be arguably the sensible choice to make) – however, for me the purpose of these pet projects to learn new things, so I took the hard way and went for an implementation of my own . A good reference I mostly based upon is [Ledoux], which focuses both on the implementation details as well as handling of the (many) degenerate cases of a DT, using an incremental flip-based algorithm.

The idea for this mesh tetrahedralization is to start with a single large tetrahedron covering the entire domain, and progressively insert all the points, one at a time. On each insertion, we determine which tetrahedron on the existing partition contains the new point, and split that into smaller tetrahedra (the actual number depends on whether the point lies in the inside of the shape or on one of its faces). The algorithm is called flip-based because after each split follows a neighbour-fixing stage where the adjacent tetrahedra are recursively tested for the empty circumsphere criterion and fixed by flipping them in a series of configurations (Figure 5) until they’re Delaunay. The process is guaranteed to converge eventually, assuming we deal robustly with numerical precision issues. Once we’re done with the Delaunay tetrahedralization, we can easily obtain its dual Voronoi Diagram without further computations.

Fig.5. Flips used to construct the tetrahedralization

The cells resulting from the 3D Voronoi Diagram are not actually flat polygons like the ones in Figure 3, but instead convex volumes made up of N-sided faces. An important lesson learnt here is how crucial is getting the DT step right and ensure each tetrahedron produced satisfies the Delaunay Condition (i.e. empty circumspheres) – even slight errors errors in this stage are magnified in the dual VD, resulting in more apparent problems such as cell overlapment or incoherent sizes which will produce horribly-looking cracks. 

To be continued…

This is the end of the first part, I hope it was an interesting read! At this point we got all the elements we need to proceed splitting the mesh. In the next post I’ll elaborate about how the fragments are actually generated from the VD cells, and what didn’t turn up that well in the process…

Downloads

I’ve implemented the algorithm as a Maya plugin, there’s a version compiled for Maya 2011 32-bit (Windows) which you can download from here. In the file you’ll find the .mll along with a sample MEL script to set up a sample scene. Please read the license and note this is a prototype to illustrate the described method, rather than a production-ready solution – no technical support is provided. Use it at your own risk (and enjoy it!).

The source code can be found in github

References

[Ledoux] Hugo Ledoux. Computing the 3D Voronoi Diagram Robustly: An Easy Explanation.

Reaction Diffusion

By playing with the parameters of an RD system, it is possible to simulate a variety of patterns ranging from spots to stripes. Unfortunately I haven’t found huge amounts of information on the subject (let alone up-to-date resources): one of the most complete texts is Greg Turk’s thesis [Turk 91], which covers both Turing-like spots and stripes. For the later, the most extended set of equations is the one proposed by [Meinhardt 82] which generates the kind of images I wanted to achieve.

The main problem in practice, however, is that Meinhardt’s description consists on 5 equations with several magic constants (which might make sense in chemical terms, but the authors seem to obviate) that lead to a large search space. Even more fun, the system is really sensitive to small variations on the value of each constant, and I spent quite a lot of time searching for a set of values which was stable and converged to something meaningful. Luckily, [Asai 99] and [Kondo 09] provide an alternative formulation of Meinhardt’s RD system, and an illustrated set of values to play with, which is how I ended up implementing it.

Below is an applet created with Processing which illustrates the process, it starts with a randomized concentration of the activator an inhibitor agents, and a small set of frozen points in the grid which help creating “vortices”. If you let it run long enough, it should converge to a nice seamless pattern with spirals (not that I’ve seen a fish with that pattern, but it would be cool!).

'Space' key to reset simulation. 'Enter' key to stop. Built with Processing


And it produces some great displacement maps!

References

[Turing 52] Turing, Alan, “The Chemical Basis of Morphogenesis,” Philosophical Transactions of the Royal Society B, Vol. 237, pp. 37–72 (August 14, 1952).

[Meinhardt 82] Meinhardt, Hans, Models of Biological Pattern Formation, Academic Press, London, 1982.

[Turk 91] Turk, Greg, “Texturing Surfaces Using Reaction-Diffusion” Ph.D. Thesis, 1981

[Asai 99] Rihito Asaia, Emiko Taguchia, Yukari Kumea, Mayumi Saitoa, Shigeru Kondoa, “Zebrafish Leopard gene as a component of the putative reaction-diffusion system”.

[Kondo 09] Shigeru Kondo, Hideaki Shirota, “Theoretical analysis of mechanisms that generate the pigmentation pattern of animals”, Division of Biological Science, Graduate School of Science, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8602, Japan

This project arose from a conversation with my friend Gino, who is currently working on Medical Visualization, about the difficulty for an artist to create growing structures such as veins -let alone animate them- and how to make them look organic, and therefore irregular and adapted to the surrounding they grow upon.

The Problem

I decided to tackle the problem, and the first approach that came to my mind were L-Systems, which have proven useful for modelling vegetation. There are two restrictions, however, which discouraged me from taking this approach:

1. We should be able to bind the growing of the structures to a surface, and therefore the growth needs to be context aware.
2. We need to be able to direct the growth by specifying the areas the structure is allowed to occupy.

Doing some research I ended up taking the approach described in the paper Modeling trees with a space colonization algorithm, listed in Algorithmic Botany, which perfectly suits the aforementioned requirements.

Space Colonization Algorithm

Space Colonization Algorithm

The idea for the Space Colonization Algorithm is brilliantly simple: instead of generating the nodes that build our structure by using a set of rules (e.g. a grammar), we begin by producing a set of attractor points by sampling a domain -the surface of an object in our case-. Then we provide one or several start locations, and we iteratively start adding nodes by looking for the set of nearest active attractors, and growing towards them. When a new node is added, it kills all the attractor points which fall within a specified range, preventing other branches to fill in the same space, and directing the growth forward.

Implementation in Maya

Given the stages of the SCA algorithm, the natural solution felt to devise them as a set of separate DAG nodes which will be linked together. This allows Maya to cache the result of each individual stage and update the graph efficiently when a given parameter is tweaked without requiring the whole solution to be re-evaluated.

• The Sampler node: takes a mesh as the input and produces a number of attractor points by sampling it. The more points we generate, the better the growing structure will fit upon the surface, but that will also tend to generate denser heavier structures.
• The Grower node: implements the SCA algorithm. It is fed with the set of attractor points, which are stored in a kd-tree for efficient nearest neighbor queries. An input locator provides the start position for the growth, and the exposed parameters of the algorithm control how the attractors are consumed to produce the resulting structure:
  • Search Radius: specifies what’s the maximum distance an attractor point can influence over a node. It is given in % of the mesh bounding box extents. A high value will link distant areas, allowing the structure to “jump” between them.
  • Kill Radius: The range within an attractor point will be killed if it is near a grower node. Large kill radius result in smoother branches, whereas lower values will generate curlier shapes.
  • Grow Distance: the distance between nodes, the lower the value, the denser the resulting geometry.
  • Max Neighbors: for the nearest neighbor searches – controls how many attractor points affect a given node at a time. Large values will smooth out the branches, and a low number of neighbors on the other hand will produce noisier more organic lines (up to a limit).
• The Trimmer node: provides a Percent Length parameter which can be used to animate the growth of the structure generated.
• The GrowerShape node: takes the (trimmed) result of the SCA and generates a polygonal mesh. It exposes a Tube Sections, and Thickness parameters, which control the smoothness and size of the branches.

The node hierarchy can be easily built with the provided command “grow” -  Just select a mesh and an optional a locator node, and type grow on the MEL console.

Directing the growth

The sampler node allows restricting which faces are used to generate the attractor points. To do so, it uses a vertex color set associated with the input mesh: the lightness of the vertex colors is used to determine the probability for each face to be sampled, hence, black faces will be excluded from sampling. See the example video below.

Download

If you feel like trying it out, I have compiled a version of the plugin for Maya 2011 32 and 64-bit. You can find it here.

Also, the source code is hosted under github

Isosurfaces

An isosurface is a shape which is described by taking a scalar field and building a surface from all the points on the field of which value matches a given threshold level we provide. To generate a polygonal surface we can use in computer graphics, our scalar field will be given by a 3D matrix -or grid- of numbers, which we will iterate in order to fetch the points that describe our shape. Unlike the mathematical scalar field, which  is defined in every point, our grid describes a discreet sampling of the field, for we only know the values of the field in the vertices of the grid cells.

Fig 1. Some of the Marching Cubes configurations

To build the surface from this discretized field, we generate triangles from each one of the cells conforming the grid by using the Marching Cubes algorithm. Given a grid cell, we mark each one of the 8 vertices as above or below the threshold level we provided. From those values, the algorithm provides a set of triangles which describe the surface inside that cell (Fig 1).

Maya Implementation

This particular implementation I provide here is integrated with Maya as a Surface Plug-in. It uses Maya’s particle systems to generate the scalar field – for each particle, we define a radius that describes the sphere of influence of that particle. The influence rapidly decays to zero as soon as we move away from the particle. Next, a grid within the particle’s bounding box is created, and sampled at a frequency given by the Grid Res attribute in the shape node. The surface node can also be plugged the particle colors and interpolate them in the resulting triangles, producing nice color blends. The higher the sampling frequency, the smaller the triangles and the more detailed surface we obtain. Bear in mind however that the number of cells grows proportionally to the cube of the grid resolution, so the surface generation quickly becomes a processor-intensive task.

Optimization

Fig 2. Fraction of the cells actually used

I started prototyping the plug-in  using Maya Python, but it soon became too slow to be practical, and I had to switch to C++. While working on it, I realized that for “splashy” particle systems such as the one depicted in the picture above, applying the triangulation to the entire grid is very wasteful, as only a small portion of the total cells (around 20% in that particular example) tend to be touched by the particles and thus will actually produce triangles. Therefore I optimized the triangle generation by performing a first pass over the grid and tagging those cells that will potentially be affected by the particles, skipping all the rest. On a second pass, the marching cubes is only applied to those tagged cells (colored in white in Fig 2.). Finally, Maya seems to like copying data around when handling nodes, and the internal cell occupancy cache can generate a lot of memory fragmentation if created/deleted every frame. To avoid this, a lazy-copy wrapper around te memory chunk is provided, so that the number of memory operations is minimized.

Downloads

Finally, if you would like to try it, I have posted a version of the plug-in compiled for Maya 2011 x64 here, as well as a sample scene. For the programmers among us, the source code can be found here.