Overview
This was one of my first procedural tools. As many beginners do I looked for a tutorial and found one from SideFX and Ben Schrijvers, principal technical artist at Guerilla games. That tutorial served as my gateway into Houdini, as well as providing an opportunity to learn the underlying math.
While certain elements of that tutorial remain unchanged in the final version of my tool, I simplified aspects and restructured the tool to leverage PDG. First I separated the tool into two components, a controller and a processor.
The controller is a series of primitive attributes attached to a curve, which the user can manipulate to change the path and properties of the river. The processor takes a river controller and generates two outputs: a riverbank to be projected into the landscape, and a river object with a surface, volume, and spline to be exported to Unreal.
Using separate HDAs to process controllers and construct the geometry respectively allows them to be easily looped over and multithread in SOP use-cases, and partitioned into separate work items in TOPs.
Pathfinding
The first part of this tool involves processing the control curve, and turning it into a path for a river. Pathfinding is broken up into three sections: winding, erosion, and cascades, which when combined create naturalistic river paths.
Winding
The function starts by sampling the height and gradient of the heightfield at the position of each point along the guide curve, then sets the position and normal point attributes of the curve respectively. To get the normal from the heightfield gradient and up vector are added together.
//sample the y position of the terrain given the x and z position of the curve point
@P = set(@P.x, volumesample(1, 0, @P), @P.z);
//get the gradient of the heightfield, and blend it with up vector
//x and z match heightfield, y faces up
@N = normalize(volumegradient(1, 0, @P)+{0,1,0});Then, with a bit of math, it gets the vector that points along the downward slope of the terrain. The cross product of two vectors returns a perpendicular vector. So the cross product of the normal and up vector returns a horizontal vector perpendicular to the normal. Since the normal is sampled from the heightfield, the cross product of the normal and horizonal vector returns a vector representing the downward slope or gradient of the heightfield. Adjusting the magnitude of the downward slope and adding it to the initial position of the point causes the curve to "slide" down the terrain.
//without changing the point position, get a vector that points sideways relative to the terrain
vector get_sideways_vector_on_terrain(vector n_vec; vector normal) {
return normalize(cross(normal, cross(normal, {0,1,0})));
}
//get user parameters
float winding_bias = chf("winding_bias");
//get vector pointing along direction of curve
int next_point = @ptnum+1;
vector next_point_pos = point(0, "P", next_point);
vector n_vec = normalize(next_point_pos-@P);
//add the positions of the downard slope and the current point position to get new point position
@P = normalize(get_sideways_vector_on_terrain(n_vec, @N)*winding_bias)+@P;The curve is then resampled using a separate winding resolution. This allows the user to control the size and density of the resulting bends. A smaller value means more, smaller bends along the river. If this value is smaller than the master resolution of the river, then the detail will ultimately be lost.
Iterating over this function causes each point in the curve to “slide” along the downward gradient of the heightfield This results in narrow straight rivers forming where the heightfield has steep narrow canyons, and wide, winding rivers to form where the heightfield is flattest.
If you followed the tutorial I included in the introduction, this will look very familiar. However, this method of pathfinding comes with a severe performance problem. Since the curve will eventually settle into a local optimum, any iterations of the function beyond the optimum yield no results. However they still carry a performance cost. In essence, the tool wastes time to accomplish nothing. This performance cost is relevant for rivers with lots of points.
Luckily, this problem is easily solved with a bit of simple calculus. Each iteration of the pathfinding function can be mapped on one axis, and the change in y-position of the point mapped on another axis. The derivative of each iteration is equal to the change in y-position relative to the previous iteration. As the derivative approaches zero, we know the point has achieved a local optimum. Once a sufficient percentage of points have achieved their local optimum, the loop is stopped. The final pathfinding script looks like this:
//get vector representing downward slope of terrain
vector get_sideways_vector_on_terrain(vector normal; vector up; vector dir) {
vector side_vec = normalize(cross(normal, cross(normal, up)));
//hold points back if they are sliding towards the next point
//helps prevent self-overlap
if (chi("limit_bends") == 0)
{
side_vec -= dot(dir, side_vec)*dir;
}
return side_vec;
}
//sample the y position of the terrain given the x and z position of the curve point
@P = set(@P.x, volumesample(1, 0, @P), @P.z);
//store initial y position to get derivative later
float init_y=@P.y;
//get the gradient of the heightfield, and blend it with up vector
//x and z match heightfield, y faces up
@N = normalize(volumegradient(1, 0, @P)+{0,1,0});
//get user parameters
float step_length = chf("winding_step_length");
vector up = chv("up");
//get vector pointing along flow direction of curve
int next_point = @ptnum+1;
vector next_point_pos = point(0, "P", next_point);
vector n_vec = normalize(next_point_pos-@P);
//slide point down the terrain
@P = @P-normalize(get_sideways_vector_on_terrain(@N, up, n_vec)*step_length);
//sample the y position of the terrain given the x and z position of the curve point
@P = set(@P.x, volumesample(1, 0, @P), @P.z);
//create derivative attribute
f@deriv = init_y-@P.y;
Another attribute wrangle uses the deriv attribute to determine is a sufficient number of points are at a local optimum. The user sets two parameters, derivative threshold and optimal threshold. Derivative threshold is the maximum derivative that a point can be considered optimal. The smaller the derivative threshold, the more iterations will be needed to find a local optimum. Optimization threshold is the percentage of total points on the curve which need to be optimal to break the loop. It is a normalized range, so a value of 1 requires the entire curve to be optimal to stop the loop.
//attribute being referenced for stop condition on loop
//by default loop continues
i@optimized=0;
//track number of points at local minima
int loc_min = 0;
//if the derivative approaches zero
for (int p = 0; p<@numpt; p++) {
if (abs(point(0, "deriv", 0)) <= ch("derivative_threshold")) {
//increment total optimal points
loc_min+=1;
}
}
//if enough points are optimized, the loop stops
if (float(loc_min)/float(npoints(0)) >= ch("optimization_threshold")) {
i@optimized = 1;
}Finally, I added a toggle to disable winding if the user wishes. So the final graph and parameter interface looks like this.
Erosion
Simulating erosion is a simple process. In an attribute wrangle node, each point is set to be equal to or lower than the point position. There are two inputs, erosion and slope. Erosion is a value between 0 and 1 which interpolates between the original point positions and their eroded positions. By default erosion is set to 1. Slope determines the steepness of the river. It is the absolute value in meters by which each point will be lowered relative to the previous point. Higher values will cause the river to flow downwards at a steeper angle.
float erosion = ch("erosion");
float slope = ch("slope");
float FLT_MIN = -99999999999999999999999;
float FLT_MAX = 1.7976931348623157e+308;
float min_height = FLT_MAX;
float new_height = 0.0;
float clamp_in = FLT_MAX;
float clamp_out = FLT_MIN;
for (int i = 0; i < @numpt; i++){
float current_x = getcomp(vector(point(0, "P", i)), 0);
float current_y = getcomp(vector(point(0, "P", i)), 1);
float current_z = getcomp(vector(point(0, "P", i)), 2);
min_height -= slope;
min_height = min(min_height, currentY);
new_height = (min_height * erosion) + (currentY * (1-erosion));
vector p = set(currentX, new_height, currentZ);
setpointattrib(0,"P",i,p);
}
Scale by Slope
This is a very simple method to scale the river depending on its slope. Ideally, a river should be wider where it is flat, and narrower where it is steep. This code in a point attribute wrangle measures the steepness of the river at any given point via the dot product of the normal and up vector, then scales it according to a user defined range, where the minimum value indicates the scale of the river at its steepest.
vector up = {0,1,0};
v@up=up;
float min_scale = chf("min_scale");
float max_scale = chf("max_scale");
//get dot product between the direction of the river and up vector
//remap dot product to alpha for lerp
float dotp = dot(up, @N);
dotp = fit(dotp, -1,0, 0,1);
//width
float width_alpha = float(@ptnum)/float(@numpt);
float width = chramp("width_ramp", width_alpha);
width = fit(width, 0,1,0,10);
float scale = lerp(min_scale, max_scale, dotp);
scale*=width;
v@pscale = set(scale, scale, scale);Cascades
In nature, rivers do not maintain an even slope across their length. Sections with hard sediment and rock erode more slowly than softer sections with more sediment. This causes cataracts and waterfalls to form at the edges of rock formations. I developed two methods to emulate cascades.
Height Threshold Method
In the first method, the user controls the minimum and maximum height of the cascades. A vex script then splits the river into segments by their total slope, flattening each segment to its lowest respective point. The end point of each cascade is omitted from the subsequent cascade, pinning the point in place. This forms a cascade at the beginning of each segment.
void cascade_by_range (int geo; float steepness, min_slope, max_slope; int seed, numpt) {
float slope_threshold_alpha = rand(seed);
float slope_threshold = lerp(min_slope, max_slope, slope_threshold_alpha);
/*
cascade_start: point number for first point in cascade
cascade_index: tracks how many cascades have been formed.
min_y: lowest point in the current cascade array
start_y: y position of the first point in the cascade
point_y: y position of the current point in the loop
*/
int cascade_start = 0;
int cascade_index = 0;
float min_y;
float start_y;
float point_y;
//array to store point number of segments of the river which will be flattened
int cascade_points[];
//for each point in curve
for (int i = 0; i < numpt; i++) {
//get the y position of the current point
point_y = getcomp(vector(point(geo, "P", i)), 1);
//get the y level of the start of the most recent cascade
start_y = getcomp(vector(point(geo, "P", cascade_start)), 1);
//if the current point y is sufficienly below the start of the current cascade
if ((start_y - point_y > slope_threshold) && (len(cascade_points)>1)) {
//make sure the cascade attribute starts at 0
if (cascade_index != 0) {
setpointattrib(geo, "cascade", cascade_points[0], cascade_index-1);
setpointattrib(geo, "cascade", cascade_points[0]-1, cascade_index-1);
}
//get tops of casacdes for particle emission
setpointattrib(geo, "cascade_emitter", i, 1);
cascade_index+=1;
/*
for each cascade point lerp between the
lowest y position in the cascade and current y position
using steepness as alpha
*/
//get lowest y position out of the cascade
min_y = getcomp(vector(point(geo, "P", i)), 1);
foreach(int num; cascade_points) {
if (num==cascade_points[0]) {
//blend first point in cascade
//control steepness of cascades
vector parent_point = point(geo, "P", cascade_start);
vector current = point(geo, "P", num);
vector p = set(
lerp(current[0], parent_point[0], steepness),
min_y,
lerp(current[2], parent_point[2], steepness));
setpointattrib(geo, "P", num, p);
} else {
vector current = point(geo, "P", num);
vector p = set(
current[0],
min_y,
current[2]);
setpointattrib(geo,"P",num,p);
}
//update point attribtes separating waterfalls and flat sections
setpointattrib(geo, "segment", num, cascade_index);
setpointattrib(geo, "cascade", num, -1);
}
//setpointattrib(0, "cascade_top", cascade_start, 1);
//clear cascade points array
cascade_start = i;
resize(cascade_points, 0);
//get new random seed and sample slope threshold range again
seed *= 2;
slope_threshold_alpha = rand(seed);
slope_threshold = lerp(min_slope, max_slope, slope_threshold_alpha);
} else {
if (i != cascade_start) {
push(cascade_points, i);
}
}
}
}Slope Method
This method is requires a slope threshold input from the user. Like the height threshold method, this tool loops through all the points in the curve, comparing the remapped dot product of the flow direction and up vector to a slope threshold input. I start by generating a flow attribute.
//current point position
vector position = @P;
//position of next point in curve
vector next_position = point(0, "P", @ptnum+1);
v@flow_dir = position-next_position;
//fix last point in curve
if (@ptnum == @numpt-1) {
v@flow_dir = normalize(point(0, "P", @ptnum-1) - @P);
}If the absolute value of the dot product of the up and flow vectors is greater than the slope threshold, then a cascade is formed. Each subsequent point that also exceeds the slope threshold will be appended to the cascade. Once the slope is sufficiently flat, again, the cascade ends. All the points in the cascade align their x and z positions, forming a cascade. Additional parameters control the steepness and position of the cascade along its length.
The benefits of the slope method over the range method is that cascades of any size may form, conforming to the input landscape. There is less need to tweak the settings or seed value to get a good result. Additionally, because it does not adjust the y-position of any of the points, it retains the slope generated by the erosion function between cascades.
void cascade_by_slope (int geo; float slope_threshold, position_alpha, steepness; vector up; int numpt) {
//variables
/*
cascade_start/end: start and end positions of the cascade respectively
cascade_position: x and z position of the cascade
current_position: initial position of the current point in the cascade loop
new_position: updated position of the current point in the cascade loop
n_1: normal of first point in cascade
n_2: normal of last point in cascade
cascade_points: array which holds all points in a cascade
cascade_index: int id for each cascade
segment_index: int id for each segment
*/
vector flow;
vector cascade_start;
vector cascade_end;
vector cascade_position;
vector current_position;
vector new_position;
vector n_1;
vector n_2;
int cascade_points[];
int cascade_index = 0;
int segment_index = 0;
for (int i = 0; i < numpt; i++) { //for each point
//flatten the flow of the river horizontally
flow = point(geo, "flow", i);
//if the point has a sufficient downward
//append it to the cascade array
if (abs(dot(flow, up)) > slope_threshold)
{
push(cascade_points, i);
}
else
{
if (len(cascade_points)>1) {
//set variables for new cascade
cascade_start = point(geo, "P", cascade_points[0]);
cascade_end = point(geo, "P", cascade_points[-1]);
cascade_position[0] = lerp(cascade_start[0], cascade_end[0], position_alpha);
cascade_position[2] = lerp(cascade_start[2], cascade_end[2], position_alpha);
n_1 = point(geo, "N", cascade_points[0]);
n_2 = point(geo, "N", cascade_points[-1]);
//cascade loop
foreach(int num; cascade_points)
{
//store cascade number
setpointattrib(geo, "cascade", num, cascade_index);
//update position
current_position = point(geo, "P", num);
new_position = set(
lerp(current_position[0], cascade_position[0], steepness),
current_position[1],
lerp(current_position[2], cascade_position[2], steepness)
);
setpointattrib(geo, "P", num, new_position);
//fix normals
setpointattrib(geo, "N", num, lerp(n_1, n_2, 0.5));
}
cascade_index +=1;
segment_index+=1;
//clear array
resize(cascade_points, 0);
} else {
setpointattrib(geo, "cascade", i, -1);
setpointattrib(geo, "segment", i, segment_index);
resize(cascade_points, 0);
}
}
}
}Skin Riverbank
Now that the guide curve is complete, generating geometry for the riverbank and riverbed is a relatively simple process. A sweep node generates a series of columns using a line node to control the horizontal resolution. The resulting columns are then separated into two groups, riverbank and riverbed based. The user controls the proportions of the river using a water inset value. A user then controls the profile of the bank and bed, as well as the blend with the landscape via an alpha point attribute. The riverbed is then projected into the landscape.
The water surface likewise follows the guide curve. Knowing the maximum width of the river and the water inset, I sweep a flat surface. Separating this process from the riverbank geometry allows finer control over the resolution of the geometry. I generate two uv sets for the water surface. One is a normalized set, the other is world aligned. Having both sets of uvs is more a convenience than a necessity. It merely allows the user to switch between different shading methods quickly instead of having to re-bake the geometry.
In the river controller, a string attribute is created containing the primitive numbers of the columns which form the bank and bed of the river respectively. This is done with an attribute wrangle, a for-loop, and a bit of simple maths. These values are based on the width resolution of the river, and a water inset value which determines the transition point between the bank and bed.
