实验实验 5：频域滤波 名称 实验目的 掌握图像进行频域滤波的方法和步骤。 1、掌握图像频域 DFT 变换和反变换的方法。 2、掌握图像频域滤波的步骤 1、灰度图像的 DFT 和 IDFT。 具体内容：利用 OpenCV 提供的 cvDFT 函数对图像进行 DFT 和 IDFT 变换 2、利用理想高通和低通滤波器对灰度图像进行频域滤波 具体内容：利用 cvDFT 函数实现 DFT，在频域上利用理想高通和低通滤波器进行滤波，并把滤波过后的图像显示在屏幕上（观察振铃现象），要求截止频率可输入。 3、利用布特沃斯高通和低通滤波器对灰度图像进行频域滤波。 具体内容：利用 cvDFT 函数实现 DFT，在频域上进行利用布特沃斯高通和低通滤波器进行滤波，并把滤波过后的图像显示在屏幕上（观察振铃现象），要求截止频率和 n 可输入。 1、 实验步骤：利用 OpenCV 提供的 cvDFT 函数对图像进行 DFT 和 IDFT 变换 核心代码如下： //DFT变换 IplImage *DFT(IplImage * src) { IplImage* fourier = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,2); int dft_H, dft_W; dft_H = src->height; dft_W = src->width; CvMat *src_Re = cvCreateMat(dft_H,dft_W, CV_64FC1); // double Re, Im; CvMat *src_Im = cvCreateMat(dft_H,dft_W, CV_64FC1); //Imaginary part CvMat *sum_src =cvCreateMat(dft_H,dft_W, CV_64FC2); //2 channels (src_Re, src_Im) CvMat *sum_dst =cvCreateMat(dft_H,dft_W, CV_64FC2); //2 channels (dst_Re, dst_Im) cvConvert(src, src_Re); cvZero(src_Im); cvMerge(src_Re, src_Im, 0, 0, sum_src); cvDFT(sum_src,sum_dst,CV_DXT_FORWARD,0); cvConvert(sum_dst, fourier); cvReleaseMat(&src_Re); cvReleaseMat(&src_Im); cvReleaseMat(&sum_src); cvReleaseMat(&sum_dst); return fourier; } 实验内容 实验完成情况
//DFT反变换 IplImage *IDFT(IplImage * fourier) { IplImage* dst = cvCreateImage(cvGetSize(fourier),IPL_DEPTH_8U,1); int dft_H, dft_W; dft_H = fourier->height; dft_W = fourier->width; CvMat *dst_Re = cvCreateMat(dft_H,dft_W, CV_64FC1); // double Re, Im; CvMat *dst_Im = cvCreateMat(dft_H,dft_W, CV_64FC1); //Imaginary part CvMat *sum_dst =cvCreateMat(dft_H,dft_W, CV_64FC2); //2 channels (dst_Re, dst_Im) CvMat *sum_src = cvCreateMat(dft_H,dft_W, CV_64FC2 ); cvConvert(fourier, sum_src); cvDFT(sum_src,sum_dst,CV_DXT_INV_SCALE,0); cvSplit(sum_dst,dst_Re,dst_Im,0,0); cvConvert(dst_Re, dst); cvReleaseMat(&dst_Re); cvReleaseMat(&dst_Im); cvReleaseMat(&sum_src); cvReleaseMat(&sum_dst); return dst; } //归一化，将灰度映射到0~255之间, 并将能量最高的四角移到中心, 生成图片频域能量图 void BuildDFTImage(IplImage *fourier, IplImage *dst) { IplImage *image_Re = 0, *image_Im = 0; image_Re = cvCreateImage(cvGetSize(fourier), IPL_DEPTH_64F, 1); image_Im = cvCreateImage(cvGetSize(fourier), IPL_DEPTH_64F, 1); //Imaginary part cvSplit(fourier, image_Re, image_Im, 0, 0 ); // Compute the magnitude of the spectrum Mag = sqrt(Re^2 + Im^2) cvPow( image_Re, image_Re, 2.0); cvPow( image_Im, image_Im, 2.0); cvAdd( image_Re, image_Im, image_Re); cvPow( image_Re, image_Re, 0.5 ); cvReleaseImage(&image_Im);
cvAddS(image_Re, cvScalar(1.0), image_Re); // 1 + Mag cvLog(image_Re, image_Re ); // log(1 + Mag) //重新安排傅里叶图像中心 // Rearrange the quadrants of Fourier image so that the origin is at // the image center double minVal = 0, maxVal = 0; cvMinMaxLoc( image_Re, &minVal, &maxVal ); // Localize minimum and maximum values CvScalar min; min.val[0] = minVal; double scale = 255 / (maxVal - minVal); cvSubS(image_Re, min, image_Re); cvConvertScale(image_Re, dst, scale); cvReleaseImage(&image_Re); // Rearrange the quadrants of Fourier image so that the origin is at // the image center int nRow, nCol, i, j, cy, cx; uchar tmp13, tmp24; nRow = fourier->height; nCol = fourier->width; cy = nRow/2; // image center cx = nCol/2; for( j = 0; j < cy; j++ ) { for( i = 0; i < cx; i++ ) { tmp13 = CV_IMAGE_ELEM( dst, uchar, j, i); CV_IMAGE_ELEM( dst, uchar, j, i) = CV_IMAGE_ELEM(dst, uchar, j+cy, i+cx); CV_IMAGE_ELEM( dst, uchar, j+cy, i+cx) = tmp13; tmp24 = CV_IMAGE_ELEM( dst, uchar, j, i+cx); CV_IMAGE_ELEM( dst, uchar, j, i+cx) = CV_IMAGE_ELEM( dst, uchar, j+cy, i); CV_IMAGE_ELEM( dst, uchar, j+cy, i) = tmp24; } } } 实验结果如图：
2、 实验步骤：利用 cvDFT 函数实现 DFT，在频域上利用理想高通和低通滤波器进行滤波，并把滤波过后的图像显示在屏幕上（观察振铃现象），截止频率可输入。 核心代码如下： void PassFilter(IplImage * fourier, int FLAG, double d0, int n1) { int i, j; int state = -1; double tempD; long width, height; width = fourier->width; height = fourier->height; long x, y; x = width / 2; y = height / 2; CvMat* H_mat;
H_mat = cvCreateMat(fourier->height,fourier->width, CV_64FC2); for(i = 0; i < height; i++){ for(j = 0; j < width; j++){ if(i > y && j > x){ state = 3; }else if(i > y){ state = 1; }else if(j > x){ state = 2; }else{ state = 0; } switch(state){ case 0: tempD = (double)sqrt(1.0*i * i + j * j);break; case 1: tempD = (double)sqrt(1.0*(height - i) * (height - i) + j * j);break; case 2: tempD = (double)sqrt(1.0*i * i + (width - j) * (width - j));break; case 3: tempD = (double)sqrt(1.0*(height - i) * (height - i) + (width - j) * (width - j));break; default: break; } switch(FLAG){ case IDEAL_LOW: if(tempD <= D0){ ((double*)(H_mat->data.ptr + H_mat->step * i))[j * 2] = 1.0; ((double*)(H_mat->data.ptr + H_mat->step * i))[j * 2 + 1] = 0.0; }else{ ((double*)(H_mat->data.ptr + H_mat->step * i))[j * 2] = 0.0; ((double*)(H_mat->data.ptr + H_mat->step * i))[j * 2 + 1] = 0.0; } break; case IDEAL_HIGH: if(tempD <= D0){ ((double*)(H_mat->data.ptr + H_mat->step * i))[j * 2] = 0.0; ((double*)(H_mat->data.ptr + H_mat->step * i))[j * 2 + 1] = 0.0; }else{