我有一个关于FFT的问题。 我已经设法在C中使用FFTW向前和向后进行FFT。现在,我想应用高通滤波器进行边缘检测,我的一些消息来源说只是将幅度的中心归零。
这是我的输入图片https://sofzh.miximages.com/c/2wnxvfl.jpg
基本上我所做的是:
这是原始幅度,处理幅度和结果
https://sofzh.miximages.com/c/aysx9s.png
有人可以帮助我,告诉我哪个部分是错的,以及如何在C中使用FFTW进行高通滤波
谢谢。
源代码:
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 maxReal ? dft[i][0] : maxReal; printf("MAX REAL : %fn",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]; } }
你的总体想法是好的。 从输出中,很难判断程序中是否存在简单的会计问题,或者这是否是预期的结果。 尝试使用更多空白空间填充源图像,并过滤出频域中较小的区域。
作为旁注,在C中执行此操作显得非常痛苦。 这是Matlab中的等效实现。 不包括绘图,它是大约10行代码。 您也可以尝试使用Numerical Python(NumPy)。
% Demonstrate frequency-domain image filtering in Matlab % Define the grid x = linspace(-1, 1, 1001); y = x; [X, Y] = meshgrid(x, y); % Make a square (source image) rect = (abs(X) < 0.1) & (abs(Y) < 0.1); % Compute the transform rect_hat = fft2(rect); % Make the high-pass filter R = sqrt(X.^2 + Y.^2); filt = (R > 0.05); % Apply the filter rect_hat_filtered = rect_hat .* ifftshift(filt); % Compute the inverse transform rect_filtered = ifft2(rect_hat_filtered); %% Plot everything figure(1) imagesc(rect); title('source'); axis square saveas(gcf, 'fig1.png'); figure(2) imagesc(abs(fftshift(rect_hat))); title('fft(source)'); axis square saveas(gcf, 'fig2.png'); figure(3) imagesc(filt); title('filter (frequency domain)'); axis square saveas(gcf, 'fig3.png'); figure(4) imagesc(fftshift(abs(rect_hat_filtered))); title('fft(source) .* filter'); axis square saveas(gcf, 'fig4.png'); figure(5) imagesc(abs(rect_filtered)) title('result'); axis square saveas(gcf, 'fig5.png');
源图片:
源图像的傅里叶变换:
filter:
将滤波器应用(乘以)源图像的傅里叶变换的结果:
采用逆变换得出最终结果:
以上就是c/c++开发分享在C中使用FFTW的高通滤波器相关内容,想了解更多C/C++开发(异常处理)及C/C++游戏开发关注计算机技术网(www.ctvol.com)!)。
本文来自网络收集,不代表计算机技术网立场,如涉及侵权请联系管理员删除。
ctvol管理联系方式QQ:251552304
本文章地址:https://www.ctvol.com/c-cdevelopment/550005.html