function spline()

knot_vector=[1,1,1,2,3,3,4,5,6,7,7,7]

precision = 0.0001


compute_nr_basis_functions = @(knot_vector,p) size(knot_vector, 2) - p - 1
mesh   = @(a,c) [a:precision*(c-a):c]


p = compute_p(knot_vector)
t = check_sanity(knot_vector,p)
nr = compute_nr_basis_functions(knot_vector,p)

x_begin = knot_vector(1)
x_end = knot_vector(size(knot_vector,2))



x=mesh(x_begin,x_end);


y=compute_spline(knot_vector,p,1,x);
plot(x,y)

hold on

for i=2:nr
  y=compute_spline(knot_vector,p,i,x);
  plot(x,y)
endfor

hold off





function p=compute_p(knot_vector)

  initial = knot_vector(1)
  kvsize = size(knot_vector,2)
  i=1

  while (i+1 < kvsize) && (initial == knot_vector(i+1))
    i=i+1
  endwhile
  
  p = i-1
  
  return
endfunction



function y=compute_spline(knot_vector,p,nr,x)
  
  fC= @(x,a,b) (x)/(b-a)-a/(b-a);
  fD= @(x,c,d) (1-x)/(d-c)+(d-1)/(d-c);
  
  
  a = knot_vector(nr);
  b = knot_vector(nr+p);
  c = knot_vector(nr+1);
  d = knot_vector(nr+p+1);

  if (p==0)
    y = 0 .* (x < a) + 1 .* (a <= x & x <= d) .+ 0 .* (x > d);
    return
  endif
  
  lp = compute_spline(knot_vector,p-1,nr,x);
  rp = compute_spline(knot_vector,p-1,nr+1,x);
  
  if (a==b)
    y1 = 0 .* (x < a) + 1 .* (a <= x & x <= b) + 0 .* (x > b);
  else
    y1 = 0 .* (x < a) + fC(x,a,b) .* (a <= x & x <= b) + 0 .* (x > b);
  endif
  
  if (c==d)
    y2 = 0 .* (x < c) + 1 .* (c < x & x <= d) + 0 .* (d < x);
  else
    y2 = 0 .* (x < c) + fD(x,c,d) .* (c < x & x <= d) + 0 .* (d < x);
  endif
  
  y = lp .* y1 .+ rp .* y2;
  return
  
endfunction


function t=check_sanity(knot_vector,p)

  initial = knot_vector(1)
  kvsize = size(knot_vector,2)

  t = true
  counter = 1

  for i=1:p+1
    if (initial != knot_vector(i))
      t = false
      return
    endif
  endfor


  for i=p+2:kvsize-p-1
    if (initial == knot_vector(i))
      counter = counter + 1
      if (counter > p)
        t = false  
        return
      endif
    else
      initial = knot_vector(i)
      counter = 1
    endif
  endfor

  initial = knot_vector(kvsize)

  for i=kvsize-p:kvsize
    if (initial != knot_vector(i))
      t = false
      return
    endif
  endfor
  
  for i=1:kvsize-1
    if (knot_vector(i)>knot_vector(i+1))
      t = false
    endif
  endfor

  return

endfunction



endfunction