On-Line Geometric Modeling Notes
BICUBIC UNIFORM B-SPLINE
SURFACE REFINEMENT

Overview

Subdivision surfaces are based upon the binary subdivision of the uniform B-spline surface. In general, they are defined by a initial polygonal mesh, along with a subdivision (or refinement) operation which, given a polygonal mesh, will generate a new mesh that has a greater number of polygonal elements, and is hopefully ``closer'' to some resulting surface. By repetitively applying the subdivision procedure to the initial mesh, we generate a sequence of meshes that (hopefully) converges to a resulting surface. As it turns out, this is a well known process when the mesh has a ``rectangular'' structure and the subdivision procedure is an extension of binary subdivision for uniform B-spline surfaces.

Ed Catmull and Jim Clark, while still graduate students at the University of Utah, sought to extend Doo and Sabin's algorithm to bicubic surfaces. Since the methods of Doo-Sabin is based upon the binary subdivision of the uniform bi-quadratic B-spline patches, Catmull and Clark believed that study of the cubic case would lead to a better subdivision surface generation scheme.

For a pdf version of these notes look here.

A Matrix Equation for the Bicubic Uniform Spline Surfaces

Consider the bicubic uniform B-spline surface ${\bf P} (u,v)$ defined by the $4 \times 4$ array of control points

$$P \: = \: \left[ \begin{array}{cccc} {\bf P} _{0,0} & {\bf P} _{0... ... _{3,0} & {\bf P} _{3,1} & {\bf P} _{3,2} & {\bf P} _{3,3} \end{array}\right]$$

where

$${\bf P} (u,v) \: = \: \left[ \begin{array}{cccc} 1 & u & u^2 & u... ...ht] M P M^T \left[ \begin{array}{c} 1 \\ v \\ v^2 \\ v^3 \end{array}\right]$$

where $M$ is the $4 \times 4$ matrix

$$M \: = \: \frac{1}{6} \left[ \begin{array}{rrrr} 1 & 4 & 1 & 0 \\ -3 & 0 & 3 & 0 \\ 3 & -6 & 3 & 0 \\ -1 & 3 & -3 & 1 \end{array}\right]$$

The control point mesh for such a patch is shown in the figure below.

Subdividing the Bicubic Patch

The bicubic uniform B-spline patch can be subdivided into four subpatches, which can be generated from 25 unique sub-control points. We focus on the subdivision scheme for only one of the four (the subpatch corresponding to $0 \leq u,v \leq \frac{1}{2}$), as the others will follow by symmetry. The following figure illustrates the 25 points produced by subdividing into four subpatches. We have outlined the initial subpatch that we consider below. It should be noted that the nine ``interior'' control points are utilized by each of the four subpatch components of the subdivision.

This subpatch can be generated by reparameterizing the surface by the variables $u'$ and $v'$, where $u'=\frac{u}{2}$ and $v'=\frac{v}{2}$. Substituting these into the equation, we obtain

$${\bf P} (u',v') = {\bf P} ( \frac{u}{2}, \frac{v}{2} )$$

$$= \left[ \begin{array}{cccc} 1 & \frac{u}{2} & \left(\frac{u}{2... ...}{2}\right)^2 \\ \left(\frac{v}{2}\right)^3 \end{array}\right]$$

$$= \left[ \begin{array}{cccc} 1 & u & u^2 & u^3 \end{array}\righ... ...\left[ \begin{array}{c} 1 \\ v \\ v^2 \\ v^3 \end{array}\right]$$

$$= \left[ \begin{array}{cccc} 1 & u & u^2 & u^3 \end{array}\righ... ...\left[ \begin{array}{c} 1 \\ v \\ v^2 \\ v^3 \end{array}\right]$$

$$= \left[ \begin{array}{cccc} 1 & u & u^2 & u^3 \end{array}\righ... ...\left[ \begin{array}{c} 1 \\ v \\ v^2 \\ v^3 \end{array}\right]$$

$$= \left[ \begin{array}{cccc} 1 & u & u^2 & u^3 \end{array}\righ... ...\left[ \begin{array}{c} 1 \\ v \\ v^2 \\ v^3 \end{array}\right]$$

where $P' = S P S^T$ and

$$S \: = \: M^{-1} \left[ \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 ... ... 0 & 0 & \frac{1}{4} & 0 \\ 0 & 0 & 0 & \frac{1}{8} \end{array}\right] M$$

Through this process, we have written the surface ${\bf P} '(u,v)$ as

$${\bf P} '(u,v) \: = \: \left[ \begin{array}{cccc} 1 & u & u^2 & ... ...] M P' M^T \left[ \begin{array}{c} 1 \\ v \\ v^2 \\ v^3 \end{array}\right]$$

for some $4 \times 4$ control point array $P'$. This implies that ${\bf P} '(u,v)$ is a uniform bicubic B-spline patch. The matrix $S$ is typically called the ``splitting matrix'', and is straightforward to calculate. It is given by

$$H_1 = \frac{1}{8} \left[ \begin{array}{cccc} 4 & 4 & 0 & 0 \\ 1 & 6 & 1 & 0\\ 0 & 4 & 4 & 0 \\ 0 & 1 & 6 & 1 \end{array}\right]$$

By carrying out the algebra, we can calculate the control point array $P'$ by

$$P' = \frac{1}{8} \left[ \begin{array}{cccc} 4 & 4 & 0 & 0 \\ 1 & ... ... 1 & 0\\ 0 & 4 & 4 & 0 \\ 0 & 1 & 6 & 1 \end{array}\right] ^T$$

$$= \frac{1}{8} \left[ \begin{array}{cccc} 4 {\bf P} _{0,0} + 4 {... ...6 & 4 & 1\\ 0 & 1 & 4 & 6 \\ 0 & 0 & 6 & 1 \end{array}\right]$$

and we obtain

$${\bf P} _{0,0}' = \frac{ {\bf P} _{0,0}+ {\bf P} _{1,0}+ {\bf P} _{0,1}+ {\bf P} _{1,1}}{4}$$

$${\bf P} _{0,1}' = \frac{ {\bf P} _{0,0}+ {\bf P} _{1,0}+6( {\bf P} _{0,1}+ {\bf P} _{1,1})+ {\bf P} _{0,2}+ {\bf P} _{1,2}}{16}$$

$${\bf P} _{0,2}' = \frac{ {\bf P} _{0,1}+ {\bf P} _{1,1}+ {\bf P} _{0,2}+ {\bf P} _{1,2}}{4}$$

$${\bf P} _{0,3}' = \frac{ {\bf P} _{0,1}+ {\bf P} _{1,1}+6( {\bf P} _{0,2}+ {\bf P} _{1,2})+ {\bf P} _{0,3}+ {\bf P} _{1,3}}{16}$$

$${\bf P} _{1,0}' = \frac{ {\bf P} _{0,0}+ {\bf P} _{0,1}+6( {\bf P} _{1,0}+ {\bf P} _{1,1})+ {\bf P} _{2,0}+ {\bf P} _{2,1}}{16}$$

$${\bf P} _{1,1}' = \frac{ {\bf P} _{0,0}+6 {\bf P} _{1,0}+ {\bf P} _{2,0}+6( {\bf P}... ...} _{1,1}+ {\bf P} _{2,1})+ {\bf P} _{0,2}+6 {\bf P} _{1,2}+ {\bf P} _{2,2}}{64}$$

$${\bf P} _{1,2}' = \frac{ {\bf P} _{0,1}+ {\bf P} _{0,2}+6( {\bf P} _{1,1}+ {\bf P} _{1,2})+ {\bf P} _{2,1}+ {\bf P} _{2,2}}{16}$$

$${\bf P} _{1,3}' = \frac{ {\bf P} _{0,1}+6 {\bf P} _{1,1}+ {\bf P} _{2,1}+6( {\bf P}... ...} _{1,2}+ {\bf P} _{2,2})+ {\bf P} _{0,3}+6 {\bf P} _{1,3}+ {\bf P} _{2,3}}{64}$$

$${\bf P} _{2,0}' = \frac{ {\bf P} _{1,0}+ {\bf P} _{2,0}+ {\bf P} _{1,1}+ {\bf P} _{2,1}}{4}$$

$${\bf P} _{2,1}' = \frac{ {\bf P} _{1,0}+ {\bf P} _{2,0}+6( {\bf P} _{1,1}+ {\bf P} _{2,1})+ {\bf P} _{1,2}+ {\bf P} _{2,2}}{16}$$

$${\bf P} _{2,2}' = \frac{ {\bf P} _{1,1}+ {\bf P} _{2,1}+ {\bf P} _{1,2}+ {\bf P} _{2,2}}{4}$$

$${\bf P} _{2,3}' = \frac{ {\bf P} _{1,1}+ {\bf P} _{2,1}+6( {\bf P} _{1,2}+ {\bf P} _{2,2})+ {\bf P} _{1,3}+ {\bf P} _{2,3}}{16}$$

$${\bf P} _{3,0}' = \frac{ {\bf P} _{1,0}+ {\bf P} _{1,1}+6( {\bf P} _{2,0}+ {\bf P} _{2,1})+ {\bf P} _{3,0}+ {\bf P} _{3,1}}{16}$$

$${\bf P} _{3,1}' = \frac{ {\bf P} _{1,0}+6 {\bf P} _{2,0}+ {\bf P} _{3,0}+6( {\bf P}... ...} _{2,1}+ {\bf P} _{3,1})+ {\bf P} _{1,2}+6 {\bf P} _{2,2}+ {\bf P} _{3,2}}{64}$$

$${\bf P} _{3,2}' = \frac{ {\bf P} _{1,1}+ {\bf P} _{1,2}+6( {\bf P} _{2,1}+ {\bf P} _{2,2})+ {\bf P} _{3,1}+ {\bf P} _{3,2}}{16}$$

$${\bf P} _{3,3}' = \frac{ {\bf P} _{1,1}+6 {\bf P} _{2,1}+ {\bf P} _{3,1}+6( {\bf P}... ...} _{2,2}+ {\bf P} _{3,2})+ {\bf P} _{1,3}+6 {\bf P} _{2,3}+ {\bf P} _{3,3}}{64}$$

Each of these points can be classified into three categories - face points, edge points and vertex points - depending on each points relationship to the original control point mesh. The points ${\bf P} _{0,0}'$, ${\bf P} _{0,2}'$, ${\bf P} _{2,0}'$ and ${\bf P} _{2,2}'$, which are shown in the figure below

are called ``face'' points, and are calculated as the average of the four points that bound the respective face. If we define the face point $F_{i,j}$ to be the average of the points $P_{i,j}$, $P_{i+1,j}$, $P_{i,j+1}$ and $P_{i+1,j+1}$, then we can rewrite the above equations with these face points substituted on the right-hand side, and obtain

$${\bf P} _{0,0}' = {\bf F} _{0,0}$$

$${\bf P} _{0,1}' = \frac{4 {\bf F} _{0,0}+4 {\bf F} _{0,1}+4 {\bf P} _{0,1}+4 {\bf P} _{1,1}}{16}$$

$${\bf P} _{0,2}' = {\bf F} _{0,1}$$

$${\bf P} _{0,3}' = \frac{4 {\bf F} _{0,1}+4 {\bf F} _{0,2}+4 {\bf P} _{0,2}+4 {\bf P} _{1,2}}{16}$$

$${\bf P} _{1,0}' = \frac{4 {\bf F} _{0,0}+4 {\bf F} _{1,0}+4 {\bf P} _{1,0}+4 {\bf P} _{1,1}}{16}$$

$${\bf P} _{1,1}' = \frac{4 {\bf F} _{0,0}+4 {\bf F} _{0,1}+4 {\bf F} _{1,0}+4 {\bf F... ...{1,0}+4 {\bf P} _{0,1}+32 {\bf P} _{1,1}+4 {\bf P} _{2,1}+4 {\bf P} _{1,2}}{64}$$

$${\bf P} _{1,2}' = \frac{4 {\bf F} _{0,1}+4 {\bf F} _{1,1}+4 {\bf P} _{1,1}+4 {\bf P} _{1,2}}{16}$$

$${\bf P} _{1,3}' = \frac{4 {\bf F} _{0,1}+4 {\bf F} _{0,2}+4 {\bf F} _{1,1}+4 {\bf F... ...{1,1}+4 {\bf P} _{0,2}+32 {\bf P} _{1,2}+4 {\bf P} _{2,2}+4 {\bf P} _{1,3}}{64}$$

$${\bf P} _{2,0}' = {\bf F} _{1,0}$$

$${\bf P} _{2,1}' = \frac{4 {\bf F} _{1,0}+4 {\bf F} _{1,1}+4 {\bf P} _{1,1}+4 {\bf P} _{2,1}}{16}$$

$${\bf P} _{2,2}' = {\bf F} _{1,1}$$

$${\bf P} _{2,3}' = \frac{4 {\bf F} _{1,1}+4 {\bf F} _{1,2}+4 {\bf P} _{1,2}+4 {\bf P} _{2,2}}{16}$$

$${\bf P} _{3,0}' = \frac{4 {\bf F} _{1,0}+4 {\bf F} _{2,0}+4 {\bf P} _{2,0}+4 {\bf P} _{2,1}}{16}$$

$${\bf P} _{3,1}' = \frac{4 {\bf F} _{1,0}+4 {\bf F} _{2,0}+4 {\bf F} _{1,1}+4 {\bf F... ...{2,0}+4 {\bf P} _{1,1}+32 {\bf P} _{2,1}+4 {\bf P} _{3,1}+4 {\bf P} _{2,2}}{64}$$

$${\bf P} _{3,2}' = \frac{4 {\bf F} _{1,1}+4 {\bf F} _{2,1}+4 {\bf P} _{2,1}+4 {\bf P} _{2,2}}{16}$$

$${\bf P} _{3,3}' = \frac{4 {\bf F} _{1,1}+4 {\bf F} _{2,1}+4 {\bf F} _{1,2}+4 {\bf F... ...{2,1}+4 {\bf P} _{1,2}+32 {\bf P} _{2,2}+4 {\bf P} _{3,2}+4 {\bf P} _{2,3}}{64}$$

Simplifying these equations, we obtain

$${\bf P} _{0,0}' = {\bf F} _{0,0}$$

$${\bf P} _{0,1}' = \frac{ {\bf F} _{0,0}+ {\bf F} _{0,1}+ {\bf P} _{0,1}+ {\bf P} _{1,1}}{4}$$

$${\bf P} _{0,2}' = {\bf F} _{0,1}$$

$${\bf P} _{0,3}' = \frac{ {\bf F} _{0,1}+ {\bf F} _{0,2}+ {\bf P} _{0,2}+ {\bf P} _{1,2}}{4}$$

$${\bf P} _{1,0}' = \frac{ {\bf F} _{0,0}+ {\bf F} _{1,0}+ {\bf P} _{1,0}+ {\bf P} _{1,1}}{4}$$

$${\bf P} _{1,1}' = \frac{ {\bf F} _{0,0}+ {\bf F} _{0,1}+ {\bf F} _{1,0}+ {\bf F} _{... ...P} _{1,0}+ {\bf P} _{0,1}+8 {\bf P} _{1,1}+ {\bf P} _{2,1}+ {\bf P} _{1,2}}{16}$$

$${\bf P} _{1,2}' = \frac{ {\bf F} _{0,1}+ {\bf F} _{1,1}+ {\bf P} _{1,1}+ {\bf P} _{1,2}}{4}$$

$${\bf P} _{1,3}' = \frac{ {\bf F} _{0,1}+ {\bf F} _{0,2}+ {\bf F} _{1,1}+ {\bf F} _{... ...P} _{1,1}+ {\bf P} _{0,2}+8 {\bf P} _{1,2}+ {\bf P} _{2,2}+ {\bf P} _{1,3}}{16}$$

$${\bf P} _{2,0}' = {\bf F} _{1,0}$$

$${\bf P} _{2,1}' = \frac{ {\bf F} _{1,0}+ {\bf F} _{1,1}+ {\bf P} _{1,1}+ {\bf P} _{2,1}}{4}$$

$${\bf P} _{2,2}' = {\bf F} _{1,1}$$

$${\bf P} _{2,3}' = \frac{ {\bf F} _{1,1}+ {\bf F} _{1,2}+ {\bf P} _{1,2}+ {\bf P} _{2,2}}{4}$$

$${\bf P} _{3,0}' = \frac{ {\bf F} _{1,0}+ {\bf F} _{2,0}+ {\bf P} _{2,0}+ {\bf P} _{2,1}}{4}$$

$${\bf P} _{3,1}' = \frac{ {\bf F} _{1,0}+ {\bf F} _{2,0}+ {\bf F} _{1,1}+ {\bf F} _{... ...P} _{2,0}+ {\bf P} _{1,1}+8 {\bf P} _{2,1}+ {\bf P} _{3,1}+ {\bf P} _{2,2}}{16}$$

$${\bf P} _{3,2}' = \frac{ {\bf F} _{1,1}+ {\bf F} _{2,1}+ {\bf P} _{2,1}+ {\bf P} _{2,2}}{4}$$

$${\bf P} _{3,3}' = \frac{ {\bf F} _{1,1}+ {\bf F} _{2,1}+ {\bf F} _{1,2}+ {\bf F} _{... ...P} _{2,1}+ {\bf P} _{1,2}+8 {\bf P} _{2,2}+ {\bf P} _{3,2}+ {\bf P} _{2,3}}{16}$$

In examining these equations, we see that the points ${\bf P} _{0,1}'$, ${\bf P} _{0,3}'$, ${\bf P} _{1,0}'$, ${\bf P} _{1,2}'$, ${\bf P} _{2,1}'$, ${\bf P} _{2,3}'$, ${\bf P} _{3,0}'$ and ${\bf P} _{3,2}'$, which are called ``edge'' points, are given as the average of four points - the two points that define the original edge and the two two new face points of the faces sharing the edge. This is shown in the following figure.

These edge points ${\bf E} _{i,j}$ can be calculated either as

$${\bf E} _{i,j} \: = \: \frac{ {\bf F} _{i,j-1} + {\bf F} _{i,j} + {\bf P} _{i,j} + {\bf P} _{i+1,j}}{4}$$

or

$${\bf E} _{i,j} \: = \: \frac{ {\bf F} _{i-1,j} + {\bf F} _{i,j} + {\bf P} _{i,j} + {\bf P} _{i,j+1}}{4}$$

depending on which side of the edge point the two faces lie. If we replace these edge points on the right-hand side of the equations above, we obtain

$${\bf P} _{0,0}' = {\bf F} _{0,0}$$

$${\bf P} _{0,1}' = {\bf E} _{0,1}$$

$${\bf P} _{0,2}' = {\bf F} _{0,1}$$

$${\bf P} _{0,3}' = {\bf E} _{0,2}$$

$${\bf P} _{1,0}' = {\bf E} _{1,0}$$

$${\bf P} _{1,1}' = \frac{ {\bf F} _{0,0}+ {\bf F} _{0,1}+ {\bf F} _{1,0}+ {\bf F} _{... ...P} _{1,0}+ {\bf P} _{0,1}+8 {\bf P} _{1,1}+ {\bf P} _{2,1}+ {\bf P} _{1,2}}{16}$$

$${\bf P} _{1,2}' = {\bf E} _{1,2}$$

$${\bf P} _{1,3}' = \frac{ {\bf F} _{0,1}+ {\bf F} _{0,2}+ {\bf F} _{1,1}+ {\bf F} _{... ...P} _{1,1}+ {\bf P} _{0,2}+8 {\bf P} _{1,2}+ {\bf P} _{2,2}+ {\bf P} _{1,3}}{16}$$

$${\bf P} _{2,0}' = {\bf F} _{1,0}$$

$${\bf P} _{2,1}' = {\bf E} _{2,1}$$

$${\bf P} _{2,2}' = {\bf F} _{1,1}$$

$${\bf P} _{2,3}' = {\bf E} _{2,2}$$

$${\bf P} _{3,0}' = {\bf E} _{3,0}$$

$${\bf P} _{3,1}' = \frac{ {\bf F} _{1,0}+ {\bf F} _{2,0}+ {\bf F} _{1,1}+ {\bf F} _{... ...P} _{2,0}+ {\bf P} _{1,1}+8 {\bf P} _{2,1}+ {\bf P} _{3,1}+ {\bf P} _{2,2}}{16}$$

$${\bf P} _{3,2}' = {\bf E} _{3,2}$$

$${\bf P} _{3,3}' = \frac{ {\bf F} _{1,1}+ {\bf F} _{2,1}+ {\bf F} _{1,2}+ {\bf F} _{... ...P} _{2,1}+ {\bf P} _{1,2}+8 {\bf P} _{2,2}+ {\bf P} _{3,2}+ {\bf P} _{2,3}}{16}$$

The remaining four points, ${\bf P} _{1,1}'$, ${\bf P} _{1,3}'$, ${\bf P} _{3,1}'$ and ${\bf P} _{3,3}'$, as shown in the figure below

are called ``vertex points''. These points, as can be seen above, are somewhat complex, but after some reduction it can be seen that

$${\bf P} _{i,j}' \: = \: \frac{ {\bf Q} + 2 {\bf R} + {\bf S} }{4}$$

where ${\bf Q}$ is the average of the face points of the faces adjacent to the vertex point, ${\bf R}$ is the average of the midpoints of the edges adjacent to the vertex point and ${\bf S}$ is the corresponding vertex from the original mesh. For example, if we consider the point ${\bf P} _{1,3}'$, then

$${\bf Q} = \frac{ {\bf F} _{0,1}+ {\bf F} _{0,2}+ {\bf F} _{1,1}+ {\bf F} _{1,2}}{4}$$

$${\bf R} = \frac{ \frac{ {\bf P} _{0,2}+ {\bf P} _{1,2}}{2} + \frac{ {\bf P}... ...f P} _{1,3}+ {\bf P} _{1,2}}{2} + \frac{ {\bf P} _{2,2}+ {\bf P} _{1,2}}{2}}{4}$$

$${\bf S} = {\bf P} _{1,2}$$

All sixteen points of the subdivision have now been characterized in terms of face points, edge points and vertex points, and a geometric method has been developed to calculate these points.

Extending this Subdivision Procedure to the Entire Patch

We note that all 25 of the points can actually be calculated in this manner, as for example ${\bf P} _{4,4}'$ is a face point and can be calculated as the average of the four points bounding the face. In general, we call the mesh generated by the 25 points as a refinement of the original mesh. In this case, we can state the following rules to generate the points for the refinement of the surface:

• For each face in the original mesh, generate the new face points - which are the average of all the original points defining the face.
• For each internal edge of the original mesh (i.e. an edge not on the boundary), generate the new edge points - which are calculated as the average of four points: the two points which define the edge, and the two new face points for the faces that are adjacent to the edge.
• For each internal vertex of the original mesh (i.e. a vertex not on the boundary of the mesh), generate the new vertex points - which are calculated as the average of ${\bf Q}$, $2 {\bf R}$ and ${\bf S}$, where ${\bf Q}$ is the average of the new face points of all faces adjacent to the original vertex point, ${\bf R}$ is the average of the midpoints of all original edges incident on the original vertex point, and ${\bf S}$ is the original vertex point.

The process generated from these rules actually extends to arbitrary rectangular meshes, so we can perform this process again on our refined mesh of 25 elements, producing a second refinement of the original mesh. In this case, we know that this represents yet another subdivision and that eventually, if we keep refining, this ``limit mesh'' will converge to the original uniform B-spline surface.

Thus, this process gives us a sequence of meshes, each of which is a refinement of the mesh directly above, and which converges to the surface in the limit. For the purposes of rendering such a surface we can simply let the refinement process go until we have a mesh that is ``sufficiently close'' to the actual surface and then utilize the mesh for rendering purposes.

Summary

The subdivision of the bicubic uniform B-spline surface produces a simple procedure based upon face points, edge points and vertex points, and can be extended to be a refinement procedure for an extended mesh based upon a rectangular topology. Catmull and Clark were able to take this procedure and produce a refinement strategy that works on a mesh of arbitrary topology.

Bibliography

1

CATMULL, E., AND CLARK, J.
Recursively generated B-spline surfaces on arbitrary topological meshes.
Computer-Aided Design 10 (Sept. 1978), 350-355.