1. 程式人生 > >在OpenCV環境下寫的灰度影象二維傅立葉換,幅值計算,頻譜平移和將數值歸一化到0到255區間的四個函式

在OpenCV環境下寫的灰度影象二維傅立葉換,幅值計算,頻譜平移和將數值歸一化到0到255區間的四個函式

影象處理開發資料、影象處理開發需求、影象處理接私活掙零花錢,可以搜尋公眾號"qxsf321",並關注!
影象處理開發資料、影象處理開發需求、影象處理接私活掙零花錢,可以搜尋公眾號"qxsf321",並關注!
影象處理開發資料、影象處理開發需求、影象處理接私活掙零花錢,可以搜尋公眾號"qxsf321",並關注!

灰度影象的二維傅立葉變換(cv_gray_fft2[我用C語言寫成]和DFT[用C++寫成]),二維傅立葉變換結果的幅值計算(cv_abs),頻譜平移(cv_gray_fft2shift)【頻譜平移的作用詳見我的博文http://blog.csdn.net/wenhao_ir/article/details/51689960

】,將數值歸一化到0到255區間(cv_range_0to255)是非常常用的四個功能!所以寫成四個函式,方便將來呼叫。

傅立葉變換的意義這裡就不多說了,大家在高等數學、訊號與系統中應該都學過,不清楚的可以去查閱相關資料!這裡我簡單說兩點!①影象的頻率是表徵影象中灰度變化劇烈程式的指標;②影象低頻分量表示背景和緩慢變化區域,高頻分量表示影象的邊緣、細節及相關噪聲。

原始碼如下

#include <opencv2/opencv.hpp>  
#include <opencv2/legacy/compat.hpp> 
#include <fstream>
using namespace std;  
#pragma comment(linker, "/subsystem:\"windows\" /entry:\"mainCRTStartup\"")  

void cv_gray_fft2(IplImage *src, IplImage *dst)  //注意:此函式僅適用灰度影象,src因為是灰度影象,所以通道數要求為1,dst的屬性為:IPL_DEPTH_64F,2(兩通道)
{
IplImage *image_Re = 0, *image_Im = 0, *Fourier = 0;
image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);  //實部
image_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);  //虛部
Fourier = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 2);//2 channels,分別儲存image_Re、image_Im

// Real part conversion from u8 to 64f (double)
cvConvertScale(src, image_Re);

// Imaginary part (zeros)
cvZero(image_Im);

// Join real and imaginary parts and stock them in Fourier image
cvMerge(image_Re, image_Im, 0, 0, Fourier);

// Application of the forward Fourier transform
cvDFT(Fourier, dst, CV_DXT_FORWARD);
cvReleaseImage(&image_Re);
cvReleaseImage(&image_Im);
cvReleaseImage(&Fourier);
}

void cv_abs(IplImage *src,IplImage *fft2_out, IplImage *dst)//src就是灰度影象原始資料,fft2_out要求是cv_gray_fft2的dst資料型別,即"64F, 2";dst為"64F, 1"
{
IplImage *image_Re = 0, *image_Im = 0;
image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);  //實部
image_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);  //虛部
cvSplit(fft2_out,image_Re,image_Im,0,0);  
cvPow(image_Re,image_Re,2);                 
    cvPow(image_Im,image_Im,2); 
cvAdd(image_Re,image_Im,image_Re,NULL);
cvPow(image_Re,dst,0.5);

cvReleaseImage(&image_Re);
cvReleaseImage(&image_Im);
}

void cv_gray_fft2shift(IplImage *src)//引數是"64F, 2"
{
int nRow, nCol, i, j, cy, cx;  
double tmp13, tmp24;

//Rearrange the quadrants of Fourier image so that the origin is at the image center  
    nRow = src->height;  
    nCol = src->width;  
    cy = nRow/2; // image center  
    cx = nCol/2;  
    //CV_IMAGE_ELEM為OpenCV定義的巨集,用來讀取影象的畫素值,這一部分就是進行中心變換  
    for( j = 0; j < cy; j++ ){  
        for( i = 0; i < cx; i++ ){  
            //中心化,將整體份成四塊進行對角交換  
            tmp13 = CV_IMAGE_ELEM( src, double, j, i);  
            CV_IMAGE_ELEM( src, double, j, i) = CV_IMAGE_ELEM(  
                src, double, j+cy, i+cx);  
            CV_IMAGE_ELEM( src, double, j+cy, i+cx) = tmp13;  
  
            tmp24 = CV_IMAGE_ELEM( src, double, j, i+cx);  
            CV_IMAGE_ELEM( src, double, j, i+cx) =  
                CV_IMAGE_ELEM( src, double, j+cy, i);  
            CV_IMAGE_ELEM( src, double, j+cy, i) = tmp24;  
        }  
    }  
}

void cv_range_0to255(IplImage *src, IplImage *dst)//兩個引數都是"64F, 2"
{
double m,M; 
double scale;  
double shift;  

cvMinMaxLoc(src,&m,&M,NULL,NULL); 
scale = 255/(M - m);  
shift = -m * scale;  
cvConvertScale(src,dst,scale,shift);  
}

int main()
{
int i=0;//迴圈變數
//從檔案中載入原圖  
IplImage *pSrcImage = cvLoadImage("coins.png", CV_LOAD_IMAGE_UNCHANGED);  

//建立輸出的影象
IplImage *pOutImage = cvCreateImage(cvGetSize(pSrcImage), IPL_DEPTH_64F,1);

//建立中間結果影象
IplImage *cv_gray_fft2_out = cvCreateImage(cvGetSize(pSrcImage), IPL_DEPTH_64F,2);//注意這裡是二通道!
IplImage *cv_abs_out = cvCreateImage(cvGetSize(pSrcImage),  IPL_DEPTH_64F, 1);

cv_gray_fft2(pSrcImage,cv_gray_fft2_out);

double watch_cv_gray_fft2_out_Re[100];//利用除錯時的區域性變數視窗觀察cvDFTOut實部的前100個值,看與MATLAB執行的結果否相同
double watch_cv_gray_fft2_out_Im[100];//利用除錯時的區域性變數視窗觀察cvDFTOut虛部的前100個值,看與MATLAB執行的結果否相同
for(i=0;i<10;i++)//一般情況下觀察10個值就夠了
{
watch_cv_gray_fft2_out_Re[i]=cvGet2D(cv_gray_fft2_out,0,i).val[0];  
}
for(i=0;i<10;i++)//一般情況下觀察10個值就夠了
{
watch_cv_gray_fft2_out_Im[i]=cvGet2D(cv_gray_fft2_out,0,i).val[1];  
}

cv_abs(pSrcImage,cv_gray_fft2_out,cv_abs_out);
double watch_cv_abs_out[100];//利用除錯時的區域性變數視窗觀察cv_abs_out的前100個值,看與MATLAB執行的結果否相同
for(i=0;i<10;i++)//一般情況下觀察10個值就夠了
{
watch_cv_abs_out[i]=cvGet2D(cv_abs_out,0,i).val[0];  
}

cv_gray_fft2shift(cv_abs_out);
double watch_cv_gray_fft2shift_out[100];//利用除錯時的區域性變數視窗觀察cv_gray_fft2shift_out的前100個值,看與MATLAB執行的結果否相同
for(i=0;i<10;i++)//一般情況下觀察10個值就夠了
{
watch_cv_gray_fft2shift_out[i]=cvGet2D(cv_abs_out,0,i).val[0];  
}

cv_range_0to255(cv_abs_out,pOutImage);
double watch_cv_range_0to255_out[100];//利用除錯時的區域性變數視窗觀察cv_range_0to255_out的前100個值,看與MATLAB執行的結果否相同
for(i=0;i<10;i++)//一般情況下觀察10個值就夠了
{
watch_cv_range_0to255_out[i]=cvGet2D(pOutImage,0,i).val[0];  
}
return 0;
}

上面四個函式與MATLAB執行的結果一致,MATLAB程式如下:

clear all;
clc;
I=imread('coins.png');
F=fft2(I);
MatlabFFt2Out=F;
MatlabAbsOut=abs(F);
MatlabFftshiftOut=fftshift(MatlabAbsOut);
T=MatlabFftshiftOut;
max_T=max(max(T));
min_T=min(min(T));
shift_T=-min_T*255/(max_T-min_T);
MatlabRange0to255=T.*255/(max_T-min_T)+shift_T;

執行結果截圖如下

cv_gray_fft2

cv_abs

cv_gray_fft2shift

cv_range_0to255

下面再給出OpenCV下的C++程式碼實現,程式碼中的DFT函式已經包含了上述函式的功能

//OpenCV版本2.4.9    
//交流QQ2487872782   
  
# include <opencv2/core/core.hpp>  
# include <opencv2/highgui/highgui.hpp>  
# include <opencv2/imgproc/imgproc.hpp>  
#include <iostream>  
using namespace cv;  
using namespace std;  
cv::Mat DFT(cv::Mat srcImage)  
{  
       cv::Mat srcGray;  
       cvtColor(srcImage,srcGray,CV_RGB2GRAY);  
       // 將輸入影象延擴到最佳的尺寸(快速傅立葉變換是基於影象尺寸是2、3或5的倍數完成的,因此對於輸入源影象,首先應將其轉換成DFTsize,OpenCV中提供了函式getOptimalDFTSize來實現尺寸轉換)  
       int nRows = getOptimalDFTSize(srcGray.rows);  
       int nCols = getOptimalDFTSize(srcGray.cols);  
       cv::Mat resultImage;  
       // 把灰度影象放在左上角,在右邊和下邊擴充套件影象,  
       // 新增的畫素初始化為0  
       copyMakeBorder(srcGray, resultImage, 0,  
            nRows - srcGray.rows,   
            0, nCols - srcGray.cols,   
            BORDER_CONSTANT, Scalar::all(0));  
       // 為傅立葉變換的結果(實部和虛部)分配儲存空間  
       cv::Mat planes[] = { cv::Mat_<float>(resultImage),   
            cv::Mat::zeros(resultImage.size(), CV_32F)};//這裡實際上是建立MAT陣列,陣列有兩個成員:  
														//第一個就是sizeConvMat這個物件(只是資料型別轉換成了float型別)  
                                                        //第二個是全0的型別為32F的物件  
       Mat completeI;  
       // 為延擴後的影象增添一個初始化為0的通道  
       merge(planes,2,completeI); //把groupMats的第0和第1個物件合併到mergeMat,通過這個操作mergeMat是雙通道的資料陣列了  
       // 進行離散傅立葉變換  
       dft(completeI,completeI);  
       // 將複數轉換為幅度  
       split(completeI,planes);  
       magnitude(planes[0],planes[1],planes[0]);//0中存的是實部值,1中存的是虛部值  
       cv::Mat dftResultImage = planes[0];   
       // 對數尺度(logarithmic scale)縮放  
       dftResultImage += 1;//陣列加1作對數變換,以擴大頻域動態顯示範圍  
       log(dftResultImage,dftResultImage);//作對數變換  
       // 剪下和重分佈幅度圖象限  
       dftResultImage= dftResultImage(Rect(0,  
            0,srcGray.cols,srcGray.rows));  
       // 歸一化影象  
       normalize(dftResultImage,dftResultImage,  
            0,1,CV_MINMAX);  
       int cx = dftResultImage.cols/2;  
       int cy = dftResultImage.rows/2;   
       Mat tmp;  
       // Top-Left - 為每一個象限建立ROI  
       Mat q0(dftResultImage,Rect(0,0,cx,cy));  
       // Top-Right  
       Mat q1(dftResultImage,Rect(cx,0,cx,cy));  
       // Bottom-Left  
       Mat q2(dftResultImage,Rect(0,cy,cx,cy));  
       // Bottom-Right  
       Mat q3(dftResultImage,Rect(cx,cy,cx,cy));   
       // 交換象限 (Top-Left with Bottom-Right)      
       q0.copyTo(tmp);  
       q3.copyTo(q0);  
       tmp.copyTo(q3);  
       // 交換象限 (Top-Right with Bottom-Left)  
       q1.copyTo(tmp);  
       q2.copyTo(q1);  
       tmp.copyTo(q2);  
       return dftResultImage;  
}  
int main()  
{  
      cv::Mat srcImage = imread("coins.png");  
      if(srcImage.empty())  
      return-1;  
      imshow("srcImage", srcImage);  
      cv::Mat resultImage = DFT(srcImage);  
      imshow("resultImage", resultImage);  
      cv::waitKey(0);  
      return 0;  
}  

執行結果如下圖所示