function  nucs(  )
%Non-uniform cubic B-splines  NON CUMULATIVES
%b. Vandeportaele
%d'après http://graphics.cs.cmu.edu/nsp/course/15-464/Fall05/papers/kimKimShin.pdf
close all
clear all

%%%%IMPORTANT%%%%%%%%%
%pour spline linéaire
%k=2
%pour spline quadratique
%k=3
%pour spline cubique
k=4;% pour avoir 3 fonctions de base non nulles entre t=3 et t=4. Il s'agit de B0 à B3
%j'affiche les fonctions de base en surnombre, mais en fait il suffit de
%regarder les 4 fonctions de base sur l'intervale t=3 à t=4. et de
%substituer les indices pour les P et les t
for numerotest=1:2
    figure
    hold on
    axis equal
    if numerotest==1
        %pour distribution uniforme
        title('pour distribution uniforme');
        P0=2; t0=0;
        P1=3; t1=1;
        P2=4; t2=2;
        P3=5; t3=3;
        %more points
        P4=6; t4=4;
        P5=7; t5=5;
        P6=8; t6=6;
        P7=9; t7=7;
        P8=10; t8=8;
        P9=11; t9=9;
    else
        %pour distribution non uniforme
        title('pour distribution non uniforme');
        t1=1.1;
        t2=2.9;
    end
    %rangement dans des vecteurs
    listt=[t0;t1;t2;t3;t4;t5;t6;t7;t8;t9];
    listp=[P0;P1;P2;P3;P4;P5;P6;P7;P8;P9];
    %affichage des points de contrôle
    plot(listt, listp,'ko');
    Base=zeros(4,1);
    for t=0:0.01:9
        pt=0;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %pour indices de 0 à +3: plus simple à comprendre mais ne gère pas
        %bien les 3 premiers points de la spline; ce sera la version à
        %implémenter lorsque l'on changera les indices au fûr et à mesure
        %for i=0:5
        %    Base(i +1) = B( i,k,t,listt);
        %    pt=pt+listp(i +1)*Base(i +1);
        %end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %pour indices de -2 à +5: permet une visu correcte de la courbe
        %pour qu'elle passe par les points de contrôle
        %je suis obligé de ranger les fonctions de bases en décalé...
        %je n'affiche pas ces 2 fonctions de base B-2 et B-1 sur la figure
        %mais elles sont utilisées pour le calcul de pt
        if k==2
            indiceBaseMin=-1;
        elseif k==3
            indiceBaseMin=-2;  %ce n'est pas correct, il y a un offset... avec -1 aussi... ce n'est pas forcément grave car on s'intéresse aux splines cubiques
        else %k==4
            indiceBaseMin=-2;
        end
        for i=indiceBaseMin:5
            Base(-indiceBaseMin+ i +1) = B( i,k,t,listt);
            pt=pt+listp(-indiceBaseMin+ i +1)*Base(-indiceBaseMin+ i +1);
        end
        %Base
        %sum(Base)
        plot(t,Base(1),'gx');
        plot(t,Base(2),'bx');
        plot(t,Base(3),'kx');
        plot(t,Base(4),'mx');
        %affiche les fonctions de bases suivantes...
        plot(t,Base(5),'yx');
        plot(t,Base(6),'cx');
        %la courbe interpolée
        plot(t,pt,'r+');
    end
    legend('PI','B0','B1','B2','B3','B4','B5','p(t)');
    ylabel('p');
    xlabel('t');
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ val ] = B( i,k,t,listt)
%for first index=1 because of matlab
if (i>=0)
    ti=listt(i           + 1);
else
    ti=listt(1);
end


if (i+1>=0)
    tip1=listt(i+1       + 1);
else
    tip1=listt(1);
end
if (i+k>=0)
    tipk=listt(i+k       + 1);
else
    tipk=listt(1);
end
if (i+k-1>=0)
    tipkm1=listt(i+k-1   + 1);
else
    tipkm1=listt(1);
end
if k==1
    if  (ti<= t) && (t<tip1)
        val =  1 ;
    else
        val = 0 ;
    end
    return;
else
    if (tipkm1 -ti) ~=0
        val1= ((t- ti) / (tipkm1 -ti)) * B( i,k-1,t,listt);
    else
        val1= 0;
    end
    if (tipk - tip1)~=0
        %version http://graphics.cs.cmu.edu/nsp/course/15-464/Fall05/papers/kimKimShin.pdf
        val2= ((tipk - t) / (tipk - tip1)) * B( i+1,k-1,t,listt);
        %version http://www.inf.ed.ac.uk/teaching/courses/cg/d3/nonbspline.html
        %il y a une erreur!
        %val2= ((tipk - t) / (tipkm1 - ti)) * B( i+1,k-1,t,listt);
        %autre version foireuse en (4.7 de https://theses.lib.vt.edu/theses/available/etd-100699-171723/unrestricted/chapter4.PDF)
        %val2= ((tip1 - t) / (tipk - tip1)) * B( i+1,k-1,t,listt);
    else
        val2= 0;
    end
    val=val1+val2;
    return;
end
end