Scholars Research, Summer 2006: Brandon Butch

Faculty Advisor: Professor Dinh, Computer Science Department

Subdivision applied to an initial coarse vector field

Introduction
The purpose of this research project was to integrate several related fields of computer graphics in order to acheive the ultimate goal of simulating the realistic erosion of surfaces. Many factors govern the erosion of an arbitrary surface, but arguably the most imporant is the means by which the surface is being eroded - Is the surface being eroded by wind, water or rock particles? In which direction are these particles traveling? At what point on the surface are they hitting? While factors such as the hardness and shape of the target surface will, for obvious reasons, be an important part in any such simulation, one of the most prominent aspects will be the external conditions of the system under investigation. For my research, I investigated and implemented methods that allow a user to specify these important conditions and how they might change over time without placing upon the user the burden of providing explicit and extremely detailed information. This is made possible by a refinement technique called subdivision, which interpolates unspecified informtion based on the data that is provided. Since many of the external forces involved in erosion (such as wind, water particles) can be described by the laws of fluid dynamics, using a subdivision scheme based on the laws of fluids would yield an interpolation method that holds true to the user's original intent. Thus, the user need only give some initial boundary conditions that must be adhered to and our subdivision scheme will take care of all the behind the scenes work. Futhermore, it is possible to alter the subdivision schemes to include dynamic elements, not only allowing the user the specify some initial bounday conditions, but how the resulting forces might change over the course of the simulation. It is natural to think of these external conditions as a vector field in which the target object is placed. This is perfect, as our subdivision allows for the refinement of coarse vector fields.

Subdivision Schemes
Subdivision schemes are generally algorithms used to iteratively calculate refined approximations from some initial, coarse shape. Each successive iteration calculates the next refinement by using the previous approximation (or the initial specification, in the case of the first calculated approximation), thereby refining the initial shape to an arbitrary degree. In other words, given some coarse shape, a subdivision scheme defines a relation of the form U_k = (S)U_k-1, where U₀ is the initial shape to which the subdivision scheme is applied. In such a relation, the S term is a predictor (dependent on the shape with which the scheme is being run as well as the initial conditions) which, when applied to the previous approximation, U_k-1, returns a more refined shape. This S term is derived by using appropriate boundary conditions for the problem at hand (see below). Theoretically, taking the limit of iterations to infinity would yield the continuous shape sought after. For this project, I was concerned mainly with subdivision schemes for discrete vector fields. Henrik Weimer and Joe Warren, in their paper "Subdivision Schemes for Fluid Flow," outline such a scheme, using the laws of fluid dynamics (specifically the Navier-Stokes equations) to dervie the predictor for the subdivision relation.

Single Variable Subdivision - Polynomial Splines
I began by investigating subdivision schemes for curve representations of polynomial splines of a single variable, so as to get a better understanding of subdivision schemes before proceeding to the more complicated multi-variable, vector field case. In the polynomial spline case, a coarse curve is specified by a generating function in which the coefficients of the function correspond with points over the initial, coarse domain. After each iteration of subdivision, the curve is smoothed by having additional points interpolated in the spaces of the coarse domain, i.e. the domain is mapped from one consisting of N points to one consisting of 2N-1 points while remaining within the same bounds. Let u(x) represent the initial generating function of order m to which we wish to apply the subdivision scheme. Then the relation to be used is defined as: u_k(x) = s(x)u_k-1(x²) where s(x) is the subdivision mask given by: s(x) = [1/2^m-1][(1 + x)^m/x^m/2]. In this case, the subdivision mask, s(x), was derived based on the condition that for a polynomial spline, u(x), of order m, the following differential equation must be satisfied: D(x)^mu(x) = 0, where D(x)^m is the differential operator applied m times.

For example, consider applying the subdivision scheme to an initial generating function of: 5 + 2x + 10x² + 7x³ + 2x⁴ over the initial coarse domain D₀ = {0, 1, 2, 3, 4}. The first round of subdivision will map this domain to the new, more refined domain D₁ = {0, 1/2, 1, 3/2, 2, 5/2, 3, 7/2, 4}, while interpolating the new coefficients for this domain. Each successive round will continue this process. Shown below are several rounds of subdivision applied to this function. In each frame. the initial points are displayed in blue, with the current refinement in green.
Click any frame for a larger picture


I implemented this subdivision scheme in Matlab such that it will refine an initial set of points with the above specifications. Each of these plots was created by that program.

Multivariable Subdivision - Two Dimensional Vector Fields
The two dimensional vector field case is very similiar to the polynomial spline case in that much of the underlying framework can be easily extended. The coarse initial conditions for a two dimensional vector field are specified by a pair of generating functions, u(x,y) and v(x,y). Since we are concerned with fluid flow the coefficients of u(x,y) and v(x,y) specify the u and v components of the velocity of the fluid at any given point (x,y) that lies in the domain. The same subdivision relation holds - that is, each new approximation is calculated by applying a subdivision mask to the previous result, and on each round of subdivision the domain is taken from an NxN two dimensional domain to a 2N-1x2N-1 refined domain. The initial discussion of subdivision masks still holds, only now different boundary conditions must be considered. Weimer and Warren derived the subdivision mask by considering the following differential equation which governs slow flow: u^(2,1)(x,y) + u^(0,3)(x,y) - v^(3,0)(x,y) - v^(1,2)(x,y) = 0. The (i,j) superscript represents the i^th partial derivative with respect to x and the j^th partial derivative with respect to y. By manipulating this boundary equation they produced a matrix relation for the subdivision scheme associated with two dimensional slow flow: where the matrix of generating functions representing the subdivision mask can be calculated by numerically solving the linear program representing the following: In this relation, the notation D[x²] represents the discrete differential operator with x² substituted into the operator and L[x,y] = D[x] + D[y], the discrete Laplacian operator. Weimer and Warren made publicly available the pre-calculated coefficients for each of the generating functions in the subdivision mask and it was with these values that I implemented the subdivision scheme. The following sequence of frames displays the subdivision scheme, again implemented in Matlab, as applied to an initial coarse vector field. The subdivision scheme refines the vector field, interpolating the values for which nothing was specified, over an increasingly dense domain. In this particular example, the subdivision was applied initially to a 3x3 vector field where all the values were specified to be zero, with the exception of the two vectors lying on the center line of equal and opposite value.

Click any frame for a larger picture


It is not hard to see from this example that the subdivision scheme refines the vector field in a way that is intuitively pleasing. Given the initial vector field from the first frame, the subsequent frames define refined fields that, in the context of fluid flow, make sense from a physically intuitive point of view. Below is another example, this one being perhaps defined by an initial field whose conditions are more like those one would expect from fluid flow.

Click any frame for a larger picture

Subdivision Over Irregular Domains

It is immediately clear from the above two examples why a subdivision scheme such as this can be a great benefit to a user wishing to specify some initial conditions for simulation. With minimal effort, the user can specify a coarse vector field that represents the forces of the system in which an object is being eroded. The subdivision scheme will then interpolate the remaining data and yield a dense field that accurately represents external forces governed by the laws of fluid dynamics. Obviously, this greatly reduces the burden of the user. While this is a great advantage, it would be even more helpful if the user were able to specify additional information or place constraints on this field. For example, the user may wish to specify a dynamic field - one that varies with time, or one with boundaries that confine the field to an irregular domain. Such conditions deviate from the standard subdivision scheme for vector fields and thus require the subdivision scheme to be altered. Therefore, I explored ways to extend the standard subdivision scheme already implemented to include such conditions.

In their derivation of the subdivision scheme for fluid flow, Weimer and Warren dealt only with refining coarse vector fields that were initially specified over domains of equal dimension - an NxN grid in the two dimensional case. It seems that it would be beneficial to extend this framework to a more general setting, one in which the domain were not restricted to an NxN grid. To achieve this, I altered the original subdivision scheme put forth by Weimer and Warren. Their scheme implicitly assumed that every point in the domain should be given a refined value after each iteration. By restricting the scheme to eliminate this assumption it is possible to refine coarse vector fields over vastly varied domains. Particularly, this can be incorporated into the subdivision scheme by resetting the values of certain sections of the domain after each iteration. For example, if one were studying the flow of fluid through an area which contained an obstruction, the scheme would consider the entire area (including the obstruction), but after each iteration the points in the area containing the obstruction would be reset to zero, as no fluid would flow through the obstruction. (Remember, the values of each vector represent the fluid velocity, so this modification ensures that the fluid velocity is always zero inside the obstruction.) The following screen shots display this type of subdivision.

Click any frame for a larger picture


Dynamic Elements
I attempted to incorporate an element of dynamics into the subdivision scheme by allowing an object to be specified in the way described above (nothing more than a set of boundary conditions that allow no flow in a certain area - analagous to an object in the field), along with a set of directions on how the object will move over some time step. This was done in the following way: set the boundary for the object's initial position, refine the coarse field around the object just as in the irregular domain case, set the new boundary conditions for the object's position after movement, subdivide the old position of the object (the hole left after movement) to fill that spot in. While this approach did provide results that might agree somewhat with intuition, there were some problems. There seems to be some discrepency at the borders of the filled-in hole and the subdivided field. This happens because the hole is filled in without any regard to the overall conditions of the field - it doesn't know that it should have to match up with the external field. There seem to be two fixes to this problem: --Extend the area being subdivided when filling in the hole in order to include more of the external field. This would hopefully force the field inside the hole to correspond with the external one - at the cost of a more time-expensive subdivision. --Implement the flow tiles as described by Stephen Chenney. The flow tiles enforce certain conditions at the boundaries such as conservation of velocity and flux. Ensuring that these quantities are conserved may help resolve the issue of inconsistent boundary conditions. In my opinion, the flow tiles method would result in more physically accurate dynamics, as the author of the flow tiles paper claims that it may be used for such a purpose. However, given the nature of the current implementation, it would be much easier to revise the code to simply include more in the domain under consideration. Hence, the first option might not produce results as good as the second, but would be easier to implement given the time constraints. Implementing the first method described above did little to alleviate the problem. As the inner domain being refined still had no knowledge of the outer boundary conditions, it failed to force the vector field to correctly match up at the edges. It seems that the only way to fix this problem would be to devise a subdivision scheme in which the inner domain being filled in were forced to match the outer boundary conditions. See below for screen shots of the problem.

Dynamic Elements

I attempted to incorporate an element of dynamics into the subdivision scheme by allowing an object to be specified in the way described above (nothing more than a set of boundary conditions that allow no flow in a certain area - analagous to an object in the field), along with a set of directions on how the object will move over some time step. This was done in the following way: set the boundary for the object's initial position, refine the coarse field around the object just as in the irregular domain case, set the new boundary conditions for the object's position after movement, subdivide the old position of the object (the hole left after movement) to fill that spot in. While this approach did provide results that might agree somewhat with intuition, there were some problems. There seems to be some discrepency at the borders of the filled-in hole and the subdivided field. This happens because the hole is filled in without any regard to the overall conditions of the field - it doesn't know that it should have to match up with the external field.

There seem to be two fixes to this problem:
--Extend the area being subdivided when filling in the hole in order to include more of the external field. This would hopefully force the field inside the hole to correspond with the external one - at the cost of a more time-expensive subdivision.
--Implement the flow tiles as described by Stephen Chenney. The flow tiles enforce certain conditions at the boundaries such as conservation of velocity and flux. Ensuring that these quantities are conserved may help resolve the issue of inconsistent boundary conditions.

In my opinion, the flow tiles method would result in more physically accurate dynamics, as the author of the flow tiles paper claims that it may be used for such a purpose. However, given the nature of the current implementation, it would be much easier to revise the code to simply include more in the domain under consideration. Hence, the first option might not produce results as good as the second, but would be easier to implement given the time constraints.

Implementing the first method described above did little to alleviate the problem. As the inner domain being refined still had no knowledge of the outer boundary conditions, it failed to force the vector field to correctly match up at the edges. It seems that the only way to fix this problem would be to devise a subdivision scheme in which the inner domain being filled in were forced to match the outer boundary conditions. See below for screen shots of the problem.

Porting to C++
The only solution that provided results which seem to agree with one's physical intuition were obtained by resubdividing the entire domain. This, however, is undesirable because it greatly increases the time complexity of the problem at hand.

Weimer and Warren have already implemented their subdivision scheme in C++. The code is freely available on their website. I modified their code, writing a more suitable front end for our purposes and deleting much of their code that had to do with various OpenGL issues that were already taken care of in Doug and Philip's code. The modified code is broken up as follows: --VectorField3d.cpp is the main data structure that stores any vector field (coarse or fine) used in the subdivision scheme. The underlying stucture is a dynamically allocated flat array that stores the x,y,z components of the field's velocity at that point by looping over the x,y,z points in the field with z as the fastest index and x as the slowest. Also, the first three values in the array store the x,y,z dimensions of the field. So, for example, if the array consists of 3 3 3 1 1 1 2 2 2 ... the vector field is a 3x3x3 field with an x component of 1 at the point (0,0,0), a y component of 1 at the point (0,0,0), a z component of 1 at the point (0,0,0), an x component of 2 at the point (0,0,1) and so forth. Also included in VectorField.cpp are the associated functions such as trimField(), which returns the center 2n-1 block of values as needed. --FlowSubdiv3D.cpp contains the actual subdivision code, which is mostly just the implementation of discrete convolution. The code in FlowSubdiv3D.cpp takes the coarse vector field and the mask and appropriately combines the two, resulting in the refined vector field. --Mask3D.cpp is much like VectorField3D.cpp, except that it is used to store the mask being used. As with VectorField3D, the mask is stored in a single flat array in the same format. Also included are the associated mask functions. --subdivide.cpp is the front end I wrote to accomodate Philip's code. The arguments provided to subdivide.cpp define the number of times to run the subdivision scheme, the size of the mask to use, the coarse vector field and the space allocated for the refined vector field. To use apply the subdivision scheme, simply include subdivide.h and call subdivide() with the appropriate arguments. Open the source code and see the comments for more information on how to use it.