Rendering Popsicles

  • by Shiwei Song (shiweis at and Loren Yu (lorenyu at

  • Stanford CS348B Spring 2009


We were inspired by the image above and wanted to create realistic looking popsicles. Popsicles come in different shapes, sizes, flavors, and textures. Their outer appearance is a strong indicator of their taste. Our goal is not just tailoring an implementation to reproduce the image above, but to create an easily parameterizable material that allows creation of realistic looking popsicles of different flavors and shapes.

To tackle this problem, we analyzed many images of popsicles we found on the web to try to come up with reasons why they look the way they are. We then broke down the problem into external surface appearance and internal appearance, each described below.

Link to our original proposal.


We didn't have access to 3D modeling software, so we first tried modeling in Google Sketchup, which worked well. However, it was hard to export the models to pbrt since we didn't have the pro version of Google Sketchup, which costs around $500. We ended up modeling the popsicle procedurally using a python script to generate the popsicle model based on parameters of length, dome radius (radius of top of popsicle), base width, and base corner radius (for rounded rectangle shape of the base).



Surface Appearance (Icyness and Frost)

The surface appearance of a popsicle consists of how icy it looks (how hard or soft it looks) and also the amount of frost on it. Surface appearance nicely fits into the textures and material plugins. Bump mapping is used to give effect of icy surface and frost. An ice texture is created to procedurally generate the bump map. An icepop material is created to take in the bump map and return the correct BRDF depending on whether a point is a frosty or non-frosty region.


In the icepop material, there is an icyness parameter that controls how hard or soft a popsicle surface looks. As experimented in [5], an icy looking surface can be created by bump mapping with FBm noise [1]. A hard/icy surface has ripples like the surface of water with relatively low frequency noise. A soft surface has many small bumps/pores consisting of high frequency noise. As described in [6], FBm noise consists of layering geometric progression of scaled Perlin noise. Since this is a solid texture computable at any point in space, it is perfect here as it works for arbitrary geometry.

We used the general bump mapping function of iscale*FBm(P*ifactor, dpdx, dpdy, iomega, ioctaves) where: iscale is scaling factor controlling how tall the bumps are, P is the point to evaluate, ifactor scales the point and controls the frequency of the noise at the highest level, iomega is the weight of each octave of noise, and ioctaves is the number of octaves used.

The specific values of the parameters above are interpolated from the single icyness parameter, which ranges from 0 to 1 (from not icy to icy). An icy surface will have higher frequency noise and weight of each octave of noise is fairly small. On the other hand, a soft surface has a lot of high frequency noise rather than low frequency noise. From experimentation, we arrived at the settings of iomega = 0.4, ioctaves = 3, iscale = 0.1, ifactor = 1 for an icy surface and iomega = 0.6, ioctaves = 3, iscale = 0.05, ifactor = 10 for a soft surface. Because of time constraints, we did not come up with a complex interpolation function between the two but instead did a linear interpolation as follow:

iomega = 0.6 - 0.2*icyness;
ioctaves = 3;
iscale = 0.05 + 0.05*icyness;
ifactor = 10-9*icyness;

Bump mapping alone does not give an appearance of icy surface. Another aspect of icy surface is specularity, and is implemented in icepop material. When the surface is icy, it is very specular; when it is soft, it is fairly matte because of pores. The icepop material closely relates to plastic as it has both specular and matte components. The specific amount is again controlled by the icyness parameter.

Through experimentation, we found that an icy surface will be highly specular with specular highlight color approaching material color (diffuse color). A soft surface will be matte and a specular highlight color of black. Again the specularity and highlight color is set using a linear interpolation based on icyness. The images below shows a comparison between an icy and soft surface.


Another important surface appearance on a popsicle is frost. We thought about building frost using a physically based approach as described in [4]. However, from looking at images of frost on popsicles, we found that frost on popsicles occur mostly where the wrapper touches popsicle. The frost appearance on popsicles tend to be randomly scattered specks along the contours where it touches the wrapper. So we decided to again utilize solid textures to generate the frost.

As part of our goal to make the system flexible, we want users to be able to input where the frost should be and how it appears. We extended our shape class to take in an array of points defining frost lines (where frost should grow) and floats controlling appearance of each line. For each frost line, the user can control the rate of frost falloff as a point deviates further from the line as well as how dense the frost should be.

For a given point P, first the distance from it to the closest frost line is found. The distance is used to calculate the appropriate frost fall off at that point and determines whether that point contains frost at all. Frost is modeled as a weighted combination of a low frequency Perlin noise (which controls the general overall distribution of the frost) and a high frequency Perlin noise (which controls the size and closeness of the frost specks). The actual function is

dw*((5-mw)*Noise(25/sw*P) + (2+mw)*Noise(1.25/sw*P))

where dw is the fall off weight as a function of distance to frost line and parameter setting, mw adjusts the combination of low and high frequency noise, sw adjusts the size of the noise overall. All these settings are derived from the simple settings passed in.

We thought about building up frost by physically adding geometry to the scene. However, we found that frost does not have a dramatic impact on actual popsicle geometry and bump mapping would be a much faster/easier way to achieve frost. Frost bump map texture is also built into the ice texture. The ice texture can be viewed in a sense as a merge of two different bump map textures. One for the frost, which is calculated by the function above. The other is the surface icyness described in earlier section. If the value of the function above falls below a certain threshold, then that point is considered to be not covered by the frost and the regular icyness bump map value is returned, otherwise the frost value is returned.

Frost also has a different BRDF function than the non-frost regions. To handle this, the ice texture also returns a boolean for whether a point was frost or just regular surface. The icepop material in response will return a different BRDF if the point was frost. Frost is modeled as a fairly matte material that is approximately white. However, it also has a hint of the underlying material color. The final color of the frost consists of 60% underlying color and 40% very light grey. Below is a comparison image of actual frost and an early test view of rendered frost lines (with different settings).

Internal Appearance with Subsurface Scattering

We believed the translucent appearance of the popsicles was due to subsurface scattering, so we set out to implement subsurface scattering using the model described in Jenson et al's "A Practical Model for Subsurface Light Transport" [2]. The image above shows a comparison of plastic popsicles vs. subsurface scattered popsicles. Plastic creates a very solid looking popsicle with harsh edges between triangles of different orientations.

After many hours of tweaking, we realized that the popsicles in our reference image looked the way they look in a large part because they were translucent/transparent, so it was impossible to generate an image that looked like that with the subsurface scattering model, which was better suited for diffuse effects due to multiple scattering. For example, if you look closely at the reference image, you can see a little bit of the popsicle stick inside the popsicles, which implies that the popsicle material is partially transparent. To achieve this effect, we could use volume rendering techniques, which is potential future work.

Due to the impossibility of generating translucent looking ice (which is common in cheap popsicles that are mostly water and food coloring), we decided to render denser looking popsicles like the ones you get from freezing smoothies.

Dipole Approximation

We implemented the dipole approximation using the two-pass rendering approach described in Jensen's "A Rapid Hierarchical Rendering Technique for Translucent Materials" [3]. We implemented this technique in SubsurfaceIntegrator which was a subclass of the SurfaceIntegrator class.

In the first pass, we computed an irradiance estimate at a set of sample points. For each translucent object in the scene (our sampling method is described below) we computed the irradiance from direct lighting by importance sampling all the lights. We used uniform sampling for each light. This first pass was done in the subsurface integrator's Preprocess method, which is called at the beginning of the rendering. Since we only compute irradiance from direct lighting, our implementation does not account for global illumination.

With the irradiance values cached, we can now compute the subsurface scattering approximation. We based the SubsurfaceIntegrator's Li method (which computes the radiance along a ray) off of pbrt's DirectLighting's Li method. Thus, when a ray intersects a non-translucent object, we simply compute the radiance the same way the DirectLighting integrator does.

When a ray intersects a translucent object, we compute the radiant exitance at the point of intersection in the direction of the ray by the following formula:

  • Lo = sumirradiance samples {Dipole diffusion approximation} * {Area associated with sample} * {Irradiance of sample}

where the diffuse reflectance coefficient is

  • {Dipole diffusion approximation} = integralwi BSSRDF(wo , xo , xi , wi) dwi

where xo is the point of intersection between the ray and the object, and xi is the sample position.

First we need the material parameters that dictate how much the material scatters. The absorption coefficient sa is a measure of how much the material absorbs a given wavelength, and the reduced scattering coefficient ss' is a measure of how much the material scatters light. ss' is a measure of the effective amount of scattering, not the actual amount of scattering. The actual amount of scattering ss is related to ss' by ss' = (1 - g)ss , where if g is positive the material is forward scattering and if g is negative the material is backwards scattering. Since g is positive for most materials, it is usually the case that ss' < ss so that's why it is called the reduced scattering coefficient. These two coefficients dictate the scattering properties of the material, and from these two you can compute the reduced extinction coefficient st' = ss' + sa , and the reduced albedo alpha' = ss' / st'. Another useful quantity is the mean free path lu = st' , which gives the average distance a photon travels in the material before it scatters. Finally, the effective transport coefficient str = sqrt(3 sa st').

In our implementation, the diffuse reflectance coefficient is equal to the multiple scattering term of the dipole diffusion approximation in Jenson's paper (we omitted the single scattering term). To compute the dipole diffusion approximation, we first compute where the virtual and real light sources are. The real light source is at a distance of lu below xi. So the position of the real light source is xi - lu n, where ni is the normal at xi. The virtual light source is at xi + lu(1 + 4A/3) n, where A = (1+Fdr)/(1-Fdr) and Fdr is the diffuse Fresnel reflectance. The diffuse Fresnel reflectance is computed from the relative refractive index eta using the following approximation:

  • Fdr = -1.440/eta2 + 0.710/eta + 0.668 + 0.0636eta.

Now, given the distance dr from xo to the real light source and the distance dz from xo to the virtual light source, we now the dipole diffusion approximation can be computed as:

  • {Dipole diffusion approximation} = 0.25 pi alpha' {lu (str + 1/dr) exp( -dr str ) / dr2 + lu(1 + 4A/3)(str + 1/dv) exp( -dv str ) / dv2}.

We approximated the area associated with each sample to be the total surface area of the object's triangle mesh divided by the number of samples. This approximation works well if the samples are reasonably uniformly distributed, and our goal was in fact to retain some of the low frequency variation that comes about from non-uniformly distributed samples which gives the popsicle the appearance of not having a perfectly uniform density.

Note that since ss' and sa are wavelength dependent, all of the quantities we described which are derived from ss' and sa are also wavelength dependent. Thus, our subsurface scattering integrator requires the material properties ss' and sa for each of the R,G, and B frequencies. An extension of pbrt that uses a different Spectrum class could still use our subsurface scattering integrator as long as the appropriate parameters are provided.

Although the refractive index eta is also wavelength dependent, our implementation chose to keep eta constant across wavelengths since the data we were able to obtain from the PBRT textbook provided 1.31 as the refractive index of ice, independent of wavelength. It would be simple to extend our implementation to allow for wavelength dependent indices of refraction. Finally, although lu is also wavelength dependent, we chose to keep it constant across wavelengths since it is easier to use the mean free path as a material parameter representing the amount of "transluceny" of the material than to use the scattering and absorption coefficients directly. We discuss this in more detail below.


The dipole approximation of subsurface scattering needs cached irradiance values over the surface of the object. The paper uses the values by weighting each sample point by the surface area associated with that sample. To reduce noise, there must be enough sample points and the calculated area of the points must be fairly accurate. We approximated the area associated with a sample point to total surface area/total number of samples. This approximation works well for uniformly distributed samples, and in our case it helped to retain some of the low frequency variation that comes about from non-uniformly distributed samples which gives the popsicle the appearance of not having a perfectly uniform density. One way to ensure a more uniform looking surface is to implement Turk's "Re-Tiling Polygonal Surfaces" [7] that guarantees points to be uniformly spaced, which is what Jenson does. A better approach is to find a way a more accurate estimate of the surface area for a point. This way you can have more samples in regions of high curvature or corners where the irradiance varies rapidly across the surface.

First we tried to simply randomly assign points over the surface. This is achieved by giving each triangle in the triangle mesh a probability of getting a point proportional to its area, and then randomly pick uv values to get the actual random point. We found that this method results in clusters of points, which introduces noise.

The second method we tried was to use stratified sampling over each triangle. This is done by first figuring out how many points a triangle should have based on its area and then pick that many samples using triangle stratified sampling. This, however, gives problems when triangles are small and only contain fraction of points. One method to solve this is to round up to a whole number. This causes too many samples in low polygon area. Another method is to give the triangle an additional point with probability of the fractional count. However, this produces too few samples in large number of small triangles. The problems are illustrated below. The left image shows too few points in small triangle and right shows too many points in small triangles.

The final method we tried and used was a hybrid method. We give each triangle a point with probability proportional to its area like the random sampling method. Once we know how many points a triangle have, we stratify sample that many points with some jitter. Comparison of the results are shown below. Note that the points are more evenly distributed than random, which reduces noise. Some noise is still present but we did not feel the need to remove it completely as that would result in a popsicle that looks too perfect to be real.


Choosing parameters for the material to achieve a desired effect is extremely difficult since the relationship between the diffuse color of an object and it's scattering and absorption coefficients is highly nonlinear. The following graph shows the relationship between diffuse reflectance and the material's albedo which is the relative amount the material scatters. The equation represented in the graph is from Jenson's paper and is given by:

  • {diffuse reflection coefficient} = 0.5 alpha' {1 + exp((-4A/3) sqrt(3(1 - alpha'))} exp( -sqrt(3(1 - alpha')) )

In the graph below, we kept the refractive index eta constant at 1.31.

Most materials have an albedo close to 1, and as you can see, a small change of the albedo from 0.95 to 1.00 changes the diffuse reflection coefficient from ~0.4 to 1.0. Thus, in order to reliably reproduce colors that we wanted, we reparameterized the subsurface scattering model in the manner described in Jenson's paper. In particular, given the refractive index of the material and a desired diffuse color, we first computed the reduced albedo for each frequency using a simple binary search. Jenson's paper used the secant method to compute the reduced albedo, but due to the extreme slope of the diffuse reflectance function when the albedo is close to 1.0, this leads to slopes that are close to infinity and numerical roundoff errors. Since the function is monotonically increasing the binary search method is more robust and still extremely fast.

Now, with the compute albedo and the desired mean free path lu, we compute ss' and sa by

  • st' = 1 / lu

  • ss' = alpha' st'

  • sa, = st' - ss'

Problems and Solutions

The sections above mentioned some of the problems we encountered and how we solved them. They will be summarized here along with other problems we came across.

  • Generating models that pbrt can import was a problem. Google SketchUp provided the easiest way for us to model, but it can only export to kmz format. We tried importing this into Maya, export as obj file, and then convert obj to pbrt file. This works but takes too many steps and is not very flexible. We ended up generating the models using a script. The script can export in both obj and pbrt format so we can use available obj visualizers to look at our wireframes.

  • Generating points to cache for subsurface scattering required a fairly uniform sampling. Random sampling created too many cluster of points while stratified sampling per triangle did not work well for small triangles. We solved this by using a hybrid approach.
  • Not all popsicles can be modeled with subsurface scattering. Our reference image, for example, is fairly transparent. We ended up rendering a denser looking popsicle, which is less transparent.

Results and Conclusion

One of our more colorful renderings:

Final rendering:

Work Breakdown

The work was shared pretty evenly. At first we split up the work according to surface appearance and subsurface scattering, where Shiwei worked on surface appearance and Loren worked on subsurface scattering. Later we realized the subsurface scattering had a lot of components to it, which included obtaining the subsurface scattering parameters for the ice material, appropriately sampling the surface of the popsicle geometry for irradiance caching, estimating the resulting radiant exitance using the irradiance samples and the dipole approximation, and combining the diffuse reflection component computed from subsurface scattering with the results from adding the frost. Ideas were discussed together, but Shiwei implemented the sampling method for uniformly sampling a triangle mesh whereas Loren implemented the irradiance computation and computation of the dipole approximation. Loren implemented a simple script to calculate subsurface scattering parameters from the desired diffuse color, but the majority of the work in obtaining and tweaking the parameters for the final image was shared equally.



  1. Elias. "Perlin Noise".

  2. Jenson, et al. "A Practical Model for Subsurface Light Transport".

  3. Jensen and Buhler. "A Rapid Hierarchical Rendering Technique for Translucent Materials".

  4. Kim, et al. "A Hybrid Algorithm for Modeling Ice Formation".

  5. Pan and Mok. "Rendering an Ice Cube". Fall 2006.

  6. Pharr and Humphreys. Physically Based Rendering. Chapter 11 Pg 553.

  7. Turk. "Re-Tiling Polygonal Surfaces".