void FFT(vector<vector<Intensity> > &in,fftw_complex *fftin1,fftw_complex *fftout1)
{
int x,y;
double t0,k1,k2,k3;
fftw_plan p1;
for(y = 0;y<bmp_w;y++)
{
for(x = 0;x<bmp_h;x++)
{
t0 = (double)in[y][x].Y;
fftin1[y*bmp_h+x][0] = t0;//*pow(-1,(double)(x*y));
fftin1[y*bmp_h+x][1] = 0.0;
}
}
p1 = fftw_plan_dft_2d(bmp_h,bmp_w,fftin1,fftout1,FFTW_FORWARD,FFTW_ESTIMATE);
fftw_execute(p1);
for(y=0;y<bmp_w;y++)
{
for(x=0;x<bmp_h;x++)
{
d1[y*bmp_h+x] = ((double)sqrt((double)(y*y)+(double)(x*x)));
k1 = exp((-c*pow(d1[y*bmp_h+x],2))/(d0*d0));
H[y*bmp_h+x] = (rh-rl) * (1-k1) + rl;
fftout1[y*bmp_h+x][0] *= H[y*bmp_h+x];
fftout1[y*bmp_h+x][1] *= H[y*bmp_h+x];
}
}
}