About
At last year’s SGP conference,
researchers presented an interesting method
for meshing NURBS surfaces. This piqued my interest, since I had written a
NURBS mesher earlier that year (as part of
Foxtrot, a STEP file viewer).
On a long train ride, I spent some time trying
to understand the paper. I quickly realized that I’d need to go upstream: the
NURBS meshing paper was the the latest in a series of related papers,
and heavily relied on previous works for geometric infrastructure.
The lineage is roughly
- Curve/Surface Intersection Problem by Means of Matrix Representations
- Implicit matrix representations of rational Bézier curves and surfaces
- A Line/Trimmed NURBS Surface Intersection Algorithm Using Matrix Representations
- Delaunay Meshing and Repairing of NURBS Models
This year, as part of my traditional Winter Break Graphics Hackathon, I dug into
the first three papers. At the end, I could successfully raytrace the
Utah Teapot
directly from its Bézier patches:
I gave running commentary in a Twitter thread;
this writeup goes into more depth, and should be helpful for others seeking to
understand this line of research.
The implicit matrix representation
Let’s start in the middle.
These papers describe a way of converting algebraic curves and surfaces
into a function which converts a point $(x, y, z)$ into a matrix
$M(x, y, z)$.
The matrix $M(x, y, z)$ is a linear combination of constant matrixes,
e.g.
$$
M(x, y, z) = M_0 + x M_X + y M_Y + z M_z
$$
This matrix has an interesting property: its rank drops by one
when you’re exactly on the curve or surface!
Consider the cubic Bézier curve with sample points $left[[0,0], [1,2], [2,1], [3,3]right]$:
Through a mysterious process (to be explained later), we can calculate the
$M$ matrices. They’re not particularly pretty:
$$
M_0 = left[begin{array}{cccccc}
0 & 0 & 0 \
-0.47130973& -0.69441025& 0.27461067 \
0.50379062& -0.58034736& -0.5515659 \
end{array}right]
$$
$$
M_X = left[begin{array}{cccccc}
-0.25848205& 0.25319592& -0.65863269\
-0.07003399& 0.27575858& -0.11849897 \
-0.45427448& 0.088577 & -0.05392415\
end{array}right]
$$
$$
M_Y = left[begin{array}{cccccc}
0.28634427& 0.10487212& 0.23777945\
0.28634427& 0.10487212& 0.23777945 \
0.28634427& 0.10487212& 0.23777945\
end{array}right]
$$
(since this is 2D, we’re ignoring $M_Z$)
With these matrices, we can evaluate $M(3, 3)$, which should be on the curve,
then look at its singular values (using the SVD):
$$
sigma(3, 3) = [1.6725, 0.76772, 2.77times10^{-16}]
$$
This seems to work: the last singular value is very close to zero. To confirm,
we’ll check a point that’s off the curve:
$$
sigma(3, 4) = [1.6067, 1.183, 0.153]
$$
In this case, the singular values are all non-zero, meaning the matrix has full
rank. This is what we expect for a point that’s not on the curve!
To go from singular values to a single pseudo-distance, we multiply
them together.
$$
begin{aligned}
prod_i sigma_i(3, 3) &= 3.57 times 10^{-16} \
prod_i sigma_i(3, 4) &= 0.29
end{aligned}
$$
As expected, this pseudo-distance goes to 0 on the curve.
Here’s what pseudo-distance looks like around the curve (shown in white):
That the values go to zero is very obvious when plotted in log scale:
Building the matrix representation
Now, let’s take a step back: where does $M(x,y,z)$ come from?
Remember that we’re dealing with parametric curves and surfaces.
In the general case, we’re taking some parameter $s$
(where $s$ is one-dimensional for a curve and $s = (s_1, s_2)$ for a
surface), and generating a rational curve or surface
$$
left(
frac{f_1(s)}{f_0(s)},
frac{f_2(s)}{f_0(s)},
frac{f_3(s)}{f_0(s)}
right)
$$
In our example above, we’re generating a 2D Bézier curve, which means that
we have a 1D parameterization with
$f_0(s) = 1$ and $f_3(s) = 0$.
Consider a second set of functions $g_{0,1,2,3}(s)$ with the property
$$
sum_{i=0}^3 g_i(s)f_i(s) = 0
$$
Right now, the $g_i(s)$ functions are underdefined; we don’t know exactly
what they are, just that the above property must hold for a given set of
$f_i(s)$ functions.
Now, let’s use the $g_i(s)$ functions to build a new function, which is
parameterized by $s$ and evaluated at some point in space
$left(x, y, zright)$
$$
L(s; x, y, z) = g_0(s) + xg_1(s) + yg_2(s) + zg_3(s)
$$
Again, $L$ is not a specific function; it’s just a particular way of
combining a set of $g_i(s)$ functions. However, it already has an
interesting property: what happens if we evaluated it at a point on the
curve or surface?
$$
begin{aligned}
Lleft(s; frac{f_1(s)}{f_0(s)}, frac{f_2(s)}{f_0(s)}, frac{f_3(s)}{f_0(s)}right)
&= g_0(s) + frac{f_1(s)}{f_0(s)}g_1(s) + frac{f_2(s)}{f_0(s)}g_2(s) + frac{f_3(s)}{f_0(s)}g_3(s)
\
&= frac{g_0(s)f_0(s) + f_1(s)g_1(s) + f_2(s)g_2(s) + f_3(s)g_3(s)}{f_0(s)}
\
&= frac{0}{f_0(s)}
\
&= 0
end{aligned}
$$
That’s neat: $L(s; x, y, z)$ is 0 at the point on
the surface given by $s$!
The paper gives more intuition for the behavior of $L$ by noting that once
we’ve locked down a particular value for $s$, $L_s(x, y, z)$
defines a plane in 3D. From our observation above, that plane must pass
through the point on the surface given by $s$.
Let’s jump ahead very briefly, so as to break up the wall of math with a picture.
In our Bézier curve example, we end up with three $L(s; x, y)$ functions which
each define a line in 2D for a particular value of $s$.
Here’s what they look like as you sweep the
value of the curve parameter $s$, with the point on the curve shown as
a black dot:
The three lines always intersect at the point on the curve given by
a particular $s$!
Okay, back to math! We’ve discussed a set of four functions $g_{0,1,2,3}(s)$,
given then a particular property, and combined them into a new function
$L(s; x, y, z)$.
Now, we’ll lock things down a little more.
Consider a set of basis functions $left(psi_1(s), psi_2(s), …, psi_{m_nu}(s)right)$,
where $m_nu$ is the number of basis functions of degree $leq nu$,
and $nu$ is an arbitrarly chosen degree.
We can construct $g_i(s)$ as a weighted sum of these basis functions,
e.g.
$$
g_0(s) = sum_{i=1}^{m_nu}lambda_{0,i}psi_i(s)
$$
Because all of the $g_i(s)$ functions are using the same basis functions,
we can represent their combination as
$$
L(s; x, y, z) = sum_{i=1}^{m_nu}left(lambda_{0,i} + xlambda_{1,i} + ylambda_{2,i} + zlambda_{3,i} right) psi_i(s)
$$
We will collapse the weighted sum into its own function
$$
Lambda_i(x, y, z) = lambda_{0,i} + xlambda_{1,i} + ylambda_{2,i} + zlambda_{3,i}
$$
(Did you know that $Lambda$ is a capital $lambda$?
I didn’t until just now; it’s all Greek to me!)
Now, we can rewrite our previous function:
$$
L(s; x, y, z) = sum_{i=1}^{m_nu}Lambda_i(x, y, z) psi_i(s)
$$
For reasons that will soon become clear, we can also write this as
a matrix multiplication (albeit one where it’s a single row multiplied
by a single column):
$$
left[begin{array}{cccc}
psi_1(s) & psi_2(s) & … & psi_{m_nu}(s)
end{array}
right]
left[begin{array}{c}
Lambda_1(x, y, z) \
Lambda_2(x, y, z) \
… \
Lambda_{m_nu}(x, y, z) \
end{array}right]
= left[L(s; x, y, z)right]
$$
Let’s take a step back: we’ve assumed the existance of a quartet of functions
$g_{0,1,2,3}(s)$, given them certain properties, and built this
infrastructure on top of them.
However, we haven’t constrained them enough that there’s only a single set of
$g$ functions for a given shape. You can imagine generating multiple
sets of $g$ functions, e.g.
$$
begin{aligned}
&left(g_{0,1}(s), g_{1,1}(s), g_{2,1}(s), g_{3,1}(s)right) \
&left(g_{0,2}(s), g_{1,2}(s), g_{2,2}(s), g_{3,2}(s)right) \
&… \
&left(g_{0,r_nu}(s), g_{1,r_nu}(s), g_{2,r_nu}(s), g_{3,r_nu}(s)right)
end{aligned}
$$
Here, we’re generating $r_nu$ sets of $g$ functions, $r_nu$ being
another mystery value that we’ll explain later.
Now that we’ve got multiple sets of $g$ functions, we can use
the rest of our infrastructure and build a big matrix of $Lambda$ values.
As shorthand, let’s just call it $M(x, y, z)$:
$$
M(x, y, z) = left[begin{array}{cccc}
Lambda_{1,1}(x,y,z) &Lambda_{1,2}(x,y,z) & … & Lambda_{1,r_nu}(x,y,z) \
Lambda_{2,1}(x,y,z) &Lambda_{2,2}(x,y,z) & … & Lambda_{2,r_nu}(x,y,z) \
… & … & … & … \
Lambda_{m_nu,1}(x,y,z) &Lambda_{m_nu,2}(x,y,z) & … & Lambda_{m_nu,r_nu}(x,y,z) \
end{array}right]
$$
But wait! That’s a confusing choice of variable, because we used $M$ as
our “magic matrix with drop-of-rank property” earlier.
Well, surprise! This is our magic matrix with drop-of-rank property!
Let’s see why. As before, we can multiply by the basis functions for $g$
to generate a family of functions $L(s; x, y, z)$:
$$
begin{aligned}
&left[begin{array}{cccc}
psi_1(s) & … & psi_{m_nu}(s)
end{array}
right]
M(x, y, z) \ =
&left[begin{array}{cccc}
psi_1(s) & … & psi_{m_nu}(s)
end{array}
right]
left[begin{array}{cccc}
Lambda_{1,1}(x,y,z) & … & Lambda_{1,r_nu}(x,y,z) \
… & … & … \
Lambda_{m_nu,1}(x,y,z) & … & Lambda_{m_nu,r_nu}(x,y,z) \
end{array}right] \ = &left[begin{array}{ccc}L_1(s; x, y, z)& … & L_{m_nu}(s; x, y, z)end{array}right]
end{aligned}
$$
Now we can see where the drop-of-rank property comes from. Consider evaluating
$M$ at a point on the surface given by the parameter $s$:
$$
Mleft(frac{f_1(s)}{f_0(s)}, frac{f_2(s)}{f_0(s)}, frac{f_3(s)}{f_0(s)}right)
$$
From all that has come before, we know that
$$
begin{aligned}
&
