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
+ #include " map"
8
+ #include " unordered_map"
9
+ using namespace std ;
10
+ using namespace cv ;
11
+
12
+ Mat RGB2GRAY (const Mat &src){
13
+ if (!src.data || src.channels ()!=3 ){
14
+ exit (1 );
15
+ }
16
+ int row = src.rows ;
17
+ int col = src.cols ;
18
+ Mat gray = Mat (row, col, CV_8UC1);
19
+ for (int i = 0 ; i < row; i++){
20
+ for (int j = 0 ; j < col; j++){
21
+ gray.at <uchar>(i, j) = (uchar)(0.114 * src.at <Vec3b>(i, j)[0 ] + 0.587 * src.at <Vec3b>(i, j)[1 ] + 0.299 * src.at <Vec3b>(i, j)[2 ]);
22
+ }
23
+ }
24
+ return gray;
25
+ }
26
+
27
+ void SobelGradDirction (Mat &src, Mat &sobelX, Mat &sobelY){
28
+ int row = src.rows ;
29
+ int col = src.cols ;
30
+ sobelX = Mat::zeros (src.size (), CV_32SC1);
31
+ sobelY = Mat::zeros (src.size (), CV_32SC1);
32
+
33
+ for (int i = 0 ; i < row-1 ; i++){
34
+ for (int j = 0 ; j < col-1 ; j++){
35
+ double gradY = src.at <uchar>(i+1 , j-1 ) + src.at <uchar>(i+1 , j) * 2 + src.at <uchar>(i+1 , j+1 ) - src.at <uchar>(i-1 , j-1 ) - src.at <uchar>(i-1 , j) * 2 - src.at <uchar>(i-1 , j+1 );
36
+ sobelY.at <float >(i, j) = abs (gradY);
37
+ double gradX = src.at <uchar>(i-1 , j+1 ) + src.at <uchar>(i, j+1 ) * 2 + src.at <uchar>(i+1 , j+1 ) - src.at <uchar>(i-1 , j-1 ) - src.at <uchar>(i, j-1 ) * 2 - src.at <uchar>(i+1 , j-1 );
38
+ sobelX.at <float >(i, j) = abs (gradX);
39
+ }
40
+ }
41
+ // 将梯度数组转换为8位无符号整形
42
+ convertScaleAbs (sobelX, sobelX);
43
+ convertScaleAbs (sobelY, sobelY);
44
+ }
45
+
46
+ Mat SobelXX (const Mat src){
47
+ int row = src.rows ;
48
+ int col = src.cols ;
49
+ Mat_<float > dst (row, col, CV_32FC1);
50
+ for (int i = 0 ; i < row; i++){
51
+ for (int j = 0 ; j < col; j++){
52
+ dst.at <float >(i, j) = src.at <uchar>(i, j) * src.at <uchar>(i, j);
53
+ }
54
+ }
55
+ return dst;
56
+ }
57
+
58
+ Mat SobelYY (const Mat src){
59
+ int row = src.rows ;
60
+ int col = src.cols ;
61
+ Mat_<float > dst (row, col, CV_32FC1);
62
+ for (int i = 0 ; i < row; i++){
63
+ for (int j = 0 ; j < col; j++){
64
+ dst.at <float >(i, j) = src.at <uchar>(i, j) * src.at <uchar>(i, j);
65
+ }
66
+ }
67
+ return dst;
68
+ }
69
+
70
+ Mat SobelXY (const Mat src1, const Mat &src2){
71
+ int row = src1.rows ;
72
+ int col = src1.cols ;
73
+ Mat_<float > dst (row, col, CV_32FC1);
74
+ for (int i = 0 ; i < row; i++){
75
+ for (int j = 0 ; j < col; j++){
76
+ dst.at <float >(i, j) = src1.at <uchar>(i, j) * src2.at <uchar>(i, j);
77
+ }
78
+ }
79
+ return dst;
80
+ }
81
+
82
+ void separateGaussianFilter (Mat_<float > &src, Mat_<float > &dst, int ksize, double sigma){
83
+ CV_Assert (src.channels ()==1 || src.channels () == 3 ); // 只处理单通道或者三通道图像
84
+ // 生成一维的
85
+ double *matrix = new double [ksize];
86
+ double sum = 0 ;
87
+ int origin = ksize / 2 ;
88
+ for (int i = 0 ; i < ksize; i++){
89
+ double g = exp (-(i-origin) * (i-origin) / (2 * sigma * sigma));
90
+ sum += g;
91
+ matrix[i] = g;
92
+ }
93
+ for (int i = 0 ; i < ksize; i++) matrix[i] /= sum;
94
+ int border = ksize / 2 ;
95
+ copyMakeBorder (src, dst, border, border, border, border, BORDER_CONSTANT);
96
+ int channels = dst.channels ();
97
+ int rows = dst.rows - border;
98
+ int cols = dst.cols - border;
99
+ // 水平方向
100
+ for (int i = border; i < rows - border; i++){
101
+ for (int j = border; j < cols - border; j++){
102
+ float sum[3 ] = {0 };
103
+ for (int k = -border; k<=border; k++){
104
+ if (channels == 1 ){
105
+ sum[0 ] += matrix[border + k] * dst.at <float >(i, j+k);
106
+ }else if (channels == 3 ){
107
+ Vec3f rgb = dst.at <Vec3f>(i, j+k);
108
+ sum[0 ] += matrix[border+k] * rgb[0 ];
109
+ sum[1 ] += matrix[border+k] * rgb[1 ];
110
+ sum[2 ] += matrix[border+k] * rgb[2 ];
111
+ }
112
+ }
113
+ for (int k = 0 ; k < channels; k++){
114
+ if (sum[k] < 0 ) sum[k] = 0 ;
115
+ else if (sum[k] > 255 ) sum[k] = 255 ;
116
+ }
117
+ if (channels == 1 )
118
+ dst.at <float >(i, j) = static_cast <float >(sum[0 ]);
119
+ else if (channels == 3 ){
120
+ Vec3f rgb = {static_cast <float >(sum[0 ]), static_cast <float >(sum[1 ]), static_cast <float >(sum[2 ])};
121
+ dst.at <Vec3f>(i, j) = rgb;
122
+ }
123
+ }
124
+ }
125
+ // 竖直方向
126
+ for (int i = border; i < rows - border; i++){
127
+ for (int j = border; j < cols - border; j++){
128
+ float sum[3 ] = {0 };
129
+ for (int k = -border; k<=border; k++){
130
+ if (channels == 1 ){
131
+ sum[0 ] += matrix[border + k] * dst.at <float >(i+k, j);
132
+ }else if (channels == 3 ){
133
+ Vec3f rgb = dst.at <Vec3f>(i+k, j);
134
+ sum[0 ] += matrix[border+k] * rgb[0 ];
135
+ sum[1 ] += matrix[border+k] * rgb[1 ];
136
+ sum[2 ] += matrix[border+k] * rgb[2 ];
137
+ }
138
+ }
139
+ for (int k = 0 ; k < channels; k++){
140
+ if (sum[k] < 0 ) sum[k] = 0 ;
141
+ else if (sum[k] > 255 ) sum[k] = 255 ;
142
+ }
143
+ if (channels == 1 )
144
+ dst.at <float >(i, j) = static_cast <float >(sum[0 ]);
145
+ else if (channels == 3 ){
146
+ Vec3f rgb = {static_cast <float >(sum[0 ]), static_cast <float >(sum[1 ]), static_cast <float >(sum[2 ])};
147
+ dst.at <Vec3f>(i, j) = rgb;
148
+ }
149
+ }
150
+ }
151
+ delete [] matrix;
152
+ }
153
+
154
+ Mat harrisResponse (Mat_<float > GaussXX, Mat_<float > GaussYY, Mat_<float > GaussXY, float k){
155
+ int row = GaussXX.rows ;
156
+ int col = GaussXX.cols ;
157
+ Mat_<float > dst (row, col, CV_32FC1);
158
+ for (int i = 0 ; i < row; i++){
159
+ for (int j = 0 ; j < col; j++){
160
+ float a = GaussXX.at <float >(i, j);
161
+ float b = GaussYY.at <float >(i, j);
162
+ float c = GaussXY.at <float >(i, j);
163
+ dst.at <float >(i, j) = a * b - c * c - k * (a + b) * (a + b);
164
+ }
165
+ }
166
+ return dst;
167
+ }
168
+
169
+ int dir[8 ][2 ] = {0 , 1 , 0 , -1 , 1 , 0 , -1 , 0 , 1 , 1 , 1 , -1 , -1 , 1 , -1 , -1 };
170
+
171
+ void MaxLocalValue (Mat_<float >&resultData, Mat srcGray, Mat &resultImage, int kSize ){
172
+ int r = kSize / 2 ;
173
+ resultImage = srcGray.clone ();
174
+ int row = resultImage.rows ;
175
+ int col = resultImage.cols ;
176
+ for (int i = r; i < row - r; i++){
177
+ for (int j = r; j < col - r; j++){
178
+ bool flag = true ;
179
+ for (int k = 0 ; k < 8 ; k++){
180
+ int tx = i + dir[k][0 ];
181
+ int ty = j + dir[k][1 ];
182
+ if (resultData.at <float >(i, j) < resultData.at <float >(tx, ty)){
183
+ flag = false ;
184
+ }
185
+ }
186
+ if (flag && (int )resultData.at <float >(i, j) > 18000 ){
187
+ circle (resultImage, Point (i, j), 5 , Scalar (0 , 0 , 255 ), 2 , 8 , 0 );
188
+ }
189
+ }
190
+ }
191
+ }
192
+
193
+ int main (){
194
+ Mat src = imread (" ../lena.jpg" );
195
+ Mat srcGray = RGB2GRAY (src);
196
+ Mat imageSobelX;
197
+ Mat imageSobelY;
198
+ Mat resultImage;
199
+ Mat_<float > imageSobelXX, imageSobelYY, imageSobelXY;
200
+ Mat_<float > GaussianXX, GaussianXY, GaussianYY, HarrisResponse;
201
+ SobelGradDirction (srcGray, imageSobelX, imageSobelY);
202
+ imageSobelXX = SobelXX (imageSobelX);
203
+ imageSobelYY = SobelYY (imageSobelY);
204
+ imageSobelXY = SobelXY (imageSobelX, imageSobelY);
205
+ separateGaussianFilter (imageSobelXX, GaussianXX, 3 , 1.0 );
206
+ separateGaussianFilter (imageSobelYY, GaussianYY, 3 , 1.0 );
207
+ separateGaussianFilter (imageSobelXY, GaussianXY, 3 , 1.0 );
208
+ HarrisResponse = harrisResponse (GaussianXX, GaussianYY, GaussianXY, 0.05 );
209
+ MaxLocalValue (HarrisResponse, srcGray, resultImage, 3 );
210
+ imshow (" origin" , src);
211
+ imshow (" gray" , srcGray);
212
+ imshow (" result" , resultImage);
213
+ imwrite (" ../result.jpg" , resultImage);
214
+ waitKey (0 );
215
+ return 0 ;
216
+ }
0 commit comments