步骤
- 用高斯滤波器平滑图像;
- 用一阶偏导有限差分计算梯度幅值和方向;
- 对梯度幅值应用非极大值抑制;
- 用双阈值算法检测和连接边缘。
一、openCV在C++中的应用
首先,在mac的Xcode上安装配置openCV库,参考一下链接(科学上网访问)https://medium.com/@jaskaranvirdi/setting-up-opencv-and-c-development-environment-in-xcode-b6027728003
using namespace cv; // 使用命名空间cv
如此可以减少输入,例如 cv :: Mat 就可省略为 Mat
Mat
Mat的优点是不再需要手动分配其内存,并在不需要它时立即发布它。在执行此操作仍然是可能的情况下,大多数OpenCV功能将自动分配其输出数据。
Mat作为一个类,包含
- 矩阵头(包含矩阵的大小,用于存储的方法,存储在哪个地址的信息等等)
- 指向包含像素值(取决于所选存储方法的任何维度)
从文件中加载图像:
Mat img = imread(filename);
如果需要加载灰度图:
Mat img = imread(filename, IMREAD_GRAYSCALE);
显示图像:
namedWindow("图片"); //打开名为“图片”的窗口
imshow("图片", img); //显示图像
openCV加载图像显示
using namespace cv; // 使用命名空间cv
int main()
{
Mat img = imread("building.jpg", IMREAD_GRAYSCALE); //从文件中加载灰度图像
//显示图像
namedWindow("原图");
imshow("原图", img);
waitKey(); //等待键值输入
return 0;
}
如何访问图像每一个像素点
利用指针访问,调用 Mat::ptr(i) 来获取第i行的首地址,通过循环进行访问。
// 按行遍历所有点(单通道)
for (int j = 0; j < nr; j++) {
for (int i = 0; i < nc; i++) {
//每个点为 img.ptr<uchar>(j)[i]
}
}
二、用高斯滤波器平滑图像
高斯滤波器(openCV)
openCV自带的高斯滤波器:cv :: GaussianBlur
void cv::GaussianBlur | ( | InputArray | src, |
OutputArray | dst, | ||
Size | ksize, | ||
double | sigmaX, | ||
double | sigmaY = 0 , | ||
int | borderType = BORDER_DEFAULT | ||
) |
#include <opencv2/imgproc.hpp>
使用高斯滤镜模糊图像。
该函数将源图像与指定的高斯内核进行卷积。支持就地过滤。
参量
src | 输入图像;图像可以具有任意数量的通道,这些通道可以独立处理,但深度应为CV_8U,CV_16U,CV_16S,CV_32F或CV_64F。 |
dst | 输出与src大小和类型相同的图像。 |
size | 高斯核大小。ksize.width和ksize.height可以不同,但它们都必须为正数和奇数。或者,它们可以为零,然后根据sigma计算得出。 |
sigmaX | X方向上的高斯核标准偏差。 |
sigmaY | Y方向的高斯核标准差;如果sigmaY为零,则将其设置为等于sigmaX;如果两个sigmas为零,则分别从ksize.width和ksize.height计算得出(有关详细信息,请参见getGaussianKernel);为了完全控制结果,而不管将来可能对所有这些语义进行的修改,建议指定所有ksize,sigmaX和sigmaY。 |
borderType | 像素外推方法,请参见BorderTypes |
高斯滤波器的C++实现
- 对图像使用一维高斯卷积模版,在一个方向上进行滤波(例如水平方向);
- 转置图像;
- 对转置后的图像使用同一个高斯卷积模版,在同样的方向上进行滤波;
- 将图像转置回原来的位置,得到二维高斯滤波的图像。
一维高斯卷积模版可以由二项式展开的系数来模拟,如3*3模版: 1/4 * [1 2 1]
/**
高斯滤波器,利用3*3的高斯模版进行高斯卷积
img 输入原图像
dst 高斯滤波后的输出图像
*/
void gaussianFilter(Mat &img, Mat &dst) {
// 对水平方向进行滤波
Mat dst1 = img.clone();
gaussianConvolution(img, dst1);
//图像矩阵转置
Mat dst2;
transpose(dst1, dst2);
// 对垂直方向进行滤波
Mat dst3 = dst2.clone();
gaussianConvolution(dst2, dst3);
// 再次转置
transpose(dst3, dst);
}
/**
一维高斯卷积,对每行进行高斯卷积
img 输入原图像
dst 一维高斯卷积后的输出图像
*/
void gaussianConvolution(Mat &img, Mat &dst) {
int nr = img.rows;
int nc = img.cols;
int templates[3] = {1, 2, 1};
// 按行遍历除每行边缘点的所有点
for (int j = 0; j < nr; j++) {
uchar* data= img.ptr<uchar>(j); //提取该行地址
for (int i = 1; i < nc-1; i++) {
int sum = 0;
for (int n = 0; n < 3; n++) {
sum += data[i-1+n] * templates[n]; //相称累加
}
sum /= 4;
dst.ptr<uchar>(j)[i] = sum;
}
}
}
高斯滤波前后的图像:
![](https://img-blog.csdnimg.cn/20191108122357494.jpg?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0hVU1RFUl9HeQ==,size_16,color_FFFFFF,t_70)
三、用一阶偏导有限差分计算梯度幅值和方向
用一阶偏导有限差分计算偏导数的两个阵列P与Q
再由P和Q算出梯度幅值和方向角
/**
用一阶偏导有限差分计算梯度幅值和方向
img 输入原图像
gradXY 输出的梯度幅值
theta 输出的梯度方向
*/
void getGrandient (Mat &img, Mat &gradXY, Mat &theta) {
gradXY = Mat::zeros(img.size(), CV_8U);
theta = Mat::zeros(img.size(), CV_8U);
for (int i = 0; i < img.rows-1; i++) {
for (int j = 0; j < img.cols-1; j++) {
double p = (img.ptr<uchar>(j)[i+1] - img.ptr<uchar>(j)[i] + img.ptr<uchar>(j+1)[i+1] - img.ptr<uchar>(j+1)[i])/2;
double q = (img.ptr<uchar>(j+1)[i] - img.ptr<uchar>(j)[i] + img.ptr<uchar>(j+1)[i+1] - img.ptr<uchar>(j)[i+1])/2;
gradXY.ptr<uchar>(j)[i] = sqrt(p*p + q*q); //计算梯度
theta.ptr<uchar>(j)[i] = atan(q/p);
}
}
}
此时输入输出图像为:
![](https://img-blog.csdnimg.cn/20191108121702859.jpg?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0hVU1RFUl9HeQ==,size_16,color_FFFFFF,t_70)
可以看出,二维计算梯度只区分出了部分边界,边界损失过大,于是采用三维算法计算梯度((y,x)为a11)。
a00 | a01 | a02 |
a10 | a11 | a12 |
a20 | a21 | a22 |
double gradX = double(a02 + 2 * a12 + a22 - a00 - 2 * a10 - a20);
double gradY = double(a00 + 2 * a01 + a02 - a20 - 2 * a21 - a22);
/**
用一阶偏导有限差分计算梯度幅值和方向
img 输入原图像
gradXY 输出的梯度幅值
theta 输出的梯度方向
*/
void getGrandient (Mat &img, Mat &gradXY, Mat &theta) {
gradXY = Mat::zeros(img.size(), CV_8U);
theta = Mat::zeros(img.size(), CV_8U);
for (int j = 1; j < img.rows-1; j++) {
for (int i = 1; i < img.cols-1; i++) {
double gradY = double(img.ptr<uchar>(j-1)[i-1] + 2 * img.ptr<uchar>(j-1)[i] + img.ptr<uchar>(j-1)[i+1] - img.ptr<uchar>(j+1)[i-1] - 2 * img.ptr<uchar>(j+1)[i] - img.ptr<uchar>(j+1)[i+1]);
double gradX = double(img.ptr<uchar>(j-1)[i+1] + 2 * img.ptr<uchar>(j)[i+1] + img.ptr<uchar>(j+1)[i+1] - img.ptr<uchar>(j-1)[i-1] - 2 * img.ptr<uchar>(j)[i-1] - img.ptr<uchar>(j+1)[i-1]);
gradXY.ptr<uchar>(j)[i] = sqrt(gradX*gradX + gradY*gradY); //计算梯度
theta.ptr<uchar>(j)[i] = atan(gradY/gradX); //计算梯度方向
}
}
}
三维梯度算法的输入输出图像:
![](https://img-blog.csdnimg.cn/20191108132415101.jpg?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0hVU1RFUl9HeQ==,size_16,color_FFFFFF,t_70)
四、对梯度幅值应用非极大值抑制
仅仅得到全局梯度并不足以确定边缘,保留局部梯度最大的点,而抑制非极大点。
- 将梯度角的变化范围减小到圆周的四个扇区之一;
- 四个扇区的标号为0到3,对应3*3领域的四种可能组合方向;
- 每一个点上领域的中心像素M与沿着梯度线的两个像素比较;
- 如果M梯度值不比沿梯度线的两个相邻像素梯度值大,则令M=0。
由 atan() 得到的角度在 范围内,将此范围均分为四个等份。
/**
局部非极大值抑制
gradXY 输入的梯度幅值
theta 输入的梯度方向
dst 输出的经局部非极大值抑制后的图像
*/
void nonLocalMaxValue (Mat &gradXY, Mat &theta, Mat &dst) {
dst = gradXY.clone();
for (int j = 1; j < gradXY.rows-1; j++) {
for (int i = 1; i < gradXY.cols-1; i++) {
double t = double(theta.ptr<uchar>(j)[i]);
double g = double(dst.ptr<uchar>(j)[i]);
if (g == 0.0) {
continue;
}
double g0, g1;
if ((t >= -(3*M_PI/8)) && (t < -(M_PI/8))) {
g0 = double(dst.ptr<uchar>(j-1)[i-1]);
g1 = double(dst.ptr<uchar>(j+1)[i+1]);
}
else if ((t >= -(M_PI/8)) && (t < M_PI/8)) {
g0 = double(dst.ptr<uchar>(j)[i-1]);
g1 = double(dst.ptr<uchar>(j)[i+1]);
}
else if ((t >= M_PI/8) && (t < 3*M_PI/8)) {
g0 = double(dst.ptr<uchar>(j-1)[i+1]);
g1 = double(dst.ptr<uchar>(j+1)[i-1]);
}
else {
g0 = double(dst.ptr<uchar>(j-1)[i]);
g1 = double(dst.ptr<uchar>(j+1)[i]);
}
if (g <= g0 || g <= g1) {
dst.ptr<uchar>(j)[i] = 0.0;
}
}
}
}
输入的经梯度计算后的图像和输出的局部非极大值抑制后的图像:
![](https://img-blog.csdnimg.cn/20191108174023525.jpg?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0hVU1RFUl9HeQ==,size_16,color_FFFFFF,t_70)
五、用双阈值算法检测和连接边缘
1、Canny算法采用双阈值,高阈值一般是低阈值的两倍,遍历所有像素点:
X < 低阈值 ,像素点置0,被抑制掉;
低阈值 < X <高阈值,像素点为弱边缘点,像素点值先不变;
X > 高阈值,像素点为强边缘点,置255。
2、弱边缘点补充连接强边缘点:
如果弱边缘点的8邻点域存在强边缘点,则将此点置255,用以连接强边缘点;如果不存在强边缘点,则这是一个孤立的弱边缘点,此点置0。
/**
用双阈值算法检测和连接边缘
low 输入的低阈值
high 输入的高阈值
img 输入的原图像
dst 输出的用双阈值算法检测和连接边缘后的图像
*/
void doubleThreshold (double low, double high, Mat &img, Mat &dst) {
dst = img.clone();
// 区分出弱边缘点和强边缘点
for (int j = 0; j < img.rows-1; j++) {
for (int i = 0; i < img.cols-1; i++) {
double x = double(dst.ptr<uchar>(j)[i]);
// 像素点为强边缘点,置255
if (x > high) {
dst.ptr<uchar>(j)[i] = 255;
}
// 像素点置0,被抑制掉
else if (x < low) {
dst.ptr<uchar>(j)[i] = 0;
}
}
}
// 弱边缘点补充连接强边缘点
doubleThresholdLink(dst);
}
/**
弱边缘点补充连接强边缘点
img 弱边缘点补充连接强边缘点的输入和输出图像
*/
void doubleThresholdLink (Mat &img) {
// 循环找到强边缘点,把其领域内的弱边缘点变为强边缘点
for (int j = 1; j < img.rows-2; j++) {
for (int i = 1; i < img.cols-2; i++) {
// 如果该点是强边缘点
if (img.ptr<uchar>(j)[i] == 255) {
// 遍历该强边缘点领域
for (int m = -1; m < 1; m++) {
for (int n = -1; n < 1; n++) {
// 该点为弱边缘点(不是强边缘点,也不是被抑制的0点)
if (img.ptr<uchar>(j+m)[i+n] != 0 && img.ptr<uchar>(j+m)[i+n] != 255) {
img.ptr<uchar>(j+m)[i+n] = 255; //该弱边缘点补充为强边缘点
}
}
}
}
}
}
for (int j = 0; j < img.rows-1; j++) {
for (int i = 0; i < img.cols-1; i++) {
// 如果该点依旧是弱边缘点,及此点是孤立边缘点
if (img.ptr<uchar>(j)[i] != 255 && img.ptr<uchar>(j)[i] != 255) {
img.ptr<uchar>(j)[i] = 0; //该孤立弱边缘点抑制
}
}
}
}
双阈值算法前后的输入输出图像 :
![](https://img-blog.csdnimg.cn/20191109012221209.jpg?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0hVU1RFUl9HeQ==,size_16,color_FFFFFF,t_70)
Canny边缘检测代码
#include <opencv2/opencv.hpp>
#include <math.h>
#define _USE_MATH_DEFINES
using namespace cv;// 使用命名空间cv
/**
将两个图像拼接,以便在同一个窗口显示
dst 输出的拼接后的图像
src1 拼接的第一幅图
src2 拼接的第二幅图
*/
void mergeImg(Mat & dst,Mat &src1,Mat &src2) {
int rows = src1.rows;
int cols = src1.cols+5+src2.cols;
CV_Assert(src1.type () == src2.type ());
dst.create (rows,cols,src1.type ());
src1.copyTo (dst(Rect(0,0,src1.cols,src1.rows)));
src2.copyTo (dst(Rect(src1.cols+5,0,src2.cols,src2.rows)));
}
/**
一维高斯卷积,对每行进行高斯卷积
img 输入原图像
dst 一维高斯卷积后的输出图像
*/
void gaussianConvolution(Mat &img, Mat &dst) {
int nr = img.rows;
int nc = img.cols;
int templates[3] = {1, 2, 1};
// 按行遍历除每行边缘点的所有点
for (int j = 0; j < nr; j++) {
uchar* data= img.ptr<uchar>(j); //提取该行地址
for (int i = 1; i < nc-1; i++) {
int sum = 0;
for (int n = 0; n < 3; n++) {
sum += data[i-1+n] * templates[n]; //相称累加
}
sum /= 4;
dst.ptr<uchar>(j)[i] = sum;
}
}
}
/**
高斯滤波器,利用3*3的高斯模版进行高斯卷积
img 输入原图像
dst 高斯滤波后的输出图像
*/
void gaussianFilter(Mat &img, Mat &dst) {
// 对水平方向进行滤波
Mat dst1 = img.clone();
gaussianConvolution(img, dst1);
//图像矩阵转置
Mat dst2;
transpose(dst1, dst2);
// 对垂直方向进行滤波
Mat dst3 = dst2.clone();
gaussianConvolution(dst2, dst3);
// 再次转置
transpose(dst3, dst);
}
/**
用一阶偏导有限差分计算梯度幅值和方向
img 输入原图像
gradXY 输出的梯度幅值
theta 输出的梯度方向
*/
void getGrandient (Mat &img, Mat &gradXY, Mat &theta) {
gradXY = Mat::zeros(img.size(), CV_8U);
theta = Mat::zeros(img.size(), CV_8U);
for (int j = 1; j < img.rows-1; j++) {
for (int i = 1; i < img.cols-1; i++) {
double gradY = double(img.ptr<uchar>(j-1)[i-1] + 2 * img.ptr<uchar>(j-1)[i] + img.ptr<uchar>(j-1)[i+1] - img.ptr<uchar>(j+1)[i-1] - 2 * img.ptr<uchar>(j+1)[i] - img.ptr<uchar>(j+1)[i+1]);
double gradX = double(img.ptr<uchar>(j-1)[i+1] + 2 * img.ptr<uchar>(j)[i+1] + img.ptr<uchar>(j+1)[i+1] - img.ptr<uchar>(j-1)[i-1] - 2 * img.ptr<uchar>(j)[i-1] - img.ptr<uchar>(j+1)[i-1]);
gradXY.ptr<uchar>(j)[i] = sqrt(gradX*gradX + gradY*gradY); //计算梯度
theta.ptr<uchar>(j)[i] = atan(gradY/gradX); //计算梯度方向
}
}
}
/**
局部非极大值抑制
gradXY 输入的梯度幅值
theta 输入的梯度方向
dst 输出的经局部非极大值抑制后的图像
*/
void nonLocalMaxValue (Mat &gradXY, Mat &theta, Mat &dst) {
dst = gradXY.clone();
for (int j = 1; j < gradXY.rows-1; j++) {
for (int i = 1; i < gradXY.cols-1; i++) {
double t = double(theta.ptr<uchar>(j)[i]);
double g = double(dst.ptr<uchar>(j)[i]);
if (g == 0.0) {
continue;
}
double g0, g1;
if ((t >= -(3*M_PI/8)) && (t < -(M_PI/8))) {
g0 = double(dst.ptr<uchar>(j-1)[i-1]);
g1 = double(dst.ptr<uchar>(j+1)[i+1]);
}
else if ((t >= -(M_PI/8)) && (t < M_PI/8)) {
g0 = double(dst.ptr<uchar>(j)[i-1]);
g1 = double(dst.ptr<uchar>(j)[i+1]);
}
else if ((t >= M_PI/8) && (t < 3*M_PI/8)) {
g0 = double(dst.ptr<uchar>(j-1)[i+1]);
g1 = double(dst.ptr<uchar>(j+1)[i-1]);
}
else {
g0 = double(dst.ptr<uchar>(j-1)[i]);
g1 = double(dst.ptr<uchar>(j+1)[i]);
}
if (g <= g0 || g <= g1) {
dst.ptr<uchar>(j)[i] = 0.0;
}
}
}
}
/**
弱边缘点补充连接强边缘点
img 弱边缘点补充连接强边缘点的输入和输出图像
*/
void doubleThresholdLink (Mat &img) {
// 循环找到强边缘点,把其领域内的弱边缘点变为强边缘点
for (int j = 1; j < img.rows-2; j++) {
for (int i = 1; i < img.cols-2; i++) {
// 如果该点是强边缘点
if (img.ptr<uchar>(j)[i] == 255) {
// 遍历该强边缘点领域
for (int m = -1; m < 1; m++) {
for (int n = -1; n < 1; n++) {
// 该点为弱边缘点(不是强边缘点,也不是被抑制的0点)
if (img.ptr<uchar>(j+m)[i+n] != 0 && img.ptr<uchar>(j+m)[i+n] != 255) {
img.ptr<uchar>(j+m)[i+n] = 255; //该弱边缘点补充为强边缘点
}
}
}
}
}
}
for (int j = 0; j < img.rows-1; j++) {
for (int i = 0; i < img.cols-1; i++) {
// 如果该点依旧是弱边缘点,及此点是孤立边缘点
if (img.ptr<uchar>(j)[i] != 255 && img.ptr<uchar>(j)[i] != 255) {
img.ptr<uchar>(j)[i] = 0; //该孤立弱边缘点抑制
}
}
}
}
/**
用双阈值算法检测和连接边缘
low 输入的低阈值
high 输入的高阈值
img 输入的原图像
dst 输出的用双阈值算法检测和连接边缘后的图像
*/
void doubleThreshold (double low, double high, Mat &img, Mat &dst) {
dst = img.clone();
// 区分出弱边缘点和强边缘点
for (int j = 0; j < img.rows-1; j++) {
for (int i = 0; i < img.cols-1; i++) {
double x = double(dst.ptr<uchar>(j)[i]);
// 像素点为强边缘点,置255
if (x > high) {
dst.ptr<uchar>(j)[i] = 255;
}
// 像素点置0,被抑制掉
else if (x < low) {
dst.ptr<uchar>(j)[i] = 0;
}
}
}
// 弱边缘点补充连接强边缘点
doubleThresholdLink(dst);
}
int main () {
Mat img = imread("woman.jpg", IMREAD_GRAYSCALE); //从文件中加载灰度图像
// 读取图片失败,则停止
if (img.empty()) {
printf("读取图像文件失败");
system("pause");
return 0;
}
// 高斯滤波
Mat gauss_img;
gaussianFilter(img, gauss_img); //高斯滤波器
// 用一阶偏导有限差分计算梯度幅值和方向
Mat gradXY, theta;
getGrandient(gauss_img, gradXY, theta);
// 局部非极大值抑制
Mat local_img;
nonLocalMaxValue(gradXY, theta, local_img);
// 用双阈值算法检测和连接边缘
Mat dst;
doubleThreshold(40, 80, local_img, dst);
// 图像显示
Mat outImg;
mergeImg (outImg,img,dst); //图像拼接
namedWindow("img");
imshow("img",outImg);// 图像显示
imwrite("canny算法.jpg", outImg);
waitKey(); //等待键值输入
return 0;
}
Canny边缘检测的前后图像
![](https://img-blog.csdnimg.cn/20191109012754920.jpg?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0hVU1RFUl9HeQ==,size_16,color_FFFFFF,t_70)
参考资源:
【1】Setting up OpenCV and C++ development environment in Xcode for Computer Vision projects,地址:https://medium.com/@jaskaranvirdi/setting-up-opencv-and-c-development-environment-in-xcode-b6027728003
【2】OpenCV Tutorials,地址:https://docs.opencv.org/master/d9/df8/tutorial_root.html
【3】OpenCV教程,地址:https://www.w3cschool.cn/opencv/opencv-2gnx28u3.html