B-Spline - Problème aux extrémités

Bonjour,

Je suis entrain d'écrire un bout de code pour calculer / dessiner des courbes splines et j'ai un soucis dans le cas particulier où la multiplicité des premiers et derniers noeuds est égale au degré+1 (i.e. la condition pour que la courbe passe par le premier et le dernier point de contrôle).

Si U est le vecteur de noeuds et CP le nombre de points de contrôle, la courbe est définie sur l'intervalle [U_degré,U_CP] (sur lequel, d'ailleurs, la somme des valeurs des fonctions de base vaut 1 ... ce sont les maths).

Tout se passe bien (1ere et 2e image, lorsque degre=2 et CP=4) lorsque la multiplicité des premiers de derniers noeuds est <= degré. La somme des fonctions de la base est en particulier égale à 1 à chaque noeud. Par contre, lorsque la multiplicité est égale à degré +1, la somme n'est plus égale à 1 pour U_CP (et on observe un décalage des valeurs de U pour lesquelles la somme vaut 1 vers les index plus petits). Cette somme vaut en fait 0, ce qui signifie pratiquement que la courbe ne passe pas par le dernier point de contrôle.

Toute la littérature sur le sujet indique qu'il faut une multiplicité égale à degré+1 pour que la courbe passe par les deux points de contrôles des extrémités. Mais en fait, quand on y regarde de plus près, la valeur de U_CP vaut 1.0 et toutes les fonctions de base sont nulles pour 1.0 ... donc c'est impossible ...

Y a-t-il quelque chose que je ne capte pas ?

Christian

Le problème ne serait-il pas juste que tu n'a pas assez élargi les bornes de ton intervalle de calcul ?

Plus le degré de la spline est élevé, plus il faut de noeuds pour calculer sa valeur, et donc plus il faut repousser les bornes inf/sup de l'intervalle.

Merci pour ta réponse pseudocode ... en fait, les bornes de l'intervalle sont définies mathématiquement par [U_degre,U_CP] où CP est le nombre de points de contrôle. Comme "degre+CP=nombre noeuds-1", il y a donc "degré" noeuds après U_CP ... qui sont tous égaux à U_CP (soit 1.0). Christian

Citation:

---

Envoyé par coberle

Merci pour ta réponse pseudocode ... en fait, les bornes de l'intervalle sont définies mathématiquement par [U_degre,U_CP] où CP est le nombre de points de contrôle. Comme "degre+CP=nombre noeuds-1", il y a donc "degré" noeuds après U_CP ... qui sont tous égaux à U_CP (soit 1.0). Christian

---

J'ai du mal a suivre tes notations, par rapport a celles de wikipedia (http://fr.wikipedia.org/wiki/B-spline).

Une chose est sure, la somme des fonctions de base fait toujours 1, que ce soit aux extremités ou au milieu de ton intervalle d'évaluation.

Par exemple, pour une spline cubique (degré 3) avec "t" dans l'intervalle [0,3] et partant du point "1" pour aller au point "100", j'obtiens cela:

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14

3 points de controle à 1 3 points de controle à 100 --------------------------- --------------------------------- Spline(0,00) = 1,000*1 + 0,000*1 + 0,000*1 + 0,000*100 + 0,000*100 + 0,000*100 = 1,000000 Spline(0,30) = 0,343*1 + 0,542*1 + 0,110*1 + 0,005*100 + 0,000*100 + 0,000*100 = 1,445500 Spline(0,60) = 0,064*1 + 0,558*1 + 0,342*1 + 0,036*100 + 0,000*100 + 0,000*100 = 4,564000 Spline(0,90) = 0,001*1 + 0,331*1 + 0,547*1 + 0,122*100 + 0,000*100 + 0,000*100 = 13,028500 Spline(1,20) = 0,000*1 + 0,128*1 + 0,588*1 + 0,282*100 + 0,002*100 + 0,000*100 = 29,116000 Spline(1,50) = 0,000*1 + 0,031*1 + 0,469*1 + 0,469*100 + 0,031*100 + 0,000*100 = 50,500000 Spline(1,80) = 0,000*1 + 0,002*1 + 0,282*1 + 0,588*100 + 0,128*100 + 0,000*100 = 71,884000 Spline(2,10) = 0,000*1 + 0,000*1 + 0,121*1 + 0,547*100 + 0,331*100 + 0,001*100 = 87,971500 Spline(2,40) = 0,000*1 + 0,000*1 + 0,036*1 + 0,342*100 + 0,558*100 + 0,064*100 = 96,436000 Spline(2,70) = 0,000*1 + 0,000*1 + 0,004*1 + 0,110*100 + 0,542*100 + 0,343*100 = 99,554500 Spline(3,00) = 0,000*1 + 0,000*1 + 0,000*1 + 0,000*100 + 0,000*100 + 1,000*100 = 100,000000

---

Pseudocode, voici un exemple:
- le vecteur nodal est U[]={0,0,0,0,0.2,0.4,0.6,0.8,1,1,1,1};
- m=11 (12 noeuds)
- n=3 (degré)
- m-n-1=7

La courbe B-Spline est donc (pour t dans [0,1]) S(t)=sigma (i=0-7) [b_i_3(t)*Pi].

Lorsque t=0, b_0_3(0)=1.0 (et tous les autres b_i sont nulles). Par contre, lorsque t=1.0, b_7_3(1.0)=0.0 (et tous les autres également) ... j'imagine que tu dois obtenir b_7_3(1.0)=1.0 ... Christian

D'après la définition, les fonctions Bi,0 valent 1 si t(i) <= t < t(i+1)
L'inégalité "stricte" à droite implique que Bi,0 est nul lorsque t=t(i+1).

Il faut donc que l'intervalle contenant les noeuds t(i) soit plus grand que l'intervalle d'évaluation de la spline. Sinon, pour le dernier noeud, la spline est une somme de 0.

Citation:

---

Envoyé par coberle

Pseudocode, voici un exemple:
- le vecteur nodal est U[]={0,0,0,0,0.2,0.4,0.6,0.8,1,1,1,1};
- m=11 (12 noeuds)
- n=3 (degré)
- m-n-1=7

---

Soit c'est moi, soit tu mélanges points de contrôle et nœuds.

P[] = {0,0,0,0,0.2,0.4,0.6,0.8,1,1,1,1} // 12 points de contrôle

La Spline est de degré 3.

n=3

Le noeuds représentent "l'emplacement" des fonctions de base.
Il faut donc m = 12+n+1 = 16 noeuds.
Par exemple, uniformément répartis (spline uniforme):

T[] = { -3,-2,-1,0,1,2,3,4,5,6,7,8,9,10,11,12 }

il faut 4 noeuds (= n+1) pour évaluer un point de la spline.

Dans notre exemple, l'intervalle d'évaluation de la spline est donc de T[3] à T[12], c'est a dire 0<=t<=9

On peut faire une règle de 3 pour se ramener à [0,1].

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49

public class Spline {   double[] P = { 0,0,0,0,0.2,0.4,0.6,0.8,1,1,1,1 }; int n=3;   double[] T = { -3,-2,-1,0,1,2,3,4,5,6,7,8,9,10,11,12 }; int m=15; // T.length - 1 = P.length + n   double BEGIN = T[n], END = T[T.length-1-n];   double B(int j, int n, double t) { if (n==0) { if (t>=T[j] && t<T[j+1]) return 1; return 0; }   double Bjnt = 0;   double d0 = T[j+n]-T[j]; if (d0!=0) Bjnt += ( (t-T[j])/d0 ) * B(j,n-1,t);   double d1 = T[j+n+1]-T[j+1]; if (d1!=0) Bjnt += ( (T[j+n+1]-t)/d1 ) * B(j+1,n-1,t);   return Bjnt; }   double S(double t) { double S=0; for(int i = 0 ; i<=m-n-1; i++) S+=P[i]*B(i,n,t); return S; }   void test() { for(double t=0;t<=1.0;t+=0.1) { double u = BEGIN + (END-BEGIN)*t; double s = S(u); System.out.printf("t=%.2f --> Spline(%.2f) = %f\n",t,u,s); } }   public static void main(String[] args) { new Spline().test(); }   }

---

Citation:

---

t=0,00 --> Spline(0,00) = 0,000000
t=0,10 --> Spline(0,90) = 0,000000
t=0,20 --> Spline(1,80) = 0,017067
t=0,30 --> Spline(2,70) = 0,140900
t=0,40 --> Spline(3,60) = 0,320000
t=0,50 --> Spline(4,50) = 0,500000
t=0,60 --> Spline(5,40) = 0,680000
t=0,70 --> Spline(6,30) = 0,859100
t=0,80 --> Spline(7,20) = 0,982933
t=0,90 --> Spline(8,10) = 1,000000
t=1,00 --> Spline(9,00) = 1,000000

---

Pseudocode,
Encore merci de prendre le temps de m'expliquer ... c'est vraiment sympa. Du coup, tu va peut-être me trouver retord et être déçu de ma réponse. J'espère cependant que nous pourrons continuer la discussion, car ce n'est pas super clair.
Dans ta réponse, tu confirmes que pour le dernier noeud, la somme est 0, ce qui est bien ce que j'avais observé (ouf !)
Là où je ne suis pas certain de te suivre est sur l'aspect confusion entre point de contrôle et noeuds (et du coup, l'exemple que tu donnes me confuse pas mal, car tu donnes une multiplicité aux points de contrôle !). Voir par exemple le tableau en bas de page ici (http://pagesperso-orange.fr/reboli/b....htm#coxdeboor). Pour ma part, les points de contrôle sont des points avec des coordonnées X,Y,Z et si on veut que le courbe passe par le premier et dernier de ces points, on "joue" avec les noeuds mais pas avec ces points ...

Pour reprendre ta notation : si on a 9 points de contrôles P0-P8 distincts (avec, disons, Pi=2i), comment choisis-tu ton vecteur nodal pour que la spline passe par les deux extrémités P0=0 et P8=16 ?

Encore merci pour ton aide.
Bien cordialement,
Christian

Citation:

---

Envoyé par coberle

Là où je ne suis pas certain de te suivre est sur l'aspect confusion entre point de contrôle et noeuds (et du coup, l'exemple que tu donnes me confuse pas mal, car tu donnes une multiplicité aux points de contrôle !). Voir par exemple le tableau en bas de page ici (http://pagesperso-orange.fr/reboli/b....htm#coxdeboor).

---

C'est une variante de la définition de la spline. On peut soit :

* imposer que le noeuds ne soient pas confondus:

t0 < t1 < ... < tm-1 (:arrow: http://fr.wikipedia.org/wiki/B-spline)

Auquel cas, il faut jouer sur la multiplicité des poles (points de controle) pour assurer les points de passages


* Ou alors permettre que les noeuds soient confondus:

t0 <= t1 <= ... <= tm-1 (:arrow: http://en.wikipedia.org/wiki/B-spline)

Auquel cas, il faut jouer sur la multiplicité des noeuds pour assurer les points de passages

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52

public class Spline {   double[] P = { 0,0.2,0.4,0.6,0.8,1 }; // 6 poles int n=3;   double[] T = { 0,0,0,0,1,2,3,3,3,3 }; // 10 noeuds int m=9; // T.length - 1 = P.length + n   /*** Le reste du code ne change pas ***/   double BEGIN = T[n], END = T[T.length-1-n];   double B(int j, int n, double t) { if (n==0) { if (t>=T[j] && t<T[j+1]) return 1; return 0; }   double Bjnt = 0;   double d0 = T[j+n]-T[j]; if (d0!=0) Bjnt += ( (t-T[j])/d0 ) * B(j,n-1,t);   double d1 = T[j+n+1]-T[j+1]; if (d1!=0) Bjnt += ( (T[j+n+1]-t)/d1 ) * B(j+1,n-1,t);   return Bjnt; }   double S(double t) { double S=0; for(int i = 0 ; i<=m-n-1; i++) { S+=P[i]*B(i,n,t); } return S; }   void test() { for(double t=0;t<=1.0;t+=0.1) { double u = BEGIN + (END-BEGIN)*t; double s = S(u); System.out.printf("t=%.2f --> Spline(%.2f) = %f\n",t,u,s); } }   public static void main(String[] args) { new Spline().test(); }   }

---

Citation:

---

t=0,00 --> Spline(0,00) = 0,000000
t=0,10 --> Spline(0,30) = 0,155250
t=0,20 --> Spline(0,60) = 0,270000
t=0,30 --> Spline(0,90) = 0,357750
t=0,40 --> Spline(1,20) = 0,431600
t=0,50 --> Spline(1,50) = 0,500000
t=0,60 --> Spline(1,80) = 0,568400
t=0,70 --> Spline(2,10) = 0,642250
t=0,80 --> Spline(2,40) = 0,730000
t=0,90 --> Spline(2,70) = 0,844750
t=1,00 --> Spline(3,00) = 0,000000 /* <-- ton problème actuel :aie: */

---

On notera que les valeurs obtenues ne sont pas les mêmes, surtout aux extrémités (du fait que dans ce cas la spline n'est pas uniforme)

Citation:

---

Envoyé par coberle

Pour reprendre ta notation : si on a 9 points de contrôles P0-P8 distincts (avec, disons, Pi=2i), comment choisis-tu ton vecteur nodal pour que la spline passe par les deux extrémités P0=0 et P8=16 ?

---

P = {0,2,4,6,8,10,12,14,16} // 9 poles

Soit "n" le degré de la spline

Il faut une multiplicité mininum de "n" sur un point de controle pour que la spline passe par ce point.

Code:

---

1 2	n=3 P = {0,0,0,2,4,6,8,10,12,14,16,16,16} // 13 poles

---

Il nous faut alors P.length+n+1 = 17 noeuds

Code:

---

1 2	T = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17} m = 16 // T.length - 1 = P.length + n

---

Et evaluer la spline entre T[n]=4 et T[T.length-1-n]=14

Citation:

---

Spline( 4,00) = 0,000000
Spline( 4,50) = 0,041667
Spline( 5,00) = 0,333333
Spline( 5,50) = 1,041667
Spline( 6,00) = 2,000000
Spline( 6,50) = 3,000000
Spline( 7,00) = 4,000000
Spline( 7,50) = 5,000000
Spline( 8,00) = 6,000000
Spline( 8,50) = 7,000000
Spline( 9,00) = 8,000000
Spline( 9,50) = 9,000000
Spline(10,00) = 10,000000
Spline(10,50) = 11,000000
Spline(11,00) = 12,000000
Spline(11,50) = 13,000000
Spline(12,00) = 14,000000
Spline(12,50) = 14,958333
Spline(13,00) = 15,666667
Spline(13,50) = 15,958333
Spline(14,00) = 16,000000

---

Pseudocode,
J'avais regardé en fin de matinée et tu arrivais bien à faire passer la courbe par ses deux extrémités avec une multiplicité sur les noeuds. Du coup, j'étais assez content. Puis, tu as apporté une modification à tu as maintenant 0 (là où tu avais précédemment 1), i.e. la courbe ne passe plus par l'une de ses extrémités ??? ... mais ton commentaire montre que tu as bien saisi mon problème (que tu résouts au passage en prenant une multiplicité sur les pôles). En ne considérant que la multiplicité sur les noeuds, as-tu modifié quelque chose à ton code pour que tu obtiennes 0 (là où tu avais 1 avant) ? Christian

Citation:

---

Envoyé par coberle

Pseudocode,
J'avais regardé en fin de matinée et tu arrivais bien à faire passer la courbe par ses deux extrémités avec une multiplicité sur les noeuds. Du coup, j'étais assez content. Puis, tu as apporté une modification à tu as maintenant 0 (là où tu avais précédemment 1), i.e. la courbe ne passe plus par l'une de ses extrémités ??? ... mais ton commentaire montre que tu as bien saisi mon problème (que tu résouts au passage en prenant une multiplicité sur les pôles). En ne considérant que la multiplicité sur les noeuds, as-tu modifié quelque chose à ton code pour que tu obtiennes 0 (là où tu avais 1 avant) ? Christian

---

Mon programme affiche "1" a cause des arrondis de calculs sur les "double" (la dernière valeur de "t" n'est pas 1.0 mais 0.999...). La dernière valeur de la Spline vaut donc bien mathématiquement zéro. C'est pourquoi j'ai édité pour ne pas te laisser de faux espoirs. :aie:

Cela dit, si c'est juste cette derniere valeur qui te pose problème, tu peux simplement ajouter un EPSILON sur le dernier noeud :

double EPSILON = 1E-10;
double[] T = { 0,0,0,0,1,2,3,3,3,3+EPSILON }; // 10 noeuds


Ou alors utiliser l'autre définition de la spline (uniforme, avec point de controle redondant) dans lequel ce problème ne se pose pas

Merci beaucoup pour ton aide. Après avoir passé de (trop) nombreuses heures sur ce problème (qui en fait n'en est pas un mais est complètement occulté par l'ensemble de la littérature sur le sujet) je me sens moins seul maintenant dans la mesure où tu confirmes ce que j'ai aussi observé. Ouf ... je vais pouvoir dormir tranquille (j'ai cherché comme un dingue une erreur dans mon code) ! Encore merci pour ton aide ! Christian :D