Abstract

This article introduces the basic principles of double cubic sampling and its application in three-dimensional terrain generation, along with a code implementation.

References

[1] wikipedia: bicubic interpolation
[2] wikipedia: perlin noise
[3] Finite Difference Approximations of Derivatives

Terrain Generation

I have recently been researching three-dimensional terrain generation, which should be one of the most common applications in procedural generation. The essence of terrain generation is generating heightmaps, and generally, more complex schemes will synthesize multiple noise maps[2], allowing each noise map to control different details. One troublesome issue is that if the generated noise maps have different resolutions, the smaller noise map needs to be upsampled, or the larger map needs to be downsampled. The common mathematical method for upsampling is interpolation, and the two most intuitive interpolation methods are nearest-neighbor interpolation and bilinear interpolation. Bicubic sampling is not as intuitive, but it is the only one of the three methods that guarantees continuity, meaning visually, you will not see abrupt changes.

Taking the two figures above as an example, when the number of interpolation points far exceeds the number of sampling points, bilinear interpolation will produce noticeable segmentation lines, while bicubic interpolation can completely avoid this problem. However, bilinear interpolation can also mask this issue when the number of interpolation points is comparable to that of the sampling points.

Mathematical Principles

Each interpolation point is influenced by 16 surrounding sample points in the grid (the 4 closest to $(x, y)$ and the outer ring of 12). The final evaluation formula uses the definition of a cubic function in two dimensions, which is detailed on Wikipedia, so it will not be expanded here. The interpolation evaluation formula is as follows:

$$ p(x,y) = \sum_{i=0}^3 \sum_{i=0}^3 a_{ij} x^i y^j $$

To determine $ a_{ij} $, we must also list the following continuity constraints: for the four points inside, their interpolation $p(x, y)$ and sampled values $f(x, y)$ must have the same derivatives in the x, y, and the mixed x-y directions. With these constraints, all the $a_{ij}$ can be expressed using four sampled values and their derivatives in the three directions.

@x

In the discrete case:

$$ f_x(x, y) = \frac{f(x - 1, y) + f(x + 1, y)}{2} $$

$$ f_y(x, y) = \frac{f(x, y - 1) + f(x, y + 1)}{2} $$

$$ f_{xy}(x, y) = \frac{f(x + 1, y + 1) - f(x - 1, y + 1) - f(x + 1, y - 1) + f(x - 1, y - 1)}{4} $$

The value of $f_{xy}$ is not very intuitive; however, it can be derived from the approximation of second-order discrete derivatives [3]. The final interpolation formula is:

@y

C# Implementation

In engineering implementation, all three derivatives of the sampling points should first be calculated and stored, and then the interpolation positions should be traversed to reduce computation. For sampling points outside the effective boundary, the nearest effective sampling points can be used as substitutes. Since interpolation calculations can be represented using matrix operations, if there are packages that can call SIMD instructions (such as Unity.Mathematics here), the computation can be greatly accelerated.

using Unity.Mathematics;
using UnityEngine;

public class BicubicSampler {
    public float[] data;
    public float4x4[] B;
    public int width;
    public int height;

    public BicubicSampler(float[] data, int width) {
        this.data = data;
        this.width = width;
        height = data.Length / width;
        B = new float4x4[(height - 1) * (width - 1)];
        SetUp();
    }

    public float F(int x, int y) {
        int _x = (x < 0) ? 0 : ((x >= width) ? width - 1 : x);
        int _y = (y < 0) ? 0 : ((y >= height) ? height - 1 : y);
        return data[_y * width + _x];
    }

    public float Fx(int x, int y) {
        return (F(x + 1, y) - F(x - 1, y)) * 0.5f;
    }

    public float Fy(int x, int y) {
        return (F(x, y + 1) - F(x, y - 1)) * 0.5f;
    }

    public float Fxy(int x, int y) {
        return (F(x + 1, y + 1) - F(x - 1, y + 1) - F(x + 1, y - 1) + F(x - 1, y - 1)) * 0.25f;
    }

    public void SetUp() {
        let left = new float4x4(1, 0, 0, 0, 0, 0, 1, 0, -3, 3, -2, -1, 2, -2, 1, 1);
        let right = new float4x4(1, 0, -3, 2, 0, 0, 3, -2, 0, 1, -2, 1, 0, 0, -1, 1);

        for (int y = 0; y < height - 1; y++) {
            for (int x = 0; x < width - 1; x++) {
                let sample = new float4x4(
                    F(x - 1, y - 1), F(x - 1, y), Fy(x - 1, y - 1), Fy(x - 1, y),
                    F(x, y - 1), F(x, y), Fy(x, y - 1), Fy(x, y),
                    Fx(x - 1, y - 1), Fx(x - 1, y), Fxy(x - 1, y - 1), Fxy(x - 1, y),
                    Fx(x, y - 1), Fx(x, y), Fxy(x, y - 1), Fxy(x, y)
                );
                B[y * (width - 1) + x] = math.mul(left, math.mul(sample, right));
            }
        }
    }

    public float Interpolate(Vector2 pos) {
        float _x = pos.x * (width - 1);
        float _y = pos.y * (height - 1);
        int idX = (int)_x;
        int idY = (int)_y;
        float u = _x - idX;
        float v = _y - idY;
        let left = new float4(1, u, u * u, u * u * u);
        let right = new float4(1, v, v * v, v * v * v);
        let res = math.mul(left, math.mul(B[idY * (width - 1) + idX], right));
        return res;
    }
}