TP3 : Bsplines, De Boor's algorithm
Code
Update your local repo
git pull
or clone everything if needed
git clone https://github.com/bbrrck/geonum2016.git
As usual,
cd TP3/
mkdir build
cd build
cmake ..
make
./geonum_TP3 simple
We’re still using Eigen so keep the quick reference open and ready.
Python script for rendering is included in the plots/
folder.
You can now pass the data name as an argument
python ../plots/plot.py spiral
If you want to use gnuplot, you will need to modify the script from TP1 yourself.
Bsplines
Given the degree $k$, $n+1$ control points $\mathbf d_0,\dots,\mathbf d_n$, and the knot vector $t_0 \leq t_1 \leq \dots \leq t_m$ with $m = n+k+1$, the Bspline curve $S(t)$ is defined as
The $N_j^k$ are the recursivelydefined basis functions (hence the name Bspline)
Looks complicated? Don’t worry if you cannot get your head around all the indices and whatnot; De Boor is here to help you!
Bspline basis functions $N^k_j$ up to degree 5 for the knot sequence $(0,1,2,3,4,5,6,7)$.
De Boor’s algorithm
… also called the De BoorCox algorithm. It can be seen as the generalization of the de Casteljau. (A Bézier curve is in fact a special case of Bspline.)
 input :
degree $k$
$n+1$ control points $\mathbf{d_{0}},\dots,\mathbf{d_{n}}$
knot vector $t_0 \leq t_1 \leq \dots \leq t_m$ with $m = n+k+1$
parameter $t \in [t_i, t_{i+1}) \quad k \leq i \leq n$  output : The point $\mathbf S(t) = \mathbf d_j^k$ on the curve.
 compute : For $j=ik, \dots, i,$ set $\mathbf d_j^0 = \mathbf d_j$. Then compute the points \begin{align} \mathbf d_{j}^{r} &= (1w_{j,k(r1)}) \mathbf d_{j1}^{r1} + w_{j,k(r1)} \mathbf d_{j }^{r1} \end{align} for \begin{align} \quad r = 1,\dots,k, \quad \quad j = ik+r,\dots,i \end{align} with \begin{align} w_{j,k(r1)} &= \frac{ t  t_j }{ t_{j+k(r1)}  t_j }. \end{align}
Be careful with the indices! Here I’ve expressed the point at depth $r$ in terms of the points at depth $r1$; that is why there is the $r1$ everywhere in the formula. (It becomes much more elegant if we express $r+1$ in terms of $r$.) This might be a bit annoying, but I think it’s also more practical for the implementation.
A cubic Bspline with 16 segments and endpoint interpolation.
ToDo
 Implement the De Boor’s algorithm. (
Vec2 DeBoor
)  Implement the evaluation of a Bspline curve. (
evaluateBSpline
)  Compute the Bspline curve for the
simple
dataset. Then, modify the knot vector and recompute. What happened?  Compute the Bspline curve for the
spiral
dataset. Try using the knot vector0 0 0 0 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5
. What changed?  Compute the Bspline curve for the
camel
dataset. Move the front leg by changing the xcoordinate of the very last control point to1.5
. Which segments of the curve have changed? Why?
NURBS
If you’ve done the tests with the circle.bspline
dataset, you might have been disappointed,
as the resulting curve is far from being a nice circle.
The truth is, it’s mathematically impossible to represent the unit circle as a Bspline.
However, it is possible (and not so hard) using a further generalization of the Bspline curves: the nonuniform rational Bsplines or NURBS.
If you examine the circle.nurbs
file, you’ll find it’s not so different from the circle.bspline
.
The only thing that’s changed is the additional number after each control point.
It almost looks like the control points don’t live in 2D, but in 3D.
If that thought occured to you, you’re not far from the truth:
together, those three numbers define the homogeneous (or projective) coordinates of a control point.
The third number is often called a weight of the control point.
Here’s the good news: even with the homogeneous coordinates, we can apply exactly the same De Boor’s algorithm without any modifications. So, if you’re up for it, here’s a secret recipe for transforming your 2D Bsplines into 2D NURBS:

Use a matrix of type
MatX3
to store the control points (the first and second column are the x and y coordinates, respectively; the third column hold the weights). For reading thecircle.nurbs
file, you can still use the methodreadBSpline
, no change here; just make sure the second argument you pass is a threecolumn matrix. Let’s call this matrixControlPoints3
. 
DeBoor3:
you can copypaste everything from the 2D version, just be sure to change the types toVec3
andMatX3
. 
Before you begin the evaluation, you need to multiply the x and y coordinates (first two columns of the
ControlPoints3
matrix) by the weights (third column). You can do this withControlPoints3.col(i) = ControlPoints3.col(i).array() * ControlPoints3.col(2).array();
fori=0,1.
The.array()
tells the Eigen to perform the operation elementwise. 
Proceed by evaluating the
SplinePoints3
with three coordinates. 
The last step is to return to the plane coordinates. To do that, you need to divide by the weights. Much like before, this is done with
SplinePoints2.col(i) = SplinePoints3.col(i).array() / SplinePoints3.col(2).array();
Bonus ToDo
 Implement the NURBS curve by extending your Bspline methods.
 Compute and visualise the circle as a NURBS curve.
Resources
 Bspline and De Boor’s algorithm
 1.4.2 Bspline curve and 1.4.3 Algorithms for Bspline curves, online chapters from the book Shape Interrogation for Computer Aided Design and Manufacturing by N. Patrikalakis, T. Maekawa & W. Cho
 NURBS on wikipedia (includes the circle example)
 homepage of Prof. de Boor