canny实现边缘检测
主要流程见主函数
主要参考：https://blog.csdn.net/dcrmg/article/details/52344902
https://www.cnblogs.com/love6tao/p/5152020.html
不过有些细节有问题，我改过来了
#include <iostream>
#include<opencv2/opencv.hpp>
#include <string>using namespace std;
using namespace cv;#define PI 3.141592653589793void conv(char *G, Mat pic, Mat res)
{
int i, j, i1, j1;int tmp;for (i = 1; i < pic.rows-1; i++){
for (j = 1; j < pic.cols - 1; j++){
tmp = 0;for (i1 = -1; i1 <= 1; i1++){
for (j1 = -1; j1 <= 1; j1++){
tmp += pic.at<uchar>(i + i1, j + j1)*G[(j1 + 1) * 3 + i1 + 1];}}res.at<short>(i, j) = tmp;}}
}void merge1(Mat resX, Mat resY, Mat &res, Mat &theta)
{
int tmp;double theTmp;for (int i = 0; i < res.rows; i++){
for (int j = 0; j < res.cols; j++){
tmp = abs(resX.at<short>(i, j)) + abs(resY.at<short>(i, j));//tmp = sqrt(resX.at<short>(i, j)*resX.at<short>(i, j) + resY.at<short>(i, j)*resY.at<short>(i, j));if (tmp > 255)tmp = 255;res.at<uchar>(i,j) = tmp;theTmp = atan2(resY.at<short>(i, j), resX.at<short>(i, j))*180/PI;if (theTmp < 0)theTmp += 360;theta.at<float>(i, j) = theTmp;}}
}void mySobel(Mat pic, Mat &src, Mat &theta)
{
char Gx[3 * 3] = {
-1, 0, 1,-2, 0, 2,-1, 0, 1};char Gy[3 * 3] = {
1, 2, 1,0, 0, 0,-1, -2, -1};Mat resX = Mat::zeros(pic.rows, pic.cols, CV_16SC1);Mat resY = Mat::zeros(pic.rows, pic.cols, CV_16SC1);conv(Gx, pic, resX);conv(Gy, pic, resY);merge1(resX, resY, src, theta);}void BGR2Gray(Mat src, Mat &grayPic)
{
for (int i = 0; i < src.rows; i++){
for (int j = 0; j < src.cols; j++){
//Gray = R*0.299 + G*0.587 + B*0.114grayPic.at<uchar>(i, j) = src.at<Vec3b>(i, j)[0] * 0.114 + src.at<Vec3b>(i, j)[1] * 0.587 + src.at<Vec3b>(i, j)[2] * 0.299;}}
}void getGaussKer(double **K, int size, double sigma)
{
int center = size / 2;double sum = 0;for (int i = 0; i < size; i++){
for (int j = 0; j < size; j++){
K[i][j] = exp(((i - center)*(i - center) + (j - center)*(j - center)) / (-2 * sigma*sigma));///最后还是需要归一化加不加系数无所谓sum += (K[i][j]);}}for (int i = 0; i < size; i++){
for (int j = 0; j < size; j++){
K[i][j] /= sum;}}}void GaussFilter(Mat pic, Mat &src, double **K, int size)
{
int center = size / 2;for (int i = 0; i < pic.rows; i++){
for (int j = 0; j < pic.cols; j++){
double tmp = 0;int nowRow;int nowCol;for (int i1 = -size / 2; i1 <= size / 2; i1++){
for (int j1 = -size / 2; j1 <= size / 2; j1++){
nowRow = i + i1;if (nowRow < 0|| nowRow>=pic.rows)continue;nowCol = j + j1;if (nowCol < 0 || nowCol >= pic.cols)continue;tmp += K[center + i1][center + j1] * pic.at<uchar>(nowRow, nowCol);}}if (tmp < 0)tmp = 0;if (tmp > 255)tmp = 255;src.at<uchar>(i, j) = tmp;}}
}void locMaxValue(Mat amp, Mat theta, Mat &src)
{
Mat output = amp.clone();深拷贝for (int i = 1; i < amp.rows - 1; i++){
for (int j = 1; j < amp.cols - 1; j++)///判断每个像素点，最外一圈8领域不完整就不管了{
float thetaNow = theta.at<float>(i, j);float k=abs(tan(thetaNow));float g1, g2, g3, g4;float dTmp1, dTmp2;float nowGray = amp.at<uchar>(i, j);if ((thetaNow >= 90 && thetaNow < 135) || (thetaNow >= 270 && thetaNow < 315)){
g1 = output.at<uchar>(i - 1, j-1);g2 = output.at<uchar>(i - 1, j);g3 = output.at<uchar>(i + 1, j);g4 = output.at<uchar>(i + 1, j + 1);dTmp1 = g1 / k + (1 - 1 / k)*g2;dTmp2 = g4 / k + (1 - 1 / k)*g3;}else if ((thetaNow >= 135 && thetaNow < 180) || (thetaNow >= 315 && thetaNow < 360)){
g1 = output.at<uchar>(i - 1, j - 1);g2 = output.at<uchar>(i, j - 1);g3 = output.at<uchar>(i, j + 1);g4 = output.at<uchar>(i + 1, j + 1);dTmp1 = g1*k + (1 - k)*g2;dTmp2 = g4*k + (1 - k)*g3;}else if ((thetaNow >= 45 && thetaNow < 90) || (thetaNow >= 225 && thetaNow < 270)){
g1 = output.at<uchar>(i - 1, j );g2 = output.at<uchar>(i-1, j + 1);g3 = output.at<uchar>(i-1, j );g4 = output.at<uchar>(i + 1, j + 1);dTmp1 = g2 / k + (1 - 1 / k)*g1;dTmp2 = g3 / k + (1 - 1 / k)*g4;}else if (thetaNow >= 0 && thetaNow < 45 || (thetaNow >= 180 && thetaNow < 225)){
g1 = output.at<uchar>(i - 1, j+1);g2 = output.at<uchar>(i , j + 1);g3 = output.at<uchar>(i , j-1);g4 = output.at<uchar>(i + 1, j - 1);dTmp1 = g1*k + (1 - k)*g2;dTmp2 = g4*k + (1 - k)*g3;}if (nowGray < dTmp1 || nowGray < dTmp2){
output.at<uchar>(i, j) = 0;}}}src = output.clone();}void traceLink(Mat &src, int minValue, int i, int j)
{
int xn[8] = {
1,1,0,-1,-1,-1,0,1 };int yn[8] = {
0, 1, 1, 1, 0, -1, -1, -1 };int xx = 0;int yy = 0;bool change = true;while (change){
change = false;for (int k = 0; k < 8; k++){
xx = i + xn[k];yy = j + yn[k];if (src.at<uchar>(xx,yy)!=255&&src.at<uchar>(xx, yy) > minValue){
change = true;src.at<uchar>(xx, yy) = 255;i = xx;j = yy;break;}}}}void doubleThreshold(Mat &src, int maxValue, int minValue)
{
for (int i = 1; i < src.rows-1; i++){
for (int j = 1; j < src.cols-1; j++){
int value = src.at<uchar>(i, j);//int t=value;if (src.at<uchar>(i, j) >= maxValue){
src.at<uchar>(i, j) = 255;traceLink(src, minValue, i, j);}}}for (int i = 1; i < src.rows-1; i++){
for (int j = 1; j < src.cols-1; j++){
if (src.at<uchar>(i, j) != 255){
src.at<uchar>(i, j) = 0;}}}
}int main()
{
int size = 3;double sigma = 1;double **k = new double*[size];for (int i = 0; i < size; i++){
k[i] = new double[size];}Mat pic = imread("C:\\Users\\Thinkpad\\Desktop\\c++work\\pic\\lena.jpg", 1);Mat grayPic=Mat::zeros(pic.rows, pic.cols, CV_8UC1); Mat gaussPic=Mat::zeros(pic.rows, pic.cols, CV_8UC1);Mat amplitude = Mat::zeros(pic.rows, pic.cols, CV_8UC1);Mat theta = Mat::zeros(pic.rows, pic.cols, CV_32FC1);Mat src;BGR2Gray(pic, grayPic);转成灰度图getGaussKer(k, size, sigma);///计算高斯核GaussFilter(grayPic, gaussPic, k, 3);///计算高斯模糊后的图mySobel(gaussPic, amplitude, theta);///用sobel算子计算幅值和角度locMaxValue(amplitude, theta, src);非极大抑制，梯度方向不是极大值的进行抑制doubleThreshold(src, 100, 60);双阈值链接imshow("1", src);imshow("0", amplitude);waitKey();/*for (int i = 0; i < 3; i++){for (int j = 0; j < 3; j++){cout << k[i][j] << ",";}cout << endl;}*/for (int i = 0; i < size; i++){
delete[]k[i];}delete[]k;return 0;
}
这里简单介绍下梯度算子的方向：通常OpenCV打开的的图像坐标原点在左上角，x轴方向向下，y轴方向向右，但是在选择算子进行梯度计算时比如sobel算子，其“X”方向的计算是右边减去左边，“Y”方向的计算是上边减去下边，这种梯度的计算方式和我们正常左手系坐标系的计算方式相同，也就是x轴方向向右，y轴方向向上，而dx,dy的计算与坐标系选取无关，但是由dx,dy计算得到的角度是与“X”轴正方向的夹角，所以无论dx,dy是基于哪个坐标系算的，sobel算子计算的梯度方向都是以右为x轴正方向的左手系坐标系的夹角，下图是基于此，我进行非极大抑制的计算过程草图
