Méthode Level Set

Version imprimable

17/02/2009, 22h04
menakikou

Méthode Level Set

bonjour à Tous,
Cela fait deux jours que je tourne en rond et je n'arrive pas à comprendre la méthode des level Set, j'ai lu pleins de Doc la dessus mais je comprend plus rien quand ca parle que l'on rajoute une nouvelle dimension.

EN fait , j'ai déja fait une étude sur les méthodes contour actif ( les snakes) que j'ai bien compris mais la la méthode level set je suis perdu.

Moi je voudrai savoir par exemple , j'ai une image 2D ou il y a deux cercles, comment ca se passe le level set , alors apparament il y a un niveau zéro que représente - il ? ensuite il faut mettre des 0 et -1 je ne sais ou et ensuite distance euclidienne.
SVP si quelqu"un pouvait m'aider ca serai trés trés gentil , je mets beaucoup d'effort pour comprendre ce concept mais j'y n'arrive pas :cry:.

Merci d'avance .
17/02/2009, 23h14
pseudocode

DISCLAMER: Attention, ce que je vais dire est totalement faux d'un point de vue strictement, fondamentalement et purement mathématique (topologie, géodésie, tout ca...). Donc inutile de m'insulter :aie:. C'est juste pour "aider" à la compréhension.

Si on considère que l'image est un terrain en 3D (carte d'élévation), alors le snake est une courbe 3D fermée qui se déplace sur la surface 3D jusqu'a trouver un équilibre. Le terrain 3D ne bouge pas, ce sont les points du snake qui bougent dans les 3 dimensions (x,y,z).

Pour le level-set c'est un peu le contraire. Les points du snake restent à l'altitude z=0 et ce sont les points du terrain dont on modifie l'altitude. Par effet de conséquence, les points du snake bougent aussi mais seulement sur (x,y) car leur z reste à 0.

Par la magie des level-set, TOUS les points qui sont à l'altitude z=0 font partie du snake ! Cela peut donc créer des courbes complexes (fermées, ouvertes, multiples, croisées, ...).

Est-ce plus clair ?
18/02/2009, 02h00
menakikou

Oui c plus ou moins clair, j'ai vu pour segmenter une image 3D il faut ajouter une dimension donc ca serai 4 !!!!
Donc notre cas , c'est une image 2D , donc on rajoute une coordonnés z pour notre modele deformable donc ,Le snake ne bouge pas au niveau du Z il reste à zero, seulement son (x,y) change, et comment peux se faire la segmentation d'un objet dans cette image 2D alors ? je vois pas trop comment on peut faire de la segmentation avec les level set, on me parle de distance euclidienne , oui mais une distance entre koi et koi ?
parce que au final, on fait bouger le snake dans une scene en 3D mais l'objet a segmenter se trouve dans une image 2D ?
Bizarre , franchement je dois énerver certains:aie: mais je vois pas de segmentation !
18/02/2009, 10h10
pseudocode

Citation:

Envoyé par menakikou

Donc notre cas , c'est une image 2D , donc on rajoute une coordonnés z pour notre modele deformable donc ,Le snake ne bouge pas au niveau du Z il reste à zero, seulement son (x,y) change, et comment peux se faire la segmentation d'un objet dans cette image 2D alors ? je vois pas trop comment on peut faire de la segmentation avec les level set, on me parle de distance euclidienne , oui mais une distance entre koi et koi ?
parce que au final, on fait bouger le snake dans une scene en 3D mais l'objet a segmenter se trouve dans une image 2D ?
Bizarre , franchement je dois énerver certains:aie: mais je vois pas de segmentation !

Les snake, comme les level-set, séparent (segmentent) une image en 2 parties:

- Pour les snakes, il y a les points à l'intérieur de la courbe fermée, et les points à l'extérieur. La frontière est constituée des points qui sont sur la courbe.

- Pour les level-set, il y a les points avec une altitude négative, et les points avec une altitude positive. La frontière est constituée des points qui sont à l'altitude 0.
02/04/2009, 22h52
hk0006

bonjour,

un peut d'explication sur cet méthode svp

Méthode Level Set

merci
03/04/2009, 09h23
cjuliard

Quelques détails supplémentaires

Basiquement, la méthode des level sets permet de faire évoluer une courbe ou une surface fermée de manière quelconque (avec une vitesse "locale" quelconque). Mais au lieu de modéliser explicitement cette courbe/surface par un ensemble de points connectés (courbe) ou par un maillage 3d (surface), la méthode level set travaille sur une représentation "implicite" de ces objets.Cette représentation implicite consiste généralement en une fonction "distance" à la courbe/surface considérée.
Par exemple, pour faire évoluer une courbe fermée par level sets, on va faire évoluer la fonction distance d(x,y) à cette courbe, définie comme une image 2D, plutôt que les points de la courbe elle-même. Pour chaque pixel (x,y), la valeur de l'image d(x,y) correspond tout bêtement à la distance euclidienne du point (x,y) à la courbe considérée.

Pour récupérer explicitement les points de la courbe à un instant t de l'évolution, il suffit d'extraire les points (x,y) vérifiant d(x,y) = 0 (isophote 0 de la fonction distance). Cette conversion "implicite->explicite" se fait en général par une méthode de type 'marching squares' (pour une courbe) ou 'marching cubes' (pour une surface).

Pourquoi s'embêter la vie avec ces fonctions implicites ? Parce que il est en réalité beaucoup plus facile de faire évoluer une fonction implicite, que de faire évoluer explicitement les points d'une courbe ou d'une surface, avec tous les problèmes que cela pose : changement de topologie, nécéssité de remaillage, etc..

A noter que l'utilisation de la fonction distance comme fonction implicite de représentation de la courbe/surface n'est pas obligatoire. N'importe quelle représentation telle que la courbe soit l'isophote de valeur 0 peut-être utilisée. Néanmoins, la fonction distance a de bonne propriétés géométriques et son utilisation permet de convertir facilement la vitesse d'évolution de la courbe/surface d'une forme explicite à une forme implicite, et permet aussi de calculer aisément la courbure locale K de la courbe/surface, qui est une caractéristique qui apparait très souvent dans les termes de vitesse d'évolutions de courbes/surfaces (en segmentation notamment).

D'un point de vue logiciel, beaucoup de bibliothèques implémentent la méthode des level sets. Pour ma part, je conseille fortement d'aller fouiller du côté de CImg, notamment avec les sources suivantes :

http://cimg.cvs.sourceforge.net/view..._levelsets.cpp

(évolution d'une courbe)

et

http://cimg.cvs.sourceforge.net/view...evelsets3d.cpp

(évolution d'une surface)

qui permet de bien comprendre comment ça marche (code court, visualisation de l'évolution possible).

-- Cécile.--
29/04/2009, 13h06
level_set

segmentation (level set)

Salut

Cela fait des jours que je tourne en rond et je n'arrive pas à modifier ce source : http://cimg.cvs.sourceforge.net/view..._levelsets.cpp pour appliquer le level set sur une image.

je voudrai savoir comment segmenter une image (je travail sur des images satellitaire) en utilisant ce source.

Merci.

IMPLEMENTATION DE LEVEL SET AVEC C++ WITH Qt

voila l'implementation de level set dans c++ et qt
modeles de mumford
le fichier "levelset.h"

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
 
#ifndef head_levelset
#define head_levelset
 
#include <QtGui/QImage>
#include <QtGui/QColor>
 
namespace levelset
{
	double H(double x);
	double Dirac(int H);
	double dphi(double **phi,int x,int y);
	double C(double **phi,int **intensite,int width,int height,int test);
	double lenght(double **phi, int width,int height);
	double area(double **phi, int width,int height);
	double div(double **phi,int x,int y);
	double F2(double **phi,int intensite,int x,int y,double c1,double c2,double landa1,double landa2,double mu,double nu);
	double F1(int intensite,double lenght,double area,double c1,double c2,double landa1,double landa2,double mu,double nu);
	double F3(double **phi,int x,int y,int u0);
	void tab_intensite(QImage *image,int **intensite);
	void level_zero(int x,int y,int ray,double **phi,int width,int height);
	void essai(QImage *image,int x,int y,int ray,int iter,int nbr,double landa1,double landa2,double mu,double nu,double dt);
	void evoluer_phi(double **phinew,double **phiold,int **intensite, int width,int height,int nbr,double landa1,double landa2,double mu,double nu,double dt);   
}
#endif

le fichier "level set.cpp"

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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
 
#include "levelset.h"
#include <stdio.h>
#include <math.h>
 
/*#define landa1	1
#define landa2	1
#define dt		1
#define mu		0.01*(255*255)
#define nu		0*/
#define ep		1
 
double max(double x,double y)
{
	if(x>y) return x;
	else return y;
}
double min(double x,double y)
{
	if(x>y) return y;
	else return x;
}
 
/////////////LA CONVERSTION VERS LE NIVEAU DE GRIS
void levelset::tab_intensite(QImage *image,int **intensite)
{
	for(int w=0;w<image->width();w++)
		for(int h=0;h<image->height();h++)
		{
			QColor coleur(image->pixel(w,h));
			intensite[w][h]=coleur.red();
		}
}
 
 
//////////////////////LE LEVEL  ZERO
void levelset::level_zero(int x,int y,int ray,double **phi,int width,int height)
{
	for(int w=0;w<width;w++)
		for(int h=0;h<height;h++)
		{
			float result=sqrt((float)((w-x)*(w-x)+(h-y)*(h-y)));
			if((int)(result)==ray)  phi[w][h]=0;
			else phi[w][h]=-(result-ray);			
		}
}
 
 
/////////////fonction HEAVSIDE
double levelset::H(double x)
{
	return 0.5*(1+(2/3.14 )*atan( x/ ep ) ); 
}
 
///////////////fonction DIRAC
double levelset::Dirac(int H)
{	
	return  ( 1/3.14 )*( ep/( ep*ep + H*H)); 
}
 
///////////////fonction DELTA  PHI
double levelset::dphi(double **phi,int x,int y)
{
	double inter;
	double part1=phi[x+1][y]-phi[x-1][y];
	double part2=phi[x][y+1]-phi[x][y-1];
	inter=0,5*0,5*(part1*part1+part2*part2);
	return sqrt(inter);
}
 
//////////////////////LENGHT
double levelset::lenght(double **phi, int width,int height)
{
	double len=0;
	for(int w=1;w<width-1;w++)
		for(int h=1;h<height-1;h++)
		{
			double max_phi = max(max(phi[w][h], phi[w + 1][h]), max(phi[w][h + 1], phi[w + 1][h + 1]));
			double min_phi = min(min(phi[w][h], phi[w + 1][h]), min(phi[w][h + 1], phi[w + 1][h + 1]));
			if(max_phi > 0.0 && min_phi < 0.0||phi[w][h]==0)  len+=Dirac(H(phi[w][h]))*dphi(phi,w,h);
		}
	return len;
}
 
 
 
////////////////SURFACE
double levelset::area(double **phi, int width,int height)
{
	double area=0;
	for(int w=1;w<width-1;w++)
		for(int h=1;h<height-1;h++)
		if(phi[w][h]>=0)  area+=H(phi[w][h]);
	return area;
}
 
///////////////fonction de CALCULE DE C1 et C2
double levelset::C(double **phi,int **intensite,int width,int height,int test)
{
	double nbr=0;
	double c=0;
	for(int w=0;w<width;w++)
		for(int h=0;h<height;h++)
		{
			if(phi[w][h]<=0 && test==1)			
			{
				c+=intensite[w][h]*(1-H(phi[w][h]));
				nbr+=(1-H(phi[w][h]));
			}
			else if(phi[w][h]>0 && test==2)
			{
				c+=intensite[w][h]*(H(phi[w][h]));
				nbr+=(H(phi[w][h]));
			}
		}
	return c/nbr;
}
 
 
 
///////////////////LA FONCTION DIV
double levelset::div(double **phi,int i,int j)
{
	double phi_x=(phi[i][j+1]-phi[i][j-1])/2.0;
	double phi_y=(phi[i+1][j]-phi[i-1][j])/2.0 ;
	double phi_xx=(phi[i][j+1]+phi[i][j-1]-2*phi[i][j])/1.0;
	double phi_yy=(phi[i+1][j]+phi[i-1][j]-2*phi[i][j])/1.0;
	double phi_xy =(phi[i+1][j+1]+phi[i-1][j-1]-phi[i-1][j+1]-phi[i+1][j-1])/4.0 ;
 
	double k=0;
	if(phi_x==0 && phi_y==0) k=0;
	else k=(phi_xx*phi_y*phi_y-2*phi_x*phi_y*phi_xy+phi_yy*phi_x*phi_x)/ pow((phi_x*phi_x+phi_y*phi_y),3.0/2.0);
 
	return k;
}
 
////////////////////LES TERME DES VITESSE
 
///////////////1
double levelset::F1(int u0,double lenght,double area,double c1,double c2,double landa1,double landa2,double mu,double nu)
{
	double part1=abs (u0-c1);
	double part2=abs (u0-c2);
	double inter= mu*lenght + nu*area + landa1*part1 - landa2*part2;
	return inter;
}
 
 
///////////////2
double levelset::F2(double **phi,int u0,int x,int y,double c1,double c2,double landa1,double landa2,double mu,double nu)
{
	double dirac = Dirac ( H(phi[x][y]) );
	double k = div(phi,x,y);
	double part1 = abs (u0-c1);
	double part2 = abs (u0-c2);
	double part3 = mu *  k  - nu + landa1 * part1*part1 - landa2 * part2*part2;
	double inter = dirac*part3;
	return inter;
}
 
/////////////3
double levelset::F3(double **phi,int x,int y,int u0)
{
	double k = div ( phi ,x ,y);
	double result = ( 1 + ( ep * k ) ) / ( 1 + u0 );
	return result;
}
 
/////////////////////LA FONCTION LEVEL SET 
void levelset::evoluer_phi(double **phinew,double **phiold,int **intensite, int width,int height,int nbr,double landa1,double landa2,double mu,double nu,double dt)
{
	int c1 = C(phiold, intensite, width, height, 1);
	int c2 = C(phiold, intensite, width, height, 2);
	double ar = area( phiold, width, height);
	double la = lenght( phiold, width, height);
 
	for(int w=1;w<width-1;w++)
		for(int h=1;h<height-1;h++)
		{
			if(nbr==1)		phinew[w][h]=phiold[w][h]+F1(intensite[w][h],la,ar,c1,c2,landa1,landa2,mu,nu)*dt;
			else if(nbr==2)	phinew[w][h]=phiold[w][h]+F2(phiold,intensite[w][h],w,h,c1,c2,landa1,landa2,mu,nu)*dt;
			//phinew[w][h] = phiold[w][h] + F3(phiold, w, h, intensite[w][h]) * dt * dphi(phiold,w,h);
		}
}
 
//////////////////ESSAI
void levelset::essai(QImage *image,int x,int y,int ray,int iter,int nbr,double landa1,double landa2,double mu,double nu,double dt)
{
	double cour[150];
 
	int wid=image->width();
	int hei=image->height();
 
	double **phiold=new double *[wid];
	for(int i=0;i<wid;i++)		phiold[i]=new double [hei];
 
	double **phinew=new double *[wid];
	for(int i=0;i<wid;i++)		phinew[i]=new double [hei];
 
	int **intensite=new int *[wid];
	for(int i=0;i<wid;i++)		intensite[i]=new int [hei];
 
	levelset::tab_intensite(image, intensite);
	level_zero(x, y, ray, phiold, wid, hei);
 
	int k=0;
	for(int i=0;i<iter;i++)
	{
		evoluer_phi(phinew,phiold,intensite,wid,hei,nbr,landa1,landa2,mu,nu,dt);
		for(int w=0;w<wid;w++)
			for(int h=0;h<hei;h++)
				phiold[w][h]=phinew[w][h];
 
		double somme=0;
		for(int w=0;w<wid;w++)
			for(int h=0;h<hei;h++)
				somme+=div(phiold,x,h);
		cour[k]=somme/wid*hei;
		k++;
	}
 
	//QImage *ima=new QImage("inter.bmp");
	FILE *out;
        out=fopen("rabah2.txt","w");
	for(int i=0;i<k;i++)
	{
		//image->setPixel(50+j,150+cour[j],qRgb(255,0,0));
		fprintf(out,"%lf\n",cour[i]);
	}
	fclose(out);	
 
	//ima->save("inter.bmp");
 
	for(int w=1;w<wid-1;w++)
		for(int h=1;h<hei-1;h++)
		{
			double max_phi=max(max(phiold[w][h],phiold[w+1][h]),max(phiold[w][h+1],phiold[w+1][h+1]));
			double min_phi=min(min(phiold[w][h],phiold[w+1][h]),min(phiold[w][h+1],phiold[w+1][h+1]));
 
                        if(max_phi>0.0&&min_phi<0.0||phiold[w][h]==0)   image->setPixel(w,h,qRgb(255,0,0));
		}
}

les fonction F1,F2,F3 sont les fonction de vitesse
pour F3 y'a une erreur.

ce programme est testé et il est correct......................

Citation:

Envoyé par bibouh123

voila l'implementation de level set dans c++ et qt
modeles de mumford
le fichier "levelset.h"

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
 
#ifndef head_levelset
#define head_levelset
 
#include <QtGui/QImage>
#include <QtGui/QColor>
 
namespace levelset
{
	double H(double x);
	double Dirac(int H);
	double dphi(double **phi,int x,int y);
	double C(double **phi,int **intensite,int width,int height,int test);
	double lenght(double **phi, int width,int height);
	double area(double **phi, int width,int height);
	double div(double **phi,int x,int y);
	double F2(double **phi,int intensite,int x,int y,double c1,double c2,double landa1,double landa2,double mu,double nu);
	double F1(int intensite,double lenght,double area,double c1,double c2,double landa1,double landa2,double mu,double nu);
	double F3(double **phi,int x,int y,int u0);
	void tab_intensite(QImage *image,int **intensite);
	void level_zero(int x,int y,int ray,double **phi,int width,int height);
	void essai(QImage *image,int x,int y,int ray,int iter,int nbr,double landa1,double landa2,double mu,double nu,double dt);
	void evoluer_phi(double **phinew,double **phiold,int **intensite, int width,int height,int nbr,double landa1,double landa2,double mu,double nu,double dt);   
}
#endif

le fichier "level set.cpp"

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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
 
#include "levelset.h"
#include <stdio.h>
#include <math.h>
 
/*#define landa1	1
#define landa2	1
#define dt		1
#define mu		0.01*(255*255)
#define nu		0*/
#define ep		1
 
double max(double x,double y)
{
	if(x>y) return x;
	else return y;
}
double min(double x,double y)
{
	if(x>y) return y;
	else return x;
}
 
/////////////LA CONVERSTION VERS LE NIVEAU DE GRIS
void levelset::tab_intensite(QImage *image,int **intensite)
{
	for(int w=0;w<image->width();w++)
		for(int h=0;h<image->height();h++)
		{
			QColor coleur(image->pixel(w,h));
			intensite[w][h]=coleur.red();
		}
}
 
 
//////////////////////LE LEVEL  ZERO
void levelset::level_zero(int x,int y,int ray,double **phi,int width,int height)
{
	for(int w=0;w<width;w++)
		for(int h=0;h<height;h++)
		{
			float result=sqrt((float)((w-x)*(w-x)+(h-y)*(h-y)));
			if((int)(result)==ray)  phi[w][h]=0;
			else phi[w][h]=-(result-ray);			
		}
}
 
 
/////////////fonction HEAVSIDE
double levelset::H(double x)
{
	return 0.5*(1+(2/3.14 )*atan( x/ ep ) ); 
}
 
///////////////fonction DIRAC
double levelset::Dirac(int H)
{	
	return  ( 1/3.14 )*( ep/( ep*ep + H*H)); 
}
 
///////////////fonction DELTA  PHI
double levelset::dphi(double **phi,int x,int y)
{
	double inter;
	double part1=phi[x+1][y]-phi[x-1][y];
	double part2=phi[x][y+1]-phi[x][y-1];
	inter=0,5*0,5*(part1*part1+part2*part2);
	return sqrt(inter);
}
 
//////////////////////LENGHT
double levelset::lenght(double **phi, int width,int height)
{
	double len=0;
	for(int w=1;w<width-1;w++)
		for(int h=1;h<height-1;h++)
		{
			double max_phi = max(max(phi[w][h], phi[w + 1][h]), max(phi[w][h + 1], phi[w + 1][h + 1]));
			double min_phi = min(min(phi[w][h], phi[w + 1][h]), min(phi[w][h + 1], phi[w + 1][h + 1]));
			if(max_phi > 0.0 && min_phi < 0.0||phi[w][h]==0)  len+=Dirac(H(phi[w][h]))*dphi(phi,w,h);
		}
	return len;
}
 
 
 
////////////////SURFACE
double levelset::area(double **phi, int width,int height)
{
	double area=0;
	for(int w=1;w<width-1;w++)
		for(int h=1;h<height-1;h++)
		if(phi[w][h]>=0)  area+=H(phi[w][h]);
	return area;
}
 
///////////////fonction de CALCULE DE C1 et C2
double levelset::C(double **phi,int **intensite,int width,int height,int test)
{
	double nbr=0;
	double c=0;
	for(int w=0;w<width;w++)
		for(int h=0;h<height;h++)
		{
			if(phi[w][h]<=0 && test==1)			
			{
				c+=intensite[w][h]*(1-H(phi[w][h]));
				nbr+=(1-H(phi[w][h]));
			}
			else if(phi[w][h]>0 && test==2)
			{
				c+=intensite[w][h]*(H(phi[w][h]));
				nbr+=(H(phi[w][h]));
			}
		}
	return c/nbr;
}
 
 
 
///////////////////LA FONCTION DIV
double levelset::div(double **phi,int i,int j)
{
	double phi_x=(phi[i][j+1]-phi[i][j-1])/2.0;
	double phi_y=(phi[i+1][j]-phi[i-1][j])/2.0 ;
	double phi_xx=(phi[i][j+1]+phi[i][j-1]-2*phi[i][j])/1.0;
	double phi_yy=(phi[i+1][j]+phi[i-1][j]-2*phi[i][j])/1.0;
	double phi_xy =(phi[i+1][j+1]+phi[i-1][j-1]-phi[i-1][j+1]-phi[i+1][j-1])/4.0 ;
 
	double k=0;
	if(phi_x==0 && phi_y==0) k=0;
	else k=(phi_xx*phi_y*phi_y-2*phi_x*phi_y*phi_xy+phi_yy*phi_x*phi_x)/ pow((phi_x*phi_x+phi_y*phi_y),3.0/2.0);
 
	return k;
}
 
////////////////////LES TERME DES VITESSE
 
///////////////1
double levelset::F1(int u0,double lenght,double area,double c1,double c2,double landa1,double landa2,double mu,double nu)
{
	double part1=abs (u0-c1);
	double part2=abs (u0-c2);
	double inter= mu*lenght + nu*area + landa1*part1 - landa2*part2;
	return inter;
}
 
 
///////////////2
double levelset::F2(double **phi,int u0,int x,int y,double c1,double c2,double landa1,double landa2,double mu,double nu)
{
	double dirac = Dirac ( H(phi[x][y]) );
	double k = div(phi,x,y);
	double part1 = abs (u0-c1);
	double part2 = abs (u0-c2);
	double part3 = mu *  k  - nu + landa1 * part1*part1 - landa2 * part2*part2;
	double inter = dirac*part3;
	return inter;
}
 
/////////////3
double levelset::F3(double **phi,int x,int y,int u0)
{
	double k = div ( phi ,x ,y);
	double result = ( 1 + ( ep * k ) ) / ( 1 + u0 );
	return result;
}
 
/////////////////////LA FONCTION LEVEL SET 
void levelset::evoluer_phi(double **phinew,double **phiold,int **intensite, int width,int height,int nbr,double landa1,double landa2,double mu,double nu,double dt)
{
	int c1 = C(phiold, intensite, width, height, 1);
	int c2 = C(phiold, intensite, width, height, 2);
	double ar = area( phiold, width, height);
	double la = lenght( phiold, width, height);
 
	for(int w=1;w<width-1;w++)
		for(int h=1;h<height-1;h++)
		{
			if(nbr==1)		phinew[w][h]=phiold[w][h]+F1(intensite[w][h],la,ar,c1,c2,landa1,landa2,mu,nu)*dt;
			else if(nbr==2)	phinew[w][h]=phiold[w][h]+F2(phiold,intensite[w][h],w,h,c1,c2,landa1,landa2,mu,nu)*dt;
			//phinew[w][h] = phiold[w][h] + F3(phiold, w, h, intensite[w][h]) * dt * dphi(phiold,w,h);
		}
}
 
//////////////////ESSAI
void levelset::essai(QImage *image,int x,int y,int ray,int iter,int nbr,double landa1,double landa2,double mu,double nu,double dt)
{
	double cour[150];
 
	int wid=image->width();
	int hei=image->height();
 
	double **phiold=new double *[wid];
	for(int i=0;i<wid;i++)		phiold[i]=new double [hei];
 
	double **phinew=new double *[wid];
	for(int i=0;i<wid;i++)		phinew[i]=new double [hei];
 
	int **intensite=new int *[wid];
	for(int i=0;i<wid;i++)		intensite[i]=new int [hei];
 
	levelset::tab_intensite(image, intensite);
	level_zero(x, y, ray, phiold, wid, hei);
 
	int k=0;
	for(int i=0;i<iter;i++)
	{
		evoluer_phi(phinew,phiold,intensite,wid,hei,nbr,landa1,landa2,mu,nu,dt);
		for(int w=0;w<wid;w++)
			for(int h=0;h<hei;h++)
				phiold[w][h]=phinew[w][h];
 
		double somme=0;
		for(int w=0;w<wid;w++)
			for(int h=0;h<hei;h++)
				somme+=div(phiold,x,h);
		cour[k]=somme/wid*hei;
		k++;
	}
 
	//QImage *ima=new QImage("inter.bmp");
	FILE *out;
        out=fopen("rabah2.txt","w");
	for(int i=0;i<k;i++)
	{
		//image->setPixel(50+j,150+cour[j],qRgb(255,0,0));
		fprintf(out,"%lf\n",cour[i]);
	}
	fclose(out);	
 
	//ima->save("inter.bmp");
 
	for(int w=1;w<wid-1;w++)
		for(int h=1;h<hei-1;h++)
		{
			double max_phi=max(max(phiold[w][h],phiold[w+1][h]),max(phiold[w][h+1],phiold[w+1][h+1]));
			double min_phi=min(min(phiold[w][h],phiold[w+1][h]),min(phiold[w][h+1],phiold[w+1][h+1]));
 
                        if(max_phi>0.0&&min_phi<0.0||phiold[w][h]==0)   image->setPixel(w,h,qRgb(255,0,0));
		}
}

les fonction F1,F2,F3 sont les fonction de vitesse
pour F3 y'a une erreur.

ce programme est testé et il est correct......................

Bonjour bibouh123,

J'espère que tu vas bien.

J'aimerai te poser quelques questions, notamment pour la méthode Level Set.

J'ai vu que t'as donné un code et que t'as mentionné qu'il était testé et approuvé.

Je l'ai bien entendu pris et essayer, mais en lisant le code j'ai trouvé des anomalies, comme:

-La matrice phi non initialisée (phiold) ?
-Les paramètres Ray, x, y de la fonction void levelset::level_zero(int x,int y,int ray,double **phi,int width,int height) représente quoi au juste? Il faut les donner en entrées, et on connait rien d'eux.

-Si on désire détecter le contour d'un objet dans une image qui contient d'autres objets, comment procéder?

Salutations.

24/03/2011, 11h47
pseudocode

Citation:

Envoyé par NightQuest

Je l'ai bien entendu pris et essayer, mais en lisant le code j'ai trouvé des anomalies, comme:

-La matrice phi non initialisée (phiold) ?

Initialisé dans la fonction levelset::essai()

Citation:

-Les paramètres Ray, x, y de la fonction void levelset::level_zero(int x,int y,int ray,double **phi,int width,int height) représente quoi au juste? Il faut les donner en entrées, et on connait rien d'eux.

Vu le code de la fonction, ca semble initialiser la fonction phi avec un cône.

x,y = le centre du cône
ray = le rayon du cône a l'altitude 0.

Citation:

-Si on désire détecter le contour d'un objet dans une image qui contient d'autres objets, comment procéder?

Et bien, c'est le principe du LevelSet de converger vers la meilleure segmentation possible (suivant les energies choisies). Le mieux étant tout de meme d'initialiser le zero de la fonction phi (-> choisir x,y,ray) pour être au plus proche de l'objet cherché.
14/04/2011, 11h02
bibouh123

Code:

-La matrice phi non initialisée (phiold) ?

ligne 193-194 dans le fichier .cpp (la fonction essai)

Code:

-Les paramètres Ray, x, y de la fonction void levelset::level_zero(int x,int y,int ray,double **phi,int width,int height) représente quoi au juste? Il faut les donner en entrées, et on connait rien d'eux.

la courbe initiale(un cercle) avec comme centre "p(x,y)" et un rayon "ray"

Code:

-Si on désire détecter le contour d'un objet dans une image qui contient d'autres objets, comment procéder?

la méthode level sets permet de détecter tous les objet dans l'image

et pour les appelles
1-"level_zero" la fonction initiale
2-"essai " c'est la fonction principale

Salutations.

Humm

Citation:

Envoyé par bibouh123

voila l'implementation de level set dans c++ et qt
modeles de mumford
le fichier "levelset.h"

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
 
#ifndef head_levelset
#define head_levelset
 
#include <QtGui/QImage>
#include <QtGui/QColor>
 
namespace levelset
{
	double H(double x);
	double Dirac(int H);
	double dphi(double **phi,int x,int y);
	double C(double **phi,int **intensite,int width,int height,int test);
	double lenght(double **phi, int width,int height);
	double area(double **phi, int width,int height);
	double div(double **phi,int x,int y);
	double F2(double **phi,int intensite,int x,int y,double c1,double c2,double landa1,double landa2,double mu,double nu);
	double F1(int intensite,double lenght,double area,double c1,double c2,double landa1,double landa2,double mu,double nu);
	double F3(double **phi,int x,int y,int u0);
	void tab_intensite(QImage *image,int **intensite);
	void level_zero(int x,int y,int ray,double **phi,int width,int height);
	void essai(QImage *image,int x,int y,int ray,int iter,int nbr,double landa1,double landa2,double mu,double nu,double dt);
	void evoluer_phi(double **phinew,double **phiold,int **intensite, int width,int height,int nbr,double landa1,double landa2,double mu,double nu,double dt);   
}
#endif

le fichier "level set.cpp"

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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
 
#include "levelset.h"
#include <stdio.h>
#include <math.h>
 
/*#define landa1	1
#define landa2	1
#define dt		1
#define mu		0.01*(255*255)
#define nu		0*/
#define ep		1
 
double max(double x,double y)
{
	if(x>y) return x;
	else return y;
}
double min(double x,double y)
{
	if(x>y) return y;
	else return x;
}
 
/////////////LA CONVERSTION VERS LE NIVEAU DE GRIS
void levelset::tab_intensite(QImage *image,int **intensite)
{
	for(int w=0;w<image->width();w++)
		for(int h=0;h<image->height();h++)
		{
			QColor coleur(image->pixel(w,h));
			intensite[w][h]=coleur.red();
		}
}
 
 
//////////////////////LE LEVEL  ZERO
void levelset::level_zero(int x,int y,int ray,double **phi,int width,int height)
{
	for(int w=0;w<width;w++)
		for(int h=0;h<height;h++)
		{
			float result=sqrt((float)((w-x)*(w-x)+(h-y)*(h-y)));
			if((int)(result)==ray)  phi[w][h]=0;
			else phi[w][h]=-(result-ray);			
		}
}
 
 
/////////////fonction HEAVSIDE
double levelset::H(double x)
{
	return 0.5*(1+(2/3.14 )*atan( x/ ep ) ); 
}
 
///////////////fonction DIRAC
double levelset::Dirac(int H)
{	
	return  ( 1/3.14 )*( ep/( ep*ep + H*H)); 
}
 
///////////////fonction DELTA  PHI
double levelset::dphi(double **phi,int x,int y)
{
	double inter;
	double part1=phi[x+1][y]-phi[x-1][y];
	double part2=phi[x][y+1]-phi[x][y-1];
	inter=0,5*0,5*(part1*part1+part2*part2);
	return sqrt(inter);
}
 
//////////////////////LENGHT
double levelset::lenght(double **phi, int width,int height)
{
	double len=0;
	for(int w=1;w<width-1;w++)
		for(int h=1;h<height-1;h++)
		{
			double max_phi = max(max(phi[w][h], phi[w + 1][h]), max(phi[w][h + 1], phi[w + 1][h + 1]));
			double min_phi = min(min(phi[w][h], phi[w + 1][h]), min(phi[w][h + 1], phi[w + 1][h + 1]));
			if(max_phi > 0.0 && min_phi < 0.0||phi[w][h]==0)  len+=Dirac(H(phi[w][h]))*dphi(phi,w,h);
		}
	return len;
}
 
 
 
////////////////SURFACE
double levelset::area(double **phi, int width,int height)
{
	double area=0;
	for(int w=1;w<width-1;w++)
		for(int h=1;h<height-1;h++)
		if(phi[w][h]>=0)  area+=H(phi[w][h]);
	return area;
}
 
///////////////fonction de CALCULE DE C1 et C2
double levelset::C(double **phi,int **intensite,int width,int height,int test)
{
	double nbr=0;
	double c=0;
	for(int w=0;w<width;w++)
		for(int h=0;h<height;h++)
		{
			if(phi[w][h]<=0 && test==1)			
			{
				c+=intensite[w][h]*(1-H(phi[w][h]));
				nbr+=(1-H(phi[w][h]));
			}
			else if(phi[w][h]>0 && test==2)
			{
				c+=intensite[w][h]*(H(phi[w][h]));
				nbr+=(H(phi[w][h]));
			}
		}
	return c/nbr;
}
 
 
 
///////////////////LA FONCTION DIV
double levelset::div(double **phi,int i,int j)
{
	double phi_x=(phi[i][j+1]-phi[i][j-1])/2.0;
	double phi_y=(phi[i+1][j]-phi[i-1][j])/2.0 ;
	double phi_xx=(phi[i][j+1]+phi[i][j-1]-2*phi[i][j])/1.0;
	double phi_yy=(phi[i+1][j]+phi[i-1][j]-2*phi[i][j])/1.0;
	double phi_xy =(phi[i+1][j+1]+phi[i-1][j-1]-phi[i-1][j+1]-phi[i+1][j-1])/4.0 ;
 
	double k=0;
	if(phi_x==0 && phi_y==0) k=0;
	else k=(phi_xx*phi_y*phi_y-2*phi_x*phi_y*phi_xy+phi_yy*phi_x*phi_x)/ pow((phi_x*phi_x+phi_y*phi_y),3.0/2.0);
 
	return k;
}
 
////////////////////LES TERME DES VITESSE
 
///////////////1
double levelset::F1(int u0,double lenght,double area,double c1,double c2,double landa1,double landa2,double mu,double nu)
{
	double part1=abs (u0-c1);
	double part2=abs (u0-c2);
	double inter= mu*lenght + nu*area + landa1*part1 - landa2*part2;
	return inter;
}
 
 
///////////////2
double levelset::F2(double **phi,int u0,int x,int y,double c1,double c2,double landa1,double landa2,double mu,double nu)
{
	double dirac = Dirac ( H(phi[x][y]) );
	double k = div(phi,x,y);
	double part1 = abs (u0-c1);
	double part2 = abs (u0-c2);
	double part3 = mu *  k  - nu + landa1 * part1*part1 - landa2 * part2*part2;
	double inter = dirac*part3;
	return inter;
}
 
/////////////3
double levelset::F3(double **phi,int x,int y,int u0)
{
	double k = div ( phi ,x ,y);
	double result = ( 1 + ( ep * k ) ) / ( 1 + u0 );
	return result;
}
 
/////////////////////LA FONCTION LEVEL SET 
void levelset::evoluer_phi(double **phinew,double **phiold,int **intensite, int width,int height,int nbr,double landa1,double landa2,double mu,double nu,double dt)
{
	int c1 = C(phiold, intensite, width, height, 1);
	int c2 = C(phiold, intensite, width, height, 2);
	double ar = area( phiold, width, height);
	double la = lenght( phiold, width, height);
 
	for(int w=1;w<width-1;w++)
		for(int h=1;h<height-1;h++)
		{
			if(nbr==1)		phinew[w][h]=phiold[w][h]+F1(intensite[w][h],la,ar,c1,c2,landa1,landa2,mu,nu)*dt;
			else if(nbr==2)	phinew[w][h]=phiold[w][h]+F2(phiold,intensite[w][h],w,h,c1,c2,landa1,landa2,mu,nu)*dt;
			//phinew[w][h] = phiold[w][h] + F3(phiold, w, h, intensite[w][h]) * dt * dphi(phiold,w,h);
		}
}
 
//////////////////ESSAI
void levelset::essai(QImage *image,int x,int y,int ray,int iter,int nbr,double landa1,double landa2,double mu,double nu,double dt)
{
	double cour[150];
 
	int wid=image->width();
	int hei=image->height();
 
	double **phiold=new double *[wid];
	for(int i=0;i<wid;i++)		phiold[i]=new double [hei];
 
	double **phinew=new double *[wid];
	for(int i=0;i<wid;i++)		phinew[i]=new double [hei];
 
	int **intensite=new int *[wid];
	for(int i=0;i<wid;i++)		intensite[i]=new int [hei];
 
	levelset::tab_intensite(image, intensite);
	level_zero(x, y, ray, phiold, wid, hei);
 
	int k=0;
	for(int i=0;i<iter;i++)
	{
		evoluer_phi(phinew,phiold,intensite,wid,hei,nbr,landa1,landa2,mu,nu,dt);
		for(int w=0;w<wid;w++)
			for(int h=0;h<hei;h++)
				phiold[w][h]=phinew[w][h];
 
		double somme=0;
		for(int w=0;w<wid;w++)
			for(int h=0;h<hei;h++)
				somme+=div(phiold,x,h);
		cour[k]=somme/wid*hei;
		k++;
	}
 
	//QImage *ima=new QImage("inter.bmp");
	FILE *out;
        out=fopen("rabah2.txt","w");
	for(int i=0;i<k;i++)
	{
		//image->setPixel(50+j,150+cour[j],qRgb(255,0,0));
		fprintf(out,"%lf\n",cour[i]);
	}
	fclose(out);	
 
	//ima->save("inter.bmp");
 
	for(int w=1;w<wid-1;w++)
		for(int h=1;h<hei-1;h++)
		{
			double max_phi=max(max(phiold[w][h],phiold[w+1][h]),max(phiold[w][h+1],phiold[w+1][h+1]));
			double min_phi=min(min(phiold[w][h],phiold[w+1][h]),min(phiold[w][h+1],phiold[w+1][h+1]));
 
                        if(max_phi>0.0&&min_phi<0.0||phiold[w][h]==0)   image->setPixel(w,h,qRgb(255,0,0));
		}
}

les fonction F1,F2,F3 sont les fonction de vitesse
pour F3 y'a une erreur.

ce programme est testé et il est correct......................

Est-ce que la procédure calculant la longueur ("lenght") ne serait pas incorrecte (je ne parle pas de la faute d'orthographe) ? Plus précisément,
ne serait-ce pas plutôt

Code:

1
2
 
len+=Dirac(phi[w][h])*dphi(phi,w,h);

et d'autre part ce choix de fonction "Dirac" est non-standard à ma connaissance.
A+

Citation:

Envoyé par bibouh123

voila l'implementation de level set dans c++ et qt
modeles de mumford
le fichier "levelset.h"

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
 
#ifndef head_levelset
#define head_levelset
 
#include <QtGui/QImage>
#include <QtGui/QColor>
 
namespace levelset
{
	double H(double x);
	double Dirac(int H);
	double dphi(double **phi,int x,int y);
	double C(double **phi,int **intensite,int width,int height,int test);
	double lenght(double **phi, int width,int height);
	double area(double **phi, int width,int height);
	double div(double **phi,int x,int y);
	double F2(double **phi,int intensite,int x,int y,double c1,double c2,double landa1,double landa2,double mu,double nu);
	double F1(int intensite,double lenght,double area,double c1,double c2,double landa1,double landa2,double mu,double nu);
	double F3(double **phi,int x,int y,int u0);
	void tab_intensite(QImage *image,int **intensite);
	void level_zero(int x,int y,int ray,double **phi,int width,int height);
	void essai(QImage *image,int x,int y,int ray,int iter,int nbr,double landa1,double landa2,double mu,double nu,double dt);
	void evoluer_phi(double **phinew,double **phiold,int **intensite, int width,int height,int nbr,double landa1,double landa2,double mu,double nu,double dt);   
}
#endif

le fichier "level set.cpp"

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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
 
#include "levelset.h"
#include <stdio.h>
#include <math.h>
 
/*#define landa1	1
#define landa2	1
#define dt		1
#define mu		0.01*(255*255)
#define nu		0*/
#define ep		1
 
double max(double x,double y)
{
	if(x>y) return x;
	else return y;
}
double min(double x,double y)
{
	if(x>y) return y;
	else return x;
}
 
/////////////LA CONVERSTION VERS LE NIVEAU DE GRIS
void levelset::tab_intensite(QImage *image,int **intensite)
{
	for(int w=0;w<image->width();w++)
		for(int h=0;h<image->height();h++)
		{
			QColor coleur(image->pixel(w,h));
			intensite[w][h]=coleur.red();
		}
}
 
 
//////////////////////LE LEVEL  ZERO
void levelset::level_zero(int x,int y,int ray,double **phi,int width,int height)
{
	for(int w=0;w<width;w++)
		for(int h=0;h<height;h++)
		{
			float result=sqrt((float)((w-x)*(w-x)+(h-y)*(h-y)));
			if((int)(result)==ray)  phi[w][h]=0;
			else phi[w][h]=-(result-ray);			
		}
}
 
 
/////////////fonction HEAVSIDE
double levelset::H(double x)
{
	return 0.5*(1+(2/3.14 )*atan( x/ ep ) ); 
}
 
///////////////fonction DIRAC
double levelset::Dirac(int H)
{	
	return  ( 1/3.14 )*( ep/( ep*ep + H*H)); 
}
 
///////////////fonction DELTA  PHI
double levelset::dphi(double **phi,int x,int y)
{
	double inter;
	double part1=phi[x+1][y]-phi[x-1][y];
	double part2=phi[x][y+1]-phi[x][y-1];
	inter=0,5*0,5*(part1*part1+part2*part2);
	return sqrt(inter);
}
 
//////////////////////LENGHT
double levelset::lenght(double **phi, int width,int height)
{
	double len=0;
	for(int w=1;w<width-1;w++)
		for(int h=1;h<height-1;h++)
		{
			double max_phi = max(max(phi[w][h], phi[w + 1][h]), max(phi[w][h + 1], phi[w + 1][h + 1]));
			double min_phi = min(min(phi[w][h], phi[w + 1][h]), min(phi[w][h + 1], phi[w + 1][h + 1]));
			if(max_phi > 0.0 && min_phi < 0.0||phi[w][h]==0)  len+=Dirac(H(phi[w][h]))*dphi(phi,w,h);
		}
	return len;
}
 
 
 
////////////////SURFACE
double levelset::area(double **phi, int width,int height)
{
	double area=0;
	for(int w=1;w<width-1;w++)
		for(int h=1;h<height-1;h++)
		if(phi[w][h]>=0)  area+=H(phi[w][h]);
	return area;
}
 
///////////////fonction de CALCULE DE C1 et C2
double levelset::C(double **phi,int **intensite,int width,int height,int test)
{
	double nbr=0;
	double c=0;
	for(int w=0;w<width;w++)
		for(int h=0;h<height;h++)
		{
			if(phi[w][h]<=0 && test==1)			
			{
				c+=intensite[w][h]*(1-H(phi[w][h]));
				nbr+=(1-H(phi[w][h]));
			}
			else if(phi[w][h]>0 && test==2)
			{
				c+=intensite[w][h]*(H(phi[w][h]));
				nbr+=(H(phi[w][h]));
			}
		}
	return c/nbr;
}
 
 
 
///////////////////LA FONCTION DIV
double levelset::div(double **phi,int i,int j)
{
	double phi_x=(phi[i][j+1]-phi[i][j-1])/2.0;
	double phi_y=(phi[i+1][j]-phi[i-1][j])/2.0 ;
	double phi_xx=(phi[i][j+1]+phi[i][j-1]-2*phi[i][j])/1.0;
	double phi_yy=(phi[i+1][j]+phi[i-1][j]-2*phi[i][j])/1.0;
	double phi_xy =(phi[i+1][j+1]+phi[i-1][j-1]-phi[i-1][j+1]-phi[i+1][j-1])/4.0 ;
 
	double k=0;
	if(phi_x==0 && phi_y==0) k=0;
	else k=(phi_xx*phi_y*phi_y-2*phi_x*phi_y*phi_xy+phi_yy*phi_x*phi_x)/ pow((phi_x*phi_x+phi_y*phi_y),3.0/2.0);
 
	return k;
}
 
////////////////////LES TERME DES VITESSE
 
///////////////1
double levelset::F1(int u0,double lenght,double area,double c1,double c2,double landa1,double landa2,double mu,double nu)
{
	double part1=abs (u0-c1);
	double part2=abs (u0-c2);
	double inter= mu*lenght + nu*area + landa1*part1 - landa2*part2;
	return inter;
}
 
 
///////////////2
double levelset::F2(double **phi,int u0,int x,int y,double c1,double c2,double landa1,double landa2,double mu,double nu)
{
	double dirac = Dirac ( H(phi[x][y]) );
	double k = div(phi,x,y);
	double part1 = abs (u0-c1);
	double part2 = abs (u0-c2);
	double part3 = mu *  k  - nu + landa1 * part1*part1 - landa2 * part2*part2;
	double inter = dirac*part3;
	return inter;
}
 
/////////////3
double levelset::F3(double **phi,int x,int y,int u0)
{
	double k = div ( phi ,x ,y);
	double result = ( 1 + ( ep * k ) ) / ( 1 + u0 );
	return result;
}
 
/////////////////////LA FONCTION LEVEL SET 
void levelset::evoluer_phi(double **phinew,double **phiold,int **intensite, int width,int height,int nbr,double landa1,double landa2,double mu,double nu,double dt)
{
	int c1 = C(phiold, intensite, width, height, 1);
	int c2 = C(phiold, intensite, width, height, 2);
	double ar = area( phiold, width, height);
	double la = lenght( phiold, width, height);
 
	for(int w=1;w<width-1;w++)
		for(int h=1;h<height-1;h++)
		{
			if(nbr==1)		phinew[w][h]=phiold[w][h]+F1(intensite[w][h],la,ar,c1,c2,landa1,landa2,mu,nu)*dt;
			else if(nbr==2)	phinew[w][h]=phiold[w][h]+F2(phiold,intensite[w][h],w,h,c1,c2,landa1,landa2,mu,nu)*dt;
			//phinew[w][h] = phiold[w][h] + F3(phiold, w, h, intensite[w][h]) * dt * dphi(phiold,w,h);
		}
}
 
//////////////////ESSAI
void levelset::essai(QImage *image,int x,int y,int ray,int iter,int nbr,double landa1,double landa2,double mu,double nu,double dt)
{
	double cour[150];
 
	int wid=image->width();
	int hei=image->height();
 
	double **phiold=new double *[wid];
	for(int i=0;i<wid;i++)		phiold[i]=new double [hei];
 
	double **phinew=new double *[wid];
	for(int i=0;i<wid;i++)		phinew[i]=new double [hei];
 
	int **intensite=new int *[wid];
	for(int i=0;i<wid;i++)		intensite[i]=new int [hei];
 
	levelset::tab_intensite(image, intensite);
	level_zero(x, y, ray, phiold, wid, hei);
 
	int k=0;
	for(int i=0;i<iter;i++)
	{
		evoluer_phi(phinew,phiold,intensite,wid,hei,nbr,landa1,landa2,mu,nu,dt);
		for(int w=0;w<wid;w++)
			for(int h=0;h<hei;h++)
				phiold[w][h]=phinew[w][h];
 
		double somme=0;
		for(int w=0;w<wid;w++)
			for(int h=0;h<hei;h++)
				somme+=div(phiold,x,h);
		cour[k]=somme/wid*hei;
		k++;
	}
 
	//QImage *ima=new QImage("inter.bmp");
	FILE *out;
        out=fopen("rabah2.txt","w");
	for(int i=0;i<k;i++)
	{
		//image->setPixel(50+j,150+cour[j],qRgb(255,0,0));
		fprintf(out,"%lf\n",cour[i]);
	}
	fclose(out);	
 
	//ima->save("inter.bmp");
 
	for(int w=1;w<wid-1;w++)
		for(int h=1;h<hei-1;h++)
		{
			double max_phi=max(max(phiold[w][h],phiold[w+1][h]),max(phiold[w][h+1],phiold[w+1][h+1]));
			double min_phi=min(min(phiold[w][h],phiold[w+1][h]),min(phiold[w][h+1],phiold[w+1][h+1]));
 
                        if(max_phi>0.0&&min_phi<0.0||phiold[w][h]==0)   image->setPixel(w,h,qRgb(255,0,0));
		}
}

les fonction F1,F2,F3 sont les fonction de vitesse
pour F3 y'a une erreur.

ce programme est testé et il est correct......................

bonjour,
j'aimerai bien savoir pourquoi quand je le tape ce programme ,il m'affiche beaucoup d'erreur

void levelset::tab_intensite(QImage *image ,int **intensite)

for(int w=0;w<wid;w++) déclaration attendue

wid*hei; identification wid et hei non définis
out=fopen("rabah2.txt","w"); cette déclaration out n'a pas de classe de stockage

merci d'avance

01/02/2012, 09h36
leo5th

bonjour,

Je ne comprends pas pourquoi tu mets la valeur de variable inter comme ca?
(Dans la fonction de calcul delta de phi)

inter=0,5*0,5*(part1*part1+part2*part2);

La valeur de inter est toujours à 0. Ce n'est pas:

inter=0.5*0.5*(part1*part1+part2*part2);

J'ai essayé avec 0.5. Mais le résultat n'est pas bon.