OpenCV实现傅里叶描述子(下): 形状特征提取

简述

在前面的更新中 OpenCV实现傅里叶描述子(上): 边界重建 有简单介绍了一下傅里叶描述子的内容，并利用边界的傅里叶描述子对其进行重建，使边界变得更加的平滑。但傅里叶描述子还是以其作为图像中形状的特征点得到广泛应用，例如手势识别，字符识别等。

常见的形状描述子有链码，傅里叶描述子和Hu不变距等，前二是基于形状轮廓的，Hu不变距是基于形状区域的，通常的形状特征描述子需要具有平移、缩放、旋转不变性，但是直接傅里叶变换得到的描述子与形状的尺度，方向，和曲线的起点位置是有关的，因此需要对傅里叶描述子进行处理以满足形状特征描述子不变性的要求。根据文献中推导并化简的公式为：
$$d(k)=\frac{||a(k)||}{||a(1)||}$$
上式$a(k)$为直接傅里叶变换得到的描述子，$d(k)$称为归一化傅里叶描述子，此描述子具有形状的旋转、平移、尺度不变性，并且和起始位置无关。

---

算法实现要点：

• 一般选择前10至15个描述子即可达到物体形状的描述作用，这10至15指的是低频系数，由于在进行傅里叶变换前，有对图像(一列)进行中心化，所以即为在图像中心$a(0)$两侧对称取10-15个归一化描述子系数。
• $d(1)==1$，所以$d(1)$项不用考虑，如果取12个描述子作为形状特征，在$a(0)$两侧各取6个，并去除$d(1)$，实际的描述子数量为11个。
• 对原图进行旋转、缩放、平移并使用傅里叶描述子提取特征，利用相关系数计算与原图的相似度。

代码实现

仅仅输出一个小demo，根据以上步骤做一个效果展示

代码如下：

#include<opencv2/core.hpp>
#include<opencv2/imgproc.hpp>
#include<opencv2/highgui.hpp>

#include<iostream>
#include<string>
#include<algorithm>



//傅里叶描述子特征提取
// contour:	形状轮廓
// number:	描述子数量
std::vector<double> _calc_FourierDescriptors(std::vector<cv::Point>& contour, int number = 12);

//相关系数计算
static double _calc_correlation(const std::vector<double>& lhs, const std::vector<double>& rhs);


int main()
{
	std::string path = "F:\\NoteImage\\moseme.tif";
	std::string path1 = "F:\\NoteImage\\descriptor.tif";

	cv::Mat src = cv::imread(path, cv::IMREAD_GRAYSCALE);
	cv::Mat othersrc = cv::imread(path1, cv::IMREAD_GRAYSCALE);
	if (!src.data) {
		std::cout << "Could not open or find the image" << std::endl;
		return -1;
	}
	
    /*********************测试代码*********************/
	cv::Mat _Template_testSrc;

	cv::threshold(src, _Template_testSrc, 50, 255, cv::THRESH_BINARY);
    //cv::threshold(othersrc, _Template_testSrc, 50, 255, cv::THRESH_BINARY);
    
	std::vector<std::vector<cv::Point>> contours;
	cv::findContours(_Template_testSrc, contours, cv::RETR_EXTERNAL, cv::CHAIN_APPROX_NONE, cv::Point());
	std::vector<cv::Point> conPoints = contours[0];			//先默认图片中就只有一条轮廓

	std::cout << "获取模板傅里叶描述子特征" << std::endl;
	std::vector<double> _Template = _calc_FourierDescriptors(conPoints);
	for (int i = 0; i < _Template.size(); i++)
	{
		printf("%d:%5.3f  ", (i + 1), _Template[i]);
	}
	std::cout << std::endl;



	std::cout << "获取变换图像后的傅里叶描述子特征" << std::endl;
	cv::Mat src1, src2, src3;
	cv::rotate(src, src1, 0);
	cv::rotate(src, src2, 1);
	cv::rotate(src, src3, 2);

	cv::Mat src4;
	cv::resize(src2, src4, cv::Size(), 2, 2);

	std::vector<cv::Mat> testMats = { src1, src2, src3, src4 };

	for (int k = 0; k < testMats.size(); ++k)
	{
		cv::Mat testSrc;

		cv::threshold(testMats[k], testSrc, 50, 255, cv::THRESH_BINARY);

		std::vector<std::vector<cv::Point>> contours;
		cv::findContours(testSrc, contours, cv::RETR_EXTERNAL, cv::CHAIN_APPROX_NONE, cv::Point());
		std::vector<cv::Point> conPoints = contours[0];			//先默认图片中就只有一条轮廓

		
		std::vector<double> descriptors = _calc_FourierDescriptors(conPoints);
		double corr = _calc_correlation(_Template, descriptors);

		for (int i = 0; i < descriptors.size(); i++)
		{
			printf("%d:%5.3f  ", (i + 1), descriptors[i]);
		}
		printf("  corr = %5.3f\n", corr);
	}

	cv::waitKey(0);
	return 0;
}


//傅里叶描述子特征提取
std::vector<double> _calc_FourierDescriptors(std::vector<cv::Point>& contour, int number)
{
	CV_Assert(contour.size() > (3 * number));
	number = (number & -2);							//取偶整数
	number = (number < 12 ? 12 : number);			//最小为12个点

	//将轮廓变为复数图像
	cv::Mat Point_Complex(contour.size(), 1, CV_32FC2);

	for (int i = 0; i < Point_Complex.rows; ++i)
	{
		cv::Vec2f* p = Point_Complex.ptr<cv::Vec2f>(i);
		p[0][0] = (float)contour[i].x;
		p[0][1] = (float)contour[i].y;
	}
	//将图像变换到适合傅里叶变换的最佳尺寸
	cv::Mat paddedMat;
	int m = cv::getOptimalDFTSize(Point_Complex.rows);
	int n = cv::getOptimalDFTSize(Point_Complex.cols);
	//验证填充后的轮廓数量是否为偶数(如果不为偶数，可以自行改为偶数，也不会耗费太多时长)
	CV_Assert(m % 2 == 0);

	cv::copyMakeBorder(Point_Complex, paddedMat, 0, m - Point_Complex.rows, 0, n - Point_Complex.cols,
		cv::BORDER_CONSTANT, cv::Scalar::all(0));

	//中心化(由于只有一列数据，x = 0)
	for (int y = 0; y < paddedMat.rows; ++y)
	{
		if ((y % 2) != 0)
		{
			paddedMat.ptr<cv::Vec2f>(y)[0][0] = -paddedMat.ptr<cv::Vec2f>(y)[0][0];
			paddedMat.ptr<cv::Vec2f>(y)[0][1] = -paddedMat.ptr<cv::Vec2f>(y)[0][1];
		}
	}

	//对轮廓图像进行傅里叶变换
	cv::Mat Point_complex_dft;
	cv::dft(paddedMat, Point_complex_dft, cv::DFT_COMPLEX_INPUT | cv::DFT_COMPLEX_OUTPUT);

	//取左右两侧的点
	int center = m / 2;			//z0
	int half_number = number / 2;
	std::vector<double> descriptors;
	double z1 = 0.0;

	for (int i = center - half_number; i <= center + half_number; ++i)
	{
		if (i != center)
		{
			cv::Vec2f* p = Point_complex_dft.ptr<cv::Vec2f>(i);
			double real = (double)p[0][0];
			double imag = (double)p[0][1];
			//求模
			double mag = std::sqrt(real * real + imag * imag);
			descriptors.push_back(mag);

			if (mag > z1)
				z1 = mag;
		}
	}

	CV_Assert(descriptors.size() == number);
	for (int i = 0; i < descriptors.size(); ++i)
		descriptors[i] /= z1;
	// 找到d1，并删除d1
	auto index = std::find(descriptors.begin(), descriptors.end(), 1.0);
	descriptors.erase(index);

	return descriptors;
}



//相关系数计算
static double _calc_correlation(const std::vector<double>& lhs, const std::vector<double>& rhs)
{
	double sum1 = 0, sum2 = 0;

	for (int i = 0; i < lhs.size(); ++i)
	{
		sum1 += lhs[i];
		sum2 += rhs[i];
	}

	double m1 = sum1 / lhs.size();
	double m2 = sum2 / rhs.size();

	double var1 = 0, var2 = 0;
	double numerator = 0;

	for (int i = 0; i < lhs.size(); ++i)
	{
		var1 += (lhs[i] - m1) * (lhs[i] - m1);
		var2 += (rhs[i] - m2) * (rhs[i] - m2);

		numerator += (lhs[i] - m1) * (rhs[i] - m2);
	}

	double _sqrt_var1var2 = std::sqrt(var1 * var2);

	return numerator / _sqrt_var1var2;
}

结果分析

还是使用染色体图片为例，对其进行旋转、坐标轴翻转并提取傅里叶描述子特征，利用相关系数方法进行相似度计算.

试验用图：

傅里叶描述子特征及其与原图特征的相关度：

当使用以下图作为模板与染色体的图片进行相关测试：

结果表明，不同形状轮廓的相关度大幅下降，表明无法匹配。
此外：本实验测试的图片量不是很多，暂时还未发现bug，如果真有网友根据此文去尝试并发现bug，还请私信与我。

---

OpenCV实现傅里叶描述子(上): 边界重建