基於Gabor濾波器的紋理分割

Gabor變換是一種短時傅里葉變換方法，其實質是在傅里葉變換中加入一個窗函數，經過窗函數來實現信號的時頻分析。當選取高斯函數做爲窗函數時，短時傅里葉變換稱爲Gabor變換。在一維狀況中，Gabor變換表明着時頻分析的最優方法，二維狀況中則是空間頻域分析的方法。對於圖像來講，窗函數決定了它在空域的局部性，因此能夠經過移動窗口的中心來得到不一樣位置的空間域信息。因爲髙斯函數在通過傅里葉變換之後仍然是高斯函數，這就使得Gabor變換在頻域上仍然是局部的。簡直完美！

與傳統的傅立葉變換相比，Gabor小波變換具備良好的時頻局部化特性。即很是容易地調整Gabor濾波器的方向、基頻帶寬及中心頻率從而可以最好的兼顧信號在時空域和頻域中的分辨能力；Gabor小波變換具備多分辨率特性即變焦能力；採用多通道濾波技術，將一組具備不一樣時頻域特性的Gabor小波應用於圖像變換，每一個通道都可以獲得輸入圖像的某種局部特性，這樣能夠根據須要在不一樣粗細粒度上分析圖像。

Gabor濾波器紋理分析方法是一種重要的基於變換的紋理特徵提取方法，經過選用某一特定的Gabor函數，設計Gabor濾波器，來實現圖像的多尺度、多方向的特徵提取。

gabor濾波器的參數不少，看着就讓人頭暈，搞清楚各個參數的物理含義及相互之間的關係後，公式就好理解了。後面使用的話設好參數就能夠直接用。

標準差σx、σy分別爲高斯函數沿 x 軸和 y 軸的標方差， 控制高斯窗口做用的範圍，σx、σy越大，高斯窗做用的範圍也越大。變換σx、σy可使gabor濾波器匹配圖像中不一樣空間尺度的結構。一般取σx=σy。

這裏的縱橫比y/x=1。關於縱橫比說一下，縱橫比控制高斯函數橫截面的橢圓形狀，爲1時截面就是圓。

(U，V)是中心頻率，其頻率調製的方向角φ=arctan(V/U)。

對高斯基本函數g(x，y)進行尺度變換和旋轉變換獲得一簇自類似函數，即Gabor小波，這樣gabor濾波器就能夠在不一樣尺度、方向上分析紋理，粗糙的紋理用低分辨率分析，細緻的紋理用高分辨率去分析。其中尺度變換和旋轉函數以下

θ 表示空域座標繞x軸逆時針旋轉θ 角，對θ 離散化， 令θ=n*π /N ,n= 0,1, N-1，N 爲整數.則可用n定義濾波器的方向。改變θ，濾波器就能夠在不一樣方向上檢測紋理。a爲伸縮因子，a>1，通常取a=1.414，m 定義濾波器的尺度，m=1,2,...M-1。

一般狀況下，取Gabor函數中高斯函數部分的方向θ與複數調製部分函數的輻角φ相等，即令θ=φ。

令

則gabor函數可化簡爲

定義波長（λ）：λ是COS調製因子的波長，它的值以像素爲單位指定，一般大於等於2.但不能大於輸入圖像尺寸的五分之一。 λ和F有以下關係：

                                                           λ=2*π/F

要實現多尺度檢測的話，要麼改變波長，要麼改變中心頻率，濾波器的中心頻率越小（即波長越長），提取的紋理特徵的尺度越大。當濾波器紋理與圖像做用時，濾波器覆蓋下的局部紋理頻率與濾波器的頻率越接近響應就越大，反之越小。

φ是cos調製因子的相位偏移值，φ決定了Gabor函數的對稱性，好比在φ=0,π的時候，Gabor函數是中心（對於（U，V）來講）對稱的，而當 φ=-π/2，π/2的時候，Gabor函數是中心反對稱的，而且全部其餘狀況都是這兩種狀況的組合。

gabor函數能夠表示爲一個橢圓高斯函數與一個複平面函數的乘積形式。複數形式的二維 gabor函數包含兩個正交投影： 實數的偶對稱餘弦部分和虛數的奇對稱正弦部分。通常根據須要提取實數部分或虛數部分或複數幅值。

#ifndef __GABOR__
#define __GABOR__


class gabor
{
//#define PI 3.1415926
	typedef unsigned char BYTE;
	const double phase=0;//相位The phase offset φ in the 
	//argument of the cosine factor of the Gabor function is specified in degrees
	//const double gama=1.8;//縱橫比
	const double OriBandwidth = PI / 6;//Orientation Bandwidth of a Gabor filter
	const double FreBandwidth = 1;
	//const int kernel_size = 21;
	const double bandwidth = 1;
	const int theta_bins = 8;const int img_size = 256;
private:
	//void resize256();//若是最長邊大於256，則將圖像按比例縮放到256*，爲了減少運算量
	void get_kernel(double*&kernel,const double theta,const double lamda,int &kernel_w,int&kernel_h);
	//void low_pass_filter();
	//double calculateSigma(double waveLength, double bandwidth);
	void filter2D(double*kernel, const int kernel_w, const int kernel_h, 
		BYTE*img,const int img_w,const int img_h,double**&out,int cnt);
public:
	gabor();
	~gabor();
	void apply_filter(BYTE*grey_img,const int imgH,const int imgW,double**&output);

};

#endif

#include "stdafx.h"
#include "gabor.h"
#include<cmath>

gabor::gabor()
{
}

gabor::~gabor()
{
	
}

//double gabor::gabor_func(BYTE*Image, const int imgH, const int imgW, const double oriBand,
//	const double theta, const double lambda, const double phi, const double freBand)
//{
//	//determine Sx and Sy by Oritation Frequency and Bandwidth Bandwidth
//	double Sx = lambda / PI*sqrt(log(2) / 2)*((pow(freBand,2) + 1) / (pow(freBand,2) - 1));
//	double Sy = lambda / PI* sqrt(log(2) / 2) / tan(oriBand / 2);
//
//}


void gabor::apply_filter(BYTE*grey_img, const int imgH, const int imgW, double**&output)
{
	int lamda_n = (log(img_size / 8.0) + 0.01) / log(2.0);
	double*lamda=new double[lamda_n];//The wavelength of the cosine factor of the Gabor filter
	for (int i = 0; i < lamda_n; i++)
	{
		lamda[i] = 2 * pow(sqrt(2.0), double(i));
		printf("%lf    ", lamda[i]);
	}
	double*kernel = NULL;
	double*theta = new double[theta_bins];//The orientation of Gabor filter
	for (int i = 0; i < theta_bins; i++)
		theta[i] = double(i) / theta_bins*PI;
	
	output = new double*[imgH * imgW];
	for (int i = 0; i < imgH*imgW; i++)
		{
			output[i] = new double[theta_bins * lamda_n];
		}
	for (int i = 0; i < lamda_n; i++)
		for (int j = 0; j < theta_bins; j++)
		{
			int kernel_w=0, kernel_h=0;
			get_kernel(kernel,theta[j], lamda[i], kernel_w, kernel_h);
			filter2D(kernel, kernel_w, kernel_h, grey_img, img_size, img_size, output, i * 6 + j);
			delete[]kernel;
			kernel = NULL;
		}
	delete[]lamda,theta;
}


void gabor::get_kernel(double*&kernel, const double theta, const double lamda, int &kernel_w, int&kernel_h)
{
	double Sx = lamda / PI*sqrt(log(2) / 2)*((pow(2, FreBandwidth) + 1) / (pow(2 ,FreBandwidth) - 1));
	double Sy = lamda / PI* sqrt(log(2) / 2) / tan(OriBandwidth / 2);

	double Sigma_max = Sx; if (Sy>Sx)Sigma_max=Sy;

	//Gabor Mask Size
	kernel_w = ceil(Sigma_max * 3) * 2 + 1;
	kernel_h = ceil(Sigma_max * 3) * 2 + 1;

	if (kernel)
		delete[]kernel;
	kernel = new double[kernel_w*kernel_w];

	for (int m = 0; m < kernel_h;m++)
		for (int n = 0; n < kernel_w; n++)
		{
			int x = m - ceil(kernel_w / 2);
			int y = n - ceil(kernel_h / 2);
			double xPrime = x * cos(theta) + y * sin(theta);
			double yPrime = y * cos(theta) - x * sin(theta);
			kernel[m*kernel_w + n] = 1 / (2 * PI*Sx*Sy)*exp(-0.5*(pow(xPrime / Sx, 2.0) 
				+ (yPrime / Sy)*(yPrime / Sy)))*cos(2 * PI*xPrime / lamda + phase);
			//printf("%lf    ", kernel[m*kernel_w + n]);
		}
	//printf("\n");
}

void gabor::filter2D(double*kernel, const int kernel_w, const int kernel_h, BYTE*img,
	const int img_w, const int img_h, double**&out, int chanel)
{
	BYTE*img_data = new BYTE[img_h*img_w];

	for (int i = 0; i < img_h; i++)
		for (int j = 0; j < img_w; j++)
		{
			double temp = 0;
			int cnt = 0;
			for (int m = -kernel_h/2; m <= kernel_h/2; m++)
				for (int n = -kernel_w/2; n <= kernel_w/2; n++)
				{
					double img_value = (i + m >= 0 && i + m < img_h
						&&j + n >= 0 && j + n < img_w) ? img[(i + m)*img_w + j + n] : 0;
					temp += img_value*kernel[cnt++];
				}
			out[i*img_w + j][chanel] = abs(temp);
			img_data[i*img_w + j] = abs(50*temp);
		}
	
	/*save_gray("11.bmp", img_data, img_h, img_w);
	CxImage cx;
	cx.Load("11.bmp", CXIMAGE_FORMAT_BMP);
	int templates[25] = { 1, 4, 7, 4, 1,
		4, 16, 26, 16, 4,
		7, 26, 41, 26, 7,
		4, 16, 26, 16, 4,
		1, 4, 7, 4, 1 };
	cx.Filter(templates, 25, 273, 0);*/
	delete[]img_data;
	return ;
}

opencv裏的gabor核函數生成代碼

/*
 Gabor filters and such. To be greatly extended to have full texture analysis.
 For the formulas and the explanation of the parameters see:
 http://en.wikipedia.org/wiki/Gabor_filter
*/

cv::Mat cv::getGaborKernel( Size ksize, double sigma, double theta,
                            double lambd, double gamma, double psi, int ktype )
{
    double sigma_x = sigma;
    double sigma_y = sigma/gamma;
    int nstds = 3;
    int xmin, xmax, ymin, ymax;
    double c = cos(theta), s = sin(theta);

    if( ksize.width > 0 )
        xmax = ksize.width/2;
    else
        xmax = cvRound(std::max(fabs(nstds*sigma_x*c), fabs(nstds*sigma_y*s)));

    if( ksize.height > 0 )
        ymax = ksize.height/2;
    else
        ymax = cvRound(std::max(fabs(nstds*sigma_x*s), fabs(nstds*sigma_y*c)));

    xmin = -xmax;
    ymin = -ymax;

    CV_Assert( ktype == CV_32F || ktype == CV_64F );

    Mat kernel(ymax - ymin + 1, xmax - xmin + 1, ktype);
    double scale = 1;
    double ex = -0.5/(sigma_x*sigma_x);
    double ey = -0.5/(sigma_y*sigma_y);
    double cscale = CV_PI*2/lambd;

    for( int y = ymin; y <= ymax; y++ )
        for( int x = xmin; x <= xmax; x++ )
        {
            double xr = x*c + y*s;
            double yr = -x*s + y*c;

            double v = scale*exp(ex*xr*xr + ey*yr*yr)*cos(cscale*xr + psi);
            if( ktype == CV_32F )
                kernel.at<float>(ymax - y, xmax - x) = (float)v;
            else
                kernel.at<double>(ymax - y, xmax - x) = v;
        }

    return kernel;
}

        

分割步驟

對待分割圖像F(m，n)進行多尺度、多方向的Gabor變換。

對獲得的各子帶採用非線性變換函數進行變換，並經過低通濾波來下降相同紋理區域內特徵的變換，增長不一樣區域的區別。

將上述處理得的各子帶對應的像素組成向量

採用分類算法對特徵向量進行分類，實現紋理圖像的分割。

Log Gabor