saliences.com : Tone Preservation with Ice Crystal Growth

Table of Contents [Hide]

Introduction

This was originally my final research project for a course taught by Professor Craig S. Kaplan
The solidification of condensed water on a cold surface has the property that it creates complex and varying patterns on the surface thanks to the inherent crystallization process and the imperfections that are present on the surface. Typically, this solidification process produces patterns that vaguely resemble the results of a fractal recursion. Whilst there is a level of control over this naturally occurring crystallization process, it is generally insufficient for creating actual art, and it is highly error prone (although examples countering this exist). By using a computer to simulate this natural phenomenon, a high level of directorial control can be given to the end user for creating artistic renderings that use crystals as the “pencil”. It is therefore my hope to be able to create a system wherein the computer can preserve the tone of an input image, all the while “interpreting” the image as a ice crystal formation. This write-up describes the rendering system in full detail and presents its strengths and shortcomings. A successful implementation of this idea is one that successfully captures the tone and edge properties of the input image, and yet still looks “icy” to a degree. Unfortunately, because of the difficulty of modelling these properties in real life, there is no existing metric for determining whether or not the resulting images are “good” or “bad”. The natural ice growth phenomenon does a spectacular job of mimicking the tonal properties of foliage, so attempting to duplicate those results in simulation serves as a good starting point. The system should be general enough that arbitrary images should be rendered faithfully as well.

Previous Works

To my knowledge, the only similar body of work that has been done in relation to ice crystals in a Computer Graphics setting is the paper that inspired this project: Kim and Lin’s 2003 EuroGraphics paper. Their primary interest when developing their system was in a) physically accurate renderings and b) edge preservation [KL1]. This is evidenced by the direction that their research takes with the 2004 and 2006 follow-ups [KL2], [KL3]. Kobiyashi’s crystal growth simulation PDE serves as the workhorse for both this project and the Kim and Lin papers [RK]. But aside from these, there is not much else to report. Certainly crystal growth simulation is a very active field of research, but aside from the preceding papers, no other body of work has extended the idea beyond the crystal’s natural beauty.

Natural Phenomenon

In a natural setting, the formation of the ice frost on a window is readily controlled by placing slightly objects on the frosting surface. Wherever the object was placed will not frost over. This can be performed in layers to get interesting results. Alternatively, placing slightly heated objects on the surface allows for control over where the condensation gathers, giving a different kind of directorial control over the formations. Although the methods are known, there is not many examples of individuals doing art in this fashion, perhaps because of the high error rates associated with the techniques. Figures 1a and 1b provide examples of “art” created using the techniques described above.

Figure 1a
By placing a constraint on the windowframe, ice growth is influenced.
Figure 1b
Dust particles create branching patterns on the windowsill.

The Kobiyashi Formulation

The work horse behind this application is the pair of PDEs described by Kobiyashi which reasonably depict the growth of some forms of ice crystals. Figure 2 shows the PDE pair:

Figure 2 - The Kobiyashi Formulation

The results section contains some solutions to the Kobiyashi Formulation using different initial conditions. The implementation of a solver for the PDE pair is troublesome because the equations are not well suited for higher order explicit differencing stencils such as Rung-Kutta, so the timestep must be very small [KL1], [KL2]. Kim and Lin successfully applied a Lax-Wendroff and Crank-Nicholson algorithm to increase the timestep by a factor of four, but for this implementation I did not want to deal with the added complexity of an implicit method. A basic FTCS technique was used to solve the equation. A short attempt was made at increasing the timestep by replacing the Laplacian terms with a Noyd-Deghgan second order stencil without any success [ND]. This was caused by a lack of understanding of how to adapt the Laplacian terms in the equation to the stencil.

The cross terms in equation (1) can be dropped. This makes the results of the simulations a little smoother, but since the postprocessing step is going to reintroduce details anyways, it is a worthwhile simplification to make [KL2]

The Kim and Lin Pipeline

Kim and Lin describe an elaborate pipeline for turning the solutions of the Kobiyashi formulation into something visually appealing. The primary goal of the pipeline is to reintroduce edges and details into the “blob” that the Kobiyashi formulation creates. In brief, the pipeline consists of the following steps:

  1. Extract Border and Medial Axis images from the Phase field.
  2. Using a Depth First Search algorithm, extract edges segments from the Border and Medial Axis images.
  3. Create a Distance Function map using the Border and Medial Axis images.
  4. Convert the edges into a connected mesh.
  5. Using the Distance Function map, the mesh, and the segments create a subdivision surface version of mesh that will contain the creases and edges desired.

The pipeline is explained in some detail in the source paper, so there is no need to describe it further here. The only difference between this implementation and Kim and Lin’s is that instead of creating a subdivision surface, I generate a Normal Map and render the mesh directly on the GPU, instead of feeding the resulting mesh into another program. Because the paper did not describe the generation of the normal map, I had to improvise the generation of the normal map.

Plate 1 - The top image is an evolved image of a specific initial condition. The middle image is the generated normal map from that image. The bottom plate is the final result of the rendering process.

The Distance Function is generated by iterating over all pixels in the field and determining their distance to the closest Medial Axis edge and the closest Border edge. Then the height values of the Medial Axis edge and the Border edge are interpolated to calculate the height value at the specific pixel. To speed up the generation of the Distance Function map, a KD Tree was used to speed up the search for a nearest edge segment to a point. The gradient of the Distance Function map is used to determine the direction of the normal vector at the specific pixel for the Normal Map. This gives reasonable creasing results as demonstrated by Plate 1.

j and Θ0 mapping

Plate 2 - The top image used constant values for j and θ0. The bottom image used a lookup table for the different values. Notice how much more random the bottom image is.

Equation (2) from defines constants j and Θ0 which determine the number of lobes and the prefered direction of growth in the Kobiyashi Formulation [KL1]. Having these two terms constant across the entire Phase and Temperature fields results in crystals that grow too uniformily in one direction. This in turn reduces the suspension of disbelief of the viewer because that uniformity does not exist in nature. To remedy this problem, a simple solution was taken were the values for j and θ0 are texture mapped. Based on the location on the phase field, different values for j and θ0 are taken. The texture maps used are very low frequency Perlin Noise textures which are then mapped to the applicable ranges for j and θ0. For j the texture is mapped from 2 to 8. θ0 is mapped between 0 and 2π. Plate 2 shows the results of this mapping on the simulation.

Implementation

The entire application was implemented in C++, using several third party libraries which are described in the API Links section of this document. In the interest of saving time, the user is responsible for providing a partially preprocessed PNG image to the program which is then used to set the initial conditions of the joint PDEs. The user is responsible for making sure that the input image has “power of 2” dimensions and that it is in black and white. Additionally, a slight Gaussian Blur tends to give better results. The user also has the option of providing two input images that explicitly define the starting conditions of the PDEs, but minor changes are required in the code for this to work.

The following subsections deal with implementation details that are part of each component of the rendering system.

Edge Preservation

Plate 3 - The top image used constant values for the diffuse colour. The bottom image used a perlin noise texture for the diffuse colour. Note how there is more internal detail in the bottom image.

Edge Preservation is achieved by reimplementing the Kim and Lin pipeline in full. As mentioned previously, the difference between this project and the source implementation is that I use normal mapping on the GPU to achieve the creasing phenomenon that Kim and Lin describe. To add a little more randomness into the final image, high frequency perlin noise was added as a texture to the diffuse colour of the rendered mesh. This adds slightly more realism than the results that Kim and Lin provided as demonstrated by Plate 3. Results from this rendering technique are present in the Results section of this document.

Tone Preservation

Two different attempts were made at creating Tone Preservation with Ice Crystals. The first method, the 3D Kobiyashi Formulation produced results that were far too erratic and did not preserve tone at all. The second method, Stacked 2D Kobiyashi Formulation gives better results, but still significantly falls short in preserving the tone of the input image. Below is a description of the different techniques.

3D Formulation

The first method attempted for tone preservation was to simply extend the Kobiyashi equations into the third dimension and simulate over phase and temperature volumes. The two PDEs are readily adapted by simply adding the z partial derivative into the equation into the expected places. The one part of the equation that is difficult to adapt is equation (2) from Figure 2, which is explictly locked down to two dimensions by Kim and Lin. On further inspection however, equation (2) is an extended Polar Rose equation, which Wejchurt successfully extended into 3D with very interesting results [JW]. Using his formulations of the Spherical Rose equation in place of the Polar Rose equation gave the desired results of extending the Formulation into 3D. Unfortunately, either because of poor coding or poor formulation, the resulting PDE never managed to give adequate results in terms of shapes that resembled ice crystals. No amount of tweaking of the simulation constants resulted in a solution that visually looked like ice similar to the 2D case. As a result, this method was abandoned in favour of the Stacked 2D Formulation described below.

Stacked 2D Formulation

In this method, multiple phase and temperature fields are stacked on top of each other with the hopes that they will approximate the tonal properties of the input image. Each layer is rendered using an adjustable alpha value which causes the different levels to merge with each other into a multitoned image. The key to this method is the choice of initial conditions for each field slice. The initial conditions for each slice were calculated from the input image as follows:

foreach pixel P in temperature field
	I = input_intensity
	if ( I > 1 - slice_number / slice_count )
		P = 1.0
	else
		P = 0.0
	end if
next
foreach pixel P in phase field
	I = input_intensity
	if ( I > 1 - slice_number / slice_count )
		P = 0.0
	else
		P = 1.0
	end if
next
	

This is in effect a screening process whereby each slice is a different thresholded version of the input image. The phase field is then further processed by extracting the border image from the field and using the border image as the phase field for the specific slice level. This pairing of temperature and phase initial conditions ensures that the crystals only ever grow inside the contours defined by the screened image, and not outside. This ensures a strong preservation of the edges of the specific levels. In general 8 slices gives the best results. This is the technique that was used in the Results section of this document. Unfortunately, as demonstrated by the results, this does a very poor job of simulating tone simply because there is insufficient overlap between the different layers to produce the desired mixing results.

GPU-based PDE Simuation

To speed up the simulation of the Kobiyashi PDE, the PDE solver code was migrated onto the GPU using the GLSL shading language. As the PDE simulation is highly paralellizable, the GPU is strongly suited for the task since it is in essense an SIMD device. Kim and Lin achieved a five fold speed increase with the GPUs of the time, and GPU technology has progressed at a blistering pace since 2004. When initially implementing the GLSL version of the PDE solver, the expectation was to get a similar five fold increase bearing in mind that time constraints on the project do not allow for time to heavily optimize the routines. However, the first implementation of the GPU based solver far surpassed expectations by giving up to a 350 fold increase in speed over the CPU version of the solver in some circumstances. Table 1 gives the benchmark results of the simulations for one iteration of the simple FTCS technique used to solve the equations.

Table 1 - PDE Solver Comparisons
Field Dimensions (WxHxD)CPU (avg ms)GPU (avg ms)Relative Increase
256 x 256 x 110701N/A
256 x 256 x 83501350x
256 x 256 x 1612341488x
512 x 512 x 13212314x
512 x 512 x 828393874x
512 x 512 x 1652266580x
1024 x 1024 x 169057283232x

All tests were performed on a Athlon 64 3200+ CPU with 1 GB of RAM. The GPU was a AGP ATI Radeon X1650 PRO with 256 MB of RAM (Catalyst 8.4 Driver). The Operating System is Windows XP SP2. The compiler was Microsoft Visual C++ 2005 with full optimizations.

Results are the averages across multiple iterations

1 - Less than the resolution of the timer used

2 - There was graphical corruption on the results. Possibly driver related. However, the benchmark scales with the previous benchmarks so I still consider it valid.

The one potential concern with these test results is that the OpenGL driver may be executing the commands asyncronously, causing the instrumentation code to return immediately. However, the same instrumentation techniques were used in the past with accurate results [SY1], [SY2]. In addition, the simulation is visibly sped up, so I am inclined to believe the results. In addition it should be noted that because of the “banded optimization” that has been implemented in the CPU version, the CPU version gets progressively slower as the simulation proceeds. The GPU version does not suffer from this fault and produces a uniform speed regardless of the field's contents. As demonstrated by Table 1, it does not even make sense to compare the two different implementations because the GPU version is on average two orders of magnitude faster than the CPU version. When a problem domain is as highly paralellizable as this, the first platform target should be the GPU.

Future Works

Although the results from the tone preservation section of this project leave something to be desired. I believe that the basic technique employed here is sound in that it will be able to produce better results with further tweaking. Longer term, the inclusion of ice melts in to the equation will open up new possibilities for extensions to the basic idea of rendering icy images. If tone preservation is successfully implemented, then the next step of temporally cohesive icy video renderings is the natural extension to this project. Doing temporal cohesion requires ice melts and some speed up to the rendering techniques as they can take up to 10 minutes as it stands, so rendering a video would be prohibitive right now. Adding ice melts should be possible by implementing lifetimes to each pixel in the image. Once the lifetime is up (after n iterations), the reverse of the Kobiyashi process could be applied to the pixel to melt the ice at a physically plausible rate.

The Abstract Framework of the project can give rise to other interesting media as long as the simulation step gives some form of directorial control over the initial conditions. There are numerous PDEs and Fractals that have inherent beauty and complexity in them, and if there was a degree of directorial control over them, some new art forms would emerge out of them. As a contrived example: A rendering of my dog as if it was a Mandelbrot Fractal. The important thing is that the pipeline that is outlined is flexible enough that different simulation parameters can give drastically different results. Speculatively, a mixture of different simulation types stacked on top of each may give the most visually interesting results if done correctly.

Conclusion

This project successfully implements Kim and Lin’s 2003 paper: Visual Simulation of Ice Crystal Growth. The project successfully delineates the edges of an input image and renders them in a style that is similar to frost formations. The source algorithm is successfully refined to remove some of the perceived imperfections in the original paper using j and θ0 mapping, and perlin noise texturing. The algorithm is further extended to account for image tone. However this extension to the pipeline produces unsatisfactory results, failing to successfully preserve tone in all but the most simple of images. Essentially, the hypothesis of stacking ice crystals on top of each other did not pan out for this project and resulted in the substandard output demonstrated in this document. I still believe that it is possible to successfully apply the stacking idea to create renderings that preserve tone, but for the time frame of this project, those goals were not successfully met. Finally, the iterative simulation of the Kobiyashi Formulation was moved onto the GPU using the GLSL shading language for a drastic speed improvement over the CPU version of the same iterative simulator.

Results

Note that the Simulation and Rendering methods were tweaked for each class of images to give optimal results.

Crystal Growth

The following are the results of the Kobiyashi Simulation with specific initial conditions and different simulation constants.

j = 2
j = 4
j = 6
j = 8
Style:

Initial Constraints

The Temperature Map used to constrain the crystal growth.
The resulting crystal.

The following video demonstrates the simulation running with all of the different stages of the Kim and Lin pipeline being shown: ice-tone.avi (2.1 MB)

Edge Preservation

The following are images that preserve the edges of the image that they are rendering.

Style:
Notice how the ice grows along the edges of the grass blades andd crevices of the rock formation. The simulation resolution is a little low, causing the ice growth to be a little too blobby.
Style:
Nicely delineated edges are present in this image, making it suitable for an ice rendering.
Style:
Compare this image to the left image from Figure 1. Notice how both are suspiciously similar in their ice growth patterns.

Tone Preservation

The following are results from the Tone Preservation component of the project. As can be seen, the results are sub-optimal. These two images are representative of the quality of the renderings that the program returns. Other images had similarly bad tone preservation properties.

Style:
Some tone variations can be seen around the lips and neck of the subject. Aside from that the results do not preserve the tone of the image anywhere else.
Style:
The result is essentially disasterous. No tone from the original image is preserved. The best the renderer does is preserve the edges of the girl’s features and some semblance of the flowers in the foreground. But even that is a bit of a stretch.

Bibliography

API Links