Discussion :

[coberle]

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

[pseudocode]

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.

[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

[pseudocode]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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

[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

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

[pseudocode]

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.

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 java :

Sélectionner tout - Visualiser dans une fenêtre à part

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();
	}
 
}

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

[coberle]

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

[pseudocode]

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 ( 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 ( http://en.wikipedia.org/wiki/B-spline)

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

Code java :

Sélectionner tout - Visualiser dans une fenêtre à part

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();
	}
 
}

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 */

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)

[pseudocode]

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 java :

Sélectionner tout - Visualiser dans une fenêtre à part

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 java :

Sélectionner tout - Visualiser dans une fenêtre à part

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

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

[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

[pseudocode]

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.

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

[coberle]

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