在做物体检测时，由于成本和应用场合的限制，不能够一味地增加相机的分辨率，或者已经用了分辨率很高的相机，但是视野范围很大，仍然无法实现很高的精度，这时就要考虑亚像素技术，亚像素技术就是在两个像素点之间进行进一步的细分，从而得到亚像素级别的边缘点的坐标（也就是float类型的坐标），一般来说，现有的技术可以做到2细分、4细分，甚至很牛的能做到更高，通过亚像素边缘检测技术的使用，可以节约成本，提高识别精度。

    常见的亚像素技术包括灰度矩、Hu矩、空间矩、Zernike矩等，通过查阅多篇论文，综合比较精度、时间和实现难度，选择Zernike矩作为项目中的亚像素边缘检测算法。基于Zernike矩的边缘检测原理，下一篇文章详细再总结一下，这里大体把步骤说一下，并使用OpenCV实现。

    大体步骤就是，首先确定使用的模板大小，然后根据Zernike矩的公式求出各个模板系数，例如我选择的是7*7模板（模板越大精度越高，但是运算时间越长），要计算M00，M11R，M11I，M20，M31R，M31I，M40七个Zernike模板。第二步是对待检测图像进行预处理，包括滤波二值化等，也可以在进行一次Canny边缘检测。第三步是将预处理的图像跟7个Zernike矩模板分别进行卷积，得到7个Zernike矩。第四步是把上一步的矩乘上角度校正系数（这是因为利用Zernike的旋转不变性，Zernike模型把边缘都简化成了x=n的直线，这里要调整回来）。第五步是计算距离参数l和灰度差参数k，根据k和l判断该点是否为边缘点。以下是基于OpenCV的实现。

#include <math.h>
#include <opencv2/opencv.hpp>

using namespace std;
using namespace cv;

const double PI = 3.14159265358979323846;
const int g_N = 7;

Mat M00 = (Mat_<double>(7, 7) <<
	0, 0.0287, 0.0686, 0.0807, 0.0686, 0.0287, 0,
	0.0287, 0.0815, 0.0816, 0.0816, 0.0816, 0.0815, 0.0287,
	0.0686, 0.0816, 0.0816, 0.0816, 0.0816, 0.0816, 0.0686,
	0.0807, 0.0816, 0.0816, 0.0816, 0.0816, 0.0816, 0.0807,
	0.0686, 0.0816, 0.0816, 0.0816, 0.0816, 0.0816, 0.0686,
	0.0287, 0.0815, 0.0816, 0.0816, 0.0816, 0.0815, 0.0287,
	0, 0.0287, 0.0686, 0.0807, 0.0686, 0.0287, 0);

Mat M11R = (Mat_<double>(7, 7) <<
	0, -0.015, -0.019, 0, 0.019, 0.015, 0,
	-0.0224, -0.0466, -0.0233, 0, 0.0233, 0.0466, 0.0224,
	-0.0573, -0.0466, -0.0233, 0, 0.0233, 0.0466, 0.0573,
	-0.069, -0.0466, -0.0233, 0, 0.0233, 0.0466, 0.069,
	-0.0573, -0.0466, -0.0233, 0, 0.0233, 0.0466, 0.0573,
	-0.0224, -0.0466, -0.0233, 0, 0.0233, 0.0466, 0.0224,
	0, -0.015, -0.019, 0, 0.019, 0.015, 0);

Mat M11I = (Mat_<double>(7, 7) <<
	0, -0.0224, -0.0573, -0.069, -0.0573, -0.0224, 0,
	-0.015, -0.0466, -0.0466, -0.0466, -0.0466, -0.0466, -0.015,
	-0.019, -0.0233, -0.0233, -0.0233, -0.0233, -0.0233, -0.019,
	0, 0, 0, 0, 0, 0, 0,
	0.019, 0.0233, 0.0233, 0.0233, 0.0233, 0.0233, 0.019,
	0.015, 0.0466, 0.0466, 0.0466, 0.0466, 0.0466, 0.015,
	0, 0.0224, 0.0573, 0.069, 0.0573, 0.0224, 0);

Mat M20 = (Mat_<double>(7, 7) <<
	0, 0.0225, 0.0394, 0.0396, 0.0394, 0.0225, 0,
	0.0225, 0.0271, -0.0128, -0.0261, -0.0128, 0.0271, 0.0225,
	0.0394, -0.0128, -0.0528, -0.0661, -0.0528, -0.0128, 0.0394,
	0.0396, -0.0261, -0.0661, -0.0794, -0.0661, -0.0261, 0.0396,
	0.0394, -0.0128, -0.0528, -0.0661, -0.0528, -0.0128, 0.0394,
	0.0225, 0.0271, -0.0128, -0.0261, -0.0128, 0.0271, 0.0225,
	0, 0.0225, 0.0394, 0.0396, 0.0394, 0.0225, 0);

Mat M31R = (Mat_<double>(7, 7) <<
	0, -0.0103, -0.0073, 0, 0.0073, 0.0103, 0,
	-0.0153, -0.0018, 0.0162, 0, -0.0162, 0.0018, 0.0153,
	-0.0223, 0.0324, 0.0333, 0, -0.0333, -0.0324, 0.0223,
	-0.0190, 0.0438, 0.0390, 0, -0.0390, -0.0438, 0.0190,
	-0.0223, 0.0324, 0.0333, 0, -0.0333, -0.0324, 0.0223,
	-0.0153, -0.0018, 0.0162, 0, -0.0162, 0.0018, 0.0153,
	0, -0.0103, -0.0073, 0, 0.0073, 0.0103, 0);

Mat M31I = (Mat_<double>(7, 7) <<
	0, -0.0153, -0.0223, -0.019, -0.0223, -0.0153, 0,
	-0.0103, -0.0018, 0.0324, 0.0438, 0.0324, -0.0018, -0.0103,
	-0.0073, 0.0162, 0.0333, 0.039, 0.0333, 0.0162, -0.0073,
	0, 0, 0, 0, 0, 0, 0,
	0.0073, -0.0162, -0.0333, -0.039, -0.0333, -0.0162, 0.0073,
	0.0103, 0.0018, -0.0324, -0.0438, -0.0324, 0.0018, 0.0103,
	0, 0.0153, 0.0223, 0.0190, 0.0223, 0.0153, 0);

Mat M40 = (Mat_<double>(7, 7) <<
	0, 0.013, 0.0056, -0.0018, 0.0056, 0.013, 0,
	0.0130, -0.0186, -0.0323, -0.0239, -0.0323, -0.0186, 0.0130,
	0.0056, -0.0323, 0.0125, 0.0406, 0.0125, -0.0323, 0.0056,
	-0.0018, -0.0239, 0.0406, 0.0751, 0.0406, -0.0239, -0.0018,
	0.0056, -0.0323, 0.0125, 0.0406, 0.0125, -0.0323, 0.0056,
	0.0130, -0.0186, -0.0323, -0.0239, -0.0323, -0.0186, 0.0130,
	0, 0.013, 0.0056, -0.0018, 0.0056, 0.013, 0);



int main()
{

	Mat OriginalImage = imread("original.png", 0);
	Mat SubImage = OriginalImage;

	Mat NewSmoothImage;
	medianBlur(OriginalImage, NewSmoothImage, 13);

	Mat NewAdaThresImage;
	adaptiveThreshold(NewSmoothImage, NewAdaThresImage, 255, ADAPTIVE_THRESH_GAUSSIAN_C, THRESH_BINARY_INV, 7, 4);

	vector<Point2d> SubEdgePoints;

	Mat ZerImageM00;
	filter2D(NewAdaThresImage, ZerImageM00, CV_64F, M00, Point(-1, -1), 0, BORDER_DEFAULT);
	//filter2D( cannyImage, zerImageM00, CV_64F, M00, Point(-1,-1), 0, BORDER_DEFAULT);

	Mat ZerImageM11R;
	filter2D(NewAdaThresImage, ZerImageM11R, CV_64F, M11R, Point(-1, -1), 0, BORDER_DEFAULT);
	//filter2D( cannyImage, zerImageM11R, CV_64F, M11R, Point(-1, -1), 0, BORDER_DEFAULT);

	Mat ZerImageM11I;
	filter2D(NewAdaThresImage, ZerImageM11I, CV_64F, M11I, Point(-1, -1), 0, BORDER_DEFAULT);
	//filter2D( cannyImage, zerImageM11I, CV_64F, M11I, Point(-1, -1), 0, BORDER_DEFAULT);

	Mat ZerImageM20;
	filter2D(NewAdaThresImage, ZerImageM20, CV_64F, M20, Point(-1, -1), 0, BORDER_DEFAULT);
	//filter2D( cannyImage, zerImageM20, CV_64F, M20, Point(-1, -1), 0, BORDER_DEFAULT);

	Mat ZerImageM31R;
	filter2D(NewAdaThresImage, ZerImageM31R, CV_64F, M31R, Point(-1, -1), 0, BORDER_DEFAULT);
	//filter2D(cannyImage, zerImageM31R, CV_64F, M31R, Point(-1, -1), 0, BORDER_DEFAULT);

	Mat ZerImageM31I;
	filter2D(NewAdaThresImage, ZerImageM31I, CV_64F, M31I, Point(-1, -1), 0, BORDER_DEFAULT);
	//filter2D(cannyImage, zerImageM31I, CV_64F, M31I, Point(-1, -1), 0, BORDER_DEFAULT);

	Mat ZerImageM40;
	filter2D(NewAdaThresImage, ZerImageM40, CV_64F, M40, Point(-1, -1), 0, BORDER_DEFAULT);
	//filter2D(cannyImage, zerImageM40, CV_64F, M40, Point(-1, -1), 0, BORDER_DEFAULT);



	int row_number = NewAdaThresImage.rows;
	int col_number = NewAdaThresImage.cols;
	//使用7个的7*7的Zernike模板（其本质是个矩阵）M00、M11R、M11I、M20、M31R、M31I、M40，分别与图像进行卷积，得到每个像素点的七个Zernike矩Z00、Z11R、Z11I、Z20、Z31R、Z31I、Z40
	//对于每个点，根据它的七个Zernike矩，求得距离参数l和灰度差参数k，当l和k都满足设定的条件时，则判读该点为边缘点，并进一步利用上述7个Zernike矩求出该点的亚像素级坐标
	//如果l或k不满足设定的条件，则该点不是边缘点，转到下一个点求解距离参数l和灰度差参数k
	for (int i = 0; i < row_number; i++)
	{
		for (int j = 0; j <col_number; j++)
		{
			if (ZerImageM00.at<double>(i, j) == 0)
			{
				continue;
			}

			//compute theta
			//vector<vector<double> > theta2(0);
			double theta_temporary = atan2(ZerImageM31I.at<double>(i, j), ZerImageM31R.at<double>(i, j));
			//theta2[i].push_back(thetaTem);

			//compute Z11'/Z31'
			double rotated_z11 = 0.0;
			rotated_z11 = sin(theta_temporary)*(ZerImageM11I.at<double>(i, j)) + cos(theta_temporary)*(ZerImageM11R.at<double>(i, j));
			double rotated_z31 = 0.0;
			rotated_z31 = sin(theta_temporary)*(ZerImageM31I.at<double>(i, j)) + cos(theta_temporary)*(ZerImageM31R.at<double>(i, j));

			//compute l
			double l_method1 = sqrt((5 * ZerImageM40.at<double>(i, j) + 3 * ZerImageM20.at<double>(i, j)) / (8 * ZerImageM20.at<double>(i, j)));
			double l_method2 = sqrt((5 * rotated_z31 + rotated_z11) / (6 * rotated_z11));
			double l = (l_method1 + l_method2) / 2;
			//compute k/h
			double k, h;
		
			k = 3 * rotated_z11 / 2 / pow((1 - l_method2*l_method2), 1.5);
			h = (ZerImageM00.at<double>(i, j) - k*PI / 2 + k*asin(l_method2) + k*l_method2*sqrt(1 - l_method2*l_method2)) / PI;

			//judge the edge
			double k_value = 20.0;
		
			double l_value = sqrt(2) / g_N;
		
			double absl = abs(l_method2 - l_method1);
			if (k >= k_value && absl <= l_value)
			{
				Point2d point_temporary;
				point_temporary.x = j + g_N*l*cos(theta_temporary) / 2;
				point_temporary.y = i + g_N*l*sin(theta_temporary) / 2;
				SubEdgePoints.push_back(point_temporary);
			}
			else
			{
				continue;
			}
		}
	}
	//显示所检测到的亚像素边缘
	for (size_t i = 0; i < SubEdgePoints.size(); i++)
	{
	    Point center_forshow(cvRound(SubEdgePoints[i].x), cvRound(SubEdgePoints[i].y));
		circle(OriginalImage, center_forshow, 1, Scalar(0, 97, 255), 1, 8, 0);
	}
	imshow("亚像素边缘", OriginalImage);
	
	waitKey(0);
	return 0;
}

---------------------------------------7.28更新分割线-------------------------------------------------

1.因为留言要源代码的朋友比较多，所以我把上面的代码补完整了，新的代码能够实现基本的过程，但是针对不同的应用，肯定不是最优的效果，可能需要根据不同的应用场合进行调整；

2.评论里提到的双层轮廓，我能想到的是因为对线进行边缘检测操作，因为线的两侧相当于都是前景和背景的边缘，所以两侧各会有一层边缘，不管是Canny（事实上不少Zernike论文都建议先用canny再用Zernike）还是二值化，只要送入Zernike的是线而不是区域，就会有双层轮廓，其他的原因我暂时还没想到，欢迎大家赐教；

3.双层轮廓并不一定是坏事，比如做圆检测，最后是要用这些边缘点来拟合圆，那么两层轮廓拟合出来的圆放大了说应该是处于两层轮廓中间，这样的结果其实更符合实际情况，但是如果不是规则轮廓，可能会对边缘点的表达和拟合有些问题；

声明：本篇仅仅是分享网上的开源项目，算法非本人原创。转载文章：

《A Sub-Pixel Edge Detector: an Implementation of the Canny/Devernay Algorithm》_Naruto_Q的博客-CSDN博客

整个项目，用opencv包装了，项目地址：

https://github.com/fcqing/sub-pixel-edge-detect

https://github.com/CsCsongor/subPixelEdgeDetect

GitHub - songyuncen/EdgesSubPix: Sub-Pixel Edge Detection Using OpenCV

论文《A Sub-Pixel Edge Detector: an Implementation of theCanny/Devernay Algorithm》提出了一种基于canny算法进行亚像素级边缘提取的算法，并给出了网页测试demo，可以在上面载入图像进行边缘提取测试，效果不错。这里简要翻译一下论文的方法：

Canny算法

针对边缘检测算法的研究有很多，其中Canny算法在这方面具有开创性的贡献，我在之前的博客里总结过canny的原理，这里简要说明一下：

Canny算子求边缘点具体算法步骤如下：
1. 用高斯滤波器平滑图像．
2. 计算梯度幅值和方向

3. 对梯度幅值进行非极大值抑制

非极大值抑制（Non-Maxima Suppression method）主要是为了更精确的定位某种特征，比如用梯度变化表征边缘时，梯度变化较大的区域通常比较宽，所以利用x和y方向的梯度确定一个法向arctan(y/x)，然后在法向上判断当前梯度测量是否是一个峰值（或局部极大值），如果是就保留，不是极大值就抑制（如设置为0）。这样的话就能将边缘定位在1-2像素宽（相邻像素有时候求极大值的方向恰好互不干扰）。

4. 用双阈值算法检测和连接边缘．

5.通过滞后跟踪边缘，抑制所有弱而不强的边缘。

Devernay亚像素校正算法

Canny检测出的边缘是像素级别的，而有时需要精度低于一个像素，即亚像素边缘提取，Devernay【3】在Canny算法的基础上进行了拓展：将新的边缘点定义为相邻几个梯度模值差值的最大值，可以通过计算梯度方向上相邻的三点处的梯度模值的二次函数插补。如下图所示：

上图中，||g(A,B,C)||为三个垂直于边缘方向上的梯度模值，canny算法会选择模值最大的B点作为边缘点。但是，也可能在A和C之间存在点η ，在η 点处的梯度模值大于B点及AC之间的其他点，此时，η 点能更好地表示边缘。论文[3]给出了估计η 点以及计算η 点偏移量的方法，并给出η 的亚像素位置。

为了减少计算量，论文中只用了三点来进行二次方程拟合，并求出相对于BC向量的补偿 η ：

Devernay亚像素算法准确性

影响Devernay算法准确性的四个误差来源：

1.假设梯度满足二次曲线规律带来的误差；

2.使用有限项计算梯度值带来的误差；

3.插补计算时有些点的梯度并没有通过canny运算计算出来所带来的误差；（个人理解这里应该是说在非水平和垂直方向上）

4.数据误差；

以高斯边缘为例，公式：，其中是高斯核函数，Y是阶梯函数，即：

   (1)

α 是阶梯函数的系数，下图中左是带亚像素偏移 γ 的高斯边缘特性，右图是使用抛物线估计的边缘点位置，例如途中三点abc,:(−1, a), (0, b) 和(1, c) ，边缘估计点的位置为：

          (2)

把带偏移 γ 的高斯边缘函数带入，可以得到：

   (3)

带入（2）式得到：

            (4)

右图中，高斯梯度与二次方程梯度极大值点偏差为：

 ,

文中给出了该误差的分布图5：

当γ=0时，误差为0，因为此时边缘在网格中间位置；当γ=0.5时，误差也为0，因为此时a=b或b=c，二次曲线插补值可以给出正确的偏移；σ越小，边缘越锐利，而σ越大，边缘越模糊，误差反而越小。

公式（4）也说明了第二种误差来源，即有限差分计算带来的误差。

亚像素边缘算法改进

在上面说的四种误差来源中，其中简单地使用水平和垂直方向插补可以有效去除第三种误差。Deverary使用三个点中间的插补值作为边缘， 这里进行修正：当点（x，y）处梯度满足以下条件，此时该点为水平方向的边缘点：

 , 

 ,

类似的，垂直方向的边缘点也是可以定义出来的。当x和y方向的梯度相等时，默认为还是水平边缘点，这个影响不大。修正方法的误差分布见下图6：

对比图5，修正方法的误差有了明显的减弱。另外，论文还对比了边缘是倾斜角度时的提取结果，见下图7：

边缘点链（Edge Point Chaining）

此前检测边缘点是一个独立过程，每个边缘点并不相关。后续需要将属于同一个边缘轮廓的点集归到一起--链，形成轮廓。每个边缘点对应的应是一个像素点，而不是一个5*5的区域。

首先，被归为同一链路的像素点应该具有近似的梯度方向，以像素点A和B为例，它们的夹角应小于90度，数学表示为：g(A).g(B)>0。 另外，图像轮廓可以将亮区域和暗区域分开，所以连续的链路需要将暗区域划分到曲线的同一侧，一个简单的方法是验证从边缘点A到点B的向量与点A两个可能的梯度方向中的一个是否近似正交。有两个概念需要介绍一下：前向链路（forward chaining）和后向链路(back chaining)。当 ,其中表示向量v旋转90度，此时成为A-B的前向链路，记为:，反方向称为反向链路。

一个边缘点至多有一个前向链路，也至多有一个后向链路，而不用有类似 Y 形状的链路。所以在确定AB是否成链前，需要检查A是否已经有前向链路和B是否已经有后向链路。此时可以通过距离来判断，如果新的链路比之前的距离短，那么之前链路就会被斩断，而增加这新的链路，经过这个步骤，边缘点的链接关系就能确定下来了。以下图8为例：

A,B,C,D四点的距离关系满足：||A-B||>||B-C||>||C-D||,假定链路搜索顺序为：A-B-C,你会发现当我们想时,存在，而且根据CB距离小于AB，所以会被切断；同样地，由于D的存在，也会被切断，所以最后的结果只剩下：。 再假定链路搜索顺序为C-B-A，经过判断过程，剩下这条链路。 所以，最后的链路结果为：和。

算法实现

Canny/Devernay算法主要分为四部分：图像求梯度，计算边缘点，边缘点链路编码，应用canny双阈值等方法，如下图9所示：

计算图像梯度的方法为算法2，如下图10所示，主要在水平和垂直两个方向上求取梯度：

求取边缘点及亚像素校正过程，算法3的伪代码如下图11所示：

可以发现，当某一点的水平方向和垂直方向梯度值相等时，该点可以被归为任意方向上的边缘，不影响结果。最后，边缘点链路编码过程算法4，如下图12所示：

算法4中的是指边缘点的邻域集合，取5*5，括号里的数字2应该是说邻域集中上下左右不超过2个像素，最后的边缘点取其邻域水平或垂直方向上的局部最大值。当然，过程中包含前向链路和后向链路，算是的子集。

最后，算法5描述了canny算法的应用，如下图13所示：

测试结果

算法中的输入参数除了被检测图像外，还有高斯函数的标准差S，canny算法的双阈值H和L，在线测试网址请点击

可见，当S取不同值是，边缘的粗细间隔是不一样的。S越大，边缘越粗些。

后记：第一次把英文文献翻译在CSDN上面，感慨大略看懂到能够翻译出来还是有一段距离的。

整个项目，我重新用opencv包装了，项目地址：

GitHub - fcqing/sub-pixel-edge-detect: this project is for searching an image's edge

GitHub - CsCsongor/subPixelEdgeDetect: An image edge detector is described which produces chained edge points with sub-pixel accuracy.The method incorporates the main ideas of the classic Canny and Devernay algorithms. Theanalysis shows that a slight modification to the original formulation improves the sub-pixelaccuracy of the edge points.

https://github.com/songyuncen/EdgesSubPix

参考

https://en.wikipedia.org/wiki/Canny_edge_detector

http://blog.csdn.net/qq_14845119/article/details/52064928

《A non-maxima suppression method for edge detection with sub-pixel accuracy》

http://ipolcore.ipol.im/demo/clientApp/demo.html?id=216