1
+ #include " opencv2/opencv.hpp"
2
+ #include " opencv2/imgproc/imgproc.hpp"
3
+ #include " iostream"
4
+ #include " algorithm"
5
+ #include " vector"
6
+ #include " stdio.h"
7
+ using namespace std ;
8
+ using namespace cv ;
9
+
10
+ const double PI = 3.1415926 ;
11
+
12
+ double getSum (Mat src){
13
+ double sum = 0 ;
14
+ for (int i = 0 ; i < src.rows ; i++){
15
+ for (int j = 0 ; j < src.cols ; j++){
16
+ sum += (double )src.at <double >(i, j);
17
+ }
18
+ }
19
+ return sum;
20
+ }
21
+
22
+ Mat CannyEdgeDetection (cv::Mat src, int kSize , double hightThres, double lowThres){
23
+ // if(src.channels() == 3){
24
+ // cvtColor(src, src, CV_BGR2GRAY);
25
+ // }
26
+ cv::Rect rect;
27
+ Mat gaussResult;
28
+ int row = src.rows ;
29
+ int col = src.cols ;
30
+ printf (" %d %d\n " , row, col);
31
+ Mat filterImg = Mat::zeros (row, col, CV_64FC1);
32
+ src.convertTo (src, CV_64FC1);
33
+ Mat dst = Mat::zeros (row, col, CV_64FC1);
34
+ int gaussCenter = kSize / 2 ;
35
+ double sigma = 1 ;
36
+ Mat guassKernel = Mat::zeros (kSize , kSize , CV_64FC1);
37
+ for (int i = 0 ; i < kSize ; i++){
38
+ for (int j = 0 ; j < kSize ; j++){
39
+ guassKernel.at <double >(i, j) = (1.0 / (2.0 * PI * sigma * sigma)) * (double )exp (-(((double )pow ((i - (gaussCenter+1 )), 2 ) + (double )pow ((j-(gaussCenter+1 )), 2 )) / (2.0 *sigma*sigma)));
40
+ }
41
+ }
42
+ Scalar sumValueScalar = cv::sum (guassKernel);
43
+ double sum = sumValueScalar.val [0 ];
44
+ guassKernel = guassKernel / sum;
45
+ for (int i = gaussCenter; i < row - gaussCenter; i++){
46
+ for (int j = gaussCenter; j < col - gaussCenter; j++){
47
+ rect.x = j - gaussCenter;
48
+ rect.y = i - gaussCenter;
49
+ rect.width = kSize ;
50
+ rect.height = kSize ;
51
+ // printf("%d %d\n", i, j);
52
+ // printf("%d %d %.5f\n", i, j, cv::sum(guassKernel.mul(src(rect))).val[0]);
53
+ filterImg.at <double >(i, j) = cv::sum (guassKernel.mul (src (rect))).val [0 ];
54
+ }
55
+ }
56
+ Mat gradX = Mat::zeros (row, col, CV_64FC1); // 水平梯度
57
+ Mat gradY = Mat::zeros (row, col, CV_64FC1); // 垂直梯度
58
+ Mat grad = Mat::zeros (row, col, CV_64FC1); // 梯度幅值
59
+ Mat thead = Mat::zeros (row, col, CV_64FC1); // 梯度角度
60
+ Mat whereGrad = Mat::zeros (row, col, CV_64FC1);// 区域
61
+ // x方向的Sobel算子
62
+ Mat Sx = (cv::Mat_<double >(3 , 3 ) << -1 , 0 , 1 , -2 , 0 , 2 , -1 , 0 , 1 );
63
+ // y方向的Sobel算子
64
+ Mat Sy = (cv::Mat_<double >(3 , 3 ) << 1 , 2 , 1 , 0 , 0 , 0 , -1 , -2 , -1 );
65
+ // 计算梯度赋值和角度
66
+ for (int i=1 ; i < row-1 ; i++){
67
+ for (int j=1 ; j < col-1 ; j++){
68
+ rect.x = j-1 ;
69
+ rect.y = i-1 ;
70
+ rect.width = 3 ;
71
+ rect.height = 3 ;
72
+ Mat rectImg = Mat::zeros (3 , 3 , CV_64FC1);
73
+ filterImg (rect).copyTo (rectImg);
74
+ // 梯度和角度
75
+ gradX.at <double >(i, j) += cv::sum (rectImg.mul (Sx)).val [0 ];
76
+ gradY.at <double >(i, j) += cv::sum (rectImg.mul (Sy)).val [0 ];
77
+ grad.at <double >(i, j) = sqrt (pow (gradX.at <double >(i, j), 2 ) + pow (gradY.at <double >(i, j), 2 ));
78
+ thead.at <double >(i, j) = atan (gradY.at <double >(i, j)/gradX.at <double >(i, j));
79
+ if (0 <= thead.at <double >(i, j) <= (PI/4.0 )){
80
+ whereGrad.at <double >(i, j) = 0 ;
81
+ }else if (PI/4.0 < thead.at <double >(i, j) <= (PI/2.0 )){
82
+ whereGrad.at <double >(i, j) = 1 ;
83
+ }else if (-PI/2.0 <= thead.at <double >(i, j) <= (-PI/4.0 )){
84
+ whereGrad.at <double >(i, j) = 2 ;
85
+ }else if (-PI/4.0 < thead.at <double >(i, j) < 0 ){
86
+ whereGrad.at <double >(i, j) = 3 ;
87
+ }
88
+ }
89
+ }
90
+ // printf("success\n");
91
+ // 梯度归一化
92
+ double gradMax;
93
+ cv::minMaxLoc (grad, &gradMax);
94
+ if (gradMax != 0 ){
95
+ grad = grad / gradMax;
96
+ }
97
+ // 双阈值确定
98
+ cv::Mat caculateValue = cv::Mat::zeros (row, col, CV_64FC1); // grad变成一维
99
+ resize (grad, caculateValue, Size (1 , grad.rows * grad.cols ));
100
+ cv::sort (caculateValue, caculateValue, CV_SORT_EVERY_COLUMN + CV_SORT_ASCENDING);// 升序
101
+ long long highIndex= row * col * hightThres;
102
+ double highValue = caculateValue.at <double >(highIndex, 0 ); // 最大阈值
103
+ double lowValue = highValue * lowThres;
104
+ // NMS
105
+ for (int i = 1 ; i < row-1 ; i++ ){
106
+ for ( int j = 1 ; j<col-1 ; j++){
107
+ // 八个方位
108
+ double N = grad.at <double >(i-1 , j);
109
+ double NE = grad.at <double >(i-1 , j+1 );
110
+ double E = grad.at <double >(i, j+1 );
111
+ double SE = grad.at <double >(i+1 , j+1 );
112
+ double S = grad.at <double >(i+1 , j);
113
+ double SW = grad.at <double >(i-1 , j-1 );
114
+ double W = grad.at <double >(i, j-1 );
115
+ double NW = grad.at <double >(i -1 , j -1 ); // 区域判断,线性插值处理
116
+ double tanThead; // tan角度
117
+ double Gp1; // 两个方向的梯度强度
118
+ double Gp2; // 求角度,绝对值
119
+ tanThead = abs (tan (thead.at <double >(i,j)));
120
+ switch ((int )whereGrad.at <double >(i,j)) {
121
+ case 0 : Gp1 = (1 - tanThead) * E + tanThead * NE; Gp2 = (1 - tanThead) * W + tanThead * SW; break ;
122
+ case 1 : Gp1 = (1 - tanThead) * N + tanThead * NE; Gp2 = (1 - tanThead) * S + tanThead * SW; break ;
123
+ case 2 : Gp1 = (1 - tanThead) * N + tanThead * NW; Gp2 = (1 - tanThead) * S + tanThead * SE; break ;
124
+ case 3 : Gp1 = (1 - tanThead) * W + tanThead *NW; Gp2 = (1 - tanThead) * E + tanThead *SE; break ;
125
+ default : break ;
126
+ }
127
+ // NMS -非极大值抑制和双阈值检测
128
+ if (grad.at <double >(i, j) >= Gp1 && grad.at <double >(i, j) >= Gp2){
129
+ // 双阈值检测
130
+ if (grad.at <double >(i, j) >= highValue){
131
+ grad.at <double >(i, j) = highValue;
132
+ dst.at <double >(i, j) = 255 ;
133
+ } else if (grad.at <double >(i, j) < lowValue){
134
+ grad.at <double >(i, j) = 0 ;
135
+ } else {
136
+ grad.at <double >(i, j) = lowValue;
137
+ }
138
+ } else {
139
+ grad.at <double >(i, j) = 0 ;
140
+ }
141
+ }
142
+ }
143
+ // 抑制孤立低阈值点 3*3. 找到高阈值就255
144
+ for (int i = 1 ; i < row - 1 ; i++){
145
+ for (int j = 1 ; j < col - 1 ; j++){
146
+ if (grad.at <double >(i, j) == lowValue){
147
+ // 3*3 区域强度
148
+ rect.x = i-1 ;
149
+ rect.y = j-1 ;
150
+ rect.width = 3 ;
151
+ rect.height = 3 ;
152
+ for (int x = 0 ; x < 3 ; x++){
153
+ for (int y = 0 ; y < 3 ; y++){
154
+ if (grad (rect).at <double >(x, y)==highValue){
155
+ dst.at <double >(i, j) = 255 ;
156
+ break ;
157
+ }
158
+ }
159
+ }
160
+ }
161
+ }
162
+ }
163
+ return dst;
164
+ }
165
+
166
+ int main (){
167
+ Mat src = imread (" ../lena.jpg" );
168
+ Mat gray;
169
+ cvtColor (src, gray, CV_BGR2GRAY);
170
+ Mat dst = CannyEdgeDetection (gray, 3 , 0.8 , 0.5 );
171
+ imshow (" origin" , src);
172
+ imshow (" gray" , gray);
173
+ imshow (" result" , dst);
174
+ waitKey (0 );
175
+ imwrite (" ../result.jpg" , gray);
176
+ imwrite (" ../result2.jpg" , dst);
177
+ return 0 ;
178
+ }
0 commit comments