User:Raul Miller/NURBS

I am trying to work up to dealing with non-uniform b-splines. I am struggling with terminology and my inability to find good examples that reflect the terms I am trying to understand.

[Uniform] B-Splines

First, b-splines are based on convolution. Convolution is a math concept where we combine the results of two functions: an impulse function (g) and an impulse response function (h). Ideally, J would supply a convolution operator, but no one has implemented one yet. So... hypothetically speaking, an impulse represents some function which depends on time (or whatever), and impulse result can be described the same way. The difference is that the impulse function is a "cause" (it happens first) and impulse response is an "effect".

Here's an approximation of what convolution looks like

${\begin{aligned}f(t)=\int \!g(s)h(t-s)\,\mathrm {d} s\end{aligned}}$

Here, parenthesis show that the enclosed value is an argument of the function whose name appears on the left of the parenthesis. t and s are both representations of time (or whatever) with t being a free variable and s being a variable of integration (where we are adding up the results of this expression for all possible values of s), and ds is the size of the distance between s values (in the ideal case it's so close to zero that we can't tell the difference). And, of course, integral means "sum".

Meanwhile, the impulse response function for b-splines is:

   h=: 0&<: * 1&>

If C is a convolution operator, the impulse functions for b-splines would be: h`(h C h)`(h C h C h)`(h C h C h C h)`...

So, this means that we can get an idea of what b-spline basis functions look like using:

   require'plot'
   dt=: 0.01
   t=: _1 +  dt * i.9%dt
   plot t;B0=: h t
   plot t;B1=: +/ B0 * dt * h t-~/t
   plot t;B2=: +/ B1 * dt * h t-~/t
   plot t;B3=: +/ B2 * dt * h t-~/t

Or, all together:

   plot t;B0,B1,B2,:B3

And, since the area under the curve B0 is a unit square, we expect all of these to have a total area of 1:

   +/B0*dt
1
   +/B1*dt
1
   +/B2*dt
1
   +/B3*dt
1

But these were just approximations -- ideally, we want more exact numbers, and we do not want to have to go through thousands of steps to get them. In other words, take that original function h and integrate it symbolically. Unfortunately, J currently does not include the mechanics we need to do this directly on h:

   (h d. _1) 0.5
|domain error

But we can use piecewise polynomials here. In other words, we will use one polynomial for the argument range 0..1, another for the argument 1..2, another for the argument range 2..3 and so on. And when we are exactly on an integer value, we define the polynomials to be equal -- each one starts where the previous one leaves off.

It would be convenient, here, if J's indexing function would give us fills (zeros) when we ask for an out of range value -- that would give us the "zeros everywhere we do not have a definition for" character of h and its convolutions. But that's not how J works. So let's define this function first, and while we are at it, let's make it ignore any fractional part in the indices, so we can use it directly on our argument values:

   from =: (i.@#@] i. <.@[) { ({.~ (1+#))@]
   2 3.14 5 7 from 1 2 3 4
3 4 0 0

The left argument to 'from' is indices (but we only use the integer part) and we do some extra work so that otherwise undefined indices select 0s.

With this, we can define a piecewise polynomial function, based on a list of polynomials (which will be our left argument):

   f=: from~ p. 1 | ]

Now all we need is our convolutions of h. I do not have to define general purpose convolution here, I can get away with:

   cnv=: [: ((0, p.&1) 0}"_1 ,&0 - 0&,) 0&p..

In other words: 0&p.. integrates each of our polynomials (using 0 for the constant of integration, since we do not know that yet). Then we subtract neighboring polynomials (this corresponds to the unit length bounds of our integration). We also evaluate each of our integrated polynomials at the value 1 to find our constants of integration.

Now we just need the list of polynomials corresponding to our initial basis function:

   B0p=: ,. 1

Here, B0p corresponds to B0, above, but it's a "sequence" of piecewise "polynomials". (It's just one polynomial and it's just the constant 1.)

And, we can define a routine to get the nth list of polynomials:

   b=: cnv@]^:[&B0p

Here's an example:

   b 3
       0    0   0  0.166667
0.166667  0.5 0.5      _0.5
0.666667    0  _1       0.5
0.166667 _0.5 0.5 _0.166667

This example has four polynomials, the first polynomial is one sixth of x cubed (where x is variable of the polynomial), and the last polynomial (used for values between 3 and 4) is: ((one sixth) plus (negative half of x) plus (half of x squared) plus (negative one sixth of x cubed)).

and we can plot an arbitrary sequence of basis functions:

   plot t;(b i.8) f"2 1 t

Non-Uniform B-Splines

(I am currently working out what the terminology about "knots" has to do with these basis functions.)

The b-splines I described above can be characterized as uniform b-spline basis functions. Basis functions can be combined with a series of control points using inner product. Additionally, the uniform b-splines can be generalized into non-uniform b-splines by replacing the integer boundary points between polynomials with an arbitrary non-decreasing knot vector. Each value in the knot vector represents a transition point between polynomial curves.

A knot almost corresponds to the boundary between piecewise polynomials. "Knots" are a sorted list of numbers (non-decreasing) which mark the "discontinuities" where we switch from one curve to its neighbor, with some rules about repeated knot values: In the above example where I showed the result of (b 3) we had knots at the values 0, 1, 2, 3 and 4. But knots are defined in terms of a the Cox deBoor algorithm rather than in terms of convolution.

Here's an implementation:

basis=:2 :0
:
  NB. x: index
  NB. m: knots
  NB. n: degree
  NB. y: parameter
  if. 1>n do.
    ((x{m)<:/y) *. ((x+1) { m,{:m)>/y
  else.
    't0 tn t1 tn1'=. 0>.(_1+#m)<.((,1+])0,n)+/x
    b0=. x m basis (n-1) y
    b1=. (1+x) m basis (n-1) y
    (b0 * (t0-~/y)%tn-t0) + b1 * (tn1-/y)%tn1-t1
  end.
)

NURB=:2 :0
  NB. m: knots
  NB. n: control points
  NB. y: parameter
  n+/ .*(i.#n) m basis (_1+m-&#n) y
)

Note that the order of the b-spline basis functions used to implement the NURB is determined by the length of the knot and control vectors. For example:

plot (;0 1 2 3 4 5 NURB 1 0) 0.01*i.400

represents a degree 3 basis function:

Warning: this rest of this page is hasty and incomplete -- it's something I wrote years ago, and am not yet prepared to delete.

NURBS stands for "Non-Uniform Rational B-Splines".

Some j wiki pages use NURBS but they do not include any documentation: Studio/OpenGL/Teapot and Studio/OpenGL/BraidKnot

To understand NURBS I should first understand Bezier curves and B-Splines.

Here's an example of using a one-dimensional (rank zero coordinates) Bezier curve

   Bz=: [ +/@:* (i. ! <:)@#@[ * ] (^~/~ * -.@[ ^~/~ |.@]) i.@#@[
   require'plot'
   _10 20 _20 10 Bz 0.01*i.101

The domain of Bezier curves is real numbers 0 through 1. The range of Bezier curves over that domain is {.x through {:x. Look at some plots for further insight.

TODO: Bezier curves can be subdivided into a sequence of Bezier curves. This requires a treatment of Bezier curves with arbitrary end points.

   Bez=:2 :'[ Bz ((n-m) %~ m -~ ])'

The domain of m Bez n is m through n.

TODO: implement b-splines in one dimension (NURBS seem to imply 3 or 4 dimensional coordinates with the fourth value in each coordinate being 1 in some expressions), then go back and do all of these for arbitrary dimensions (rank 1 coordinates)

B-Splines have a basis from the Cox de Boor formula:

Cox_deBoor=:4 : 0"1 0
  NB. x: knots
  N=.,:(x <: y) *. 1 |. y < x
  for_p. }.i.(#x)-I.{.N do.
     Ni=. {:N
     xp1=. (p+1)|. x
     N=.N,(Ni*(x-y)%x-p|.x) + (}.Ni,0)*(xp1-y)%xp1-1|.x
  end.
)

NB. example use:
   1 2 3 4 5 Cox_deBoor 1.5
   1 2 3 4 5 Cox_deBoor 3.5

Informally, the values in this matrix represent interpolation of the previous rows values for the given value of y.

That's a rather awkward expression, and I should clean it up (for example, note that values rotated around the end by |. do not matter here as they will always be multiplied by 0).

Uniform B-Splines have uniform spacing between the knots. In other words, x would be replaced by x0 + (x1-x0) * i.n where x0 and x1 are the first two knot values.

Cox_deBoorU=:4 :0"1 0
  'x0 x1 n'=. x
  d=.x1-x0
  (x0+d*i.n+<.(y-x0)%d) Cox_deBoor y
)
   1 2 7 Cox_deBoorU 6.5 8.5 NB. seven row Cox de Boor bases

This result looks slightly odd, and seems to conflict with my informal description, above.

   1 2 7 Cox_deBoorU 5

However, when compared with nearby values, that result seems appropriate:

   1 2 7 Cox_deBoorU 4.9 5 5.1

See also: Catmull Clark subdivision surface

---