Log-Gabor滤波器构造,opencv实现

算法 同时被 2 个专栏收录
5 篇文章 0 订阅
3 篇文章 0 订阅

参照链接:https://github.com/carandraug/PeterKovesiImage/blob/master/PhaseCongruency/gaborconvolve.m点击打开链接

毕设要用Log-Gabor滤波器来实现视网膜血管增强,这东西真是折腾了我好一段时间。。。。Matlab代码转到C++代码没动什么脑子,基本就是按逻辑翻译过来的

Log-Gabor滤波器的概念本文暂不表述了,待后期添加

首先是一些与滤波器无关的matlab转C++用 的代码:

void circshift(Mat& out, const Point& delta)
{
	Size sz = out.size();


	// 错误检查
	assert(sz.height > 0 && sz.width > 0);
	// 此种情况不需要移动
	if ((sz.height == 1 && sz.width == 1) || (delta.x == 0 && delta.y == 0))
		return;


	// 需要移动参数的变换,这样就能输入各种整数了
	int x = delta.x;
	int y = delta.y;
	if (x > 0) x = x % sz.width;
	if (y > 0) y = y % sz.height;
	if (x < 0) x = x % sz.width + sz.width;
	if (y < 0) y = y % sz.height + sz.height;




	// 多维的Mat也能移动
	vector<Mat> planes;
	split(out, planes);


	for (size_t i = 0; i < planes.size(); i++)
	{
		// 竖直方向移动
		Mat tmp0, tmp1, tmp2, tmp3;
		Mat q0(planes[i], Rect(0, 0, sz.width, sz.height - y));
		Mat q1(planes[i], Rect(0, sz.height - y, sz.width, y));
		q0.copyTo(tmp0);
		q1.copyTo(tmp1);
		tmp0.copyTo(planes[i](Rect(0, y, sz.width, sz.height - y)));
		tmp1.copyTo(planes[i](Rect(0, 0, sz.width, y)));


		// 水平方向移动
		Mat q2(planes[i], Rect(0, 0, sz.width - x, sz.height));
		Mat q3(planes[i], Rect(sz.width - x, 0, x, sz.height));
		q2.copyTo(tmp2);
		q3.copyTo(tmp3);
		tmp2.copyTo(planes[i](Rect(x, 0, sz.width - x, sz.height)));
		tmp3.copyTo(planes[i](Rect(0, 0, x, sz.height)));
	}


	merge(planes, out);
}


void fftshift(Mat& out)
{
	Size sz = out.size();
	Point pt(0, 0);
	pt.x = (int)floor(sz.width / 2.0);
	pt.y = (int)floor(sz.height / 2.0);
	circshift(out, pt);
}


void ifftshift(Mat& out)
{
	Size sz = out.size();
	Point pt(0, 0);
	pt.x = (int)ceil(sz.width / 2.0);
	pt.y = (int)ceil(sz.height / 2.0);
	circshift(out, pt);
}


void meshgrid(double xStart, double xEnd, double yStart, double yEnd,
              Mat& X, Mat& Y)
{
	std::vector<double> t_x, t_y;
	while (xStart <= xEnd)
	{
		t_x.push_back(xStart);


		++xStart;
	}
	while (yStart <= yEnd)
	{
		t_y.push_back(yStart);


		++yStart;
	}


	repeat(Mat(t_x).t(), t_y.size(), 1, X);
	repeat(Mat(t_y), 1, t_x.size(), Y);
}

下文需要用到的是meshgrid函数、fftshift、ifftshift函数,这两个函数是干嘛用的自行百度

mathlab里[x,y] = meshgrid(xrange, yrange); 等价于上述的 meshgrid(xrange.lower,xrange.upper,yrange.lower,yrange.upper)//大意伪代码

然后是Log-Gabor核构造需要的低通滤波器

Mat lowpassFilter(Size size, double cutOff, int n)
{
	int rows, cols;


	double xStart, xEnd, yStart, yEnd;


	if (cutOff < 0 || cutOff > 0.5)
		throw std::exception("cutoff frequency must be between 0 and 0.5");

	if (size.width == 1)
	{
		rows = size.width;
		cols = size.width;
	}
	else
	{
		rows = size.height;
		cols = size.width;
	}


	Mat x(cols, rows, CV_64FC1);
	Mat y(cols, rows, CV_64FC1);

	bool isRowsOdd = (rows % 2 == 1);
	bool isColsOdd = (cols % 2 == 1);
	if (isRowsOdd)
	{
		yStart = -(rows - 1) / 2.0;
		yEnd = (rows - 1) / 2.0;
	}
	else
	{
		yStart = -(rows / 2.0);
		yEnd = (rows / 2.0 - 1);
	}

	if (isColsOdd)
	{
		xStart = -(cols - 1) / 2.0;
		xEnd = (cols - 1) / 2.0;
	}
	else
	{
		xStart = -(cols / 2.0);
		xEnd = (cols / 2.0 - 1);
	}


	meshgrid(xStart, xEnd, yStart, yEnd, x, y);
	if (isColsOdd)
		x /= (cols - 1);
	else
		x /= cols;

	if (isRowsOdd)
		y /= (rows - 1);
	else
		y /= rows;


	Mat radius;
	Mat x2;
	Mat y2;


	pow(x, 2, x2);
	pow(y, 2, y2);
	sqrt(x2 + y2, radius);
	Mat f;
	Mat temp;
	pow((radius / cutOff), 2 * n, temp);
	divide(Mat::ones(radius.rows, radius.cols, CV_64F), 1.0 + temp, f);
	ifftshift(f);
	return f;
}

返回的mat即为构造的低通滤波器核

 

然后是一些矩阵运算操作

void complexMul(Mat left[2], Mat right[2], Mat* result)
{
	Mat real1;
	Mat real2;
	Mat imag1;
	Mat imag2;

	real1 = left[0].mul(right[0]);
	real2 = (left[1].mul(right[1]));
	imag1 = left[1].mul(right[0]);
	imag2 = left[0].mul(right[1]);

	result[0] = real1 - real2;
	result[1] = imag1 + imag2;
}

void complexDiv(Mat left[2], Mat right[2], Mat* result)
{
	Mat negRight[2];
	negRight[0] = right[0];
	negRight[1] = -right[1];

	Mat upper[2];
	Mat lower[2];
	complexMul(left, negRight, upper);
	complexMul(right, negRight, lower);


	divide(upper[0], lower[0], result[0]);
	divide(upper[1], lower[0], result[1]);
}

void complexAbs(Mat mat[2], Mat& result)
{
	Mat left, right;
	cv::pow(mat[0], 2, left);
	pow(mat[1], 2, right);
	sqrt(left + right, result);
}

注意,之所以是Mat[2]是因为 传入的应是复数,mat[0] [1]分别为实部和虚部

上述三个函数分别是计算 × ÷ 和  强度的(平方和开根)

最后是返回结果的结构体的定义

struct GaborConvolveResult
{

	vector<vector<Mat>>EO;
	vector<Mat>BP;
};

EO为 EO[nOrient,nScale],分别代表了不同角度和不同尺度的处理结果

BP为BP[nScale] ,分别表示了不同尺度下的处理结果(方向已融合)

下面就是Log-Gabor核的具体构造了

GaborConvolveResult garborConvolve(const Mat& mat, int nScale, int nOrient, double minWaveLength, double mult, double sigmaOnf,
                    double dThetaSigma, int Lnorm = 0, double feedback = 0)
{
	//mat应已确认是灰度图,并且是double
	int rows = mat.rows;
	int cols = mat.cols;
	

	Mat matDft;
	toFre(mat, matDft);
	
	

	vector<vector<Mat>>EO(nOrient, vector<Mat>(nScale, Mat(cols, rows, CV_64F,Scalar(0))));
	vector<Mat>BP(nScale, Mat(cols, rows, CV_64F, Scalar(0)));

	Mat x;
	Mat y;
	double xStart, xEnd, yStart, yEnd;
	bool isRowsOdd = (rows % 2 == 1);
	bool isColsOdd = (cols % 2 == 1);
	if (isRowsOdd)
	{
		yStart = -(rows - 1) / 2.0;
		yEnd = (rows - 1) / 2.0;
	}
	else
	{
		yStart = -(rows / 2.0);
		yEnd = (rows / 2.0 - 1);
	}

	if (isColsOdd)
	{
		xStart = -(cols - 1) / 2.0;
		xEnd = (cols - 1) / 2.0;
	}
	else
	{
		xStart = -(cols / 2.0);
		xEnd = (cols / 2.0 - 1);
	}




	meshgrid(xStart, xEnd, yStart, yEnd, x, y);

	if (isColsOdd)
		x /= (cols - 1);
	else
		x /= cols;

	if (isRowsOdd)
		y /= (rows - 1);
	else
		y /= rows;

	Mat radius;
	Mat x2;
	Mat y2;


	pow(x, 2, x2);
	pow(y, 2, y2);
	sqrt(x2 + y2, radius);

	Mat theta(rows, cols, CV_64FC1);

	for (int i = 0; i < rows; ++i)
		for (int j = 0; j < cols; ++j)
			theta.at<double>(i, j) = atan2(y.at<double>(i, j), x.at<double>(i, j));

	ifftshift(radius);
	ifftshift(theta);

	radius.at<double>(0, 0) = 1;

	Mat sinTheta(rows, cols, CV_64FC1);
	Mat cosTheta(rows, cols, CV_64FC1);

	for (int i = 0; i < rows; ++i)
		for (int j = 0; j < cols; ++j)
			sinTheta.at<double>(i, j) = sin(theta.at<double>(i, j));

	for (int i = 0; i < rows; ++i)
		for (int j = 0; j < cols; ++j)
			cosTheta.at<double>(i, j) = cos(theta.at<double>(i, j));

	Mat lp = lowpassFilter(Size(cols,rows), 0.45, 15);

	

	//vector<Mat>logGabor(nScale, Mat(rows, cols, CV_64F,Scalar(0,0)));
	vector<Mat>logGabor;
	for(int  s = 0;s<nScale;++s)
	{
		logGabor.push_back(Mat(rows, cols, CV_64F));
		double waveLength = minWaveLength* pow(mult, s );
		double fo = 1.0 / waveLength;

		

		Mat tempUpper;

		log(radius / fo, tempUpper);
		pow(tempUpper, 2, tempUpper);

		double tempLower = pow(log(sigmaOnf), 2);

		double factory = -1 / 2.0;
		tempUpper = tempUpper / tempLower * factory;



		exp(tempUpper, logGabor[s]);

		logGabor[s] = logGabor[s].mul(lp);
		logGabor[s].at<double>(0, 0) = 0;

		double L = 0;
		switch (Lnorm)
		{
		case 0:
			L = 1;
			break;
		case 1:
			{
			Mat planes[2] = { Mat(rows,cols,CV_64F),Mat(rows,cols,CV_64F) };
			Mat complex;
			idft(logGabor[s], complex, DFT_COMPLEX_OUTPUT + DFT_SCALE);
			split(complex, planes);
			Mat realPart = planes[0];

			L = sum(abs(realPart))[0];
			break;
			}
			
		case 2:
		{
			Mat temp;
			pow(logGabor[s], 2, temp);

			L = sqrt(sum(temp)[0]);

		}

			break;
		default:
			break;
		}

		logGabor[s] = logGabor[s] / L;
		//cout << logGabor[s] << endl;
		//cout << curLogGabor;
		

		Mat complex;
		Mat planes[2] = { Mat(rows,cols,CV_64F),Mat(rows,cols,CV_64F) };
		split(matDft, planes);

		planes[0] = planes[0].mul(logGabor[s]);
		planes[1] = planes[1].mul(logGabor[s]);

		//idft(complex, EO, DFT_COMPLEX_OUTPUT + DFT_SCALE);
		

		//cout << temp.depth() << "  " << temp.channels();
		//split(temp, planes);
		
		Mat complexd;
		merge(planes, 2, complexd);
		 idft(complexd,BP[s], DFT_COMPLEX_OUTPUT + DFT_SCALE);

	}
	//cout << logGabor[0] << endl;

	for(int o = 0 ; o < nOrient;++o)
	{
		double angl = o*CV_PI / nOrient;
		double waveLength = minWaveLength;


		Mat ds = sinTheta * cos(angl) - cosTheta * sin(angl);
	Mat dc = cosTheta * cos(angl) + sinTheta * sin(angl);

	Mat dTheta(rows, cols, CV_64F);
	for (int i = 0; i < rows; ++i)
		for (int j = 0; j < cols; ++j)
			dTheta.at<double>(i, j) = abs(atan2(ds.at<double>(i, j), dc.at<double>(i, j)));

	Mat temp;
	pow(dTheta, 2, temp);
	temp = -temp;
	Mat spread;
	double thetaSigma = CV_PI / nOrient / dThetaSigma;
	exp(temp / (2 * pow(thetaSigma, 2)), spread);

		for(int s =0 ;s<nScale;++s)
		{
			Mat filter = spread.mul(logGabor[s]);
			double L = 0;
			switch (Lnorm)
			{
			case 0: L = 1;
				break;
			case 1:
				{
				Mat planes[2] = { Mat(rows,cols,CV_64F),Mat(rows,cols,CV_64F) };
				Mat complex;
				idft(filter, complex, DFT_COMPLEX_OUTPUT + DFT_SCALE);
				split(complex, planes);
				Mat realPart = planes[0];
				L = sum(abs(realPart))[0];
				}
				break;
			case 2:
				{
				
				
				Mat planes[2] = { Mat(rows,cols,CV_64F),Mat(rows,cols,CV_64F) };
				
				split(temp, planes);
				Mat realPart = planes[0];
				
					
				Mat imagPart = planes[1];
				pow(realPart, 2, realPart);
				pow(imagPart, 2, imagPart);
				

				L = sqrt(sum(realPart)[0] + sum(imagPart)[0]);
				}
				break;
			default:
				break;
			}
			filter = filter / L;

			Mat complex;
			Mat planes[2] = { Mat(rows,cols,CV_64F),Mat(rows,cols,CV_64F) };
			cv::split(matDft, planes);
		
			planes[0] = planes[0].mul(filter);
			planes[1] = planes[1].mul(filter);
		
			merge(planes,2, complex);

			//here
			//Mat  multed = matDft.mul(filter);
			//cout << filter << endl;
			idft(complex, EO[o][s], DFT_COMPLEX_OUTPUT + DFT_SCALE);
			//cout << EO[s][o].cols << " " << EO[s][o].rows << EO[s][o].channels() << " " << EO[s][o].depth() << endl;

			Mat EOPlanes[2] = { Mat(rows,cols,CV_64F),Mat(rows,cols,CV_64F) };
			split(EO[o][s], EOPlanes);
		
			
			waveLength = waveLength *mult;
		}

		
	}



	GaborConvolveResult result;
	result.BP = BP;
	result.EO = EO;
	return result;
}

\

 

 

 

 

  • 0
    点赞
  • 5
    评论
  • 4
    收藏
  • 打赏
    打赏
  • 扫一扫,分享海报

©️2022 CSDN 皮肤主题:编程工作室 设计师:CSDN官方博客 返回首页

打赏作者

MaloryVer9

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值