Hello everybody,
I am trying to filter a image and do some edge detection using the Marr
Hildreth method. I have (hopefully!) two problems:
1) I read the image and stores it in a array. Then, I create a laplacian
of a gaussian using the same size of the imagem. I transform both using
FFTW r2c_2d, do the pointwise complex multiplication, and then transform
back using c2r_2d. On the returned array I do the normalization, dividing
the transformed array by width*height. But the returned array is never
normalized: depending the sigma of the gaussian, the returned array can be
very large or very small. So here is the first problem: where, and how can
I normalize the returned array? I tried to use also 1/sqrt(width*height),
but with no success. I think that maybe this is related to the
normalization of the laplacian of the gaussian, because the original
signal (the image) is on the intensity scale, 0 for white and 1 for black.
I saw some posts here at DSPrelated but with no luck.
2) I don't know if it is related to the above problem, but the laplacian
of a gaussian, that should return positive and negative values, is
returning only positive or negative values for several sigmas, and return
positive and negative ones with just one sigma, 0.4. I think that this is
odd, because the sigma should control the smoothing and the scale of the
transform...
I have written the gaussian in the four corners of the image because the
periodical characteristic of the fft of a gaussian.
Here is some code:
**** Initializations goes hre
float sigma = 0.42;
for(i=0; i < altura; i++)
{
for(j=0; j < largura; j++)
{
in_sinal[i*largura + j] = (double)sinal[i*largura + j];
in_gauss[i*largura + j]
= -(1/sqrt(largura*altura))*
((1 - ( (i-altura)*(i-altura) + (j)*(j) )/(sigma*sigma))
*exp(-( (i-altura)*(i-altura) + (j*j) ) /(2*sigma*sigma))
+ (1-( (i)*(i) + (j-largura)*(j-largura) )/(sigma*sigma))
*exp(-( i*i + (j-largura)*(j-largura) ) /(2*sigma*sigma))
+ (1-( (i-altura)*(i-altura) + (j-largura)*(j-largura)
)/(sigma*sigma))
*exp(-( (i-altura)*(i-altura) + (j-largura)*(j-largura) )
/(2*sigma*sigma))
+ (1-( i*i + j*j )/(sigma*sigma))*exp(-( (i*i + j*j))
/(2*sigma*sigma)));
}
}
p = fftw_plan_dft_r2c_2d(altura, largura,in_sinal,
out_sinal,FFTW_ESTIMATE);
fftw_execute(p);
p_gauss = fftw_plan_dft_r2c_2d(altura, largura,in_gauss,
out_gauss,FFTW_ESTIMATE);
fftw_execute(p_gauss);
for(l=0; l < altura/2 +1; l++)
{
for(m=0; m < largura; m++)
{
out_sinal[l*largura + m] =
((creal(out_sinal[l*largura + m])*creal(out_gauss[l*largura + m])
- cimag(out_sinal[l*largura + m])*cimag(out_gauss[l*largura+m]))
+ I*(creal(out_sinal[l*largura + m])*cimag(out_gauss[l*largura + m])
+ cimag(out_sinal[l*largura + m])*creal(out_gauss[l*largura
+m])));
}
}
p = fftw_plan_dft_c2r_2d(altura, largura, out_sinal,
in_sinal,FFTW_ESTIMATE);
fftw_execute(p);
return in_sinal;
*** end of it
Thank you so much!