1 ORB test实验代码

来源：互联网发布：压缩感知凸优化编辑：程序博客网时间：2024/06/14 06:00

目录(?)[+]

提取左右图特征
做穷举式的最近邻检索
绘图
估计单应矩阵计算重投影误差
结果分析
完整代码如下

转：http://www.cvchina.info/2011/09/25/orb-test/

之前介绍了ORB，一种具备旋转不变形的局部特征描述子。OpenCV2.3中提供了实现，但是缺少使用例程。下面是一个简单的样例程序。

随便拍了两张图片作为测试图像。

下面上下两图分别为模板图像和查询图像：

提取左右图特征：

?

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

Mat
img1 = imread(image_filename1, 0);

Mat
img2 = imread(image_filename2, 0);

//GaussianBlur(img1,
img1, Size(5, 5), 0);

//GaussianBlur(img2,
img2, Size(5, 5), 0);

 

<a
href="http://www.cvchina.info/tag/orb/"

class="st_tag
internal_tag" 
rel="tag"

title="标签
orb 下的日志">ORB</a>
orb1(3000, <a href="http://www.cvchina.info/tag/orb/"

class="st_tag
internal_tag" 
rel="tag"

title="标签
orb 下的日志">ORB</a>::CommonParams(1.2,
8));

<a
href="http://www.cvchina.info/tag/orb/"

class="st_tag
internal_tag" 
rel="tag"

title="标签
orb 下的日志">ORB</a>
orb2(100, <a href="http://www.cvchina.info/tag/orb/"

class="st_tag
internal_tag" 
rel="tag"

title="标签
orb 下的日志">ORB</a>::CommonParams(1.2,
1));

 

vector
keys1, keys2;

Mat
descriptors1, descriptors2;

 

<a
href="http://www.cvchina.info/tag/orb/"

class="st_tag
internal_tag" 
rel="tag"

title="标签
orb 下的日志">orb</a>1(img1,
Mat(), keys1, descriptors1, false);

printf("tem
feat num: %d\n",
keys1.size());

 

int64
st, et;

st
= cvGetTickCount();

<a
href="http://www.cvchina.info/tag/orb/"

class="st_tag
internal_tag" 
rel="tag"

title="标签
orb 下的日志">orb</a>2(img2,
Mat(), keys2, descriptors2, false);

et
= cvGetTickCount();

printf("<a
href="http://www.cvchina.info/tag/orb/"
class="st_tag internal_tag" rel="tag" title="标签 orb 下的日志">orb</a>2 extraction time: %f\n", (et-st)/(double)cvGetTickFrequency()/1000.);

printf("query
feat num: %d\n",
keys2.size());

注：模板图像在多尺度提取特征，查询图像只在提取原始尺度上的特征。

做穷举式的最近邻检索：

?

1

2

3

4

5

6

7

8

9

10

11

//
find matches

vector
matches;

st
= cvGetTickCount();

//for(int
i = 0; i < 10; i++){

naive_nn_search2(keys1,
descriptors1, keys2, descriptors2, matches);

//}

et
= cvGetTickCount();

printf("match
time: %f\n",
(et-st)/(double)cvGetTickFrequency()/1000.);

printf("matchs
num: %d\n",
matches.size());

hamming距离测算通过查找表实现：

?

1

2

3

4

5

6

7

8

unsigned
int

hamdist2(unsigned char*
a, unsigned char*
b, size_t

size)

{

HammingLUT
lut;

unsigned
int

result;

result
= lut((a), (b), size);

return

result;

}

绘图：

?

1

2

3

4

5

6

Mat
showImg;

drawMatches(img2,
keys2, img1, keys1, matches, showImg, CV_RGB(0, 255, 0), CV_RGB(0, 0, 255));

string
winName = "Matches";

namedWindow(
winName, WINDOW_AUTOSIZE );

imshow(
winName, showImg );

waitKey();

估计单应矩阵，计算重投影误差：

?

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

Mat
homo;

st
= cvGetTickCount();

homo
= findHomography(pt1, pt2, Mat(), CV_RANSAC, 5);

et
= cvGetTickCount();

printf("ransac
time: %f\n",
(et-st)/(double)cvGetTickFrequency()/1000.);

printf("homo\n"

"%f
%f %f\n"

"%f
%f %f\n"

"%f
%f %f\n",

homo.at(0,0),
homo.at(0,1), homo.at(0,2),

homo.at(1,0),
homo.at(1,1), homo.at(1,2),

homo.at(2,0),homo.at(2,1),homo.at(2,2));

vector

reproj;

reproj.resize(pt1.size());

perspectiveTransform(pt1,
reproj, homo);

Mat
diff;

diff
= Mat(reproj) - Mat(pt2);

int

inlier = 0;

double

err_sum = 0;

for(int

i = 0; i < diff.rows; i++){

float*
ptr = diff.ptr(i);

float

err = ptr[0]*ptr[0] + ptr[1]*ptr[1];

if(err
< 25.f){

inlier++;

err_sum
+= sqrt(err);

}

}

printf("inlier
num: %d\n",
inlier);

printf("ratio
%f\n",
inlier / (float)(diff.rows));

printf("mean
reprojection error: %f\n",
err_sum / inlier);

结果分析：

tem feat num: 743
orb2 extraction time: 1.672435
query feat num: 100
match time: 3.698276
matchs num: 8
ransac time: 143.570586
homo
0.974942 0.410833 4.426035
-0.182418 0.828115 52.742661
0.001191 0.000144 1.000000
inlier num: 8
ratio 1.000000
mean reprojection error: 0.976777

可见最近邻检索是系统的瓶颈，（进行了743*100次hamming距离（32bytes）计算。）一个简单的优化如下，分段计算hamming距离，先计算前16byte的hamming距离，如超过某一阈值，则直接认为非候选，如小于某阈值，则继续进行后一半16bytes的距离计算。（粗略估计可以减少30%+的最近邻查询时间）。更复杂的办法是使用LSH，此处按下不提，有空再续。

完整代码如下：

?

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

#include
"opencv2/objdetect/objdetect.hpp"

#include
"opencv2/features2d/features2d.hpp"

#include
"opencv2/highgui/highgui.hpp"

#include
"opencv2/calib3d/calib3d.hpp"

#include
"opencv2/imgproc/imgproc_c.h"

#include
"opencv2/imgproc/imgproc.hpp"

#include

#include

#include

using

namespace 
std;

using

namespace 
cv;

char*
image_filename1 = "apple_vinegar_0.png";

char*
image_filename2 = "apple_vinegar_2.png";

unsigned
int

hamdist(unsigned int

x, unsigned int

y)

{

unsigned
int

dist = 0, val = x ^ y;

//
Count the number of set bits

while(val)

{

++dist;

val
&= val - 1;

}

return

dist;

}

unsigned
int

hamdist2(unsigned char*
a, unsigned char*
b, size_t

size)

{

HammingLUT
lut;

unsigned
int

result;

result
= lut((a), (b), size);

return

result;

}

void

naive_nn_search(vector& keys1, Mat& descp1,

vector&
keys2, Mat& descp2,

vector&
matches)

{

for(
int

i = 0; i < (int)keys2.size();
i++){

unsigned
int

min_dist = INT_MAX;

int

min_idx = -1;

unsigned
char*
query_feat = descp2.ptr(i);

for(
int

j = 0; j < (int)keys1.size();
j++){

unsigned
char*
train_feat = descp1.ptr(j);

unsigned
int

dist =  hamdist2(query_feat, train_feat, 32);

if(dist
< min_dist){

min_dist
= dist;

min_idx
= j;

}

}

//if(min_dist
<= (unsigned int)(second_dist * 0.8)){

if(min_dist
<= 50){

matches.push_back(DMatch(i,
min_idx, 0, (float)min_dist));

}

}

}

void

naive_nn_search2(vector& keys1, Mat& descp1,

vector&
keys2, Mat& descp2,

vector&
matches)

{

for(
int

i = 0; i < (int)keys2.size();
i++){

unsigned
int

min_dist = INT_MAX;

unsigned
int

sec_dist = INT_MAX;

int

min_idx = -1, sec_idx = -1;

unsigned
char*
query_feat = descp2.ptr(i);

for(
int

j = 0; j < (int)keys1.size();
j++){

unsigned
char*
train_feat = descp1.ptr(j);

unsigned
int

dist =  hamdist2(query_feat, train_feat, 32);

if(dist
< min_dist){

sec_dist
= min_dist;

sec_idx
= min_idx;

min_dist
= dist;

min_idx
= j;

}else

if(dist
< sec_dist){

sec_dist
= dist;

sec_idx
= j;

}

}

if(min_dist
<= (unsigned int)(sec_dist
* 0.8) && min_dist <=50){

//if(min_dist
<= 50){

matches.push_back(DMatch(i,
min_idx, 0, (float)min_dist));

}

}

}

int

main(int

argc, char*
argv[])

{

Mat
img1 = imread(image_filename1, 0);

Mat
img2 = imread(image_filename2, 0);

//GaussianBlur(img1,
img1, Size(5, 5), 0);

//GaussianBlur(img2,
img2, Size(5, 5), 0);

ORB
orb1(3000, ORB::CommonParams(1.2, 8));

ORB
orb2(100, ORB::CommonParams(1.2, 1));

vector
keys1, keys2;

Mat
descriptors1, descriptors2;

orb1(img1,
Mat(), keys1, descriptors1, false);

printf("tem
feat num: %d\n",
keys1.size());

int64
st, et;

st
= cvGetTickCount();

orb2(img2,
Mat(), keys2, descriptors2, false);

et
= cvGetTickCount();

printf("orb2
extraction time: %f\n",
(et-st)/(double)cvGetTickFrequency()/1000.);

printf("query
feat num: %d\n",
keys2.size());

//
find matches

vector
matches;

st
= cvGetTickCount();

//for(int
i = 0; i < 10; i++){

naive_nn_search2(keys1,
descriptors1, keys2, descriptors2, matches);

//}

et
= cvGetTickCount();

printf("match
time: %f\n",
(et-st)/(double)cvGetTickFrequency()/1000.);

printf("matchs
num: %d\n",
matches.size());

Mat
showImg;

drawMatches(img2,
keys2, img1, keys1, matches, showImg, CV_RGB(0, 255, 0), CV_RGB(0, 0, 255));

string
winName = "Matches";

namedWindow(
winName, WINDOW_AUTOSIZE );

imshow(
winName, showImg );

waitKey();

vector

pt1;

vector

pt2;

for(int

i = 0; i < (int)matches.size();
i++){

pt1.push_back(Point2f(keys2[matches[i].queryIdx].pt.x,
keys2[matches[i].queryIdx].pt.y));

pt2.push_back(Point2f(keys1[matches[i].trainIdx].pt.x,
keys1[matches[i].trainIdx].pt.y));

}

Mat
homo;

st
= cvGetTickCount();

homo
= findHomography(pt1, pt2, Mat(), CV_RANSAC, 5);

et
= cvGetTickCount();

printf("ransac
time: %f\n",
(et-st)/(double)cvGetTickFrequency()/1000.);

printf("homo\n"

"%f
%f %f\n"

"%f
%f %f\n"

"%f
%f %f\n",

homo.at(0,0),
homo.at(0,1), homo.at(0,2),

homo.at(1,0),
homo.at(1,1), homo.at(1,2),

homo.at(2,0),homo.at(2,1),homo.at(2,2));

vector

reproj;

reproj.resize(pt1.size());

perspectiveTransform(pt1,
reproj, homo);

Mat
diff;

diff
= Mat(reproj) - Mat(pt2);

int

inlier = 0;

double

err_sum = 0;

for(int

i = 0; i < diff.rows; i++){

float*
ptr = diff.ptr(i);

float

err = ptr[0]*ptr[0] + ptr[1]*ptr[1];

if(err
< 25.f){

inlier++;

err_sum
+= sqrt(err);

}

}

printf("inlier
num: %d\n",
inlier);

printf("ratio
%f\n",
inlier / (float)(diff.rows));

printf("mean
reprojection error: %f\n",
err_sum / inlier);

return

0;

}

0 0