Triangles on Web Ch9 Bézier Curve, B-spline Curve and Surface

Previously, we took a polygonal approach for shapes. That is, we defined the shape with straight shapes, like straight lines and triangles. However, sometimes we may need inherently curve shapes, like in vector graphics design. In this chapter, we will learn about Bézier curve and B-spline curve, which are two popular curve representations. Each two kind of two curves combined together, we get a surface.

Curve Math

In a mathematical sense, a curve is a continuous mapping from a segment of the real points in space. That is,

$$Curve(\sigma): \mathbb{R} \rightarrow \mathbb{R}^n$$

Or a more accessible form,

$$Curve(\sigma) = (x(\sigma), y(\sigma), z(\sigma)) \space \text{where} \space \sigma \in [\sigma_{min}, \sigma_{max}]$$

Where $\sigma$ is the parameter that defines the curve. For example, in a straight line, $\sigma$ can be the distance from the starting point. In a circle, $\sigma$ can be the angle.

Due to the different nature of the $\sigma$, we usually re-parameterize the curve. Re-parameterization is the process of changing the parameter $\sigma$ to another parameter $\tau$. The re-parameterization is usually a function $\tau = rep(\sigma)$.

Then the new curve is,

$$Curve(\tau) = (x(rep(\tau)), y(rep(\tau)), z(rep(\tau))) \space \text{where} \space \tau \in [\tau_{min}, \tau_{max}]$$

We often take a unified parameter called the curve-length. In an intuitive sense, the curve-length is basically how long a point have traveled on the curve from an arbitrary starting point. The curve-length is denoted as $s$.

In a more rigorous sense, the curve-length is defined as,

$$s(\sigma) = \int_{\sigma_{min}}^{\sigma} \sqrt{\left(\frac{dx}{d\sigma}\right)^2 + \left(\frac{dy}{d\sigma}\right)^2 + \left(\frac{dz}{d\sigma}\right)^2} d\sigma$$

So a curve can be represented as,

$$Curve(s) = (x(s), y(s), z(s)) \space \text{where} \space s \in [s_{min}, s_{max}]$$

Often, we prefer a normalized form. With the curve-length parameter, the unit of the parameter is unified, but it is not normalized. We can simply use a linear transformation from the curve-length to a new parameter $t$.

$$t = \frac{s - s_{min}}{s_{max} - s_{min}}$$

In all, for any given curve with any given parameter $\sigma$, we can always re-parameterize it to a normalized curve-length parameter $t$.

Bézier Curve

Then we need to consider how we can define a curve- you can't just let the user input the curve function- that is just, stupid. Currently, two popular methods are Bézier curve and B-spline curve. We will first talk about Bézier curve.

The idea of Bézier curve is simple- instead of asking user the curve function, we ask the user to define a few control points- the points that controls the shape of the curve, then we define a good curve that roughly fits the straight lines between the control points.

Let's suppose the control points are $CP_i \space \text{where} i = 0, 1, ..., n$, then Bézier curve is defined as,

$$Curve(t) = \sum_{i=0}^{n} CP_i \cdot B_{n}^{i}(t)$$

Where $B_{n}^{i+1}$ is called the Bernstein polynomial. It is defined as,

$$B_{n}^{i}(t) = \binom{n}{i} t^i (1 - t)^{n - i}$$

Where $\binom{n}{i}$ is the binomial coefficient. It is defined as,

$$\binom{n}{i} = \frac{n!}{i! (n - i)!}$$

The idea behind the Bézier curve is that, for a certain point at parameter $t$, the further the control point is, the less influence it has on the point. So we can write down,

$$Curve(t) = \sum_{i=0}^{n} CP_i \cdot Weight_i(t)$$

A natural assumption we can make is that, we want the sum of the weight stay the same throughout the curve. If we force the sum of the weight to be 1, we can write down,

$$1 = \sum_{i=0}^{n} Weight_i(t)$$

For Bézier curve, we can introduce the Bernstein polynomial in the following way,

$$1 = (1 - t + t)^n = \sum_{i=0}^{n} B_{n}^{i}(t)$$

So the weight is defined as Bernstein polynomial.

B-spline Curve

The B-spline is also based on the weight-control-point idea. In Bézier curve, since at any given point, the Bernstein weight is not zero, except for the starting and ending point, so if you make any adjustment to any control point, the whole curve will be affected. This is not always desirable.

B-spline Curve address the issue by introducing a knot vector. The knot vector is a vector that defines the weight of the control points. The B-spline curve is defined as,

$$Curve(t) = \sum_{i=0}^{n-1} CP_i \cdot N_{i, p}(t)$$

Where $N_{i, p}(t)$ is the B-spline basis function. It is defined as,

$$N_{i, p}(t) = \frac{t - t_i}{t_{i+p} - t_i} N_{i, p-1}(t) + \frac{t_{i+p+1} - t}{t_{i+p+1} - t_{i+1}} N_{i+1, p-1}(t)$$

Woa, that's scary. But don't be afraid, it is actually quite simple. The basic idea is that the instead of using the Bernstein polynomial, the B-spline function only spans a few control points, and evaluated zero outside the span. The number of points in the span is commensurate degree of the B-spline function, that is the $p$ in the $N_{i, p}(t)$.

Let's suppose we have some $CP_i$, and we want $CP_i$ to have the greatest influence at $t_{i+1}$. Considering the weight restriction, we make the following restriction on the function we construct, $N_{i, p}(t)$,

1. should be zero outside the span of $CP_i$.
2. should be one at the corresponding $t_{i+1}$ for a given $CP_i$.
3. sum of all $N_{i, p}(t)$ over $i$ should be one.

In the following, when encountering zero divided by zero, we can just take the limit as its value- it will always have one. However, actually, the values are always zero for such cases.

For the zero degree, we only consider the sole point $CP_i$. It is easy to define the weights as, one for a certain segment of $t$ and zero for the rest.

$$N_{i, 0}(t) = \begin{cases} 1 & \text{if} \space t_i \leq t < t_{i+1} \\ 0 & \text{otherwise} \end{cases}$$

So this is just a bunch of control points.

Then, we increase the degree- the span. That is, we make one control point $CP_i$ influence the next one control point, $P_{i+1}$. So we have,

$$N_{i, 1}(t) = \begin{cases} \frac{t - t_i}{t_{i+1} - t_i} & \text{if} \space t_i \leq t < t_{i+1} \\ \frac{t_{i+1} - t}{t_{i+2} - t_{i+1}} & \text{if} \space t_{i+1} \leq t < t_{i+2} \\ 0 & \text{otherwise} \end{cases}$$

This function will form a triangle spike at $t_{i+1}$, spanning from $t_{i}$ to $t_{i+2}$.

The one-zero restriction is satisfied, we check the sum now.

$$N_{i, 1}(t) + N_{i + 1, 1}(t) \\ =\frac{ t - t_{i + 1}}{t_{i + 2} - t_{i + 1}} + \frac{t_{i + 2} - t}{t_{i + 2} - t_{i + 1}} = 1$$

This sum is done for $t_{i+1} < t < t_{i+2}$.

Then, we increase the degree again. We make the control point $CP_i$ influence the next two control points. This time, instead of interpolating the $t$ value, we will interpolate the degree function. Like the following,

$$N_{i, 2}(t) = \frac{t - t_i}{t_{i+1} - t_i} N_{i, 1}(t) + \frac{t_{i+2} - t}{t_{i+2} - t_{i+1}} N_{i+1, 1}(t)$$

We know that for $N_{i, 1}(t)$, is is one at $t_{i+1}$ and zero at $t_i$ and $t_{i+2}$. Similarly, $N_{i+1, 1}(t)$ is one at $t_{i+2}$ and zero at $t_{i+1}$ and $t_{i+3}$.

So for $t_{i + 1}$, the former term is always one, while the latter stays zero.

In addition, if we sum all the $N_{i, p}(t)$ over $i$, we only have three terms to consider,

$$N_{i, 2}(t) = \frac{t - t_i}{t_{i+1} - t_i} N_{i, 1}(t) + \frac{t_{i+2} - t}{t_{i+2} - t_{i+1}} N_{i+1, 1}(t) \\ N_{i + 1, 2}(t) = \frac{t - t_{i+1}}{t_{i+2} - t_{i+1}} N_{i+1, 1}(t) + \frac{t_{i+3} - t}{t_{i+3} - t_{i+2}} N_{i+2, 1}(t) \\ N_{i + 2, 2}(t) = \frac{t - t_{i+2}}{t_{i+3} - t_{i+2}} N_{i+2, 1}(t) + \frac{t_{i+4} - t}{t_{i+4} - t_{i+3}} N_{i+3, 1}(t) \\ $$

Adding up the three terms, we have,

$$N_{i, 2}(t) + N_{i + 1, 2}(t) + N_{i + 2, 2}(t) = 1$$

The sum is done for $t_{i+2} < t < t_{i+3}$. So the remaining terms are zero.

Note that, in this calculation, the only assumption is that, $N_{i, 1} + N_{i + 1, 1} = 1$.

If we take,

$$N_{i, p} = \frac{t - t_i}{t_{i+p} - t_i} N_{i, p-1} + \frac{t_{i+p+1} - t}{t_{i+p+1} - t_{i+1}} N_{i+1, p-1}$$

The same calculation can be done for any $p$, it is just more terms to consider.

Here is a simple illustration for that. We can take $N_{i, p}$ and $N_{i+1, p}$, they will produce $N_{i, p+1}$ and $N_{i+1, p+1}$. If we iterate through, we simply get all the $N_{i, p}$. So we can conclude that the weight is always one.

The one-zero restriction is more obviously, we just skip the calculation here.

Producing Bézier Curve

This is a programming tutorial, not purely math. So we need to implement Bezier and B-spline curve. Since their only difference is the weight function, we will only implement Bézier curve for reference.

However, up until now, WebGPU does not directly support Bézier curve nor B-spline curve- since polygonal shapes are more common in computer graphics.

However, we can perform CPU calculation and gather the lines, then draw the lines with WebGPU.

That's enough rambling, let's implement the Bézier curve.

Let's review how to set up WebGPU again- this will be, sadly, the last time we do this in the series.

First, request an adapter, from which we then request a device. The adapter is obtained via navigator.gpu.requestAdapter().

const requestDevice = async (): Promise<[GPUAdapter, GPUDevice] | null> => {
    const adapter = await navigator.gpu.requestAdapter();
    if (!adapter) {
        console.error('WebGPU not supported');
        return null;
    }
    const device = await adapter.requestDevice();
    console.log(device);
    return [adapter, device];
}

Then, we set the canvas and the context. The context is obtained via canvas.getContext('webgpu').

const getContext = async (device: GPUDevice): Promise<[GPUCanvasContext, GPUTextureFormat]> => {
    const canvas = document.getElementById('app') as HTMLCanvasElement;
    const context = canvas.getContext("webgpu")!;
    const canvasFormat = navigator.gpu.getPreferredCanvasFormat();
    context.configure({
        device: device,
        format: canvasFormat,
    });
    return [context, canvasFormat];
}

Then, we set up the shader module. The vertex shader now takes a t parameter, and gives the position of the point on the curve. The fragment shader is the same as before.

We suppose we will have a bézier curve with,

$$Curve(t) = (-0.5, -0.5) B_4^0(t) + (0.5, -0.5) B_4^1(t) + (0.5, 0.5) B_4^2(t) + (-0.5, 0.5) B_4^3(t) + (-0.5, -0.5) B_4^4(t)$$

So we have the following code,

const getShaderModule = async (device: GPUDevice): Promise<GPUShaderModule> => {
    return device.createShaderModule({
        label: "shader",
        code: `
            fn factorial(n: i32) -> i32 {
                var result: i32 = 1;
                for (var i: i32 = 1; i <= n; i = i + 1) {
                    result = result * i;
                }
                return result;
            }

            fn comb(n: i32, k: i32) -> i32 {
                return factorial(n) / (factorial(k) * factorial(n - k));
            }

            @vertex
            fn vertexMain(@location(0) t: f32) -> @builtin(position) vec4<f32> {
                // calculate all the bezier functions
                var b40 = f32(comb(4, 0)) * pow(1.0 - t, 4 - 0) * pow(t, 0);
                var b41 = f32(comb(4, 1)) * pow(1.0 - t, 4 - 1) * pow(t, 1);
                var b42 = f32(comb(4, 2)) * pow(1.0 - t, 4 - 2) * pow(t, 2);
                var b43 = f32(comb(4, 3)) * pow(1.0 - t, 4 - 3) * pow(t, 3);
                var b44 = f32(comb(4, 4)) * pow(1.0 - t, 4 - 4) * pow(t, 4);
                // control points are (-0.5, -0.5), (-0.5, 0.5), (0.5, 0.5), (0.5, -0.5), (-0.5, -0.5)
                var p0 = vec2<f32>(-0.5, -0.5);
                var p1 = vec2<f32>(-0.5, 0.5);
                var p2 = vec2<f32>(0.5, 0.5);
                var p3 = vec2<f32>(0.5, -0.5);
                var p4 = vec2<f32>(-0.5, -0.5);
                var res = vec2<f32>(0.0, 0.0);
                res = res + b40 * p0;
                res = res + b41 * p1;
                res = res + b42 * p2;
                res = res + b43 * p3;
                res = res + b44 * p4;
                return vec4<f32>(res, 0.0, 1.0);
            }

            @fragment
            fn fragmentMain() -> @location(0) vec4<f32> {
                return vec4<f32>(1.0, 1.0, 1.0, 1.0);
            }
        `
    });
}

Then the vertex buffer,

const getVertexBuffer = async (device: GPUDevice): Promise<[GPUBuffer, GPUVertexBufferLayout]> => {
    const totalPoints = 100;
    const t = new Float32Array(
        Array.from(
            new Array(totalPoints - 1),
            (_, i) => i
        ).flatMap((i) => [i / totalPoints, (i + 1) / totalPoints])
    );
    const tBufferLayout: GPUVertexBufferLayout = {
        attributes: [
            {
                shaderLocation: 0,
                offset: 0,
                format: "float32"
            }
        ],
        arrayStride: 4,
        stepMode: "vertex"
    };
    const tBuffer = device.createBuffer({
        size: t.byteLength,
        usage: GPUBufferUsage.VERTEX,
        mappedAtCreation: true
    });
    new Float32Array(tBuffer.getMappedRange()).set(t);
    tBuffer.unmap();
    return [tBuffer, tBufferLayout];
}

Then, the encoder, render pass, pipeline- everything,

const main = async () => {
    const [_, device] = (await requestDevice())!;
    const [context, format] = await getContext(device);
    const shaderModule = await getShaderModule(device);

    const [tBuffer, tBufferLayout] = await getVertexBuffer(device);
    const encoder = device.createCommandEncoder();
    const renderPassDescriptor: GPURenderPassDescriptor = {
        colorAttachments: [{
            view: context.getCurrentTexture().createView(),
            clearValue: { r: 0.1, g: 0.1, b: 0.1, a: 1.0 },
            storeOp: "store",
            loadOp: "clear",
        }],
    }
    const passEncoder = encoder.beginRenderPass(renderPassDescriptor);
    const pipeline = device.createRenderPipeline({
        vertex: {
            module: shaderModule,
            entryPoint: "vertexMain",
            buffers: [tBufferLayout]
        },
        fragment: {
            module: shaderModule,
            entryPoint: "fragmentMain",
            targets: [{
                format: format
            }]
        },
        primitive: {
            topology: "line-list",
        },
        layout: "auto"
    });
    passEncoder.setPipeline(pipeline);
    passEncoder.setVertexBuffer(0, tBuffer);
    passEncoder.draw(200);
    passEncoder.end();
    device.queue.submit([encoder.finish()]);
}

Now, you can see the Bézier curve on the screen. The curve is defined by the control points, and the curve is evaluated at the given parameter $t$.

Surface

We have understood two kinds of curves- now the surface is just a combination of two curves. The surface is defined as,

$$Surface(u, v) = \sum_{i=0}^{n} \sum_{j=0}^{m} CP_{i, j} \cdot Weight_{i, j}(u, v)$$

The weight can simply be defined as the product of the weights of two base functions. That is, if you like, you can write that as,

$$Surface(u, v) = \sum_{i=0}^{n} Weight_{i}(u) \sum_{j=0}^{m} CP_{i, j} \cdot Weight_{j}(v)$$

Please note that the result of the following part

$$\sum_{j=0}^{m} CP_{i, j} \cdot Weight_{j}(v)$$

is a point. So if we fix a $v$, we get a curve. Since $u$ and $v$ are interchangeable, we can conclude that if we fix a $u$, we also get a curve.

If we choose both Bézier curve for $u$ and $v$, we have a Bézier surface. If we choose B-spline curve for $u$ and $v$, we have a B-spline surface.

So, surface is just a combination of two curves- we have drawn the curves, we can draw the surface. But since drawing 3D things requires us to, ugh, again, deal with all the perspective projection, lighting, and so on, we will not draw the surface with code in this tutorial.

Conclusion

Okay, so this whole series ends. We have learned about the basic concepts of computer graphics, WebGPU, and some advanced topics like Bézier curve and B-spline curve. We have also learned about the shader language, WGSL, and how to load 3D models.

Hope it is helpful for you.