1 | /* |
---|
2 | * Speeded-Up Robust Features (SURF) |
---|
3 | * http://people.ee.ethz.ch/~surf |
---|
4 | * |
---|
5 | * Authors: Herbert Bay, Andreas Ess, Geert Willems |
---|
6 | * Windows port by Stefan Saur |
---|
7 | * |
---|
8 | * Copyright (2006): ETH Zurich, Switzerland |
---|
9 | * Katholieke Universiteit Leuven, Belgium |
---|
10 | * All rights reserved. |
---|
11 | * |
---|
12 | * For details, see the paper: |
---|
13 | * Herbert Bay, Tinne Tuytelaars, Luc Van Gool, |
---|
14 | * "SURF: Speeded Up Robust Features" |
---|
15 | * Proceedings of the ninth European Conference on Computer Vision, May 2006 |
---|
16 | * |
---|
17 | * Permission to use, copy, modify, and distribute this software and |
---|
18 | * its documentation for educational, research, and non-commercial |
---|
19 | * purposes, without fee and without a signed licensing agreement, is |
---|
20 | * hereby granted, provided that the above copyright notice and this |
---|
21 | * paragraph appear in all copies modifications, and distributions. |
---|
22 | * |
---|
23 | * Any commercial use or any redistribution of this software |
---|
24 | * requires a license from one of the above mentioned establishments. |
---|
25 | * |
---|
26 | * For further details, contact Andreas Ess (aess@vision.ee.ethz.ch). |
---|
27 | */ |
---|
28 | |
---|
29 | #include <iostream> |
---|
30 | #include <fstream> |
---|
31 | #include <vector> |
---|
32 | #include <cmath> |
---|
33 | #include <string.h> |
---|
34 | |
---|
35 | #ifdef WIN32 |
---|
36 | #include "surfWINDLL.h" |
---|
37 | #endif |
---|
38 | |
---|
39 | #include "imload.h" |
---|
40 | #include "surflib.h" |
---|
41 | |
---|
42 | #include "os_mapping.h" |
---|
43 | |
---|
44 | using namespace std; |
---|
45 | using namespace surf; |
---|
46 | |
---|
47 | // Length of the descriptor vector |
---|
48 | int VLength; |
---|
49 | |
---|
50 | // Forward declaration of the functions to load/save the SURF points |
---|
51 | void loadIpoints(string fn, vector< Ipoint >& keys, bool bVerbose = false); |
---|
52 | void saveIpoints(string fn, const vector< Ipoint >& keys, bool bVerbose = false, bool bLaplacian = true); |
---|
53 | |
---|
54 | int main (int argc, char **argv) |
---|
55 | { |
---|
56 | // Initial sampling step (default 2) |
---|
57 | int samplingStep = 2; |
---|
58 | // Number of analysed octaves (default 4) |
---|
59 | int octaves = 4; |
---|
60 | // Blob response treshold |
---|
61 | double thres = 4.0; |
---|
62 | // Set this flag "true" to double the image size |
---|
63 | bool doubleImageSize = false; |
---|
64 | // Initial lobe size, default 3 and 5 (with double image size) |
---|
65 | int initLobe = 3; |
---|
66 | // Upright SURF or rotation invaraiant |
---|
67 | bool upright = false; |
---|
68 | // If the extended flag is turned on, SURF 128 is used |
---|
69 | bool extended = false; |
---|
70 | // Spatial size of the descriptor window (default 4) |
---|
71 | int indexSize = 4; |
---|
72 | // Variables for the timing measure |
---|
73 | osmapping::os_TIME tim1, tim2; //STS |
---|
74 | // verbose output |
---|
75 | bool bVerbose = true; |
---|
76 | // skip sign of laplacian |
---|
77 | bool bLaplacian = true; |
---|
78 | |
---|
79 | bool bLoadRegions = false; |
---|
80 | string sRegionFile = ""; |
---|
81 | |
---|
82 | // Print command line help |
---|
83 | if (argc==1) { |
---|
84 | cerr << "./surf -i img.pgm -o img.surf [options]\n" |
---|
85 | << " blob response threshold -thres 1000\n" |
---|
86 | << " double image size: -d\n" |
---|
87 | << " custom lobe size: -ms 3\n" |
---|
88 | << " initial sampling step: -ss 2\n" |
---|
89 | << " number of octaves: -oc 4\n" |
---|
90 | << " U-SURF (not rotation invariant): -u\n" |
---|
91 | << " extended descriptor (SURF-128): -e\n" |
---|
92 | << " descriptor size: -in 4\n" |
---|
93 | << " input regions: -p1 <file>\n" |
---|
94 | << " verbose output: -v\n" |
---|
95 | << " don't write laplacian: -nl\n" |
---|
96 | << " quiet mode: -q\n"; |
---|
97 | return(0); |
---|
98 | } |
---|
99 | |
---|
100 | // Read the arguments |
---|
101 | ImLoad ImageLoader; |
---|
102 | int arg = 0; |
---|
103 | string fn = "out.surf"; |
---|
104 | Image *im=NULL; |
---|
105 | while (++arg < argc) { |
---|
106 | if (! strcmp(argv[arg], "-i")) |
---|
107 | im = ImageLoader.readImage(argv[++arg]); |
---|
108 | if (! strcmp(argv[arg], "-o")) |
---|
109 | fn = argv[++arg]; |
---|
110 | if (! strcmp(argv[arg], "-thres")) |
---|
111 | thres = (atof(argv[++arg]))/10000; |
---|
112 | if (! strcmp(argv[arg], "-d")) |
---|
113 | doubleImageSize = true; |
---|
114 | if (! strcmp(argv[arg], "-ms")) |
---|
115 | initLobe = atoi(argv[++arg]); |
---|
116 | if (! strcmp(argv[arg], "-oc")) |
---|
117 | octaves = atoi(argv[++arg]); |
---|
118 | if (! strcmp(argv[arg], "-ss")) |
---|
119 | samplingStep = atoi(argv[++arg]); |
---|
120 | if (! strcmp(argv[arg], "-u")) |
---|
121 | upright = true; |
---|
122 | if (! strcmp(argv[arg], "-e")) |
---|
123 | extended = true; |
---|
124 | if (! strcmp(argv[arg], "-in")) |
---|
125 | indexSize = atoi(argv[++arg]); |
---|
126 | if (! strcmp(argv[arg], "-p1")) { |
---|
127 | bLoadRegions = true; |
---|
128 | sRegionFile = argv[++arg]; |
---|
129 | } |
---|
130 | if (! strcmp(argv[arg], "-v")) |
---|
131 | bVerbose = true; |
---|
132 | if (! strcmp(argv[arg], "-nl")) |
---|
133 | bLaplacian = false; |
---|
134 | if (! strcmp(argv[arg], "-q")) |
---|
135 | bVerbose = false; |
---|
136 | } |
---|
137 | |
---|
138 | |
---|
139 | // Start measuring the time |
---|
140 | osmapping::os_GetTime(&tim1); |
---|
141 | |
---|
142 | // Create the integral image |
---|
143 | Image iimage(im, doubleImageSize); |
---|
144 | |
---|
145 | // Start finding the SURF points |
---|
146 | if( bVerbose ) |
---|
147 | cout << "Finding SURFs...\n"; |
---|
148 | |
---|
149 | // These are the interest points |
---|
150 | vector< Ipoint > ipts; |
---|
151 | ipts.reserve(1000); |
---|
152 | |
---|
153 | // Extract interest points with Fast-Hessian |
---|
154 | FastHessian fh(&iimage, /* pointer to integral image */ |
---|
155 | ipts, |
---|
156 | thres, /* blob response threshold */ |
---|
157 | doubleImageSize, /* double image size flag */ |
---|
158 | initLobe * 3 /* 3 times lobe size equals the mask size */, |
---|
159 | samplingStep, /* subsample the blob response map */ |
---|
160 | octaves /* number of octaves to be analysed */); |
---|
161 | |
---|
162 | if( bLoadRegions ) { |
---|
163 | // Load the interest points from disk |
---|
164 | loadIpoints( sRegionFile, ipts, bVerbose ); |
---|
165 | } else { |
---|
166 | // Extract them and get their pointer |
---|
167 | fh.getInterestPoints(); |
---|
168 | } |
---|
169 | |
---|
170 | // Initialise the SURF descriptor |
---|
171 | Surf des(&iimage, /* pointer to integral image */ |
---|
172 | doubleImageSize, /* double image size flag */ |
---|
173 | upright, /* rotation invariance or upright */ |
---|
174 | extended, /* use the extended descriptor */ |
---|
175 | indexSize /* square size of the descriptor window (default 4x4)*/); |
---|
176 | |
---|
177 | // Get the length of the descriptor vector resulting from the parameters |
---|
178 | VLength = des.getVectLength(); |
---|
179 | |
---|
180 | if(upright) |
---|
181 | cout << "Using upright = " << upright << ", with descriptor length " << VLength << ", threshold = " << thres << endl; |
---|
182 | // If the extended flag is turned on, SURF 128 is used |
---|
183 | |
---|
184 | // Compute the orientation and the descriptor for every interest point |
---|
185 | for (unsigned n=0; n<ipts.size(); n++){ |
---|
186 | //for (Ipoint *k = ipts; k != NULL; k = k->next){ |
---|
187 | // set the current interest point |
---|
188 | des.setIpoint(&ipts[n]); |
---|
189 | // assign reproducible orientation |
---|
190 | des.assignOrientation(); |
---|
191 | // make the SURF descriptor |
---|
192 | des.makeDescriptor(); |
---|
193 | } |
---|
194 | |
---|
195 | // stop measuring the time, we're all done |
---|
196 | osmapping::os_GetTime(&tim2); |
---|
197 | |
---|
198 | // save the interest points in the output file |
---|
199 | saveIpoints(fn, ipts, bVerbose, bLaplacian); |
---|
200 | |
---|
201 | // print some nice information on the command prompt |
---|
202 | if( bVerbose ) |
---|
203 | cout << "Detection time: " << osmapping::os_TimeDiff(&tim2, &tim1) << " ms" << endl; |
---|
204 | |
---|
205 | delete im; |
---|
206 | |
---|
207 | return 0; |
---|
208 | } |
---|
209 | |
---|
210 | // Save the interest points to a regular ASCII file |
---|
211 | void saveIpoints(string sFileName, const vector< Ipoint >& ipts, bool bVerbose, bool bLaplacian) |
---|
212 | { |
---|
213 | ofstream ipfile(sFileName.c_str()); |
---|
214 | if( !ipfile ) { |
---|
215 | cerr << "ERROR in loadIpoints(): " |
---|
216 | << "Couldn't open file '" << sFileName.c_str() << "'!" << endl; //STS |
---|
217 | return; |
---|
218 | } |
---|
219 | |
---|
220 | double sc; |
---|
221 | unsigned count = ipts.size(); |
---|
222 | |
---|
223 | // Write the file header |
---|
224 | if (bLaplacian) |
---|
225 | ipfile << VLength + 1 << endl << count << endl; |
---|
226 | else |
---|
227 | ipfile << VLength << endl << count << endl; |
---|
228 | |
---|
229 | // In order to just save the interest points without descriptor, comment |
---|
230 | // the above and uncomment the following command. |
---|
231 | // ipfile << 1.0 << endl << count << endl; |
---|
232 | // Save interest point with descriptor in the format of Krystian Mikolajczyk |
---|
233 | // for reasons of comparison with other descriptors. As our interest points |
---|
234 | // are circular in any case, we use the second component of the ellipse to |
---|
235 | // provide some information about the strength of the interest point. This is |
---|
236 | // important for 3D reconstruction as only the strongest interest points are |
---|
237 | // considered. Replace the strength with 0.0 in order to perform Krystian's |
---|
238 | // comparisons. |
---|
239 | for (unsigned n=0; n<ipts.size(); n++){ |
---|
240 | // circular regions with diameter 5 x scale |
---|
241 | sc = 2.5 * ipts[n].scale; sc*=sc; |
---|
242 | ipfile << ipts[n].x /* x-location of the interest point */ |
---|
243 | << " " << ipts[n].y /* y-location of the interest point */ |
---|
244 | << " " << 1.0/sc /* 1/r^2 */ |
---|
245 | << " " << 0.0 //(*ipts)[n]->strength /* 0.0 */ |
---|
246 | << " " << 1.0/sc; /* 1/r^2 */ |
---|
247 | |
---|
248 | if (bLaplacian) |
---|
249 | ipfile << " " << ipts[n].laplace; |
---|
250 | |
---|
251 | // Here comes the descriptor |
---|
252 | for (int i = 0; i < VLength; i++) { |
---|
253 | ipfile << " " << ipts[n].ivec[i]; |
---|
254 | } |
---|
255 | ipfile << endl; |
---|
256 | } |
---|
257 | |
---|
258 | // Write message to terminal. |
---|
259 | if( bVerbose ) |
---|
260 | cout << count << " interest points found" << endl; |
---|
261 | } |
---|
262 | |
---|
263 | |
---|
264 | // Load the interest points from a regular ASCII file |
---|
265 | void loadIpoints(string sFileName, vector< Ipoint >& ipts, bool bVerbose) |
---|
266 | { |
---|
267 | ifstream ipfile(sFileName.c_str()); |
---|
268 | if( !ipfile ) { |
---|
269 | cerr << "ERROR in loadIpoints(): " |
---|
270 | << "Couldn't open file '" << sFileName.c_str() << "'!" << endl; //STS |
---|
271 | return; |
---|
272 | } |
---|
273 | |
---|
274 | // Load the file header |
---|
275 | float dummy; |
---|
276 | unsigned count; |
---|
277 | ipfile >> dummy >> count; |
---|
278 | |
---|
279 | // create a new interest point vector |
---|
280 | ipts.clear(); |
---|
281 | ipts.resize(count); |
---|
282 | |
---|
283 | // Load the interest points in Mikolajczyk's format |
---|
284 | for (unsigned n=0; n<count; n++){ |
---|
285 | // circular regions with diameter 5 x scale |
---|
286 | float x, y, a, b, c; |
---|
287 | ipfile >> x >> y >> a >> b >> c; |
---|
288 | |
---|
289 | float det = sqrt((a-c)*(a-c) + 4.0*b*b); |
---|
290 | float e1 = 0.5*(a+c + det); |
---|
291 | float e2 = 0.5*(a+c - det); |
---|
292 | float l1 = (1.0/sqrt(e1)); |
---|
293 | float l2 = (1.0/sqrt(e2)); |
---|
294 | float sc = sqrt( l1*l2 ); |
---|
295 | |
---|
296 | ipts[n].x = x; |
---|
297 | ipts[n].y = y; |
---|
298 | ipts[n].scale = sc/2.5; |
---|
299 | } |
---|
300 | |
---|
301 | // close the interest point file again |
---|
302 | ipfile.close(); |
---|
303 | |
---|
304 | // Write message to terminal. |
---|
305 | if( bVerbose ) |
---|
306 | cout << "read in " << count << " interest points." << endl; |
---|
307 | } |
---|