Andrew Corrigan - Modeling Implicit Surfaces using Compactly-Supported Radial Basis Functions


Description:

In this project I am working on optimizing the computation involved in representing surfaces implicitly using radial basis functions.


Background / Motivation:

Implicit representation of a surface is often intuitive and elegant. Unfortunately, the interactive manipulation of implicit surfaces is not currently feasible. To represent a surface implicitly, one can choose amongst a variety of basis functions. These basis functions are composed to represent this surface. In order to determine how to properly compose these basis functions to represent a surface, a large system of equations must be solved. This operation is not trivial and is a bottleneck in interactive manipulation of implicit surfaces. This project looks to remove that bottleneck.


Sources:

1. Bolz, Farmer, Grinspun, Shroder: Sparse Matrix Solvers on the GPU: Conjugate Gradients and Multigrid Download
2. Morse, Yoo, Chen, Rheingans, Subramanian: Interpolating Implicit Surfaces from Scattered Surface Data Using Compactly Supported Radial Basis Functions Download Abstract
3. David Eberly: Magic Software (sparse system solver) Web Page
4. Isosurface extraction code: Jules Bloomenthal: "An Implicit Surface Polygonizer", jbloom@beauty.gmu.edu in "Graphics Gems IV", Academic Press, 1994
5. Ply model file code: Greg Turk
6. Shewchuk, Jonathan: An Introduction to the Conjugate Gradient Method Without the Agonizing Pain Download



Theory:

1. Compactly Supported Radial Basis Functions

Basis functions are used to create implicit surfaces from scattered samples of the surface. The surfaces are represented as a weighted sum of these basis functions. In [2] the authors provided evidence that compactly supported radial basis functions provide qualitatively comparable results to non-compactly supported basis functions. This is a good thing. The weights of the basis functions are found by solving a large system of equations must be solved. Since radial basis functions are only affected by local constraints the resulting system is very sparse. In general basis functions may be affected globally by constraints, resulting in a dense system. In addition to being sparse, the system is symmetric and positive definite. This results in a significant saving in storing the system, and a faster computation of the system’s solution, by using the conjugate gradients method. Figure 3 shows the effect of a small radius of support, which emphasizes the local influence the constraints have on the radial basis functions.


2. Solving Linear Systems on the GPU

In [1] the authors mapped the Conjugate Gradients method to solve systems of equations onto the GPU. The GPU, is the processor on a computers graphics card. Usually the GPU aids in operations such as transformations, lighting, and scan line conversion. Fragment shaders allow for a programmer to alter the previously fixed rendering pipeline. Fragment shaders are also given access to texture memory. If a system of equations is stored in texture memory, one can implement a custom fragment shader to perform operations such as matrix-vector multiplication very rapidly. The authors of [2] did this, mapping all of the conjugate gradients method to the GPU and reported gains in efficiency.


3. The Conjugate Gradients Method

In [6] the author describes the conjugate gradients method, which can be used to solve sparse, positive definite systems of equations. Given that the system is positive definite, convergence is guaranteed. In addition much less memory is used, since no factorization is required, as in LU decomposition. Intuitively this algorithm works by iteratively minimizing a function, which corresponds to solving the system of equations.


What has been implemented?

1. A hardcoded implicit function has been succesfully sampled and displayed. (Figure 1)
2. Systems have been succesfully setup and solved. (File 1)
3. Interpolation of a hard coded cube with compactly supported radial basis functions. (Figure 2)
4. Same as 3 but with general models.
5. The algorithm has been optimized to use a sparse symmetric representation.
6. Many tests have been performed.


What has yet to be implemented?

1. Performing sparse matrix-vector multiplication on the GPU.
2. Implementation of a reduction operator on the GPU.
3. Mapping the Conjugate Gradients method onto the GPU, using 1 and 2.


Results:

Figure 1. This is x² + y² + z² = 1 evaluated using the Marching Cubes algorithm.



Figure 2. This is a polygonal representation of a model created by interpolating radial basis functions between the corners of a cube.



Figure 3. This is an implicit surface composed of radial basis functions with a relatively small radius of support. This should give an intuition of how radial basis functions are only affected by local constraints.



Figure 4. Two implicit surfaces represented by radial basis functions with different radii of support. The system for the sphere on the left was solved 0.01 seconds faster, and was 0.05 seconds faster to resample.



Figure 4. Two implicit surfaces represented by radial basis functions with different radii of support. The system for the torus on the left was solved 0.441 seconds faster, and was 0.11 seconds faster to resample.



Figure 4. Two implicit surfaces represented by radial basis functions with different radii of support. The system for the bunny on the left was solved 14.041 seconds faster, and was 7.6 seconds faster to resample.




Discussion:

A spreadsheet of my tests on these algorithms (Notice that there is another tab with analysis on these results).

My results are partially a confirmation of those of [2]. Choosing a larger radius of support generally did not signficantly increase the quality of the image visually, at the same time resulting in a large increase in computational cost. The amount of storage was greatly reduced as well for smaller radii of support. In terms of aesthetics there was one issue with using a smaller radius of support. Small protrusions such as ears were troublesome for smaller radii of support. As indicated on the table, the ears on the rabbit, did not receive enough weight, and while the rest of the surface was reasonably smooth, the ears were incomplete. To properly display them, the marching cube size would need to be reduced, but this would increase the computational cost of marching cubes, so one must make a tradeoff either way. For the sphere, once the radius got large, a hole appeared in the surface. This could be due to the constraints on the radial basis functions, cancelling each other out. Using a sparse representation, verse a non sparse representation, produces the exact same result, at a lower cost computationally and also using much less memory. On average building a sparse system required slightly more computation, but compared to the enormous savings in storage, this cost was insignificant.


Downloads:

This program implemented accepts a PLY model as input. It represents the surface using radial basis functions in a non-sparse and sparse representation. The code then uses iso-surface extraction to re-evaluate the resulting surface and outputs it to a PLY file. For details on how to run this program refer to the provided documentation.

Windows version. This includes source, binaries, and visual studio 6 workspace for sparse and non-sparse solver, required WildMagic code, binary for viewer, ply file handling code, and isosurface extraction code.
*nix version. This includes source with makefile for sparse and non-sparse solver, WildMagic code, no viewer, ply file handling code, and isosurface extraction code.

The following will be available later:

Blender PLY Importer. A plug in to import PLY models into Blender.
Blender PLY Exporter. A plug in to export PLY models from Blender.
Sparse Matrix Vector Multiply Sparse matrix vector multiply implemented in Nvidia's Cg
Reduction Operator Reduction operator implemented in Nvidia's Cg
Windows Version. The same project as the one currently available, with conjugate gradient implemented in Nvidia's Cg
*nix Version. The same project as the one currently available, with conjugate gradient implemented in Nvidia's Cg