利用频率陷波滤波处理图像

最新推荐文章于 2024-06-27 00:40:16 发布

Cross_Entropy

最新推荐文章于 2024-06-27 00:40:16 发布

阅读量1k

点赞数 1

分类专栏：数字视频处理文章标签：计算机视觉视频处理

本文链接：https://blog.csdn.net/Cross_Entropy/article/details/105585078

版权

数字视频处理专栏收录该内容

5 篇文章 0 订阅

订阅专栏

编程实现2-d DFT正变换和反变换
测试图像：moon.yuv （宽464，高538，4：4：4取样）
要求：(1) 调用1-D FFT模块实现2-D FFT
(2) 封装2-D FFT函数，通过调用2-D FFT实现2-D IDFT
频率域陷波滤波
测试图像：moon.yuv （宽464 高538）
实验平台：Visual C++
实验方法：频率域陷波滤波。具体要求如下：
- 陷波滤波：
  $f(x)=\left\{ \begin{aligned} x & = 0, (u,v) = (\frac{M}{2}, \frac{N}{2}) \\ y & = 1, else\\ \end{aligned} \right.$
  在滤波器中加入一个常量，以避免F(0,0)被完全消除,即：
  $H_e(u,v)=kH(u,v)+c$ ,为了方便，取k=1。
- 只需要对Y分量进行频率域高通滤波。Cb、Cr分量直接置为128即可。
- 频域滤波的具体步骤：
  - 用 $1)^{x+y}$ 乘以输入图像进行频率中心变换
  - 计算1中的DFT F(u, v)
  - 用滤波器函数H(u, v)乘以F(u, v)
  - 计算（3）中结果的反DFT
  - 得到（4）中结果的实部
  - 用 $1)^{x+y}$ 乘以5中的结果，取消输入图像的乘数

代码实现：

fft.h

void compute_W(int n, double *W_re, double *W_im); 
void permute_bitrev(int n, double *A_re, double *A_im); 
int  bitrev(int inp, int numbits); 
int  log_2(int n);  
void fft_butterfly(int n, double *A_re, double *A_im, double *W_re, double *W_im);
void fft_1d(int N, double *A_re, double *A_im);
void ifft_2d(double* Real, double* Imag, int width, int height);
void fft_2d(double* Real, double* Imag, int width, int height);

fft.cpp

#include <stdio.h> 
#include <stdlib.h>
#include <math.h>   
#include <assert.h>

#include "fft.h"

#define  M_PI 3.1415926


/*********************************************************************************************
对N点序列进行fft变换：
1. A_re为该序列的实部，A_im为该序列的虚部；
2. 运算结果仍然存放在A_re和A_im 
3. 输入输出都是自然顺序

Note：
N必须为2的整数次幂
**********************************************************************************************/
void fft_1d(int N, double *A_re, double *A_im) 
{ 
	double *W_re, *W_im;

	W_re = (double*)malloc(sizeof(double)*N/2); 
	W_im = (double*)malloc(sizeof(double)*N/2); 
	assert(W_re != NULL && W_im != NULL); 

	compute_W(N, W_re, W_im);     
	fft_butterfly(N, A_re, A_im, W_re, W_im);
	permute_bitrev(N, A_re, A_im); 

	free(W_re); 
	free(W_im);

}



/* W will contain roots of unity so that W[bitrev(i,log2n-1)] = e^(2*pi*i/n)
* n should be a power of 2
* Note: W is bit-reversal permuted because fft(..) goes faster if this is done.
*       see that function for more details on why we treat 'i' as a (log2n-1) bit number.
*/
void compute_W(int n, double *W_re, double *W_im)
{
	int i, br;
	int log2n = log_2(n);

	for (i=0; i<(n/2); i++)
	{
		br = bitrev(i,log2n-1); 
		W_re[br] = cos(((double)i*2.0*M_PI)/((double)n));  
		W_im[br] = -sin(((double)i*2.0*M_PI)/((double)n));  
	}

}


/* permutes the array using a bit-reversal permutation */ 
void permute_bitrev(int n, double *A_re, double *A_im) 
{ 
	int i, bri, log2n;
	double t_re, t_im;

	log2n = log_2(n); 

	for (i=0; i<n; i++)
	{
		bri = bitrev(i, log2n);

		/* skip already swapped elements */
		if (bri <= i) continue;

		t_re = A_re[i];
		t_im = A_im[i];
		A_re[i]= A_re[bri];
		A_im[i]= A_im[bri];
		A_re[bri]= t_re;
		A_im[bri]= t_im;
	}  
} 


/* treats inp as a numbits number and bitreverses it. 
* inp < 2^(numbits) for meaningful bit-reversal
*/ 
int bitrev(int inp, int numbits)
{
	int i, rev=0;
	for (i=0; i < numbits; i++)
	{
		rev = (rev << 1) | (inp & 1);
		inp >>= 1;
	}
	return rev;
}


/* returns log n (to the base 2), if n is positive and power of 2 */ 
int log_2(int n) 
{
	int res; 
	for (res=0; n >= 2; res++) 
		n = n >> 1; 
	return res; 
}

/* fft on a set of n points given by A_re and A_im. Bit-reversal permuted roots-of-unity lookup table
* is given by W_re and W_im. More specifically,  W is the array of first n/2 nth roots of unity stored
* in a permuted bitreversal order.
*
* FFT - Decimation In Time FFT with input array in correct order and output array in bit-reversed order.
*
* REQ: n should be a power of 2 to work. 
*
* Note: - See www.cs.berkeley.edu/~randit for her thesis on VIRAM FFTs and other details about VHALF section of the algo
*         (thesis link - http://www.cs.berkeley.edu/~randit/papers/csd-00-1106.pdf)
*       - See the foll. CS267 website for details of the Decimation In Time FFT implemented here.
*         (www.cs.berkeley.edu/~demmel/cs267/lecture24/lecture24.html)
*       - Also, look "Cormen Leicester Rivest [CLR] - Introduction to Algorithms" book for another variant of Iterative-FFT
*/
void fft_butterfly(int n, double *A_re, double *A_im, double *W_re, double *W_im) 
{
	double w_re, w_im, u_re, u_im, t_re, t_im;
	int m, g, b;
	int mt, k;

	/* for each stage */  
	for (m=n; m>=2; m=m>>1) 
	{
		/* m = n/2^s; mt = m/2; */
		mt = m >> 1;

		/* for each group of butterfly */ 
		for (g=0,k=0; g<n; g+=m,k++) 
		{
			/* each butterfly group uses only one root of unity. actually, it is the bitrev of this group's number k.
			* BUT 'bitrev' it as a log2n-1 bit number because we are using a lookup array of nth root of unity and
			* using cancellation lemma to scale nth root to n/2, n/4,... th root.
			*
			* It turns out like the foll.
			*   w.re = W[bitrev(k, log2n-1)].re;
			*   w.im = W[bitrev(k, log2n-1)].im;
			* Still, we just use k, because the lookup array itself is bit-reversal permuted. 
			*/
			w_re = W_re[k];
			w_im = W_im[k];

			/* for each butterfly */ 
			for (b=g; b<(g+mt); b++) 
			{
				/* t = w * A[b+mt] */
				t_re = w_re * A_re[b+mt] - w_im * A_im[b+mt];
				t_im = w_re * A_im[b+mt] + w_im * A_re[b+mt];

				/* u = A[b]; in[b] = u + t; in[b+mt] = u - t; */
				u_re = A_re[b];
				u_im = A_im[b];
				A_re[b] = u_re + t_re;
				A_im[b] = u_im + t_im;
				A_re[b+mt] = u_re - t_re;
				A_im[b+mt] = u_im - t_im;

			}
		}
	}
}

void fft_2d(double* Real, double* Imag, int width, int height)
{
	double* reLine, * imLine, * reRow, * imRow;
	reRow = new double[height]; imRow = new double[height];
	for (int i = 0; i < height; ++i)
		fft_1d(width, Real + i * width, Imag + i * width);

	for (int num = 0; num < width; ++num)
	{
		for (int i = 0; i < height; ++i)
		{
			reRow[i] = Real[i * width + num];
			imRow[i] = Imag[i * width + num];
		}
		fft_1d(height, reRow, imRow);
		for (int i = 0; i < height; ++i)
		{
			Real[i * width + num] = reRow[i];
			Imag[i * width + num] = imRow[i];
		}
	};
}
void ifft_2d(double* Real, double* Imag, int width, int height)
{
	for (int i = 0; i < width * height; ++i) Imag[i] = -Imag[i];
	fft_2d(Real, Imag, width, height);
}

main.cpp

#include <cstdio> 
#include <cstdlib>
#include <cmath>   
#include <cassert>
#include <iostream>
#include <fstream>
#include "fft.h"
#include <cstring>
using namespace std;

unsigned char filter(double x)
{
	if (x > 255) x = 255;
	if (x < 0) x = 0;
	return x;
}
int main(int argc, char* argv[])
{
	const string inPath = "moon.yuv";
	const string outPath = "monn_filted.yuv";
	int width = 464, height = 538;
	ifstream in( inPath, ios::binary );
	unsigned char* yBuffer, * uBuffer, * vBuffer;
	int size = width * height;
	for (auto i : { &yBuffer, &uBuffer, &vBuffer }) *i = new unsigned char[size];
	for (auto i : { &yBuffer, &uBuffer, &vBuffer }) in.read((char*)(*i), size);
	
	double* Real, * Imag;
	int newWidth = 1 << (log_2(width) + 1), newHeight = 1 << (log_2(height) + 1);
	int newSize = newWidth * newHeight;
	Real = new double[newSize]; Imag = new double[newSize];
	for (int i = 0; i < newSize; ++i) Imag[i] = Real[i] = 0;
	for (int i = 0; i < height; ++i)
		for (int j = 0; j < width; ++j)
			Real[i * newWidth + j] =  ((i+j)&1 ? -1.0 : 1) * yBuffer[i * width + j];
	
	fft_2d(Real, Imag, newWidth, newHeight);
	for (int i = 0; i < newSize; ++i)
	{
		Real[i] *= 1 + 0.5;
	}
	Real[newWidth / 2] = Imag[newHeight / 2] = 0;
	ifft_2d(Real, Imag, newWidth, newHeight);
	for (int i = 0; i < height; ++i)
	{
		for (int j = 0; j < width; ++j)
			yBuffer[i * width + j]  = filter(((i + j) & 1 ? -1.0 : 1) * Real[i * newWidth + j] / newSize  );
	}
	ofstream out(outPath, ios::binary);
	for (auto i : { &yBuffer, &uBuffer, &vBuffer }) out.write((char*)(*i), size);
}