Animating Moving Soil or Sand
Notes prepared by Hui Gong and Howard Hamilton
Introduction
We describe animation techniques suitable for illustrating the interaction of machinery and soil (dirt, sand, etc.) on the ground. This is of particular interest when creating animated simulations for the purpose of training people to use specialized construction equipment.
When implementing a virtual world where other models can interact with the ground, it is useful to converta portion of, or the entire ground surface, into many particles, which represent dirt, sand, clay and similar substances. Large particles that we can distinctly perceive are called granular material[5]. Particles that represent granular material exist in a broad range of contexts. Granular material is the second most manipulated material in industrial contexts [6].
Granular material is different from other materials, such as liquids or solid objects, because of the variety of ways in which it can behave. Multiple particles frequently interact with each other. When they are scattered in the air, they can break into pieces. When they are pushed, they can flow like a fluid. When they are lying on the ground, they act like part of the surface.To simulate all of these behaviors, it is worthwhile to find an efficient method to simulate soil treated as either ground or particles, depending on the context.
Thecritical slope angleis the maximum slope such that soil does not slip. This angle depends on the cohesion of the soil, the size of the particles, and the moisture level [7, 12].
Height Maps
A height mapor heightfieldis a 2D grid of heightswith indices corresponding to individual locations. In graphics and animation, it is typically stored as an 8-bit greyscaleimage, where each individual pixel gives the height at the corresponding (X, Y) location.If the size of a 2Dheight map is nx n,then the height map resolution is n. A 2D height map with resolution ndefines a grid of (n – 1) x (n – 1) cells (i.e., spaces between the grid points); here the n – 1 is the grid size. Height map methods are generally appropriate when the terrain is of significant size and does not have any overhangs or caves. These methods represent granular material as 2D terrain [23].
Height Map Approach to Animating Moving Sand
In the Height Map Approach, terrain is based on a height map. Granular material that is not moving is represented in the height map, and granular material that is moving is represented separately. Sandin motion is simulated using many particles, which are represented as cubes or other shapes. Piles of sand are formed by moving the particlesaround. We ensure that the surface is sloped at no more than the critical slope angle θafter the sand comes to rest. Resting sandabove the critical slope angle ismovedto nearby locations and the height of the surface is changed accordingly.
Deformation of the terrain where sand falls
Figure 3. Adding height distribution
The goal is to change the height of the terrain as sand falls. If sand falls at a single point, the goal is to form a round pile with a slope of θ. The pile should be roughly cone shaped.
Since our approach is based on a height map, the height of the terrain is represented as if it consisted of the heights of poles at individual grid coordinates, as shown in Figure 3. Imagine
Suppose sand with height H0 is added at point P0in 3D space, which is vertically aligned with coordinates (nx, ny) in the grid space. The variable wx represents the distance from the left boundary of the grid cell to the coordinates (nx, ny) in the x direction of the terrain, and wy is the distance from the bottom boundary of this grid cell to the coordinates (nx, ny) in the y direction of the terrain. Since the increased height can only be applied to the grid coordinates on the height map, we distribute the added height from position P0 to the nearest four grid coordinates, here called P00, P01, P10, and P11. To add the small height bfrom each sand particle to the height map, such that the height added to each coordinate is related to its distance from the point of impact, we use the formulas below.
convert nx from world coordinates to height map coordinate
convert ny from world coordinates to height map coordinate
wx = nx – floor(nx);
wy = ny – floor(ny);
// lower-left coordinate P00 = ( floor(nx), floor(ny) )
H00 = b * ( 1 - wx ) * ( 1 - wy );
P00.height += H00;
// upper-right coordinate P01 = ( floor(nx) + 1, floor(ny) )
H01 = b * wx * ( 1 – wy );
P01.height += H01;
// lower-right coordinate ( floor(nx), floor(ny) + 1 )
H10 = b * ( 1 - wx ) * wy;
P10.height += H10;
// upper-right coordinate ( floor(nx) + 1, floor(ny) + 1 )
H11 = b * wx * wy;
P11.height += H11;
To ensure that the sand is stable, we check that it does not exceed the critical slope angle. So after distributing the increase in height at the four coordinates as described above, we check the slope angle at each of these coordinates against the critical angle. Recall that the critical slope angle θ is the maximum slope angle for which the sand does not slip.We check the slope angle using the four height map coordinates instead of P0 because they are available in the representation.It easier to deal with the one special case of sand being added at P0 and then all subsequent cases are of moving sand between poles.
Consider the pole at position Pijwith height Hijto which we have added some sand.Suppose thatH1, ..., Hnare the heights of n nearby poles at grid coordinatesP1, ..., Pn. The slope angle ak between any pair of grid coordinates (e.g. Pij and Pk) can be calculated from the difference Hk - Hijin the height of the poles as follows:
ak = arctan,
where k = 1, 2, …, n.
Next the method repeatedly chooses the nearby point Pk for which the angle ak is greatest. If the angle ak is greater than θ, we move sand from Pij to Pk according to the amount that will fit there.If the angle ak is not greater than θ, then we stop moving sand. The pseudocode for determining slippage is as follows:
// move sand until all points are below critical angle
// or no sand remains
sand = some amount of sand to move
do {
Calculate all angles a1, a2, …, an
k = indexOfMax( /* a1, a2, …, an */ );
// slip occurs?
if(akθ) {
/* move someamount d of sandfrom Pij to Pk */
sand -= d;
}
} while(sand > 0 & akθ);
Overall, to deform terrain when sand falls, the above methodadjusts the height of the surface at each grid point in a circle around the location where the sand falls. The effect is tochange theheights of some poles near (nx, ny).
In practice, to avoid calculating the arctan repeatedly, we can approximate the task of finding the maximum angle by finding the maximum distance from the expected angle at the distance. This approximation is described shortly. As well, to ensure more uniform response time, the “while” statement in the above can be replaced by an “if” statement. If there is some sand that is not moved to the first candidate neighbor selected above, it is left at its current position.
(a) Calculation of slope angles at Pij
(b) Offset heights of neighbors
Figure 4. Calculation of slope angles based on the height map
We now describe the approach to adding sand to terrain in more detail. In Section 3.1.1, we explain the preprocessing steps and in Section 3.1.2, we give the detailed algorithm.
3.1.1. Preprocessing
1) Get the height map 2D coordinates (nx, ny) from the 3D location P0 = (x, y, z) where the sand is falling;
nx = (P0.x / WORLD_WIDTH) * (heightmap.Width- 1);
ny = (P0.z / WORLD_LENGTH) // 3D z --> 2D y
* (heightmap.Length- 1);
2) Set the fractional parts of the coordinates (nx, ny) as weight parameters wx and wy;
wx = nx - floor(nx);
wy = ny - floor(ny);
3) Let b be a variable representing the total height of all sand particles deposited at location P0in a given frame.Set bto be a fixed value related to a particle’s volume(for details see Section 4.4) when initializing the terrain to ensure that the terrain changes exactly by the same amount when each sand particle falls. Note that if b is too small, there will be no change to the height of the height map.
float b = 1 * WORLD_WIDTH/(heightmap.Width- 1);
4) Add height H00, H01, H10, and H11 to the four grid points in the height map.
3.1.2. Sand placement
Figure 5. Pile of soil with r = 2.3
To ensure that the method applies to any shape of terrain, we first distribute the height to the nearest grid points. Then we distribute it again from those points and so on. By using these local distributions, we avoid the need to consider the overall shape of the terrain.
To form a (nearly) round pile, we examine all grid points in a circle of radius r from H0, for example r = 4. The most appropriate value for r depends on the height map resolution and the detail we want to achieve. Values of less than tend to result in square piles. Let n represent the number of neighbors located in a circle of radius r. The best value for r can then be determined by solving the Gauss circle problem [29].
Figure 6. Gauss Circle Problem
As Figure 6 shows, Gauss proved that the number of integer coordinates N(r) located in a circle of radius r have the following relationship:
N(r)=r2 + E(r),
where |E(r)|2r. We define n = N(r) – 1, since we do not count the center as a neighbor of itself.
For all neighbors, we determine the angle between the top of the neighbor’s pole and the pole at the current point Pij. If any angle is greater than the critical slope angle, then move one unit height of sand toward one of the directions where the angle is greatest. Then repeat the process until no angle is greater than the critical slope angle. An alternative method would be to equally distribute the amount amongall poles where the angle is greatest.
If a sufficient number of neighbors are used, then the result will be a pile shaped roughly like a cone. Table 1 shows the number of neighbors and the resulting pile shape for a range of possible r. For example, if r is 4, then n is 46 and the shape has 20 edges.
Table 1. Radius r vs. pile shape
r / n / Base Shape/ 4 / Square
/ 8 / Square
/ 12 / Square
/ 20 / Octagon
/ 24 / Square
/ 28 / 16-gon
(16-edge Polygon)
/ 36 / Octagon
/ 42 / Octagon
/ 46 / 20-gon
(20-edge Polygon)
The algorithm for creating cone-shapedsand piles, given as function CreatePile, is as follows. The inputs include:
heightmap: a 2D height map,
(x, y): the pole position where sandhas been added (the source pole)
height: the height of the source pole
r : the radius of the circle of neighbors
1) Initialize some useful constants based on the In the code below, HORIZONTAL_CONVERSION is the conversion factor used to convert between world (terrain) coordinates (in say meters) and height maparray indices. We assume that this conversion factor is the same for both horizontal dimensions. CRITICAL_SLOPE_ANGLE is the critical slope angle θ (in radians).
WORLD_WIDTH = 200; // in meters
CRITICAL_SLOPE_ANGLE =(30.0 / 180.0 * PI);
HORIZONTAL_CONVERSION = WORLD_WIDTH / (heightmap.Width – 1);
SLOPE = tan(CRITICAL_SLOPE_ANGLE) * HORIZONTAL_CONVERSION;
2) Generate a roughly circular shape of radius r to cover all nearby integer coordinates. The poles at these locations are called the neighbors. Consider the offset height of all neighbors, where the offset height is the amount by which the height of the neighboris less than theexpected (cone-line) height of a pole at this distance from the source pole if the source pole were the center of a cone with slope θ, given that the neighbor is at its current distance from the source pole (see Figure 4(b)). Selectthe neighbor with the maximum offset height and name it the best neighbor. We will store best neighbor’s position and height in a structured variable named bestNeighbor for temporary use. The fields of the structure are named x, y, and offset. In the loop below, chr is the cone height reduction, i.e., the reduction in height that would be expected at this distance from the source pole if it were the center of a cone with slope θ.
// initialize neighbor at which maximum offset height occurs
bestNeighbor.x = x; // dummy value
bestNeighbor.y = y; // dummy value
bestNeighbor.offset = 0;
ny= floor(r);
// process all points across length of circle
for(j = y -ny; j <= y +ny; j++)
{
nx = floor( sqrt(r*r - (j-y)*(j-y)) );
// process all points across the width of circle
// at this y level of the circle
for(i = x – nx;i <= x + nx; i++)
{
//ignore any point (i,j) that is at the
// center of the circle or beyond the edge of
// the terrainor beyond the edge of the circle
if((j == y & i == x) ||
j < 0 ||
j > heightmap.Length – 1 ||
i < 0 ||
iheightmap.Width - 1||
(j - y) * (j - y) + (i - x) * (i - x) > r * r)
continue;
// assuming a cone centered at (x, y) should exist,
// determine the cone height reductionfor
// the neighbor at (i, j)based on the distance
// between it and the sourcepole at (x, y)
distance = sqrt((i-x)*(i-x) +(j-y)*(j-y) );
chr = SLOPE * distance;
expectedHeight = height – chr;
// the offset of this neighbor is the amount its
// current height is less than its expected height
curHeight = heightmap.Data[j, i] * VERTICAL_CONVERSION;
offset = expectedHeight – curHeight;
// if this offset is greater than the highest known
// offset, then update the selectionof the
// best neighbor
if(offsetbestNeighbor.offset)
{
bestNeighbor.x = i;
bestNeighbor.y = j;
bestNeighbor.offset = offset;
}
}
}
3) Determine whether or not to let the sandat this coordinate flow to the neighbor found above or not. If so, do each of the steps below using this neighbor. In the following code, c is the offset height of the neighbor found above.
assert(y >= 0 & y <= heightmap.Length – 1);
assert(x >= 0 & x <= heightmap.Width - 1);
deltaHeight = height – heightmap.GetHeight(x, y);
// positive offset?
if(bestNeighbor.offset > 0)
{
// Determine additional height to add to best neighbor.
// Only give a neighbor enough to bring it up to its
// expected height.
if (bestNeighbor.offsetdeltaHeight)
extraHeight = bestNeighbor.offset;
else
extraHeight = deltaHeight;
// increase height of pileat best neighbor
CreatePile(
bestNeighbor.x,
bestNeighbor.y,
heightmap.GetHeight(bestNeighbor.x, bestNeighbor.y)
+ extraHeight);
// any remaining height goes to (x, y) in height map
heightmap.Data[y,x] =
(height - extraHeight) / VERTICAL_CONVERSION;
}
else
{
heightmapData[y,x] = height / VERTICAL_CONVERSION;
}
4)
Texturing of the terrain where sand falls
We should also changethe textureof the terrain when sand is falling on the ground. When the height at a neighboring location is above that allowed according the maximum slope, we assume that the location isalready solid. In this case, we do not change the textures on it.
if(heightmap.GetHeight(x, y) <
heightmap.GetHeight(bestNeighbor.x, bestNeighbor.y) + chr)
// process points (i,j) around current point (x,y)
for(j = -1; j <= 0; j++)
for(i = -1; i <= 0; i++)
for(k = 0; k < heightmap.alphaMapLayers; k++)
{
// is sand falling outside of the map?
if(y + j < 0 ||
x + i < 0 ||
y + j >= heightmap.Length – 1 ||
x + i >= heightmap.Width – 1)
continue;
// change texture layer number to the maximum,
// if it represents the texture of sand;
// otherwise, change it to the minimum
if(k == 0)
splatMapData[y + j, x + i, k] = 1;
else
splatMapData[y + j, x + i, k] = 0;
}
3.3.4. PouringSand
To achieve realistic effects, some parameters of the sandparticles are adjustedbefore the particles hit the ground. When particles fall from the edge of the bucket, the speed in the XOZ plane is set to a smaller value, producing a visual result that seems more like sand than ice balls. At the moment that a falling cube hits the ground, the height of its highest point is at most, wherebis the side length of the cube. To avoid recalculating the square root, is represented as 0.86 in the code.The movement of the sand particle is implemented in the FixedUpdate function since a smaller time step is needed to increase accuracy of the movements. The code follows below.
// NOTATION IN THE FOLLOWING PARTS OF THIS PAPER HAS
// NOT YET BEEN MADE CONSISTENT WITH THE FIRST PART
#define LEVEL_VEL_CONS 1.1f
#define DOWNWARD_VEL_CONS -10.0f
void drop()
{
// process each sand particle
foreach(GameObject d in mylist.ToList())
{
d.rigidbody.velocity = new Vector3
(d.rigidbody.velocity.x / LEVEL_VEL_CONS,
d.rigidbody.velocity.y,
d.rigidbody.velocity.z / LEVEL_VEL_CONS);
// if particle is moving upwards, change it to moving
// downwards
if(d.rigidbody.velocity.y > 0)
d.rigidbody.velocity =
new Vector3(0, DOWNWARD_VEL_CONS, 0);
//convert the 3D world position of object d to
// 2D terrain height map coordinates
Vector3 p = d.transform.position;
intnx = wtonx(p);
intny = wtony(p);
// point outside of terrain?
if(nxMAPGRIDSIZE || nyMAPGRIDSIZE)
continue;
// find the height of the terrain at the neighboring
// grid point that has the greatest height and the
// height of the lowest part of the bucket
float max = 0; // maximum height
float min = MAX_TERRAIN_HEIGHT; // minimum height
for(int j = 0; j <= 1; j++)
for(inti = 0; i <= 1; i++)
{
tmp = heightmapData[ny + j, nx + i] *
td.heightmapScale.y;
if (tmp > max) // highest neighboring point
max = tmp;
tmp = low[ny + j, nx + i]; // lowest part of bucket
if (tmp < min)
min = tmp;
}
// if the sandparticle is sufficiently far under the bucket
// (i.e., it is no longer touching the bucket) and it is
// touching the terrain, then it becomes part of the
// terrain. The maximum height of a cube with side length
// b stood on its corner is sqrt(3)/2 * b.
const float SQRT_3_OVER_2 = 0.86;
if ((p.y= min –SQRT_3_OVER_2 * b)
(p.y <= max + SQRT_3_OVER_2 * b))
{
hit(p);
Destroy(d);
mylist.Remove(d);
}
}
}
3.3.5. Hitting
To simulate the action of fallen soil becoming a part of the ground, we remove the particles representing this soil from the simulation, as described at the end of the previous section,and adjust the terrain height map at the exact positions of the particles. Here we describe the details of the hit function that performs the adjustment.
We apply a hit function to every position where a sand particle lands on the ground.The hit functiondeforms the ground slightly for each particle. It implements the approach described in Section 3.1. The only parameter used is the (x,y,z) position of the hit.
void hit(Vector3 pos)
{
// compute standard local positions and offsets
Vector3 terrainLocalPos = pos - terr.transform.position;
float nx = Mathf.InverseLerp(0, td.size.x, terrainLocalPos.x)
* (td.HeightMapWidth - 1);
float ny = Mathf.InverseLerp(0, td.size.z, terrainLocalPos.z)
* (td.HeightMapLength - 1);
float fx = Mathf.Floor(nx);
float fy = Mathf.Floor(ny);
float wx = nx - fx;
float wy = ny - fy;
// deform each neighboring grid point in the height map
for(int j = 0; j <= 1; j++)
for(inti = 0; i <= 1; i++)
{
// compute increased height resulting from impact
float addHeight = (float)((1 - i) + (2 * i - 1) * wx) *
((1 - j) + (2 * j - 1) * wy) * b;
float oldHeight =
heightmapData [(int)fy+j, (int)fx+i] * td.heightmapScale.y;
// add increased height to current point
float[,] newHeightData =
new float[1, 1] { { oldHeight + addHeight } };
CreatePile((int)fx+i , (int)fy+j, newHeightData[0, 0]);
}
}
4.4. Volume conversion
When the sand is dug up, a unit of sand is converted from terrain to sand particles. Since the change of terrain is calculated by changing the height map, we can get the total changed volume by accumulating all the changes from every 2D-coordinate of the height map. For each change in a coordinate, we obtain the volume change by considering the nearest six coordinates as the terrain height map using a tree structure to represent the change.