Gabor

Gabor滤波器改

float W_line = 0.56;/使用自己写的Gabor滤波器时，脊线宽度所占整个周期（脊线加谷线）的比例

float PI = 3.141592653589793;

double modPI(double x)将角度x转成一个周期内的表达形式
{
	if (x >= 0 && x < (2*PI))
		return x;
	int t = (int)(x /(2* PI));
	return (x - t*2*PI);
}

double myGaborCos(double x)
{
	double res;
	if (x < 0)
		x = -x;
	if (x - 6.3153 < 0.001 && x>6.3)
		x = x;
	double y = modPI(x);///y在0-2PI之间
	
	if (y >= 0 && y < PI*W_line)
	{
		res = cos(x / (2 * W_line));
	}
	else if (y >= PI*W_line && y < (2 * PI - PI*W_line))
	{
		res = -cos((x - PI) / (2 - 2 * W_line));
	}
	else if (y >= (2 * PI - PI*W_line) && y < 2 * PI)
	{
		res = cos((x - 2 * PI) / (2 * W_line));
	}

	return res;
	
	
}

//拷贝OpenCV源码的并且稍微改了下
Mat MygetGaborKernel(Size ksize, double sigma, double theta,
	double lambd, double gamma, double psi, int ktype = CV_64F)
{
	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;

			// Gabor滤波器的公式，将余弦值的计算换成实际需要的
			double v = scale*std::exp(ex*xr*xr + ey*yr*yr)*myGaborCos(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;
}