边缘检测算法的C++实现
基于sobel算子的边缘检测算法的C++实现及改进
实验内容
- 对已有的sobel算子边缘检测算法的学习
- 对已有算法进行C++实现
- 对已有算法进行改进
已有的边缘检测算法流程
- 对图像进行灰度化得到灰度图像
- 对灰度图像进行高斯模糊去除部分噪声
- 利用sobel算子作为卷积核对图像进行卷积,计算梯度
- 用非极大值抑制法进行边缘细化
- 利用双阈值进行噪声去除
算法实现
图像灰度化
cv::Mat Graying(cv::Mat img) {
int row = img.rows;
int col = img.cols;
// cv::imshow("灰度图", img);
cv::Mat grayImg = Mat(row, col, CV_8U);
cout << row << '\n' << col;
for(register int i = 0; i < row; ++i) {
Vec3b *p = img.ptr<Vec3b>(i);
uchar *p2 = grayImg.ptr<uchar>(i);
for(register int j = 0; j <=col; ++j) {
*p2 = (*p)[0] * 0.114 + (*p)[1] * 0.587 + (*p)[2] * 0.299;
//Gray = 0.2989*R+0.5870*G+0.1140*B 注意opencv的Mat类型是BGR顺序
p2++;
p++;
}
}
return grayImg;
}
灰度化实现效果如下:
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191217072013outimg.jpg)
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191217072000outgrayimg.jpg)
边缘检测主体实现
高斯模糊
GaussianBlur(grayImg, grayImg, Size(5, 5), BORDER_DEFAULT); //进行5*5的高斯模糊(这个是最后懒了,觉得代码只是繁琐就不想手写,调的库
sobel算子滤波
sobel算子分为x方向和y方向,在某个像素点处分别用两个算子卷积求出x方向和y方向的导数值,再通过加权平均可以求出在某个像素点处的梯度值,而由于图像边缘就是图像灰度值发生剧烈变化的地方,梯度值大的地方是边缘的可能就大,因此sobel算子可以检测出边缘
\(G_x = \left[ \begin{matrix} -1 & 0 & 1\\-2 & 0 & 2\\-1 & 0 & 1 \end{matrix} \right] * A\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,G_y = \left[ \begin{matrix} -1 & -2 & -1\\0 & 0 & 0\\1 & 2 & 1 \end{matrix} \right] * A\\\\G = \sqrt{{G_x}^2+{G_y}^2}\)
cv::Mat filterlengthways(cv::Mat img1) {
cv::Mat img;
img1.copyTo(img);
int row = img.rows;
int col = img.cols;
vector<int> temp[row];
for (int i = 0; i < row; ++i) {
uchar *p = img.ptr<uchar>(i);
for (int j = 0; j < col; ++j) {
temp[i].push_back((int) (*p));
p++;
}
}
cv::Mat zero = Mat(cv::Size(col, row), CV_8UC1, Scalar(0));
zero.copyTo(img);
for (register int i = 1; i < row - 1; ++i) { //竖直方向滤波,算子为
for (register int j = 1; j < col - 1; ++j) { //-1 0 1
int temp1 = 0; //-2 0 2
temp1 -= temp[i - 1][j - 1]; //-1 0 1 该算子可以求出水平方向的导数
temp1 -= temp[i][j - 1] * 2;
temp1 -= temp[i + 1][j - 1];
temp1 += temp[i - 1][j + 1];
temp1 += temp[i][j + 1] * 2;
temp1 += temp[i + 1][j + 1];
if (temp1 < 0) temp1 = -temp1;
if (temp1 > 255) temp1 = 255;
img.at<uchar>(i, j) = temp1;
}
}
return img;
}
cv::Mat filtercrossways(cv::Mat img1) {
cv::Mat img;
img1.copyTo(img);
int row = img.rows;
int col = img.cols;
vector<int> temp[row];
for (int i = 0; i < row; ++i) {
uchar *p = img.ptr<uchar>(i);
for (int j = 0; j < col; ++j) {
temp[i].push_back((int) (*p));
p++;
}
}
cv::Mat zero = Mat(cv::Size(col, row), CV_8UC1, Scalar(0));
zero.copyTo(img);
for (register int i = 1; i < row - 1; ++i) {
for (register int j = 1; j < col - 1; ++j) { //水平方向滤波, 算子为
int temp1 = 0; //-1 -2 -1
temp1 -= temp[i - 1][j - 1]; // 0 0 0
temp1 -= temp[i - 1][j] * 2; // 1 2 1 求出竖直方向的导数
temp1 -= temp[i - 1][j + 1];
temp1 += temp[i + 1][j - 1];
temp1 += temp[i + 1][j] * 2;
temp1 += temp[i + 1][j + 1];
if (temp1 < 0) temp1 = -temp1;
if (temp1 > 255) temp1 = 255;
img.at<uchar>(i, j) = temp1;
}
// puts("");
}
return img;
}
vector<double> gradient;//求出梯度方向和大小,并进行梯度方向的规范化
for (int i = 0; i < row; ++i) {
for (int j = 0; j < col; ++j) {
int a = img1.at<uchar>(i, j), b = img2.at<uchar>(i, j);
int temp = sqrt(a * a + b * b);//几何平均求出梯度大小
gradient.push_back(atan2(b, a));
if (gradient[getNum(i, j, col)] < 0) gradient[getNum(i, j, col)] += 2 * pi;//使得梯度方向为0到2pi,方便处理
img.at<uchar>(i, j) = temp;
}
}
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191217072024outzong.jpg)
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191217072007outheng.jpg)
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191217071936outbeforenms.jpg)
非极大值抑制
由于sobel算子计算出的边缘过于模糊和粗大,因此需要进行边缘的细化。
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191217084155d.png)
不妨假设图像边缘部分的梯度如图,且之前的sobel算子滤波会将橙线之上的部分都当成边缘,因此边缘会粗大而模糊
于是就产生一种想法即只保留顶点所对应的像素点作为边缘,其它部分置为0,这就产生了非极大值抑制
之前sobel算子滤波时计算出了每个像素点梯度的大小和方向
因为图像中方格是离散的,因此对于梯度的方向进行区域的划分,如图所示
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191217085709wsl.png)
对于每个像素点,计算出梯度大小和梯度方向后,将其梯度大小与对应区域相邻像素的梯度大小进行比较,如果该像素点的梯度大小比任意一个对应区域相邻像素小,则该像素点置为0,反之保持原样,这样就实现了边缘的细化
const double pi = acos(-1);
const double zone[8] = {atan2(1, 3), atan2(3, 1), atan2(3, -1), atan2(1, -3), 2 * pi + atan2(-1, -3), 2 * pi + atan2(-3, -1), 2 * pi + atan2(-3, 1), 2 * pi + atan2(-1, 3)};
//进行区域划分,实际上是4个,但为了编码方便,化成8个
//注意之前求梯度方向时对atan2的结果都加上了2*pi
int getType(double t) { //返回梯度方向所属区域
for (int i = 0; i < 8; ++i) {
if (t >= zone[7] || t <= zone[0]) return 0;
if (t >= zone[0] && t <= zone[1]) return 1;
if (t >= zone[1] && t <= zone[2]) return 2;
if (t >= zone[2] && t <= zone[3]) return 3;
if (t >= zone[3] && t <= zone[4]) return 0;
if (t >= zone[4] && t <= zone[5]) return 1;
if (t >= zone[5] && t <= zone[6]) return 2;
if (t >= zone[6] && t <= zone[7]) return 3;
}
throw "error";
}
void nonMaximumSuppression(Mat img, vector<double> &gradient) {
int col = img.cols, row = img.rows;
for (int i = 0; i < row; ++i) {
for (int j = 0; j < col; ++j) {
switch (getType(gradient[getNum(i, j, col)])) {
case 0:
if (check(i, j - 1, row, col) &&
img.at<uchar>(i, j - 1) > img.at<uchar>(i, j)) {
img.at<uchar>(i, j) = 0;
}
if (check(i, j + 1, row, col) &&
img.at<uchar>(i, j + 1) > img.at<uchar>(i, j)) {
img.at<uchar>(i, j) = 0;
}
break;
case 1:
if (check(i - 1, j + 1, row, col) &&
img.at<uchar>(i - 1, j + 1) > img.at<uchar>(i, j)) {
img.at<uchar>(i, j) = 0;
}
if (check(i + 1, j - 1, row, col) &&
img.at<uchar>(i + 1, j - 1) > img.at<uchar>(i, j)) {
img.at<uchar>(i, j) = 0;
}
break;
case 2:
if (check(i - 1, j, row, col) &&
img.at<uchar>(i - 1, j) > img.at<uchar>(i, j)) {
img.at<uchar>(i, j) = 0;
}
if (check(i + 1, j, row, col) &&
img.at<uchar>(i + 1, j) > img.at<uchar>(i, j)) {
img.at<uchar>(i, j) = 0;
}
break;
case 3:
if (check(i + 1, j + 1, row, col) &&
img.at<uchar>(i + 1, j + 1) > img.at<uchar>(i, j)) {
img.at<uchar>(i, j) = 0;
}
if (check(i - 1, j - 1, row, col) &&
img.at<uchar>(i - 1, j - 1) > img.at<uchar>(i, j)) {
img.at<uchar>(i, j) = 0;
}
break;
}
}
}
}
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191217071930outbeforedtd.jpg)
双阈值判定
非极大值抑制后的图像仍存在许多噪声,可以将灰度值小于某个固定阈值的像素点置为0,可以一定程度上消除噪声
void imageBinaryzation(Mat img, int thresholdVal = 0) { //尝试对边缘进行单固定阈值,但效果不好
Mat img1;
img.copyTo(img1);
Mat_<uchar>::iterator it1 = img1.begin<uchar>();
Mat_<uchar>::iterator itend1 = img1.end<uchar>();
for (; it1 != itend1; ++it1) {
(*it1) > thresholdVal ? 0 : (*it1) = 0; //如果像素点灰度值大于阈值,灰度不变,否则灰度置为0;
}
OutPut("Binaryzation", img1);
}
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191217091052outBinaryzation.jpg)
然而固定单阈值法有许多问题
- 对于不同类型的图像,不同要求的边缘检测,阈值不同,得到效果好的阈值很麻烦
- 极容易造成边缘的中断
因此采取双阈值法可以一定程度上解决问题
我们定义高阈值和低阈值,对于灰度值高于高阈值的像素点我们将其确定为边缘,而灰度值低于低阈值的像素点我们认为它一定不是边缘,对于灰度值介于高阈值和低阈值之间的像素点,它有可能是边缘,当它相邻8个像素点存在确定的边缘时,它也会变成确定的边缘,否则它被判定为不是边缘它的灰度值将被置为0
void doubleThresholdDetection(Mat &img1) {
int col = img1.cols, row = img1.rows;
cv::Mat temp;
img1.copyTo(temp);
for (int i = 0; i < row; ++i) {
for (int j = 0; j < col; ++j) {
if (temp.at<uchar>(i, j) > highTh) {
continue;
} else {
if (temp.at<uchar>(i, j) < lowTh) {
img1.at<uchar>(i, j) = 0;
continue;
} else {
if (lowthcheck(temp, i - 1, j - 1, row, col) ||
lowthcheck(temp, i - 1, j, row, col) ||
lowthcheck(temp, i - 1, j + 1, row, col) ||
lowthcheck(temp, i, j - 1, row, col) ||
// lowthcheck(temp, i, j , row, col) ||
lowthcheck(temp, i, j + 1, row, col) ||
lowthcheck(temp, i + 1, j - 1, row, col) ||
lowthcheck(temp, i + 1, j, row, col) ||
lowthcheck(temp, i + 1, j + 1, row, col)
) {
temp.at<uchar>(i, j) = highTh + 1;
continue;
} else {
img1.at<uchar>(i, j) = 0;
}
}
}
}
}
}
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191217071950outfilter.jpg)
可以很明显的发现,双阈值法较好地解决了单阈值法会造成边缘中断的问题,并且由于介于高低阈值之间的点仍能被保留,因此检测效果较好的阈值范围变大,阈值更容易调整
改进的边缘细化方法
由于非极大值抑制对边缘进行细化时,只会保留梯度最大的一个点,因此边缘大多被细化成了一个像素点宽度的线,并且有可能使较粗的边缘变成两条线,这就使得边缘锯齿化严重结果不够美观,因此我试图用另一种思路进行边缘细化
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191218091609outautoBinaryzation6.jpg)
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191217071930outbeforedtd.jpg)
\(cntCrossWay\) 和 \(cntLengthWay\) 分别为上述两个矩阵的系数和
\(averCrossWay = {sumCrossWay \over cntCrossWay}\)
\(averLengthWay = {sumLengthWay \over cntLengthWay}\)
\(threshold Value = 0.92 *\sqrt{averLengthWay ^2+averCrossWay^2}\)
代码实现:
void autoImageBinaryzation6(Mat img) {
int row = img.rows, col = img.cols;
Mat img1;
img.copyTo(img1);
int averCrossWay, averLengthWay, sumCrossWay, sumLengthWay, cntCrossWay, cntLengthWay, aver;
for(int j = 0; j < col; ++j) {
for(int i = 0; i < row; ++i) {
sumCrossWay = sumLengthWay = cntCrossWay = cntLengthWay = 0;
for(int k = -4; k <= 4; ++k) {
for(int g = (k <= 1 ? -2 : ((k <= 3 ) ? -1 : 0)) ; g <= (k <= 1 ? 2 : ((k <=3) ? 1 : 0)); ++g) {
if(check(i + g, j + k, row, col)) {
sumLengthWay += img.at<uchar>(i + g, j + k) * (1 + (abs(k) >= 2));
cntLengthWay += (1 + (abs(k) >= 2));
}
if(check(i + k, j + g, row, col)) {
sumCrossWay += img.at<uchar>(i + k, j + g) * (1 + (abs(k) >= 2));
cntCrossWay += (1 + (abs(k) >= 2));
}
}
}
averLengthWay = sumLengthWay / cntLengthWay;
averCrossWay = sumCrossWay / cntCrossWay;
aver = sqrt(averCrossWay * averCrossWay + averLengthWay * averLengthWay);
img1.at<uchar>(i, j) > aver * 0.92 ? 0 : img1.at<uchar>(i, j) = 0;
}
}
OutPut("autoBinaryzation6", img1);
doubleThresholdDetection(img1);
OutPut("autoBinaryzation6dtd", img1);
OutPut("myfilter", img1);
}
该方法看起来时间复杂度略有增大,但averLengthWay和averCrossWay可以进行动态维护,不需要每次进行遍历求值,代码实现中直接for循环遍历求值只是为了呈现更加直观。 且只涉及整型操作,因此实际运行速度并不比涉及到浮点运算以及atan2操作的非极大值抑制有劣势。
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191218092510outmyfilter.jpg)
效果展示
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191217071936outbeforenms.jpg)
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191218091609outautoBinaryzation6.jpg)
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191217071930outbeforedtd.jpg)
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191218092510outmyfilter.jpg)
![](https://images.cnblogs.com/cnblogs_com/thhyj/1615217/o_191217071950outfilter.jpg)