Bonjour,
je cherche à utiliser la bibliotheque fftw pour faire une fft sur une image et passer un filtre gaussien.
J'ai utilisé le tuto de millie mais ca ne marche pas. Il y a un truc que je n'ai pas du comprendre.
Voici mon code :
Les fonctions utiles :
Code C : 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
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
 
void FFT(const float* imentree, float* reimsortie, float* imimsortie,unsigned int largeur, unsigned int hauteur)
{
	fftw_complex* spatial_repr;
	fftw_complex* frequency_repr;
	unsigned int i;
	unsigned int j;
	fftw_plan plan;
	int x,y;
 
	spatial_repr= (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*largeur*hauteur);
	frequency_repr= (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*largeur*hauteur);
 
 
	for(i=0;i<largeur*hauteur;i++)
	{
	   spatial_repr[i][0] = imentree[i];
	   spatial_repr[i][1] =  0.0f;
	}
 
	 /*on calcule le plan d'exécution*/
	 plan=fftw_plan_dft_2d(hauteur, largeur, spatial_repr, frequency_repr, FFTW_FORWARD, FFTW_ESTIMATE);
 
	 /*on calcule la transformée*/
	 fftw_execute(plan);
 
	for(j=0;j<hauteur;j++)
	{
		for(i=0;i<largeur;i++)
		{
	        /*on recentre l'image*/
			x=i;
			y=j;
			if (i<largeur/2 && j<hauteur/2)
			{ 
				x=i+largeur/2; 
				y=j+hauteur/2; 
			}
			if (i>=largeur/2 && j<hauteur/2)
			{
				x=i-largeur/2;
				y=j+hauteur/2;
			}
			if (i<largeur/2 && j>=hauteur/2)
			{ 
				x=i+largeur/2; 
				y=j-hauteur/2; 
			}
			if (i>=largeur/2 && j>=hauteur/2)
			{ 
				x=i-largeur/2; 
				y=j-hauteur/2; 
			}
			reimsortie[y*largeur+x]=frequency_repr[j*largeur+i][0];
			imimsortie[y*largeur+x]=frequency_repr[j*largeur+i][1];
		}
	}
 
	fftw_destroy_plan(plan);
	fftw_free(spatial_repr);
	fftw_free(frequency_repr);
 
 
}
void FFTInverse( const float* reIn, const float* imIn, float* out, unsigned int largeur, unsigned int hauteur)
{
	/* reIn : partie réel de l'image dans l'espace de Fourier
 * imIn : partie imaginaire de l'image dans l'espace de Fourier
 * out : image de sortie
 * largeur : largeur des images d'entrée et de sortie
 */
 
   fftw_complex* spatial_repr;
   fftw_complex* frequency_repr;
   unsigned int i;
   unsigned int j;
   int x,y;
   fftw_plan plan;
 
   spatial_repr=(fftw_complex*) fftw_malloc(sizeof(fftw_complex)*largeur*hauteur);
   frequency_repr= (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*largeur*hauteur);
 
   for(j=0;j<hauteur;j++)
   {
		for(i=0;i<largeur;i++)
		{
			/*on décentre*/
			x=i;
			y=j;
			if(i<largeur/2 && j<hauteur/2)
			{ 
				x=i+largeur/2;
				y=j+hauteur/2;
			}
			if(i>=largeur/2 && j<hauteur/2)
			{ 
				x=i-largeur/2; 
				y=j+hauteur/2; 
			}
			if(i<largeur/2 && j>=hauteur/2)
			{ 
				x=i+largeur/2; 
				y=j-hauteur/2; 
			}
			if(i>=largeur/2 && j>=hauteur/2)
			{ 
				x=i-largeur/2; 
				y=j-hauteur/2; 
			}
			frequency_repr[j*largeur+i][0]=reIn[y*largeur+x];
			frequency_repr[j*largeur+i][1]=imIn[y*largeur+x];
		}
	}
	plan=fftw_plan_dft_2d(hauteur, largeur, frequency_repr, spatial_repr, FFTW_BACKWARD, FFTW_ESTIMATE);
 
	fftw_execute(plan);
 
	 /*on retranscrit l'image complexe en image réelle, sans oublier de diviser par largeur*hauteur*/
	for(i=0;i<largeur*hauteur;i++)
	{
		out[i]=spatial_repr[i][0]/(largeur*hauteur);
	}
 
	fftw_destroy_plan(plan);
	fftw_free(spatial_repr);
	fftw_free(frequency_repr);
}
 
void filtreAppliquer(float* reTabImage,float* imTabImage,const float* reTabFiltre, const float* imTabFiltre,unsigned int largeur,unsigned int hauteur)
{
 
	float a,b,c,d;
	unsigned int i;
 
	for(i=0; i< largeur * hauteur; i++)
	{
		a = reTabImage[i];
		b = imTabImage[i];
		c = reTabFiltre[i];
		d = imTabFiltre[i];
 
    /*calcul du produit des complexes*/
		reTabImage[i] = a*c - b*d;
		imTabImage[i] = b*c + a*d;
 
	}
 
}
 
void filtreFlouGaussien(float* reTabImage,float* imTabImage,unsigned int largeur,unsigned int hauteur,float sigma) /*sigma du flou*/
{
	unsigned int i;
	unsigned int j;
 
	float* filtreR = (float*)malloc(sizeof(float) * largeur*hauteur);
	float* filtreI = (float*)malloc(sizeof(float) * largeur*hauteur);
 
	if((filtreR==NULL)||(filtreI ==NULL))
	{
		fprintf(stderr, "Probleme d'allocation memoire");
		exit(EXIT_FAILURE);
	}
 
	for(j=0; j<hauteur;j++)
	{
        for(i=0; i<largeur;i++)
		{
		/*on applique la formule du flou gaussien*/
			filtreR[i+largeur*j]=exp(-carre(sigma)*(carre(i-largeur/2.0)+carre(j-hauteur/2.0))/(2.0));
			filtreI[i+largeur*j]=0.0f;
    /*on a pensé à recentrer en 0, car le point (0,0) correspond au point (largeur/2, hauteur/2)*/
		}
	}
	filtreAppliquer(reTabImage, imTabImage, filtreR, filtreI, largeur, hauteur);
 
	free(filtreR);
	free(filtreI);
}
double carre(double a)
{
	return (a*a);
}
void ExtractionRGB(Image &im1,float *R,float *G,float *B)
{
	int k=0;
	for(int j=0;j<im1.hauteur;j++)
	{
		for(int i=0;i<im1.largeur*3;i+=3)
		{
			R[k]=im1.ptrLigne[j][i];
			G[k]=im1.ptrLigne[j][i+1];
			B[k]=im1.ptrLigne[j][i+2];
			k++;
		}
	}
}
void FusionRGB(float *R,float *G,float *B,Image &im1)
{
	int k=0;
	for(int j=0;j<im1.hauteur;j++)
	{
		for(int i=0;i<im1.largeur*3;i+=3)
		{
			im1.ptrLigne[j][i]=R[k];
			im1.ptrLigne[j][i+1]=G[k];
			im1.ptrLigne[j][i+2]=B[k];
			k++;
		}
	}
}

La fonction principale où im1 est mon image:
Code C : 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
 
float *R,*G,*B;
float *resortieR,*resortieG,*resortieB;
float *imsortieR,*imsortieG,*imsortieB;
int n=im1.largeur*im1.hauteur;
R=new float[n];
G=new float[n];
B=new float[n];
 
	resortieR=new float[n];
	imsortieR=new float[n];
	resortieG=new float[n];
	imsortieG=new float[n];
	resortieB=new float[n];
	imsortieB=new float[n];
 
	ExtractionRGB(im1,R,G,B);
	//FFT sur les composantes
	FFT(R,resortieR,imsortieR,im1.largeur,im1.hauteur);
	FFT(G,resortieG,imsortieG,im1.largeur,im1.hauteur);
	FFT(B,resortieB,imsortieB,im1.largeur,im1.hauteur);
	//Filtrage passe haut du tableau. 
	/*filtreFlouGaussien(resortieS,imsortieS,im1.largeur,im1.hauteur,0.1);
	filtreFlouGaussien(resortieL,imsortieL,im1.largeur,im1.hauteur,0.1);*/
	//FFT inverse des composantes 
	FFTInverse(resortieR,imsortieR,R,im1.largeur,im1.hauteur);
	FFTInverse(resortieG,imsortieG,G,im1.largeur,im1.hauteur);
	FFTInverse(resortieB,imsortieB,B,im1.largeur,im1.hauteur);
	//FusionRGB
	FusionRGB(R,G,B,im1);
Si je retire le flou gaussien, je devrais avoir comme image finale, la meme que l'image initiale, non???
Et je n'ai pas la meme.
Si j'enleve la FFT et la FFTInverse, ca marche donc mes fonctions ExtractionRGB et FusionRGB marchent bien.
Ca semble etre mon utilisation de FFT et FFTInverse qui posent probleme....

Si quelqu'un voit mon erreur......