opencv 低通滤波总结

<pre name="code" class="cpp">#include "StdAfx.h"
#include "cxcore.h"
#include "cv.h"
#include "highgui.h"


double D0=20;
 void ILPF(CvMat* src, const double D0)
{
	int i, j;
	int state = -1;
	double tempD;  
	long width, height;
	width = src->width;
	height = src->height;
    
	long x, y;
	x = width / 2;
	y = height / 2;


	CvMat* H_mat;
	H_mat = cvCreateMat(src->height,src->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)  (i * i + j * j);tempD = sqrt(tempD);break;
				case 1:
				tempD = (double)  ((height - i) * (height - i) + j * j);tempD = sqrt(tempD);break;
				case 2:
				tempD = (double)  (i * i + (width - j) * (width - j));tempD = sqrt(tempD);break;
				case 3:
				tempD = (double)  ((height - i) * (height - i) + (width - j) * (width - j));tempD = sqrt(tempD);break;
				default:
				break;
			}

	//二维高斯低通滤波器传递函数


			tempD = exp(-0.5 * pow(tempD / D0, 2));
			((double*)(H_mat->data.ptr + H_mat->step * i))[j * 2] = tempD;
			((double*)(H_mat->data.ptr + H_mat->step * i))[j * 2 + 1] = 0.0;


			//衰减系数为2的二维指数低通滤波器传递函数
			/* tempD = exp(-pow(tempD / D0, 2));   
			((double*)(H_mat->data.ptr + H_mat->step * i))[j * 2] = tempD;
			((double*)(H_mat->data.ptr + H_mat->step * i))[j * 2 + 1] = 0.0;*/


			//2阶巴特沃思低通滤波器传递函数
			/*tempD = 1 / (1 + pow(tempD / D0, 2 * 2));
			((double*)(H_mat->data.ptr + H_mat->step * i))[j * 2] = tempD;
			((double*)(H_mat->data.ptr + H_mat->step * i))[j * 2 + 1] = 0.0;*/




			//二维理想低通滤波器传递函数
			// 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;
			// }
		}
	}
	cvMulSpectrums(src, H_mat, src, CV_DXT_ROWS);
	cvReleaseMat(&H_mat);
  
}


int main(int argc, char ** argv)
{
    const char* filename = "E:\\work\\hh\\lena.jpg";
    IplImage * im;
 
    IplImage * realInput;
    IplImage * imaginaryInput;
    IplImage * complexInput;
    int dft_M, dft_N;
    CvMat* dft_A, tmp, *dft_B;
    IplImage * image_Re;
    IplImage * image_Im;
    double m, M;
 
    im = cvLoadImage( filename, CV_LOAD_IMAGE_GRAYSCALE );
    if( !im )
	{
        return -1;
	}
 
    realInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1);
    imaginaryInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1);
    complexInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 2);
 
    cvScale(im, realInput, 1.0, 0.0);
    cvZero(imaginaryInput);
    cvMerge(realInput, imaginaryInput, NULL, NULL, complexInput);
 
    dft_M = cvGetOptimalDFTSize( im->height - 1 );
    dft_N = cvGetOptimalDFTSize( im->width - 1 );
	dft_B = cvCreateMat( dft_M, dft_N, CV_64FC2 );
    dft_A = cvCreateMat( dft_M, dft_N, CV_64FC2 );
	cvZero(dft_B);


    image_Re = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);
    image_Im = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);
 
    
    cvGetSubRect( dft_A,&tmp, cvRect(0,0, im->width, im->height));
    cvCopy( complexInput, &tmp, NULL );
   
 
   
 
	cvDFT( dft_A, dft_A, CV_DXT_FORWARD, complexInput->height );

	ILPF(dft_A,  D0);
	cvDFT( dft_A, dft_A, CV_DXT_INVERSE , complexInput->height );


    cvNamedWindow("win", 0);
    cvNamedWindow("magnitude", 0);
    cvShowImage("win", im);
 
    
    cvSplit( dft_A, image_Re, image_Im, 0, 0 );
 
   
 
    cvMinMaxLoc(image_Re, &m, &M, NULL, NULL, NULL);
    cvScale(image_Re, image_Re, 1.0/(M-m), 1.0*(-m)/(M-m));




    cvShowImage("magnitude", image_Re);
 
    cvWaitKey(-1);
    return 0;
}

//next /

// Example : apply butterworth low pass filtering to input image/video
// usage: prog {<image_name> | <video_name>}

// Author : Toby Breckon, toby.breckon@cranfield.ac.uk

// Copyright (c) 2011 School of Engineering, Cranfield University
// License : LGPL - http://www.gnu.org/licenses/lgpl.html

#include <cv.h>   		// open cv general include file
#include <highgui.h>	// open cv GUI include file
#include <iostream>		// standard C++ I/O

using namespace cv; // OpenCV API is in the C++ "cv" namespace

/******************************************************************************/
// setup the cameras properly based on OS platform

// 0 in linux gives first camera for v4l
//-1 in windows gives first device or user dialog selection

#ifdef linux
	#define CAMERA_INDEX 0
#else
	#define CAMERA_INDEX -1
#endif
/******************************************************************************/
// Rearrange the quadrants of a Fourier image so that the origin is at
// the image center

void shiftDFT(Mat& fImage )
{
  	Mat tmp, q0, q1, q2, q3;

	// first crop the image, if it has an odd number of rows or columns

	fImage = fImage(Rect(0, 0, fImage.cols & -2, fImage.rows & -2));

	int cx = fImage.cols/2;
	int cy = fImage.rows/2;

	// rearrange the quadrants of Fourier image
	// so that the origin is at the image center

	q0 = fImage(Rect(0, 0, cx, cy));
	q1 = fImage(Rect(cx, 0, cx, cy));
	q2 = fImage(Rect(0, cy, cx, cy));
	q3 = fImage(Rect(cx, cy, cx, cy));

	q0.copyTo(tmp);
	q3.copyTo(q0);
	tmp.copyTo(q3);

	q1.copyTo(tmp);
	q2.copyTo(q1);
	tmp.copyTo(q2);
}

/******************************************************************************/
// return a floating point spectrum magnitude image scaled for user viewing
// complexImg- input dft (2 channel floating point, Real + Imaginary fourier image)
// rearrange - perform rearrangement of DFT quadrants if true

// return value - pointer to output spectrum magnitude image scaled for user viewing

Mat create_spectrum_magnitude_display(Mat& complexImg, bool rearrange)
{
    Mat planes[2];

    // compute magnitude spectrum (N.B. for display)
    // compute log(1 + sqrt(Re(DFT(img))**2 + Im(DFT(img))**2))

    split(complexImg, planes);
    magnitude(planes[0], planes[1], planes[0]);

    Mat mag = (planes[0]).clone();
    mag += Scalar::all(1);
    log(mag, mag);

    if (rearrange)
    {
        // re-arrange the quaderants
        shiftDFT(mag);
    }

    normalize(mag, mag, 0, 1, CV_MINMAX);

    return mag;

}
/******************************************************************************/

// create a 2-channel butterworth low-pass filter with radius D, order n
// (assumes pre-aollocated size of dft_Filter specifies dimensions)

void create_butterworth_lowpass_filter(Mat &dft_Filter, int D, int n)
{
	Mat tmp = Mat(dft_Filter.rows, dft_Filter.cols, CV_32F);

	Point centre = Point(dft_Filter.rows / 2, dft_Filter.cols / 2);
	double radius;

	// based on the forumla in the IP notes (p. 130 of 2009/10 version)
	// see also HIPR2 on-line

	for(int i = 0; i < dft_Filter.rows; i++)
	{
		for(int j = 0; j < dft_Filter.cols; j++)
		{
			radius = (double) sqrt(pow((i - centre.x), 2.0) + pow((double) (j - centre.y), 2.0));
			tmp.at<float>(i,j) = (float)
						( 1 / (1 + pow((double) (radius /  D), (double) (2 * n))));
		}
	}

    Mat toMerge[] = {tmp, tmp};
	merge(toMerge, 2, dft_Filter);
}

/******************************************************************************/

int main( int argc, char** argv )
{

  Mat img, imgGray, imgOutput;	// image object(s)
  VideoCapture cap; // capture object

  Mat padded;		// fourier image objects and arrays
  Mat complexImg, filter, filterOutput;
  Mat planes[2], mag;

  int N, M; // fourier image sizes

  int radius = 30;				// low pass filter parameter
  int order = 2;				// low pass filter parameter

  const string originalName = "Input Image (grayscale)"; // window name
  const string spectrumMagName = "Magnitude Image (log transformed)"; // window name
  const string lowPassName = "Butterworth Low Pass Filtered (grayscale)"; // window name
  const string filterName = "Filter Image"; // window nam

  bool keepProcessing = true;	// loop control flag
  int  key;						// user input
  int  EVENT_LOOP_DELAY = 40;	// delay for GUI window
                                // 40 ms equates to 1000ms/25fps = 40ms per frame

  // if command line arguments are provided try to read image/video_name
  // otherwise default to capture from attached H/W camera

    if(
	  ( argc == 2 && (!(img = imread( argv[1], CV_LOAD_IMAGE_COLOR)).empty()))||
	  ( argc == 2 && (cap.open(argv[1]) == true )) ||
	  ( argc != 2 && (cap.open(CAMERA_INDEX) == true))
	  )
    {
      // create window object (use flag=0 to allow resize, 1 to auto fix size)

      namedWindow(originalName, 0);
	  namedWindow(spectrumMagName, 0);
	  namedWindow(lowPassName, 0);
      namedWindow(filterName, 0);

        // if capture object in use (i.e. video/camera)
        // get image from capture object

		  if (cap.isOpened()) {

			  cap >> img;
			  if(img.empty()){
				if (argc == 2){
					std::cerr << "End of video file reached" << std::endl;
				} else {
					std::cerr << "ERROR: cannot get next fram from camera"
						      << std::endl;
				}
				exit(0);
			  }

		   }

      // setup the DFT image sizes

      M = getOptimalDFTSize( img.rows );
      N = getOptimalDFTSize( img.cols );

      // add adjustable trackbar for low pass filter threshold parameter

      createTrackbar("Radius", lowPassName, &radius, (min(M, N) / 2));
	  createTrackbar("Order", lowPassName, &order, 10);

	  // start main loop

	  while (keepProcessing) {

		  // if capture object in use (i.e. video/camera)
		  // get image from capture object

		  if (cap.isOpened()) {

			  cap >> img;
			  if(img.empty()){
				if (argc == 2){
					std::cerr << "End of video file reached" << std::endl;
				} else {
					std::cerr << "ERROR: cannot get next fram from camera"
						      << std::endl;
				}
				exit(0);
			  }

		  }

		  // ***

		    // convert input to grayscale

		  	cvtColor(img, imgGray, CV_BGR2GRAY);

			// setup the DFT images

	  		copyMakeBorder(imgGray, padded, 0, M - imgGray.rows, 0,
			      N - imgGray.cols, BORDER_CONSTANT, Scalar::all(0));
	  		planes[0] = Mat_<float>(padded);
		    planes[1] = Mat::zeros(padded.size(), CV_32F);

	  		merge(planes, 2, complexImg);

			// do the DFT

		    dft(complexImg, complexImg);

		    // construct the filter (same size as complex image)

		    filter = complexImg.clone();
		    create_butterworth_lowpass_filter(filter, radius, order);

		    // apply filter
		    shiftDFT(complexImg);
            mulSpectrums(complexImg, filter, complexImg, 0);
            shiftDFT(complexImg);

			// create magnitude spectrum for display

		    mag = create_spectrum_magnitude_display(complexImg, true);

            // do inverse DFT on filtered image

            idft(complexImg, complexImg);

            // split into planes and extract plane 0 as output image

            split(complexImg, planes);
            normalize(planes[0], imgOutput, 0, 1, CV_MINMAX);

            // do the same with the filter image

            split(filter, planes);
            normalize(planes[0], filterOutput, 0, 1, CV_MINMAX);

		  // ***

		  // display image in window

		  imshow(originalName, imgGray);
		  imshow(spectrumMagName, mag);
		  imshow(lowPassName, imgOutput);
		  imshow(filterName, filterOutput);

		  // start event processing loop (very important,in fact essential for GUI)
	      // 40 ms roughly equates to 1000ms/25fps = 4ms per frame

		  key = waitKey(EVENT_LOOP_DELAY);

		  if (key == 'x'){

	   		// if user presses "x" then exit

			  	std::cout << "Keyboard exit requested : exiting now - bye!"
				  		  << std::endl;
	   			keepProcessing = false;
		  }
	  }

	  // the camera will be deinitialized automatically in VideoCapture destructor

      // all OK : main returns 0

      return 0;
    }

    // not OK : main returns -1

    return -1;
}
/******************************************************************************/

matlab lowpass

function [H,D]=lpfilter(type,M,N,D0,n) 
%LPFILTER Computes frequency domain lowpass filters. 
% [H,D]=LPFILTER(TYPE,M,N,D0,n) creates the transfer function of a lowpass 
% filter, H, of the specified TYPE and size (M-by-N). To view the filter as 
% an image or mesh plot, it should be centered using H=fftshift(H). 
% 
% valid values for TYPE,D0, and n are: 
% 
% 'ideal'  Ideal lowpass filter with cutoff frequency D0. n need not be 
%          supplied. D0 must be positive. 
% 
%  'btw'   Butterworth lowpass filter of order n, and cutoff D0. The 
%          default value for n is 1.0. D0 must be positive. 
% 
%  'gaussian'  Gaussian lowpass filter with cutoff(standard deviation) D0. 
%              n need not be supplied. D0 must be positive. 
%Use funtion dftuv to set up the meshgrid arrays needed for computing the 
%required distances. 
[U,V]=dftuv(M,N); 
 
%compute the distances D(U,V) 
D=sqrt(U.^2+V.^2); 
 
%Begin filter computation. 
switch type 
    case 'ideal' 
        H=double(D<=D0); 
    case 'btw' 
        if nargin==4 
            n=1; 
        end 
        H=1./(1+(D./D0).^(2*n)); 
    case 'gaussian' 
        H=exp(-(D.^2)./(2*(D0^2))); 
    otherwise 
        error('Unknow filter type'); 
end