DSPRelated.com
Forums

FFTW multiplication, marr hildreth filter, and normalization

Started by mfactor July 4, 2007
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!




Nevermind, found the error.