Introduction

Last day, we started talking about Bezier curves and Catmull Rom Splines again.
Our starting motivation was smooth surface generation as that which might be obtained by Perlin noise terrain generation or which might be obtained by the random displacement method.
When we talked about these curves last semester we were interested in smoothly animating motion along curves.
My notes for today come from another course so for this lecture I will tend to write `vec p` for points `tilde p` using our book's notation.
We first talked about degree 3, 1-dimensional Bezier Curves and at the end said one could have 3D Bezier curves with equations like:
`\vec q(u) = B_0(u)\vec p_0 + B_1(u)\vec p_1 + B_2(u)\vec p_2 + B_3(u)\vec p_3`
where `B_i(u) = ((3),(i))u^i(1-u)^(3-i) ` are called blending functions.
Blending functions sum to 1:
`\sum_{i=0}^3B_i(u) = \sum_{i=0}^3((3),(i))u^i(1-u)^(3-i) = (u+(1-u))^3 = 1`
So `\vec q(u)` is always a weighted average of the four control points. From the definition of the `B_i` we also have `\vec q(0) = \vec p_0` and `\vec q(1) = \vec p_3`.
To do terrain generation, one problem to solve is smoothly interpolation a curve through a set of points, `vec p_0, \ldots, vec p_m`.
The solution we have for this was to use the translation of degree three curves and to make sure the end points satisfied some "smoothness" condition on the first derivatives. The resulting curve is a sequence of Bezier curves with domains `[0,1]`, `[1,2]`, ..., `[n-1, n]`, so with overall domain `[0,n]` called a Catmull-Rom spline.
The points `0`, `1`, `2`, ..., `n` are called the "knots" of this spline and are evenly spaced. In general, we could imagine splitting our domain interval `[0, n]` at non-uniformly spaced knot points `u_i` and ask that our function `vec q(u)` satisfy `vec q(u_i) = vec p_i`.
We start today with some problems evenly space knots cause, and a brief introduction on how to solve them.
To do the random displacement method, we might want to split a curve into parts quickly.
To make terrain we also want to say how to generalize things to surfaces rather than just curves.

Bessel-Overhauser Splines

example of not Catmull-Rom Splines with not evenly spaced control points

Let `vec p_i` be the points we are interpolating through.
Recall the control point of the `i`th curve which made up a Catmull-Rom spline were given by
`vec p_i`, `vec p_i^+`, `vec p_(i+1)^-`, and `vec p_(i+1)` where
`vec p_i^+ = vec p_i +1/3 vec L_i` and `vec p_i^(-) = vec p_i -1/3 vec L_i` where `vec L_i = 1/2(vec p_(i+1) - vec p_(i-1))`.
So `\vec L_i` is defined in terms of the `i-1`st and `i+1`st control points. Because of this we can only define `\vec p_i^+` and `\vec p_i^-` when `0 < i < m`, so for Catmull-Rom Splines the first and last points are there to make the equations of the intermediate piecewise Bezier curves work out, but will not be interpolated.
With Catmull-Rom splines, when knot points are not evenly spaced, you can get bad "overshoot" effects as seen above when two close control are next to widely separated control points.
To solve this problem, one can use a chord-length parametrization. That is, we choose the knots so that `u_(i+1) - u_i = ||\vec p_(i+1) - \vec p_i||`.
This will roughly make the arclength of the curve between `\vec p_i` and `\vec p_(i+1)` proportional to `u_(i+1) - u_i`.

More Bessel-Overhauser Splines

To effect this parametrization, we need to modify our definitions of `\vec p_i^-` and `\vec p_i^+` using the Bessel tangent method (aka the Overhauser method).
Suppose we have a set of knots `u_i`. We can define the average velocity of the curve between `\vec p_i` and `\vec p_(i+1)` as
`\vec v_(i+1/2) = \frac{\vec p_(i+1) - \vec p_(i)}{u_(i+1) - u_(i)}`.
If the chord length approximation was good, then these will be unit vectors, and we can approximate the velocity at `\vec p_i` as weighted average:
`\vec v_i = \frac{(u_(i+1) - u_(i))\vec v_(i+1/2) + (u_(i) - u_(i-1))\vec v_(i-1/2)}{u_(i+1) - u_(i-1)}`
Finally, we can define the Bessel-Overhauser Splines, in the same way as Catmull-Rom splines but using `\vec p_i^+` and `\vec p_i^-` defined by:
`\vec p_i^(-) = vec p_i - \frac{1}{3}(u_i - u_(i-1))\vec v_i`
`\vec p_i^(+) = vec p_i + \frac{1}{3}(u_(i+1) - u_(i))\vec v_i`

Recursive Subdivision

It is much easier to draw lines than curves.
So it useful to be able to come up with techniques to approximate Bezier curves with line segments.
Splitting up curves is also useful for the random displacement method.
In order to do this, we will develop techniques to split an existing Bezier curve into two Bezier curves.
If `\vec q(u)` is given by `\vec p_0`, `\vec p_1`, `\vec p_2`, `\vec p_3`, then one way to define two curves that together give all of `\vec q(u)` is to set:
`\vec q_1(u) = \vec q(u/2)` and `\vec q_2(u) = \vec q((u + 1)/2).`
`\vec q_1(u)` and `\vec q_2(u)` are cubic because we substituted a linear factor in `\vec q(u)` which was cubic, but are they Bezier?
Yes! You can mechanically verify that `\vec q_1(u)` is the Bezier curve with control points: `\vec p_0`, `\vec r_0`, `\vec s_0`, `\vec t_0` and `\vec q_2(u)` is the Bezier curve with control points: `\vec t_0`, `\vec s_1`, `\vec r_2`, `\vec p_3` where the `r_i`, `s_i`, and `t_i` are the interpolants from de Casteljau's method:
`\vec r_i = (1-u)\vec p_i + u\vec p_(i+1).` (in this case we want to take `u=1/2`.)
`\vec s_i = (1-u)\vec r_i + u\vec r_(i+1).`
`\vec t_0 = (1-u)\vec s_0 + u \vec s_1`.
What if you wanted to make the split into two curves at point `\vec q(u_0)`? Instead, of `\vec q_1(u)` and `\vec q_2(u)` as above, you set
`\vec q_1(u) = \vec q(u_0u)` and `\vec q_2(u) = \vec q(u_0 + (1 - u_0)u).`
if we define `t_0 = q(u_0)` and define the `r_i`'s and `s_i`'s using `u_0` then `\vec q_1(u)` and `\vec q_2(u)` will be given by `\vec p_0`, `\vec r_0`, `\vec s_0`, `\vec t_0` and `\vec t_0`, `\vec s_1`, `\vec r_2`, `\vec p_3`.

Applications of Recursive Subdivision

Given a Bezier curve, one way to draw it is to recursively split the curve into two halve curves and draw each of the halves.
We keep splitting until the difference between the split Bezier curve and a straight line is less than some error `\delta`.
We then draw a straight line for the Bezier curve.
The Bezier curve is contained within the region specified by its four control points.
So we thus want the approximating line segment `\bar(p_0p_3)` to be close to the points `p_1` and `p_2`.
A quick and usually correct way to check for this is to see if `||\vec q(1/2) - 1/2(p_0 + p_3)|| < \delta`
One can check `\vec q(1/2) = 1/8 \vec p_0 + 3/8 \vec p_1 + 3/8 \vec p_2 +1/8 \vec p_3`.
Substituting this means checking that:
`||\vec p_0 - \vec p_1 - \vec p_2 - \vec p_3||^2 < (\frac{8\delta}{3})^2`.

Another Application of Recursive Subdivision

Another use of recursive subdivision is intersection testing. i.e., You might want to see check if some Bezier surface is outside the viewing frustum and so does not need to be rendered.
You can check if a ray intersects with a Bezier curve by checking first if the ray intersects with the curves convex hull. If no, return false.
If yes, split the Bezier curve into two halves and check for intersection on both halves.
The recursion continues until we get to the line segment approximation level.
Finally, if the ray and line segment intersect we actually output yes for the whole intersection test.

Bezier Surface Patches

A picture of the degree three Bezier patch and its sixteen control points

A Bezier curve is a one-dimensional parametric curve; a Bezier patch is a two-dimensional surface.
Such a patch is parametrized using two variables `u` and `v` each ranging over `[0, 1]`.
The patch is then the parametric surface `\vec q(u,v)` where `\vec q` is a vector-valued function defined on the unit square.
A degree three Bezier patch is given by a 4 x 4 array of control points `\vec p_(i,j)`, where `i`, `j` take on values `0, 1, 2, 3`.
The patch defined by these control points is given by the formula:
`\vec q(u,v) = \sum_(i=0)^3\sum_(j=0)^3B_i(u)B_j(v)\vec p_(i,j)`.

More on Bezier Patches

Notice if you hold `u` fixed in the equation for `\vec q` and let `v` vary, the resulting curve is a Bezier curve with control points
`\vec s_j = \sum_(i=0)^3B_i(u)\vec p_(i,j)`.
You also get Bezier curves if you hold `v` fixed and let `u` vary.
Since we eventually want to be able to join patches together in order to take complex surfaces, knowing the partial derivatives of a Bezier patch at its control points is useful.
The holding one parameter fixed property we just noticed allows us to figure out the partial derivatives for Bezier patches using our knowledge of the derivatives of Bezier curves.
Namely, we have:
`\frac{\partial \vec q}{\partial v}(u, 0) = \sum_(i=0)^3 3B_i(u)(\vec p_(i,1) - \vec p_(i,0))`
`\frac{\partial \vec q}{\partial v}(u, 1) = \sum_(i=0)^3 3B_i(u)(\vec p_(i,3) - \vec p_(i,2))`
`\frac{\partial \vec q}{\partial u}(0, v) = \sum_(j=0)^3 3B_j(v)(\vec p_(1,j) - \vec p_(0,j))`
`\frac{\partial \vec q}{\partial u}(1, v) = \sum_(j=0)^3 3B_j(v)(\vec p_(3,j) - \vec p_(2,j))`
To get corners of patches to join up correctly, we are also going to need to know the second-order partial derivatives. These can be calculated from the equation above as:
`\frac{\partial^2 \vec q}{\partial u \partial v}(0, 0) = 9\cdot (\vec p_(1,1) - \vec p_(0,1) - \vec p_(1, 0) - \vec p_(0,0))`
`\frac{\partial^2 \vec q}{\partial u \partial v}(0, 1) = 9\cdot (\vec p_(1,3) - \vec p_(0,3) - \vec p_(1, 2) - \vec p_(0,2))`
`\frac{\partial^2 \vec q}{\partial u \partial v}(1, 0) = 9\cdot (\vec p_(2,1) - \vec p_(2,1) - \vec p_(3, 0) - \vec p_(2,0))`
`\frac{\partial^2 \vec q}{\partial u \partial v}(1, 1) = 9\cdot (\vec p_(3,3) - \vec p_(2,3) - \vec p_(3, 2) - \vec p_(2,2))`

Joining Bezier Patches

We now want to consider gluing multiple such patches together to specify a smooth surface.
Suppose we have two Bezier Patches `\vec q_1(u,v)` and `\vec q_2(u,v)` with control points `\vec p_(i,j)` and `\vec r_(i,j)` respectively.
We will let `\vec q_2(u,v)` have the domain `[0,1]\times[0,1]` and translate `\vec q_1(u,v)`'s domain so that it is `[-1,0]\times[0,1]`.
To make the right boundary of `\vec q_1` match with the left boundary of `\vec q_2` we need to ensure that the control points along these boundaries match up. This requires: `\vec p_(3,j) = \vec r_(0,j)` for `j = 0, 1, 2, 3`.
Suppose we now have continuous boundaries and want the partial derivatives at the boundaries to be continuous.
So we are interested in continuity in the `u` direction which will be across the boundary as the `v` direction is handle by the two curves lining up on the boundary.
It turns out for `C^1`-continuity that it is sufficient that `\vec p_(3,j) - \vec p_(2,j) = \vec r_(1,j) - \vec r_(0,j)` where `j = 0, 1, 2, 3`. For `G^1`-continuity you just need these differences are proportional.

Subdividing Bezier Patches

two Bezier patches, one subdivided, one not resulting in a crack

As with Bezier curves you can use recursive subdivision to draw Bezier patches
The idea is to recursively split a given patch into smaller patches until each patch is essentially a flat quadrilateral and then draw this using two triangles.
There is now one new difficulty: since some patches may need to be subdivided further than others, it can happen that a surface is subdivided and its neighbor is not.
This can result in a problem called "cracking" illustrated by the diagram above.
One fix is to replace the common boundary of the two patches with a straight line. However, this can cause problems with shading because the surface normals along the common boundary should also line up.
Another way to fix this is to subdivide `q_2` further so that the common boundary has the same points and normals with respect to both patches.

OpenGL and Bezier Curves and Patches

The built-in Bezier functions in OpenGL can only be used for drawing not animating.
In the 1-D case the basic idea to use these functions is to create your own array with the control points of your curve.
Then you use the glMap1*(GL_MAP1_VERTEX_3, ...) function to tell OpenGL that these points will be used to draw a Bezier curve and it also tells OpenGL about how it is parameterized.
Next to enable OpenGL to use the set-up curve you have a command like: glEnable(GL_MAP1_VERTEX_3);
Finally, whenever you want to get the coordinates of `u` along the curve, you call the function: glEvalCoord1*(u);
As usual, by '*' in the above I mean either f, or i, or iv, fv (for arrays).
The process of drawing a Bezier curve is illustrated by the code on the next slide.

OpenGL Bezier Curve Example

/*-----------------------------------------------*/
void drawBezier()
/*
 PURPOSE: Drawing a Bezier curve with glMap1* functions 
 RECEIVES: nothing
 RETURNS: nothing
 REMARKS: 
 */
{
    GLfloat ctrlPts[4][3] = {{-40.0, -40.0, 0.0}, 
        {-10.0, 200.0, 0.0}, {10.0, -200.0, 0.0}, {40.0, 40.0, 0.0}};
    glMap1f (GL_MAP1_VERTEX_3, 0.0 /* uMin*/, 1.0 /* uMax*/, 
        3 /* stride*/, 4 /*num points */,
	          *ctrlPts);
    glEnable(GL_MAP1_VERTEX_3);
    GLint k;	
    // draw curve in blue
    glColor3f(0.0, 0.0, 1.0);	
    glBegin(GL_LINE_STRIP);
    for (k = 0; k <= 50; k++) {
        glEvalCoord1f(GLfloat(k)/50.0);
    }
    glEnd();
	
    // draw control points in red
    glColor3f(1.0, 0.0, 0.0);
    glPointSize(5.0);	
    glBegin(GL_POINTS);
    for (k = 0; k < 4; k++) {
        glVertex3fv(&ctrlPts[k][0]);
    }
    glEnd();
	
    glFlush();
}

Remarks

OpenGL provides a functions glMapGrid1*(n, u1, u2) and glEvalMesh1(mode, n1, n2), to generate uniform spaced points along the curve.
In particular, we could have replaced the for loop of the last example with:
```
glMapGrid1f(50, 0.0, 1.0);
glEvalMesh1(GL_LINE, 0, 50);
```
For Bezier patches, we use the functions glMap2*(GL_MAP_VERTEX_3, uMin, uMax, uStride, nuPts, vMin, vMax, vStride, nvPts, *ctrlPts); to tell OpenGL about a Bezier patch.
Then we call glEnable(GL_MAP2_VERTEX_3); to get it render to calculate points on it. We call glDisable(GL_MAP2_VERTEX_3); later if we want to shut this off.
To evaluate points on the surface we can call glEvalCoord2*(uValue, vValue);
Alternatively, we could also used: glMapGrid2*(nu, u1, u2, nv, v1, v2); and glEvalMesh2(mode, nu1, nu2, nv1, nv2);

Drawing a Bezier Patch

/*-----------------------------------------------*/
void drawBezierPatch()
/*
 PURPOSE: Drawing a Bezier Patch with glMap2* functions 
 RECEIVES: nothing
 RETURNS: nothing
 REMARKS:   
 */
{
    GLfloat ctrlPts[4][4][3] = {
        {{-1.5, -1.5, 4.0}, {-.5, -1.5, 2.0}, 
         {-.5, -1.5, -1.0}, {1.5, -1.5, 2.0} }, 
        {{-1.5, -.5, 1.0}, {-.5, -.5, 3.0},
         {.5, -.5, 0.0}, {1.5, -.5, 1.0} }, 
        {{-1.5, .5, 4.0}, {-.5, .5, 0.0}, 
         {.5, .5, 3.0}, {1.5, .5, 4.0} },
        {{-1.5, 1.5, -2.0}, {-.5, 1.5, -2.0}, 
         {.5, 1.5, 0.0}, {1.5, 1.5, -1.0} }
    };
	
    glMap2f (GL_MAP2_VERTEX_3, 0.0 /* uMin*/, 1.0 /* uMax*/, 
        3 /* u stride*/, 4 /*num points u*/,
        0.0 /* vMin*/, 1.0 /* vMax*/, 12 /*v stride*/, 
        4, &ctrlPts[0][0][0]);

    glEnable(GL_MAP2_VERTEX_3);
    GLint k, j;

    //draw patch in red
    glColor3f(1.0, 0.0, 0.0);
    glMapGrid2f(50, 0.0, 1.0, 50, 0.0, 1.0);
    glEvalMesh2(GL_FILL, 0, 50, 0, 50);

    // draw curves on patch in blue
    glColor3f(0.0, 0.0, 1.0);
	
    for (k = 0 ; k <= 8; k++) {
        glBegin(GL_LINE_STRIP);
        for(j = 0; j <= 40; j++) {
            glEvalCoord2f(GLfloat(j)/40.0, GLfloat(k)/8.0);
        }
        glEnd();
        glBegin(GL_LINE_STRIP);
        for (j = 0; j <= 40; j++) {
            glEvalCoord2f(GLfloat(k)/8.0, GLfloat(j)/40.0);
        }
        glEnd();	
    }
    glFlush();
}

More Splines, Recursive Subdivision, Patches, OpenGL

Outline

Introduction

Bessel-Overhauser Splines

More Bessel-Overhauser Splines

Recursive Subdivision

Applications of Recursive Subdivision

Another Application of Recursive Subdivision

Bezier Surface Patches

More on Bezier Patches

Joining Bezier Patches

Subdividing Bezier Patches

OpenGL and Bezier Curves and Patches

OpenGL Bezier Curve Example

Example Screenshot

Remarks

Drawing a Bezier Patch

Example Screenshot