У меня вопрос по БПФ. Мне уже удается делать БПФ вперед и назад, используя FFTW в C. Теперь я хочу применить фильтр верхних частот для обнаружения краев, некоторые из моих источников сказали, что просто обнуляют центр величины.
Это мое входное изображение http://i62.tinypic.com/2wnxvfl.jpg
В основном, что я делаю:
- Вперед БПФ
- Преобразование вывода в 2D-массив
- Выполнить сдвиг вперед БПФ
- Сделайте значение real и imag равным 0, когда расстояние от центра составляет 25% от высоты.
- Генерация величины
- Выполнить обратное смещение БПФ
- Преобразовать в одномерный массив
- Выполните обратное БПФ.
Это исходная величина, обработанная величина и результат.
http://i58.tinypic.com/aysx9s.png
может кто-нибудь помочь мне, чтобы сказать мне, какая часть неверна и как сделать фильтрацию верхних частот с помощью FFTW в C.
Спасибо.
Исходный код:
unsigned char **FFT2(int width,int height, unsigned char **pixel, char line1[100],char line2[100], char line3[100],char filename[100])
{
fftw_complex* in, * dft, * idft, * dft2;
//fftw_complex tmp1,tmp2;
fftw_plan plan_f,plan_i;
int i,j,k,w,h,N,w2,h2;
w = width;
h = height;
N = w*h;
unsigned char **pixel_out;
pixel_out = malloc(h*sizeof(unsigned char*));
for(i = 0 ; i<h;i++)
pixel_out[i]=malloc(w*sizeof(unsigned char));
in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N);
dft = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N);
dft2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N);
idft = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N);
/*run forward FFT*/
plan_f = fftw_plan_dft_2d(w,h,in,dft,FFTW_FORWARD,FFTW_ESTIMATE);
for(i = 0,k = 0 ; i < h ; i++)
{
for(j = 0 ; j < w ; j++,k++)
{
in[k][0] = pixel[i][j];
in[k][1] = 0.0;
}
}
fftw_execute(plan_f);
double maxReal = 0.0;
for(i = 0 ; i < N ; i++)
maxReal = dft[i][0] > maxReal ? dft[i][0] : maxReal;
printf("MAX REAL : %f\n",maxReal);
/*fftshift*/
//convert to 2d
double ***temp1;
temp1 = malloc(h * sizeof (double**));
for (i = 0;i < h; i++){
temp1[i] = malloc(w*sizeof (double*));
for (j = 0; j < w; j++){
temp1[i][j] = malloc(2*sizeof(double));
}
}
double ***temp2;
temp2 = malloc(h * sizeof (double**));
for (i = 0;i < h; i++){
temp2[i] = malloc(w*sizeof (double*));
for (j = 0; j < w; j++){
temp2[i][j] = malloc(2*sizeof(double));
}
}
for (i = 0;i < h; i++){
for (j = 0; j < w; j++){
temp1[i][j][0] = dft[i*w+j][0];
temp1[i][j][1] = dft[i*w+j][1];
}
}
int m2 = h/2;
int n2 = w/2;
//forward shifting
for (i = 0; i < m2; i++)
{
for (k = 0; k < n2; k++)
{
double tmp13[2] = {temp1[i][k][0],temp1[i][k][1]};
temp1[i][k][0] = temp1[i+m2][k+n2][0];
temp1[i][k][1] = temp1[i+m2][k+n2][1];
temp1[i+m2][k+n2][0] = tmp13[0];
temp1[i+m2][k+n2][1] = tmp13[1];
double tmp24[2] = {temp1[i+m2][k][0],temp1[i+m2][k][1]};
temp1[i+m2][k][0] = temp1[i][k+n2][0];
temp1[i+m2][k][1] = temp1[i][k+n2][1];
temp1[i][k+n2][0] = tmp24[0];
temp1[i][k+n2][1] = tmp24[1];
}
}
//process
for (i = 0;i < h; i++){
for (j = 0; j < w; j++){
if(distance_to_center(i,j,m2,n2) < 0.25*h)
{
temp1[i][j][0] = (double)0.0;
temp1[i][j][1] = (double)0.0;
}
}
}
/* copy for magnitude */
for (i = 0;i < h; i++){
for (j = 0; j < w; j++){
temp2[i][j][0] = temp1[i][j][0];
temp2[i][j][1] = temp1[i][j][1];
}
}
//backward shifting
for (i = 0; i < m2; i++)
{
for (k = 0; k < n2; k++)
{
double tmp13[2] = {temp1[i][k][0],temp1[i][k][1]};
temp1[i][k][0] = temp1[i+m2][k+n2][0];
temp1[i][k][1] = temp1[i+m2][k+n2][1];
temp1[i+m2][k+n2][0] = tmp13[0];
temp1[i+m2][k+n2][1] = tmp13[1];
double tmp24[2] = {temp1[i+m2][k][0],temp1[i+m2][k][1]};
temp1[i+m2][k][0] = temp1[i][k+n2][0];
temp1[i+m2][k][1] = temp1[i][k+n2][1];
temp1[i][k+n2][0] = tmp24[0];
temp1[i][k+n2][1] = tmp24[1];
}
}
//convert back to 1d
for (i = 0;i < h; i++){
for (j = 0; j < w; j++){
dft[i*w+j][0] = temp1[i][j][0];
dft[i*w+j][1] = temp1[i][j][1];
dft2[i*w+j][0] = temp2[i][j][0];
dft2[i*w+j][1] = temp2[i][j][1];
}
}
/* magnitude */
double max = 0;
double min = 0;
double mag=0;
for (i = 0, k = 1; i < h; i++){
for (j = 0; j < w; j++, k++){
mag = sqrt(pow(dft2[i*w+j][0],2) + pow(dft2[i*w+j][1],2));
if (max < mag)
max = mag;
}
}
double **magTemp;
magTemp = malloc(h * sizeof (double*));
for (i = 0;i < h; i++){
magTemp[i] = malloc(w*sizeof (double));
}
for(i = 0,k = 0 ; i < h ; i++)
{
for(j = 0 ; j < w ; j++,k++)
{
double mag = sqrt(pow(dft2[i*w+j][0],2) + pow(dft2[i*w+j][1],2));
mag = 255*(mag/max);
//magTemp[i][j] = 255-mag; //Putih
magTemp[i][j] = mag; //Item
}
}
/* brightening magnitude*/
for(i = 0,k = 0 ; i < h ; i++)
{
for(j = 0 ; j < w ; j++,k++)
{
//double temp = magTemp[i][j];
double temp = (double)(255/(log(1+255)))*log(1+magTemp[i][j]);
pixel_out[i][j] = (unsigned char)temp;
}
}
generateImage(width,height,pixel_out,line1,line2,line3,filename,"magnitude");
/* backward fft */
plan_i = fftw_plan_dft_2d(w,h,dft,idft,FFTW_BACKWARD,FFTW_ESTIMATE);
fftw_execute(plan_i);
for(i = 0,k = 0 ; i < h ; i++)
{
for(j = 0 ; j < w ; j++,k++)
{
double temp = idft[i*w+j][0]/N;
pixel_out[i][j] = (unsigned char)temp; //+ pixel[i][j];
}
}
generateImage(width,height,pixel_out,line1,line2,line3,filename,"backward");
return pixel_out;
}
ИЗМЕНИТЬ новый исходный код
Я добавляю эту часть перед переключением вперед, результат также соответствует ожиданиям.
//proses
//create filter
unsigned char **pixel_filter;
pixel_filter = malloc(h*sizeof(unsigned char*));
for(i = 0 ; i<h;i++)
pixel_filter[i]=malloc(w*sizeof(unsigned char));
for (i = 0;i < h; i++){
for (j = 0; j < w; j++){
if(distance_to_center(i,j,m2,n2) < 20)
{
pixel_filter[i][j] = 0;
}
else
{
pixel_filter[i][j] = 255;
}
}
}
generateImage(width,height,pixel_filter,line1,line2,line3,filename,"filter1");
for (i = 0; i < m2; i++)
{
for (k = 0; k < n2; k++)
{
unsigned char tmp13 = pixel_filter[i][k];
pixel_filter[i][k] = pixel_filter[i+m2][k+n2];
pixel_filter[i+m2][k+n2] = tmp13;
unsigned char tmp24 = pixel_filter[i+m2][k];
pixel_filter[i+m2][k] = pixel_filter[i][k+n2];
pixel_filter[i][k+n2] = tmp24;
}
}
generateImage(width,height,pixel_filter,line1,line2,line3,filename,"filter2");
for (i = 0;i < h; i++){
for (j = 0; j < w; j++){
temp1[i][j][0] *= pixel_filter[i][j];
temp1[i][j][1] *= pixel_filter[i][j];
}
}
