傅里叶变换用于特征匹配+相似性度量应用总结

最近经交流被问到傅里叶变换用于特征匹配的具体原理、及在解决分类问题时，为何欧式距离可以用于n-1维向量空间的相似性度量，奈何一时无语凝噎，难以用简洁通俗的语言来解释，故在此总结留念：

1.FFT用于特征匹配(即相位相关法)

图像配准的基本问题是找出一种图像转换方法,用以纠正图像的形变。造成图像形变的原因多种多样,例如对于遥感图像而言,传感器噪声、由传感器视点变化或平台不稳定造成的透视变化、被拍摄物体的移动、变形或生长等变化、闪电和大气变化,以及阴影和云层遮盖都使图像产生不同形式的形变。正是图像形变原因和形式的不同决定了多种多样的图像配准技术。

迄今已报道了多种图像配准方法,主要有互相关法、傅立叶变换法、点映射法口脚外和弹性模型法。其中傅立叶变换法基于傅立叶变换的相位匹配是利用傅立叶变换的性质而出现的一种图像配准方法。图像经过傅立叶变换,由空域变换到频率缘则两组数据在空何上的相关运算可以变为频谱的复数乘法运算,同时图像在变换域中还能获得在空域中很难获得的特征。
划重点 – 平移不影响傅氏变换的幅值（谱），对应的幅值谱和原图像是一样的。旋转在傅氏变换中是小变量。根据傅氏变换的旋转特性，旋转一幅图，在频域相当于对这幅图的傅氏变换做相同的旋转。使用频域方法的好处是计算简单，速度快（可使用MIT的fftw库）,同时傅立叶变换可以采用方法提高执行的速度。因此,傅氏变换是图像配准中常用的方法之一。下面我们就具体分析当图像发生平移、旋转和缩放时,图像信号在频域中的表现。



通过求取互功率谱的傅立叶反变换,得到一个狄拉克函数（脉冲函数）,再寻找函数峰值点对应的坐标,即可得到我们所要求得的配准点。实际上,在计算机处理中,连续域要用离散域代替,这使得狄拉克函数转化为离散时间单位冲击函数序列的形式。在实际运算中,两幅图像互功率谱相位的反变换,总是含有一个相关峰值代表两幅图像的配准点,和一些非相关峰值,相关峰值直接反映两幅图像间的一致程度。更精确的讲,相关峰的能量对应重叠区域的所占百分比,非相关峰对应非重叠区域所占百分比。由此我们可以看出,当两幅图像重叠区域较小时,采用本方法就不能检测出两幅图像的平移量。

算法流程：

调用opencv实现代码：

void PhaseCorrelation2D(const BYTE *signal,//原信号  
    const BYTE *pattern,//带配准信号  
    int &height_offset,//高的偏移量  
    int &width_offset)//宽的偏移量  
{  
  
    fftw_complex *signal_img = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*nRow*nCol);  
    fftw_complex *pattern_img = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*nRow*nCol);  
  
    for (int i = 0; i < nRow*nCol; i++)  
    {  
        signal_img[i][0] = signal[i];  
        signal_img[i][1] = 0;  
    }  
    for (int j = 0; j < nRow*nCol; j++)  
    {  
        pattern_img[j][0] = pattern[j];  
        pattern_img[j][1] = 0;  
    }  
  
    // 对两幅图像傅里叶变换  
    fftw_plan signal_forward_plan = fftw_plan_dft_2d(nRow, nCol, signal_img, signal_img,  
        FFTW_FORWARD, FFTW_ESTIMATE);  
    fftw_plan pattern_forward_plan = fftw_plan_dft_2d(nRow, nCol, pattern_img, pattern_img,  
        FFTW_FORWARD, FFTW_ESTIMATE);  
    fftw_execute(signal_forward_plan);  
    fftw_execute(pattern_forward_plan);  
  
    // 求互功率谱  
    fftw_complex *cross_img = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*nRow*nCol);  
    double temp;  
    for (int i = 0; i < nRow*nCol; i++)  
    {  
        cross_img[i][0] = (signal_img[i][0] * pattern_img[i][0]) -   
            (signal_img[i][1] * (-1.0*pattern_img[i][1]));  
        cross_img[i][1] = (signal_img[i][0] * (-1.0*pattern_img[i][1])) +   
            (signal_img[i][1] * pattern_img[i][0]);  
        temp = sqrt(cross_img[i][0] * cross_img[i][0] + cross_img[i][1] * cross_img[i][1]) + 0.001;  
        cross_img[i][0] /= temp;  
        cross_img[i][1] /= temp;  
    }  
  
    //对互功率谱求反变换  
    fftw_plan cross_backward_plan = fftw_plan_dft_2d(nRow, nCol, cross_img, cross_img,  
        FFTW_BACKWARD, FFTW_ESTIMATE);  
    fftw_execute(cross_backward_plan);  
  
    // 释放内存  
    fftw_destroy_plan(signal_forward_plan);  
    fftw_destroy_plan(pattern_forward_plan);  
    fftw_destroy_plan(cross_backward_plan);  
    fftw_free(signal_img);  
    fftw_free(pattern_img);  
  
    double *cross_real=new double[nRow*nCol];  
    for (int i = 0; i < nRow*nCol; i++)  
        cross_real[i] = cross_img[i][0];  
  
    int max_loc = 0;//准备存放最大值的位置坐标（注意，只有一个值）  
    double max_vlaue = 0.0;  
    for (int i = 0; i < nRow*nCol; i++)  
    {  
        if (max_vlaue<cross_real[i])  
        {  
            max_vlaue = cross_real[i];  
            max_loc = i;  
        }  
    }  
  
    height_offset = floor(((int)max_loc) / nCol);  
    width_offset = (int)max_loc - nCol*height_offset;  
  
    if (height_offset > 0.5*nRow)  
        height_offset = height_offset - nRow;  
    if (width_offset > 0.5*nCol)  
        width_offset = width_offset - nCol;  
  
    delete[] cross_real;  
    cross_real = NULL;  
}  
int main(int argc, char** argv)
{
	Mat srcImage11 = imread("Brain13x17y.bmp");
	Mat srcImage21 = imread("BrainP.bmp");
	Mat dst1, dst2;
	Mat srcImage1, srcImage2;
	cvtColor(srcImage11, srcImage1, CV_BGR2GRAY);
	cvtColor(srcImage21, srcImage2, CV_BGR2GRAY);


	srcImage1.convertTo(dst1, CV_32FC1);
	srcImage2.convertTo(dst2, CV_32FC1);

	QueryPerformanceFrequency(&freq);
	QueryPerformanceCounter(&startTime);//开始计时
	Point2d phase_shift;
	phase_shift = phaseCorrelate(dst1, dst2);
	cout << "OpneCV库函数实现 :" << endl << "\tX方向的偏移 : " << phase_shift.x << ",\tY方向的偏移 : " << phase_shift.y << endl;
	QueryPerformanceCounter(&stopTime);
	time = 1e3*(stopTime.QuadPart - startTime.QuadPart) / freq.QuadPart;

	cout << "你电脑的频率为：" << freq.QuadPart << endl;
	cout << "opencv库函数消耗时间为：" << time << "毫秒" << endl;
	
	QueryPerformanceCounter(&startTime);//开始计时
	int height_offset = 0;
	int width_offset = 0;
	PhaseCorrelation2D(srcImage1, srcImage2, height_offset, width_offset);//宽的偏移量  
	cout << "Phase Correlation自实现 :" << endl << "高度偏移量" << -height_offset << ",   宽度偏移量" << -width_offset << endl;
	QueryPerformanceCounter(&stopTime);//结束计时
	time = 1e3*(stopTime.QuadPart - startTime.QuadPart) / freq.QuadPart;
	cout << "自实现相位相关函数消耗时间为：" << time << "毫秒" << endl;
	waitKey(0);
	getchar();
	return 0;
}

参考：
数字图像处理-相位相关算法解决图像的刚性平移问题

2.相似性度量总结

Distance Classification
Distance
欧氏距离(Euclidean Distance)
闵可夫斯基距离（Minkowski distance）
曼哈顿距离(Manhattan distance)
切比雪夫距离 ( Chebyshev distance )
标准化欧氏距离(Standardized Euclidean distance )
马氏距离(Mahalanobis Distance)
汉明距离(Hamming distance)
余弦相似度（Cosine Similarity）
杰卡德相似系数(Jaccardsimilarity coefficient)
参考：
相似性度量的各种距离（Distance）计算归类详解及应用

3.支持向量机SVM原理分析

参考：
支持向量机 SVM