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

简述

傅里叶描述子在数字图像处理中是属于特征提取的内容，其最主要的作用之一是应用在形状识别中，如字符识别。其描述子组成的特征向量经过适当的处理可以具有旋转、平移和尺度不变性。

这次博客的内容主要是边界重建，形状识别放在以后的更新当中，下面是关于一些学习内容的总结。

内容总结及算法步骤:

• 将由K个点组成的边界(轮廓点)的坐标对$(x_0,y_0),(x_1,y_1),\cdots ,(x_{k-1},y_{k-1})$，视为一个个复数，于是有：$s(k)=x(k)+\text{j}y(k)$，他将二维描述问题转化为一维描述问题。
• 对$s(k)$进行离散傅里叶变换，获得复数系数$a(u)$称为边界的傅里叶描述子：
  $$a(u)=\sum _{K=0}^{K-1}s(k){\exp }^{-j2\pi uk/K}$$
• 重建过程就发生在这一步的反变换过程中，如果用了所有的傅里叶级数，那么反变换就等于原输入，但假设只使用前$P$个系数，大于$P$令为0，即$a(u>p)=0$，此时的结果对反变换的$s(k)$有如下近似：
  $$\hat{s}(k)=\frac{1}{K}\sum _{u=0}^{p-1}a(u){\exp }^{j2\pi uk/K}$$
  参数$k$的值域仍然是0到K-1，即近似边界中存在相同数量的点，但是未在每个点的重建中使用那么多项。

---

其实重建过程和频率域低通滤波是一样的，都是屏蔽高频系数，保留低频然后在反变换，所以在傅里叶变换之前要先进行中心化，反变换回来再次去中心化，由于DFT的对称性，边界上的点数及其反变换必须是偶数，意味着计算反变换之前，删除的系数(设置为0)数量必须是偶数，因为有中心化，所以删除系数是将系数两端的一半数量设为0。高频决定细节，低频决定整体形状，设置为0的数量越多，边界细节丢失越多。

代码实现

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

代码如下：

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

#include<iostream>
#include<string>


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

	cv::Mat src = cv::imread(path, cv::IMREAD_GRAYSCALE);
	if (!src.data) {
		std::cout << "Could not open or find the image" << std::endl;
		return -1;
	}

	cv::threshold(src, src, 50, 255, cv::THRESH_BINARY);

	std::vector<std::vector<cv::Point>> contours;
	cv::findContours(src, contours, cv::RETR_EXTERNAL, cv::CHAIN_APPROX_NONE, cv::Point());
	//定义一个画板
	cv::Mat drawmat = cv::Mat::zeros(src.size(), CV_8UC3);
	
	std::vector<cv::Point> conPoints = contours[0];			//先默认图片中就只有一条轮廓
	//将轮廓变为复数图像
	cv::Mat Point_Complex(conPoints.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)conPoints[i].x;
		p[0][1] = (float)conPoints[i].y;
	}
	//将图像变换到适合傅里叶变换的最佳尺寸
	cv::Mat paddedMat;
	int m = cv::getOptimalDFTSize(Point_Complex.rows);
	int n = cv::getOptimalDFTSize(Point_Complex.cols);

	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);

	//确定需要用来进行边界重构的描述子数目
	float Descriptor_UsingRatio = 0.0063;								//傅里叶描述子的使用比例
	int DescriptorNumber = cvRound(Descriptor_UsingRatio * m) & -2;	//傅里叶描述子的使用数，如果计算出来为奇数要变为偶数
	int half_RemoveNumber = (m - DescriptorNumber) / 2;				//每一端需要去除的描述子数量

	Point_complex_dft.rowRange(0, half_RemoveNumber) = cv::Scalar::all(0);	//删除系数(设置为0)
	Point_complex_dft.rowRange(Point_complex_dft.rows - half_RemoveNumber, Point_complex_dft.rows) = cv::Scalar::all(0);
	
	cv::Mat Point_complex_idft;
	cv::idft(Point_complex_dft, Point_complex_idft, cv::DFT_COMPLEX_INPUT | cv::DFT_COMPLEX_OUTPUT | cv::DFT_SCALE);

	//对反变换回来的图像去中心化(由于只有一列数据，x = 0)
	for (int y = 0; y < Point_complex_idft.rows; ++y)
	{
		if ((y % 2) != 0)
		{
			Point_complex_idft.ptr<cv::Vec2f>(y)[0][0] = -Point_complex_idft.ptr<cv::Vec2f>(y)[0][0];
			Point_complex_idft.ptr<cv::Vec2f>(y)[0][1] = -Point_complex_idft.ptr<cv::Vec2f>(y)[0][1];
		}
	}
	
	//裁剪回原来的图像大小得到结果
	cv::Mat result = Point_complex_idft.rowRange(0, Point_Complex.rows).clone();
	result.convertTo(result, CV_32SC2);

	std::vector<std::vector<cv::Point>> fixed_contours;
	std::vector<cv::Point> fixed_points;
	//将图像变回轮廓点并在画板上进行输出
	for (int y = 0; y < result.rows; ++y)
	{
		cv::Point dot;
		dot.x = result.ptr<cv::Vec2i>(y)[0][0];
		dot.y = result.ptr<cv::Vec2i>(y)[0][1];
		fixed_points.push_back(dot);
	}

	fixed_contours.push_back(fixed_points);
	cv::drawContours(drawmat, fixed_contours, -1, cv::Scalar(0, 255, 0), 1);

	cv::imshow("result", drawmat);
	cv::waitKey(0);
	return 0;
}

试验结果

通过调整傅里叶描述子使用比例Descriptor_UsingRatio来对比不同的边界重建结果。

原图：

重建效果图(分别对应0.63%，1.25%)：

描述子使用比例越少，损失的细节就越多，重建边界越平滑。
更多的细节请参考：《数字图像处理》

频率域滤波参考链接：OpenCV实现图像基础频率域滤波
形状特征提取：OpenCV实现傅里叶描述子(下)：形状特征提取