（五）OpenCV | 斑点中心检测（图像矩）

---

1. 图像矩

前一文介绍了斑点的检测，本文将介绍斑点中心的检测，主要用到了OpenCV中图像矩的概念。图像矩不仅可以描述图像的全局特征，还可以提供大量如图像大小、位置、方向和形状等图像信息。OpenCV中图像矩的计算定义在类Moments中，其部分源码为：

Moments(double m00, double m10, double m01, double m20, double m11,
		double m02, double m30, double m21, double m12, double m03);
// spatial moments
CV_PROP_RW double m00, m10, m01, m20, m11, m02, m30, m21, m12, m03;
// central moments
CV_PROP_RW double mu20, mu11, mu02, mu30, mu21, mu12, mu03;
// central normalized moments
CV_PROP_RW double nu20, nu11, nu02, nu30, nu21, nu12, nu03;

首先是空间矩。其中，$i+j$ 阶矩的计算公式如下：$${\mathrm{m}}_{ij}=\sum _{x,y}array(x,y)\cdot x^i\cdot y^j\text{\hspace{1em}(1)}$$

这里，对于二值图像来说，$array(x,y)$ 为 $0$ 或 $1$。则零阶矩 ${\mathrm{m}}_{00}$ 表示该二值图像的非零像素点数量，易得该值与轮廓所在位置无关；一阶矩 ${\mathrm{m}}_{01}$ 和 ${\mathrm{m}}_{10}$ 分别表示该二值图像的非零像素点横坐标和纵坐标之和，可以反映图像的形状特征。则此时该二值图像对应的质心坐标可以通过下式得到：$$\bar{x}=\frac{{\mathrm{m}}_{10}}{{\mathrm{m}}_{00}},\bar{y}=\frac{{\mathrm{m}}_{01}}{{\mathrm{m}}_{00}}\text{\hspace{1em}(2)}$$

同理，可以得到二阶矩和三阶矩，而二阶矩和三阶矩可以进一步得到图像的中心矩和归一化中心矩，统称为不变矩。不变矩可以反映图像的统计特征，其具有旋转、平移、缩放的不变特征。先来看中心矩，其中 $i+j$ 阶中心矩的计算公式如下：$${\mathrm{m}\mathrm{u}}_{ij}=\sum _{x,y}array(x,y)\cdot (x-\bar{x}{)}^i\cdot (y-\bar{y}{)}^j\text{\hspace{1em}(3)}$$

通过上式，易得零阶中心矩等于零阶矩 ${\mathrm{m}}_{00}$ ，一阶中心距为零。对于更高阶的矩，由于是相对物体质量而得到的，所以它们也不随物体在图像中的绝对位置改变而改变，即中心矩的平移不变性。再来看归一化中心矩，其中 $i+j$ 阶归一化中心矩的计算公式如下：$${\mathrm{n}\mathrm{u}}_{ij}=\frac{{\mathrm{m}\mathrm{u}}_{ij}}{{\mathrm{m}}_{00}^{(\mathrm{i}+\mathrm{j})/2+1}}\text{\hspace{1em}(4)}$$

归一化中心矩就是在中心矩加上一个归一化因子，这个因子是二值化图像中的非零像素点数量的幂（对于更高阶的矩，这个幂会变大）。中心矩通过减去均值而获得平移不变性，归一化中心矩通过除以物体的总尺寸而获得缩放不变性（当物体缩放时，上式的分子和分母以同等比例缩放）。同时，归一化中心矩也具有平移不变性。

最后，为了叙述的完整性，介绍Hu不变矩。它是归一化中心矩的线性组合。通过组合不同的归一化中心矩，我们可以得到一个反映图像不同特征的不同函数，这个函数不随平移、缩放、旋转等的变化而变化。
$$\begin{array}{rl}{\mathrm{h}\mathrm{u}}_1& ={\mathrm{n}\mathrm{u}}_{20}+{\mathrm{n}\mathrm{u}}_{02}\\ {\mathrm{h}\mathrm{u}}_2& =({\mathrm{n}\mathrm{u}}_{20}+{\mathrm{n}\mathrm{u}}_{02}{)}^2+4{\mathrm{n}\mathrm{u}}_{11}^2\\ {\mathrm{h}\mathrm{u}}_3& =({\mathrm{n}\mathrm{u}}_{30}-3{\mathrm{n}\mathrm{u}}_{12}{)}^2+3({\mathrm{n}\mathrm{u}}_{21}-{\mathrm{n}\mathrm{u}}_{03}{)}^2\\ {\mathrm{h}\mathrm{u}}_4& =({\mathrm{n}\mathrm{u}}_{30}+{\mathrm{n}\mathrm{u}}_{12}{)}^2+({\mathrm{n}\mathrm{u}}_{21}+{\mathrm{n}\mathrm{u}}_{03}{)}^2\\ {\mathrm{h}\mathrm{u}}_5& =({\mathrm{n}\mathrm{u}}_{30}-3{\mathrm{n}\mathrm{u}}_{12})({\mathrm{n}\mathrm{u}}_{30}+{\mathrm{n}\mathrm{u}}_{12})\left[({\mathrm{n}\mathrm{u}}_{30}+{\mathrm{n}\mathrm{u}}_{12}{)}^2-3({\mathrm{n}\mathrm{u}}_{21}+{\mathrm{n}\mathrm{u}}_{03}{)}^2\right]\\ & +(3{\mathrm{n}\mathrm{u}}_{21}-{\mathrm{n}\mathrm{u}}_{03})({\mathrm{n}\mathrm{u}}_{21}+{\mathrm{n}\mathrm{u}}_{03})\left[3({\mathrm{n}\mathrm{u}}_{30}+{\mathrm{n}\mathrm{u}}_{12})-({\mathrm{n}\mathrm{u}}_{21}+{\mathrm{n}\mathrm{u}}_{03})\right]\\ {\mathrm{h}\mathrm{u}}_6& =({\mathrm{n}\mathrm{u}}_{20}-{\mathrm{n}\mathrm{u}}_{02})\left[({\mathrm{n}\mathrm{u}}_{30}+{\mathrm{n}\mathrm{u}}_{12}{)}^2-({\mathrm{n}\mathrm{u}}_{21}+{\mathrm{n}\mathrm{u}}_{03}{)}^2\right]+4{\mathrm{n}\mathrm{u}}_{11}({\mathrm{n}\mathrm{u}}_{30}+{\mathrm{n}\mathrm{u}}_{12})({\mathrm{n}\mathrm{u}}_{21}+{\mathrm{n}\mathrm{u}}_{03}{)}^2\\ {\mathrm{h}\mathrm{u}}_7& =(3{\mathrm{n}\mathrm{u}}_{21}-{\mathrm{n}\mathrm{u}}_{03})({\mathrm{n}\mathrm{u}}_{30}+{\mathrm{n}\mathrm{u}}_{12})\left[({\mathrm{n}\mathrm{u}}_{30}+{\mathrm{n}\mathrm{u}}_{12}{)}^2-3({\mathrm{n}\mathrm{u}}_{21}+{\mathrm{n}\mathrm{u}}_{03}{)}^2\right]\\ & -({\mathrm{n}\mathrm{u}}_{30}-3{\mathrm{n}\mathrm{u}}_{12})({\mathrm{n}\mathrm{u}}_{21}+{\mathrm{n}\mathrm{u}}_{03})\left[3({\mathrm{n}\mathrm{u}}_{30}+{\mathrm{n}\mathrm{u}}_{12}{)}^2-({\mathrm{n}\mathrm{u}}_{21}+{\mathrm{n}\mathrm{u}}_{03}{)}^2\right]\end{array}\text{\hspace{1em}(5)}$$

这里，直接使用一个例子来说明Hu不变矩的平移、旋转、缩放等不变性特征。

#include<iostream>
#include<opencv2/opencv.hpp>

using namespace cv;
using namespace std;

int main(int argc, char** argv) {
	// 读入图像
	Mat gray;
	Mat image = imread("figures/moments_1.png");
	// 转换为灰度图
	cvtColor(image, gray, CV_BGR2GRAY);
	// 计算Hu不变矩
	Moments m = moments(gray);
	double h[7];
	HuMoments(m, h);
	// 打印
	for (int i = 0; i < 7; ++i) {
		cout << h[i] << endl;
	}
	return 0;
}

各输入图像：

第二幅图像是第一幅图像平移后的结果；第三幅图像是第一幅图像放大后的结果；第四幅图像是第一幅图像旋转后的结果。

	${\mathrm{h}\mathrm{u}}_1$	${\mathrm{h}\mathrm{u}}_2$	${\mathrm{h}\mathrm{u}}_3$	${\mathrm{h}\mathrm{u}}_4$	${\mathrm{h}\mathrm{u}}_5$	${\mathrm{h}\mathrm{u}}_6$	${\mathrm{h}\mathrm{u}}_7$
1	0.00130764	5.4248e-08	2.80186e-11	2.15264e-12	1.58235e-23	2.33778e-16	-5.39496e-24
2	0.00130711	5.43252e-08	2.75262e-11	2.14713e-12	1.55648e-23	2.29674e-16	-5.49622e-24
3	0.00130433	5.58127e-08	2.67299e-11	1.98447e-12	1.36036e-23	2.17734e-16	-4.88229e-24
4	0.00130697	5.41004e-08	2.76969e-11	2.03061e-12	1.44933e-23	2.253e-16	-4.67461e-24

由表可知，在误差范围内，Hu不变矩满足平移、旋转、缩放等特征的不变性。

最后，根据以上知识来求斑点的中心。首先是单个斑点的中心检测：

#include<stdio.h>
#include<stdlib.h>
#include<iostream>
#include<opencv2/highgui/highgui.hpp>
#include<opencv2/imgproc/imgproc.hpp>

using namespace cv;
using namespace std;

int main(int argc, char** argv) {
	// 定义变量存放读入图像、转换后的图像以及阈值
	Mat src, gray, thr;
	// 读入图像，第二个参数flag=1表示读入彩色图像
	src = imread("circle.png", 1);
	// 转换为灰度图
	cvtColor(src, gray, COLOR_BGR2GRAY);
	// 将像素值二值化，如果当前位置像素值大于100，则为255；否则为0
	threshold(gray, thr, 100, 255, THRESH_BINARY);
	// 得到二值图像的矩，第二个参数为true表示输入为二值图像，即所有非零像素点值为1
	Moments m = moments(thr, true);
	// 得到中心坐标，即质心坐标
	Point p(m.m10 / m.m00, m.m01 / m.m00);
	// 输出中心坐标
	cout << Mat(p) << endl;
	// 绘制中心坐标
	circle(src, p, 5, Scalar(128, 0, 0), -1);
	// 显示图像
	imshow("Image with Center", src);
	waitKey(0);

	return 0;
}

输入图像及输出结果：

---

2. Canny边缘检测

Canny是经典的边缘检测算法，其主要流程分为以下几步：

2.1 高斯平滑

由于边缘检测易受到图像中噪点的影响，Canny边缘检测第一步首先使用高斯平滑对图像进行降噪处理。大小为$(2k+1)\times (2k+1)$的高斯核为：$$H_{ij}=\frac{1}{2\pi {\sigma }^2}\exp \left(-\frac{(i-(k+1){)}^2+(j-(k+1){)}^2}{2{\sigma }^2}\right);1\le i,j\le (2k+1)\text{\hspace{1em}(6)}$$

这里，选择高斯核的大小将会影响Canny检测器的性能。一般情况下，高斯核的尺寸越大，检测器对噪点的敏感度越低；但随着高斯核尺寸的增大，边缘检测的误差会增大。

2.2 图像梯度

图像中的边缘可能指向任意方向，Canny检测算法使用四个过滤器来检测水平、垂直和对角线方向的边缘。使用边缘检测算子（OpenCV中使用Sobel算子）得到水平方向和垂直方向的梯度，分别为${\mathbf{G}}_x$和${\mathbf{G}}_y$。由此可以确定边缘的梯度和方向：$$\mathbf{G}=\sqrt{{\mathbf{G}}_x^2+{\mathbf{G}}_y^2}\text{\hspace{1em}(7)}$$

$$\Theta =\arctan ({\mathbf{G}}_y,{\mathbf{G}}_x)\text{\hspace{1em}(8)}$$

其中，各参数的含义如下：$${\mathbf{G}}_x=\left[\begin{array}{ccc}-1& 0& +1\\ -2& 0& +2\\ -1& 0& +1\end{array}\right]\ast \mathbf{I},{\mathbf{G}}_y=\left[\begin{array}{ccc}+1& +2& +1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right]\ast \mathbf{I}\text{\hspace{1em}(9)}$$

这里，$\Theta$取整为0°、45°、90°或135°，如 [0°, 22.5°] 和 [157.5°, 180°] 映射为0°。

2.3 非极大值抑制

仅仅基于梯度值提取的边缘仍然很模糊，非极大值抑制是一种边缘稀疏技术。具体地，将当前像素的梯度与沿正负梯度方向上的两个像素进行对比，如果当前像素的梯度比另外两个像素的梯度大，则保留该像素点为边缘点；否则将其抑制。如果某点像素$P$的梯度用红线表示：

为了简便，在最初的非极大值抑制中，采用量化的方式确定某个像素点的梯度：

• $\Theta \in [0,22.5)\cup (-22.5,0]\cup (157.5,180]\cup (-180,157.5]$，量化为水平方向0°的梯度
• $\Theta \in [22.5,67.5)\cup (-157.5,-112.5]$，量化为对角线方向45°的梯度
• $\Theta \in [67.5,112.5)\cup (-112.5,-67.5]$，量化为垂直方向方向90°的梯度
• $\Theta \in [112.5,157.5)\cup (-67.5,-22.5]$，量化为对角线方向135°的梯度

后来，为了得到更加精确的计算结果，在梯度跨越方向的两个相邻像素之间使用插值的方式来得到待比较像素的梯度。伪代码为：

tan_theta = G_y / G_x
# 正方向梯度
G_1 = (1 - tan_theta) * G_0 + tan_theta * G_45
# 负方向梯度
G_2 = (1 - tan_theta) * G_180 + tan_theta * G_225

2.4 双阈值

为了进一步强化边缘，Canny边缘检测引入双阈值法。具体地，如果该像素值的梯度大于高阈值，则视其为强边缘；如果该像素值介于低阈值与高阈值之间，则视其为弱边缘；否则将其抑制。

2.5 边缘跟踪

基于上一步的结果，对于强边缘的像素点将其视为最后的边缘。其次，如果弱边缘与强边缘有连接，也将其视为最后的边缘；否则将其抑制。这样就得到了最后物体的边缘，即Canny边缘检测整体流程。

基于Canny边缘检测提取多斑点中心：

#include<stdio.h>
#include<stdlib.h>
#include<iostream>
#include<opencv2/highgui/highgui.hpp>
#include<opencv2/imgproc/imgproc.hpp>

using namespace cv;
using namespace std;

RNG rng(12345);

void find_moments(Mat src);

int main(int argc, char** argv) {
	// 读入图像
	Mat src, gray;
	src = imread("figures/multiple-blob.png", 1);
	// 转换为灰度图
	cvtColor(src, gray, COLOR_BGR2GRAY);
	// 显示原图
	namedWindow("Source", WINDOW_AUTOSIZE);
	imshow("Source", src);
	// 调用函数寻找多个斑点的中心
	find_moments(gray);
	waitKey(0);

	return 0;
}

void find_moments(Mat gray) {
	// 用于存放Canny边缘检测后得到的边缘信息
	Mat canny_output;
	// contours用于存放所有集合，contours.size()为轮廓数量
	vector<vector<Point>> contours;
	// Vec4i表示四维int向量，hierarchy用于存放后一个轮廓、前一个轮廓、父轮廓、内嵌轮廓的索引
	// 如果不存在上述轮廓，则值为-1
	vector<Vec4i> hierarchy;
	// Canny边缘检测，第三/四个参数表示阈值，第五个参数表示Sobel算子大小
	Canny(gray, canny_output, 50, 150, 3);
	// 寻找轮廓，RETR_TREE提取所有轮廓并重建网状轮廓结构
	// CHAIN_APPROX_SIMPLE压缩水平方向、垂直方向、对角线方向的元素，只保留重点坐标，如矩阵只需要四个点
	// Point(0, 0)表示偏移坐标
	findContours(canny_output, contours, hierarchy, RETR_TREE, CHAIN_APPROX_SIMPLE, Point(0, 0));
	// 将每个轮廓定义为一个Moments对象，利用单个斑点中心检测方法检测所有斑点的中心
	vector<Moments> mu(contours.size());
	for (int i = 0; i < contours.size(); ++i) {
		mu[i] = moments(contours[i], false);
	}
	// 存储所有斑点的中心（质心坐标）
	vector<Point2f> mc(contours.size());
	for (int i = 0; i < contours.size(); ++i) {
		mc[i] = Point2f(mu[i].m10 / mu[i].m00, mu[i].m01 / mu[i].m00);
	}
	// 绘制
	Mat drawing(canny_output.size(), CV_8UC3, Scalar(255, 255, 255));
	for (int i = 0; i < contours.size(); ++i) {
		Scalar color = Scalar(167, 151, 0);
		drawContours(drawing, contours, i, color, 2, 8, hierarchy, 0, Point());
		circle(drawing, mc[i], 4, color, -1, 7, 0);
	}
	// 显示
	namedWindow("Contours", WINDOW_AUTOSIZE);
	imshow("Contours", drawing);
	waitKey(0);
}

输入图形及输出结果：

---

参考

https://github.com/spmallick/learnopencv/tree/master/CenterofBlob

---

完

---