Fast Rational Bézier Surface Evaluation Problem

Question

I am currently writing a NURBS ray tracer. What I do is convert the NURBS into rational Bézier patches and then perform the intersection test using Newton's method. To do this fast (the ray tracer should be later implemented using GPGPU programming and should run in real time) a point on the surface and first derivatives at that point need to be found. I have found two algorithms that should be fast.

First one is by S. Mann and T. DeRose ("Computing values and derivatives of Bézier and B-spline tensor products" (1995)), the problem here is that there is something incorrect about the weights, so round shapes look a bit square. See here.

The second one is by T. W. Sederberg ("Point and tangent computation of tensor product rational Bézier surfaces" (1995)), description of the algorithm that is same as in the paper can be found here, p. 175

For $0 \leq s, t \leq 0.5$ the evaluation is same as in the paper by Sederberg:

$$ {\mathbf{P}^{kl}(s,t)\over(1-s)^{m-1}(1-t)^{n-1}}=\hat{\mathbf{P}}^{kl}(u,v) = \sum_{i=k}^{m+k-l} \sum_{j=l}^{n+l-1} \hat{\mathbf{P}}^{kl}_{ij} u^i v^j; \quad k,l = 0, 1 $$

with $u = {s\over1-s}$, $v = {t\over1-t}$ and $\hat{\mathbf{P}}^{kl}_{ij} = {m-1 \choose i-k} {n-1\choose j-l} \mathbf{P}_{ij}$.

I tried to derive the equations for $0.5 \leq s, t \leq 1$, because there are problems when $s$ or $t$ is close to 1 and came up with:

$$ {\mathbf{P}^{kl}(s,t)\over s^{m-1}t^{n-1}}=\hat{\mathbf{P}}^{kl}(u,v) = \sum_{i=k}^{m+k-l} \sum_{j=l}^{n+l-1} \hat{\mathbf{P}}^{kl}_{m+k-i, n+l-j} u^i v^j; \quad k,l = 0, 1 $$

where $u = {1-s\over s}$, $v = {1-t\over t}$ and $\hat{\mathbf{P}}^{kl}_{ij} = {m-1 \choose i-k} {n-1\choose j-l} \mathbf{P}_{ij}$. I apply the post--multiplication by $s^{m-1}t^{n-1}$ just to be safe, but Sederberg mentions in the curve case that it isn't necessary so I don't really know.

Also I guess the derivatives need the correct magnitude for the Newton's method to work, Sederberg mentions that for a curve the derviative is obtained as:

$$ \mathbf{q}(t) - \mathbf{r}(t) \equiv \Pi(\mathbf{Q}(t) - \mathbf{R}(t)) $$

with correct-magnitude-derivative being equal to:

$$ {d\mathbf{p}(t)\over dt} = n {\mathbf{R}_w(t) \mathbf{Q}_w(t)\over((1 - t)\mathbf{Q}_w(t)+t\mathbf{R}_w(t))^2}\left[ \mathbf{r}(t) - \mathbf{q}(t) \right] $$

Now for the surface derivatives. The curve derivative

$$ \Pi(\mathbf{Q}(t) - \mathbf{R}(t)) $$

looks similar to the surface derivative equations

$$ \mathbf{p}_s(s,t) = \Pi(\underbrace{(1-t) \mathbf{P}^{00}(s,t) + t \mathbf{P}^{01}(s,t))}_{\mathbf{Q}_s(s,t)} - \Pi(\underbrace{(1-t) \mathbf{P}^{10}(s,t) + t \mathbf{P}^{11}(s,t)}_{\mathbf{R}_s(s,t)}) $$

$$ \mathbf{p}_t(s,t) = \Pi(\underbrace{(1-s) \mathbf{P}^{00}(s,t) + s \mathbf{P}^{10}(s,t))}_{\mathbf{Q}_t(s,t)} - \Pi(\underbrace{(1-s) \mathbf{P}^{01}(s,t) + s \mathbf{P}^{11}(s,t)}_{\mathbf{R}_t(s,t)}) $$

So I plugged $\mathbf{Q}_s(s,t)$, $\mathbf{R}_s(s,t)$ and $\mathbf{Q}_t(s,t)$, $\mathbf{R}_t(s,t)$ into the equation for curve derivative but it doesn't seem to work as the Newton's method doesn't converge at all. What am I missing?

Answer 1

So after some time I got both of the algorithms working. In scenes with surfaces of lower degree, algorithm by Mann and DeRose seems to be faster, in other scenes (about 80% of scenes I test my ray tracer with) the algorithm by Sederberg is faster and requires a constant amount of memory so it is more suitable for GPU (no dynamic allocation in OpenCL 1.2). Actually, both of the algorithms have small errors in their versions on Elsevier, so I will write the corrected versions together with the correct magnitude first derivatives and working on the whole parametric domain (Sederberg's algorithm needs to be edited to work when params are close to 1).

First the algorithm by Mann and DeRose, note that from now on $m$ is degree in $U$ direction and $n$ is degree in $V$ direction. Capital letters mean that the point is in the homogeneous space and needs to be projected back into euclidean space by operator $\Pi$, $P = (wx, wy, wz, w)$, $\Pi(P) = p = (P_x/w, P_y/w, P_z/w)$. The derivative in $U$ direction is $p_u$ and the derivative in $V$ direction is $p_v$. Also i drop the $(u,v)$, shouldn't be too confusing. The algorithm as it is presented overwrites the points, so you should work with their copy if you would like to keep the data.

$$ \sum_{i=0}^m \sum_{k=1}^n \sum_{j=0}^{n-k}P_{i,j}=(1-v) \cdot P_{i,j}+v \cdot P_{i,j+1} $$ $$ \sum_{j=0}^1 \sum_{k=1}^m \sum_{i=0}^{m-k}P_{i,j}=(1-u) \cdot P_{i,j}+u \cdot P_{i+1,j} $$ (think of the sums as if they were for loops, outer sum is outer loop etc.) $$ L_{u,0} = (1-u)\cdot P_{0,0} + u \cdot P_{1,0} $$ $$ L_{u,1} = (1-u)\cdot P_{0,1} + u \cdot P_{1,1} $$ $$ L_{v,0} = (1-v)\cdot P_{0,0} + v \cdot P_{0,1} $$ $$ L_{v,1} = (1-v)\cdot P_{1,0} + v \cdot P_{1,1} $$ $$ P = (1-v)\cdot L_{u,0} + v \cdot L_{u,1} $$ $$ p = \Pi(P) $$ $$ P_u = L_{v,1} - L_{v,0} $$ $$ p_u = m \frac{{P_u} \cdot P_w - P \cdot {P_u}_w}{P_w^2} $$ $$ P_v = L_{u,1} - L_{u,0} $$ $$ p_v = m \frac{{P_v} \cdot P_w - P \cdot {P_v}_w}{P_w^2} $$

Now the algorithm by Sederberg, we need to obtain the four points $P^{00}$, $P^{01}$, $P^{10}$, $P^{11}$. To keep the notation from Sederberg's paper the input parameters are $s$, $t$ instead of $u$, $v$ from the last algorithm, sorry for the confusion.

For $s$ and $t$ between $0$ and $0.5$ the equation is: $$ {{P}^{kl}(s,t)}=\hat{{P}}^{kl}(u,v) = \sum_{i=k}^{m+k-l} \sum_{j=l}^{n+l-1} \hat{{P}}^{kl}_{i,j} u^{i-k} v^{j-l}; \quad k,l = 0, 1 $$ where $u = \frac{s}{1-s}$ and $v = \frac{t}{1-t}$ and $\hat{{P}}^{kl}_{i,j} = {m-1\choose i-k} {n-1\choose j-l} {P}_{i,j}$

For $s$ and $t$ between $0.5$ and $1$ the equation is: $$ {{P}^{kl}(s,t)}=\hat{{P}}^{kl}(u,v) = \sum_{i=k}^{m+k-l} \sum_{j=l}^{n+l-1} \hat{{P}}^{kl}_{m + 2k - i - 1,j} u^{i-k} v^{j-l}; \quad k,l = 0, 1 $$ where $u = \frac{1-s}{s}$ and $v = \frac{1-t}{t}$ and $\hat{{P}}^{kl}_{i,j}$ is same as above.

Note that you need to combine between these two equations when for example $s$ is greater than $0.5$ and $t$ is less than $0.5$.

$$ P = (1-s)(1-t)P^{00} + s(1-t)P^{10} + (1-s)tP^{01} + stP^{11} $$

If $s > 0.5$: $$ P_t = (1-t)P^{00} + tP^{01} $$ $$ p_t = -\frac{m}{s}\frac{{P_t}_w}{P_w} \cdot (p_t - p) $$ else $$ P_t = (1-t)P^{10} + tP^{11} $$ $$ p_t = \frac{m}{1-s}\frac{{P_t}_w}{P_w} \cdot (p_t - p) $$

If $t > 0.5$: $$ P_s = (1-s)P^{00} + sP^{10} $$ $$ p_s = -\frac{n}{t}\frac{{P_s}_w}{P_w} \cdot (p_s - p) $$ else $$ P_s = (1-s)P^{01} + sP^{11} $$ $$ p_s = \frac{n}{1-t}\frac{{P_s}_w}{P_w} \cdot (p_s - p) $$

Answer 2

NURBs are vastly over-rated, in my view. Since they're rational, computing derivatives is a mess. Computing integrals is even worse. If you have the option, I'd strongly recommend sticking with polynomials. People will tell you that polynomials can't represent conics and quadrics exactly. This is true, but we're using computers, so nothing is exact. And polynomial approximations of conics and quadrics are perfectly adequate for any application that I'm aware of. They're certainly good enough for graphics

I don't know if Newton's method is well suited to GPU programming. It has too many special cases. It might be better to tesselate finely and then shoot rays at triangles, which is easy.

Seems like there are two somewhat independent pieces of your work. To figure out what's going wrong, I'd recommend that you test them independently. First, test the Newton-Raphson code using simple slow Bezier position/derivative functions that are known to work correctly. The NURBs Book by Tiller and Piegl has decent code that you might be able to use. When that works properly, you can start playing with faster (but more complex) position/derivative functions. You will have a benchmark against which you can test them.

Another suggestion: compute derivatives using finite-difference approximations. These will be plenty good enough for graphics applications. Even if you don't end up using these functions, they will give you a way to test your more complex position/derivative functions.